Context Navigation

eosbn2.F90 @ 9383

Last change on this file since 9383 was 9383, checked in by andmirek, 6 years ago
#2050 fixes and changes
File size: 74.6 KB

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

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/UKMO/test_moci_test_suite_namelist_read/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90 @ 9383

Download in other formats: