Line data Source code
1 : module ndrop
2 :
3 : !---------------------------------------------------------------------------------
4 : ! Purpose:
5 : ! CAM Interface for droplet activation by modal aerosols
6 : !
7 : ! ***N.B.*** This module is currently hardcoded to recognize only the modes that
8 : ! affect the climate calculation. This is implemented by using list
9 : ! index 0 in all the calls to rad_constituent interfaces.
10 : !---------------------------------------------------------------------------------
11 :
12 : use shr_kind_mod, only: r8 => shr_kind_r8, shr_kind_cs
13 : use ppgrid, only: pcols, pver
14 : use physconst, only: pi, rhoh2o, mwh2o, r_universal, rh2o, &
15 : gravit, latvap, cpair, rair
16 : use constituents, only: pcnst, cnst_get_ind, cnst_name, cnst_spec_class_gas, cnst_species_class
17 : use physics_types, only: physics_state, physics_ptend, physics_ptend_init
18 : use physics_buffer, only: physics_buffer_desc, pbuf_get_index, pbuf_get_field
19 :
20 : use wv_saturation, only: qsat
21 : use phys_control, only: phys_getopts
22 : use ref_pres, only: top_lev => trop_cloud_top_lev
23 : use shr_spfn_mod, only: erf => shr_spfn_erf
24 : use cam_history, only: addfld, add_default, horiz_only, fieldname_len, outfld
25 : use cam_abortutils, only: endrun
26 : use cam_logfile, only: iulog
27 :
28 : use aerosol_properties_mod, only: aerosol_properties
29 : use aerosol_state_mod, only: aerosol_state, ptr2d_t
30 :
31 : implicit none
32 : private
33 : save
34 :
35 : public ndrop_init, dropmixnuc, activate_aerosol
36 :
37 : ! mathematical constants
38 : real(r8), parameter :: zero = 0._r8
39 : real(r8), parameter :: third = 1._r8/3._r8
40 : real(r8), parameter :: twothird = 2._r8*third
41 : real(r8), parameter :: sixth = 1._r8/6._r8
42 : real(r8), parameter :: sq2 = sqrt(2._r8)
43 : real(r8), parameter :: sq2pi = sqrt(2._r8*pi)
44 : real(r8), parameter :: sqpi = sqrt(pi)
45 : real(r8), parameter :: surften = 0.076_r8
46 : real(r8), parameter :: tmelt = 273._r8
47 :
48 : real(r8) :: aten
49 :
50 : ! CCN diagnostic fields
51 : integer, parameter :: psat=6 ! number of supersaturations to calc ccn concentration
52 : real(r8), parameter :: supersat(psat)= & ! supersaturation (%) to determine ccn concentration
53 : (/ 0.02_r8, 0.05_r8, 0.1_r8, 0.2_r8, 0.5_r8, 1.0_r8 /)
54 : character(len=8) :: ccn_name(psat)= &
55 : (/'CCN1','CCN2','CCN3','CCN4','CCN5','CCN6'/)
56 :
57 : ! indices in state and pbuf structures
58 : integer :: numliq_idx = -1
59 : integer :: kvh_idx = -1
60 :
61 : logical :: history_aerosol ! Output the aerosol tendencies
62 : character(len=fieldname_len), allocatable :: fieldname(:) ! names for drop nuc tendency output fields
63 : character(len=fieldname_len), allocatable :: fieldname_cw(:) ! names for drop nuc tendency output fields
64 :
65 : ! Indices for aerosol species in the ptend%q array.
66 : integer, allocatable :: aer_cnst_idx(:,:)
67 :
68 : logical :: lq(pcnst) = .false. ! set flags true for constituents with non-zero tendencies
69 : ! in the ptend object
70 :
71 : !===============================================================================
72 : contains
73 : !===============================================================================
74 :
75 1536 : subroutine ndrop_init(aero_props)
76 :
77 : class(aerosol_properties), intent(in) :: aero_props
78 :
79 : integer :: l, m, mm
80 : integer :: idxtmp = -1
81 : character(len=32) :: tmpname
82 : character(len=32) :: tmpname_cw
83 : character(len=128) :: long_name
84 : character(len=8) :: unit
85 : logical :: history_amwg ! output the variables used by the AMWG diag package
86 :
87 : !-------------------------------------------------------------------------------
88 :
89 : ! get indices into state%q and pbuf structures
90 1536 : call cnst_get_ind('NUMLIQ', numliq_idx)
91 :
92 1536 : kvh_idx = pbuf_get_index('kvh')
93 :
94 1536 : aten = 2._r8*mwh2o*surften/(r_universal*tmelt*rhoh2o)
95 :
96 : allocate( &
97 0 : aer_cnst_idx(aero_props%nbins(),0:maxval(aero_props%nmasses())), &
98 0 : fieldname(aero_props%ncnst_tot()), &
99 21504 : fieldname_cw(aero_props%ncnst_tot()) )
100 :
101 : ! Add dropmixnuc tendencies for all modal aerosol species
102 :
103 : call phys_getopts(history_amwg_out = history_amwg, &
104 1536 : history_aerosol_out = history_aerosol)
105 :
106 :
107 9216 : do m = 1, aero_props%nbins()
108 87552 : do l = 0, aero_props%nmasses(m)
109 :
110 78336 : mm = aero_props%indexer(m,l)
111 :
112 78336 : unit = 'kg/m2/s'
113 78336 : if (l == 0) then ! number
114 7680 : unit = '#/m2/s'
115 : end if
116 :
117 78336 : if (l == 0) then ! number
118 7680 : call aero_props%num_names( m, tmpname, tmpname_cw)
119 : else
120 70656 : call aero_props%mmr_names( m,l, tmpname, tmpname_cw)
121 : end if
122 :
123 78336 : fieldname(mm) = trim(tmpname) // '_mixnuc1'
124 78336 : fieldname_cw(mm) = trim(tmpname_cw) // '_mixnuc1'
125 :
126 : ! To set tendencies in the ptend object need to get the constituent indices
127 : ! for the prognostic species
128 :
129 78336 : call cnst_get_ind(tmpname, idxtmp, abort=.false.)
130 78336 : aer_cnst_idx(m,l) = idxtmp
131 :
132 78336 : if (idxtmp>0) then
133 78336 : lq(idxtmp) = .true.
134 : end if
135 :
136 : ! Add tendency fields to the history only when prognostic MAM is enabled.
137 78336 : long_name = trim(tmpname) // ' dropmixnuc mixnuc column tendency'
138 78336 : call addfld(fieldname(mm), horiz_only, 'A', unit, long_name, sampled_on_subcycle=.true.)
139 :
140 78336 : long_name = trim(tmpname_cw) // ' dropmixnuc mixnuc column tendency'
141 78336 : call addfld(fieldname_cw(mm), horiz_only, 'A', unit, long_name, sampled_on_subcycle=.true.)
142 :
143 86016 : if (history_aerosol) then
144 0 : call add_default(fieldname(mm), 1, ' ')
145 0 : call add_default(fieldname_cw(mm), 1, ' ')
146 : end if
147 :
148 : end do
149 : end do
150 :
151 3072 : call addfld('CCN1',(/ 'lev' /), 'A','#/cm3','CCN concentration at S=0.02%', sampled_on_subcycle=.true.)
152 3072 : call addfld('CCN2',(/ 'lev' /), 'A','#/cm3','CCN concentration at S=0.05%', sampled_on_subcycle=.true.)
153 3072 : call addfld('CCN3',(/ 'lev' /), 'A','#/cm3','CCN concentration at S=0.1%', sampled_on_subcycle=.true.)
154 3072 : call addfld('CCN4',(/ 'lev' /), 'A','#/cm3','CCN concentration at S=0.2%', sampled_on_subcycle=.true.)
155 3072 : call addfld('CCN5',(/ 'lev' /), 'A','#/cm3','CCN concentration at S=0.5%', sampled_on_subcycle=.true.)
156 3072 : call addfld('CCN6',(/ 'lev' /), 'A','#/cm3','CCN concentration at S=1.0%', sampled_on_subcycle=.true.)
157 :
158 :
159 3072 : call addfld('WTKE', (/ 'lev' /), 'A', 'm/s', 'Standard deviation of updraft velocity', sampled_on_subcycle=.true.)
160 3072 : call addfld('NDROPMIX', (/ 'lev' /), 'A', '#/kg/s', 'Droplet number mixing', sampled_on_subcycle=.true.)
161 3072 : call addfld('NDROPSRC', (/ 'lev' /), 'A', '#/kg/s', 'Droplet number source', sampled_on_subcycle=.true.)
162 3072 : call addfld('NDROPSNK', (/ 'lev' /), 'A', '#/kg/s', 'Droplet number loss by microphysics', sampled_on_subcycle=.true.)
163 1536 : call addfld('NDROPCOL', horiz_only, 'A', '#/m2', 'Column droplet number', sampled_on_subcycle=.true.)
164 :
165 : ! set the add_default fields
166 1536 : if (history_amwg) then
167 1536 : call add_default('CCN3', 1, ' ')
168 : endif
169 :
170 3072 : end subroutine ndrop_init
171 :
172 : !===============================================================================
173 :
174 58302720 : subroutine dropmixnuc( aero_props, aero_state, &
175 : state, ptend, dtmicro, pbuf, wsub, wmixmin, &
176 241920 : cldn, cldo, cldliqf, tendnd, factnum)
177 :
178 : ! vertical diffusion and nucleation of cloud droplets
179 : ! assume cloud presence controlled by cloud fraction
180 : ! doesn't distinguish between warm, cold clouds
181 :
182 : ! arguments
183 : type(physics_state), target, intent(in) :: state
184 : type(physics_ptend), intent(out) :: ptend
185 : real(r8), intent(in) :: dtmicro ! time step for microphysics (s)
186 : real(r8), intent(in) :: wmixmin ! minimum turbulence vertical velocity (m/s)
187 :
188 : type(physics_buffer_desc), pointer :: pbuf(:)
189 :
190 : class(aerosol_properties), intent(in) :: aero_props
191 : class(aerosol_state), intent(in) :: aero_state
192 :
193 : ! arguments
194 : real(r8), intent(in) :: wsub(pcols,pver) ! subgrid vertical velocity
195 : real(r8), intent(in) :: cldn(pcols,pver) ! cloud fraction
196 : real(r8), intent(in) :: cldo(pcols,pver) ! cloud fraction on previous time step
197 : real(r8), intent(in) :: cldliqf(pcols,pver) ! liquid cloud fraction (liquid / (liquid + ice))
198 :
199 : ! output arguments
200 : real(r8), intent(out) :: tendnd(pcols,pver) ! change in droplet number concentration (#/kg/s)
201 : real(r8), intent(out) :: factnum(:,:,:) ! activation fraction for aerosol number
202 : !--------------------Local storage-------------------------------------
203 :
204 : integer :: lchnk ! chunk identifier
205 : integer :: ncol ! number of columns
206 : integer :: nbin ! number of modes/bins
207 : integer :: nele_tot ! total number of aerosol elements
208 :
209 241920 : real(r8), pointer :: ncldwtr(:,:) ! droplet number concentration (#/kg)
210 241920 : real(r8), pointer :: temp(:,:) ! temperature (K)
211 241920 : real(r8), pointer :: pmid(:,:) ! mid-level pressure (Pa)
212 241920 : real(r8), pointer :: pint(:,:) ! pressure at layer interfaces (Pa)
213 241920 : real(r8), pointer :: pdel(:,:) ! pressure thickess of layer (Pa)
214 241920 : real(r8), pointer :: rpdel(:,:) ! inverse of pressure thickess of layer (/Pa)
215 241920 : real(r8), pointer :: zm(:,:) ! geopotential height of level (m)
216 :
217 241920 : real(r8), pointer :: kvh(:,:) ! vertical diffusivity (m2/s)
218 :
219 241920 : type(ptr2d_t), allocatable :: raer(:) ! aerosol mass, number mixing ratios
220 241920 : type(ptr2d_t), allocatable :: qqcw(:)
221 : real(r8) :: raertend(pver) ! tendency of aerosol mass, number mixing ratios
222 : real(r8) :: qqcwtend(pver) ! tendency of cloudborne aerosol mass, number mixing ratios
223 :
224 : real(r8), parameter :: zkmin = 0.01_r8, zkmax = 100._r8
225 : integer :: i, k, l, m, mm, n
226 : integer :: km1, kp1
227 : integer :: nnew, nsav, ntemp
228 : integer :: lptr
229 : integer :: nsubmix, nsubmix_bnd
230 : integer, save :: count_submix(100)
231 : integer :: phase ! phase of aerosol
232 :
233 : real(r8) :: arg
234 : real(r8) :: dtinv
235 : real(r8) :: dtmin, tinv, dtt
236 : real(r8) :: lcldn(pcols,pver)
237 : real(r8) :: lcldo(pcols,pver)
238 :
239 : real(r8) :: zs(pver) ! inverse of distance between levels (m)
240 : real(r8) :: qcld(pver) ! cloud droplet number mixing ratio (#/kg)
241 : real(r8) :: qncld(pver) ! droplet number nucleated on cloud boundaries
242 : real(r8) :: srcn(pver) ! droplet source rate (/s)
243 : real(r8) :: cs(pcols,pver) ! air density (kg/m3)
244 : real(r8) :: csbot(pver) ! air density at bottom (interface) of layer (kg/m3)
245 : real(r8) :: csbot_cscen(pver) ! csbot(i)/cs(i,k)
246 : real(r8) :: dz(pcols,pver) ! geometric thickness of layers (m)
247 :
248 : real(r8) :: wtke(pcols,pver) ! turbulent vertical velocity at base of layer k (m/s)
249 : real(r8) :: wtke_cen(pcols,pver) ! turbulent vertical velocity at center of layer k (m/s)
250 : real(r8) :: wbar, wmix, wmin, wmax
251 :
252 : real(r8) :: zn(pver) ! g/pdel (m2/g) for layer
253 : real(r8) :: flxconv ! convergence of flux into lowest layer
254 :
255 : real(r8) :: wdiab ! diabatic vertical velocity
256 : real(r8) :: ekd(pver) ! diffusivity for droplets (m2/s)
257 : real(r8) :: ekk(0:pver) ! density*diffusivity for droplets (kg/m3 m2/s)
258 : real(r8) :: ekkp(pver) ! zn*zs*density*diffusivity
259 : real(r8) :: ekkm(pver) ! zn*zs*density*diffusivity
260 :
261 : real(r8) :: dum, dumc
262 : real(r8) :: tmpa
263 : real(r8) :: dact
264 : real(r8) :: fluxntot ! (#/cm2/s)
265 : real(r8) :: dtmix
266 : real(r8) :: alogarg
267 : real(r8) :: overlapp(pver), overlapm(pver) ! cloud overlap
268 :
269 : real(r8) :: nsource(pcols,pver) ! droplet number source (#/kg/s)
270 : real(r8) :: ndropmix(pcols,pver) ! droplet number mixing (#/kg/s)
271 : real(r8) :: ndropcol(pcols) ! column droplet number (#/m2)
272 : real(r8) :: cldo_tmp, cldn_tmp
273 : real(r8) :: tau_cld_regenerate
274 : real(r8) :: taumix_internal_pver_inv ! 1/(internal mixing time scale for k=pver) (1/s)
275 :
276 :
277 241920 : real(r8), allocatable :: nact(:,:) ! fractional aero. number activation rate (/s)
278 241920 : real(r8), allocatable :: mact(:,:) ! fractional aero. mass activation rate (/s)
279 :
280 241920 : real(r8), allocatable :: raercol(:,:,:) ! single column of aerosol mass, number mixing ratios
281 241920 : real(r8), allocatable :: raercol_cw(:,:,:) ! same as raercol but for cloud-borne phase
282 :
283 :
284 : real(r8) :: na(pcols), va(pcols), hy(pcols)
285 241920 : real(r8), allocatable :: naermod(:) ! (1/m3)
286 241920 : real(r8), allocatable :: hygro(:) ! hygroscopicity of aerosol mode
287 241920 : real(r8), allocatable :: vaerosol(:) ! interstit+activated aerosol volume conc (cm3/cm3)
288 :
289 : real(r8) :: source(pver)
290 :
291 241920 : real(r8), allocatable :: fn(:) ! activation fraction for aerosol number
292 241920 : real(r8), allocatable :: fm(:) ! activation fraction for aerosol mass
293 :
294 241920 : real(r8), allocatable :: fluxn(:) ! number activation fraction flux (cm/s)
295 241920 : real(r8), allocatable :: fluxm(:) ! mass activation fraction flux (cm/s)
296 : real(r8) :: flux_fullact(pver) ! 100% activation fraction flux (cm/s)
297 : ! note: activation fraction fluxes are defined as
298 : ! fluxn = [flux of activated aero. number into cloud (#/cm2/s)]
299 : ! / [aero. number conc. in updraft, just below cloudbase (#/cm3)]
300 :
301 :
302 241920 : real(r8), allocatable :: coltend(:,:) ! column tendency for diagnostic output
303 241920 : real(r8), allocatable :: coltend_cw(:,:) ! column tendency
304 : real(r8) :: ccn(pcols,pver,psat) ! number conc of aerosols activated at supersat
305 :
306 : !for gas species turbulent mixing
307 : real(r8), pointer :: rgas(:, :, :)
308 : real(r8), allocatable :: rgascol(:, :, :)
309 : real(r8), allocatable :: coltendgas(:)
310 : real(r8) :: zerogas(pver)
311 : character*200 fieldnamegas
312 :
313 : integer :: errnum
314 : character(len=shr_kind_cs) :: errstr
315 : !-------------------------------------------------------------------------------
316 :
317 241920 : lchnk = state%lchnk
318 241920 : ncol = state%ncol
319 241920 : nbin = aero_props%nbins()
320 241920 : nele_tot = aero_props%ncnst_tot()
321 :
322 241920 : ncldwtr => state%q(:,:,numliq_idx)
323 241920 : temp => state%t
324 241920 : pmid => state%pmid
325 241920 : pint => state%pint
326 241920 : pdel => state%pdel
327 241920 : rpdel => state%rpdel
328 241920 : zm => state%zm
329 :
330 241920 : call pbuf_get_field(pbuf, kvh_idx, kvh)
331 :
332 : ! Create the liquid weighted cloud fractions that were passsed in
333 : ! before. This doesn't seem like the best variable, since the cloud could
334 : ! have liquid condensate, but the part of it that is changing could be the
335 : ! ice portion; however, this is what was done before.
336 119460096 : lcldo(:ncol,:) = cldo(:ncol,:) * cldliqf(:ncol,:)
337 119460096 : lcldn(:ncol,:) = cldn(:ncol,:) * cldliqf(:ncol,:)
338 :
339 :
340 241920 : arg = 1.0_r8
341 241920 : if (abs(0.8427_r8 - erf(arg))/0.8427_r8 > 0.001_r8) then
342 0 : write(iulog,*) 'erf(1.0) = ',ERF(arg)
343 0 : call endrun('dropmixnuc: Error function error')
344 : endif
345 241920 : arg = 0.0_r8
346 241920 : if (erf(arg) /= 0.0_r8) then
347 0 : write(iulog,*) 'erf(0.0) = ',erf(arg)
348 0 : write(iulog,*) 'dropmixnuc: Error function error'
349 0 : call endrun('dropmixnuc: Error function error')
350 : endif
351 :
352 241920 : dtinv = 1._r8/dtmicro
353 :
354 : allocate( &
355 : nact(pver,nbin), &
356 : mact(pver,nbin), &
357 0 : raer(nele_tot), &
358 0 : qqcw(nele_tot), &
359 : raercol(pver,nele_tot,2), &
360 : raercol_cw(pver,nele_tot,2), &
361 : coltend(pcols,nele_tot), &
362 : coltend_cw(pcols,nele_tot), &
363 0 : naermod(nbin), &
364 0 : hygro(nbin), &
365 0 : vaerosol(nbin), &
366 0 : fn(nbin), &
367 0 : fm(nbin), &
368 0 : fluxn(nbin), &
369 5322240 : fluxm(nbin) )
370 :
371 : ! Init pointers to mode number and specie mass mixing ratios in
372 : ! intersitial and cloud borne phases.
373 241920 : call aero_state%get_states( aero_props, raer, qqcw )
374 :
375 659473920 : factnum = 0._r8
376 241920 : wtke = 0._r8
377 241920 : nsource = 0._r8
378 241920 : ndropmix = 0._r8
379 241920 : ndropcol = 0._r8
380 241920 : tendnd = 0._r8
381 :
382 : ! initialize aerosol tendencies
383 241920 : call physics_ptend_init(ptend, state%psetcols, 'ndrop', lq=lq)
384 :
385 : ! overall_main_i_loop
386 3725568 : do i = 1, ncol
387 :
388 111476736 : do k = top_lev, pver-1
389 111476736 : zs(k) = 1._r8/(zm(i,k) - zm(i,k+1))
390 : end do
391 3483648 : zs(pver) = zs(pver-1)
392 :
393 : ! load number nucleated into qcld on cloud boundaries
394 :
395 114960384 : do k = top_lev, pver
396 :
397 111476736 : qcld(k) = ncldwtr(i,k)
398 111476736 : qncld(k) = 0._r8
399 111476736 : srcn(k) = 0._r8
400 111476736 : cs(i,k) = pmid(i,k)/(rair*temp(i,k)) ! air density (kg/m3)
401 111476736 : dz(i,k) = 1._r8/(cs(i,k)*gravit*rpdel(i,k)) ! layer thickness in m
402 :
403 668860416 : do m = 1, nbin
404 557383680 : nact(k,m) = 0._r8
405 668860416 : mact(k,m) = 0._r8
406 : end do
407 :
408 111476736 : zn(k) = gravit*rpdel(i,k)
409 :
410 111476736 : if (k < pver) then
411 107993088 : ekd(k) = kvh(i,k+1)
412 107993088 : ekd(k) = max(ekd(k), zkmin)
413 107993088 : ekd(k) = min(ekd(k), zkmax)
414 107993088 : csbot(k) = 2.0_r8*pint(i,k+1)/(rair*(temp(i,k) + temp(i,k+1)))
415 107993088 : csbot_cscen(k) = csbot(k)/cs(i,k)
416 : else
417 3483648 : ekd(k) = 0._r8
418 3483648 : csbot(k) = cs(i,k)
419 3483648 : csbot_cscen(k) = 1.0_r8
420 : end if
421 :
422 : ! rce-comment - define wtke at layer centers for new-cloud activation
423 : ! and at layer boundaries for old-cloud activation
424 111476736 : wtke_cen(i,k) = wsub(i,k)
425 111476736 : wtke(i,k) = wsub(i,k)
426 111476736 : wtke_cen(i,k) = max(wtke_cen(i,k), wmixmin)
427 111476736 : wtke(i,k) = max(wtke(i,k), wmixmin)
428 :
429 114960384 : nsource(i,k) = 0._r8
430 :
431 : end do
432 :
433 181149696 : nsav = 1
434 181149696 : nnew = 2
435 :
436 181149696 : do mm = 1,nele_tot
437 5862979584 : raercol_cw(:,mm,nsav) = 0.0_r8
438 5862979584 : raercol(:,mm,nsav) = 0.0_r8
439 5862979584 : raercol_cw(top_lev:pver,mm,nsav) = qqcw(mm)%fld(i,top_lev:pver)
440 5866463232 : raercol(top_lev:pver,mm,nsav) = raer(mm)%fld(i,top_lev:pver)
441 : end do
442 :
443 : ! droplet nucleation/aerosol activation
444 :
445 : ! tau_cld_regenerate = time scale for regeneration of cloudy air
446 : ! by (horizontal) exchange with clear air
447 3483648 : tau_cld_regenerate = 3600.0_r8 * 3.0_r8
448 :
449 : ! k-loop for growing/shrinking cloud calcs .............................
450 : ! grow_shrink_main_k_loop: &
451 114960384 : do k = top_lev, pver
452 :
453 : ! This code was designed for liquid clouds, but the cloudbourne
454 : ! aerosol can be either from liquid or ice clouds. For the ice clouds,
455 : ! we do not do regeneration, but as cloud fraction decreases the
456 : ! aerosols should be returned interstitial. The lack of a liquid cloud
457 : ! should not mean that all of the aerosol is realease. Therefor a
458 : ! section has been added for shrinking ice clouds and checks were added
459 : ! to protect ice cloudbourne aerosols from being released when no
460 : ! liquid cloud is present.
461 :
462 : ! shrinking ice cloud ......................................................
463 111476736 : cldo_tmp = cldo(i,k) * (1._r8 - cldliqf(i,k))
464 111476736 : cldn_tmp = cldn(i,k) * (1._r8 - cldliqf(i,k))
465 :
466 111476736 : if (cldn_tmp < cldo_tmp) then
467 :
468 : ! convert activated aerosol to interstitial in decaying cloud
469 :
470 10256574 : dumc = (cldn_tmp - cldo_tmp)/cldo_tmp * (1._r8 - cldliqf(i,k))
471 533341848 : do mm = 1,nele_tot
472 523085274 : dact = raercol_cw(k,mm,nsav)*dumc
473 523085274 : raercol_cw(k,mm,nsav) = raercol_cw(k,mm,nsav) + dact ! cloud-borne aerosol
474 533341848 : raercol(k,mm,nsav) = raercol(k,mm,nsav) - dact
475 : end do
476 :
477 : end if
478 :
479 : ! shrinking liquid cloud ......................................................
480 : ! treat the reduction of cloud fraction from when cldn(i,k) < cldo(i,k)
481 : ! and also dissipate the portion of the cloud that will be regenerated
482 111476736 : cldo_tmp = lcldo(i,k)
483 111476736 : cldn_tmp = lcldn(i,k) * exp( -dtmicro/tau_cld_regenerate )
484 : ! alternate formulation
485 : ! cldn_tmp = cldn(i,k) * max( 0.0_r8, (1.0_r8-dtmicro/tau_cld_regenerate) )
486 :
487 : ! fraction is also provided.
488 111476736 : if (cldn_tmp < cldo_tmp) then
489 : ! droplet loss in decaying cloud
490 : !++ sungsup
491 13775976 : nsource(i,k) = nsource(i,k) + qcld(k)*(cldn_tmp - cldo_tmp)/cldo_tmp*cldliqf(i,k)*dtinv
492 13775976 : qcld(k) = qcld(k)*(1._r8 + (cldn_tmp - cldo_tmp)/cldo_tmp)
493 : !-- sungsup
494 :
495 : ! convert activated aerosol to interstitial in decaying cloud
496 :
497 13775976 : dumc = (cldn_tmp - cldo_tmp)/cldo_tmp * cldliqf(i,k)
498 716350752 : do mm = 1,nele_tot
499 702574776 : dact = raercol_cw(k,mm,nsav)*dumc
500 702574776 : raercol_cw(k,mm,nsav) = raercol_cw(k,mm,nsav) + dact ! cloud-borne aerosol
501 716350752 : raercol(k,mm,nsav) = raercol(k,mm,nsav) - dact
502 : end do
503 :
504 : end if
505 :
506 : ! growing liquid cloud ......................................................
507 : ! treat the increase of cloud fraction from when cldn(i,k) > cldo(i,k)
508 : ! and also regenerate part of the cloud
509 111476736 : cldo_tmp = cldn_tmp
510 111476736 : cldn_tmp = lcldn(i,k)
511 :
512 114960384 : if (cldn_tmp-cldo_tmp > 0.01_r8) then
513 :
514 : ! rce-comment - use wtke at layer centers for new-cloud activation
515 7507451 : wbar = wtke_cen(i,k)
516 7507451 : wmix = 0._r8
517 7507451 : wmin = 0._r8
518 7507451 : wmax = 10._r8
519 7507451 : wdiab = 0._r8
520 :
521 : ! load aerosol properties, assuming external mixtures
522 :
523 7507451 : phase = 1 ! interstitial
524 45044706 : do m = 1, nbin
525 : call aero_state%loadaer( aero_props, &
526 : i, i, k, &
527 : m, cs, phase, na, va, &
528 37537255 : hy, errnum, errstr)
529 37537255 : if (errnum/=0) then
530 0 : call endrun('dropmixnuc : '//trim(errstr))
531 : end if
532 37537255 : naermod(m) = na(i)
533 37537255 : vaerosol(m) = va(i)
534 45044706 : hygro(m) = hy(i)
535 : end do
536 :
537 : call activate_aerosol( &
538 : wbar, wmix, wdiab, wmin, wmax, &
539 7507451 : temp(i,k), cs(i,k), naermod, nbin, &
540 : vaerosol, hygro, aero_props, fn, fm, fluxn, &
541 15014902 : fluxm,flux_fullact(k))
542 :
543 45044706 : factnum(i,k,:) = fn
544 :
545 7507451 : dumc = (cldn_tmp - cldo_tmp)
546 :
547 45044706 : do m = 1, nbin
548 37537255 : mm = aero_props%indexer(m,0)
549 37537255 : dact = dumc*fn(m)*raer(mm)%fld(i,k) ! interstitial only
550 37537255 : qcld(k) = qcld(k) + dact
551 37537255 : nsource(i,k) = nsource(i,k) + dact*dtinv
552 37537255 : raercol_cw(k,mm,nsav) = raercol_cw(k,mm,nsav) + dact ! cloud-borne aerosol
553 37537255 : raercol(k,mm,nsav) = raercol(k,mm,nsav) - dact
554 37537255 : dum = dumc*fm(m)
555 390387452 : do l = 1,aero_props%nmasses(m)
556 345342746 : mm = aero_props%indexer(m,l)
557 345342746 : dact = dum*raer(mm)%fld(i,k) ! interstitial only
558 345342746 : raercol_cw(k,mm,nsav) = raercol_cw(k,mm,nsav) + dact ! cloud-borne aerosol
559 382880001 : raercol(k,mm,nsav) = raercol(k,mm,nsav) - dact
560 : enddo
561 : enddo
562 :
563 : endif
564 :
565 : enddo ! grow_shrink_main_k_loop
566 : ! end of k-loop for growing/shrinking cloud calcs ......................
567 :
568 : ! ......................................................................
569 : ! start of k-loop for calc of old cloud activation tendencies ..........
570 : !
571 : ! rce-comment
572 : ! changed this part of code to use current cloud fraction (cldn) exclusively
573 : ! consider case of cldo(:)=0, cldn(k)=1, cldn(k+1)=0
574 : ! previous code (which used cldo below here) would have no cloud-base activation
575 : ! into layer k. however, activated particles in k mix out to k+1,
576 : ! so they are incorrectly depleted with no replacement
577 :
578 : ! old_cloud_main_k_loop
579 114960384 : do k = top_lev, pver
580 111476736 : kp1 = min0(k+1, pver)
581 111476736 : taumix_internal_pver_inv = 0.0_r8
582 :
583 114960384 : if (lcldn(i,k) > 0.01_r8) then
584 :
585 10024290 : wdiab = 0._r8
586 10024290 : wmix = 0._r8 ! single updraft
587 10024290 : wbar = wtke(i,k) ! single updraft
588 10024290 : if (k == pver) wbar = wtke_cen(i,k) ! single updraft
589 10024290 : wmax = 10._r8
590 10024290 : wmin = 0._r8
591 :
592 10024290 : if (lcldn(i,k) - lcldn(i,kp1) > 0.01_r8 .or. k == pver) then
593 :
594 : ! cloud base
595 :
596 : ! ekd(k) = wtke(i,k)*dz(i,k)/sq2pi
597 : ! rce-comments
598 : ! first, should probably have 1/zs(k) here rather than dz(i,k) because
599 : ! the turbulent flux is proportional to ekd(k)*zs(k),
600 : ! while the dz(i,k) is used to get flux divergences
601 : ! and mixing ratio tendency/change
602 : ! second and more importantly, using a single updraft velocity here
603 : ! means having monodisperse turbulent updraft and downdrafts.
604 : ! The sq2pi factor assumes a normal draft spectrum.
605 : ! The fluxn/fluxm from activate must be consistent with the
606 : ! fluxes calculated in explmix.
607 4998466 : ekd(k) = wbar/zs(k)
608 :
609 4998466 : alogarg = max(1.e-20_r8, 1._r8/lcldn(i,k) - 1._r8)
610 4998466 : wmin = wbar + wmix*0.25_r8*sq2pi*log(alogarg)
611 4998466 : phase = 1 ! interstitial
612 :
613 29990796 : do m = 1, nbin
614 : ! rce-comment - use kp1 here as old-cloud activation involves
615 : ! aerosol from layer below
616 : call aero_state%loadaer( aero_props, &
617 : i, i, kp1, &
618 : m, cs, phase, na, va, &
619 24992330 : hy, errnum, errstr)
620 24992330 : if (errnum/=0) then
621 0 : call endrun('dropmixnuc : '//trim(errstr))
622 : end if
623 24992330 : naermod(m) = na(i)
624 24992330 : vaerosol(m) = va(i)
625 29990796 : hygro(m) = hy(i)
626 : end do
627 :
628 : call activate_aerosol( &
629 : wbar, wmix, wdiab, wmin, wmax, &
630 4998466 : temp(i,k), cs(i,k), naermod, nbin, &
631 : vaerosol, hygro, aero_props, fn, fm, fluxn, &
632 9996932 : fluxm, flux_fullact(k))
633 :
634 29990796 : factnum(i,k,:) = fn
635 :
636 4998466 : if (k < pver) then
637 4502329 : dumc = lcldn(i,k) - lcldn(i,kp1)
638 : else
639 496137 : dumc = lcldn(i,k)
640 : endif
641 :
642 4998466 : fluxntot = 0
643 :
644 : ! rce-comment 1
645 : ! flux of activated mass into layer k (in kg/m2/s)
646 : ! = "actmassflux" = dumc*fluxm*raercol(kp1,lmass)*csbot(k)
647 : ! source of activated mass (in kg/kg/s) = flux divergence
648 : ! = actmassflux/(cs(i,k)*dz(i,k))
649 : ! so need factor of csbot_cscen = csbot(k)/cs(i,k)
650 : ! dum=1./(dz(i,k))
651 4998466 : dum=csbot_cscen(k)/(dz(i,k))
652 :
653 : ! rce-comment 2
654 : ! code for k=pver was changed to use the following conceptual model
655 : ! in k=pver, there can be no cloud-base activation unless one considers
656 : ! a scenario such as the layer being partially cloudy,
657 : ! with clear air at bottom and cloudy air at top
658 : ! assume this scenario, and that the clear/cloudy portions mix with
659 : ! a timescale taumix_internal = dz(i,pver)/wtke_cen(i,pver)
660 : ! in the absence of other sources/sinks, qact (the activated particle
661 : ! mixratio) attains a steady state value given by
662 : ! qact_ss = fcloud*fact*qtot
663 : ! where fcloud is cloud fraction, fact is activation fraction,
664 : ! qtot=qact+qint, qint is interstitial particle mixratio
665 : ! the activation rate (from mixing within the layer) can now be
666 : ! written as
667 : ! d(qact)/dt = (qact_ss - qact)/taumix_internal
668 : ! = qtot*(fcloud*fact*wtke/dz) - qact*(wtke/dz)
669 : ! note that (fcloud*fact*wtke/dz) is equal to the nact/mact
670 : ! also, d(qact)/dt can be negative. in the code below
671 : ! it is forced to be >= 0
672 : !
673 : ! steve --
674 : ! you will likely want to change this. i did not really understand
675 : ! what was previously being done in k=pver
676 : ! in the cam3_5_3 code, wtke(i,pver) appears to be equal to the
677 : ! droplet deposition velocity which is quite small
678 : ! in the cam3_5_37 version, wtke is done differently and is much
679 : ! larger in k=pver, so the activation is stronger there
680 : !
681 4998466 : if (k == pver) then
682 496137 : taumix_internal_pver_inv = flux_fullact(k)/dz(i,k)
683 : end if
684 :
685 29990796 : do m = 1, nbin
686 24992330 : mm = aero_props%indexer(m,0)
687 24992330 : fluxn(m) = fluxn(m)*dumc
688 24992330 : fluxm(m) = fluxm(m)*dumc
689 24992330 : nact(k,m) = nact(k,m) + fluxn(m)*dum
690 24992330 : mact(k,m) = mact(k,m) + fluxm(m)*dum
691 29990796 : if (k < pver) then
692 : ! note that kp1 is used here
693 : fluxntot = fluxntot &
694 22511645 : + fluxn(m)*raercol(kp1,mm,nsav)*cs(i,k)
695 : else
696 4961370 : tmpa = raercol(kp1,mm,nsav)*fluxn(m) &
697 : + raercol_cw(kp1,mm,nsav)*(fluxn(m) &
698 7442055 : - taumix_internal_pver_inv*dz(i,k))
699 2480685 : fluxntot = fluxntot + max(0.0_r8, tmpa)*cs(i,k)
700 : end if
701 : end do
702 4998466 : srcn(k) = srcn(k) + fluxntot/(cs(i,k)*dz(i,k))
703 4998466 : nsource(i,k) = nsource(i,k) + fluxntot/(cs(i,k)*dz(i,k))
704 :
705 : endif ! (cldn(i,k) - cldn(i,kp1) > 0.01 .or. k == pver)
706 :
707 : else
708 :
709 : ! no liquid cloud
710 101452446 : nsource(i,k) = nsource(i,k) - qcld(k)*dtinv
711 101452446 : qcld(k) = 0
712 :
713 101452446 : if (cldn(i,k) < 0.01_r8) then
714 : ! no ice cloud either
715 :
716 : ! convert activated aerosol to interstitial in decaying cloud
717 :
718 4675332792 : do mm = 1,nele_tot
719 4585422546 : raercol(k,mm,nsav) = raercol(k,mm,nsav) + raercol_cw(k,mm,nsav) ! cloud-borne aerosol
720 4675332792 : raercol_cw(k,mm,nsav) = 0._r8
721 : end do
722 :
723 : end if
724 : end if
725 :
726 : end do ! old_cloud_main_k_loop
727 :
728 : ! switch nsav, nnew so that nnew is the updated aerosol
729 3483648 : ntemp = nsav
730 3483648 : nsav = nnew
731 3483648 : nnew = ntemp
732 :
733 : ! load new droplets in layers above, below clouds
734 :
735 3483648 : dtmin = dtmicro
736 3483648 : ekk(top_lev-1) = 0.0_r8
737 3483648 : ekk(pver) = 0.0_r8
738 111476736 : do k = top_lev, pver-1
739 : ! rce-comment -- ekd(k) is eddy-diffusivity at k/k+1 interface
740 : ! want ekk(k) = ekd(k) * (density at k/k+1 interface)
741 : ! so use pint(i,k+1) as pint is 1:pverp
742 : ! ekk(k)=ekd(k)*2.*pint(i,k)/(rair*(temp(i,k)+temp(i,k+1)))
743 : ! ekk(k)=ekd(k)*2.*pint(i,k+1)/(rair*(temp(i,k)+temp(i,k+1)))
744 111476736 : ekk(k) = ekd(k)*csbot(k)
745 : end do
746 :
747 114960384 : do k = top_lev, pver
748 111476736 : km1 = max0(k-1, top_lev)
749 111476736 : ekkp(k) = zn(k)*ekk(k)*zs(k)
750 111476736 : ekkm(k) = zn(k)*ekk(k-1)*zs(km1)
751 111476736 : tinv = ekkp(k) + ekkm(k)
752 :
753 : ! rce-comment -- tinv is the sum of all first-order-loss-rates
754 : ! for the layer. for most layers, the activation loss rate
755 : ! (for interstitial particles) is accounted for by the loss by
756 : ! turb-transfer to the layer above.
757 : ! k=pver is special, and the loss rate for activation within
758 : ! the layer must be added to tinv. if not, the time step
759 : ! can be too big, and explmix can produce negative values.
760 : ! the negative values are reset to zero, resulting in an
761 : ! artificial source.
762 111476736 : if (k == pver) tinv = tinv + taumix_internal_pver_inv
763 :
764 114960384 : if (tinv .gt. 1.e-6_r8) then
765 26201922 : dtt = 1._r8/tinv
766 26201922 : dtmin = min(dtmin, dtt)
767 : end if
768 : end do
769 :
770 3483648 : dtmix = 0.9_r8*dtmin
771 3483648 : nsubmix = int(dtmicro/dtmix) + 1
772 3483648 : if (nsubmix > 100) then
773 : nsubmix_bnd = 100
774 : else
775 : nsubmix_bnd = nsubmix
776 : end if
777 3483648 : count_submix(nsubmix_bnd) = count_submix(nsubmix_bnd) + 1
778 3483648 : dtmix = dtmicro/nsubmix
779 :
780 114960384 : do k = top_lev, pver
781 111476736 : kp1 = min(k+1, pver)
782 111476736 : km1 = max(k-1, top_lev)
783 : ! maximum overlap assumption
784 111476736 : if (cldn(i,kp1) > 1.e-10_r8) then
785 37085337 : overlapp(k) = min(cldn(i,k)/cldn(i,kp1), 1._r8)
786 : else
787 74391399 : overlapp(k) = 1._r8
788 : end if
789 114960384 : if (cldn(i,km1) > 1.e-10_r8) then
790 34280979 : overlapm(k) = min(cldn(i,k)/cldn(i,km1), 1._r8)
791 : else
792 77195757 : overlapm(k) = 1._r8
793 : end if
794 : end do
795 :
796 :
797 : ! rce-comment
798 : ! the activation source(k) = mact(k,m)*raercol(kp1,lmass)
799 : ! should not exceed the rate of transfer of unactivated particles
800 : ! from kp1 to k which = ekkp(k)*raercol(kp1,lmass)
801 : ! however it might if things are not "just right" in subr activate
802 : ! the following is a safety measure to avoid negatives in explmix
803 111476736 : do k = top_lev, pver-1
804 651442176 : do m = 1, nbin
805 539965440 : nact(k,m) = min( nact(k,m), ekkp(k) )
806 647958528 : mact(k,m) = min( mact(k,m), ekkp(k) )
807 : end do
808 : end do
809 :
810 :
811 : ! old_cloud_nsubmix_loop
812 12903360 : do n = 1, nsubmix
813 9419712 : qncld(:) = qcld(:)
814 : ! switch nsav, nnew so that nsav is the updated aerosol
815 9419712 : ntemp = nsav
816 9419712 : nsav = nnew
817 9419712 : nnew = ntemp
818 9419712 : srcn(:) = 0.0_r8
819 :
820 56518272 : do m = 1, nbin
821 47098560 : mm = aero_props%indexer(m,0)
822 :
823 : ! update droplet source
824 : ! rce-comment- activation source in layer k involves particles from k+1
825 : ! srcn(:)=srcn(:)+nact(:,m)*(raercol(:,mm,nsav))
826 1507153920 : srcn(top_lev:pver-1) = srcn(top_lev:pver-1) + nact(top_lev:pver-1,m)*(raercol(top_lev+1:pver,mm,nsav))
827 :
828 : ! rce-comment- new formulation for k=pver
829 : ! srcn( pver )=srcn( pver )+nact( pver ,m)*(raercol( pver,mm,nsav))
830 : tmpa = raercol(pver,mm,nsav)*nact(pver,m) &
831 47098560 : + raercol_cw(pver,mm,nsav)*(nact(pver,m) - taumix_internal_pver_inv)
832 56518272 : srcn(pver) = srcn(pver) + max(0.0_r8,tmpa)
833 : end do
834 : call explmix( &
835 : qcld, srcn, ekkp, ekkm, overlapp, &
836 : overlapm, qncld, zero, zero, pver, &
837 9419712 : dtmix, .false.)
838 :
839 : ! rce-comment
840 : ! the interstitial particle mixratio is different in clear/cloudy portions
841 : ! of a layer, and generally higher in the clear portion. (we have/had
842 : ! a method for diagnosing the the clear/cloudy mixratios.) the activation
843 : ! source terms involve clear air (from below) moving into cloudy air (above).
844 : ! in theory, the clear-portion mixratio should be used when calculating
845 : ! source terms
846 60001920 : do m = 1, nbin
847 47098560 : mm = aero_props%indexer(m,0)
848 : ! rce-comment - activation source in layer k involves particles from k+1
849 : ! source(:)= nact(:,m)*(raercol(:,mm,nsav))
850 1507153920 : source(top_lev:pver-1) = nact(top_lev:pver-1,m)*(raercol(top_lev+1:pver,mm,nsav))
851 : ! rce-comment - new formulation for k=pver
852 : ! source( pver )= nact( pver, m)*(raercol( pver,mm,nsav))
853 : tmpa = raercol(pver,mm,nsav)*nact(pver,m) &
854 47098560 : + raercol_cw(pver,mm,nsav)*(nact(pver,m) - taumix_internal_pver_inv)
855 47098560 : source(pver) = max(0.0_r8, tmpa)
856 47098560 : flxconv = 0._r8
857 :
858 : call explmix( &
859 : raercol_cw(:,mm,nnew), source, ekkp, ekkm, overlapp, &
860 : overlapm, raercol_cw(:,mm,nsav), zero, zero, pver, &
861 47098560 : dtmix, .false.)
862 :
863 : call explmix( &
864 : raercol(:,mm,nnew), source, ekkp, ekkm, overlapp, &
865 : overlapm, raercol(:,mm,nsav), zero, flxconv, pver, &
866 47098560 : dtmix, .true., raercol_cw(:,mm,nsav))
867 :
868 489825024 : do l = 1,aero_props%nmasses(m)
869 433306752 : mm = aero_props%indexer(m,l)
870 : ! rce-comment - activation source in layer k involves particles from k+1
871 : ! source(:)= mact(:,m)*(raercol(:,mm,nsav))
872 13865816064 : source(top_lev:pver-1) = mact(top_lev:pver-1,m)*(raercol(top_lev+1:pver,mm,nsav))
873 : ! rce-comment- new formulation for k=pver
874 : ! source( pver )= mact( pver ,m)*(raercol( pver,mm,nsav))
875 : tmpa = raercol(pver,mm,nsav)*mact(pver,m) &
876 433306752 : + raercol_cw(pver,mm,nsav)*(mact(pver,m) - taumix_internal_pver_inv)
877 433306752 : source(pver) = max(0.0_r8, tmpa)
878 433306752 : flxconv = 0._r8
879 :
880 : call explmix( &
881 : raercol_cw(:,mm,nnew), source, ekkp, ekkm, overlapp, &
882 : overlapm, raercol_cw(:,mm,nsav), zero, zero, pver, &
883 433306752 : dtmix, .false.)
884 :
885 : call explmix( &
886 : raercol(:,mm,nnew), source, ekkp, ekkm, overlapp, &
887 : overlapm, raercol(:,mm,nsav), zero, flxconv, pver, &
888 480405312 : dtmix, .true., raercol_cw(:,mm,nsav))
889 :
890 : end do
891 : end do
892 :
893 : end do ! old_cloud_nsubmix_loop
894 :
895 : ! evaporate particles again if no cloud (either ice or liquid)
896 :
897 114960384 : do k = top_lev, pver
898 114960384 : if (cldn(i,k) == 0._r8) then
899 : ! no ice or liquid cloud
900 75793578 : qcld(k)=0._r8
901 :
902 : ! convert activated aerosol to interstitial in decaying cloud
903 3941266056 : do mm = 1,nele_tot
904 3865472478 : raercol(k,mm,nnew) = raercol(k,mm,nnew) + raercol_cw(k,mm,nnew)
905 3941266056 : raercol_cw(k,mm,nnew) = 0._r8
906 : end do
907 :
908 : end if
909 : end do
910 :
911 : ! droplet number
912 :
913 3483648 : ndropcol(i) = 0._r8
914 114960384 : do k = top_lev, pver
915 111476736 : ndropmix(i,k) = (qcld(k) - ncldwtr(i,k))*dtinv - nsource(i,k)
916 111476736 : tendnd(i,k) = (max(qcld(k), 1.e-6_r8) - ncldwtr(i,k))*dtinv
917 114960384 : ndropcol(i) = ndropcol(i) + ncldwtr(i,k)*pdel(i,k)
918 : end do
919 3483648 : ndropcol(i) = ndropcol(i)/gravit
920 :
921 3483648 : raertend = 0._r8
922 3483648 : qqcwtend = 0._r8
923 :
924 21143808 : do m = 1, nbin
925 198567936 : do l = 0, aero_props%nmasses(m)
926 :
927 177666048 : mm = aero_props%indexer(m,l)
928 177666048 : lptr = aer_cnst_idx(m,l)
929 :
930 5862979584 : raertend(top_lev:pver) = (raercol(top_lev:pver,mm,nnew) - raer(mm)%fld(i,top_lev:pver))*dtinv
931 5862979584 : qqcwtend(top_lev:pver) = (raercol_cw(top_lev:pver,mm,nnew) - qqcw(mm)%fld(i,top_lev:pver))*dtinv
932 :
933 5862979584 : coltend(i,mm) = sum( pdel(i,:)*raertend )/gravit
934 5862979584 : coltend_cw(i,mm) = sum( pdel(i,:)*qqcwtend )/gravit
935 :
936 : ! check for advected aerosol constituents
937 177666048 : if (lptr>0) then ! advected aerosol parts
938 5862979584 : ptend%q(i,:,lptr) = 0.0_r8
939 5862979584 : ptend%q(i,top_lev:pver,lptr) = raertend(top_lev:pver) ! set tendencies for interstitial aerosol
940 : else
941 0 : raer(mm)%fld(i,:) = 0.0_r8
942 0 : raer(mm)%fld(i,top_lev:pver) = raercol(top_lev:pver,mm,nnew) ! update non-advected interstitial aerosol (pbuf)
943 : end if
944 :
945 5862979584 : qqcw(mm)%fld(i,:) = 0.0_r8
946 5880397824 : qqcw(mm)%fld(i,top_lev:pver) = raercol_cw(top_lev:pver,mm,nnew) ! update cloud-borne aerosol
947 :
948 : end do
949 : end do
950 :
951 : end do ! overall_main_i_loop
952 : ! end of main loop over i/longitude ....................................
953 :
954 241920 : call outfld('NDROPCOL', ndropcol, pcols, lchnk)
955 241920 : call outfld('NDROPSRC', nsource, pcols, lchnk)
956 241920 : call outfld('NDROPMIX', ndropmix, pcols, lchnk)
957 241920 : call outfld('WTKE ', wtke, pcols, lchnk)
958 :
959 241920 : call ccncalc(aero_state, aero_props, state, cs, ccn)
960 1693440 : do l = 1, psat
961 1693440 : call outfld(ccn_name(l), ccn(1,1,l), pcols, lchnk)
962 : enddo
963 :
964 : ! do column tendencies
965 1451520 : do m = 1, nbin
966 13789440 : do l = 0,aero_props%nmasses(m)
967 12337920 : mm = aero_props%indexer(m,l)
968 12337920 : call outfld(fieldname(mm), coltend(:,mm), pcols, lchnk)
969 13547520 : call outfld(fieldname_cw(mm), coltend_cw(:,mm), pcols, lchnk)
970 : end do
971 : end do
972 :
973 0 : deallocate( &
974 : nact, &
975 0 : mact, &
976 0 : raer, &
977 0 : qqcw, &
978 0 : raercol, &
979 0 : raercol_cw, &
980 0 : coltend, &
981 0 : coltend_cw, &
982 0 : naermod, &
983 0 : hygro, &
984 0 : vaerosol, &
985 0 : fn, &
986 0 : fm, &
987 0 : fluxn, &
988 241920 : fluxm )
989 :
990 483840 : end subroutine dropmixnuc
991 :
992 : !===============================================================================
993 :
994 1940460672 : subroutine explmix( q, src, ekkp, ekkm, overlapp, overlapm, &
995 480405312 : qold, surfrate, flxconv, pver, dt, is_unact, qactold )
996 :
997 : ! explicit integration of droplet/aerosol mixing
998 : ! with source due to activation/nucleation
999 :
1000 :
1001 : integer, intent(in) :: pver ! number of levels
1002 : real(r8), intent(out) :: q(pver) ! mixing ratio to be updated
1003 : real(r8), intent(in) :: qold(pver) ! mixing ratio from previous time step
1004 : real(r8), intent(in) :: src(pver) ! source due to activation/nucleation (/s)
1005 : real(r8), intent(in) :: ekkp(pver) ! zn*zs*density*diffusivity (kg/m3 m2/s) at interface
1006 : ! below layer k (k,k+1 interface)
1007 : real(r8), intent(in) :: ekkm(pver) ! zn*zs*density*diffusivity (kg/m3 m2/s) at interface
1008 : ! above layer k (k,k+1 interface)
1009 : real(r8), intent(in) :: overlapp(pver) ! cloud overlap below
1010 : real(r8), intent(in) :: overlapm(pver) ! cloud overlap above
1011 : real(r8), intent(in) :: surfrate ! surface exchange rate (/s)
1012 : real(r8), intent(in) :: flxconv ! convergence of flux from surface
1013 : real(r8), intent(in) :: dt ! time step (s)
1014 : logical, intent(in) :: is_unact ! true if this is an unactivated species
1015 : real(r8), intent(in),optional :: qactold(pver)
1016 : ! mixing ratio of ACTIVATED species from previous step
1017 : ! *** this should only be present
1018 : ! if the current species is unactivated number/sfc/mass
1019 :
1020 : integer k,kp1,km1
1021 :
1022 970230336 : if ( is_unact ) then
1023 : ! the qactold*(1-overlap) terms are resuspension of activated material
1024 15853375296 : do k=top_lev,pver
1025 15372969984 : kp1=min(k+1,pver)
1026 15372969984 : km1=max(k-1,top_lev)
1027 30745939968 : q(k) = qold(k) + dt*( - src(k) + ekkp(k)*(qold(kp1) - qold(k) + &
1028 15372969984 : qactold(kp1)*(1.0_r8-overlapp(k))) &
1029 15372969984 : + ekkm(k)*(qold(km1) - qold(k) + &
1030 61491879936 : qactold(km1)*(1.0_r8-overlapm(k))) )
1031 : ! force to non-negative
1032 : ! if(q(k)<-1.e-30)then
1033 : ! write(iulog,*)'q=',q(k),' in explmix'
1034 15853375296 : q(k)=max(q(k),0._r8)
1035 : ! endif
1036 : end do
1037 :
1038 : ! diffusion loss at base of lowest layer
1039 480405312 : q(pver)=q(pver)-surfrate*qold(pver)*dt+flxconv*dt
1040 : ! force to non-negative
1041 : ! if(q(pver)<-1.e-30)then
1042 : ! write(iulog,*)'q=',q(pver),' in explmix'
1043 480405312 : q(pver)=max(q(pver),0._r8)
1044 : ! endif
1045 : else
1046 16164225792 : do k=top_lev,pver
1047 15674400768 : kp1=min(k+1,pver)
1048 15674400768 : km1=max(k-1,top_lev)
1049 31348801536 : q(k) = qold(k) + dt*(src(k) + ekkp(k)*(overlapp(k)*qold(kp1)-qold(k)) + &
1050 47023202304 : ekkm(k)*(overlapm(k)*qold(km1)-qold(k)) )
1051 : ! force to non-negative
1052 : ! if(q(k)<-1.e-30)then
1053 : ! write(iulog,*)'q=',q(k),' in explmix'
1054 16164225792 : q(k)=max(q(k),0._r8)
1055 : ! endif
1056 : end do
1057 : ! diffusion loss at base of lowest layer
1058 489825024 : q(pver)=q(pver)-surfrate*qold(pver)*dt+flxconv*dt
1059 : ! force to non-negative
1060 : ! if(q(pver)<-1.e-30)then
1061 : ! write(iulog,*)'q=',q(pver),' in explmix'
1062 489825024 : q(pver)=max(q(pver),0._r8)
1063 :
1064 : end if
1065 :
1066 970230336 : end subroutine explmix
1067 :
1068 : !===============================================================================
1069 :
1070 14062087 : subroutine activate_aerosol(wbar, sigw, wdiab, wminf, wmaxf, tair, rhoair, &
1071 14062087 : na, nbins, volume, hygro, aero_props, &
1072 14062087 : fn, fm, fluxn, fluxm, flux_fullact, smax_prescribed, in_cloud_in, smax_f)
1073 :
1074 : ! calculates number, surface, and mass fraction of aerosols activated as CCN
1075 : ! calculates flux of cloud droplets, surface area, and aerosol mass into cloud
1076 : ! assumes an internal mixture within each of up to nbin multiple aerosol bins
1077 : ! a gaussiam spectrum of updrafts can be treated.
1078 :
1079 : ! mks units
1080 :
1081 : ! Abdul-Razzak and Ghan, A parameterization of aerosol activation.
1082 : ! 2. Multiple aerosol types. J. Geophys. Res., 105, 6837-6844.
1083 :
1084 : ! input
1085 :
1086 : real(r8), intent(in) :: wbar ! grid cell mean vertical velocity (m/s)
1087 : real(r8), intent(in) :: sigw ! subgrid standard deviation of vertical vel (m/s)
1088 : real(r8), intent(in) :: wdiab ! diabatic vertical velocity (0 if adiabatic)
1089 : real(r8), intent(in) :: wminf ! minimum updraft velocity for integration (m/s)
1090 : real(r8), intent(in) :: wmaxf ! maximum updraft velocity for integration (m/s)
1091 : real(r8), intent(in) :: tair ! air temperature (K)
1092 : real(r8), intent(in) :: rhoair ! air density (kg/m3)
1093 : real(r8), intent(in) :: na(:) ! aerosol number concentration (/m3)
1094 : integer, intent(in) :: nbins ! number of aerosol bins
1095 : real(r8), intent(in) :: volume(:) ! aerosol volume concentration (m3/m3)
1096 : real(r8), intent(in) :: hygro(:) ! hygroscopicity of aerosol mode
1097 :
1098 : class(aerosol_properties), intent(in) :: aero_props
1099 :
1100 : ! output
1101 :
1102 : real(r8), intent(out) :: fn(:) ! number fraction of aerosols activated
1103 : real(r8), intent(out) :: fm(:) ! mass fraction of aerosols activated
1104 : real(r8), intent(out) :: fluxn(:) ! flux of activated aerosol number fraction into cloud (cm/s)
1105 : real(r8), intent(out) :: fluxm(:) ! flux of activated aerosol mass fraction into cloud (cm/s)
1106 : real(r8), intent(out) :: flux_fullact ! flux of activated aerosol fraction assuming 100% activation (cm/s)
1107 : ! rce-comment
1108 : ! used for consistency check -- this should match (ekd(k)*zs(k))
1109 : ! also, fluxm/flux_fullact gives fraction of aerosol mass flux
1110 : ! that is activated
1111 :
1112 : ! optional
1113 : real(r8), optional, intent(in) :: smax_prescribed ! prescribed max. supersaturation for secondary activation
1114 : logical, optional, intent(in) :: in_cloud_in ! switch to modify calculations when above cloud base
1115 : real(r8), optional, intent(in) :: smax_f ! droplet and rain size distr factor in the smax calculation
1116 : ! used when in_cloud=.true.
1117 :
1118 : ! local
1119 :
1120 : integer, parameter:: nx=200
1121 : real(r8) integ,integf
1122 : real(r8), parameter :: p0 = 1013.25e2_r8 ! reference pressure (Pa)
1123 : real(r8) pres ! pressure (Pa)
1124 : real(r8) diff0,conduct0
1125 : real(r8) es ! saturation vapor pressure
1126 : real(r8) qs ! water vapor saturation mixing ratio
1127 : real(r8) dqsdt ! change in qs with temperature
1128 : real(r8) g ! thermodynamic function (m2/s)
1129 28124174 : real(r8) zeta(nbins), eta(nbins)
1130 : real(r8) alpha
1131 : real(r8) gamma
1132 : real(r8) beta
1133 : real(r8) sqrtg
1134 28124174 : real(r8) :: amcube(nbins) ! cube of dry bin radius (m)
1135 28124174 : real(r8) smc(nbins) ! critical supersaturation for number bin radius
1136 : real(r8) sumflx_fullact
1137 28124174 : real(r8) sumflxn(nbins)
1138 28124174 : real(r8) sumflxm(nbins)
1139 28124174 : real(r8) sumfn(nbins)
1140 28124174 : real(r8) sumfm(nbins)
1141 28124174 : real(r8) fnold(nbins) ! number fraction activated
1142 28124174 : real(r8) fmold(nbins) ! mass fraction activated
1143 : real(r8) wold,gold
1144 : real(r8) wmin,wmax,w,dw,dwmax,dwmin,wnuc,dwnew,wb
1145 : real(r8) dfmin,dfmax,fnew,fold,fnmin,fnbar,fmbar
1146 : real(r8) alw,sqrtalw
1147 : real(r8) smax
1148 : real(r8) z,z1,z2,wf1,wf2,zf1,zf2,gf1,gf2,gf
1149 14062087 : real(r8) etafactor1,etafactor2(nbins),etafactor2max
1150 : real(r8) grow
1151 : character(len=*), parameter :: subname='activate_aerosol'
1152 :
1153 : logical :: in_cloud
1154 : integer m,n
1155 : ! numerical integration parameters
1156 : real(r8), parameter :: eps=0.3_r8,fmax=0.99_r8,sds=3._r8
1157 :
1158 : real(r8), parameter :: namin=1.e6_r8 ! minimum aerosol number concentration (/m3)
1159 :
1160 : integer ndist(nx) ! accumulates frequency distribution of integration bins required
1161 : data ndist/nx*0/
1162 : save ndist
1163 :
1164 14062087 : if (present(in_cloud_in)) then
1165 0 : if (.not. present(smax_f)) call endrun(subname//' error: smax_f must be supplied when in_cloud is used')
1166 0 : in_cloud = in_cloud_in
1167 : else
1168 : in_cloud = .false.
1169 : end if
1170 :
1171 84372522 : fn(:)=0._r8
1172 84372522 : fm(:)=0._r8
1173 84372522 : fluxn(:)=0._r8
1174 84372522 : fluxm(:)=0._r8
1175 14062087 : flux_fullact=0._r8
1176 :
1177 14062087 : if(nbins.eq.1.and.na(1).lt.1.e-20_r8)return
1178 :
1179 14062087 : if(sigw.le.1.e-5_r8.and.wbar.le.0._r8)return
1180 :
1181 14062087 : if ( present( smax_prescribed ) ) then
1182 1299743 : if (smax_prescribed <= 0.0_r8) return
1183 : end if
1184 :
1185 14062087 : pres=rair*rhoair*tair
1186 14062087 : diff0=0.211e-4_r8*(p0/pres)*(tair/tmelt)**1.94_r8
1187 14062087 : conduct0=(5.69_r8+0.017_r8*(tair-tmelt))*4.186e2_r8*1.e-5_r8 ! convert to J/m/s/deg
1188 14062087 : call qsat(tair, pres, es, qs)
1189 14062087 : dqsdt=latvap/(rh2o*tair*tair)*qs
1190 14062087 : alpha=gravit*(latvap/(cpair*rh2o*tair*tair)-1._r8/(rair*tair))
1191 14062087 : gamma=(1.0_r8+latvap/cpair*dqsdt)/(rhoair*qs)
1192 14062087 : etafactor2max=1.e10_r8/(alpha*wmaxf)**1.5_r8 ! this should make eta big if na is very small.
1193 :
1194 : grow = 1._r8/(rhoh2o/(diff0*rhoair*qs) &
1195 14062087 : + latvap*rhoh2o/(conduct0*tair)*(latvap/(rh2o*tair) - 1._r8))
1196 14062087 : sqrtg = sqrt(grow)
1197 14062087 : beta = 2._r8*pi*rhoh2o*grow*gamma
1198 :
1199 84372522 : do m=1,nbins
1200 :
1201 84372522 : if(volume(m).gt.1.e-39_r8.and.na(m).gt.1.e-39_r8)then
1202 : ! number mode radius (m)
1203 69654993 : amcube(m)=aero_props%amcube(m, volume(m),na(m))
1204 : ! growth coefficent Abdul-Razzak & Ghan 1998 eqn 16
1205 : ! should depend on mean radius of mode to account for gas kinetic effects
1206 : ! see Fountoukis and Nenes, JGR2005 and Meskhidze et al., JGR2006
1207 : ! for approriate size to use for effective diffusivity.
1208 69654993 : etafactor2(m)=1._r8/(na(m)*beta*sqrtg)
1209 69654993 : if(hygro(m).gt.1.e-10_r8)then
1210 69654993 : smc(m)=2._r8*aten*sqrt(aten/(27._r8*hygro(m)*amcube(m))) ! only if variable size dist
1211 : else
1212 0 : smc(m)=100._r8
1213 : endif
1214 : else
1215 655442 : smc(m)=1._r8
1216 655442 : etafactor2(m)=etafactor2max ! this should make eta big if na is very small.
1217 : endif
1218 :
1219 : enddo
1220 :
1221 14062087 : if(sigw.gt.1.e-5_r8)then ! spectrum of updrafts
1222 :
1223 0 : wmax=min(wmaxf,wbar+sds*sigw)
1224 0 : wmin=max(wminf,-wdiab)
1225 0 : wmin=max(wmin,wbar-sds*sigw)
1226 0 : w=wmin
1227 0 : dwmax=eps*sigw
1228 0 : dw=dwmax
1229 0 : dfmax=0.2_r8
1230 0 : dfmin=0.1_r8
1231 0 : if (wmax <= w) return
1232 0 : do m=1,nbins
1233 0 : sumflxn(m)=0._r8
1234 0 : sumfn(m)=0._r8
1235 0 : fnold(m)=0._r8
1236 0 : sumflxm(m)=0._r8
1237 0 : sumfm(m)=0._r8
1238 0 : fmold(m)=0._r8
1239 : enddo
1240 0 : sumflx_fullact=0._r8
1241 :
1242 0 : fold=0._r8
1243 0 : wold=0._r8
1244 0 : gold=0._r8
1245 :
1246 0 : dwmin = min( dwmax, 0.01_r8 )
1247 0 : do n = 1, nx
1248 :
1249 0 : 100 wnuc=w+wdiab
1250 : ! write(iulog,*)'wnuc=',wnuc
1251 0 : alw=alpha*wnuc
1252 0 : sqrtalw=sqrt(alw)
1253 0 : etafactor1=alw*sqrtalw
1254 :
1255 0 : do m=1,nbins
1256 0 : eta(m)=etafactor1*etafactor2(m)
1257 0 : zeta(m)=twothird*sqrtalw*aten/sqrtg
1258 : enddo
1259 :
1260 0 : if ( present( smax_prescribed ) ) then
1261 0 : smax = smax_prescribed
1262 : else
1263 0 : smax = aero_props%maxsat(zeta,eta,smc)
1264 : endif
1265 :
1266 0 : call aero_props%actfracs( nbins, smc(nbins), smax, fnew, fm(nbins) )
1267 :
1268 0 : dwnew = dw
1269 0 : if(fnew-fold.gt.dfmax.and.n.gt.1)then
1270 : ! reduce updraft increment for greater accuracy in integration
1271 0 : if (dw .gt. 1.01_r8*dwmin) then
1272 0 : dw=0.7_r8*dw
1273 0 : dw=max(dw,dwmin)
1274 0 : w=wold+dw
1275 0 : go to 100
1276 : else
1277 : dwnew = dwmin
1278 : endif
1279 : endif
1280 :
1281 0 : if(fnew-fold.lt.dfmin)then
1282 : ! increase updraft increment to accelerate integration
1283 0 : dwnew=min(1.5_r8*dw,dwmax)
1284 : endif
1285 0 : fold=fnew
1286 :
1287 0 : z=(w-wbar)/(sigw*sq2)
1288 0 : g=exp(-z*z)
1289 0 : fnmin=1._r8
1290 :
1291 0 : do m=1,nbins
1292 : ! modal
1293 0 : call aero_props%actfracs( m, smc(m), smax, fn(m), fm(m) )
1294 0 : fnmin=min(fn(m),fnmin)
1295 : ! integration is second order accurate
1296 : ! assumes linear variation of f*g with w
1297 0 : fnbar=(fn(m)*g+fnold(m)*gold)
1298 0 : fmbar=(fm(m)*g+fmold(m)*gold)
1299 0 : wb=(w+wold)
1300 0 : if(w.gt.0._r8)then
1301 : sumflxn(m)=sumflxn(m)+sixth*(wb*fnbar &
1302 0 : +(fn(m)*g*w+fnold(m)*gold*wold))*dw
1303 : sumflxm(m)=sumflxm(m)+sixth*(wb*fmbar &
1304 0 : +(fm(m)*g*w+fmold(m)*gold*wold))*dw
1305 : endif
1306 0 : sumfn(m)=sumfn(m)+0.5_r8*fnbar*dw
1307 0 : fnold(m)=fn(m)
1308 0 : sumfm(m)=sumfm(m)+0.5_r8*fmbar*dw
1309 0 : fmold(m)=fm(m)
1310 : enddo
1311 : ! same form as sumflxm but replace the fm with 1.0
1312 : sumflx_fullact = sumflx_fullact &
1313 0 : + sixth*(wb*(g+gold) + (g*w+gold*wold))*dw
1314 0 : gold=g
1315 0 : wold=w
1316 0 : dw=dwnew
1317 0 : if (n > 1 .and. (w > wmax .or. fnmin > fmax)) exit
1318 0 : w=w+dw
1319 0 : if (n == nx) then
1320 0 : write(iulog,*)'do loop is too short in activate'
1321 0 : write(iulog,*)'wmin=',wmin,' w=',w,' wmax=',wmax,' dw=',dw
1322 0 : write(iulog,*)'wbar=',wbar,' sigw=',sigw,' wdiab=',wdiab
1323 0 : write(iulog,*)'wnuc=',wnuc
1324 0 : write(iulog,*)'na=',(na(m),m=1,nbins)
1325 0 : write(iulog,*)'fn=',(fn(m),m=1,nbins)
1326 : ! dump all subr parameters to allow testing with standalone code
1327 : ! (build a driver that will read input and call activate)
1328 0 : write(iulog,*)'wbar,sigw,wdiab,tair,rhoair,nbins='
1329 0 : write(iulog,*) wbar,sigw,wdiab,tair,rhoair,nbins
1330 0 : write(iulog,*)'na=',na
1331 0 : write(iulog,*)'volume=', (volume(m),m=1,nbins)
1332 0 : write(iulog,*)'hygro='
1333 0 : write(iulog,*) hygro
1334 0 : call endrun(subname)
1335 : end if
1336 :
1337 : enddo
1338 :
1339 0 : ndist(n)=ndist(n)+1
1340 0 : if(w.lt.wmaxf)then
1341 :
1342 : ! contribution from all updrafts stronger than wmax
1343 : ! assuming constant f (close to fmax)
1344 0 : wnuc=w+wdiab
1345 :
1346 0 : z1=(w-wbar)/(sigw*sq2)
1347 0 : z2=(wmaxf-wbar)/(sigw*sq2)
1348 0 : g=exp(-z1*z1)
1349 0 : integ=sigw*0.5_r8*sq2*sqpi*(erf(z2)-erf(z1))
1350 : ! consider only upward flow into cloud base when estimating flux
1351 0 : wf1=max(w,zero)
1352 0 : zf1=(wf1-wbar)/(sigw*sq2)
1353 0 : gf1=exp(-zf1*zf1)
1354 0 : wf2=max(wmaxf,zero)
1355 0 : zf2=(wf2-wbar)/(sigw*sq2)
1356 0 : gf2=exp(-zf2*zf2)
1357 0 : gf=(gf1-gf2)
1358 0 : integf=wbar*sigw*0.5_r8*sq2*sqpi*(erf(zf2)-erf(zf1))+sigw*sigw*gf
1359 :
1360 0 : do m=1,nbins
1361 0 : sumflxn(m)=sumflxn(m)+integf*fn(m)
1362 0 : sumfn(m)=sumfn(m)+fn(m)*integ
1363 0 : sumflxm(m)=sumflxm(m)+integf*fm(m)
1364 0 : sumfm(m)=sumfm(m)+fm(m)*integ
1365 : enddo
1366 : ! same form as sumflxm but replace the fm with 1.0
1367 0 : sumflx_fullact = sumflx_fullact + integf
1368 : ! sumg=sumg+integ
1369 : endif
1370 :
1371 :
1372 0 : do m=1,nbins
1373 0 : fn(m)=sumfn(m)/(sq2*sqpi*sigw)
1374 : ! fn(m)=sumfn(m)/(sumg)
1375 0 : if(fn(m).gt.1.01_r8)then
1376 0 : write(iulog,*)'fn=',fn(m),' > 1 in activate'
1377 0 : write(iulog,*)'w,m,na,amcube=',w,m,na(m),amcube(m)
1378 0 : write(iulog,*)'integ,sumfn,sigw=',integ,sumfn(m),sigw
1379 0 : call endrun('activate')
1380 : endif
1381 0 : fluxn(m)=sumflxn(m)/(sq2*sqpi*sigw)
1382 0 : fm(m)=sumfm(m)/(sq2*sqpi*sigw)
1383 : ! fm(m)=sumfm(m)/(sumg)
1384 0 : if(fm(m).gt.1.01_r8)then
1385 0 : write(iulog,*)'fm=',fm(m),' > 1 in activate'
1386 : endif
1387 0 : fluxm(m)=sumflxm(m)/(sq2*sqpi*sigw)
1388 : enddo
1389 : ! same form as fluxm
1390 0 : flux_fullact = sumflx_fullact/(sq2*sqpi*sigw)
1391 :
1392 : else
1393 :
1394 : ! single updraft
1395 14062087 : wnuc=wbar+wdiab
1396 :
1397 14062087 : if(wnuc.gt.0._r8)then
1398 :
1399 14062087 : w=wbar
1400 :
1401 14062087 : if(in_cloud) then
1402 :
1403 0 : if (smax_f > 0._r8) then
1404 0 : smax = alpha*w/(2.0_r8*pi*rhoh2o*grow*gamma*smax_f)
1405 : else
1406 0 : smax = 1.e-20_r8
1407 : end if
1408 :
1409 : else ! at cloud base
1410 14062087 : alw = alpha*wnuc
1411 14062087 : sqrtalw = sqrt(alw)
1412 14062087 : etafactor1 = alw*sqrtalw
1413 :
1414 84372522 : do m = 1, nbins
1415 70310435 : eta(m) = etafactor1*etafactor2(m)
1416 84372522 : zeta(m) = twothird*sqrtalw*aten/sqrtg
1417 : end do
1418 14062087 : if ( present(smax_prescribed) ) then
1419 1299743 : smax = smax_prescribed
1420 : else
1421 12762344 : smax = aero_props%maxsat(zeta,eta,smc)
1422 : end if
1423 : end if
1424 :
1425 84372522 : do m=1,nbins
1426 :
1427 70310435 : call aero_props%actfracs( m, smc(m), smax, fn(m), fm(m) )
1428 :
1429 84372522 : if(wbar.gt.0._r8)then
1430 70310435 : fluxn(m)=fn(m)*w
1431 70310435 : fluxm(m)=fm(m)*w
1432 : endif
1433 : enddo
1434 14062087 : flux_fullact = w
1435 : endif
1436 :
1437 : endif
1438 :
1439 : end subroutine activate_aerosol
1440 :
1441 : !===============================================================================
1442 :
1443 241920 : subroutine ccncalc(aero_state, aero_props, state, cs, ccn)
1444 :
1445 : ! calculates number concentration of aerosols activated as CCN at
1446 : ! supersaturation supersat.
1447 : ! assumes an internal mixture of a multiple externally-mixed aerosol modes
1448 : ! cgs units
1449 :
1450 : ! Ghan et al., Atmos. Res., 1993, 198-221.
1451 :
1452 : ! arguments
1453 : class(aerosol_state), intent(in) :: aero_state
1454 : class(aerosol_properties), intent(in) :: aero_props
1455 :
1456 : type(physics_state), target, intent(in) :: state
1457 :
1458 : real(r8), intent(in) :: cs(pcols,pver) ! air density (kg/m3)
1459 : real(r8), intent(out) :: ccn(pcols,pver,psat) ! number conc of aerosols activated at supersat (#/m3)
1460 :
1461 : ! local
1462 :
1463 : integer :: ncol ! number of columns
1464 : integer :: nbin ! number of bins
1465 241920 : real(r8), pointer :: tair(:,:) ! air temperature (K)
1466 :
1467 : real(r8) naerosol(pcols) ! interstit+activated aerosol number conc (/m3)
1468 : real(r8) vaerosol(pcols) ! interstit+activated aerosol volume conc (m3/m3)
1469 :
1470 : real(r8) amcube(pcols)
1471 241920 : real(r8), allocatable :: argfactor(:)
1472 : real(r8) surften_coef
1473 : real(r8) a(pcols) ! surface tension parameter
1474 : real(r8) hygro(pcols) ! aerosol hygroscopicity
1475 : real(r8) sm(pcols) ! critical supersaturation at mode radius
1476 : real(r8) arg(pcols)
1477 : integer l,m,i,k, astat
1478 : real(r8) smcoef(pcols)
1479 : integer phase ! phase of aerosol
1480 :
1481 : integer :: errnum
1482 : character(len=shr_kind_cs) :: errstr
1483 :
1484 : ! mathematical constants
1485 : real(r8), parameter :: super(psat) = supersat(:psat)*0.01_r8
1486 : real(r8), parameter :: smcoefcoef = 2._r8/sqrt(27._r8)
1487 :
1488 : !-------------------------------------------------------------------------------
1489 :
1490 483840 : nbin = aero_props%nbins()
1491 241920 : ncol = state%ncol
1492 241920 : tair => state%t
1493 :
1494 725760 : allocate( argfactor(nbin), stat=astat )
1495 241920 : if (astat/=0) then
1496 0 : call endrun('ndrop::ccncalc : not able to allocate argfactor')
1497 : end if
1498 :
1499 241920 : surften_coef=2._r8*mwh2o*surften/(r_universal*rhoh2o)
1500 :
1501 1451520 : do m=1,nbin
1502 1451520 : argfactor(m)=twothird/(sq2*aero_props%alogsig(m))
1503 : end do
1504 :
1505 241920 : ccn = 0._r8
1506 7983360 : do k=top_lev,pver
1507 :
1508 119218176 : do i=1,ncol
1509 111476736 : a(i)=surften_coef/tair(i,k)
1510 119218176 : smcoef(i)=smcoefcoef*a(i)*sqrt(a(i))
1511 : end do
1512 :
1513 46690560 : do m=1,nbin
1514 :
1515 38707200 : phase=3 ! interstitial+cloudborne
1516 :
1517 : call aero_state%loadaer( aero_props, &
1518 : 1, ncol, k, &
1519 : m, cs, phase, naerosol, vaerosol, &
1520 38707200 : hygro, errnum, errstr)
1521 38707200 : if (errnum/=0) then
1522 0 : call endrun('ccncalc : '//trim(errstr))
1523 : end if
1524 :
1525 2903040000 : where(naerosol(:ncol)>1.e-3_r8 .and. hygro(:ncol)>1.e-10_r8)
1526 38707200 : amcube(:ncol)=aero_props%amcube(m, vaerosol(:ncol), naerosol(:ncol) )
1527 38707200 : sm(:ncol)=smcoef(:ncol)/sqrt(hygro(:ncol)*amcube(:ncol)) ! critical supersaturation
1528 : elsewhere
1529 : sm(:ncol)=1._r8 ! value shouldn't matter much since naerosol is small
1530 : endwhere
1531 278691840 : do l=1,psat
1532 3615252480 : do i=1,ncol
1533 3344302080 : arg(i)=argfactor(m)*log(sm(i)/super(l))
1534 3576545280 : ccn(i,k,l)=ccn(i,k,l)+naerosol(i)*0.5_r8*(1._r8-erf(arg(i)))
1535 : enddo
1536 : enddo
1537 : enddo
1538 : enddo
1539 717002496 : ccn(:ncol,:,:)=ccn(:ncol,:,:)*1.e-6_r8 ! convert from #/m3 to #/cm3
1540 :
1541 241920 : deallocate( argfactor )
1542 :
1543 241920 : end subroutine ccncalc
1544 :
1545 : !===============================================================================
1546 : end module ndrop
|
