• API
• FAQ
• Tools
• Archive
daily pastebin goal
57%
SHARE
TWEET

# Untitled

a guest Dec 13th, 2018 59 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. import numpy as np
2. def K0_Weiss(S, TC):
3.     TK=TC+273.15 #Temperatura de Celsius a Kelvin
4.     lnK0 = 9345.17/TK - 60.2409 + 23.3585 * np.log(TK/100) + S * (0.023517 - 0.00023656 * TK + 4.7036e-07 *
5.         TK * TK)
6.     K0=np.exp(lnK0)
7.     return K0
8.
9.
10. """
11. Created on Mon Sep  3 11:51:00 2018
12.
13. @author: fergomez
14. Esta funcion calcula K1 y K2 del sistema CO2-H20 segun
15. Roy et al., (1993), The dissociation constant of carbonic acid in sewater
16. at salinities of 5 to 45 and temperatures of 0 to 45 desgrees Celsius.
17. Marine Chemistry 44(2-4), 249-267.
19. K1, K2=mol*kg-soln-1
20. Ejemplo de uso:
21. Input
22. K1, K2=K1_K2_Roy(35, 25)
23. print(K1, K2)
24. Output:
25. 1.3921075396202872e-06 1.1887254858040348e-09
26. """
27.
28.
29. def K1_K2_Roy(S, TC):
30.     TK=TC+273.15 #Temperatura de Celsius a Kelvin
31.     #K1 usando Roy et al., 1993
32.     tmp1 = 2.83655 - 2307.1266/TK - 1.5529413 * np.log(TK)
33.     tmp2 = -(0.207608410 + 4.0484/TK) * np.sqrt(S)
34.     tmp3 = 0.08468345 * S - 0.00654208 * S**1.5
35.     tmp4 = np.log(1 - 0.001005 * S)
36.     lnK1roy = tmp1 + tmp2 + tmp3 + tmp4
37.     K1 = np.exp(lnK1roy)
38.     #K2 usando Roy et al., 1993
39.     tmp1 = -9.226508 - 3351.6106/TK - 0.2005743 * np.log(TK)
40.     tmp2 = (-0.106901773 - 23.9722/TK) * np.sqrt(S)
41.     tmp3 = 0.1130822 * S - 0.00846934 * S**1.5 + np.log(1 -
42.         0.001005 * S)
43.     lnK2roy = tmp1 + tmp2 + tmp3
44.     K2 = np.exp(lnK2roy)
45.     return K1, K2
46.
47. def K1_K2_SBY(S):
48.     pK1=6.1568-0.00352*S
49.     pK2=8.5503-0.0080*S
50.     K1=10**(-pK1)
51.     K2=10**(-pK2)
52.     return K1, K2
53. """
54. Created on Mon Sep  3 11:53:53 2018
55.
56. @author: fergomez
57. Esta funcion calcula K0 del sistema CO2-H20 segun recomienda Dickson et al 2007
58. Guide for best practices for ocean CO2 measurements. PICES SP 3, 191 pp.
59. Conforme a Mucci, A., 1983. The solubility of calcite and aragonite in seawater
60. at various salinities and temperatures and one atmosphere of total pressure.
61. American journal of Science, 283, 780-799.
62.
63. Ejemplo de uso:
64. Input
65. Ksp_a=Kspa_Mucci(35, 25)
66. print(Ksp_a)
67. Output:
68. 6.48175906801198e-07
69. """
70.
71. def Kspa_Mucci(S, TC):
72.     TK=TC+273.15 #Temperatura de Celsius a Kelvin
73.     tmp1 = -171.945 - 0.077993 * TK + 2903.293/TK + 71.595 * np.log10(TK)
74.     tmp2 = +(-0.068393 + 0.0017276 * TK + 88.135/TK) * np.sqrt(S)
75.     tmp3 = -0.10018 * S + 0.0059415 * S**1.5
76.     log10Kspa = tmp1 + tmp2 + tmp3
77.     Ksp_a=10**(log10Kspa)
78.     return Ksp_a
79.
80.
81. """
82. Created on Mon Sep  3 11:54:42 2018
83.
84. @author: fergomez
85.
86. Esta funcion calcula K0 del sistema CO2-H20 segun recomienda Dickson et al 2007
87. Guide for best practices for ocean CO2 measurements. PICES SP 3, 191 pp.
88. Conforme a Mucci, A., 1983. The solubility of calcite and aragonite in seawater
89. at various salinities and temperatures and one atmosphere of total pressure.
90. American journal of Science, 283, 780-799.
91.
92. Ejemplo de uso:
93. Input
94. Ksp_c=Kspc_Mucci(35, 25)
95. print(Ksp_c)
96. Output:
97. 4.2723509278626e-07
98. """
99.
100. def Kspc_Mucci(S, TC):
101.     TK=TC+273.15 #Temperatura de Celsius a Kelvin
102.     tmp1 = -171.9065 - 0.077993 * TK + 2839.319/TK + 71.595 * np.log10(TK)
103.     tmp2 = +(-0.77712 + 0.0028426 * TK + 178.34/TK) * np.sqrt(S)
104.     tmp3 = -0.07711 * S + 0.0041249 * S**1.5
105.     log10Kspc = tmp1 + tmp2 + tmp3
106.     Ksp_c=10**(log10Kspc)
107.     return Ksp_c
108.
109. # -*- coding: utf-8 -*-
110. """
111. Created on Mon Sep  3 11:55:56 2018
112. Calcula el factor de correccion de la presion para las constantes de
113. equilibrio Ki usando Millero, F.J., (1995). Thermodynamics of the carbon
114. dioxide system in the oceans. Geochimica et Cosmochimica Acta 59:661-677.
115. @author: fergomez
116.
117. Ejemplo de uso:
118. A presion atmosferica Kspa=4.27 e-07
119. Input (para aragonita):
120. fcorr=P_corr_Millero(-45.96, 0.5304, 0.0, -11.76e-03, 0.3692e-03, 0.0, 25, 300)
121. print(fcorr)
122. Output:
123. 1.4787577790538065
124.
125. entonces Kspa(presion atmosferica)*fcorr=Kspa_corregido
126. 4.27 e-07*1.4787577790538065= 9.58 e-07
127. """
128.
129.
130. #pagina 186 Andreas Hoffman PhD Thesis
131. def P_corr_Millero(a0, a1, a2, b0, b1, b2, TC, P):
132.     TK=TC+273.15 #Temperatura de Celsius a Kelvin
133.     R=83.131
134.     delta_V= a0+a1*TC+a2*TC*TC
135.     delta_k= (b0+b1*TC+b2*TC*TC)
136.     lnKp=-(delta_V/(R * TK))*P+0.5*(delta_k/(R*TK))*P*P
137.     fcorr=np.exp(lnKp)
138.     return fcorr
139.
140. #Llamamos las librerias que utilizaremos
141. import numpy as np
142. #import matplotlib.pyplot as plt
143. import math
144.
145. from numpy import exp,arange
146. from pylab import meshgrid,cm,imshow,contour, clabel,colorbar,axis,title,show, contourf
147.
148. #from constantes import K0_Weiss, K1_K2_Roy, Kspc_Mucci, Kspa_Mucci, P_corr_Millero
149. import matplotlib.pyplot as plt
150.
151. #============================================================================
152. #Inicio de seccion de INPUT==================================================
153. #============================================================================
154.
155. S=6.88 #ppt
156. TC=6.42 #Celsius
157. P=1.0 #atmosferas
158.
159.
160.
161. K0=K0_Weiss(S, TC)
162. K1, K2=K1_K2_Roy(S, TC)
163. Ksp_c=Kspc_Mucci(S, TC)
164. Ksp_a=Kspa_Mucci(S, TC)
165.
166. #============================================================================
167. #Fin de seccion de INPUT======================================================
168. #============================================================================
169.
170. #=============================================================================
171. #Estimacion de Calcio de Tyrrell et al., 2008 y valores de ejemplo agua de mar
172. #============================================================================
173. #Ca=0.0103 #valores agua de mar
174. Ca=(0.331*S+0.392)*1e-03 #Mar Baltico
175. #Ca=(0.375*S+0.0368)*1e-03 # Bothnian Bay (parte norte del mar Baltico)
176. #=============================================================================
177.
178. #def co2_t(pco2, K0):
179. #    co2=pco2*K0
180. #    return co2
181. #=============================================================================
182. #vco2_t = np.vectorize(co2_t)
183. #Pasamos de pco2 a co2
184. #co2=co2_t(pco2, K0)
185.
186.
187. #=============================================================================
188. #=============================================================================
189. #Calculo del factor de correccion por presion
190. #=============================================================================
191. fcorr_K1=P_corr_Millero(-25.50, 0.1271, 0.0, -3.08e-03, 0.0877e-03, 0.0, TC, P)
192. fcorr_K2=P_corr_Millero(-15.82, -0.219, 0.0, 1.136e-03, -0.1475e-03, 0.0, TC, P)
193. fcorr_Kspc=P_corr_Millero(-48.76, 0.5304, 0.0, -11.76e-03, 0.3692e-03, 0.0, TC, P)
194. fcorr_Kspa=P_corr_Millero(-45.96, 0.5304, 0.0, -11.76e-03, 0.3692e-03, 0.0, TC, P)
195. #=============================================================================
196.
197. #=============================================================================
198. #Correccion de ctes de equilibrio por presion
199. #=============================================================================
200. K1=K1*fcorr_K1
201. K2=K2*fcorr_K2
202. Ksp_c=Ksp_c*fcorr_Kspc
203. Ksp_a=Ksp_a*fcorr_Kspa
204.
205. #=============================================================================
206. #Funcion que realiza calculos del sistema carbonato y omega
207. #=============================================================================
208.
209. def carbEq(co2, dic, Ca):
210.     #-----------------------------------------------
211.     # Resolvemos para obtener H+ (cf. a Zeebe and Wolf-Gladrow, 2000)
212.     a1=co2-dic
213.     a2=K1*co2
214.     a3=K1*K2*co2
215.
216.     p= [a1, a2, a3]
217.     r = np.roots(p)
218.     h= max(np.real(r)) # Esto es para seleccionar la raiz real mas grande
219.     #
220.     # Calculamos HCO3, CO3 and CO2aq, usando DIC, AlK y H+
221.     hco3=dic/(1+h/K1+K2/h)
222.     co3=dic/(1+h/K2+h*h/(K1*K2))
223.     #co2=dic/(1+K1/h+K1*K2/(h*h))
224.     fco2=co2 / K0
225.     pH=-math.log10(h)
226.     alk=2*co3+hco3
227.
228.     #Saturacion de Calcita y Aragonita
229.     omega_ar=Ca*co3/Ksp_a
230.     omega_cal=Ca*co3/Ksp_c
231.
232.     return fco2, pH, co2, hco3, co3, h, dic, alk, omega_ar, omega_cal
233.
234. vcarbEq = np.vectorize(carbEq)
235. #=============================================================================
236. #Uso de la funcion - ejemplo
237. #=============================================================================
238. #co2=co2*1000000
239. #dic=dic*1000000
240.
241. pco2=np.linspace(0.000090, 0.001100, 100)
242. co2=pco2*K0
243. dic=np.linspace(0.001526, 0.002100, 100) #moles (para pasar de micromoles a moles, 1587*e-06)
244.
245. #pco2=pco2*1000
246. #dic=dic*1000
247. co2, dic=meshgrid(co2*1000000, dic*1000000)
248.
249. fco2, pH, co2, hco3, co3, h, dic, alk, omega_ar, omega_cal=vcarbEq(co2=co2, dic=dic, Ca=Ca)
250.
251.
252. #Ejemplo con valores de referencia de uso para agua de mar normal
253. #fco2, pH, co2, hco3, co3, h, dic, alk, omega_ar, omega_cal=carbEq(co2=0.00001032997, dic=0.002108, Ca=Ca)
254. #para este usar calcio=0.0103
255. #=============================================================================
257. #=============================================================================
258. print(co2)
259. print('====================================')
260. print(dic)
261. print('====================================')
262. print(pH)
263. #print(omega_cal)
264. #fig, ax = plt.subplots()
265. im = imshow(pH, cmap=cm.RdBu) # dibujo la funcion
266. plt.xlabel('CO2')
267. plt.ylabel('DIC')
268. plt.ylim(0, 100)
269.
270. plt.scatter(40, 40, color='k')
271. #ax.set_xlim(0.00000525, 0.000065)
272. #ax.set_ylim(0.001, 0.0021)
273. #plt.xlim(min(co2), max(co2))
274. # agrego lineas de contorno y rotulos
275. cset = contour(pH, arange(6.0,9.0,0.2),linewidths=2, cmap=cm.Set2)
276. clabel(cset,inline=True,fmt='%1.1f',fontsize=10)
277. colorbar(im) # agrego la barra de colores al costado
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.

Top