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 particle heating.
6 : !!
7 : !! The net energy absorbed by each particle is calculatated as <qrad>, and
8 : !! this heating rate is included in the caclulation of <dmdt> in growevapl. The
9 : !! particle temperature perturbation realtive to atmospheric temperature <dtpart>
10 : !! and the radiative heating of the atmosphere by particles <partheat>
11 : !! are also calculated.
12 : !!
13 : !! This algorithm is based upon the model described in the appendix of
14 : !! Toon et al., J. Geophys. Res., 94, 11359-11380, 1989.
15 : !!
16 : !! This routine assumes that the following variable/tables have already been
17 : !! set up:
18 : !!
19 : !! <radint> intensity of incoming radiance (solar+ir) [erg/cm2/sr/s/cm]
20 : !! <wave> wavelengths used for integration [cm]
21 : !! <dwave> width of wavelength bands for integration [cm]
22 : !! <do_wave_emit> whether planck emission should be doen for the band
23 : !! <qext> extinction [cm2]
24 : !! <ssa> single scattering albedo
25 : !!
26 : !! @author Chuck Bardeen
27 : !! @version Jan-2010
28 33704986635 : subroutine pheat(carma, cstate, iz, igroup, iepart, ibin, igas, dmdt, rc)
29 :
30 : ! types
31 : use carma_precision_mod
32 : use carma_enums_mod
33 : use carma_constants_mod
34 : use carma_types_mod
35 : use carmastate_mod
36 : use carma_mod
37 :
38 : use planck, only : planckIntensity, planckBandIntensity, planckBandIntensityWidger1976, planckBandIntensityConley2011
39 :
40 : implicit none
41 :
42 :
43 : type(carma_type), intent(in) :: carma !! the carma object
44 : type(carmastate_type), intent(inout) :: cstate !! the carma state object
45 : integer, intent(in) :: iz !! vertical index
46 : integer, intent(in) :: igroup !! group index
47 : integer, intent(in) :: iepart !! group's concentration element index
48 : integer, intent(in) :: ibin !! bin index
49 : integer, intent(in) :: igas !! gas index
50 : real(kind=f), intent(out) :: dmdt !! particle growth rate (g/s)
51 : integer, intent(inout) :: rc !! return code, negative indicates failure
52 :
53 : ! Local declarations
54 : integer, parameter :: MAX_ITER = 10 ! Maximum number of iterations
55 : real(kind=f), parameter :: DDTP_LIMIT = 0.01_f ! Convergence criteria for iteration.
56 :
57 : integer :: iter ! iteration
58 : integer :: iwvl ! wavelength band index
59 67409973270 : integer :: ieother(NELEM)
60 : integer :: nother
61 : integer :: ieoth_rel
62 : integer :: ieoth_abs
63 : integer :: jother
64 : integer :: isol
65 67409973270 : real(kind=f) :: otherm(NELEM)
66 : real(kind=f) :: argsol
67 : real(kind=f) :: othermtot
68 : real(kind=f) :: othervtot
69 : real(kind=f) :: condm
70 : real(kind=f) :: condv
71 : real(kind=f) :: volfrc
72 : real(kind=f) :: akas
73 : real(kind=f) :: expon
74 : real(kind=f) :: g0
75 : real(kind=f) :: g1
76 : real(kind=f) :: g2
77 : real(kind=f) :: ss
78 : real(kind=f) :: pvap
79 : real(kind=f) :: qrad ! particle net radiation (erg/s)
80 : ! real(kind=f) :: qrad0 ! particle net radiation (Tp=Ta) (erg/s)
81 : real(kind=f) :: rlh ! latent heat (erg/g)
82 : real(kind=f) :: tp ! particle temperature (K)
83 : real(kind=f) :: dtp ! change in particle temperature (K)
84 : real(kind=f) :: dtpl ! last change in particle temperature (K)
85 : real(kind=f) :: ddtp ! change in particle temperature in last iteration (K)
86 : real(kind=f) :: plkint ! planck intensity
87 :
88 : ! <akas> is combined kelvin (curvature) and solute factors.
89 : !
90 : ! Ignore solute factor for ice particles.
91 33704986635 : if( is_grp_ice(igroup) )then
92 0 : expon = akelvini(iz,igas) / rup_wet(iz,ibin,igroup)
93 :
94 : ! Ice can't be neutralized, so set the volume fraction to 0.
95 0 : volfrc = 0._f
96 : else
97 :
98 33704986635 : argsol = 0._f
99 :
100 : ! Consider growth of average particle at radius <rup(ibin,igroup)>.
101 : !
102 : ! Treat solute effect first: <asol> is solute factor.
103 : !
104 : ! Only need to treat solute effect if <nelemg(igroup)> > 1
105 33704986635 : if( nelemg(igroup) .gt. 1 )then
106 :
107 : ! <condm> is mass concentration of condensed gas <igas> in particle.
108 : ! <nother> is number of other elements in group having mass.
109 : ! <otherm> are mass concentrations of other elements in particle group.
110 : ! <othermtot> is total mass concentrations of other elements in particle.
111 : nother = 0
112 : othermtot = 0._f
113 : othervtot = 0._f
114 :
115 : ! <ieoth_rel> is relative element number of other element in group.
116 >15625*10^7 : do ieoth_rel = 2,nelemg(igroup)
117 :
118 : ! <ieoth_abs> is absolute element number of other element.
119 >13021*10^7 : ieoth_abs = iepart + ieoth_rel - 1
120 :
121 >15625*10^7 : if( itype(ieoth_abs) .eq. I_COREMASS )then
122 >13021*10^7 : nother = nother + 1
123 >13021*10^7 : ieother(nother) = ieoth_abs
124 >13021*10^7 : otherm(nother) = pc(iz,ibin,ieoth_abs)
125 >13021*10^7 : othermtot = othermtot + otherm(nother)
126 >13021*10^7 : othervtot = othervtot + otherm(nother) / pc(iz,ibin,iepart) / rhoelem(ibin,ieoth_abs)
127 : endif
128 : enddo
129 :
130 26043237471 : condm = rmass(ibin,igroup) * pc(iz,ibin,iepart) - othermtot
131 26043237471 : condv = min(0._f, (rmass(ibin,igroup) / rhoelem(ibin,iepart)) - othervtot)
132 :
133 26043237471 : if( condm .le. 0._f )then
134 :
135 : ! Zero mass for the condensate -- <asol> is a small value << 1
136 : argsol = 1e6_f
137 :
138 : ! If there is no condensed mass, then the volume fraction of core is 1.
139 : volfrc = 1._f
140 : else
141 :
142 : ! Sum over masses of other elements in group for argument of solute factor.
143 >15580*10^7 : do jother = 1,nother
144 >12983*10^7 : isol = isolelem(ieother(jother))
145 :
146 : ! Some elements aren't soluble, so skip them.
147 >15580*10^7 : if(isol .gt. 0 ) argsol = argsol + sol_ions(isol)*otherm(jother)/solwtmol(isol)
148 : enddo
149 :
150 25967279674 : argsol = argsol*gwtmol(igas)/condm
151 :
152 25967279674 : volfrc = othervtot / (othervtot + condv)
153 : endif
154 : endif ! nelemg(igroup) > 1
155 33704986635 : expon = akelvin(iz,igas) / rup_wet(iz,ibin,igroup) - argsol
156 : endif
157 :
158 33704986635 : expon = max(-POWMAX, expon)
159 33704986635 : akas = exp( expon )
160 :
161 : ! Trick for removing haze droplets from droplet bins:
162 : ! allows haze droplets to exist under supersaturated conditions;
163 : ! when below supersaturation, haze droplets will evaporate.
164 : ! if( (.not. is_grp_ice(igroup)) .and. (akas .lt. 1._f) .and. &
165 : ! (supsatl(iz,igas) .lt. 0._f) ) akas = 1._f
166 :
167 : ! <dmdt> is growth rate in mass space [g/s].
168 33704986635 : g0 = gro(iz,ibin+1,igroup)
169 33704986635 : g1 = gro1(iz,ibin+1,igroup)
170 33704986635 : g2 = gro2(iz,igroup)
171 :
172 33704986635 : if( is_grp_ice(igroup) )then
173 0 : ss = supsati(iz,igas)
174 0 : pvap = pvapi(iz,igas)
175 : else
176 33704986635 : ss = supsatl(iz,igas)
177 33704986635 : pvap = pvapl(iz,igas)
178 : endif
179 :
180 :
181 : ! If particle heating is being considered, then determine qrad and tpart to
182 : ! determine dmdt.
183 : !
184 : ! NOTE: If no optical properties, then can't do the particle heating calculation.
185 33704986635 : if ((.not. do_pheat) .or. (.not. do_mie(igroup))) then
186 :
187 : ! Ignore the qrad term.
188 33704986635 : dmdt = pvap * ( ss + 1._f - akas ) * g0 / ( 1._f + g0 * g1 * pvap )
189 :
190 : ! Is neutralization set up for the group?
191 33704986635 : if (neutral_volfrc(igroup) > 0._f) then
192 :
193 : ! When the particle is less than fully neutralized, calculate a new
194 : ! dmdt based upon assuming that the saturation vapor pressure (pvap)
195 : ! is 0.
196 0 : if (volfrc >= neutral_volfrc(igroup)) then
197 0 : dmdt = max((pvap * (ss + 1._f)) * g0, dmdt)
198 : else
199 :
200 : ! You can only lose sulfuric acid (condensed) mass until the volume fraction
201 : ! for neutralization is reached. At that point the particle is fully
202 : ! neutralized and the vapor pressure goes to 0. The volume of condensed gas
203 : ! in excess of full neutralization is:
204 : !
205 : ! condv - othervtot * ((1 - neutral_volfrc) / neutral_volfrc)
206 : !
207 : ! NOTE: Limit the growth rate so that the neutralized volume fraction is
208 : ! not overshot. Test have shown that this requires reducing the rate by a
209 : ! factor of 2; although, other values probably work too.
210 : dmdt = max(-(condv - othervtot * ((1._f - neutral_volfrc(igroup)) / neutral_volfrc(igroup))) &
211 0 : * rhoelem(ibin,iepart) / 2._f / dtime, &
212 0 : dmdt)
213 : end if
214 33704986635 : elseif(neutral_volfrc(igroup) .eq. -1._f)then
215 : ! The particle does not evaporate due to neutralization.
216 : ! for example, in trop_strat model, the sulfate does not
217 : ! evaporate in the mixed particles, because it is assumed
218 : ! that it is neutralized by ammonia.
219 26043237471 : dmdt = max((pvap * (ss + 1._f)) * g0, dmdt)
220 : !write(*,*) "igroup",igroup," in pheat.F90:neutral_volfrc"
221 : end if
222 : else
223 :
224 : ! Latent heat of condensing gas
225 0 : if( is_grp_ice(igroup) )then
226 0 : rlh = rlhe(iz,igas) + rlhm(iz,igas)
227 : else
228 0 : rlh = rlhe(iz,igas)
229 : endif
230 :
231 : ! The particle temperature must be solved for by iterating, with an
232 : ! initial guess that the particle temperature is the ambient temperature.
233 : !
234 : ! NOTE: We could also try a guest that is based upon an equilibrium
235 : ! between upwelling IR and collisonal heating, which was identified by
236 : ! Jensen [1989] as the dominant terms.
237 : !
238 : ! radp = 0.d0
239 : !
240 : ! do iwvl = 1, Nwave
241 : ! radp = radp + (4.0d0*PI * absk(iwvl,ibin+1,igroup) *
242 : ! $ radint3(ixyz,iwvl) * dwave(iwvl))
243 : ! end do
244 : !
245 : ! dtp2 = radp /
246 : ! $ (4.d0*PI*rlow(ibin+1,igroup)*thcondnc(iz)*ft(iz,ibin+1,igroup))
247 0 : tp = t(iz)
248 0 : dtp = 0._f
249 0 : dtpl = 0._f
250 :
251 0 : do iter = 1, MAX_ITER
252 :
253 : ! Calculate the net radiative flux on the particle, which requires
254 : ! integrating the incoming and outgoing flux over the spectral
255 : ! interval.
256 0 : qrad = 0._f
257 :
258 0 : do iwvl = 1, NWAVE
259 :
260 : ! There may be overlap between bands, so only do the emission
261 : ! for each range of wavelengths once.
262 0 : if (do_wave_emit(iwvl)) then
263 :
264 : ! Get an integral across the entire band. There are several
265 : ! techniques for doing this that vary in accuracy and
266 : ! performance. Comments below are based on the CAM RRTMG
267 : ! band structure.
268 :
269 : ! Just use the band center.
270 : !
271 : ! NOTE: This generates about a 20% error, but is the fastest
272 : ! plkint = planckIntensity(wave(iwvl), tp)
273 :
274 : ! Brute Force integral
275 : !
276 : ! The slowest technique, and not as accurate as either Widger
277 : ! and Woodall or Conley, even at 100 iterations.
278 : ! plkint = planckBandIntensity(wave(iwvl), dwave(iwvl), tp, 60)
279 :
280 : ! Integral using Widger and Woodall, 1976.
281 : !
282 : ! NOTE: One of the fastest technique at 2 iterations, but yields errors
283 : ! of about 2%. Can handle wide rage of band sizes.
284 : ! plkint = planckBandIntensityWidger1976(wave(iwvl), dwave(iwvl), tp, 2)
285 :
286 : ! Using method developed by Andrew Conley.
287 : !
288 : ! This is similar in performance to Widger and Woodall, but is more
289 : ! accurate with errors of about 0.3%. It had trouble with SW bands that
290 : ! are very large, but the latest version has improved performance and
291 : ! it does work with the RRTMG band structure.
292 0 : plkint = planckBandIntensityConley2011(wave(iwvl), dwave(iwvl), tp, 1)
293 :
294 : else
295 : plkint = 0._f
296 : end if
297 :
298 0 : qrad = qrad + 4.0_f * PI * (1._f - ssa(iwvl,ibin+1,igroup)) * &
299 0 : qext(iwvl,ibin+1,igroup) * PI * (rlow_wet(iz,ibin+1,igroup) ** 2) &
300 0 : * arat(ibin+1,igroup) * (radint(iz,iwvl) - plkint) * dwave(iwvl)
301 : end do
302 :
303 : ! Save of the Qrad association with the ambient air temperature.
304 : ! if (iter == 0) then
305 : ! qrad0 = qrad
306 : ! end if
307 :
308 : ! Calculate the change in mass using eq. A3 from Toon et al. [1989].
309 : dmdt = pvap * ( ss + 1._f - akas * (1._f + qrad * g1 * g2 )) * &
310 0 : g0 / ( 1._f + g0 * g1 * pvap )
311 :
312 : ! Calculate a new particle temperature based upon the loss of mass and
313 : ! energy being absorbed.
314 0 : if ((dmdt * dtime) .le. (- rmass(ibin+1, igroup))) then
315 : dtp = ((rlh * (- rmass(ibin+1, igroup) / dtime)) + qrad) / &
316 0 : (4._f * PI * rlow_wet(iz,ibin+1,igroup) * thcondnc(iz,ibin+1,igroup) * ft(iz,ibin+1,igroup))
317 : else
318 : dtp = ((rlh * dmdt) + qrad) / &
319 0 : (4._f * PI * rlow_wet(iz,ibin+1,igroup) * thcondnc(iz,ibin+1,igroup) * ft(iz,ibin+1,igroup))
320 : end if
321 :
322 0 : tp = t(iz) + dtp
323 :
324 0 : ddtp = dtp - dtpl
325 0 : dtpl = dtp
326 :
327 0 : if (abs(ddtp) .le. DDTP_LIMIT) then
328 : exit
329 : end if
330 :
331 0 : if ((iter .gt. 1) .and. (ddtp .gt. dtpl)) then
332 : exit
333 : end if
334 : end do
335 :
336 0 : dtpart(iz,ibin,igroup) = dtp
337 :
338 : ! Calculate the contribution of this bin to the heating of the atmosphere. CARMA does
339 : ! not actually apply this heating to change the temperature.
340 : !
341 : ! From Pruppacher & Klett [2000], eq. 13-19, the heat transfer to
342 : ! one particle is:
343 : !
344 : ! dq/dt = 4*pi*r*thcondnc*Ft(r)*(T - Tp(r))
345 : !
346 : ! so the total heating rate of the air by the particle is:
347 : !
348 : ! dT/dt = -Sum((4*pi*r*thcondnc*Ft(r)*(T-Tp(r))*pc(r))) / (Cp,air*arho)
349 : !
350 : ! or
351 : !
352 : ! dT/dt = Sum((4*pi*r*thcondnc*Ft(r)*dtp*pc(r))) / (Cp,air*arho)
353 : !
354 : ! where dtp = Tp(r) - T
355 : !
356 : ! NOTE: Using these terms will cause the model parent model to go out of
357 : ! energy balance, since qrad difference is not being communicated to the
358 : ! other layers.
359 0 : if (do_pheatatm) then
360 :
361 : ! NOTE: If the particle is going to evaporate entirely during the timestep,
362 : ! then assume that there is no particle heating.
363 0 : if ((dmdt * dtime) .gt. (- rmass(ibin+1, igroup))) then
364 :
365 : ! If the particles are radiatively active, then the parent model's radiation
366 : ! code is calculated based upon Ta, not Tp. Adjust for this error in Qrad.
367 : ! phprod = phprod + (qrad - qrad0) * pc(iz,ibin+1,iepart) / CP / rhoa(iz)
368 :
369 : ! Now add in the heating from thermal conduction.
370 0 : phprod = phprod + 4._f * PI * rlow_wet(iz,ibin+1,igroup) * &
371 0 : thcondnc(iz,ibin+1,igroup) * ft(iz,ibin+1,igroup) * dtp * &
372 0 : pc(iz,ibin+1,iepart) / (CP * rhoa(iz))
373 : end if
374 : end if
375 : end if
376 :
377 : ! Return to caller with particle loss rates for growth and evaporation
378 : ! evaluated.
379 33704986635 : return
380 33704986635 : end
|
