Home

Context Navigation

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

Last change on this file since 1350 was 1321, checked in by raasch, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 10.2 KB

Rev	Line
[96]	1	MODULE eqn_state_seawater_mod
	2
[1320]	3	!------------------------------------------------------------------------------!
[1036]	4	! This file is part of PALM.
	5	!
	6	! PALM is free software: you can redistribute it and/or modify it under the terms
	7	! of the GNU General Public License as published by the Free Software Foundation,
	8	! either version 3 of the License, or (at your option) any later version.
	9	!
	10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
	11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
	12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
	13	!
	14	! You should have received a copy of the GNU General Public License along with
	15	! PALM. If not, see <http://www.gnu.org/licenses/>.
	16	!
[1310]	17	! Copyright 1997-2014 Leibniz Universitaet Hannover
[1320]	18	!------------------------------------------------------------------------------!
[1036]	19	!
[484]	20	! Current revisions:
[96]	21	! -----------------
[1321]	22	!
	23	!
	24	! Former revisions:
	25	! -----------------
	26	! $Id: eqn_state_seawater.f90 1321 2014-03-20 09:40:40Z maronga $
	27	!
	28	! 1320 2014-03-20 08:40:49Z raasch
[1320]	29	! ONLY-attribute added to USE-statements,
	30	! kind-parameters added to all INTEGER and REAL declaration statements,
	31	! kinds are defined in new module kinds,
	32	! revision history before 2012 removed,
	33	! comment fields (!:) to be used for variable explanations added to
	34	! all variable declaration statements
[96]	35	!
[1037]	36	! 1036 2012-10-22 13:43:42Z raasch
	37	! code put under GPL (PALM 3.9)
	38	!
[98]	39	! 97 2007-06-21 08:23:15Z raasch
[96]	40	! Initial revision
	41	!
	42	!
	43	! Description:
	44	! ------------
	45	! Equation of state for seawater as a function of potential temperature,
	46	! salinity, and pressure.
	47	! For coefficients see Jackett et al., 2006: J. Atm. Ocean Tech.
[97]	48	! eqn_state_seawater calculates the potential density referred at hyp(0).
	49	! eqn_state_seawater_func calculates density.
[96]	50	!------------------------------------------------------------------------------!
[1320]	51
	52	USE kinds
[96]	53
	54	IMPLICIT NONE
	55
	56	PRIVATE
	57	PUBLIC eqn_state_seawater, eqn_state_seawater_func
	58
[1320]	59	REAL(wp), DIMENSION(12), PARAMETER :: nom = &
	60	(/ 9.9984085444849347D2, 7.3471625860981584D0, &
	61	-5.3211231792841769D-2, 3.6492439109814549D-4, &
	62	2.5880571023991390D0, -6.7168282786692354D-3, &
	63	1.9203202055760151D-3, 1.1798263740430364D-2, &
	64	9.8920219266399117D-8, 4.6996642771754730D-6, &
	65	-2.5862187075154352D-8, -3.2921414007960662D-12 /)
	66	!:
[96]	67
[1320]	68	REAL(wp), DIMENSION(13), PARAMETER :: den = &
	69	(/ 1.0D0, 7.2815210113327091D-3, &
	70	-4.4787265461983921D-5, 3.3851002965802430D-7, &
	71	1.3651202389758572D-10, 1.7632126669040377D-3, &
	72	-8.8066583251206474D-6, -1.8832689434804897D-10, &
	73	5.7463776745432097D-6, 1.4716275472242334D-9, &
	74	6.7103246285651894D-6, -2.4461698007024582D-17, &
	75	-9.1534417604289062D-18 /)
	76	!:
[96]	77
	78	INTERFACE eqn_state_seawater
	79	MODULE PROCEDURE eqn_state_seawater
	80	MODULE PROCEDURE eqn_state_seawater_ij
	81	END INTERFACE eqn_state_seawater
	82
	83	INTERFACE eqn_state_seawater_func
	84	MODULE PROCEDURE eqn_state_seawater_func
	85	END INTERFACE eqn_state_seawater_func
	86
	87	CONTAINS
	88
	89
	90	!------------------------------------------------------------------------------!
	91	! Call for all grid points
	92	!------------------------------------------------------------------------------!
	93	SUBROUTINE eqn_state_seawater
	94
[1320]	95	USE arrays_3d, &
	96	ONLY: hyp, prho, pt_p, rho, sa_p
	97	USE indices, &
	98	ONLY: nxl, nxr, nyn, nys, nzb_s_inner, nzt
[96]	99
	100	IMPLICIT NONE
	101
[1320]	102	INTEGER(iwp) :: i !:
	103	INTEGER(iwp) :: j !:
	104	INTEGER(iwp) :: k !:
[96]	105
[1320]	106	REAL(wp) :: pden !:
	107	REAL(wp) :: pnom !:
	108	REAL(wp) :: p1 !:
	109	REAL(wp) :: p2 !:
	110	REAL(wp) :: p3 !:
	111	REAL(wp) :: pt1 !:
	112	REAL(wp) :: pt2 !:
	113	REAL(wp) :: pt3 !:
	114	REAL(wp) :: pt4 !:
	115	REAL(wp) :: sa1 !:
	116	REAL(wp) :: sa15 !:
	117	REAL(wp) :: sa2 !:
	118
	119
[96]	120
	121	DO i = nxl, nxr
	122	DO j = nys, nyn
[97]	123	DO k = nzb_s_inner(j,i)+1, nzt
[96]	124	!
	125	!-- Pressure is needed in dbar
[336]	126	p1 = hyp(k) * 1E-4
[96]	127	p2 = p1 * p1
	128	p3 = p2 * p1
	129
	130	!
	131	!-- Temperature needed in degree Celsius
	132	pt1 = pt_p(k,j,i) - 273.15
	133	pt2 = pt1 * pt1
	134	pt3 = pt1 * pt2
	135	pt4 = pt2 * pt2
	136
	137	sa1 = sa_p(k,j,i)
	138	sa15 = sa1 * SQRT( sa1 )
	139	sa2 = sa1 * sa1
	140
[1320]	141	pnom = nom(1) + nom(2)pt1 + nom(3)pt2 + &
	142	nom(4)pt3 + nom(5)sa1 + nom(6)sa1pt1 + &
[388]	143	nom(7)*sa2
[96]	144
[1320]	145	pden = den(1) + den(2)pt1 + den(3)pt2 + &
	146	den(4)pt3 + den(5)pt4 + den(6)*sa1 + &
	147	den(7)sa1pt1 + den(8)sa1pt3 + den(9)*sa15 + &
[388]	148	den(10)sa15pt2
	149
	150	!
	151	!-- Potential density (without pressure terms)
	152	prho(k,j,i) = pnom / pden
	153
[1320]	154	pnom = pnom + nom(8)p1 + nom(9)p1*pt2 + &
[388]	155	nom(10)p1sa1 + nom(11)p2 + nom(12)p2*pt2
	156
[1320]	157	pden = pden + den(11)p1 + den(12)p2*pt3 + &
[388]	158	den(13)p3pt1
	159
	160	!
	161	!-- In-situ density
	162	rho(k,j,i) = pnom / pden
	163
[96]	164	ENDDO
[97]	165	!
	166	!-- Neumann conditions are assumed at bottom and top boundary
[388]	167	prho(nzt+1,j,i) = prho(nzt,j,i)
	168	prho(nzb_s_inner(j,i),j,i) = prho(nzb_s_inner(j,i)+1,j,i)
	169	rho(nzt+1,j,i) = rho(nzt,j,i)
	170	rho(nzb_s_inner(j,i),j,i) = rho(nzb_s_inner(j,i)+1,j,i)
	171
[96]	172	ENDDO
	173	ENDDO
	174
	175	END SUBROUTINE eqn_state_seawater
	176
	177
	178	!------------------------------------------------------------------------------!
	179	! Call for grid point i,j
	180	!------------------------------------------------------------------------------!
	181	SUBROUTINE eqn_state_seawater_ij( i, j )
	182
[1320]	183	USE arrays_3d, &
	184	ONLY: hyp, prho, pt_p, rho, sa_p
	185
	186	USE indices, &
	187	ONLY: nzb_s_inner, nzt
[96]	188
	189	IMPLICIT NONE
	190
[1320]	191	INTEGER(iwp) :: i, j, k
[96]	192
[1320]	193	REAL(wp) :: pden, pnom, p1, p2, p3, pt1, pt2, pt3, pt4, sa1, sa15, &
	194	sa2
[96]	195
[97]	196	DO k = nzb_s_inner(j,i)+1, nzt
[96]	197	!
	198	!-- Pressure is needed in dbar
[336]	199	p1 = hyp(k) * 1E-4
[96]	200	p2 = p1 * p1
	201	p3 = p2 * p1
	202
	203	!
	204	!-- Temperature needed in degree Celsius
	205	pt1 = pt_p(k,j,i) - 273.15
	206	pt2 = pt1 * pt1
	207	pt3 = pt1 * pt2
	208	pt4 = pt2 * pt2
	209
	210	sa1 = sa_p(k,j,i)
	211	sa15 = sa1 * SQRT( sa1 )
	212	sa2 = sa1 * sa1
	213
[1320]	214	pnom = nom(1) + nom(2)pt1 + nom(3)pt2 + &
	215	nom(4)pt3 + nom(5)sa1 + nom(6)sa1pt1 + &
[388]	216	nom(7)*sa2
	217
[1320]	218	pden = den(1) + den(2)pt1 + den(3)pt2 + &
	219	den(4)pt3 + den(5)pt4 + den(6)*sa1 + &
	220	den(7)sa1pt1 + den(8)sa1pt3 + den(9)*sa15 + &
[388]	221	den(10)sa15pt2
	222
	223	!
	224	!-- Potential density (without pressure terms)
	225	prho(k,j,i) = pnom / pden
	226
[1320]	227	pnom = pnom + nom(8)p1 + nom(9)p1*pt2 + &
[388]	228	nom(10)p1sa1 + nom(11)p2 + nom(12)p2*pt2
[1320]	229	pden = pden + den(11)p1 + den(12)p2*pt3 + &
[388]	230	den(13)p3pt1
	231
	232	!
	233	!-- In-situ density
	234	rho(k,j,i) = pnom / pden
	235
	236
[96]	237	ENDDO
[388]	238
[97]	239	!
	240	!-- Neumann conditions are assumed at bottom and top boundary
[388]	241	prho(nzt+1,j,i) = prho(nzt,j,i)
	242	prho(nzb_s_inner(j,i),j,i) = prho(nzb_s_inner(j,i)+1,j,i)
	243	rho(nzt+1,j,i) = rho(nzt,j,i)
	244	rho(nzb_s_inner(j,i),j,i) = rho(nzb_s_inner(j,i)+1,j,i)
[96]	245
	246	END SUBROUTINE eqn_state_seawater_ij
	247
	248
	249	!------------------------------------------------------------------------------!
	250	! Equation of state as a function
	251	!------------------------------------------------------------------------------!
[1320]	252	REAL(wp) FUNCTION eqn_state_seawater_func( p, pt, sa )
[96]	253
	254	IMPLICIT NONE
	255
[1320]	256	REAL(wp) :: p !:
	257	REAL(wp) :: p1 !:
	258	REAL(wp) :: p2 !:
	259	REAL(wp) :: p3 !:
	260	REAL(wp) :: pt !:
	261	REAL(wp) :: pt1 !:
	262	REAL(wp) :: pt2 !:
	263	REAL(wp) :: pt3 !:
	264	REAL(wp) :: pt4 !:
	265	REAL(wp) :: sa !:
	266	REAL(wp) :: sa15 !:
	267	REAL(wp) :: sa2 !:
[96]	268
	269	!
	270	!-- Pressure is needed in dbar
	271	p1 = p * 1E-4
	272	p2 = p1 * p1
	273	p3 = p2 * p1
	274
	275	!
	276	!-- Temperature needed in degree Celsius
	277	pt1 = pt - 273.15
	278	pt2 = pt1 * pt1
	279	pt3 = pt1 * pt2
	280	pt4 = pt2 * pt2
	281
	282	sa15 = sa * SQRT( sa )
	283	sa2 = sa * sa
	284
	285
	286	eqn_state_seawater_func = &
	287	( nom(1) + nom(2)pt1 + nom(3)pt2 + nom(4)*pt3 + &
	288	nom(5)sa + nom(6)sapt1 + nom(7)sa2 + nom(8)*p1 + &
	289	nom(9)p1pt2 + nom(10)p1sa + nom(11)p2 + nom(12)p2*pt2 &
	290	) / &
	291	( den(1) + den(2)pt1 + den(3)pt2 + den(4)*pt3 + &
	292	den(5)pt4 + den(6)sa + den(7)sapt1 + den(8)sapt3 + &
	293	den(9)sa15 + den(10)sa15pt2 + den(11)p1 + den(12)p2pt3 + &
	294	den(13)p3pt1 &
	295	)
	296
	297
	298	END FUNCTION eqn_state_seawater_func
	299
	300	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 |