Context Navigation

eosbn2_tam.F90 @ 2587

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

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

Context Navigation

source: branches/TAM_V3_0/NEMOTAM/OPATAM_SRC/eosbn2_tam.F90 @ 2587