Context Navigation

eosbn2.F90 @ 12377

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

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

Property svn:keywords set to Id

File size: 71.8 KB

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

Download in other formats: