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