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/TAM_V3_0/NEMOTAM/OPATAM_SRC – NEMO

Context Navigation

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

Last change on this file since 1885 was 1885, checked in by rblod, 14 years ago
add TAM sources
File size: 184.1 KB

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

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

Download in other formats: