-
Notifications
You must be signed in to change notification settings - Fork 0
/
Aero_water_mod.f90
executable file
·374 lines (297 loc) · 11.5 KB
/
Aero_water_mod.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
! <Aero_water_mod.f90 - A component of the EMEP MSC-W Unified Eulerian
! Chemical transport Model>
!*****************************************************************************!
!*
!* Copyright (C) 2007-2011 met.no
!*
!* Contact information:
!* Norwegian Meteorological Institute
!* Box 43 Blindern
!* 0313 OSLO
!* NORWAY
!* email: [email protected]
!* http://www.emep.int
!*
!* This program is free software: you can redistribute it and/or modify
!* it under the terms of the GNU General Public License as published by
!* the Free Software Foundation, either version 3 of the License, or
!* (at your option) any later version.
!*
!* This program is distributed in the hope that it will be useful,
!* but WITHOUT ANY WARRANTY; without even the implied warranty of
!* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!* GNU General Public License for more details.
!*
!* You should have received a copy of the GNU General Public License
!* along with this program. If not, see <http://www.gnu.org/licenses/>.
!*****************************************************************************!
module Aero_water_mod
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
! DESCRIPTION
!
! Purpose : Calculates aerosols' liquid water content
!
! Subroutine: Awater
! Input: - relative humidity: relh
! - number of micromoles/(m^3 of air) for sulfate,
! ammonium, and nitrate: mso4, mnh4, mno3
! Output: - water amount in micrograms/(m^3 of air): wh2o
!
! Author : Dr. Francis S. Binkowski, 4/8/96
! : modified slightly for EMEP model
!
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
implicit none
private
! subroutines:
public :: Awater
! functions
private :: poly4, poly6
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
! coefficients for polynomials (function poly4) to be defined
! at start of routine:
! Define saved variables:
real, private,parameter, dimension(4) :: & !(for x = 0, 1, 1.5 and 2)
C0 = (/ 0.798079, -1.574367, 2.536686, -1.735297 /),&
C1 = (/ 0.9995178, -0.7952896, 0.99683673, -1.143874 /),&
C15= (/ 1.697092, -4.045936, 5.833688, -3.463783 /)!,&
!C2 = (/ 2.085067, -6.024139, 8.967967, -5.002934 /)
real, private,parameter, dimension(6) :: &
KNO3 = (/ 0.2906, 6.83665, -26.9093, 46.6983, -38.803, 11.8837/),&
KSO4 = (/ 2.27515, -11.147, 36.3369, -64.2134, 56.8341 ,-20.0953/)
! Set molecular weights:
real, private, parameter :: &
MWSO4 = 96.0 &
,MWNH4 = 18.0 &
,MWNO3 = 62.0 &
,MW2 = MWSO4 + 2.0 * MWNH4 &
,MWANO3 = MWNO3 + MWNH4
contains
subroutine Awater(relh,mso4,mnh4,mno3,wh2o)
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
! This routine uses polynomials rather than tables, and uses empirical
! polynomials for the mass fraction of solute (mfs) as a function of
! water activity where:
!
! mfs = ms / ( ms + mw)
! ms - the mass of solute
! mw - the mass of water
!
! Define y = mw / ms
!
! then mfs = 1 / (1 + y)
!
! y can then be obtained from the values of mfs as
!
! y = (1 - mfs) / mfs
!
! The aerosol is assumed to be in a metastable state if the relative
! humidity (rh) is is below the rh of deliquescence, but above the
! rh of crystallization.
!
! The Zdanovskii-Stokes-Robinson relation ('ZSR') is used for sulfates
! with x (molar ratio of ammonium to sulfate) in the range 0 <= x <= 2,
! subdivided into four sections:
!
! section 1: 0 <= x < 1
! section 2: 1 <= x < 1.5
! section 3: 1.5 <= x < 2.0
! section 4: 2 <= x
!
! In sections 1 through 3 only the sulfates can affect the amount
! of water on the particles. In section 4 we have fully neutralized
! sulfate, and extra ammonium which allows more nitrate to be present.
! Thus, the ammount of water is calculated using ZSR for ammonium
! sulfate and ammonium nitrate. Crystallization is assumed to occur
! in sections 2, 3, and 4 (see detailed discussion below).
!
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!
! Description of variables used in the routine:
!
!
! mso4, mnh4, and mno3 : number of micromoles/(cubic meter of air) for
! sulfate, ammonium, and nitrate, respectively
! relh : relative humidity (%)
! wh2o : returned water amount in micrograms /(cubic meter of air)
! x : molar ratio of ammonium to sulfate
! y0, y1, y15, y2 : water contents in mass of water/mass of solute
! for pure aqueous solutions with x equal to 0, 1,
! 1.5, and 2, respectively
! y3 : the value of the mass ratio of water to solute for a pure
! ammonium nitrate solution.
!
!
! C1, C15, C2: (defined at start of routine!)
! The polynomials use data for relh as a function of mfs from Tang
! and Munkelwitz, JGR. 99: 18801-18808, 1994.
! The polynomials were fit to Tang's values of water activity as a
! function of mfs.
! C1, C15, C2 are the coefficients of polynomials fit to Tang and
! Munkelwitz data giving mfs as a function of water activity.
!
! C0: (defined at start of routine!)
! fit to data from
! Nair and Vohra J. Aerosol Sci., 6: 265-271, 1975
! Giaque et al. J. Am. Chem. Soc., 82: 62-70, 1960
! Zeleznik J. Phys. Chem. ref. data, 20: 157-1200
!
! KNO3, KSO4: (defined at start of routine!)
! The polynomials for ammonium nitrate and ammonium sulfate are from:
! Chan et al.1992, atmospheric environment (26a): 1661-1673.
!
!
! tso4, tnh4, tno3: mole concentrations used in the calculations
! ( tso4 = max(mso4,0.), tnh4 = max(mnh4,0.), tno3 = max(mno3,0.) )
!
! aw : relative humidity used in the calculations (.01<aw<.95 required)
! awc : relative humidity along the crystallization curve
! u : help variable for rh used for interpolation (u=40%)
! x : molar ratio (x = tnh4 / tso4 for non-zero sulfate, and =10 else)
!
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!-- in
real, intent(in) :: relh, mso4, mnh4, mno3
!-- out
real, intent(out) :: wh2o
!-- local
real :: tso4, tnh4, tno3
real :: x, awc, aw, u &
, mfs0, mfs1, mfs15 &
, mfsso4, mfsno3 &
, y, y0, y1, y15, y2, y3 &
, y40, y140, y1540, yc
! Check range of per cent relative humidity:
! aw - water activity = fractional relative humidity
aw = relh
! Set aw between 0.01 and 0.95:
aw = max(0.01,aw)
aw = min(aw,0.95)
tso4 = max(mso4,0.)
tnh4 = max(mnh4,0.)
tno3 = max(mno3,0.)
! Initialise molar ratio:
x = 0.
! If there is non-zero sulfate calculate the molar ratio:
if (tso4 > 0.) then
x = tnh4 / tso4
else
! ... otherwise check for non-zero nitrate and ammonium
if (tno3 > 0. .and. tnh4 > 0.) x = 10.
end if
! Begin screen on x for calculating wh2o:
if ( x < 1. ) then
mfs0 = poly4(C0,aw)
mfs1 = poly4(C1,aw)
y0 = ( 1. - mfs0 ) / mfs0
y1 = ( 1. - mfs1 ) / mfs1
y = ( 1. - x ) * y0 + x * y1
else if ( x < 1.5) then
if ( aw >= 0.40 ) then
mfs1 = poly4(C1,aw)
mfs15 = poly4(C15,aw)
y1 = (1. - mfs1 ) / mfs1
y15 = (1. - mfs15) / mfs15
y = 2. * ( y1 * (1.5 - x) + y15 *( x - 1.) )
else
! Setup for crystalization:
! Crystallization is done as follows:
!
! for 1.5 <= x : crystallization is assumed to occur at rh = 0.4
! for x <= 1.0 : crystallization is assumed to occur at an rh < 0.01
!
! and since the code does not allow rh < 0.01, crystallization is
! assumed not to occur in this range.
!
! for 1.0 <= x <= 1.5 : the crystallization curve is a straight line
! from a value of y15 at rh = 0.4 to a value of
! zero at y1. From point b to point a in the
! diagram the algorithm does a double inter-
! polation to calculate the amount of water.
!
! y1(0.40) y15(0.40)
! + + point b
!
!
!
! +---------------------+
! x=1 x=1.5
! point a
!
awc = 0.80 * (x - 1.0) ! rh along the crystallization curve.
y = 0.0
u=0.40
if ( aw >= awc ) then
! Interpolate using crystalization curve:
mfs1 = poly4(C1,u)
mfs15 = poly4(C15,u)
y140 = (1.0 - mfs1 ) / mfs1
y1540 = (1.0 - mfs15) / mfs15
y40 = 2.0 * ( y140 * (1.5 - x) + y1540 *( x - 1.0) )
yc = 2.0 * y1540 * (x -1.0) ! y along crystallization curve
y = y40 - (y40 - yc) * (u - aw) / (u - awc)
end if ! end of "if ( aw >= awc ) then"
end if ! end of "if ( aw >= 0.40 ) then"
else if ( x < 1.9999) then
y= 0.0
if (aw >= 0.40) then
mfs15 = poly4(C15,aw)
!mfs2 = poly4(C2,aw)
!y2 = (1.0 - mfs2) / mfs2
MFSSO4 = POLY6( KSO4, AW ) ! PW-18/12/2014 after Changed 05/30/2000 by FSB
Y2 = ( 1.0d0 - MFSSO4 ) / MFSSO4
y15 = (1.0 - mfs15) / mfs15
y = 2.0 * (y15 * (2.0 - x) + y2 * (x - 1.5) )
end if ! end of check for crystallization
else
! i.e. x >= 1.9999
!
! Regime where ammonium sulfate and ammonium nitrate are in solution!
!
! Following cf&s for both ammonium sulfate and ammonium nitrate
! Check for crystallization here. Their data indicate a 40% value
! is appropriate.
y2 = 0.0
y3 = 0.0
if (aw >= 0.40) then
mfsso4 = poly6(KSO4,aw)
mfsno3 = poly6(KNO3,aw)
y2 = (1.0 - mfsso4) / mfsso4
y3 = (1.0 - mfsno3) / mfsno3
end if
end if ! end of 'if ( x < 1. ) then'
! Now set up output of wh2o
!
! wh2o units are micrograms(liquid water) / (cubic meter of air)
if ( x < 1.9999) then
wh2o = y * (tso4 * MWSO4 + tnh4 * MWNH4 )
else
! This is the case when all the sulfate is ammonium sulfate
! and the excess ammonium forms ammonum nitrate
wh2o = y2 * tso4 * MW2 + y3 * tno3 * MWANO3
end if
end subroutine Awater
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
real function poly4 (a,x)
! Calculates the polynomial based on 4 coefficients a(1:4):
!-- arguments
real, dimension(4), intent(in) :: a
real, intent(in) :: x
poly4 = a(1) + x * ( a(2) + x * ( a(3) + x * ( a(4) ) ) )
end function poly4
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
real function poly6(a,x)
! Calculates the polynomial based on 6 coefficients a(1:6):
!-- arguments
real, dimension(6), intent(in) :: a
real, intent(in) :: x
poly6 = a(1) + x * ( a(2) + x * ( a(3) + x * ( a(4) + &
x * ( a(5) + x * (a(6) ) ) ) ) )
end function poly6
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
end module Aero_water_mod