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 @ 3660

Last change on this file since 3660 was 3640, checked in by pabouttier, 12 years ago
Missing allocation/deallocation in TAM routines - See ticket #1013
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	CALL wrk_dealloc( jpi, jpj, zws )
772	!
773	END SUBROUTINE eos_insitu_2d_tan
774
775	SUBROUTINE eos_insitu_adj(pts, pts_ad, prd_ad)
776	!!-----------------------------------------------------------------------
777	!!
778	!! * ROUTINE eos_insitu_tan : Adjoint OF ROUTINE eos_insitu *
779	!!
780	!! ** Purpose of direct routine : Compute the in situ density
781	!! (ratio rho/rau0) from potential temperature and salinity
782	!! using an equation of state defined through the namelist
783	!! parameter nneos.
784	!!
785	!! ** Method of direct routine : 3 cases:
786	!! nn_eos = 0 : Jackett and McDougall (1994) equation of state.
787	!! the in situ density is computed directly as a function of
788	!! potential temperature relative to the surface (the opa t
789	!! variable), salt and pressure (assuming no pressure variation
790	!! along geopotential surfaces, i.e. the pressure p in decibars
791	!! is approximated by the depth in meters.
792	!! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0
793	!! with pressure p decibars
794	!! potential temperature t deg celsius
795	!! salinity s psu
796	!! reference volumic mass rau0 kg/m**3
797	!! in situ volumic mass rho kg/m**3
798	!! in situ density anomalie prd no units
799	!! Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar,
800	!! t = 40 deg celcius, s=40 psu
801	!! nn_eos = 1 : linear equation of state function of temperature only
802	!! prd(t) = 0.0285 - rn_alpha * t
803	!! nn_eos = 2 : linear equation of state function of temperature and
804	!! salinity
805	!! prd(t,s) = rn_beta * s - rn_alpha * tn - 1.
806	!! Note that no boundary condition problem occurs in this routine
807	!! as (ptem,psal) are defined over the whole domain.
808	!!
809	!! ** Comments on Adjoint Routine :
810	!! Care has been taken to avoid division by zero when computing
811	!! the inverse of the square root of salinity at masked salinity
812	!! points.
813	!!
814	!! ** Action :
815	!!
816	!! References :
817	!!
818	!! History :
819	!! 8.2 ! 05-03 ((F. Van den Berghe, A. Weaver, N. Daget) - eostan.F
820	!! 9.0 ! 08-08 (A. Vidard) 9.0 version
821	!!-----------------------------------------------------------------------
822	!! * Modules used
823	!! * Arguments
824	REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pts ! 1 : potential temperature [Celcius]
825	! ! 2 : salinity [psu]
826	REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pts_ad ! 1 : TL of potential temperature [Celsius]
827	! 2 : TL of salinity [psu]
828	REAL(wp), DIMENSION(:,:,:), INTENT( inout ) :: &
829	& prd_ad ! TL of potential density (surface referenced)
830
831	!! * Local declarations
832	INTEGER :: ji, jj, jk ! dummy loop indices
833	REAL(wp) :: & ! temporary scalars
834	zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw, &
835	zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, &
836	ztad, zsad, zhad, zsrad, zr1ad, zr2ad, zr3ad, &
837	zr4ad, zrhopad, zead, zbwad, &
838	zbad, zdad, zcad, zawad, zaad, zb1ad, za1ad, &
839	zkwad, zk0ad, zpes, zrdc1, zrdc2, zeps, &
840	zmask, zrau0r
841	REAL(wp), POINTER, DIMENSION(:,:,:) :: zws
842	!!----------------------------------------------------------------------
843	IF( nn_timing == 1 ) CALL timing_start('eos_adj')
844	!
845	CALL wrk_alloc( jpi, jpj, jpk, zws )
846	!
847	! initialization of adjoint variables
848	ztad = 0.0_wp
849	zsad = 0.0_wp
850	zhad = 0.0_wp
851	zsrad = 0.0_wp
852	zr1ad = 0.0_wp
853	zr2ad = 0.0_wp
854	zr3ad = 0.0_wp
855	zr4ad = 0.0_wp
856	zrhopad = 0.0_wp
857	zead = 0.0_wp
858	zbwad = 0.0_wp
859	zbad = 0.0_wp
860	zdad = 0.0_wp
861	zcad = 0.0_wp
862	zawad = 0.0_wp
863	zaad = 0.0_wp
864	zb1ad = 0.0_wp
865	za1ad = 0.0_wp
866	zkwad = 0.0_wp
867	zk0ad = 0.0_wp
868	SELECT CASE ( nn_eos )
869
870	CASE ( 0 ) !== Jackett and McDougall (1994) formulation ==!
871	zrau0r = 1.e0 / rau0
872	#ifdef key_sp
873	zeps = 1.e-7
874	#else
875	zeps = 1.e-14
876	#endif
877	!CDIR NOVERRCHK
878	zws(:,:,:) = SQRT( ABS( pts(:,:,:, jp_sal) ) )
879	DO jk = jpkm1, 1, -1
880	DO jj = jpj, 1, -1
881	DO ji = jpi, 1, -1
882	zt = pts(ji,jj,jk, jp_tem)
883	zs = pts(ji,jj,jk, jp_sal)
884	zh = fsdept(ji,jj,jk) ! depth
885	zsr= zws(ji,jj,jk) ! square root salinity
886	! compute volumic mass pure water at atm pressure
887	zr1= ( ( ( ( 6.536332e-9_wpzt-1.120083e-6_wp )zt+1.001685e-4_wp)*zt &
888	-9.095290e-3_wp )zt+6.793952e-2_wp )zt+999.842594_wp
889	! seawater volumic mass atm pressure
890	zr2= ( ( ( 5.3875e-9_wpzt-8.2467e-7_wp ) zt+7.6438e-5_wp ) *zt &
891	-4.0899e-3_wp ) *zt+0.824493_wp
892	zr3= ( -1.6546e-6_wpzt+1.0227e-4_wp ) zt-5.72466e-3_wp
893	zr4= 4.8314e-4_wp
894
895	! potential volumic mass (reference to the surface)
896	zrhop= ( zr4zs + zr3zsr + zr2 ) *zs + zr1
897
898	! add the compression terms
899	ze = ( -3.508914e-8_wpzt-1.248266e-8_wp ) zt-2.595994e-6_wp
900	zbw= ( 1.296821e-6_wpzt-5.782165e-9_wp ) zt+1.045941e-4_wp
901	zb = zbw + ze * zs
902
903	zd = -2.042967e-2_wp
904	zc = (-7.267926e-5_wpzt+2.598241e-3_wp ) zt+0.1571896_wp
905	zaw= ( ( 5.939910e-6_wpzt+2.512549e-3_wp ) zt-0.1028859_wp ) *zt - 4.721788_wp
906	za = ( zdzsr + zc ) zs + zaw
907
908	zb1= (-0.1909078_wpzt+7.390729_wp ) zt-55.87545_wp
909	za1= ( ( 2.326469e-3_wpzt+1.553190_wp)zt-65.00517_wp ) *zt+1044.077_wp
910	zkw= ( ( (-1.361629e-4_wpzt-1.852732e-2_wp ) zt-30.41638_wp ) zt + 2098.925_wp ) zt+190925.6_wp
911	zk0= ( zb1zsr + za1 )zs + zkw
912
913	zrdc1 = 1.0 / ( zk0 - zh * ( za - zh * zb ) )
914	zrdc2 = 1.0 / ( 1.0 - zh * zrdc1 )
915	! ============
916	! Adjoint part
917	! ============
918
919	! Masked in situ density anomaly
920
921	zrhopad = zrhopad + prd_ad(ji,jj,jk) * tmask(ji,jj,jk) &
922	& * zrdc2 * zrau0r
923	zk0ad = zk0ad - prd_ad(ji,jj,jk) * tmask(ji,jj,jk) &
924	& * zrdc2 * zrdc2 * zh &
925	& * zrdc1*2 zrhop &
926	& * zrau0r
927	zaad = zaad + prd_ad(ji,jj,jk) * tmask(ji,jj,jk) &
928	& * zrdc2 * zrdc2 * zh &
929	& * zrdc1*2 zrhop &
930	& * zh * zrau0r
931	zbad = zbad - prd_ad(ji,jj,jk) * tmask(ji,jj,jk) &
932	& * zrdc2 * zrdc2 * zh &
933	& * zrdc1*2 zrhop &
934	& * zh * zh * zrau0r
935	prd_ad(ji,jj,jk) = 0.0_wp
936
937	zkwad = zkwad + zk0ad
938	zsrad = zsrad + zk0ad * zb1 * zs
939	zb1ad = zb1ad + zk0ad * zs * zsr
940	za1ad = za1ad + zk0ad * zs
941	zsad = zsad + zk0ad * ( zb1 * zsr + za1 )
942	zk0ad = 0.0_wp
943
944	ztad = ztad + zkwad * ( ( (-4.1.361629e-4_wp zt &
945	& -3.1.852732e-2_wp ) zt &
946	& -2.30.41638_wp ) zt &
947	& + 2098.925_wp )
948	zkwad = 0.0_wp
949
950	ztad = ztad + za1ad * ( ( 3.2.326469e-3_wp zt &
951	& +2.1.553190_wp ) zt &
952	& - 65.00517_wp )
953	za1ad = 0.0_wp
954
955	ztad = ztad + zb1ad * (-2.0.1909078_wp zt &
956	& + 7.390729_wp )
957	zb1ad = 0.0_wp
958
959	zawad = zawad + zaad
960	zsrad = zsrad + zaad * zd * zs
961	zcad = zcad + zaad * zs
962	zsad = zsad + zaad * ( zd * zsr + zc )
963	zaad = 0.0_wp
964
965	ztad = ztad + zawad * ( ( 3.5.939910e-6_wp zt &
966	& +2.2.512549e-3_wp ) zt &
967	& - 0.1028859_wp )
968	zawad = 0.0_wp
969
970	ztad = ztad + zcad * (-2.7.267926e-5_wp zt &
971	& + 2.598241e-3_wp )
972	zcad = 0.0_wp
973
974	zbwad = zbwad + zbad
975	zead = zead + zbad * zs
976	zsad = zsad + zbad * ze
977	zbad = 0.0_wp
978
979	ztad = ztad + zbwad * ( 2.1.296821e-6_wp zt &
980	& - 5.782165e-9_wp )
981	zbwad = 0.0_wp
982
983	ztad = ztad + zead * (-2.3.508914e-8_wp zt &
984	& - 1.248266e-8_wp )
985	zead = 0.0_wp
986
987	zr1ad = zr1ad + zrhopad
988	zr2ad = zr2ad + zrhopad * zs
989	zr3ad = zr3ad + zrhopad * zsr * zs
990	zsrad = zsrad + zrhopad * zr3 * zs
991	zsad = zsad + zrhopad * ( 2. * zr4 * zs + zr2 &
992	& + zr3 * zsr )
993	zrhopad = 0.0_wp
994
995	ztad = ztad + zr3ad * (-2.1.6546e-6_wp zt &
996	& + 1.0227e-4_wp )
997	zr3ad = 0.0_wp
998
999	ztad = ztad + zr2ad * ( ( ( 4.5.3875e-9_wp zt &
1000	& -3.8.2467e-7_wp ) zt &
1001	& +2.7.6438e-5_wp ) zt &
1002	& - 4.0899e-3_wp )
1003	zr2ad = 0.0_wp
1004
1005	ztad = ztad + zr1ad * ( ( ( ( 5.6.536332e-9_wp zt &
1006	& -4.1.120083e-6_wp ) zt &
1007	& +3.1.001685e-4_wp ) zt &
1008	& -2.9.095290e-3_wp ) zt &
1009	& + 6.793952e-2_wp )
1010	zr1ad = 0.0_wp
1011
1012	zsad = zsad + zsrad * ( 1.0 / MAX( 2.*zsr, zeps ) ) &
1013	& * tmask(ji,jj,jk)
1014	zsrad = 0.0_wp
1015
1016	pts_ad(ji,jj,jk, jp_sal) = pts_ad(ji,jj,jk, jp_sal) + zsad
1017	pts_ad(ji,jj,jk,jp_tem) = pts_ad(ji,jj,jk, jp_tem) + ztad
1018	ztad = 0.0_wp
1019	zsad = 0.0_wp
1020	END DO
1021	END DO
1022	END DO
1023	!
1024	CASE ( 1 ) !== Linear formulation function of temperature only ==!
1025	DO jk = jpkm1, 1, -1
1026	pts_ad(:,:,jk,jp_tem) = pts_ad(:,:,jk,jp_tem) - rn_alpha * prd_ad(:,:,jk) * tmask(:,:,jk)
1027	prd_ad(:,:,jk) = 0.0_wp
1028	END DO
1029	!
1030	CASE ( 2 ) !== Linear formulation function of temperature and salinity ==!
1031	DO jk = jpkm1, 1, -1
1032	pts_ad(:,:,jk,jp_tem) = pts_ad(:,:,jk,jp_tem) - rn_alpha * prd_ad(:,:,jk) * tmask(:,:,jk)
1033	pts_ad(:,:,jk,jp_sal) = pts_ad(:,:,jk,jp_sal) + rn_beta * prd_ad( :,:,jk) * tmask(:,:,jk)
1034	prd_ad( :,:,jk) = 0.0_wp
1035	END DO
1036	!
1037	END SELECT
1038	!
1039	CALL wrk_dealloc( jpi, jpj, jpk, zws )
1040	!
1041	IF( nn_timing == 1 ) CALL timing_stop('eos_adj')
1042	!
1043	END SUBROUTINE eos_insitu_adj
1044
1045	SUBROUTINE eos_insitu_pot_adj ( pts, pts_ad, prd_ad, prhop_ad )
1046	!!----------------------------------------------------------------------
1047	!! * ROUTINE eos_insitu_pot_adj *
1048	!!
1049	!! ** Purpose or the direct routine:
1050	!! Compute the in situ density (ratio rho/rau0) and the
1051	!! potential volumic mass (Kg/m3) from potential temperature and
1052	!! salinity fields using an equation of state defined through the
1053	!! namelist parameter nn_eos.
1054	!!
1055	!! ** Method :
1056	!! nn_eos = 0 : Jackett and McDougall (1994) equation of state.
1057	!! the in situ density is computed directly as a function of
1058	!! potential temperature relative to the surface (the opa t
1059	!! variable), salt and pressure (assuming no pressure variation
1060	!! along geopotential surfaces, i.e. the pressure p in decibars
1061	!! is approximated by the depth in meters.
1062	!! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0
1063	!! rhop(t,s) = rho(t,s,0)
1064	!! with pressure p decibars
1065	!! potential temperature t deg celsius
1066	!! salinity s psu
1067	!! reference volumic mass rau0 kg/m**3
1068	!! in situ volumic mass rho kg/m**3
1069	!! in situ density anomalie prd no units
1070	!!
1071	!! Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar,
1072	!! t = 40 deg celcius, s=40 psu
1073	!!
1074	!! neos = 1 : linear equation of state function of temperature only
1075	!! prd(t) = ( rho(t) - rau0 ) / rau0 = 0.028 - ralpha * t
1076	!! rhop(t,s) = rho(t,s)
1077	!!
1078	!! nn_eos = 2 : linear equation of state function of temperature and
1079	!! salinity
1080	!! prd(t,s) = ( rho(t,s) - rau0 ) / rau0
1081	!! = rn_beta * s - rn_alpha * tn - 1.
1082	!! rhop(t,s) = rho(t,s)
1083	!! Note that no boundary condition problem occurs in this routine
1084	!! as (tn,sn) or (ta,sa) are defined over the whole domain.
1085	!!
1086	!! ** Action : - prd , the in situ density (no units)
1087	!! - prhop, the potential volumic mass (Kg/m3)
1088	!!
1089	!! References :
1090	!! Jackett, D.R., and T.J. McDougall. J. Atmos. Ocean. Tech., 1994
1091	!! Brown, J. A. and K. A. Campana. Mon. Weather Rev., 1978
1092	!!
1093	!! History of the adjoint routine:
1094	!! 9.0 ! 08-06 (A. Vidard) Initial version
1095	!!----------------------------------------------------------------------
1096	!! * Arguments
1097
1098	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! 1 : potential temperature/salinity [Celcius/psu]
1099	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(inout) :: pts_ad ! 1 : potential temperature/salinity [Celcius/psu]
1100	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(inout) :: prd_ad ! TL of in_situ density [-]
1101	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(inout) :: prhop_ad ! TL of potential density (surface referenced)
1102	!! * Local declarations
1103	INTEGER :: ji, jj, jk ! dummy loop indices
1104	REAL(wp) :: & ! temporary scalars
1105	zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw, &
1106	zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, &
1107	ztad, zsad, zhad, zsrad, zr1ad, zr2ad, zr3ad, &
1108	zr4ad, zrhopad, zead, zbwad, &
1109	zbad, zdad, zcad, zawad, zaad, zb1ad, za1ad, &
1110	zkwad, zk0ad, zpes, zrdc1, zrdc2, zeps, &
1111	zmask, zrau0r
1112	REAL(wp), POINTER, DIMENSION(:,:,:) :: zws
1113	!!----------------------------------------------------------------------
1114	!
1115	IF( nn_timing == 1 ) CALL timing_start('eos-p_adj')
1116	!
1117	CALL wrk_alloc( jpi, jpj, jpk, zws )
1118	!
1119	! initialization of adjoint variables
1120	ztad = 0.0_wp
1121	zsad = 0.0_wp
1122	zhad = 0.0_wp
1123	zsrad = 0.0_wp
1124	zr1ad = 0.0_wp
1125	zr2ad = 0.0_wp
1126	zr3ad = 0.0_wp
1127	zr4ad = 0.0_wp
1128	zrhopad = 0.0_wp
1129	zead = 0.0_wp
1130	zbwad = 0.0_wp
1131	zbad = 0.0_wp
1132	zdad = 0.0_wp
1133	zcad = 0.0_wp
1134	zawad = 0.0_wp
1135	zaad = 0.0_wp
1136	zb1ad = 0.0_wp
1137	za1ad = 0.0_wp
1138	zkwad = 0.0_wp
1139	zk0ad = 0.0_wp
1140
1141	SELECT CASE ( nn_eos )
1142
1143	CASE ( 0 ) !== Jackett and McDougall (1994) formulation ==!
1144	zrau0r = 1.e0 / rau0
1145	#ifdef key_sp
1146	zeps = 1.e-7
1147	#else
1148	zeps = 1.e-14
1149	#endif
1150	!CDIR NOVERRCHK
1151	zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) )
1152	!
1153	DO jk = jpkm1, 1, -1
1154	DO jj = jpj, 1, -1
1155	DO ji = jpi, 1, -1
1156	! direct recomputing
1157	zt = pts(ji,jj,jk,jp_tem)
1158	zs = pts(ji,jj,jk,jp_sal)
1159	zh = fsdept(ji,jj,jk) ! depth
1160	zsr = zws(ji,jj,jk) ! square root salinity
1161	! compute volumic mass pure water at atm pressure
1162	zr1 = ( ( ( ( 6.536332e-9_wp * zt - 1.120083e-6_wp ) * zt &
1163	& + 1.001685e-4_wp ) * zt - 9.095290e-3_wp ) * zt &
1164	& + 6.793952e-2_wp ) * zt + 999.842594_wp
1165	! seawater volumic mass atm pressure
1166	zr2 = ( ( ( 5.3875e-9_wp * zt - 8.2467e-7_wp ) * zt &
1167	& + 7.6438e-5_wp ) * zt - 4.0899e-3_wp ) * zt + 0.824493_wp
1168	zr3 = ( -1.6546e-6_wp * zt + 1.0227e-4_wp ) * zt - 5.72466e-3_wp
1169	zr4 = 4.8314e-4_wp
1170	! potential volumic mass (reference to the surface)
1171	zrhop = ( zr4 * zs + zr3zsr + zr2 ) zs + zr1
1172	! add the compression terms
1173	ze = ( -3.508914e-8_wp * zt - 1.248266e-8_wp ) * zt - 2.595994e-6_wp
1174	zbw = ( 1.296821e-6_wp * zt - 5.782165e-9_wp ) * zt + 1.045941e-4_wp
1175	zb = zbw + ze * zs
1176
1177	zd = -2.042967e-2_wp
1178	zc = (-7.267926e-5_wp * zt + 2.598241e-3_wp ) * zt + 0.1571896_wp
1179	zaw= ( ( 5.939910e-6_wp * zt + 2.512549e-3_wp ) * zt - 0.1028859_wp &
1180	& ) * zt - 4.721788
1181	za = ( zd * zsr + zc ) * zs + zaw
1182
1183	zb1= (-0.1909078_wp * zt + 7.390729_wp ) * zt - 55.87545_wp
1184	za1= ( ( 2.326469e-3_wp * zt + 1.553190_wp ) * zt - 65.00517_wp &
1185	& ) * zt + 1044.077_wp
1186	zkw= ( ( (-1.361629e-4_wp * zt - 1.852732e-2_wp ) * zt - 30.41638_wp &
1187	& ) * zt + 2098.925_wp ) * zt + 190925.6_wp
1188	zk0= ( zb1 * zsr + za1 ) * zs + zkw
1189
1190
1191	zrdc1 = 1.0 / ( zk0 - zh * ( za - zh * zb ) )
1192	zrdc2 = 1.0 / ( 1.0 - zh * zrdc1 )
1193
1194	! ============
1195	! Adjoint part
1196	! ============
1197
1198	! Masked in situ density anomaly
1199
1200	zrhopad = zrhopad + prd_ad(ji,jj,jk) * tmask(ji,jj,jk) &
1201	& * zrdc2 * zrau0r
1202	zk0ad = zk0ad - prd_ad(ji,jj,jk) * tmask(ji,jj,jk) &
1203	& * zrdc2 * zrdc2 * zh &
1204	& * zrdc1*2 zrhop &
1205	& * zrau0r
1206	zaad = zaad + prd_ad(ji,jj,jk) * tmask(ji,jj,jk) &
1207	& * zrdc2 * zrdc2 * zh &
1208	& * zrdc1*2 zrhop &
1209	& * zh * zrau0r
1210	zbad = zbad - prd_ad(ji,jj,jk) * tmask(ji,jj,jk) &
1211	& * zrdc2 * zrdc2 * zh &
1212	& * zrdc1*2 zrhop &
1213	& * zh * zh * zrau0r
1214	prd_ad(ji,jj,jk) = 0.0_wp
1215
1216	zkwad = zkwad + zk0ad
1217	zsrad = zsrad + zk0ad * zb1 * zs
1218	zb1ad = zb1ad + zk0ad * zs * zsr
1219	za1ad = za1ad + zk0ad * zs
1220	zsad = zsad + zk0ad * ( zb1 * zsr + za1 )
1221	zk0ad = 0.0_wp
1222
1223	ztad = ztad + zkwad * ( ( (-4.1.361629e-4_wp zt &
1224	& -3.1.852732e-2_wp ) zt &
1225	& -2.30.41638_wp ) zt &
1226	& + 2098.925_wp )
1227	zkwad = 0.0_wp
1228
1229	ztad = ztad + za1ad * ( ( 3.2.326469e-3_wp zt &
1230	& +2.1.553190_wp ) zt &
1231	& - 65.00517_wp )
1232	za1ad = 0.0_wp
1233
1234	ztad = ztad + zb1ad * (-2.0.1909078_wp zt &
1235	& + 7.390729_wp )
1236	zb1ad = 0.0_wp
1237
1238	zawad = zawad + zaad
1239	zsrad = zsrad + zaad * zd * zs
1240	zcad = zcad + zaad * zs
1241	zsad = zsad + zaad * ( zd * zsr + zc )
1242	zaad = 0.0_wp
1243
1244	ztad = ztad + zawad * ( ( 3.5.939910e-6_wp zt &
1245	& +2.2.512549e-3_wp ) zt &
1246	& - 0.1028859_wp )
1247	zawad = 0.0_wp
1248
1249	ztad = ztad + zcad * (-2.7.267926e-5_wp zt &
1250	& + 2.598241e-3_wp )
1251	zcad = 0.0_wp
1252
1253
1254	zsad = zsad + zbad * ze
1255	zead = zead + zbad * zs
1256	zbwad = zbwad + zbad
1257	zbad = 0.0_wp
1258
1259	ztad = ztad + zbwad * ( 2.1.296821e-6_wp zt &
1260	& - 5.782165e-9_wp )
1261	zbwad = 0.0_wp
1262
1263	ztad = ztad + zead * (-2.3.508914e-8_wp zt &
1264	& - 1.248266e-8_wp )
1265	zead = 0.0_wp
1266
1267	zrhopad = zrhopad + prhop_ad(ji,jj,jk) * tmask(ji,jj,jk)
1268	prhop_ad(ji,jj,jk) = 0.0_wp
1269
1270	zr1ad = zr1ad + zrhopad
1271	zr2ad = zr2ad + zrhopad * zs
1272	zr3ad = zr3ad + zrhopad * zsr * zs
1273	zsrad = zsrad + zrhopad * zr3 * zs
1274	zsad = zsad + zrhopad * ( 2. * zr4 * zs + zr2 &
1275	& + zr3 * zsr )
1276	zrhopad = 0.0_wp
1277
1278	ztad = ztad + zr3ad * (-2.1.6546e-6_wp zt &
1279	& + 1.0227e-4_wp )
1280	zr3ad = 0.0_wp
1281
1282	ztad = ztad + zr2ad * ( ( ( 4.5.3875e-9_wp zt &
1283	& -3.8.2467e-7_wp ) zt &
1284	& +2.7.6438e-5_wp ) zt &
1285	& - 4.0899e-3_wp )
1286	zr2ad = 0.0_wp
1287
1288	ztad = ztad + zr1ad * ( ( ( ( 5.6.536332e-9_wp zt &
1289	& -4.1.120083e-6_wp ) zt &
1290	& +3.1.001685e-4_wp ) zt &
1291	& -2.9.095290e-3_wp ) zt &
1292	& + 6.793952e-2_wp )
1293	zr1ad = 0.0_wp
1294
1295	zsad = zsad + zsrad * ( 1.0 / MAX( 2.*zsr, zeps ) ) &
1296	& * tmask(ji,jj,jk)
1297	zsrad = 0.0_wp
1298
1299	pts_ad(ji,jj,jk,jp_sal) = pts_ad(ji,jj,jk,jp_sal) + zsad
1300	pts_ad(ji,jj,jk,jp_tem) = pts_ad(ji,jj,jk,jp_tem) + ztad
1301	ztad = 0.0_wp
1302	zsad = 0.0_wp
1303	END DO
1304	END DO
1305	END DO
1306	!
1307	CASE ( 1 ) !== Linear formulation = F( temperature ) ==!
1308	DO jk = jpkm1, 1, -1
1309	prd_ad(:,:,jk) = prd_ad(:,:,jk) + rau0 * prhop_ad(:,:,jk) * tmask(:,:,jk)
1310	prhop_ad(:,:,jk) = 0.0_wp
1311	pts_ad(:,:,jk,jp_tem) = pts_ad(:,:,jk,jp_tem) - rn_alpha * prd_ad(:,:,jk) * tmask(:,:,jk)
1312	prd_ad(:,:,jk) = 0.0_wp
1313	END DO
1314	!
1315	CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==!
1316	DO jk = jpkm1, 1, -1
1317	prd_ad( :,:,jk) = prd_ad(:,:,jk) + rau0 * prhop_ad(:,:,jk) * tmask(:,:,jk)
1318	prhop_ad(:,:,jk) = 0.0_wp
1319	pts_ad( :,:,jk,jp_tem) = pts_ad(:,:,jk,jp_tem) - rn_alpha * prd_ad(:,:,jk) * tmask(:,:,jk)
1320	pts_ad( :,:,jk,jp_sal) = pts_ad(:,:,jk,jp_sal) + rn_beta * prd_ad(:,:,jk) * tmask(:,:,jk)
1321	prd_ad( :,:,jk) = 0.0_wp
1322	END DO
1323	!
1324	END SELECT
1325	CALL wrk_dealloc( jpi, jpj, jpk, zws )
1326	!
1327	IF( nn_timing == 1 ) CALL timing_stop('eos-p_adj')
1328	!
1329	END SUBROUTINE eos_insitu_pot_adj
1330
1331	SUBROUTINE eos_insitu_2d_adj( pts, pdep, pts_ad, prd_ad )
1332	!!-----------------------------------------------------------------------
1333	!!
1334	!! * ROUTINE eos_insitu_2d_adj : adj OF ROUTINE eos_insitu_2d *
1335	!!
1336	!! ** Purpose of direct routine : Compute the in situ density
1337	!! (ratio rho/rau0) from potential temperature and salinity
1338	!! using an equation of state defined through the namelist
1339	!! parameter nn_eos. * 2D field case
1340	!!
1341	!! ** Method of direct routine : 3 cases:
1342	!! nn_eos = 0 : Jackett and McDougall (1994) equation of state.
1343	!! the in situ density is computed directly as a function of
1344	!! potential temperature relative to the surface (the opa t
1345	!! variable), salt and pressure (assuming no pressure variation
1346	!! along geopotential surfaces, i.e. the pressure p in decibars
1347	!! is approximated by the depth in meters.
1348	!! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0
1349	!! with pressure p decibars
1350	!! potential temperature t deg celsius
1351	!! salinity s psu
1352	!! reference volumic mass rau0 kg/m**3
1353	!! in situ volumic mass rho kg/m**3
1354	!! in situ density anomalie prd no units
1355	!! Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar,
1356	!! t = 40 deg celcius, s=40 psu
1357	!! nn_eos = 1 : linear equation of state function of temperature only
1358	!! prd(t) = 0.0285 - rn_alpha * t
1359	!! nn_eos = 2 : linear equation of state function of temperature and
1360	!! salinity
1361	!! prd(t,s) = rn_beta * s - rn_alpha * tn - 1.
1362	!! Note that no boundary condition problem occurs in this routine
1363	!! as (ptem,psal) are defined over the whole domain.
1364	!!
1365	!! ** Comments on Adjoint Routine :
1366	!! Care has been taken to avoid division by zero when computing
1367	!! the inverse of the square root of salinity at masked salinity
1368	!! points.
1369	!!
1370	!! ** Action :
1371	!!
1372	!! References :
1373	!!
1374	!! History :
1375	!! 8.2 ! 05-03 ((F. Van den Berghe, A. Weaver, N. Daget) - eosadj.F
1376	!! 9.0 ! 08-07 (A. Vidard) first version based on eosadj
1377	!!-----------------------------------------------------------------------
1378	!! * Modules used
1379	!! * Arguments
1380	REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in ) :: pts ! 1 : potential temperature [Celcius]
1381	! ! 2 : salinity [psu]
1382	REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(inout) :: pts_ad ! 1 : TL of potential temperature [Celcius]
1383	! ! 2 : TL of salinity [psu]
1384	REAL(wp), DIMENSION(jpi,jpj ), INTENT( inout) :: prd_ad ! TL of in_situ density [-]
1385	REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pdep ! depth [m]
1386	!
1387
1388	INTEGER :: ji, jj ! dummy loop indices
1389	REAL(wp) :: & ! temporary scalars
1390	zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw, &
1391	zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, &
1392	ztad, zsad, zhad, zsrad, zr1ad, zr2ad, zr3ad, &
1393	zr4ad, zrhopad, zead, zbwad, &
1394	zbad, zdad, zcad, zawad, zaad, zb1ad, za1ad, &
1395	zkwad, zk0ad, zpes, zrdc1, zrdc2, zeps, &
1396	zmask
1397	REAL(wp), POINTER, DIMENSION(:,:) :: zws
1398	!!----------------------------------------------------------------------
1399	!
1400	IF( nn_timing == 1 ) CALL timing_start('eos2d_adj')
1401	!
1402	CALL wrk_alloc( jpi, jpj, zws )
1403	!
1404	! initialization of adjoint variables
1405	ztad = 0.0_wp
1406	zsad = 0.0_wp
1407	zhad = 0.0_wp
1408	zsrad = 0.0_wp
1409	zr1ad = 0.0_wp
1410	zr2ad = 0.0_wp
1411	zr3ad = 0.0_wp
1412	zr4ad = 0.0_wp
1413	zrhopad = 0.0_wp
1414	zead = 0.0_wp
1415	zbwad = 0.0_wp
1416	zbad = 0.0_wp
1417	zdad = 0.0_wp
1418	zcad = 0.0_wp
1419	zawad = 0.0_wp
1420	zaad = 0.0_wp
1421	zb1ad = 0.0_wp
1422	za1ad = 0.0_wp
1423	zkwad = 0.0_wp
1424	zk0ad = 0.0_wp
1425	SELECT CASE ( nn_eos )
1426	CASE ( 0 ) !== Jackett and McDougall (1994) formulation ==!
1427
1428	#ifdef key_sp
1429	zeps = 1.e-7
1430	#else
1431	zeps = 1.e-14
1432	#endif
1433	!CDIR NOVERRCHK
1434	DO jj = jpjm1, 1, -1
1435	!CDIR NOVERRCHK
1436	DO ji = fs_jpim1, 1, -1 ! vector opt.
1437	zws(ji,jj) = SQRT( ABS( pts(ji,jj,jp_sal) ) )
1438	END DO
1439	END DO
1440	!
1441	DO jj = jpjm1, 1, -1
1442	DO ji = fs_jpim1, 1, -1 ! vector opt.
1443	zmask = tmask(ji,jj,1) ! land/sea bottom mask = surf. mask
1444	zt = pts (ji,jj,jp_tem) ! interpolated T
1445	zs = pts (ji,jj,jp_sal) ! interpolated S
1446	zsr= zws(ji,jj) ! square root of interpolated S
1447	zh = pdep(ji,jj) ! depth at the partial step level
1448	! compute volumic mass pure water at atm pressure
1449	zr1= ( ( ( ( 6.536332e-9_wpzt-1.120083e-6_wp )zt+1.001685e-4_wp)*zt &
1450	& -9.095290e-3_wp )zt+6.793952e-2_wp )zt+999.842594_wp
1451	! seawater volumic mass atm pressure
1452	zr2= ( ( ( 5.3875e-9_wpzt-8.2467e-7_wp ) zt+7.6438e-5_wp ) *zt &
1453	& -4.0899e-3_wp ) *zt+0.824493_wp
1454	zr3= ( -1.6546e-6_wpzt+1.0227e-4_wp ) zt-5.72466e-3_wp
1455	zr4= 4.8314e-4_wp
1456
1457	! potential volumic mass (reference to the surface)
1458	zrhop= ( zr4zs + zr3zsr + zr2 ) *zs + zr1
1459
1460	! add the compression terms
1461	ze = ( -3.508914e-8_wpzt-1.248266e-8_wp ) zt-2.595994e-6_wp
1462	zbw= ( 1.296821e-6_wpzt-5.782165e-9_wp ) zt+1.045941e-4_wp
1463	zb = zbw + ze * zs
1464
1465	zd = -2.042967e-2_wp
1466	zc = (-7.267926e-5_wpzt+2.598241e-3_wp ) zt+0.1571896_wp
1467	zaw= ( ( 5.939910e-6_wpzt+2.512549e-3_wp ) zt-0.1028859_wp ) *zt - 4.721788_wp
1468	za = ( zdzsr + zc ) zs + zaw
1469
1470	zb1= (-0.1909078_wpzt+7.390729_wp ) zt-55.87545_wp
1471	za1= ( ( 2.326469e-3_wpzt+1.553190_wp)zt-65.00517_wp ) *zt+1044.077_wp
1472	zkw= ( ( (-1.361629e-4_wpzt-1.852732e-2_wp ) zt-30.41638_wp ) zt + 2098.925_wp ) zt+190925.6_wp
1473	zk0= ( zb1zsr + za1 )zs + zkw
1474
1475	zrdc1 = 1.0 / ( zk0 - zh * ( za - zh * zb ) )
1476	zrdc2 = 1.0 / ( 1.0 - zh * zrdc1 )
1477	! ============
1478	! Adjoint part
1479	! ============
1480
1481	! Masked in situ density anomaly
1482
1483	zrhopad = zrhopad + prd_ad(ji,jj) * tmask(ji,jj,1) &
1484	& * zrdc2 / rau0
1485	zk0ad = zk0ad - prd_ad(ji,jj) * tmask(ji,jj,1) &
1486	& * zrdc2 * zrdc2 * zh &
1487	& * zrdc1*2 zrhop &
1488	& / rau0
1489	zaad = zaad + prd_ad(ji,jj) * tmask(ji,jj,1) &
1490	& * zrdc2 * zrdc2 * zh &
1491	& * zrdc1*2 zrhop &
1492	& * zh / rau0
1493	zbad = zbad - prd_ad(ji,jj) * tmask(ji,jj,1) &
1494	& * zrdc2 * zrdc2 * zh &
1495	& * zrdc1*2 zrhop &
1496	& * zh * zh / rau0
1497	prd_ad(ji,jj) = 0.0_wp
1498
1499	zkwad = zkwad + zk0ad
1500	zsrad = zsrad + zk0ad * zb1 * zs
1501	zb1ad = zb1ad + zk0ad * zs * zsr
1502	za1ad = za1ad + zk0ad * zs
1503	zsad = zsad + zk0ad * ( zb1 * zsr + za1 )
1504	zk0ad = 0.0_wp
1505
1506	ztad = ztad + zkwad * ( ( (-4.1.361629e-4_wp zt &
1507	& -3.1.852732e-2_wp ) zt &
1508	& -2.30.41638_wp ) zt &
1509	& + 2098.925_wp )
1510	zkwad = 0.0_wp
1511
1512	ztad = ztad + za1ad * ( ( 3.2.326469e-3_wp zt &
1513	& +2.1.553190_wp ) zt &
1514	& - 65.00517_wp )
1515	za1ad = 0.0_wp
1516
1517	ztad = ztad + zb1ad * (-2.0.1909078_wp zt &
1518	& + 7.390729_wp )
1519	zb1ad = 0.0_wp
1520
1521	zawad = zawad + zaad
1522	zsrad = zsrad + zaad * zd * zs
1523	zcad = zcad + zaad * zs
1524	zsad = zsad + zaad * ( zd * zsr + zc )
1525	zaad = 0.0_wp
1526
1527	ztad = ztad + zawad * ( ( 3.5.939910e-6_wp zt &
1528	& +2.2.512549e-3_wp ) zt &
1529	& - 0.1028859_wp )
1530	zawad = 0.0_wp
1531
1532	ztad = ztad + zcad * (-2.7.267926e-5_wp zt &
1533	& + 2.598241e-3_wp )
1534	zcad = 0.0_wp
1535
1536	zbwad = zbwad + zbad
1537	zead = zead + zbad * zs
1538	zsad = zsad + zbad * ze
1539	zbad = 0.0_wp
1540
1541	ztad = ztad + zbwad * ( 2.1.296821e-6_wp zt &
1542	& - 5.782165e-9_wp )
1543	zbwad = 0.0_wp
1544
1545	ztad = ztad + zead * (-2.3.508914e-8_wp zt &
1546	& - 1.248266e-8_wp )
1547	zead = 0.0_wp
1548
1549	zr1ad = zr1ad + zrhopad
1550	zr2ad = zr2ad + zrhopad * zs
1551	zr3ad = zr3ad + zrhopad * zsr * zs
1552	zsrad = zsrad + zrhopad * zr3 * zs
1553	zsad = zsad + zrhopad * ( 2. * zr4 * zs + zr2 &
1554	& + zr3 * zsr )
1555	zrhopad = 0.0_wp
1556
1557	ztad = ztad + zr3ad * (-2.1.6546e-6_wp zt &
1558	& + 1.0227e-4_wp )
1559	zr3ad = 0.0_wp
1560
1561	ztad = ztad + zr2ad * ( ( ( 4.5.3875e-9_wp zt &
1562	& -3.8.2467e-7_wp ) zt &
1563	& +2.7.6438e-5_wp ) zt &
1564	& - 4.0899e-3_wp )
1565	zr2ad = 0.0_wp
1566
1567	ztad = ztad + zr1ad * ( ( ( ( 5.6.536332e-9_wp zt &
1568	& -4.1.120083e-6_wp ) zt &
1569	& +3.1.001685e-4_wp ) zt &
1570	& -2.9.095290e-3_wp ) zt &
1571	& + 6.793952e-2_wp )
1572	zr1ad = 0.0_wp
1573
1574	zsad = zsad + zsrad * ( 1.0 / MAX( 2.*zsr, zeps ) ) &
1575	& * tmask(ji,jj, 1)
1576	zsrad = 0.0_wp
1577
1578	pts_ad(ji,jj,jp_sal) = pts_ad(ji,jj,jp_sal) + zsad
1579	pts_ad(ji,jj,jp_tem) = pts_ad(ji,jj,jp_tem) + ztad
1580	ztad = 0.0_wp
1581	zsad = 0.0_wp
1582	END DO
1583	END DO
1584	!
1585	CASE ( 1 ) !== Linear formulation = F( temperature ) ==!
1586	DO jj = jpjm1, 1, -1
1587	DO ji = fs_jpim1, 1, -1 ! vector opt.
1588	pts_ad(ji,jj,jp_tem) = pts_ad(ji,jj,jp_tem) - prd_ad(ji,jj) * rn_alpha * tmask(ji,jj,1)
1589	prd_ad(ji,jj) = 0.0_wp
1590	END DO
1591	END DO
1592	!
1593	CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==!
1594	DO jj = jpjm1, 1, -1
1595	DO ji = fs_jpim1, 1, -1 ! vector opt.
1596	pts_ad(ji,jj,jp_tem) = pts_ad(ji,jj,jp_tem) - prd_ad(ji,jj) * rn_alpha * tmask(ji,jj,1)
1597	pts_ad(ji,jj,jp_sal) = pts_ad(ji,jj,jp_sal) + prd_ad(ji,jj) * rn_beta * tmask(ji,jj,1)
1598	prd_ad (ji,jj) = 0.0_wp
1599	END DO
1600	END DO
1601	!
1602	END SELECT
1603	!
1604	CALL wrk_dealloc( jpi, jpj, zws )
1605	!
1606	IF( nn_timing == 1 ) CALL timing_stop('eos2d_adj')
1607	!
1608	END SUBROUTINE eos_insitu_2d_adj
1609
1610	SUBROUTINE eos_bn2_tan ( pts, pts_tl, pn2_tl )
1611	!!----------------------------------------------------------------------
1612	!! * ROUTINE eos_bn2_tan *
1613	!!
1614	!! ** Purpose of the direct routine: Compute the local
1615	!! Brunt-Vaisala frequency at the time-step of the input arguments
1616	!!
1617	!! ** Method of the direct routine:
1618	!! * nn_eos = 0 : UNESCO sea water properties
1619	!! The brunt-vaisala frequency is computed using the polynomial
1620	!! polynomial expression of McDougall (1987):
1621	!! N^2 = grav * beta * ( alpha/beta*dk[ t ] - dk[ s ] )/e3w
1622	!! If lk_zdfddm=T, the heat/salt buoyancy flux ratio Rrau is
1623	!! computed and used in zdfddm module :
1624	!! Rrau = alpha/beta * ( dk[ t ] / dk[ s ] )
1625	!! * nn_eos = 1 : linear equation of state (temperature only)
1626	!! N^2 = grav * rn_alpha * dk[ t ]/e3w
1627	!! * nn_eos = 2 : linear equation of state (temperature & salinity)
1628	!! N^2 = grav * (rn_alpha * dk[ t ] - rn_beta * dk[ s ] ) / e3w
1629	!! The use of potential density to compute N^2 introduces e r r o r
1630	!! in the sign of N^2 at great depths. We recommand the use of
1631	!! nn_eos = 0, except for academical studies.
1632	!! Macro-tasked on horizontal slab (jk-loop)
1633	!! N.B. N^2 is set to zero at the first level (JK=1) in inidtr
1634	!! and is never used at this level.
1635	!!
1636	!! ** Action : - pn2 : the brunt-vaisala frequency
1637	!!
1638	!! References :
1639	!! McDougall, T. J., J. Phys. Oceanogr., 17, 1950-1964, 1987.
1640	!!
1641	!! History:
1642	!! ! 08-07 (A. Vidard) First version
1643	!!----------------------------------------------------------------------
1644	!! * Arguments
1645	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts, & ! 1 : potential temperature [Celcius]
1646	! ! 2 : salinity [psu]
1647	& pts_tl ! 1 : TL of potential temperature [Celsius]
1648	! 2 : TL of salinity [psu]
1649	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out ) :: &
1650	& pn2_tl ! TL of potential density (surface referenced)
1651	!! * Local declarations
1652	INTEGER :: ji, jj, jk ! dummy loop indices
1653	REAL(wp) :: &
1654	zgde3w, zt, zs, zh, & ! temporary scalars
1655	zalbet, zbeta ! " "
1656	REAL(wp) :: &
1657	zttl, zstl, & ! temporary scalars
1658	zalbettl, zbetatl ! " "
1659	#if defined key_zdfddm
1660	REAL(wp) :: zds, zdstl ! temporary scalars
1661	#endif
1662
1663	! pn2_tl : interior points only (2=< jk =< jpkm1 )
1664	! --------------------------
1665	SELECT CASE ( nn_eos )
1666
1667	CASE ( 0 ) !== Jackett and McDougall (1994) formulation ==!
1668	DO jk = 2, jpkm1
1669	DO jj = 1, jpj
1670	DO ji = 1, jpi
1671	zgde3w = grav / fse3w(ji,jj,jk)
1672	zt = 0.5 * ( pts(ji,jj,jk,jp_tem) + pts(ji,jj,jk-1,jp_tem) ) ! potential temperature at w-point
1673	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
1674	zh = fsdepw(ji,jj,jk) ! depth in meters at w-point
1675
1676	zalbet = ( ( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt & ! ratio alpha/beta
1677	& - 0.203814e-03_wp ) * zt &
1678	& + 0.170907e-01_wp ) * zt &
1679	& + 0.665157e-01_wp &
1680	& + ( - 0.678662e-05_wp * zs &
1681	& - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs &
1682	& + ( ( - 0.302285e-13_wp * zh &
1683	& - 0.251520e-11_wp * zs &
1684	& + 0.512857e-12_wp * zt * zt ) * zh &
1685	& - 0.164759e-06_wp * zs &
1686	& +( 0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt &
1687	& + 0.380374e-04_wp ) * zh
1688
1689	zbeta = ( ( -0.415613e-09_wp * zt + 0.555579e-07_wp ) * zt & ! beta
1690	& - 0.301985e-05_wp ) * zt &
1691	& + 0.785567e-03_wp &
1692	& + ( 0.515032e-08_wp * zs &
1693	& + 0.788212e-08_wp * zt - 0.356603e-06_wp ) * zs &
1694	& +( ( 0.121551e-17_wp * zh &
1695	& - 0.602281e-15_wp * zs &
1696	& - 0.175379e-14_wp * zt + 0.176621e-12_wp ) * zh &
1697	& + 0.408195e-10_wp * zs &
1698	& + ( - 0.213127e-11_wp * zt + 0.192867e-09_wp ) * zt &
1699	& - 0.121555e-07_wp ) * zh
1700
1701	!! tangent part
1702	zttl = 0.5 * ( pts_tl(ji,jj,jk,jp_tem) + pts_tl(ji,jj,jk-1,jp_tem) ) ! potential temperature at w-point
1703	zstl = 0.5 * ( pts_tl(ji,jj,jk,jp_sal) + pts_tl(ji,jj,jk-1,jp_sal) ) ! salinity anomaly at w-point
1704	zalbettl = ( ( ( -4.0.255019e-07_wp zt & ! ratio alpha/beta
1705	& +3.0.298357e-05_wp ) zt &
1706	& -2.0.203814e-03_wp ) zt &
1707	& + 0.170907e-01_wp &
1708	& - 0.846960e-04_wp * zs &
1709	& - ( 0.933746e-06_wp &
1710	& - ( 2.*0.791325e-08_wp &
1711	& +2.0.512857e-12_wp zh ) * zt ) * zh ) * zttl &
1712	& + ( - 2.0.678662e-05_wp zs &
1713	& - 0.846960e-04_wp * zt &
1714	& + 0.378110e-02_wp &
1715	& + ( - 0.164759e-06_wp &
1716	& - 0.251520e-11_wp * zh ) * zh ) * zstl
1717
1718	zbetatl = ( ( -3.0.415613e-09_wp zt &
1719	& +2.0.555579e-07_wp ) zt &
1720	& - 0.301985e-05_wp &
1721	& + 0.788212e-08_wp * zs &
1722	& + ( -2.0.213127e-11_wp zt &
1723	& - 0.175379e-14_wp * zh &
1724	& + 0.192867e-09_wp ) * zh ) * zttl &
1725	& + ( 2.0.515032e-08_wp zs &
1726	& + 0.788212e-08_wp * zt &
1727	& - 0.356603e-06_wp &
1728	& + ( - 0.602281e-15_wp * zh &
1729	& + 0.408195e-10_wp ) * zh ) * zstl
1730
1731	pn2_tl(ji,jj,jk) = zgde3w * tmask(ji,jj,jk) * ( &
1732	& zbeta * ( zalbet &
1733	& * ( pts_tl(ji,jj,jk-1,jp_tem) - pts_tl(ji,jj,jk,jp_tem) ) &
1734	& + zalbettl &
1735	& * ( pts (ji,jj,jk-1,jp_tem) - pts (ji,jj,jk,jp_tem) ) &
1736	& - ( pts_tl(ji,jj,jk-1,jp_sal) - pts_tl(ji,jj,jk,jp_sal) ) ) &
1737	& + zbetatl * ( zalbet &
1738	& * ( pts (ji,jj,jk-1,jp_tem) - pts (ji,jj,jk,jp_tem) ) &
1739	& - ( pts (ji,jj,jk-1,jp_sal) - pts (ji,jj,jk,jp_sal) ) ) )
1740	#if defined key_zdfddm
1741	zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )
1742	zdstl = ( pts_tl(ji,jj,jk-1,jp_sal) - pts_tl(ji,jj,jk,jp_sal) )
1743	IF ( ABS( zds) <= 1.e-20 ) THEN
1744	zds = 1.e-20
1745	rrau_tl(ji,jj,jk) = zalbettl * &
1746	& ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds &
1747	& + zalbet * &
1748	& ( ( pts_tl(ji,jj,jk-1,jp_tem) - pts_tl(ji,jj,jk,jp_tem) ) / zds )
1749	ELSE
1750	rrau_tl(ji,jj,jk) = zalbettl * &
1751	& ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds &
1752	& + zalbet * &
1753	& ( ( pts_tl(ji,jj,jk-1,jp_tem) - pts_tl(ji,jj,jk,jp_tem) ) / zds &
1754	& - ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) * zdstl/ zds**2 )
1755	ENDIF
1756	#endif
1757	END DO
1758	END DO
1759	END DO
1760	!
1761	CASE ( 1 ) !== Linear formulation = F( temperature ) ==!
1762	DO jk = 2, jpkm1
1763	pn2_tl(:,:,jk) = rn_alpha * ( pts_tl(:,:,jk-1,jp_tem) - pts_tl(:,:,jk,jp_tem) ) * grav / fse3w(:,:,jk) * tmask(:,:,jk)
1764	END DO
1765	!
1766	CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==!
1767	DO jk = 2, jpkm1
1768	pn2_tl(:,:,jk) = ( rn_alpha * ( pts_tl(:,:,jk-1,jp_tem) - pts_tl(:,:,jk,jp_tem) ) &
1769	& - rn_beta * ( pts_tl(:,:,jk-1,jp_sal) - pts_tl(:,:,jk,jp_sal) ) ) &
1770	& * grav / fse3w(:,:,jk) * tmask(:,:,jk)
1771	END DO
1772	#if defined key_zdfddm
1773	DO jk = 2, jpkm1
1774	DO jj = 1, jpj
1775	DO ji = 1, jpi
1776	zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )
1777	zdstl = pts_tl(ji,jj,jk-1,jp_sal) - pts_tl(ji,jj,jk,jp_sal)
1778	IF ( ABS( zds) <= 1.e-20 ) THEN
1779	zds = 1.e-20
1780	rrau_tl(ji,jj,jk) = ralpbet * &
1781	& ( ( pts_tl(ji,jj,jk-1,jp_tem) - pts_tl(ji,jj,jk,jp_tem) ) / zds )
1782	ELSE
1783	rrau_tl(ji,jj,jk) = ralpbet * &
1784	& ( ( pts_tl(ji,jj,jk-1,jp_tem) - pts_tl(ji,jj,jk,jp_tem) ) / zds &
1785	& - ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) * zdstl / zds**2 )
1786	ENDIF
1787	rrau(ji,jj,jk) = ralpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds
1788	END DO
1789	END DO
1790	END DO
1791	#endif
1792	END SELECT
1793	END SUBROUTINE eos_bn2_tan
1794
1795	SUBROUTINE eos_bn2_adj ( pts, pts_ad, pn2_ad )
1796	!!----------------------------------------------------------------------
1797	!! * ROUTINE eos_bn2_adj *
1798	!!
1799	!! ** Purpose of the direct routine: Compute the local
1800	!! Brunt-Vaisala frequency at the time-step of the input arguments
1801	!!
1802	!! ** Method of the direct routine:
1803	!! * nn_eos = 0 : UNESCO sea water properties
1804	!! The brunt-vaisala frequency is computed using the polynomial
1805	!! polynomial expression of McDougall (1987):
1806	!! N^2 = grav * beta * ( alpha/beta*dk[ t ] - dk[ s ] )/e3w
1807	!! If lk_zdfddm=T, the heat/salt buoyancy flux ratio Rrau is
1808	!! computed and used in zdfddm module :
1809	!! Rrau = alpha/beta * ( dk[ t ] / dk[ s ] )
1810	!! * nn_eos = 1 : linear equation of state (temperature only)
1811	!! N^2 = grav * rn_alpha * dk[ t ]/e3w
1812	!! * nn_eos = 2 : linear equation of state (temperature & salinity)
1813	!! N^2 = grav * (rn_alpha * dk[ t ] - rn_beta * dk[ s ] ) / e3w
1814	!! The use of potential density to compute N^2 introduces e r r o r
1815	!! in the sign of N^2 at great depths. We recommand the use of
1816	!! nn_eos = 0, except for academical studies.
1817	!! Macro-tasked on horizontal slab (jk-loop)
1818	!! N.B. N^2 is set to zero at the first level (JK=1) in inidtr
1819	!! and is never used at this level.
1820	!!
1821	!! ** Action : - pn2 : the brunt-vaisala frequency
1822	!!
1823	!! References :
1824	!! McDougall, T. J., J. Phys. Oceanogr., 17, 1950-1964, 1987.
1825	!!
1826	!! History:
1827	!! ! 08-07 (A. Vidard) First version
1828	!!----------------------------------------------------------------------
1829	!! * Arguments
1830	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! 1 : potential temperature [Celcius]
1831	! ! 2 : salinity [psu]
1832	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(inout ) :: pts_ad ! 1 : Adjoint of potential temperature [Celsius]
1833	! 2 : Adjoint of salinity [psu]
1834	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: &
1835	& pn2_ad ! Adjoint of potential density (surface referenced)
1836	!
1837	!! * Local declarations
1838	INTEGER :: ji, jj, jk ! dummy loop indices
1839	REAL(wp) :: &
1840	zgde3w, zt, zs, zh, & ! temporary scalars
1841	zalbet, zbeta ! " "
1842	REAL(wp) :: &
1843	ztad, zsad, & ! temporary scalars
1844	zalbetad, zbetaad ! " "
1845	#if defined key_zdfddm
1846	REAL(wp) :: zds, zdsad ! temporary scalars
1847	#endif
1848
1849	! pn2_tl : interior points only (2=< jk =< jpkm1 )
1850	! --------------------------
1851	zalbetad = 0.0_wp
1852	zbetaad = 0.0_wp
1853	ztad = 0.0_wp
1854	zsad = 0.0_wp
1855	#if defined key_zdfddm
1856	zdsad = 0.0_wp
1857	#endif
1858
1859	SELECT CASE ( nn_eos )
1860
1861	CASE ( 0 ) !== Jackett and McDougall (1994) formulation ==!
1862	DO jk = jpkm1, 2, -1
1863	DO jj = jpj, 1, -1
1864	DO ji = jpi, 1, -1
1865	zgde3w = grav / fse3w(ji,jj,jk)
1866	zt = 0.5 * ( pts(ji,jj,jk,jp_tem) + pts(ji,jj,jk-1,jp_tem) ) ! potential temperature at w-point
1867	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
1868	zh = fsdepw(ji,jj,jk) ! depth in meters at w-point
1869
1870	zalbet = ( ( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt & ! ratio alpha/beta
1871	& - 0.203814e-03_wp ) * zt &
1872	& + 0.170907e-01_wp ) * zt &
1873	& + 0.665157e-01_wp &
1874	& + ( - 0.678662e-05_wp * zs &
1875	& - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs &
1876	& + ( ( - 0.302285e-13_wp * zh &
1877	& - 0.251520e-11_wp * zs &
1878	& + 0.512857e-12_wp * zt * zt ) * zh &
1879	& - 0.164759e-06_wp * zs &
1880	& +( 0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt &
1881	& + 0.380374e-04_wp ) * zh
1882
1883	zbeta = ( ( -0.415613e-09_wp * zt + 0.555579e-07_wp ) * zt & ! beta
1884	& - 0.301985e-05_wp ) * zt &
1885	& + 0.785567e-03_wp &
1886	& + ( 0.515032e-08_wp * zs &
1887	& + 0.788212e-08_wp * zt - 0.356603e-06_wp ) * zs &
1888	& +( ( 0.121551e-17_wp * zh &
1889	& - 0.602281e-15_wp * zs &
1890	& - 0.175379e-14_wp * zt + 0.176621e-12_wp ) * zh &
1891	& + 0.408195e-10_wp * zs &
1892	& + ( - 0.213127e-11_wp * zt + 0.192867e-09_wp ) * zt &
1893	& - 0.121555e-07_wp ) * zh
1894
1895	#if defined key_zdfddm
1896
1897	zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )
1898	IF ( ABS( zds) <= 1.e-20 ) THEN
1899	zds = 1.e-20
1900	zdsad = 0.0_wp
1901	ELSE
1902	zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )
1903	zdsad = rrau_ad(ji,jj,jk) * zalbet ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds*2
1904	zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )
1905	ENDIF
1906	pts_ad(ji,jj,jk-1,jp_tem) = pts_ad(ji,jj,jk-1,jp_tem) + rrau_ad(ji,jj,jk) * zalbet / zds
1907	pts_ad(ji,jj,jk,jp_tem ) = pts_ad(ji,jj,jk,jp_tem ) - rrau_ad(ji,jj,jk) * zalbet / zds
1908	zalbetad = zalbetad + rrau_ad(ji,jj,jk) * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds
1909	rrau_ad(ji,jj,jk) = 0._wp
1910	pts_ad(ji,jj,jk-1,jp_sal) = pts_ad(ji,jj,jk-1,jp_sal) + zdsad
1911	pts_ad(ji,jj,jk,jp_sal ) = pts_ad(ji,jj,jk,jp_sal ) - zdsad
1912	zdsad = 0._wp
1913	#endif
1914	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)
1915	pts_ad(ji,jj,jk,jp_tem ) = pts_ad(ji,jj,jk,jp_tem ) - zalbetzbetazgde3wtmask(ji,jj,jk)pn2_ad(ji,jj,jk)
1916	zalbetad = zalbetad + zbetazgde3wtmask(ji,jj,jk) &
1917	& ( pts (ji,jj,jk-1,jp_tem) - pts (ji,jj,jk,jp_tem) ) pn2_ad(ji,jj,jk)
1918	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)
1919	pts_ad(ji,jj,jk,jp_sal ) = pts_ad(ji,jj,jk,jp_sal ) + zbetatmask(ji,jj,jk)zgde3w*pn2_ad(ji,jj,jk)
1920	zbetaad = zbetaad &
1921	& + zgde3w tmask(ji,jj,jk) ( zalbet * ( pts (ji,jj,jk-1,jp_tem) - pts (ji,jj,jk,jp_tem) ) &
1922	& - ( pts (ji,jj,jk-1,jp_sal) - pts (ji,jj,jk,jp_sal) ) )*pn2_ad(ji,jj,jk)
1923
1924	pn2_ad(ji,jj,jk) = 0.0_wp
1925
1926	ztad = ztad + ( ( -3.0.415613e-09_wp zt &
1927	& +2.0.555579e-07_wp ) zt &
1928	& - 0.301985e-05_wp &
1929	& + 0.788212e-08_wp * zs &
1930	& + ( -2.0.213127e-11_wp zt &
1931	& - 0.175379e-14_wp * zh &
1932	& + 0.192867e-09_wp ) * zh ) *zbetaad
1933
1934	zsad = zsad + ( 2.0.515032e-08_wp zs &
1935	& + 0.788212e-08_wp * zt &
1936	& - 0.356603e-06_wp &
1937	& + ( - 0.602281e-15_wp * zh &
1938	& + 0.408195e-10_wp ) * zh ) * zbetaad
1939
1940	zbetaad = 0.0_wp
1941
1942	ztad = ztad + ( ( ( -4.0.255019e-07_wp zt &! ratio alpha/beta
1943	& +3.0.298357e-05_wp ) zt &
1944	& -2.0.203814e-03_wp ) zt &
1945	& + 0.170907e-01_wp &
1946	& - 0.846960e-04_wp * zs &
1947	& - ( 0.933746e-06_wp &
1948	& - ( 2.*0.791325e-08_wp &
1949	& +2.0.512857e-12_wp zh ) * zt ) * zh &
1950	& ) *zalbetad
1951
1952	zsad = zsad + ( - 2.0.678662e-05_wp zs &
1953	& - 0.846960e-04_wp * zt &
1954	& + 0.378110e-02_wp &
1955	& + ( - 0.164759e-06_wp &
1956	& - 0.251520e-11_wp * zh ) * zh &
1957	& ) *zalbetad
1958
1959	zalbetad = 0.0_wp
1960
1961
1962	pts_ad(ji,jj,jk,jp_sal) = pts_ad(ji,jj,jk,jp_sal) + 0.5 * zsad
1963	pts_ad(ji,jj,jk-1,jp_sal) = pts_ad(ji,jj,jk-1,jp_sal) + 0.5 * zsad
1964	zsad = 0.0_wp
1965
1966	pts_ad(ji,jj,jk,jp_tem) = pts_ad(ji,jj,jk,jp_tem) + 0.5 * ztad
1967	pts_ad(ji,jj,jk-1,jp_tem) = pts_ad(ji,jj,jk-1,jp_tem) + 0.5 * ztad
1968	ztad = 0.0_wp
1969
1970	END DO
1971	END DO
1972	END DO
1973	!
1974	CASE ( 1 ) !== Linear formulation = F( temperature ) ==!
1975	DO jk = jpkm1, 2, -1
1976	pts_ad(:,:,jk-1,jp_tem) = pts_ad(:,:,jk-1,jp_tem) + rn_alpha * pn2_ad(:,:,jk) &
1977	& * grav / fse3w(:,:,jk) * tmask(:,:,jk)
1978	pts_ad(:,:,jk,jp_tem ) = pts_ad(:,:,jk,jp_tem ) - rn_alpha * pn2_ad(:,:,jk) &
1979	& * grav / fse3w(:,:,jk) * tmask(:,:,jk)
1980	pn2_ad(:,:,jk) = 0.0_wp
1981	END DO
1982	!
1983	CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==!
1984	#if defined key_zdfddm
1985	DO jk = jpkm1, 2, -1
1986	DO jj = jpj, 1, -1
1987	DO ji = jpi, 1, -1
1988	zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )
1989	IF ( ABS( zds) <= 1.e-20 ) THEN
1990	zds = 1.e-20
1991	zdsad = 0.0_wp
1992	ELSE
1993	zdsad = zdsad - rrau_ad(ji,jj,jk) * ralpbet &
1994	& * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds**2
1995	ENDIF
1996	rrau(ji,jj,jk) = ralpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds
1997	pts_ad(ji,jj,jk-1,jp_tem) = pts_ad(ji,jj,jk-1,jp_tem) &
1998	& + rrau_ad(ji,jj,jk) * ralpbet / zds
1999	pts_ad(ji,jj,jk,jp_tem ) = pts_ad(ji,jj,jk,jp_tem ) &
2000	& - rrau_ad(ji,jj,jk) * ralpbet / zds
2001	rrau_ad(ji,jj,jk) = 0._wp
2002
2003	pts_ad(ji,jj,jk-1,jp_sal) = pts_ad(ji,jj,jk-1,jp_sal) + zdsad
2004	pts_ad(ji,jj,jk,jp_sal ) = pts_ad(ji,jj,jk,jp_sal ) - zdsad
2005	zdsad = 0._wp
2006	END DO
2007	END DO
2008	END DO
2009	#endif
2010	DO jk = jpkm1, 2, -1
2011	pts_ad(:,:,jk-1,jp_tem) = pts_ad(:,:,jk-1,jp_tem) + rn_alpha * pn2_ad(:,:,jk) &
2012	& * grav / fse3w(:,:,jk) * tmask(:,:,jk)
2013	pts_ad(:,:,jk,jp_tem ) = pts_ad(:,:,jk,jp_tem ) - rn_alpha * pn2_ad(:,:,jk) &
2014	& * grav / fse3w(:,:,jk) * tmask(:,:,jk)
2015	pts_ad(:,:,jk-1,jp_sal) = pts_ad(:,:,jk-1,jp_sal) - rn_beta * pn2_ad(:,:,jk) &
2016	& * grav / fse3w(:,:,jk) * tmask(:,:,jk)
2017	pts_ad(:,:,jk,jp_sal ) = pts_ad(:,:,jk,jp_sal ) + rn_beta * pn2_ad(:,:,jk) &
2018	& * grav / fse3w(:,:,jk) * tmask(:,:,jk)
2019	pn2_ad(:,:,jk) = 0.0_wp
2020	END DO
2021	END SELECT
2022	END SUBROUTINE eos_bn2_adj
2023
2024	SUBROUTINE eos_alpbet_tan( pts, pts_tl, palpbet_tl, beta0_tl )
2025	!!----------------------------------------------------------------------
2026	!! * ROUTINE eos_alpbet_tan *
2027	!!
2028	!! ** Purpose of the direct routine :
2029	!! Calculates the in situ thermal/haline expansion ratio at T-points
2030	!!
2031	!! ** Method of the direct routine :
2032	!! calculates alpha / beta ratio at T-points
2033	!! * nn_eos = 0 : UNESCO sea water properties
2034	!! The alpha/beta ratio is returned as 3-D array palpbet using the polynomial
2035	!! polynomial expression of McDougall (1987).
2036	!! Scalar beta0 is returned = 1.
2037	!! * nn_eos = 1 : linear equation of state (temperature only)
2038	!! The ratio is undefined, so we return alpha as palpbet
2039	!! Scalar beta0 is returned = 0.
2040	!! * nn_eos = 2 : linear equation of state (temperature & salinity)
2041	!! The alpha/beta ratio is returned as ralpbet
2042	!! Scalar beta0 is returned = 1.
2043	!!
2044	!! ** Action : - palpbet : thermal/haline expansion ratio at T-points
2045	!! : beta0 : 1. or 0.
2046	!!----------------------------------------------------------------------
2047	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts_tl, & ! linear tangent of pot. temperature & salinity
2048	& pts ! pot. temperature & salinity
2049	REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: palpbet_tl ! thermal/haline expansion ratio
2050	REAL(wp), INTENT( out) :: beta0_tl ! set = 1 except with case 1 eos, rho=rho(T)
2051	!!
2052	INTEGER :: ji, jj, jk ! dummy loop indices
2053	REAL(wp) :: zt, zs, zh, & ! local scalars
2054	& zttl, zstl
2055	!!----------------------------------------------------------------------
2056	!
2057	IF( nn_timing == 1 ) CALL timing_start('eos_alpbet_tan')
2058	!
2059	SELECT CASE ( nn_eos )
2060	!
2061	CASE ( 0 ) ! Jackett and McDougall (1994) formulation
2062	DO jk = 1, jpk
2063	DO jj = 1, jpj
2064	DO ji = 1, jpi
2065	zt = pts(ji,jj,jk,jp_tem) ! potential temperature
2066	zs = pts(ji,jj,jk,jp_sal) - 35._wp ! salinity anomaly (s-35)
2067	zh = fsdept(ji,jj,jk) ! depth in meters
2068	!! Tangent part
2069	zttl = pts_tl(ji,jj,jk,jp_tem) ! potential temperature
2070	zstl = pts_tl(ji,jj,jk,jp_sal) ! salinity anomaly (s-35)
2071	palpbet_tl(ji,jj,jk) = &
2072	& ( ( ( ( - 4. * 0.255019e-07_wp * zt &
2073	& + 3. * 0.298357e-05_wp ) * zt &
2074	& - 2. * 0.203814e-03_wp ) * zt &
2075	& + 0.170907e-01_wp * zt ) &
2076	& - 0.846960e-04_wp * zs &
2077	& + ( ( 2. * 0.512857e-12_wp * zt ) * zh &
2078	& + ( 2. * 0.791325e-08_wp * zt &
2079	& - 0.933746e-06_wp ) ) * zh ) * zttl &
2080	& + ( ( - 2. * 0.678662e-05_wp * zs &
2081	& - 0.846960e-04_wp * zt &
2082	& + 0.378110e-02_wp ) &
2083	& + ( - 0.251520e-11_wp * zh &
2084	& - 0.164759e-06_wp ) * zh ) * zstl
2085	END DO
2086	END DO
2087	END DO
2088	beta0_tl = 0._wp
2089	!
2090	CASE ( 1 ) !== Linear formulation = F( temperature ) ==!
2091	palpbet_tl(:,:,:) = 0._wp
2092	beta0_tl = 0._wp
2093	!
2094	CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==!
2095	palpbet_tl(:,:,:) = 0._wp
2096	beta0_tl = 0._wp
2097	!
2098	CASE DEFAULT
2099	IF(lwp) WRITE(numout,cform_err)
2100	IF(lwp) WRITE(numout,*) ' bad flag value for nn_eos = ', nn_eos
2101	nstop = nstop + 1
2102	!
2103	END SELECT
2104	!
2105	IF( nn_timing == 1 ) CALL timing_stop('eos_alpbet_tan')
2106	!
2107	END SUBROUTINE eos_alpbet_tan
2108
2109	SUBROUTINE eos_alpbet_adj( pts, pts_ad, palpbet_ad, beta0_ad )
2110	!!----------------------------------------------------------------------
2111	!! * ROUTINE eos_alpbet_adj *
2112	!!
2113	!! ** Purpose of the direct routine :
2114	!! Calculates the in situ thermal/haline expansion ratio at T-points
2115	!!
2116	!! ** Method of the direct routine :
2117	!! calculates alpha / beta ratio at T-points
2118	!! * nn_eos = 0 : UNESCO sea water properties
2119	!! The alpha/beta ratio is returned as 3-D array palpbet using the polynomial
2120	!! polynomial expression of McDougall (1987).
2121	!! Scalar beta0 is returned = 1.
2122	!! * nn_eos = 1 : linear equation of state (temperature only)
2123	!! The ratio is undefined, so we return alpha as palpbet
2124	!! Scalar beta0 is returned = 0.
2125	!! * nn_eos = 2 : linear equation of state (temperature & salinity)
2126	!! The alpha/beta ratio is returned as ralpbet
2127	!! Scalar beta0 is returned = 1.
2128	!!
2129	!! ** Action : - palpbet : thermal/haline expansion ratio at T-points
2130	!! : beta0 : 1. or 0.
2131	!!----------------------------------------------------------------------
2132	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(inout) :: pts_ad ! linear tangent of pot. temperature & salinity
2133	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in) :: pts ! pot. temperature & salinity
2134	REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(inout) :: palpbet_ad ! thermal/haline expansion ratio
2135	REAL(wp), INTENT(inout) :: beta0_ad ! set = 1 except with case 1 eos, rho=rho(T)
2136	!!
2137	INTEGER :: ji, jj, jk ! dummy loop indices
2138	REAL(wp) :: zt, zs, zh, & ! local scalars
2139	& ztad, zsad
2140	!!----------------------------------------------------------------------
2141	!
2142	ztad = 0.0_wp
2143	zsad = 0.0_wp
2144	IF( nn_timing == 1 ) CALL timing_start('eos_alpbet_adj')
2145	!
2146	SELECT CASE ( nn_eos )
2147	!
2148	CASE ( 0 ) ! Jackett and McDougall (1994) formulation
2149	DO jk = jpk, 1, -1
2150	DO jj = jpj, 1, -1
2151	DO ji = jpi, 1, -1
2152	zt = pts(ji,jj,jk,jp_tem) ! potential temperature
2153	zs = pts(ji,jj,jk,jp_sal) - 35._wp ! salinity anomaly (s-35)
2154	zh = fsdept(ji,jj,jk) ! depth in meters
2155	!! Adjoint part
2156	ztad = ztad + ( ( ( ( - 4. * 0.255019e-07_wp * zt &
2157	& + 3. * 0.298357e-05_wp ) * zt &
2158	& - 2. * 0.203814e-03_wp ) * zt &
2159	& + 0.170907e-01_wp * zt ) &
2160	& - 0.846960e-04_wp * zs &
2161	& + ( ( 2. * 0.512857e-12_wp * zt ) * zh &
2162	& + ( 2. * 0.791325e-08_wp * zt &
2163	& - 0.933746e-06_wp ) ) * zh ) * palpbet_ad(ji,jj,jk)
2164	zsad = zsad + ( ( - 2. * 0.678662e-05_wp * zs &
2165	& - 0.846960e-04_wp * zt &
2166	& + 0.378110e-02_wp ) &
2167	& + ( - 0.251520e-11_wp * zh &
2168	& - 0.164759e-06_wp ) * zh ) * palpbet_ad(ji,jj,jk)
2169	palpbet_ad(ji,jj,jk) = 0.0_wp
2170	pts_ad(ji,jj,jk,jp_tem) = pts_ad(ji,jj,jk,jp_tem) + ztad
2171	pts_ad(ji,jj,jk,jp_sal) = pts_ad(ji,jj,jk,jp_sal) + zsad
2172	ztad = 0.0_wp
2173	zsad = 0.0_wp
2174	END DO
2175	END DO
2176	END DO
2177	beta0_ad = 0._wp
2178	!
2179	CASE ( 1 ) !== Linear formulation = F( temperature ) ==!
2180	palpbet_ad(:,:,:) = 0._wp
2181	beta0_ad = 0._wp
2182	!
2183	CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==!
2184	palpbet_ad(:,:,:) = 0._wp
2185	beta0_ad = 0._wp
2186	!
2187	CASE DEFAULT
2188	IF(lwp) WRITE(numout,cform_err)
2189	IF(lwp) WRITE(numout,*) ' bad flag value for nn_eos = ', nn_eos
2190	nstop = nstop + 1
2191	!
2192	END SELECT
2193	!
2194	IF( nn_timing == 1 ) CALL timing_stop('eos_alpbet_adj')
2195	!
2196	END SUBROUTINE eos_alpbet_adj
2197
2198	#if defined key_tam
2199	SUBROUTINE eos_insitu_adj_tst( kumadt )
2200	!!-----------------------------------------------------------------------
2201	!!
2202	!! * ROUTINE eos_adj_tst *
2203	!!
2204	!! ** Purpose : Test the adjoint routine.
2205	!!
2206	!! ** Method : Verify the scalar product
2207	!!
2208	!! ( L dx )^T W dy = dx^T L^T W dy
2209	!!
2210	!! where L = tangent routine
2211	!! L^T = adjoint routine
2212	!! W = diagonal matrix of scale factors
2213	!! dx = input perturbation (random field)
2214	!! dy = L dx
2215	!!
2216	!!
2217	!! History :
2218	!! ! 08-07 (A. Vidard)
2219	!!-----------------------------------------------------------------------
2220	!! * Modules used
2221
2222	!! * Arguments
2223	INTEGER, INTENT(IN) :: &
2224	& kumadt ! Output unit
2225	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
2226	zts ! potential temperature
2227	! salinity
2228	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
2229	& zts_adout ! potential temperature
2230	! salinity
2231	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2232	& zrd_adin ! anomaly density
2233	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
2234	& zts_tlin ! potential temperature
2235	! salinity
2236	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
2237	& znts ! potential temperature
2238	! salinity
2239	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2240	& zrd_tlout ! anomaly density
2241	REAL(KIND=wp) :: &
2242	& zsp1, & ! scalar product involving the tangent routine
2243	& zsp2 ! scalar product involving the adjoint routine
2244	INTEGER :: &
2245	& ji, &
2246	& jj, &
2247	& jk, &
2248	& jn, &
2249	& jeos
2250	CHARACTER(LEN=14) :: cl_name
2251	ALLOCATE( &
2252	& zts( jpi, jpj, jpk, jpts ), &
2253	& znts( jpi, jpj, jpk, jpts ), &
2254	& zts_adout( jpi, jpj, jpk,jpts ), &
2255	& zrd_adin( jpi, jpj, jpk ), &
2256	& zts_tlin( jpi, jpj, jpk,jpts ), &
2257	& zrd_tlout(jpi, jpj, jpk ) )
2258	! Initialize the reference state
2259	zts = tsn
2260	! store initial nn_eos
2261	jeos = nn_eos
2262	DO jn = 0, 2
2263	nn_eos = jn
2264	!=============================================================
2265	! 1) dx = ( T ) and dy = ( T )
2266	!=============================================================
2267
2268	!--------------------------------------------------------------------
2269	! Reset the tangent and adjoint variables
2270	!--------------------------------------------------------------------
2271	zts_tlin(:,:,:,:) = 0.0_wp
2272	zrd_tlout(:,:,:) = 0.0_wp
2273	zts_adout(:,:,:,:) = 0.0_wp
2274	zrd_adin(:,:,:) = 0.0_wp
2275
2276	!--------------------------------------------------------------------
2277	! Initialize the tangent input with random noise: dx
2278	!--------------------------------------------------------------------
2279	CALL grid_random( znts(:,:,:,jp_tem), 'T', 0.0_wp, stdt )
2280	CALL grid_random( znts(:,:,:,jp_sal), 'T', 0.0_wp, stds )
2281
2282	DO jk = 1, jpk
2283	DO jj = nldj, nlej
2284	DO ji = nldi, nlei
2285	zts_tlin(ji,jj,jk,jp_tem) = znts(ji,jj,jk,jp_tem)
2286	zts_tlin(ji,jj,jk,jp_sal) = znts(ji,jj,jk,jp_sal)
2287	END DO
2288	END DO
2289	END DO
2290	CALL eos_insitu_tan(zts, zts_tlin, zrd_tlout)
2291
2292	DO jk = 1, jpk
2293	DO jj = nldj, nlej
2294	DO ji = nldi, nlei
2295	zrd_adin(ji,jj,jk) = zrd_tlout(ji,jj,jk) &
2296	& * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)&
2297	& * tmask(ji,jj,jk)
2298	END DO
2299	END DO
2300	END DO
2301
2302	!--------------------------------------------------------------------
2303	! Compute the scalar product: ( L dx )^T W dy
2304	!--------------------------------------------------------------------
2305
2306	zsp1 = DOT_PRODUCT( zrd_tlout, zrd_adin )
2307
2308	!--------------------------------------------------------------------
2309	! Call the adjoint routine: dx^* = L^T dy^*
2310	!--------------------------------------------------------------------
2311	CALL eos_insitu_adj(zts, zts_adout, zrd_adin)
2312	zsp2 = DOT_PRODUCT( zts_tlin(:,:,:,jp_tem), zts_adout(:,:,:,jp_tem) ) + &
2313	& DOT_PRODUCT( zts_tlin(:,:,:,jp_sal), zts_adout(:,:,:,jp_sal) )
2314
2315	! Compare the scalar products
2316
2317	! Compare the scalar products
2318	! 14 char:'12345678901234'
2319	SELECT CASE( jn )
2320	CASE (0) ; cl_name = 'eos_adj ins T1'
2321	CASE (1) ; cl_name = 'eos_adj ins T2'
2322	CASE (2) ; cl_name = 'eos_adj ins T3'
2323	END SELECT
2324	CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
2325	ENDDO
2326	! restore initial nn_eos
2327	nn_eos = jeos
2328
2329	! Deallocate memory
2330
2331	DEALLOCATE( &
2332	& zts, &
2333	& zts_adout, &
2334	& zrd_adin, &
2335	& zts_tlin, &
2336	& zrd_tlout, &
2337	& znts &
2338	& )
2339	END SUBROUTINE eos_insitu_adj_tst
2340
2341	SUBROUTINE eos_insitu_pot_adj_tst( kumadt )
2342	!!-----------------------------------------------------------------------
2343	!!
2344	!! * ROUTINE eos_adj_tst *
2345	!!
2346	!! ** Purpose : Test the adjoint routine.
2347	!!
2348	!! ** Method : Verify the scalar product
2349	!!
2350	!! ( L dx )^T W dy = dx^T L^T W dy
2351	!!
2352	!! where L = tangent routine
2353	!! L^T = adjoint routine
2354	!! W = diagonal matrix of scale factors
2355	!! dx = input perturbation (random field)
2356	!! dy = L dx
2357	!!
2358	!! ** Action : Separate tests are applied for the following dx and dy:
2359	!!
2360	!! 1) dx = ( SSH ) and dy = ( SSH )
2361	!!
2362	!! History :
2363	!! ! 08-07 (A. Vidard)
2364	!!-----------------------------------------------------------------------
2365	!! * Modules used
2366
2367	!! * Arguments
2368	INTEGER, INTENT(IN) :: &
2369	& kumadt ! Output unit
2370	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
2371	zts ! potential temperature
2372	! salinity
2373	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
2374	& zts_adout ! potential temperature
2375	! salinity
2376	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2377	& zrd_adin ! anomaly density
2378	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2379	& zrhop_adin ! volume mass
2380	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
2381	& zts_tlin ! potential temperature
2382	! salinity
2383	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
2384	& znts ! potential temperature
2385	! salinity
2386	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2387	& zrd_tlout ! anomaly density
2388	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2389	& zrhop_tlout ! volume mass
2390	REAL(KIND=wp) :: &
2391	& zsp1, & ! scalar product involving the tangent routine
2392	& zsp2 ! scalar product involving the adjoint routine
2393	INTEGER :: &
2394	& ji, &
2395	& jj, &
2396	& jk, &
2397	& jn, &
2398	& jeos
2399	CHARACTER(LEN=14) :: cl_name
2400
2401	! Allocate memory
2402	ALLOCATE( &
2403	& zts( jpi, jpj, jpk, jpts ), &
2404	& zts_adout( jpi, jpj, jpk, jpts ), &
2405	& zrhop_adin( jpi, jpj, jpk ), &
2406	& zrd_adin( jpi, jpj, jpk ), &
2407	& zts_tlin( jpi, jpj, jpk, jpts ), &
2408	& znts( jpi, jpj, jpk, jpts ), &
2409	& zrd_tlout(jpi, jpj, jpk ), &
2410	& zrhop_tlout(jpi, jpj, jpk ) )
2411
2412	! Initialize random field standard deviationsthe reference state
2413	zts = tsn
2414
2415	! store initial nn_eos
2416	jeos = nn_eos
2417	DO jn = 0, 2
2418	nn_eos = jn
2419	!=============================================
2420	! testing of eos_insitu_pot
2421	!=============================================
2422
2423	!=============================================================
2424	! 1) dx = ( T ) and dy = ( T )
2425	!=============================================================
2426
2427	!--------------------------------------------------------------------
2428	! Reset the tangent and adjoint variables
2429	!--------------------------------------------------------------------
2430	zts_tlin(:,:,:,:) = 0.0_wp
2431	zrd_tlout(:,:,:) = 0.0_wp
2432	zrhop_tlout(:,:,:) = 0.0_wp
2433	zts_adout(:,:,:,:) = 0.0_wp
2434	zrhop_adin(:,:,:) = 0.0_wp
2435	zrd_adin(:,:,:) = 0.0_wp
2436
2437	!--------------------------------------------------------------------
2438	! Initialize the tangent input with random noise: dx
2439	!--------------------------------------------------------------------
2440	CALL grid_random( znts(:,:,:,jp_tem), 'T', 0.0_wp, stdt )
2441	CALL grid_random( znts(:,:,:,jp_sal), 'T', 0.0_wp, stds )
2442	DO jk = 1, jpk
2443	DO jj = nldj, nlej
2444	DO ji = nldi, nlei
2445	zts_tlin(ji,jj,jk,:) = znts(ji,jj,jk,:)
2446	END DO
2447	END DO
2448	END DO
2449	!--------------------------------------------------------------------
2450	! Call the tangent routine: dy = L dx
2451	!--------------------------------------------------------------------
2452
2453	call eos_insitu_pot_tan ( zts, zts_tlin, zrd_tlout, zrhop_tlout )
2454
2455	!--------------------------------------------------------------------
2456	! Initialize the adjoint variables: dy^* = W dy
2457	!--------------------------------------------------------------------
2458	DO jk = 1, jpk
2459	DO jj = nldj, nlej
2460	DO ji = nldi, nlei
2461	zrd_adin(ji,jj,jk) = zrd_tlout(ji,jj,jk) &
2462	& * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)&
2463	& * tmask(ji,jj,jk)
2464	zrhop_adin(ji,jj,jk) = zrhop_tlout(ji,jj,jk) &
2465	& * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)&
2466	& * tmask(ji,jj,jk)
2467	END DO
2468	END DO
2469	END DO
2470
2471	!--------------------------------------------------------------------
2472	! Compute the scalar product: ( L dx )^T W dy
2473	!--------------------------------------------------------------------
2474
2475	zsp1 = DOT_PRODUCT( zrd_tlout , zrd_adin ) &
2476	& + DOT_PRODUCT( zrhop_tlout, zrhop_adin )
2477	!--------------------------------------------------------------------
2478	! Call the adjoint routine: dx^* = L^T dy^*
2479	!--------------------------------------------------------------------
2480
2481	CALL eos_insitu_pot_adj( zts, zts_adout, zrd_adin, zrhop_adin )
2482	!--------------------------------------------------------------------
2483	! Compute the scalar product: dx^T L^T W dy
2484	!--------------------------------------------------------------------
2485
2486	zsp2 = DOT_PRODUCT( zts_tlin(:,:,:,jp_tem), zts_adout(:,:,:,jp_tem) ) + &
2487	& DOT_PRODUCT( zts_tlin(:,:,:,jp_sal), zts_adout(:,:,:,jp_sal) )
2488	! Compare the scalar products
2489
2490	! 14 char:'12345678901234'
2491	SELECT CASE( jn )
2492	CASE (0) ; cl_name = 'eos_adj pot T1'
2493	CASE (1) ; cl_name = 'eos_adj pot T2'
2494	CASE (2) ; cl_name = 'eos_adj pot T3'
2495	END SELECT
2496	CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
2497
2498	ENDDO
2499
2500	! restore initial nn_eos
2501	nn_eos = jeos
2502
2503	! Deallocate memory
2504	DEALLOCATE( &
2505	& zts, &
2506	& zts_adout, &
2507	& zrd_adin, &
2508	& zrhop_adin, &
2509	& zts_tlin, &
2510	& zrd_tlout, &
2511	& zrhop_tlout,&
2512	& znts )
2513	END SUBROUTINE eos_insitu_pot_adj_tst
2514
2515	SUBROUTINE eos_alpbet_adj_tst(kumadt)
2516	!!-----------------------------------------------------------------------
2517	!!
2518	!! * ROUTINE eos_adj_tst *
2519	!!
2520	!! ** Purpose : Test the adjoint routine.
2521	!!
2522	!! ** Method : Verify the scalar product
2523	!!
2524	!! ( L dx )^T W dy = dx^T L^T W dy
2525	!!
2526	!! where L = tangent routine
2527	!! L^T = adjoint routine
2528	!! W = diagonal matrix of scale factors
2529	!! dx = input perturbation (random field)
2530	!! dy = L dx
2531	!!
2532	!! ** Action : Separate tests are applied for the following dx and dy:
2533	!!
2534	!! 1) dx = ( SSH ) and dy = ( SSH )
2535	!!
2536	!! History :
2537	!! ! 05-12 (P.-A. Bouttier)
2538	!!-----------------------------------------------------------------------
2539	!! * Modules used
2540
2541	!! * Arguments
2542	INTEGER, INTENT(IN) :: &
2543	& kumadt ! Output unit
2544	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
2545	& zts ! potential temperature
2546	! salinity
2547	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
2548	& zts_adout ! potential temperature
2549	! salinity
2550	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2551	& zpalpbet_adin
2552	REAL(wp) :: &
2553	& zbeta0_adin
2554	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
2555	& zts_tlin ! potential temperature
2556	! salinity
2557	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
2558	& znts ! potential temperature
2559	! salinity
2560	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2561	& zpalpbet_tlout
2562	REAL(wp) :: &
2563	& zbeta0_tlout
2564	REAL(KIND=wp) :: &
2565	& zsp1, & ! scalar product involving the tangent routine
2566	& zsp2 ! scalar product involving the adjoint routine
2567	INTEGER :: &
2568	& ji, &
2569	& jj, &
2570	& jk, &
2571	& jn, &
2572	& jeos
2573	CHARACTER(LEN=14) :: cl_name
2574
2575	! Allocate memory
2576	ALLOCATE( &
2577	& zts( jpi, jpj, jpk, jpts ), &
2578	& zts_adout( jpi, jpj, jpk, jpts ), &
2579	& zpalpbet_adin( jpi, jpj, jpk ), &
2580	& zts_tlin( jpi, jpj, jpk, jpts ), &
2581	& znts( jpi, jpj, jpk, jpts ), &
2582	& zpalpbet_tlout(jpi, jpj, jpk ) )
2583
2584
2585	! Initialize random field standard deviationsthe reference state
2586	zts = tsn
2587
2588	! store initial nn_eos
2589	jeos = nn_eos
2590	DO jn = 0, 2
2591	nn_eos = jn
2592	!=============================================
2593	! testing of eos_insitu_pot
2594	!=============================================
2595
2596	!=============================================================
2597	! 1) dx = ( T ) and dy = ( T )
2598	!=============================================================
2599
2600	!--------------------------------------------------------------------
2601	! Reset the tangent and adjoint variables
2602	!--------------------------------------------------------------------
2603	zts_tlin(:,:,:,:) = 0.0_wp
2604	zpalpbet_tlout(:,:,:) = 0.0_wp
2605	zbeta0_tlout = 0.0_wp
2606	zts_adout(:,:,:,:) = 0.0_wp
2607	zbeta0_adin = 0.0_wp
2608	zpalpbet_adin(:,:,:) = 0.0_wp
2609
2610	!--------------------------------------------------------------------
2611	! Initialize the tangent input with random noise: dx
2612	!--------------------------------------------------------------------
2613	CALL grid_random( znts(:,:,:,jp_tem), 'T', 0.0_wp, stdt )
2614	CALL grid_random( znts(:,:,:,jp_sal), 'T', 0.0_wp, stds )
2615	DO jk = 1, jpk
2616	DO jj = nldj, nlej
2617	DO ji = nldi, nlei
2618	zts_tlin(ji,jj,jk,:) = znts(ji,jj,jk,:)
2619	END DO
2620	END DO
2621	END DO
2622	!--------------------------------------------------------------------
2623	! Call the tangent routine: dy = L dx
2624	!--------------------------------------------------------------------
2625
2626	call eos_alpbet_tan ( zts, zts_tlin, zpalpbet_tlout, zbeta0_tlout )
2627
2628	!--------------------------------------------------------------------
2629	! Initialize the adjoint variables: dy^* = W dy
2630	!--------------------------------------------------------------------
2631	DO jk = 1, jpk
2632	DO jj = nldj, nlej
2633	DO ji = nldi, nlei
2634	zpalpbet_adin(ji,jj,jk) = zpalpbet_tlout(ji,jj,jk) &
2635	& * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)&
2636	& * tmask(ji,jj,jk)
2637	END DO
2638	END DO
2639	END DO
2640	zbeta0_adin = zbeta0_tlout
2641
2642	!--------------------------------------------------------------------
2643	! Compute the scalar product: ( L dx )^T W dy
2644	!--------------------------------------------------------------------
2645
2646	zsp1 = DOT_PRODUCT( zpalpbet_tlout , zpalpbet_adin ) &
2647	& + zbeta0_tlout * zbeta0_adin
2648	!--------------------------------------------------------------------
2649	! Call the adjoint routine: dx^* = L^T dy^*
2650	!--------------------------------------------------------------------
2651
2652	CALL eos_alpbet_adj( zts, zts_adout, zpalpbet_adin, zbeta0_adin )
2653
2654	!--------------------------------------------------------------------
2655	! Compute the scalar product: dx^T L^T W dy
2656	!--------------------------------------------------------------------
2657
2658	zsp2 = DOT_PRODUCT( zts_tlin(:,:,:,jp_tem), zts_adout(:,:,:,jp_tem) ) + &
2659	& DOT_PRODUCT( zts_tlin(:,:,:,jp_sal), zts_adout(:,:,:,jp_sal) )
2660	! Compare the scalar products
2661
2662	! 14 char:'12345678901234'
2663	SELECT CASE( jn )
2664	CASE (0) ; cl_name = 'eos_adj ab T1'
2665	CASE (1) ; cl_name = 'eos_adj ab T2'
2666	CASE (2) ; cl_name = 'eos_adj ab T3'
2667	END SELECT
2668	CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
2669
2670	ENDDO
2671
2672	! restore initial nn_eos
2673	nn_eos = jeos
2674
2675	! Deallocate memory
2676	DEALLOCATE( &
2677	& zts, &
2678	& zts_adout, &
2679	& zpalpbet_adin, &
2680	& zts_tlin, &
2681	& zpalpbet_tlout, &
2682	& znts )
2683
2684	END SUBROUTINE eos_alpbet_adj_tst
2685
2686	SUBROUTINE eos_insitu_2d_adj_tst( kumadt )
2687	!!-----------------------------------------------------------------------
2688	!!
2689	!! * ROUTINE eos_adj_tst *
2690	!!
2691	!! ** Purpose : Test the adjoint routine.
2692	!!
2693	!! ** Method : Verify the scalar product
2694	!!
2695	!! ( L dx )^T W dy = dx^T L^T W dy
2696	!!
2697	!! where L = tangent routine
2698	!! L^T = adjoint routine
2699	!! W = diagonal matrix of scale factors
2700	!! dx = input perturbation (random field)
2701	!! dy = L dx
2702	!!
2703	!! ** Action : Separate tests are applied for the following dx and dy:
2704	!!
2705	!! 1) dx = ( SSH ) and dy = ( SSH )
2706	!!
2707	!! History :
2708	!! ! 08-07 (A. Vidard)
2709	!!-----------------------------------------------------------------------
2710	!! * Modules used
2711
2712	!! * Arguments
2713	INTEGER, INTENT(IN) :: &
2714	& kumadt ! Output unit
2715	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
2716	zdep ! depth
2717	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2718	zts ! potential temperature
2719	! salinity
2720	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2721	& zts_adout ! potential temperature
2722	! salinity
2723	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
2724	& zrd_adin ! anomaly density
2725	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2726	& zts_tlin ! potential temperature
2727	! salinity
2728	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2729	& znts ! potential temperature
2730	! salinity
2731	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
2732	& zrd_tlout ! anomaly density
2733	REAL(KIND=wp) :: &
2734	& zsp1, & ! scalar product involving the tangent routine
2735	& zsp2 ! scalar product involving the adjoint routine
2736	INTEGER :: &
2737	& ji, &
2738	& jj, &
2739	& jn, &
2740	& jeos
2741	CHARACTER(LEN=14) :: cl_name
2742	! Allocate memory
2743
2744	ALLOCATE( &
2745	& zdep( jpi, jpj), &
2746	& zts( jpi, jpj, jpts ), &
2747	& znts( jpi, jpj, jpts ), &
2748	& zts_adout( jpi, jpj, jpts ), &
2749	& zrd_adin( jpi, jpj ), &
2750	& zts_tlin( jpi, jpj, jpts ), &
2751	& zrd_tlout(jpi, jpj ) )
2752
2753	! Initialize the reference state
2754	zts(:,:,:) = tsn(:,:,2,:)
2755	zdep(:,:) = fsdept(:,:,2)
2756
2757	! store initial nn_eos
2758	jeos = nn_eos
2759	DO jn = 0, 2
2760	nn_eos = jn
2761	!=============================================================
2762	! 1) dx = ( T ) and dy = ( T )
2763	!=============================================================
2764
2765	!--------------------------------------------------------------------
2766	! Reset the tangent and adjoint variables
2767	!--------------------------------------------------------------------
2768	zts_tlin(:,:,:) = 0.0_wp
2769	zrd_tlout(:,:) = 0.0_wp
2770	zts_adout(:,:,:) = 0.0_wp
2771	zrd_adin(:,:) = 0.0_wp
2772
2773	!--------------------------------------------------------------------
2774	! Initialize the tangent input with random noise: dx
2775	!--------------------------------------------------------------------
2776	CALL grid_random( znts(:,:,jp_tem), 'T', 0.0_wp, stdt )
2777	CALL grid_random( znts(:,:,jp_sal), 'T', 0.0_wp, stds )
2778	DO jj = nldj, nlej
2779	DO ji = nldi, nlei
2780	zts_tlin(ji,jj,:) = znts(ji,jj,:)
2781	END DO
2782	END DO
2783
2784	CALL eos_insitu_2d_tan(zts, zdep, zts_tlin, zrd_tlout)
2785
2786	DO jj = nldj, nlej
2787	DO ji = nldi, nlei
2788	zrd_adin(ji,jj) = zrd_tlout(ji,jj) &
2789	& * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,2)&
2790	& * tmask(ji,jj,2)
2791	END DO
2792	END DO
2793
2794	!--------------------------------------------------------------------
2795	! Compute the scalar product: ( L dx )^T W dy
2796	!--------------------------------------------------------------------
2797
2798	zsp1 = DOT_PRODUCT( zrd_tlout, zrd_adin )
2799
2800	!--------------------------------------------------------------------
2801	! Call the adjoint routine: dx^* = L^T dy^*
2802	!--------------------------------------------------------------------
2803	CALL eos_insitu_2d_adj(zts, zdep, zts_adout, zrd_adin)
2804	zsp2 = DOT_PRODUCT( zts_tlin(:,:,jp_tem), zts_adout(:,:,jp_tem) ) + &
2805	& DOT_PRODUCT( zts_tlin(:,:,jp_sal), zts_adout(:,:,jp_sal) )
2806
2807	! Compare the scalar products
2808
2809	! 14 char:'12345678901234'
2810	SELECT CASE( jn )
2811	CASE (0) ; cl_name = 'eos_adj 2d T1'
2812	CASE (1) ; cl_name = 'eos_adj 2d T2'
2813	CASE (2) ; cl_name = 'eos_adj 2d T3'
2814	END SELECT
2815	CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
2816
2817	ENDDO
2818
2819	! restore initial nn_eos
2820	nn_eos = jeos
2821
2822	! Deallocate memory
2823
2824	DEALLOCATE( &
2825	& zdep, &
2826	& zts, &
2827	& zts_adout, &
2828	& zrd_adin, &
2829	& zts_tlin, &
2830	& zrd_tlout, &
2831	& znts )
2832	END SUBROUTINE eos_insitu_2d_adj_tst
2833
2834	SUBROUTINE eos_adj_tst( kumadt )
2835	!!-----------------------------------------------------------------------
2836	!!
2837	!! * ROUTINE eos_adj_tst *
2838	!!
2839	!! ** Purpose : Test the adjoint routine.
2840	!!
2841	!! History :
2842	!! ! 08-07 (A. Vidard)
2843	!!-----------------------------------------------------------------------
2844	!! * Arguments
2845	INTEGER, INTENT(IN) :: &
2846	& kumadt ! Output unit
2847
2848	CALL eos_insitu_adj_tst( kumadt )
2849	CALL eos_insitu_pot_adj_tst( kumadt )
2850	CALL eos_insitu_2d_adj_tst( kumadt )
2851	CALL eos_alpbet_adj_tst( kumadt )
2852
2853	END SUBROUTINE eos_adj_tst
2854
2855	SUBROUTINE bn2_adj_tst( kumadt )
2856	!!-----------------------------------------------------------------------
2857	!!
2858	!! * ROUTINE bn2_adj_tst *
2859	!!
2860	!! ** Purpose : Test the adjoint routine.
2861	!!
2862	!! ** Method : Verify the scalar product
2863	!!
2864	!! ( L dx )^T W dy = dx^T L^T W dy
2865	!!
2866	!! where L = tangent routine
2867	!! L^T = adjoint routine
2868	!! W = diagonal matrix of scale factors
2869	!! dx = input perturbation (random field)
2870	!! dy = L dx
2871	!!
2872	!! ** Action : Separate tests are applied for the following dx and dy:
2873	!!
2874	!! 1) dx = ( SSH ) and dy = ( SSH )
2875	!!
2876	!! History :
2877	!! ! 08-07 (A. Vidard)
2878	!!-----------------------------------------------------------------------
2879	!! * Modules used
2880
2881	!! * Arguments
2882	INTEGER, INTENT(IN) :: &
2883	& kumadt ! Output unit
2884	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
2885	zts ! potential temperature
2886	! salinity
2887	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
2888	& zts_adout ! potential temperature
2889	! salinity
2890	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2891	& zrd_adin, & ! potential density (surface referenced)
2892	& zrd_adout ! potential density (surface referenced)
2893	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
2894	& zts_tlin, & ! potential temperature
2895	! salinity
2896	& zts_tlout ! potential temperature
2897	! salinity
2898	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: &
2899	& zrd_tlout ! potential density (surface referenced)
2900	REAL(KIND=wp) :: &
2901	& zsp1, & ! scalar product involving the tangent routine
2902	& zsp2 ! scalar product involving the adjoint routine
2903	REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: &
2904	& znts ! potential temperature
2905	! salinity
2906	INTEGER :: &
2907	& ji, &
2908	& jj, &
2909	& jk, &
2910	& jn, &
2911	& jeos
2912	CHARACTER(LEN=14) :: cl_name
2913
2914	! Allocate memory
2915	ALLOCATE( &
2916	& zts( jpi, jpj, jpk, jpts ), &
2917	& zts_adout( jpi, jpj, jpk, jpts ), &
2918	& zrd_adin( jpi, jpj, jpk ), &
2919	& zrd_adout(jpi, jpj, jpk ), &
2920	& zts_tlin( jpi, jpj, jpk, jpts ), &
2921	& znts( jpi, jpj, jpk, jpts ), &
2922	& zts_tlout( jpi, jpj, jpk, jpts ), &
2923	& zrd_tlout(jpi, jpj, jpk ) )
2924
2925	! Initialize random field standard deviationsthe reference state
2926	zts = tsn
2927	! store initial nn_eos
2928	jeos = nn_eos
2929	DO jn = 0, 2
2930	nn_eos = jn
2931	!=============================================================
2932	! 1) dx = ( T ) and dy = ( T )
2933	!=============================================================
2934
2935	!--------------------------------------------------------------------
2936	! Reset the tangent and adjoint variables
2937	!--------------------------------------------------------------------
2938	zts_tlin(:,:,:,:) = 0.0_wp
2939	zts_tlout(:,:,:,:) = 0.0_wp
2940	zrd_tlout(:,:,:) = 0.0_wp
2941	zts_adout(:,:,:,:) = 0.0_wp
2942	zrd_adin(:,:,:) = 0.0_wp
2943	zrd_adout(:,:,:) = 0.0_wp
2944
2945	!--------------------------------------------------------------------
2946	! Initialize the tangent input with random noise: dx
2947	!--------------------------------------------------------------------
2948	CALL grid_random( znts(:,:,:,jp_tem), 'T', 0.0_wp, stdt )
2949	CALL grid_random( znts(:,:,:,jp_sal), 'T', 0.0_wp, stds )
2950	DO jk = 1, jpk
2951	DO jj = nldj, nlej
2952	DO ji = nldi, nlei
2953	zts_tlin(ji,jj,jk,:) = znts(ji,jj,jk,:)
2954	END DO
2955	END DO
2956	END DO
2957	!--------------------------------------------------------------------
2958	! Call the tangent routine: dy = L dx
2959	!--------------------------------------------------------------------
2960	zts_tlout(:,:,:,:) = zts_tlin
2961	CALL eos_bn2_tan( zts, zts_tlout, zrd_tlout )
2962	!--------------------------------------------------------------------
2963	! Initialize the adjoint variables: dy^* = W dy
2964	!--------------------------------------------------------------------
2965	DO jk = 1, jpk
2966	DO jj = nldj, nlej
2967	DO ji = nldi, nlei
2968	zrd_adin(ji,jj,jk) = zrd_tlout(ji,jj,jk) &
2969	& * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
2970	& * tmask(ji,jj,jk)
2971	END DO
2972	END DO
2973	END DO
2974
2975	!--------------------------------------------------------------------
2976	! Compute the scalar product: ( L dx )^T W dy
2977	!--------------------------------------------------------------------
2978
2979	zsp1 = DOT_PRODUCT( zrd_tlout, zrd_adin )
2980
2981	!--------------------------------------------------------------------
2982	! Call the adjoint routine: dx^* = L^T dy^*
2983	!--------------------------------------------------------------------
2984
2985	zrd_adout(:,:,:) = zrd_adin(:,:,:)
2986	CALL eos_bn2_adj( zts, zts_adout, zrd_adout )
2987
2988	!--------------------------------------------------------------------
2989	! Compute the scalar product: dx^T L^T W dy
2990	!--------------------------------------------------------------------
2991	zsp2 = DOT_PRODUCT( zts_tlin(:,:,:,jp_tem), zts_adout(:,:,:,jp_tem )) + &
2992	& DOT_PRODUCT( zts_tlin(:,:,:,jp_sal), zts_adout(:,:,:,jp_sal ))
2993
2994	! Compare the scalar products
2995	! 14 char:'12345678901234'
2996	SELECT CASE( jn )
2997	CASE (0) ; cl_name = 'bn2_adj T1'
2998	CASE (1) ; cl_name = 'bn2_adj T2'
2999	CASE (2) ; cl_name = 'bn2_adj T3'
3000	END SELECT
3001	CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
3002	ENDDO
3003	! restore initial nn_eos
3004	nn_eos = jeos
3005
3006	! Deallocate memory
3007	DEALLOCATE( &
3008	& zts, &
3009	& zts_adout, &
3010	& zrd_adin, &
3011	& zrd_adout, &
3012	& zts_tlin, &
3013	& zts_tlout, &
3014	& zrd_tlout, &
3015	& znts )
3016	END SUBROUTINE bn2_adj_tst
3017	#endif
3018	!!======================================================================
3019	END MODULE eosbn2_tam

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

Download in other formats: