Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

ice_gas_solu.f @ 24

Last change on this file since 24 was 4, checked in by vancop, 8 years ago
initial import /Users/ioulianikolskaia/Boulot/CODES/LIM1D/ARCHIVE/TMP/LIM1D_v3.20/
File size: 11.7 KB

Line
1	SUBROUTINE ice_gas_solu
2
3	!------------------------------------------------------------------------------!
4	! * ice_gas_solu *
5	!
6	! This routine computes for Argon, Oxygen, CO2 and Nitrogen
7	! - Saturation concentration (csat_gas, mmol/m3)
8	! - Solubility (sol_gas, mmol/m3/atm)
9	! in each vertical layer
10	!
11	! Refs
12	! Argon: Hamme and Emerson (2004)
13	! Oxygen: Garcia and Gordon (1992)
14	! CO2: Sarmiento and Gruber (2006), Weiss (1974)
15	! N2: Sarmiento and Gruber (2006) and references therein
16	!
17	! Original code (in ice_gas.f): S. Moreau, M. Vancoppenolle, 2012-2015
18	!
19	! Rewriting: M. Vancopppenolle, Nov 2015
20	!
21	! Present status :
22	! - code for Ar, not for N2 and O2
23	! - should put mixing ratios in the namelist
24	!------------------------------------------------------------------------------!
25
26	USE lib_fortran
27
28	INCLUDE 'type.com'
29	INCLUDE 'para.com'
30	INCLUDE 'const.com'
31	INCLUDE 'ice.com'
32	INCLUDE 'thermo.com'
33	INCLUDE 'bio.com'
34
35	REAL(8), DIMENSION(nlay_bio) ::
36	& ztc, ! Celsius temperature
37	& ztlog ! Re-scaled temperature
38
39	LOGICAL ::
40	& ln_write_gas = .TRUE.
41
42	!==============================================================================!
43
44	!
45	!------------------------------------------------------------------------------!
46	! X) Starting the routine
47	!------------------------------------------------------------------------------!
48	!
49
50	IF ( ln_write_gas ) THEN
51
52	WRITE(numout,*)
53	WRITE(numout,) ' ** ice_gas_solu : '
54	WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~'
55	WRITE(numout,*)
56
57	ENDIF
58
59	ji = 1
60	zpatm = psbqb(ji) / 101325. ! atmospheric pressure in atmospheres
61
62	CALL ice_brine
63
64	IF ( ln_write_gas ) THEN
65	WRITE(numout,*) ' ztc : ', ztc(1:nlay_bio)
66	WRITE(numout,*) ' sbr_bio : ', sbr_bio(1:nlay_bio)
67	WRITE(numout,*) ' rhobr_bio: ', rhobr_bio(1:nlay_bio)
68	WRITE(numout,*)
69	ENDIF
70
71	DO jn = 1, ntra_bio
72
73	IF ( flag_active(jn) .AND. ( biotr_i_typ(jn) .EQ. 'gas' ) )
74	& THEN
75
76	!
77	!------------------------------------------------------------------------------!
78	! X) Celsius temperature and rescaled temperature
79	!------------------------------------------------------------------------------!
80	!
81
82	DO layer = 1, nlay_bio
83
84	ztc(layer) = tc_bio(layer) ! temperature in celsius
85
86	ztlog(layer) = LOG( ( 298.15 - ztc(layer) ) / ! rescaled log temperature
87	& ( 273.15 + ztc(layer) ) )
88
89	END DO
90
91	!
92	!------------------------------------------------------------------------------!
93	! X) Saturation concentrations
94	!------------------------------------------------------------------------------!
95	! !-------
96	IF ( biotr_i_nam(jn) .EQ. 'Arg' ) THEN ! Argon
97	!-------
98
99	IF ( ln_write_gas ) THEN
100	WRITE(numout,*) ' --- Argon ...'
101	ENDIF
102
103	! Hamme and Emerson, 2004
104	! form ln C = AO + A1 * Ts + A2 * Ts^2 + A3 * Ts^3
105	! + S ( B0 + B1 * Ts + B2 * Ts^2 )
106	! Ts = f(Tc), Tc in Celsius, S in pss
107
108	!--- Polynomial coefficients (mumol/kg, Table 4)
109	za0 = 2.79150 ; za1 = 3.17609 ; za2 = 4.13116
110	za3 = 4.90379
111	zb0 = -6.96233e-3 ; zb1 = -7.66670e-3 ; zb2 = -1.16888e-2
112
113	DO layer = 1, nlay_bio
114
115	!--- Rescaled log temperature powers
116	zt = ztlog(layer)
117	zt2 = zt * zt
118	zt3 = zt * zt2
119
120	!--- Saturation concentration at equilibrium with a moist atmosphere at 1-atm pressure
121	zln_csat = za0 + za1zt + za2zt2 + za3*zt3 +
122	& ( zb0 + zb1zt + zb2zt2 ) * sbr_bio(layer)
123
124	csat_gas(jn,layer) = EXP( zln_csat ) ! in mumol/kg
125
126	csat_gas(jn,layer) = csat_gas(jn,layer) * ! in mmol/m3
127	& rhobr_bio(layer) / 1000.
128
129	END DO
130
131	ENDIF ! 'Arg'
132
133	!-----------
134	IF ( biotr_i_nam(jn) .EQ. 'Oxy' ) THEN ! di-Oxygen
135	!-----------
136
137	IF ( ln_write_gas ) THEN
138	WRITE(numout,*) ' --- Oxygen ...'
139	ENDIF
140
141	! Garcia and Gordon, L&O 1992
142	! form ln C = f(S,T), eq 8 page 1310 (! small mistake in the equation, extra-term)
143
144	!--- Polynomial coefficients (mumol/kg, Table 1, last column)
145	za0_sat_oxy = 5.80818; za1_sat_oxy = 3.20684
146	za2_sat_oxy = 4.11890; za3_sat_oxy = 4.93845
147	za4_sat_oxy = 1.01567; za5_sat_oxy = 1.41575
148	zb0_sat_oxy = -7.01211e-3; zb1_sat_oxy = -7.25958e-3
149	zb2_sat_oxy = -7.933343e-3; zb3_sat_oxy = -5.54491e-3
150	zc0_sat_oxy = -1.32412e-7
151
152	DO layer = 1, nlay_bio
153
154	!--- Rescaled log temperature powers
155	zt = ztlog(layer)
156	zt2 = zt * zt
157	zt3 = zt2 * zt
158	zt4 = zt2 * zt2
159	zt5 = zt2 * zt3
160
161	!--- Saturation concentration at equilibrium with a moist atmosphere at 1-atm pressure
162	zln_csat = ( za0_sat_oxy + za1_sat_oxy * zt +
163	& za2_sat_oxy * zt2 +
164	& za3_sat_oxy * zt3 +
165	& za4_sat_oxy * zt4 +
166	& za5_sat_oxy * zt5 +
167	& sbr_bio(layer) *
168	& ( zb0_sat_oxy + zb1_sat_oxy * zt +
169	& zb2_sat_oxy * zt2 +
170	& zb3_sat_oxy * zt3 ) +
171	& zc0_sat_oxy *
172	& sbr_bio(layer) * sbr_bio(layer) )
173
174	csat_gas(jn,layer) = EXP(zln_csat) ! mumol/kg
175
176	csat_gas(jn,layer) = csat_gas(jn,layer) * ! in mmol/m3
177	& rhobr_bio(layer) / 1000.
178
179	END DO
180
181	ENDIF ! 'Oxy'
182
183	!-------------
184	IF ( biotr_i_nam(jn) .EQ. 'Nit' ) THEN ! di-Nitrogen
185	!-------------
186
187	IF ( ln_write_gas ) THEN
188	WRITE(numout,*) ' --- Nitrogen ...'
189	ENDIF
190
191	! Sarmiento and Gruber (2006, page 74) and references therein
192	! form ...
193
194	!--- Polynomial coefficients (mumol/kg)
195	za0_sat_nit = 6.42931
196	za1_sat_nit = 2.92704
197	za2_sat_nit = 4.32531
198	za3_sat_nit = 4.69149
199	zb0_sat_nit = -0.00744129
200	zb1_sat_nit = -0.00802566
201	zb2_sat_nit = -0.0146775
202
203	DO layer = 1, nlay_bio
204
205	!--- Rescaled log temperature powers
206	zts = ztlog(layer)
207	zts2 = zts * zts
208	zts3 = zts2 * zts
209
210	zln_csat = za0_sat_nit + za1_sat_nit * zts +
211	& za2_sat_nit * zts2 +
212	& za3_sat_nit * zts3
213	& + sbr_bio(layer) *
214	& ( zb0_sat_nit + zb1_sat_nit * zts +
215	& zb2_sat_nit * zts2 )
216
217	csat_gas(jn,layer) = EXP(zln_csat) ! mumol/kg
218
219	csat_gas(jn,layer) = csat_gas(jn,layer) * ! in mmol/m3
220	& rhobr_bio(layer) / 1000.
221
222	END DO ! layer
223
224	ENDIF ! 'Nit'
225
226	!-------------
227	IF ( biotr_i_nam(jn) .EQ. 'CO2' ) THEN ! CO2
228	!-------------
229
230	IF ( ln_write_gas ) THEN
231	WRITE(numout,*) ' --- CO2 ...'
232	ENDIF
233
234	! Sarmiento and Gruber, 2006 (Weiss, 1974 with total ph scale)
235	! form C = mix_rat * exp[ f(S,T) ], Table 3.2.2
236
237	za0_sat_CO2 = -160.7333; za1_sat_CO2 = 215.4152
238	za2_sat_CO2 = 89.8920; za3_sat_CO2 = -1.47759
239
240	zb0_sat_CO2 = 0.029941; zb1_sat_CO2 = -0.027455
241	zb2_sat_CO2 = 0.0053407
242
243	DO layer = 1, nlay_bio
244
245	!--- temperature powers
246	zt = t_i_bio(layer) ! Temperature (K)
247	z100_t = 100. / zt
248	zt_100 = zt / 100.
249	zt_100_2 = zt_100 * zt_100
250
251	zsat_CO2 = mixr_gas(jn) *
252	& EXP ( za0_sat_CO2 +
253	& za1_sat_CO2 * z100_t + ! zsat_CO2 is the CO2 concentration at saturation
254	& za2_sat_CO2 * log (zt_100) + ! zsat_CO2 is in mol/kg
255	& za3_sat_CO2 * zt_100_2 +
256	& sbr_bio(layer) *
257	& ( zb0_sat_CO2 +
258	& zb1_sat_CO2 * zt_100 +
259	& zb2_sat_CO2 * zt_100_2 ) )
260
261	csat_gas(jn,layer) = zsat_CO2 * rhobr_bio(layer) * 1000. ! mmol/m3
262
263	! WRITE(numout,*) '--------------------------- '
264	! WRITE(numout,*) ' jn, layer: ', jn,layer
265	! WRITE(numout,*) ' mixr_gas : ', mixr_gas(jn)
266	! WRITE(numout,*) ' rhobr_bio: ', rhobr_bio(layer)
267	! WRITE(numout,*) ' zsat_CO2 : ', zsat_CO2
268	! WRITE(numout,*) ' zt : ', zt
269	! WRITE(numout,*) ' z100_t : ', z100_t
270	! WRITE(numout,*) ' zt_100 : ', zt_100
271	! WRITE(numout,*) ' zt_100_2 : ', zt_100_2
272	! WRITE(numout,*) ' csat_gas : ', csat_gas(jn,layer)
273	! WRITE(numout,*) '--------------------------- '
274
275	END DO ! layer
276
277	ENDIF ! 'CO2'
278
279	WRITE(numout,*) ' ---> csat_gas: ', csat_gas(jn,:)
280	WRITE(numout,*)
281
282	!------------------------------------------------------------------------------!
283
284	ENDIF ! flags
285	END DO ! jn
286
287	!
288	!------------------------------------------------------------------------------!
289	! X) Solubility and partial pressure
290	!------------------------------------------------------------------------------!
291	!
292	DO jn = 1, ntra_bio
293
294	sol_gas(jn,:) = 0.
295
296	IF ( flag_active(jn) .AND. ( biotr_i_typ(jn) .EQ. 'gas' ) )
297	& THEN
298
299	sol_gas(jn,:) = csat_gas(jn,:) / mixr_gas(jn) ! solubility, ./m3/atm
300
301	IF ( ln_write_gas ) THEN
302	WRITE(numout,*) ' Tracer number : ', jn
303	WRITE(numout,*) ' --- ', biotr_i_nam(jn)
304	WRITE(numout,*) ' sol_gas(jn,:) : ', sol_gas(jn,:)
305	ENDIF
306
307	ENDIF
308
309	END DO
310	!
311	!------------------------------------------------------------------------------!
312	! X) End of the routine
313	!------------------------------------------------------------------------------!
314	!
315	WRITE(numout,*)
316	WRITE(numout,*) ' End of ice_gas_solu '
317	WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
318	WRITE(numout,*)
319
320	RETURN
321
322	END

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

Context Navigation

source: branches/2016/dev_v3.20_2016_gravity_drainage/SOURCES/source_3.20/ice_gas_solu.f @ 24

Download in other formats: