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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

eosbn2_tam.F90 in branches/2012/dev_r3604_LEGI8_TAM/NEMOGCM/NEMO/OPATAM_SRC/TRA – NEMO

Context Navigation

source: branches/2012/dev_r3604_LEGI8_TAM/NEMOGCM/NEMO/OPATAM_SRC/TRA/eosbn2_tam.F90 @ 3627

Last change on this file since 3627 was 3627, checked in by rblod, 11 years ago
Correct long lines in TAM 4.3 see ticket #1010
Property svn:executable set to ``*
File size: 140.3 KB

Line
1	MODULE eosbn2_tam
2	!!==============================================================================
3	!! * MODULE eosbn2_tam *
4	!! Ocean diagnostic variable : equation of state - in situ and potential density
5	!! - Brunt-Vaisala frequency
6	!! Tangent and Adjoint module
7	!!==============================================================================
8	!! History : OPA ! 1989-03 (O. Marti) Original code
9	!! 6.0 ! 1994-07 (G. Madec, M. Imbard) add bn2
10	!! 6.0 ! 1994-08 (G. Madec) Add Jackett & McDougall eos
11	!! 7.0 ! 1996-01 (G. Madec) statement function for e3
12	!! 8.1 ! 1997-07 (G. Madec) density instead of volumic mass
13	!! - ! 1999-02 (G. Madec, N. Grima) semi-implicit pressure gradient
14	!! 8.2 ! 2001-09 (M. Ben Jelloul) bugfix on linear eos
15	!! NEMO 1.0 ! 2002-10 (G. Madec) add eos_init
16	!! - ! 2002-11 (G. Madec, A. Bozec) partial step, eos_insitu_2d
17	!! - ! 2003-08 (G. Madec) F90, free form
18	!! 3.0 ! 2006-08 (G. Madec) add tfreez function
19	!! 3.3 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA
20	!! - ! 2010-10 (G. Nurser, G. Madec) add eos_alpbet used in ldfslp
21	!! History of the TAM Module :
22	!! 8.2 ! 2005-03 (F. Van den Berghe, A. Weaver, N. Daget)
23	!! - eostan.F
24	!! 9.0 ! 2007-07 (K. Mogensen) Initial version based on eostan.F
25	!! ! 2008-07 (A. Vidard) bug fix in computation of prd_tl if neos=1
26	!! ! 2008-11 (A. Vidard) TAM of the 06-08 version
27	!! NEMO 3.2 ! 2010-04 (F. Vigilant) version 3.2
28	!! 3.4 ! 2012-04 (P.-A. Bouttier) version 3.4
29	!!----------------------------------------------------------------------
30
31	!!----------------------------------------------------------------------
32	!! Direct subroutines
33	!! eos : generic interface of the equation of state
34	!! eos_insitu : Compute the in situ density
35	!! eos_insitu_pot : Compute the insitu and surface referenced potential
36	!! volumic mass
37	!! eos_insitu_2d : Compute the in situ density for 2d fields
38	!! eos_bn2 : Compute the Brunt-Vaisala frequency
39	!! eos_alpbet : calculates the in situ thermal/haline expansion ratio
40	!! tfreez : Compute the surface freezing temperature
41	!! eos_init : set eos parameters (namelist)
42	!!----------------------------------------------------------------------
43	!! * Modules used
44	#if defined key_zdfddm
45	USE oce_tam
46	#endif
47	USE dom_oce
48	USE par_kind
49	USE par_oce
50	USE oce
51	USE phycst
52	USE in_out_manager
53	#if defined key_zdfddm
54	USE zdfddm ! vertical physics: double diffusion
55	#endif
56	USE eosbn2
57	USE gridrandom
58	USE dotprodfld
59	USE tstool_tam
60	USE wrk_nemo
61	USE timing
62	USE lib_mpp
63	IMPLICIT NONE
64	PRIVATE
65	!
66	!! * Interface
67	INTERFACE eos_tan
68	MODULE PROCEDURE eos_insitu_tan, eos_insitu_pot_tan, eos_insitu_2d_tan, &
69	eos_alpbet_tan
70	END INTERFACE
71	INTERFACE eos_adj
72	MODULE PROCEDURE eos_insitu_adj, eos_insitu_pot_adj, eos_insitu_2d_adj, &
73	eos_alpbet_adj
74	END INTERFACE
75	INTERFACE bn2_tan
76	MODULE PROCEDURE eos_bn2_tan
77	END INTERFACE
78	INTERFACE bn2_adj
79	MODULE PROCEDURE eos_bn2_adj
80	END INTERFACE
81	!
82	!! * Routine accessibility
83	PUBLIC eos_tan ! called by step.F90, inidtr.F90, tranpc.F90 and intgrd.F90
84	PUBLIC bn2_tan ! called by step.F90
85	PUBLIC eos_adj ! called by step.F90, inidtr.F90, tranpc.F90 and intgrd.F90
86	PUBLIC bn2_adj ! called by step.F90
87	#if defined key_tam
88	PUBLIC eos_adj_tst
89	PUBLIC bn2_adj_tst
90	#endif
91	!! * Substitutions
92	# include "domzgr_substitute.h90"
93	# include "vectopt_loop_substitute.h90"
94	!!----------------------------------------------------------------------
95	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
96	!! $Id$
97	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
98	!!----------------------------------------------------------------------
99	CONTAINS
100	SUBROUTINE eos_insitu_tan( pts, pts_tl, prd_tl )
101	!!-----------------------------------------------------------------------
102	!!
103	!! * ROUTINE eos_insitu_tan : TL OF ROUTINE eos_insitu *
104	!!
105	!! ** Purpose of direct routine : Compute the in situ density
106	!! (ratio rho/rau0) from potential temperature and salinity
107	!! using an equation of state defined through the namelist
108	!! parameter nn_eos.
109	!!
110	!! ** Method of direct routine : 3 cases:
111	!! nn_eos = 0 : Jackett and McDougall (1994) equation of state.
112	!! the in situ density is computed directly as a function of
113	!! potential temperature relative to the surface (the opa t
114	!! variable), salt and pressure (assuming no pressure variation
115	!! along geopotential surfaces, i.e. the pressure p in decibars
116	!! is approximated by the depth in meters.
117	!! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0
118	!! with pressure p decibars
119	!! potential temperature t deg celsius
120	!! salinity s psu
121	!! reference volumic mass rau0 kg/m**3
122	!! in situ volumic mass rho kg/m**3
123	!! in situ density anomalie prd no units
124	!! Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar,
125	!! t = 40 deg celcius, s=40 psu
126	!! nn_eos = 1 : linear equation of state function of temperature only
127	!! prd(t) = 0.0285 - rn_alpha * t
128	!! nn_eos = 2 : linear equation of state function of temperature and
129	!! salinity
130	!! prd(t,s) = rn_beta * s - rn_alpha * tn - 1.
131	!! Note that no boundary condition problem occurs in this routine
132	!! as (ptem,psal) are defined over the whole domain.
133	!!
134	!! ** Comments on Adjoint Routine :
135	!! Care has been taken to avoid division by zero when computing
136	!! the inverse of the square root of salinity at masked salinity
137	!! points.
138	!!
139	!! * Arguments
140	REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pts, & ! 1 : potential temperature [Celcius]
141	! ! 2 : salinity [psu]
142	& pts_tl ! 1 : TL of potential temperature [Celsius]
143	! 2 : TL of salinity [psu]
144	REAL(wp), DIMENSION(:,:,:), INTENT( out ) :: &
145	& prd_tl ! TL of potential density (surface referenced)
146	!! * Local declarations
147	INTEGER :: ji, jj, jk ! dummy loop indices
148	REAL(wp) :: & ! temporary scalars
149	zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw, &
150	zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, &
151	zttl, zstl, zhtl, zsrtl, zr1tl, zr2tl, zr3tl, &
152	zr4tl, zrhoptl, zetl, zbwtl, &
153	zbtl, zdtl, zctl, zawtl, zatl, zb1tl, za1tl, &
154	zkwtl, zk0tl, zpes, zrdc1, zrdc2, zeps, &
155	zmask, zrau0r
156	REAL(wp), POINTER, DIMENSION(:,:,:) :: zws
157	!!----------------------------------------------------------------------
158	!
159	IF( nn_timing == 1 ) CALL timing_start('eos_tan')
160	!
161	CALL wrk_alloc( jpi, jpj, jpk, zws )
162	!
163	SELECT CASE ( nn_eos )
164
165	CASE ( 0 ) !== Jackett and McDougall (1994) formulation ==!
166	zrau0r = 1._wp / rau0
167	#ifdef key_sp
168	zeps = 1.e-7_wp
169	#else
170	zeps = 1.e-14_wp
171	#endif
172	!CDIR NOVERRCHK
173	zws(:,:,:) = SQRT( ABS( pts(:,:,:, jp_sal) ) )
174	!
175	DO jk = 1, jpkm1
176	DO jj = 1, jpj
177	DO ji = 1, jpi
178	zt = pts(ji,jj,jk,jp_tem)
179	zs = pts(ji,jj,jk,jp_sal)
180	zh = fsdept(ji,jj,jk) ! depth
181	zsr= zws(ji,jj,jk) ! square root salinity
182	! compute volumic mass pure water at atm pressure
183	zr1= ( ( ( ( 6.536332e-9_wpzt-1.120083e-6_wp )zt+1.001685e-4_wp)*zt &
184	-9.095290e-3_wp )zt+6.793952e-2_wp )zt+999.842594_wp
185	! seawater volumic mass atm pressure
186	zr2= ( ( ( 5.3875e-9_wpzt-8.2467e-7_wp ) zt+7.6438e-5_wp ) *zt &
187	-4.0899e-3_wp ) *zt+0.824493_wp
188	zr3= ( -1.6546e-6_wpzt+1.0227e-4_wp ) zt-5.72466e-3_wp
189	zr4= 4.8314e-4_wp
190	!
191	! potential volumic mass (reference to the surface)
192	zrhop= ( zr4zs + zr3zsr + zr2 ) *zs + zr1
193	!
194	! add the compression terms
195	ze = ( -3.508914e-8_wpzt-1.248266e-8_wp ) zt-2.595994e-6_wp
196	zbw= ( 1.296821e-6_wpzt-5.782165e-9_wp ) zt+1.045941e-4_wp
197	zb = zbw + ze * zs
198
199	zd = -2.042967e-2_wp
200	zc = (-7.267926e-5_wpzt+2.598241e-3_wp ) zt+0.1571896_wp
201	zaw= ( ( 5.939910e-6_wpzt+2.512549e-3_wp ) zt-0.1028859_wp ) *zt - 4.721788_wp
202	za = ( zdzsr + zc ) zs + zaw
203
204	zb1= (-0.1909078_wpzt+7.390729_wp ) zt-55.87545_wp
205	za1= ( ( 2.326469e-3_wpzt+1.553190_wp)zt-65.00517_wp ) *zt+1044.077_wp
206	zkw= ( ( (-1.361629e-4_wpzt-1.852732e-2_wp ) zt-30.41638_wp ) zt + 2098.925_wp ) zt+190925.6_wp
207	zk0= ( zb1zsr + za1 )zs + zkw
208
209	! Tangent linear part
210
211	zttl = pts_tl(ji,jj,jk, jp_tem)
212	zstl = pts_tl(ji,jj,jk, jp_sal)
213
214	zsrtl= ( 1.0 / MAX( 2.*zsr, zeps ) ) &
215	& * tmask(ji,jj,jk) * zstl
216
217	zr1tl= ( ( ( ( 5.6.536332e-9_wp zt &
218	& -4.1.120083e-6_wp ) zt &
219	& +3.1.001685e-4_wp ) zt &
220	& -2.9.095290e-3_wp ) zt &
221	& + 6.793952e-2_wp ) * zttl
222
223	zr2tl= ( ( ( 4.5.3875e-9_wp zt &
224	& -3.8.2467e-7_wp ) zt &
225	& +2.7.6438e-5_wp ) zt &
226	& - 4.0899e-3_wp ) * zttl
227
228	zr3tl= ( -2.1.6546e-6_wp zt &
229	& + 1.0227e-4_wp ) * zttl
230
231	zrhoptl= zr1tl &
232	& + zs * zr2tl &
233	& + zsr * zs * zr3tl &
234	& + zr3 * zs * zsrtl &
235	& + ( 2. * zr4 * zs + zr2 &
236	& + zr3 * zsr ) * zstl
237
238	zetl = ( -2.3.508914e-8_wp zt &
239	& - 1.248266e-8_wp ) * zttl
240
241	zbwtl= ( 2.1.296821e-6_wp zt &
242	& - 5.782165e-9_wp ) * zttl
243
244	zbtl = zbwtl &
245	& + zs * zetl &
246	& + ze * zstl
247
248	zctl = ( -2.7.267926e-5_wp zt &
249	& + 2.598241e-3_wp ) * zttl
250
251	zawtl= ( ( 3.5.939910e-6_wp zt &
252	& +2.2.512549e-3_wp ) zt &
253	& - 0.1028859_wp ) * zttl
254
255	zatl = zawtl &
256	& + zd * zs * zsrtl &
257	& + zs * zctl &
258	& + ( zd * zsr + zc ) * zstl
259
260	zb1tl= ( -2.0.1909078_wp zt &
261	& + 7.390729_wp ) * zttl
262
263	za1tl= ( ( 3.2.326469e-3_wp zt &
264	& +2.1.553190_wp ) zt &
265	& - 65.00517_wp ) * zttl
266
267	zkwtl= ( ( ( -4.1.361629e-4_wp zt &
268	& -3.1.852732e-2_wp ) zt &
269	& -2.30.41638_wp ) zt &
270	& + 2098.925_wp ) * zttl
271
272	zk0tl= zkwtl &
273	& + zb1 * zs * zsrtl &
274	& + zs * zsr * zb1tl &
275	& + zs * za1tl &
276	& + ( zb1 * zsr + za1 ) * zstl
277
278	! Masked in situ density anomaly
279
280	zrdc1 = 1.0 / ( zk0 - zh * ( za - zh * zb ) )
281	zrdc2 = 1.0 / ( 1.0 - zh * zrdc1 )
282
283	prd_tl(ji,jj,jk) = tmask(ji,jj,jk) * zrdc2 * &
284	& ( zrhoptl &
285	& - zrdc2 * zh * zrdc1*2 zrhop &
286	& * ( zk0tl &
287	& - zh * ( zatl &
288	& - zh * zbtl ) ) )&
289	& * zrau0r
290	END DO
291	END DO
292	END DO
293	!
294	CASE ( 1 ) !== Linear formulation function of temperature only ==!
295	DO jk = 1, jpkm1
296	prd_tl(:,:,jk) = ( - rn_alpha * pts_tl(:,:,jk,jp_tem) ) * tmask(:,:,jk)
297	END DO
298	!
299	CASE ( 2 ) !== Linear formulation function of temperature and salinity ==!
300
301	! ! ===============
302	DO jk = 1, jpkm1
303	prd_tl(:,:,jk) = ( rn_beta * pts_tl(:,:,jk,jp_sal) - rn_alpha * pts_tl(:,:,jk,jp_tem ) ) * tmask(:,:,jk)
304	END DO
305	!
306	END SELECT
307	!
308	CALL wrk_dealloc( jpi, jpj, jpk, zws )
309	!
310	IF( nn_timing == 1 ) CALL timing_stop('eos_tan')
311	!
312	END SUBROUTINE eos_insitu_tan
313
314	SUBROUTINE eos_insitu_pot_tan( pts, pts_tl, prd_tl, prhop_tl)
315	!!----------------------------------------------------------------------
316	!! * ROUTINE eos_insitu_pot_tan *
317	!!
318	!! ** Purpose or the direct routine:
319	!! Compute the in situ density (ratio rho/rau0) and the
320	!! potential volumic mass (Kg/m3) from potential temperature and
321	!! salinity fields using an equation of state defined through the
322	!! namelist parameter nn_eos.
323	!!
324	!! ** Method :
325	!! nn_eos = 0 : Jackett and McDougall (1994) equation of state.
326	!! the in situ density is computed directly as a function of
327	!! potential temperature relative to the surface (the opa t
328	!! variable), salt and pressure (assuming no pressure variation
329	!! along geopotential surfaces, i.e. the pressure p in decibars
330	!! is approximated by the depth in meters.
331	!! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0
332	!! rhop(t,s) = rho(t,s,0)
333	!! with pressure p decibars
334	!! potential temperature t deg celsius
335	!! salinity s psu
336	!! reference volumic mass rau0 kg/m**3
337	!! in situ volumic mass rho kg/m**3
338	!! in situ density anomalie prd no units
339	!!
340	!! Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar,
341	!! t = 40 deg celcius, s=40 psu
342	!!
343	!! nn_eos = 1 : linear equation of state function of temperature only
344	!! prd(t) = ( rho(t) - rau0 ) / rau0 = 0.028 - rn_alpha * t
345	!! rhop(t,s) = rho(t,s)
346	!!
347	!! nn_eos = 2 : linear equation of state function of temperature and
348	!! salinity
349	!! prd(t,s) = ( rho(t,s) - rau0 ) / rau0
350	!! = rn_beta * s - rn_alpha * tn - 1.
351	!! rhop(t,s) = rho(t,s)
352	!! Note that no boundary condition problem occurs in this routine
353	!! as (tn,sn) or (ta,sa) are defined over the whole domain.
354	!!
355	!! ** Action : - prd , the in situ density (no units)
356	!! - prhop, the potential volumic mass (Kg/m3)
357	!!
358	!! References :
359	!! Jackett, D.R., and T.J. McDougall. J. Atmos. Ocean. Tech., 1994
360	!! Brown, J. A. and K. A. Campana. Mon. Weather Rev., 1978
361	!!
362	!!----------------------------------------------------------------------
363	!! * Arguments
364	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts, & ! 1 : potential temperature [Celcius]
365	! ! 2 : salinity [psu]
366	& pts_tl ! 1 : TL of potential temperature [Celcius]
367	! ! 2 : TL of salinity [psu]
368	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT( out) :: prd_tl ! TL of in_situ density [-]
369	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT( out) :: prhop_tl ! TL of potential density (surface referenced)
370	!! * Local declarations
371	INTEGER :: ji, jj, jk ! dummy loop indices
372	REAL(wp) :: & ! temporary scalars
373	zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw, &
374	zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, &
375	zttl, zstl, zhtl, zsrtl, zr1tl, zr2tl, zr3tl, &
376	zr4tl, zrhoptl, zetl, zbwtl, &
377	zbtl, zdtl, zctl, zawtl, zatl, zb1tl, za1tl, &
378	zkwtl, zk0tl, zpes, zrdc1, zrdc2, zeps, &
379	zmask, zrau0r
380	REAL(wp), POINTER, DIMENSION(:,:,:) :: zws
381	!!----------------------------------------------------------------------
382	!
383	IF( nn_timing == 1 ) CALL timing_start('eos-p_tan')
384	!
385	CALL wrk_alloc( jpi, jpj, jpk, zws )
386	!
387	SELECT CASE ( nn_eos )
388	!
389	CASE ( 0 ) !== Jackett and McDougall (1994) formulation ==!
390	zrau0r = 1.e0 / rau0
391	#ifdef key_sp
392	zeps = 1.e-7
393	#else
394	zeps = 1.e-14
395	#endif
396	!CDIR NOVERRCHK
397	zws(:,:,:) = SQRT( ABS( pts(:,:,:, jp_sal) ) )
398	!
399	DO jk = 1, jpkm1
400	DO jj = 1, jpj
401	DO ji = 1, jpi
402	zt = pts(ji,jj,jk, jp_tem)
403	zs = pts(ji,jj,jk, jp_sal)
404	zh = fsdept(ji,jj,jk) ! depth
405	zsr = zws(ji,jj,jk) ! square root salinity
406	! compute volumic mass pure water at atm pressure
407	zr1 = ( ( ( ( 6.536332e-9_wp * zt - 1.120083e-6_wp ) * zt &
408	& + 1.001685e-4_wp ) * zt - 9.095290e-3_wp ) * zt &
409	& + 6.793952e-2_wp ) * zt + 999.842594_wp
410	! seawater volumic mass atm pressure
411	zr2 = ( ( ( 5.3875e-9_wp * zt - 8.2467e-7_wp ) * zt &
412	& + 7.6438e-5_wp ) * zt - 4.0899e-3_wp ) * zt &
413	& + 0.824493_wp
414	zr3 = ( -1.6546e-6_wp * zt + 1.0227e-4_wp ) * zt - 5.72466e-3_wp
415	zr4 = 4.8314e-4_wp
416
417	! potential volumic mass (reference to the surface)
418	zrhop= ( zr4 * zs + zr3 * zsr + zr2 ) * zs + zr1
419
420	! add the compression terms
421	ze = ( -3.508914e-8_wp * zt - 1.248266e-8_wp ) * zt - 2.595994e-6_wp
422	zbw = ( 1.296821e-6_wp * zt - 5.782165e-9_wp ) * zt + 1.045941e-4_wp
423	zb = zbw + ze * zs
424
425	zd = -2.042967e-2_wp
426	zc = (-7.267926e-5_wp * zt + 2.598241e-3_wp ) * zt + 0.1571896_wp
427	zaw= ( ( 5.939910e-6_wp * zt + 2.512549e-3_wp ) * zt - 0.1028859_wp ) * zt - 4.721788_wp
428	za = ( zd * zsr + zc ) * zs + zaw
429
430	zb1 = (-0.1909078_wp * zt + 7.390729_wp ) * zt - 55.87545_wp
431	za1 = ( ( 2.326469e-3_wp * zt + 1.553190_wp ) * zt - 65.00517_wp &
432	& ) * zt + 1044.077_wp
433	zkw = ( ( (-1.361629e-4_wp * zt - 1.852732e-2_wp ) * zt - 30.41638_wp &
434	& ) * zt + 2098.925_wp ) * zt + 190925.6_wp
435	zk0 = ( zb1 * zsr + za1 ) * zs + zkw
436
437
438	zrdc1 = 1.0 / ( zk0 - zh * ( za - zh * zb ) )
439	zrdc2 = 1.0 / ( 1.0 - zh * zrdc1 )
440
441	! Tangent linear part
442
443	zttl = pts_tl(ji,jj,jk, jp_tem)
444	zstl = pts_tl(ji,jj,jk, jp_sal)
445
446	zsrtl= ( 1.0 / MAX( 2.*zsr, zeps ) ) &
447	& * tmask(ji,jj,jk) * zstl
448
449	zr1tl= ( ( ( ( 5.6.536332e-9_wp zt &
450	& -4.1.120083e-6_wp ) zt &
451	& +3.1.001685e-4_wp ) zt &
452	& -2.9.095290e-3_wp ) zt &
453	& + 6.793952e-2_wp ) * zttl
454
455	zr2tl= ( ( ( 4.5.3875e-9_wp zt &
456	& -3.8.2467e-7_wp ) zt &
457	& +2.7.6438e-5_wp ) zt &
458	& - 4.0899e-3_wp ) * zttl
459
460	zr3tl= ( -2.1.6546e-6_wp zt &
461	& + 1.0227e-4_wp ) * zttl
462
463	zrhoptl= zr1tl &
464	& + zs * zr2tl &
465	& + zsr * zs * zr3tl &
466	& + zr3 * zs * zsrtl &
467	& + ( 2. * zr4 * zs + zr2 &
468	& + zr3 * zsr ) * zstl
469
470	prhop_tl(ji,jj,jk) = zrhoptl * tmask(ji,jj,jk)
471
472	zetl = ( -2.3.508914e-8_wp zt &
473	& - 1.248266e-8_wp ) * zttl
474
475	zbwtl= ( 2.1.296821e-6_wp zt &
476	& - 5.782165e-9_wp ) * zttl
477
478	zbtl = zbwtl &
479	& + zs * zetl &
480	& + ze * zstl
481
482	zctl = ( -2.7.267926e-5_wp zt &
483	& + 2.598241e-3_wp ) * zttl
484
485	zawtl= ( ( 3.5.939910e-6_wp zt &
486	& +2.2.512549e-3_wp ) zt &
487	& - 0.1028859_wp ) * zttl
488
489	zatl = zawtl &
490	& + zd * zs * zsrtl &
491	& + zs * zctl &
492	& + ( zd * zsr + zc ) * zstl
493
494	zb1tl= ( -2.0.1909078_wp zt &
495	& + 7.390729_wp ) * zttl
496
497	za1tl= ( ( 3.2.326469e-3_wp zt &
498	& +2.1.553190_wp ) zt &
499	& - 65.00517_wp ) * zttl
500
501	zkwtl= ( ( ( -4.1.361629e-4_wp zt &
502	& -3.1.852732e-2_wp ) zt &
503	& -2.30.41638_wp ) zt &
504	& + 2098.925_wp ) * zttl
505
506	zk0tl= zkwtl &
507	& + zb1 * zs * zsrtl &
508	& + zs * zsr * zb1tl &
509	& + zs * za1tl &
510	& + ( zb1 * zsr + za1 ) * zstl
511
512	! Masked in situ density anomaly
513
514	prd_tl(ji,jj,jk) = tmask(ji,jj,jk) * zrdc2 * &
515	& ( zrhoptl &
516	& - zrdc2 * zh * zrdc1*2 zrhop &
517	& * ( zk0tl &
518	& - zh * ( zatl &
519	& - zh * zbtl ) ) )&
520	& * zrau0r
521	END DO
522	END DO
523	END DO
524	!
525	CASE ( 1 ) !== Linear formulation = F( temperature ) ==!
526	DO jk = 1, jpkm1
527	prd_tl (:,:,jk) = ( - rn_alpha * pts_tl(:,:,jk, jp_tem) ) * tmask(:,:,jk)
528	prhop_tl(:,:,jk) = ( rau0 * prd_tl(:,:,jk) ) * tmask(:,:,jk)
529	END DO
530	!
531	CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==!
532	DO jk = 1, jpkm1
533	prd_tl(:,:,jk) = ( rn_beta * pts_tl(:,:,jk, jp_sal) - rn_alpha * pts_tl(:,:,jk, jp_tem) ) * tmask(:,:,jk)
534	prhop_tl(:,:,jk) = ( rau0 * prd_tl(:,:,jk) ) * tmask(:,:,jk)
535	END DO
536	!
537	END SELECT
538	!
539	CALL wrk_dealloc( jpi, jpj, jpk, zws )
540	!
541	IF( nn_timing == 1 ) CALL timing_stop('eos-p_tan')
542	!
543	END SUBROUTINE eos_insitu_pot_tan
544	SUBROUTINE eos_insitu_2d_tan( pts, pdep, pts_tl, prd_tl )
545	!!-----------------------------------------------------------------------
546	!!
547	!! * ROUTINE eos_insitu_2d_tan : TL OF ROUTINE eos_insitu_2d *
548	!!
549	!! ** Purpose of direct routine : Compute the in situ density
550	!! (ratio rho/rau0) from potential temperature and salinity
551	!! using an equation of state defined through the namelist
552	!! parameter nn_eos. * 2D field case
553	!!
554	!! ** Method of direct routine : 3 cases:
555	!! nn_eos = 0 : Jackett and McDougall (1994) equation of state.
556	!! the in situ density is computed directly as a function of
557	!! potential temperature relative to the surface (the opa t
558	!! variable), salt and pressure (assuming no pressure variation
559	!! along geopotential surfaces, i.e. the pressure p in decibars
560	!! is approximated by the depth in meters.
561	!! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0
562	!! with pressure p decibars
563	!! potential temperature t deg celsius
564	!! salinity s psu
565	!! reference volumic mass rau0 kg/m**3
566	!! in situ volumic mass rho kg/m**3
567	!! in situ density anomalie prd no units
568	!! Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar,
569	!! t = 40 deg celcius, s=40 psu
570	!! nn_eos = 1 : linear equation of state function of temperature only
571	!! prd(t) = 0.0285 - ralpha * t
572	!! nn_eos = 2 : linear equation of state function of temperature and
573	!! salinity
574	!! prd(t,s) = rn_beta * s - rn_alpha * tn - 1.
575	!! Note that no boundary condition problem occurs in this routine
576	!! as (ptem,psal) are defined over the whole domain.
577	!!
578	!! ** Comments on Adjoint Routine :
579	!! Care has been taken to avoid division by zero when computing
580	!! the inverse of the square root of salinity at masked salinity
581	!! points.
582	!!
583	!! ** Action :
584	!!
585	!! References :
586	!!
587	!! History :
588	!! 8.2 ! 05-03 ((F. Van den Berghe, A. Weaver, N. Daget) - eostan.F
589	!! 9.0 ! 07-07 (K. Mogensen) Initial version based on eostan.F
590	!! ! 08-07 (A. Vidard) bug fix in computation of prd_tl if neos=1
591	!!-----------------------------------------------------------------------
592	!! * Modules used
593	!! * Arguments
594	REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in ) :: pts, & ! 1 : potential temperature [Celcius]
595	! ! 2 : salinity [psu]
596	& pts_tl ! 1 : TL of potential temperature [Celcius]
597	! ! 2 : TL of salinity [psu]
598	REAL(wp), DIMENSION(jpi,jpj ), INTENT( out) :: prd_tl ! TL of in_situ density [-]
599	REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pdep ! depth [m]
600	!
601	!! * Local declarations
602	INTEGER :: ji, jj, jk ! dummy loop indices
603	REAL(wp) :: & ! temporary scalars
604	zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw, &
605	zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, &
606	zttl, zstl, zhtl, zsrtl, zr1tl, zr2tl, zr3tl, &
607	zr4tl, zrhoptl, zetl, zbwtl, &
608	zbtl, zdtl, zctl, zawtl, zatl, zb1tl, za1tl, &
609	zkwtl, zk0tl, zpes, zrdc1, zrdc2, zeps, &
610	zmask
611	REAL(wp), POINTER, DIMENSION(:,:) :: zws
612	!!----------------------------------------------------------------------
613	!
614	IF( nn_timing == 1 ) CALL timing_start('eos2d')
615	!
616	CALL wrk_alloc( jpi, jpj, zws )
617	!
618	SELECT CASE ( nn_eos )
619	!
620	CASE ( 0 ) !== Jackett and McDougall (1994) formulation ==!
621	!
622	#ifdef key_sp
623	zeps = 1.e-7
624	#else
625	zeps = 1.e-14
626	#endif
627	!CDIR NOVERRCHK
628	DO jj = 1, jpjm1
629	!CDIR NOVERRCHK
630	DO ji = 1, fs_jpim1 ! vector opt.
631	zws(ji,jj) = SQRT( ABS( pts(ji,jj, jp_sal) ) )
632	END DO
633	END DO
634	DO jj = 1, jpjm1
635	DO ji = 1, fs_jpim1 ! vector opt.
636
637	zmask = tmask(ji,jj,1) ! land/sea bottom mask = surf. mask
638
639	zt = pts (ji,jj, jp_tem) ! interpolated T
640	zs = pts (ji,jj, jp_sal) ! interpolated S
641	zsr= zws(ji,jj) ! square root of interpolated S
642	zh = pdep(ji,jj) ! depth at the partial step level
643	! compute volumic mass pure water at atm pressure
644	zr1= ( ( ( ( 6.536332e-9_wpzt-1.120083e-6_wp )zt+1.001685e-4_wp)*zt &
645	& -9.095290e-3_wp )zt+6.793952e-2_wp )zt+999.842594_wp
646	! seawater volumic mass atm pressure
647	zr2= ( ( ( 5.3875e-9_wpzt-8.2467e-7_wp ) zt+7.6438e-5_wp ) *zt &
648	& -4.0899e-3_wp ) *zt+0.824493_wp
649	zr3= ( -1.6546e-6_wpzt+1.0227e-4_wp ) zt-5.72466e-3_wp
650	zr4= 4.8314e-4_wp
651
652	! potential volumic mass (reference to the surface)
653	zrhop= ( zr4zs + zr3zsr + zr2 ) *zs + zr1
654
655	! add the compression terms
656	ze = ( -3.508914e-8_wpzt-1.248266e-8_wp ) zt-2.595994e-6_wp
657	zbw= ( 1.296821e-6zt-5.782165e-9_wp ) zt+1.045941e-4_wp
658	zb = zbw + ze * zs
659
660	zd = -2.042967e-2_wp
661	zc = (-7.267926e-5_wpzt+2.598241e-3_wp ) zt+0.1571896_wp
662	zaw= ( ( 5.939910e-6_wpzt+2.512549e-3_wp ) zt-0.1028859_wp ) *zt - 4.721788_wp
663	za = ( zdzsr + zc ) zs + zaw
664
665	zb1= (-0.1909078_wpzt+7.390729_wp ) zt-55.87545_wp
666	za1= ( ( 2.326469e-3_wpzt+1.553190_wp)zt-65.00517_wp ) *zt+1044.077_wp
667	zkw= ( ( (-1.361629e-4_wpzt-1.852732e-2_wp ) zt-30.41638_wp ) zt + 2098.925_wp ) zt+190925.6_wp
668	zk0= ( zb1zsr + za1 )zs + zkw
669
670	! Tangent linear part
671
672	zttl = pts_tl(ji,jj, jp_tem)
673	zstl = pts_tl(ji,jj, jp_sal)
674
675	zsrtl= ( 1.0 / MAX( 2.*zsr, zeps ) ) &
676	& * tmask(ji,jj,1) * zstl
677
678	zr1tl= ( ( ( ( 5.6.536332e-9_wp zt &
679	& -4.1.120083e-6_wp ) zt &
680	& +3.1.001685e-4_wp ) zt &
681	& -2.9.095290e-3_wp ) zt &
682	& + 6.793952e-2_wp ) * zttl
683
684	zr2tl= ( ( ( 4.5.3875e-9_wp zt &
685	& -3.8.2467e-7_wp ) zt &
686	& +2.7.6438e-5_wp ) zt &
687	& - 4.0899e-3_wp ) * zttl
688
689	zr3tl= ( -2.1.6546e-6_wp zt &
690	& + 1.0227e-4_wp ) * zttl
691
692	zrhoptl= zr1tl &
693	& + zs * zr2tl &
694	& + zsr * zs * zr3tl &
695	& + zr3 * zs * zsrtl &
696	& + ( 2. * zr4 * zs + zr2 &
697	& + zr3 * zsr ) * zstl
698
699	zetl = ( -2.3.508914e-8_wp zt &
700	& - 1.248266e-8_wp ) * zttl
701
702	zbwtl= ( 2.1.296821e-6_wp zt &
703	& - 5.782165e-9_wp ) * zttl
704
705	zbtl = zbwtl &
706	& + zs * zetl &
707	& + ze * zstl
708
709	zctl = ( -2.7.267926e-5_wp zt &
710	& + 2.598241e-3_wp ) * zttl
711
712	zawtl= ( ( 3.5.939910e-6_wp zt &
713	& +2.2.512549e-3_wp ) zt &
714	& - 0.1028859_wp ) * zttl
715
716	zatl = zawtl &
717	& + zd * zs * zsrtl &
718	& + zs * zctl &
719	& + ( zd * zsr + zc ) * zstl
720
721	zb1tl= ( -2.0.1909078_wp zt &
722	& + 7.390729_wp ) * zttl
723
724	za1tl= ( ( 3.2.326469e-3_wp zt &
725	& +2.1.553190_wp ) zt &
726	& - 65.00517_wp ) * zttl
727
728	zkwtl= ( ( ( -4.1.361629e-4_wp zt &
729	& -3.1.852732e-2_wp ) zt &
730	& -2.30.41638_wp ) zt &
731	& + 2098.925_wp ) * zttl
732
733	zk0tl= zkwtl &
734	& + zb1 * zs * zsrtl &
735	& + zs * zsr * zb1tl &
736	& + zs * za1tl &
737	& + ( zb1 * zsr + za1 ) * zstl
738
739	! Masked in situ density anomaly
740
741	zrdc1 = 1.0 / ( zk0 - zh * ( za - zh * zb ) )
742	zrdc2 = 1.0 / ( 1.0 - zh * zrdc1 )
743
744	prd_tl(ji,jj) = tmask(ji,jj,1) * zrdc2 * &
745	& ( zrhoptl &
746	& - zrdc2 * zh * zrdc1*2 zrhop &
747	& * ( zk0tl &
748	& - zh * ( zatl &
749	& - zh * zbtl ) ) )&
750	& / rau0
751
752
753	END DO
754	END DO
755	!
756	CASE ( 1 ) !== Linear formulation = F( temperature ) ==!
757	DO jj = 1, jpjm1
758	DO ji = 1, fs_jpim1 ! vector opt.
759	prd_tl(ji,jj) = ( - rn_alpha * pts_tl(ji,jj,jp_tem) ) * tmask(ji,jj,1)
760	END DO
761	END DO
762	!
763	CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==!
764	DO jj = 1, jpjm1
765	DO ji = 1, fs_jpim1 ! vector opt.
766	prd_tl (ji,jj) = ( rn_beta * pts_tl(ji,jj, jp_sal) - rn_alpha * pts_tl(ji,jj, jp_tem) ) * tmask(ji,jj,1)
767	END DO
768	END DO
769	!
770	END SELECT
771	!
772	END SUBROUTINE eos_insitu_2d_tan
773
774	SUBROUTINE eos_insitu_adj(pts, pts_ad, prd_ad)
775	!!-----------------------------------------------------------------------
776	!!
777	!! * ROUTINE eos_insitu_tan : Adjoint OF ROUTINE eos_insitu *
778	!!
779	!! ** Purpose of direct routine : Compute the in situ density
780	!! (ratio rho/rau0) from potential temperature and salinity
781	!! using an equation of state defined through the namelist
782	!! parameter nneos.
783	!!
784	!! ** Method of direct routine : 3 cases:
785	!! nn_eos = 0 : Jackett and McDougall (1994) equation of state.
786	!! the in situ density is computed directly as a function of
787	!! potential temperature relative to the surface (the opa t
788	!! variable), salt and pressure (assuming no pressure variation
789	!! along geopotential surfaces, i.e. the pressure p in decibars
790	!! is approximated by the depth in meters.
791	!! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0
792	!! with pressure p decibars
793	!! potential temperature t deg celsius
794	!! salinity s psu
795	!! reference volumic mass rau0 kg/m**3
796	!! in situ volumic mass rho kg/m**3
797	!! in situ density anomalie prd no units
798	!! Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar,
799	!! t = 40 deg celcius, s=40 psu
800	!! nn_eos = 1 : linear equation of state function of temperature only
801	!! prd(t) = 0.0285 - rn_alpha * t
802	!! nn_eos = 2 : linear equation of state function of temperature and
803	!! salinity
804	!! prd(t,s) = rn_beta * s - rn_alpha * tn - 1.
805	!! Note that no boundary condition problem occurs in this routine
806	!! as (ptem,psal) are defined over the whole domain.
807	!!
808	!! ** Comments on Adjoint Routine :
809	!! Care has been taken to avoid division by zero when computing
810	!! the inverse of the square root of salinity at masked salinity
811	!! points.
812	!!
813	!! ** Action :
814	!!
815	!! References :
816	!!
817	!! History :
818	!! 8.2 ! 05-03 ((F. Van den Berghe, A. Weaver, N. Daget) - eostan.F
819	!! 9.0 ! 08-08 (A. Vidard) 9.0 version
820	!!-----------------------------------------------------------------------
821	!! * Modules used
822	!! * Arguments
823	REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pts ! 1 : potential temperature [Celcius]
824	! ! 2 : salinity [psu]
825	REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pts_ad ! 1 : TL of potential temperature [Celsius]
826	! 2 : TL of salinity [psu]
827	REAL(wp), DIMENSION(:,:,:), INTENT( inout ) :: &
828	& prd_ad ! TL of potential density (surface referenced)
829
830	!! * Local declarations
831	INTEGER :: ji, jj, jk ! dummy loop indices
832	REAL(wp) :: & ! temporary scalars
833	zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw, &
834	zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, &
835	ztad, zsad, zhad, zsrad, zr1ad, zr2ad, zr3ad, &
836	zr4ad, zrhopad, zead, zbwad, &
837	zbad, zdad, zcad, zawad, zaad, zb1ad, za1ad, &
838	zkwad, zk0ad, zpes, zrdc1, zrdc2, zeps, &
839	zmask, zrau0r
840	REAL(wp), POINTER, DIMENSION(:,:,:) :: zws
841	!!----------------------------------------------------------------------
842	IF( nn_timing == 1 ) CALL timing_start('eos_adj')
843	!
844	CALL wrk_alloc( jpi, jpj, jpk, zws )
845	!
846	! initialization of adjoint variables
847	ztad = 0.0_wp
848	zsad = 0.0_wp
849	zhad = 0.0_wp
850	zsrad = 0.0_wp
851	zr1ad = 0.0_wp
852	zr2ad = 0.0_wp
853	zr3ad = 0.0_wp
854	zr4ad = 0.0_wp
855	zrhopad = 0.0_wp
856	zead = 0.0_wp
857	zbwad = 0.0_wp
858	zbad = 0.0_wp
859	zdad = 0.0_wp
860	zcad = 0.0_wp
861	zawad = 0.0_wp
862	zaad = 0.0_wp
863	zb1ad = 0.0_wp
864	za1ad = 0.0_wp
865	zkwad = 0.0_wp
866	zk0ad = 0.0_wp
867	SELECT CASE ( nn_eos )
868
869	CASE ( 0 ) !== Jackett and McDougall (1994) formulation ==!
870	zrau0r = 1.e0 / rau0
871	#ifdef key_sp
872	zeps = 1.e-7
873	#else
874	zeps = 1.e-14
875	#endif
876	!CDIR NOVERRCHK
877	zws(:,:,:) = SQRT( ABS( pts(:,:,:, jp_sal) ) )
878	DO jk = jpkm1, 1, -1
879	DO jj = jpj, 1, -1
880	DO ji = jpi, 1, -1
881	zt = pts(ji,jj,jk, jp_tem)
882	zs = pts(ji,jj,jk, jp_sal)
883	zh = fsdept(ji,jj,jk) ! depth
884	zsr= zws(ji,jj,jk) ! square root salinity
885	! compute volumic mass pure water at atm pressure
886	zr1= ( ( ( ( 6.536332e-9_wpzt-1.120083e-6_wp )zt+1.001685e-4_wp)*zt &
887	-9.095290e-3_wp )zt+6.793952e-2_wp )zt+999.842594_wp
888	! seawater volumic mass atm pressure
889	zr2= ( ( ( 5.3875e-9_wpzt-8.2467e-7_wp ) zt+7.6438e-5_wp ) *zt &
890	-4.0899e-3_wp ) *zt+0.824493_wp
891	zr3= ( -1.6546e-6_wpzt+1.0227e-4_wp ) zt-5.72466e-3_wp
892	zr4= 4.8314e-4_wp
893
894	! potential volumic mass (reference to the surface)
895	zrhop= ( zr4zs + zr3zsr + zr2 ) *zs + zr1
896
897	! add the compression terms
898	ze = ( -3.508914e-8_wpzt-1.248266e-8_wp ) zt-2.595994e-6_wp
899	zbw= ( 1.296821e-6_wpzt-5.782165e-9_wp ) zt+1.045941e-4_wp
900	zb = zbw + ze * zs
901
902	zd = -2.042967e-2_wp
903	zc = (-7.267926e-5_wpzt+2.598241e-3_wp ) zt+0.1571896_wp
904	zaw= ( ( 5.939910e-6_wpzt+2.512549e-3_wp ) zt-0.1028859_wp ) *zt - 4.721788_wp
905	za = ( zdzsr + zc ) zs + zaw
906
907	zb1= (-0.1909078_wpzt+7.390729_wp ) zt-55.87545_wp
908	za1= ( ( 2.326469e-3_wpzt+1.553190_wp)zt-65.00517_wp ) *zt+1044.077_wp
909	zkw= ( ( (-1.361629e-4_wpzt-1.852732e-2_wp ) zt-30.41638_wp ) zt + 2098.925_wp ) zt+190925.6_wp
910	zk0= ( zb1zsr + za1 )zs + zkw
911
912	zrdc1 = 1.0 / ( zk0 - zh * ( za - zh * zb ) )
913	zrdc2 = 1.0 / ( 1.0 - zh * zrdc1 )
914	! ============
915	! Adjoint part
916	! ============
917
918	! Masked in situ density anomaly
919
920	zrhopad = zrhopad + prd_ad(ji,jj,jk) * tmask(ji,jj,jk) &
921	& * zrdc2 * zrau0r
922	zk0ad = zk0ad - prd_ad(ji,jj,jk) * tmask(ji,jj,jk) &
923	& * zrdc2 * zrdc2 * zh &
924	& * zrdc1*2 zrhop &
925	& * zrau0r
926	zaad = zaad + prd_ad(ji,jj,jk) * tmask(ji,jj,jk) &
927	& * zrdc2 * zrdc2 * zh &
928	& * zrdc1*2 zrhop &
929	& * zh * zrau0r
930	zbad = zbad - prd_ad(ji,jj,jk) * tmask(ji,jj,jk) &
931	& * zrdc2 * zrdc2 * zh &
932	& * zrdc1*2 zrhop &
933	& * zh * zh * zrau0r
934	prd_ad(ji,jj,jk) = 0.0_wp
935
936	zkwad = zkwad + zk0ad
937	zsrad = zsrad + zk0ad * zb1 * zs
938	zb1ad = zb1ad + zk0ad * zs * zsr
939	za1ad = za1ad + zk0ad * zs
940	zsad = zsad + zk0ad * ( zb1 * zsr + za1 )
941	zk0ad = 0.0_wp
942
943	ztad = ztad + zkwad * ( ( (-4.1.361629e-4_wp zt &
944	& -3.1.852732e-2_wp ) zt &
945	& -2.30.41638_wp ) zt &
946	& + 2098.925_wp )
947	zkwad = 0.0_wp
948
949	ztad = ztad + za1ad * ( ( 3.2.326469e-3_wp zt &
950	& +2.1.553190_wp ) zt &
951	& - 65.00517_wp )
952	za1ad = 0.0_wp
953
954	ztad = ztad + zb1ad * (-2.0.1909078_wp zt &
955	& + 7.390729_wp )
956	zb1ad = 0.0_wp
957
958	zawad = zawad + zaad
959	zsrad = zsrad + zaad * zd * zs
960	zcad = zcad + zaad * zs
961	zsad = zsad + zaad * ( zd * zsr + zc )
962	zaad = 0.0_wp
963
964	ztad = ztad + zawad * ( ( 3.5.939910e-6_wp zt &
965	& +2.2.512549e-3_wp ) zt &
966	& - 0.1028859_wp )
967	zawad = 0.0_wp
968
969	ztad = ztad + zcad * (-2.7.267926e-5_wp zt &
970	& + 2.598241e-3_wp )
971	zcad = 0.0_wp
972
973	zbwad = zbwad + zbad
974	zead = zead + zbad * zs
975	zsad = zsad + zbad * ze
976	zbad = 0.0_wp
977
978	ztad = ztad + zbwad * ( 2.1.296821e-6_wp zt &
979	& - 5.782165e-9_wp )
980	zbwad = 0.0_wp
981
982	ztad = ztad + zead * (-2.3.508914e-8_wp zt &
983	& - 1.248266e-8_wp )
984	zead = 0.0_wp
985
986	zr1ad = zr1ad + zrhopad
987	zr2ad = zr2ad + zrhopad * zs
988	zr3ad = zr3ad + zrhopad * zsr * zs
989	zsrad = zsrad + zrhopad * zr3 * zs
990	zsad = zsad + zrhopad * ( 2. * zr4 * zs + zr2 &
991	& + zr3 * zsr )
992	zrhopad = 0.0_wp
993
994	ztad = ztad + zr3ad * (-2.1.6546e-6_wp zt &
995	& + 1.0227e-4_wp )
996	zr3ad = 0.0_wp
997
998	ztad = ztad + zr2ad * ( ( ( 4.5.3875e-9_wp zt &
999	& -3.8.2467e-7_wp ) zt &
1000	& +2.7.6438e-5_wp ) zt &
1001	& - 4.0899e-3_wp )
1002	zr2ad = 0.0_wp
1003
1004	ztad = ztad + zr1ad * ( ( ( ( 5.6.536332e-9_wp zt &
1005	& -4.1.120083e-6_wp ) zt &
1006	& +3.1.001685e-4_wp ) zt &
1007	& -2.9.095290e-3_wp ) zt &
1008	& + 6.793952e-2_wp )
1009	zr1ad = 0.0_wp
1010
1011	zsad = zsad + zsrad * ( 1.0 / MAX( 2.*zsr, zeps ) ) &
1012	& * tmask(ji,jj,jk)
1013	zsrad = 0.0_wp
1014
1015	pts_ad(ji,jj,jk, jp_sal) = pts_ad(ji,jj,jk, jp_sal) + zsad
1016	pts_ad(ji,jj,jk,jp_tem) = pts_ad(ji,jj,jk, jp_tem) + ztad
1017	ztad = 0.0_wp
1018	zsad = 0.0_wp
1019	END DO
1020	END DO
1021	END DO
1022	!
1023	CASE ( 1 ) !== Linear formulation function of temperature only ==!
1024	DO jk = jpkm1, 1, -1
1025	pts_ad(:,:,jk,jp_tem) = pts_ad(:,:,jk,jp_tem) - rn_alpha * prd_ad(:,:,jk) * tmask(:,:,jk)
1026	prd_ad(:,:,jk) = 0.0_wp
1027	END DO
1028	!
1029	CASE ( 2 ) !== Linear formulation function of temperature and salinity ==!
1030	DO jk = jpkm1, 1, -1
1031	pts_ad(:,:,jk,jp_tem) = pts_ad(:,:,jk,jp_tem) - rn_alpha * prd_ad(:,:,jk) * tmask(:,:,jk)
1032	pts_ad(:,:,jk,jp_sal) = pts_ad(:,:,jk,jp_sal) + rn_beta * prd_ad( :,:,jk) * tmask(:,:,jk)
1033	prd_ad( :,:,jk) = 0.0_wp
1034	END DO
1035	!
1036	END SELECT
1037	!
1038	CALL wrk_dealloc( jpi, jpj, jpk, zws )
1039	!
1040	IF( nn_timing == 1 ) CALL timing_stop('eos_adj')
1041	!
1042	END SUBROUTINE eos_insitu_adj
1043
1044	SUBROUTINE eos_insitu_pot_adj ( pts, pts_ad, prd_ad, prhop_ad )
1045	!!----------------------------------------------------------------------
1046	!! * ROUTINE eos_insitu_pot_adj *
1047	!!
1048	!! ** Purpose or the direct routine:
1049	!! Compute the in situ density (ratio rho/rau0) and the
1050	!! potential volumic mass (Kg/m3) from potential temperature and
1051	!! salinity fields using an equation of state defined through the
1052	!! namelist parameter nn_eos.
1053	!!
1054	!! ** Method :
1055	!! nn_eos = 0 : Jackett and McDougall (1994) equation of state.
1056	!! the in situ density is computed directly as a function of
1057	!! potential temperature relative to the surface (the opa t
1058	!! variable), salt and pressure (assuming no pressure variation
1059	!! along geopotential surfaces, i.e. the pressure p in decibars
1060	!! is approximated by the depth in meters.
1061	!! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0
1062	!! rhop(t,s) = rho(t,s,0)
1063	!! with pressure p decibars
1064	!! potential temperature t deg celsius
1065	!! salinity s psu
1066	!! reference volumic mass rau0 kg/m**3
1067	!! in situ volumic mass rho kg/m**3
1068	!! in situ density anomalie prd no units
1069	!!
1070	!! Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar,
1071	!! t = 40 deg celcius, s=40 psu
1072	!!
1073	!! neos = 1 : linear equation of state function of temperature only
1074	!! prd(t) = ( rho(t) - rau0 ) / rau0 = 0.028 - ralpha * t
1075	!! rhop(t,s) = rho(t,s)
1076	!!
1077	!! nn_eos = 2 : linear equation of state function of temperature and
1078	!! salinity
1079	!! prd(t,s) = ( rho(t,s) - rau0 ) / rau0
1080	!! = rn_beta * s - rn_alpha * tn - 1.
1081	!! rhop(t,s) = rho(t,s)
1082	!! Note that no boundary condition problem occurs in this routine
1083	!! as (tn,sn) or (ta,sa) are defined over the whole domain.
1084	!!
1085	!! ** Action : - prd , the in situ density (no units)
1086	!! - prhop, the potential volumic mass (Kg/m3)
1087	!!
1088	!! References :
1089	!! Jackett, D.R., and T.J. McDougall. J. Atmos. Ocean. Tech., 1994
1090	!! Brown, J. A. and K. A. Campana. Mon. Weather Rev., 1978
1091	!!
1092	!! History of the adjoint routine:
1093	!! 9.0 ! 08-06 (A. Vidard) Initial version
1094	!!----------------------------------------------------------------------
1095	!! * Arguments
1096
1097	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! 1 : potential temperature/salinity [Celcius/psu]
1098	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(inout) :: pts_ad ! 1 : potential temperature/salinity [Celcius/psu]
1099	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(inout) :: prd_ad ! TL of in_situ density [-]
1100	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(inout) :: prhop_ad ! TL of potential density (surface referenced)
1101	!! * Local declarations
1102	INTEGER :: ji, jj, jk ! dummy loop indices
1103	REAL(wp) :: & ! temporary scalars
1104	zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw, &
1105	zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, &
1106	ztad, zsad, zhad, zsrad, zr1ad, zr2ad, zr3ad, &
1107	zr4ad, zrhopad, zead, zbwad, &
1108	zbad, zdad, zcad, zawad, zaad, zb1ad, za1ad, &
1109	zkwad, zk0ad, zpes, zrdc1, zrdc2, zeps, &
1110	zmask, zrau0r
1111	REAL(wp), POINTER, DIMENSION(:,:,:) :: zws
1112	!!----------------------------------------------------------------------
1113	!
1114	IF( nn_timing == 1 ) CALL timing_start('eos-p_adj')
1115	!
1116	CALL wrk_alloc( jpi, jpj, jpk, zws )
1117	!
1118	! initialization of adjoint variables
1119	ztad = 0.0_wp
1120	zsad = 0.0_wp
1121	zhad = 0.0_wp
1122	zsrad = 0.0_wp
1123	zr1ad = 0.0_wp
1124	zr2ad = 0.0_wp
1125	zr3ad = 0.0_wp
1126	zr4ad = 0.0_wp
1127	zrhopad = 0.0_wp
1128	zead = 0.0_wp
1129	zbwad = 0.0_wp
1130	zbad = 0.0_wp
1131	zdad = 0.0_wp
1132	zcad = 0.0_wp
1133	zawad = 0.0_wp
1134	zaad = 0.0_wp
1135	zb1ad = 0.0_wp
1136	za1ad = 0.0_wp
1137	zkwad = 0.0_wp
1138	zk0ad = 0.0_wp
1139
1140	SELECT CASE ( nn_eos )
1141
1142	CASE ( 0 ) !== Jackett and McDougall (1994) formulation ==!
1143	zrau0r = 1.e0 / rau0
1144	#ifdef key_sp
1145	zeps = 1.e-7
1146	#else
1147	zeps = 1.e-14
1148	#endif
1149	!CDIR NOVERRCHK
1150	zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) )
1151	!
1152	DO jk = jpkm1, 1, -1
1153	DO jj = jpj, 1, -1
1154	DO ji = jpi, 1, -1
1155	! direct recomputing
1156	zt = pts(ji,jj,jk,jp_tem)
1157	zs = pts(ji,jj,jk,jp_sal)
1158	zh = fsdept(ji,jj,jk) ! depth
1159	zsr = zws(ji,jj,jk) ! square root salinity
1160	! compute volumic mass pure water at atm pressure
1161	zr1 = ( ( ( ( 6.536332e-9_wp * zt - 1.120083e-6_wp ) * zt &
1162	& + 1.001685e-4_wp ) * zt - 9.095290e-3_wp ) * zt &
1163	& + 6.793952e-2_wp ) * zt + 999.842594_wp
1164	! seawater volumic mass atm pressure
1165	zr2 = ( ( ( 5.3875e-9_wp * zt - 8.2467e-7_wp ) * zt &
1166	& + 7.6438e-5_wp ) * zt - 4.0899e-3_wp ) * zt + 0.824493_wp
1167	zr3 = ( -1.6546e-6_wp * zt + 1.0227e-4_wp ) * zt - 5.72466e-3_wp
1168	zr4 = 4.8314e-4_wp
1169	! potential volumic mass (reference to the surface)
1170	zrhop = ( zr4 * zs + zr3zsr + zr2 ) zs + zr1
1171	! add the compression terms
1172	ze = ( -3.508914e-8_wp * zt - 1.248266e-8_wp ) * zt - 2.595994e-6_wp
1173	zbw = ( 1.296821e-6_wp * zt - 5.782165e-9_wp ) * zt + 1.045941e-4_wp
1174	zb = zbw + ze * zs
1175
1176	zd = -2.042967e-2_wp
1177	zc = (-7.267926e-5_wp * zt + 2.598241e-3_wp ) * zt + 0.1571896_wp
1178	zaw= ( ( 5.939910e-6_wp * zt + 2.512549e-3_wp ) * zt - 0.1028859_wp &
1179	& ) * zt - 4.721788
1180	za = ( zd * zsr + zc ) * zs + zaw
1181
1182	zb1= (-0.1909078_wp * zt + 7.390729_wp ) * zt - 55.87545_wp
1183	za1= ( ( 2.326469e-3_wp * zt + 1.553190_wp ) * zt - 65.00517_wp &
1184	& ) * zt + 1044.077_wp
1185	zkw= ( ( (-1.361629e-4_wp * zt - 1.852732e-2_wp ) * zt - 30.41638_wp &
1186	& ) * zt + 2098.925_wp ) * zt + 190925.6_wp
1187	zk0= ( zb1 * zsr + za1 ) * zs + zkw
1188
1189
1190	zrdc1 = 1.0 / ( zk0 - zh * ( za - zh * zb ) )
1191	zrdc2 = 1.0 / ( 1.0 - zh * zrdc1 )
1192
1193	! ============
1194	! Adjoint part
1195	! ============
1196
1197	! Masked in situ density anomaly
1198
1199	zrhopad = zrhopad + prd_ad(ji,jj,jk) * tmask(ji,jj,jk) &
1200	& * zrdc2 * zrau0r
1201	zk0ad = zk0ad - prd_ad(ji,jj,jk) * tmask(ji,jj,jk) &
1202	& * zrdc2 * zrdc2 * zh &
1203	& * zrdc1*2 zrhop &
1204	& * zrau0r
1205	zaad = zaad + prd_ad(ji,jj,jk) * tmask(ji,jj,jk) &
1206	& * zrdc2 * zrdc2 * zh &
1207	& * zrdc1*2 zrhop &
1208	& * zh * zrau0r
1209	zbad = zbad - prd_ad(ji,jj,jk) * tmask(ji,jj,jk) &
1210	& * zrdc2 * zrdc2 * zh &
1211	& * zrdc1*2 zrhop &
1212	& * zh * zh * zrau0r
1213	prd_ad(ji,jj,jk) = 0.0_wp
1214
1215	zkwad = zkwad + zk0ad
1216	zsrad = zsrad + zk0ad * zb1 * zs
1217	zb1ad = zb1ad + zk0ad * zs * zsr
1218	za1ad = za1ad + zk0ad * zs
1219	zsad = zsad + zk0ad * ( zb1 * zsr + za1 )
1220	zk0ad = 0.0_wp
1221
1222	ztad = ztad + zkwad * ( ( (-4.1.361629e-4_wp zt &
1223	& -3.1.852732e-2_wp ) zt &
1224	& -2.30.41638_wp ) zt &
1225	& + 2098.925_wp )
1226	zkwad = 0.0_wp
1227
1228	ztad = ztad + za1ad * ( ( 3.2.326469e-3_wp zt &
1229	& +2.1.553190_wp ) zt &
1230	& - 65.00517_wp )
1231	za1ad = 0.0_wp
1232
1233	ztad = ztad + zb1ad * (-2.0.1909078_wp zt &
1234	& + 7.390729_wp )
1235	zb1ad = 0.0_wp
1236
1237	zawad = zawad + zaad
1238	zsrad = zsrad + zaad * zd * zs
1239	zcad = zcad + zaad * zs
1240	zsad = zsad + zaad * ( zd * zsr + zc )
1241	zaad = 0.0_wp
1242
1243	ztad = ztad + zawad * ( ( 3.5.939910e-6_wp zt &
1244	& +2.2.512549e-3_wp ) zt &
1245	& - 0.1028859_wp )
1246	zawad = 0.0_wp
1247
1248	ztad = ztad + zcad * (-2.7.267926e-5_wp zt &
1249	& + 2.598241e-3_wp )
1250	zcad = 0.0_wp
1251
1252
1253	zsad = zsad + zbad * ze
1254	zead = zead + zbad * zs
1255	zbwad = zbwad + zbad
1256	zbad = 0.0_wp
1257
1258	ztad = ztad + zbwad * ( 2.1.296821e-6_wp zt &
1259	& - 5.782165e-9_wp )
1260	zbwad = 0.0_wp
1261
1262	ztad = ztad + zead * (-2.3.508914e-8_wp zt &
1263	& - 1.248266e-8_wp )
1264	zead = 0.0_wp
1265
1266	zrhopad = zrhopad + prhop_ad(ji,jj,jk) * tmask(ji,jj,jk)
1267	prhop_ad(ji,jj,jk) = 0.0_wp
1268
1269	zr1ad = zr1ad + zrhopad
1270	zr2ad = zr2ad + zrhopad * zs
1271	zr3ad = zr3ad + zrhopad * zsr * zs
1272	zsrad = zsrad + zrhopad * zr3 * zs
1273	zsad = zsad + zrhopad * ( 2. * zr4 * zs + zr2 &
1274	& + zr3 * zsr )
1275	zrhopad = 0.0_wp
1276
1277	ztad = ztad + zr3ad * (-2.1.6546e-6_wp zt &
1278	& + 1.0227e-4_wp )
1279	zr3ad = 0.0_wp
1280
1281	ztad = ztad + zr2ad * ( ( ( 4.5.3875e-9_wp zt &
1282	& -3.8.2467e-7_wp ) zt &
1283	& +2.7.6438e-5_wp ) zt &
1284	& - 4.0899e-3_wp )
1285	zr2ad = 0.0_wp
1286
1287	ztad = ztad + zr1ad * ( ( ( ( 5.6.536332e-9_wp zt &
1288	& -4.1.120083e-6_wp ) zt &
1289	& +3.1.001685e-4_wp ) zt &
1290	& -2.9.095290e-3_wp ) zt &
1291	& + 6.793952e-2_wp )
1292	zr1ad = 0.0_wp
1293
1294	zsad = zsad + zsrad * ( 1.0 / MAX( 2.*zsr, zeps ) ) &
1295	& * tmask(ji,jj,jk)
1296	zsrad = 0.0_wp
1297
1298	pts_ad(ji,jj,jk,jp_sal) = pts_ad(ji,jj,jk,jp_sal) + zsad
1299	pts_ad(ji,jj,jk,jp_tem) = pts_ad(ji,jj,jk,jp_tem) + ztad
1300	ztad = 0.0_wp
1301	zsad = 0.0_wp
1302	END DO
1303	END DO
1304	END DO
1305	!
1306	CASE ( 1 ) !== Linear formulation = F( temperature ) ==!
1307	DO jk = jpkm1, 1, -1
1308	prd_ad(:,:,jk) = prd_ad(:,:,jk) + rau0 * prhop_ad(:,:,jk) * tmask(:,:,jk)
1309	prhop_ad(:,:,jk) = 0.0_wp
1310	pts_ad(:,:,jk,jp_tem) = pts_ad(:,:,jk,jp_tem) - rn_alpha * prd_ad(:,:,jk) * tmask(:,:,jk)
1311	prd_ad(:,:,jk) = 0.0_wp
1312	END DO
1313	!
1314	CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==!
1315	DO jk = jpkm1, 1, -1
1316	prd_ad( :,:,jk) = prd_ad(:,:,jk) + rau0 * prhop_ad(:,:,jk) * tmask(:,:,jk)
1317	prhop_ad(:,:,jk) = 0.0_wp
1318	pts_ad( :,:,jk,jp_tem) = pts_ad(:,:,jk,jp_tem) - rn_alpha * prd_ad(:,:,jk) * tmask(:,:,jk)
1319	pts_ad( :,:,jk,jp_sal) = pts_ad(:,:,jk,jp_sal) + rn_beta * prd_ad(:,:,jk) * tmask(:,:,jk)
1320	prd_ad( :,:,jk) = 0.0_wp
1321	END DO
1322	!
1323	END SELECT
1324	CALL wrk_dealloc( jpi, jpj, jpk, zws )
1325	!
1326	IF( nn_timing == 1 ) CALL timing_stop('eos-p_adj')
1327	!
1328	END SUBROUTINE eos_insitu_pot_adj
1329
1330	SUBROUTINE eos_insitu_2d_adj( pts, pdep, pts_ad, prd_ad )
1331	!!-----------------------------------------------------------------------
1332	!!
1333	!! * ROUTINE eos_insitu_2d_adj : adj OF ROUTINE eos_insitu_2d *
1334	!!
1335	!! ** Purpose of direct routine : Compute the in situ density
1336	!! (ratio rho/rau0) from potential temperature and salinity
1337	!! using an equation of state defined through the namelist
1338	!! parameter nn_eos. * 2D field case
1339	!!
1340	!! ** Method of direct routine : 3 cases:
1341	!! nn_eos = 0 : Jackett and McDougall (1994) equation of state.
1342	!! the in situ density is computed directly as a function of
1343	!! potential temperature relative to the surface (the opa t
1344	!! variable), salt and pressure (assuming no pressure variation
1345	!! along geopotential surfaces, i.e. the pressure p in decibars
1346	!! is approximated by the depth in meters.
1347	!! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0
1348	!! with pressure p decibars
1349	!! potential temperature t deg celsius
1350	!! salinity s psu
1351	!! reference volumic mass rau0 kg/m**3
1352	!! in situ volumic mass rho kg/m**3
1353	!! in situ density anomalie prd no units
1354	!! Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar,
1355	!! t = 40 deg celcius, s=40 psu
1356	!! nn_eos = 1 : linear equation of state function of temperature only
1357	!! prd(t) = 0.0285 - rn_alpha * t
1358	!! nn_eos = 2 : linear equation of state function of temperature and
1359	!! salinity
1360	!! prd(t,s) = rn_beta * s - rn_alpha * tn - 1.
1361	!! Note that no boundary condition problem occurs in this routine
1362	!! as (ptem,psal) are defined over the whole domain.
1363	!!
1364	!! ** Comments on Adjoint Routine :
1365	!! Care has been taken to avoid division by zero when computing
1366	!! the inverse of the square root of salinity at masked salinity
1367	!! points.
1368	!!
1369	!! ** Action :
1370	!!
1371	!! References :
1372	!!
1373	!! History :
1374	!! 8.2 ! 05-03 ((F. Van den Berghe, A. Weaver, N. Daget) - eosadj.F
1375	!! 9.0 ! 08-07 (A. Vidard) first version based on eosadj
1376	!!-----------------------------------------------------------------------
1377	!! * Modules used
1378	!! * Arguments
1379	REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in ) :: pts ! 1 : potential temperature [Celcius]
1380	! ! 2 : salinity [psu]
1381	REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(inout) :: pts_ad ! 1 : TL of potential temperature [Celcius]
1382	! ! 2 : TL of salinity [psu]
1383	REAL(wp), DIMENSION(jpi,jpj ), INTENT( inout) :: prd_ad ! TL of in_situ density [-]
1384	REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pdep ! depth [m]
1385	!
1386
1387	INTEGER :: ji, jj ! dummy loop indices
1388	REAL(wp) :: & ! temporary scalars
1389	zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw, &
1390	zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, &
1391	ztad, zsad, zhad, zsrad, zr1ad, zr2ad, zr3ad, &
1392	zr4ad, zrhopad, zead, zbwad, &
1393	zbad, zdad, zcad, zawad, zaad, zb1ad, za1ad, &
1394	zkwad, zk0ad, zpes, zrdc1, zrdc2, zeps, &
1395	zmask
1396	REAL(wp), POINTER, DIMENSION(:,:) :: zws
1397	!!----------------------------------------------------------------------
1398	!
1399	IF( nn_timing == 1 ) CALL timing_start('eos2d_adj')
1400	!
1401	CALL wrk_alloc( jpi, jpj, zws )
1402	!
1403	! initialization of adjoint variables
1404	ztad = 0.0_wp
1405	zsad = 0.0_wp
1406	zhad = 0.0_wp
1407	zsrad = 0.0_wp
1408	zr1ad = 0.0_wp
1409	zr2ad = 0.0_wp
1410	zr3ad = 0.0_wp
1411	zr4ad = 0.0_wp
1412	zrhopad = 0.0_wp
1413	zead = 0.0_wp
1414	zbwad = 0.0_wp
1415	zbad = 0.0_wp
1416	zdad = 0.0_wp
1417	zcad = 0.0_wp
1418	zawad = 0.0_wp
1419	zaad = 0.0_wp
1420	zb1ad = 0.0_wp
1421	za1ad = 0.0_wp
1422	zkwad = 0.0_wp
1423	zk0ad = 0.0_wp
1424	SELECT CASE ( nn_eos )
1425	CASE ( 0 ) !== Jackett and McDougall (1994) formulation ==!
1426
1427	#ifdef key_sp
1428	zeps = 1.e-7
1429	#else
1430	zeps = 1.e-14
1431	#endif
1432	!CDIR NOVERRCHK
1433	DO jj = jpjm1, 1, -1
1434	!CDIR NOVERRCHK
1435	DO ji = fs_jpim1, 1, -1 ! vector opt.
1436	zws(ji,jj) = SQRT( ABS( pts(ji,jj,jp_sal) ) )
1437	END DO
1438	END DO
1439	!
1440	DO jj = jpjm1, 1, -1
1441	DO ji = fs_jpim1, 1, -1 ! vector opt.
1442	zmask = tmask(ji,jj,1) ! land/sea bottom mask = surf. mask
1443	zt = pts (ji,jj,jp_tem) ! interpolated T
1444	zs = pts (ji,jj,jp_sal) ! interpolated S
1445	zsr= zws(ji,jj) ! square root of interpolated S
1446	zh = pdep(ji,jj) ! depth at the partial step level
1447	! compute volumic mass pure water at atm pressure
1448	zr1= ( ( ( ( 6.536332e-9_wpzt-1.120083e-6_wp )zt+1.001685e-4_wp)*zt &
1449	& -9.095290e-3_wp )zt+6.793952e-2_wp )zt+999.842594_wp
1450	! seawater volumic mass atm pressure
1451	zr2= ( ( ( 5.3875e-9_wpzt-8.2467e-7_wp ) zt+7.6438e-5_wp ) *zt &
1452	& -4.0899e-3_wp ) *zt+0.824493_wp
1453	zr3= ( -1.6546e-6_wpzt+1.0227e-4_wp ) zt-5.72466e-3_wp
1454	zr4= 4.8314e-4_wp
1455
1456	! potential volumic mass (reference to the surface)
1457	zrhop= ( zr4zs + zr3zsr + zr2 ) *zs + zr1
1458
1459	! add the compression terms
1460	ze = ( -3.508914e-8_wpzt-1.248266e-8_wp ) zt-2.595994e-6_wp
1461	zbw= ( 1.296821e-6_wpzt-5.782165e-9_wp ) zt+1.045941e-4_wp
1462	zb = zbw + ze * zs
1463
1464	zd = -2.042967e-2_wp
1465	zc = (-7.267926e-5_wpzt+2.598241e-3_wp ) zt+0.1571896_wp
1466	zaw= ( ( 5.939910e-6_wpzt+2.512549e-3_wp ) zt-0.1028859_wp ) *zt - 4.721788_wp
1467	za = ( zdzsr + zc ) zs + zaw
1468
1469	zb1= (-0.1909078_wpzt+7.390729_wp ) zt-55.87545_wp
1470	za1= ( ( 2.326469e-3_wpzt+1.553190_wp)zt-65.00517_wp ) *zt+1044.077_wp
1471	zkw= ( ( (-1.361629e-4_wpzt-1.852732e-2_wp ) zt-30.41638_wp ) zt + 2098.925_wp ) zt+190925.6_wp
1472	zk0= ( zb1zsr + za1 )zs + zkw
1473
1474	zrdc1 = 1.0 / ( zk0 - zh * ( za - zh * zb ) )
1475	zrdc2 = 1.0 / ( 1.0 - zh * zrdc1 )
1476	! ============
1477	! Adjoint part
1478	! ============
1479
1480	! Masked in situ density anomaly
1481
1482	zrhopad = zrhopad + prd_ad(ji,jj) * tmask(ji,jj,1) &
1483	& * zrdc2 / rau0
1484	zk0ad = zk0ad - prd_ad(ji,jj) * tmask(ji,jj,1) &
1485	& * zrdc2 * zrdc2 * zh &
1486	& * zrdc1*2 zrhop &
1487	& / rau0
1488	zaad = zaad + prd_ad(ji,jj) * tmask(ji,jj,1) &
1489	& * zrdc2 * zrdc2 * zh &
1490	& * zrdc1*2 zrhop &
1491	& * zh / rau0
1492	zbad = zbad - prd_ad(ji,jj) * tmask(ji,jj,1) &
1493	& * zrdc2 * zrdc2 * zh &
1494	& * zrdc1*2 zrhop &
1495	& * zh * zh / rau0
1496	prd_ad(ji,jj) = 0.0_wp
1497
1498	zkwad = zkwad + zk0ad
1499	zsrad = zsrad + zk0ad * zb1 * zs
1500	zb1ad = zb1ad + zk0ad * zs * zsr
1501	za1ad = za1ad + zk0ad * zs
1502	zsad = zsad + zk0ad * ( zb1 * zsr + za1 )
1503	zk0ad = 0.0_wp
1504
1505	ztad = ztad + zkwad * ( ( (-4.1.361629e-4_wp zt &
1506	& -3.1.852732e-2_wp ) zt &
1507	& -2.30.41638_wp ) zt &
1508	& + 2098.925_wp )
1509	zkwad = 0.0_wp
1510
1511	ztad = ztad + za1ad * ( ( 3.2.326469e-3_wp zt &
1512	& +2.1.553190_wp ) zt &
1513	& - 65.00517_wp )
1514	za1ad = 0.0_wp
1515
1516	ztad = ztad + zb1ad * (-2.0.1909078_wp zt &
1517	& + 7.390729_wp )
1518	zb1ad = 0.0_wp
1519
1520	zawad = zawad + zaad
1521	zsrad = zsrad + zaad * zd * zs
1522	zcad = zcad + zaad * zs
1523	zsad = zsad + zaad * ( zd * zsr + zc )
1524	zaad = 0.0_wp
1525
1526	ztad = ztad + zawad * ( ( 3.5.939910e-6_wp zt &
1527	& +2.2.512549e-3_wp ) zt &
1528	& - 0.1028859_wp )
1529	zawad = 0.0_wp
1530
1531	ztad = ztad + zcad * (-2.7.267926e-5_wp zt &
1532	& + 2.598241e-3_wp )
1533	zcad = 0.0_wp
1534
1535	zbwad = zbwad + zbad
1536	zead = zead + zbad * zs
1537	zsad = zsad + zbad * ze
1538	zbad = 0.0_wp
1539
1540	ztad = ztad + zbwad * ( 2.1.296821e-6_wp zt &
1541	& - 5.782165e-9_wp )
1542	zbwad = 0.0_wp
1543
1544	ztad = ztad + zead * (-2.3.508914e-8_wp zt &
1545	& - 1.248266e-8_wp )
1546	zead = 0.0_wp
1547
1548	zr1ad = zr1ad + zrhopad
1549	zr2ad = zr2ad + zrhopad * zs
1550	zr3ad = zr3ad + zrhopad * zsr * zs
1551	zsrad = zsrad + zrhopad * zr3 * zs
1552	zsad = zsad + zrhopad * ( 2. * zr4 * zs + zr2 &
1553	& + zr3 * zsr )
1554	zrhopad = 0.0_wp
1555
1556	ztad = ztad + zr3ad * (-2.1.6546e-6_wp zt &
1557	& + 1.0227e-4_wp )
1558	zr3ad = 0.0_wp
1559
1560	ztad = ztad + zr2ad * ( ( ( 4.5.3875e-9_wp zt &
1561	& -3.8.2467e-7_wp ) zt &
1562	& +2.7.6438e-5_wp ) zt &
1563	& - 4.0899e-3_wp )
1564	zr2ad = 0.0_wp
1565
1566	ztad = ztad + zr1ad * ( ( ( ( 5.6.536332e-9_wp zt &
1567	& -4.1.120083e-6_wp ) zt &
1568	& +3.1.001685e-4_wp ) zt &
1569	& -2.9.095290e-3_wp ) zt &
1570	& + 6.793952e-2_wp )
1571	zr1ad = 0.0_wp
1572
1573	zsad = zsad + zsrad * ( 1.0 / MAX( 2.*zsr, zeps ) ) &
1574	& * tmask(ji,jj, 1)
1575	zsrad = 0.0_wp
1576
1577	pts_ad(ji,jj,jp_sal) = pts_ad(ji,jj,jp_sal) + zsad
1578	pts_ad(ji,jj,jp_tem) = pts_ad(ji,jj,jp_tem) + ztad
1579	ztad = 0.0_wp
1580	zsad = 0.0_wp
1581	END DO
1582	END DO
1583	!
1584	CASE ( 1 ) !== Linear formulation = F( temperature ) ==!
1585	DO jj = jpjm1, 1, -1
1586	DO ji = fs_jpim1, 1, -1 ! vector opt.
1587	pts_ad(ji,jj,jp_tem) = pts_ad(ji,jj,jp_tem) - prd_ad(ji,jj) * rn_alpha * tmask(ji,jj,1)
1588	prd_ad(ji,jj) = 0.0_wp
1589	END DO
1590	END DO
1591	!
1592	CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==!
1593	DO jj = jpjm1, 1, -1
1594	DO ji = fs_jpim1, 1, -1 ! vector opt.
1595	pts_ad(ji,jj,jp_tem) = pts_ad(ji,jj,jp_tem) - prd_ad(ji,jj) * rn_alpha * tmask(ji,jj,1)
1596	pts_ad(ji,jj,jp_sal) = pts_ad(ji,jj,jp_sal) + prd_ad(ji,jj) * rn_beta * tmask(ji,jj,1)
1597	prd_ad (ji,jj) = 0.0_wp
1598	END DO
1599	END DO
1600	!
1601	END SELECT
1602	!
1603	CALL wrk_dealloc( jpi, jpj, zws )
1604	!
1605	IF( nn_timing == 1 ) CALL timing_stop('eos2d_adj')
1606	!
1607	END SUBROUTINE eos_insitu_2d_adj
1608
1609	SUBROUTINE eos_bn2_tan ( pts, pts_tl, pn2_tl )
1610	!!----------------------------------------------------------------------
1611	!! * ROUTINE eos_bn2_tan *
1612	!!
1613	!! ** Purpose of the direct routine: Compute the local
1614	!! Brunt-Vaisala frequency at the time-step of the input arguments
1615	!!
1616	!! ** Method of the direct routine:
1617	!! * nn_eos = 0 : UNESCO sea water properties
1618	!! The brunt-vaisala frequency is computed using the polynomial
1619	!! polynomial expression of McDougall (1987):
1620	!! N^2 = grav * beta * ( alpha/beta*dk[ t ] - dk[ s ] )/e3w
1621	!! If lk_zdfddm=T, the heat/salt buoyancy flux ratio Rrau is
1622	!! computed and used in zdfddm module :
1623	!! Rrau = alpha/beta * ( dk[ t ] / dk[ s ] )
1624	!! * nn_eos = 1 : linear equation of state (temperature only)
1625	!! N^2 = grav * rn_alpha * dk[ t ]/e3w
1626	!! * nn_eos = 2 : linear equation of state (temperature & salinity)
1627	!! N^2 = grav * (rn_alpha * dk[ t ] - rn_beta * dk[ s ] ) / e3w
1628	!! The use of potential density to compute N^2 introduces e r r o r
1629	!! in the sign of N^2 at great depths. We recommand the use of
1630	!! nn_eos = 0, except for academical studies.
1631	!! Macro-tasked on horizontal slab (jk-loop)
1632	!! N.B. N^2 is set to zero at the first level (JK=1) in inidtr
1633	!! and is never used at this level.
1634	!!
1635	!! ** Action : - pn2 : the brunt-vaisala frequency
1636	!!
1637	!! References :
1638	!! McDougall, T. J., J. Phys. Oceanogr., 17, 1950-1964, 1987.
1639	!!
1640	!! History:
1641	!! ! 08-07 (A. Vidard) First version
1642	!!----------------------------------------------------------------------
1643	!! * Arguments
1644	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts, & ! 1 : potential temperature [Celcius]
1645	! ! 2 : salinity [psu]
1646	& pts_tl ! 1 : TL of potential temperature [Celsius]
1647	! 2 : TL of salinity [psu]
1648	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out ) :: &
1649	& pn2_tl ! TL of potential density (surface referenced)
1650	!! * Local declarations
1651	INTEGER :: ji, jj, jk ! dummy loop indices
1652	REAL(wp) :: &
1653	zgde3w, zt, zs, zh, & ! temporary scalars
1654	zalbet, zbeta ! " "
1655	REAL(wp) :: &
1656	zttl, zstl, & ! temporary scalars
1657	zalbettl, zbetatl ! " "
1658	#if defined key_zdfddm
1659	REAL(wp) :: zds, zdstl ! temporary scalars
1660	#endif
1661
1662	! pn2_tl : interior points only (2=< jk =< jpkm1 )
1663	! --------------------------
1664	SELECT CASE ( nn_eos )
1665
1666	CASE ( 0 ) !== Jackett and McDougall (1994) formulation ==!
1667	DO jk = 2, jpkm1
1668	DO jj = 1, jpj
1669	DO ji = 1, jpi
1670	zgde3w = grav / fse3w(ji,jj,jk)
1671	zt = 0.5 * ( pts(ji,jj,jk,jp_tem) + pts(ji,jj,jk-1,jp_tem) ) ! potential temperature at w-point
1672	zs = 0.5 * ( pts(ji,jj,jk,jp_sal) + pts(ji,jj,jk-1,jp_sal) ) - 35.0 ! salinity anomaly (s-35) at w-point
1673	zh = fsdepw(ji,jj,jk) ! depth in meters at w-point
1674
1675	zalbet = ( ( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt & ! ratio alpha/beta
1676	& - 0.203814e-03_wp ) * zt &
1677	& + 0.170907e-01_wp ) * zt &
1678	& + 0.665157e-01_wp &
1679	& + ( - 0.678662e-05_wp * zs &
1680	& - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs &
1681	& + ( ( - 0.302285e-13_wp * zh &
1682	& - 0.251520e-11_wp * zs &
1683	& + 0.512857e-12_wp * zt * zt ) * zh &
1684	& - 0.164759e-06_wp * zs &
1685	& +( 0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt &
1686	& + 0.380374e-04_wp ) * zh
1687
1688	zbeta = ( ( -0.415613e-09_wp * zt + 0.555579e-07_wp ) * zt & ! beta
1689	& - 0.301985e-05_wp ) * zt &
1690	& + 0.785567e-03_wp &
1691	& + ( 0.515032e-08_wp * zs &
1692	& + 0.788212e-08_wp * zt - 0.356603e-06_wp ) * zs &
1693	& +( ( 0.121551e-17_wp * zh &
1694	& - 0.602281e-15_wp * zs &
1695	& - 0.175379e-14_wp * zt + 0.176621e-12_wp ) * zh &
1696	& + 0.408195e-10_wp * zs &
1697	& + ( - 0.213127e-11_wp * zt + 0.192867e-09_wp ) * zt &
1698	& - 0.121555e-07_wp ) * zh
1699
1700	!! tangent part
1701	zttl = 0.5 * ( pts_tl(ji,jj,jk,jp_tem) + pts_tl(ji,jj,jk-1,jp_tem) ) ! potential temperature at w-point
1702	zstl = 0.5 * ( pts_tl(ji,jj,jk,jp_sal) + pts_tl(ji,jj,jk-1,jp_sal) ) ! salinity anomaly at w-point
1703	zalbettl = ( ( ( -4.0.255019e-07_wp zt & ! ratio alpha/beta
1704	& +3.0.298357e-05_wp ) zt &
1705	& -2.0.203814e-03_wp ) zt &
1706	& + 0.170907e-01_wp &
1707	& - 0.846960e-04_wp * zs &
1708	& - ( 0.933746e-06_wp &
1709	& - ( 2.*0.791325e-08_wp &
1710	& +2.0.512857e-12_wp zh ) * zt ) * zh ) * zttl &
1711	& + ( - 2.0.678662e-05_wp zs &
1712	& - 0.846960e-04_wp * zt &
1713	& + 0.378110e-02_wp &
1714	& + ( - 0.164759e-06_wp &
1715	& - 0.251520e-11_wp * zh ) * zh ) * zstl
1716
1717	zbetatl = ( ( -3.0.415613e-09_wp zt &
1718	& +2.0.555579e-07_wp ) zt &
1719	& - 0.301985e-05_wp &
1720	& + 0.788212e-08_wp * zs &
1721	& + ( -2.0.213127e-11_wp zt &
1722	& - 0.175379e-14_wp * zh &
1723	& + 0.192867e-09_wp ) * zh ) * zttl &
1724	& + ( 2.0.515032e-08_wp zs &
1725	& + 0.788212e-08_wp * zt &
1726	& - 0.356603e-06_wp &
1727	& + ( - 0.602281e-15_wp * zh &
1728	& + 0.408195e-10_wp ) * zh ) * zstl
1729
1730	pn2_tl(ji,jj,jk) = zgde3w * tmask(ji,jj,jk) * ( &
1731	& zbeta * ( zalbet &
1732	& * ( pts_tl(ji,jj,jk-1,jp_tem) - pts_tl(ji,jj,jk,jp_tem) ) &
1733	& + zalbettl &
1734	& * ( pts (ji,jj,jk-1,jp_tem) - pts (ji,jj,jk,jp_tem) ) &
1735	& - ( pts_tl(ji,jj,jk-1,jp_sal) - pts_tl(ji,jj,jk,jp_sal) ) ) &
1736	& + zbetatl * ( zalbet &
1737	& * ( pts (ji,jj,jk-1,jp_tem) - pts (ji,jj,jk,jp_tem) ) &
1738	& - ( pts (ji,jj,jk-1,jp_sal) - pts (ji,jj,jk,jp_sal) ) ) )
1739	#if defined key_zdfddm
1740	zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )
1741	zdstl = ( pts_tl(ji,jj,jk-1,jp_sal) - pts_tl(ji,jj,jk,jp_sal) )
1742	IF ( ABS( zds) <= 1.e-20 ) THEN
1743	zds = 1.e-20
1744	rrau_tl(ji,jj,jk) = zalbettl * &
1745	& ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds &
1746	& + zalbet * &
1747	& ( ( pts_tl(ji,jj,jk-1,jp_tem) - pts_tl(ji,jj,jk,jp_tem) ) / zds )
1748	ELSE
1749	rrau_tl(ji,jj,jk) = zalbettl * &
1750	& ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds &
1751	& + zalbet * &
1752	& ( ( pts_tl(ji,jj,jk-1,jp_tem) - pts_tl(ji,jj,jk,jp_tem) ) / zds &
1753	& - ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) * zdstl/ zds**2 )
1754	ENDIF
1755	#endif
1756	END DO
1757	END DO
1758	END DO
1759	!
1760	CASE ( 1 ) !== Linear formulation = F( temperature ) ==!
1761	DO jk = 2, jpkm1
1762	pn2_tl(:,:,jk) = rn_alpha * ( pts_tl(:,:,jk-1,jp_tem) - pts_tl(:,:,jk,jp_tem) ) * grav / fse3w(:,:,jk) * tmask(:,:,jk)
1763	END DO
1764	!
1765	CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==!
1766	DO jk = 2, jpkm1
1767	pn2_tl(:,:,jk) = ( rn_alpha * ( pts_tl(:,:,jk-1,jp_tem) - pts_tl(:,:,jk,jp_tem) ) &
1768	& - rn_beta * ( pts_tl(:,:,jk-1,jp_sal) - pts_tl(:,:,jk,jp_sal) ) ) &
1769	& * grav / fse3w(:,:,jk) * tmask(:,:,jk)
1770	END DO
1771	#if defined key_zdfddm
1772	DO jk = 2, jpkm1
1773	DO jj = 1, jpj
1774	DO ji = 1, jpi
1775	zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )
1776	zdstl = pts_tl(ji,jj,jk-1,jp_sal) - pts_tl(ji,jj,jk,jp_sal)
1777	IF ( ABS( zds) <= 1.e-20 ) THEN
1778	zds = 1.e-20
1779	rrau_tl(ji,jj,jk) = ralpbet * &
1780	& ( ( pts_tl(ji,jj,jk-1,jp_tem) - pts_tl(ji,jj,jk,jp_tem) ) / zds )
1781	ELSE
1782	rrau_tl(ji,jj,jk) = ralpbet * &
1783	& ( ( pts_tl(ji,jj,jk-1,jp_tem) - pts_tl(ji,jj,jk,jp_tem) ) / zds &
1784	& - ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) * zdstl / zds**2 )
1785	ENDIF
1786	rrau(ji,jj,jk) = ralpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds
1787	END DO
1788	END DO
1789	END DO
1790	#endif
1791	END SELECT
1792	END SUBROUTINE eos_bn2_tan
1793
1794	SUBROUTINE eos_bn2_adj ( pts, pts_ad, pn2_ad )
1795	!!----------------------------------------------------------------------
1796	!! * ROUTINE eos_bn2_adj *
1797	!!
1798	!! ** Purpose of the direct routine: Compute the local
1799	!! Brunt-Vaisala frequency at the time-step of the input arguments
1800	!!
1801	!! ** Method of the direct routine:
1802	!! * nn_eos = 0 : UNESCO sea water properties
1803	!! The brunt-vaisala frequency is computed using the polynomial
1804	!! polynomial expression of McDougall (1987):
1805	!! N^2 = grav * beta * ( alpha/beta*dk[ t ] - dk[ s ] )/e3w
1806	!! If lk_zdfddm=T, the heat/salt buoyancy flux ratio Rrau is
1807	!! computed and used in zdfddm module :
1808	!! Rrau = alpha/beta * ( dk[ t ] / dk[ s ] )
1809	!! * nn_eos = 1 : linear equation of state (temperature only)
1810	!! N^2 = grav * rn_alpha * dk[ t ]/e3w
1811	!! * nn_eos = 2 : linear equation of state (temperature & salinity)
1812	!! N^2 = grav * (rn_alpha * dk[ t ] - rn_beta * dk[ s ] ) / e3w
1813	!! The use of potential density to compute N^2 introduces e r r o r
1814	!! in the sign of N^2 at great depths. We recommand the use of
1815	!! nn_eos = 0, except for academical studies.
1816	!! Macro-tasked on horizontal slab (jk-loop)
1817	!! N.B. N^2 is set to zero at the first level (JK=1) in inidtr
1818	!! and is never used at this level.
1819	!!
1820	!! ** Action : - pn2 : the brunt-vaisala frequency
1821	!!
1822	!! References :
1823	!! McDougall, T. J., J. Phys. Oceanogr., 17, 1950-1964, 1987.
1824	!!
1825	!! History:
1826	!! ! 08-07 (A. Vidard) First version
1827	!!----------------------------------------------------------------------
1828	!! * Arguments
1829	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! 1 : potential temperature [Celcius]
1830	! ! 2 : salinity [psu]
1831	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(inout ) :: pts_ad ! 1 : Adjoint of potential temperature [Celsius]
1832	! 2 : Adjoint of salinity [psu]
1833	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: &
1834	& pn2_ad ! Adjoint of potential density (surface referenced)
1835	!
1836	!! * Local declarations
1837	INTEGER :: ji, jj, jk ! dummy loop indices
1838	REAL(wp) :: &
1839	zgde3w, zt, zs, zh, & ! temporary scalars
1840	zalbet, zbeta ! " "
1841	REAL(wp) :: &
1842	ztad, zsad, & ! temporary scalars
1843	zalbetad, zbetaad ! " "
1844	#if defined key_zdfddm
1845	REAL(wp) :: zds, zdsad ! temporary scalars
1846	#endif
1847
1848	! pn2_tl : interior points only (2=< jk =< jpkm1 )
1849	! --------------------------
1850	zalbetad = 0.0_wp
1851	zbetaad = 0.0_wp
1852	ztad = 0.0_wp
1853	zsad = 0.0_wp
1854	#if defined key_zdfddm
1855	zdsad = 0.0_wp
1856	#endif
1857
1858	SELECT CASE ( nn_eos )
1859
1860	CASE ( 0 ) !== Jackett and McDougall (1994) formulation ==!
1861	DO jk = jpkm1, 2, -1
1862	DO jj = jpj, 1, -1
1863	DO ji = jpi, 1, -1
1864	zgde3w = grav / fse3w(ji,jj,jk)
1865	zt = 0.5 * ( pts(ji,jj,jk,jp_tem) + pts(ji,jj,jk-1,jp_tem) ) ! potential temperature at w-point
1866	zs = 0.5 * ( pts(ji,jj,jk,jp_sal) + pts(ji,jj,jk-1,jp_sal) ) - 35.0 ! salinity anomaly (s-35) at w-point
1867	zh = fsdepw(ji,jj,jk) ! depth in meters at w-point
1868
1869	zalbet = ( ( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt & ! ratio alpha/beta
1870	& - 0.203814e-03_wp ) * zt &
1871	& + 0.170907e-01_wp ) * zt &
1872	& + 0.665157e-01_wp &
1873	& + ( - 0.678662e-05_wp * zs &
1874	& - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs &
1875	& + ( ( - 0.302285e-13_wp * zh &
1876	& - 0.251520e-11_wp * zs &
1877	& + 0.512857e-12_wp * zt * zt ) * zh &
1878	& - 0.164759e-06_wp * zs &
1879	& +( 0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt &
1880	& + 0.380374e-04_wp ) * zh
1881
1882	zbeta = ( ( -0.415613e-09_wp * zt + 0.555579e-07_wp ) * zt & ! beta
1883	& - 0.301985e-05_wp ) * zt &
1884	& + 0.785567e-03_wp &
1885	& + ( 0.515032e-08_wp * zs &
1886	& + 0.788212e-08_wp * zt - 0.356603e-06_wp ) * zs &
1887	& +( ( 0.121551e-17_wp * zh &
1888	& - 0.602281e-15_wp * zs &
1889	& - 0.175379e-14_wp * zt + 0.176621e-12_wp ) * zh &
1890	& + 0.408195e-10_wp * zs &
1891	& + ( - 0.213127e-11_wp * zt + 0.192867e-09_wp ) * zt &
1892	& - 0.121555e-07_wp ) * zh
1893
1894	#if defined key_zdfddm
1895
1896	zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )
1897	IF ( ABS( zds) <= 1.e-20 ) THEN
1898	zds = 1.e-20
1899	zdsad = 0.0_wp
1900	ELSE
1901	zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )
1902	zdsad = rrau_ad(ji,jj,jk) * zalbet ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds*2
1903	zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )
1904	ENDIF
1905	pts_ad(ji,jj,jk-1,jp_tem) = pts_ad(ji,jj,jk-1,jp_tem) + rrau_ad(ji,jj,jk) * zalbet / zds
1906	pts_ad(ji,jj,jk,jp_tem ) = pts_ad(ji,jj,jk,jp_tem ) - rrau_ad(ji,jj,jk) * zalbet / zds
1907	zalbetad = zalbetad + rrau_ad(ji,jj,jk) * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds
1908	rrau_ad(ji,jj,jk) = 0._wp
1909	pts_ad(ji,jj,jk-1,jp_sal) = pts_ad(ji,jj,jk-1,jp_sal) + zdsad
1910	pts_ad(ji,jj,jk,jp_sal ) = pts_ad(ji,jj,jk,jp_sal ) - zdsad
1911	zdsad = 0._wp
1912	#endif
1913	pts_ad(ji,jj,jk-1,jp_tem) = pts_ad(ji,jj,jk-1,jp_tem) + zalbetzbetazgde3wtmask(ji,jj,jk)pn2_ad(ji,jj,jk)
1914	pts_ad(ji,jj,jk,jp_tem ) = pts_ad(ji,jj,jk,jp_tem ) - zalbetzbetazgde3wtmask(ji,jj,jk)pn2_ad(ji,jj,jk)
1915	zalbetad = zalbetad + zbetazgde3wtmask(ji,jj,jk) &
1916	& ( pts (ji,jj,jk-1,jp_tem) - pts (ji,jj,jk,jp_tem) ) pn2_ad(ji,jj,jk)
1917	pts_ad(ji,jj,jk-1,jp_sal) = pts_ad(ji,jj,jk-1,jp_sal) - zbetatmask(ji,jj,jk)zgde3w*pn2_ad(ji,jj,jk)
1918	pts_ad(ji,jj,jk,jp_sal ) = pts_ad(ji,jj,jk,jp_sal ) + zbetatmask(ji,jj,jk)zgde3w*pn2_ad(ji,jj,jk)
1919	zbetaad = zbetaad &
1920	& + zgde3w tmask(ji,jj,jk) ( zalbet * ( pts (ji,jj,jk-1,jp_tem) - pts (ji,jj,jk,jp_tem) ) &
1921	& - ( pts (ji,jj,jk-1,jp_sal) - pts (ji,jj,jk,jp_sal) ) )*pn2_ad(ji,jj,jk)
1922
1923	pn2_ad(ji,jj,jk) = 0.0_wp
1924
1925	ztad = ztad + ( ( -3.0.415613e-09_wp zt &
1926	& +2.0.555579e-07_wp ) zt &
1927	& - 0.301985e-05_wp &
1928	& + 0.788212e-08_wp * zs &
1929	& + ( -2.0.213127e-11_wp zt &
1930	& - 0.175379e-14_wp * zh &
1931	& + 0.192867e-09_wp ) * zh ) *zbetaad
1932
1933	zsad = zsad + ( 2.0.515032e-08_wp zs &
1934	& + 0.788212e-08_wp * zt &
1935	& - 0.356603e-06_wp &
1936	& + ( - 0.602281e-15_wp * zh &
1937	& + 0.408195e-10_wp ) * zh ) * zbetaad
1938
1939	zbetaad = 0.0_wp
1940
1941	ztad = ztad + ( ( ( -4.0.255019e-07_wp zt &! ratio alpha/beta
1942	& +3.0.298357e-05_wp ) zt &
1943	& -2.0.203814e-03_wp ) zt &
1944	& + 0.170907e-01_wp &
1945	& - 0.846960e-04_wp * zs &
1946	& - ( 0.933746e-06_wp &
1947	& - ( 2.*0.791325e-08_wp &
1948	& +2.0.512857e-12_wp zh ) * zt ) * zh &
1949	& ) *zalbetad
1950
1951	zsad = zsad + ( - 2.0.678662e-05_wp zs &
1952	& - 0.846960e-04_wp * zt &
1953	& + 0.378110e-02_wp &
1954	& + ( - 0.164759e-06_wp &
1955	& - 0.251520e-11_wp * zh ) * zh &
1956	& ) *zalbetad
1957
1958	zalbetad = 0.0_wp
1959
1960
1961	pts_ad(ji,jj,jk,jp_sal) = pts_ad(ji,jj,jk,jp_sal) + 0.5 * zsad
1962	pts_ad(ji,jj,jk-1,jp_sal) = pts_ad(ji,jj,jk-1,jp_sal) + 0.5 * zsad
1963	zsad = 0.0_wp
1964
1965	pts_ad(ji,jj,jk,jp_tem) = pts_ad(ji,jj,jk,jp_tem) + 0.5 * ztad
1966	pts_ad(ji,jj,jk-1,jp_tem) = pts_ad(ji,jj,jk-1,jp_tem) + 0.5 * ztad
1967	ztad = 0.0_wp
1968
1969	END DO
1970	END DO
1971	END DO
1972	!
1973	CASE ( 1 ) !== Linear formulation = F( temperature ) ==!
1974	DO jk = jpkm1, 2, -1
1975	pts_ad(:,:,jk-1,jp_tem) = pts_ad(:,:,jk-1,jp_tem) + rn_alpha * pn2_ad(:,:,jk) &
1976	& * grav / fse3w(:,:,jk) * tmask(:,:,jk)
1977	pts_ad(:,:,jk,jp_tem ) = pts_ad(:,:,jk,jp_tem ) - rn_alpha * pn2_ad(:,:,jk) &
1978	& * grav / fse3w(:,:,jk) * tmask(:,:,jk)
1979	pn2_ad(:,:,jk) = 0.0_wp
1980	END DO
1981	!
1982	CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==!
1983	#if defined key_zdfddm
1984	DO jk = jpkm1, 2, -1
1985	DO jj = jpj, 1, -1
1986	DO ji = jpi, 1, -1
1987	zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )
1988	IF ( ABS( zds) <= 1.e-20 ) THEN
1989	zds = 1.e-20
1990	zdsad = 0.0_wp
1991	ELSE
1992	zdsad = zdsad - rrau_ad(ji,jj,jk) * ralpbet &
1993	& * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds**2
1994	ENDIF
1995	rrau(ji,jj,jk) = ralpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds
1996	pts_ad(ji,jj,jk-1,jp_tem) = pts_ad(ji,jj,jk-1,jp_tem) &
1997	& + rrau_ad(ji,jj,jk) * ralpbet / zds
1998	pts_ad(ji,jj,jk,jp_tem ) = pts_ad(ji,jj,jk,jp_tem ) &
1999	& - rrau_ad(ji,jj,jk) * ralpbet / zds
2000	rrau_ad(ji,jj,jk) = 0._wp
2001
2002	pts_ad(ji,jj,jk-1,jp_sal) = pts_ad(ji,jj,jk-1,jp_sal) + zdsad
2003	pts_ad(ji,jj,jk,jp_sal ) = pts_ad(ji,jj,jk,jp_sal ) - zdsad
2004	zdsad = 0._wp
2005	END DO
2006	END DO
2007	END DO
2008	#endif
2009	DO jk = jpkm1, 2, -1
2010	pts_ad(:,:,jk-1,jp_tem) = pts_ad(:,:,jk-1,jp_tem) + rn_alpha * pn2_ad(:,:,jk) &
2011	& * grav / fse3w(:,:,jk) * tmask(:,:,jk)
2012	pts_ad(:,:,jk,jp_tem ) = pts_ad(:,:,jk,jp_tem ) - rn_alpha * pn2_ad(:,:,jk) &
2013	& * grav / fse3w(:,:,jk) * tmask(:,:,jk)
2014	pts_ad(:,:,jk-1,jp_sal) = pts_ad(:,:,jk-1,jp_sal) - rn_beta * pn2_ad(:,:,jk) &
2015	& * grav / fse3w(:,:,jk) * tmask(:,:,jk)
2016	pts_ad(:,:,jk,jp_sal ) = pts_ad(:,:,jk,jp_sal ) + rn_beta * pn2_ad(:,:,jk) &
2017	& * grav / fse3w(:,:,jk) * tmask(:,:,jk)
2018	pn2_ad(:,:,jk) = 0.0_wp
2019	END DO
2020	END SELECT
2021	END SUBROUTINE eos_bn2_adj
2022
2023	SUBROUTINE eos_alpbet_tan( pts, pts_tl, palpbet_tl, beta0_tl )
2024	!!----------------------------------------------------------------------
2025	!! * ROUTINE eos_alpbet_tan *
2026	!!
2027	!! ** Purpose of the direct routine :
2028	!! Calculates the in situ thermal/haline expansion ratio at T-points
2029	!!
2030	!! ** Method of the direct routine :
2031	!! calculates alpha / beta ratio at T-points
2032	!! * nn_eos = 0 : UNESCO sea water properties
2033	!! The alpha/beta ratio is returned as 3-D array palpbet using the polynomial
2034	!! polynomial expression of McDougall (1987).
2035	!! Scalar beta0 is returned = 1.
2036	!! * nn_eos = 1 : linear equation of state (temperature only)
2037	!! The ratio is undefined, so we return alpha as palpbet
2038	!! Scalar beta0 is returned = 0.
2039	!! * nn_eos = 2 : linear equation of state (temperature & salinity)
2040	!! The alpha/beta ratio is returned as ralpbet
2041	!! Scalar beta0 is returned = 1.
2042	!!
2043	!! ** Action : - palpbet : thermal/haline expansion ratio at T-points
2044	!! : beta0 : 1. or 0.
2045	!!----------------------------------------------------------------------
2046	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts_tl, & ! linear tangent of pot. temperature & salinity
2047	& pts ! pot. temperature & salinity
2048	REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: palpbet_tl ! thermal/haline expansion ratio
2049	REAL(wp), INTENT( out) :: beta0_tl ! set = 1 except with case 1 eos, rho=rho(T)
2050	!!
2051	INTEGER :: ji, jj, jk ! dummy loop indices
2052	REAL(wp) :: zt, zs, zh, & ! local scalars
2053	& zttl, zstl
2054	!!----------------------------------------------------------------------
2055	!
2056	IF( nn_timing == 1 ) CALL timing_start('eos_alpbet_tan')
2057	!
2058	SELECT CASE ( nn_eos )
2059	!
2060	CASE ( 0 ) ! Jackett and McDougall (1994) formulation
2061	DO jk = 1, jpk
2062	DO jj = 1, jpj
2063	DO ji = 1, jpi
2064	zt = pts(ji,jj,jk,jp_tem) ! potential temperature
2065	zs = pts(ji,jj,jk,jp_sal) - 35._wp ! salinity anomaly (s-35)
2066	zh = fsdept(ji,jj,jk) ! depth in meters
2067	!! Tangent part
2068	zttl = pts_tl(ji,jj,jk,jp_tem) ! potential temperature
2069	zstl = pts_tl(ji,jj,jk,jp_sal) ! salinity anomaly (s-35)
2070	palpbet_tl(ji,jj,jk) = &
2071	& ( ( ( ( - 4. * 0.255019e-07_wp * zt &
2072	& + 3. * 0.298357e-05_wp ) * zt &
2073	& - 2. * 0.203814e-03_wp ) * zt &
2074	& + 0.170907e-01_wp * zt ) &
2075	& - 0.846960e-04_wp * zs &
2076	& + ( ( 2. * 0.512857e-12_wp * zt ) * zh &
2077	& + ( 2. * 0.791325e-08_wp * zt &
2078	& - 0.933746e-06_wp ) ) * zh ) * zttl &
2079	& + ( ( - 2. * 0.678662e-05_wp * zs &
2080	& - 0.846960e-04_wp * zt &
2081	& + 0.378110e-02_wp ) &
2082	& + ( - 0.251520e-11_wp * zh &
2083	& - 0.164759e-06_wp ) * zh ) * zstl
2084	END DO
2085	END DO
2086	END DO
2087	beta0_tl = 0._wp
2088	!
2089	CASE ( 1 ) !== Linear formulation = F( temperature ) ==!
2090	palpbet_tl(:,:,:) = 0._wp
2091	beta0_tl = 0._wp
2092	!
2093	CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==!
2094	palpbet_tl(:,:,:) = 0._wp
2095	beta0_tl = 0._wp
2096	!
2097	CASE DEFAULT
2098	IF(lwp) WRITE(numout,cform_err)
2099	IF(lwp) WRITE(numout,*) ' bad flag value for nn_eos = ', nn_eos
2100	nstop = nstop + 1
2101	!
2102	END SELECT
2103	!
2104	IF( nn_timing == 1 ) CALL timing_stop('eos_alpbet_tan')
2105	!
2106	END SUBROUTINE eos_alpbet_tan
2107
2108	SUBROUTINE eos_alpbet_adj( pts, pts_ad, palpbet_ad, beta0_ad )
2109	!!----------------------------------------------------------------------
2110	!! * ROUTINE eos_alpbet_adj *
2111	!!
2112	!! ** Purpose of the direct routine :
2113	!! Calculates the in situ thermal/haline expansion ratio at T-points
2114	!!
2115	!! ** Method of the direct routine :
2116	!! calculates alpha / beta ratio at T-points
2117	!! * nn_eos = 0 : UNESCO sea water properties
2118	!! The alpha/beta ratio is returned as 3-D array palpbet using the polynomial
2119	!! polynomial expression of McDougall (1987).
2120	!! Scalar beta0 is returned = 1.
2121	!! * nn_eos = 1 : linear equation of state (temperature only)
2122	!! The ratio is undefined, so we return alpha as palpbet
2123	!! Scalar beta0 is returned = 0.
2124	!! * nn_eos = 2 : linear equation of state (temperature & salinity)
2125	!! The alpha/beta ratio is returned as ralpbet
2126	!! Scalar beta0 is returned = 1.
2127	!!
2128	!! ** Action : - palpbet : thermal/haline expansion ratio at T-points
2129	!! : beta0 : 1. or 0.
2130	!!----------------------------------------------------------------------
2131	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(inout) :: pts_ad ! linear tangent of pot. temperature & salinity
2132	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in) :: pts ! pot. temperature & salinity
2133	REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(inout) :: palpbet_ad ! thermal/haline expansion ratio
2134	REAL(wp), INTENT(inout) :: beta0_ad ! set = 1 except with case 1 eos, rho=rho(T)
2135	!!
2136	INTEGER :: ji, jj, jk ! dummy loop indices
2137	REAL(wp) :: zt, zs, zh, & ! local scalars
2138	& ztad, zsad
2139	!!----------------------------------------------------------------------
2140	!
2141	ztad = 0.0_wp
2142	zsad = 0.0_wp
2143	IF( nn_timing == 1 ) CALL timing_start('eos_alpbet_adj')
2144	!
2145	SELECT CASE ( nn_eos )
2146	!
2147	CASE ( 0 ) ! Jackett and McDougall (1994) formulation
2148	DO jk = jpk, 1, -1
2149	DO jj = jpj, 1, -1
2150	DO ji = jpi, 1, -1
2151	zt = pts(ji,jj,jk,jp_tem) ! potential temperature
2152	zs = pts(ji,jj,jk,jp_sal) - 35._wp ! salinity anomaly (s-35)
2153	zh = fsdept(ji,jj,jk) ! depth in meters
2154	!! Adjoint part
2155	ztad = ztad + ( ( ( ( - 4. * 0.255019e-07_wp * zt &
2156	& + 3. * 0.298357e-05_wp ) * zt &
2157	& - 2. * 0.203814e-03_wp ) * zt &
2158	& + 0.170907e-01_wp * zt ) &
2159	& - 0.846960e-04_wp * zs &
2160	& + ( ( 2. * 0.512857e-12_wp * zt ) * zh &
2161	& + ( 2. * 0.791325e-08_wp * zt &
2162	& - 0.933746e-06_wp ) ) * zh ) * palpbet_ad(ji,jj,jk)
2163	zsad = zsad + ( ( - 2. * 0.678662e-05_wp * zs &
2164	& - 0.846960e-04_wp * zt &
2165	& + 0.378110e-02_wp ) &
2166	& + ( - 0.251520e-11_wp * zh &
2167	& - 0.164759e-06_wp ) * zh ) * palpbet_ad(ji,jj,jk)
2168	palpbet_ad(ji,jj,jk) = 0.0_wp
2169	pts_ad(ji,jj,jk,jp_tem) = pts_ad(ji,jj,jk,jp_tem) + ztad
2170	pts_ad(ji,jj,jk,jp_sal) = pts_ad(ji,jj,jk,jp_sal) + zsad
2171	ztad = 0.0_wp
2172	zsad = 0.0_wp
2173	END DO
2174	END DO
2175	END DO
2176	beta0_ad = 0._wp
2177	!
2178	CASE ( 1 ) !== Linear formulation = F( temperature ) ==!
2179	palpbet_ad(:,:,:) = 0._wp
2180	beta0_ad = 0._wp
2181	!
2182	CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==!
2183	palpbet_ad(:,:,:) = 0._wp
2184	beta0_ad = 0._wp
2185	!
2186	CASE DEFAULT
2187	IF(lwp) WRITE(numout,cform_err)
2188	IF(lwp) WRITE(numout,*) ' bad flag value for nn_eos = ', nn_eos
2189	nstop = nstop + 1
2190	!
2191	END SELECT
2192	!
2193	IF( nn_timing == 1 ) CALL timing_stop('eos_alpbet_adj')
2194	!
2195	END SUBROUTINE eos_alpbet_adj
2196
2197	#if defined key_tam
2198	SUBROUTINE eos_insitu_adj_tst( kumadt )
2199	!!-----------------------------------------------------------------------
2200	!!
2201	!! * ROUTINE eos_adj_tst *
2202	!!
2203	!! ** Purpose : Test the adjoint routine.
2204	!!
2205	!! ** Method : Verify the scalar product
2206	!!
2207	!! ( L dx )^T W dy = dx^T L^T W dy
2208	!!
2209	!! where L = tangent routine
2210	!! L^T = adjoint routine
2211	!! W = diagonal matrix of scale factors
2212	!! dx = input perturbation (random field)
2213	!! dy = L dx
2214	!!
2215	!!
2216	!! History :
2217	!! ! 08-07 (A. Vidard)
2218	!!-----------------------------------------------------------------------
2219	!! * Modules used
2220
2221	!! * Arguments
2222	INTEGER, INTENT(IN) :: &
2223	& kumadt ! Output unit
2224	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
2225	zts ! potential temperature
2226	! salinity
2227	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
2228	& zts_adout ! potential temperature
2229	! salinity
2230	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2231	& zrd_adin ! anomaly density
2232	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
2233	& zts_tlin ! potential temperature
2234	! salinity
2235	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
2236	& znts ! potential temperature
2237	! salinity
2238	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2239	& zrd_tlout ! anomaly density
2240	REAL(KIND=wp) :: &
2241	& zsp1, & ! scalar product involving the tangent routine
2242	& zsp2 ! scalar product involving the adjoint routine
2243	INTEGER :: &
2244	& ji, &
2245	& jj, &
2246	& jk, &
2247	& jn, &
2248	& jeos
2249	CHARACTER(LEN=14) :: cl_name
2250	ALLOCATE( &
2251	& zts( jpi, jpj, jpk, jpts ), &
2252	& znts( jpi, jpj, jpk, jpts ), &
2253	& zts_adout( jpi, jpj, jpk,jpts ), &
2254	& zrd_adin( jpi, jpj, jpk ), &
2255	& zts_tlin( jpi, jpj, jpk,jpts ), &
2256	& zrd_tlout(jpi, jpj, jpk ) )
2257	! Initialize the reference state
2258	zts = tsn
2259	! store initial nn_eos
2260	jeos = nn_eos
2261	DO jn = 0, 2
2262	nn_eos = jn
2263	!=============================================================
2264	! 1) dx = ( T ) and dy = ( T )
2265	!=============================================================
2266
2267	!--------------------------------------------------------------------
2268	! Reset the tangent and adjoint variables
2269	!--------------------------------------------------------------------
2270	zts_tlin(:,:,:,:) = 0.0_wp
2271	zrd_tlout(:,:,:) = 0.0_wp
2272	zts_adout(:,:,:,:) = 0.0_wp
2273	zrd_adin(:,:,:) = 0.0_wp
2274
2275	!--------------------------------------------------------------------
2276	! Initialize the tangent input with random noise: dx
2277	!--------------------------------------------------------------------
2278	CALL grid_random( znts(:,:,:,jp_tem), 'T', 0.0_wp, stdt )
2279	CALL grid_random( znts(:,:,:,jp_sal), 'T', 0.0_wp, stds )
2280
2281	DO jk = 1, jpk
2282	DO jj = nldj, nlej
2283	DO ji = nldi, nlei
2284	zts_tlin(ji,jj,jk,jp_tem) = znts(ji,jj,jk,jp_tem)
2285	zts_tlin(ji,jj,jk,jp_sal) = znts(ji,jj,jk,jp_sal)
2286	END DO
2287	END DO
2288	END DO
2289	CALL eos_insitu_tan(zts, zts_tlin, zrd_tlout)
2290
2291	DO jk = 1, jpk
2292	DO jj = nldj, nlej
2293	DO ji = nldi, nlei
2294	zrd_adin(ji,jj,jk) = zrd_tlout(ji,jj,jk) &
2295	& * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)&
2296	& * tmask(ji,jj,jk)
2297	END DO
2298	END DO
2299	END DO
2300
2301	!--------------------------------------------------------------------
2302	! Compute the scalar product: ( L dx )^T W dy
2303	!--------------------------------------------------------------------
2304
2305	zsp1 = DOT_PRODUCT( zrd_tlout, zrd_adin )
2306
2307	!--------------------------------------------------------------------
2308	! Call the adjoint routine: dx^* = L^T dy^*
2309	!--------------------------------------------------------------------
2310	CALL eos_insitu_adj(zts, zts_adout, zrd_adin)
2311	zsp2 = DOT_PRODUCT( zts_tlin(:,:,:,jp_tem), zts_adout(:,:,:,jp_tem) ) + &
2312	& DOT_PRODUCT( zts_tlin(:,:,:,jp_sal), zts_adout(:,:,:,jp_sal) )
2313
2314	! Compare the scalar products
2315
2316	! Compare the scalar products
2317	! 14 char:'12345678901234'
2318	SELECT CASE( jn )
2319	CASE (0) ; cl_name = 'eos_adj ins T1'
2320	CASE (1) ; cl_name = 'eos_adj ins T2'
2321	CASE (2) ; cl_name = 'eos_adj ins T3'
2322	END SELECT
2323	CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
2324	ENDDO
2325	! restore initial nn_eos
2326	nn_eos = jeos
2327
2328	! Deallocate memory
2329
2330	DEALLOCATE( &
2331	& zts, &
2332	& zts_adout, &
2333	& zrd_adin, &
2334	& zts_tlin, &
2335	& zrd_tlout, &
2336	& znts &
2337	& )
2338	END SUBROUTINE eos_insitu_adj_tst
2339
2340	SUBROUTINE eos_insitu_pot_adj_tst( kumadt )
2341	!!-----------------------------------------------------------------------
2342	!!
2343	!! * ROUTINE eos_adj_tst *
2344	!!
2345	!! ** Purpose : Test the adjoint routine.
2346	!!
2347	!! ** Method : Verify the scalar product
2348	!!
2349	!! ( L dx )^T W dy = dx^T L^T W dy
2350	!!
2351	!! where L = tangent routine
2352	!! L^T = adjoint routine
2353	!! W = diagonal matrix of scale factors
2354	!! dx = input perturbation (random field)
2355	!! dy = L dx
2356	!!
2357	!! ** Action : Separate tests are applied for the following dx and dy:
2358	!!
2359	!! 1) dx = ( SSH ) and dy = ( SSH )
2360	!!
2361	!! History :
2362	!! ! 08-07 (A. Vidard)
2363	!!-----------------------------------------------------------------------
2364	!! * Modules used
2365
2366	!! * Arguments
2367	INTEGER, INTENT(IN) :: &
2368	& kumadt ! Output unit
2369	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
2370	zts ! potential temperature
2371	! salinity
2372	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
2373	& zts_adout ! potential temperature
2374	! salinity
2375	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2376	& zrd_adin ! anomaly density
2377	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2378	& zrhop_adin ! volume mass
2379	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
2380	& zts_tlin ! potential temperature
2381	! salinity
2382	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
2383	& znts ! potential temperature
2384	! salinity
2385	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2386	& zrd_tlout ! anomaly density
2387	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2388	& zrhop_tlout ! volume mass
2389	REAL(KIND=wp) :: &
2390	& zsp1, & ! scalar product involving the tangent routine
2391	& zsp2 ! scalar product involving the adjoint routine
2392	INTEGER :: &
2393	& ji, &
2394	& jj, &
2395	& jk, &
2396	& jn, &
2397	& jeos
2398	CHARACTER(LEN=14) :: cl_name
2399
2400	! Allocate memory
2401	ALLOCATE( &
2402	& zts( jpi, jpj, jpk, jpts ), &
2403	& zts_adout( jpi, jpj, jpk, jpts ), &
2404	& zrhop_adin( jpi, jpj, jpk ), &
2405	& zrd_adin( jpi, jpj, jpk ), &
2406	& zts_tlin( jpi, jpj, jpk, jpts ), &
2407	& znts( jpi, jpj, jpk, jpts ), &
2408	& zrd_tlout(jpi, jpj, jpk ), &
2409	& zrhop_tlout(jpi, jpj, jpk ) )
2410
2411	! Initialize random field standard deviationsthe reference state
2412	zts = tsn
2413
2414	! store initial nn_eos
2415	jeos = nn_eos
2416	DO jn = 0, 2
2417	nn_eos = jn
2418	!=============================================
2419	! testing of eos_insitu_pot
2420	!=============================================
2421
2422	!=============================================================
2423	! 1) dx = ( T ) and dy = ( T )
2424	!=============================================================
2425
2426	!--------------------------------------------------------------------
2427	! Reset the tangent and adjoint variables
2428	!--------------------------------------------------------------------
2429	zts_tlin(:,:,:,:) = 0.0_wp
2430	zrd_tlout(:,:,:) = 0.0_wp
2431	zrhop_tlout(:,:,:) = 0.0_wp
2432	zts_adout(:,:,:,:) = 0.0_wp
2433	zrhop_adin(:,:,:) = 0.0_wp
2434	zrd_adin(:,:,:) = 0.0_wp
2435
2436	!--------------------------------------------------------------------
2437	! Initialize the tangent input with random noise: dx
2438	!--------------------------------------------------------------------
2439	CALL grid_random( znts(:,:,:,jp_tem), 'T', 0.0_wp, stdt )
2440	CALL grid_random( znts(:,:,:,jp_sal), 'T', 0.0_wp, stds )
2441	DO jk = 1, jpk
2442	DO jj = nldj, nlej
2443	DO ji = nldi, nlei
2444	zts_tlin(ji,jj,jk,:) = znts(ji,jj,jk,:)
2445	END DO
2446	END DO
2447	END DO
2448	!--------------------------------------------------------------------
2449	! Call the tangent routine: dy = L dx
2450	!--------------------------------------------------------------------
2451
2452	call eos_insitu_pot_tan ( zts, zts_tlin, zrd_tlout, zrhop_tlout )
2453
2454	!--------------------------------------------------------------------
2455	! Initialize the adjoint variables: dy^* = W dy
2456	!--------------------------------------------------------------------
2457	DO jk = 1, jpk
2458	DO jj = nldj, nlej
2459	DO ji = nldi, nlei
2460	zrd_adin(ji,jj,jk) = zrd_tlout(ji,jj,jk) &
2461	& * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)&
2462	& * tmask(ji,jj,jk)
2463	zrhop_adin(ji,jj,jk) = zrhop_tlout(ji,jj,jk) &
2464	& * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)&
2465	& * tmask(ji,jj,jk)
2466	END DO
2467	END DO
2468	END DO
2469
2470	!--------------------------------------------------------------------
2471	! Compute the scalar product: ( L dx )^T W dy
2472	!--------------------------------------------------------------------
2473
2474	zsp1 = DOT_PRODUCT( zrd_tlout , zrd_adin ) &
2475	& + DOT_PRODUCT( zrhop_tlout, zrhop_adin )
2476	!--------------------------------------------------------------------
2477	! Call the adjoint routine: dx^* = L^T dy^*
2478	!--------------------------------------------------------------------
2479
2480	CALL eos_insitu_pot_adj( zts, zts_adout, zrd_adin, zrhop_adin )
2481	!--------------------------------------------------------------------
2482	! Compute the scalar product: dx^T L^T W dy
2483	!--------------------------------------------------------------------
2484
2485	zsp2 = DOT_PRODUCT( zts_tlin(:,:,:,jp_tem), zts_adout(:,:,:,jp_tem) ) + &
2486	& DOT_PRODUCT( zts_tlin(:,:,:,jp_sal), zts_adout(:,:,:,jp_sal) )
2487	! Compare the scalar products
2488
2489	! 14 char:'12345678901234'
2490	SELECT CASE( jn )
2491	CASE (0) ; cl_name = 'eos_adj pot T1'
2492	CASE (1) ; cl_name = 'eos_adj pot T2'
2493	CASE (2) ; cl_name = 'eos_adj pot T3'
2494	END SELECT
2495	CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
2496
2497	ENDDO
2498
2499	! restore initial nn_eos
2500	nn_eos = jeos
2501
2502	! Deallocate memory
2503	DEALLOCATE( &
2504	& zts, &
2505	& zts_adout, &
2506	& zrd_adin, &
2507	& zrhop_adin, &
2508	& zts_tlin, &
2509	& zrd_tlout, &
2510	& zrhop_tlout,&
2511	& znts )
2512	END SUBROUTINE eos_insitu_pot_adj_tst
2513
2514	SUBROUTINE eos_alpbet_adj_tst(kumadt)
2515	!!-----------------------------------------------------------------------
2516	!!
2517	!! * ROUTINE eos_adj_tst *
2518	!!
2519	!! ** Purpose : Test the adjoint routine.
2520	!!
2521	!! ** Method : Verify the scalar product
2522	!!
2523	!! ( L dx )^T W dy = dx^T L^T W dy
2524	!!
2525	!! where L = tangent routine
2526	!! L^T = adjoint routine
2527	!! W = diagonal matrix of scale factors
2528	!! dx = input perturbation (random field)
2529	!! dy = L dx
2530	!!
2531	!! ** Action : Separate tests are applied for the following dx and dy:
2532	!!
2533	!! 1) dx = ( SSH ) and dy = ( SSH )
2534	!!
2535	!! History :
2536	!! ! 05-12 (P.-A. Bouttier)
2537	!!-----------------------------------------------------------------------
2538	!! * Modules used
2539
2540	!! * Arguments
2541	INTEGER, INTENT(IN) :: &
2542	& kumadt ! Output unit
2543	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
2544	& zts ! potential temperature
2545	! salinity
2546	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
2547	& zts_adout ! potential temperature
2548	! salinity
2549	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2550	& zpalpbet_adin
2551	REAL(wp) :: &
2552	& zbeta0_adin
2553	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
2554	& zts_tlin ! potential temperature
2555	! salinity
2556	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
2557	& znts ! potential temperature
2558	! salinity
2559	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2560	& zpalpbet_tlout
2561	REAL(wp) :: &
2562	& zbeta0_tlout
2563	REAL(KIND=wp) :: &
2564	& zsp1, & ! scalar product involving the tangent routine
2565	& zsp2 ! scalar product involving the adjoint routine
2566	INTEGER :: &
2567	& ji, &
2568	& jj, &
2569	& jk, &
2570	& jn, &
2571	& jeos
2572	CHARACTER(LEN=14) :: cl_name
2573
2574	! Allocate memory
2575	ALLOCATE( &
2576	& zts( jpi, jpj, jpk, jpts ), &
2577	& zts_adout( jpi, jpj, jpk, jpts ), &
2578	& zpalpbet_adin( jpi, jpj, jpk ), &
2579	& zts_tlin( jpi, jpj, jpk, jpts ), &
2580	& znts( jpi, jpj, jpk, jpts ), &
2581	& zpalpbet_tlout(jpi, jpj, jpk ) )
2582
2583
2584	! Initialize random field standard deviationsthe reference state
2585	zts = tsn
2586
2587	! store initial nn_eos
2588	jeos = nn_eos
2589	DO jn = 0, 2
2590	nn_eos = jn
2591	!=============================================
2592	! testing of eos_insitu_pot
2593	!=============================================
2594
2595	!=============================================================
2596	! 1) dx = ( T ) and dy = ( T )
2597	!=============================================================
2598
2599	!--------------------------------------------------------------------
2600	! Reset the tangent and adjoint variables
2601	!--------------------------------------------------------------------
2602	zts_tlin(:,:,:,:) = 0.0_wp
2603	zpalpbet_tlout(:,:,:) = 0.0_wp
2604	zbeta0_tlout = 0.0_wp
2605	zts_adout(:,:,:,:) = 0.0_wp
2606	zbeta0_adin = 0.0_wp
2607	zpalpbet_adin(:,:,:) = 0.0_wp
2608
2609	!--------------------------------------------------------------------
2610	! Initialize the tangent input with random noise: dx
2611	!--------------------------------------------------------------------
2612	CALL grid_random( znts(:,:,:,jp_tem), 'T', 0.0_wp, stdt )
2613	CALL grid_random( znts(:,:,:,jp_sal), 'T', 0.0_wp, stds )
2614	DO jk = 1, jpk
2615	DO jj = nldj, nlej
2616	DO ji = nldi, nlei
2617	zts_tlin(ji,jj,jk,:) = znts(ji,jj,jk,:)
2618	END DO
2619	END DO
2620	END DO
2621	!--------------------------------------------------------------------
2622	! Call the tangent routine: dy = L dx
2623	!--------------------------------------------------------------------
2624
2625	call eos_alpbet_tan ( zts, zts_tlin, zpalpbet_tlout, zbeta0_tlout )
2626
2627	!--------------------------------------------------------------------
2628	! Initialize the adjoint variables: dy^* = W dy
2629	!--------------------------------------------------------------------
2630	DO jk = 1, jpk
2631	DO jj = nldj, nlej
2632	DO ji = nldi, nlei
2633	zpalpbet_adin(ji,jj,jk) = zpalpbet_tlout(ji,jj,jk) &
2634	& * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)&
2635	& * tmask(ji,jj,jk)
2636	END DO
2637	END DO
2638	END DO
2639	zbeta0_adin = zbeta0_tlout
2640
2641	!--------------------------------------------------------------------
2642	! Compute the scalar product: ( L dx )^T W dy
2643	!--------------------------------------------------------------------
2644
2645	zsp1 = DOT_PRODUCT( zpalpbet_tlout , zpalpbet_adin ) &
2646	& + zbeta0_tlout * zbeta0_adin
2647	!--------------------------------------------------------------------
2648	! Call the adjoint routine: dx^* = L^T dy^*
2649	!--------------------------------------------------------------------
2650
2651	CALL eos_alpbet_adj( zts, zts_adout, zpalpbet_adin, zbeta0_adin )
2652
2653	!--------------------------------------------------------------------
2654	! Compute the scalar product: dx^T L^T W dy
2655	!--------------------------------------------------------------------
2656
2657	zsp2 = DOT_PRODUCT( zts_tlin(:,:,:,jp_tem), zts_adout(:,:,:,jp_tem) ) + &
2658	& DOT_PRODUCT( zts_tlin(:,:,:,jp_sal), zts_adout(:,:,:,jp_sal) )
2659	! Compare the scalar products
2660
2661	! 14 char:'12345678901234'
2662	SELECT CASE( jn )
2663	CASE (0) ; cl_name = 'eos_adj ab T1'
2664	CASE (1) ; cl_name = 'eos_adj ab T2'
2665	CASE (2) ; cl_name = 'eos_adj ab T3'
2666	END SELECT
2667	CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
2668
2669	ENDDO
2670
2671	! restore initial nn_eos
2672	nn_eos = jeos
2673
2674	! Deallocate memory
2675	DEALLOCATE( &
2676	& zts, &
2677	& zts_adout, &
2678	& zpalpbet_adin, &
2679	& zts_tlin, &
2680	& zpalpbet_tlout, &
2681	& znts )
2682
2683	END SUBROUTINE eos_alpbet_adj_tst
2684
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	zts ! potential temperature
2718	! salinity
2719	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2720	& zts_adout ! potential temperature
2721	! salinity
2722	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
2723	& zrd_adin ! anomaly density
2724	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2725	& zts_tlin ! potential temperature
2726	! salinity
2727	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2728	& znts ! potential temperature
2729	! 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	INTEGER :: &
2736	& ji, &
2737	& jj, &
2738	& jn, &
2739	& jeos
2740	CHARACTER(LEN=14) :: cl_name
2741	! Allocate memory
2742
2743	ALLOCATE( &
2744	& zdep( jpi, jpj), &
2745	& zts( jpi, jpj, jpts ), &
2746	& znts( jpi, jpj, jpts ), &
2747	& zts_adout( jpi, jpj, jpts ), &
2748	& zrd_adin( jpi, jpj ), &
2749	& zts_tlin( jpi, jpj, jpts ), &
2750	& zrd_tlout(jpi, jpj ) )
2751
2752	! Initialize the reference state
2753	zts(:,:,:) = tsn(:,:,2,:)
2754	zdep(:,:) = fsdept(:,:,2)
2755
2756	! store initial nn_eos
2757	jeos = nn_eos
2758	DO jn = 0, 2
2759	nn_eos = jn
2760	!=============================================================
2761	! 1) dx = ( T ) and dy = ( T )
2762	!=============================================================
2763
2764	!--------------------------------------------------------------------
2765	! Reset the tangent and adjoint variables
2766	!--------------------------------------------------------------------
2767	zts_tlin(:,:,:) = 0.0_wp
2768	zrd_tlout(:,:) = 0.0_wp
2769	zts_adout(:,:,:) = 0.0_wp
2770	zrd_adin(:,:) = 0.0_wp
2771
2772	!--------------------------------------------------------------------
2773	! Initialize the tangent input with random noise: dx
2774	!--------------------------------------------------------------------
2775	CALL grid_random( znts(:,:,jp_tem), 'T', 0.0_wp, stdt )
2776	CALL grid_random( znts(:,:,jp_sal), 'T', 0.0_wp, stds )
2777	DO jj = nldj, nlej
2778	DO ji = nldi, nlei
2779	zts_tlin(ji,jj,:) = znts(ji,jj,:)
2780	END DO
2781	END DO
2782
2783	CALL eos_insitu_2d_tan(zts, zdep, zts_tlin, zrd_tlout)
2784
2785	DO jj = nldj, nlej
2786	DO ji = nldi, nlei
2787	zrd_adin(ji,jj) = zrd_tlout(ji,jj) &
2788	& * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,2)&
2789	& * tmask(ji,jj,2)
2790	END DO
2791	END DO
2792
2793	!--------------------------------------------------------------------
2794	! Compute the scalar product: ( L dx )^T W dy
2795	!--------------------------------------------------------------------
2796
2797	zsp1 = DOT_PRODUCT( zrd_tlout, zrd_adin )
2798
2799	!--------------------------------------------------------------------
2800	! Call the adjoint routine: dx^* = L^T dy^*
2801	!--------------------------------------------------------------------
2802	CALL eos_insitu_2d_adj(zts, zdep, zts_adout, zrd_adin)
2803	zsp2 = DOT_PRODUCT( zts_tlin(:,:,jp_tem), zts_adout(:,:,jp_tem) ) + &
2804	& DOT_PRODUCT( zts_tlin(:,:,jp_sal), zts_adout(:,:,jp_sal) )
2805
2806	! Compare the scalar products
2807
2808	! 14 char:'12345678901234'
2809	SELECT CASE( jn )
2810	CASE (0) ; cl_name = 'eos_adj 2d T1'
2811	CASE (1) ; cl_name = 'eos_adj 2d T2'
2812	CASE (2) ; cl_name = 'eos_adj 2d T3'
2813	END SELECT
2814	CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
2815
2816	ENDDO
2817
2818	! restore initial nn_eos
2819	nn_eos = jeos
2820
2821	! Deallocate memory
2822
2823	DEALLOCATE( &
2824	& zdep, &
2825	& zts, &
2826	& zts_adout, &
2827	& zrd_adin, &
2828	& zts_tlin, &
2829	& zrd_tlout, &
2830	& znts )
2831	END SUBROUTINE eos_insitu_2d_adj_tst
2832
2833	SUBROUTINE eos_adj_tst( kumadt )
2834	!!-----------------------------------------------------------------------
2835	!!
2836	!! * ROUTINE eos_adj_tst *
2837	!!
2838	!! ** Purpose : Test the adjoint routine.
2839	!!
2840	!! History :
2841	!! ! 08-07 (A. Vidard)
2842	!!-----------------------------------------------------------------------
2843	!! * Arguments
2844	INTEGER, INTENT(IN) :: &
2845	& kumadt ! Output unit
2846
2847	CALL eos_insitu_adj_tst( kumadt )
2848	CALL eos_insitu_pot_adj_tst( kumadt )
2849	CALL eos_insitu_2d_adj_tst( kumadt )
2850	CALL eos_alpbet_adj_tst( kumadt )
2851
2852	END SUBROUTINE eos_adj_tst
2853
2854	SUBROUTINE bn2_adj_tst( kumadt )
2855	!!-----------------------------------------------------------------------
2856	!!
2857	!! * ROUTINE bn2_adj_tst *
2858	!!
2859	!! ** Purpose : Test the adjoint routine.
2860	!!
2861	!! ** Method : Verify the scalar product
2862	!!
2863	!! ( L dx )^T W dy = dx^T L^T W dy
2864	!!
2865	!! where L = tangent routine
2866	!! L^T = adjoint routine
2867	!! W = diagonal matrix of scale factors
2868	!! dx = input perturbation (random field)
2869	!! dy = L dx
2870	!!
2871	!! ** Action : Separate tests are applied for the following dx and dy:
2872	!!
2873	!! 1) dx = ( SSH ) and dy = ( SSH )
2874	!!
2875	!! History :
2876	!! ! 08-07 (A. Vidard)
2877	!!-----------------------------------------------------------------------
2878	!! * Modules used
2879
2880	!! * Arguments
2881	INTEGER, INTENT(IN) :: &
2882	& kumadt ! Output unit
2883	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
2884	zts ! potential temperature
2885	! salinity
2886	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
2887	& zts_adout ! potential temperature
2888	! salinity
2889	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2890	& zrd_adin, & ! potential density (surface referenced)
2891	& zrd_adout ! potential density (surface referenced)
2892	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
2893	& zts_tlin, & ! potential temperature
2894	! salinity
2895	& zts_tlout ! potential temperature
2896	! salinity
2897	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2898	& zrd_tlout ! potential density (surface referenced)
2899	REAL(KIND=wp) :: &
2900	& zsp1, & ! scalar product involving the tangent routine
2901	& zsp2 ! scalar product involving the adjoint routine
2902	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
2903	& znts ! potential temperature
2904	! salinity
2905	INTEGER :: &
2906	& ji, &
2907	& jj, &
2908	& jk, &
2909	& jn, &
2910	& jeos
2911	CHARACTER(LEN=14) :: cl_name
2912
2913	! Allocate memory
2914	ALLOCATE( &
2915	& zts( jpi, jpj, jpk, jpts ), &
2916	& zts_adout( jpi, jpj, jpk, jpts ), &
2917	& zrd_adin( jpi, jpj, jpk ), &
2918	& zrd_adout(jpi, jpj, jpk ), &
2919	& zts_tlin( jpi, jpj, jpk, jpts ), &
2920	& znts( jpi, jpj, jpk, jpts ), &
2921	& zts_tlout( jpi, jpj, jpk, jpts ), &
2922	& zrd_tlout(jpi, jpj, jpk ) )
2923
2924	! Initialize random field standard deviationsthe reference state
2925	zts = tsn
2926	! store initial nn_eos
2927	jeos = nn_eos
2928	DO jn = 0, 2
2929	nn_eos = jn
2930	!=============================================================
2931	! 1) dx = ( T ) and dy = ( T )
2932	!=============================================================
2933
2934	!--------------------------------------------------------------------
2935	! Reset the tangent and adjoint variables
2936	!--------------------------------------------------------------------
2937	zts_tlin(:,:,:,:) = 0.0_wp
2938	zts_tlout(:,:,:,:) = 0.0_wp
2939	zrd_tlout(:,:,:) = 0.0_wp
2940	zts_adout(:,:,:,:) = 0.0_wp
2941	zrd_adin(:,:,:) = 0.0_wp
2942	zrd_adout(:,:,:) = 0.0_wp
2943
2944	!--------------------------------------------------------------------
2945	! Initialize the tangent input with random noise: dx
2946	!--------------------------------------------------------------------
2947	CALL grid_random( znts(:,:,:,jp_tem), 'T', 0.0_wp, stdt )
2948	CALL grid_random( znts(:,:,:,jp_sal), 'T', 0.0_wp, stds )
2949	DO jk = 1, jpk
2950	DO jj = nldj, nlej
2951	DO ji = nldi, nlei
2952	zts_tlin(ji,jj,jk,:) = znts(ji,jj,jk,:)
2953	END DO
2954	END DO
2955	END DO
2956	!--------------------------------------------------------------------
2957	! Call the tangent routine: dy = L dx
2958	!--------------------------------------------------------------------
2959	zts_tlout(:,:,:,:) = zts_tlin
2960	CALL eos_bn2_tan( zts, zts_tlout, zrd_tlout )
2961	!--------------------------------------------------------------------
2962	! Initialize the adjoint variables: dy^* = W dy
2963	!--------------------------------------------------------------------
2964	DO jk = 1, jpk
2965	DO jj = nldj, nlej
2966	DO ji = nldi, nlei
2967	zrd_adin(ji,jj,jk) = zrd_tlout(ji,jj,jk) &
2968	& * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
2969	& * tmask(ji,jj,jk)
2970	END DO
2971	END DO
2972	END DO
2973
2974	!--------------------------------------------------------------------
2975	! Compute the scalar product: ( L dx )^T W dy
2976	!--------------------------------------------------------------------
2977
2978	zsp1 = DOT_PRODUCT( zrd_tlout, zrd_adin )
2979
2980	!--------------------------------------------------------------------
2981	! Call the adjoint routine: dx^* = L^T dy^*
2982	!--------------------------------------------------------------------
2983
2984	zrd_adout(:,:,:) = zrd_adin(:,:,:)
2985	CALL eos_bn2_adj( zts, zts_adout, zrd_adout )
2986
2987	!--------------------------------------------------------------------
2988	! Compute the scalar product: dx^T L^T W dy
2989	!--------------------------------------------------------------------
2990	zsp2 = DOT_PRODUCT( zts_tlin(:,:,:,jp_tem), zts_adout(:,:,:,jp_tem )) + &
2991	& DOT_PRODUCT( zts_tlin(:,:,:,jp_sal), zts_adout(:,:,:,jp_sal ))
2992
2993	! Compare the scalar products
2994	! 14 char:'12345678901234'
2995	SELECT CASE( jn )
2996	CASE (0) ; cl_name = 'bn2_adj T1'
2997	CASE (1) ; cl_name = 'bn2_adj T2'
2998	CASE (2) ; cl_name = 'bn2_adj T3'
2999	END SELECT
3000	CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
3001	ENDDO
3002	! restore initial nn_eos
3003	nn_eos = jeos
3004
3005	! Deallocate memory
3006	DEALLOCATE( &
3007	& zts, &
3008	& zts_adout, &
3009	& zrd_adin, &
3010	& zrd_adout, &
3011	& zts_tlin, &
3012	& zts_tlout, &
3013	& zrd_tlout, &
3014	& znts )
3015	END SUBROUTINE bn2_adj_tst
3016	#endif
3017	!!======================================================================
3018	END MODULE eosbn2_tam

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

Download in other formats: