Home

Context Navigation

source: palm/trunk/SOURCE/eqn_state_seawater.f90 @ 755

Last change on this file since 755 was 484, checked in by raasch, 15 years ago
typo in file headers removed
Property svn:keywords set to `Id`
File size: 7.8 KB

Rev	Line
[96]	1	MODULE eqn_state_seawater_mod
	2
	3	!------------------------------------------------------------------------------!
[484]	4	! Current revisions:
[96]	5	! -----------------
[392]	6	!
[96]	7	!
	8	! Former revisions:
	9	! -----------------
[97]	10	! $Id: eqn_state_seawater.f90 484 2010-02-05 07:36:54Z witha $
[96]	11	!
[392]	12	! 388 2009-09-23 09:40:33Z raasch
	13	! Potential density is additionally calculated in eqn_state_seawater,
	14	! first constant in array den also defined as type double.
	15	!
[98]	16	! 97 2007-06-21 08:23:15Z raasch
[96]	17	! Initial revision
	18	!
	19	!
	20	! Description:
	21	! ------------
	22	! Equation of state for seawater as a function of potential temperature,
	23	! salinity, and pressure.
	24	! For coefficients see Jackett et al., 2006: J. Atm. Ocean Tech.
[97]	25	! eqn_state_seawater calculates the potential density referred at hyp(0).
	26	! eqn_state_seawater_func calculates density.
[96]	27	!------------------------------------------------------------------------------!
	28
	29	IMPLICIT NONE
	30
	31	PRIVATE
	32	PUBLIC eqn_state_seawater, eqn_state_seawater_func
	33
	34	REAL, DIMENSION(12), PARAMETER :: nom = &
	35	(/ 9.9984085444849347D2, 7.3471625860981584D0, &
	36	-5.3211231792841769D-2, 3.6492439109814549D-4, &
	37	2.5880571023991390D0, -6.7168282786692354D-3, &
	38	1.9203202055760151D-3, 1.1798263740430364D-2, &
	39	9.8920219266399117D-8, 4.6996642771754730D-6, &
	40	-2.5862187075154352D-8, -3.2921414007960662D-12 /)
	41
	42	REAL, DIMENSION(13), PARAMETER :: den = &
[231]	43	(/ 1.0D0, 7.2815210113327091D-3, &
[96]	44	-4.4787265461983921D-5, 3.3851002965802430D-7, &
	45	1.3651202389758572D-10, 1.7632126669040377D-3, &
	46	-8.8066583251206474D-6, -1.8832689434804897D-10, &
	47	5.7463776745432097D-6, 1.4716275472242334D-9, &
	48	6.7103246285651894D-6, -2.4461698007024582D-17, &
	49	-9.1534417604289062D-18 /)
	50
	51	INTERFACE eqn_state_seawater
	52	MODULE PROCEDURE eqn_state_seawater
	53	MODULE PROCEDURE eqn_state_seawater_ij
	54	END INTERFACE eqn_state_seawater
	55
	56	INTERFACE eqn_state_seawater_func
	57	MODULE PROCEDURE eqn_state_seawater_func
	58	END INTERFACE eqn_state_seawater_func
	59
	60	CONTAINS
	61
	62
	63	!------------------------------------------------------------------------------!
	64	! Call for all grid points
	65	!------------------------------------------------------------------------------!
	66	SUBROUTINE eqn_state_seawater
	67
	68	USE arrays_3d
	69	USE indices
	70
	71	IMPLICIT NONE
	72
	73	INTEGER :: i, j, k
	74
[388]	75	REAL :: pden, pnom, p1, p2, p3, pt1, pt2, pt3, pt4, sa1, sa15, sa2
[96]	76
	77	DO i = nxl, nxr
	78	DO j = nys, nyn
[97]	79	DO k = nzb_s_inner(j,i)+1, nzt
[96]	80	!
	81	!-- Pressure is needed in dbar
[336]	82	p1 = hyp(k) * 1E-4
[96]	83	p2 = p1 * p1
	84	p3 = p2 * p1
	85
	86	!
	87	!-- Temperature needed in degree Celsius
	88	pt1 = pt_p(k,j,i) - 273.15
	89	pt2 = pt1 * pt1
	90	pt3 = pt1 * pt2
	91	pt4 = pt2 * pt2
	92
	93	sa1 = sa_p(k,j,i)
	94	sa15 = sa1 * SQRT( sa1 )
	95	sa2 = sa1 * sa1
	96
[388]	97	pnom = nom(1) + nom(2)pt1 + nom(3)pt2 + &
	98	nom(4)pt3 + nom(5)sa1 + nom(6)sa1pt1 + &
	99	nom(7)*sa2
[96]	100
[388]	101	pden = den(1) + den(2)pt1 + den(3)pt2 + &
	102	den(4)pt3 + den(5)pt4 + den(6)*sa1 + &
	103	den(7)sa1pt1 + den(8)sa1pt3 + den(9)*sa15 + &
	104	den(10)sa15pt2
	105
	106	!
	107	!-- Potential density (without pressure terms)
	108	prho(k,j,i) = pnom / pden
	109
	110	pnom = pnom + nom(8)p1 + nom(9)p1*pt2 + &
	111	nom(10)p1sa1 + nom(11)p2 + nom(12)p2*pt2
	112
	113	pden = pden + den(11)p1 + den(12)p2*pt3 + &
	114	den(13)p3pt1
	115
	116	!
	117	!-- In-situ density
	118	rho(k,j,i) = pnom / pden
	119
[96]	120	ENDDO
[97]	121	!
	122	!-- Neumann conditions are assumed at bottom and top boundary
[388]	123	prho(nzt+1,j,i) = prho(nzt,j,i)
	124	prho(nzb_s_inner(j,i),j,i) = prho(nzb_s_inner(j,i)+1,j,i)
	125	rho(nzt+1,j,i) = rho(nzt,j,i)
	126	rho(nzb_s_inner(j,i),j,i) = rho(nzb_s_inner(j,i)+1,j,i)
	127
[96]	128	ENDDO
	129	ENDDO
	130
	131	END SUBROUTINE eqn_state_seawater
	132
	133
	134	!------------------------------------------------------------------------------!
	135	! Call for grid point i,j
	136	!------------------------------------------------------------------------------!
	137	SUBROUTINE eqn_state_seawater_ij( i, j )
	138
	139	USE arrays_3d
	140	USE indices
	141
	142	IMPLICIT NONE
	143
	144	INTEGER :: i, j, k
	145
[388]	146	REAL :: pden, pnom, p1, p2, p3, pt1, pt2, pt3, pt4, sa1, sa15, sa2
[96]	147
[97]	148	DO k = nzb_s_inner(j,i)+1, nzt
[96]	149	!
	150	!-- Pressure is needed in dbar
[336]	151	p1 = hyp(k) * 1E-4
[96]	152	p2 = p1 * p1
	153	p3 = p2 * p1
	154
	155	!
	156	!-- Temperature needed in degree Celsius
	157	pt1 = pt_p(k,j,i) - 273.15
	158	pt2 = pt1 * pt1
	159	pt3 = pt1 * pt2
	160	pt4 = pt2 * pt2
	161
	162	sa1 = sa_p(k,j,i)
	163	sa15 = sa1 * SQRT( sa1 )
	164	sa2 = sa1 * sa1
	165
[388]	166	pnom = nom(1) + nom(2)pt1 + nom(3)pt2 + &
	167	nom(4)pt3 + nom(5)sa1 + nom(6)sa1pt1 + &
	168	nom(7)*sa2
	169
	170	pden = den(1) + den(2)pt1 + den(3)pt2 + &
	171	den(4)pt3 + den(5)pt4 + den(6)*sa1 + &
	172	den(7)sa1pt1 + den(8)sa1pt3 + den(9)*sa15 + &
	173	den(10)sa15pt2
	174
	175	!
	176	!-- Potential density (without pressure terms)
	177	prho(k,j,i) = pnom / pden
	178
	179	pnom = pnom + nom(8)p1 + nom(9)p1*pt2 + &
	180	nom(10)p1sa1 + nom(11)p2 + nom(12)p2*pt2
	181	pden = pden + den(11)p1 + den(12)p2*pt3 + &
	182	den(13)p3pt1
	183
	184	!
	185	!-- In-situ density
	186	rho(k,j,i) = pnom / pden
	187
	188
[96]	189	ENDDO
[388]	190
[97]	191	!
	192	!-- Neumann conditions are assumed at bottom and top boundary
[388]	193	prho(nzt+1,j,i) = prho(nzt,j,i)
	194	prho(nzb_s_inner(j,i),j,i) = prho(nzb_s_inner(j,i)+1,j,i)
	195	rho(nzt+1,j,i) = rho(nzt,j,i)
	196	rho(nzb_s_inner(j,i),j,i) = rho(nzb_s_inner(j,i)+1,j,i)
[96]	197
	198	END SUBROUTINE eqn_state_seawater_ij
	199
	200
	201	!------------------------------------------------------------------------------!
	202	! Equation of state as a function
	203	!------------------------------------------------------------------------------!
	204	REAL FUNCTION eqn_state_seawater_func( p, pt, sa )
	205
	206	IMPLICIT NONE
	207
	208	REAL :: p, p1, p2, p3, pt, pt1, pt2, pt3, pt4, sa, sa15, sa2
	209
	210	!
	211	!-- Pressure is needed in dbar
	212	p1 = p * 1E-4
	213	p2 = p1 * p1
	214	p3 = p2 * p1
	215
	216	!
	217	!-- Temperature needed in degree Celsius
	218	pt1 = pt - 273.15
	219	pt2 = pt1 * pt1
	220	pt3 = pt1 * pt2
	221	pt4 = pt2 * pt2
	222
	223	sa15 = sa * SQRT( sa )
	224	sa2 = sa * sa
	225
	226
	227	eqn_state_seawater_func = &
	228	( nom(1) + nom(2)pt1 + nom(3)pt2 + nom(4)*pt3 + &
	229	nom(5)sa + nom(6)sapt1 + nom(7)sa2 + nom(8)*p1 + &
	230	nom(9)p1pt2 + nom(10)p1sa + nom(11)p2 + nom(12)p2*pt2 &
	231	) / &
	232	( den(1) + den(2)pt1 + den(3)pt2 + den(4)*pt3 + &
	233	den(5)pt4 + den(6)sa + den(7)sapt1 + den(8)sapt3 + &
	234	den(9)sa15 + den(10)sa15pt2 + den(11)p1 + den(12)p2pt3 + &
	235	den(13)p3pt1 &
	236	)
	237
	238
	239	END FUNCTION eqn_state_seawater_func
	240
	241	END MODULE eqn_state_seawater_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |