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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

eosbn2.F90 in branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90 @ 4666

Last change on this file since 4666 was 4666, checked in by mathiot, 10 years ago
#1331 : add ISOMIP config files + ice shelf code
Property svn:keywords set to `Id`
File size: 40.6 KB

Rev	Line
[3]	1	MODULE eosbn2
	2	!!==============================================================================
	3	!! * MODULE eosbn2 *
	4	!! Ocean diagnostic variable : equation of state - in situ and potential density
[3294]	5	!! - Brunt-Vaisala frequency
[3]	6	!!==============================================================================
[1559]	7	!! History : OPA ! 1989-03 (O. Marti) Original code
	8	!! 6.0 ! 1994-07 (G. Madec, M. Imbard) add bn2
	9	!! 6.0 ! 1994-08 (G. Madec) Add Jackett & McDougall eos
	10	!! 7.0 ! 1996-01 (G. Madec) statement function for e3
	11	!! 8.1 ! 1997-07 (G. Madec) density instead of volumic mass
	12	!! - ! 1999-02 (G. Madec, N. Grima) semi-implicit pressure gradient
	13	!! 8.2 ! 2001-09 (M. Ben Jelloul) bugfix on linear eos
	14	!! NEMO 1.0 ! 2002-10 (G. Madec) add eos_init
	15	!! - ! 2002-11 (G. Madec, A. Bozec) partial step, eos_insitu_2d
	16	!! - ! 2003-08 (G. Madec) F90, free form
	17	!! 3.0 ! 2006-08 (G. Madec) add tfreez function
[2528]	18	!! 3.3 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA
	19	!! - ! 2010-10 (G. Nurser, G. Madec) add eos_alpbet used in ldfslp
[888]	20	!!----------------------------------------------------------------------
[3]	21
	22	!!----------------------------------------------------------------------
	23	!! eos : generic interface of the equation of state
	24	!! eos_insitu : Compute the in situ density
	25	!! eos_insitu_pot : Compute the insitu and surface referenced potential
	26	!! volumic mass
	27	!! eos_insitu_2d : Compute the in situ density for 2d fields
	28	!! eos_bn2 : Compute the Brunt-Vaisala frequency
[3294]	29	!! eos_alpbet : calculates the in situ thermal/haline expansion ratio
[888]	30	!! tfreez : Compute the surface freezing temperature
[3]	31	!! eos_init : set eos parameters (namelist)
	32	!!----------------------------------------------------------------------
	33	USE dom_oce ! ocean space and time domain
	34	USE phycst ! physical constants
[2715]	35	USE zdfddm ! vertical physics: double diffusion
[3]	36	USE in_out_manager ! I/O manager
[2715]	37	USE lib_mpp ! MPP library
[258]	38	USE prtctl ! Print control
[3294]	39	USE wrk_nemo ! Memory Allocation
	40	USE timing ! Timing
[3]	41
	42	IMPLICIT NONE
	43	PRIVATE
	44
[3294]	45	! !! * Interface
[3]	46	INTERFACE eos
	47	MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d
[3294]	48	END INTERFACE
[3]	49	INTERFACE bn2
	50	MODULE PROCEDURE eos_bn2
[3294]	51	END INTERFACE
[3]	52
[888]	53	PUBLIC eos ! called by step, istate, tranpc and zpsgrd modules
[988]	54	PUBLIC eos_init ! called by istate module
[888]	55	PUBLIC bn2 ! called by step module
[2528]	56	PUBLIC eos_alpbet ! called by ldfslp module
[4666]	57	PUBLIC tfreez ! called by sbcice_... modules and sbcisf module
	58	PUBLIC tfreez1D ! called by trasbc modules
[3]	59
[4147]	60	! !!* Namelist (nameos) *
	61	INTEGER , PUBLIC :: nn_eos !: = 0/1/2 type of eq. of state and Brunt-Vaisala frequ.
	62	REAL(wp), PUBLIC :: rn_alpha !: thermal expension coeff. (linear equation of state)
	63	REAL(wp), PUBLIC :: rn_beta !: saline expension coeff. (linear equation of state)
[1601]	64
[2528]	65	REAL(wp), PUBLIC :: ralpbet !: alpha / beta ratio
[3294]	66
[3]	67	!! * Substitutions
	68	# include "domzgr_substitute.h90"
	69	# include "vectopt_loop_substitute.h90"
	70	!!----------------------------------------------------------------------
[2528]	71	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[888]	72	!! $Id$
[2528]	73	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	74	!!----------------------------------------------------------------------
	75	CONTAINS
	76
[4292]	77	SUBROUTINE eos_insitu( pts, prd, pdep )
[3]	78	!!----------------------------------------------------------------------
	79	!! * ROUTINE eos_insitu *
[3294]	80	!!
	81	!! ** Purpose : Compute the in situ density (ratio rho/rau0) from
[3]	82	!! potential temperature and salinity using an equation of state
[1601]	83	!! defined through the namelist parameter nn_eos.
[3]	84	!!
	85	!! ** Method : 3 cases:
[1601]	86	!! nn_eos = 0 : Jackett and McDougall (1994) equation of state.
[3]	87	!! the in situ density is computed directly as a function of
	88	!! potential temperature relative to the surface (the opa t
	89	!! variable), salt and pressure (assuming no pressure variation
	90	!! along geopotential surfaces, i.e. the pressure p in decibars
	91	!! is approximated by the depth in meters.
	92	!! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0
	93	!! with pressure p decibars
	94	!! potential temperature t deg celsius
	95	!! salinity s psu
	96	!! reference volumic mass rau0 kg/m**3
	97	!! in situ volumic mass rho kg/m**3
	98	!! in situ density anomalie prd no units
	99	!! Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar,
	100	!! t = 40 deg celcius, s=40 psu
[1601]	101	!! nn_eos = 1 : linear equation of state function of temperature only
	102	!! prd(t) = 0.0285 - rn_alpha * t
	103	!! nn_eos = 2 : linear equation of state function of temperature and
[3]	104	!! salinity
[1601]	105	!! prd(t,s) = rn_beta * s - rn_alpha * tn - 1.
[3]	106	!! Note that no boundary condition problem occurs in this routine
[2528]	107	!! as pts are defined over the whole domain.
[3]	108	!!
	109	!! ** Action : compute prd , the in situ density (no units)
	110	!!
[888]	111	!! References : Jackett and McDougall, J. Atmos. Ocean. Tech., 1994
	112	!!----------------------------------------------------------------------
[719]	113	!!
[2715]	114	REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pts ! 1 : potential temperature [Celcius]
	115	! ! 2 : salinity [psu]
	116	REAL(wp), DIMENSION(:,:,:) , INTENT( out) :: prd ! in situ density [-]
[4292]	117	REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pdep ! depth [m]
[2715]	118	!!
[1559]	119	INTEGER :: ji, jj, jk ! dummy loop indices
[2715]	120	REAL(wp) :: zt , zs , zh , zsr ! local scalars
	121	REAL(wp) :: zr1, zr2, zr3, zr4 ! - -
	122	REAL(wp) :: zrhop, ze, zbw, zb ! - -
	123	REAL(wp) :: zd , zc , zaw, za ! - -
	124	REAL(wp) :: zb1, za1, zkw, zk0 ! - -
[3294]	125	REAL(wp), POINTER, DIMENSION(:,:,:) :: zws
[3]	126	!!----------------------------------------------------------------------
	127
[3294]	128	!
	129	IF( nn_timing == 1 ) CALL timing_start('eos')
	130	!
	131	CALL wrk_alloc( jpi, jpj, jpk, zws )
	132	!
[1601]	133	SELECT CASE( nn_eos )
[888]	134	!
[1559]	135	CASE( 0 ) !== Jackett and McDougall (1994) formulation ==!
[3]	136	!CDIR NOVERRCHK
[2528]	137	zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) )
[3294]	138	!
[1559]	139	DO jk = 1, jpkm1
[3]	140	DO jj = 1, jpj
	141	DO ji = 1, jpi
[2528]	142	zt = pts (ji,jj,jk,jp_tem)
	143	zs = pts (ji,jj,jk,jp_sal)
[4292]	144	zh = pdep(ji,jj,jk) ! depth
[1559]	145	zsr= zws (ji,jj,jk) ! square root salinity
	146	!
[3]	147	! compute volumic mass pure water at atm pressure
[2528]	148	zr1= ( ( ( ( 6.536332e-9_wp zt - 1.120083e-6_wp )zt + 1.001685e-4_wp )*zt &
	149	& -9.095290e-3_wp )zt + 6.793952e-2_wp )zt + 999.842594_wp
[3]	150	! seawater volumic mass atm pressure
[2528]	151	zr2= ( ( ( 5.3875e-9_wpzt-8.2467e-7_wp ) zt+7.6438e-5_wp ) *zt &
	152	& -4.0899e-3_wp ) *zt+0.824493_wp
	153	zr3= ( -1.6546e-6_wpzt+1.0227e-4_wp ) zt-5.72466e-3_wp
	154	zr4= 4.8314e-4_wp
[1559]	155	!
[3]	156	! potential volumic mass (reference to the surface)
	157	zrhop= ( zr4zs + zr3zsr + zr2 ) *zs + zr1
[1559]	158	!
[3]	159	! add the compression terms
[2528]	160	ze = ( -3.508914e-8_wpzt-1.248266e-8_wp ) zt-2.595994e-6_wp
	161	zbw= ( 1.296821e-6_wpzt-5.782165e-9_wp ) zt+1.045941e-4_wp
[3]	162	zb = zbw + ze * zs
[1559]	163	!
[2528]	164	zd = -2.042967e-2_wp
	165	zc = (-7.267926e-5_wpzt+2.598241e-3_wp ) zt+0.1571896_wp
	166	zaw= ( ( 5.939910e-6_wpzt+2.512549e-3_wp ) zt-0.1028859_wp ) *zt - 4.721788_wp
[3]	167	za = ( zdzsr + zc ) zs + zaw
[1559]	168	!
[2528]	169	zb1= (-0.1909078_wpzt+7.390729_wp ) zt-55.87545_wp
	170	za1= ( ( 2.326469e-3_wpzt+1.553190_wp) zt-65.00517_wp ) *zt+1044.077_wp
	171	zkw= ( ( (-1.361629e-4_wpzt-1.852732e-2_wp ) zt-30.41638_wp ) zt + 2098.925_wp ) zt+190925.6_wp
[3]	172	zk0= ( zb1zsr + za1 )zs + zkw
[1559]	173	!
[3]	174	! masked in situ density anomaly
[2528]	175	prd(ji,jj,jk) = ( zrhop / ( 1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) ) ) &
[3625]	176	& - rau0 ) * r1_rau0 * tmask(ji,jj,jk)
[3]	177	END DO
	178	END DO
[1559]	179	END DO
[888]	180	!
[1559]	181	CASE( 1 ) !== Linear formulation function of temperature only ==!
	182	DO jk = 1, jpkm1
[2528]	183	prd(:,:,jk) = ( 0.0285_wp - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask(:,:,jk)
[1559]	184	END DO
[888]	185	!
[1559]	186	CASE( 2 ) !== Linear formulation function of temperature and salinity ==!
	187	DO jk = 1, jpkm1
[2528]	188	prd(:,:,jk) = ( rn_beta * pts(:,:,jk,jp_sal) - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask(:,:,jk)
[1559]	189	END DO
[888]	190	!
[3]	191	END SELECT
[888]	192	!
[1559]	193	IF(ln_ctl) CALL prt_ctl( tab3d_1=prd, clinfo1=' eos : ', ovlap=1, kdim=jpk )
[888]	194	!
[3294]	195	CALL wrk_dealloc( jpi, jpj, jpk, zws )
[2715]	196	!
[3294]	197	IF( nn_timing == 1 ) CALL timing_stop('eos')
	198	!
[3]	199	END SUBROUTINE eos_insitu
	200
	201
[4292]	202	SUBROUTINE eos_insitu_pot( pts, prd, prhop, pdep )
[3]	203	!!----------------------------------------------------------------------
	204	!! * ROUTINE eos_insitu_pot *
[3294]	205	!!
[3]	206	!! ** Purpose : Compute the in situ density (ratio rho/rau0) and the
	207	!! potential volumic mass (Kg/m3) from potential temperature and
[3294]	208	!! salinity fields using an equation of state defined through the
[1601]	209	!! namelist parameter nn_eos.
[3]	210	!!
	211	!! ** Method :
[1601]	212	!! nn_eos = 0 : Jackett and McDougall (1994) equation of state.
[3]	213	!! the in situ density is computed directly as a function of
	214	!! potential temperature relative to the surface (the opa t
	215	!! variable), salt and pressure (assuming no pressure variation
	216	!! along geopotential surfaces, i.e. the pressure p in decibars
	217	!! is approximated by the depth in meters.
	218	!! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0
	219	!! rhop(t,s) = rho(t,s,0)
	220	!! with pressure p decibars
	221	!! potential temperature t deg celsius
	222	!! salinity s psu
	223	!! reference volumic mass rau0 kg/m**3
	224	!! in situ volumic mass rho kg/m**3
	225	!! in situ density anomalie prd no units
	226	!!
	227	!! Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar,
	228	!! t = 40 deg celcius, s=40 psu
	229	!!
[1601]	230	!! nn_eos = 1 : linear equation of state function of temperature only
	231	!! prd(t) = ( rho(t) - rau0 ) / rau0 = 0.028 - rn_alpha * t
[3]	232	!! rhop(t,s) = rho(t,s)
	233	!!
[1601]	234	!! nn_eos = 2 : linear equation of state function of temperature and
[3]	235	!! salinity
[3294]	236	!! prd(t,s) = ( rho(t,s) - rau0 ) / rau0
[1601]	237	!! = rn_beta * s - rn_alpha * tn - 1.
[3]	238	!! rhop(t,s) = rho(t,s)
	239	!! Note that no boundary condition problem occurs in this routine
	240	!! as (tn,sn) or (ta,sa) are defined over the whole domain.
	241	!!
	242	!! ** Action : - prd , the in situ density (no units)
	243	!! - prhop, the potential volumic mass (Kg/m3)
	244	!!
[888]	245	!! References : Jackett and McDougall, J. Atmos. Ocean. Tech., 1994
	246	!! Brown and Campana, Mon. Weather Rev., 1978
[3]	247	!!----------------------------------------------------------------------
[2715]	248	!!
	249	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! 1 : potential temperature [Celcius]
	250	! ! 2 : salinity [psu]
	251	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT( out) :: prd ! in situ density [-]
[2528]	252	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT( out) :: prhop ! potential density (surface referenced)
[4292]	253	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pdep ! depth [m]
[2715]	254	!
[1559]	255	INTEGER :: ji, jj, jk ! dummy loop indices
[2715]	256	REAL(wp) :: zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw ! local scalars
[3625]	257	REAL(wp) :: zb, zd, zc, zaw, za, zb1, za1, zkw, zk0 ! - -
[3294]	258	REAL(wp), POINTER, DIMENSION(:,:,:) :: zws
[3]	259	!!----------------------------------------------------------------------
[3294]	260	!
	261	IF( nn_timing == 1 ) CALL timing_start('eos-p')
	262	!
	263	CALL wrk_alloc( jpi, jpj, jpk, zws )
	264	!
[1601]	265	SELECT CASE ( nn_eos )
[888]	266	!
[1559]	267	CASE( 0 ) !== Jackett and McDougall (1994) formulation ==!
[3]	268	!CDIR NOVERRCHK
[2528]	269	zws(:,:,:) = SQRT( ABS( pts(:,:,:,jp_sal) ) )
[3294]	270	!
[1559]	271	DO jk = 1, jpkm1
[3]	272	DO jj = 1, jpj
	273	DO ji = 1, jpi
[2528]	274	zt = pts (ji,jj,jk,jp_tem)
	275	zs = pts (ji,jj,jk,jp_sal)
[4292]	276	zh = pdep(ji,jj,jk) ! depth
[1559]	277	zsr= zws (ji,jj,jk) ! square root salinity
	278	!
[3]	279	! compute volumic mass pure water at atm pressure
[2528]	280	zr1= ( ( ( ( 6.536332e-9_wpzt-1.120083e-6_wp )zt+1.001685e-4_wp )*zt &
	281	& -9.095290e-3_wp )zt+6.793952e-2_wp )zt+999.842594_wp
[3]	282	! seawater volumic mass atm pressure
[2528]	283	zr2= ( ( ( 5.3875e-9_wpzt-8.2467e-7_wp ) zt+7.6438e-5_wp ) *zt &
	284	& -4.0899e-3_wp ) *zt+0.824493_wp
	285	zr3= ( -1.6546e-6_wpzt+1.0227e-4_wp ) zt-5.72466e-3_wp
	286	zr4= 4.8314e-4_wp
[1559]	287	!
[3]	288	! potential volumic mass (reference to the surface)
	289	zrhop= ( zr4zs + zr3zsr + zr2 ) *zs + zr1
[1559]	290	!
[3]	291	! save potential volumic mass
	292	prhop(ji,jj,jk) = zrhop * tmask(ji,jj,jk)
[1559]	293	!
[3]	294	! add the compression terms
[2528]	295	ze = ( -3.508914e-8_wpzt-1.248266e-8_wp ) zt-2.595994e-6_wp
	296	zbw= ( 1.296821e-6_wpzt-5.782165e-9_wp ) zt+1.045941e-4_wp
[3]	297	zb = zbw + ze * zs
[1559]	298	!
[2528]	299	zd = -2.042967e-2_wp
	300	zc = (-7.267926e-5_wpzt+2.598241e-3_wp ) zt+0.1571896_wp
	301	zaw= ( ( 5.939910e-6_wpzt+2.512549e-3_wp ) zt-0.1028859_wp ) *zt - 4.721788_wp
[3]	302	za = ( zdzsr + zc ) zs + zaw
[1559]	303	!
[2528]	304	zb1= ( -0.1909078_wp zt+7.390729_wp ) zt-55.87545_wp
	305	za1= ( ( 2.326469e-3_wpzt+1.553190_wp ) zt-65.00517_wp ) *zt + 1044.077_wp
	306	zkw= ( ( (-1.361629e-4_wpzt-1.852732e-2_wp ) zt-30.41638_wp ) zt + 2098.925_wp ) zt+190925.6_wp
[3]	307	zk0= ( zb1zsr + za1 )zs + zkw
[1559]	308	!
[3]	309	! masked in situ density anomaly
[2528]	310	prd(ji,jj,jk) = ( zrhop / ( 1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) ) ) &
[3625]	311	& - rau0 ) * r1_rau0 * tmask(ji,jj,jk)
[3]	312	END DO
	313	END DO
[1559]	314	END DO
[888]	315	!
[1559]	316	CASE( 1 ) !== Linear formulation = F( temperature ) ==!
	317	DO jk = 1, jpkm1
[2528]	318	prd (:,:,jk) = ( 0.0285_wp - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask(:,:,jk)
	319	prhop(:,:,jk) = ( 1.e0_wp + prd (:,:,jk) ) * rau0 * tmask(:,:,jk)
[1559]	320	END DO
[888]	321	!
[1559]	322	CASE( 2 ) !== Linear formulation = F( temperature , salinity ) ==!
	323	DO jk = 1, jpkm1
[2528]	324	prd (:,:,jk) = ( rn_beta * pts(:,:,jk,jp_sal) - rn_alpha * pts(:,:,jk,jp_tem) ) * tmask(:,:,jk)
	325	prhop(:,:,jk) = ( 1.e0_wp + prd (:,:,jk) ) * rau0 * tmask(:,:,jk)
[1559]	326	END DO
[888]	327	!
[3]	328	END SELECT
[888]	329	!
	330	IF(ln_ctl) CALL prt_ctl( tab3d_1=prd, clinfo1=' eos-p: ', tab3d_2=prhop, clinfo2=' pot : ', ovlap=1, kdim=jpk )
	331	!
[3294]	332	CALL wrk_dealloc( jpi, jpj, jpk, zws )
[2715]	333	!
[3294]	334	IF( nn_timing == 1 ) CALL timing_stop('eos-p')
	335	!
[888]	336	END SUBROUTINE eos_insitu_pot
[719]	337
	338
[2528]	339	SUBROUTINE eos_insitu_2d( pts, pdep, prd )
[3]	340	!!----------------------------------------------------------------------
	341	!! * ROUTINE eos_insitu_2d *
	342	!!
[3294]	343	!! ** Purpose : Compute the in situ density (ratio rho/rau0) from
[3]	344	!! potential temperature and salinity using an equation of state
[1601]	345	!! defined through the namelist parameter nn_eos. * 2D field case
[3]	346	!!
	347	!! ** Method :
[1601]	348	!! nn_eos = 0 : Jackett and McDougall (1994) equation of state.
[3]	349	!! the in situ density is computed directly as a function of
	350	!! potential temperature relative to the surface (the opa t
	351	!! variable), salt and pressure (assuming no pressure variation
	352	!! along geopotential surfaces, i.e. the pressure p in decibars
	353	!! is approximated by the depth in meters.
	354	!! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0
	355	!! with pressure p decibars
	356	!! potential temperature t deg celsius
	357	!! salinity s psu
	358	!! reference volumic mass rau0 kg/m**3
	359	!! in situ volumic mass rho kg/m**3
	360	!! in situ density anomalie prd no units
	361	!! Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar,
	362	!! t = 40 deg celcius, s=40 psu
[1601]	363	!! nn_eos = 1 : linear equation of state function of temperature only
	364	!! prd(t) = 0.0285 - rn_alpha * t
	365	!! nn_eos = 2 : linear equation of state function of temperature and
[3]	366	!! salinity
[1601]	367	!! prd(t,s) = rn_beta * s - rn_alpha * tn - 1.
[3]	368	!! Note that no boundary condition problem occurs in this routine
[2528]	369	!! as pts are defined over the whole domain.
[3]	370	!!
	371	!! ** Action : - prd , the in situ density (no units)
	372	!!
[888]	373	!! References : Jackett and McDougall, J. Atmos. Ocean. Tech., 1994
	374	!!----------------------------------------------------------------------
[2715]	375	!!
[2528]	376	REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in ) :: pts ! 1 : potential temperature [Celcius]
	377	! ! 2 : salinity [psu]
	378	REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pdep ! depth [m]
[3294]	379	REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: prd ! in situ density
[719]	380	!!
[1559]	381	INTEGER :: ji, jj ! dummy loop indices
	382	REAL(wp) :: zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw ! temporary scalars
	383	REAL(wp) :: zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zmask ! - -
[3294]	384	REAL(wp), POINTER, DIMENSION(:,:) :: zws
[3]	385	!!----------------------------------------------------------------------
[3294]	386	!
	387	IF( nn_timing == 1 ) CALL timing_start('eos2d')
	388	!
	389	CALL wrk_alloc( jpi, jpj, zws )
	390	!
[3]	391
[2715]	392	prd(:,:) = 0._wp
	393
[1601]	394	SELECT CASE( nn_eos )
[888]	395	!
[1559]	396	CASE( 0 ) !== Jackett and McDougall (1994) formulation ==!
[888]	397	!
[3]	398	!CDIR NOVERRCHK
[4666]	399	DO jj = 1, jpj
[3]	400	!CDIR NOVERRCHK
[4666]	401	DO ji = 1, jpi ! vector opt.
[2528]	402	zws(ji,jj) = SQRT( ABS( pts(ji,jj,jp_sal) ) )
[3]	403	END DO
	404	END DO
[4666]	405	DO jj = 1, jpj
	406	DO ji = 1, jpi ! vector opt.
	407	zmask = lmask(ji,jj) ! land/sea bottom mask = surf. mask
[2528]	408	zt = pts (ji,jj,jp_tem) ! interpolated T
	409	zs = pts (ji,jj,jp_sal) ! interpolated S
[1559]	410	zsr = zws (ji,jj) ! square root of interpolated S
	411	zh = pdep (ji,jj) ! depth at the partial step level
	412	!
[3]	413	! compute volumic mass pure water at atm pressure
[2528]	414	zr1 = ( ( ( ( 6.536332e-9_wpzt-1.120083e-6_wp )zt+1.001685e-4_wp )*zt &
	415	& -9.095290e-3_wp )zt+6.793952e-2_wp )zt+999.842594_wp
[3]	416	! seawater volumic mass atm pressure
[2528]	417	zr2 = ( ( ( 5.3875e-9_wpzt-8.2467e-7_wp )zt+7.6438e-5_wp ) *zt &
	418	& -4.0899e-3_wp ) *zt+0.824493_wp
	419	zr3 = ( -1.6546e-6_wpzt+1.0227e-4_wp ) zt-5.72466e-3_wp
	420	zr4 = 4.8314e-4_wp
[1559]	421	!
[3]	422	! potential volumic mass (reference to the surface)
	423	zrhop= ( zr4zs + zr3zsr + zr2 ) *zs + zr1
[1559]	424	!
[3]	425	! add the compression terms
[2528]	426	ze = ( -3.508914e-8_wpzt-1.248266e-8_wp ) zt-2.595994e-6_wp
	427	zbw= ( 1.296821e-6_wpzt-5.782165e-9_wp ) zt+1.045941e-4_wp
[3]	428	zb = zbw + ze * zs
[1559]	429	!
[2528]	430	zd = -2.042967e-2_wp
	431	zc = (-7.267926e-5_wpzt+2.598241e-3_wp ) zt+0.1571896_wp
	432	zaw= ( ( 5.939910e-6_wpzt+2.512549e-3_wp ) zt-0.1028859_wp ) *zt -4.721788_wp
[3]	433	za = ( zdzsr + zc ) zs + zaw
[1559]	434	!
[2528]	435	zb1= (-0.1909078_wp zt+7.390729_wp ) zt-55.87545_wp
	436	za1= ( ( 2.326469e-3_wpzt+1.553190_wp ) zt-65.00517_wp ) *zt+1044.077_wp
	437	zkw= ( ( (-1.361629e-4_wpzt-1.852732e-2_wp ) zt-30.41638_wp ) *zt &
	438	& +2098.925_wp ) *zt+190925.6_wp
[3]	439	zk0= ( zb1zsr + za1 )zs + zkw
[1559]	440	!
[3]	441	! masked in situ density anomaly
[2528]	442	prd(ji,jj) = ( zrhop / ( 1.0_wp - zh / ( zk0 - zh * ( za - zh * zb ) ) ) - rau0 ) / rau0 * zmask
[3]	443	END DO
[1559]	444	END DO
[888]	445	!
[1559]	446	CASE( 1 ) !== Linear formulation = F( temperature ) ==!
[4666]	447	DO jj = 1, jpj
	448	DO ji = 1, jpi ! vector opt.
	449	prd(ji,jj) = ( 0.0285_wp - rn_alpha * pts(ji,jj,jp_tem) ) * lmask(ji,jj)
[3]	450	END DO
[1559]	451	END DO
[888]	452	!
[1559]	453	CASE( 2 ) !== Linear formulation = F( temperature , salinity ) ==!
[4666]	454	DO jj = 1, jpj
	455	DO ji = 1, jpi ! vector opt.
	456	prd(ji,jj) = ( rn_beta * pts(ji,jj,jp_sal) - rn_alpha * pts(ji,jj,jp_tem) ) * lmask(ji,jj)
[3]	457	END DO
[1559]	458	END DO
[888]	459	!
[3]	460	END SELECT
	461
[888]	462	IF(ln_ctl) CALL prt_ctl( tab2d_1=prd, clinfo1=' eos2d: ' )
	463	!
[3294]	464	CALL wrk_dealloc( jpi, jpj, zws )
[2715]	465	!
[3294]	466	IF( nn_timing == 1 ) CALL timing_stop('eos2d')
	467	!
[3]	468	END SUBROUTINE eos_insitu_2d
	469
	470
[2528]	471	SUBROUTINE eos_bn2( pts, pn2 )
[3]	472	!!----------------------------------------------------------------------
	473	!! * ROUTINE eos_bn2 *
	474	!!
	475	!! ** Purpose : Compute the local Brunt-Vaisala frequency at the time-
	476	!! step of the input arguments
[3294]	477	!!
[3]	478	!! ** Method :
[1601]	479	!! * nn_eos = 0 : UNESCO sea water properties
[3]	480	!! The brunt-vaisala frequency is computed using the polynomial
	481	!! polynomial expression of McDougall (1987):
[15]	482	!! N^2 = grav * beta * ( alpha/beta*dk[ t ] - dk[ s ] )/e3w
[3]	483	!! If lk_zdfddm=T, the heat/salt buoyancy flux ratio Rrau is
	484	!! computed and used in zdfddm module :
	485	!! Rrau = alpha/beta * ( dk[ t ] / dk[ s ] )
[1601]	486	!! * nn_eos = 1 : linear equation of state (temperature only)
	487	!! N^2 = grav * rn_alpha * dk[ t ]/e3w
	488	!! * nn_eos = 2 : linear equation of state (temperature & salinity)
	489	!! N^2 = grav * (rn_alpha * dk[ t ] - rn_beta * dk[ s ] ) / e3w
[3]	490	!! The use of potential density to compute N^2 introduces e r r o r
[3294]	491	!! in the sign of N^2 at great depths. We recommand the use of
[1601]	492	!! nn_eos = 0, except for academical studies.
[3]	493	!! Macro-tasked on horizontal slab (jk-loop)
	494	!! N.B. N^2 is set to zero at the first level (JK=1) in inidtr
	495	!! and is never used at this level.
	496	!!
	497	!! ** Action : - pn2 : the brunt-vaisala frequency
	498	!!
[888]	499	!! References : McDougall, J. Phys. Oceanogr., 17, 1950-1964, 1987.
[3]	500	!!----------------------------------------------------------------------
[2528]	501	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! 1 : potential temperature [Celcius]
	502	! ! 2 : salinity [psu]
[2715]	503	REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: pn2 ! Brunt-Vaisala frequency [s-1]
[1559]	504	!!
[3]	505	INTEGER :: ji, jj, jk ! dummy loop indices
[3294]	506	REAL(wp) :: zgde3w, zt, zs, zh, zalbet, zbeta ! local scalars
[15]	507	#if defined key_zdfddm
[2715]	508	REAL(wp) :: zds ! local scalars
[15]	509	#endif
[3]	510	!!----------------------------------------------------------------------
	511
[3294]	512	!
	513	IF( nn_timing == 1 ) CALL timing_start('bn2')
	514	!
[3]	515	! pn2 : interior points only (2=< jk =< jpkm1 )
[3294]	516	! --------------------------
[888]	517	!
[1601]	518	SELECT CASE( nn_eos )
[1559]	519	!
	520	CASE( 0 ) !== Jackett and McDougall (1994) formulation ==!
	521	DO jk = 2, jpkm1
[3]	522	DO jj = 1, jpj
	523	DO ji = 1, jpi
[15]	524	zgde3w = grav / fse3w(ji,jj,jk)
[4666]	525	zt = 0.5_wp * ( pts(ji,jj,jk,jp_tem) + pts(ji,jj,jk-1,jp_tem) ) ! potential temperature at w-pt
	526	zs = 0.5_wp * ( pts(ji,jj,jk,jp_sal) + pts(ji,jj,jk-1,jp_sal) ) - 35.0_wp ! salinity anomaly (s-35) at w-pt
[2715]	527	zh = fsdepw(ji,jj,jk) ! depth in meters at w-point
[1559]	528	!
[2528]	529	zalbet = ( ( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt & ! ratio alpha/beta
	530	& - 0.203814e-03_wp ) * zt &
	531	& + 0.170907e-01_wp ) * zt &
	532	& + 0.665157e-01_wp &
	533	& + ( - 0.678662e-05_wp * zs &
	534	& - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs &
	535	& + ( ( - 0.302285e-13_wp * zh &
	536	& - 0.251520e-11_wp * zs &
	537	& + 0.512857e-12_wp * zt * zt ) * zh &
	538	& - 0.164759e-06_wp * zs &
	539	& +( 0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt &
	540	& + 0.380374e-04_wp ) * zh
[1559]	541	!
[2528]	542	zbeta = ( ( -0.415613e-09_wp * zt + 0.555579e-07_wp ) * zt & ! beta
	543	& - 0.301985e-05_wp ) * zt &
	544	& + 0.785567e-03_wp &
	545	& + ( 0.515032e-08_wp * zs &
	546	& + 0.788212e-08_wp * zt - 0.356603e-06_wp ) * zs &
	547	& + ( ( 0.121551e-17_wp * zh &
	548	& - 0.602281e-15_wp * zs &
	549	& - 0.175379e-14_wp * zt + 0.176621e-12_wp ) * zh &
	550	& + 0.408195e-10_wp * zs &
	551	& + ( - 0.213127e-11_wp * zt + 0.192867e-09_wp ) * zt &
	552	& - 0.121555e-07_wp ) * zh
[1559]	553	!
[4666]	554	pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) & ! N^2
[2528]	555	& * ( zalbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) &
	556	& - ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) )
[3]	557	#if defined key_zdfddm
	558	! !!bug **** caution a traiter zds=dk[S]= 0 !!!!
[2528]	559	zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) ) ! Rrau = (alpha / beta) (dk[t] / dk[s])
	560	IF ( ABS( zds) <= 1.e-20_wp ) zds = 1.e-20_wp
	561	rrau(ji,jj,jk) = zalbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds
[3]	562	#endif
	563	END DO
	564	END DO
[1559]	565	END DO
[888]	566	!
[1559]	567	CASE( 1 ) !== Linear formulation = F( temperature ) ==!
	568	DO jk = 2, jpkm1
[4666]	569	pn2(:,:,jk) = grav * rn_alpha * ( pts(:,:,jk-1,jp_tem) - pts(:,:,jk,jp_tem) ) &
	570	& / fse3w(:,:,jk) * tmask(:,:,jk) * tmask(:,:,jk-1)
[1559]	571	END DO
[888]	572	!
[1559]	573	CASE( 2 ) !== Linear formulation = F( temperature , salinity ) ==!
	574	DO jk = 2, jpkm1
[2528]	575	pn2(:,:,jk) = grav * ( rn_alpha * ( pts(:,:,jk-1,jp_tem) - pts(:,:,jk,jp_tem) ) &
	576	& - rn_beta * ( pts(:,:,jk-1,jp_sal) - pts(:,:,jk,jp_sal) ) ) &
[4666]	577	& / fse3w(:,:,jk) * tmask(:,:,jk) * tmask(:,:,jk-1)
[3294]	578	END DO
[3]	579	#if defined key_zdfddm
[1559]	580	DO jk = 2, jpkm1 ! Rrau = (alpha / beta) (dk[t] / dk[s])
[3]	581	DO jj = 1, jpj
	582	DO ji = 1, jpi
[3294]	583	zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )
[2528]	584	IF ( ABS( zds ) <= 1.e-20_wp ) zds = 1.e-20_wp
	585	rrau(ji,jj,jk) = ralpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds
[3]	586	END DO
	587	END DO
[1559]	588	END DO
[3]	589	#endif
	590	END SELECT
	591
[1559]	592	IF(ln_ctl) CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2 : ', ovlap=1, kdim=jpk )
[49]	593	#if defined key_zdfddm
[1559]	594	IF(ln_ctl) CALL prt_ctl( tab3d_1=rrau, clinfo1=' rrau : ', ovlap=1, kdim=jpk )
[49]	595	#endif
[888]	596	!
[3294]	597	IF( nn_timing == 1 ) CALL timing_stop('bn2')
	598	!
[3]	599	END SUBROUTINE eos_bn2
	600
	601
[3294]	602	SUBROUTINE eos_alpbet( pts, palpbet, beta0 )
[2528]	603	!!----------------------------------------------------------------------
[3294]	604	!! * ROUTINE eos_alpbet *
[2528]	605	!!
[3294]	606	!! ** Purpose : Calculates the in situ thermal/haline expansion ratio at T-points
[2528]	607	!!
[3294]	608	!! ** Method : calculates alpha / beta ratio at T-points
[2528]	609	!! * nn_eos = 0 : UNESCO sea water properties
[3294]	610	!! The alpha/beta ratio is returned as 3-D array palpbet using the polynomial
	611	!! polynomial expression of McDougall (1987).
	612	!! Scalar beta0 is returned = 1.
[2528]	613	!! * nn_eos = 1 : linear equation of state (temperature only)
[3294]	614	!! The ratio is undefined, so we return alpha as palpbet
	615	!! Scalar beta0 is returned = 0.
[2528]	616	!! * nn_eos = 2 : linear equation of state (temperature & salinity)
[3294]	617	!! The alpha/beta ratio is returned as ralpbet
	618	!! Scalar beta0 is returned = 1.
[2528]	619	!!
[3294]	620	!! ** Action : - palpbet : thermal/haline expansion ratio at T-points
	621	!! : beta0 : 1. or 0.
[2528]	622	!!----------------------------------------------------------------------
[3294]	623	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature & salinity
	624	REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: palpbet ! thermal/haline expansion ratio
	625	REAL(wp), INTENT( out) :: beta0 ! set = 1 except with case 1 eos, rho=rho(T)
	626	!!
[2528]	627	INTEGER :: ji, jj, jk ! dummy loop indices
[3294]	628	REAL(wp) :: zt, zs, zh ! local scalars
[2528]	629	!!----------------------------------------------------------------------
	630	!
[3294]	631	IF( nn_timing == 1 ) CALL timing_start('eos_alpbet')
	632	!
[2528]	633	SELECT CASE ( nn_eos )
	634	!
	635	CASE ( 0 ) ! Jackett and McDougall (1994) formulation
	636	DO jk = 1, jpk
	637	DO jj = 1, jpj
	638	DO ji = 1, jpi
	639	zt = pts(ji,jj,jk,jp_tem) ! potential temperature
	640	zs = pts(ji,jj,jk,jp_sal) - 35._wp ! salinity anomaly (s-35)
[3294]	641	zh = fsdept(ji,jj,jk) ! depth in meters
[2528]	642	!
[3294]	643	palpbet(ji,jj,jk) = &
	644	& ( ( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt &
	645	& - 0.203814e-03_wp ) * zt &
	646	& + 0.170907e-01_wp ) * zt &
	647	& + 0.665157e-01_wp &
	648	& + ( - 0.678662e-05_wp * zs &
	649	& - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs &
	650	& + ( ( - 0.302285e-13_wp * zh &
	651	& - 0.251520e-11_wp * zs &
	652	& + 0.512857e-12_wp * zt * zt ) * zh &
	653	& - 0.164759e-06_wp * zs &
	654	& +( 0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt &
	655	& + 0.380374e-04_wp ) * zh
[2528]	656	END DO
	657	END DO
	658	END DO
[3294]	659	beta0 = 1._wp
[2528]	660	!
[3294]	661	CASE ( 1 ) !== Linear formulation = F( temperature ) ==!
	662	palpbet(:,:,:) = rn_alpha
	663	beta0 = 0._wp
[2528]	664	!
[3294]	665	CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==!
	666	palpbet(:,:,:) = ralpbet
	667	beta0 = 1._wp
[2528]	668	!
	669	CASE DEFAULT
	670	IF(lwp) WRITE(numout,cform_err)
	671	IF(lwp) WRITE(numout,*) ' bad flag value for nn_eos = ', nn_eos
	672	nstop = nstop + 1
	673	!
	674	END SELECT
	675	!
[3294]	676	IF( nn_timing == 1 ) CALL timing_stop('eos_alpbet')
	677	!
[2528]	678	END SUBROUTINE eos_alpbet
	679
	680
[4162]	681	FUNCTION tfreez( psal, pdep ) RESULT( ptf )
[3]	682	!!----------------------------------------------------------------------
	683	!! * ROUTINE eos_init *
	684	!!
[888]	685	!! ** Purpose : Compute the sea surface freezing temperature [Celcius]
[3]	686	!!
[888]	687	!! ** Method : UNESCO freezing point at the surface (pressure = 0???)
	688	!! freezing point [Celcius]=(-.0575+1.710523e-3sqrt(abs(s))-2.154996e-4s)s-7.53e-4p
	689	!! checkvalue: tf= -2.588567 Celsius for s=40.0psu, p=500. decibars
[3]	690	!!
[888]	691	!! Reference : UNESCO tech. papers in the marine science no. 28. 1978
[3]	692	!!----------------------------------------------------------------------
[888]	693	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: psal ! salinity [psu]
[4162]	694	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ), OPTIONAL :: pdep ! depth [decibars]
[2715]	695	! Leave result array automatic rather than making explicitly allocated
[888]	696	REAL(wp), DIMENSION(jpi,jpj) :: ptf ! freezing temperature [Celcius]
[3]	697	!!----------------------------------------------------------------------
[1559]	698	!
[2528]	699	ptf(:,:) = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal(:,:) ) &
	700	& - 2.154996e-4_wp * psal(:,:) ) * psal(:,:)
[4162]	701	IF ( PRESENT( pdep ) ) THEN
	702	ptf(:,:) = ptf(:,:) - 7.53e-4_wp * pdep(:,:)
	703	ENDIF
[1559]	704	!
[888]	705	END FUNCTION tfreez
	706
[4666]	707	FUNCTION tfreez1D( psal, pdep ) RESULT( ptf )
	708	!!----------------------------------------------------------------------
	709	!! * ROUTINE eos_init *
	710	!!
	711	!! ** Purpose : Compute the sea surface freezing temperature [Celcius]
	712	!!
	713	!! ** Method : UNESCO freezing point at the surface (pressure = 0???)
	714	!! freezing point [Celcius]=(-.0575+1.710523e-3sqrt(abs(s))-2.154996e-4s)s-7.53e-4p
	715	!! checkvalue: tf= -2.588567 Celsius for s=40.0psu, p=500. decibars
	716	!!
	717	!! Reference : UNESCO tech. papers in the marine science no. 28. 1978
	718	!!----------------------------------------------------------------------
	719	REAL(wp), INTENT(in ) :: psal ! salinity [psu]
	720	REAL(wp), INTENT(in ), OPTIONAL :: pdep ! pressure [dBar]
	721	! Leave result array automatic rather than making explicitly allocated
	722	REAL(wp) :: ptf ! freezing temperature [Celcius]
	723	!!----------------------------------------------------------------------
	724	!
	725	ptf = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal ) &
	726	& - 2.154996e-4_wp * psal ) * psal
	727	IF ( PRESENT( pdep ) ) THEN
	728	ptf = ptf - 7.53e-4_wp * pdep
	729	ENDIF
	730	!
	731	END FUNCTION tfreez1D
[888]	732
[4666]	733
	734
[888]	735	SUBROUTINE eos_init
[719]	736	!!----------------------------------------------------------------------
[888]	737	!! * ROUTINE eos_init *
	738	!!
	739	!! ** Purpose : initializations for the equation of state
	740	!!
	741	!! ** Method : Read the namelist nameos and control the parameters
	742	!!----------------------------------------------------------------------
[1601]	743	NAMELIST/nameos/ nn_eos, rn_alpha, rn_beta
[1559]	744	!!----------------------------------------------------------------------
[4147]	745	INTEGER :: ios
[1559]	746	!
[4147]	747	REWIND( numnam_ref ) ! Namelist nameos in reference namelist : equation of state
	748	READ ( numnam_ref, nameos, IOSTAT = ios, ERR = 901 )
	749	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nameos in reference namelist', lwp )
	750
	751	REWIND( numnam_cfg ) ! Namelist nameos in configuration namelist : equation of state
	752	READ ( numnam_cfg, nameos, IOSTAT = ios, ERR = 902 )
	753	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nameos in configuration namelist', lwp )
[4624]	754	IF(lwm) WRITE( numond, nameos )
[1559]	755	!
	756	IF(lwp) THEN ! Control print
[3]	757	WRITE(numout,*)
	758	WRITE(numout,*) 'eos_init : equation of state'
	759	WRITE(numout,*) '~~~~~~~~'
	760	WRITE(numout,*) ' Namelist nameos : set eos parameters'
[1601]	761	WRITE(numout,*) ' flag for eq. of state and N^2 nn_eos = ', nn_eos
	762	WRITE(numout,*) ' thermal exp. coef. (linear) rn_alpha = ', rn_alpha
	763	WRITE(numout,*) ' saline exp. coef. (linear) rn_beta = ', rn_beta
[3]	764	ENDIF
[1559]	765	!
[1601]	766	SELECT CASE( nn_eos ) ! check option
[1559]	767	!
[1601]	768	CASE( 0 ) !== Jackett and McDougall (1994) formulation ==!
[888]	769	IF(lwp) WRITE(numout,*)
[3]	770	IF(lwp) WRITE(numout,*) ' use of Jackett & McDougall (1994) equation of state and'
	771	IF(lwp) WRITE(numout,*) ' McDougall (1987) Brunt-Vaisala frequency'
[888]	772	!
[1601]	773	CASE( 1 ) !== Linear formulation = F( temperature ) ==!
[888]	774	IF(lwp) WRITE(numout,*)
[1601]	775	IF(lwp) WRITE(numout,) ' use of linear eos rho(T) = rau0 ( 1.0285 - rn_alpha * T )'
[474]	776	IF( lk_zdfddm ) CALL ctl_stop( ' double diffusive mixing parameterization requires', &
	777	& ' that T and S are used as state variables' )
[888]	778	!
[1601]	779	CASE( 2 ) !== Linear formulation = F( temperature , salinity ) ==!
	780	ralpbet = rn_alpha / rn_beta
[888]	781	IF(lwp) WRITE(numout,*)
[1601]	782	IF(lwp) WRITE(numout,) ' use of linear eos rho(T,S) = rau0 ( rn_beta * S - rn_alpha * T )'
[888]	783	!
[1601]	784	CASE DEFAULT !== ERROR in nn_eos ==!
	785	WRITE(ctmp1,*) ' bad flag value for nn_eos = ', nn_eos
[474]	786	CALL ctl_stop( ctmp1 )
[1559]	787	!
[3]	788	END SELECT
[1559]	789	!
[3]	790	END SUBROUTINE eos_init
	791
	792	!!======================================================================
[3294]	793	END MODULE eosbn2

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

Download in other formats: