Line data Source code
1 : ! Include shortname defintions, so that the F77 code does not have to be modified to
2 : ! reference the CARMA structure.
3 : #include "carma_globaer.h"
4 :
5 : module sulfate_utils
6 : use carma_precision_mod
7 : use carma_enums_mod
8 : use carma_constants_mod
9 : use carma_types_mod
10 : use carmastate_mod
11 : use carma_mod
12 :
13 : implicit none
14 :
15 : ! Declare the public methods.
16 : public wtpct_tabaz
17 : public sulfate_density
18 : public sulfate_surf_tens
19 :
20 : real(kind=f), public:: dnwtp(46), dnc0(46), dnc1(46)
21 :
22 : data dnwtp / 0._f, 1._f, 5._f, 10._f, 20._f, 25._f, 30._f, 35._f, 40._f, &
23 : 41._f, 45._f, 50._f, 53._f, 55._f, 56._f, 60._f, 65._f, 66._f, 70._f, &
24 : 72._f, 73._f, 74._f, 75._f, 76._f, 78._f, 79._f, 80._f, 81._f, 82._f, &
25 : 83._f, 84._f, 85._f, 86._f, 87._f, 88._f, 89._f, 90._f, 91._f, 92._f, &
26 : 93._f, 94._f, 95._f, 96._f, 97._f, 98._f, 100._f /
27 :
28 : data dnc0 / 1._f, 1.13185_f, 1.17171_f, 1.22164_f, 1.3219_f, 1.37209_f, &
29 : 1.42185_f, 1.4705_f, 1.51767_f, 1.52731_f, 1.56584_f, 1.61834_f, 1.65191_f, &
30 : 1.6752_f, 1.68708_f, 1.7356_f, 1.7997_f, 1.81271_f, 1.86696_f, 1.89491_f, &
31 : 1.9092_f, 1.92395_f, 1.93904_f, 1.95438_f, 1.98574_f, 2.00151_f, 2.01703_f, &
32 : 2.03234_f, 2.04716_f, 2.06082_f, 2.07363_f, 2.08461_f, 2.09386_f, 2.10143_f,&
33 : 2.10764_f, 2.11283_f, 2.11671_f, 2.11938_f, 2.12125_f, 2.1219_f, 2.12723_f, &
34 : 2.12654_f, 2.12621_f, 2.12561_f, 2.12494_f, 2.12093_f /
35 :
36 : data dnc1 / 0._f, -0.000435022_f, -0.000479481_f, -0.000531558_f, -0.000622448_f,&
37 : -0.000660866_f, -0.000693492_f, -0.000718251_f, -0.000732869_f, -0.000735755_f, &
38 : -0.000744294_f, -0.000761493_f, -0.000774238_f, -0.00078392_f, -0.000788939_f, &
39 : -0.00080946_f, -0.000839848_f, -0.000845825_f, -0.000874337_f, -0.000890074_f, &
40 : -0.00089873_f, -0.000908778_f, -0.000920012_f, -0.000932184_f, -0.000959514_f, &
41 : -0.000974043_f, -0.000988264_f, -0.00100258_f, -0.00101634_f, -0.00102762_f, &
42 : -0.00103757_f, -0.00104337_f, -0.00104563_f, -0.00104458_f, -0.00104144_f, &
43 : -0.00103719_f, -0.00103089_f, -0.00102262_f, -0.00101355_f, -0.00100249_f, &
44 : -0.00100934_f, -0.000998299_f, -0.000990961_f, -0.000985845_f, -0.000984529_f, &
45 : -0.000989315_f /
46 : contains
47 :
48 : !! This function calculates the weight % H2SO4 composition of
49 : !! sulfate aerosol, using Tabazadeh et. al. (GRL, 1931, 1997).
50 : !! Rated for T=185-260K, activity=0.01-1.0
51 : !!
52 : !! Argument list input:
53 : !! temp = temperature (K)
54 : !! h2o_mass = water vapor mass concentration (g/cm3)
55 : !! h2o_vp = water eq. vaper pressure (dynes/cm2)
56 : !!
57 : !! Output:
58 : !! wtpct_tabaz = weight % H2SO4 in H2O/H2SO4 particle (0-100)
59 : !!
60 : !! Include global constants and variables (BK=Boltzman constant,
61 : !! AVG=Avogadro's constant)
62 : !!
63 : !! @author Jason English
64 : !! @ version Apr-2010
65 9414463709 : function wtpct_tabaz(carma, temp, h2o_mass, h2o_vp, rc)
66 :
67 : real(kind=f) :: wtpct_tabaz
68 : type(carma_type), intent(in) :: carma !! the carma object
69 : real(kind=f), intent(in) :: temp !! temperature [K]
70 : real(kind=f), intent(in) :: h2o_mass !! water vapor mass concentration (g/cm3)
71 : real(kind=f), intent(in) :: h2o_vp !! water eq. vaper pressure (dynes/cm2)
72 : integer, intent(inout) :: rc !! return code, negative indicates failure
73 :
74 : ! Declare variables for this routine only
75 : real(kind=f) :: atab1,btab1,ctab1,dtab1,atab2,btab2,ctab2,dtab2
76 : real(kind=f) :: h2o_num, p_h2o, vp_h2o
77 : real(kind=f) :: contl, conth, contt, conwtp
78 : real(kind=f) :: activ
79 :
80 : ! Get number density of water (/cm3) from mass concentration (g/cm3)
81 9414463709 : h2o_num=h2o_mass*AVG/gwtmol(1)
82 :
83 : ! Get partial pressure of water (dynes/cm2) from concentration (/cm3)
84 : ! Ideal gas law: P=nkT
85 9414463709 : p_h2o=h2o_num*bk*temp
86 :
87 : ! Convert from dynes/cm2 to mb (hPa)
88 9414463709 : p_h2o=p_h2o/1000.0_f ! partial pressure
89 9414463709 : vp_h2o=h2o_vp/1000.0_f ! eq. vp
90 :
91 : ! Activity = water pp in mb / water eq. vp over pure water in mb
92 9414463709 : activ = p_h2o/vp_h2o
93 :
94 9414463709 : if (activ.lt.0.05_f) then
95 4950544760 : activ = max(activ,1.e-32_f) ! restrict minimum activity
96 4950544760 : atab1 = 12.37208932_f
97 4950544760 : btab1 = -0.16125516114_f
98 4950544760 : ctab1 = -30.490657554_f
99 4950544760 : dtab1 = -2.1133114241_f
100 4950544760 : atab2 = 13.455394705_f
101 4950544760 : btab2 = -0.1921312255_f
102 4950544760 : ctab2 = -34.285174607_f
103 4950544760 : dtab2 = -1.7620073078_f
104 4463918949 : elseif (activ.ge.0.05_f.and.activ.le.0.85_f) then
105 : atab1 = 11.820654354_f
106 : btab1 = -0.20786404244_f
107 : ctab1 = -4.807306373_f
108 : dtab1 = -5.1727540348_f
109 : atab2 = 12.891938068_f
110 : btab2 = -0.23233847708_f
111 : ctab2 = -6.4261237757_f
112 : dtab2 = -4.9005471319_f
113 593230431 : elseif (activ.gt.0.85_f) then
114 593230431 : activ = min(activ,1._f) ! restrict maximum activity
115 593230431 : atab1 = -180.06541028_f
116 593230431 : btab1 = -0.38601102592_f
117 593230431 : ctab1 = -93.317846778_f
118 593230431 : dtab1 = 273.88132245_f
119 593230431 : atab2 = -176.95814097_f
120 593230431 : btab2 = -0.36257048154_f
121 593230431 : ctab2 = -90.469744201_f
122 593230431 : dtab2 = 267.45509988_f
123 : else
124 0 : if (do_print) write(LUNOPRT,*) 'invalid activity: activity,pp,vp=',activ, p_h2o
125 0 : rc = RC_ERROR
126 0 : return
127 : endif
128 :
129 9414463709 : contl = atab1*(activ**btab1)+ctab1*activ+dtab1
130 9414463709 : conth = atab2*(activ**btab2)+ctab2*activ+dtab2
131 :
132 9414463709 : contt = contl + (conth-contl) * ((temp -190._f)/70._f)
133 9414463709 : conwtp = (contt*98._f) + 1000._f
134 :
135 9414463709 : wtpct_tabaz = (100._f*contt*98._f)/conwtp
136 9414463709 : wtpct_tabaz = min(max(wtpct_tabaz,1._f),100._f) ! restrict between 1 and 100 %
137 :
138 9414463709 : return
139 : end function wtpct_tabaz
140 :
141 : !! Calculates specific gravity (g/cm3) of sulfate of
142 : !! different compositions as a linear function of temperature,
143 : !! based of measurements of H2SO4/H2O solution densities made
144 : !! at 0 to 100C tabulated in the International Critical Tables
145 : !! (Washburn, ed., NRC, 1928). Measurements have confirmed that
146 : !! this data may be linearly extrapolated to stratospheric
147 : !! temperatures (180-380K) with excellent accuracy
148 : !! (Beyer, Ravishankara, & Lovejoy, JGR, 1996).
149 : !!
150 : !! Argument list input:
151 : !! wtp = aerosol composition in weight % H2SO4 (0-100)
152 : !! temp = temperature in Kelvin
153 : !!
154 : !! Output:
155 : !! sulfate_density (g/cm3) [function name]
156 : !!
157 : !! This function requires setup_sulfate_density to be run
158 : !! first to read in the density coefficients DNC0 and DNC1
159 : !! and the tabulated weight percents DNWTP.
160 : !!
161 : !! @author Mike Mills
162 : !! @version Mar-2013
163 9150481871 : function sulfate_density(carma, wtp, temp, rc)
164 :
165 : !! Include global constants and variables
166 :
167 : real(kind=f) :: sulfate_density
168 : type(carma_type), intent(in) :: carma !! the carma object
169 : real(kind=f), intent(in) :: wtp !! weight percent
170 : real(kind=f), intent(in) :: temp !! temperature
171 : integer, intent(inout) :: rc !! return code, negative indicates failure
172 :
173 : ! Local declarations
174 : integer :: i
175 : real(kind=f) :: den1, den2
176 : real(kind=f) :: frac, temp_loc
177 :
178 9150481871 : if (wtp .lt. 0.0_f .or. wtp .gt. 100.0_f) then
179 0 : if (do_print) write(LUNOPRT,*)'sulfate_density: Illegal value for wtp:',wtp
180 0 : rc = RC_ERROR
181 0 : return
182 : endif
183 :
184 : ! limit temperature to bounds of extrapolation
185 9150481871 : temp_loc=min(temp, 380.0_f)
186 9150481871 : temp_loc=max(temp_loc, 180.0_f)
187 :
188 9150481871 : i=1
189 :
190 >20349*10^7 : do while (wtp .gt. dnwtp(i))
191 >19434*10^7 : i=i+1
192 : end do
193 :
194 9150481871 : den2=dnc0(i)+dnc1(i)*temp_loc
195 :
196 9150481871 : if (i.eq.1 .or. wtp.eq.dnwtp(i)) then
197 9150481871 : sulfate_density=den2
198 : return
199 : endif
200 :
201 9106589908 : den1=dnc0(i-1)+dnc1(i-1)*temp_loc
202 9106589908 : frac=(dnwtp(i)-wtp)/(dnwtp(i)-dnwtp(i-1))
203 9106589908 : sulfate_density=den1*frac+den2*(1.0_f-frac)
204 :
205 9106589908 : return
206 : end function sulfate_density
207 :
208 : !! Calculates surface tension (erg/cm2 = dyne/cm) of sulfate of
209 : !! different compositions as a linear function of temperature,
210 : !! as described in Mills (Ph.D. Thesis, 1996), derived from
211 : !! the measurements of Sabinina and Terpugow (1935).
212 : !!
213 : !! Argument list input:
214 : !! WTP = aerosol composition in weight % H2SO4 (0-100)
215 : !! TEMP = temperature in Kelvin
216 : !!
217 : !! Output:
218 : !! sulfate_surf_tens (erg/cm2) [function name]
219 : !!
220 : !! This function requires setup_sulfate_density to be run
221 : !! first to read in the density coefficients DNC0 and DNC1
222 : !! and the tabulated weight percents DNWTP.
223 : !!
224 : !! @author Mike Mills
225 : !! @version Mar-2013
226 244727514 : function sulfate_surf_tens(carma, wtp, temp, rc)
227 :
228 : real(kind=f) :: sulfate_surf_tens
229 : type(carma_type), intent(in) :: carma !! the carma object
230 : real(kind=f), intent(in) :: wtp !! weight percent
231 : real(kind=f), intent(in) :: temp !! temperature
232 : integer, intent(inout) :: rc !! return code, negative indicates failure
233 :
234 : ! Local declarations
235 : integer :: i
236 : real(kind=f) :: sig1, sig2
237 : real(kind=f) :: frac, temp_loc
238 : real(kind=f) :: stwtp(15), stc0(15), stc1(15)
239 :
240 : data stwtp/0._f, 23.8141_f, 38.0279_f, 40.6856_f, 45.335_f, 52.9305_f, 56.2735_f, &
241 : & 59.8557_f, 66.2364_f, 73.103_f, 79.432_f, 85.9195_f, 91.7444_f, 97.6687_f, 100._f/
242 :
243 : data stc0/117.564_f, 103.303_f, 101.796_f, 100.42_f, 98.4993_f, 91.8866_f, &
244 : & 88.3033_f, 86.5546_f, 84.471_f, 81.2939_f, 79.3556_f, 75.608_f, 70.0777_f, &
245 : & 63.7412_f, 61.4591_f /
246 :
247 : data stc1/-0.153641_f, -0.0982007_f, -0.0872379_f, -0.0818509_f, &
248 : & -0.0746702_f, -0.0522399_f, -0.0407773_f, -0.0357946_f, -0.0317062_f, &
249 : & -0.025825_f, -0.0267212_f, -0.0269204_f, -0.0276187_f, -0.0302094_f, &
250 : & -0.0303081_f /
251 :
252 : ! limit temperature to reasonable bounds of extrapolation
253 244727514 : temp_loc=min(temp, 380.0_f)
254 244727514 : temp_loc=max(temp_loc, 180.0_f)
255 :
256 244727514 : if (wtp .lt. 0.0_f .OR. wtp .gt. 100.0_f) then
257 0 : if (do_print) write(LUNOPRT,*)'sulfate_surf_tens: Illegal value for wtp:',wtp
258 0 : if (do_print) write(LUNOPRT,*)'sulfate_surf_tens: temp=',temp
259 0 : rc = RC_ERROR
260 0 : return
261 : endif
262 :
263 : i=1
264 :
265 2019407037 : do while (wtp.gt.stwtp(i))
266 1774679523 : i=i+1
267 : end do
268 :
269 244727514 : sig2=stc0(i)+stc1(i)*temp_loc
270 :
271 244727514 : if (i.eq.1 .or. wtp.eq.stwtp(i)) then
272 244727514 : sulfate_surf_tens=sig2
273 : return
274 : end if
275 :
276 244727514 : sig1=stc0(i-1)+stc1(i-1)*temp_loc
277 244727514 : frac=(stwtp(i)-wtp)/(stwtp(i)-stwtp(i-1))
278 244727514 : sulfate_surf_tens=sig1*frac+sig2*(1.0_f-frac)
279 :
280 244727514 : return
281 : end function sulfate_surf_tens
282 :
283 : end module sulfate_utils
|
