Context Navigation

source: trunk/SpectralModelF90/PROJECTS/trivial/equation_of_state.f @ 6

Last change on this file since 6 was 6, checked in by xlvlod, 18 years ago
import initial from SVN_BASE_TRUNK
File size: 6.5 KB

Line
1	subroutine equation_of_state(
2	* scalars,s1,s2,s1_scale, s2_scale,
3	* nx, ny, locnz, rho_0, DGRAD,
4	* myid, numprocs, Lx, Ly, Lz, pd)
5
6	! Routine to calculate d'less pd from values of d'less s1 and s2
7	! and their reference profiles s1_bar(x,z) and s2_bar(x,z)
8	! The "perturbation" density is defined relative to the reference
9	! profile rho_bar(x,z).
10	!======================================================================
11	! scalars='TS'
12	! EOS from MILLERO et.al. (1980) DEEP-SEA RES.
13	! The scalar s1=temparature in deg C and
14	! the scalar s2 is salinity in psu (i.e. ~ 36)
15
16	implicit none
17	character*2 scalars
18	integer nx, ny, locnz, myid, numprocs, nyplanes
19	real Lx, Ly, Lz
20	real s1(nx+2,ny+1,locnz+1),s2(nx+2,ny+1,locnz+1),pd(nx+2,ny+1,locnz+1)
21	real s1_scale, s2_scale, rho_0, DGRAD
22
23	! local variables
24	! S,T,P all DOUBLE precision arguments of function rho: P in MPa, T in deg C, S in conc units
25	double precision Temp,Salinity
26	double precision P_Mpa,density
27	double precision rho_bar,total_density
28	real z,depth,s1_bar,s2_bar,dx,dy,dz
29	real junk,xprime,yprime,zprime,xx,yy
30	integer i, j, k, kend
31
32	if(ny .eq. 0) then
33	nyplanes = 1
34	else
35	nyplanes = ny
36	endif
37
38	P_Mpa = 0.0 ! compute potential density, not in situ density
39
40	if( scalars .eq. 'TS' ) then
41	xx=1.0 ! use this to keep Salinity
42	yy=1.0 ! use this to keep Temperature
43	elseif( scalars .eq. 'SP' .or. scalars .eq. 'S0' ) then
44	xx=1.0 ! use this to keep Salinity
45	yy=0.0 ! use this to neglect Temperature
46	elseif( scalars .eq. 'TP' .or. scalars .eq. 'T0' ) then
47	yy=1.0 ! use this to keep Temperature
48	xx=0.0 ! use this to neglect Salinity
49	endif
50
51	kend=locnz
52	if( myid .eq. numprocs-1) kend=locnz+1
53
54	dx=Lx/float(nx) ! [m]
55	dy=Ly/float(nyplanes) ! [m]
56	dz=Lz/float(locnz*numprocs) ! [m]
57
58	if( scalars .eq. 'TS'
59	* .or. scalars .eq. 'TP' .or. scalars .eq. 'SP'
60	* .or. scalars .eq. 'T0' .or. scalars .eq. 'S0' ) then
61
62	do k=1,kend
63	z = (myidlocnz+k-1)dz ![m]
64	zprime=z/Lz ! d'less z
65	do j=1,nyplanes
66	yprime=(j-1.)*dy/Lz
67	do i=1,nx
68	xprime=(i-1)*dx/Lz
69	call s1_bar_func(xprime,yprime,zprime,s1_bar,
70	* junk,junk,junk,junk,
71	* Lz,s1_scale)
72	call s2_bar_func(xprime,yprime,zprime,s2_bar,
73	* junk,junk,junk,junk,
74	* Lz,s2_scale)
75	Temp = yy(s1_bars1_scale) ! [deg C]
76	Salinity = (s2_bar*s2_scale) ! [psu]
77	Salinity=xx*Salinity/1000. ! fortran routine rho uses concentration units
78	call rho(Salinity,Temp,P_mpa,density)
79	rho_bar = density ! ambient density [kg/m3]
80
81	! now repeat the calculation for total s1 & s2
82	Temp = yy*((s1_bar
83	* + s1(i,j,k))*s1_scale) ! dimensional quantity: total Temp.
84	Salinity = ((s2_bar
85	* + s2(i,j,k))*s2_scale) ! dimensional quantity: total Salinity
86	Salinity=xx*Salinity/1000. ! fortran routine rho uses concentration units
87
88	call rho(Salinity,Temp,P_Mpa,density) ! total density [kg/m3]
89	total_density = density ! [kg/m3]
90
91	pd(i,j,k) = total_density - rho_bar ! perturbation in kg/m3
92	pd(i,j,k) = pd(i,j,k)/(Lz*DGRAD) ! dimensionless perturbation density
93
94	enddo
95	enddo
96	enddo
97
98	elseif( scalars .eq. 'RP' .or. scalars .eq. 'R0' ) then ! s1 is pd
99
100	do k=1,kend
101	do j=1,nyplanes
102	do i=1,nx
103	pd(i,j,k)=s1(i,j,k)
104	enddo
105	enddo
106	enddo
107
108	elseif( scalars .eq. 'PP' ) then ! both scalars passive, set pd=0
109
110	do k=1,kend
111	do j=1,nyplanes
112	do i=1,nx
113	pd(i,j,k)=0.0
114	enddo
115	enddo
116	enddo
117
118	endif
119
120	end subroutine equation_of_state
121
122
123	! "rho.f"
124
125	! Parameter : equation of state
126	! Source : MILLERO et.al. (1980) DEEP-SEA RES.
127	! Vol 27A, pp 255-264
128	! Units : Kg / m**3
129	! Input : s (concentration units), t (deg C), p (MPa)
130	! Range : s (0 to 0.040), t (-2 to 40), p (MPa)
131	! Ck value : K (0, 0, 100) = 2.297721e4
132	! K (0.035, 0, 100) = 2.499200e4
133	! K (0, 25, 100) = 2.54051e4
134	! K (0.035, 25, 100) = 2.710895e4
135	! Coded : Ngoc Dang, July 1982
136
137	subroutine rho (s, t, p, density)
138
139	double precision s, t, p
140	double precision a, b, c, d, e, k, t1, t2, t3, t4, ta, tb
141	double precision k0, k0w, alpha0, alpha, rho0, aw, bw, density
142
143
144	! convert input units
145
146	s = 1000.0d0 * s
147	p = 10.0d0 * p
148
149	! compute alpha0 = alpha (s, t, 0)
150
151	t1 = 9.99841594d-1 + t * (6.793952d-5 + t * (-9.095290d-6 + t *
152	& (1.001685d-7 + t * (-1.120083d-9 + t * 6.536332d-12))))
153	t2 = 8.25917d-4 + t * (-4.4490d-6 + t * (1.0485d-7 + t *
154	& (-1.258d-9 + t * 3.315d-12)))
155	t3 = -6.33761d-6 + t * (2.8441d-7 + t * (-1.6871d-8 + t *
156	& 2.83258d-10))
157	t4 = 5.4705d-7 + t * (-1.97975d-8 + t * (1.6641d-9 - t *
158	& 3.1203d-11))
159	rho0 = t1 + t2 * s + t3 * s*(1.5) + t4 s**2
160	alpha0 = 1.0d0 / rho0
161
162	! compute the pressure coefficient
163
164	! the first term
165	k0w = 1.965221d4 + t * (1.484206d2 + t * (-2.327105d0 + t *
166	& (1.360477d-2 - t * 5.155288d-5)))
167	ta = 5.46746d1 + t * (-6.03459d-1 + t * (1.09987d-2 - t *
168	& 6.167d-5))
169	tb = 7.944d-2 + t * (1.6483d-2 - t * 5.3009d-4)
170	k0 = k0w + ta * s + tb * s**(1.5)
171
172	! the second term
173	aw = 3.239908d0 + t * (1.43713d-3 + t * (1.16092d-4 - t *
174	& 5.77905d-7))
175	c = 2.2838d-3 + t * (-1.0981d-5 - t * 1.6078d-6)
176	d = 1.9107d-4
177	a = aw + c * s + d * s**(1.5)
178
179	! the third term
180	bw = 8.50935d-5 + t * (-6.12293d-6 + t * 5.2787d-8)
181	e = -9.9348d-7 + t * (2.0816d-8 + t * 9.1697d-10)
182	b = bw + e * s
183
184	k = k0 + a * p + b * p**2
185
186	! compute density
187	alpha = alpha0 * (1.0d0 - p / k)
188	density = 1.0d0 / alpha
189
190	! convert to kg / m**3
191	density = 1000.0d0 * density
192	s = s / 1000.0d0
193	p = p / 10.0d0
194
195	return
196	end

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

Download in other formats: