Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

eosbn2_tam.F90 @ 2578

Last change on this file since 2578 was 2578, checked in by rblod, 13 years ago
first import of NEMOTAM 3.2.2
File size: 195.3 KB

Line
1	MODULE eosbn2_tam
2	!!==============================================================================
3	!! * MODULE eosbn2_tam *
4	!! Ocean diagnostic variable : equation of state - in situ and potential density
5	!! - Brunt-Vaisala frequency
6	!! Tangent and Adjoint Module
7	!!===========================================================================
8	!! History of the direct Module :
9	!! OPA ! 1989-03 (O. Marti) Original code
10	!! 6.0 ! 1994-07 (G. Madec, M. Imbard) add bn2
11	!! 6.0 ! 1994-08 (G. Madec) Add Jackett & McDougall eos
12	!! 7.0 ! 1996-01 (G. Madec) statement function for e3
13	!! 8.1 ! 1997-07 (G. Madec) density instead of volumic mass
14	!! - ! 1999-02 (G. Madec, N. Grima) semi-implicit pressure gradient
15	!! 8.2 ! 2001-09 (M. Ben Jelloul) bugfix on linear eos
16	!! NEMO 1.0 ! 2002-10 (G. Madec) add eos_init
17	!! - ! 2002-11 (G. Madec, A. Bozec) partial step, eos_insitu_2d
18	!! - ! 2003-08 (G. Madec) F90, free form
19	!! 3.0 ! 2006-08 (G. Madec) add tfreez function
20	!! ! 2003-09 (M. Balmaseda) compute refrence rho prof
21	!! History of the TAM Module :
22	!! 8.2 ! 2005-03 (F. Van den Berghe, A. Weaver, N. Daget)
23	!! - eostan.F
24	!! 9.0 ! 2007-07 (K. Mogensen) Initial version based on eostan.F
25	!! ! 2008-07 (A. Vidard) bug fix in computation of prd_tl if neos=1
26	!! ! 2008-11 (A. Vidard) TAM of the 06-08 version
27	!! NEMO 3.2 ! 2010-04 (F. Vigilant) version 3.2
28	!!----------------------------------------------------------------------
29
30	!!----------------------------------------------------------------------
31	!! Direct subroutines
32	!! eos : generic interface of the equation of state
33	!! eos_insitu : Compute the in situ density
34	!! eos_insitu_pot : Compute the insitu and surface referenced potential
35	!! volumic mass
36	!! eos_insitu_2d : Compute the in situ density for 2d fields
37	!! eos_insitu_pot_1pt : Compute the in situ density for a single point
38	!! eos_bn2 : Compute the Brunt-Vaisala frequency
39	!! tfreez : Compute the surface freezing temperature (NOT IN TAM)
40	!! eos_init : set eos parameters (namelist)
41	!! eos_rprof : Compute the in situ density of a reference profile
42	!!----------------------------------------------------------------------
43	!! * Modules used
44	#if defined key_zdfddm
45	USE oce_tam , ONLY: &
46	& rrau_tl, &
47	& rrau_ad
48	#endif
49	USE dom_oce , ONLY: & ! ocean space and time domain
50	& tmask, &
51	& e1t, &
52	& e2t, &
53	#if defined key_zco
54	& e3t_0, &
55	& e3w_0, &
56	& gdept_0, &
57	& gdepw_0, &
58	#else
59	& e3t, &
60	& e3w, &
61	& gdept, &
62	& gdepw, &
63	#endif
64	& mig, &
65	& mjg, &
66	& nperio, &
67	& nldi, &
68	& nldj, &
69	& nlei, &
70	& nlej
71	USE par_kind , ONLY: &
72	& wp
73	USE par_oce , ONLY: &
74	& jpi, &
75	& jpj, &
76	& jpk, &
77	& jpim1, &
78	& jpjm1, &
79	& jpkm1, &
80	& jpiglo
81	USE oce , ONLY: &
82	& tn, &
83	& sn
84	USE phycst , ONLY: & ! physical constants
85	& rau0, &
86	& grav
87	USE in_out_manager, ONLY: & ! I/O manager
88	& lwp, &
89	& ctmp1, &
90	& numout, &
91	& ctl_stop
92	#if defined key_zdfddm
93	USE zdfddm ! vertical physics: double diffusion
94	#endif
95	USE eosbn2 , ONLY: &
96	& eos_init, &
97	& neos_init, &
98	& nn_eos, &
99	& rn_alpha, &
100	& rn_beta, &
101	& ralpbet
102	USE gridrandom , ONLY: & ! Random Gaussian noise on grids
103	& grid_random
104	USE dotprodfld , ONLY: & ! Computes dot product for 3D and 2D fields
105	& dot_product
106	USE tstool_tam , ONLY: &
107	& prntst_adj, &
108	& prntst_tlm, & !
109	& stds, &
110	& stdt
111	IMPLICIT NONE
112	PRIVATE
113
114	!! * Interface
115	INTERFACE eos_tan
116	MODULE PROCEDURE eos_insitu_tan, eos_insitu_pot_tan, eos_pot_1pt_tan, eos_insitu_2d_tan
117	END INTERFACE
118	INTERFACE eos_adj
119	MODULE PROCEDURE eos_insitu_adj, eos_insitu_pot_adj, eos_pot_1pt_adj, eos_insitu_2d_adj
120	END INTERFACE
121	INTERFACE bn2_tan
122	MODULE PROCEDURE eos_bn2_tan
123	END INTERFACE
124	INTERFACE bn2_adj
125	MODULE PROCEDURE eos_bn2_adj
126	END INTERFACE
127
128	!! * Routine accessibility
129	PUBLIC eos_tan ! called by step.F90, inidtr.F90, tranpc.F90 and intgrd.F90
130	PUBLIC bn2_tan ! called by step.F90
131	PUBLIC eos_adj ! called by step.F90, inidtr.F90, tranpc.F90 and intgrd.F90
132	PUBLIC bn2_adj ! called by step.F90
133	#if defined key_tam
134	PUBLIC eos_adj_tst
135	PUBLIC bn2_adj_tst
136	#if defined key_tst_tlm
137	PUBLIC eos_tlm_tst
138	PUBLIC bn2_tlm_tst
139	#endif
140	#endif
141
142
143	!! * Substitutions
144	# include "domzgr_substitute.h90"
145	# include "vectopt_loop_substitute.h90"
146
147	CONTAINS
148
149	SUBROUTINE eos_insitu_tan( ptem, psal, ptem_tl, psal_tl, prd_tl )
150	!!-----------------------------------------------------------------------
151	!!
152	!! * ROUTINE eos_insitu_tan : TL OF ROUTINE eos_insitu *
153	!!
154	!! ** Purpose of direct routine : Compute the in situ density
155	!! (ratio rho/rau0) from potential temperature and salinity
156	!! using an equation of state defined through the namelist
157	!! parameter nn_eos.
158	!!
159	!! ** Method of direct routine : 3 cases:
160	!! nn_eos = 0 : Jackett and McDougall (1994) equation of state.
161	!! the in situ density is computed directly as a function of
162	!! potential temperature relative to the surface (the opa t
163	!! variable), salt and pressure (assuming no pressure variation
164	!! along geopotential surfaces, i.e. the pressure p in decibars
165	!! is approximated by the depth in meters.
166	!! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0
167	!! with pressure p decibars
168	!! potential temperature t deg celsius
169	!! salinity s psu
170	!! reference volumic mass rau0 kg/m**3
171	!! in situ volumic mass rho kg/m**3
172	!! in situ density anomalie prd no units
173	!! Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar,
174	!! t = 40 deg celcius, s=40 psu
175	!! nn_eos = 1 : linear equation of state function of temperature only
176	!! prd(t) = 0.0285 - rn_alpha * t
177	!! nn_eos = 2 : linear equation of state function of temperature and
178	!! salinity
179	!! prd(t,s) = rn_beta * s - rn_alpha * tn - 1.
180	!! Note that no boundary condition problem occurs in this routine
181	!! as (ptem,psal) are defined over the whole domain.
182	!!
183	!! ** Comments on Adjoint Routine :
184	!! Care has been taken to avoid division by zero when computing
185	!! the inverse of the square root of salinity at masked salinity
186	!! points.
187	!!
188	!! * Arguments
189	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: &
190	& ptem, & ! potential temperature
191	& psal, & ! salinity
192	& ptem_tl, & ! TL of potential temperature
193	& psal_tl ! TL of salinity
194	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out ) :: &
195	& prd_tl ! TL of potential density (surface referenced)
196	!! * Local declarations
197	INTEGER :: ji, jj, jk ! dummy loop indices
198	REAL(wp) :: & ! temporary scalars
199	zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw, &
200	zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, &
201	zttl, zstl, zhtl, zsrtl, zr1tl, zr2tl, zr3tl, &
202	zr4tl, zrhoptl, zetl, zbwtl, &
203	zbtl, zdtl, zctl, zawtl, zatl, zb1tl, za1tl, &
204	zkwtl, zk0tl, zpes, zrdc1, zrdc2, zeps, &
205	zmask, zrau0r
206	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zws
207	!!----------------------------------------------------------------------
208
209	SELECT CASE ( nn_eos )
210
211	CASE ( 0 ) !== Jackett and McDougall (1994) formulation ==!
212	zrau0r = 1.e0 / rau0
213	#ifdef key_sp
214	zeps = 1.e-7
215	#else
216	zeps = 1.e-14
217	#endif
218
219	!CDIR NOVERRCHK
220	zws(:,:,:) = SQRT( ABS( psal(:,:,:) ) )
221	!
222	DO jk = 1, jpkm1
223	DO jj = 1, jpj
224	DO ji = 1, jpi
225	zt = ptem(ji,jj,jk)
226	zs = psal(ji,jj,jk)
227	zh = fsdept(ji,jj,jk) ! depth
228	zsr= zws(ji,jj,jk) ! square root salinity
229	! compute volumic mass pure water at atm pressure
230	zr1= ( ( ( ( 6.536332e-9zt-1.120083e-6 )zt+1.001685e-4)*zt &
231	-9.095290e-3 )zt+6.793952e-2 )zt+999.842594
232	! seawater volumic mass atm pressure
233	zr2= ( ( ( 5.3875e-9zt-8.2467e-7 ) zt+7.6438e-5 ) *zt &
234	-4.0899e-3 ) *zt+0.824493
235	zr3= ( -1.6546e-6zt+1.0227e-4 ) zt-5.72466e-3
236	zr4= 4.8314e-4
237
238	! potential volumic mass (reference to the surface)
239	zrhop= ( zr4zs + zr3zsr + zr2 ) *zs + zr1
240
241	! add the compression terms
242	ze = ( -3.508914e-8zt-1.248266e-8 ) zt-2.595994e-6
243	zbw= ( 1.296821e-6zt-5.782165e-9 ) zt+1.045941e-4
244	zb = zbw + ze * zs
245
246	zd = -2.042967e-2
247	zc = (-7.267926e-5zt+2.598241e-3 ) zt+0.1571896
248	zaw= ( ( 5.939910e-6zt+2.512549e-3 ) zt-0.1028859 ) *zt - 4.721788
249	za = ( zdzsr + zc ) zs + zaw
250
251	zb1= (-0.1909078zt+7.390729 ) zt-55.87545
252	za1= ( ( 2.326469e-3zt+1.553190)zt-65.00517 ) *zt+1044.077
253	zkw= ( ( (-1.361629e-4zt-1.852732e-2 ) zt-30.41638 ) zt + 2098.925 ) zt+190925.6
254	zk0= ( zb1zsr + za1 )zs + zkw
255
256	! Tangent linear part
257
258	zttl = ptem_tl(ji,jj,jk)
259	zstl = psal_tl(ji,jj,jk)
260
261	zsrtl= ( 1.0 / MAX( 2.*zsr, zeps ) ) &
262	& * tmask(ji,jj,jk) * zstl
263
264	zr1tl= ( ( ( ( 5.6.536332e-9 zt &
265	& -4.1.120083e-6 ) zt &
266	& +3.1.001685e-4 ) zt &
267	& -2.9.095290e-3 ) zt &
268	& + 6.793952e-2 ) * zttl
269
270	zr2tl= ( ( ( 4.5.3875e-9 zt &
271	& -3.8.2467e-7 ) zt &
272	& +2.7.6438e-5 ) zt &
273	& - 4.0899e-3 ) * zttl
274
275	zr3tl= ( -2.1.6546e-6 zt &
276	& + 1.0227e-4 ) * zttl
277
278	zrhoptl= zr1tl &
279	& + zs * zr2tl &
280	& + zsr * zs * zr3tl &
281	& + zr3 * zs * zsrtl &
282	& + ( 2. * zr4 * zs + zr2 &
283	& + zr3 * zsr ) * zstl
284
285	zetl = ( -2.3.508914e-8 zt &
286	& - 1.248266e-8 ) * zttl
287
288	zbwtl= ( 2.1.296821e-6 zt &
289	& - 5.782165e-9 ) * zttl
290
291	zbtl = zbwtl &
292	& + zs * zetl &
293	& + ze * zstl
294
295	zctl = ( -2.7.267926e-5 zt &
296	& + 2.598241e-3 ) * zttl
297
298	zawtl= ( ( 3.5.939910e-6 zt &
299	& +2.2.512549e-3 ) zt &
300	& - 0.1028859 ) * zttl
301
302	zatl = zawtl &
303	& + zd * zs * zsrtl &
304	& + zs * zctl &
305	& + ( zd * zsr + zc ) * zstl
306
307	zb1tl= ( -2.0.1909078 zt &
308	& + 7.390729 ) * zttl
309
310	za1tl= ( ( 3.2.326469e-3 zt &
311	& +2.1.553190 ) zt &
312	& - 65.00517 ) * zttl
313
314	zkwtl= ( ( ( -4.1.361629e-4 zt &
315	& -3.1.852732e-2 ) zt &
316	& -2.30.41638 ) zt &
317	& + 2098.925 ) * zttl
318
319	zk0tl= zkwtl &
320	& + zb1 * zs * zsrtl &
321	& + zs * zsr * zb1tl &
322	& + zs * za1tl &
323	& + ( zb1 * zsr + za1 ) * zstl
324
325	! Masked in situ density anomaly
326
327	zrdc1 = 1.0 / ( zk0 - zh * ( za - zh * zb ) )
328	zrdc2 = 1.0 / ( 1.0 - zh * zrdc1 )
329
330	prd_tl(ji,jj,jk) = tmask(ji,jj,jk) * zrdc2 * &
331	& ( zrhoptl &
332	& - zrdc2 * zh * zrdc1*2 zrhop &
333	& * ( zk0tl &
334	& - zh * ( zatl &
335	& - zh * zbtl ) ) )&
336	& * zrau0r
337	END DO
338	END DO
339	END DO
340	!
341	CASE ( 1 ) !== Linear formulation function of temperature only ==!
342	DO jk = 1, jpkm1
343	prd_tl(:,:,jk) = ( - rn_alpha * ptem_tl(:,:,jk) ) * tmask(:,:,jk)
344	END DO
345	!
346	CASE ( 2 ) !== Linear formulation function of temperature and salinity ==!
347
348	! ! ===============
349	DO jk = 1, jpkm1
350	prd_tl(:,:,jk) = ( rn_beta * psal_tl(:,:,jk) - rn_alpha * ptem_tl(:,:,jk ) ) * tmask(:,:,jk)
351	END DO
352	!
353	END SELECT
354
355	END SUBROUTINE eos_insitu_tan
356
357	SUBROUTINE eos_insitu_pot_tan( ptem, psal, ptem_tl, psal_tl, prd_tl, prhop_tl)
358	!!----------------------------------------------------------------------
359	!! * ROUTINE eos_insitu_pot_tan *
360	!!
361	!! ** Purpose or the direct routine:
362	!! Compute the in situ density (ratio rho/rau0) and the
363	!! potential volumic mass (Kg/m3) from potential temperature and
364	!! salinity fields using an equation of state defined through the
365	!! namelist parameter nn_eos.
366	!!
367	!! ** Method :
368	!! nn_eos = 0 : Jackett and McDougall (1994) equation of state.
369	!! the in situ density is computed directly as a function of
370	!! potential temperature relative to the surface (the opa t
371	!! variable), salt and pressure (assuming no pressure variation
372	!! along geopotential surfaces, i.e. the pressure p in decibars
373	!! is approximated by the depth in meters.
374	!! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0
375	!! rhop(t,s) = rho(t,s,0)
376	!! with pressure p decibars
377	!! potential temperature t deg celsius
378	!! salinity s psu
379	!! reference volumic mass rau0 kg/m**3
380	!! in situ volumic mass rho kg/m**3
381	!! in situ density anomalie prd no units
382	!!
383	!! Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar,
384	!! t = 40 deg celcius, s=40 psu
385	!!
386	!! nn_eos = 1 : linear equation of state function of temperature only
387	!! prd(t) = ( rho(t) - rau0 ) / rau0 = 0.028 - rn_alpha * t
388	!! rhop(t,s) = rho(t,s)
389	!!
390	!! nn_eos = 2 : linear equation of state function of temperature and
391	!! salinity
392	!! prd(t,s) = ( rho(t,s) - rau0 ) / rau0
393	!! = rn_beta * s - rn_alpha * tn - 1.
394	!! rhop(t,s) = rho(t,s)
395	!! Note that no boundary condition problem occurs in this routine
396	!! as (tn,sn) or (ta,sa) are defined over the whole domain.
397	!!
398	!! ** Action : - prd , the in situ density (no units)
399	!! - prhop, the potential volumic mass (Kg/m3)
400	!!
401	!! References :
402	!! Jackett, D.R., and T.J. McDougall. J. Atmos. Ocean. Tech., 1994
403	!! Brown, J. A. and K. A. Campana. Mon. Weather Rev., 1978
404	!!
405	!!----------------------------------------------------------------------
406	!! * Arguments
407	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: &
408	ptem, & ! potential temperature
409	psal ! salinity
410	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: &
411	ptem_tl,& ! potential temperature
412	psal_tl ! salinity
413	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out ) :: &
414	prd_tl, & ! potential density (surface referenced)
415	prhop_tl ! potential density (surface referenced)
416	!! * Local declarations
417	INTEGER :: ji, jj, jk ! dummy loop indices
418	REAL(wp) :: & ! temporary scalars
419	zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw, &
420	zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, &
421	zttl, zstl, zhtl, zsrtl, zr1tl, zr2tl, zr3tl, &
422	zr4tl, zrhoptl, zetl, zbwtl, &
423	zbtl, zdtl, zctl, zawtl, zatl, zb1tl, za1tl, &
424	zkwtl, zk0tl, zpes, zrdc1, zrdc2, zeps, &
425	zmask, zrau0r
426	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zws
427	!!----------------------------------------------------------------------
428
429	SELECT CASE ( nn_eos )
430
431	CASE ( 0 ) !== Jackett and McDougall (1994) formulation ==!
432	zrau0r = 1.e0 / rau0
433	#ifdef key_sp
434	zeps = 1.e-7
435	#else
436	zeps = 1.e-14
437	#endif
438
439	!CDIR NOVERRCHK
440	zws(:,:,:) = SQRT( ABS( psal(:,:,:) ) )
441	!
442	DO jk = 1, jpkm1
443	DO jj = 1, jpj
444	DO ji = 1, jpi
445	zt = ptem(ji,jj,jk)
446	zs = psal(ji,jj,jk)
447	zh = fsdept(ji,jj,jk) ! depth
448	zsr = zws(ji,jj,jk) ! square root salinity
449	! compute volumic mass pure water at atm pressure
450	zr1 = ( ( ( ( 6.536332e-9 * zt - 1.120083e-6 ) * zt &
451	& + 1.001685e-4 ) * zt - 9.095290e-3 ) * zt &
452	& + 6.793952e-2 ) * zt + 999.842594
453	! seawater volumic mass atm pressure
454	zr2 = ( ( ( 5.3875e-9 * zt - 8.2467e-7 ) * zt &
455	& + 7.6438e-5 ) * zt - 4.0899e-3 ) * zt &
456	& + 0.824493
457	zr3 = ( -1.6546e-6 * zt + 1.0227e-4 ) * zt - 5.72466e-3
458	zr4 = 4.8314e-4
459
460	! potential volumic mass (reference to the surface)
461	zrhop= ( zr4 * zs + zr3 * zsr + zr2 ) * zs + zr1
462
463	! add the compression terms
464	ze = ( -3.508914e-8 * zt - 1.248266e-8 ) * zt - 2.595994e-6
465	zbw = ( 1.296821e-6 * zt - 5.782165e-9 ) * zt + 1.045941e-4
466	zb = zbw + ze * zs
467
468	zd = -2.042967e-2
469	zc = (-7.267926e-5 * zt + 2.598241e-3 ) * zt + 0.1571896
470	zaw= ( ( 5.939910e-6 * zt + 2.512549e-3 ) * zt - 0.1028859 ) * zt - 4.721788
471	za = ( zd * zsr + zc ) * zs + zaw
472
473	zb1 = (-0.1909078 * zt + 7.390729 ) * zt - 55.87545
474	za1 = ( ( 2.326469e-3 * zt + 1.553190 ) * zt - 65.00517 &
475	& ) * zt + 1044.077
476	zkw = ( ( (-1.361629e-4 * zt - 1.852732e-2 ) * zt - 30.41638 &
477	& ) * zt + 2098.925 ) * zt + 190925.6
478	zk0 = ( zb1 * zsr + za1 ) * zs + zkw
479
480
481	zrdc1 = 1.0 / ( zk0 - zh * ( za - zh * zb ) )
482	zrdc2 = 1.0 / ( 1.0 - zh * zrdc1 )
483
484	! Tangent linear part
485
486	zttl = ptem_tl(ji,jj,jk)
487	zstl = psal_tl(ji,jj,jk)
488
489	zsrtl= ( 1.0 / MAX( 2.*zsr, zeps ) ) &
490	& * tmask(ji,jj,jk) * zstl
491
492	zr1tl= ( ( ( ( 5.6.536332e-9 zt &
493	& -4.1.120083e-6 ) zt &
494	& +3.1.001685e-4 ) zt &
495	& -2.9.095290e-3 ) zt &
496	& + 6.793952e-2 ) * zttl
497
498	zr2tl= ( ( ( 4.5.3875e-9 zt &
499	& -3.8.2467e-7 ) zt &
500	& +2.7.6438e-5 ) zt &
501	& - 4.0899e-3 ) * zttl
502
503	zr3tl= ( -2.1.6546e-6 zt &
504	& + 1.0227e-4 ) * zttl
505
506	zrhoptl= zr1tl &
507	& + zs * zr2tl &
508	& + zsr * zs * zr3tl &
509	& + zr3 * zs * zsrtl &
510	& + ( 2. * zr4 * zs + zr2 &
511	& + zr3 * zsr ) * zstl
512
513	prhop_tl(ji,jj,jk) = zrhoptl * tmask(ji,jj,jk)
514
515	zetl = ( -2.3.508914e-8 zt &
516	& - 1.248266e-8 ) * zttl
517
518	zbwtl= ( 2.1.296821e-6 zt &
519	& - 5.782165e-9 ) * zttl
520
521	zbtl = zbwtl &
522	& + zs * zetl &
523	& + ze * zstl
524
525	zctl = ( -2.7.267926e-5 zt &
526	& + 2.598241e-3 ) * zttl
527
528	zawtl= ( ( 3.5.939910e-6 zt &
529	& +2.2.512549e-3 ) zt &
530	& - 0.1028859 ) * zttl
531
532	zatl = zawtl &
533	& + zd * zs * zsrtl &
534	& + zs * zctl &
535	& + ( zd * zsr + zc ) * zstl
536
537	zb1tl= ( -2.0.1909078 zt &
538	& + 7.390729 ) * zttl
539
540	za1tl= ( ( 3.2.326469e-3 zt &
541	& +2.1.553190 ) zt &
542	& - 65.00517 ) * zttl
543
544	zkwtl= ( ( ( -4.1.361629e-4 zt &
545	& -3.1.852732e-2 ) zt &
546	& -2.30.41638 ) zt &
547	& + 2098.925 ) * zttl
548
549	zk0tl= zkwtl &
550	& + zb1 * zs * zsrtl &
551	& + zs * zsr * zb1tl &
552	& + zs * za1tl &
553	& + ( zb1 * zsr + za1 ) * zstl
554
555	! Masked in situ density anomaly
556
557	prd_tl(ji,jj,jk) = tmask(ji,jj,jk) * zrdc2 * &
558	& ( zrhoptl &
559	& - zrdc2 * zh * zrdc1*2 zrhop &
560	& * ( zk0tl &
561	& - zh * ( zatl &
562	& - zh * zbtl ) ) )&
563	& * zrau0r
564	END DO
565	END DO
566	END DO
567	!
568	CASE ( 1 ) !== Linear formulation = F( temperature ) ==!
569	DO jk = 1, jpkm1
570	prd_tl (:,:,jk) = ( - rn_alpha * ptem_tl(:,:,jk) ) * tmask(:,:,jk)
571	prhop_tl(:,:,jk) = ( rau0 * prd_tl(:,:,jk) ) * tmask(:,:,jk)
572	END DO
573	!
574	CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==!
575	DO jk = 1, jpkm1
576	prd_tl(:,:,jk) = ( rn_beta * psal_tl(:,:,jk) - rn_alpha * ptem_tl(:,:,jk) ) * tmask(:,:,jk)
577	prhop_tl(:,:,jk) = ( rau0 * prd_tl(:,:,jk) ) * tmask(:,:,jk)
578	END DO
579	!
580	END SELECT
581
582	END SUBROUTINE eos_insitu_pot_tan
583	SUBROUTINE eos_pot_1pt_tan( ptem, psal, ptem_tl, psal_tl, prhop_tl)
584	!!----------------------------------------------------------------------
585	!! * ROUTINE eos_pot_1pt_tan *
586	!!
587	!! ** Purpose of the direct routine:
588	!! Compute the
589	!! potential volumic mass (Kg/m3) from potential temperature and
590	!! salinity fields using an equation of state defined through the
591	!! namelist parameter neos.
592	!!
593	!! ** Method :
594	!! nn_eos = 0 : Jackett and McDougall (1994) equation of state.
595	!! the in situ density is computed directly as a function of
596	!! potential temperature relative to the surface (the opa t
597	!! variable), salt and pressure (assuming no pressure variation
598	!! along geopotential surfaces, i.e. the pressure p in decibars
599	!! is approximated by the depth in meters.
600	!! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0
601	!! rhop(t,s) = rho(t,s,0)
602	!! with pressure p decibars
603	!! potential temperature t deg celsius
604	!! salinity s psu
605	!! reference volumic mass rau0 kg/m**3
606	!! in situ volumic mass rho kg/m**3
607	!! in situ density anomalie prd no units
608	!!
609	!! Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar,
610	!! t = 40 deg celcius, s=40 psu
611	!!
612	!! nn_eos = 1 : linear equation of state function of temperature only
613	!! prd(t) = ( rho(t) - rau0 ) / rau0 = 0.028 - ralpha * t
614	!! rhop(t,s) = rho(t,s)
615	!!
616	!! nn_eos = 2 : linear equation of state function of temperature and
617	!! salinity
618	!! prd(t,s) = ( rho(t,s) - rau0 ) / rau0
619	!! = rn_beta * s - rn_alpha * tn - 1.
620	!! rhop(t,s) = rho(t,s)
621	!! Note that no boundary condition problem occurs in this routine
622	!! as (tn,sn) or (ta,sa) are defined over the whole domain.
623	!!
624	!! ** Action : - prd , the in situ density (no units)
625	!! - prhop, the potential volumic mass (Kg/m3)
626	!!
627	!! References :
628	!! Jackett, D.R., and T.J. McDougall. J. Atmos. Ocean. Tech., 1994
629	!! Brown, J. A. and K. A. Campana. Mon. Weather Rev., 1978
630	!!
631	!! History of the direct routine:
632	!! 4.0 ! 89-03 (O. Marti)
633	!! ! 94-08 (G. Madec)
634	!! ! 96-01 (G. Madec) statement function for e3
635	!! ! 97-07 (G. Madec) introduction of neos, OPA8.1
636	!! ! 97-07 (G. Madec) density instead of volumic mass
637	!! ! 99-02 (G. Madec, N. Grima) semi-implicit pressure gradient
638	!! ! 01-09 (M. Ben Jelloul) bugfix
639	!! 9.0 ! 03-08 (G. Madec) F90, free form
640	!! ! 06-11 (G. Smith) single point case
641	!! History of the tangent routine:
642	!! 9.0 ! 08-07 (A. Vidard) Initial version
643	!!----------------------------------------------------------------------
644	!! * Arguments
645	REAL(wp), INTENT( in ) :: &
646	ptem, & ! potential temperature
647	psal ! salinity
648	REAL(wp), INTENT( in ) :: &
649	ptem_tl, & ! potential temperature
650	psal_tl ! salinity
651	REAL(wp), INTENT( out ) :: &
652	prhop_tl ! potential density (surface referenced)
653
654	!! * Local declarations
655	INTEGER :: ji, jj, jk ! dummy loop indices
656	REAL(wp) :: & ! temporary scalars
657	& zt, zs, zsr, zr1, zr2, zr3, zr4, zrhop, &
658	& zttl, zstl, zsrtl, zr1tl, zr2tl, zr3tl, zr4tl, zrhoptl
659	REAL(wp) :: zws, zwstl, zeps
660	!!----------------------------------------------------------------------
661
662	#ifdef key_sp
663	zeps = 1.e-7
664	#else
665	zeps = 1.e-14
666	#endif
667
668	SELECT CASE ( nn_eos )
669
670	CASE ( 0 ) ! Jackett and McDougall (1994) formulation
671	zws = SQRT( ABS( psal ) )
672
673	zt = ptem
674	zs = psal
675	! square root salinity
676	zsr= zws
677	! compute volumic mass pure water at atm pressure
678	zr1= ( ( ( ( 6.536332e-9zt-1.120083e-6 )zt+1.001685e-4)*zt &
679	-9.095290e-3 )zt+6.793952e-2 )zt+999.842594
680	! seawater volumic mass atm pressure
681	zr2= ( ( ( 5.3875e-9zt-8.2467e-7 ) zt+7.6438e-5 ) *zt &
682	-4.0899e-3 ) *zt+0.824493
683	zr3= ( -1.6546e-6zt+1.0227e-4 ) zt-5.72466e-3
684	zr4= 4.8314e-4
685
686	! ==================
687	! tangent linear part
688	! ==================
689
690	zwstl = 1/MAX( 2zws , zeps ) psal_tl
691	zttl = ptem_tl
692	zstl = psal_tl
693	! square root salinity
694	zsrtl= zwstl
695	! compute volumic mass pure water at atm pressure
696	zr1tl= ( ( ( ( 5.6.536332e-9 zt &
697	& -4.1.120083e-6 ) zt &
698	& +3.1.001685e-4 ) zt &
699	& -2.9.095290e-3 ) zt &
700	& + 6.793952e-2 ) * zttl
701	! seawater volumic mass atm pressure
702	zr2tl= ( ( ( 4.5.3875e-9 zt &
703	& -3.8.2467e-7 ) zt &
704	& +2.7.6438e-5 ) zt &
705	& - 4.0899e-3 ) * zttl
706
707	zr3tl= ( -2.1.6546e-6 zt &
708	& + 1.0227e-4 ) * zttl
709
710	! potential volumic mass (reference to the surface)
711	zrhoptl= zr1tl &
712	& + zs * zr2tl &
713	& + zsr * zs * zr3tl &
714	& + zr3 * zs * zsrtl &
715	& + ( 2. * zr4 * zs + zr2 &
716	& + zr3 * zsr ) * zstl
717
718	! save potential volumic mass
719	prhop_tl = zrhoptl
720
721	CASE ( 1 ) ! Linear formulation function of temperature only
722
723	zttl = ptem_tl
724	! ... potential volumic mass
725	prhop_tl = ( rau0 * ( - rn_alpha * zttl ) )
726
727	CASE ( 2 ) ! Linear formulation function of temperature and salinity
728
729	zttl = ptem_tl
730	zstl = psal_tl
731	! ... potential volumic mass
732	prhop_tl = rau0 * ( rn_beta * zstl - rn_alpha * zttl )
733
734	CASE DEFAULT
735
736	WRITE(ctmp1,*) ' bad flag value for nn_eos = ', nn_eos
737	CALL ctl_stop( ctmp1 )
738
739	END SELECT
740
741	END SUBROUTINE eos_pot_1pt_tan
742	SUBROUTINE eos_insitu_2d_tan( ptem, psal, pdep, ptem_tl, psal_tl, prd_tl )
743	!!-----------------------------------------------------------------------
744	!!
745	!! * ROUTINE eos_insitu_2d_tan : TL OF ROUTINE eos_insitu_2d *
746	!!
747	!! ** Purpose of direct routine : Compute the in situ density
748	!! (ratio rho/rau0) from potential temperature and salinity
749	!! using an equation of state defined through the namelist
750	!! parameter nn_eos. * 2D field case
751	!!
752	!! ** Method of direct routine : 3 cases:
753	!! nn_eos = 0 : Jackett and McDougall (1994) equation of state.
754	!! the in situ density is computed directly as a function of
755	!! potential temperature relative to the surface (the opa t
756	!! variable), salt and pressure (assuming no pressure variation
757	!! along geopotential surfaces, i.e. the pressure p in decibars
758	!! is approximated by the depth in meters.
759	!! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0
760	!! with pressure p decibars
761	!! potential temperature t deg celsius
762	!! salinity s psu
763	!! reference volumic mass rau0 kg/m**3
764	!! in situ volumic mass rho kg/m**3
765	!! in situ density anomalie prd no units
766	!! Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar,
767	!! t = 40 deg celcius, s=40 psu
768	!! nn_eos = 1 : linear equation of state function of temperature only
769	!! prd(t) = 0.0285 - ralpha * t
770	!! nn_eos = 2 : linear equation of state function of temperature and
771	!! salinity
772	!! prd(t,s) = rn_beta * s - rn_alpha * tn - 1.
773	!! Note that no boundary condition problem occurs in this routine
774	!! as (ptem,psal) are defined over the whole domain.
775	!!
776	!! ** Comments on Adjoint Routine :
777	!! Care has been taken to avoid division by zero when computing
778	!! the inverse of the square root of salinity at masked salinity
779	!! points.
780	!!
781	!! ** Action :
782	!!
783	!! References :
784	!!
785	!! History :
786	!! 8.2 ! 05-03 ((F. Van den Berghe, A. Weaver, N. Daget) - eostan.F
787	!! 9.0 ! 07-07 (K. Mogensen) Initial version based on eostan.F
788	!! ! 08-07 (A. Vidard) bug fix in computation of prd_tl if neos=1
789	!!-----------------------------------------------------------------------
790	!! * Modules used
791	!! * Arguments
792	REAL(wp), DIMENSION(jpi,jpj), INTENT( in ) :: &
793	& ptem, & ! potential temperature
794	& psal, & ! salinity
795	& pdep, & ! depth
796	& ptem_tl, & ! TL of potential temperature
797	& psal_tl ! TL of salinity
798	REAL(wp), DIMENSION(jpi,jpj), INTENT( out ) :: &
799	& prd_tl ! TL of potential density (surface referenced)
800
801	!! * Local declarations
802	INTEGER :: ji, jj, jk ! dummy loop indices
803	REAL(wp) :: & ! temporary scalars
804	zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw, &
805	zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, &
806	zttl, zstl, zhtl, zsrtl, zr1tl, zr2tl, zr3tl, &
807	zr4tl, zrhoptl, zetl, zbwtl, &
808	zbtl, zdtl, zctl, zawtl, zatl, zb1tl, za1tl, &
809	zkwtl, zk0tl, zpes, zrdc1, zrdc2, zeps, &
810	zmask
811	REAL(wp), DIMENSION(jpi,jpj) :: zws
812	!!----------------------------------------------------------------------
813
814	SELECT CASE ( nn_eos )
815
816	CASE ( 0 ) !== Jackett and McDougall (1994) formulation ==!
817
818	#ifdef key_sp
819	zeps = 1.e-7
820	#else
821	zeps = 1.e-14
822	#endif
823
824	!CDIR NOVERRCHK
825	DO jj = 1, jpjm1
826	!CDIR NOVERRCHK
827	DO ji = 1, fs_jpim1 ! vector opt.
828	zws(ji,jj) = SQRT( ABS( psal(ji,jj) ) )
829	END DO
830	END DO
831	DO jj = 1, jpjm1
832	DO ji = 1, fs_jpim1 ! vector opt.
833
834	zmask = tmask(ji,jj,1) ! land/sea bottom mask = surf. mask
835
836	zt = ptem (ji,jj) ! interpolated T
837	zs = psal (ji,jj) ! interpolated S
838	zsr= zws(ji,jj) ! square root of interpolated S
839	zh = pdep(ji,jj) ! depth at the partial step level
840	! compute volumic mass pure water at atm pressure
841	zr1= ( ( ( ( 6.536332e-9zt-1.120083e-6 )zt+1.001685e-4)*zt &
842	& -9.095290e-3 )zt+6.793952e-2 )zt+999.842594
843	! seawater volumic mass atm pressure
844	zr2= ( ( ( 5.3875e-9zt-8.2467e-7 ) zt+7.6438e-5 ) *zt &
845	& -4.0899e-3 ) *zt+0.824493
846	zr3= ( -1.6546e-6zt+1.0227e-4 ) zt-5.72466e-3
847	zr4= 4.8314e-4
848
849	! potential volumic mass (reference to the surface)
850	zrhop= ( zr4zs + zr3zsr + zr2 ) *zs + zr1
851
852	! add the compression terms
853	ze = ( -3.508914e-8zt-1.248266e-8 ) zt-2.595994e-6
854	zbw= ( 1.296821e-6zt-5.782165e-9 ) zt+1.045941e-4
855	zb = zbw + ze * zs
856
857	zd = -2.042967e-2
858	zc = (-7.267926e-5zt+2.598241e-3 ) zt+0.1571896
859	zaw= ( ( 5.939910e-6zt+2.512549e-3 ) zt-0.1028859 ) *zt - 4.721788
860	za = ( zdzsr + zc ) zs + zaw
861
862	zb1= (-0.1909078zt+7.390729 ) zt-55.87545
863	za1= ( ( 2.326469e-3zt+1.553190)zt-65.00517 ) *zt+1044.077
864	zkw= ( ( (-1.361629e-4zt-1.852732e-2 ) zt-30.41638 ) zt + 2098.925 ) zt+190925.6
865	zk0= ( zb1zsr + za1 )zs + zkw
866
867	! Tangent linear part
868
869	zttl = ptem_tl(ji,jj)
870	zstl = psal_tl(ji,jj)
871
872	zsrtl= ( 1.0 / MAX( 2.*zsr, zeps ) ) &
873	& * tmask(ji,jj,1) * zstl
874
875	zr1tl= ( ( ( ( 5.6.536332e-9 zt &
876	& -4.1.120083e-6 ) zt &
877	& +3.1.001685e-4 ) zt &
878	& -2.9.095290e-3 ) zt &
879	& + 6.793952e-2 ) * zttl
880
881	zr2tl= ( ( ( 4.5.3875e-9 zt &
882	& -3.8.2467e-7 ) zt &
883	& +2.7.6438e-5 ) zt &
884	& - 4.0899e-3 ) * zttl
885
886	zr3tl= ( -2.1.6546e-6 zt &
887	& + 1.0227e-4 ) * zttl
888
889	zrhoptl= zr1tl &
890	& + zs * zr2tl &
891	& + zsr * zs * zr3tl &
892	& + zr3 * zs * zsrtl &
893	& + ( 2. * zr4 * zs + zr2 &
894	& + zr3 * zsr ) * zstl
895
896	zetl = ( -2.3.508914e-8 zt &
897	& - 1.248266e-8 ) * zttl
898
899	zbwtl= ( 2.1.296821e-6 zt &
900	& - 5.782165e-9 ) * zttl
901
902	zbtl = zbwtl &
903	& + zs * zetl &
904	& + ze * zstl
905
906	zctl = ( -2.7.267926e-5 zt &
907	& + 2.598241e-3 ) * zttl
908
909	zawtl= ( ( 3.5.939910e-6 zt &
910	& +2.2.512549e-3 ) zt &
911	& - 0.1028859 ) * zttl
912
913	zatl = zawtl &
914	& + zd * zs * zsrtl &
915	& + zs * zctl &
916	& + ( zd * zsr + zc ) * zstl
917
918	zb1tl= ( -2.0.1909078 zt &
919	& + 7.390729 ) * zttl
920
921	za1tl= ( ( 3.2.326469e-3 zt &
922	& +2.1.553190 ) zt &
923	& - 65.00517 ) * zttl
924
925	zkwtl= ( ( ( -4.1.361629e-4 zt &
926	& -3.1.852732e-2 ) zt &
927	& -2.30.41638 ) zt &
928	& + 2098.925 ) * zttl
929
930	zk0tl= zkwtl &
931	& + zb1 * zs * zsrtl &
932	& + zs * zsr * zb1tl &
933	& + zs * za1tl &
934	& + ( zb1 * zsr + za1 ) * zstl
935
936	! Masked in situ density anomaly
937
938	zrdc1 = 1.0 / ( zk0 - zh * ( za - zh * zb ) )
939	zrdc2 = 1.0 / ( 1.0 - zh * zrdc1 )
940
941	prd_tl(ji,jj) = tmask(ji,jj,1) * zrdc2 * &
942	& ( zrhoptl &
943	& - zrdc2 * zh * zrdc1*2 zrhop &
944	& * ( zk0tl &
945	& - zh * ( zatl &
946	& - zh * zbtl ) ) )&
947	& / rau0
948
949
950	END DO
951	END DO
952	!
953	CASE ( 1 ) !== Linear formulation = F( temperature ) ==!
954	DO jj = 1, jpjm1
955	DO ji = 1, fs_jpim1 ! vector opt.
956	prd_tl(ji,jj) = ( - rn_alpha * ptem_tl(ji,jj) ) * tmask(ji,jj,1)
957	END DO
958	END DO
959	!
960	CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==!
961	DO jj = 1, jpjm1
962	DO ji = 1, fs_jpim1 ! vector opt.
963	prd_tl (ji,jj) = ( rn_beta * psal_tl(ji,jj) - rn_alpha * ptem_tl(ji,jj) ) * tmask(ji,jj,1)
964	END DO
965	END DO
966	!
967	END SELECT
968
969
970	END SUBROUTINE eos_insitu_2d_tan
971
972	SUBROUTINE eos_insitu_adj(ptem, psal, ptem_ad, psal_ad, prd_ad)
973	!!-----------------------------------------------------------------------
974	!!
975	!! * ROUTINE eos_insitu_tan : Adjoint OF ROUTINE eos_insitu *
976	!!
977	!! ** Purpose of direct routine : Compute the in situ density
978	!! (ratio rho/rau0) from potential temperature and salinity
979	!! using an equation of state defined through the namelist
980	!! parameter nneos.
981	!!
982	!! ** Method of direct routine : 3 cases:
983	!! nn_eos = 0 : Jackett and McDougall (1994) equation of state.
984	!! the in situ density is computed directly as a function of
985	!! potential temperature relative to the surface (the opa t
986	!! variable), salt and pressure (assuming no pressure variation
987	!! along geopotential surfaces, i.e. the pressure p in decibars
988	!! is approximated by the depth in meters.
989	!! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0
990	!! with pressure p decibars
991	!! potential temperature t deg celsius
992	!! salinity s psu
993	!! reference volumic mass rau0 kg/m**3
994	!! in situ volumic mass rho kg/m**3
995	!! in situ density anomalie prd no units
996	!! Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar,
997	!! t = 40 deg celcius, s=40 psu
998	!! nn_eos = 1 : linear equation of state function of temperature only
999	!! prd(t) = 0.0285 - rn_alpha * t
1000	!! nn_eos = 2 : linear equation of state function of temperature and
1001	!! salinity
1002	!! prd(t,s) = rn_beta * s - rn_alpha * tn - 1.
1003	!! Note that no boundary condition problem occurs in this routine
1004	!! as (ptem,psal) are defined over the whole domain.
1005	!!
1006	!! ** Comments on Adjoint Routine :
1007	!! Care has been taken to avoid division by zero when computing
1008	!! the inverse of the square root of salinity at masked salinity
1009	!! points.
1010	!!
1011	!! ** Action :
1012	!!
1013	!! References :
1014	!!
1015	!! History :
1016	!! 8.2 ! 05-03 ((F. Van den Berghe, A. Weaver, N. Daget) - eostan.F
1017	!! 9.0 ! 08-08 (A. Vidard) 9.0 version
1018	!!-----------------------------------------------------------------------
1019	!! * Modules used
1020	!! * Arguments
1021	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: &
1022	ptem, & ! potential temperature
1023	psal ! salinity
1024	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: &
1025	ptem_ad, & ! potential temperature
1026	psal_ad ! salinity
1027	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: &
1028	prd_ad ! potential density (surface referenced)
1029	!! * Local declarations
1030	INTEGER :: ji, jj, jk ! dummy loop indices
1031	REAL(wp) :: & ! temporary scalars
1032	zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw, &
1033	zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, &
1034	ztad, zsad, zhad, zsrad, zr1ad, zr2ad, zr3ad, &
1035	zr4ad, zrhopad, zead, zbwad, &
1036	zbad, zdad, zcad, zawad, zaad, zb1ad, za1ad, &
1037	zkwad, zk0ad, zpes, zrdc1, zrdc2, zeps, &
1038	zmask, zrau0r
1039	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zws
1040	!!----------------------------------------------------------------------
1041
1042	! initialization of adjoint variables
1043	ztad = 0.0_wp
1044	zsad = 0.0_wp
1045	zhad = 0.0_wp
1046	zsrad = 0.0_wp
1047	zr1ad = 0.0_wp
1048	zr2ad = 0.0_wp
1049	zr3ad = 0.0_wp
1050	zr4ad = 0.0_wp
1051	zrhopad = 0.0_wp
1052	zead = 0.0_wp
1053	zbwad = 0.0_wp
1054	zbad = 0.0_wp
1055	zdad = 0.0_wp
1056	zcad = 0.0_wp
1057	zawad = 0.0_wp
1058	zaad = 0.0_wp
1059	zb1ad = 0.0_wp
1060	za1ad = 0.0_wp
1061	zkwad = 0.0_wp
1062	zk0ad = 0.0_wp
1063
1064	SELECT CASE ( nn_eos )
1065
1066	CASE ( 0 ) !== Jackett and McDougall (1994) formulation ==!
1067	zrau0r = 1.e0 / rau0
1068	#ifdef key_sp
1069	zeps = 1.e-7
1070	#else
1071	zeps = 1.e-14
1072	#endif
1073
1074	!CDIR NOVERRCHK
1075	zws(:,:,:) = SQRT( ABS( psal(:,:,:) ) )
1076	DO jk = 1, jpkm1
1077	DO jj = 1, jpj
1078	DO ji = 1, jpi
1079	zt = ptem(ji,jj,jk)
1080	zs = psal(ji,jj,jk)
1081	zh = fsdept(ji,jj,jk) ! depth
1082	zsr= zws(ji,jj,jk) ! square root salinity
1083	! compute volumic mass pure water at atm pressure
1084	zr1= ( ( ( ( 6.536332e-9zt-1.120083e-6 )zt+1.001685e-4)*zt &
1085	-9.095290e-3 )zt+6.793952e-2 )zt+999.842594
1086	! seawater volumic mass atm pressure
1087	zr2= ( ( ( 5.3875e-9zt-8.2467e-7 ) zt+7.6438e-5 ) *zt &
1088	-4.0899e-3 ) *zt+0.824493
1089	zr3= ( -1.6546e-6zt+1.0227e-4 ) zt-5.72466e-3
1090	zr4= 4.8314e-4
1091
1092	! potential volumic mass (reference to the surface)
1093	zrhop= ( zr4zs + zr3zsr + zr2 ) *zs + zr1
1094
1095	! add the compression terms
1096	ze = ( -3.508914e-8zt-1.248266e-8 ) zt-2.595994e-6
1097	zbw= ( 1.296821e-6zt-5.782165e-9 ) zt+1.045941e-4
1098	zb = zbw + ze * zs
1099
1100	zd = -2.042967e-2
1101	zc = (-7.267926e-5zt+2.598241e-3 ) zt+0.1571896
1102	zaw= ( ( 5.939910e-6zt+2.512549e-3 ) zt-0.1028859 ) *zt - 4.721788
1103	za = ( zdzsr + zc ) zs + zaw
1104
1105	zb1= (-0.1909078zt+7.390729 ) zt-55.87545
1106	za1= ( ( 2.326469e-3zt+1.553190)zt-65.00517 ) *zt+1044.077
1107	zkw= ( ( (-1.361629e-4zt-1.852732e-2 ) zt-30.41638 ) zt + 2098.925 ) zt+190925.6
1108	zk0= ( zb1zsr + za1 )zs + zkw
1109
1110	zrdc1 = 1.0 / ( zk0 - zh * ( za - zh * zb ) )
1111	zrdc2 = 1.0 / ( 1.0 - zh * zrdc1 )
1112	! ============
1113	! Adjoint part
1114	! ============
1115
1116	! Masked in situ density anomaly
1117
1118	zrhopad = zrhopad + prd_ad(ji,jj,jk) * tmask(ji,jj,jk) &
1119	& * zrdc2 * zrau0r
1120	zk0ad = zk0ad - prd_ad(ji,jj,jk) * tmask(ji,jj,jk) &
1121	& * zrdc2 * zrdc2 * zh &
1122	& * zrdc1*2 zrhop &
1123	& * zrau0r
1124	zaad = zaad + prd_ad(ji,jj,jk) * tmask(ji,jj,jk) &
1125	& * zrdc2 * zrdc2 * zh &
1126	& * zrdc1*2 zrhop &
1127	& * zh * zrau0r
1128	zbad = zbad - prd_ad(ji,jj,jk) * tmask(ji,jj,jk) &
1129	& * zrdc2 * zrdc2 * zh &
1130	& * zrdc1*2 zrhop &
1131	& * zh * zh * zrau0r
1132	prd_ad(ji,jj,jk) = 0.0_wp
1133
1134	zkwad = zkwad + zk0ad
1135	zsrad = zsrad + zk0ad * zb1 * zs
1136	zb1ad = zb1ad + zk0ad * zs * zsr
1137	za1ad = za1ad + zk0ad * zs
1138	zsad = zsad + zk0ad * ( zb1 * zsr + za1 )
1139	zk0ad = 0.0_wp
1140
1141	ztad = ztad + zkwad * ( ( (-4.1.361629e-4 zt &
1142	& -3.1.852732e-2 ) zt &
1143	& -2.30.41638 ) zt &
1144	& + 2098.925 )
1145	zkwad = 0.0_wp
1146
1147	ztad = ztad + za1ad * ( ( 3.2.326469e-3 zt &
1148	& +2.1.553190 ) zt &
1149	& - 65.00517 )
1150	za1ad = 0.0_wp
1151
1152	ztad = ztad + zb1ad * (-2.0.1909078 zt &
1153	& + 7.390729 )
1154	zb1ad = 0.0_wp
1155
1156	zawad = zawad + zaad
1157	zsrad = zsrad + zaad * zd * zs
1158	zcad = zcad + zaad * zs
1159	zsad = zsad + zaad * ( zd * zsr + zc )
1160	zaad = 0.0_wp
1161
1162	ztad = ztad + zawad * ( ( 3.5.939910e-6 zt &
1163	& +2.2.512549e-3 ) zt &
1164	& - 0.1028859 )
1165	zawad = 0.0_wp
1166
1167	ztad = ztad + zcad * (-2.7.267926e-5 zt &
1168	& + 2.598241e-3 )
1169	zcad = 0.0_wp
1170
1171	zbwad = zbwad + zbad
1172	zead = zead + zbad * zs
1173	zsad = zsad + zbad * ze
1174	zbad = 0.0_wp
1175
1176	ztad = ztad + zbwad * ( 2.1.296821e-6 zt &
1177	& - 5.782165e-9 )
1178	zbwad = 0.0_wp
1179
1180	ztad = ztad + zead * (-2.3.508914e-8 zt &
1181	& - 1.248266e-8 )
1182	zead = 0.0_wp
1183
1184	zr1ad = zr1ad + zrhopad
1185	zr2ad = zr2ad + zrhopad * zs
1186	zr3ad = zr3ad + zrhopad * zsr * zs
1187	zsrad = zsrad + zrhopad * zr3 * zs
1188	zsad = zsad + zrhopad * ( 2. * zr4 * zs + zr2 &
1189	& + zr3 * zsr )
1190	zrhopad = 0.0_wp
1191
1192	ztad = ztad + zr3ad * (-2.1.6546e-6 zt &
1193	& + 1.0227e-4 )
1194	zr3ad = 0.0_wp
1195
1196	ztad = ztad + zr2ad * ( ( ( 4.5.3875e-9 zt &
1197	& -3.8.2467e-7 ) zt &
1198	& +2.7.6438e-5 ) zt &
1199	& - 4.0899e-3 )
1200	zr2ad = 0.0_wp
1201
1202	ztad = ztad + zr1ad * ( ( ( ( 5.6.536332e-9 zt &
1203	& -4.1.120083e-6 ) zt &
1204	& +3.1.001685e-4 ) zt &
1205	& -2.9.095290e-3 ) zt &
1206	& + 6.793952e-2 )
1207	zr1ad = 0.0_wp
1208
1209	zsad = zsad + zsrad * ( 1.0 / MAX( 2.*zsr, zeps ) ) &
1210	& * tmask(ji,jj,jk)
1211	zsrad = 0.0_wp
1212
1213	psal_ad(ji,jj,jk) = psal_ad(ji,jj,jk) + zsad
1214	ptem_ad(ji,jj,jk) = ptem_ad(ji,jj,jk) + ztad
1215	ztad = 0.0_wp
1216	zsad = 0.0_wp
1217	END DO
1218	END DO
1219	END DO
1220	!
1221	CASE ( 1 ) !== Linear formulation function of temperature only ==!
1222	DO jk = 1, jpkm1
1223	ptem_ad(:,:,jk) = ptem_ad(:,:,jk) - rn_alpha * prd_ad(:,:,jk) * tmask(:,:,jk)
1224	prd_ad(:,:,jk) = 0.0_wp
1225	END DO
1226	!
1227	CASE ( 2 ) !== Linear formulation function of temperature and salinity ==!
1228	DO jk = 1, jpkm1
1229	ptem_ad(:,:,jk) = ptem_ad(:,:,jk) - rn_alpha * prd_ad(:,:,jk) * tmask(:,:,jk)
1230	psal_ad(:,:,jk) = psal_ad(:,:,jk) + rn_beta * prd_ad( :,:,jk) * tmask(:,:,jk)
1231	prd_ad( :,:,jk) = 0.0_wp
1232	END DO
1233	!
1234	END SELECT
1235	END SUBROUTINE eos_insitu_adj
1236
1237	SUBROUTINE eos_insitu_pot_adj ( ptem, psal, ptem_ad, psal_ad, prd_ad, prhop_ad )
1238	!!----------------------------------------------------------------------
1239	!! * ROUTINE eos_insitu_pot_adj *
1240	!!
1241	!! ** Purpose or the direct routine:
1242	!! Compute the in situ density (ratio rho/rau0) and the
1243	!! potential volumic mass (Kg/m3) from potential temperature and
1244	!! salinity fields using an equation of state defined through the
1245	!! namelist parameter nn_eos.
1246	!!
1247	!! ** Method :
1248	!! nn_eos = 0 : Jackett and McDougall (1994) equation of state.
1249	!! the in situ density is computed directly as a function of
1250	!! potential temperature relative to the surface (the opa t
1251	!! variable), salt and pressure (assuming no pressure variation
1252	!! along geopotential surfaces, i.e. the pressure p in decibars
1253	!! is approximated by the depth in meters.
1254	!! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0
1255	!! rhop(t,s) = rho(t,s,0)
1256	!! with pressure p decibars
1257	!! potential temperature t deg celsius
1258	!! salinity s psu
1259	!! reference volumic mass rau0 kg/m**3
1260	!! in situ volumic mass rho kg/m**3
1261	!! in situ density anomalie prd no units
1262	!!
1263	!! Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar,
1264	!! t = 40 deg celcius, s=40 psu
1265	!!
1266	!! neos = 1 : linear equation of state function of temperature only
1267	!! prd(t) = ( rho(t) - rau0 ) / rau0 = 0.028 - ralpha * t
1268	!! rhop(t,s) = rho(t,s)
1269	!!
1270	!! nn_eos = 2 : linear equation of state function of temperature and
1271	!! salinity
1272	!! prd(t,s) = ( rho(t,s) - rau0 ) / rau0
1273	!! = rn_beta * s - rn_alpha * tn - 1.
1274	!! rhop(t,s) = rho(t,s)
1275	!! Note that no boundary condition problem occurs in this routine
1276	!! as (tn,sn) or (ta,sa) are defined over the whole domain.
1277	!!
1278	!! ** Action : - prd , the in situ density (no units)
1279	!! - prhop, the potential volumic mass (Kg/m3)
1280	!!
1281	!! References :
1282	!! Jackett, D.R., and T.J. McDougall. J. Atmos. Ocean. Tech., 1994
1283	!! Brown, J. A. and K. A. Campana. Mon. Weather Rev., 1978
1284	!!
1285	!! History of the adjoint routine:
1286	!! 9.0 ! 08-06 (A. Vidard) Initial version
1287	!!----------------------------------------------------------------------
1288	!! * Arguments
1289
1290	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: &
1291	ptem, & ! potential temperature
1292	psal ! salinity
1293	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: &
1294	ptem_ad, & ! potential temperature
1295	psal_ad ! salinity
1296	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: &
1297	prd_ad, & ! potential density (surface referenced)
1298	prhop_ad
1299	!! * Local declarations
1300	INTEGER :: ji, jj, jk ! dummy loop indices
1301	REAL(wp) :: & ! temporary scalars
1302	zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw, &
1303	zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, &
1304	ztad, zsad, zhad, zsrad, zr1ad, zr2ad, zr3ad, &
1305	zr4ad, zrhopad, zead, zbwad, &
1306	zbad, zdad, zcad, zawad, zaad, zb1ad, za1ad, &
1307	zkwad, zk0ad, zpes, zrdc1, zrdc2, zeps, &
1308	zmask, zrau0r
1309	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zws
1310	!!----------------------------------------------------------------------
1311
1312	! initialization of adjoint variables
1313	ztad = 0.0_wp
1314	zsad = 0.0_wp
1315	zhad = 0.0_wp
1316	zsrad = 0.0_wp
1317	zr1ad = 0.0_wp
1318	zr2ad = 0.0_wp
1319	zr3ad = 0.0_wp
1320	zr4ad = 0.0_wp
1321	zrhopad = 0.0_wp
1322	zead = 0.0_wp
1323	zbwad = 0.0_wp
1324	zbad = 0.0_wp
1325	zdad = 0.0_wp
1326	zcad = 0.0_wp
1327	zawad = 0.0_wp
1328	zaad = 0.0_wp
1329	zb1ad = 0.0_wp
1330	za1ad = 0.0_wp
1331	zkwad = 0.0_wp
1332	zk0ad = 0.0_wp
1333
1334	SELECT CASE ( nn_eos )
1335
1336	CASE ( 0 ) !== Jackett and McDougall (1994) formulation ==!
1337	zrau0r = 1.e0 / rau0
1338	#ifdef key_sp
1339	zeps = 1.e-7
1340	#else
1341	zeps = 1.e-14
1342	#endif
1343
1344	!CDIR NOVERRCHK
1345	zws(:,:,:) = SQRT( ABS( psal(:,:,:) ) )
1346	!
1347	DO jk = jpkm1, 1, -1
1348	DO jj = jpj, 1, -1
1349	DO ji = jpi, 1, -1
1350	! direct recomputing
1351	zt = ptem(ji,jj,jk)
1352	zs = psal(ji,jj,jk)
1353	zh = fsdept(ji,jj,jk) ! depth
1354	zsr = zws(ji,jj,jk) ! square root salinity
1355	! compute volumic mass pure water at atm pressure
1356	zr1 = ( ( ( ( 6.536332e-9 * zt - 1.120083e-6 ) * zt &
1357	& + 1.001685e-4 ) * zt - 9.095290e-3 ) * zt &
1358	& + 6.793952e-2 ) * zt + 999.842594
1359	! seawater volumic mass atm pressure
1360	zr2 = ( ( ( 5.3875e-9 * zt - 8.2467e-7 ) * zt &
1361	& + 7.6438e-5 ) * zt - 4.0899e-3 ) * zt + 0.824493
1362	zr3 = ( -1.6546e-6 * zt + 1.0227e-4 ) * zt - 5.72466e-3
1363	zr4 = 4.8314e-4
1364	! potential volumic mass (reference to the surface)
1365	zrhop = ( zr4 * zs + zr3zsr + zr2 ) zs + zr1
1366	! add the compression terms
1367	ze = ( -3.508914e-8 * zt - 1.248266e-8 ) * zt - 2.595994e-6
1368	zbw = ( 1.296821e-6 * zt - 5.782165e-9 ) * zt + 1.045941e-4
1369	zb = zbw + ze * zs
1370
1371	zd = -2.042967e-2
1372	zc = (-7.267926e-5 * zt + 2.598241e-3 ) * zt + 0.1571896
1373	zaw= ( ( 5.939910e-6 * zt + 2.512549e-3 ) * zt - 0.1028859 &
1374	& ) * zt - 4.721788
1375	za = ( zd * zsr + zc ) * zs + zaw
1376
1377	zb1= (-0.1909078 * zt + 7.390729 ) * zt - 55.87545
1378	za1= ( ( 2.326469e-3 * zt + 1.553190 ) * zt - 65.00517 &
1379	& ) * zt + 1044.077
1380	zkw= ( ( (-1.361629e-4 * zt - 1.852732e-2 ) * zt - 30.41638 &
1381	& ) * zt + 2098.925 ) * zt + 190925.6
1382	zk0= ( zb1 * zsr + za1 ) * zs + zkw
1383
1384
1385	zrdc1 = 1.0 / ( zk0 - zh * ( za - zh * zb ) )
1386	zrdc2 = 1.0 / ( 1.0 - zh * zrdc1 )
1387
1388	! ============
1389	! Adjoint part
1390	! ============
1391
1392	! Masked in situ density anomaly
1393
1394	zrhopad = zrhopad + prd_ad(ji,jj,jk) * tmask(ji,jj,jk) &
1395	& * zrdc2 * zrau0r
1396	zk0ad = zk0ad - prd_ad(ji,jj,jk) * tmask(ji,jj,jk) &
1397	& * zrdc2 * zrdc2 * zh &
1398	& * zrdc1*2 zrhop &
1399	& * zrau0r
1400	zaad = zaad + prd_ad(ji,jj,jk) * tmask(ji,jj,jk) &
1401	& * zrdc2 * zrdc2 * zh &
1402	& * zrdc1*2 zrhop &
1403	& * zh * zrau0r
1404	zbad = zbad - prd_ad(ji,jj,jk) * tmask(ji,jj,jk) &
1405	& * zrdc2 * zrdc2 * zh &
1406	& * zrdc1*2 zrhop &
1407	& * zh * zh * zrau0r
1408	prd_ad(ji,jj,jk) = 0.0_wp
1409
1410	zkwad = zkwad + zk0ad
1411	zsrad = zsrad + zk0ad * zb1 * zs
1412	zb1ad = zb1ad + zk0ad * zs * zsr
1413	za1ad = za1ad + zk0ad * zs
1414	zsad = zsad + zk0ad * ( zb1 * zsr + za1 )
1415	zk0ad = 0.0_wp
1416
1417	ztad = ztad + zkwad * ( ( (-4.1.361629e-4 zt &
1418	& -3.1.852732e-2 ) zt &
1419	& -2.30.41638 ) zt &
1420	& + 2098.925 )
1421	zkwad = 0.0_wp
1422
1423	ztad = ztad + za1ad * ( ( 3.2.326469e-3 zt &
1424	& +2.1.553190 ) zt &
1425	& - 65.00517 )
1426	za1ad = 0.0_wp
1427
1428	ztad = ztad + zb1ad * (-2.0.1909078 zt &
1429	& + 7.390729 )
1430	zb1ad = 0.0_wp
1431
1432	zawad = zawad + zaad
1433	zsrad = zsrad + zaad * zd * zs
1434	zcad = zcad + zaad * zs
1435	zsad = zsad + zaad * ( zd * zsr + zc )
1436	zaad = 0.0_wp
1437
1438	ztad = ztad + zawad * ( ( 3.5.939910e-6 zt &
1439	& +2.2.512549e-3 ) zt &
1440	& - 0.1028859 )
1441	zawad = 0.0_wp
1442
1443	ztad = ztad + zcad * (-2.7.267926e-5 zt &
1444	& + 2.598241e-3 )
1445	zcad = 0.0_wp
1446
1447
1448	zsad = zsad + zbad * ze
1449	zead = zead + zbad * zs
1450	zbwad = zbwad + zbad
1451	zbad = 0.0_wp
1452
1453	ztad = ztad + zbwad * ( 2.1.296821e-6 zt &
1454	& - 5.782165e-9 )
1455	zbwad = 0.0_wp
1456
1457	ztad = ztad + zead * (-2.3.508914e-8 zt &
1458	& - 1.248266e-8 )
1459	zead = 0.0_wp
1460
1461	zrhopad = zrhopad + prhop_ad(ji,jj,jk) * tmask(ji,jj,jk)
1462	prhop_ad(ji,jj,jk) = 0.0_wp
1463
1464	zr1ad = zr1ad + zrhopad
1465	zr2ad = zr2ad + zrhopad * zs
1466	zr3ad = zr3ad + zrhopad * zsr * zs
1467	zsrad = zsrad + zrhopad * zr3 * zs
1468	zsad = zsad + zrhopad * ( 2. * zr4 * zs + zr2 &
1469	& + zr3 * zsr )
1470	zrhopad = 0.0_wp
1471
1472	ztad = ztad + zr3ad * (-2.1.6546e-6 zt &
1473	& + 1.0227e-4 )
1474	zr3ad = 0.0_wp
1475
1476	ztad = ztad + zr2ad * ( ( ( 4.5.3875e-9 zt &
1477	& -3.8.2467e-7 ) zt &
1478	& +2.7.6438e-5 ) zt &
1479	& - 4.0899e-3 )
1480	zr2ad = 0.0_wp
1481
1482	ztad = ztad + zr1ad * ( ( ( ( 5.6.536332e-9 zt &
1483	& -4.1.120083e-6 ) zt &
1484	& +3.1.001685e-4 ) zt &
1485	& -2.9.095290e-3 ) zt &
1486	& + 6.793952e-2 )
1487	zr1ad = 0.0_wp
1488
1489	zsad = zsad + zsrad * ( 1.0 / MAX( 2.*zsr, zeps ) ) &
1490	& * tmask(ji,jj,jk)
1491	zsrad = 0.0_wp
1492
1493	psal_ad(ji,jj,jk) = psal_ad(ji,jj,jk) + zsad
1494	ptem_ad(ji,jj,jk) = ptem_ad(ji,jj,jk) + ztad
1495	ztad = 0.0_wp
1496	zsad = 0.0_wp
1497	END DO
1498	END DO
1499	END DO
1500	!
1501	CASE ( 1 ) !== Linear formulation = F( temperature ) ==!
1502	DO jk = jpkm1, 1, -1
1503	prd_ad(:,:,jk) = prd_ad(:,:,jk) + rau0 * prhop_ad(:,:,jk) * tmask(:,:,jk)
1504	prhop_ad(:,:,jk) = 0.0_wp
1505	ptem_ad(:,:,jk) = ptem_ad(:,:,jk) - rn_alpha * prd_ad(:,:,jk) * tmask(:,:,jk)
1506	prd_ad(:,:,jk) = 0.0_wp
1507	END DO
1508	!
1509	CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==!
1510	DO jk = 1, jpkm1
1511	prd_ad( :,:,jk) = prd_ad(:,:,jk) + rau0 * prhop_ad(:,:,jk) * tmask(:,:,jk)
1512	prhop_ad(:,:,jk) = 0.0_wp
1513	ptem_ad( :,:,jk) = ptem_ad(:,:,jk) - rn_alpha * prd_ad(:,:,jk) * tmask(:,:,jk)
1514	psal_ad( :,:,jk) = psal_ad(:,:,jk) + rn_beta * prd_ad(:,:,jk) * tmask(:,:,jk)
1515	prd_ad( :,:,jk) = 0.0_wp
1516	END DO
1517	!
1518	END SELECT
1519
1520	END SUBROUTINE eos_insitu_pot_adj
1521
1522	SUBROUTINE eos_pot_1pt_adj( ptem, psal, ptem_ad, psal_ad, prhop_ad )
1523	REAL(wp), INTENT( in ) :: &
1524	ptem, & ! potential temperature
1525	psal ! salinity
1526	REAL(wp), INTENT( inout ) :: &
1527	ptem_ad, & ! potential temperature
1528	psal_ad ! salinity
1529	REAL(wp), INTENT( inout ) :: &
1530	prhop_ad ! potential density (surface referenced)
1531
1532	!! * Local declarations
1533	INTEGER :: ji, jj, jk ! dummy loop indices
1534	REAL(wp) :: & ! temporary scalars
1535	& zt, zs, zsr, zr1, zr2, zr3, zr4, zrhop, &
1536	& ztad, zsad, zsrad, zr1ad, zr2ad, zr3ad, zr4ad, zrhopad
1537	REAL(wp) :: zws, zwsad, zeps
1538	!!----------------------------------------------------------------------
1539
1540	#ifdef key_sp
1541	zeps = 1.e-7
1542	#else
1543	zeps = 1.e-14
1544	#endif
1545	zwsad = 0.0_wp
1546	ztad = 0.0_wp
1547	zsad = 0.0_wp
1548	zsrad = 0.0_wp
1549	zr1ad = 0.0_wp
1550	zr2ad = 0.0_wp
1551	zr3ad = 0.0_wp
1552	zr4ad = 0.0_wp
1553	zrhopad = 0.0_wp
1554	SELECT CASE ( nn_eos )
1555
1556	CASE ( 0 ) ! Jackett and McDougall (1994) formulation
1557	zws = SQRT( ABS( psal ) )
1558
1559	zt = ptem
1560	zs = psal
1561	! square root salinity
1562	zsr= zws
1563	! compute volumic mass pure water at atm pressure
1564	zr1= ( ( ( ( 6.536332e-9zt-1.120083e-6 )zt+1.001685e-4)*zt &
1565	-9.095290e-3 )zt+6.793952e-2 )zt+999.842594
1566	! seawater volumic mass atm pressure
1567	zr2= ( ( ( 5.3875e-9zt-8.2467e-7 ) zt+7.6438e-5 ) *zt &
1568	-4.0899e-3 ) *zt+0.824493
1569	zr3= ( -1.6546e-6zt+1.0227e-4 ) zt-5.72466e-3
1570	zr4= 4.8314e-4
1571
1572	! ============
1573	! Adjoint part
1574	! ============
1575	! save potential volumic mass
1576	zrhopad = zrhopad + prhop_ad
1577	prhop_ad = 0.0_wp
1578	! potential volumic mass (reference to the surface)
1579	zr1ad = zr1ad + zrhopad
1580	zr2ad = zr2ad + zrhopad * zs
1581	zr3ad = zr3ad + zrhopad * zsr * zs
1582	zsrad = zsrad + zrhopad * zr3 * zs
1583	zsad = zsad + zrhopad * ( 2. * zr4 * zs + zr2 &
1584	& + zr3 * zsr )
1585	zrhopad= 0.0_wp
1586
1587	! seawater volumic mass atm pressure
1588	ztad = ztad + zr3ad * ( -2.1.6546e-6 zt &
1589	& + 1.0227e-4 )
1590	zr3ad = 0.0_wp
1591	ztad = ztad + zr2ad * ( ( ( 4.5.3875e-9 zt &
1592	& -3.8.2467e-7 ) zt &
1593	& +2.7.6438e-5 ) zt &
1594	& - 4.0899e-3 )
1595	zr2ad = 0.0_wp
1596	! compute volumic mass pure water at atm pressure
1597	ztad = ztad + zr1ad * ( ( ( ( 5.6.536332e-9 zt &
1598	& -4.1.120083e-6 ) zt &
1599	& +3.1.001685e-4 ) zt &
1600	& -2.9.095290e-3 ) zt &
1601	& + 6.793952e-2 )
1602	zr1ad = 0.0_wp
1603	! square root salinity
1604	zwsad = zwsad + zsrad
1605	zsrad = 0.0_wp
1606
1607	ptem_ad = ptem_ad + ztad
1608	psal_ad = psal_ad + zsad
1609	ztad = 0.0_wp
1610	zsad = 0.0_wp
1611
1612	psal_ad = psal_ad + 1/MAX( 2zws , zeps ) zwsad
1613	zwsad = 0.0_wp
1614
1615	CASE ( 1 ) ! Linear formulation function of temperature only
1616
1617	! ... potential volumic mass
1618	ptem_ad = ptem_ad - prhop_ad * rau0 * rn_alpha
1619	prhop_ad = 0.0_wp
1620
1621	CASE ( 2 ) ! Linear formulation function of temperature and salinity
1622
1623	! ... potential volumic mass
1624	ptem_ad = ptem_ad - prhop_ad * rau0 * rn_alpha
1625	psal_ad = psal_ad + prhop_ad * rau0 * rn_beta
1626	prhop_ad = 0.0_wp
1627
1628	CASE DEFAULT
1629
1630	WRITE(ctmp1,*) ' bad flag value for nn_eos = ', nn_eos
1631	CALL ctl_stop( ctmp1 )
1632
1633	END SELECT
1634
1635	END SUBROUTINE eos_pot_1pt_adj
1636	SUBROUTINE eos_insitu_2d_adj( ptem, psal, pdep, ptem_ad, psal_ad, prd_ad )
1637	!!-----------------------------------------------------------------------
1638	!!
1639	!! * ROUTINE eos_insitu_2d_adj : adj OF ROUTINE eos_insitu_2d *
1640	!!
1641	!! ** Purpose of direct routine : Compute the in situ density
1642	!! (ratio rho/rau0) from potential temperature and salinity
1643	!! using an equation of state defined through the namelist
1644	!! parameter nn_eos. * 2D field case
1645	!!
1646	!! ** Method of direct routine : 3 cases:
1647	!! nn_eos = 0 : Jackett and McDougall (1994) equation of state.
1648	!! the in situ density is computed directly as a function of
1649	!! potential temperature relative to the surface (the opa t
1650	!! variable), salt and pressure (assuming no pressure variation
1651	!! along geopotential surfaces, i.e. the pressure p in decibars
1652	!! is approximated by the depth in meters.
1653	!! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0
1654	!! with pressure p decibars
1655	!! potential temperature t deg celsius
1656	!! salinity s psu
1657	!! reference volumic mass rau0 kg/m**3
1658	!! in situ volumic mass rho kg/m**3
1659	!! in situ density anomalie prd no units
1660	!! Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar,
1661	!! t = 40 deg celcius, s=40 psu
1662	!! nn_eos = 1 : linear equation of state function of temperature only
1663	!! prd(t) = 0.0285 - rn_alpha * t
1664	!! nn_eos = 2 : linear equation of state function of temperature and
1665	!! salinity
1666	!! prd(t,s) = rn_beta * s - rn_alpha * tn - 1.
1667	!! Note that no boundary condition problem occurs in this routine
1668	!! as (ptem,psal) are defined over the whole domain.
1669	!!
1670	!! ** Comments on Adjoint Routine :
1671	!! Care has been taken to avoid division by zero when computing
1672	!! the inverse of the square root of salinity at masked salinity
1673	!! points.
1674	!!
1675	!! ** Action :
1676	!!
1677	!! References :
1678	!!
1679	!! History :
1680	!! 8.2 ! 05-03 ((F. Van den Berghe, A. Weaver, N. Daget) - eosadj.F
1681	!! 9.0 ! 08-07 (A. Vidard) first version based on eosadj
1682	!!-----------------------------------------------------------------------
1683	!! * Modules used
1684	!! * Arguments
1685	REAL(wp), DIMENSION(jpi,jpj), INTENT( in ) :: &
1686	& ptem, & ! potential temperature
1687	& psal, & ! salinity
1688	& pdep ! depth
1689	REAL(wp), DIMENSION(jpi,jpj), INTENT( inout ) :: &
1690	& ptem_ad, & ! TL of potential temperature
1691	& psal_ad ! TL of salinity
1692	REAL(wp), DIMENSION(jpi,jpj), INTENT( inout ) :: &
1693	& prd_ad ! TL of potential density (surface referenced)
1694	INTEGER :: ji, jj ! dummy loop indices
1695	REAL(wp) :: & ! temporary scalars
1696	zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw, &
1697	zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, &
1698	ztad, zsad, zhad, zsrad, zr1ad, zr2ad, zr3ad, &
1699	zr4ad, zrhopad, zead, zbwad, &
1700	zbad, zdad, zcad, zawad, zaad, zb1ad, za1ad, &
1701	zkwad, zk0ad, zpes, zrdc1, zrdc2, zeps, &
1702	zmask
1703	REAL(wp), DIMENSION(jpi,jpj) :: zws
1704	!!----------------------------------------------------------------------
1705
1706	! initialization of adjoint variables
1707	ztad = 0.0_wp
1708	zsad = 0.0_wp
1709	zhad = 0.0_wp
1710	zsrad = 0.0_wp
1711	zr1ad = 0.0_wp
1712	zr2ad = 0.0_wp
1713	zr3ad = 0.0_wp
1714	zr4ad = 0.0_wp
1715	zrhopad = 0.0_wp
1716	zead = 0.0_wp
1717	zbwad = 0.0_wp
1718	zbad = 0.0_wp
1719	zdad = 0.0_wp
1720	zcad = 0.0_wp
1721	zawad = 0.0_wp
1722	zaad = 0.0_wp
1723	zb1ad = 0.0_wp
1724	za1ad = 0.0_wp
1725	zkwad = 0.0_wp
1726	zk0ad = 0.0_wp
1727
1728
1729	SELECT CASE ( nn_eos )
1730
1731	CASE ( 0 ) !== Jackett and McDougall (1994) formulation ==!
1732
1733	#ifdef key_sp
1734	zeps = 1.e-7
1735	#else
1736	zeps = 1.e-14
1737	#endif
1738
1739	!CDIR NOVERRCHK
1740	DO jj = 1, jpjm1
1741	!CDIR NOVERRCHK
1742	DO ji = 1, fs_jpim1 ! vector opt.
1743	zws(ji,jj) = SQRT( ABS( psal(ji,jj) ) )
1744	END DO
1745	END DO
1746	!
1747	DO jj = 1, jpjm1
1748	DO ji = 1, fs_jpim1 ! vector opt.
1749	zmask = tmask(ji,jj,1) ! land/sea bottom mask = surf. mask
1750	zt = ptem (ji,jj) ! interpolated T
1751	zs = psal (ji,jj) ! interpolated S
1752	zsr= zws(ji,jj) ! square root of interpolated S
1753	zh = pdep(ji,jj) ! depth at the partial step level
1754	! compute volumic mass pure water at atm pressure
1755	zr1= ( ( ( ( 6.536332e-9zt-1.120083e-6 )zt+1.001685e-4)*zt &
1756	& -9.095290e-3 )zt+6.793952e-2 )zt+999.842594
1757	! seawater volumic mass atm pressure
1758	zr2= ( ( ( 5.3875e-9zt-8.2467e-7 ) zt+7.6438e-5 ) *zt &
1759	& -4.0899e-3 ) *zt+0.824493
1760	zr3= ( -1.6546e-6zt+1.0227e-4 ) zt-5.72466e-3
1761	zr4= 4.8314e-4
1762
1763	! potential volumic mass (reference to the surface)
1764	zrhop= ( zr4zs + zr3zsr + zr2 ) *zs + zr1
1765
1766	! add the compression terms
1767	ze = ( -3.508914e-8zt-1.248266e-8 ) zt-2.595994e-6
1768	zbw= ( 1.296821e-6zt-5.782165e-9 ) zt+1.045941e-4
1769	zb = zbw + ze * zs
1770
1771	zd = -2.042967e-2
1772	zc = (-7.267926e-5zt+2.598241e-3 ) zt+0.1571896
1773	zaw= ( ( 5.939910e-6zt+2.512549e-3 ) zt-0.1028859 ) *zt - 4.721788
1774	za = ( zdzsr + zc ) zs + zaw
1775
1776	zb1= (-0.1909078zt+7.390729 ) zt-55.87545
1777	za1= ( ( 2.326469e-3zt+1.553190)zt-65.00517 ) *zt+1044.077
1778	zkw= ( ( (-1.361629e-4zt-1.852732e-2 ) zt-30.41638 ) zt + 2098.925 ) zt+190925.6
1779	zk0= ( zb1zsr + za1 )zs + zkw
1780
1781	zrdc1 = 1.0 / ( zk0 - zh * ( za - zh * zb ) )
1782	zrdc2 = 1.0 / ( 1.0 - zh * zrdc1 )
1783	! ============
1784	! Adjoint part
1785	! ============
1786
1787	! Masked in situ density anomaly
1788
1789	zrhopad = zrhopad + prd_ad(ji,jj) * tmask(ji,jj,1) &
1790	& * zrdc2 / rau0
1791	zk0ad = zk0ad - prd_ad(ji,jj) * tmask(ji,jj,1) &
1792	& * zrdc2 * zrdc2 * zh &
1793	& * zrdc1*2 zrhop &
1794	& / rau0
1795	zaad = zaad + prd_ad(ji,jj) * tmask(ji,jj,1) &
1796	& * zrdc2 * zrdc2 * zh &
1797	& * zrdc1*2 zrhop &
1798	& * zh / rau0
1799	zbad = zbad - prd_ad(ji,jj) * tmask(ji,jj,1) &
1800	& * zrdc2 * zrdc2 * zh &
1801	& * zrdc1*2 zrhop &
1802	& * zh * zh / rau0
1803	prd_ad(ji,jj) = 0.0_wp
1804
1805	zkwad = zkwad + zk0ad
1806	zsrad = zsrad + zk0ad * zb1 * zs
1807	zb1ad = zb1ad + zk0ad * zs * zsr
1808	za1ad = za1ad + zk0ad * zs
1809	zsad = zsad + zk0ad * ( zb1 * zsr + za1 )
1810	zk0ad = 0.0_wp
1811
1812	ztad = ztad + zkwad * ( ( (-4.1.361629e-4 zt &
1813	& -3.1.852732e-2 ) zt &
1814	& -2.30.41638 ) zt &
1815	& + 2098.925 )
1816	zkwad = 0.0_wp
1817
1818	ztad = ztad + za1ad * ( ( 3.2.326469e-3 zt &
1819	& +2.1.553190 ) zt &
1820	& - 65.00517 )
1821	za1ad = 0.0_wp
1822
1823	ztad = ztad + zb1ad * (-2.0.1909078 zt &
1824	& + 7.390729 )
1825	zb1ad = 0.0_wp
1826
1827	zawad = zawad + zaad
1828	zsrad = zsrad + zaad * zd * zs
1829	zcad = zcad + zaad * zs
1830	zsad = zsad + zaad * ( zd * zsr + zc )
1831	zaad = 0.0_wp
1832
1833	ztad = ztad + zawad * ( ( 3.5.939910e-6 zt &
1834	& +2.2.512549e-3 ) zt &
1835	& - 0.1028859 )
1836	zawad = 0.0_wp
1837
1838	ztad = ztad + zcad * (-2.7.267926e-5 zt &
1839	& + 2.598241e-3 )
1840	zcad = 0.0_wp
1841
1842	zbwad = zbwad + zbad
1843	zead = zead + zbad * zs
1844	zsad = zsad + zbad * ze
1845	zbad = 0.0_wp
1846
1847	ztad = ztad + zbwad * ( 2.1.296821e-6 zt &
1848	& - 5.782165e-9 )
1849	zbwad = 0.0_wp
1850
1851	ztad = ztad + zead * (-2.3.508914e-8 zt &
1852	& - 1.248266e-8 )
1853	zead = 0.0_wp
1854
1855	zr1ad = zr1ad + zrhopad
1856	zr2ad = zr2ad + zrhopad * zs
1857	zr3ad = zr3ad + zrhopad * zsr * zs
1858	zsrad = zsrad + zrhopad * zr3 * zs
1859	zsad = zsad + zrhopad * ( 2. * zr4 * zs + zr2 &
1860	& + zr3 * zsr )
1861	zrhopad = 0.0_wp
1862
1863	ztad = ztad + zr3ad * (-2.1.6546e-6 zt &
1864	& + 1.0227e-4 )
1865	zr3ad = 0.0_wp
1866
1867	ztad = ztad + zr2ad * ( ( ( 4.5.3875e-9 zt &
1868	& -3.8.2467e-7 ) zt &
1869	& +2.7.6438e-5 ) zt &
1870	& - 4.0899e-3 )
1871	zr2ad = 0.0_wp
1872
1873	ztad = ztad + zr1ad * ( ( ( ( 5.6.536332e-9 zt &
1874	& -4.1.120083e-6 ) zt &
1875	& +3.1.001685e-4 ) zt &
1876	& -2.9.095290e-3 ) zt &
1877	& + 6.793952e-2 )
1878	zr1ad = 0.0_wp
1879
1880	zsad = zsad + zsrad * ( 1.0 / MAX( 2.*zsr, zeps ) ) &
1881	& * tmask(ji,jj, 1)
1882	zsrad = 0.0_wp
1883
1884	psal_ad(ji,jj) = psal_ad(ji,jj) + zsad
1885	ptem_ad(ji,jj) = ptem_ad(ji,jj) + ztad
1886	ztad = 0.0_wp
1887	zsad = 0.0_wp
1888	END DO
1889	END DO
1890	!
1891	CASE ( 1 ) !== Linear formulation = F( temperature ) ==!
1892	DO jj = 1, jpjm1
1893	DO ji = 1, fs_jpim1 ! vector opt.
1894	ptem_ad(ji,jj) = ptem_ad(ji,jj) - prd_ad(ji,jj) * rn_alpha * tmask(ji,jj,1)
1895	prd_ad(ji,jj) = 0.0_wp
1896	END DO
1897	END DO
1898	!
1899	CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==!
1900	DO jj = 1, jpjm1
1901	DO ji = 1, fs_jpim1 ! vector opt.
1902	ptem_ad(ji,jj) = ptem_ad(ji,jj) - prd_ad(ji,jj) * rn_alpha * tmask(ji,jj,1)
1903	psal_ad(ji,jj) = psal_ad(ji,jj) + prd_ad(ji,jj) * rn_beta * tmask(ji,jj,1)
1904	prd_ad (ji,jj) = 0.0_wp
1905	END DO
1906	END DO
1907	!
1908	END SELECT
1909	END SUBROUTINE eos_insitu_2d_adj
1910
1911	SUBROUTINE eos_bn2_tan ( ptem, psal, ptem_tl, psal_tl, pn2_tl )
1912	!!----------------------------------------------------------------------
1913	!! * ROUTINE eos_bn2_tan *
1914	!!
1915	!! ** Purpose of the direct routine: Compute the local
1916	!! Brunt-Vaisala frequency at the time-step of the input arguments
1917	!!
1918	!! ** Method of the direct routine:
1919	!! * nn_eos = 0 : UNESCO sea water properties
1920	!! The brunt-vaisala frequency is computed using the polynomial
1921	!! polynomial expression of McDougall (1987):
1922	!! N^2 = grav * beta * ( alpha/beta*dk[ t ] - dk[ s ] )/e3w
1923	!! If lk_zdfddm=T, the heat/salt buoyancy flux ratio Rrau is
1924	!! computed and used in zdfddm module :
1925	!! Rrau = alpha/beta * ( dk[ t ] / dk[ s ] )
1926	!! * nn_eos = 1 : linear equation of state (temperature only)
1927	!! N^2 = grav * rn_alpha * dk[ t ]/e3w
1928	!! * nn_eos = 2 : linear equation of state (temperature & salinity)
1929	!! N^2 = grav * (rn_alpha * dk[ t ] - rn_beta * dk[ s ] ) / e3w
1930	!! The use of potential density to compute N^2 introduces e r r o r
1931	!! in the sign of N^2 at great depths. We recommand the use of
1932	!! nn_eos = 0, except for academical studies.
1933	!! Macro-tasked on horizontal slab (jk-loop)
1934	!! N.B. N^2 is set to zero at the first level (JK=1) in inidtr
1935	!! and is never used at this level.
1936	!!
1937	!! ** Action : - pn2 : the brunt-vaisala frequency
1938	!!
1939	!! References :
1940	!! McDougall, T. J., J. Phys. Oceanogr., 17, 1950-1964, 1987.
1941	!!
1942	!! History:
1943	!! ! 08-07 (A. Vidard) First version
1944	!!----------------------------------------------------------------------
1945	!! * Arguments
1946
1947	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: &
1948	ptem, & ! potential temperature
1949	psal ! salinity
1950	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: &
1951	ptem_tl, & ! potential temperature
1952	psal_tl ! salinity
1953	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out ) :: &
1954	pn2_tl ! Brunt-Vaisala frequency
1955
1956	!! * Local declarations
1957	INTEGER :: ji, jj, jk ! dummy loop indices
1958	REAL(wp) :: &
1959	zgde3w, zt, zs, zh, & ! temporary scalars
1960	zalbet, zbeta ! " "
1961	REAL(wp) :: &
1962	zttl, zstl, & ! temporary scalars
1963	zalbettl, zbetatl ! " "
1964	#if defined key_zdfddm
1965	REAL(wp) :: zds, zdstl ! temporary scalars
1966	#endif
1967
1968	! pn2_tl : interior points only (2=< jk =< jpkm1 )
1969	! --------------------------
1970
1971	SELECT CASE ( nn_eos )
1972
1973	CASE ( 0 ) !== Jackett and McDougall (1994) formulation ==!
1974	DO jk = 2, jpkm1
1975	DO jj = 1, jpj
1976	DO ji = 1, jpi
1977	zgde3w = grav / fse3w(ji,jj,jk)
1978	zt = 0.5 * ( ptem(ji,jj,jk) + ptem(ji,jj,jk-1) ) ! potential temperature at w-point
1979	zs = 0.5 * ( psal(ji,jj,jk) + psal(ji,jj,jk-1) ) - 35.0 ! salinity anomaly (s-35) at w-point
1980	zh = fsdepw(ji,jj,jk) ! depth in meters at w-point
1981
1982	zalbet = ( ( ( - 0.255019e-07 * zt + 0.298357e-05 ) * zt & ! ratio alpha/beta
1983	& - 0.203814e-03 ) * zt &
1984	& + 0.170907e-01 ) * zt &
1985	& + 0.665157e-01 &
1986	& + ( - 0.678662e-05 * zs &
1987	& - 0.846960e-04 * zt + 0.378110e-02 ) * zs &
1988	& + ( ( - 0.302285e-13 * zh &
1989	& - 0.251520e-11 * zs &
1990	& + 0.512857e-12 * zt * zt ) * zh &
1991	& - 0.164759e-06 * zs &
1992	& +( 0.791325e-08 * zt - 0.933746e-06 ) * zt &
1993	& + 0.380374e-04 ) * zh
1994
1995	zbeta = ( ( -0.415613e-09 * zt + 0.555579e-07 ) * zt & ! beta
1996	& - 0.301985e-05 ) * zt &
1997	& + 0.785567e-03 &
1998	& + ( 0.515032e-08 * zs &
1999	& + 0.788212e-08 * zt - 0.356603e-06 ) * zs &
2000	& +( ( 0.121551e-17 * zh &
2001	& - 0.602281e-15 * zs &
2002	& - 0.175379e-14 * zt + 0.176621e-12 ) * zh &
2003	& + 0.408195e-10 * zs &
2004	& + ( - 0.213127e-11 * zt + 0.192867e-09 ) * zt &
2005	& - 0.121555e-07 ) * zh
2006
2007
2008	!! tangent part
2009	zttl = 0.5 * ( ptem_tl(ji,jj,jk) + ptem_tl(ji,jj,jk-1) ) ! potential temperature at w-point
2010	zstl = 0.5 * ( psal_tl(ji,jj,jk) + psal_tl(ji,jj,jk-1) ) ! salinity anomaly at w-point
2011	zalbettl = ( ( ( -4.0.255019e-07 zt &! ratio alpha/beta
2012	& +3.0.298357e-05 ) zt &
2013	& -2.0.203814e-03 ) zt &
2014	& + 0.170907e-01 &
2015	& - 0.846960e-04 * zs &
2016	& - ( 0.933746e-06 &
2017	& - ( 2.*0.791325e-08 &
2018	& +2.0.512857e-12 zh ) * zt ) * zh ) * zttl &
2019	& + ( - 2.0.678662e-05 zs &
2020	& - 0.846960e-04 * zt &
2021	& + 0.378110e-02 &
2022	& + ( - 0.164759e-06 &
2023	& - 0.251520e-11 * zh ) * zh ) * zstl
2024
2025	zbetatl = ( ( -3.0.415613e-09 zt &
2026	& +2.0.555579e-07 ) zt &
2027	& - 0.301985e-05 &
2028	& + 0.788212e-08 * zs &
2029	& + ( -2.0.213127e-11 zt &
2030	& - 0.175379e-14 * zh &
2031	& + 0.192867e-09 ) * zh ) * zttl &
2032	& + ( 2.0.515032e-08 zs &
2033	& + 0.788212e-08 * zt &
2034	& - 0.356603e-06 &
2035	& + ( - 0.602281e-15 * zh &
2036	& + 0.408195e-10 ) * zh ) * zstl
2037
2038	pn2_tl(ji,jj,jk) = zgde3w * tmask(ji,jj,jk) * ( &
2039	& zbeta * ( zalbet &
2040	& * ( ptem_tl(ji,jj,jk-1) - ptem_tl(ji,jj,jk) ) &
2041	& + zalbettl &
2042	& * ( ptem (ji,jj,jk-1) - ptem (ji,jj,jk) ) &
2043	& - ( psal_tl(ji,jj,jk-1) - psal_tl(ji,jj,jk) ) ) &
2044	& + zbetatl * ( zalbet &
2045	& * ( ptem (ji,jj,jk-1) - ptem (ji,jj,jk) ) &
2046	& - ( psal (ji,jj,jk-1) - psal (ji,jj,jk) ) ) )
2047	#if defined key_zdfddm
2048	zds = ( psal(ji,jj,jk-1) - psal(ji,jj,jk) )
2049	zdstl = ( psal_tl(ji,jj,jk-1) - psal_tl(ji,jj,jk) )
2050	IF ( ABS( zds) <= 1.e-20 ) THEN
2051	zds = 1.e-20
2052	rrau_tl(ji,jj,jk) = zalbettl * &
2053	& ( ptem(ji,jj,jk-1) - ptem(ji,jj,jk) ) / zds &
2054	& + zalbet * &
2055	& ( ( ptem_tl(ji,jj,jk-1) - ptem_tl(ji,jj,jk) ) / zds )
2056	ELSE
2057	rrau_tl(ji,jj,jk) = zalbettl * &
2058	& ( ptem(ji,jj,jk-1) - ptem(ji,jj,jk) ) / zds &
2059	& + zalbet * &
2060	& ( ( ptem_tl(ji,jj,jk-1) - ptem_tl(ji,jj,jk) ) / zds &
2061	& - ( ptem(ji,jj,jk-1) - ptem(ji,jj,jk) ) * zdstl/ zds**2 )
2062	ENDIF
2063	#endif
2064	END DO
2065	END DO
2066	END DO
2067	!
2068	CASE ( 1 ) !== Linear formulation = F( temperature ) ==!
2069	DO jk = 2, jpkm1
2070	pn2_tl(:,:,jk) = rn_alpha * ( ptem_tl(:,:,jk-1) - ptem_tl(:,:,jk) ) * grav / fse3w(:,:,jk) * tmask(:,:,jk)
2071	END DO
2072	!
2073	CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==!
2074	DO jk = 2, jpkm1
2075	pn2_tl(:,:,jk) = ( rn_alpha * ( ptem_tl(:,:,jk-1) - ptem_tl(:,:,jk) ) &
2076	& - rn_beta * ( psal_tl(:,:,jk-1) - psal_tl(:,:,jk) ) ) &
2077	& * grav / fse3w(:,:,jk) * tmask(:,:,jk)
2078	END DO
2079	#if defined key_zdfddm
2080	DO jk = 2, jpkm1
2081	DO jj = 1, jpj
2082	DO ji = 1, jpi
2083	zds = ( psal(ji,jj,jk-1) - psal(ji,jj,jk) )
2084	zdstl = psal_tl(ji,jj,jk-1) - psal_tl(ji,jj,jk)
2085	IF ( ABS( zds) <= 1.e-20 ) THEN
2086	zds = 1.e-20
2087	rrau_tl(ji,jj,jk) = ralpbet * &
2088	& ( ( ptem_tl(ji,jj,jk-1) - ptem_tl(ji,jj,jk) ) / zds )
2089	ELSE
2090	rrau_tl(ji,jj,jk) = ralpbet * &
2091	& ( ( ptem_tl(ji,jj,jk-1) - ptem_tl(ji,jj,jk) ) / zds &
2092	& - ( ptem(ji,jj,jk-1) - ptem(ji,jj,jk) ) * zdstl / zds**2 )
2093	ENDIF
2094	rrau(ji,jj,jk) = ralpbet * ( ptem(ji,jj,jk-1) - ptem(ji,jj,jk) ) / zds
2095	END DO
2096	END DO
2097	END DO
2098	#endif
2099	END SELECT
2100
2101	END SUBROUTINE eos_bn2_tan
2102
2103	SUBROUTINE eos_bn2_adj ( ptem, psal, ptem_ad, psal_ad, pn2_ad )
2104	!!----------------------------------------------------------------------
2105	!! * ROUTINE eos_bn2_adj *
2106	!!
2107	!! ** Purpose of the direct routine: Compute the local
2108	!! Brunt-Vaisala frequency at the time-step of the input arguments
2109	!!
2110	!! ** Method of the direct routine:
2111	!! * nn_eos = 0 : UNESCO sea water properties
2112	!! The brunt-vaisala frequency is computed using the polynomial
2113	!! polynomial expression of McDougall (1987):
2114	!! N^2 = grav * beta * ( alpha/beta*dk[ t ] - dk[ s ] )/e3w
2115	!! If lk_zdfddm=T, the heat/salt buoyancy flux ratio Rrau is
2116	!! computed and used in zdfddm module :
2117	!! Rrau = alpha/beta * ( dk[ t ] / dk[ s ] )
2118	!! * nn_eos = 1 : linear equation of state (temperature only)
2119	!! N^2 = grav * rn_alpha * dk[ t ]/e3w
2120	!! * nn_eos = 2 : linear equation of state (temperature & salinity)
2121	!! N^2 = grav * (rn_alpha * dk[ t ] - rn_beta * dk[ s ] ) / e3w
2122	!! The use of potential density to compute N^2 introduces e r r o r
2123	!! in the sign of N^2 at great depths. We recommand the use of
2124	!! nn_eos = 0, except for academical studies.
2125	!! Macro-tasked on horizontal slab (jk-loop)
2126	!! N.B. N^2 is set to zero at the first level (JK=1) in inidtr
2127	!! and is never used at this level.
2128	!!
2129	!! ** Action : - pn2 : the brunt-vaisala frequency
2130	!!
2131	!! References :
2132	!! McDougall, T. J., J. Phys. Oceanogr., 17, 1950-1964, 1987.
2133	!!
2134	!! History:
2135	!! ! 08-07 (A. Vidard) First version
2136	!!----------------------------------------------------------------------
2137	!! * Arguments
2138
2139	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: &
2140	ptem, & ! potential temperature
2141	psal ! salinity
2142	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: &
2143	ptem_ad, & ! adjoint potential temperature
2144	psal_ad ! adjoint salinity
2145	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: &
2146	pn2_ad ! adjoint Brunt-Vaisala frequency
2147
2148	!! * Local declarations
2149	INTEGER :: ji, jj, jk ! dummy loop indices
2150	REAL(wp) :: &
2151	zgde3w, zt, zs, zh, & ! temporary scalars
2152	zalbet, zbeta ! " "
2153	REAL(wp) :: &
2154	ztad, zsad, & ! temporary scalars
2155	zalbetad, zbetaad ! " "
2156	#if defined key_zdfddm
2157	REAL(wp) :: zds, zdsad ! temporary scalars
2158	#endif
2159
2160	! pn2_tl : interior points only (2=< jk =< jpkm1 )
2161	! --------------------------
2162	zalbetad = 0.0_wp
2163	zbetaad = 0.0_wp
2164	ztad = 0.0_wp
2165	zsad = 0.0_wp
2166	#if defined key_zdfddm
2167	zdsad = 0.0_wp
2168	#endif
2169
2170	SELECT CASE ( nn_eos )
2171
2172	CASE ( 0 ) !== Jackett and McDougall (1994) formulation ==!
2173	DO jk = jpkm1, 2, -1
2174	DO jj = jpj, 1, -1
2175	DO ji = jpi, 1, -1
2176	zgde3w = grav / fse3w(ji,jj,jk)
2177	zt = 0.5 * ( ptem(ji,jj,jk) + ptem(ji,jj,jk-1) ) ! potential temperature at w-point
2178	zs = 0.5 * ( psal(ji,jj,jk) + psal(ji,jj,jk-1) ) - 35.0 ! salinity anomaly (s-35) at w-point
2179	zh = fsdepw(ji,jj,jk) ! depth in meters at w-point
2180
2181	zalbet = ( ( ( - 0.255019e-07 * zt + 0.298357e-05 ) * zt & ! ratio alpha/beta
2182	& - 0.203814e-03 ) * zt &
2183	& + 0.170907e-01 ) * zt &
2184	& + 0.665157e-01 &
2185	& + ( - 0.678662e-05 * zs &
2186	& - 0.846960e-04 * zt + 0.378110e-02 ) * zs &
2187	& + ( ( - 0.302285e-13 * zh &
2188	& - 0.251520e-11 * zs &
2189	& + 0.512857e-12 * zt * zt ) * zh &
2190	& - 0.164759e-06 * zs &
2191	& +( 0.791325e-08 * zt - 0.933746e-06 ) * zt &
2192	& + 0.380374e-04 ) * zh
2193
2194	zbeta = ( ( -0.415613e-09 * zt + 0.555579e-07 ) * zt & ! beta
2195	& - 0.301985e-05 ) * zt &
2196	& + 0.785567e-03 &
2197	& + ( 0.515032e-08 * zs &
2198	& + 0.788212e-08 * zt - 0.356603e-06 ) * zs &
2199	& +( ( 0.121551e-17 * zh &
2200	& - 0.602281e-15 * zs &
2201	& - 0.175379e-14 * zt + 0.176621e-12 ) * zh &
2202	& + 0.408195e-10 * zs &
2203	& + ( - 0.213127e-11 * zt + 0.192867e-09 ) * zt &
2204	& - 0.121555e-07 ) * zh
2205
2206
2207	#if defined key_zdfddm
2208
2209	zds = ( psal(ji,jj,jk-1) - psal(ji,jj,jk) )
2210	IF ( ABS( zds) <= 1.e-20 ) THEN
2211	zds = 1.e-20
2212	zdsad = 0.0_wp
2213	ELSE
2214	zdsad = rrau_ad(ji,jj,jk) * zalbet ( ptem(ji,jj,jk-1) - ptem(ji,jj,jk) ) / zds*2
2215	ENDIF
2216	ptem_ad(ji,jj,jk-1) = ptem_ad(ji,jj,jk-1) + rrau_ad(ji,jj,jk) * zalbet / zds
2217	ptem_ad(ji,jj,jk ) = ptem_ad(ji,jj,jk ) - rrau_ad(ji,jj,jk) * zalbet / zds
2218	zalbetad = zalbetad + rrau_ad(ji,jj,jk) * ( ptem(ji,jj,jk-1) - ptem(ji,jj,jk) ) / zds
2219	rrau_ad(ji,jj,jk) = 0._wp
2220
2221	psal_ad(ji,jj,jk-1) = psal_ad(ji,jj,jk-1) + zdsad
2222	psal_ad(ji,jj,jk ) = psal_ad(ji,jj,jk ) - zdsad
2223	zdsad = 0._wp
2224	#endif
2225	ptem_ad(ji,jj,jk-1) = ptem_ad(ji,jj,jk-1) + zalbetzbetazgde3wtmask(ji,jj,jk)pn2_ad(ji,jj,jk)
2226	ptem_ad(ji,jj,jk) = ptem_ad(ji,jj,jk ) - zalbetzbetazgde3wtmask(ji,jj,jk)pn2_ad(ji,jj,jk)
2227	zalbetad = zalbetad + zbetazgde3wtmask(ji,jj,jk)( ptem (ji,jj,jk-1) - ptem (ji,jj,jk) ) pn2_ad(ji,jj,jk)
2228	psal_ad(ji,jj,jk-1) = psal_ad(ji,jj,jk-1) - zbetatmask(ji,jj,jk)zgde3w*pn2_ad(ji,jj,jk)
2229	psal_ad(ji,jj,jk ) = psal_ad(ji,jj,jk ) + zbetatmask(ji,jj,jk)zgde3w*pn2_ad(ji,jj,jk)
2230	zbetaad = zbetaad &
2231	& + zgde3w tmask(ji,jj,jk) ( zalbet * ( ptem (ji,jj,jk-1) - ptem (ji,jj,jk) ) &
2232	& - ( psal (ji,jj,jk-1) - psal (ji,jj,jk) ) )*pn2_ad(ji,jj,jk)
2233
2234	pn2_ad(ji,jj,jk) = 0.0_wp
2235
2236	ztad = ztad + ( ( -3.0.415613e-09 zt &
2237	& +2.0.555579e-07 ) zt &
2238	& - 0.301985e-05 &
2239	& + 0.788212e-08 * zs &
2240	& + ( -2.0.213127e-11 zt &
2241	& - 0.175379e-14 * zh &
2242	& + 0.192867e-09 ) * zh ) *zbetaad
2243
2244	zsad = zsad + ( 2.0.515032e-08 zs &
2245	& + 0.788212e-08 * zt &
2246	& - 0.356603e-06 &
2247	& + ( - 0.602281e-15 * zh &
2248	& + 0.408195e-10 ) * zh ) * zbetaad
2249
2250	zbetaad = 0.0_wp
2251
2252	ztad = ztad + ( ( ( -4.0.255019e-07 zt &! ratio alpha/beta
2253	& +3.0.298357e-05 ) zt &
2254	& -2.0.203814e-03 ) zt &
2255	& + 0.170907e-01 &
2256	& - 0.846960e-04 * zs &
2257	& - ( 0.933746e-06 &
2258	& - ( 2.*0.791325e-08 &
2259	& +2.0.512857e-12 zh ) * zt ) * zh &
2260	& ) *zalbetad
2261
2262	zsad = zsad + ( - 2.0.678662e-05 zs &
2263	& - 0.846960e-04 * zt &
2264	& + 0.378110e-02 &
2265	& + ( - 0.164759e-06 &
2266	& - 0.251520e-11 * zh ) * zh &
2267	& ) *zalbetad
2268
2269	zalbetad = 0.0_wp
2270
2271
2272	psal_ad(ji,jj,jk) = psal_ad(ji,jj,jk) + 0.5 * zsad
2273	psal_ad(ji,jj,jk-1) = psal_ad(ji,jj,jk-1) + 0.5 * zsad
2274	zsad = 0.0_wp
2275
2276	ptem_ad(ji,jj,jk) = ptem_ad(ji,jj,jk) + 0.5 * ztad
2277	ptem_ad(ji,jj,jk-1) = ptem_ad(ji,jj,jk-1) + 0.5 * ztad
2278	ztad = 0.0_wp
2279
2280	END DO
2281	END DO
2282	END DO
2283	!
2284	CASE ( 1 ) !== Linear formulation = F( temperature ) ==!
2285	DO jk = jpkm1, 2, -1
2286	ptem_ad(:,:,jk-1) = ptem_ad(:,:,jk-1) + rn_alpha * pn2_ad(:,:,jk) &
2287	& * grav / fse3w(:,:,jk) * tmask(:,:,jk)
2288	ptem_ad(:,:,jk ) = ptem_ad(:,:,jk ) - rn_alpha * pn2_ad(:,:,jk) &
2289	& * grav / fse3w(:,:,jk) * tmask(:,:,jk)
2290	pn2_ad(:,:,jk) = 0.0_wp
2291	END DO
2292	!
2293	CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==!
2294	#if defined key_zdfddm
2295	DO jk = jpkm1, 2, -1
2296	DO jj = jpj, 1, -1
2297	DO ji = jpi, 1, -1
2298	zds = ( psal(ji,jj,jk-1) - psal(ji,jj,jk) )
2299	IF ( ABS( zds) <= 1.e-20 ) THEN
2300	zds = 1.e-20
2301	zdsad = 0.0_wp
2302	ELSE
2303	zdsad = zdsad - rrau_ad(ji,jj,jk) * ralpbet &
2304	& * ( ptem(ji,jj,jk-1) - ptem(ji,jj,jk) ) / zds**2
2305	ENDIF
2306	rrau(ji,jj,jk) = ralpbet * ( ptem(ji,jj,jk-1) - ptem(ji,jj,jk) ) / zds
2307	ptem_ad(ji,jj,jk-1) = ptem_ad(ji,jj,jk-1) &
2308	& + rrau_ad(ji,jj,jk) * ralpbet / zds
2309	ptem_ad(ji,jj,jk ) = ptem_ad(ji,jj,jk ) &
2310	& - rrau_ad(ji,jj,jk) * ralpbet / zds
2311	rrau_ad(ji,jj,jk) = 0._wp
2312
2313	psal_ad(ji,jj,jk-1) = psal_ad(ji,jj,jk-1) + zdsad
2314	psal_ad(ji,jj,jk ) = psal_ad(ji,jj,jk ) - zdsad
2315	zdsad = 0._wp
2316	END DO
2317	END DO
2318	END DO
2319	#endif
2320	DO jk = jpkm1, 2, -1
2321	ptem_ad(:,:,jk-1) = ptem_ad(:,:,jk-1) + rn_alpha * pn2_ad(:,:,jk) &
2322	& * grav / fse3w(:,:,jk) * tmask(:,:,jk)
2323	ptem_ad(:,:,jk ) = ptem_ad(:,:,jk ) - rn_alpha * pn2_ad(:,:,jk) &
2324	& * grav / fse3w(:,:,jk) * tmask(:,:,jk)
2325	psal_ad(:,:,jk-1) = psal_ad(:,:,jk-1) - rn_beta * pn2_ad(:,:,jk) &
2326	& * grav / fse3w(:,:,jk) * tmask(:,:,jk)
2327	psal_ad(:,:,jk ) = psal_ad(:,:,jk ) + rn_beta * pn2_ad(:,:,jk) &
2328	& * grav / fse3w(:,:,jk) * tmask(:,:,jk)
2329	pn2_ad(:,:,jk) = 0.0_wp
2330	END DO
2331	END SELECT
2332
2333	END SUBROUTINE eos_bn2_adj
2334
2335	#if defined key_tam
2336	SUBROUTINE eos_insitu_adj_tst( kumadt )
2337	!!-----------------------------------------------------------------------
2338	!!
2339	!! * ROUTINE eos_adj_tst *
2340	!!
2341	!! ** Purpose : Test the adjoint routine.
2342	!!
2343	!! ** Method : Verify the scalar product
2344	!!
2345	!! ( L dx )^T W dy = dx^T L^T W dy
2346	!!
2347	!! where L = tangent routine
2348	!! L^T = adjoint routine
2349	!! W = diagonal matrix of scale factors
2350	!! dx = input perturbation (random field)
2351	!! dy = L dx
2352	!!
2353	!!
2354	!! History :
2355	!! ! 08-07 (A. Vidard)
2356	!!-----------------------------------------------------------------------
2357	!! * Modules used
2358
2359	!! * Arguments
2360	INTEGER, INTENT(IN) :: &
2361	& kumadt ! Output unit
2362	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2363	ztem, & ! potential temperature
2364	zsal ! salinity
2365	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2366	& zt_adout, & ! potential temperature
2367	& zs_adout ! salinity
2368	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2369	& zrd_adin ! anomaly density
2370	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2371	& zt_tlin, & ! potential temperature
2372	& zs_tlin ! salinity
2373	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2374	& znt, & ! potential temperature
2375	& zns ! salinity
2376	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2377	& zrd_tlout ! anomaly density
2378	REAL(KIND=wp) :: &
2379	& zsp1, & ! scalar product involving the tangent routine
2380	& zsp2, & ! scalar product involving the adjoint routine
2381	& zsp2_1, & ! scalar product involving the adjoint routine
2382	& zsp2_2 ! scalar product involving the adjoint routine
2383	INTEGER, DIMENSION(jpi,jpj) :: &
2384	& iseed_2d ! 2D seed for the random number generator
2385	INTEGER :: &
2386	& ji, &
2387	& jj, &
2388	& jk, &
2389	& jn, &
2390	& jeos
2391
2392	CHARACTER(LEN=14) :: cl_name
2393
2394	ALLOCATE( &
2395	& ztem( jpi, jpj, jpk ), &
2396	& zsal( jpi, jpj, jpk ), &
2397	& znt( jpi, jpj, jpk ), &
2398	& zns( jpi, jpj, jpk ), &
2399	& zt_adout( jpi, jpj, jpk ), &
2400	& zs_adout( jpi, jpj, jpk ), &
2401	& zrd_adin( jpi, jpj, jpk ), &
2402	& zs_tlin( jpi, jpj, jpk ), &
2403	& zt_tlin( jpi, jpj, jpk ), &
2404	& zrd_tlout(jpi, jpj, jpk ) )
2405	! Initialize the reference state
2406	ztem(:,:,:) = tn(:,:,:)
2407	zsal(:,:,:) = sn(:,:,:)
2408
2409	! store initial nn_eos
2410	jeos = nn_eos
2411	DO jn = 0, 2
2412	nn_eos = jn
2413	!=============================================================
2414	! 1) dx = ( T ) and dy = ( T )
2415	!=============================================================
2416
2417	!--------------------------------------------------------------------
2418	! Reset the tangent and adjoint variables
2419	!--------------------------------------------------------------------
2420	zt_tlin(:,:,:) = 0.0_wp
2421	zs_tlin(:,:,:) = 0.0_wp
2422	zrd_tlout(:,:,:) = 0.0_wp
2423	zt_adout(:,:,:) = 0.0_wp
2424	zs_adout(:,:,:) = 0.0_wp
2425	zrd_adin(:,:,:) = 0.0_wp
2426
2427	!--------------------------------------------------------------------
2428	! Initialize the tangent input with random noise: dx
2429	!--------------------------------------------------------------------
2430	DO jj = 1, jpj
2431	DO ji = 1, jpi
2432	iseed_2d(ji,jj) = - ( 456953 + &
2433	& mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
2434	END DO
2435	END DO
2436	CALL grid_random( iseed_2d, znt, 'T', 0.0_wp, stdt )
2437	DO jj = 1, jpj
2438	DO ji = 1, jpi
2439	iseed_2d(ji,jj) = - ( 395703 + &
2440	& mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
2441	END DO
2442	END DO
2443	CALL grid_random( iseed_2d, zns, 'T', 0.0_wp, stds )
2444
2445	DO jk = 1, jpk
2446	DO jj = nldj, nlej
2447	DO ji = nldi, nlei
2448	zt_tlin(ji,jj,jk) = znt(ji,jj,jk)
2449	zs_tlin(ji,jj,jk) = zns(ji,jj,jk)
2450	END DO
2451	END DO
2452	END DO
2453
2454	CALL eos_insitu_tan(ztem, zsal, zt_tlin, zs_tlin, zrd_tlout)
2455
2456	DO jk = 1, jpk
2457	DO jj = nldj, nlej
2458	DO ji = nldi, nlei
2459	zrd_adin(ji,jj,jk) = zrd_tlout(ji,jj,jk) &
2460	& * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)&
2461	& * tmask(ji,jj,jk)
2462	END DO
2463	END DO
2464	END DO
2465
2466	!--------------------------------------------------------------------
2467	! Compute the scalar product: ( L dx )^T W dy
2468	!--------------------------------------------------------------------
2469
2470	zsp1 = DOT_PRODUCT( zrd_tlout, zrd_adin )
2471
2472	!--------------------------------------------------------------------
2473	! Call the adjoint routine: dx^* = L^T dy^*
2474	!--------------------------------------------------------------------
2475
2476	CALL eos_insitu_adj(ztem, zsal, zt_adout, zs_adout, zrd_adin)
2477
2478	zsp2_1 = DOT_PRODUCT( zt_tlin, zt_adout )
2479	zsp2_2 = DOT_PRODUCT( zs_tlin, zs_adout )
2480	zsp2 = zsp2_1 + zsp2_2
2481
2482	! Compare the scalar products
2483
2484	! Compare the scalar products
2485	! 14 char:'12345678901234'
2486	SELECT CASE( jn )
2487	CASE (0) ; cl_name = 'eos_adj ins T1'
2488	CASE (1) ; cl_name = 'eos_adj ins T2'
2489	CASE (2) ; cl_name = 'eos_adj ins T3'
2490	END SELECT
2491	CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
2492	ENDDO
2493	! restore initial nn_eos
2494	nn_eos = jeos
2495
2496	! Deallocate memory
2497
2498	DEALLOCATE( &
2499	& ztem, &
2500	& zsal, &
2501	& zt_adout, &
2502	& zs_adout, &
2503	& zrd_adin, &
2504	& zt_tlin, &
2505	& zs_tlin, &
2506	& zrd_tlout, &
2507	& znt, &
2508	& zns &
2509	& )
2510
2511
2512	END SUBROUTINE eos_insitu_adj_tst
2513	SUBROUTINE eos_insitu_pot_adj_tst( kumadt )
2514	!!-----------------------------------------------------------------------
2515	!!
2516	!! * ROUTINE eos_adj_tst *
2517	!!
2518	!! ** Purpose : Test the adjoint routine.
2519	!!
2520	!! ** Method : Verify the scalar product
2521	!!
2522	!! ( L dx )^T W dy = dx^T L^T W dy
2523	!!
2524	!! where L = tangent routine
2525	!! L^T = adjoint routine
2526	!! W = diagonal matrix of scale factors
2527	!! dx = input perturbation (random field)
2528	!! dy = L dx
2529	!!
2530	!! ** Action : Separate tests are applied for the following dx and dy:
2531	!!
2532	!! 1) dx = ( SSH ) and dy = ( SSH )
2533	!!
2534	!! History :
2535	!! ! 08-07 (A. Vidard)
2536	!!-----------------------------------------------------------------------
2537	!! * Modules used
2538
2539	!! * Arguments
2540	INTEGER, INTENT(IN) :: &
2541	& kumadt ! Output unit
2542	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2543	ztem, & ! potential temperature
2544	zsal ! salinity
2545	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2546	& zt_adout, & ! potential temperature
2547	& zs_adout ! salinity
2548	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2549	& zrd_adin ! anomaly density
2550	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2551	& zrhop_adin ! volume mass
2552	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2553	& zt_tlin, & ! potential temperature
2554	& zs_tlin ! salinity
2555	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2556	& znt, & ! potential temperature
2557	& zns ! salinity
2558	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2559	& zrd_tlout ! anomaly density
2560	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2561	& zrhop_tlout ! volume mass
2562	REAL(KIND=wp) :: &
2563	& zsp1, & ! scalar product involving the tangent routine
2564	& zsp2, & ! scalar product involving the adjoint routine
2565	& zsp2_1, & ! scalar product involving the adjoint routine
2566	& zsp2_2 ! scalar product involving the adjoint routine
2567	INTEGER, DIMENSION(jpi,jpj) :: &
2568	& iseed_2d ! 2D seed for the random number generator
2569	INTEGER :: &
2570	& ji, &
2571	& jj, &
2572	& jk, &
2573	& jn, &
2574	& jeos
2575	CHARACTER(LEN=14) :: cl_name
2576
2577	! Allocate memory
2578	ALLOCATE( &
2579	& ztem( jpi, jpj, jpk ), &
2580	& zsal( jpi, jpj, jpk ), &
2581	& zt_adout( jpi, jpj, jpk ), &
2582	& zs_adout( jpi, jpj, jpk ), &
2583	& zrhop_adin( jpi, jpj, jpk ), &
2584	& zrd_adin( jpi, jpj, jpk ), &
2585	& zs_tlin( jpi, jpj, jpk ), &
2586	& zt_tlin( jpi, jpj, jpk ), &
2587	& zns( jpi, jpj, jpk ), &
2588	& znt( jpi, jpj, jpk ), &
2589	& zrd_tlout(jpi, jpj, jpk ), &
2590	& zrhop_tlout(jpi, jpj, jpk ) )
2591
2592	! Initialize random field standard deviationsthe reference state
2593	ztem = tn
2594	zsal = sn
2595
2596	! store initial nn_eos
2597	jeos = nn_eos
2598	DO jn = 0, 2
2599	nn_eos = jn
2600	!=============================================
2601	! testing of eos_insitu_pot
2602	!=============================================
2603
2604	!=============================================================
2605	! 1) dx = ( T ) and dy = ( T )
2606	!=============================================================
2607
2608	!--------------------------------------------------------------------
2609	! Reset the tangent and adjoint variables
2610	!--------------------------------------------------------------------
2611	zt_tlin(:,:,:) = 0.0_wp
2612	zs_tlin(:,:,:) = 0.0_wp
2613	zrd_tlout(:,:,:) = 0.0_wp
2614	zrhop_tlout(:,:,:) = 0.0_wp
2615	zt_adout(:,:,:) = 0.0_wp
2616	zs_adout(:,:,:) = 0.0_wp
2617	zrhop_adin(:,:,:) = 0.0_wp
2618	zrd_adin(:,:,:) = 0.0_wp
2619
2620	!--------------------------------------------------------------------
2621	! Initialize the tangent input with random noise: dx
2622	!--------------------------------------------------------------------
2623	DO jj = 1, jpj
2624	DO ji = 1, jpi
2625	iseed_2d(ji,jj) = - ( 456953 + &
2626	& mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
2627	END DO
2628	END DO
2629	CALL grid_random( iseed_2d, znt, 'T', 0.0_wp, stdt )
2630	DO jj = 1, jpj
2631	DO ji = 1, jpi
2632	iseed_2d(ji,jj) = - ( 456953 + &
2633	& mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
2634	END DO
2635	END DO
2636	CALL grid_random( iseed_2d, zns, 'T', 0.0_wp, stds )
2637	DO jk = 1, jpk
2638	DO jj = nldj, nlej
2639	DO ji = nldi, nlei
2640	zt_tlin(ji,jj,jk) = znt(ji,jj,jk)
2641	zs_tlin(ji,jj,jk) = zns(ji,jj,jk)
2642	END DO
2643	END DO
2644	END DO
2645	!--------------------------------------------------------------------
2646	! Call the tangent routine: dy = L dx
2647	!--------------------------------------------------------------------
2648
2649	call eos_insitu_pot_tan ( ztem, zsal, zt_tlin, zs_tlin, zrd_tlout, zrhop_tlout )
2650
2651	!--------------------------------------------------------------------
2652	! Initialize the adjoint variables: dy^* = W dy
2653	!--------------------------------------------------------------------
2654	DO jk = 1, jpk
2655	DO jj = nldj, nlej
2656	DO ji = nldi, nlei
2657	zrd_adin(ji,jj,jk) = zrd_tlout(ji,jj,jk) &
2658	& * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)&
2659	& * tmask(ji,jj,jk)
2660	zrhop_adin(ji,jj,jk) = zrhop_tlout(ji,jj,jk) &
2661	& * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)&
2662	& * tmask(ji,jj,jk)
2663	END DO
2664	END DO
2665	END DO
2666
2667	!--------------------------------------------------------------------
2668	! Compute the scalar product: ( L dx )^T W dy
2669	!--------------------------------------------------------------------
2670
2671	zsp1 = DOT_PRODUCT( zrd_tlout , zrd_adin ) &
2672	& + DOT_PRODUCT( zrhop_tlout, zrhop_adin )
2673	!--------------------------------------------------------------------
2674	! Call the adjoint routine: dx^* = L^T dy^*
2675	!--------------------------------------------------------------------
2676
2677	CALL eos_insitu_pot_adj( ztem, zsal, zt_adout, zs_adout, zrd_adin, zrhop_adin )
2678	!--------------------------------------------------------------------
2679	! Compute the scalar product: dx^T L^T W dy
2680	!--------------------------------------------------------------------
2681
2682	zsp2_1 = DOT_PRODUCT( zt_tlin, zt_adout )
2683	zsp2_2 = DOT_PRODUCT( zs_tlin, zs_adout )
2684	zsp2 = zsp2_1 + zsp2_2
2685	! Compare the scalar products
2686
2687	! 14 char:'12345678901234'
2688	SELECT CASE( jn )
2689	CASE (0) ; cl_name = 'eos_adj pot T1'
2690	CASE (1) ; cl_name = 'eos_adj pot T2'
2691	CASE (2) ; cl_name = 'eos_adj pot T3'
2692	END SELECT
2693	CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
2694
2695	ENDDO
2696
2697	! restore initial nn_eos
2698	nn_eos = jeos
2699
2700	! Deallocate memory
2701	DEALLOCATE( &
2702	& ztem, &
2703	& zsal, &
2704	& zt_adout, &
2705	& zs_adout, &
2706	& zrd_adin, &
2707	& zrhop_adin, &
2708	& zt_tlin, &
2709	& zs_tlin, &
2710	& zrd_tlout, &
2711	& zrhop_tlout,&
2712	& zns, znt )
2713
2714	END SUBROUTINE eos_insitu_pot_adj_tst
2715	SUBROUTINE eos_insitu_2d_adj_tst( kumadt )
2716	!!-----------------------------------------------------------------------
2717	!!
2718	!! * ROUTINE eos_adj_tst *
2719	!!
2720	!! ** Purpose : Test the adjoint routine.
2721	!!
2722	!! ** Method : Verify the scalar product
2723	!!
2724	!! ( L dx )^T W dy = dx^T L^T W dy
2725	!!
2726	!! where L = tangent routine
2727	!! L^T = adjoint routine
2728	!! W = diagonal matrix of scale factors
2729	!! dx = input perturbation (random field)
2730	!! dy = L dx
2731	!!
2732	!! ** Action : Separate tests are applied for the following dx and dy:
2733	!!
2734	!! 1) dx = ( SSH ) and dy = ( SSH )
2735	!!
2736	!! History :
2737	!! ! 08-07 (A. Vidard)
2738	!!-----------------------------------------------------------------------
2739	!! * Modules used
2740
2741	!! * Arguments
2742	INTEGER, INTENT(IN) :: &
2743	& kumadt ! Output unit
2744	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
2745	zdep ! depth
2746	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
2747	ztem, & ! potential temperature
2748	zsal ! salinity
2749	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
2750	& zt_adout, & ! potential temperature
2751	& zs_adout ! salinity
2752	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
2753	& zrd_adin ! anomaly density
2754	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
2755	& zt_tlin, & ! potential temperature
2756	& zs_tlin ! salinity
2757	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
2758	& znt, & ! potential temperature
2759	& zns ! salinity
2760	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
2761	& zrd_tlout ! anomaly density
2762	REAL(KIND=wp) :: &
2763	& zsp1, & ! scalar product involving the tangent routine
2764	& zsp2, & ! scalar product involving the adjoint routine
2765	& zsp2_1, & ! scalar product involving the adjoint routine
2766	& zsp2_2 ! scalar product involving the adjoint routine
2767	INTEGER, DIMENSION(jpi,jpj) :: &
2768	& iseed_2d ! 2D seed for the random number generator
2769	INTEGER :: &
2770	& ji, &
2771	& jj, &
2772	& jn, &
2773	& jeos
2774	CHARACTER(LEN=14) :: cl_name
2775
2776	! Allocate memory
2777
2778	ALLOCATE( &
2779	& zdep( jpi, jpj ), &
2780	& ztem( jpi, jpj ), &
2781	& zsal( jpi, jpj ), &
2782	& znt( jpi, jpj ), &
2783	& zns( jpi, jpj ), &
2784	& zt_adout( jpi, jpj ), &
2785	& zs_adout( jpi, jpj ), &
2786	& zrd_adin( jpi, jpj ), &
2787	& zs_tlin( jpi, jpj ), &
2788	& zt_tlin( jpi, jpj ), &
2789	& zrd_tlout(jpi, jpj ) )
2790
2791	! Initialize the reference state
2792	ztem(:,:) = tn(:,:,2)
2793	zsal(:,:) = sn(:,:,2)
2794	zdep(:,:) = fsdept(:,:,2)
2795
2796	! store initial nn_eos
2797	jeos = nn_eos
2798	DO jn = 0, 2
2799	nn_eos = jn
2800	!=============================================================
2801	! 1) dx = ( T ) and dy = ( T )
2802	!=============================================================
2803
2804	!--------------------------------------------------------------------
2805	! Reset the tangent and adjoint variables
2806	!--------------------------------------------------------------------
2807	zt_tlin(:,:) = 0.0_wp
2808	zs_tlin(:,:) = 0.0_wp
2809	zrd_tlout(:,:) = 0.0_wp
2810	zt_adout(:,:) = 0.0_wp
2811	zs_adout(:,:) = 0.0_wp
2812	zrd_adin(:,:) = 0.0_wp
2813
2814	!--------------------------------------------------------------------
2815	! Initialize the tangent input with random noise: dx
2816	!--------------------------------------------------------------------
2817	DO jj = 1, jpj
2818	DO ji = 1, jpi
2819	iseed_2d(ji,jj) = - ( 456953 + &
2820	& mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
2821	END DO
2822	END DO
2823	CALL grid_random( iseed_2d, znt, 'T', 0.0_wp, stdt )
2824	DO jj = 1, jpj
2825	DO ji = 1, jpi
2826	iseed_2d(ji,jj) = - ( 456953 + &
2827	& mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
2828	END DO
2829	END DO
2830	CALL grid_random( iseed_2d, zns, 'T', 0.0_wp, stds )
2831	DO jj = nldj, nlej
2832	DO ji = nldi, nlei
2833	zt_tlin(ji,jj) = znt(ji,jj)
2834	zs_tlin(ji,jj) = zns(ji,jj)
2835	END DO
2836	END DO
2837
2838	CALL eos_insitu_2d_tan(ztem, zsal, zdep, zt_tlin, zs_tlin, zrd_tlout)
2839
2840	DO jj = nldj, nlej
2841	DO ji = nldi, nlei
2842	zrd_adin(ji,jj) = zrd_tlout(ji,jj) &
2843	& * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,2)&
2844	& * tmask(ji,jj,2)
2845	END DO
2846	END DO
2847
2848	!--------------------------------------------------------------------
2849	! Compute the scalar product: ( L dx )^T W dy
2850	!--------------------------------------------------------------------
2851
2852	zsp1 = DOT_PRODUCT( zrd_tlout, zrd_adin )
2853
2854	!--------------------------------------------------------------------
2855	! Call the adjoint routine: dx^* = L^T dy^*
2856	!--------------------------------------------------------------------
2857
2858	CALL eos_insitu_2d_adj(ztem, zsal, zdep, zt_adout, zs_adout, zrd_adin)
2859
2860	zsp2_1 = DOT_PRODUCT( zt_tlin, zt_adout )
2861	zsp2_2 = DOT_PRODUCT( zs_tlin, zs_adout )
2862	zsp2 = zsp2_1 + zsp2_2
2863
2864	! Compare the scalar products
2865
2866	! 14 char:'12345678901234'
2867	SELECT CASE( jn )
2868	CASE (0) ; cl_name = 'eos_adj 2d T1'
2869	CASE (1) ; cl_name = 'eos_adj 2d T2'
2870	CASE (2) ; cl_name = 'eos_adj 2d T3'
2871	END SELECT
2872	CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
2873
2874	ENDDO
2875
2876	! restore initial nn_eos
2877	nn_eos = jeos
2878
2879	! Deallocate memory
2880
2881	DEALLOCATE( &
2882	& zdep, &
2883	& ztem, &
2884	& zsal, &
2885	& zt_adout, &
2886	& zs_adout, &
2887	& zrd_adin, &
2888	& zt_tlin, &
2889	& zs_tlin, &
2890	& zrd_tlout, &
2891	& zns, znt )
2892
2893
2894	END SUBROUTINE eos_insitu_2d_adj_tst
2895	SUBROUTINE eos_pot_1pt_adj_tst( kumadt )
2896	!!-----------------------------------------------------------------------
2897	!!
2898	!! * ROUTINE eos_adj_tst *
2899	!!
2900	!! ** Purpose : Test the adjoint routine.
2901	!!
2902	!! ** Method : Verify the scalar product
2903	!!
2904	!! ( L dx )^T W dy = dx^T L^T W dy
2905	!!
2906	!! where L = tangent routine
2907	!! L^T = adjoint routine
2908	!! W = diagonal matrix of scale factors
2909	!! dx = input perturbation (random field)
2910	!! dy = L dx
2911	!!
2912	!! ** Action : Separate tests are applied for the following dx and dy:
2913	!!
2914	!! 1) dx = ( SSH ) and dy = ( SSH )
2915	!!
2916	!! History :
2917	!! ! 08-07 (A. Vidard)
2918	!!-----------------------------------------------------------------------
2919	!! * Modules used
2920
2921	!! * Arguments
2922	INTEGER, INTENT(IN) :: &
2923	& kumadt ! Output unit
2924
2925	!! * Local declarations
2926	REAL(KIND=wp) :: &
2927	& ztem, &
2928	& zsal, &
2929	& zt_tlin, &
2930	& zs_tlin, &
2931	& zrhop_tlout, &
2932	& zt_adout, &
2933	& zs_adout, &
2934	& zrhop_adin
2935	REAL(KIND=wp) :: &
2936	& zsp1, & ! scalar product involving the tangent routine
2937	& zsp2 ! scalar product involving the adjoint routine
2938	CHARACTER(LEN=14) :: cl_name
2939
2940	INTEGER :: &
2941	& jn, &
2942	& jeos
2943
2944	! Initialize the reference state
2945	ztem = 23.7
2946	zsal = 30.1
2947	! store initial nn_eos
2948	jeos = nn_eos
2949	DO jn = 0, 2
2950	nn_eos = jn
2951	!==================================================================
2952	! 1) dx = ( un_tl, vn_tl, hdivn_tl ) and
2953	! dy = ( hdivb_tl, hdivn_tl )
2954	!==================================================================
2955
2956	!--------------------------------------------------------------------
2957	! Reset the tangent and adjoint variables
2958	!--------------------------------------------------------------------
2959
2960	zt_tlin = 1.12_wp
2961	zs_tlin = 0.123_wp
2962	zrhop_tlout = 0.0_wp
2963	zt_adout = 0.0_wp
2964	zs_adout = 0.0_wp
2965	zrhop_adin = 0.0_wp
2966
2967	CALL eos_pot_1pt_tan( ztem, zsal, zt_tlin, zs_tlin, zrhop_tlout )
2968
2969	!--------------------------------------------------------------------
2970	! Initialize the adjoint variables: dy^* = W dy
2971	!--------------------------------------------------------------------
2972
2973	zrhop_adin = zrhop_tlout * e1t(1,1) * e2t(1,1) * fse3t(1,1,1)
2974
2975	!--------------------------------------------------------------------
2976	! Compute the scalar product: ( L dx )^T W dy
2977	!--------------------------------------------------------------------
2978
2979	zsp1 = zrhop_adin * zrhop_tlout
2980
2981	CALL eos_pot_1pt_adj( ztem, zsal, zt_adout, zs_adout, zrhop_adin )
2982
2983	zsp2 = zt_tlin * zt_adout + zs_tlin * zs_adout
2984
2985	! 14 char:'12345678901234'
2986	SELECT CASE( jn )
2987	CASE (0) ; cl_name = 'eos_adj 1pt T1'
2988	CASE (1) ; cl_name = 'eos_adj 1pt T2'
2989	CASE (2) ; cl_name = 'eos_adj 1pt T3'
2990	END SELECT
2991	CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
2992	ENDDO
2993	! restore initial nn_eos
2994	nn_eos = jeos
2995
2996	END SUBROUTINE eos_pot_1pt_adj_tst
2997	SUBROUTINE eos_adj_tst( kumadt )
2998	!!-----------------------------------------------------------------------
2999	!!
3000	!! * ROUTINE eos_adj_tst *
3001	!!
3002	!! ** Purpose : Test the adjoint routine.
3003	!!
3004	!! History :
3005	!! ! 08-07 (A. Vidard)
3006	!!-----------------------------------------------------------------------
3007	!! * Arguments
3008	INTEGER, INTENT(IN) :: &
3009	& kumadt ! Output unit
3010
3011	CALL eos_insitu_adj_tst( kumadt )
3012
3013	CALL eos_insitu_pot_adj_tst( kumadt )
3014
3015	CALL eos_insitu_2d_adj_tst( kumadt )
3016
3017	CALL eos_pot_1pt_adj_tst( kumadt )
3018
3019	END SUBROUTINE eos_adj_tst
3020	SUBROUTINE bn2_adj_tst( kumadt )
3021	!!-----------------------------------------------------------------------
3022	!!
3023	!! * ROUTINE bn2_adj_tst *
3024	!!
3025	!! ** Purpose : Test the adjoint routine.
3026	!!
3027	!! ** Method : Verify the scalar product
3028	!!
3029	!! ( L dx )^T W dy = dx^T L^T W dy
3030	!!
3031	!! where L = tangent routine
3032	!! L^T = adjoint routine
3033	!! W = diagonal matrix of scale factors
3034	!! dx = input perturbation (random field)
3035	!! dy = L dx
3036	!!
3037	!! ** Action : Separate tests are applied for the following dx and dy:
3038	!!
3039	!! 1) dx = ( SSH ) and dy = ( SSH )
3040	!!
3041	!! History :
3042	!! ! 08-07 (A. Vidard)
3043	!!-----------------------------------------------------------------------
3044	!! * Modules used
3045
3046	!! * Arguments
3047	INTEGER, INTENT(IN) :: &
3048	& kumadt ! Output unit
3049	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
3050	ztem, & ! potential temperature
3051	zsal ! salinity
3052	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
3053	& zt_adout, & ! potential temperature
3054	& zs_adout ! salinity
3055	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
3056	& zrd_adin, & ! potential density (surface referenced)
3057	& zrd_adout ! potential density (surface referenced)
3058	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
3059	& zt_tlin, & ! potential temperature
3060	& zs_tlin, & ! salinity
3061	& zt_tlout, & ! potential temperature
3062	& zs_tlout ! salinity
3063	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
3064	& zrd_tlout ! potential density (surface referenced)
3065	REAL(KIND=wp) :: &
3066	& zsp1, & ! scalar product involving the tangent routine
3067	& zsp2, & ! scalar product involving the adjoint routine
3068	& zsp2_1, & ! scalar product involving the adjoint routine
3069	& zsp2_2 ! scalar product involving the adjoint routine
3070	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
3071	& znt, & ! potential temperature
3072	& zns ! salinity
3073	INTEGER, DIMENSION(jpi,jpj) :: &
3074	& iseed_2d ! 2D seed for the random number generator
3075	INTEGER :: &
3076	& iseed, &
3077	& ji, &
3078	& jj, &
3079	& jk, &
3080	& jn, &
3081	& jeos
3082	CHARACTER(LEN=14) :: cl_name
3083
3084	! Allocate memory
3085	ALLOCATE( &
3086	& ztem( jpi, jpj, jpk ), &
3087	& zsal( jpi, jpj, jpk ), &
3088	& zt_adout( jpi, jpj, jpk ), &
3089	& zs_adout( jpi, jpj, jpk ), &
3090	& zrd_adin( jpi, jpj, jpk ), &
3091	& zrd_adout(jpi, jpj, jpk ), &
3092	& zs_tlin( jpi, jpj, jpk ), &
3093	& zt_tlin( jpi, jpj, jpk ), &
3094	& zns( jpi, jpj, jpk ), &
3095	& znt( jpi, jpj, jpk ), &
3096	& zt_tlout( jpi, jpj, jpk ), &
3097	& zs_tlout( jpi, jpj, jpk ), &
3098	& zrd_tlout(jpi, jpj, jpk ) )
3099
3100	! Initialize random field standard deviationsthe reference state
3101	ztem = tn
3102	zsal = sn
3103	! store initial nn_eos
3104	jeos = nn_eos
3105	DO jn = 0, 2
3106	nn_eos = jn
3107	!=============================================================
3108	! 1) dx = ( T ) and dy = ( T )
3109	!=============================================================
3110
3111	!--------------------------------------------------------------------
3112	! Reset the tangent and adjoint variables
3113	!--------------------------------------------------------------------
3114	zt_tlin(:,:,:) = 0.0_wp
3115	zs_tlin(:,:,:) = 0.0_wp
3116	zt_tlout(:,:,:) = 0.0_wp
3117	zs_tlout(:,:,:) = 0.0_wp
3118	zrd_tlout(:,:,:) = 0.0_wp
3119	zt_adout(:,:,:) = 0.0_wp
3120	zs_adout(:,:,:) = 0.0_wp
3121	zrd_adin(:,:,:) = 0.0_wp
3122	zrd_adout(:,:,:) = 0.0_wp
3123
3124	!--------------------------------------------------------------------
3125	! Initialize the tangent input with random noise: dx
3126	!--------------------------------------------------------------------
3127	DO jj = 1, jpj
3128	DO ji = 1, jpi
3129	iseed_2d(ji,jj) = - ( 456953 + &
3130	& mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
3131	END DO
3132	END DO
3133	CALL grid_random( iseed_2d, znt, 'T', 0.0_wp, stdt )
3134	DO jj = 1, jpj
3135	DO ji = 1, jpi
3136	iseed_2d(ji,jj) = - ( 456953 + &
3137	& mig(ji) + ( mjg(jj) - 1 ) * jpiglo )
3138	END DO
3139	END DO
3140	CALL grid_random( iseed_2d, zns, 'T', 0.0_wp, stds )
3141	DO jk = 1, jpk
3142	DO jj = nldj, nlej
3143	DO ji = nldi, nlei
3144	zt_tlin(ji,jj,jk) = znt(ji,jj,jk)
3145	zs_tlin(ji,jj,jk) = zns(ji,jj,jk)
3146	END DO
3147	END DO
3148	END DO
3149	!--------------------------------------------------------------------
3150	! Call the tangent routine: dy = L dx
3151	!--------------------------------------------------------------------
3152	zt_tlout(:,:,:) = zt_tlin
3153	zs_tlout(:,:,:) = zs_tlin
3154
3155	CALL eos_bn2_tan( ztem, zsal, zt_tlout, zs_tlout, zrd_tlout )
3156	!--------------------------------------------------------------------
3157	! Initialize the adjoint variables: dy^* = W dy
3158	!--------------------------------------------------------------------
3159	DO jk = 1, jpk
3160	DO jj = nldj, nlej
3161	DO ji = nldi, nlei
3162	zrd_adin(ji,jj,jk) = zrd_tlout(ji,jj,jk) &
3163	& * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
3164	& * tmask(ji,jj,jk)
3165	END DO
3166	END DO
3167	END DO
3168
3169	!--------------------------------------------------------------------
3170	! Compute the scalar product: ( L dx )^T W dy
3171	!--------------------------------------------------------------------
3172
3173	zsp1 = DOT_PRODUCT( zrd_tlout, zrd_adin )
3174
3175	!--------------------------------------------------------------------
3176	! Call the adjoint routine: dx^* = L^T dy^*
3177	!--------------------------------------------------------------------
3178
3179	zrd_adout(:,:,:) = zrd_adin(:,:,:)
3180
3181	CALL eos_bn2_adj( ztem, zsal, zt_adout, zs_adout, zrd_adout )
3182
3183	!--------------------------------------------------------------------
3184	! Compute the scalar product: dx^T L^T W dy
3185	!--------------------------------------------------------------------
3186
3187	zsp2_1 = DOT_PRODUCT( zt_tlin, zt_adout )
3188	zsp2_2 = DOT_PRODUCT( zs_tlin, zs_adout )
3189	zsp2 = zsp2_1 + zsp2_2
3190
3191	! Compare the scalar products
3192
3193	! 14 char:'12345678901234'
3194	SELECT CASE( jn )
3195	CASE (0) ; cl_name = 'bn2_adj T1'
3196	CASE (1) ; cl_name = 'bn2_adj T2'
3197	CASE (2) ; cl_name = 'bn2_adj T3'
3198	END SELECT
3199	CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
3200	ENDDO
3201	! restore initial nn_eos
3202	nn_eos = jeos
3203
3204	! Deallocate memory
3205
3206	DEALLOCATE( &
3207	& ztem, &
3208	& zsal, &
3209	& zt_adout, &
3210	& zs_adout, &
3211	& zrd_adin, &
3212	& zrd_adout, &
3213	& zt_tlin, &
3214	& zs_tlin, &
3215	& zt_tlout, &
3216	& zs_tlout, &
3217	& zrd_tlout, &
3218	& zns, znt )
3219
3220
3221	END SUBROUTINE bn2_adj_tst
3222	#if defined key_tst_tlm
3223	SUBROUTINE eos_insitu_tlm_tst( kumadt )
3224	!!-----------------------------------------------------------------------
3225	!!
3226	!! * ROUTINE eos_insitu_tlm_tst *
3227	!!
3228	!! ** Purpose : Test the tangent routine.
3229	!!
3230	!! ** Method : Verify the tangent with Taylor expansion
3231	!!
3232	!! M(x+hdx) = M(x) + L(hdx) + O(h^2)
3233	!!
3234	!! where L = tangent routine
3235	!! M = direct routine
3236	!! dx = input perturbation (random field)
3237	!! h = ration on perturbation
3238	!!
3239	!! In the tangent test we verify that:
3240	!! M(x+h*dx) - M(x)
3241	!! g(h) = ------------------ ---> 1 as h ---> 0
3242	!! L(h*dx)
3243	!! and
3244	!! g(h) - 1
3245	!! f(h) = ---------- ---> k (costant) as h ---> 0
3246	!! p
3247	!!
3248	!! History :
3249	!! ! 09-08 (A. Vigilant)
3250	!!-----------------------------------------------------------------------
3251	!! * Modules used
3252	USE eosbn2, ONLY: & ! horizontal & vertical advective trend
3253	& eos
3254	USE tamtrj ! writing out state trajectory
3255	USE par_tlm, ONLY: &
3256	& tlm_bch, &
3257	& cur_loop, &
3258	& h_ratio
3259	USE istate_mod
3260	USE gridrandom, ONLY: &
3261	& grid_rd_sd
3262	USE trj_tam
3263	USE oce , ONLY: & ! ocean dynamics and tracers variables
3264	& tn, sn, rhd, rhop
3265	USE in_out_manager, ONLY: & ! I/O manager
3266	& nitend, &
3267	& nit000
3268	USE tamctl, ONLY: & ! Control parameters
3269	& numtan, numtan_sc
3270	!! * Arguments
3271	INTEGER, INTENT(IN) :: &
3272	& kumadt ! Output unit
3273	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
3274	& zrd_out, & !
3275	& zt_tlin , & !
3276	& zs_tlin , &
3277	& zrd_tl , &
3278	& zrd_wop, &
3279	& z3r
3280	REAL(KIND=wp) :: &
3281	& zsp1, &
3282	& zsp2, &
3283	& zsp3, &
3284	& zzsp, &
3285	& gamma, &
3286	& zgsp1, &
3287	& zgsp2, &
3288	& zgsp3, &
3289	& zgsp4, &
3290	& zgsp5, &
3291	& zgsp6, &
3292	& zgsp7
3293	INTEGER :: &
3294	& ji, &
3295	& jj, &
3296	& jk
3297	CHARACTER(LEN=14) :: cl_name
3298	CHARACTER (LEN=128) :: file_out_sc, file_wop, file_out, file_xdx
3299	CHARACTER (LEN=90) :: FMT
3300	REAL(KIND=wp), DIMENSION(100):: &
3301	& zscrd, zscerrrd
3302	INTEGER, DIMENSION(100):: &
3303	& iiposrd, ijposrd, ikposrd
3304	INTEGER:: &
3305	& ii, numsctlm, &
3306	& numtlm, &
3307	& isamp=40,jsamp=40, ksamp=10
3308	REAL(KIND=wp), DIMENSION(jpi,jpj,jpk) :: &
3309	& zerrrd
3310	ALLOCATE( &
3311	& zrd_out( jpi, jpj, jpk ), &
3312	& zrd_tl( jpi, jpj, jpk ), &
3313	& zs_tlin( jpi, jpj, jpk ), &
3314	& zt_tlin( jpi, jpj, jpk ), &
3315	& zrd_wop( jpi, jpj, jpk ), &
3316	& z3r ( jpi, jpj, jpk ) )
3317
3318	!--------------------------------------------------------------------
3319	! Reset the tangent and adjoint variables
3320	!--------------------------------------------------------------------
3321	zt_tlin( :,:,:) = 0.0_wp
3322	zs_tlin( :,:,:) = 0.0_wp
3323	zrd_out( :,:,:) = 0.0_wp
3324	zrd_wop( :,:,:) = 0.0_wp
3325	zscerrrd(:) = 0.0_wp
3326	zscrd(:) = 0.0_wp
3327	IF ( tlm_bch == 2 ) zrd_tl ( :,:,:) = 0.0_wp
3328	!--------------------------------------------------------------------
3329	! Output filename Xn=F(X0)
3330	!--------------------------------------------------------------------
3331	! CALL tlm_namrd
3332	gamma = h_ratio
3333	file_wop='trj_wop_eos_insitu'
3334	file_xdx='trj_xdx_eos_insitu'
3335	!--------------------------------------------------------------------
3336	! Initialize the tangent input with random noise: dx
3337	!--------------------------------------------------------------------
3338	IF (( cur_loop .NE. 0) .OR. ( gamma .NE. 0.0_wp) )THEN
3339	CALL grid_rd_sd( 596035, z3r, 'T', 0.0_wp, stdt)
3340	DO jk = 1, jpk
3341	DO jj = nldj, nlej
3342	DO ji = nldi, nlei
3343	zt_tlin(ji,jj,jk) = z3r(ji,jj,jk)
3344	END DO
3345	END DO
3346	END DO
3347	CALL grid_rd_sd( 371836, z3r, 'S', 0.0_wp, stds)
3348	DO jk = 1, jpk
3349	DO jj = nldj, nlej
3350	DO ji = nldi, nlei
3351	zs_tlin(ji,jj,jk) = z3r(ji,jj,jk)
3352	END DO
3353	END DO
3354	END DO
3355	ENDIF
3356
3357	!--------------------------------------------------------------------
3358	! Complete Init for Direct
3359	!-------------------------------------------------------------------
3360	IF ( tlm_bch /= 2 ) CALL istate_p
3361
3362	! *** initialize the reference trajectory
3363	! ------------
3364	CALL trj_rea( nit000-1, 1 )
3365	CALL trj_rea( nit000, 1 )
3366
3367	IF (( cur_loop .NE. 0) .OR. ( gamma .NE. 0.0_wp) )THEN
3368	zt_tlin(:,:,:) = gamma * zt_tlin(:,:,:)
3369	tn(:,:,:) = tn(:,:,:) + zt_tlin(:,:,:)
3370
3371	zs_tlin(:,:,:) = gamma * zs_tlin(:,:,:)
3372	sn(:,:,:) = sn(:,:,:) + zs_tlin(:,:,:)
3373	ENDIF
3374	!--------------------------------------------------------------------
3375	! Compute the direct model F(X0,t=n) = Xn
3376	!--------------------------------------------------------------------
3377	IF ( tlm_bch /= 2 ) CALL eos(tn, sn, zrd_out)
3378	rhd(:,:,:)= zrd_out(:,:,:)
3379	IF ( tlm_bch == 0 ) CALL trj_wri_spl(file_wop)
3380	IF ( tlm_bch == 1 ) CALL trj_wri_spl(file_xdx)
3381	!--------------------------------------------------------------------
3382	! Compute the Tangent
3383	!--------------------------------------------------------------------
3384	IF ( tlm_bch == 2 ) THEN
3385	!--------------------------------------------------------------------
3386	! Initialize the tangent variables: dy^* = W dy
3387	!--------------------------------------------------------------------
3388	CALL trj_rea( nit000-1, 1 )
3389	CALL trj_rea( nit000, 1 )
3390	!-----------------------------------------------------------------------
3391	! Initialization of the dynamics and tracer fields for the tangent
3392	!-----------------------------------------------------------------------
3393	CALL eos_insitu_tan(tn, sn, zt_tlin, zs_tlin, zrd_tl)
3394	!--------------------------------------------------------------------
3395	! Compute the scalar product: ( L(t0,tn) gamma dx0 ) )
3396	!--------------------------------------------------------------------
3397	zsp2 = DOT_PRODUCT( zrd_tl, zrd_tl )
3398	!--------------------------------------------------------------------
3399	! Storing data
3400	!--------------------------------------------------------------------
3401	CALL trj_rd_spl(file_wop)
3402	zrd_wop (:,:,:) = rhd (:,:,:)
3403	CALL trj_rd_spl(file_xdx)
3404	zrd_out (:,:,:) = rhd (:,:,:)
3405	!--------------------------------------------------------------------
3406	! Compute the Linearization Error
3407	! Nn = M( X0+gamma.dX0, t0,tn) - M(X0, t0,tn)
3408	! and
3409	! Compute the Linearization Error
3410	! En = Nn -TL(gamma.dX0, t0,tn)
3411	!--------------------------------------------------------------------
3412	! Warning: Here we re-use local variables z()_out and z()_wop
3413	ii=0
3414	DO jk = 1, jpk
3415	DO jj = 1, jpj
3416	DO ji = 1, jpi
3417	zrd_out (ji,jj,jk) = zrd_out (ji,jj,jk) - zrd_wop (ji,jj,jk)
3418	zrd_wop (ji,jj,jk) = zrd_out (ji,jj,jk) - zrd_tl (ji,jj,jk)
3419	IF ( zrd_tl(ji,jj,jk) .NE. 0.0_wp ) &
3420	& zerrrd(ji,jj,jk) = zrd_out(ji,jj,jk)/zrd_tl(ji,jj,jk)
3421	IF( (MOD(ji, isamp) .EQ. 0) .AND. &
3422	& (MOD(jj, jsamp) .EQ. 0) .AND. &
3423	& (MOD(jk, ksamp) .EQ. 0) ) THEN
3424	ii = ii+1
3425	iiposrd(ii) = ji
3426	ijposrd(ii) = jj
3427	ikposrd(ii) = jk
3428	IF ( INT(tmask(ji,jj,jk)) .NE. 0) THEN
3429	zscrd (ii) = zrd_wop(ji,jj,jk)
3430	zscerrrd (ii) = ( zerrrd(ji,jj,jk) - 1.0_wp ) / gamma
3431	ENDIF
3432	ENDIF
3433	END DO
3434	END DO
3435	END DO
3436
3437	zsp1 = DOT_PRODUCT( zrd_out, zrd_out )
3438	zsp3 = DOT_PRODUCT( zrd_wop, zrd_wop )
3439
3440	!--------------------------------------------------------------------
3441	! Print the linearization error En - norme 2
3442	!--------------------------------------------------------------------
3443	! 14 char:'12345678901234'
3444	cl_name = 'eos_insitu:En '
3445	zzsp = SQRT(zsp3)
3446	zgsp5 = zzsp
3447	CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio )
3448	!--------------------------------------------------------------------
3449	! Compute TLM norm2
3450	!--------------------------------------------------------------------
3451	zzsp = SQRT(zsp2)
3452	zgsp4 = zzsp
3453	cl_name = 'eos_insitu:Ln2'
3454	CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio )
3455	!--------------------------------------------------------------------
3456	! Print the linearization error Nn - norme 2
3457	!--------------------------------------------------------------------
3458	zzsp = SQRT(zsp1)
3459	cl_name = 'eosins:Mhdx-Mx'
3460	CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio )
3461
3462	zgsp3 = SQRT( zsp3/zsp2 )
3463	zgsp7 = zgsp3/gamma
3464	zgsp1 = zzsp
3465	zgsp2 = zgsp1 / zgsp4
3466	zgsp6 = (zgsp2 - 1.0_wp)/gamma
3467
3468	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)"
3469	WRITE(numtan,FMT) 'eosinsit', cur_loop, h_ratio, zgsp1, zgsp2, zgsp3, zgsp4, zgsp5, zgsp6, zgsp7
3470
3471	!--------------------------------------------------------------------
3472	! Unitary calculus
3473	!--------------------------------------------------------------------
3474	FMT = "(A8,2X,A8,2X,I4.4,2X,E6.1,2X,I4.4,2X,I4.4,2X,I4.4,2X,E20.13,1X)"
3475	cl_name = 'eosinsit'
3476	IF(lwp) THEN
3477	DO ii=1, 100, 1
3478	IF ( zscrd(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscrd ', &
3479	& cur_loop, h_ratio, iiposrd(ii), ijposrd(ii), ikposrd(ii), zscrd(ii)
3480	ENDDO
3481	DO ii=1, 100, 1
3482	IF ( zscerrrd(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscerrrd', &
3483	& cur_loop, h_ratio, iiposrd(ii), ijposrd(ii), ikposrd(ii), zscerrrd(ii)
3484	ENDDO
3485	! write separator
3486	WRITE(numtan_sc,"(A4)") '===='
3487	ENDIF
3488
3489	ENDIF
3490
3491	DEALLOCATE( &
3492	& zrd_out, zrd_tl, zrd_wop, &
3493	& zt_tlin, zs_tlin, z3r &
3494	& )
3495	END SUBROUTINE eos_insitu_tlm_tst
3496
3497	SUBROUTINE eos_insitu_pot_tlm_tst( kumadt )
3498	!!-----------------------------------------------------------------------
3499	!!
3500	!! * ROUTINE eos_insitu_pot_tlm_tst *
3501	!!
3502	!! ** Purpose : Test the tangent routine.
3503	!!
3504	!! ** Method : Verify the tangent with Taylor expansion
3505	!!
3506	!! M(x+hdx) = M(x) + L(hdx) + O(h^2)
3507	!!
3508	!! where L = tangent routine
3509	!! M = direct routine
3510	!! dx = input perturbation (random field)
3511	!! h = ration on perturbation
3512	!!
3513	!! In the tangent test we verify that:
3514	!! M(x+h*dx) - M(x)
3515	!! g(h) = ------------------ ---> 1 as h ---> 0
3516	!! L(h*dx)
3517	!! and
3518	!! g(h) - 1
3519	!! f(h) = ---------- ---> k (costant) as h ---> 0
3520	!! p
3521	!!
3522	!! History :
3523	!! ! 09-08 (A. Vigilant)
3524	!!-----------------------------------------------------------------------
3525	!! * Modules used
3526	USE eosbn2, ONLY: & ! horizontal & vertical advective trend
3527	& eos
3528	USE tamtrj ! writing out state trajectory
3529	USE par_tlm, ONLY: &
3530	& tlm_bch, &
3531	& cur_loop, &
3532	& h_ratio
3533	USE istate_mod
3534	USE gridrandom, ONLY: &
3535	& grid_rd_sd
3536	USE trj_tam
3537	USE oce , ONLY: & ! ocean dynamics and tracers variables
3538	& tn, sn, rhd, rhop
3539	USE in_out_manager, ONLY: & ! I/O manager
3540	& nitend, &
3541	& nit000
3542	USE tamctl, ONLY: & ! Control parameters
3543	& numtan, numtan_sc
3544	!! * Arguments
3545	INTEGER, INTENT(IN) :: &
3546	& kumadt ! Output unit
3547	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
3548	& zrd_out, & !
3549	& zrh_out, &
3550	& zt_tlin , & !
3551	& zs_tlin , &
3552	& zrh_tl , &
3553	& zrd_tl , &
3554	& zrd_wop , &
3555	& zrh_wop , &
3556	& z3r
3557	REAL(KIND=wp) :: &
3558	& zsp1, &
3559	& zsp2, &
3560	& zsp3, &
3561	& zsp1_Rd, Zsp1_Rh, &
3562	& zsp2_Rd, zsp2_Rh, &
3563	& zsp3_Rd, zsp3_Rh, &
3564	& zzsp, &
3565	& gamma, &
3566	& zgsp1, &
3567	& zgsp2, &
3568	& zgsp3, &
3569	& zgsp4, &
3570	& zgsp5, &
3571	& zgsp6, &
3572	& zgsp7
3573	INTEGER :: &
3574	& ji, &
3575	& jj, &
3576	& jk
3577	CHARACTER(LEN=14) :: cl_name
3578	CHARACTER (LEN=128) :: file_out, file_wop,file_xdx
3579	CHARACTER (LEN=90) :: FMT
3580	REAL(KIND=wp), DIMENSION(100):: &
3581	& zscrd, zscerrrd, &
3582	& zscrh, zscerrrh
3583	INTEGER, DIMENSION(100):: &
3584	& iiposrd, ijposrd, ikposrd, &
3585	& iiposrh, ijposrh, ikposrh
3586	INTEGER:: &
3587	& ii, numsctlm, &
3588	& isamp=40,jsamp=40, ksamp=10
3589	REAL(KIND=wp), DIMENSION(jpi,jpj,jpk) :: &
3590	& zerrrd, zerrrh
3591	ALLOCATE( &
3592	& zrd_out( jpi, jpj, jpk ), &
3593	& zrh_out( jpi, jpj, jpk ), &
3594	& zrd_tl( jpi, jpj, jpk ), &
3595	& zrh_tl( jpi, jpj, jpk ), &
3596	& zs_tlin( jpi, jpj, jpk ), &
3597	& zt_tlin( jpi, jpj, jpk ), &
3598	& zrd_wop( jpi, jpj, jpk ), &
3599	& zrh_wop( jpi, jpj, jpk ), &
3600	& z3r ( jpi, jpj, jpk ) )
3601
3602	!--------------------------------------------------------------------
3603	! Reset the tangent and adjoint variables
3604	!--------------------------------------------------------------------
3605	zt_tlin( :,:,:) = 0.0_wp
3606	zs_tlin( :,:,:) = 0.0_wp
3607	zrd_out( :,:,:) = 0.0_wp
3608	zrh_out( :,:,:) = 0.0_wp
3609	zrd_wop( :,:,:) = 0.0_wp
3610	zrh_wop( :,:,:) = 0.0_wp
3611	zscerrrd( :) = 0.0_wp
3612	zscerrrh( :) = 0.0_wp
3613	zscrd(:) = 0.0_wp
3614	zscrh(:) = 0.0_wp
3615	IF ( tlm_bch == 2 ) THEN
3616	zrd_tl ( :,:,:) = 0.0_wp
3617	zrh_tl ( :,:,:) = 0.0_wp
3618	ENDIF
3619	!--------------------------------------------------------------------
3620	! Output filename Xn=F(X0)
3621	!--------------------------------------------------------------------
3622	!! CALL tlm_namrd
3623	gamma = h_ratio
3624	file_wop='trj_wop_eos_pot'
3625	file_xdx='trj_xdx_eos_pot'
3626	!--------------------------------------------------------------------
3627	! Initialize the tangent input with random noise: dx
3628	!--------------------------------------------------------------------
3629	IF (( cur_loop .NE. 0) .OR. ( gamma .NE. 0.0_wp) )THEN
3630	CALL grid_rd_sd( 596035, z3r, 'T', 0.0_wp, stdt)
3631	DO jk = 1, jpk
3632	DO jj = nldj, nlej
3633	DO ji = nldi, nlei
3634	zt_tlin(ji,jj,jk) = z3r(ji,jj,jk)
3635	END DO
3636	END DO
3637	END DO
3638	CALL grid_rd_sd( 371836, z3r, 'S', 0.0_wp, stds)
3639	DO jk = 1, jpk
3640	DO jj = nldj, nlej
3641	DO ji = nldi, nlei
3642	zs_tlin(ji,jj,jk) = z3r(ji,jj,jk)
3643	END DO
3644	END DO
3645	END DO
3646	ENDIF
3647	!--------------------------------------------------------------------
3648	! Complete Init for Direct
3649	!-------------------------------------------------------------------
3650	IF ( tlm_bch /= 2 ) CALL istate_p
3651
3652	! *** initialize the reference trajectory
3653	! ------------
3654	CALL trj_rea( nit000-1, 1 )
3655	CALL trj_rea( nit000, 1 )
3656
3657	IF (( cur_loop .NE. 0) .OR. ( gamma .NE. 0.0_wp) )THEN
3658	zt_tlin(:,:,:) = gamma * zt_tlin(:,:,:)
3659	tn(:,:,:) = tn(:,:,:) + zt_tlin(:,:,:)
3660
3661	zs_tlin(:,:,:) = gamma * zs_tlin(:,:,:)
3662	sn(:,:,:) = sn(:,:,:) + zs_tlin(:,:,:)
3663	ENDIF
3664	!--------------------------------------------------------------------
3665	! Compute the direct model F(X0,t=n) = Xn
3666	!--------------------------------------------------------------------
3667	IF ( tlm_bch /= 2 ) CALL eos(tn, sn, zrd_out, zrh_out)
3668	rhd (:,:,:) = zrd_out(:,:,:)
3669	rhop(:,:,:) = zrh_out(:,:,:)
3670	IF ( tlm_bch == 0 ) CALL trj_wri_spl(file_wop)
3671	IF ( tlm_bch == 1 ) CALL trj_wri_spl(file_xdx)
3672	!--------------------------------------------------------------------
3673	! Compute the Tangent
3674	!--------------------------------------------------------------------
3675	IF ( tlm_bch == 2 ) THEN
3676	!--------------------------------------------------------------------
3677	! Initialize the tangent variables: dy^* = W dy
3678	!--------------------------------------------------------------------
3679	CALL trj_rea( nit000-1, 1 )
3680	CALL trj_rea( nit000, 1 )
3681	!-----------------------------------------------------------------------
3682	! Initialization of the dynamics and tracer fields for the tangent
3683	!-----------------------------------------------------------------------
3684	CALL eos_insitu_pot_tan(tn, sn, zt_tlin, zs_tlin, zrd_tl, zrh_tl)
3685	!--------------------------------------------------------------------
3686	! Compute the scalar product: ( L(t0,tn) gamma dx0 ) )
3687	!--------------------------------------------------------------------
3688
3689	zsp2_Rd = DOT_PRODUCT( zrd_tl, zrd_tl )
3690	zsp2_Rh = DOT_PRODUCT( zrh_tl, zrh_tl )
3691	zsp2 = zsp2_Rd + zsp2_Rh
3692	!--------------------------------------------------------------------
3693	! Storing data
3694	!--------------------------------------------------------------------
3695	CALL trj_rd_spl(file_wop)
3696	zrd_wop (:,:,:) = rhd (:,:,:)
3697	zrh_wop (:,:,:) = rhop (:,:,:)
3698	CALL trj_rd_spl(file_xdx)
3699	zrd_out (:,:,:) = rhd (:,:,:)
3700	zrh_out (:,:,:) = rhop (:,:,:)
3701	!--------------------------------------------------------------------
3702	! Compute the Linearization Error
3703	! Nn = M( X0+gamma.dX0, t0,tn) - M(X0, t0,tn)
3704	! and
3705	! Compute the Linearization Error
3706	! En = Nn -TL(gamma.dX0, t0,tn)
3707	!--------------------------------------------------------------------
3708	! Warning: Here we re-use local variables z()_out and z()_wop
3709	ii=0
3710	DO jk = 1, jpk
3711	DO jj = 1, jpj
3712	DO ji = 1, jpi
3713	zrd_out (ji,jj,jk) = zrd_out (ji,jj,jk) - zrd_wop (ji,jj,jk)
3714	zrd_wop (ji,jj,jk) = zrd_out (ji,jj,jk) - zrd_tl (ji,jj,jk)
3715	IF ( zrd_tl(ji,jj,jk) .NE. 0.0_wp ) &
3716	& zerrrd(ji,jj,jk) = zrd_out(ji,jj,jk)/zrd_tl(ji,jj,jk)
3717	IF( (MOD(ji, isamp) .EQ. 0) .AND. &
3718	& (MOD(jj, jsamp) .EQ. 0) .AND. &
3719	& (MOD(jk, ksamp) .EQ. 0) ) THEN
3720	ii = ii+1
3721	iiposrd(ii) = ji
3722	ijposrd(ii) = jj
3723	ikposrd(ii) = jk
3724	IF ( INT(tmask(ji,jj,jk)) .NE. 0) THEN
3725	zscrd (ii) = zrd_wop(ji,jj,jk)
3726	zscerrrd (ii) = ( zerrrd( ji,jj,jk) - 1.0_wp ) / gamma
3727	ENDIF
3728	ENDIF
3729	END DO
3730	END DO
3731	END DO
3732	ii=0
3733	DO jk = 1, jpk
3734	DO jj = 1, jpj
3735	DO ji = 1, jpi
3736	zrh_out (ji,jj,jk) = zrh_out (ji,jj,jk) - zrh_wop (ji,jj,jk)
3737	zrh_wop (ji,jj,jk) = zrh_out (ji,jj,jk) - zrh_tl (ji,jj,jk)
3738	IF ( zrh_tl(ji,jj,jk) .NE. 0.0_wp ) &
3739	& zerrrh(ji,jj,jk) = zrh_out(ji,jj,jk)/zrh_tl(ji,jj,jk)
3740	IF( (MOD(ji, isamp) .EQ. 0) .AND. &
3741	& (MOD(jj, jsamp) .EQ. 0) .AND. &
3742	& (MOD(jk, ksamp) .EQ. 0) ) THEN
3743	ii = ii+1
3744	iiposrh(ii) = ji
3745	ijposrh(ii) = jj
3746	ikposrh(ii) = jk
3747	IF ( INT(tmask(ji,jj,jk)) .NE. 0) THEN
3748	zscrh (ii) = zrh_wop(ji,jj,jk)
3749	zscerrrh (ii) = ( zerrrh( ji,jj,jk) - 1.0_wp ) /gamma
3750	ENDIF
3751	ENDIF
3752	END DO
3753	END DO
3754	END DO
3755	zsp1_Rd = DOT_PRODUCT( zrd_out, zrd_out )
3756	zsp1_Rh = DOT_PRODUCT( zrh_out, zrh_out )
3757	zsp1 = zsp1_Rd + zsp1_Rh
3758
3759	zsp3_Rd = DOT_PRODUCT( zrd_wop, zrd_wop )
3760	zsp3_Rh = DOT_PRODUCT( zrh_wop, zrh_wop )
3761	zsp3 = zsp3_Rd + zsp3_Rh
3762	!--------------------------------------------------------------------
3763	! Print the linearization error En - norme 2
3764	!--------------------------------------------------------------------
3765	! 14 char:'12345678901234'
3766	cl_name = 'eos_pot:En '
3767	zzsp = SQRT(zsp3)
3768	zgsp5 = zzsp
3769	CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio )
3770	!--------------------------------------------------------------------
3771	! Compute TLM norm2
3772	!--------------------------------------------------------------------
3773	zzsp = SQRT(zsp2)
3774	zgsp4 = zzsp
3775	cl_name = 'eos_pot:Ln2'
3776	CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio )
3777	!--------------------------------------------------------------------
3778	! Print the linearization error Nn - norme 2
3779	!--------------------------------------------------------------------
3780	zzsp = SQRT(zsp1)
3781	cl_name = 'eospot:Mhdx-Mx'
3782	CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio )
3783
3784	zgsp3 = SQRT( zsp3/zsp2 )
3785	zgsp7 = zgsp3/gamma
3786	zgsp1 = zzsp
3787	zgsp2 = zgsp1 / zgsp4
3788	zgsp6 = (zgsp2 - 1.0_wp)/gamma
3789
3790	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)"
3791	WRITE(numtan,FMT) 'eosinsp ', cur_loop, h_ratio, zgsp1, zgsp2, zgsp3, zgsp4, zgsp5, zgsp6, zgsp7
3792	!--------------------------------------------------------------------
3793	! Unitary calculus
3794	!--------------------------------------------------------------------
3795	FMT = "(A8,2X,A8,2X,I4.4,2X,E6.1,2X,I4.4,2X,I4.4,2X,I4.4,2X,E20.13,1X)"
3796	cl_name = 'eosinsp '
3797	IF(lwp) THEN
3798	DO ii=1, 100, 1
3799	IF ( zscrd(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscrd ', &
3800	& cur_loop, h_ratio, iiposrd(ii), ijposrd(ii), ikposrd(ii), zscrd(ii)
3801	ENDDO
3802	DO ii=1, 100, 1
3803	IF ( zscrh(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscrh ', &
3804	& cur_loop, h_ratio, iiposrh(ii), ijposrh(ii), ikposrh(ii), zscrh(ii)
3805	ENDDO
3806	DO ii=1, 100, 1
3807	IF ( zscerrrd(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscerrrd', &
3808	& cur_loop, h_ratio, iiposrd(ii), ijposrd(ii), ikposrd(ii), zscerrrd(ii)
3809	ENDDO
3810	DO ii=1, 100, 1
3811	IF ( zscerrrh(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscerrrh', &
3812	& cur_loop, h_ratio, iiposrh(ii), ijposrh(ii), ikposrh(ii), zscerrrh(ii)
3813	ENDDO
3814	! write separator
3815	WRITE(numtan_sc,"(A4)") '===='
3816	ENDIF
3817	ENDIF
3818
3819	DEALLOCATE( &
3820	& zrd_out, zrd_tl, zrd_wop, &
3821	& zrh_out, zrh_tl, zrh_wop, &
3822	& zt_tlin, zs_tlin, z3r &
3823	& )
3824	END SUBROUTINE eos_insitu_pot_tlm_tst
3825
3826	SUBROUTINE eos_insitu_2d_tlm_tst( kumadt )
3827	!!-----------------------------------------------------------------------
3828	!!
3829	!! * ROUTINE eos_insitu_2d_tlm_tst *
3830	!!
3831	!! ** Purpose : Test the tangent routine.
3832	!!
3833	!! ** Method : Verify the tangent with Taylor expansion
3834	!!
3835	!! M(x+hdx) = M(x) + L(hdx) + O(h^2)
3836	!!
3837	!! where L = tangent routine
3838	!! M = direct routine
3839	!! dx = input perturbation (random field)
3840	!! h = ration on perturbation
3841	!!
3842	!! In the tangent test we verify that:
3843	!! M(x+h*dx) - M(x)
3844	!! g(h) = ------------------ ---> 1 as h ---> 0
3845	!! L(h*dx)
3846	!! and
3847	!! g(h) - 1
3848	!! f(h) = ---------- ---> k (costant) as h ---> 0
3849	!! p
3850	!!
3851	!! History :
3852	!! ! 09-08 (A. Vigilant)
3853	!!-----------------------------------------------------------------------
3854	!! * Modules used
3855	USE eosbn2, ONLY: & ! horizontal & vertical advective trend
3856	& eos
3857	USE tamtrj ! writing out state trajectory
3858	USE par_tlm, ONLY: &
3859	& tlm_bch, &
3860	& cur_loop, &
3861	& h_ratio
3862	USE istate_mod
3863	USE gridrandom, ONLY: &
3864	& grid_rd_sd
3865	USE trj_tam
3866	USE oce , ONLY: & ! ocean dynamics and tracers variables
3867	& tn, sn, rhd, rhop
3868	USE in_out_manager, ONLY: & ! I/O manager
3869	& nitend, &
3870	& nit000
3871	USE tamctl, ONLY: & ! Control parameters
3872	& numtan, numtan_sc
3873	!! * Arguments
3874	INTEGER, INTENT(IN) :: &
3875	& kumadt ! Output unit
3876	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
3877	& zdep, & ! depth
3878	& ztem, & ! potential temperature
3879	& zsal ! salinity
3880	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
3881	& zrd_out, & !
3882	& zt_tlin , & !
3883	& zs_tlin , &
3884	& zrd_tl , &
3885	& zrd_wop , &
3886	& z2r
3887	REAL(KIND=wp) :: &
3888	& zsp1, &
3889	& zsp2, &
3890	& zsp3, &
3891	& zzsp, &
3892	& gamma, &
3893	& zgsp1, &
3894	& zgsp2, &
3895	& zgsp3, &
3896	& zgsp4, &
3897	& zgsp5, &
3898	& zgsp6, &
3899	& zgsp7
3900	INTEGER :: &
3901	& ji, &
3902	& jj
3903	CHARACTER(LEN=14) :: cl_name
3904	CHARACTER (LEN=128) :: file_out, file_wop, file_xdx
3905	CHARACTER (LEN=90) :: FMT
3906	REAL(KIND=wp), DIMENSION(100):: &
3907	& zscrd, zscerrrd
3908	INTEGER, DIMENSION(100):: &
3909	& iiposrd, ijposrd
3910	INTEGER:: &
3911	& ii, numsctlm, &
3912	& isamp=40,jsamp=40
3913	REAL(KIND=wp), DIMENSION(jpi,jpj) :: &
3914	& zerrrd
3915	ALLOCATE( &
3916	& ztem ( jpi, jpj ), &
3917	& zsal ( jpi, jpj ), &
3918	& zdep ( jpi, jpj ), &
3919	& zrd_out( jpi, jpj ), &
3920	& zrd_tl( jpi, jpj ), &
3921	& zs_tlin( jpi, jpj ), &
3922	& zt_tlin( jpi, jpj ), &
3923	& zrd_wop( jpi, jpj ), &
3924	& z2r ( jpi, jpj ) )
3925
3926	!--------------------------------------------------------------------
3927	! Reset the tangent and adjoint variables
3928	!--------------------------------------------------------------------
3929	ztem ( :,:) = 0.0_wp
3930	zsal ( :,:) = 0.0_wp
3931	zt_tlin( :,:) = 0.0_wp
3932	zs_tlin( :,:) = 0.0_wp
3933	zrd_out( :,:) = 0.0_wp
3934	zrd_wop( :,:) = 0.0_wp
3935	zscerrrd( :) = 0.0_wp
3936	zscrd(:) = 0.0_wp
3937	IF ( tlm_bch == 2 ) zrd_tl ( :,:) = 0.0_wp
3938	!--------------------------------------------------------------------
3939	! Output filename Xn=F(X0)
3940	!--------------------------------------------------------------------
3941	!! CALL tlm_namrd
3942	gamma = h_ratio
3943	file_wop='trj_wop_eos_2d'
3944	file_xdx='trj_xdx_eos_2d'
3945	!--------------------------------------------------------------------
3946	! Initialize the tangent input with random noise: dx
3947	!--------------------------------------------------------------------
3948	IF (( cur_loop .NE. 0) .OR. ( gamma .NE. 0.0_wp) )THEN
3949	CALL grid_rd_sd( 596035, z2r, 'T', 0.0_wp, stdt)
3950	DO jj = nldj, nlej
3951	DO ji = nldi, nlei
3952	zt_tlin(ji,jj) = z2r(ji,jj)
3953	END DO
3954	END DO
3955	CALL grid_rd_sd( 371836, z2r, 'S', 0.0_wp, stds)
3956	DO jj = nldj, nlej
3957	DO ji = nldi, nlei
3958	zs_tlin(ji,jj) = z2r(ji,jj)
3959	END DO
3960	END DO
3961	ENDIF
3962
3963	!--------------------------------------------------------------------
3964	! Complete Init for Direct
3965	!-------------------------------------------------------------------
3966	IF ( tlm_bch /= 2 ) CALL istate_p
3967	! *** initialize the reference trajectory
3968	! ------------
3969	CALL trj_rea( nit000-1, 1 )
3970	CALL trj_rea( nit000, 1 )
3971	! Initialize the reference state
3972	ztem(:,:) = tn(:,:,2)
3973	zsal(:,:) = sn(:,:,2)
3974	zdep(:,:) = fsdept(:,:,2)
3975
3976	IF (( cur_loop .NE. 0) .OR. ( gamma .NE. 0.0_wp) )THEN
3977	zt_tlin(:,:) = gamma * zt_tlin(:,:)
3978	ztem(:,:) = ztem(:,:) + zt_tlin(:,:)
3979
3980	zs_tlin(:,:) = gamma * zs_tlin(:,:)
3981	zsal(:,:) = zsal(:,:) + zs_tlin(:,:)
3982	ENDIF
3983	!--------------------------------------------------------------------
3984	! Compute the direct model F(X0,t=n) = Xn
3985	!--------------------------------------------------------------------
3986	IF ( tlm_bch /= 2 ) CALL eos(ztem, zsal, zdep, zrd_out)
3987	rhd (:,:,2) = zrd_out(:,:)
3988	IF ( tlm_bch == 0 ) CALL trj_wri_spl(file_wop)
3989	IF ( tlm_bch == 1 ) CALL trj_wri_spl(file_xdx)
3990	!--------------------------------------------------------------------
3991	! Compute the Tangent
3992	!--------------------------------------------------------------------
3993	IF ( tlm_bch == 2 ) THEN
3994	!--------------------------------------------------------------------
3995	! Initialize the tangent variables: dy^* = W dy
3996	!--------------------------------------------------------------------
3997	CALL trj_rea( nit000-1, 1 )
3998	CALL trj_rea( nit000, 1 )
3999	ztem(:,:) = tn(:,:,2)
4000	zsal(:,:) = sn(:,:,2)
4001	!-----------------------------------------------------------------------
4002	! Initialization of the dynamics and tracer fields for the tangent
4003	!-----------------------------------------------------------------------
4004	CALL eos_insitu_2d_tan(ztem, zsal, zdep, zt_tlin, zs_tlin, zrd_tl)
4005	!--------------------------------------------------------------------
4006	! Compute the scalar product: ( L(t0,tn) gamma dx0 ) )
4007	!--------------------------------------------------------------------
4008	zsp2 = DOT_PRODUCT( zrd_tl, zrd_tl )
4009	!--------------------------------------------------------------------
4010	! Storing data
4011	!--------------------------------------------------------------------
4012	CALL trj_rd_spl(file_wop)
4013	zrd_wop (:,:) = rhd (:,:,2)
4014	CALL trj_rd_spl(file_xdx)
4015	zrd_out (:,:) = rhd (:,:,2)
4016	!--------------------------------------------------------------------
4017	! Compute the Linearization Error
4018	! Nn = M( X0+gamma.dX0, t0,tn) - M(X0, t0,tn)
4019	! and
4020	! Compute the Linearization Error
4021	! En = Nn -TL(gamma.dX0, t0,tn)
4022	!--------------------------------------------------------------------
4023	! Warning: Here we re-use local variables z()_out and z()_wop
4024	ii=0
4025	DO jj = 1, jpj
4026	DO ji = 1, jpi
4027	zrd_out (ji,jj) = zrd_out (ji,jj) - zrd_wop (ji,jj)
4028	zrd_wop (ji,jj) = zrd_out (ji,jj) - zrd_tl (ji,jj)
4029	IF ( zrd_tl(ji,jj) .NE. 0.0_wp ) &
4030	& zerrrd(ji,jj) = zrd_out(ji,jj)/zrd_tl(ji,jj)
4031	IF( (MOD(ji, isamp) .EQ. 0) .AND. &
4032	& (MOD(jj, jsamp) .EQ. 0) ) THEN
4033	ii = ii+1
4034	iiposrd(ii) = ji
4035	ijposrd(ii) = jj
4036	IF ( INT(tmask(ji,jj,2)) .NE. 0) THEN
4037	zscrd (ii) = zrd_wop(ji,jj)
4038	zscerrrd (ii) = ( zerrrd( ji,jj) - 1.0_wp ) / gamma
4039	ENDIF
4040	ENDIF
4041	END DO
4042	END DO
4043
4044	zsp1 = DOT_PRODUCT( zrd_out, zrd_out )
4045
4046	zsp3 = DOT_PRODUCT( zrd_wop, zrd_wop )
4047	!--------------------------------------------------------------------
4048	! Print the linearization error En - norme 2
4049	!--------------------------------------------------------------------
4050	! 14 char:'12345678901234'
4051	cl_name = 'eos_2d:En '
4052	zzsp = SQRT(zsp3)
4053	zgsp5 = zzsp
4054	CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio )
4055	!--------------------------------------------------------------------
4056	! Compute TLM norm2
4057	!--------------------------------------------------------------------
4058	zzsp = SQRT(zsp2)
4059	zgsp4 = zzsp
4060	cl_name = 'eos_2d:Ln2'
4061	CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio )
4062	!--------------------------------------------------------------------
4063	! Print the linearization error Nn - norme 2
4064	!--------------------------------------------------------------------
4065	zzsp = SQRT(zsp1)
4066	cl_name = 'eos_2d:Mhdx-Mx'
4067	CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio )
4068
4069	zgsp3 = SQRT( zsp3/zsp2 )
4070	zgsp7 = zgsp3/gamma
4071	zgsp1 = zzsp
4072	zgsp2 = zgsp1 / zgsp4
4073	zgsp6 = (zgsp2 - 1.0_wp)/gamma
4074
4075	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)"
4076	WRITE(numtan,FMT) 'eosins2d', cur_loop, h_ratio, zgsp1, zgsp2, zgsp3, zgsp4, zgsp5, zgsp6, zgsp7
4077	!--------------------------------------------------------------------
4078	! Unitary calculus
4079	!--------------------------------------------------------------------
4080	FMT = "(A8,2X,A8,2X,I4.4,2X,E6.1,2X,I4.4,2X,I4.4,2X,I4.4,2X,E20.13,1X)"
4081	cl_name = 'eosins2d'
4082	IF(lwp) THEN
4083	DO ii=1, 100, 1
4084	IF ( zscrd(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscrd ', &
4085	& cur_loop, h_ratio, ii, iiposrd(ii), ijposrd(ii), zscrd(ii)
4086	ENDDO
4087	DO ii=1, 100, 1
4088	IF ( zscerrrd(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscerrrd', &
4089	& cur_loop, h_ratio, ii, iiposrd(ii), ijposrd(ii), zscerrrd(ii)
4090	ENDDO
4091	! write separator
4092	WRITE(numtan_sc,"(A4)") '===='
4093	ENDIF
4094
4095	ENDIF
4096
4097	DEALLOCATE( &
4098	& zrd_out, zrd_tl, zrd_wop, &
4099	& ztem, zsal, zdep, &
4100	& zt_tlin, zs_tlin, z2r &
4101	& )
4102	END SUBROUTINE eos_insitu_2d_tlm_tst
4103
4104	SUBROUTINE eos_pot_1pt_tlm_tst( kumadt )
4105	!!-----------------------------------------------------------------------
4106	!!
4107	!! * ROUTINE eos_pot_1pt_tlm_tst *
4108	!!
4109	!! ** Purpose : Test the tangent routine.
4110	!!
4111	!! ** Method : Verify the tangent with Taylor expansion
4112	!!
4113	!! M(x+hdx) = M(x) + L(hdx) + O(h^2)
4114	!!
4115	!! where L = tangent routine
4116	!! M = direct routine
4117	!! dx = input perturbation (random field)
4118	!! h = ration on perturbation
4119	!!
4120	!! In the tangent test we verify that:
4121	!! M(x+h*dx) - M(x)
4122	!! g(h) = ------------------ ---> 1 as h ---> 0
4123	!! L(h*dx)
4124	!! and
4125	!! g(h) - 1
4126	!! f(h) = ---------- ---> k (costant) as h ---> 0
4127	!! p
4128	!!
4129	!! History :
4130	!! ! 09-08 (A. Vigilant)
4131	!!-----------------------------------------------------------------------
4132	!! * Modules used
4133	USE eosbn2, ONLY: & ! horizontal & vertical advective trend
4134	& eos
4135	USE tamtrj ! writing out state trajectory
4136	USE par_tlm, ONLY: &
4137	& tlm_bch, &
4138	& cur_loop, &
4139	& h_ratio
4140	USE istate_mod
4141	USE trj_tam
4142	USE oce , ONLY: & ! ocean dynamics and tracers variables
4143	& tn, sn, rhd, rhop
4144	USE in_out_manager, ONLY: & ! I/O manager
4145	& nitend, &
4146	& nit000
4147	USE tamctl, ONLY: & ! Control parameters
4148	& numtan, numtan_sc
4149	!! * Arguments
4150	INTEGER, INTENT(IN) :: &
4151	& kumadt ! Output unit
4152	INTEGER:: &
4153	& numsctlm
4154	!! * Local declarations
4155	REAL(KIND=wp) :: &
4156	& ztem, &
4157	& zsal, &
4158	& zt_tlin, &
4159	& zs_tlin, &
4160	& zrh_out, &
4161	& zrh_tl, &
4162	& zrh_wop
4163	REAL(KIND=wp) :: &
4164	& zsp1, zsp2, &
4165	& zsp3, zzsp, &
4166	& gamma, &
4167	& zgsp1, &
4168	& zgsp2, &
4169	& zgsp3, &
4170	& zgsp4, &
4171	& zgsp5, &
4172	& zgsp6, &
4173	& zgsp7
4174	CHARACTER(LEN=14) :: cl_name
4175	CHARACTER (LEN=128) :: file_out, file_wop, file_xdx
4176	CHARACTER (LEN=90) :: FMT
4177	! Initialize the reference state
4178	ztem = 23.7
4179	zsal = 30.1
4180
4181	!--------------------------------------------------------------------
4182	! Reset the tangentvariables
4183	!--------------------------------------------------------------------
4184	zt_tlin = 1.12_wp
4185	zs_tlin = 0.123_wp
4186	zrh_out = 0.0_wp
4187	zrh_wop = 0.0_wp
4188	IF ( tlm_bch == 2 ) zrh_tl = 0.0_wp
4189	!--------------------------------------------------------------------
4190	! Output filename Xn=F(X0)
4191	!--------------------------------------------------------------------
4192	!! CALL tlm_namrd
4193	gamma = h_ratio
4194	file_wop='trj_wop_eos_1pt'
4195	file_xdx='trj_xdx_eos_1pt'
4196	!--------------------------------------------------------------------
4197	! Complete Init for Direct
4198	!-------------------------------------------------------------------
4199	IF ( tlm_bch /= 2 ) CALL istate_p
4200
4201	! *** initialize the reference trajectory
4202	! ------------
4203	CALL trj_rea( nit000-1, 1 )
4204	CALL trj_rea( nit000, 1 )
4205
4206	IF (( cur_loop .NE. 0) .OR. ( gamma .NE. 0.0_wp) )THEN
4207	zt_tlin = gamma * zt_tlin
4208	ztem = ztem + zt_tlin
4209
4210	zs_tlin = gamma * zs_tlin
4211	zsal = zsal + zs_tlin
4212	ENDIF
4213	!--------------------------------------------------------------------
4214	! Compute the direct model F(X0,t=n) = Xn
4215	!--------------------------------------------------------------------
4216	IF ( tlm_bch /= 2 ) CALL eos(ztem, zsal, zrh_out)
4217	rhop (1,1,1) = zrh_out
4218	IF ( tlm_bch == 0 ) CALL trj_wri_spl(file_wop)
4219	IF ( tlm_bch == 1 ) CALL trj_wri_spl(file_xdx)
4220	!--------------------------------------------------------------------
4221	! Compute the Tangent
4222	!--------------------------------------------------------------------
4223	IF ( tlm_bch == 2 ) THEN
4224	!--------------------------------------------------------------------
4225	! Initialize the tangent variables: dy^* = W dy
4226	!--------------------------------------------------------------------
4227	CALL trj_rea( nit000-1, 1 )
4228	CALL trj_rea( nit000, 1 )
4229	ztem = 23.7
4230	zsal = 30.1
4231	!-----------------------------------------------------------------------
4232	! Initialization of the dynamics and tracer fields for the tangent
4233	!-----------------------------------------------------------------------
4234	CALL eos_pot_1pt_tan(ztem, zsal, zt_tlin, zs_tlin, zrh_tl)
4235	!--------------------------------------------------------------------
4236	! Compute the scalar product: ( L(t0,tn) gamma dx0 ) )
4237	!--------------------------------------------------------------------
4238	zsp2 = abs(zrh_tl)
4239	!--------------------------------------------------------------------
4240	! Storing data
4241	!--------------------------------------------------------------------
4242	CALL trj_rd_spl(file_wop)
4243	zrh_wop = rhop (1,1,1)
4244	CALL trj_rd_spl(file_xdx)
4245	zrh_out = rhop (1,1,1)
4246	!--------------------------------------------------------------------
4247	! Compute the Linearization Error
4248	! Nn = M( X0+gamma.dX0, t0,tn) - M(X0, t0,tn)
4249	! and
4250	! Compute the Linearization Error
4251	! En = Nn -TL(gamma.dX0, t0,tn)
4252	!--------------------------------------------------------------------
4253	! Warning: Here we re-use local variables z()_out and z()_wop
4254	zrh_out = zrh_out - zrh_wop
4255	zrh_wop = zrh_out - zrh_tl
4256
4257	zsp1 = abs(zrh_out)
4258	zsp3 = abs(zrh_wop)
4259	!--------------------------------------------------------------------
4260	! Print the linearization error En - norme 2
4261	!--------------------------------------------------------------------
4262	! 14 char:'12345678901234'
4263	cl_name = 'eos_1pt:En '
4264	zzsp = zsp3
4265	zgsp5 = zzsp
4266	CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio )
4267	!--------------------------------------------------------------------
4268	! Compute TLM norm2
4269	!--------------------------------------------------------------------
4270	zzsp = zsp2
4271	zgsp4 = zzsp
4272	cl_name = 'eos_1pt:Ln2'
4273	CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio )
4274	!--------------------------------------------------------------------
4275	! Print the linearization error Nn - norme 2
4276	!--------------------------------------------------------------------
4277	zzsp = zsp1
4278	cl_name = 'eos_1pt:Mhdx-Mx'
4279	CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio )
4280
4281	zgsp3 = SQRT( zsp3/zsp2 )
4282	zgsp7 = zgsp3/gamma
4283	zgsp1 = zzsp
4284	zgsp2 = zgsp1 / zgsp4
4285	zgsp6 = (zgsp2 - 1.0_wp)/gamma
4286
4287	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)"
4288	WRITE(numtan,FMT) 'eos1pt ', cur_loop, h_ratio, zgsp1, zgsp2, zgsp3, zgsp4, zgsp5, zgsp6, zgsp7
4289	ENDIF
4290
4291	END SUBROUTINE eos_pot_1pt_tlm_tst
4292
4293	SUBROUTINE bn2_tlm_tst( kumadt )
4294	!!-----------------------------------------------------------------------
4295	!!
4296	!! * ROUTINE bn2_tlm_tst *
4297	!!
4298	!! ** Purpose : Test the tangent routine.
4299	!!
4300	!! ** Method : Verify the tangent with Taylor expansion
4301	!!
4302	!! M(x+hdx) = M(x) + L(hdx) + O(h^2)
4303	!!
4304	!! where L = tangent routine
4305	!! M = direct routine
4306	!! dx = input perturbation (random field)
4307	!! h = ration on perturbation
4308	!!
4309	!! In the tangent test we verify that:
4310	!! M(x+h*dx) - M(x)
4311	!! g(h) = ------------------ ---> 1 as h ---> 0
4312	!! L(h*dx)
4313	!! and
4314	!! g(h) - 1
4315	!! f(h) = ---------- ---> k (costant) as h ---> 0
4316	!! p
4317	!!
4318	!! History :
4319	!! ! 09-12 (A. Vigilant)
4320	!!-----------------------------------------------------------------------
4321	!! * Modules used
4322	USE eosbn2, ONLY: & ! horizontal & vertical advective trend
4323	& bn2
4324	USE tamtrj ! writing out state trajectory
4325	USE par_tlm, ONLY: &
4326	& tlm_bch, &
4327	& cur_loop, &
4328	& h_ratio
4329	USE istate_mod
4330	USE sshwzv ! vertical velocity
4331	USE gridrandom, ONLY: &
4332	& grid_rd_sd
4333	USE trj_tam
4334	USE oce , ONLY: & ! ocean dynamics and tracers variables
4335	& tn, sn, rn2
4336	USE oce_tam , ONLY: &
4337	& tn_tl, &
4338	& sn_tl, &
4339	& rn2_tl
4340	USE in_out_manager, ONLY: & ! I/O manager &
4341	& nit000
4342	USE tamctl, ONLY: & ! Control parameters
4343	& numtan, numtan_sc
4344	!! * Arguments
4345	INTEGER, INTENT(IN) :: &
4346	& kumadt ! Output unit
4347	!! * Local declarations
4348	INTEGER :: &
4349	& ji, & ! dummy loop indices
4350	& jj, &
4351	& jk
4352	REAL(KIND=wp) :: &
4353	& zsp1, & ! scalar product involving the tangent routine
4354	& zsp2, & ! scalar product involving the tangent routine
4355	& zsp3, & ! scalar product involving the tangent routine
4356	& zzsp, & ! scalar product involving the tangent routine
4357	& gamma, &
4358	& zgsp1, &
4359	& zgsp2, &
4360	& zgsp3, &
4361	& zgsp4, &
4362	& zgsp5, &
4363	& zgsp6, &
4364	& zgsp7
4365	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
4366	& ztn_tlin, & ! potential temperature
4367	& zsn_tlin, & ! salinity
4368	& zrn2_out, & ! Brunt-Vaisala frequency [s-1] Direct output
4369	& zrn2_wop, & ! Brunt-Vaisala frequency [s-1] Direct output w/o perturbation
4370	& z3r
4371	CHARACTER(LEN=14) :: cl_name
4372	CHARACTER (LEN=128) :: file_out, file_wop, file_xdx
4373	CHARACTER (LEN=90) :: FMT
4374	REAL(KIND=wp), DIMENSION(100):: &
4375	& zscrn2, &
4376	& zscerrrn2
4377	INTEGER, DIMENSION(100):: &
4378	& iiposrn2,ijposrn2, ikposrn2
4379	INTEGER:: &
4380	& ii, &
4381	& isamp=40, &
4382	& jsamp=40, &
4383	& ksamp=10, &
4384	& numsctlm
4385	REAL(KIND=wp), DIMENSION(jpi,jpj,jpk) :: &
4386	& zerrrn2
4387
4388	! Allocate memory
4389	ALLOCATE( &
4390	& zsn_tlin( jpi, jpj, jpk ), &
4391	& ztn_tlin( jpi, jpj, jpk ), &
4392	& zrn2_out( jpi, jpj, jpk ), &
4393	& zrn2_wop( jpi, jpj, jpk ), &
4394	& z3r( jpi, jpj, jpk ) )
4395
4396
4397	!--------------------------------------------------------------------
4398	! Output filename Xn=F(X0)
4399	!--------------------------------------------------------------------
4400	!! CALL tlm_namrd
4401	gamma = h_ratio
4402	file_wop='trj_wop_bn2'
4403	file_xdx='trj_xdx_bn2'
4404	!--------------------------------------------------------------------
4405	! Initialize the tangent input with random noise: dx
4406	!--------------------------------------------------------------------
4407	IF (( cur_loop .NE. 0) .OR. ( gamma .NE. 0.0_wp) )THEN
4408	CALL grid_rd_sd( 264940, z3r, 'T', 0.0_wp, stdt)
4409	DO jk = 1, jpk
4410	DO jj = nldj, nlej
4411	DO ji = nldi, nlei
4412	ztn_tlin(ji,jj,jk) = z3r(ji,jj,jk)
4413	END DO
4414	END DO
4415	END DO
4416	CALL grid_rd_sd( 618304, z3r, 'T', 0.0_wp, stds)
4417	DO jk = 1, jpk
4418	DO jj = nldj, nlej
4419	DO ji = nldi, nlei
4420	zsn_tlin(ji,jj,jk) = z3r(ji,jj,jk)
4421	END DO
4422	END DO
4423	END DO
4424	ENDIF
4425	!--------------------------------------------------------------------
4426	! Complete Init for Direct
4427	!-------------------------------------------------------------------
4428	IF ( tlm_bch /= 2 ) CALL istate_p
4429
4430	! *** initialize the reference trajectory
4431	! ------------
4432	CALL trj_rea( nit000-1, 1 )
4433
4434	IF (( cur_loop .NE. 0) .OR. ( gamma .NE. 0.0_wp) )THEN
4435
4436	ztn_tlin(:,:,:) = gamma * ztn_tlin(:,:,:)
4437	tn(:,:,:) = tn(:,:,:) + ztn_tlin(:,:,:)
4438
4439	zsn_tlin(:,:,:) = gamma * zsn_tlin(:,:,:)
4440	sn(:,:,:) = sn(:,:,:) + zsn_tlin(:,:,:)
4441
4442	ENDIF
4443
4444	!--------------------------------------------------------------------
4445	! Compute the direct model F(X0,t=n) = Xn
4446	!--------------------------------------------------------------------
4447	IF ( tlm_bch /= 2 ) CALL bn2(tn, sn, rn2)
4448	IF ( tlm_bch == 0 ) CALL trj_wri_spl(file_wop)
4449	IF ( tlm_bch == 1 ) CALL trj_wri_spl(file_xdx)
4450	!--------------------------------------------------------------------
4451	! Compute the Tangent
4452	!--------------------------------------------------------------------
4453	IF ( tlm_bch == 2 ) THEN
4454	!--------------------------------------------------------------------
4455	! Initialize the tangent variables: dy^* = W dy
4456	!--------------------------------------------------------------------
4457	CALL trj_rea( nit000-1, 1 )
4458	tn_tl (:,:,:) = ztn_tlin (:,:,:)
4459	sn_tl (:,:,:) = zsn_tlin (:,:,:)
4460
4461	!-----------------------------------------------------------------------
4462	! Initialization of the dynamics and tracer fields for the tangent
4463	!-----------------------------------------------------------------------
4464	CALL bn2_tan(tn, sn, ztn_tlin, zsn_tlin, rn2_tl)
4465
4466	!--------------------------------------------------------------------
4467	! Compute the scalar product: ( L(t0,tn) gamma dx0 ) )
4468	!--------------------------------------------------------------------
4469	zsp2 = DOT_PRODUCT( rn2_tl, rn2_tl )
4470
4471	!--------------------------------------------------------------------
4472	! Storing data
4473	!--------------------------------------------------------------------
4474	CALL trj_rd_spl(file_wop)
4475	zrn2_wop (:,:,:) = rn2 (:,:,:)
4476	CALL trj_rd_spl(file_xdx)
4477	zrn2_out (:,:,:) = rn2 (:,:,:)
4478	!--------------------------------------------------------------------
4479	! Compute the Linearization Error
4480	! Nn = M( X0+gamma.dX0, t0,tn) - M(X0, t0,tn)
4481	! and
4482	! Compute the Linearization Error
4483	! En = Nn -TL(gamma.dX0, t0,tn)
4484	!--------------------------------------------------------------------
4485	! Warning: Here we re-use local variables z()_out and z()_wop
4486	ii=0
4487	DO jk = 1, jpk
4488	DO jj = 1, jpj
4489	DO ji = 1, jpi
4490	zrn2_out (ji,jj,jk) = zrn2_out (ji,jj,jk) - zrn2_wop (ji,jj,jk)
4491	zrn2_wop (ji,jj,jk) = zrn2_out (ji,jj,jk) - rn2_tl (ji,jj,jk)
4492	IF ( rn2_tl(ji,jj,jk) .NE. 0.0_wp ) &
4493	& zerrrn2(ji,jj,jk) = zrn2_out(ji,jj,jk)/rn2_tl(ji,jj,jk)
4494	IF( (MOD(ji, isamp) .EQ. 0) .AND. &
4495	& (MOD(jj, jsamp) .EQ. 0) .AND. &
4496	& (MOD(jk, ksamp) .EQ. 0) ) THEN
4497	ii = ii+1
4498	iiposrn2(ii) = ji
4499	ijposrn2(ii) = jj
4500	ikposrn2(ii) = jk
4501	IF ( INT(tmask(ji,jj,jk)) .NE. 0) THEN
4502	zscrn2 (ii) = zrn2_wop(ji,jj,jk)
4503	zscerrrn2 (ii) = ( zerrrn2(ji,jj,jk) - 1.0_wp ) / gamma
4504	ENDIF
4505	ENDIF
4506	END DO
4507	END DO
4508	END DO
4509
4510	zsp1 = DOT_PRODUCT( zrn2_out, zrn2_out )
4511
4512	zsp3 = DOT_PRODUCT( zrn2_wop, zrn2_wop )
4513
4514	!--------------------------------------------------------------------
4515	! Print the linearization error En - norme 2
4516	!--------------------------------------------------------------------
4517	! 14 char:'12345678901234'
4518	cl_name = 'eosbn2_tam:En '
4519	zzsp = SQRT(zsp3)
4520	zgsp5 = zzsp
4521	CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio )
4522	!--------------------------------------------------------------------
4523	! Compute TLM norm2
4524	!--------------------------------------------------------------------
4525	zzsp = SQRT(zsp2)
4526	zgsp4 = zzsp
4527	cl_name = 'eosbn2_tam:Ln2'
4528	CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio )
4529	!--------------------------------------------------------------------
4530	! Print the linearization error Nn - norme 2
4531	!--------------------------------------------------------------------
4532	zzsp = SQRT(zsp1)
4533	cl_name = 'traadv:Mhdx-Mx'
4534	CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio )
4535
4536	zgsp3 = SQRT( zsp3/zsp2 )
4537	zgsp7 = zgsp3/gamma
4538	zgsp1 = zzsp
4539	zgsp2 = zgsp1 / zgsp4
4540	zgsp6 = (zgsp2 - 1.0_wp)/gamma
4541
4542	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)"
4543	WRITE(numtan,FMT) 'eosbn2 ', cur_loop, h_ratio, zgsp1, zgsp2, zgsp3, zgsp4, zgsp5, zgsp6, zgsp7
4544
4545	!--------------------------------------------------------------------
4546	! Unitary calculus
4547	!--------------------------------------------------------------------
4548	FMT = "(A8,2X,A8,2X,I4.4,2X,E6.1,2X,I4.4,2X,I4.4,2X,I4.4,2X,E20.13,1X)"
4549	cl_name = 'eosbn2 '
4550	IF(lwp) THEN
4551	DO ii=1, 100, 1
4552	IF ( zscrn2(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscrn2 ', &
4553	& cur_loop, h_ratio, ii, iiposrn2(ii), ijposrn2(ii), zscrn2(ii)
4554	ENDDO
4555	! write separator
4556	WRITE(numtan_sc,"(A4)") '===='
4557	ENDIF
4558
4559	ENDIF
4560
4561	DEALLOCATE( &
4562	& ztn_tlin, zsn_tlin, &
4563	& zrn2_out, zrn2_wop, &
4564	& z3r &
4565	& )
4566
4567	END SUBROUTINE bn2_tlm_tst
4568
4569	SUBROUTINE eos_tlm_tst( kumadt )
4570	!!-----------------------------------------------------------------------
4571	!!
4572	!! * ROUTINE eos_tlm_tst *
4573	!!
4574	!! ** Purpose : Test the tangent routine.
4575	!!
4576	!! History :
4577	!! ! 09-08 (F. Vigilant)
4578	!!-----------------------------------------------------------------------
4579	!! * Modules used
4580	!! * Arguments
4581	INTEGER, INTENT(IN) :: &
4582	& kumadt ! Output unit
4583
4584	CALL eos_insitu_tlm_tst( kumadt )
4585	CALL eos_insitu_pot_tlm_tst( kumadt )
4586	CALL eos_insitu_2d_tlm_tst( kumadt )
4587	CALL eos_pot_1pt_tlm_tst( kumadt )
4588
4589	END SUBROUTINE eos_tlm_tst
4590	#endif
4591	#endif
4592	!!======================================================================
4593	END MODULE eosbn2_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_2_2/NEMOTAM/OPATAM_SRC/TRA/eosbn2_tam.F90 @ 2578

Download in other formats: