Context Navigation

trazdf_imp_tam.F90 @ 2587

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

Line
1	MODULE trazdf_imp_tam
2	#ifdef key_tam
3	!!==============================================================================
4	!! * MODULE trazdf_imp_tam *
5	!! Ocean active tracers: vertical component of the tracer mixing trend
6	!! Tangent and Adjoint Module
7	!!==============================================================================
8	!! History of the direct module:
9	!! 6.0 ! 90-10 (B. Blanke) Original code
10	!! 7.0 ! 91-11 (G. Madec)
11	!! ! 92-06 (M. Imbard) correction on tracer trend loops
12	!! ! 96-01 (G. Madec) statement function for e3
13	!! ! 97-05 (G. Madec) vertical component of isopycnal
14	!! ! 97-07 (G. Madec) geopotential diffusion in s-coord
15	!! ! 00-08 (G. Madec) double diffusive mixing
16	!! 8.5 ! 02-08 (G. Madec) F90: Free form and module
17	!! 9.0 ! 06-11 (G. Madec) New step reorganisation
18	!! History of the T&A module:
19	!! ! 09-01 (A. Vidard) tam of the 06-11 version
20	!!----------------------------------------------------------------------
21	!! tra_zdf_imp_tan : Update the tracer trend with the diagonal vertical
22	!! part of the mixing tensor (tangent).
23	!! tra_zdf_imp_adj : Update the tracer trend with the diagonal vertical
24	!! part of the mixing tensor (adjoint).
25	!!----------------------------------------------------------------------
26	!! * Modules used
27	USE par_kind , ONLY: & ! Precision variables
28	& wp
29	USE par_oce , ONLY: & ! Ocean space and time domain variables
30	& jpi, &
31	& jpj, &
32	& jpk, &
33	& jpim1, &
34	& jpjm1, &
35	& jpkm1, &
36	& jpiglo
37	USE oce_tam , ONLY: & ! ocean dynamics and active tracers
38	& tb_tl, &
39	& sb_tl, &
40	& ta_tl, &
41	& sa_tl, &
42	& tb_ad, &
43	& sb_ad, &
44	& ta_ad, &
45	& sa_ad
46	USE dom_oce , ONLY: & ! ocean space and time domain
47	& e1t, &
48	& e2t, &
49	# if defined key_vvl
50	& e3t_1, &
51	# else
52	# if defined key_zco
53	& e3t_0, &
54	& e3w_0, &
55	# else
56	& e3t, &
57	& e3w, &
58	# endif
59	# endif
60	& tmask, &
61	& lk_vvl, &
62	& mig, &
63	& mjg, &
64	& nldi, &
65	& nldj, &
66	& nlei, &
67	& nlej, &
68	& rdttra
69	USE oce , ONLY: & ! ocean dynamics and tracers variables
70	& l_traldf_rot
71	USE zdf_oce , ONLY: & ! ocean vertical physics
72	& avt
73	USE ldftra_oce , ONLY: & ! ocean active tracers: lateral physics
74	& ahtw, &
75	& aht0
76	#if defined key_ldfslp
77	USE ldfslp , ONLY: & ! lateral physics: slope of diffusion
78	& wslpi, & !: i_slope at W-points
79	& wslpj !: j-slope at W-points
80	#endif
81	#if defined key_zdfddm
82	USE zdfddm , ONLY: &
83	& avs
84	#endif
85	USE traldf_tam
86	USE in_out_manager, ONLY: & ! I/O manager
87	& lwp, &
88	& numout, &
89	& nitend, &
90	& nit000
91	USE gridrandom , ONLY: & ! Random Gaussian noise on grids
92	& grid_random
93	USE dotprodfld , ONLY: & ! Computes dot product for 3D and 2D fields
94	& dot_product
95	USE tstool_tam , ONLY: &
96	& prntst_adj, & !
97	& prntst_tlm, & !
98	& stdt, & ! stdev for u-velocity
99	& stds ! v-velocity
100	USE paresp , ONLY: & ! Weights for an energy-type scalar product
101	& wesp_t, &
102	& wesp_s
103
104	IMPLICIT NONE
105	PRIVATE
106
107	!! * Routine accessibility
108	PUBLIC tra_zdf_imp_tan ! routine called by tra_zdf_tan.F90
109	PUBLIC tra_zdf_imp_adj ! routine called by tra_zdf_adj.F90
110	PUBLIC tra_zdf_imp_adj_tst ! routine called by tst.F90
111	#if defined key_tst_tlm
112	PUBLIC tra_zdf_imp_tlm_tst ! routine called by tamtst.F90
113	#endif
114
115	!! * Substitutions
116	# include "domzgr_substitute.h90"
117	# include "ldftra_substitute.h90"
118	# include "zdfddm_substitute.h90"
119	# include "vectopt_loop_substitute.h90"
120	!!----------------------------------------------------------------------
121	!!----------------------------------------------------------------------
122	!! OPA 9.0 , LOCEAN-IPSL (2005)
123	!! $Id: trazdf_imp.F90 1156 2008-06-26 16:06:45Z rblod $
124	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
125	!!----------------------------------------------------------------------
126	CONTAINS
127
128	SUBROUTINE tra_zdf_imp_tan( kt, p2dt )
129	!!----------------------------------------------------------------------
130	!! * ROUTINE tra_zdf_imp_tan *
131	!!
132	!! ** Purpose of the direct routine:
133	!! Compute the trend due to the vertical tracer diffusion
134	!! including the vertical component of lateral mixing (only for 2nd
135	!! order operator, for fourth order it is already computed and add
136	!! to the general trend in traldf.F) and add it to the general trend
137	!! of the tracer equations.
138	!!
139	!! ** Method of the direct routine :
140	!! The vertical component of the lateral diffusive trends
141	!! is provided by a 2nd order operator rotated along neutral or geo-
142	!! potential surfaces to which an eddy induced advection can be
143	!! added. It is computed using before fields (forward in time) and
144	!! isopycnal or geopotential slopes computed in routine ldfslp.
145	!!
146	!! Second part: vertical trend associated with the vertical physics
147	!! =========== (including the vertical flux proportional to dk[t]
148	!! associated with the lateral mixing, through the
149	!! update of avt)
150	!! The vertical diffusion of tracers (t & s) is given by:
151	!! difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) )
152	!! It is computed using a backward time scheme (t=ta).
153	!! Surface and bottom boundary conditions: no diffusive flux on
154	!! both tracers (bottom, applied through the masked field avt).
155	!! Add this trend to the general trend ta,sa :
156	!! ta = ta + dz( avt dz(t) )
157	!! (sa = sa + dz( avs dz(t) ) if lk_zdfddm=T )
158	!!
159	!! Third part: recover avt resulting from the vertical physics
160	!! ========== alone, for further diagnostics (for example to
161	!! compute the turbocline depth in zdfmxl.F90).
162	!! avt = zavt
163	!! (avs = zavs if lk_zdfddm=T )
164	!!
165	!! ** Remarks on the tangent routine : - key_vvl is not available in tangent yet.
166	!! Once it will be this routine wil need to be rewritten
167	!! - simplified version, slopes (wslp[ij])
168	!! assumed to be constant (read from the trajectory). same for av[ts]
169	!!
170	!!---------------------------------------------------------------------
171	!! * Modules used
172	USE oce , ONLY : zwd => ua, & ! ua used as workspace
173	zws => va ! va " "
174	!! * Arguments
175	INTEGER, INTENT( in ) :: kt ! ocean time-step index
176	REAL(wp), DIMENSION(jpk), INTENT( in ) :: &
177	p2dt ! vertical profile of tracer time-step
178
179	!! * Local declarations
180	INTEGER :: ji, jj, jk ! dummy loop indices
181	REAL(wp) :: zavi, zrhstl, znvvl, & ! temporary scalars
182	ze3tb, ze3tn, ze3ta, zvsfvvl ! variable vertical scale factors
183	REAL(wp), DIMENSION(jpi,jpj,jpk) :: &
184	zwi, zwt, zavsi ! workspace arrays
185	!!---------------------------------------------------------------------
186
187	IF( kt == nit000 ) THEN
188	IF(lwp)WRITE(numout,*)
189	IF(lwp)WRITE(numout,*) 'tra_zdf_imp : implicit vertical mixing'
190	IF(lwp)WRITE(numout,*) '~~~~~~~~~~~ '
191	zavi = 0._wp ! avoid warning at compilation phase when lk_ldfslp=F
192	ENDIF
193
194	! I. Local initialization
195	! -----------------------
196	zwd (1,:, : ) = 0._wp ; zwd (jpi,:,:) = 0._wp
197	zws (1,:, : ) = 0._wp ; zws (jpi,:,:) = 0._wp
198	zwi (1,:, : ) = 0._wp ; zwi (jpi,:,:) = 0._wp
199	zwt (1,:, : ) = 0._wp ; zwt (jpi,:,:) = 0._wp
200	zavsi(1,:, : ) = 0._wp ; zavsi(jpi,:,:) = 0._wp
201	zwt (:,:,jpk) = 0._wp ; zwt ( : ,:,1) = 0._wp
202	zavsi(:,:,jpk) = 0._wp ; zavsi( : ,:,1) = 0._wp
203
204	! I.1 Variable volume : to take into account vertical variable vertical scale factors
205	! -------------------
206	! ... not available in tangent yet
207	! II. Vertical trend associated with the vertical physics
208	! =======================================================
209	! (including the vertical flux proportional to dk[t] associated
210	! with the lateral mixing, through the avt update)
211	! dk[ avt dk[ (t,s) ] ] diffusive trends
212
213	! II.0 Matrix construction
214	! ------------------------
215	#if defined key_ldfslp
216	! update and save of avt (and avs if double diffusive mixing)
217	IF( l_traldf_rot ) THEN
218	DO jk = 2, jpkm1
219	DO jj = 2, jpjm1
220	DO ji = fs_2, fs_jpim1 ! vector opt.
221	zavi = fsahtw(ji,jj,jk) & ! vertical mixing coef. due to lateral mixing
222	& * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk) &
223	& + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) )
224	zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi ! zwt=avt+zavi (total vertical mixing coef. on temperature)
225	# if defined key_zdfddm
226	zavsi(ji,jj,jk) = fsavs(ji,jj,jk) + zavi ! dd mixing: zavsi = total vertical mixing coef. on salinity
227	# endif
228	END DO
229	END DO
230	END DO
231	ENDIF
232	#else
233	! No isopycnal diffusion
234	zwt(:,:,:) = avt(:,:,:)
235	# if defined key_zdfddm
236	zavsi(:,:,:) = avs(:,:,:)
237	# endif
238
239	#endif
240
241	! Diagonal, inferior, superior (including the bottom boundary condition via avt masked)
242	DO jk = 1, jpkm1
243	DO jj = 2, jpjm1
244	DO ji = fs_2, fs_jpim1 ! vector opt.
245	ze3ta = 1._wp ! after scale factor at T-point
246	ze3tn = fse3t(ji,jj,jk) ! now scale factor at T-point
247	zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk ) / ( ze3tn * fse3w(ji,jj,jk ) )
248	zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
249	zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
250	END DO
251	END DO
252	END DO
253
254	! Surface boudary conditions
255	DO jj = 2, jpjm1
256	DO ji = fs_2, fs_jpim1 ! vector opt.
257	ze3ta = 1._wp ! after scale factor at T-point
258	zwi(ji,jj,1) = 0._wp
259	zwd(ji,jj,1) = ze3ta - zws(ji,jj,1)
260	END DO
261	END DO
262
263	! II.1. Vertical diffusion on t
264	! ---------------------------
265
266	!! Matrix inversion from the first level
267	!!----------------------------------------------------------------------
268	! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
269	!
270	! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
271	! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
272	! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
273	! ( ... )( ... ) ( ... )
274	! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
275	!
276	! m is decomposed in the product of an upper and lower triangular matrix
277	! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
278	! The second member is in 2d array zwy
279	! The solution is in 2d array zwx
280	! The 3d arry zwt is a work space array
281	! zwy is used and then used as a work space array : its value is modified!
282
283	! first recurrence: Tk = Dk - Ik Sk-1 / Tk-1 (increasing k)
284	DO jj = 2, jpjm1
285	DO ji = fs_2, fs_jpim1
286	zwt(ji,jj,1) = zwd(ji,jj,1)
287	END DO
288	END DO
289	DO jk = 2, jpkm1
290	DO jj = 2, jpjm1
291	DO ji = fs_2, fs_jpim1
292	zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1)
293	END DO
294	END DO
295	END DO
296
297	! second recurrence: Zk = Yk - Ik / Tk-1 Zk-1
298	DO jj = 2, jpjm1
299	DO ji = fs_2, fs_jpim1
300	ze3tb = 1._wp
301	ze3tn = 1._wp
302	ta_tl(ji,jj,1) = ze3tb * tb_tl(ji,jj,1) + p2dt(1) * ze3tn * ta_tl(ji,jj,1)
303	END DO
304	END DO
305
306	DO jk = 2, jpkm1
307	DO jj = 2, jpjm1
308	DO ji = fs_2, fs_jpim1
309	ze3tb = 1._wp
310	ze3tn = 1._wp
311	zrhstl = ze3tb * tb_tl(ji,jj,jk) + p2dt(jk) * ze3tn * ta_tl(ji,jj,jk) ! zrhs=right hand side
312	ta_tl(ji,jj,jk) = zrhstl - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *ta_tl(ji,jj,jk-1)
313	END DO
314	END DO
315	END DO
316
317	! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
318	! Save the masked temperature after in ta
319	! (c a u t i o n: temperature not its trend, Leap-frog scheme done it will not be done in tranxt)
320	DO jj = 2, jpjm1
321	DO ji = fs_2, fs_jpim1
322	ta_tl(ji,jj,jpkm1) = ta_tl(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
323	END DO
324	END DO
325	DO jk = jpk-2, 1, -1
326	DO jj = 2, jpjm1
327	DO ji = fs_2, fs_jpim1
328	ta_tl(ji,jj,jk) = ( ta_tl(ji,jj,jk) - zws(ji,jj,jk) * ta_tl(ji,jj,jk+1) ) / zwt(ji,jj,jk) * tmask(ji,jj,jk)
329	END DO
330	END DO
331	END DO
332
333	! II.2 Vertical diffusion on salinity
334	! -----------------------------------
335
336	#if defined key_zdfddm
337	! Rebuild the Matrix as avt /= avs
338
339	! Diagonal, inferior, superior (including the bottom boundary condition via avs masked)
340	DO jk = 1, jpkm1
341	DO jj = 2, jpjm1
342	DO ji = fs_2, fs_jpim1 ! vector opt.
343	ze3ta = 1._wp ! after scale factor at T-point
344	ze3tn = fse3t(ji,jj,jk) ! now scale factor at T-point
345	zwi(ji,jj,jk) = - p2dt(jk) * zavsi(ji,jj,jk ) / ( ze3tn * fse3w(ji,jj,jk ) )
346	zws(ji,jj,jk) = - p2dt(jk) * zavsi(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
347	zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
348	END DO
349	END DO
350	END DO
351
352	! Surface boudary conditions
353	DO jj = 2, jpjm1
354	DO ji = fs_2, fs_jpim1 ! vector opt.
355	ze3ta = 1._wp ! after scale factor at T-point
356	zwi(ji,jj,1) = 0._wp
357	zwd(ji,jj,1) = ze3ta - zws(ji,jj,1)
358	END DO
359	END DO
360	#endif
361
362
363	!! Matrix inversion from the first level
364	!!----------------------------------------------------------------------
365	! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
366	!
367	! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
368	! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
369	! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
370	! ( ... )( ... ) ( ... )
371	! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
372	!
373	! m is decomposed in the product of an upper and lower triangular
374	! matrix
375	! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
376	! The second member is in 2d array zwy
377	! The solution is in 2d array zwx
378	! The 3d arry zwt is a work space array
379	! zwy is used and then used as a work space array : its value is modified!
380
381	! first recurrence: Tk = Dk - Ik Sk-1 / Tk-1 (increasing k)
382	DO jj = 2, jpjm1
383	DO ji = fs_2, fs_jpim1
384	zwt(ji,jj,1) = zwd(ji,jj,1)
385	END DO
386	END DO
387	DO jk = 2, jpkm1
388	DO jj = 2, jpjm1
389	DO ji = fs_2, fs_jpim1
390	zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1)
391	END DO
392	END DO
393	END DO
394
395	! second recurrence: Zk = Yk - Ik / Tk-1 Zk-1
396	DO jj = 2, jpjm1
397	DO ji = fs_2, fs_jpim1
398	ze3tb = 1.0_wp ! before scale factor at T-point
399	ze3tn = 1.0_wp ! now scale factor at T-point
400	sa_tl(ji,jj,1) = ze3tb * sb_tl(ji,jj,1) + p2dt(1) * ze3tn * sa_tl(ji,jj,1)
401	END DO
402	END DO
403	DO jk = 2, jpkm1
404	DO jj = 2, jpjm1
405	DO ji = fs_2, fs_jpim1
406	ze3tb = 1.0_wp ! before scale factor at T-point
407	ze3tn = 1.0_wp ! now scale factor at T-point
408	zrhstl = ze3tb * sb_tl(ji,jj,jk) + p2dt(jk) * ze3tn * sa_tl(ji,jj,jk) ! zrhs=right hand side
409	sa_tl(ji,jj,jk) = zrhstl - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *sa_tl(ji,jj,jk-1)
410	END DO
411	END DO
412	END DO
413
414	! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
415	! Save the masked temperature after in ta
416	! (c a u t i o n: temperature not its trend, Leap-frog scheme done it will not be done in tranxt)
417	DO jj = 2, jpjm1
418	DO ji = fs_2, fs_jpim1
419	sa_tl(ji,jj,jpkm1) = sa_tl(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
420	END DO
421	END DO
422	DO jk = jpk-2, 1, -1
423	DO jj = 2, jpjm1
424	DO ji = fs_2, fs_jpim1
425	sa_tl(ji,jj,jk) = ( sa_tl(ji,jj,jk) - zws(ji,jj,jk) * sa_tl(ji,jj,jk+1) ) / zwt(ji,jj,jk) * tmask(ji,jj,jk)
426	END DO
427	END DO
428	END DO
429
430	END SUBROUTINE tra_zdf_imp_tan
431	SUBROUTINE tra_zdf_imp_adj( kt, p2dt )
432	!!----------------------------------------------------------------------
433	!! * ROUTINE tra_zdf_imp_adj *
434	!!
435	!! ** Purpose of the direct routine:
436	!! Compute the trend due to the vertical tracer diffusion
437	!! including the vertical component of lateral mixing (only for 2nd
438	!! order operator, for fourth order it is already computed and add
439	!! to the general trend in traldf.F) and add it to the general trend
440	!! of the tracer equations.
441	!!
442	!! ** Method of the direct routine :
443	!! The vertical component of the lateral diffusive trends
444	!! is provided by a 2nd order operator rotated along neutral or geo-
445	!! potential surfaces to which an eddy induced advection can be
446	!! added. It is computed using before fields (forward in time) and
447	!! isopycnal or geopotential slopes computed in routine ldfslp.
448	!!
449	!! Second part: vertical trend associated with the vertical physics
450	!! =========== (including the vertical flux proportional to dk[t]
451	!! associated with the lateral mixing, through the
452	!! update of avt)
453	!! The vertical diffusion of tracers (t & s) is given by:
454	!! difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) )
455	!! It is computed using a backward time scheme (t=ta).
456	!! Surface and bottom boundary conditions: no diffusive flux on
457	!! both tracers (bottom, applied through the masked field avt).
458	!! Add this trend to the general trend ta,sa :
459	!! ta = ta + dz( avt dz(t) )
460	!! (sa = sa + dz( avs dz(t) ) if lk_zdfddm=T )
461	!!
462	!! Third part: recover avt resulting from the vertical physics
463	!! ========== alone, for further diagnostics (for example to
464	!! compute the turbocline depth in zdfmxl.F90).
465	!! avt = zavt
466	!! (avs = zavs if lk_zdfddm=T )
467	!!
468	!! ** Remarks on the adjoint routine : - key_vvl is not available in adjoint yet.
469	!! Once it will be this routine wil need to be rewritten
470	!! - simplified version, slopes (wslp[ij])
471	!! assumed to be constant (read from the trajectory). same for av[ts]
472	!!
473	!!---------------------------------------------------------------------
474	!! * Modules used
475	USE oce , ONLY : zwd => ua, & ! ua used as workspace
476	zws => va ! va " "
477	!! * Arguments
478	INTEGER, INTENT( in ) :: kt ! ocean time-step index
479	REAL(wp), DIMENSION(jpk), INTENT( in ) :: &
480	p2dt ! vertical profile of tracer time-step
481
482	!! * Local declarations
483	INTEGER :: ji, jj, jk ! dummy loop indices
484	REAL(wp) :: zavi, zrhsad, znvvl, & ! temporary scalars
485	ze3tb, ze3tn, ze3ta, zvsfvvl ! variable vertical scale factors
486	REAL(wp), DIMENSION(jpi,jpj,jpk) :: &
487	zwi, zwt, zavsi ! workspace arrays
488	!!---------------------------------------------------------------------
489
490	IF( kt == nitend ) THEN
491	IF(lwp)WRITE(numout,*)
492	IF(lwp)WRITE(numout,*) 'tra_zdf_imp_adj : implicit vertical mixing'
493	IF(lwp)WRITE(numout,*) '~~~~~~~~~~~~~~~ '
494	CALL ldf_ctl_tam ! init of l_traldf_rot
495	zavi = 0._wp ! avoid warning at compilation phase when lk_ldfslp=F
496	ENDIF
497
498	! I. Local initialization
499	! -----------------------
500	zrhsad = 0.0_wp
501	zwd (1,:, : ) = 0._wp ; zwd (jpi,:,:) = 0._wp
502	zws (1,:, : ) = 0._wp ; zws (jpi,:,:) = 0._wp
503	zwi (1,:, : ) = 0._wp ; zwi (jpi,:,:) = 0._wp
504	zwt (1,:, : ) = 0._wp ; zwt (jpi,:,:) = 0._wp
505	zavsi(1,:, : ) = 0._wp ; zavsi(jpi,:,:) = 0._wp
506	zwt (:,:,jpk) = 0._wp ; zwt ( : ,:,1) = 0._wp
507	zavsi(:,:,jpk) = 0._wp ; zavsi( : ,:,1) = 0._wp
508
509	! I.1 Variable volume : to take into account vertical variable vertical scale factors
510	! -------------------
511	! ... not available in adjoint yet
512	! II. Vertical trend associated with the vertical physics
513	! =======================================================
514	! (including the vertical flux proportional to dk[t] associated
515	! with the lateral mixing, through the avt update)
516	! dk[ avt dk[ (t,s) ] ] diffusive trends
517
518
519	! II.0 Matrix construction
520	! ------------------------
521
522	#if defined key_ldfslp
523	! update and save of avt (and avs if double diffusive mixing)
524	IF( l_traldf_rot ) THEN
525	DO jk = 2, jpkm1
526	DO jj = 2, jpjm1
527	DO ji = fs_2, fs_jpim1 ! vector opt.
528	zavi = fsahtw(ji,jj,jk) & ! vertical mixing coef. due to lateral mixing
529	& * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk) &
530	& + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) )
531	zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi ! zwt=avt+zavi (total vertical mixing coef. on temperature)
532	# if defined key_zdfddm
533	zavsi(ji,jj,jk) = fsavs(ji,jj,jk) + zavi ! dd mixing: zavsi = total vertical mixing coef. on salinity
534	# endif
535	END DO
536	END DO
537	END DO
538	ENDIF
539	#else
540	! No isopycnal diffusion
541	zwt(:,:,:) = avt(:,:,:)
542	# if defined key_zdfddm
543	zavsi(:,:,:) = avs(:,:,:)
544	# endif
545
546	#endif
547
548	! Diagonal, inferior, superior (including the bottom boundary condition via avt masked)
549	DO jk = 1, jpkm1
550	DO jj = 2, jpjm1
551	DO ji = fs_2, fs_jpim1 ! vector opt.
552	ze3ta = 1._wp ! after scale factor at T-point
553	ze3tn = fse3t(ji,jj,jk) ! now scale factor at T-point
554	zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk ) / ( ze3tn * fse3w(ji,jj,jk ) )
555	zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
556	zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
557	END DO
558	END DO
559	END DO
560
561	! Surface boudary conditions
562	DO jj = 2, jpjm1
563	DO ji = fs_2, fs_jpim1 ! vector opt.
564	ze3ta = 1._wp ! after scale factor at T-point
565	zwi(ji,jj,1) = 0._wp
566	zwd(ji,jj,1) = ze3ta - zws(ji,jj,1)
567	END DO
568	END DO
569
570
571	! II.1. Vertical diffusion on t
572	! ---------------------------
573
574	!! Matrix inversion from the first level
575	!!----------------------------------------------------------------------
576	! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
577	!
578	! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
579	! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
580	! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
581	! ( ... )( ... ) ( ... )
582	! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
583	!
584	! m is decomposed in the product of an upper and lower triangular matrix
585	! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
586	! The second member is in 2d array zwy
587	! The solution is in 2d array zwx
588	! The 3d arry zwt is a work space array
589	! zwy is used and then used as a work space array : its value is modified!
590
591	! first recurrence: Tk = Dk - Ik Sk-1 / Tk-1 (increasing k)
592	DO jj = 2, jpjm1
593	DO ji = fs_2, fs_jpim1
594	zwt(ji,jj,1) = zwd(ji,jj,1)
595	END DO
596	END DO
597	DO jk = 2, jpkm1
598	DO jj = 2, jpjm1
599	DO ji = fs_2, fs_jpim1
600	zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1)
601	END DO
602	END DO
603	END DO
604	! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
605	! Save the masked temperature after in ta
606	! (c a u t i o n: temperature not its trend, Leap-frog scheme done it will not be done in tranxt)
607	DO jk = 1, jpk-2
608	DO jj = 2, jpjm1
609	DO ji = fs_2, fs_jpim1
610	ta_ad(ji,jj,jk+1) = ta_ad(ji,jj,jk+1) - zws(ji,jj,jk) * ta_ad(ji,jj,jk) / zwt(ji,jj,jk) * tmask(ji,jj,jk)
611	ta_ad(ji,jj,jk) = ta_ad(ji,jj,jk) / zwt(ji,jj,jk) * tmask(ji,jj,jk)
612	END DO
613	END DO
614	END DO
615	DO jj = 2, jpjm1
616	DO ji = fs_2, fs_jpim1
617	ta_ad(ji,jj,jpkm1) = ta_ad(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
618	END DO
619	END DO
620	! second recurrence: Zk = Yk - Ik / Tk-1 Zk-1
621	DO jk = jpkm1, 2, -1
622	DO jj = 2, jpjm1
623	DO ji = fs_2, fs_jpim1
624	ze3tb = 1._wp
625	ze3tn = 1._wp
626	zrhsad = zrhsad + ta_ad(ji,jj,jk)
627	ta_ad(ji,jj,jk-1) = ta_ad(ji,jj,jk-1) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * ta_ad(ji,jj,jk)
628	ta_ad(ji,jj,jk) = 0.0_wp
629	tb_ad(ji,jj,jk) = tb_ad(ji,jj,jk) + ze3tb * zrhsad
630	ta_ad(ji,jj,jk) = ta_ad(ji,jj,jk) + p2dt(jk) * ze3tn * zrhsad
631	zrhsad = 0.0_wp
632	END DO
633	END DO
634	END DO
635	DO jj = 2, jpjm1
636	DO ji = fs_2, fs_jpim1
637	ze3tb = 1._wp
638	ze3tn = 1._wp
639	tb_ad(ji,jj,1) = tb_ad(ji,jj,1) + ze3tb * ta_ad(ji,jj,1)
640	ta_ad(ji,jj,1) = ta_ad(ji,jj,1) * p2dt(1) * ze3tn
641	END DO
642	END DO
643	! II.2 Vertical diffusion on salinity
644	! -----------------------------------
645
646	#if defined key_zdfddm
647	! Rebuild the Matrix as avt /= avs
648
649	! Diagonal, inferior, superior (including the bottom boundary condition via avs masked)
650	DO jk = jpkm1, 1, -1
651	DO jj = 2, jpjm1
652	DO ji = fs_2, fs_jpim1 ! vector opt.
653	ze3ta = 1._wp ! after scale factor at T-point
654	ze3tn = fse3t(ji,jj,jk) ! now scale factor at T-point
655	zwi(ji,jj,jk) = - p2dt(jk) * zavsi(ji,jj,jk ) / ( ze3tn * fse3w(ji,jj,jk ) )
656	zws(ji,jj,jk) = - p2dt(jk) * zavsi(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
657	zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
658	END DO
659	END DO
660	END DO
661
662	! Surface boudary conditions
663	DO jj = 2, jpjm1
664	DO ji = fs_2, fs_jpim1 ! vector opt.
665	ze3ta = 1._wp ! after scale factor at T-point
666	zwi(ji,jj,1) = 0._wp
667	zwd(ji,jj,1) = ze3ta - zws(ji,jj,1)
668	END DO
669	END DO
670	#endif
671
672
673	!! Matrix inversion from the first level
674	!!----------------------------------------------------------------------
675	! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
676	!
677	! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
678	! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
679	! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
680	! ( ... )( ... ) ( ... )
681	! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
682	!
683	! m is decomposed in the product of an upper and lower triangular
684	! matrix
685	! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
686	! The second member is in 2d array zwy
687	! The solution is in 2d array zwx
688	! The 3d arry zwt is a work space array
689	! zwy is used and then used as a work space array : its value is modified!
690
691	! first recurrence: Tk = Dk - Ik Sk-1 / Tk-1 (increasing k)
692	DO jj = 2, jpjm1
693	DO ji = fs_2, fs_jpim1
694	zwt(ji,jj,1) = zwd(ji,jj,1)
695	END DO
696	END DO
697	DO jk = 2, jpkm1
698	DO jj = 2, jpjm1
699	DO ji = fs_2, fs_jpim1
700	zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1)
701	END DO
702	END DO
703	END DO
704	! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
705	! Save the masked temperature after in ta
706	! (c a u t i o n: temperature not its trend, Leap-frog scheme done it will not be done in tranxt)
707	DO jk = 1, jpk-2
708	DO jj = 2, jpjm1
709	DO ji = fs_2, fs_jpim1
710	sa_ad(ji,jj,jk+1) = sa_ad(ji,jj,jk+1) - zws(ji,jj,jk) * sa_ad(ji,jj,jk) / zwt(ji,jj,jk) * tmask(ji,jj,jk)
711	sa_ad(ji,jj,jk ) = sa_ad(ji,jj,jk ) / zwt(ji,jj,jk) * tmask(ji,jj,jk)
712	END DO
713	END DO
714	END DO
715	DO jj = 2, jpjm1
716	DO ji = fs_2, fs_jpim1
717	sa_ad(ji,jj,jpkm1) = sa_ad(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
718	END DO
719	END DO
720
721	! second recurrence: Zk = Yk - Ik / Tk-1 Zk-1
722	DO jk = jpkm1, 2, -1
723	DO jj = 2, jpjm1
724	DO ji = fs_2, fs_jpim1
725	ze3tb = 1.0_wp ! before scale factor at T-point
726	ze3tn = 1.0_wp ! now scale factor at T-point
727	zrhsad = zrhsad + sa_ad(ji,jj,jk)
728	sa_ad(ji,jj,jk-1) = sa_ad(ji,jj,jk-1) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * sa_ad(ji,jj,jk)
729	sa_ad(ji,jj,jk) = 0.0_wp
730	sb_ad(ji,jj,jk) = sb_ad(ji,jj,jk) + ze3tb * zrhsad
731	sa_ad(ji,jj,jk) = sa_ad(ji,jj,jk) + p2dt(jk) * ze3tn * zrhsad
732	zrhsad = 0.0_wp
733	END DO
734	END DO
735	END DO
736	DO jj = 2, jpjm1
737	DO ji = fs_2, fs_jpim1
738	ze3tb = 1.0_wp ! before scale factor at T-point
739	ze3tn = 1.0_wp ! now scale factor at T-point
740	sb_ad(ji,jj,1) = sb_ad(ji,jj,1) + ze3tb * sa_ad(ji,jj,1)
741	sa_ad(ji,jj,1) = p2dt(1) * ze3tn * sa_ad(ji,jj,1)
742	END DO
743	END DO
744	END SUBROUTINE tra_zdf_imp_adj
745	SUBROUTINE tra_zdf_imp_adj_tst( kumadt )
746	!!-----------------------------------------------------------------------
747	!!
748	!! * ROUTINE tra_zdf_imp_adj_tst *
749	!!
750	!! ** Purpose : Test the adjoint routine.
751	!!
752	!! ** Method : Verify the scalar product
753	!!
754	!! ( L dx )^T W dy = dx^T L^T W dy
755	!!
756	!! where L = tangent routine
757	!! L^T = adjoint routine
758	!! W = diagonal matrix of scale factors
759	!! dx = input perturbation (random field)
760	!! dy = L dx
761	!!
762	!!
763	!! History :
764	!! ! 08-08 (A. Vidard)
765	!!-----------------------------------------------------------------------
766	!! * Modules used
767
768	!! * Arguments
769	INTEGER, INTENT(IN) :: &
770	& kumadt ! Output unit
771
772	!! * Local declarations
773	INTEGER :: &
774	& istp, &
775	& jstp, &
776	& ji, & ! dummy loop indices
777	& jj, &
778	& jk
779	INTEGER, DIMENSION(jpi,jpj) :: &
780	& iseed_2d ! 2D seed for the random number generator
781	REAL(KIND=wp) :: &
782	& zsp1, & ! scalar product involving the tangent routine
783	& zsp2 ! scalar product involving the adjoint routine
784	REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: &
785	& zta_tlin , & ! Tangent input
786	& ztb_tlin , & ! Tangent input
787	& zsa_tlin , & ! Tangent input
788	& zsb_tlin , & ! Tangent input
789	& zta_tlout, & ! Tangent output
790	& zsa_tlout, & ! Tangent output
791	& zta_adin , & ! Adjoint input
792	& zsa_adin , & ! Adjoint input
793	& zta_adout, & ! Adjoint output
794	& ztb_adout, & ! Adjoint output
795	& zsa_adout, & ! Adjoint output
796	& zsb_adout, & ! Adjoint output
797	& zr ! 3D random field
798	CHARACTER(LEN=14) :: cl_name
799	! Allocate memory
800
801	ALLOCATE( &
802	& zta_tlin( jpi,jpj,jpk), &
803	& zsa_tlin( jpi,jpj,jpk), &
804	& ztb_tlin( jpi,jpj,jpk), &
805	& zsb_tlin( jpi,jpj,jpk), &
806	& zta_tlout(jpi,jpj,jpk), &
807	& zsa_tlout(jpi,jpj,jpk), &
808	& zta_adin( jpi,jpj,jpk), &
809	& zsa_adin( jpi,jpj,jpk), &
810	& zta_adout(jpi,jpj,jpk), &
811	& zsa_adout(jpi,jpj,jpk), &
812	& ztb_adout(jpi,jpj,jpk), &
813	& zsb_adout(jpi,jpj,jpk), &
814	& zr( jpi,jpj,jpk) &
815	& )
816	!==================================================================
817	! 1) dx = ( un_tl, vn_tl, hdivn_tl ) and
818	! dy = ( hdivb_tl, hdivn_tl )
819	!==================================================================
820
821	! initialization (normally done in traldf)
822	l_traldf_rot = .TRUE.
823
824	! Test for time steps nit000 and nit000 + 1 (the matrix changes)
825
826	DO jstp = nit000, nit000 + 2
827	istp = jstp
828	IF ( jstp == nit000+2 ) istp = nitend
829
830	!--------------------------------------------------------------------
831	! Reset the tangent and adjoint variables
832	!--------------------------------------------------------------------
833	zta_tlin( :,:,:) = 0.0_wp
834	ztb_tlin( :,:,:) = 0.0_wp
835	zsa_tlin( :,:,:) = 0.0_wp
836	zsb_tlin( :,:,:) = 0.0_wp
837	zta_tlout(:,:,:) = 0.0_wp
838	zsa_tlout(:,:,:) = 0.0_wp
839	zta_adin( :,:,:) = 0.0_wp
840	zsa_adin( :,:,:) = 0.0_wp
841	zta_adout(:,:,:) = 0.0_wp
842	zsa_adout(:,:,:) = 0.0_wp
843	ztb_adout(:,:,:) = 0.0_wp
844	zsb_adout(:,:,:) = 0.0_wp
845	zr( :,:,:) = 0.0_wp
846	tb_ad(:,:,:) = 0.0_wp
847	sb_ad(:,:,:) = 0.0_wp
848	!--------------------------------------------------------------------
849	! Initialize the tangent input with random noise: dx
850	!--------------------------------------------------------------------
851
852	DO jj = 1, jpj
853	DO ji = 1, jpi
854	iseed_2d(ji,jj) = - ( 596035 + &
855	& mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
856	END DO
857	END DO
858	CALL grid_random( iseed_2d, zr, 'T', 0.0_wp, stdt )
859	DO jk = 1, jpk
860	DO jj = nldj, nlej
861	DO ji = nldi, nlei
862	zta_tlin(ji,jj,jk) = zr(ji,jj,jk)
863	END DO
864	END DO
865	END DO
866
867	DO jj = 1, jpj
868	DO ji = 1, jpi
869	iseed_2d(ji,jj) = - ( 352791 + &
870	& mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
871	END DO
872	END DO
873	CALL grid_random( iseed_2d, zr, 'T', 0.0_wp, stdt )
874	DO jk = 1, jpk
875	DO jj = nldj, nlej
876	DO ji = nldi, nlei
877	ztb_tlin(ji,jj,jk) = zr(ji,jj,jk)
878	END DO
879	END DO
880	END DO
881
882	DO jj = 1, jpj
883	DO ji = 1, jpi
884	iseed_2d(ji,jj) = - ( 142746 + &
885	& mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
886	END DO
887	END DO
888	CALL grid_random( iseed_2d, zr, 'T', 0.0_wp, stds )
889	DO jk = 1, jpk
890	DO jj = nldj, nlej
891	DO ji = nldi, nlei
892	zsa_tlin(ji,jj,jk) = zr(ji,jj,jk)
893	END DO
894	END DO
895	END DO
896
897	DO jj = 1, jpj
898	DO ji = 1, jpi
899	iseed_2d(ji,jj) = - ( 214934 + &
900	& mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
901	END DO
902	END DO
903	CALL grid_random( iseed_2d, zr, 'T', 0.0_wp, stds )
904	DO jk = 1, jpk
905	DO jj = nldj, nlej
906	DO ji = nldi, nlei
907	zsb_tlin(ji,jj,jk) = zr(ji,jj,jk)
908	END DO
909	END DO
910	END DO
911
912
913	ta_tl(:,:,:) = zta_tlin(:,:,:)
914	sa_tl(:,:,:) = zsa_tlin(:,:,:)
915	tb_tl(:,:,:) = ztb_tlin(:,:,:)
916	sb_tl(:,:,:) = zsb_tlin(:,:,:)
917	CALL tra_zdf_imp_tan ( istp, rdttra )
918	zta_tlout(:,:,:) = ta_tl(:,:,:)
919	zsa_tlout(:,:,:) = sa_tl(:,:,:)
920
921	!--------------------------------------------------------------------
922	! Initialize the adjoint variables: dy^* = W dy
923	!--------------------------------------------------------------------
924
925	DO jk = 1, jpk
926	DO jj = nldj, nlej
927	DO ji = nldi, nlei
928	zta_adin(ji,jj,jk) = zta_tlout(ji,jj,jk) &
929	& * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
930	& * tmask(ji,jj,jk) * wesp_t(jk)
931	zsa_adin(ji,jj,jk) = zsa_tlout(ji,jj,jk) &
932	& * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
933	& * tmask(ji,jj,jk) * wesp_s(jk)
934	END DO
935	END DO
936	END DO
937	!--------------------------------------------------------------------
938	! Compute the scalar product: ( L dx )^T W dy
939	!--------------------------------------------------------------------
940
941	zsp1 = DOT_PRODUCT( zta_tlout, zta_adin ) &
942	& + DOT_PRODUCT( zsa_tlout, zsa_adin )
943
944	!--------------------------------------------------------------------
945	! Call the adjoint routine: dx^* = L^T dy^*
946	!--------------------------------------------------------------------
947
948	ta_ad(:,:,:) = zta_adin(:,:,:)
949	sa_ad(:,:,:) = zsa_adin(:,:,:)
950
951	CALL tra_zdf_imp_adj ( istp, rdttra )
952
953	zta_adout(:,:,:) = ta_ad(:,:,:)
954	zsa_adout(:,:,:) = sa_ad(:,:,:)
955	ztb_adout(:,:,:) = tb_ad(:,:,:)
956	zsb_adout(:,:,:) = sb_ad(:,:,:)
957	zsp2 = DOT_PRODUCT( zta_tlin, zta_adout ) &
958	& + DOT_PRODUCT( zsa_tlin, zsa_adout ) &
959	& + DOT_PRODUCT( ztb_tlin, ztb_adout ) &
960	& + DOT_PRODUCT( zsb_tlin, zsb_adout )
961
962	! 14 char:'12345678901234'
963	IF ( istp == nit000 ) THEN
964	cl_name = 'trazdfimpadjT1'
965	ELSEIF ( istp == nit000 +1 ) THEN
966	cl_name = 'trazdfimpadjT2'
967	ELSEIF ( istp == nitend ) THEN
968	cl_name = 'trazdfimpadjT3'
969	END IF
970	CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
971
972	END DO
973
974	DEALLOCATE( &
975	& zta_tlin, &
976	& ztb_tlin, &
977	& zsa_tlin, &
978	& zsb_tlin, &
979	& zta_tlout, &
980	& zsa_tlout, &
981	& zta_adin, &
982	& zsa_adin, &
983	& zta_adout, &
984	& ztb_adout, &
985	& zsa_adout, &
986	& zsb_adout, &
987	& zr &
988	& )
989
990
991
992	END SUBROUTINE tra_zdf_imp_adj_tst
993	#if defined key_tst_tlm
994	SUBROUTINE tra_zdf_imp_tlm_tst( kumadt )
995	!!-----------------------------------------------------------------------
996	!!
997	!! * ROUTINE tra_adv_zdf_imp_tlm_tst *
998	!!
999	!! ** Purpose : Test the adjoint routine.
1000	!!
1001	!! ** Method : Verify the tangent with Taylor expansion
1002	!!
1003	!! M(x+hdx) = M(x) + L(hdx) + O(h^2)
1004	!!
1005	!! where L = tangent routine
1006	!! M = direct routine
1007	!! dx = input perturbation (random field)
1008	!! h = ration on perturbation
1009	!!
1010	!! In the tangent test we verify that:
1011	!! M(x+h*dx) - M(x)
1012	!! g(h) = ------------------ ---> 1 as h ---> 0
1013	!! L(h*dx)
1014	!! and
1015	!! g(h) - 1
1016	!! f(h) = ---------- ---> k (costant) as h ---> 0
1017	!! p
1018	!!
1019	!! History :
1020	!! ! 09-07 (A. Vigilant)
1021	!!-----------------------------------------------------------------------
1022	!! * Modules used
1023	USE oce , ONLY: & ! ocean dynamics and active tracers
1024	& tb, sb, ta, sa
1025	USE tamtrj ! writing out state trajectory
1026	USE par_tlm, ONLY: &
1027	& tlm_bch, &
1028	& cur_loop, &
1029	& h_ratio
1030	USE istate_mod
1031	USE wzvmod ! vertical velocity
1032	USE gridrandom, ONLY: &
1033	& grid_rd_sd
1034	USE trj_tam
1035	USE oce , ONLY: & ! ocean dynamics and tracers variables
1036	& tb, sb, ta, sa
1037	USE trazdf_imp ! vertical component of the tracer mixing trend
1038	USE opatam_tst_ini, ONLY: &
1039	& tlm_namrd
1040	USE tamctl, ONLY: & ! Control parameters
1041	& numtan, numtan_sc
1042	!! * Arguments
1043	INTEGER, INTENT(IN) :: &
1044	& kumadt ! Output unit
1045
1046	!! * Local declarations
1047	INTEGER :: &
1048	& istp, &
1049	& jstp, &
1050	& ji, & ! dummy loop indices
1051	& jj, &
1052	& jk
1053	REAL(KIND=wp) :: &
1054	& zsp1, & ! scalar product involving the tangent routine
1055	& zsp1_Ta, &
1056	& zsp1_Sa, &
1057	& zsp2, & ! scalar product involving the tangent routine
1058	& zsp2_Ta, &
1059	& zsp2_Sa, &
1060	& zsp3, & ! scalar product involving the tangent routine
1061	& zsp3_Ta, &
1062	& zsp3_Sa, &
1063	& zzsp, & ! scalar product involving the tangent routine
1064	& zzsp_Ta, &
1065	& zzsp_Sa, &
1066	& gamma, &
1067	& zgsp1, &
1068	& zgsp2, &
1069	& zgsp3, &
1070	& zgsp4, &
1071	& zgsp5, &
1072	& zgsp6, &
1073	& zgsp7
1074	REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: &
1075	& ztb_tlin, & ! Tangent input
1076	& zsb_tlin, & ! Tangent input
1077	& zta_tlin, & ! Tangent input
1078	& zsa_tlin, & ! Tangent input
1079	& zta_out , & ! Direct output
1080	& zsa_out , & ! Direct output
1081	& zta_wop , & ! Direct output w/o perturbation
1082	& zsa_wop , & ! Direct output w/o perturbation
1083	& zr ! 3D random field
1084	CHARACTER(LEN=14) ::&
1085	& cl_name
1086	CHARACTER (LEN=128) :: file_out, file_wop, file_xdx
1087	CHARACTER (LEN=90) :: &
1088	& FMT
1089	REAL(KIND=wp), DIMENSION(100):: &
1090	& zscta, zscsa, &
1091	& zscerrta, &
1092	& zscerrsa
1093	INTEGER, DIMENSION(100):: &
1094	& iipostb, iipossb, &
1095	& iiposta, iipossa, &
1096	& ijpostb, ijpossb, &
1097	& ijposta, ijpossa, &
1098	& ikpostb, ikpossb, &
1099	& ikposta, ikpossa
1100	INTEGER:: &
1101	& ii, &
1102	& isamp=40, &
1103	& jsamp=40, &
1104	& ksamp=10, &
1105	& numsctlm
1106	REAL(KIND=wp), DIMENSION(jpi,jpj,jpk) :: &
1107	& zerrta, zerrsa
1108	! Allocate memory
1109	ALLOCATE( &
1110	& ztb_tlin( jpi,jpj,jpk), &
1111	& zsb_tlin( jpi,jpj,jpk), &
1112	& zta_tlin( jpi,jpj,jpk), &
1113	& zsa_tlin( jpi,jpj,jpk), &
1114	& zta_out ( jpi,jpj,jpk), &
1115	& zsa_out ( jpi,jpj,jpk), &
1116	& zta_wop ( jpi,jpj,jpk), &
1117	& zsa_wop ( jpi,jpj,jpk), &
1118	& zr( jpi,jpj,jpk) &
1119	& )
1120	!--------------------------------------------------------------------
1121	! Reset variables
1122	!--------------------------------------------------------------------
1123	ztb_tlin( :,:,:) = 0.0_wp
1124	zsb_tlin( :,:,:) = 0.0_wp
1125	zta_tlin( :,:,:) = 0.0_wp
1126	zsa_tlin( :,:,:) = 0.0_wp
1127	zta_out ( :,:,:) = 0.0_wp
1128	zsa_out ( :,:,:) = 0.0_wp
1129	zta_wop ( :,:,:) = 0.0_wp
1130	zsa_wop ( :,:,:) = 0.0_wp
1131	zr( :,:,:) = 0.0_wp
1132	!--------------------------------------------------------------------
1133	! Output filename Xn=F(X0)
1134	!--------------------------------------------------------------------
1135	CALL tlm_namrd
1136	gamma = h_ratio
1137	file_wop='trj_wop_trazdf_imp'
1138	file_xdx='trj_xdx_trazdf_imp'
1139	!--------------------------------------------------------------------
1140	! Initialize the tangent input with random noise: dx
1141	!--------------------------------------------------------------------
1142	IF (( cur_loop .NE. 0) .OR. ( gamma .NE. 0.0_wp) )THEN
1143	CALL grid_rd_sd( 596035, zr, 'T', 0.0_wp, stdt)
1144	DO jk = 1, jpk
1145	DO jj = nldj, nlej
1146	DO ji = nldi, nlei
1147	ztb_tlin(ji,jj,jk) = zr(ji,jj,jk)
1148	END DO
1149	END DO
1150	END DO
1151	CALL grid_rd_sd( 352791, zr, 'T', 0.0_wp, stds)
1152	DO jk = 1, jpk
1153	DO jj = nldj, nlej
1154	DO ji = nldi, nlei
1155	zsb_tlin(ji,jj,jk) = zr(ji,jj,jk)
1156	END DO
1157	END DO
1158	END DO
1159	CALL grid_rd_sd( 142746, zr, 'T', 0.0_wp, stdt)
1160	DO jk = 1, jpk
1161	DO jj = nldj, nlej
1162	DO ji = nldi, nlei
1163	zta_tlin(ji,jj,jk) = zr(ji,jj,jk)
1164	END DO
1165	END DO
1166	END DO
1167	CALL grid_rd_sd( 214934, zr, 'T', 0.0_wp, stds)
1168	DO jk = 1, jpk
1169	DO jj = nldj, nlej
1170	DO ji = nldi, nlei
1171	zsa_tlin(ji,jj,jk) = zr(ji,jj,jk)
1172	END DO
1173	END DO
1174	END DO
1175	ENDIF
1176	!--------------------------------------------------------------------
1177	! Complete Init for Direct
1178	!-------------------------------------------------------------------
1179	IF ( tlm_bch /= 2 ) CALL istate_p
1180
1181	! *** initialize the reference trajectory
1182	! ------------
1183	CALL trj_rea( nit000-1, 1 )
1184	CALL trj_rea( nit000, 1 )
1185
1186	IF (( cur_loop .NE. 0) .OR. ( gamma .NE. 0.0_wp) )THEN
1187	ztb_tlin(:,:,:) = gamma * ztb_tlin(:,:,:)
1188	tb(:,:,:) = tb(:,:,:) + ztb_tlin(:,:,:)
1189
1190	zsb_tlin(:,:,:) = gamma * zsb_tlin(:,:,:)
1191	sb(:,:,:) = sb(:,:,:) + zsb_tlin(:,:,:)
1192
1193	zta_tlin(:,:,:) = gamma * zta_tlin(:,:,:)
1194	ta(:,:,:) = ta(:,:,:) + zta_tlin(:,:,:)
1195
1196	zsa_tlin(:,:,:) = gamma * zsa_tlin(:,:,:)
1197	sa(:,:,:) = sa(:,:,:) + zsa_tlin(:,:,:)
1198	ENDIF
1199	!--------------------------------------------------------------------
1200	! Compute the direct model F(X0,t=n) = Xn
1201	!--------------------------------------------------------------------
1202	IF ( tlm_bch /= 2 ) CALL tra_zdf_imp(nit000, rdttra)
1203	IF ( tlm_bch == 0 ) CALL trj_wri_spl(file_wop)
1204	IF ( tlm_bch == 1 ) CALL trj_wri_spl(file_xdx)
1205	!--------------------------------------------------------------------
1206	! Compute the Tangent
1207	!--------------------------------------------------------------------
1208	IF ( tlm_bch == 2 ) THEN
1209	!--------------------------------------------------------------------
1210	! Initialize the tangent variables: dy^* = W dy
1211	!--------------------------------------------------------------------
1212	CALL trj_rea( nit000-1, 1 )
1213	CALL trj_rea( nit000, 1 )
1214	tb_tl (:,:,:) = ztb_tlin (:,:,:)
1215	sb_tl (:,:,:) = zsb_tlin (:,:,:)
1216	ta_tl (:,:,:) = zta_tlin (:,:,:)
1217	sa_tl (:,:,:) = zsa_tlin (:,:,:)
1218
1219	!-----------------------------------------------------------------------
1220	! Initialization of the dynamics and tracer fields for the tangent
1221	!-----------------------------------------------------------------------
1222	CALL tra_zdf_imp_tan(nit000, rdttra)
1223
1224	!--------------------------------------------------------------------
1225	! Compute the scalar product: ( L(t0,tn) gamma dx0 ) )
1226	!--------------------------------------------------------------------
1227
1228	zsp2_Ta = DOT_PRODUCT( ta_tl, ta_tl )
1229	zsp2_Sa = DOT_PRODUCT( sa_tl, sa_tl )
1230
1231	zsp2 = zsp2_Ta + zsp2_Sa
1232
1233	!--------------------------------------------------------------------
1234	! Storing data
1235	!--------------------------------------------------------------------
1236	CALL trj_rd_spl(file_wop)
1237	zta_wop (:,:,:) = ta (:,:,:)
1238	zsa_wop (:,:,:) = sa (:,:,:)
1239	CALL trj_rd_spl(file_xdx)
1240	zta_out (:,:,:) = ta (:,:,:)
1241	zsa_out (:,:,:) = sa (:,:,:)
1242	!--------------------------------------------------------------------
1243	! Compute the Linearization Error
1244	! Nn = M( X0+gamma.dX0, t0,tn) - M(X0, t0,tn)
1245	! and
1246	! Compute the Linearization Error
1247	! En = Nn -TL(gamma.dX0, t0,tn)
1248	!--------------------------------------------------------------------
1249	! Warning: Here we re-use local variables z()_out and z()_wop
1250	ii=0
1251	DO jk = 1, jpk
1252	DO jj = 1, jpj
1253	DO ji = 1, jpi
1254	zta_out (ji,jj,jk) = zta_out (ji,jj,jk) - zta_wop (ji,jj,jk)
1255	zta_wop (ji,jj,jk) = zta_out (ji,jj,jk) - ta_tl (ji,jj,jk)
1256	IF ( ta_tl(ji,jj,jk) .NE. 0.0_wp ) &
1257	& zerrta(ji,jj,jk) = zta_out(ji,jj,jk)/ta_tl(ji,jj,jk)
1258	IF( (MOD(ji, isamp) .EQ. 0) .AND. &
1259	& (MOD(jj, jsamp) .EQ. 0) .AND. &
1260	& (MOD(jk, ksamp) .EQ. 0) ) THEN
1261	ii = ii+1
1262	iiposta(ii) = ji
1263	ijposta(ii) = jj
1264	ikposta(ii) = jk
1265	IF ( INT(tmask(ji,jj,jk)) .NE. 0) THEN
1266	zscta (ii) = zta_wop(ji,jj,jk)
1267	zscerrta (ii) = ( zerrta(ji,jj,jk) - 1.0_wp ) / gamma
1268	ENDIF
1269	ENDIF
1270	END DO
1271	END DO
1272	END DO
1273	ii=0
1274	DO jk = 1, jpk
1275	DO jj = 1, jpj
1276	DO ji = 1, jpi
1277	zsa_out (ji,jj,jk) = zsa_out (ji,jj,jk) - zsa_wop (ji,jj,jk)
1278	zsa_wop (ji,jj,jk) = zsa_out (ji,jj,jk) - sa_tl (ji,jj,jk)
1279	IF ( sa_tl(ji,jj,jk) .NE. 0.0_wp ) &
1280	& zerrsa(ji,jj,jk) = zsa_out(ji,jj,jk)/sa_tl(ji,jj,jk)
1281	IF( (MOD(ji, isamp) .EQ. 0) .AND. &
1282	& (MOD(jj, jsamp) .EQ. 0) .AND. &
1283	& (MOD(jk, ksamp) .EQ. 0) ) THEN
1284	ii = ii+1
1285	iipossa(ii) = ji
1286	ijpossa(ii) = jj
1287	ikpossa(ii) = jk
1288	IF ( INT(tmask(ji,jj,jk)) .NE. 0) THEN
1289	zscsa (ii) = zsa_wop(ji,jj,jk)
1290	zscerrsa (ii) = ( zerrsa(ji,jj,jk) - 1.0_wp ) / gamma
1291	ENDIF
1292	ENDIF
1293	END DO
1294	END DO
1295	END DO
1296
1297	zsp1_Ta = DOT_PRODUCT( zta_out, zta_out )
1298	zsp1_Sa = DOT_PRODUCT( zsa_out, zsa_out )
1299
1300	zsp1 = zsp1_Ta + zsp1_Sa
1301
1302	zsp3_Ta = DOT_PRODUCT( zta_wop, zta_wop )
1303	zsp3_Sa = DOT_PRODUCT( zsa_wop, zsa_wop )
1304
1305	zsp3 = zsp3_Ta + zsp3_Sa
1306
1307	!--------------------------------------------------------------------
1308	! Print the linearization error En - norme 2
1309	!--------------------------------------------------------------------
1310	! 14 char:'12345678901234'
1311	cl_name = 'trazdf_tam:En '
1312	zzsp = SQRT(zsp3)
1313	zzsp_Ta = SQRT(zsp3_Ta)
1314	zzsp_Sa = SQRT(zsp3_Sa)
1315	zgsp5 = zzsp
1316	CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio )
1317
1318	!--------------------------------------------------------------------
1319	! Compute TLM norm2
1320	!--------------------------------------------------------------------
1321	zzsp = SQRT(zsp2)
1322	zzsp_Ta = SQRT(zsp2_Ta)
1323	zzsp_Sa = SQRT(zsp2_Sa)
1324	zgsp4 = zzsp
1325	cl_name = 'trazdf_tam:Ln2'
1326	CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio )
1327
1328	!--------------------------------------------------------------------
1329	! Print the linearization error Nn - norme 2
1330	!--------------------------------------------------------------------
1331	zzsp = SQRT(zsp1)
1332	zzsp_Ta = SQRT(zsp1_Ta)
1333	zzsp_Sa = SQRT(zsp1_Sa)
1334
1335	cl_name = 'trazdf:Mhdx-Mx'
1336	CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio )
1337
1338	zgsp3 = SQRT( zsp3/zsp2 )
1339	zgsp7 = zgsp3/gamma
1340	zgsp1 = zzsp
1341	zgsp2 = zgsp1 / zgsp4
1342	zgsp6 = (zgsp2 - 1.0_wp)/gamma
1343
1344	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)"
1345	WRITE(numtan,FMT) 'tzdfimp ', cur_loop, h_ratio, zgsp1, zgsp2, zgsp3, zgsp4, zgsp5, zgsp6, zgsp7
1346
1347	!--------------------------------------------------------------------
1348	! Unitary calculus
1349	!--------------------------------------------------------------------
1350	FMT = "(A8,2X,A8,2X,I4.4,2X,E6.1,2X,I4.4,2X,I4.4,2X,I4.4,2X,E20.13,1X)"
1351	cl_name = 'tzdfimp '
1352	IF(lwp) THEN
1353	DO ii=1, 100, 1
1354	IF ( zscta(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscta ', &
1355	& cur_loop, h_ratio, ii, iiposta(ii), ijposta(ii), zscta(ii)
1356	ENDDO
1357	DO ii=1, 100, 1
1358	IF ( zscsa(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscsa ', &
1359	& cur_loop, h_ratio, ii, iipossa(ii), ijpossa(ii), zscsa(ii)
1360	ENDDO
1361	DO ii=1, 100, 1
1362	IF ( zscerrta(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscerrta ', &
1363	& cur_loop, h_ratio, ii, iiposta(ii), ijposta(ii), zscerrta(ii)
1364	ENDDO
1365	DO ii=1, 100, 1
1366	IF ( zscerrsa(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscerrsa ', &
1367	& cur_loop, h_ratio, ii, iipossa(ii), ijpossa(ii), zscerrsa(ii)
1368	ENDDO
1369	! write separator
1370	WRITE(numtan_sc, "(A4)") '===='
1371	ENDIF
1372
1373	ENDIF
1374
1375	DEALLOCATE( &
1376	& ztb_tlin, zsb_tlin, &
1377	& zta_tlin, zsa_tlin, &
1378	& zta_out, zsa_out, &
1379	& zta_wop, zsa_wop, &
1380	& zr &
1381	& )
1382	END SUBROUTINE tra_zdf_imp_tlm_tst
1383	#endif
1384	#endif
1385	END MODULE trazdf_imp_tam

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

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

Context Navigation

source: branches/TAM_V3_0/NEMOTAM/OPATAM_SRC/TRA/trazdf_imp_tam.F90 @ 2587

Download in other formats: