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_v3_4_STABLE_2012/NEMOGCM/NEMO/OPATAM_SRC/TRA – NEMO

Context Navigation

source: branches/2012/dev_v3_4_STABLE_2012/NEMOGCM/NEMO/OPATAM_SRC/TRA/eosbn2_tam.F90 @ 4590

Last change on this file since 4590 was 4590, checked in by pabouttier, 10 years ago
Add a call to timing_stop at the end of eos_insistu_2d_tan, see Ticket #1281
Property svn:executable set to ``*
File size: 140.4 KB

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

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

Download in other formats: