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

Context Navigation

source: trunk/NEMO/OPA_SRC/eosbn2.F90 @ 719

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

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

Download in other formats: