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.

trazdf_imp_tam.F90 in branches/TAM_V3_0/NEMOTAM/OPATAM_SRC/TRA – NEMO

Context Navigation

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

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

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

Download in other formats: