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

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90 @ 3389

Last change on this file since 3389 was 3294, checked in by rblod, 12 years ago
Merge of 3.4beta into the trunk
Property svn:keywords set to `Id`
File size: 38.4 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
[888]	57	PUBLIC tfreez ! called by sbcice_... modules
[3]	58
[2528]	59	! !!* Namelist (nameos) *
	60	INTEGER , PUBLIC :: nn_eos = 0 !: = 0/1/2 type of eq. of state and Brunt-Vaisala frequ.
	61	REAL(wp), PUBLIC :: rn_alpha = 2.0e-4_wp !: thermal expension coeff. (linear equation of state)
	62	REAL(wp), PUBLIC :: rn_beta = 7.7e-4_wp !: saline expension coeff. (linear equation of state)
[1601]	63
[2528]	64	REAL(wp), PUBLIC :: ralpbet !: alpha / beta ratio
[3294]	65
[3]	66	!! * Substitutions
	67	# include "domzgr_substitute.h90"
	68	# include "vectopt_loop_substitute.h90"
	69	!!----------------------------------------------------------------------
[2528]	70	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[888]	71	!! $Id$
[2528]	72	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	73	!!----------------------------------------------------------------------
	74	CONTAINS
	75
[2528]	76	SUBROUTINE eos_insitu( pts, prd )
[3]	77	!!----------------------------------------------------------------------
	78	!! * ROUTINE eos_insitu *
[3294]	79	!!
	80	!! ** Purpose : Compute the in situ density (ratio rho/rau0) from
[3]	81	!! potential temperature and salinity using an equation of state
[1601]	82	!! defined through the namelist parameter nn_eos.
[3]	83	!!
	84	!! ** Method : 3 cases:
[1601]	85	!! nn_eos = 0 : Jackett and McDougall (1994) equation of state.
[3]	86	!! the in situ density is computed directly as a function of
	87	!! potential temperature relative to the surface (the opa t
	88	!! variable), salt and pressure (assuming no pressure variation
	89	!! along geopotential surfaces, i.e. the pressure p in decibars
	90	!! is approximated by the depth in meters.
	91	!! prd(t,s,p) = ( rho(t,s,p) - rau0 ) / rau0
	92	!! with pressure p decibars
	93	!! potential temperature t deg celsius
	94	!! salinity s psu
	95	!! reference volumic mass rau0 kg/m**3
	96	!! in situ volumic mass rho kg/m**3
	97	!! in situ density anomalie prd no units
	98	!! Check value: rho = 1060.93298 kg/m**3 for p=10000 dbar,
	99	!! t = 40 deg celcius, s=40 psu
[1601]	100	!! nn_eos = 1 : linear equation of state function of temperature only
	101	!! prd(t) = 0.0285 - rn_alpha * t
	102	!! nn_eos = 2 : linear equation of state function of temperature and
[3]	103	!! salinity
[1601]	104	!! prd(t,s) = rn_beta * s - rn_alpha * tn - 1.
[3]	105	!! Note that no boundary condition problem occurs in this routine
[2528]	106	!! as pts are defined over the whole domain.
[3]	107	!!
	108	!! ** Action : compute prd , the in situ density (no units)
	109	!!
[888]	110	!! References : Jackett and McDougall, J. Atmos. Ocean. Tech., 1994
	111	!!----------------------------------------------------------------------
[719]	112	!!
[2715]	113	REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pts ! 1 : potential temperature [Celcius]
	114	! ! 2 : salinity [psu]
	115	REAL(wp), DIMENSION(:,:,:) , INTENT( out) :: prd ! in situ density [-]
	116	!!
[1559]	117	INTEGER :: ji, jj, jk ! dummy loop indices
[2715]	118	REAL(wp) :: zt , zs , zh , zsr ! local scalars
	119	REAL(wp) :: zr1, zr2, zr3, zr4 ! - -
	120	REAL(wp) :: zrhop, ze, zbw, zb ! - -
	121	REAL(wp) :: zd , zc , zaw, za ! - -
	122	REAL(wp) :: zb1, za1, zkw, zk0 ! - -
	123	REAL(wp) :: zrau0r ! - -
[3294]	124	REAL(wp), POINTER, DIMENSION(:,:,:) :: zws
[3]	125	!!----------------------------------------------------------------------
	126
[3294]	127	!
	128	IF( nn_timing == 1 ) CALL timing_start('eos')
	129	!
	130	CALL wrk_alloc( jpi, jpj, jpk, zws )
	131	!
[1601]	132	SELECT CASE( nn_eos )
[888]	133	!
[1559]	134	CASE( 0 ) !== Jackett and McDougall (1994) formulation ==!
	135	zrau0r = 1.e0 / rau0
[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)
[1559]	144	zh = fsdept(ji,jj,jk) ! depth
	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 ) ) ) &
[1559]	176	& - rau0 ) * zrau0r * 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
[2528]	202	SUBROUTINE eos_insitu_pot( pts, prd, prhop )
[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)
[2715]	253	!
[1559]	254	INTEGER :: ji, jj, jk ! dummy loop indices
[2715]	255	REAL(wp) :: zt, zs, zh, zsr, zr1, zr2, zr3, zr4, zrhop, ze, zbw ! local scalars
	256	REAL(wp) :: zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zrau0r ! - -
[3294]	257	REAL(wp), POINTER, DIMENSION(:,:,:) :: zws
[3]	258	!!----------------------------------------------------------------------
[3294]	259	!
	260	IF( nn_timing == 1 ) CALL timing_start('eos-p')
	261	!
	262	CALL wrk_alloc( jpi, jpj, jpk, zws )
	263	!
[1601]	264	SELECT CASE ( nn_eos )
[888]	265	!
[1559]	266	CASE( 0 ) !== Jackett and McDougall (1994) formulation ==!
	267	zrau0r = 1.e0 / rau0
[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)
[1559]	276	zh = fsdept(ji,jj,jk) ! depth
	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 ) ) ) &
[1559]	311	& - rau0 ) * zrau0r * 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
	399	DO jj = 1, jpjm1
	400	!CDIR NOVERRCHK
	401	DO ji = 1, fs_jpim1 ! vector opt.
[2528]	402	zws(ji,jj) = SQRT( ABS( pts(ji,jj,jp_sal) ) )
[3]	403	END DO
	404	END DO
[1559]	405	DO jj = 1, jpjm1
[3]	406	DO ji = 1, fs_jpim1 ! vector opt.
[1559]	407	zmask = tmask(ji,jj,1) ! 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 ) ==!
	447	DO jj = 1, jpjm1
[3]	448	DO ji = 1, fs_jpim1 ! vector opt.
[2528]	449	prd(ji,jj) = ( 0.0285_wp - rn_alpha * pts(ji,jj,jp_tem) ) * tmask(ji,jj,1)
[3]	450	END DO
[1559]	451	END DO
[888]	452	!
[1559]	453	CASE( 2 ) !== Linear formulation = F( temperature , salinity ) ==!
	454	DO jj = 1, jpjm1
[3]	455	DO ji = 1, fs_jpim1 ! vector opt.
[3294]	456	prd(ji,jj) = ( rn_beta * pts(ji,jj,jp_sal) - rn_alpha * pts(ji,jj,jp_tem) ) * tmask(ji,jj,1)
[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)
[2715]	525	zt = 0.5 * ( pts(ji,jj,jk,jp_tem) + pts(ji,jj,jk-1,jp_tem) ) ! potential temperature at w-pt
	526	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
	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	!
[3294]	554	pn2(ji,jj,jk) = zgde3w * zbeta * tmask(ji,jj,jk) & ! 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
[2528]	569	pn2(:,:,jk) = grav * rn_alpha * ( pts(:,:,jk-1,jp_tem) - pts(:,:,jk,jp_tem) ) / fse3w(:,:,jk) * tmask(:,:,jk)
[1559]	570	END DO
[888]	571	!
[1559]	572	CASE( 2 ) !== Linear formulation = F( temperature , salinity ) ==!
	573	DO jk = 2, jpkm1
[2528]	574	pn2(:,:,jk) = grav * ( rn_alpha * ( pts(:,:,jk-1,jp_tem) - pts(:,:,jk,jp_tem) ) &
	575	& - rn_beta * ( pts(:,:,jk-1,jp_sal) - pts(:,:,jk,jp_sal) ) ) &
[1559]	576	& / fse3w(:,:,jk) * tmask(:,:,jk)
[3294]	577	END DO
[3]	578	#if defined key_zdfddm
[1559]	579	DO jk = 2, jpkm1 ! Rrau = (alpha / beta) (dk[t] / dk[s])
[3]	580	DO jj = 1, jpj
	581	DO ji = 1, jpi
[3294]	582	zds = ( pts(ji,jj,jk-1,jp_sal) - pts(ji,jj,jk,jp_sal) )
[2528]	583	IF ( ABS( zds ) <= 1.e-20_wp ) zds = 1.e-20_wp
	584	rrau(ji,jj,jk) = ralpbet * ( pts(ji,jj,jk-1,jp_tem) - pts(ji,jj,jk,jp_tem) ) / zds
[3]	585	END DO
	586	END DO
[1559]	587	END DO
[3]	588	#endif
	589	END SELECT
	590
[1559]	591	IF(ln_ctl) CALL prt_ctl( tab3d_1=pn2, clinfo1=' bn2 : ', ovlap=1, kdim=jpk )
[49]	592	#if defined key_zdfddm
[1559]	593	IF(ln_ctl) CALL prt_ctl( tab3d_1=rrau, clinfo1=' rrau : ', ovlap=1, kdim=jpk )
[49]	594	#endif
[888]	595	!
[3294]	596	IF( nn_timing == 1 ) CALL timing_stop('bn2')
	597	!
[3]	598	END SUBROUTINE eos_bn2
	599
	600
[3294]	601	SUBROUTINE eos_alpbet( pts, palpbet, beta0 )
[2528]	602	!!----------------------------------------------------------------------
[3294]	603	!! * ROUTINE eos_alpbet *
[2528]	604	!!
[3294]	605	!! ** Purpose : Calculates the in situ thermal/haline expansion ratio at T-points
[2528]	606	!!
[3294]	607	!! ** Method : calculates alpha / beta ratio at T-points
[2528]	608	!! * nn_eos = 0 : UNESCO sea water properties
[3294]	609	!! The alpha/beta ratio is returned as 3-D array palpbet using the polynomial
	610	!! polynomial expression of McDougall (1987).
	611	!! Scalar beta0 is returned = 1.
[2528]	612	!! * nn_eos = 1 : linear equation of state (temperature only)
[3294]	613	!! The ratio is undefined, so we return alpha as palpbet
	614	!! Scalar beta0 is returned = 0.
[2528]	615	!! * nn_eos = 2 : linear equation of state (temperature & salinity)
[3294]	616	!! The alpha/beta ratio is returned as ralpbet
	617	!! Scalar beta0 is returned = 1.
[2528]	618	!!
[3294]	619	!! ** Action : - palpbet : thermal/haline expansion ratio at T-points
	620	!! : beta0 : 1. or 0.
[2528]	621	!!----------------------------------------------------------------------
[3294]	622	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(in ) :: pts ! pot. temperature & salinity
	623	REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT( out) :: palpbet ! thermal/haline expansion ratio
	624	REAL(wp), INTENT( out) :: beta0 ! set = 1 except with case 1 eos, rho=rho(T)
	625	!!
[2528]	626	INTEGER :: ji, jj, jk ! dummy loop indices
[3294]	627	REAL(wp) :: zt, zs, zh ! local scalars
[2528]	628	!!----------------------------------------------------------------------
	629	!
[3294]	630	IF( nn_timing == 1 ) CALL timing_start('eos_alpbet')
	631	!
[2528]	632	SELECT CASE ( nn_eos )
	633	!
	634	CASE ( 0 ) ! Jackett and McDougall (1994) formulation
	635	DO jk = 1, jpk
	636	DO jj = 1, jpj
	637	DO ji = 1, jpi
	638	zt = pts(ji,jj,jk,jp_tem) ! potential temperature
	639	zs = pts(ji,jj,jk,jp_sal) - 35._wp ! salinity anomaly (s-35)
[3294]	640	zh = fsdept(ji,jj,jk) ! depth in meters
[2528]	641	!
[3294]	642	palpbet(ji,jj,jk) = &
	643	& ( ( ( - 0.255019e-07_wp * zt + 0.298357e-05_wp ) * zt &
	644	& - 0.203814e-03_wp ) * zt &
	645	& + 0.170907e-01_wp ) * zt &
	646	& + 0.665157e-01_wp &
	647	& + ( - 0.678662e-05_wp * zs &
	648	& - 0.846960e-04_wp * zt + 0.378110e-02_wp ) * zs &
	649	& + ( ( - 0.302285e-13_wp * zh &
	650	& - 0.251520e-11_wp * zs &
	651	& + 0.512857e-12_wp * zt * zt ) * zh &
	652	& - 0.164759e-06_wp * zs &
	653	& +( 0.791325e-08_wp * zt - 0.933746e-06_wp ) * zt &
	654	& + 0.380374e-04_wp ) * zh
[2528]	655	END DO
	656	END DO
	657	END DO
[3294]	658	beta0 = 1._wp
[2528]	659	!
[3294]	660	CASE ( 1 ) !== Linear formulation = F( temperature ) ==!
	661	palpbet(:,:,:) = rn_alpha
	662	beta0 = 0._wp
[2528]	663	!
[3294]	664	CASE ( 2 ) !== Linear formulation = F( temperature , salinity ) ==!
	665	palpbet(:,:,:) = ralpbet
	666	beta0 = 1._wp
[2528]	667	!
	668	CASE DEFAULT
	669	IF(lwp) WRITE(numout,cform_err)
	670	IF(lwp) WRITE(numout,*) ' bad flag value for nn_eos = ', nn_eos
	671	nstop = nstop + 1
	672	!
	673	END SELECT
	674	!
[3294]	675	IF( nn_timing == 1 ) CALL timing_stop('eos_alpbet')
	676	!
[2528]	677	END SUBROUTINE eos_alpbet
	678
	679
[888]	680	FUNCTION tfreez( psal ) RESULT( ptf )
[3]	681	!!----------------------------------------------------------------------
	682	!! * ROUTINE eos_init *
	683	!!
[888]	684	!! ** Purpose : Compute the sea surface freezing temperature [Celcius]
[3]	685	!!
[888]	686	!! ** Method : UNESCO freezing point at the surface (pressure = 0???)
	687	!! freezing point [Celcius]=(-.0575+1.710523e-3sqrt(abs(s))-2.154996e-4s)s-7.53e-4p
	688	!! checkvalue: tf= -2.588567 Celsius for s=40.0psu, p=500. decibars
[3]	689	!!
[888]	690	!! Reference : UNESCO tech. papers in the marine science no. 28. 1978
[3]	691	!!----------------------------------------------------------------------
[888]	692	REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: psal ! salinity [psu]
[2715]	693	! Leave result array automatic rather than making explicitly allocated
[888]	694	REAL(wp), DIMENSION(jpi,jpj) :: ptf ! freezing temperature [Celcius]
[3]	695	!!----------------------------------------------------------------------
[1559]	696	!
[2528]	697	ptf(:,:) = ( - 0.0575_wp + 1.710523e-3_wp * SQRT( psal(:,:) ) &
	698	& - 2.154996e-4_wp * psal(:,:) ) * psal(:,:)
[1559]	699	!
[888]	700	END FUNCTION tfreez
	701
	702
	703	SUBROUTINE eos_init
[719]	704	!!----------------------------------------------------------------------
[888]	705	!! * ROUTINE eos_init *
	706	!!
	707	!! ** Purpose : initializations for the equation of state
	708	!!
	709	!! ** Method : Read the namelist nameos and control the parameters
	710	!!----------------------------------------------------------------------
[1601]	711	NAMELIST/nameos/ nn_eos, rn_alpha, rn_beta
[1559]	712	!!----------------------------------------------------------------------
	713	!
[888]	714	REWIND( numnam ) ! Read Namelist nameos : equation of state
[3]	715	READ ( numnam, nameos )
[1559]	716	!
	717	IF(lwp) THEN ! Control print
[3]	718	WRITE(numout,*)
	719	WRITE(numout,*) 'eos_init : equation of state'
	720	WRITE(numout,*) '~~~~~~~~'
	721	WRITE(numout,*) ' Namelist nameos : set eos parameters'
[1601]	722	WRITE(numout,*) ' flag for eq. of state and N^2 nn_eos = ', nn_eos
	723	WRITE(numout,*) ' thermal exp. coef. (linear) rn_alpha = ', rn_alpha
	724	WRITE(numout,*) ' saline exp. coef. (linear) rn_beta = ', rn_beta
[3]	725	ENDIF
[1559]	726	!
[1601]	727	SELECT CASE( nn_eos ) ! check option
[1559]	728	!
[1601]	729	CASE( 0 ) !== Jackett and McDougall (1994) formulation ==!
[888]	730	IF(lwp) WRITE(numout,*)
[3]	731	IF(lwp) WRITE(numout,*) ' use of Jackett & McDougall (1994) equation of state and'
	732	IF(lwp) WRITE(numout,*) ' McDougall (1987) Brunt-Vaisala frequency'
[888]	733	!
[1601]	734	CASE( 1 ) !== Linear formulation = F( temperature ) ==!
[888]	735	IF(lwp) WRITE(numout,*)
[1601]	736	IF(lwp) WRITE(numout,) ' use of linear eos rho(T) = rau0 ( 1.0285 - rn_alpha * T )'
[474]	737	IF( lk_zdfddm ) CALL ctl_stop( ' double diffusive mixing parameterization requires', &
	738	& ' that T and S are used as state variables' )
[888]	739	!
[1601]	740	CASE( 2 ) !== Linear formulation = F( temperature , salinity ) ==!
	741	ralpbet = rn_alpha / rn_beta
[888]	742	IF(lwp) WRITE(numout,*)
[1601]	743	IF(lwp) WRITE(numout,) ' use of linear eos rho(T,S) = rau0 ( rn_beta * S - rn_alpha * T )'
[888]	744	!
[1601]	745	CASE DEFAULT !== ERROR in nn_eos ==!
	746	WRITE(ctmp1,*) ' bad flag value for nn_eos = ', nn_eos
[474]	747	CALL ctl_stop( ctmp1 )
[1559]	748	!
[3]	749	END SELECT
[1559]	750	!
[3]	751	END SUBROUTINE eos_init
	752
	753	!!======================================================================
[3294]	754	END MODULE eosbn2

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

Download in other formats: