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.F90 in trunk/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90 @ 7698

Last change on this file since 7698 was 7698, checked in by mocavero, 7 years ago
update trunk with OpenMP parallelization
Property svn:keywords set to `Id`
File size: 77.6 KB

Line
1	MODULE eosbn2
2	!!==============================================================================
3	!! * MODULE eosbn2 *
4	!! Equation Of Seawater : in situ density - Brunt-Vaisala frequency
5	!!==============================================================================
6	!! History : OPA ! 1989-03 (O. Marti) Original code
7	!! 6.0 ! 1994-07 (G. Madec, M. Imbard) add bn2
8	!! 6.0 ! 1994-08 (G. Madec) Add Jackett & McDougall eos
9	!! 7.0 ! 1996-01 (G. Madec) statement function for e3
10	!! 8.1 ! 1997-07 (G. Madec) density instead of volumic mass
11	!! - ! 1999-02 (G. Madec, N. Grima) semi-implicit pressure gradient
12	!! 8.2 ! 2001-09 (M. Ben Jelloul) bugfix on linear eos
13	!! NEMO 1.0 ! 2002-10 (G. Madec) add eos_init
14	!! - ! 2002-11 (G. Madec, A. Bozec) partial step, eos_insitu_2d
15	!! - ! 2003-08 (G. Madec) F90, free form
16	!! 3.0 ! 2006-08 (G. Madec) add tfreez function (now eos_fzp function)
17	!! 3.3 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA
18	!! - ! 2010-10 (G. Nurser, G. Madec) add alpha/beta used in ldfslp
19	!! 3.7 ! 2012-03 (F. Roquet, G. Madec) add primitive of alpha and beta used in PE computation
20	!! - ! 2012-05 (F. Roquet) add Vallis and original JM95 equation of state
21	!! - ! 2013-04 (F. Roquet, G. Madec) add eos_rab, change bn2 computation and reorganize the module
22	!! - ! 2014-09 (F. Roquet) add TEOS-10, S-EOS, and modify EOS-80
23	!! - ! 2015-06 (P.A. Bouttier) eos_fzp functions changed to subroutines for AGRIF
24	!!----------------------------------------------------------------------
25
26	!!----------------------------------------------------------------------
27	!! eos : generic interface of the equation of state
28	!! eos_insitu : Compute the in situ density
29	!! eos_insitu_pot: Compute the insitu and surface referenced potential volumic mass
30	!! eos_insitu_2d : Compute the in situ density for 2d fields
31	!! bn2 : Compute the Brunt-Vaisala frequency
32	!! eos_rab : generic interface of in situ thermal/haline expansion ratio
33	!! eos_rab_3d : compute in situ thermal/haline expansion ratio
34	!! eos_rab_2d : compute in situ thermal/haline expansion ratio for 2d fields
35	!! eos_fzp_2d : freezing temperature for 2d fields
36	!! eos_fzp_0d : freezing temperature for scalar
37	!! eos_init : set eos parameters (namelist)
38	!!----------------------------------------------------------------------
39	USE dom_oce ! ocean space and time domain
40	USE phycst ! physical constants
41	USE stopar ! Stochastic T/S fluctuations
42	USE stopts ! Stochastic T/S fluctuations
43	!
44	USE in_out_manager ! I/O manager
45	USE lib_mpp ! MPP library
46	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
47	USE prtctl ! Print control
48	USE wrk_nemo ! Memory Allocation
49	USE lbclnk ! ocean lateral boundary conditions
50	USE timing ! Timing
51
52	IMPLICIT NONE
53	PRIVATE
54
55	! !! * Interface
56	INTERFACE eos
57	MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d
58	END INTERFACE
59	!
60	INTERFACE eos_rab
61	MODULE PROCEDURE rab_3d, rab_2d, rab_0d
62	END INTERFACE
63	!
64	INTERFACE eos_fzp
65	MODULE PROCEDURE eos_fzp_2d, eos_fzp_0d
66	END INTERFACE
67	!
68	PUBLIC eos ! called by step, istate, tranpc and zpsgrd modules
69	PUBLIC bn2 ! called by step module
70	PUBLIC eos_rab ! called by ldfslp, zdfddm, trabbl
71	PUBLIC eos_pt_from_ct ! called by sbcssm
72	PUBLIC eos_fzp ! called by traadv_cen2 and sbcice_... modules
73	PUBLIC eos_pen ! used for pe diagnostics in trdpen module
74	PUBLIC eos_init ! called by istate module
75
76	! !! Namelist nameos
77	LOGICAL , PUBLIC :: ln_TEOS10 ! determine if eos_pt_from_ct is used to compute sst_m
78	LOGICAL , PUBLIC :: ln_EOS80 ! determine if eos_pt_from_ct is used to compute sst_m
79	LOGICAL , PUBLIC :: ln_SEOS ! determine if eos_pt_from_ct is used to compute sst_m
80
81	! Parameters
82	LOGICAL , PUBLIC :: l_useCT ! =T in ln_TEOS10=T (i.e. use eos_pt_from_ct to compute sst_m), =F otherwise
83	INTEGER , PUBLIC :: neos ! Identifier for equation of state used
84
85	INTEGER , PARAMETER :: np_teos10 = -1 ! parameter for using TEOS10
86	INTEGER , PARAMETER :: np_eos80 = 0 ! parameter for using EOS80
87	INTEGER , PARAMETER :: np_seos = 1 ! parameter for using Simplified Equation of state
88
89	! !!! simplified eos coefficients (default value: Vallis 2006)
90	REAL(wp) :: rn_a0 = 1.6550e-1_wp ! thermal expansion coeff.
91	REAL(wp) :: rn_b0 = 7.6554e-1_wp ! saline expansion coeff.
92	REAL(wp) :: rn_lambda1 = 5.9520e-2_wp ! cabbeling coeff. in T^2
93	REAL(wp) :: rn_lambda2 = 5.4914e-4_wp ! cabbeling coeff. in S^2
94	REAL(wp) :: rn_mu1 = 1.4970e-4_wp ! thermobaric coeff. in T
95	REAL(wp) :: rn_mu2 = 1.1090e-5_wp ! thermobaric coeff. in S
96	REAL(wp) :: rn_nu = 2.4341e-3_wp ! cabbeling coeff. in theta*salt
97
98	! TEOS10/EOS80 parameters
99	REAL(wp) :: r1_S0, r1_T0, r1_Z0, rdeltaS
100
101	! EOS parameters
102	REAL(wp) :: EOS000 , EOS100 , EOS200 , EOS300 , EOS400 , EOS500 , EOS600
103	REAL(wp) :: EOS010 , EOS110 , EOS210 , EOS310 , EOS410 , EOS510
104	REAL(wp) :: EOS020 , EOS120 , EOS220 , EOS320 , EOS420
105	REAL(wp) :: EOS030 , EOS130 , EOS230 , EOS330
106	REAL(wp) :: EOS040 , EOS140 , EOS240
107	REAL(wp) :: EOS050 , EOS150
108	REAL(wp) :: EOS060
109	REAL(wp) :: EOS001 , EOS101 , EOS201 , EOS301 , EOS401
110	REAL(wp) :: EOS011 , EOS111 , EOS211 , EOS311
111	REAL(wp) :: EOS021 , EOS121 , EOS221
112	REAL(wp) :: EOS031 , EOS131
113	REAL(wp) :: EOS041
114	REAL(wp) :: EOS002 , EOS102 , EOS202
115	REAL(wp) :: EOS012 , EOS112
116	REAL(wp) :: EOS022
117	REAL(wp) :: EOS003 , EOS103
118	REAL(wp) :: EOS013
119
120	! ALPHA parameters
121	REAL(wp) :: ALP000 , ALP100 , ALP200 , ALP300 , ALP400 , ALP500
122	REAL(wp) :: ALP010 , ALP110 , ALP210 , ALP310 , ALP410
123	REAL(wp) :: ALP020 , ALP120 , ALP220 , ALP320
124	REAL(wp) :: ALP030 , ALP130 , ALP230
125	REAL(wp) :: ALP040 , ALP140
126	REAL(wp) :: ALP050
127	REAL(wp) :: ALP001 , ALP101 , ALP201 , ALP301
128	REAL(wp) :: ALP011 , ALP111 , ALP211
129	REAL(wp) :: ALP021 , ALP121
130	REAL(wp) :: ALP031
131	REAL(wp) :: ALP002 , ALP102
132	REAL(wp) :: ALP012
133	REAL(wp) :: ALP003
134
135	! BETA parameters
136	REAL(wp) :: BET000 , BET100 , BET200 , BET300 , BET400 , BET500
137	REAL(wp) :: BET010 , BET110 , BET210 , BET310 , BET410
138	REAL(wp) :: BET020 , BET120 , BET220 , BET320
139	REAL(wp) :: BET030 , BET130 , BET230
140	REAL(wp) :: BET040 , BET140
141	REAL(wp) :: BET050
142	REAL(wp) :: BET001 , BET101 , BET201 , BET301
143	REAL(wp) :: BET011 , BET111 , BET211
144	REAL(wp) :: BET021 , BET121
145	REAL(wp) :: BET031
146	REAL(wp) :: BET002 , BET102
147	REAL(wp) :: BET012
148	REAL(wp) :: BET003
149
150	! PEN parameters
151	REAL(wp) :: PEN000 , PEN100 , PEN200 , PEN300 , PEN400
152	REAL(wp) :: PEN010 , PEN110 , PEN210 , PEN310
153	REAL(wp) :: PEN020 , PEN120 , PEN220
154	REAL(wp) :: PEN030 , PEN130
155	REAL(wp) :: PEN040
156	REAL(wp) :: PEN001 , PEN101 , PEN201
157	REAL(wp) :: PEN011 , PEN111
158	REAL(wp) :: PEN021
159	REAL(wp) :: PEN002 , PEN102
160	REAL(wp) :: PEN012
161
162	! ALPHA_PEN parameters
163	REAL(wp) :: APE000 , APE100 , APE200 , APE300
164	REAL(wp) :: APE010 , APE110 , APE210
165	REAL(wp) :: APE020 , APE120
166	REAL(wp) :: APE030
167	REAL(wp) :: APE001 , APE101
168	REAL(wp) :: APE011
169	REAL(wp) :: APE002
170
171	! BETA_PEN parameters
172	REAL(wp) :: BPE000 , BPE100 , BPE200 , BPE300
173	REAL(wp) :: BPE010 , BPE110 , BPE210
174	REAL(wp) :: BPE020 , BPE120
175	REAL(wp) :: BPE030
176	REAL(wp) :: BPE001 , BPE101
177	REAL(wp) :: BPE011
178	REAL(wp) :: BPE002
179
180	!! * Substitutions
181	# include "vectopt_loop_substitute.h90"
182	!!----------------------------------------------------------------------
183	!! NEMO/OPA 3.7 , NEMO Consortium (2014)
184	!! $Id$
185	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
186	!!----------------------------------------------------------------------
187	CONTAINS
188
189	SUBROUTINE eos_insitu( pts, prd, pdep )
190	!!----------------------------------------------------------------------
191	!! * ROUTINE eos_insitu *
192	!!
193	!! ** Purpose : Compute the in situ density (ratio rho/rau0) from
194	!! potential temperature and salinity using an equation of state
195	!! selected in the nameos namelist
196	!!
197	!! ** Method : prd(t,s,z) = ( rho(t,s,z) - rau0 ) / rau0
198	!! with prd in situ density anomaly no units
199	!! t TEOS10: CT or EOS80: PT Celsius
200	!! s TEOS10: SA or EOS80: SP TEOS10: g/kg or EOS80: psu
201	!! z depth meters
202	!! rho in situ density kg/m^3
203	!! rau0 reference density kg/m^3
204	!!
205	!! ln_teos10 : polynomial TEOS-10 equation of state is used for rho(t,s,z).
206	!! Check value: rho = 1028.21993233072 kg/m^3 for z=3000 dbar, ct=3 Celsius, sa=35.5 g/kg
207	!!
208	!! ln_eos80 : polynomial EOS-80 equation of state is used for rho(t,s,z).
209	!! Check value: rho = 1028.35011066567 kg/m^3 for z=3000 dbar, pt=3 Celsius, sp=35.5 psu
210	!!
211	!! ln_seos : simplified equation of state
212	!! prd(t,s,z) = ( -a0(1+lambda/2(T-T0)+muz+nu(S-S0))(T-T0) + b0(S-S0) ) / rau0
213	!! linear case function of T only: rn_alpha<>0, other coefficients = 0
214	!! linear eos function of T and S: rn_alpha and rn_beta<>0, other coefficients=0
215	!! Vallis like equation: use default values of coefficients
216	!!
217	!! ** Action : compute prd , the in situ density (no units)
218	!!
219	!! References : Roquet et al, Ocean Modelling, in preparation (2014)
220	!! Vallis, Atmospheric and Oceanic Fluid Dynamics, 2006
221	!! TEOS-10 Manual, 2010
222	!!----------------------------------------------------------------------
223	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! 1 : potential temperature [Celsius]
224	! ! 2 : salinity [psu]
225	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT( out) :: prd ! in situ density [-]
226	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pdep ! depth [m]
227	!
228	INTEGER :: ji, jj, jk ! dummy loop indices
229	REAL(wp) :: zt , zh , zs , ztm ! local scalars
230	REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - -
231	!!----------------------------------------------------------------------
232	!
233	IF( nn_timing == 1 ) CALL timing_start('eos-insitu')
234	!
235	SELECT CASE( neos )
236	!
237	CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==!
238	!
239	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zh, zt, zs, ztm, zn3, zn2, zn1, zn0, zn)
240	DO jk = 1, jpkm1
241	DO jj = 1, jpj
242	DO ji = 1, jpi
243	!
244	zh = pdep(ji,jj,jk) * r1_Z0 ! depth
245	zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature
246	zs = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity
247	ztm = tmask(ji,jj,jk) ! tmask
248	!
249	zn3 = EOS013*zt &
250	& + EOS103*zs+EOS003
251	!
252	zn2 = (EOS022*zt &
253	& + EOS112zs+EOS012)zt &
254	& + (EOS202zs+EOS102)zs+EOS002
255	!
256	zn1 = (((EOS041*zt &
257	& + EOS131zs+EOS031)zt &
258	& + (EOS221zs+EOS121)zs+EOS021)*zt &
259	& + ((EOS311zs+EOS211)zs+EOS111)zs+EOS011)zt &
260	& + (((EOS401zs+EOS301)zs+EOS201)zs+EOS101)zs+EOS001
261	!
262	zn0 = (((((EOS060*zt &
263	& + EOS150zs+EOS050)zt &
264	& + (EOS240zs+EOS140)zs+EOS040)*zt &
265	& + ((EOS330zs+EOS230)zs+EOS130)zs+EOS030)zt &
266	& + (((EOS420zs+EOS320)zs+EOS220)zs+EOS120)zs+EOS020)*zt &
267	& + ((((EOS510zs+EOS410)zs+EOS310)zs+EOS210)zs+EOS110)zs+EOS010)zt &
268	& + (((((EOS600zs+EOS500)zs+EOS400)zs+EOS300)zs+EOS200)zs+EOS100)zs+EOS000
269	!
270	zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
271	!
272	prd(ji,jj,jk) = ( zn * r1_rau0 - 1._wp ) * ztm ! density anomaly (masked)
273	!
274	END DO
275	END DO
276	END DO
277	!
278	CASE( np_seos ) !== simplified EOS ==!
279	!
280	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zh, zt, zs, ztm, zn)
281	DO jk = 1, jpkm1
282	DO jj = 1, jpj
283	DO ji = 1, jpi
284	zt = pts (ji,jj,jk,jp_tem) - 10._wp
285	zs = pts (ji,jj,jk,jp_sal) - 35._wp
286	zh = pdep (ji,jj,jk)
287	ztm = tmask(ji,jj,jk)
288	!
289	zn = - rn_a0 * ( 1._wp + 0.5_wprn_lambda1zt + rn_mu1zh ) zt &
290	& + rn_b0 * ( 1._wp - 0.5_wprn_lambda2zs - rn_mu2zh ) zs &
291	& - rn_nu * zt * zs
292	!
293	prd(ji,jj,jk) = zn * r1_rau0 * ztm ! density anomaly (masked)
294	END DO
295	END DO
296	END DO
297	!
298	END SELECT
299	!
300	IF(ln_ctl) CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-insitu : ', ovlap=1, kdim=jpk )
301	!
302	IF( nn_timing == 1 ) CALL timing_stop('eos-insitu')
303	!
304	END SUBROUTINE eos_insitu
305
306
307	SUBROUTINE eos_insitu_pot( pts, prd, prhop, pdep )
308	!!----------------------------------------------------------------------
309	!! * ROUTINE eos_insitu_pot *
310	!!
311	!! ** Purpose : Compute the in situ density (ratio rho/rau0) and the
312	!! potential volumic mass (Kg/m3) from potential temperature and
313	!! salinity fields using an equation of state selected in the
314	!! namelist.
315	!!
316	!! ** Action : - prd , the in situ density (no units)
317	!! - prhop, the potential volumic mass (Kg/m3)
318	!!
319	!!----------------------------------------------------------------------
320	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! 1 : potential temperature [Celsius]
321	! ! 2 : salinity [psu]
322	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT( out) :: prd ! in situ density [-]
323	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT( out) :: prhop ! potential density (surface referenced)
324	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pdep ! depth [m]
325	!
326	INTEGER :: ji, jj, jk, jsmp ! dummy loop indices
327	INTEGER :: jdof
328	REAL(wp) :: zt , zh , zstemp, zs , ztm ! local scalars
329	REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - -
330	REAL(wp), DIMENSION(:), ALLOCATABLE :: zn0_sto, zn_sto, zsign ! local vectors
331	!!----------------------------------------------------------------------
332	!
333	IF( nn_timing == 1 ) CALL timing_start('eos-pot')
334	!
335	SELECT CASE ( neos )
336	!
337	CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==!
338	!
339	! Stochastic equation of state
340	IF ( ln_sto_eos ) THEN
341	ALLOCATE(zn0_sto(1:2*nn_sto_eos))
342	ALLOCATE(zn_sto(1:2*nn_sto_eos))
343	ALLOCATE(zsign(1:2*nn_sto_eos))
344	DO jsmp = 1, 2*nn_sto_eos, 2
345	zsign(jsmp) = 1._wp
346	zsign(jsmp+1) = -1._wp
347	END DO
348	!
349	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, jsmp, jdof, zh, zt, zstemp, zs, ztm, zn3, zn2, zn1)
350	DO jk = 1, jpkm1
351	DO jj = 1, jpj
352	DO ji = 1, jpi
353	!
354	! compute density (2*nn_sto_eos) times:
355	! (1) for t+dt, s+ds (with the random TS fluctutation computed in sto_pts)
356	! (2) for t-dt, s-ds (with the opposite fluctuation)
357	DO jsmp = 1, nn_sto_eos*2
358	jdof = (jsmp + 1) / 2
359	zh = pdep(ji,jj,jk) * r1_Z0 ! depth
360	zt = (pts (ji,jj,jk,jp_tem) + pts_ran(ji,jj,jk,jp_tem,jdof) * zsign(jsmp)) * r1_T0 ! temperature
361	zstemp = pts (ji,jj,jk,jp_sal) + pts_ran(ji,jj,jk,jp_sal,jdof) * zsign(jsmp)
362	zs = SQRT( ABS( zstemp + rdeltaS ) * r1_S0 ) ! square root salinity
363	ztm = tmask(ji,jj,jk) ! tmask
364	!
365	zn3 = EOS013*zt &
366	& + EOS103*zs+EOS003
367	!
368	zn2 = (EOS022*zt &
369	& + EOS112zs+EOS012)zt &
370	& + (EOS202zs+EOS102)zs+EOS002
371	!
372	zn1 = (((EOS041*zt &
373	& + EOS131zs+EOS031)zt &
374	& + (EOS221zs+EOS121)zs+EOS021)*zt &
375	& + ((EOS311zs+EOS211)zs+EOS111)zs+EOS011)zt &
376	& + (((EOS401zs+EOS301)zs+EOS201)zs+EOS101)zs+EOS001
377	!
378	zn0_sto(jsmp) = (((((EOS060*zt &
379	& + EOS150zs+EOS050)zt &
380	& + (EOS240zs+EOS140)zs+EOS040)*zt &
381	& + ((EOS330zs+EOS230)zs+EOS130)zs+EOS030)zt &
382	& + (((EOS420zs+EOS320)zs+EOS220)zs+EOS120)zs+EOS020)*zt &
383	& + ((((EOS510zs+EOS410)zs+EOS310)zs+EOS210)zs+EOS110)zs+EOS010)zt &
384	& + (((((EOS600zs+EOS500)zs+EOS400)zs+EOS300)zs+EOS200)zs+EOS100)zs+EOS000
385	!
386	zn_sto(jsmp) = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0_sto(jsmp)
387	END DO
388	!
389	! compute stochastic density as the mean of the (2*nn_sto_eos) densities
390	prhop(ji,jj,jk) = 0._wp ; prd(ji,jj,jk) = 0._wp
391	DO jsmp = 1, nn_sto_eos*2
392	prhop(ji,jj,jk) = prhop(ji,jj,jk) + zn0_sto(jsmp) ! potential density referenced at the surface
393	!
394	prd(ji,jj,jk) = prd(ji,jj,jk) + ( zn_sto(jsmp) * r1_rau0 - 1._wp ) ! density anomaly (masked)
395	END DO
396	prhop(ji,jj,jk) = 0.5_wp * prhop(ji,jj,jk) * ztm / nn_sto_eos
397	prd (ji,jj,jk) = 0.5_wp * prd (ji,jj,jk) * ztm / nn_sto_eos
398	END DO
399	END DO
400	END DO
401	DEALLOCATE(zn0_sto,zn_sto,zsign)
402	! Non-stochastic equation of state
403	ELSE
404	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zh, zt, zs, ztm, zn3, zn2, zn1, zn0, zn)
405	DO jk = 1, jpkm1
406	DO jj = 1, jpj
407	DO ji = 1, jpi
408	!
409	zh = pdep(ji,jj,jk) * r1_Z0 ! depth
410	zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature
411	zs = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity
412	ztm = tmask(ji,jj,jk) ! tmask
413	!
414	zn3 = EOS013*zt &
415	& + EOS103*zs+EOS003
416	!
417	zn2 = (EOS022*zt &
418	& + EOS112zs+EOS012)zt &
419	& + (EOS202zs+EOS102)zs+EOS002
420	!
421	zn1 = (((EOS041*zt &
422	& + EOS131zs+EOS031)zt &
423	& + (EOS221zs+EOS121)zs+EOS021)*zt &
424	& + ((EOS311zs+EOS211)zs+EOS111)zs+EOS011)zt &
425	& + (((EOS401zs+EOS301)zs+EOS201)zs+EOS101)zs+EOS001
426	!
427	zn0 = (((((EOS060*zt &
428	& + EOS150zs+EOS050)zt &
429	& + (EOS240zs+EOS140)zs+EOS040)*zt &
430	& + ((EOS330zs+EOS230)zs+EOS130)zs+EOS030)zt &
431	& + (((EOS420zs+EOS320)zs+EOS220)zs+EOS120)zs+EOS020)*zt &
432	& + ((((EOS510zs+EOS410)zs+EOS310)zs+EOS210)zs+EOS110)zs+EOS010)zt &
433	& + (((((EOS600zs+EOS500)zs+EOS400)zs+EOS300)zs+EOS200)zs+EOS100)zs+EOS000
434	!
435	zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
436	!
437	prhop(ji,jj,jk) = zn0 * ztm ! potential density referenced at the surface
438	!
439	prd(ji,jj,jk) = ( zn * r1_rau0 - 1._wp ) * ztm ! density anomaly (masked)
440	END DO
441	END DO
442	END DO
443	ENDIF
444
445	CASE( np_seos ) !== simplified EOS ==!
446	!
447	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zh, zt, zs, ztm, zn)
448	DO jk = 1, jpkm1
449	DO jj = 1, jpj
450	DO ji = 1, jpi
451	zt = pts (ji,jj,jk,jp_tem) - 10._wp
452	zs = pts (ji,jj,jk,jp_sal) - 35._wp
453	zh = pdep (ji,jj,jk)
454	ztm = tmask(ji,jj,jk)
455	! ! potential density referenced at the surface
456	zn = - rn_a0 * ( 1._wp + 0.5_wprn_lambda1zt ) * zt &
457	& + rn_b0 * ( 1._wp - 0.5_wprn_lambda2zs ) * zs &
458	& - rn_nu * zt * zs
459	prhop(ji,jj,jk) = ( rau0 + zn ) * ztm
460	! ! density anomaly (masked)
461	zn = zn - ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zh
462	prd(ji,jj,jk) = zn * r1_rau0 * ztm
463	!
464	END DO
465	END DO
466	END DO
467	!
468	END SELECT
469	!
470	IF(ln_ctl) CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-pot: ', tab3d_2=prhop, clinfo2=' pot : ', ovlap=1, kdim=jpk )
471	!
472	IF( nn_timing == 1 ) CALL timing_stop('eos-pot')
473	!
474	END SUBROUTINE eos_insitu_pot
475
476
477	SUBROUTINE eos_insitu_2d( pts, pdep, prd )
478	!!----------------------------------------------------------------------
479	!! * ROUTINE eos_insitu_2d *
480	!!
481	!! ** Purpose : Compute the in situ density (ratio rho/rau0) from
482	!! potential temperature and salinity using an equation of state
483	!! selected in the nameos namelist. * 2D field case
484	!!
485	!! ** Action : - prd , the in situ density (no units) (unmasked)
486	!!
487	!!----------------------------------------------------------------------
488	REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in ) :: pts ! 1 : potential temperature [Celsius]
489	! ! 2 : salinity [psu]
490	REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pdep ! depth [m]
491	REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: prd ! in situ density
492	!
493	INTEGER :: ji, jj, jk ! dummy loop indices
494	REAL(wp) :: zt , zh , zs ! local scalars
495	REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - -
496	!!----------------------------------------------------------------------
497	!
498	IF( nn_timing == 1 ) CALL timing_start('eos2d')
499	!
500	!$OMP PARALLEL DO schedule(static) private(jj, ji)
501	DO jj = 1, jpj
502	DO ji = 1, jpi
503	prd(ji,jj) = 0._wp
504	END DO
505	END DO
506	!
507	SELECT CASE( neos )
508	!
509	CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==!
510	!
511	!$OMP PARALLEL DO schedule(static) private(jj, ji, zh, zt, zs, zn3, zn2, zn1, zn0, zn)
512	DO jj = 1, jpjm1
513	DO ji = 1, fs_jpim1 ! vector opt.
514	!
515	zh = pdep(ji,jj) * r1_Z0 ! depth
516	zt = pts (ji,jj,jp_tem) * r1_T0 ! temperature
517	zs = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity
518	!
519	zn3 = EOS013*zt &
520	& + EOS103*zs+EOS003
521	!
522	zn2 = (EOS022*zt &
523	& + EOS112zs+EOS012)zt &
524	& + (EOS202zs+EOS102)zs+EOS002
525	!
526	zn1 = (((EOS041*zt &
527	& + EOS131zs+EOS031)zt &
528	& + (EOS221zs+EOS121)zs+EOS021)*zt &
529	& + ((EOS311zs+EOS211)zs+EOS111)zs+EOS011)zt &
530	& + (((EOS401zs+EOS301)zs+EOS201)zs+EOS101)zs+EOS001
531	!
532	zn0 = (((((EOS060*zt &
533	& + EOS150zs+EOS050)zt &
534	& + (EOS240zs+EOS140)zs+EOS040)*zt &
535	& + ((EOS330zs+EOS230)zs+EOS130)zs+EOS030)zt &
536	& + (((EOS420zs+EOS320)zs+EOS220)zs+EOS120)zs+EOS020)*zt &
537	& + ((((EOS510zs+EOS410)zs+EOS310)zs+EOS210)zs+EOS110)zs+EOS010)zt &
538	& + (((((EOS600zs+EOS500)zs+EOS400)zs+EOS300)zs+EOS200)zs+EOS100)zs+EOS000
539	!
540	zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
541	!
542	prd(ji,jj) = zn * r1_rau0 - 1._wp ! unmasked in situ density anomaly
543	!
544	END DO
545	END DO
546	!
547	CALL lbc_lnk( prd, 'T', 1. ) ! Lateral boundary conditions
548	!
549	CASE( np_seos ) !== simplified EOS ==!
550	!
551	!$OMP PARALLEL DO schedule(static) private(jj, ji, zh, zt, zs, zn)
552	DO jj = 1, jpjm1
553	DO ji = 1, fs_jpim1 ! vector opt.
554	!
555	zt = pts (ji,jj,jp_tem) - 10._wp
556	zs = pts (ji,jj,jp_sal) - 35._wp
557	zh = pdep (ji,jj) ! depth at the partial step level
558	!
559	zn = - rn_a0 * ( 1._wp + 0.5_wprn_lambda1zt + rn_mu1zh ) zt &
560	& + rn_b0 * ( 1._wp - 0.5_wprn_lambda2zs - rn_mu2zh ) zs &
561	& - rn_nu * zt * zs
562	!
563	prd(ji,jj) = zn * r1_rau0 ! unmasked in situ density anomaly
564	!
565	END DO
566	END DO
567	!
568	CALL lbc_lnk( prd, 'T', 1. ) ! Lateral boundary conditions
569	!
570	END SELECT
571	!
572	IF(ln_ctl) CALL prt_ctl( tab2d_1=prd, clinfo1=' eos2d: ' )
573	!
574	IF( nn_timing == 1 ) CALL timing_stop('eos2d')
575	!
576	END SUBROUTINE eos_insitu_2d
577
578
579	SUBROUTINE rab_3d( pts, pab )
580	!!----------------------------------------------------------------------
581	!! * ROUTINE rab_3d *
582	!!
583	!! ** Purpose : Calculates thermal/haline expansion ratio at T-points
584	!!
585	!! ** Method : calculates alpha / beta at T-points
586	!!
587	!! ** Action : - pab : thermal/haline expansion ratio at T-points
588	!!----------------------------------------------------------------------
589	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature & salinity
590	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT( out) :: pab ! thermal/haline expansion ratio
591	!
592	INTEGER :: ji, jj, jk ! dummy loop indices
593	REAL(wp) :: zt , zh , zs , ztm ! local scalars
594	REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - -
595	!!----------------------------------------------------------------------
596	!
597	IF( nn_timing == 1 ) CALL timing_start('rab_3d')
598	!
599	SELECT CASE ( neos )
600	!
601	CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==!
602	!
603	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zh, zt, zs, ztm, zn3, zn2, zn1, zn0, zn)
604	DO jk = 1, jpkm1
605	DO jj = 1, jpj
606	DO ji = 1, jpi
607	!
608	zh = gdept_n(ji,jj,jk) * r1_Z0 ! depth
609	zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature
610	zs = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity
611	ztm = tmask(ji,jj,jk) ! tmask
612	!
613	! alpha
614	zn3 = ALP003
615	!
616	zn2 = ALP012zt + ALP102zs+ALP002
617	!
618	zn1 = ((ALP031*zt &
619	& + ALP121zs+ALP021)zt &
620	& + (ALP211zs+ALP111)zs+ALP011)*zt &
621	& + ((ALP301zs+ALP201)zs+ALP101)*zs+ALP001
622	!
623	zn0 = ((((ALP050*zt &
624	& + ALP140zs+ALP040)zt &
625	& + (ALP230zs+ALP130)zs+ALP030)*zt &
626	& + ((ALP320zs+ALP220)zs+ALP120)zs+ALP020)zt &
627	& + (((ALP410zs+ALP310)zs+ALP210)zs+ALP110)zs+ALP010)*zt &
628	& + ((((ALP500zs+ALP400)zs+ALP300)zs+ALP200)zs+ALP100)*zs+ALP000
629	!
630	zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
631	!
632	pab(ji,jj,jk,jp_tem) = zn * r1_rau0 * ztm
633	!
634	! beta
635	zn3 = BET003
636	!
637	zn2 = BET012zt + BET102zs+BET002
638	!
639	zn1 = ((BET031*zt &
640	& + BET121zs+BET021)zt &
641	& + (BET211zs+BET111)zs+BET011)*zt &
642	& + ((BET301zs+BET201)zs+BET101)*zs+BET001
643	!
644	zn0 = ((((BET050*zt &
645	& + BET140zs+BET040)zt &
646	& + (BET230zs+BET130)zs+BET030)*zt &
647	& + ((BET320zs+BET220)zs+BET120)zs+BET020)zt &
648	& + (((BET410zs+BET310)zs+BET210)zs+BET110)zs+BET010)*zt &
649	& + ((((BET500zs+BET400)zs+BET300)zs+BET200)zs+BET100)*zs+BET000
650	!
651	zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
652	!
653	pab(ji,jj,jk,jp_sal) = zn / zs * r1_rau0 * ztm
654	!
655	END DO
656	END DO
657	END DO
658	!
659	CASE( np_seos ) !== simplified EOS ==!
660	!
661	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zh, zt, zs, ztm, zn)
662	DO jk = 1, jpkm1
663	DO jj = 1, jpj
664	DO ji = 1, jpi
665	zt = pts (ji,jj,jk,jp_tem) - 10._wp ! pot. temperature anomaly (t-T0)
666	zs = pts (ji,jj,jk,jp_sal) - 35._wp ! abs. salinity anomaly (s-S0)
667	zh = gdept_n(ji,jj,jk) ! depth in meters at t-point
668	ztm = tmask(ji,jj,jk) ! land/sea bottom mask = surf. mask
669	!
670	zn = rn_a0 * ( 1._wp + rn_lambda1zt + rn_mu1zh ) + rn_nu*zs
671	pab(ji,jj,jk,jp_tem) = zn * r1_rau0 * ztm ! alpha
672	!
673	zn = rn_b0 * ( 1._wp - rn_lambda2zs - rn_mu2zh ) - rn_nu*zt
674	pab(ji,jj,jk,jp_sal) = zn * r1_rau0 * ztm ! beta
675	!
676	END DO
677	END DO
678	END DO
679	!
680	CASE DEFAULT
681	IF(lwp) WRITE(numout,cform_err)
682	IF(lwp) WRITE(numout,*) ' bad flag value for neos = ', neos
683	nstop = nstop + 1
684	!
685	END SELECT
686	!
687	IF(ln_ctl) CALL prt_ctl( tab3d_1=pab(:,:,:,jp_tem), clinfo1=' rab_3d_t: ', &
688	& tab3d_2=pab(:,:,:,jp_sal), clinfo2=' rab_3d_s : ', ovlap=1, kdim=jpk )
689	!
690	IF( nn_timing == 1 ) CALL timing_stop('rab_3d')
691	!
692	END SUBROUTINE rab_3d
693
694
695	SUBROUTINE rab_2d( pts, pdep, pab )
696	!!----------------------------------------------------------------------
697	!! * ROUTINE rab_2d *
698	!!
699	!! ** Purpose : Calculates thermal/haline expansion ratio for a 2d field (unmasked)
700	!!
701	!! ** Action : - pab : thermal/haline expansion ratio at T-points
702	!!----------------------------------------------------------------------
703	REAL(wp), DIMENSION(jpi,jpj,jpts) , INTENT(in ) :: pts ! pot. temperature & salinity
704	REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pdep ! depth [m]
705	REAL(wp), DIMENSION(jpi,jpj,jpts) , INTENT( out) :: pab ! thermal/haline expansion ratio
706	!
707	INTEGER :: ji, jj, jk ! dummy loop indices
708	REAL(wp) :: zt , zh , zs ! local scalars
709	REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - -
710	!!----------------------------------------------------------------------
711	!
712	IF( nn_timing == 1 ) CALL timing_start('rab_2d')
713	!
714	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
715	DO jk = 1, jpts
716	DO jj = 1, jpj
717	DO ji = 1, jpi
718	pab(ji,jj,jk) = 0._wp
719	END DO
720	END DO
721	END DO
722	!
723	SELECT CASE ( neos )
724	!
725	CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==!
726	!
727	!$OMP PARALLEL DO schedule(static) private(jj, ji, zh, zt, zs, zn3, zn2, zn1, zn0, zn)
728	DO jj = 1, jpjm1
729	DO ji = 1, fs_jpim1 ! vector opt.
730	!
731	zh = pdep(ji,jj) * r1_Z0 ! depth
732	zt = pts (ji,jj,jp_tem) * r1_T0 ! temperature
733	zs = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity
734	!
735	! alpha
736	zn3 = ALP003
737	!
738	zn2 = ALP012zt + ALP102zs+ALP002
739	!
740	zn1 = ((ALP031*zt &
741	& + ALP121zs+ALP021)zt &
742	& + (ALP211zs+ALP111)zs+ALP011)*zt &
743	& + ((ALP301zs+ALP201)zs+ALP101)*zs+ALP001
744	!
745	zn0 = ((((ALP050*zt &
746	& + ALP140zs+ALP040)zt &
747	& + (ALP230zs+ALP130)zs+ALP030)*zt &
748	& + ((ALP320zs+ALP220)zs+ALP120)zs+ALP020)zt &
749	& + (((ALP410zs+ALP310)zs+ALP210)zs+ALP110)zs+ALP010)*zt &
750	& + ((((ALP500zs+ALP400)zs+ALP300)zs+ALP200)zs+ALP100)*zs+ALP000
751	!
752	zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
753	!
754	pab(ji,jj,jp_tem) = zn * r1_rau0
755	!
756	! beta
757	zn3 = BET003
758	!
759	zn2 = BET012zt + BET102zs+BET002
760	!
761	zn1 = ((BET031*zt &
762	& + BET121zs+BET021)zt &
763	& + (BET211zs+BET111)zs+BET011)*zt &
764	& + ((BET301zs+BET201)zs+BET101)*zs+BET001
765	!
766	zn0 = ((((BET050*zt &
767	& + BET140zs+BET040)zt &
768	& + (BET230zs+BET130)zs+BET030)*zt &
769	& + ((BET320zs+BET220)zs+BET120)zs+BET020)zt &
770	& + (((BET410zs+BET310)zs+BET210)zs+BET110)zs+BET010)*zt &
771	& + ((((BET500zs+BET400)zs+BET300)zs+BET200)zs+BET100)*zs+BET000
772	!
773	zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
774	!
775	pab(ji,jj,jp_sal) = zn / zs * r1_rau0
776	!
777	!
778	END DO
779	END DO
780	!
781	CALL lbc_lnk( pab(:,:,jp_tem), 'T', 1. ) ! Lateral boundary conditions
782	CALL lbc_lnk( pab(:,:,jp_sal), 'T', 1. )
783	!
784	CASE( np_seos ) !== simplified EOS ==!
785	!
786	!$OMP PARALLEL DO schedule(static) private(jj, ji, zh, zt, zs, zn)
787	DO jj = 1, jpjm1
788	DO ji = 1, fs_jpim1 ! vector opt.
789	!
790	zt = pts (ji,jj,jp_tem) - 10._wp ! pot. temperature anomaly (t-T0)
791	zs = pts (ji,jj,jp_sal) - 35._wp ! abs. salinity anomaly (s-S0)
792	zh = pdep (ji,jj) ! depth at the partial step level
793	!
794	zn = rn_a0 * ( 1._wp + rn_lambda1zt + rn_mu1zh ) + rn_nu*zs
795	pab(ji,jj,jp_tem) = zn * r1_rau0 ! alpha
796	!
797	zn = rn_b0 * ( 1._wp - rn_lambda2zs - rn_mu2zh ) - rn_nu*zt
798	pab(ji,jj,jp_sal) = zn * r1_rau0 ! beta
799	!
800	END DO
801	END DO
802	!
803	CALL lbc_lnk( pab(:,:,jp_tem), 'T', 1. ) ! Lateral boundary conditions
804	CALL lbc_lnk( pab(:,:,jp_sal), 'T', 1. )
805	!
806	CASE DEFAULT
807	IF(lwp) WRITE(numout,cform_err)
808	IF(lwp) WRITE(numout,*) ' bad flag value for neos = ', neos
809	nstop = nstop + 1
810	!
811	END SELECT
812	!
813	IF(ln_ctl) CALL prt_ctl( tab2d_1=pab(:,:,jp_tem), clinfo1=' rab_2d_t: ', &
814	& tab2d_2=pab(:,:,jp_sal), clinfo2=' rab_2d_s : ' )
815	!
816	IF( nn_timing == 1 ) CALL timing_stop('rab_2d')
817	!
818	END SUBROUTINE rab_2d
819
820
821	SUBROUTINE rab_0d( pts, pdep, pab )
822	!!----------------------------------------------------------------------
823	!! * ROUTINE rab_0d *
824	!!
825	!! ** Purpose : Calculates thermal/haline expansion ratio for a 2d field (unmasked)
826	!!
827	!! ** Action : - pab : thermal/haline expansion ratio at T-points
828	!!----------------------------------------------------------------------
829	REAL(wp), DIMENSION(jpts) , INTENT(in ) :: pts ! pot. temperature & salinity
830	REAL(wp), INTENT(in ) :: pdep ! depth [m]
831	REAL(wp), DIMENSION(jpts) , INTENT( out) :: pab ! thermal/haline expansion ratio
832	!
833	REAL(wp) :: zt , zh , zs ! local scalars
834	REAL(wp) :: zn , zn0, zn1, zn2, zn3 ! - -
835	!!----------------------------------------------------------------------
836	!
837	IF( nn_timing == 1 ) CALL timing_start('rab_2d')
838	!
839	pab(:) = 0._wp
840	!
841	SELECT CASE ( neos )
842	!
843	CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==!
844	!
845	!
846	zh = pdep * r1_Z0 ! depth
847	zt = pts (jp_tem) * r1_T0 ! temperature
848	zs = SQRT( ABS( pts(jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity
849	!
850	! alpha
851	zn3 = ALP003
852	!
853	zn2 = ALP012zt + ALP102zs+ALP002
854	!
855	zn1 = ((ALP031*zt &
856	& + ALP121zs+ALP021)zt &
857	& + (ALP211zs+ALP111)zs+ALP011)*zt &
858	& + ((ALP301zs+ALP201)zs+ALP101)*zs+ALP001
859	!
860	zn0 = ((((ALP050*zt &
861	& + ALP140zs+ALP040)zt &
862	& + (ALP230zs+ALP130)zs+ALP030)*zt &
863	& + ((ALP320zs+ALP220)zs+ALP120)zs+ALP020)zt &
864	& + (((ALP410zs+ALP310)zs+ALP210)zs+ALP110)zs+ALP010)*zt &
865	& + ((((ALP500zs+ALP400)zs+ALP300)zs+ALP200)zs+ALP100)*zs+ALP000
866	!
867	zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
868	!
869	pab(jp_tem) = zn * r1_rau0
870	!
871	! beta
872	zn3 = BET003
873	!
874	zn2 = BET012zt + BET102zs+BET002
875	!
876	zn1 = ((BET031*zt &
877	& + BET121zs+BET021)zt &
878	& + (BET211zs+BET111)zs+BET011)*zt &
879	& + ((BET301zs+BET201)zs+BET101)*zs+BET001
880	!
881	zn0 = ((((BET050*zt &
882	& + BET140zs+BET040)zt &
883	& + (BET230zs+BET130)zs+BET030)*zt &
884	& + ((BET320zs+BET220)zs+BET120)zs+BET020)zt &
885	& + (((BET410zs+BET310)zs+BET210)zs+BET110)zs+BET010)*zt &
886	& + ((((BET500zs+BET400)zs+BET300)zs+BET200)zs+BET100)*zs+BET000
887	!
888	zn = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0
889	!
890	pab(jp_sal) = zn / zs * r1_rau0
891	!
892	!
893	!
894	CASE( np_seos ) !== simplified EOS ==!
895	!
896	zt = pts(jp_tem) - 10._wp ! pot. temperature anomaly (t-T0)
897	zs = pts(jp_sal) - 35._wp ! abs. salinity anomaly (s-S0)
898	zh = pdep ! depth at the partial step level
899	!
900	zn = rn_a0 * ( 1._wp + rn_lambda1zt + rn_mu1zh ) + rn_nu*zs
901	pab(jp_tem) = zn * r1_rau0 ! alpha
902	!
903	zn = rn_b0 * ( 1._wp - rn_lambda2zs - rn_mu2zh ) - rn_nu*zt
904	pab(jp_sal) = zn * r1_rau0 ! beta
905	!
906	CASE DEFAULT
907	IF(lwp) WRITE(numout,cform_err)
908	IF(lwp) WRITE(numout,*) ' bad flag value for neos = ', neos
909	nstop = nstop + 1
910	!
911	END SELECT
912	!
913	IF( nn_timing == 1 ) CALL timing_stop('rab_2d')
914	!
915	END SUBROUTINE rab_0d
916
917
918	SUBROUTINE bn2( pts, pab, pn2 )
919	!!----------------------------------------------------------------------
920	!! * ROUTINE bn2 *
921	!!
922	!! ** Purpose : Compute the local Brunt-Vaisala frequency at the
923	!! time-step of the input arguments
924	!!
925	!! ** Method : pn2 = grav * (alpha dk[T] + beta dk[S] ) / e3w
926	!! where alpha and beta are given in pab, and computed on T-points.
927	!! N.B. N^2 is set one for all to zero at jk=1 in istate module.
928	!!
929	!! ** Action : pn2 : square of the brunt-vaisala frequency at w-point
930	!!
931	!!----------------------------------------------------------------------
932	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature and salinity [Celsius,psu]
933	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pab ! thermal/haline expansion coef. [Celsius-1,psu-1]
934	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT( out) :: pn2 ! Brunt-Vaisala frequency squared [1/s^2]
935	!
936	INTEGER :: ji, jj, jk ! dummy loop indices
937	REAL(wp) :: zaw, zbw, zrw ! local scalars
938	!!----------------------------------------------------------------------
939	!
940	IF( nn_timing == 1 ) CALL timing_start('bn2')
941	!
942	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zrw, zaw, zbw)
943	DO jk = 2, jpkm1 ! interior points only (2=< jk =< jpkm1 )
944	DO jj = 1, jpj ! surface and bottom value set to zero one for all in istate.F90
945	DO ji = 1, jpi
946	zrw = ( gdepw_n(ji,jj,jk ) - gdept_n(ji,jj,jk) ) &
947	& / ( gdept_n(ji,jj,jk-1) - gdept_n(ji,jj,jk) )
948	!
949	zaw = pab(ji,jj,jk,jp_tem) * (1. - zrw) + pab(ji,jj,jk-1,jp_tem) * zrw
950	zbw = pab(ji,jj,jk,jp_sal) * (1. - zrw) + pab(ji,jj,jk-1,jp_sal) * zrw
951	!
952	pn2(ji,jj,jk) = grav * ( zaw * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) &
953	& - zbw * ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ) &
954	& / e3w_n(ji,jj,jk) * tmask(ji,jj,jk)
955	END DO
956	END DO
957	END DO
958	!
959	IF(ln_ctl) CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2 : ', ovlap=1, kdim=jpk )
960	!
961	IF( nn_timing == 1 ) CALL timing_stop('bn2')
962	!
963	END SUBROUTINE bn2
964
965
966	FUNCTION eos_pt_from_ct( ctmp, psal ) RESULT( ptmp )
967	!!----------------------------------------------------------------------
968	!! * ROUTINE eos_pt_from_ct *
969	!!
970	!! ** Purpose : Compute pot.temp. from cons. temp. [Celsius]
971	!!
972	!! ** Method : rational approximation (5/3th order) of TEOS-10 algorithm
973	!! checkvalue: pt=20.02391895 Celsius for sa=35.7g/kg, ct=20degC
974	!!
975	!! Reference : TEOS-10, UNESCO
976	!! Rational approximation to TEOS10 algorithm (rms error on WOA13 values: 4.0e-5 degC)
977	!!----------------------------------------------------------------------
978	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: ctmp ! Cons. Temp [Celsius]
979	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: psal ! salinity [psu]
980	! Leave result array automatic rather than making explicitly allocated
981	REAL(wp), DIMENSION(jpi,jpj) :: ptmp ! potential temperature [Celsius]
982	!
983	INTEGER :: ji, jj ! dummy loop indices
984	REAL(wp) :: zt , zs , ztm ! local scalars
985	REAL(wp) :: zn , zd ! local scalars
986	REAL(wp) :: zdeltaS , z1_S0 , z1_T0
987	!!----------------------------------------------------------------------
988	!
989	IF ( nn_timing == 1 ) CALL timing_start('eos_pt_from_ct')
990	!
991	zdeltaS = 5._wp
992	z1_S0 = 0.875_wp/35.16504_wp
993	z1_T0 = 1._wp/40._wp
994	!
995	!$OMP PARALLEL DO schedule(static) private(jj, ji, zt, zs, ztm, zn,zd)
996	DO jj = 1, jpj
997	DO ji = 1, jpi
998	!
999	zt = ctmp (ji,jj) * z1_T0
1000	zs = SQRT( ABS( psal(ji,jj) + zdeltaS ) * r1_S0 )
1001	ztm = tmask(ji,jj,1)
1002	!
1003	zn = ((((-2.1385727895e-01_wp*zt &
1004	& - 2.7674419971e-01_wpzs+1.0728094330_wp)zt &
1005	& + (2.6366564313_wpzs+3.3546960647_wp)zs-7.8012209473_wp)*zt &
1006	& + ((1.8835586562_wpzs+7.3949191679_wp)zs-3.3937395875_wp)zs-5.6414948432_wp)zt &
1007	& + (((3.5737370589_wpzs-1.5512427389e+01_wp)zs+2.4625741105e+01_wp)*zs &
1008	& +1.9912291000e+01_wp)zs-3.2191146312e+01_wp)zt &
1009	& + ((((5.7153204649e-01_wpzs-3.0943149543_wp)zs+9.3052495181_wp)*zs &
1010	& -9.4528934807_wp)zs+3.1066408996_wp)zs-4.3504021262e-01_wp
1011	!
1012	zd = (2.0035003456_wp*zt &
1013	& -3.4570358592e-01_wpzs+5.6471810638_wp)zt &
1014	& + (1.5393993508_wpzs-6.9394762624_wp)zs+1.2750522650e+01_wp
1015	!
1016	ptmp(ji,jj) = ( zt / z1_T0 + zn / zd ) * ztm
1017	!
1018	END DO
1019	END DO
1020	!
1021	IF( nn_timing == 1 ) CALL timing_stop('eos_pt_from_ct')
1022	!
1023	END FUNCTION eos_pt_from_ct
1024
1025
1026	SUBROUTINE eos_fzp_2d( psal, ptf, pdep )
1027	!!----------------------------------------------------------------------
1028	!! * ROUTINE eos_fzp *
1029	!!
1030	!! ** Purpose : Compute the freezing point temperature [Celsius]
1031	!!
1032	!! ** Method : UNESCO freezing point (ptf) in Celsius is given by
1033	!! ptf(t,z) = (-.0575+1.710523e-3sqrt(abs(s))-2.154996e-4s)s - 7.53e-4z
1034	!! checkvalue: tf=-2.588567 Celsius for s=40psu, z=500m
1035	!!
1036	!! Reference : UNESCO tech. papers in the marine science no. 28. 1978
1037	!!----------------------------------------------------------------------
1038	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: psal ! salinity [psu]
1039	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ), OPTIONAL :: pdep ! depth [m]
1040	REAL(wp), DIMENSION(jpi,jpj), INTENT(out ) :: ptf ! freezing temperature [Celsius]
1041	!
1042	INTEGER :: ji, jj ! dummy loop indices
1043	REAL(wp) :: zt, zs, z1_S0 ! local scalars
1044	!!----------------------------------------------------------------------
1045	!
1046	SELECT CASE ( neos )
1047	!
1048	CASE ( np_teos10, np_seos ) !== CT,SA (TEOS-10 and S-EOS formulations) ==!
1049	!
1050	z1_S0 = 1._wp / 35.16504_wp
1051	!$OMP PARALLEL
1052	!$OMP DO schedule(static) private(jj, ji, zs)
1053	DO jj = 1, jpj
1054	DO ji = 1, jpi
1055	zs= SQRT( ABS( psal(ji,jj) ) * z1_S0 ) ! square root salinity
1056	ptf(ji,jj) = ((((1.46873e-03_wpzs-9.64972e-03_wp)zs+2.28348e-02_wp)*zs &
1057	& - 3.12775e-02_wp)zs+2.07679e-02_wp)zs-5.87701e-02_wp
1058	END DO
1059	END DO
1060	!$OMP DO schedule(static) private(jj, ji)
1061	DO jj = 1, jpj
1062	DO ji = 1, jpi
1063	ptf(ji,jj) = ptf(ji,jj) * psal(ji,jj)
1064	END DO
1065	END DO
1066	!$OMP END PARALLEL
1067	!
1068	IF( PRESENT( pdep ) ) THEN
1069	!$OMP PARALLEL DO schedule(static) private(jj, ji)
1070	DO jj = 1, jpj
1071	DO ji = 1, jpi
1072	ptf(ji,jj) = ptf(ji,jj) - 7.53e-4 * pdep(ji,jj)
1073	END DO
1074	END DO
1075	END IF
1076	!
1077	CASE ( np_eos80 ) !== PT,SP (UNESCO formulation) ==!
1078	!
1079	!$OMP PARALLEL DO schedule(static) private(jj, ji)
1080	DO jj = 1, jpj
1081	DO ji = 1, jpi
1082	ptf(ji,jj) = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal(ji,jj) ) &
1083	& - 2.154996e-4_wp * psal(ji,jj) ) * psal(ji,jj)
1084	END DO
1085	END DO
1086	!
1087	IF( PRESENT( pdep ) ) THEN
1088	!$OMP PARALLEL DO schedule(static) private(jj, ji)
1089	DO jj = 1, jpj
1090	DO ji = 1, jpi
1091	ptf(ji,jj) = ptf(ji,jj) - 7.53e-4 * pdep(ji,jj)
1092	END DO
1093	END DO
1094	END IF
1095	!
1096	CASE DEFAULT
1097	IF(lwp) WRITE(numout,cform_err)
1098	IF(lwp) WRITE(numout,*) ' bad flag value for neos = ', neos
1099	nstop = nstop + 1
1100	!
1101	END SELECT
1102	!
1103	END SUBROUTINE eos_fzp_2d
1104
1105
1106	SUBROUTINE eos_fzp_0d( psal, ptf, pdep )
1107	!!----------------------------------------------------------------------
1108	!! * ROUTINE eos_fzp *
1109	!!
1110	!! ** Purpose : Compute the freezing point temperature [Celsius]
1111	!!
1112	!! ** Method : UNESCO freezing point (ptf) in Celsius is given by
1113	!! ptf(t,z) = (-.0575+1.710523e-3sqrt(abs(s))-2.154996e-4s)s - 7.53e-4z
1114	!! checkvalue: tf=-2.588567 Celsius for s=40psu, z=500m
1115	!!
1116	!! Reference : UNESCO tech. papers in the marine science no. 28. 1978
1117	!!----------------------------------------------------------------------
1118	REAL(wp), INTENT(in ) :: psal ! salinity [psu]
1119	REAL(wp), INTENT(in ), OPTIONAL :: pdep ! depth [m]
1120	REAL(wp), INTENT(out) :: ptf ! freezing temperature [Celsius]
1121	!
1122	REAL(wp) :: zs ! local scalars
1123	!!----------------------------------------------------------------------
1124	!
1125	SELECT CASE ( neos )
1126	!
1127	CASE ( np_teos10, np_seos ) !== CT,SA (TEOS-10 and S-EOS formulations) ==!
1128	!
1129	zs = SQRT( ABS( psal ) / 35.16504_wp ) ! square root salinity
1130	ptf = ((((1.46873e-03_wpzs-9.64972e-03_wp)zs+2.28348e-02_wp)*zs &
1131	& - 3.12775e-02_wp)zs+2.07679e-02_wp)zs-5.87701e-02_wp
1132	ptf = ptf * psal
1133	!
1134	IF( PRESENT( pdep ) ) ptf = ptf - 7.53e-4 * pdep
1135	!
1136	CASE ( np_eos80 ) !== PT,SP (UNESCO formulation) ==!
1137	!
1138	ptf = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal ) &
1139	& - 2.154996e-4_wp * psal ) * psal
1140	!
1141	IF( PRESENT( pdep ) ) ptf = ptf - 7.53e-4 * pdep
1142	!
1143	CASE DEFAULT
1144	IF(lwp) WRITE(numout,cform_err)
1145	IF(lwp) WRITE(numout,*) ' bad flag value for neos = ', neos
1146	nstop = nstop + 1
1147	!
1148	END SELECT
1149	!
1150	END SUBROUTINE eos_fzp_0d
1151
1152
1153	SUBROUTINE eos_pen( pts, pab_pe, ppen )
1154	!!----------------------------------------------------------------------
1155	!! * ROUTINE eos_pen *
1156	!!
1157	!! ** Purpose : Calculates nonlinear anomalies of alpha_PE, beta_PE and PE at T-points
1158	!!
1159	!! ** Method : PE is defined analytically as the vertical
1160	!! primitive of EOS times -g integrated between 0 and z>0.
1161	!! pen is the nonlinear bsq-PE anomaly: pen = ( PE - rau0 gz ) / rau0 gz - rd
1162	!! = 1/z * /int_0^z rd dz - rd
1163	!! where rd is the density anomaly (see eos_rhd function)
1164	!! ab_pe are partial derivatives of PE anomaly with respect to T and S:
1165	!! ab_pe(1) = - 1/(rau0 gz) * dPE/dT + drd/dT = - d(pen)/dT
1166	!! ab_pe(2) = 1/(rau0 gz) * dPE/dS + drd/dS = d(pen)/dS
1167	!!
1168	!! ** Action : - pen : PE anomaly given at T-points
1169	!! : - pab_pe : given at T-points
1170	!! pab_pe(:,:,:,jp_tem) is alpha_pe
1171	!! pab_pe(:,:,:,jp_sal) is beta_pe
1172	!!----------------------------------------------------------------------
1173	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature & salinity
1174	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT( out) :: pab_pe ! alpha_pe and beta_pe
1175	REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: ppen ! potential energy anomaly
1176	!
1177	INTEGER :: ji, jj, jk ! dummy loop indices
1178	REAL(wp) :: zt , zh , zs , ztm ! local scalars
1179	REAL(wp) :: zn , zn0, zn1, zn2 ! - -
1180	!!----------------------------------------------------------------------
1181	!
1182	IF( nn_timing == 1 ) CALL timing_start('eos_pen')
1183	!
1184	SELECT CASE ( neos )
1185	!
1186	CASE( np_teos10, np_eos80 ) !== polynomial TEOS-10 / EOS-80 ==!
1187	!
1188	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zh, zt, zs, ztm, zn2, zn1, zn0, zn)
1189	DO jk = 1, jpkm1
1190	DO jj = 1, jpj
1191	DO ji = 1, jpi
1192	!
1193	zh = gdept_n(ji,jj,jk) * r1_Z0 ! depth
1194	zt = pts (ji,jj,jk,jp_tem) * r1_T0 ! temperature
1195	zs = SQRT( ABS( pts(ji,jj,jk,jp_sal) + rdeltaS ) * r1_S0 ) ! square root salinity
1196	ztm = tmask(ji,jj,jk) ! tmask
1197	!
1198	! potential energy non-linear anomaly
1199	zn2 = (PEN012)*zt &
1200	& + PEN102*zs+PEN002
1201	!
1202	zn1 = ((PEN021)*zt &
1203	& + PEN111zs+PEN011)zt &
1204	& + (PEN201zs+PEN101)zs+PEN001
1205	!
1206	zn0 = ((((PEN040)*zt &
1207	& + PEN130zs+PEN030)zt &
1208	& + (PEN220zs+PEN120)zs+PEN020)*zt &
1209	& + ((PEN310zs+PEN210)zs+PEN110)zs+PEN010)zt &
1210	& + (((PEN400zs+PEN300)zs+PEN200)zs+PEN100)zs+PEN000
1211	!
1212	zn = ( zn2 * zh + zn1 ) * zh + zn0
1213	!
1214	ppen(ji,jj,jk) = zn * zh * r1_rau0 * ztm
1215	!
1216	! alphaPE non-linear anomaly
1217	zn2 = APE002
1218	!
1219	zn1 = (APE011)*zt &
1220	& + APE101*zs+APE001
1221	!
1222	zn0 = (((APE030)*zt &
1223	& + APE120zs+APE020)zt &
1224	& + (APE210zs+APE110)zs+APE010)*zt &
1225	& + ((APE300zs+APE200)zs+APE100)*zs+APE000
1226	!
1227	zn = ( zn2 * zh + zn1 ) * zh + zn0
1228	!
1229	pab_pe(ji,jj,jk,jp_tem) = zn * zh * r1_rau0 * ztm
1230	!
1231	! betaPE non-linear anomaly
1232	zn2 = BPE002
1233	!
1234	zn1 = (BPE011)*zt &
1235	& + BPE101*zs+BPE001
1236	!
1237	zn0 = (((BPE030)*zt &
1238	& + BPE120zs+BPE020)zt &
1239	& + (BPE210zs+BPE110)zs+BPE010)*zt &
1240	& + ((BPE300zs+BPE200)zs+BPE100)*zs+BPE000
1241	!
1242	zn = ( zn2 * zh + zn1 ) * zh + zn0
1243	!
1244	pab_pe(ji,jj,jk,jp_sal) = zn / zs * zh * r1_rau0 * ztm
1245	!
1246	END DO
1247	END DO
1248	END DO
1249	!
1250	CASE( np_seos ) !== Vallis (2006) simplified EOS ==!
1251	!
1252	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zh, zt, zs, ztm, zn)
1253	DO jk = 1, jpkm1
1254	DO jj = 1, jpj
1255	DO ji = 1, jpi
1256	zt = pts(ji,jj,jk,jp_tem) - 10._wp ! temperature anomaly (t-T0)
1257	zs = pts (ji,jj,jk,jp_sal) - 35._wp ! abs. salinity anomaly (s-S0)
1258	zh = gdept_n(ji,jj,jk) ! depth in meters at t-point
1259	ztm = tmask(ji,jj,jk) ! tmask
1260	zn = 0.5_wp * zh * r1_rau0 * ztm
1261	! ! Potential Energy
1262	ppen(ji,jj,jk) = ( rn_a0 * rn_mu1 * zt + rn_b0 * rn_mu2 * zs ) * zn
1263	! ! alphaPE
1264	pab_pe(ji,jj,jk,jp_tem) = - rn_a0 * rn_mu1 * zn
1265	pab_pe(ji,jj,jk,jp_sal) = rn_b0 * rn_mu2 * zn
1266	!
1267	END DO
1268	END DO
1269	END DO
1270	!
1271	CASE DEFAULT
1272	IF(lwp) WRITE(numout,cform_err)
1273	IF(lwp) WRITE(numout,*) ' bad flag value for neos = ', neos
1274	nstop = nstop + 1
1275	!
1276	END SELECT
1277	!
1278	IF( nn_timing == 1 ) CALL timing_stop('eos_pen')
1279	!
1280	END SUBROUTINE eos_pen
1281
1282
1283	SUBROUTINE eos_init
1284	!!----------------------------------------------------------------------
1285	!! * ROUTINE eos_init *
1286	!!
1287	!! ** Purpose : initializations for the equation of state
1288	!!
1289	!! ** Method : Read the namelist nameos and control the parameters
1290	!!----------------------------------------------------------------------
1291	INTEGER :: ios ! local integer
1292	INTEGER :: ioptio ! local integer
1293	!!
1294	NAMELIST/nameos/ ln_TEOS10, ln_EOS80, ln_SEOS, rn_a0, rn_b0, rn_lambda1, rn_mu1, &
1295	& rn_lambda2, rn_mu2, rn_nu
1296	!!----------------------------------------------------------------------
1297	!
1298	REWIND( numnam_ref ) ! Namelist nameos in reference namelist : equation of state
1299	READ ( numnam_ref, nameos, IOSTAT = ios, ERR = 901 )
1300	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nameos in reference namelist', lwp )
1301	!
1302	REWIND( numnam_cfg ) ! Namelist nameos in configuration namelist : equation of state
1303	READ ( numnam_cfg, nameos, IOSTAT = ios, ERR = 902 )
1304	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nameos in configuration namelist', lwp )
1305	IF(lwm) WRITE( numond, nameos )
1306	!
1307	rau0 = 1026._wp !: volumic mass of reference [kg/m3]
1308	rcp = 3991.86795711963_wp !: heat capacity [J/K]
1309	!
1310	IF(lwp) THEN ! Control print
1311	WRITE(numout,*)
1312	WRITE(numout,*) 'eos_init : equation of state'
1313	WRITE(numout,*) '~~~~~~~~'
1314	WRITE(numout,*) ' Namelist nameos : Chosen the Equation Of Seawater (EOS)'
1315	WRITE(numout,*) ' TEOS-10 : rho=F(Conservative Temperature, Absolute Salinity, depth) ln_TEOS10 = ', ln_TEOS10
1316	WRITE(numout,*) ' EOS-80 : rho=F(Potential Temperature, Practical Salinity, depth) ln_EOS80 = ', ln_EOS80
1317	WRITE(numout,*) ' S-EOS : rho=F(Conservative Temperature, Absolute Salinity, depth) ln_SEOS = ', ln_SEOS
1318	ENDIF
1319
1320	! Check options for equation of state & set neos based on logical flags
1321	ioptio = 0
1322	IF( ln_TEOS10 ) THEN ; ioptio = ioptio+1 ; neos = np_teos10 ; ENDIF
1323	IF( ln_EOS80 ) THEN ; ioptio = ioptio+1 ; neos = np_eos80 ; ENDIF
1324	IF( ln_SEOS ) THEN ; ioptio = ioptio+1 ; neos = np_seos ; ENDIF
1325	IF( ioptio /= 1 ) CALL ctl_stop("Exactly one equation of state option must be selected")
1326	!
1327	SELECT CASE( neos ) ! check option
1328	!
1329	CASE( np_teos10 ) !== polynomial TEOS-10 ==!
1330	IF(lwp) WRITE(numout,*)
1331	IF(lwp) WRITE(numout,*) ' use of TEOS-10 equation of state (cons. temp. and abs. salinity)'
1332	!
1333	l_useCT = .TRUE. ! model temperature is Conservative temperature
1334	!
1335	rdeltaS = 32._wp
1336	r1_S0 = 0.875_wp/35.16504_wp
1337	r1_T0 = 1._wp/40._wp
1338	r1_Z0 = 1.e-4_wp
1339	!
1340	EOS000 = 8.0189615746e+02_wp
1341	EOS100 = 8.6672408165e+02_wp
1342	EOS200 = -1.7864682637e+03_wp
1343	EOS300 = 2.0375295546e+03_wp
1344	EOS400 = -1.2849161071e+03_wp
1345	EOS500 = 4.3227585684e+02_wp
1346	EOS600 = -6.0579916612e+01_wp
1347	EOS010 = 2.6010145068e+01_wp
1348	EOS110 = -6.5281885265e+01_wp
1349	EOS210 = 8.1770425108e+01_wp
1350	EOS310 = -5.6888046321e+01_wp
1351	EOS410 = 1.7681814114e+01_wp
1352	EOS510 = -1.9193502195_wp
1353	EOS020 = -3.7074170417e+01_wp
1354	EOS120 = 6.1548258127e+01_wp
1355	EOS220 = -6.0362551501e+01_wp
1356	EOS320 = 2.9130021253e+01_wp
1357	EOS420 = -5.4723692739_wp
1358	EOS030 = 2.1661789529e+01_wp
1359	EOS130 = -3.3449108469e+01_wp
1360	EOS230 = 1.9717078466e+01_wp
1361	EOS330 = -3.1742946532_wp
1362	EOS040 = -8.3627885467_wp
1363	EOS140 = 1.1311538584e+01_wp
1364	EOS240 = -5.3563304045_wp
1365	EOS050 = 5.4048723791e-01_wp
1366	EOS150 = 4.8169980163e-01_wp
1367	EOS060 = -1.9083568888e-01_wp
1368	EOS001 = 1.9681925209e+01_wp
1369	EOS101 = -4.2549998214e+01_wp
1370	EOS201 = 5.0774768218e+01_wp
1371	EOS301 = -3.0938076334e+01_wp
1372	EOS401 = 6.6051753097_wp
1373	EOS011 = -1.3336301113e+01_wp
1374	EOS111 = -4.4870114575_wp
1375	EOS211 = 5.0042598061_wp
1376	EOS311 = -6.5399043664e-01_wp
1377	EOS021 = 6.7080479603_wp
1378	EOS121 = 3.5063081279_wp
1379	EOS221 = -1.8795372996_wp
1380	EOS031 = -2.4649669534_wp
1381	EOS131 = -5.5077101279e-01_wp
1382	EOS041 = 5.5927935970e-01_wp
1383	EOS002 = 2.0660924175_wp
1384	EOS102 = -4.9527603989_wp
1385	EOS202 = 2.5019633244_wp
1386	EOS012 = 2.0564311499_wp
1387	EOS112 = -2.1311365518e-01_wp
1388	EOS022 = -1.2419983026_wp
1389	EOS003 = -2.3342758797e-02_wp
1390	EOS103 = -1.8507636718e-02_wp
1391	EOS013 = 3.7969820455e-01_wp
1392	!
1393	ALP000 = -6.5025362670e-01_wp
1394	ALP100 = 1.6320471316_wp
1395	ALP200 = -2.0442606277_wp
1396	ALP300 = 1.4222011580_wp
1397	ALP400 = -4.4204535284e-01_wp
1398	ALP500 = 4.7983755487e-02_wp
1399	ALP010 = 1.8537085209_wp
1400	ALP110 = -3.0774129064_wp
1401	ALP210 = 3.0181275751_wp
1402	ALP310 = -1.4565010626_wp
1403	ALP410 = 2.7361846370e-01_wp
1404	ALP020 = -1.6246342147_wp
1405	ALP120 = 2.5086831352_wp
1406	ALP220 = -1.4787808849_wp
1407	ALP320 = 2.3807209899e-01_wp
1408	ALP030 = 8.3627885467e-01_wp
1409	ALP130 = -1.1311538584_wp
1410	ALP230 = 5.3563304045e-01_wp
1411	ALP040 = -6.7560904739e-02_wp
1412	ALP140 = -6.0212475204e-02_wp
1413	ALP050 = 2.8625353333e-02_wp
1414	ALP001 = 3.3340752782e-01_wp
1415	ALP101 = 1.1217528644e-01_wp
1416	ALP201 = -1.2510649515e-01_wp
1417	ALP301 = 1.6349760916e-02_wp
1418	ALP011 = -3.3540239802e-01_wp
1419	ALP111 = -1.7531540640e-01_wp
1420	ALP211 = 9.3976864981e-02_wp
1421	ALP021 = 1.8487252150e-01_wp
1422	ALP121 = 4.1307825959e-02_wp
1423	ALP031 = -5.5927935970e-02_wp
1424	ALP002 = -5.1410778748e-02_wp
1425	ALP102 = 5.3278413794e-03_wp
1426	ALP012 = 6.2099915132e-02_wp
1427	ALP003 = -9.4924551138e-03_wp
1428	!
1429	BET000 = 1.0783203594e+01_wp
1430	BET100 = -4.4452095908e+01_wp
1431	BET200 = 7.6048755820e+01_wp
1432	BET300 = -6.3944280668e+01_wp
1433	BET400 = 2.6890441098e+01_wp
1434	BET500 = -4.5221697773_wp
1435	BET010 = -8.1219372432e-01_wp
1436	BET110 = 2.0346663041_wp
1437	BET210 = -2.1232895170_wp
1438	BET310 = 8.7994140485e-01_wp
1439	BET410 = -1.1939638360e-01_wp
1440	BET020 = 7.6574242289e-01_wp
1441	BET120 = -1.5019813020_wp
1442	BET220 = 1.0872489522_wp
1443	BET320 = -2.7233429080e-01_wp
1444	BET030 = -4.1615152308e-01_wp
1445	BET130 = 4.9061350869e-01_wp
1446	BET230 = -1.1847737788e-01_wp
1447	BET040 = 1.4073062708e-01_wp
1448	BET140 = -1.3327978879e-01_wp
1449	BET050 = 5.9929880134e-03_wp
1450	BET001 = -5.2937873009e-01_wp
1451	BET101 = 1.2634116779_wp
1452	BET201 = -1.1547328025_wp
1453	BET301 = 3.2870876279e-01_wp
1454	BET011 = -5.5824407214e-02_wp
1455	BET111 = 1.2451933313e-01_wp
1456	BET211 = -2.4409539932e-02_wp
1457	BET021 = 4.3623149752e-02_wp
1458	BET121 = -4.6767901790e-02_wp
1459	BET031 = -6.8523260060e-03_wp
1460	BET002 = -6.1618945251e-02_wp
1461	BET102 = 6.2255521644e-02_wp
1462	BET012 = -2.6514181169e-03_wp
1463	BET003 = -2.3025968587e-04_wp
1464	!
1465	PEN000 = -9.8409626043_wp
1466	PEN100 = 2.1274999107e+01_wp
1467	PEN200 = -2.5387384109e+01_wp
1468	PEN300 = 1.5469038167e+01_wp
1469	PEN400 = -3.3025876549_wp
1470	PEN010 = 6.6681505563_wp
1471	PEN110 = 2.2435057288_wp
1472	PEN210 = -2.5021299030_wp
1473	PEN310 = 3.2699521832e-01_wp
1474	PEN020 = -3.3540239802_wp
1475	PEN120 = -1.7531540640_wp
1476	PEN220 = 9.3976864981e-01_wp
1477	PEN030 = 1.2324834767_wp
1478	PEN130 = 2.7538550639e-01_wp
1479	PEN040 = -2.7963967985e-01_wp
1480	PEN001 = -1.3773949450_wp
1481	PEN101 = 3.3018402659_wp
1482	PEN201 = -1.6679755496_wp
1483	PEN011 = -1.3709540999_wp
1484	PEN111 = 1.4207577012e-01_wp
1485	PEN021 = 8.2799886843e-01_wp
1486	PEN002 = 1.7507069098e-02_wp
1487	PEN102 = 1.3880727538e-02_wp
1488	PEN012 = -2.8477365341e-01_wp
1489	!
1490	APE000 = -1.6670376391e-01_wp
1491	APE100 = -5.6087643219e-02_wp
1492	APE200 = 6.2553247576e-02_wp
1493	APE300 = -8.1748804580e-03_wp
1494	APE010 = 1.6770119901e-01_wp
1495	APE110 = 8.7657703198e-02_wp
1496	APE210 = -4.6988432490e-02_wp
1497	APE020 = -9.2436260751e-02_wp
1498	APE120 = -2.0653912979e-02_wp
1499	APE030 = 2.7963967985e-02_wp
1500	APE001 = 3.4273852498e-02_wp
1501	APE101 = -3.5518942529e-03_wp
1502	APE011 = -4.1399943421e-02_wp
1503	APE002 = 7.1193413354e-03_wp
1504	!
1505	BPE000 = 2.6468936504e-01_wp
1506	BPE100 = -6.3170583896e-01_wp
1507	BPE200 = 5.7736640125e-01_wp
1508	BPE300 = -1.6435438140e-01_wp
1509	BPE010 = 2.7912203607e-02_wp
1510	BPE110 = -6.2259666565e-02_wp
1511	BPE210 = 1.2204769966e-02_wp
1512	BPE020 = -2.1811574876e-02_wp
1513	BPE120 = 2.3383950895e-02_wp
1514	BPE030 = 3.4261630030e-03_wp
1515	BPE001 = 4.1079296834e-02_wp
1516	BPE101 = -4.1503681096e-02_wp
1517	BPE011 = 1.7676120780e-03_wp
1518	BPE002 = 1.7269476440e-04_wp
1519	!
1520	CASE( np_eos80 ) !== polynomial EOS-80 formulation ==!
1521	!
1522	IF(lwp) WRITE(numout,*)
1523	IF(lwp) WRITE(numout,*) ' use of EOS-80 equation of state (pot. temp. and pract. salinity)'
1524	!
1525	l_useCT = .FALSE. ! model temperature is Potential temperature
1526	rdeltaS = 20._wp
1527	r1_S0 = 1._wp/40._wp
1528	r1_T0 = 1._wp/40._wp
1529	r1_Z0 = 1.e-4_wp
1530	!
1531	EOS000 = 9.5356891948e+02_wp
1532	EOS100 = 1.7136499189e+02_wp
1533	EOS200 = -3.7501039454e+02_wp
1534	EOS300 = 5.1856810420e+02_wp
1535	EOS400 = -3.7264470465e+02_wp
1536	EOS500 = 1.4302533998e+02_wp
1537	EOS600 = -2.2856621162e+01_wp
1538	EOS010 = 1.0087518651e+01_wp
1539	EOS110 = -1.3647741861e+01_wp
1540	EOS210 = 8.8478359933_wp
1541	EOS310 = -7.2329388377_wp
1542	EOS410 = 1.4774410611_wp
1543	EOS510 = 2.0036720553e-01_wp
1544	EOS020 = -2.5579830599e+01_wp
1545	EOS120 = 2.4043512327e+01_wp
1546	EOS220 = -1.6807503990e+01_wp
1547	EOS320 = 8.3811577084_wp
1548	EOS420 = -1.9771060192_wp
1549	EOS030 = 1.6846451198e+01_wp
1550	EOS130 = -2.1482926901e+01_wp
1551	EOS230 = 1.0108954054e+01_wp
1552	EOS330 = -6.2675951440e-01_wp
1553	EOS040 = -8.0812310102_wp
1554	EOS140 = 1.0102374985e+01_wp
1555	EOS240 = -4.8340368631_wp
1556	EOS050 = 1.2079167803_wp
1557	EOS150 = 1.1515380987e-01_wp
1558	EOS060 = -2.4520288837e-01_wp
1559	EOS001 = 1.0748601068e+01_wp
1560	EOS101 = -1.7817043500e+01_wp
1561	EOS201 = 2.2181366768e+01_wp
1562	EOS301 = -1.6750916338e+01_wp
1563	EOS401 = 4.1202230403_wp
1564	EOS011 = -1.5852644587e+01_wp
1565	EOS111 = -7.6639383522e-01_wp
1566	EOS211 = 4.1144627302_wp
1567	EOS311 = -6.6955877448e-01_wp
1568	EOS021 = 9.9994861860_wp
1569	EOS121 = -1.9467067787e-01_wp
1570	EOS221 = -1.2177554330_wp
1571	EOS031 = -3.4866102017_wp
1572	EOS131 = 2.2229155620e-01_wp
1573	EOS041 = 5.9503008642e-01_wp
1574	EOS002 = 1.0375676547_wp
1575	EOS102 = -3.4249470629_wp
1576	EOS202 = 2.0542026429_wp
1577	EOS012 = 2.1836324814_wp
1578	EOS112 = -3.4453674320e-01_wp
1579	EOS022 = -1.2548163097_wp
1580	EOS003 = 1.8729078427e-02_wp
1581	EOS103 = -5.7238495240e-02_wp
1582	EOS013 = 3.8306136687e-01_wp
1583	!
1584	ALP000 = -2.5218796628e-01_wp
1585	ALP100 = 3.4119354654e-01_wp
1586	ALP200 = -2.2119589983e-01_wp
1587	ALP300 = 1.8082347094e-01_wp
1588	ALP400 = -3.6936026529e-02_wp
1589	ALP500 = -5.0091801383e-03_wp
1590	ALP010 = 1.2789915300_wp
1591	ALP110 = -1.2021756164_wp
1592	ALP210 = 8.4037519952e-01_wp
1593	ALP310 = -4.1905788542e-01_wp
1594	ALP410 = 9.8855300959e-02_wp
1595	ALP020 = -1.2634838399_wp
1596	ALP120 = 1.6112195176_wp
1597	ALP220 = -7.5817155402e-01_wp
1598	ALP320 = 4.7006963580e-02_wp
1599	ALP030 = 8.0812310102e-01_wp
1600	ALP130 = -1.0102374985_wp
1601	ALP230 = 4.8340368631e-01_wp
1602	ALP040 = -1.5098959754e-01_wp
1603	ALP140 = -1.4394226233e-02_wp
1604	ALP050 = 3.6780433255e-02_wp
1605	ALP001 = 3.9631611467e-01_wp
1606	ALP101 = 1.9159845880e-02_wp
1607	ALP201 = -1.0286156825e-01_wp
1608	ALP301 = 1.6738969362e-02_wp
1609	ALP011 = -4.9997430930e-01_wp
1610	ALP111 = 9.7335338937e-03_wp
1611	ALP211 = 6.0887771651e-02_wp
1612	ALP021 = 2.6149576513e-01_wp
1613	ALP121 = -1.6671866715e-02_wp
1614	ALP031 = -5.9503008642e-02_wp
1615	ALP002 = -5.4590812035e-02_wp
1616	ALP102 = 8.6134185799e-03_wp
1617	ALP012 = 6.2740815484e-02_wp
1618	ALP003 = -9.5765341718e-03_wp
1619	!
1620	BET000 = 2.1420623987_wp
1621	BET100 = -9.3752598635_wp
1622	BET200 = 1.9446303907e+01_wp
1623	BET300 = -1.8632235232e+01_wp
1624	BET400 = 8.9390837485_wp
1625	BET500 = -1.7142465871_wp
1626	BET010 = -1.7059677327e-01_wp
1627	BET110 = 2.2119589983e-01_wp
1628	BET210 = -2.7123520642e-01_wp
1629	BET310 = 7.3872053057e-02_wp
1630	BET410 = 1.2522950346e-02_wp
1631	BET020 = 3.0054390409e-01_wp
1632	BET120 = -4.2018759976e-01_wp
1633	BET220 = 3.1429341406e-01_wp
1634	BET320 = -9.8855300959e-02_wp
1635	BET030 = -2.6853658626e-01_wp
1636	BET130 = 2.5272385134e-01_wp
1637	BET230 = -2.3503481790e-02_wp
1638	BET040 = 1.2627968731e-01_wp
1639	BET140 = -1.2085092158e-01_wp
1640	BET050 = 1.4394226233e-03_wp
1641	BET001 = -2.2271304375e-01_wp
1642	BET101 = 5.5453416919e-01_wp
1643	BET201 = -6.2815936268e-01_wp
1644	BET301 = 2.0601115202e-01_wp
1645	BET011 = -9.5799229402e-03_wp
1646	BET111 = 1.0286156825e-01_wp
1647	BET211 = -2.5108454043e-02_wp
1648	BET021 = -2.4333834734e-03_wp
1649	BET121 = -3.0443885826e-02_wp
1650	BET031 = 2.7786444526e-03_wp
1651	BET002 = -4.2811838287e-02_wp
1652	BET102 = 5.1355066072e-02_wp
1653	BET012 = -4.3067092900e-03_wp
1654	BET003 = -7.1548119050e-04_wp
1655	!
1656	PEN000 = -5.3743005340_wp
1657	PEN100 = 8.9085217499_wp
1658	PEN200 = -1.1090683384e+01_wp
1659	PEN300 = 8.3754581690_wp
1660	PEN400 = -2.0601115202_wp
1661	PEN010 = 7.9263222935_wp
1662	PEN110 = 3.8319691761e-01_wp
1663	PEN210 = -2.0572313651_wp
1664	PEN310 = 3.3477938724e-01_wp
1665	PEN020 = -4.9997430930_wp
1666	PEN120 = 9.7335338937e-02_wp
1667	PEN220 = 6.0887771651e-01_wp
1668	PEN030 = 1.7433051009_wp
1669	PEN130 = -1.1114577810e-01_wp
1670	PEN040 = -2.9751504321e-01_wp
1671	PEN001 = -6.9171176978e-01_wp
1672	PEN101 = 2.2832980419_wp
1673	PEN201 = -1.3694684286_wp
1674	PEN011 = -1.4557549876_wp
1675	PEN111 = 2.2969116213e-01_wp
1676	PEN021 = 8.3654420645e-01_wp
1677	PEN002 = -1.4046808820e-02_wp
1678	PEN102 = 4.2928871430e-02_wp
1679	PEN012 = -2.8729602515e-01_wp
1680	!
1681	APE000 = -1.9815805734e-01_wp
1682	APE100 = -9.5799229402e-03_wp
1683	APE200 = 5.1430784127e-02_wp
1684	APE300 = -8.3694846809e-03_wp
1685	APE010 = 2.4998715465e-01_wp
1686	APE110 = -4.8667669469e-03_wp
1687	APE210 = -3.0443885826e-02_wp
1688	APE020 = -1.3074788257e-01_wp
1689	APE120 = 8.3359333577e-03_wp
1690	APE030 = 2.9751504321e-02_wp
1691	APE001 = 3.6393874690e-02_wp
1692	APE101 = -5.7422790533e-03_wp
1693	APE011 = -4.1827210323e-02_wp
1694	APE002 = 7.1824006288e-03_wp
1695	!
1696	BPE000 = 1.1135652187e-01_wp
1697	BPE100 = -2.7726708459e-01_wp
1698	BPE200 = 3.1407968134e-01_wp
1699	BPE300 = -1.0300557601e-01_wp
1700	BPE010 = 4.7899614701e-03_wp
1701	BPE110 = -5.1430784127e-02_wp
1702	BPE210 = 1.2554227021e-02_wp
1703	BPE020 = 1.2166917367e-03_wp
1704	BPE120 = 1.5221942913e-02_wp
1705	BPE030 = -1.3893222263e-03_wp
1706	BPE001 = 2.8541225524e-02_wp
1707	BPE101 = -3.4236710714e-02_wp
1708	BPE011 = 2.8711395266e-03_wp
1709	BPE002 = 5.3661089288e-04_wp
1710	!
1711	CASE( np_seos ) !== Simplified EOS ==!
1712	IF(lwp) THEN
1713	WRITE(numout,*)
1714	WRITE(numout,*) ' use of simplified eos: rhd(dT=T-10,dS=S-35,Z) = '
1715	WRITE(numout,) ' [-a0(1+lambda1/2dT+mu1Z)dT + b0(1+lambda2/2dT+mu2Z)dS - nudT*dS]/rau0'
1716	WRITE(numout,*)
1717	WRITE(numout,*) ' thermal exp. coef. rn_a0 = ', rn_a0
1718	WRITE(numout,*) ' saline cont. coef. rn_b0 = ', rn_b0
1719	WRITE(numout,*) ' cabbeling coef. rn_lambda1 = ', rn_lambda1
1720	WRITE(numout,*) ' cabbeling coef. rn_lambda2 = ', rn_lambda2
1721	WRITE(numout,*) ' thermobar. coef. rn_mu1 = ', rn_mu1
1722	WRITE(numout,*) ' thermobar. coef. rn_mu2 = ', rn_mu2
1723	WRITE(numout,*) ' 2nd cabbel. coef. rn_nu = ', rn_nu
1724	WRITE(numout,*) ' Caution: rn_beta0=0 incompatible with ddm parameterization '
1725	ENDIF
1726	l_useCT = .TRUE. ! Use conservative temperature
1727	!
1728	CASE DEFAULT !== ERROR in neos ==!
1729	WRITE(ctmp1,*) ' bad flag value for neos = ', neos, '. You should never see this error'
1730	CALL ctl_stop( ctmp1 )
1731	!
1732	END SELECT
1733	!
1734	rau0_rcp = rau0 * rcp
1735	r1_rau0 = 1._wp / rau0
1736	r1_rcp = 1._wp / rcp
1737	r1_rau0_rcp = 1._wp / rau0_rcp
1738	!
1739	IF(lwp) THEN
1740	IF( l_useCT ) THEN
1741	WRITE(numout,*) ' model uses Conservative Temperature'
1742	WRITE(numout,*) ' Important: model must be initialized with CT and SA fields'
1743	ELSE
1744	WRITE(numout,*) ' model does not use Conservative Temperature'
1745	ENDIF
1746	ENDIF
1747	!
1748	IF(lwp) WRITE(numout,*)
1749	IF(lwp) WRITE(numout,*) ' volumic mass of reference rau0 = ', rau0 , ' kg/m^3'
1750	IF(lwp) WRITE(numout,*) ' 1. / rau0 r1_rau0 = ', r1_rau0, ' m^3/kg'
1751	IF(lwp) WRITE(numout,*) ' ocean specific heat rcp = ', rcp , ' J/Kelvin'
1752	IF(lwp) WRITE(numout,) ' rau0 rcp rau0_rcp = ', rau0_rcp
1753	IF(lwp) WRITE(numout,) ' 1. / ( rau0 rcp ) r1_rau0_rcp = ', r1_rau0_rcp
1754	!
1755	END SUBROUTINE eos_init
1756
1757	!!======================================================================
1758	END MODULE eosbn2

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

Download in other formats: