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

Last change on this file since 3611 was 3611, checked in by pabouttier, 11 years ago
Add TAM code and ORCA2_TAM configuration - see Ticket #1007
Property svn:executable set to ``*
File size: 140.3 KB

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

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

Download in other formats: