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 : !! This routine evaluate particle loss rates due to condensational
6 : !! growth and evaporation for all condensing gases.
7 : !!
8 : !! The loss rates for each group are <growlg> and <evaplg>.
9 : !!
10 : !! Units are [s^-1].
11 : !!
12 : !! @author Andy Ackerman
13 : !! @version Dec-1995
14 1418703438 : subroutine growevapl(carma, cstate, iz, rc)
15 :
16 : ! types
17 : use carma_precision_mod
18 : use carma_enums_mod
19 : use carma_constants_mod
20 : use carma_types_mod
21 : use carmastate_mod
22 : use carma_mod
23 :
24 : implicit none
25 :
26 : type(carma_type), intent(in) :: carma !! the carma object
27 : type(carmastate_type), intent(inout) :: cstate !! the carma state object
28 : integer, intent(in) :: iz !! z index
29 : integer, intent(inout) :: rc !! return code, negative indicates failure
30 :
31 : ! Local declarations
32 : integer :: igroup
33 : integer :: iepart
34 : integer :: igas
35 : integer :: ibin
36 : integer :: isol
37 : integer :: nother
38 : integer :: ieoth_rel
39 : integer :: ieoth_abs
40 : integer :: jother
41 : real(kind=f) :: argsol
42 : real(kind=f) :: othermtot
43 : real(kind=f) :: condm
44 : real(kind=f) :: akas
45 : real(kind=f) :: expon
46 : real(kind=f) :: g0
47 : real(kind=f) :: g1
48 : real(kind=f) :: g2
49 : real(kind=f) :: ss
50 : real(kind=f) :: pvap
51 : real(kind=f) :: dpc
52 : real(kind=f) :: dpc1
53 : real(kind=f) :: dpcm1
54 : real(kind=f) :: rat1
55 : real(kind=f) :: rat2
56 : real(kind=f) :: rat3
57 : real(kind=f) :: rat4
58 : real(kind=f) :: ratt1
59 : real(kind=f) :: ratt2
60 : real(kind=f) :: ratt3
61 : real(kind=f) :: den1
62 : real(kind=f) :: test1
63 : real(kind=f) :: test2
64 : real(kind=f) :: x
65 : integer :: ieother(NELEM)
66 : real(kind=f) :: otherm(NELEM)
67 2837406876 : real(kind=f) :: dela(NBIN)
68 2837406876 : real(kind=f) :: delma(NBIN)
69 2837406876 : real(kind=f) :: aju(NBIN)
70 2837406876 : real(kind=f) :: ar(NBIN)
71 2837406876 : real(kind=f) :: al(NBIN)
72 2837406876 : real(kind=f) :: a6(NBIN)
73 2837406876 : real(kind=f) :: dmdt(NBIN)
74 : real(kind=f) :: growlg_max
75 :
76 :
77 4256110314 : do igroup = 1,NGROUP
78 :
79 : ! element of particle number concentration
80 2837406876 : iepart = ienconc(igroup)
81 :
82 : ! condensing gas
83 2837406876 : igas = igrowgas(iepart)
84 :
85 4256110314 : if (igas .ne. 0) then
86 : ! Only valid for condensing liquid water and sulfric acid currently.
87 2837406876 : if ((igas /= igash2o) .and. (igas .ne. igash2so4)) then
88 0 : if (do_print) write(LUNOPRT,*) 'growevapl::ERROR - Invalid gas (', igas, ').'
89 0 : rc = -1
90 0 : return
91 : endif
92 :
93 : ! Treat condensation of gas <igas> to/from particle group <igroup>.
94 : !
95 : ! Bypass calculation if few particles are present
96 2837406876 : if( pconmax(iz,igroup) .gt. FEW_PC )then
97 35478933300 : do ibin = 1,NBIN-1
98 :
99 : ! Determine the growth rate (dmdt). This calculation may take into account
100 : ! radiative effects on the particle which can affect the growth rates.
101 35478933300 : call pheat(carma, cstate, iz, igroup, iepart, ibin, igas, dmdt(ibin), rc)
102 :
103 : enddo ! ibin = 1,NBIN-1
104 :
105 : ! Now calculate condensation/evaporation production and loss rates.
106 : ! Use Piecewise Polynomial Method [Colela and Woodard, J. Comp. Phys.,
107 : ! 54, 174-201, 1984]
108 : !
109 : ! First, use cubic fits to estimate concentration values at bin
110 : ! boundaries
111 33704986635 : do ibin = 2,NBIN-1
112 :
113 31931039970 : dpc = pc(iz,ibin,iepart) / dm(ibin,igroup)
114 31931039970 : dpc1 = pc(iz,ibin+1,iepart) / dm(ibin+1,igroup)
115 31931039970 : dpcm1 = pc(iz,ibin-1,iepart) / dm(ibin-1,igroup)
116 31931039970 : ratt1 = pratt(1,ibin,igroup)
117 31931039970 : ratt2 = pratt(2,ibin,igroup)
118 31931039970 : ratt3 = pratt(3,ibin,igroup)
119 31931039970 : dela(ibin) = ratt1 * ( ratt2*(dpc1-dpc) + ratt3*(dpc-dpcm1) )
120 31931039970 : delma(ibin) = 0._f
121 :
122 31931039970 : if( (dpc1-dpc)*(dpc-dpcm1) .gt. 0._f ) &
123 : delma(ibin) = min( abs(dela(ibin)), 2._f*abs(dpc-dpc1), &
124 32899409334 : 2._f*abs(dpc-dpcm1) ) * sign(1._f, dela(ibin))
125 :
126 : enddo ! ibin = 2,NBIN-2
127 :
128 31931039970 : do ibin = 2,NBIN-2
129 :
130 30157093305 : dpc = pc(iz,ibin,iepart) / dm(ibin,igroup)
131 30157093305 : dpc1 = pc(iz,ibin+1,iepart) / dm(ibin+1,igroup)
132 30157093305 : dpcm1 = pc(iz,ibin-1,iepart) / dm(ibin-1,igroup)
133 30157093305 : rat1 = prat(1,ibin,igroup)
134 30157093305 : rat2 = prat(2,ibin,igroup)
135 30157093305 : rat3 = prat(3,ibin,igroup)
136 30157093305 : rat4 = prat(4,ibin,igroup)
137 30157093305 : den1 = pden1(ibin,igroup)
138 :
139 : ! <aju(ibin)> is the estimate for concentration (dn/dm) at bin
140 : ! boundary <ibin>+1/2.
141 30157093305 : aju(ibin) = dpc + rat1*(dpc1-dpc) + 1._f/den1 * &
142 : ( rat2*(rat3-rat4)*(dpc1-dpc) - &
143 30157093305 : dm(ibin,igroup)*rat3*delma(ibin+1) + &
144 31931039970 : dm(ibin+1,igroup)*rat4*delma(ibin) )
145 : enddo ! ibin = 2,NBIN-2
146 :
147 : ! Now construct polynomial functions in each bin
148 30157093305 : do ibin = 3,NBIN-2
149 28383146640 : al(ibin) = aju(ibin-1)
150 30157093305 : ar(ibin) = aju(ibin)
151 : enddo
152 :
153 : ! Use linear functions in first two and last two bins
154 1773946665 : if( NBIN .gt. 1 )then
155 1773946665 : ibin = NBIN
156 :
157 1773946665 : ar(2) = aju(2)
158 1773946665 : al(2) = pc(iz,1,iepart)/dm(1,igroup) + &
159 0 : palr(1,igroup) * &
160 1773946665 : (pc(iz,2,iepart)/dm(2,igroup)- &
161 3547893330 : pc(iz,1,iepart)/dm(1,igroup))
162 1773946665 : ar(1) = al(2)
163 0 : al(1) = pc(iz,1,iepart)/dm(1,igroup) + &
164 1773946665 : palr(2,igroup) * &
165 1773946665 : (pc(iz,2,iepart)/dm(2,igroup)- &
166 1773946665 : pc(iz,1,iepart)/dm(1,igroup))
167 :
168 1773946665 : al(ibin-1) = aju(ibin-2)
169 0 : ar(ibin-1) = pc(iz,ibin-1,iepart)/dm(ibin-1,igroup) + &
170 1773946665 : palr(3,igroup) * &
171 1773946665 : (pc(iz,ibin,iepart)/dm(ibin,igroup)- &
172 1773946665 : pc(iz,ibin-1,iepart)/dm(ibin-1,igroup))
173 1773946665 : al(ibin) = ar(ibin-1)
174 0 : ar(ibin) = pc(iz,ibin-1,iepart)/dm(ibin-1,igroup) + &
175 1773946665 : palr(4,igroup) * &
176 1773946665 : (pc(iz,ibin,iepart)/dm(ibin,igroup)- &
177 3547893330 : pc(iz,ibin-1,iepart)/dm(ibin-1,igroup))
178 : endif
179 :
180 : ! Next, ensure that polynomial functions do not deviate beyond the
181 : ! range [<al(ibin)>,<ar(ibin)>]
182 37252879965 : do ibin = 1,NBIN
183 :
184 35478933300 : dpc = pc(iz,ibin,iepart) / dm(ibin,igroup)
185 :
186 35478933300 : if( (ar(ibin)-dpc)*(dpc-al(ibin)) .le. 0._f )then
187 805577301 : al(ibin) = dpc
188 805577301 : ar(ibin) = dpc
189 : endif
190 :
191 35478933300 : test1 = (ar(ibin)-al(ibin))*(dpc - 0.5_f*(al(ibin)+ar(ibin)))
192 35478933300 : test2 = 1._f/6._f*(ar(ibin)-al(ibin))**2
193 :
194 37252879965 : if( test1 .gt. test2 )then
195 15610577933 : al(ibin) = 3._f*dpc - 2._f*ar(ibin)
196 19868355367 : elseif( test1 .lt. -test2 )then
197 372143795 : ar(ibin) = 3._f*dpc - 2._f*al(ibin)
198 : endif
199 : enddo
200 :
201 : ! Lastly, calculate fluxes across each bin boundary.
202 : !
203 : ! Use upwind advection when courant number > 1.
204 37252879965 : do ibin = 1,NBIN
205 35478933300 : dpc = pc(iz,ibin,iepart) / dm(ibin,igroup)
206 35478933300 : dela(ibin) = ar(ibin) - al(ibin)
207 37252879965 : a6(ibin) = 6._f * ( dpc - 0.5_f*(ar(ibin)+al(ibin)) )
208 : enddo
209 :
210 35478933300 : do ibin = 1,NBIN-1
211 :
212 33704986635 : if( dmdt(ibin) .gt. 0._f .and. &
213 1773946665 : pc(iz,ibin,iepart) .gt. SMALL_PC )then
214 :
215 32817939583 : x = dmdt(ibin)*dtime/dm(ibin,igroup)
216 :
217 32817939583 : if( x .lt. 1._f )then
218 0 : growlg(ibin,igroup) = dmdt(ibin)/pc(iz,ibin,iepart) &
219 : * ( ar(ibin) - 0.5_f*dela(ibin)*x + &
220 32816383287 : (x/2._f - x**2/3._f)*a6(ibin) )
221 : else
222 1556296 : growlg(ibin,igroup) = dmdt(ibin) / dm(ibin,igroup)
223 : endif
224 :
225 887047052 : elseif( dmdt(ibin) .lt. 0._f .and. &
226 : pc(iz,ibin+1,iepart) .gt. SMALL_PC )then
227 :
228 609380597 : x = -dmdt(ibin)*dtime/dm(ibin+1,igroup)
229 :
230 609380597 : if( x .lt. 1._f )then
231 0 : evaplg(ibin+1,igroup) = -dmdt(ibin)/ &
232 : pc(iz,ibin+1,iepart) &
233 605688161 : * ( al(ibin+1) + 0.5_f*dela(ibin+1)*x + &
234 1211376322 : (x/2._f - (x**2)/3._f)*a6(ibin+1) )
235 : else
236 3692436 : evaplg(ibin+1,igroup) = -dmdt(ibin) / dm(ibin+1,igroup)
237 : endif
238 :
239 : ! Boundary conditions: for evaporation out of first bin (with cores),
240 : ! use evaporation rate from second bin.
241 : ! if( ibin .eq. 1 .and. ncore(igroup) .gt. 0 )then
242 609380597 : if( ibin .eq. 1)then
243 213694496 : evaplg(1,igroup) = -dmdt(1) / dm(1,igroup)
244 : endif
245 : endif
246 :
247 : enddo ! ibin = 1,NBIN-1
248 : endif ! (pconmax .gt. FEW_PC)
249 : endif ! (igas = igrowgas(ielem)) .ne. 0
250 : enddo ! igroup = 1,NGROUP
251 :
252 :
253 : ! Return to caller with particle loss rates for growth and evaporation
254 : ! evaluated.
255 : return
256 1418703438 : end
|
