| blob | 12c3f6642e789bfc9f93e7707c3f48b129f8830b |
1 #define MODAL_AERO
2 #define WRF_PORT
3 ! module_cam_mam_newnuc.F
4 ! adapted from cam3 modal_aero_newnuc.F90 by r.c.easter, august 2010
5 ! Updated to CESM1.0.3 (CAM5.1.01) by Balwinder.Singh@pnnl.gov
6 !
7 ! 2010-08-04 rce notes
8 ! > change pcnstxx to a subr argument because gas_pcnst is
9 ! not a constant/parameter in wrf-chem
10 !--------------------------------------------------------------
11 #include "MODAL_AERO_CPP_DEFINES.h"
13 ! modal_aero_newnuc.F90
16 !----------------------------------------------------------------------
17 !BOP
18 !
19 ! !MODULE: modal_aero_newnuc --- modal aerosol new-particle nucleation
20 !
21 ! !INTERFACE:
22 module modal_aero_newnuc
23 #if (defined MODAL_AERO)
25 ! !USES:
26 use shr_kind_mod, only: r8 => shr_kind_r8
27 use shr_kind_mod, only: r4 => shr_kind_r4
28 #ifndef WRF_PORT
29 use mo_constants, only: pi
30 use chem_mods, only: gas_pcnst
31 #else
32 use physconst, only: pi
33 use module_cam_support, only: gas_pcnst => gas_pcnst_modal_aero
34 #endif
36 implicit none
37 private
38 save
40 ! !PUBLIC MEMBER FUNCTIONS:
41 public modal_aero_newnuc_sub, modal_aero_newnuc_init
43 ! !PUBLIC DATA MEMBERS:
44 #ifndef WRF_PORT
45 integer, parameter :: pcnstxx = gas_pcnst
46 #endif
47 integer :: l_h2so4_sv, l_nh3_sv, lnumait_sv, lnh4ait_sv, lso4ait_sv
49 ! min h2so4 vapor for nuc calcs = 4.0e-16 mol/mol-air ~= 1.0e4 molecules/cm3,
50 real(r8), parameter :: qh2so4_cutoff = 4.0e-16_r8
52 ! !DESCRIPTION: This module implements ...
53 !
54 ! !REVISION HISTORY:
55 !
56 ! R.Easter 2007.09.14: Adapted from MIRAGE2 code
57 !
58 !EOP
59 !----------------------------------------------------------------------
60 !BOC
62 ! list private module data here
64 !EOC
65 !----------------------------------------------------------------------
68 contains
70 !----------------------------------------------------------------------
71 !----------------------------------------------------------------------
72 !BOP
73 ! !ROUTINE: modal_aero_newnuc_sub --- ...
74 !
75 ! !INTERFACE:
76 subroutine modal_aero_newnuc_sub( &
77 lchnk, ncol, nstep, &
78 loffset, deltat, &
79 latndx, lonndx, &
80 t, pmid, pdel, &
81 zm, pblh, &
82 qv, cld, &
83 q, &
84 del_h2so4_gasprod, del_h2so4_aeruptk &
85 #ifdef WRF_PORT
86 , pcnstxx &
87 #endif
88 )
91 ! !USES:
92 use modal_aero_data
93 #ifndef WRF_PORT
94 use abortutils, only: endrun
95 use cam_history, only: outfld, fieldname_len
96 use chem_mods, only: adv_mass
97 use constituents, only: pcnst, cnst_name
98 #endif
99 use physconst, only: gravit, mwdry, r_universal
100 #ifndef WRF_PORT
101 use ppgrid, only: pcols, pver
102 use spmd_utils, only: iam, masterproc
103 #else
104 use module_cam_support, only: pcnst => pcnst_runtime, &
105 pcols, pver, &
106 fieldname_len, &
107 endrun, iam, masterproc, outfld
108 use constituents, only: cnst_name
109 use module_data_cam_mam_asect, only: adv_mass => mw_q_mo_array
110 #endif
111 use wv_saturation, only: aqsat
113 implicit none
115 ! !PARAMETERS:
116 #ifdef WRF_PORT
117 integer, intent(in) :: pcnstxx ! number of species in "incoming q array"
118 #endif
119 integer, intent(in) :: lchnk ! chunk identifier
120 integer, intent(in) :: ncol ! number of columns in chunk
121 integer, intent(in) :: nstep ! model step
122 integer, intent(in) :: latndx(pcols), lonndx(pcols)
123 integer, intent(in) :: loffset ! offset applied to modal aero "pointers"
124 real(r8), intent(in) :: deltat ! model timestep (s)
126 real(r8), intent(in) :: t(pcols,pver) ! temperature (K)
127 real(r8), intent(in) :: pmid(pcols,pver) ! pressure at model levels (Pa)
128 real(r8), intent(in) :: pdel(pcols,pver) ! pressure thickness of levels (Pa)
129 real(r8), intent(in) :: zm(pcols,pver) ! midpoint height above surface (m)
130 real(r8), intent(in) :: pblh(pcols) ! pbl height (m)
131 real(r8), intent(in) :: qv(pcols,pver) ! specific humidity (kg/kg)
132 real(r8), intent(in) :: cld(ncol,pver) ! stratiform cloud fraction
133 ! *** NOTE ncol dimension
134 real(r8), intent(inout) :: q(ncol,pver,pcnstxx)
135 ! tracer mixing ratio (TMR) array
136 ! *** MUST BE mol/mol-air or #/mol-air
137 ! *** NOTE ncol & pcnstxx dimensions
138 real(r8), intent(in) :: del_h2so4_gasprod(ncol,pver)
139 ! h2so4 gas-phase production
140 ! change over deltat (mol/mol)
141 real(r8), intent(in) :: del_h2so4_aeruptk(ncol,pver)
142 ! h2so4 gas-phase loss to
143 ! aerosol over deltat (mol/mol)
145 ! !DESCRIPTION:
146 ! computes changes due to aerosol nucleation (new particle formation)
147 ! treats both nucleation and subsequent growth of new particles
148 ! to aitken mode size
149 ! uses the following parameterizations
150 ! vehkamaki et al. (2002) parameterization for binary
151 ! homogeneous nucleation (h2so4-h2o) plus
152 ! kerminen and kulmala (2002) parameterization for
153 ! new particle loss during growth to aitken size
154 !
155 ! !REVISION HISTORY:
156 ! R.Easter 2007.09.14: Adapted from MIRAGE2 code and CMAQ V4.6 code
157 !
158 !EOP
159 !----------------------------------------------------------------------
160 !BOC
162 ! local variables
163 integer :: i, itmp, k, l, lmz, lun, m, mait
164 integer :: lnumait, lso4ait, lnh4ait
165 integer :: l_h2so4, l_nh3
166 integer :: ldiagveh02
167 integer, parameter :: ldiag1=-1, ldiag2=-1, ldiag3=-1, ldiag4=-1
168 integer, parameter :: newnuc_method_flagaa = 11
169 ! integer, parameter :: newnuc_method_flagaa = 12
170 ! 1=merikanto et al (2007) ternary 2=vehkamaki et al (2002) binary
171 ! 11=merikanto ternary + first-order boundary layer
172 ! 12=merikanto ternary + second-order boundary layer
174 real(r8) :: adjust_factor
175 real(r8) :: aircon
176 real(r8) :: cldx
177 real(r8) :: dens_nh4so4a
178 real(r8) :: dmdt_ait, dmdt_aitsv1, dmdt_aitsv2, dmdt_aitsv3
179 real(r8) :: dndt_ait, dndt_aitsv1, dndt_aitsv2, dndt_aitsv3
180 real(r8) :: dnh4dt_ait, dso4dt_ait
181 real(r8) :: dpnuc
182 real(r8) :: dplom_mode(1), dphim_mode(1)
183 real(r8) :: ev_sat(pcols,pver)
184 real(r8) :: mass1p
185 real(r8) :: mass1p_aithi, mass1p_aitlo
186 real(r8) :: mw_so4a_host
187 real(r8) :: pdel_fac
188 real(r8) :: qh2so4_cur, qh2so4_avg, qh2so4_del
189 real(r8) :: qnh3_cur, qnh3_del, qnh4a_del
190 real(r8) :: qnuma_del
191 real(r8) :: qso4a_del
192 real(r8) :: qv_sat(pcols,pver)
193 real(r8) :: qvswtr
194 real(r8) :: relhum, relhumav, relhumnn
195 real(r8) :: tmpa, tmpb, tmpc
196 real(r8) :: tmp_q1, tmp_q2, tmp_q3
197 real(r8) :: tmp_frso4, tmp_uptkrate
199 integer, parameter :: nsrflx = 1 ! last dimension of qsrflx
200 real(r8) :: qsrflx(pcols,pcnst,nsrflx)
201 ! process-specific column tracer tendencies
202 ! 1 = nucleation (for aerocom)
203 real(r8) :: dqdt(ncol,pver,pcnstxx) ! TMR tendency array -- NOTE dims
204 logical :: dotend(pcnst) ! flag for doing tendency
205 logical :: do_nh3 ! flag for doing nh3/nh4
208 character(len=1) :: tmpch1, tmpch2, tmpch3
209 character(len=fieldname_len+3) :: fieldname
212 ! begin
213 lun = 6
215 !--------------------------------------------------------------------------------
216 if (ldiag1 > 0) then
217 do i = 1, ncol
218 if (lonndx(i) /= 37) cycle
219 if (latndx(i) /= 23) cycle
220 if (nstep > 3) cycle
221 write( lun, '(/a,i7,3i5,f10.2)' ) &
222 '*** modal_aero_newnuc_sub -- nstep, iam, lat, lon =', &
223 nstep, iam, latndx(i), lonndx(i)
224 end do
225 if (nstep > 3) call endrun( '*** modal_aero_newnuc_sub -- testing halt after step 3' )
226 ! if (ncol /= -999888777) return
227 end if
228 !--------------------------------------------------------------------------------
230 !-----------------------------------------------------------------------
231 l_h2so4 = l_h2so4_sv - loffset
232 l_nh3 = l_nh3_sv - loffset
233 lnumait = lnumait_sv - loffset
234 lnh4ait = lnh4ait_sv - loffset
235 lso4ait = lso4ait_sv - loffset
237 ! skip if no aitken mode OR if no h2so4 species
238 if ((l_h2so4 <= 0) .or. (lso4ait <= 0) .or. (lnumait <= 0)) return
240 dotend(:) = .false.
241 dqdt(1:ncol,:,:) = 0.0
242 qsrflx(1:ncol,:,:) = 0.0
244 ! set dotend
245 mait = modeptr_aitken
246 dotend(lnumait) = .true.
247 dotend(lso4ait) = .true.
248 dotend(l_h2so4) = .true.
250 lnh4ait = lptr_nh4_a_amode(mait) - loffset
251 if ((l_nh3 > 0) .and. (l_nh3 <= pcnst) .and. &
252 (lnh4ait > 0) .and. (lnh4ait <= pcnst)) then
253 do_nh3 = .true.
254 dotend(lnh4ait) = .true.
255 dotend(l_nh3) = .true.
256 else
257 do_nh3 = .false.
258 end if
261 ! dry-diameter limits for "grown" new particles
262 dplom_mode(1) = exp( 0.67*log(dgnumlo_amode(mait)) &
263 + 0.33*log(dgnum_amode(mait)) )
264 dphim_mode(1) = dgnumhi_amode(mait)
266 ! mass1p_... = mass (kg) of so4 & nh4 in a single particle of diameter ...
267 ! (assuming same dry density for so4 & nh4)
268 ! mass1p_aitlo - dp = dplom_mode(1)
269 ! mass1p_aithi - dp = dphim_mode(1)
270 tmpa = specdens_so4_amode*pi/6.0
271 mass1p_aitlo = tmpa*(dplom_mode(1)**3)
272 mass1p_aithi = tmpa*(dphim_mode(1)**3)
274 ! compute qv_sat = saturation specific humidity
275 call aqsat( t, pmid, ev_sat, qv_sat, pcols, ncol, pver, 1, pver )
277 ! mw_so4a_host is molec-wght of sulfate aerosol in host code
278 ! 96 when nh3/nh4 are simulated
279 ! something else when nh3/nh4 are not simulated
280 mw_so4a_host = specmw_so4_amode
283 !
284 ! loop over levels and columns to calc the renaming
285 !
286 main_k: do k = 1, pver
287 main_i: do i = 1, ncol
289 ! skip if completely cloudy,
290 ! because all h2so4 vapor should be cloud-borne
291 if (cld(i,k) >= 0.99) cycle main_i
293 ! qh2so4_cur = current qh2so4, after aeruptk
294 qh2so4_cur = q(i,k,l_h2so4)
295 ! skip if h2so4 vapor < qh2so4_cutoff
296 if (qh2so4_cur <= qh2so4_cutoff) cycle main_i
298 tmpa = max( 0.0_r8, del_h2so4_gasprod(i,k) )
299 tmp_q3 = qh2so4_cur
300 ! tmp_q2 = qh2so4 before aeruptk
301 ! (note tmp_q3, tmp_q2 both >= 0.0)
302 tmp_q2 = tmp_q3 + max( 0.0_r8, -del_h2so4_aeruptk(i,k) )
304 ! *** temporary -- in order to get more nucleation
305 ! qh2so4_cur = qh2so4_cur*1.0e1
306 ! tmp_q3 = tmp_q3*1.0e1
307 ! tmp_q2 = tmp_q2*1.0e1
308 ! tmpa = tmpa *1.0e1
310 ! tmpb = log( tmp_q2/tmp_q3 ) BUT with some checks added
311 ! tmp_uptkrate = tmpb/deltat
312 if (tmp_q2 <= tmp_q3) then
313 tmpb = 0.0
314 else
315 tmpc = tmp_q2 * exp( -20.0_r8 )
316 if (tmp_q3 <= tmpc) then
317 tmp_q3 = tmpc
318 tmpb = 20.0_r8
319 else
320 tmpb = log( tmp_q2/tmp_q3 )
321 end if
322 end if
323 ! d[ln(qh2so4)]/dt (1/s) from uptake (condensation) to aerosol
324 tmp_uptkrate = tmpb/deltat
326 ! qh2so4_avg = estimated average qh2so4
327 ! when production & loss are done simultaneously
328 if (tmpb <= 0.1_r8) then
329 qh2so4_avg = tmp_q3*(1.0_r8 + 0.5_r8*tmpb) - 0.5_r8*tmpa
330 else
331 tmpc = tmpa/tmpb
332 qh2so4_avg = (tmp_q3 - tmpc)*((exp(tmpb)-1.0_r8)/tmpb) + tmpc
333 end if
334 if (qh2so4_avg <= qh2so4_cutoff) cycle main_i
337 if ( do_nh3 ) then
338 qnh3_cur = max( 0.0_r8, q(i,k,l_nh3) )
339 else
340 qnh3_cur = 0.0_r8
341 end if
344 ! relhumav = grid average RH
345 qvswtr = qv_sat(i,k)
346 qvswtr = max( qvswtr, 1.0d-20 )
347 relhumav = qv(i,k) / qvswtr
348 relhumav = max( 0.0_r8, min( 1.0_r8, relhumav ) )
349 ! relhum = non-cloudy area RH
350 cldx = max( 0.0_r8, cld(i,k) )
351 relhum = (relhumav - cldx) / (1.0_r8 - cldx)
352 relhum = max( 0.0_r8, min( 1.0_r8, relhum ) )
353 ! limit RH to between 0.1% and 99%
354 relhumnn = relhum
355 relhumnn = max( 0.01_r8, min( 0.99_r8, relhumnn ) )
357 ! aircon = air concentration (mol-air/m3)
358 aircon = 1.0e3*pmid(i,k)/(r_universal*t(i,k))
361 ! call ... routine to get nucleation rates
362 ldiagveh02 = -1
363 if (ldiag2 > 0) then
364 if ((lonndx(i) == 37) .and. (latndx(i) == 23)) then
365 if ((k >= 24) .or. (mod(k,4) == 0)) then
366 ldiagveh02 = +1
367 write(lun,'(/a,i8,3i4,f8.2,1p,4e10.2)') &
368 'veh02 call - nstep,lat,lon,k; tk,rh,p,cair', &
369 nstep, latndx(i), lonndx(i), k, &
370 t(i,k), relhumnn, pmid(k,k), aircon
371 end if
372 end if
373 end if
374 call mer07_veh02_nuc_mosaic_1box( &
375 newnuc_method_flagaa, &
376 deltat, t(i,k), relhumnn, pmid(i,k), &
377 zm(i,k), pblh(i), &
378 qh2so4_cur, qh2so4_avg, qnh3_cur, tmp_uptkrate, &
379 mw_so4a_host, &
380 1, 1, dplom_mode, dphim_mode, &
381 itmp, qnuma_del, qso4a_del, qnh4a_del, &
382 qh2so4_del, qnh3_del, dens_nh4so4a, ldiagveh02 )
383 ! qh2so4_del, qnh3_del, dens_nh4so4a )
384 !----------------------------------------------------------------------
385 ! subr mer07_veh02_nuc_mosaic_1box( &
386 ! newnuc_method_flagaa, &
387 ! dtnuc, temp_in, rh_in, press_in, &
388 ! qh2so4_cur, qh2so4_avg, qnh3_cur, h2so4_uptkrate, &
389 ! nsize, maxd_asize, dplom_sect, dphim_sect, &
390 ! isize_nuc, qnuma_del, qso4a_del, qnh4a_del, &
391 ! qh2so4_del, qnh3_del, dens_nh4so4a )
392 !
393 !! subr arguments (in)
394 ! real(r8), intent(in) :: dtnuc ! nucleation time step (s)
395 ! real(r8), intent(in) :: temp_in ! temperature, in k
396 ! real(r8), intent(in) :: rh_in ! relative humidity, as fraction
397 ! real(r8), intent(in) :: press_in ! air pressure (pa)
398 !
399 ! real(r8), intent(in) :: qh2so4_cur, qh2so4_avg
400 ! ! gas h2so4 mixing ratios (mol/mol-air)
401 ! real(r8), intent(in) :: qnh3_cur ! gas nh3 mixing ratios (mol/mol-air)
402 ! ! qxxx_cur = current value (after gas chem and condensation)
403 ! ! qxxx_avg = estimated average value (for simultaneous source/sink calcs)
404 ! real(r8), intent(in) :: h2so4_uptkrate ! h2so4 uptake rate to aerosol (1/s)
406 !
407 ! integer, intent(in) :: nsize ! number of aerosol size bins
408 ! integer, intent(in) :: maxd_asize ! dimension for dplom_sect, ...
409 ! real(r8), intent(in) :: dplom_sect(maxd_asize) ! dry diameter at lower bnd of bin (m)
410 ! real(r8), intent(in) :: dphim_sect(maxd_asize) ! dry diameter at upper bnd of bin (m)
411 !
412 !! subr arguments (out)
413 ! integer, intent(out) :: isize_nuc ! size bin into which new particles go
414 ! real(r8), intent(out) :: qnuma_del ! change to aerosol number mixing ratio (#/mol-air)
415 ! real(r8), intent(out) :: qso4a_del ! change to aerosol so4 mixing ratio (mol/mol-air)
416 ! real(r8), intent(out) :: qnh4a_del ! change to aerosol nh4 mixing ratio (mol/mol-air)
417 ! real(r8), intent(out) :: qh2so4_del ! change to gas h2so4 mixing ratio (mol/mol-air)
418 ! real(r8), intent(out) :: qnh3_del ! change to gas nh3 mixing ratio (mol/mol-air)
419 ! ! aerosol changes are > 0; gas changes are < 0
420 ! real(r8), intent(out) :: dens_nh4so4a ! dry-density of the new nh4-so4 aerosol mass (kg/m3)
421 !----------------------------------------------------------------------
424 ! convert qnuma_del from (#/mol-air) to (#/kmol-air)
425 qnuma_del = qnuma_del*1.0e3_r8
426 ! number nuc rate (#/kmol-air/s) from number nuc amt
427 dndt_ait = qnuma_del/deltat
428 ! fraction of mass nuc going to so4
429 tmpa = qso4a_del*specmw_so4_amode
430 tmpb = tmpa + qnh4a_del*specmw_nh4_amode
431 tmp_frso4 = max( tmpa, 1.0e-35_r8 )/max( tmpb, 1.0e-35_r8 )
432 ! mass nuc rate (kg/kmol-air/s or g/mol...) hhfrom mass nuc amts
433 dmdt_ait = max( 0.0_r8, (tmpb/deltat) )
435 dndt_aitsv1 = dndt_ait
436 dmdt_aitsv1 = dmdt_ait
437 dndt_aitsv2 = 0.0
438 dmdt_aitsv2 = 0.0
439 dndt_aitsv3 = 0.0
440 dmdt_aitsv3 = 0.0
441 tmpch1 = ' '
442 tmpch2 = ' '
444 if (dndt_ait < 1.0e2) then
445 ! ignore newnuc if number rate < 100 #/kmol-air/s ~= 0.3 #/mg-air/d
446 dndt_ait = 0.0
447 dmdt_ait = 0.0
448 tmpch1 = 'A'
450 else
451 dndt_aitsv2 = dndt_ait
452 dmdt_aitsv2 = dmdt_ait
453 tmpch1 = 'B'
455 ! mirage2 code checked for complete h2so4 depletion here,
456 ! but this is now done in mer07_veh02_nuc_mosaic_1box
457 mass1p = dmdt_ait/dndt_ait
458 dndt_aitsv3 = dndt_ait
459 dmdt_aitsv3 = dmdt_ait
461 ! apply particle size constraints
462 if (mass1p < mass1p_aitlo) then
463 ! reduce dndt to increase new particle size
464 dndt_ait = dmdt_ait/mass1p_aitlo
465 tmpch1 = 'C'
466 else if (mass1p > mass1p_aithi) then
467 ! reduce dmdt to decrease new particle size
468 dmdt_ait = dndt_ait*mass1p_aithi
469 tmpch1 = 'E'
470 end if
471 end if
473 ! *** apply adjustment factor to avoid unrealistically high
474 ! aitken number concentrations in mid and upper troposphere
475 ! adjust_factor = 0.5
476 ! dndt_ait = dndt_ait * adjust_factor
477 ! dmdt_ait = dmdt_ait * adjust_factor
479 ! set tendencies
480 pdel_fac = pdel(i,k)/gravit
482 ! dso4dt_ait, dnh4dt_ait are (kmol/kmol-air/s)
483 dso4dt_ait = dmdt_ait*tmp_frso4/specmw_so4_amode
484 dnh4dt_ait = dmdt_ait*(1.0_r8 - tmp_frso4)/specmw_nh4_amode
486 dqdt(i,k,l_h2so4) = -dso4dt_ait*(1.0-cldx)
487 qsrflx(i,l_h2so4,1) = qsrflx(i,l_h2so4,1) + dqdt(i,k,l_h2so4)*pdel_fac
488 q(i,k,l_h2so4) = q(i,k,l_h2so4) + dqdt(i,k,l_h2so4)*deltat
490 dqdt(i,k,lso4ait) = dso4dt_ait*(1.0-cldx)
491 qsrflx(i,lso4ait,1) = qsrflx(i,lso4ait,1) + dqdt(i,k,lso4ait)*pdel_fac
492 q(i,k,lso4ait) = q(i,k,lso4ait) + dqdt(i,k,lso4ait)*deltat
493 if (lnumait > 0) then
494 dqdt(i,k,lnumait) = dndt_ait*(1.0-cldx)
495 qsrflx(i,lnumait,1) = qsrflx(i,lnumait,1) &
496 + dqdt(i,k,lnumait)*pdel_fac
497 q(i,k,lnumait) = q(i,k,lnumait) + dqdt(i,k,lnumait)*deltat
498 end if
500 if (( do_nh3 ) .and. (dnh4dt_ait > 0.0_r8)) then
501 dqdt(i,k,l_nh3) = -dnh4dt_ait*(1.0-cldx)
502 qsrflx(i,l_nh3,1) = qsrflx(i,l_nh3,1) + dqdt(i,k,l_nh3)*pdel_fac
503 q(i,k,l_nh3) = q(i,k,l_nh3) + dqdt(i,k,l_nh3)*deltat
505 dqdt(i,k,lnh4ait) = dnh4dt_ait*(1.0-cldx)
506 qsrflx(i,lnh4ait,1) = qsrflx(i,lnh4ait,1) + dqdt(i,k,lnh4ait)*pdel_fac
507 q(i,k,lnh4ait) = q(i,k,lnh4ait) + dqdt(i,k,lnh4ait)*deltat
508 end if
510 !! temporary diagnostic
511 ! if (ldiag3 > 0) then
512 ! if ((dndt_ait /= 0.0_r8) .or. (dmdt_ait /= 0.0_r8)) then
513 ! write(lun,'(3a,1x,i7,3i5,1p,5e12.4)') &
514 ! 'newnucxx', tmpch1, tmpch2, nstep, lchnk, i, k, &
515 ! dndt_ait, dmdt_ait, cldx
516 !! call endrun( 'modal_aero_newnuc_sub' )
517 ! end if
518 ! end if
521 ! diagnostic output start ----------------------------------------
522 if (ldiag4 > 0) then
523 if ((lonndx(i) == 37) .and. (latndx(i) == 23)) then
524 if ((k >= 24) .or. (mod(k,4) == 0)) then
525 write(lun,97010) nstep, latndx(i), lonndx(i), k, t(i,k), aircon
526 write(lun,97020) 'pmid, pdel ', &
527 pmid(i,k), pdel(i,k)
528 write(lun,97030) 'qv,qvsw, cld, rh_av, rh_clr ', &
529 qv(i,k), qvswtr, cldx, relhumav, relhum
530 write(lun,97020) 'h2so4_cur, _pre, _av, nh3_cur', &
531 qh2so4_cur, tmp_q2, qh2so4_avg, qnh3_cur
532 write(lun,97020) 'del_h2so4_gasprod, _aeruptk ', &
533 del_h2so4_gasprod(i,k), del_h2so4_aeruptk(i,k), &
534 tmp_uptkrate*3600.0
535 write(lun,97020) ' '
536 write(lun,97050) 'tmpch1, tmpch2 ', tmpch1, tmpch2
537 write(lun,97020) 'dndt_, dmdt_aitsv1 ', &
538 dndt_aitsv1, dmdt_aitsv1
539 write(lun,97020) 'dndt_, dmdt_aitsv2 ', &
540 dndt_aitsv2, dmdt_aitsv2
541 write(lun,97020) 'dndt_, dmdt_aitsv3 ', &
542 dndt_aitsv3, dmdt_aitsv3
543 write(lun,97020) 'dndt_, dmdt_ait ', &
544 dndt_ait, dmdt_ait
545 write(lun,97020) 'dso4dt_, dnh4dt_ait ', &
546 dso4dt_ait, dnh4dt_ait
547 write(lun,97020) 'qso4a_del, qh2so4_del ', &
548 qso4a_del, qh2so4_del
549 write(lun,97020) 'qnh4a_del, qnh3_del ', &
550 qnh4a_del, qnh3_del
551 write(lun,97020) 'dqdt(h2so4), (nh3) ', &
552 dqdt(i,k,l_h2so4), dqdt(i,k,l_nh3)
553 write(lun,97020) 'dqdt(so4a), (nh4a), (numa) ', &
554 dqdt(i,k,lso4ait), dqdt(i,k,lnh4ait), dqdt(i,k,lnumait)
556 dpnuc = 0.0
557 if (dndt_aitsv1 > 1.0e-5) dpnuc = (6.0*dmdt_aitsv1/ &
558 (pi*specdens_so4_amode*dndt_aitsv1))**0.3333333
559 if (dpnuc > 0.0) then
560 write(lun,97020) 'dpnuc, dp_aitlo, _aithi ', &
561 dpnuc, dplom_mode(1), dphim_mode(1)
562 write(lun,97020) 'mass1p, mass1p_aitlo, _aithi ', &
563 mass1p, mass1p_aitlo, mass1p_aithi
564 end if
566 97010 format( / 'NEWNUC nstep,lat,lon,k,tk,cair', i8, 3i4, f8.2, 1pe12.4 )
567 97020 format( a, 1p, 6e12.4 )
568 97030 format( a, 1p, 2e12.4, 0p, 5f10.6 )
569 97040 format( 29x, 1p, 6e12.4 )
570 97050 format( a, 2(3x,a) )
571 end if
572 end if
573 end if
574 ! diagnostic output end ------------------------------------------
577 end do main_i
578 end do main_k
581 ! do history file column-tendency fields
582 do l = loffset+1, pcnst
583 lmz = l - loffset
584 if ( .not. dotend(lmz) ) cycle
586 do i = 1, ncol
587 qsrflx(i,lmz,1) = qsrflx(i,lmz,1)*(adv_mass(lmz)/mwdry)
588 end do
589 fieldname = trim(cnst_name(l)) // '_sfnnuc1'
590 call outfld( fieldname, qsrflx(:,lmz,1), pcols, lchnk )
592 ! if (( masterproc ) .and. (nstep < 1)) &
593 ! write(lun,'(2(a,2x),1p,e11.3)') &
594 ! 'modal_aero_newnuc_sub outfld', fieldname, adv_mass(lmz)
595 end do ! l = ...
598 return
599 !EOC
600 end subroutine modal_aero_newnuc_sub
604 !----------------------------------------------------------------------
605 !-----------------------------------------------------------------------
606 subroutine mer07_veh02_nuc_mosaic_1box( &
607 newnuc_method_flagaa, dtnuc, temp_in, rh_in, press_in, &
608 zm_in, pblh_in, &
609 qh2so4_cur, qh2so4_avg, qnh3_cur, h2so4_uptkrate, &
610 mw_so4a_host, &
611 nsize, maxd_asize, dplom_sect, dphim_sect, &
612 isize_nuc, qnuma_del, qso4a_del, qnh4a_del, &
613 qh2so4_del, qnh3_del, dens_nh4so4a, ldiagaa )
614 ! qh2so4_del, qnh3_del, dens_nh4so4a )
615 !.......................................................................
616 !
617 ! calculates new particle production from homogeneous nucleation
618 ! over timestep dtnuc, using nucleation rates from either
619 ! merikanto et al. (2007) h2so4-nh3-h2o ternary parameterization
620 ! vehkamaki et al. (2002) h2so4-h2o binary parameterization
621 !
622 ! the new particles are "grown" to the lower-bound size of the host code's
623 ! smallest size bin. (this "growth" is somewhat ad hoc, and would not be
624 ! necessary if the host code's size bins extend down to ~1 nm.)
625 !
626 ! if the h2so4 and nh3 mass mixing ratios (mixrats) of the grown new
627 ! particles exceed the current gas mixrats, the new particle production
628 ! is reduced so that the new particle mass mixrats match the gas mixrats.
629 !
630 ! the correction of kerminen and kulmala (2002) is applied to account
631 ! for loss of the new particles by coagulation as they are
632 ! growing to the "host code mininum size"
633 !
634 ! revision history
635 ! coded by rc easter, pnnl, xx-apr-2007
636 !
637 ! key routines called: subr ternary_nuc_napari
638 !
639 ! references:
640 ! merikanto, j., i. napari, h. vehkamaki, t. anttila,
641 ! and m. kulmala, 2007, new parameterization of
642 ! sulfuric acid-ammonia-water ternary nucleation
643 ! rates at tropospheric conditions,
644 ! j. geophys. res., 112, d15207, doi:10.1029/2006jd0027977
645 !
646 ! vehkamaki, h., m. kulmala, i. napari, k.e.j. lehtinen,
647 ! c. timmreck, m. noppel and a. laaksonen, 2002,
648 ! an improved parameterization for sulfuric acid-water nucleation
649 ! rates for tropospheric and stratospheric conditions,
650 ! j. geophys. res., 107, 4622, doi:10.1029/2002jd002184
651 !
652 ! kerminen, v., and m. kulmala, 2002,
653 ! analytical formulae connecting the "real" and the "apparent"
654 ! nucleation rate and the nuclei number concentration
655 ! for atmospheric nucleation events
656 !
657 !.......................................................................
658 implicit none
660 ! subr arguments (in)
661 real(r8), intent(in) :: dtnuc ! nucleation time step (s)
662 real(r8), intent(in) :: temp_in ! temperature, in k
663 real(r8), intent(in) :: rh_in ! relative humidity, as fraction
664 real(r8), intent(in) :: press_in ! air pressure (pa)
665 real(r8), intent(in) :: zm_in ! layer midpoint height (m)
666 real(r8), intent(in) :: pblh_in ! pbl height (m)
668 real(r8), intent(in) :: qh2so4_cur, qh2so4_avg
669 ! gas h2so4 mixing ratios (mol/mol-air)
670 real(r8), intent(in) :: qnh3_cur ! gas nh3 mixing ratios (mol/mol-air)
671 ! qxxx_cur = current value (after gas chem and condensation)
672 ! qxxx_avg = estimated average value (for simultaneous source/sink calcs)
673 real(r8), intent(in) :: h2so4_uptkrate ! h2so4 uptake rate to aerosol (1/s)
674 real(r8), intent(in) :: mw_so4a_host ! mw of so4 aerosol in host code (g/mol)
676 integer, intent(in) :: newnuc_method_flagaa ! 1=merikanto et al (2007) ternary
677 ! 2=vehkamaki et al (2002) binary
678 integer, intent(in) :: nsize ! number of aerosol size bins
679 integer, intent(in) :: maxd_asize ! dimension for dplom_sect, ...
680 real(r8), intent(in) :: dplom_sect(maxd_asize) ! dry diameter at lower bnd of bin (m)
681 real(r8), intent(in) :: dphim_sect(maxd_asize) ! dry diameter at upper bnd of bin (m)
682 integer, intent(in) :: ldiagaa
684 ! subr arguments (out)
685 integer, intent(out) :: isize_nuc ! size bin into which new particles go
686 real(r8), intent(out) :: qnuma_del ! change to aerosol number mixing ratio (#/mol-air)
687 real(r8), intent(out) :: qso4a_del ! change to aerosol so4 mixing ratio (mol/mol-air)
688 real(r8), intent(out) :: qnh4a_del ! change to aerosol nh4 mixing ratio (mol/mol-air)
689 real(r8), intent(out) :: qh2so4_del ! change to gas h2so4 mixing ratio (mol/mol-air)
690 real(r8), intent(out) :: qnh3_del ! change to gas nh3 mixing ratio (mol/mol-air)
691 ! aerosol changes are > 0; gas changes are < 0
692 real(r8), intent(out) :: dens_nh4so4a ! dry-density of the new nh4-so4 aerosol mass (kg/m3)
694 ! subr arguments (out) passed via common block
695 ! these are used to duplicate the outputs of yang zhang's original test driver
696 ! they are not really needed in wrf-chem
697 real(r8) :: ratenuclt ! j = ternary nucleation rate from napari param. (cm-3 s-1)
698 real(r8) :: rateloge ! ln (j)
699 real(r8) :: cnum_h2so4 ! number of h2so4 molecules in the critical nucleus
700 real(r8) :: cnum_nh3 ! number of nh3 molecules in the critical nucleus
701 real(r8) :: cnum_tot ! total number of molecules in the critical nucleus
702 real(r8) :: radius_cluster ! the radius of cluster (nm)
705 ! local variables
706 integer :: i
707 integer :: igrow
708 integer, save :: icase = 0, icase_reldiffmax = 0
709 ! integer, parameter :: ldiagaa = -1
710 integer :: lun
711 integer :: newnuc_method_flagaa2
713 real(r8), parameter :: onethird = 1.0/3.0
714 real(r8), parameter :: avogad = 6.022e23 ! avogadro number (molecules/mol)
715 real(r8), parameter :: mw_air = 28.966 ! dry-air mean molecular weight (g/mol)
717 real(r8), parameter :: accom_coef_h2so4 = 0.65 ! accomodation coef for h2so4 conden
719 ! dry densities (kg/m3) molecular weights of aerosol
720 ! ammsulf, ammbisulf, and sulfacid (from mosaic dens_electrolyte values)
721 ! real(r8), parameter :: dens_ammsulf = 1.769e3
722 ! real(r8), parameter :: dens_ammbisulf = 1.78e3
723 ! real(r8), parameter :: dens_sulfacid = 1.841e3
724 ! use following to match cam3 modal_aero densities
725 real(r8), parameter :: dens_ammsulf = 1.770e3
726 real(r8), parameter :: dens_ammbisulf = 1.770e3
727 real(r8), parameter :: dens_sulfacid = 1.770e3
728 real(r8), parameter :: dens_water = 1.0e3
730 ! molecular weights (g/mol) of aerosol ammsulf, ammbisulf, and sulfacid
731 ! for ammbisulf and sulfacid, use 114 & 96 here rather than 115 & 98
732 ! because we don't keep track of aerosol hion mass
733 real(r8), parameter :: mw_ammsulf = 132.0
734 real(r8), parameter :: mw_ammbisulf = 114.0
735 real(r8), parameter :: mw_sulfacid = 96.0
736 ! molecular weights of aerosol sulfate and ammonium
737 real(r8), parameter :: mw_so4a = 96.0
738 real(r8), parameter :: mw_nh4a = 18.0
739 real(r8), parameter :: mw_water = 18.0
741 real(r8), save :: reldiffmax = 0.0
743 real(r8) cair ! dry-air molar density (mol/m3)
744 real(r8) cs_prime_kk ! kk2002 "cs_prime" parameter (1/m2)
745 real(r8) cs_kk ! kk2002 "cs" parameter (1/s)
746 real(r8) dens_part ! "grown" single-particle dry density (kg/m3)
747 real(r8) dfin_kk, dnuc_kk ! kk2002 final/initial new particle wet diameter (nm)
748 real(r8) dpdry_clus ! critical cluster diameter (m)
749 real(r8) dpdry_part ! "grown" single-particle dry diameter (m)
750 real(r8) tmpa, tmpb, tmpc, tmpe, tmpq
751 real(r8) tmpa1, tmpb1
752 real(r8) tmp_m1, tmp_m2, tmp_m3, tmp_n1, tmp_n2, tmp_n3
753 real(r8) tmp_spd ! h2so4 vapor molecular speed (m/s)
754 real(r8) factor_kk
755 real(r8) fogas, foso4a, fonh4a, fonuma
756 real(r8) freduce ! reduction factor applied to nucleation rate
757 ! due to limited availability of h2so4 & nh3 gases
758 real(r8) freducea, freduceb
759 real(r8) gamma_kk ! kk2002 "gamma" parameter (nm2*m2/h)
760 real(r8) gr_kk ! kk2002 "gr" parameter (nm/h)
761 real(r8) kgaero_per_moleso4a ! (kg dry aerosol)/(mol aerosol so4)
762 real(r8) mass_part ! "grown" single-particle dry mass (kg)
763 real(r8) molenh4a_per_moleso4a ! (mol aerosol nh4)/(mol aerosol so4)
764 real(r8) nh3ppt, nh3ppt_bb ! actual and bounded nh3 (ppt)
765 real(r8) nu_kk ! kk2002 "nu" parameter (nm)
766 real(r8) qmolnh4a_del_max ! max production of aerosol nh4 over dtnuc (mol/mol-air)
767 real(r8) qmolso4a_del_max ! max production of aerosol so4 over dtnuc (mol/mol-air)
768 real(r8) ratenuclt_bb ! nucleation rate (#/m3/s)
769 real(r8) ratenuclt_kk ! nucleation rate after kk2002 adjustment (#/m3/s)
770 real(r8) rh_bb ! bounded value of rh_in
771 real(r8) so4vol_in ! concentration of h2so4 for nucl. calc., molecules cm-3
772 real(r8) so4vol_bb ! bounded value of so4vol_in
773 real(r8) temp_bb ! bounded value of temp_in
774 real(r8) voldry_clus ! critical-cluster dry volume (m3)
775 real(r8) voldry_part ! "grown" single-particle dry volume (m3)
776 real(r8) wetvol_dryvol ! grown particle (wet-volume)/(dry-volume)
777 real(r8) wet_volfrac_so4a ! grown particle (dry-volume-from-so4)/(wet-volume)
781 !
782 ! if h2so4 vapor < qh2so4_cutoff
783 ! exit with new particle formation = 0
784 !
785 isize_nuc = 1
786 qnuma_del = 0.0
787 qso4a_del = 0.0
788 qnh4a_del = 0.0
789 qh2so4_del = 0.0
790 qnh3_del = 0.0
791 ! if (qh2so4_avg .le. qh2so4_cutoff) return ! this no longer needed
792 ! if (qh2so4_cur .le. qh2so4_cutoff) return ! this no longer needed
794 if ((newnuc_method_flagaa /= 1) .and. &
795 (newnuc_method_flagaa /= 2) .and. &
796 (newnuc_method_flagaa /= 11) .and. &
797 (newnuc_method_flagaa /= 12)) return
800 !
801 ! make call to parameterization routine
802 !
804 ! calc h2so4 in molecules/cm3 and nh3 in ppt
805 cair = press_in/(temp_in*8.3144)
806 so4vol_in = qh2so4_avg * cair * avogad * 1.0e-6
807 nh3ppt = qnh3_cur * 1.0e12
808 ratenuclt = 1.0e-38_r8
809 rateloge = log( ratenuclt )
811 if ( (newnuc_method_flagaa /= 2) .and. &
812 (nh3ppt >= 0.1) ) then
813 ! make call to merikanto ternary parameterization routine
814 ! (when nh3ppt < 0.1, use binary param instead)
816 if (so4vol_in >= 5.0e4) then
817 temp_bb = max( 235.0_r8, min( 295.0_r8, temp_in ) )
818 rh_bb = max( 0.05_r8, min( 0.95_r8, rh_in ) )
819 so4vol_bb = max( 5.0e4_r8, min( 1.0e9_r8, so4vol_in ) )
820 nh3ppt_bb = max( 0.1_r8, min( 1.0e3_r8, nh3ppt ) )
821 call ternary_nuc_merik2007( &
822 temp_bb, rh_bb, so4vol_bb, nh3ppt_bb, &
823 rateloge, &
824 cnum_tot, cnum_h2so4, cnum_nh3, radius_cluster )
825 end if
826 newnuc_method_flagaa2 = 1
828 else
829 ! make call to vehkamaki binary parameterization routine
831 if (so4vol_in >= 1.0e4) then
832 temp_bb = max( 230.15_r8, min( 305.15_r8, temp_in ) )
833 rh_bb = max( 1.0e-4_r8, min( 1.0_r8, rh_in ) )
834 so4vol_bb = max( 1.0e4_r8, min( 1.0e11_r8, so4vol_in ) )
835 call binary_nuc_vehk2002( &
836 temp_bb, rh_bb, so4vol_bb, &
837 ratenuclt, rateloge, &
838 cnum_h2so4, cnum_tot, radius_cluster )
839 end if
840 cnum_nh3 = 0.0
841 newnuc_method_flagaa2 = 2
843 end if
846 ! do boundary layer nuc
847 if ((newnuc_method_flagaa == 11) .or. &
848 (newnuc_method_flagaa == 12)) then
849 if ( zm_in <= max(pblh_in,100.0_r8) ) then
850 so4vol_bb = so4vol_in
851 call pbl_nuc_wang2008( so4vol_bb, &
852 newnuc_method_flagaa, newnuc_method_flagaa2, &
853 ratenuclt, rateloge, &
854 cnum_tot, cnum_h2so4, cnum_nh3, radius_cluster )
855 end if
856 end if
859 ! if nucleation rate is less than 1e-6 #/m3/s ~= 0.1 #/cm3/day,
860 ! exit with new particle formation = 0
861 if (rateloge .le. -13.82_r8) return
862 ! if (ratenuclt .le. 1.0e-6) return
863 ratenuclt = exp( rateloge )
864 ratenuclt_bb = ratenuclt*1.0e6_r8
867 ! wet/dry volume ratio - use simple kohler approx for ammsulf/ammbisulf
868 tmpa = max( 0.10_r8, min( 0.95_r8, rh_in ) )
869 wetvol_dryvol = 1.0 - 0.56/log(tmpa)
872 ! determine size bin into which the new particles go
873 ! (probably it will always be bin #1, but ...)
874 voldry_clus = ( max(cnum_h2so4,1.0_r8)*mw_so4a + cnum_nh3*mw_nh4a ) / &
875 (1.0e3*dens_sulfacid*avogad)
876 ! correction when host code sulfate is really ammonium bisulfate/sulfate
877 voldry_clus = voldry_clus * (mw_so4a_host/mw_so4a)
878 dpdry_clus = (voldry_clus*6.0/pi)**onethird
880 isize_nuc = 1
881 dpdry_part = dplom_sect(1)
882 if (dpdry_clus <= dplom_sect(1)) then
883 igrow = 1 ! need to clusters to larger size
884 else if (dpdry_clus >= dphim_sect(nsize)) then
885 igrow = 0
886 isize_nuc = nsize
887 dpdry_part = dphim_sect(nsize)
888 else
889 igrow = 0
890 do i = 1, nsize
891 if (dpdry_clus < dphim_sect(i)) then
892 isize_nuc = i
893 dpdry_part = dpdry_clus
894 dpdry_part = min( dpdry_part, dphim_sect(i) )
895 dpdry_part = max( dpdry_part, dplom_sect(i) )
896 exit
897 end if
898 end do
899 end if
900 voldry_part = (pi/6.0)*(dpdry_part**3)
903 !
904 ! determine composition and density of the "grown particles"
905 ! the grown particles are assumed to be liquid
906 ! (since critical clusters contain water)
907 ! so any (nh4/so4) molar ratio between 0 and 2 is allowed
908 ! assume that the grown particles will have
909 ! (nh4/so4 molar ratio) = min( 2, (nh3/h2so4 gas molar ratio) )
910 !
911 if (igrow .le. 0) then
912 ! no "growing" so pure sulfuric acid
913 tmp_n1 = 0.0
914 tmp_n2 = 0.0
915 tmp_n3 = 1.0
916 else if (qnh3_cur .ge. qh2so4_cur) then
917 ! combination of ammonium sulfate and ammonium bisulfate
918 ! tmp_n1 & tmp_n2 = mole fractions of the ammsulf & ammbisulf
919 tmp_n1 = (qnh3_cur/qh2so4_cur) - 1.0
920 tmp_n1 = max( 0.0_r8, min( 1.0_r8, tmp_n1 ) )
921 tmp_n2 = 1.0 - tmp_n1
922 tmp_n3 = 0.0
923 else
924 ! combination of ammonium bisulfate and sulfuric acid
925 ! tmp_n2 & tmp_n3 = mole fractions of the ammbisulf & sulfacid
926 tmp_n1 = 0.0
927 tmp_n2 = (qnh3_cur/qh2so4_cur)
928 tmp_n2 = max( 0.0_r8, min( 1.0_r8, tmp_n2 ) )
929 tmp_n3 = 1.0 - tmp_n2
930 end if
932 tmp_m1 = tmp_n1*mw_ammsulf
933 tmp_m2 = tmp_n2*mw_ammbisulf
934 tmp_m3 = tmp_n3*mw_sulfacid
935 dens_part = (tmp_m1 + tmp_m2 + tmp_m3)/ &
936 ((tmp_m1/dens_ammsulf) + (tmp_m2/dens_ammbisulf) &
937 + (tmp_m3/dens_sulfacid))
938 dens_nh4so4a = dens_part
939 mass_part = voldry_part*dens_part
940 ! (mol aerosol nh4)/(mol aerosol so4)
941 molenh4a_per_moleso4a = 2.0*tmp_n1 + tmp_n2
942 ! (kg dry aerosol)/(mol aerosol so4)
943 kgaero_per_moleso4a = 1.0e-3*(tmp_m1 + tmp_m2 + tmp_m3)
944 ! correction when host code sulfate is really ammonium bisulfate/sulfate
945 kgaero_per_moleso4a = kgaero_per_moleso4a * (mw_so4a_host/mw_so4a)
947 ! fraction of wet volume due to so4a
948 tmpb = 1.0 + molenh4a_per_moleso4a*17.0/98.0
949 wet_volfrac_so4a = 1.0 / ( wetvol_dryvol * tmpb )
952 !
953 ! calc kerminen & kulmala (2002) correction
954 !
955 if (igrow <= 0) then
956 factor_kk = 1.0
958 else
959 ! "gr" parameter (nm/h) = condensation growth rate of new particles
960 ! use kk2002 eqn 21 for h2so4 uptake, and correct for nh3 & h2o uptake
961 tmp_spd = 14.7*sqrt(temp_in) ! h2so4 molecular speed (m/s)
962 gr_kk = 3.0e-9*tmp_spd*mw_sulfacid*so4vol_in/ &
963 (dens_part*wet_volfrac_so4a)
965 ! "gamma" parameter (nm2/m2/h)
966 ! use kk2002 eqn 22
967 !
968 ! dfin_kk = wet diam (nm) of grown particle having dry dia = dpdry_part (m)
969 dfin_kk = 1.0e9 * dpdry_part * (wetvol_dryvol**onethird)
970 ! dnuc_kk = wet diam (nm) of cluster
971 dnuc_kk = 2.0*radius_cluster
972 dnuc_kk = max( dnuc_kk, 1.0_r8 )
973 ! neglect (dmean/150)**0.048 factor,
974 ! which should be very close to 1.0 because of small exponent
975 gamma_kk = 0.23 * (dnuc_kk)**0.2 &
976 * (dfin_kk/3.0)**0.075 &
977 * (dens_part*1.0e-3)**(-0.33) &
978 * (temp_in/293.0)**(-0.75)
980 ! "cs_prime parameter" (1/m2)
981 ! instead kk2002 eqn 3, use
982 ! cs_prime ~= tmpa / (4*pi*tmpb * h2so4_accom_coef)
983 ! where
984 ! tmpa = -d(ln(h2so4))/dt by conden to particles (1/h units)
985 ! tmpb = h2so4 vapor diffusivity (m2/h units)
986 ! this approx is generally within a few percent of the cs_prime
987 ! calculated directly from eqn 2,
988 ! which is acceptable, given overall uncertainties
989 ! tmpa = -d(ln(h2so4))/dt by conden to particles (1/h units)
990 tmpa = h2so4_uptkrate * 3600.0
991 tmpa1 = tmpa
992 tmpa = max( tmpa, 0.0_r8 )
993 ! tmpb = h2so4 gas diffusivity (m2/s, then m2/h)
994 tmpb = 6.7037e-6 * (temp_in**0.75) / cair
995 tmpb1 = tmpb ! m2/s
996 tmpb = tmpb*3600.0 ! m2/h
997 cs_prime_kk = tmpa/(4.0*pi*tmpb*accom_coef_h2so4)
998 cs_kk = cs_prime_kk*4.0*pi*tmpb1
1000 ! "nu" parameter (nm) -- kk2002 eqn 11
1001 nu_kk = gamma_kk*cs_prime_kk/gr_kk
1002 ! nucleation rate adjustment factor (--) -- kk2002 eqn 13
1003 factor_kk = exp( (nu_kk/dfin_kk) - (nu_kk/dnuc_kk) )
1005 end if
1006 ratenuclt_kk = ratenuclt_bb*factor_kk
1009 ! max production of aerosol dry mass (kg-aero/m3-air)
1010 tmpa = max( 0.0_r8, (ratenuclt_kk*dtnuc*mass_part) )
1011 ! max production of aerosol so4 (mol-so4a/mol-air)
1012 tmpe = tmpa/(kgaero_per_moleso4a*cair)
1013 ! max production of aerosol so4 (mol/mol-air)
1014 ! based on ratenuclt_kk and mass_part
1015 qmolso4a_del_max = tmpe
1017 ! check if max production exceeds available h2so4 vapor
1018 freducea = 1.0
1019 if (qmolso4a_del_max .gt. qh2so4_cur) then
1020 freducea = qh2so4_cur/qmolso4a_del_max
1021 end if
1023 ! check if max production exceeds available nh3 vapor
1024 freduceb = 1.0
1025 if (molenh4a_per_moleso4a .ge. 1.0e-10) then
1026 ! max production of aerosol nh4 (ppm) based on ratenuclt_kk and mass_part
1027 qmolnh4a_del_max = qmolso4a_del_max*molenh4a_per_moleso4a
1028 if (qmolnh4a_del_max .gt. qnh3_cur) then
1029 freduceb = qnh3_cur/qmolnh4a_del_max
1030 end if
1031 end if
1032 freduce = min( freducea, freduceb )
1034 ! if adjusted nucleation rate is less than 1e-12 #/m3/s ~= 0.1 #/cm3/day,
1035 ! exit with new particle formation = 0
1036 if (freduce*ratenuclt_kk .le. 1.0e-12) return
1039 ! note: suppose that at this point, freduce < 1.0 (no gas-available
1040 ! constraints) and molenh4a_per_moleso4a < 2.0
1041 ! if the gas-available constraints is do to h2so4 availability,
1042 ! then it would be possible to condense "additional" nh3 and have
1043 ! (nh3/h2so4 gas molar ratio) < (nh4/so4 aerosol molar ratio) <= 2
1044 ! one could do some additional calculations of
1045 ! dens_part & molenh4a_per_moleso4a to realize this
1046 ! however, the particle "growing" is a crude approximate way to get
1047 ! the new particles to the host code's minimum particle size,
1048 ! are such refinements worth the effort?
1051 ! changes to h2so4 & nh3 gas (in mol/mol-air), limited by amounts available
1052 tmpa = 0.9999
1053 qh2so4_del = min( tmpa*qh2so4_cur, freduce*qmolso4a_del_max )
1054 qnh3_del = min( tmpa*qnh3_cur, qh2so4_del*molenh4a_per_moleso4a )
1055 qh2so4_del = -qh2so4_del
1056 qnh3_del = -qnh3_del
1058 ! changes to so4 & nh4 aerosol (in mol/mol-air)
1059 qso4a_del = -qh2so4_del
1060 qnh4a_del = -qnh3_del
1061 ! change to aerosol number (in #/mol-air)
1062 qnuma_del = 1.0e-3*(qso4a_del*mw_so4a + qnh4a_del*mw_nh4a)/mass_part
1064 ! do the following (tmpa, tmpb, tmpc) calculations as a check
1065 ! max production of aerosol number (#/mol-air)
1066 tmpa = max( 0.0_r8, (ratenuclt_kk*dtnuc/cair) )
1067 ! adjusted production of aerosol number (#/mol-air)
1068 tmpb = tmpa*freduce
1069 ! relative difference from qnuma_del
1070 tmpc = (tmpb - qnuma_del)/max(tmpb, qnuma_del, 1.0e-35_r8)
1073 !
1074 ! diagnostic output to fort.41
1075 ! (this should be commented-out or deleted in the wrf-chem version)
1076 !
1077 if (ldiagaa <= 0) return
1079 icase = icase + 1
1080 if (abs(tmpc) .gt. abs(reldiffmax)) then
1081 reldiffmax = tmpc
1082 icase_reldiffmax = icase
1083 end if
1084 ! do lun = 41, 51, 10
1085 do lun = 6, 6
1086 ! write(lun,'(/)')
1087 write(lun,'(a,2i9,1p,e10.2)') &
1088 'vehkam bin-nuc icase, icase_rdmax =', &
1089 icase, icase_reldiffmax, reldiffmax
1090 if (freduceb .lt. freducea) then
1091 if (abs(freducea-freduceb) .gt. &
1092 3.0e-7*max(freduceb,freducea)) write(lun,'(a,1p,2e15.7)') &
1093 'freducea, b =', freducea, freduceb
1094 end if
1095 end do
1097 ! output factors so that output matches that of ternucl03
1098 ! fogas = 1.0e6 ! convert mol/mol-air to ppm
1099 ! foso4a = 1.0e9*mw_so4a/mw_air ! convert mol-so4a/mol-air to ug/kg-air
1100 ! fonh4a = 1.0e9*mw_nh4a/mw_air ! convert mol-nh4a/mol-air to ug/kg-air
1101 ! fonuma = 1.0e3/mw_air ! convert #/mol-air to #/kg-air
1102 fogas = 1.0
1103 foso4a = 1.0
1104 fonh4a = 1.0
1105 fonuma = 1.0
1107 ! do lun = 41, 51, 10
1108 do lun = 6, 6
1110 write(lun,'(a,2i5)') 'newnuc_method_flagaa/aa2', &
1111 newnuc_method_flagaa, newnuc_method_flagaa2
1113 write(lun,9210)
1114 write(lun,9201) temp_in, rh_in, &
1115 ratenuclt, 2.0*radius_cluster*1.0e-7, dpdry_part*1.0e2, &
1116 voldry_part*1.0e6, float(igrow)
1117 write(lun,9215)
1118 write(lun,9201) &
1119 qh2so4_avg*fogas, 0.0, &
1120 qh2so4_cur*fogas, qnh3_cur*fogas, &
1121 qh2so4_del*fogas, qnh3_del*fogas, &
1122 qso4a_del*foso4a, qnh4a_del*fonh4a
1124 write(lun,9220)
1125 write(lun,9201) &
1126 dtnuc, dens_nh4so4a*1.0e-3, &
1127 (qnh3_cur/qh2so4_cur), molenh4a_per_moleso4a, &
1128 qnuma_del*fonuma, tmpb*fonuma, tmpc, freduce
1130 end do
1132 ! lun = 51
1133 lun = 6
1134 write(lun,9230)
1135 write(lun,9201) &
1136 press_in, cair*1.0e-6, so4vol_in, &
1137 wet_volfrac_so4a, wetvol_dryvol, dens_part*1.0e-3
1139 if (igrow > 0) then
1140 write(lun,9240)
1141 write(lun,9201) &
1142 tmp_spd, gr_kk, dnuc_kk, dfin_kk, &
1143 gamma_kk, tmpa1, tmpb1, cs_kk
1145 write(lun,9250)
1146 write(lun,9201) &
1147 cs_prime_kk, nu_kk, factor_kk, ratenuclt, &
1148 ratenuclt_kk*1.0e-6
1149 end if
1151 9201 format ( 1p, 40e10.2 )
1152 9210 format ( &
1153 ' temp rh', &
1154 ' ratenuc dia_clus ddry_part', &
1155 ' vdry_part igrow' )
1156 9215 format ( &
1157 ' h2so4avg h2so4pre', &
1158 ' h2so4cur nh3_cur', &
1159 ' h2so4del nh3_del', &
1160 ' so4a_del nh4a_del' )
1161 9220 format ( &
1162 ' dtnuc dens_a nh/so g nh/so a', &
1163 ' numa_del numa_dl2 reldiff freduce' )
1164 9230 format ( &
1165 ' press_in cair so4_volin', &
1166 ' wet_volfr wetv_dryv dens_part' )
1167 9240 format ( &
1168 ' tmp_spd gr_kk dnuc_kk dfin_kk', &
1169 ' gamma_kk tmpa1 tmpb1 cs_kk' )
1170 9250 format ( &
1171 ' cs_pri_kk nu_kk factor_kk ratenuclt', &
1172 ' ratenu_kk' )
1175 return
1176 end subroutine mer07_veh02_nuc_mosaic_1box
1180 !-----------------------------------------------------------------------
1181 !-----------------------------------------------------------------------
1182 subroutine pbl_nuc_wang2008( so4vol, &
1183 newnuc_method_flagaa, newnuc_method_flagaa2, &
1184 ratenucl, rateloge, &
1185 cnum_tot, cnum_h2so4, cnum_nh3, radius_cluster )
1186 !
1187 ! calculates boundary nucleation nucleation rate
1188 ! using the first or second-order parameterization in
1189 ! wang, m., and j.e. penner, 2008,
1190 ! aerosol indirect forcing in a global model with particle nucleation,
1191 ! atmos. chem. phys. discuss., 8, 13943-13998
1192 !
1193 implicit none
1195 ! subr arguments (in)
1196 real(r8), intent(in) :: so4vol ! concentration of h2so4 (molecules cm-3)
1197 integer, intent(in) :: newnuc_method_flagaa
1198 ! [11,12] value selects [first,second]-order parameterization
1200 ! subr arguments (inout)
1201 integer, intent(inout) :: newnuc_method_flagaa2
1202 real(r8), intent(inout) :: ratenucl ! binary nucleation rate, j (# cm-3 s-1)
1203 real(r8), intent(inout) :: rateloge ! log( ratenucl )
1205 real(r8), intent(inout) :: cnum_tot ! total number of molecules
1206 ! in the critical nucleus
1207 real(r8), intent(inout) :: cnum_h2so4 ! number of h2so4 molecules
1208 real(r8), intent(inout) :: cnum_nh3 ! number of nh3 molecules
1209 real(r8), intent(inout) :: radius_cluster ! the radius of cluster (nm)
1212 ! local variables
1213 real(r8) :: tmp_diam, tmp_mass, tmp_volu
1214 real(r8) :: tmp_rateloge, tmp_ratenucl
1216 ! executable
1219 ! nucleation rate
1220 if (newnuc_method_flagaa == 11) then
1221 tmp_ratenucl = 1.0e-6_r8 * so4vol
1222 else if (newnuc_method_flagaa == 12) then
1223 tmp_ratenucl = 1.0e-12_r8 * (so4vol**2)
1224 else
1225 return
1226 end if
1227 tmp_rateloge = log( tmp_ratenucl )
1229 ! exit if pbl nuc rate is lower than (incoming) ternary/binary rate
1230 if (tmp_rateloge <= rateloge) return
1232 rateloge = tmp_rateloge
1233 ratenucl = tmp_ratenucl
1234 newnuc_method_flagaa2 = newnuc_method_flagaa
1236 ! following wang 2002, assume fresh nuclei are 1 nm diameter
1237 ! subsequent code will "grow" them to aitken mode size
1238 radius_cluster = 0.5_r8
1240 ! assume fresh nuclei are pure h2so4
1241 ! since aitken size >> initial size, the initial composition
1242 ! has very little impact on the results
1243 tmp_diam = radius_cluster * 2.0e-7_r8 ! diameter in cm
1244 tmp_volu = (tmp_diam**3) * (pi/6.0_r8) ! volume in cm^3
1245 tmp_mass = tmp_volu * 1.8_r8 ! mass in g
1246 cnum_h2so4 = (tmp_mass / 98.0_r8) * 6.023e23_r8 ! no. of h2so4 molec assuming pure h2so4
1247 cnum_tot = cnum_h2so4
1248 cnum_nh3 = 0.0_r8
1251 return
1252 end subroutine pbl_nuc_wang2008
1256 !-----------------------------------------------------------------------
1257 !-----------------------------------------------------------------------
1258 subroutine binary_nuc_vehk2002( temp, rh, so4vol, &
1259 ratenucl, rateloge, &
1260 cnum_h2so4, cnum_tot, radius_cluster )
1261 !
1262 ! calculates binary nucleation rate and critical cluster size
1263 ! using the parameterization in
1264 ! vehkamaki, h., m. kulmala, i. napari, k.e.j. lehtinen,
1265 ! c. timmreck, m. noppel and a. laaksonen, 2002,
1266 ! an improved parameterization for sulfuric acid-water nucleation
1267 ! rates for tropospheric and stratospheric conditions,
1268 ! j. geophys. res., 107, 4622, doi:10.1029/2002jd002184
1269 !
1270 implicit none
1272 ! subr arguments (in)
1273 real(r8), intent(in) :: temp ! temperature (k)
1274 real(r8), intent(in) :: rh ! relative humidity (0-1)
1275 real(r8), intent(in) :: so4vol ! concentration of h2so4 (molecules cm-3)
1277 ! subr arguments (out)
1278 real(r8), intent(out) :: ratenucl ! binary nucleation rate, j (# cm-3 s-1)
1279 real(r8), intent(out) :: rateloge ! log( ratenucl )
1281 real(r8), intent(out) :: cnum_h2so4 ! number of h2so4 molecules
1282 ! in the critical nucleus
1283 real(r8), intent(out) :: cnum_tot ! total number of molecules
1284 ! in the critical nucleus
1285 real(r8), intent(out) :: radius_cluster ! the radius of cluster (nm)
1288 ! local variables
1289 real(r8) :: crit_x
1290 real(r8) :: acoe, bcoe, ccoe, dcoe, ecoe, fcoe, gcoe, hcoe, icoe, jcoe
1291 real(r8) :: tmpa, tmpb
1293 ! executable
1296 ! calc sulfuric acid mole fraction in critical cluster
1297 crit_x = 0.740997 - 0.00266379 * temp &
1298 - 0.00349998 * log (so4vol) &
1299 + 0.0000504022 * temp * log (so4vol) &
1300 + 0.00201048 * log (rh) &
1301 - 0.000183289 * temp * log (rh) &
1302 + 0.00157407 * (log (rh)) ** 2.0 &
1303 - 0.0000179059 * temp * (log (rh)) ** 2.0 &
1304 + 0.000184403 * (log (rh)) ** 3.0 &
1305 - 1.50345e-6 * temp * (log (rh)) ** 3.0
1308 ! calc nucleation rate
1309 acoe = 0.14309+2.21956*temp &
1310 - 0.0273911 * temp**2.0 &
1311 + 0.0000722811 * temp**3.0 + 5.91822/crit_x
1313 bcoe = 0.117489 + 0.462532 *temp &
1314 - 0.0118059 * temp**2.0 &
1315 + 0.0000404196 * temp**3.0 + 15.7963/crit_x
1317 ccoe = -0.215554-0.0810269 * temp &
1318 + 0.00143581 * temp**2.0 &
1319 - 4.7758e-6 * temp**3.0 &
1320 - 2.91297/crit_x
1322 dcoe = -3.58856+0.049508 * temp &
1323 - 0.00021382 * temp**2.0 &
1324 + 3.10801e-7 * temp**3.0 &
1325 - 0.0293333/crit_x
1327 ecoe = 1.14598 - 0.600796 * temp &
1328 + 0.00864245 * temp**2.0 &
1329 - 0.0000228947 * temp**3.0 &
1330 - 8.44985/crit_x
1332 fcoe = 2.15855 + 0.0808121 * temp &
1333 -0.000407382 * temp**2.0 &
1334 -4.01957e-7 * temp**3.0 &
1335 + 0.721326/crit_x
1337 gcoe = 1.6241 - 0.0160106 * temp &
1338 + 0.0000377124 * temp**2.0 &
1339 + 3.21794e-8 * temp**3.0 &
1340 - 0.0113255/crit_x
1342 hcoe = 9.71682 - 0.115048 * temp &
1343 + 0.000157098 * temp**2.0 &
1344 + 4.00914e-7 * temp**3.0 &
1345 + 0.71186/crit_x
1347 icoe = -1.05611 + 0.00903378 * temp &
1348 - 0.0000198417 * temp**2.0 &
1349 + 2.46048e-8 * temp**3.0 &
1350 - 0.0579087/crit_x
1352 jcoe = -0.148712 + 0.00283508 * temp &
1353 - 9.24619e-6 * temp**2.0 &
1354 + 5.00427e-9 * temp**3.0 &
1355 - 0.0127081/crit_x
1357 tmpa = ( &
1358 acoe &
1359 + bcoe * log (rh) &
1360 + ccoe * ( log (rh))**2.0 &
1361 + dcoe * ( log (rh))**3.0 &
1362 + ecoe * log (so4vol) &
1363 + fcoe * (log (rh)) * (log (so4vol)) &
1364 + gcoe * ((log (rh) ) **2.0) &
1365 * (log (so4vol)) &
1366 + hcoe * (log (so4vol)) **2.0 &
1367 + icoe * log (rh) &
1368 * ((log (so4vol)) **2.0) &
1369 + jcoe * (log (so4vol)) **3.0 &
1370 )
1371 rateloge = tmpa
1372 tmpa = min( tmpa, log(1.0e38_r8) )
1373 ratenucl = exp ( tmpa )
1374 ! write(*,*) 'tmpa, ratenucl =', tmpa, ratenucl
1378 ! calc number of molecules in critical cluster
1379 acoe = -0.00295413 - 0.0976834*temp &
1380 + 0.00102485 * temp**2.0 &
1381 - 2.18646e-6 * temp**3.0 - 0.101717/crit_x
1383 bcoe = -0.00205064 - 0.00758504*temp &
1384 + 0.000192654 * temp**2.0 &
1385 - 6.7043e-7 * temp**3.0 - 0.255774/crit_x
1387 ccoe = +0.00322308 + 0.000852637 * temp &
1388 - 0.0000154757 * temp**2.0 &
1389 + 5.66661e-8 * temp**3.0 &
1390 + 0.0338444/crit_x
1392 dcoe = +0.0474323 - 0.000625104 * temp &
1393 + 2.65066e-6 * temp**2.0 &
1394 - 3.67471e-9 * temp**3.0 &
1395 - 0.000267251/crit_x
1397 ecoe = -0.0125211 + 0.00580655 * temp &
1398 - 0.000101674 * temp**2.0 &
1399 + 2.88195e-7 * temp**3.0 &
1400 + 0.0942243/crit_x
1402 fcoe = -0.038546 - 0.000672316 * temp &
1403 + 2.60288e-6 * temp**2.0 &
1404 + 1.19416e-8 * temp**3.0 &
1405 - 0.00851515/crit_x
1407 gcoe = -0.0183749 + 0.000172072 * temp &
1408 - 3.71766e-7 * temp**2.0 &
1409 - 5.14875e-10 * temp**3.0 &
1410 + 0.00026866/crit_x
1412 hcoe = -0.0619974 + 0.000906958 * temp &
1413 - 9.11728e-7 * temp**2.0 &
1414 - 5.36796e-9 * temp**3.0 &
1415 - 0.00774234/crit_x
1417 icoe = +0.0121827 - 0.00010665 * temp &
1418 + 2.5346e-7 * temp**2.0 &
1419 - 3.63519e-10 * temp**3.0 &
1420 + 0.000610065/crit_x
1422 jcoe = +0.000320184 - 0.0000174762 * temp &
1423 + 6.06504e-8 * temp**2.0 &
1424 - 1.4177e-11 * temp**3.0 &
1425 + 0.000135751/crit_x
1427 cnum_tot = exp ( &
1428 acoe &
1429 + bcoe * log (rh) &
1430 + ccoe * ( log (rh))**2.0 &
1431 + dcoe * ( log (rh))**3.0 &
1432 + ecoe * log (so4vol) &
1433 + fcoe * (log (rh)) * (log (so4vol)) &
1434 + gcoe * ((log (rh) ) **2.0) &
1435 * (log (so4vol)) &
1436 + hcoe * (log (so4vol)) **2.0 &
1437 + icoe * log (rh) &
1438 * ((log (so4vol)) **2.0) &
1439 + jcoe * (log (so4vol)) **3.0 &
1440 )
1442 cnum_h2so4 = cnum_tot * crit_x
1444 ! calc radius (nm) of critical cluster
1445 radius_cluster = exp( -1.6524245 + 0.42316402*crit_x &
1446 + 0.3346648*log(cnum_tot) )
1449 return
1450 end subroutine binary_nuc_vehk2002
1454 !----------------------------------------------------------------------
1455 !----------------------------------------------------------------------
1456 subroutine modal_aero_newnuc_init
1458 !-----------------------------------------------------------------------
1459 !
1460 ! Purpose:
1461 ! set do_adjust and do_aitken flags
1462 ! create history fields for column tendencies associated with
1463 ! modal_aero_calcsize
1464 !
1465 ! Author: R. Easter
1466 !
1467 !-----------------------------------------------------------------------
1469 use modal_aero_data
1470 use modal_aero_rename
1471 #ifndef WRF_PORT
1472 use abortutils, only: endrun
1473 use cam_history, only: addfld, add_default, fieldname_len, phys_decomp
1474 use constituents, only: pcnst, cnst_get_ind, cnst_name
1475 use spmd_utils, only: masterproc
1476 use phys_control, only: phys_getopts
1477 #else
1478 use module_cam_support, only: pcnst => pcnst_runtime, &
1479 pcols, pver, &
1480 fieldname_len, &
1481 endrun, iam, masterproc,addfld, add_default, &
1482 phys_decomp
1483 use constituents, only: cnst_name, cnst_get_ind
1485 #endif
1488 implicit none
1490 !-----------------------------------------------------------------------
1491 ! arguments
1493 !-----------------------------------------------------------------------
1494 ! local
1495 integer :: l_h2so4, l_nh3
1496 integer :: lnumait, lnh4ait, lso4ait
1497 integer :: l
1498 integer :: m, mait
1500 character(len=fieldname_len) :: tmpname
1501 character(len=fieldname_len+3) :: fieldname
1502 character(128) :: long_name
1503 character(8) :: unit
1505 logical :: dotend(pcnst)
1506 logical :: history_aerosol ! Output the MAM aerosol tendencies
1508 !-----------------------------------------------------------------------
1509 #ifndef WRF_PORT
1510 call phys_getopts( history_aerosol_out = history_aerosol )
1511 #else
1512 history_aerosol = .FALSE.
1513 #endif
1516 ! set these indices
1517 ! skip if no h2so4 species
1518 ! skip if no aitken mode so4 or num species
1519 l_h2so4_sv = 0
1520 l_nh3_sv = 0
1521 lnumait_sv = 0
1522 lnh4ait_sv = 0
1523 lso4ait_sv = 0
1525 call cnst_get_ind( 'H2SO4', l_h2so4, .false. )
1526 call cnst_get_ind( 'NH3', l_nh3, .false. )
1528 mait = modeptr_aitken
1529 if (mait > 0) then
1530 lnumait = numptr_amode(mait)
1531 lso4ait = lptr_so4_a_amode(mait)
1532 lnh4ait = lptr_nh4_a_amode(mait)
1533 end if
1534 if ((l_h2so4 <= 0) .or. (l_h2so4 > pcnst)) then
1535 write(*,'(/a/)') &
1536 '*** modal_aero_newnuc bypass -- l_h2so4 <= 0'
1537 return
1538 else if ((lso4ait <= 0) .or. (lso4ait > pcnst)) then
1539 write(*,'(/a/)') &
1540 '*** modal_aero_newnuc bypass -- lso4ait <= 0'
1541 return
1542 else if ((lnumait <= 0) .or. (lnumait > pcnst)) then
1543 write(*,'(/a/)') &
1544 '*** modal_aero_newnuc bypass -- lnumait <= 0'
1545 return
1546 else if ((mait <= 0) .or. (mait > ntot_amode)) then
1547 write(*,'(/a/)') &
1548 '*** modal_aero_newnuc bypass -- modeptr_aitken <= 0'
1549 return
1550 end if
1552 l_h2so4_sv = l_h2so4
1553 l_nh3_sv = l_nh3
1554 lnumait_sv = lnumait
1555 lnh4ait_sv = lnh4ait
1556 lso4ait_sv = lso4ait
1558 !
1559 ! create history file column-tendency fields
1560 !
1561 dotend(:) = .false.
1562 dotend(lnumait) = .true.
1563 dotend(lso4ait) = .true.
1564 dotend(l_h2so4) = .true.
1565 if ((l_nh3 > 0) .and. (l_nh3 <= pcnst) .and. &
1566 (lnh4ait > 0) .and. (lnh4ait <= pcnst)) then
1567 dotend(lnh4ait) = .true.
1568 dotend(l_nh3) = .true.
1569 end if
1571 do l = 1, pcnst
1572 if ( .not. dotend(l) ) cycle
1573 tmpname = cnst_name(l)
1574 unit = 'kg/m2/s'
1575 do m = 1, ntot_amode
1576 if (l == numptr_amode(m)) unit = '#/m2/s'
1577 end do
1578 fieldname = trim(tmpname) // '_sfnnuc1'
1579 long_name = trim(tmpname) // ' modal_aero new particle nucleation column tendency'
1580 call addfld( fieldname, unit, 1, 'A', long_name, phys_decomp )
1581 if ( history_aerosol ) then
1582 call add_default( fieldname, 1, ' ' )
1583 endif
1584 if ( masterproc ) write(*,'(3(a,2x))') &
1585 'modal_aero_newnuc_init addfld', fieldname, unit
1586 end do ! l = ...
1589 return
1590 end subroutine modal_aero_newnuc_init
1594 !----------------------------------------------------------------------
1595 !----------------------------------------------------------------------
1596 subroutine ternary_nuc_merik2007( t, rh, c2, c3, j_log, ntot, nacid, namm, r )
1597 !subroutine ternary_fit( t, rh, c2, c3, j_log, ntot, nacid, namm, r )
1598 ! *************************** ternary_fit.f90 ********************************
1599 ! joonas merikanto, 2006
1600 !
1601 ! fortran 90 subroutine that calculates the parameterized composition
1602 ! and nucleation rate of critical clusters in h2o-h2so4-nh3 vapor
1603 !
1604 ! warning: the fit should not be used outside its limits of validity
1605 ! (limits indicated below)
1606 !
1607 ! in:
1608 ! t: temperature (k), limits 235-295 k
1609 ! rh: relative humidity as fraction (eg. 0.5=50%) limits 0.05-0.95
1610 ! c2: sulfuric acid concentration (molecules/cm3) limits 5x10^4 - 10^9 molecules/cm3
1611 ! c3: ammonia mixing ratio (ppt) limits 0.1 - 1000 ppt
1612 !
1613 ! out:
1614 ! j_log: logarithm of nucleation rate (1/(s cm3))
1615 ! ntot: total number of molecules in the critical cluster
1616 ! nacid: number of sulfuric acid molecules in the critical cluster
1617 ! namm: number of ammonia molecules in the critical cluster
1618 ! r: radius of the critical cluster (nm)
1619 ! ****************************************************************************
1620 implicit none
1622 real(r8), intent(in) :: t, rh, c2, c3
1623 real(r8), intent(out) :: j_log, ntot, nacid, namm, r
1624 real(r8) :: j, t_onset
1626 t_onset=143.6002929064716 + 1.0178856665693992*rh + &
1627 10.196398812974294*log(c2) - &
1628 0.1849879416839113*log(c2)**2 - 17.161783213150173*log(c3) + &
1629 (109.92469248546053*log(c3))/log(c2) + &
1630 0.7734119613144357*log(c2)*log(c3) - 0.15576469879527022*log(c3)**2
1632 if(t_onset.gt.t) then
1634 j_log=-12.861848898625231 + 4.905527742256349*c3 - 358.2337705052991*rh -&
1635 0.05463019231872484*c3*t + 4.8630382337426985*rh*t + &
1636 0.00020258394697064567*c3*t**2 - 0.02175548069741675*rh*t**2 - &
1637 2.502406532869512e-7*c3*t**3 + 0.00003212869941055865*rh*t**3 - &
1638 4.39129415725234e6/log(c2)**2 + (56383.93843154586*t)/log(c2)**2 -&
1639 (239.835990963361*t**2)/log(c2)**2 + &
1640 (0.33765136625580167*t**3)/log(c2)**2 - &
1641 (629.7882041830943*rh)/(c3**3*log(c2)) + &
1642 (7.772806552631709*rh*t)/(c3**3*log(c2)) - &
1643 (0.031974053936299256*rh*t**2)/(c3**3*log(c2)) + &
1644 (0.00004383764128775082*rh*t**3)/(c3**3*log(c2)) + &
1645 1200.472096232311*log(c2) - 17.37107890065621*t*log(c2) + &
1646 0.08170681335921742*t**2*log(c2) - 0.00012534476159729881*t**3*log(c2) - &
1647 14.833042158178936*log(c2)**2 + 0.2932631303555295*t*log(c2)**2 - &
1648 0.0016497524241142845*t**2*log(c2)**2 + &
1649 2.844074805239367e-6*t**3*log(c2)**2 - 231375.56676032578*log(c3) - &
1650 100.21645273730675*rh*log(c3) + 2919.2852552424706*t*log(c3) + &
1651 0.977886555834732*rh*t*log(c3) - 12.286497122264588*t**2*log(c3) - &
1652 0.0030511783284506377*rh*t**2*log(c3) + &
1653 0.017249301826661612*t**3*log(c3) + 2.967320346100855e-6*rh*t**3*log(c3) + &
1654 (2.360931724951942e6*log(c3))/log(c2) - &
1655 (29752.130254319443*t*log(c3))/log(c2) + &
1656 (125.04965118142027*t**2*log(c3))/log(c2) - &
1657 (0.1752996881934318*t**3*log(c3))/log(c2) + &
1658 5599.912337254629*log(c2)*log(c3) - 70.70896612937771*t*log(c2)*log(c3) + &
1659 0.2978801613269466*t**2*log(c2)*log(c3) - &
1660 0.00041866525019504*t**3*log(c2)*log(c3) + 75061.15281456841*log(c3)**2 - &
1661 931.8802278173565*t*log(c3)**2 + 3.863266220840964*t**2*log(c3)**2 - &
1662 0.005349472062284983*t**3*log(c3)**2 - &
1663 (732006.8180571689*log(c3)**2)/log(c2) + &
1664 (9100.06398573816*t*log(c3)**2)/log(c2) - &
1665 (37.771091915932004*t**2*log(c3)**2)/log(c2) + &
1666 (0.05235455395566905*t**3*log(c3)**2)/log(c2) - &
1667 1911.0303773001353*log(c2)*log(c3)**2 + &
1668 23.6903969622286*t*log(c2)*log(c3)**2 - &
1669 0.09807872005428583*t**2*log(c2)*log(c3)**2 + &
1670 0.00013564560238552576*t**3*log(c2)*log(c3)**2 - &
1671 3180.5610833308*log(c3)**3 + 39.08268568672095*t*log(c3)**3 - &
1672 0.16048521066690752*t**2*log(c3)**3 + &
1673 0.00022031380023793877*t**3*log(c3)**3 + &
1674 (40751.075322248245*log(c3)**3)/log(c2) - &
1675 (501.66977622013934*t*log(c3)**3)/log(c2) + &
1676 (2.063469732254135*t**2*log(c3)**3)/log(c2) - &
1677 (0.002836873785758324*t**3*log(c3)**3)/log(c2) + &
1678 2.792313345723013*log(c2)**2*log(c3)**3 - &
1679 0.03422552111802899*t*log(c2)**2*log(c3)**3 + &
1680 0.00014019195277521142*t**2*log(c2)**2*log(c3)**3 - &
1681 1.9201227328396297e-7*t**3*log(c2)**2*log(c3)**3 - &
1682 980.923146020468*log(rh) + 10.054155220444462*t*log(rh) - &
1683 0.03306644502023841*t**2*log(rh) + 0.000034274041225891804*t**3*log(rh) + &
1684 (16597.75554295064*log(rh))/log(c2) - &
1685 (175.2365504237746*t*log(rh))/log(c2) + &
1686 (0.6033215603167458*t**2*log(rh))/log(c2) - &
1687 (0.0006731787599587544*t**3*log(rh))/log(c2) - &
1688 89.38961120336789*log(c3)*log(rh) + 1.153344219304926*t*log(c3)*log(rh) - &
1689 0.004954549700267233*t**2*log(c3)*log(rh) + &
1690 7.096309866238719e-6*t**3*log(c3)*log(rh) + &
1691 3.1712136610383244*log(c3)**3*log(rh) - &
1692 0.037822330602328806*t*log(c3)**3*log(rh) + &
1693 0.0001500555743561457*t**2*log(c3)**3*log(rh) - &
1694 1.9828365865570703e-7*t**3*log(c3)**3*log(rh)
1696 j=exp(j_log)
1698 ntot=57.40091052369212 - 0.2996341884645408*t + &
1699 0.0007395477768531926*t**2 - &
1700 5.090604835032423*log(c2) + 0.011016634044531128*t*log(c2) + &
1701 0.06750032251225707*log(c2)**2 - 0.8102831333223962*log(c3) + &
1702 0.015905081275952426*t*log(c3) - 0.2044174683159531*log(c2)*log(c3) + &
1703 0.08918159167625832*log(c3)**2 - 0.0004969033586666147*t*log(c3)**2 + &
1704 0.005704394549007816*log(c3)**3 + 3.4098703903474368*log(j) - &
1705 0.014916956508210809*t*log(j) + 0.08459090011666293*log(c3)*log(j) - &
1706 0.00014800625143907616*t*log(c3)*log(j) + 0.00503804694656905*log(j)**2
1708 r=3.2888553966535506e-10 - 3.374171768439839e-12*t + &
1709 1.8347359507774313e-14*t**2 + 2.5419844298881856e-12*log(c2) - &
1710 9.498107643050827e-14*t*log(c2) + 7.446266520834559e-13*log(c2)**2 + &
1711 2.4303397746137294e-11*log(c3) + 1.589324325956633e-14*t*log(c3) - &
1712 2.034596219775266e-12*log(c2)*log(c3) - 5.59303954457172e-13*log(c3)**2 - &
1713 4.889507104645867e-16*t*log(c3)**2 + 1.3847024107506764e-13*log(c3)**3 + &
1714 4.141077193427042e-15*log(j) - 2.6813110884009767e-14*t*log(j) + &
1715 1.2879071621313094e-12*log(c3)*log(j) - &
1716 3.80352446061867e-15*t*log(c3)*log(j) - 1.8790172502456827e-14*log(j)**2
1718 nacid=-4.7154180661803595 + 0.13436423483953885*t - &
1719 0.00047184686478816176*t**2 - &
1720 2.564010713640308*log(c2) + 0.011353312899114723*t*log(c2) + &
1721 0.0010801941974317014*log(c2)**2 + 0.5171368624197119*log(c3) - &
1722 0.0027882479896204665*t*log(c3) + 0.8066971907026886*log(c3)**2 - &
1723 0.0031849094214409335*t*log(c3)**2 - 0.09951184152927882*log(c3)**3 + &
1724 0.00040072788891745513*t*log(c3)**3 + 1.3276469271073974*log(j) - &
1725 0.006167654171986281*t*log(j) - 0.11061390967822708*log(c3)*log(j) + &
1726 0.0004367575329273496*t*log(c3)*log(j) + 0.000916366357266258*log(j)**2
1728 namm=71.20073903979772 - 0.8409600103431923*t + &
1729 0.0024803006590334922*t**2 + &
1730 2.7798606841602607*log(c2) - 0.01475023348171676*t*log(c2) + &
1731 0.012264508212031405*log(c2)**2 - 2.009926050440182*log(c3) + &
1732 0.008689123511431527*t*log(c3) - 0.009141180198955415*log(c2)*log(c3) + &
1733 0.1374122553905617*log(c3)**2 - 0.0006253227821679215*t*log(c3)**2 + &
1734 0.00009377332742098946*log(c3)**3 + 0.5202974341687757*log(j) - &
1735 0.002419872323052805*t*log(j) + 0.07916392322884074*log(c3)*log(j) - &
1736 0.0003021586030317366*t*log(c3)*log(j) + 0.0046977006608603395*log(j)**2
1738 else
1739 ! nucleation rate less that 5e-6, setting j_log arbitrary small
1740 j_log=-300.
1741 end if
1743 return
1745 end subroutine ternary_nuc_merik2007
1749 !----------------------------------------------------------------------
1750 #endif ! (defined MODAL_AERO)
1751 end module modal_aero_newnuc