Context Navigation

eosbn2.F90 @ 12340

Last change on this file since 12340 was 12340, checked in by acc, 4 years ago

Branch 2019/dev_r11943_MERGE_2019. This commit introduces basic do loop macro
substitution to the 2019 option 1, merge branch. These changes have been SETTE
tested. The only addition is the do_loop_substitute.h90 file in the OCE directory but
the macros defined therein are used throughout the code to replace identifiable, 2D-
and 3D- nested loop opening and closing statements with single-line alternatives. Code
indents are also adjusted accordingly.

The following explanation is taken from comments in the new header file:

This header file contains preprocessor definitions and macros used in the do-loop
substitutions introduced between version 4.0 and 4.2. The primary aim of these macros
is to assist in future applications of tiling to improve performance. This is expected
to be achieved by alternative versions of these macros in selected locations. The
initial introduction of these macros simply replaces all identifiable nested 2D- and
3D-loops with single line statements (and adjusts indenting accordingly). Do loops
are identifiable if they comform to either:

DO jk = ....

DO jj = .... DO jj = ...

DO ji = .... DO ji = ...
. OR .
. .

END DO END DO

END DO END DO

END DO

and white-space variants thereof.

Additionally, only loops with recognised jj and ji loops limits are treated; these are:
Lower limits of 1, 2 or fs_2
Upper limits of jpi, jpim1 or fs_jpim1 (for ji) or jpj, jpjm1 or fs_jpjm1 (for jj)

The macro naming convention takes the form: DO_2D_BT_LR where:

B is the Bottom offset from the PE's inner domain;
T is the Top offset from the PE's inner domain;
L is the Left offset from the PE's inner domain;
R is the Right offset from the PE's inner domain

So, given an inner domain of 2,jpim1 and 2,jpjm1, a typical example would replace:

DO jj = 2, jpj

DO ji = 1, jpim1
.
.

END DO

END DO

with:

DO_2D_01_10
.
.
END_2D

similar conventions apply to the 3D loops macros. jk loop limits are retained
through macro arguments and are not restricted. This includes the possibility of
strides for which an extra set of DO_3DS macros are defined.

In the example definition below the inner PE domain is defined by start indices of
(kIs, kJs) and end indices of (kIe, KJe)

#define DO_2D_00_00 DO jj = kJs, kJe ; DO ji = kIs, kIe
#define END_2D END DO ; END DO

TO DO:

Only conventional nested loops have been identified and replaced by this step. There are constructs such as:

DO jk = 2, jpkm1

z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)

END DO

which may need to be considered.

Property svn:keywords set to Id

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

Download in other formats: