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/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90 @ 3211

Last change on this file since 3211 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.

Property svn:keywords set to Id

File size: 44.4 KB

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

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

Download in other formats: