module_mixactivate.F
References to this file elsewhere.
1 !**********************************************************************************
2 ! This computer software was prepared by Battelle Memorial Institute, hereinafter
3 ! the Contractor, under Contract No. DE-AC05-76RL0 1830 with the Department of
4 ! Energy (DOE). NEITHER THE GOVERNMENT NOR THE CONTRACTOR MAKES ANY WARRANTY,
5 ! EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE.
6 !
7 ! MOSAIC module: see module_mosaic_driver.F for information and terms of use
8 !**********************************************************************************
9
10 MODULE module_mixactivate
11
12 CONTAINS
13
14
15 !----------------------------------------------------------------------
16 !----------------------------------------------------------------------
17 ! 06-nov-2005 rce - grid_id & ktau added to arg list
18 ! 25-apr-2006 rce - dens_aer is (g/cm3), NOT (kg/m3)
19 subroutine prescribe_aerosol_mixactivate ( &
20 grid_id, ktau, dtstep, naer, &
21 rho_phy, th_phy, pi_phy, w, cldfra, cldfra_old, &
22 z, dz8w, p_at_w, t_at_w, exch_h, &
23 qv, qc, qi, qndrop3d, &
24 nsource, &
25 ims,ime, jms,jme, kms,kme, &
26 its,ite, jts,jte, kts,kte, &
27 f_qc, f_qi )
28
29 ! USE module_configure
30
31 ! wrapper to call mixactivate for mosaic description of aerosol
32
33 implicit none
34
35 ! subr arguments
36 integer, intent(in) :: &
37 grid_id, ktau, &
38 ims, ime, jms, jme, kms, kme, &
39 its, ite, jts, jte, kts, kte
40
41 real, intent(in) :: dtstep
42 real, intent(inout) :: naer ! aerosol number (/kg)
43
44 real, intent(in), &
45 dimension( ims:ime, kms:kme, jms:jme ) :: &
46 rho_phy, th_phy, pi_phy, w, &
47 z, dz8w, p_at_w, t_at_w, exch_h
48
49 real, intent(inout), &
50 dimension( ims:ime, kms:kme, jms:jme ) :: cldfra, cldfra_old
51
52 real, intent(in), &
53 dimension( ims:ime, kms:kme, jms:jme ) :: &
54 qv, qc, qi
55
56 real, intent(inout), &
57 dimension( ims:ime, kms:kme, jms:jme ) :: &
58 qndrop3d
59
60 real, intent(out), &
61 dimension( ims:ime, kms:kme, jms:jme) :: nsource
62
63 LOGICAL, OPTIONAL :: f_qc, f_qi
64
65 ! local vars
66 integer maxd_aphase, maxd_atype, maxd_asize, maxd_acomp, max_chem
67 parameter (maxd_aphase=2,maxd_atype=1,maxd_asize=1,maxd_acomp=1, max_chem=10)
68 real ddvel(its:ite, jts:jte, max_chem) ! dry deposition velosity
69 real qsrflx(ims:ime, jms:jme, max_chem) ! dry deposition flux of aerosol
70 real chem(ims:ime, kms:kme, jms:jme, max_chem) ! chem array
71 integer i,j,k,l,m,n,p
72 real hygro( its:ite, kts:kte, jts:jte, maxd_asize, maxd_atype ) ! bulk
73 integer ntype_aer, nsize_aer(maxd_atype),ncomp_aer(maxd_atype), nphase_aer
74 integer massptr_aer( maxd_acomp, maxd_asize, maxd_atype, maxd_aphase ), &
75 waterptr_aer( maxd_asize, maxd_atype ), &
76 numptr_aer( maxd_asize, maxd_atype, maxd_aphase ), &
77 ai_phase, cw_phase
78 real dlo_sect( maxd_asize, maxd_atype ), & ! minimum size of section (cm)
79 dhi_sect( maxd_asize, maxd_atype ), & ! maximum size of section (cm)
80 sigmag_aer(maxd_asize, maxd_atype), & ! geometric standard deviation of aerosol size dist
81 dgnum_aer(maxd_asize, maxd_atype), & ! mean diameter (cm) of mode
82 dens_aer( maxd_acomp, maxd_atype), & ! density (g/cm3) of material
83 mw_aer( maxd_acomp, maxd_atype) ! molecular weight (g/mole)
84 real, dimension(ims:ime,kms:kme,jms:jme) :: &
85 ccn1,ccn2,ccn3,ccn4,ccn5,ccn6 ! number conc of aerosols activated at supersat
86 integer idrydep_onoff
87 real, dimension(ims:ime,kms:kme,jms:jme) :: t_phy
88 integer msectional
89
90
91 integer ptr
92 real maer
93
94 if(naer.lt.1.)then
95 naer=1000.e6 ! #/kg default value
96 endif
97 ai_phase=1
98 cw_phase=2
99 idrydep_onoff = 0
100 msectional = 0
101
102 t_phy(:,:,:)=th_phy(:,:,:)*pi_phy(:,:,:)
103
104 ntype_aer=maxd_atype
105 do n=1,ntype_aer
106 nsize_aer(n)=maxd_asize
107 ncomp_aer(n)=maxd_acomp
108 end do
109 nphase_aer=maxd_aphase
110
111 ! set properties for each type and size
112 do n=1,ntype_aer
113 do m=1,nsize_aer(n)
114 dlo_sect( m,n )=0.01e-4 ! minimum size of section (cm)
115 dhi_sect( m,n )=0.5e-4 ! maximum size of section (cm)
116 sigmag_aer(m,n)=2. ! geometric standard deviation of aerosol size dist
117 dgnum_aer(m,n)=0.1e-4 ! mean diameter (cm) of mode
118 end do
119 do l=1,ncomp_aer(n)
120 dens_aer( l, n)=1.0 ! density (g/cm3) of material
121 mw_aer( l, n)=132. ! molecular weight (g/mole)
122 end do
123 end do
124 ptr=0
125 do p=1,nphase_aer
126 do n=1,ntype_aer
127 do m=1,nsize_aer(n)
128 ptr=ptr+1
129 numptr_aer( m, n, p )=ptr
130 if(p.eq.ai_phase)then
131 chem(:,:,:,ptr)=naer
132 else
133 chem(:,:,:,ptr)=0
134 endif
135 end do ! size
136 end do ! type
137 end do ! phase
138 do p=1,maxd_aphase
139 do n=1,ntype_aer
140 do m=1,nsize_aer(n)
141 do l=1,ncomp_aer(n)
142 ptr=ptr+1
143 if(ptr.gt.max_chem)then
144 write(6,*)'ptr,max_chem=',ptr,max_chem,' in prescribe_aerosol_mixactivate'
145 call exit(1)
146 endif
147 massptr_aer(l, m, n, p)=ptr
148 ! maer is ug/kg-air; naer is #/kg-air; dgnum is cm; dens_aer is g/cm3
149 ! 1.e6 factor converts g to ug
150 maer= 1.0e6 * naer * dens_aer(l,n) * ( (3.1416/6.) * &
151 (dgnum_aer(m,n)**3) * exp( 4.5*((log(sigmag_aer(m,n)))**2) ) )
152 if(p.eq.ai_phase)then
153 chem(:,:,:,ptr)=maer
154 else
155 chem(:,:,:,ptr)=0
156 endif
157 end do
158 end do ! size
159 end do ! type
160 end do ! phase
161 do n=1,ntype_aer
162 do m=1,nsize_aer(n)
163 ptr=ptr+1
164 if(ptr.gt.max_chem)then
165 write(6,*)'ptr,max_chem=',ptr,max_chem,' in prescribe_aerosol_mixactivate'
166 call exit(1)
167 endif
168 !wig waterptr_aer(m, n)=ptr
169 waterptr_aer(m, n)=-1
170 end do ! size
171 end do ! type
172 ddvel(:,:,:)=0.
173 hygro( :,:,:,:,:) = 0.5
174
175 ! 06-nov-2005 rce - grid_id & ktau added to arg list
176 call mixactivate( msectional, &
177 chem,max_chem,qv,qc,qi,qndrop3d, &
178 t_phy, w, ddvel, idrydep_onoff, &
179 maxd_acomp, maxd_asize, maxd_atype, maxd_aphase, &
180 ncomp_aer, nsize_aer, ntype_aer, nphase_aer, &
181 numptr_aer, massptr_aer, dlo_sect, dhi_sect, sigmag_aer, dgnum_aer, &
182 dens_aer, mw_aer, &
183 waterptr_aer, hygro, ai_phase, cw_phase, &
184 ims,ime, jms,jme, &
185 kms,kme, &
186 its,ite, jts,jte, kts,kte, &
187 rho_phy, z, dz8w, p_at_w, t_at_w, exch_h, &
188 cldfra, cldfra_old, qsrflx, &
189 ccn1, ccn2, ccn3, ccn4, ccn5, ccn6, nsource, &
190 grid_id, ktau, dtstep, &
191 F_QC=f_qc, F_QI=f_qi )
192
193
194 end subroutine prescribe_aerosol_mixactivate
195
196 !----------------------------------------------------------------------
197 !----------------------------------------------------------------------
198 ! nov-04 sg ! replaced amode with aer and expanded aerosol dimension to include type and phase
199
200 ! 06-nov-2005 rce - grid_id & ktau added to arg list
201 ! 25-apr-2006 rce - dens_aer is (g/cm3), NOT (kg/m3)
202 subroutine mixactivate( msectional, &
203 chem, num_chem, qv, qc, qi, qndrop3d, &
204 temp, w, ddvel, idrydep_onoff, &
205 maxd_acomp, maxd_asize, maxd_atype, maxd_aphase, &
206 ncomp_aer, nsize_aer, ntype_aer, nphase_aer, &
207 numptr_aer, massptr_aer, dlo_sect, dhi_sect, sigmag_aer, dgnum_aer, &
208 dens_aer, mw_aer, &
209 waterptr_aer, hygro, ai_phase, cw_phase, &
210 ims,ime, jms,jme, kms,kme, &
211 its,ite, jts,jte, kts,kte, &
212 rho, zm, dz8w, p_at_w, t_at_w, kvh, &
213 cldfra, cldfra_old, qsrflx, &
214 ccn1, ccn2, ccn3, ccn4, ccn5, ccn6, nsource, &
215 grid_id, ktau, dtstep, &
216 f_qc, f_qi )
217
218
219 ! vertical diffusion and nucleation of cloud droplets
220 ! assume cloud presence controlled by cloud fraction
221 ! doesn't distinguish between warm, cold clouds
222
223 USE module_model_constants, only: g, rhowater, xlv, cp, rvovrd, r_d, r_v, mwdry, ep_2
224 USE module_radiation_driver, only: cal_cldfra
225
226 implicit none
227
228 ! input
229
230 INTEGER, intent(in) :: grid_id, ktau
231 INTEGER, intent(in) :: num_chem
232 integer, intent(in) :: ims,ime, jms,jme, kms,kme, &
233 its,ite, jts,jte, kts,kte
234
235 integer maxd_aphase, nphase_aer, maxd_atype, ntype_aer
236 integer maxd_asize, maxd_acomp, nsize_aer(maxd_atype)
237 integer, intent(in) :: &
238 ncomp_aer( maxd_atype ), &
239 massptr_aer( maxd_acomp, maxd_asize, maxd_atype, maxd_aphase ), &
240 waterptr_aer( maxd_asize, maxd_atype ), &
241 numptr_aer( maxd_asize, maxd_atype, maxd_aphase), &
242 ai_phase, cw_phase
243 integer, intent(in) :: msectional ! 1 for sectional, 0 for modal
244 integer, intent(in) :: idrydep_onoff
245 real, intent(in) :: &
246 dlo_sect( maxd_asize, maxd_atype ), & ! minimum size of section (cm)
247 dhi_sect( maxd_asize, maxd_atype ), & ! maximum size of section (cm)
248 sigmag_aer(maxd_asize, maxd_atype), & ! geometric standard deviation of aerosol size dist
249 dgnum_aer(maxd_asize, maxd_atype), & ! mean diameter (cm) of mode
250 dens_aer( maxd_acomp, maxd_atype), & ! density (g/cm3) of material
251 mw_aer( maxd_acomp, maxd_atype) ! molecular weight (g/mole)
252
253
254 REAL, intent(inout), DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ) :: &
255 chem ! aerosol molar mixing ratio (ug/kg or #/kg)
256
257 REAL, intent(in), DIMENSION( ims:ime, kms:kme, jms:jme ) :: &
258 qv, qc, qi ! water species (vapor, cloud drops, cloud ice) mixing ratio (g/g)
259
260 LOGICAL, OPTIONAL :: f_qc, f_qi
261
262 REAL, intent(inout), DIMENSION( ims:ime, kms:kme, jms:jme ) :: &
263 qndrop3d ! water species mixing ratio (g/g)
264
265 real, intent(in) :: dtstep ! time step for microphysics (s)
266 real, intent(in) :: temp(ims:ime, kms:kme, jms:jme) ! temperature (K)
267 real, intent(in) :: w(ims:ime, kms:kme, jms:jme) ! vertical velocity (m/s)
268 real, intent(in) :: rho(ims:ime, kms:kme, jms:jme) ! density at mid-level (kg/m3)
269 REAL, intent(in) :: ddvel( its:ite, jts:jte, num_chem ) ! deposition velocity (m/s)
270 real, intent(in) :: zm(ims:ime, kms:kme, jms:jme) ! geopotential height of level (m)
271 real, intent(in) :: dz8w(ims:ime, kms:kme, jms:jme) ! layer thickness (m)
272 real, intent(in) :: p_at_w(ims:ime, kms:kme, jms:jme) ! pressure at layer interface (Pa)
273 real, intent(in) :: t_at_w(ims:ime, kms:kme, jms:jme) ! temperature at layer interface (K)
274 real, intent(in) :: kvh(ims:ime, kms:kme, jms:jme) ! vertical diffusivity (m2/s)
275 real, intent(inout) :: cldfra_old(ims:ime, kms:kme, jms:jme)! cloud fraction on previous time step
276 real, intent(inout) :: cldfra(ims:ime, kms:kme, jms:jme) ! cloud fraction
277 real, intent(in) :: hygro( its:ite, kts:kte, jts:jte, maxd_asize, maxd_atype ) ! bulk hygroscopicity &
278
279 REAL, intent(out), DIMENSION( ims:ime, jms:jme, num_chem ) :: qsrflx ! dry deposition rate for aerosol
280 real, intent(out), dimension(ims:ime,kms:kme,jms:jme) :: nsource, & ! droplet number source (#/kg/s)
281 ccn1,ccn2,ccn3,ccn4,ccn5,ccn6 ! number conc of aerosols activated at supersat
282
283
284 !--------------------Local storage-------------------------------------
285 !
286 real :: qndrop(kms:kme) ! cloud droplet number mixing ratio (#/kg)
287 real :: lcldfra(kms:kme) ! liquid cloud fraction
288 real :: lcldfra_old(kms:kme) ! liquid cloud fraction for previous timestep
289 real :: wtke(kms:kme) ! turbulent vertical velocity at base of layer k (m2/s)
290 real zn(kms:kme) ! g/pdel (m2/g) for layer
291 real zs(kms:kme) ! inverse of distance between levels (m)
292 real zkmin,zkmax
293 data zkmin/0.01/,zkmax/100./
294 save zkmin,zkmax
295 real cs(kms:kme) ! air density (kg/m3)
296 real dz(kms:kme) ! geometric thickness of layers (m)
297
298 real wdiab ! diabatic vertical velocity
299 ! real, parameter :: wmixmin = 0.1 ! minimum turbulence vertical velocity (m/s)
300 real, parameter :: wmixmin = 0.2 ! minimum turbulence vertical velocity (m/s)
301 ! real, parameter :: wmixmin = 1.0 ! minimum turbulence vertical velocity (m/s)
302 real :: qndrop_new(kms:kme) ! droplet number nucleated on cloud boundaries
303 real :: ekd(kms:kme) ! diffusivity for droplets (m2/s)
304 real :: ekk(kms:kme) ! density*diffusivity for droplets (kg/m3 m2/s)
305 real :: srcn(kms:kme) ! droplet source rate (/s)
306 real sq2pi
307 data sq2pi/2.5066282746/
308 save sq2pi
309 real dtinv
310
311 logical top ! true if cloud top, false if cloud base or new cloud
312 logical first
313 save first
314 data first/.true./
315 integer km1,kp1
316 real wbar,wmix,wmin,wmax
317 real cmincld
318 data cmincld/1.e-12/
319 save cmincld
320 real dum,dumc
321 real dact
322 real fluxntot ! (#/cm2/s)
323 real fac_srflx
324 real depvel_drop
325 real :: surfrate(num_chem) ! surface exchange rate (/s)
326 real surfratemax ! max surfrate for all species treated here
327 real surfrate_drop ! surfade exchange rate for droplelts
328 real dtmin,tinv,dtt
329 integer nsubmix,nsubmix_bnd
330 integer i,j,k,m,n,nsub
331 real dtmix
332 real alogarg
333 real qcld
334 real pi
335 integer nnew,nsav,ntemp
336 real :: overlapp(kms:kme),overlapm(kms:kme) ! cloud overlap
337 real :: ekkp(kms:kme),ekkm(kms:kme) ! zn*zs*density*diffusivity
338 integer count_submix(100)
339 save count_submix
340
341 integer lnum,lnumcw,l,lmass,lmasscw,lsfc,lsfccw,ltype,lsig,lwater
342 integer :: ntype(maxd_asize)
343
344 real :: naerosol(maxd_asize, maxd_atype) ! interstitial aerosol number conc (/m3)
345 real :: naerosolcw(maxd_asize, maxd_atype) ! activated number conc (/m3)
346 real :: maerosol(maxd_acomp,maxd_asize, maxd_atype) ! interstit mass conc (kg/m3)
347 real :: maerosolcw(maxd_acomp,maxd_asize, maxd_atype) ! activated mass conc (kg/m3)
348 real :: maerosol_tot(maxd_asize, maxd_atype) ! species-total interstit mass conc (kg/m3)
349 real :: maerosol_totcw(maxd_asize, maxd_atype) ! species-total activated mass conc (kg/m3)
350 real :: vaerosol(maxd_asize, maxd_atype) ! interstit+activated aerosol volume conc (m3/m3)
351 real :: vaerosolcw(maxd_asize, maxd_atype) ! activated aerosol volume conc (m3/m3)
352 real :: raercol(kms:kme,num_chem,2) ! aerosol mass, number mixing ratios
353 real :: source(kms:kme) !
354
355 real :: fn(maxd_asize, maxd_atype) ! activation fraction for aerosol number
356 real :: fs(maxd_asize, maxd_atype) ! activation fraction for aerosol sfcarea
357 real :: fm(maxd_asize, maxd_atype) ! activation fraction for aerosol mass
358 integer :: ncomp(maxd_atype)
359
360 real :: fluxn(maxd_asize, maxd_atype) ! number activation fraction flux (m/s)
361 real :: fluxs(maxd_asize, maxd_atype) ! sfcarea activation fraction flux (m/s)
362 real :: fluxm(maxd_asize, maxd_atype) ! mass activation fraction flux (m/s)
363 ! note: activation fraction fluxes are defined as
364 ! fluxn = [flux of activated aero. number into cloud (#/cm2/s)]
365 ! / [aero. number conc. in updraft, just below cloudbase (#/cm3)]
366
367 real :: nact(kms:kme,maxd_asize, maxd_atype) ! fractional aero. number activation rate (/s)
368 real :: mact(kms:kme,maxd_asize, maxd_atype) ! fractional aero. mass activation rate (/s)
369 real :: npv(maxd_asize, maxd_atype) ! number per volume concentration (/m3)
370 real scale
371
372 real :: hygro_aer(maxd_asize, maxd_atype) ! hygroscopicity of aerosol mode
373 real :: exp45logsig ! exp(4.5*alogsig**2)
374 real :: alogsig(maxd_asize, maxd_atype) ! natl log of geometric standard dev of aerosol
375 integer psat
376 parameter (psat=6) ! number of supersaturations to calc ccn concentration
377 real ccn(kts:kte,psat) ! number conc of aerosols activated at supersat
378 real, parameter :: supersat(psat)= &! supersaturation (%) to determine ccn concentration
379 (/0.02,0.05,0.1,0.2,0.5,1.0/)
380 real super(psat) ! supersaturation
381 real surften ! surface tension of water w/respect to air (N/m)
382 data surften/0.076/
383 save surften
384 real :: ccnfact(psat,maxd_asize, maxd_atype)
385 real :: amcube(maxd_asize, maxd_atype) ! cube of dry mode radius (m)
386 real :: argfactor(maxd_asize, maxd_atype)
387 real aten ! surface tension parameter
388 real t0 ! reference temperature
389 real sm ! critical supersaturation
390 real arg
391
392 !!$#if (defined AIX)
393 !!$#define ERF erf
394 !!$#define ERFC erfc
395 !!$#else
396 !!$#define ERF erf
397 !!$ real erf
398 !!$#define ERFC erfc
399 !!$ real erfc
400 !!$#endif
401
402 character*8, parameter :: ccn_name(psat)=(/'CCN1','CCN2','CCN3','CCN4','CCN5','CCN6'/)
403 integer ids,ide, jds,jde, kds,kde
404
405 arg = 1.0
406 if (abs(0.8427-ERF_ALT(arg))/0.8427>0.001) then
407 write (6,*) 'erf_alt(1.0) = ',ERF_ALT(arg)
408 write (6,*) 'dropmixnuc: Error function error'
409 call exit
410 endif
411 arg = 0.0
412 if (ERF_ALT(arg) /= 0.0) then
413 write (6,*) 'erf_alt(0.0) = ',ERF_ALT(arg)
414 write (6,*) 'dropmixnuc: Error function error'
415 call exit
416 endif
417
418 pi = 4.*atan(1.0)
419 dtinv=1./dtstep
420
421 depvel_drop = 0.1 ! prescribed here rather than getting it from dry_dep_driver
422 if (idrydep_onoff .le. 0) depvel_drop = 0.0
423
424 do n=1,ntype_aer
425 do m=1,nsize_aer(n)
426 ncomp(n)=ncomp_aer(n)
427 ! print *,'sigmag_aer,dgnum_aer=',sigmag_aer(m,n),dgnum_aer(m,n)
428 alogsig(m,n)=alog(sigmag_aer(m,n))
429 ! used only if number is diagnosed from volume
430 npv(m,n)=6./(pi*(0.01*dgnum_aer(m,n))**3*exp(4.5*alogsig(m,n)*alogsig(m,n)))
431 end do
432 end do
433
434 t0=273.
435 aten=2.*surften/(r_v*t0*rhowater)
436 super(:)=0.01*supersat(:)
437 do n=1,ntype_aer
438 do m=1,nsize_aer(n)
439 exp45logsig=exp(4.5*alogsig(m,n)*alogsig(m,n))
440 argfactor(m,n)=2./(3.*sqrt(2.)*alogsig(m,n))
441 amcube(m,n)=3./(4.*pi*exp45logsig*npv(m,n))
442 enddo
443 enddo
444
445 IF( PRESENT(F_QC) .AND. PRESENT ( F_QI ) ) THEN
446 CALL cal_cldfra(CLDFRA,qc,qi,f_qc,f_qi, &
447 ids,ide, jds,jde, kds,kde, &
448 ims,ime, jms,jme, kms,kme, &
449 its,ite, jts,jte, kts,kte )
450 END IF
451
452 qsrflx(:,:,:) = 0.
453
454 ! start loop over columns
455
456 do 120 j=jts,jte
457 do 100 i=its,ite
458
459 raercol(:,:,:) = 0. !~ wig: added, but should not be necessary
460 fluxn(:,:) = 0. !~
461 fluxs(:,:) = 0. !~
462 fluxm(:,:) = 0. !~
463 fn(:,:) = 0. !~
464 fs(:,:) = 0. !~
465 fm(:,:) = 0. !~
466
467 ! load number nucleated into qndrop on cloud boundaries
468
469 ! initialization for current i .........................................
470
471 do k=kts+1,kte-1
472 zs(k)=1./(zm(i,k,j)-zm(i,k-1,j))
473 enddo
474 zs(kts)=zs(kts+1)
475 zs(kte)=0.
476
477 do k=kts,kte-1
478 !!$ if(qndrop3d(i,k,j).lt.-10.e6.or.qndrop3d(i,k,j).gt.1.E20)then
479 !!$! call exit(1)
480 !!$ endif
481 if(f_qi)then
482 qcld=qc(i,k,j)+qi(i,k,j)
483 else
484 qcld=qc(i,k,j)
485 endif
486 if(qcld.lt.-1..or.qcld.gt.1.)then
487 write(6,'(a,g12.2,a,3i5)')'qcld=',qcld,' for i,k,j=',i,k,j
488 call exit(1)
489 endif
490 if(qcld.gt.1.e-20)then
491 lcldfra(k)=cldfra(i,k,j)*qc(i,k,j)/qcld
492 lcldfra_old(k)=cldfra_old(i,k,j)*qc(i,k,j)/qcld
493 else
494 lcldfra(k)=0.
495 lcldfra_old(k)=0.
496 endif
497 qndrop(k)=qndrop3d(i,k,j)
498 ! qndrop(k)=1.e5
499 cs(k)=rho(i,k,j) ! air density (kg/m3)
500 dz(k)=dz8w(i,k,j)
501 do n=1,ntype_aer
502 do m=1,nsize_aer(n)
503 nact(k,m,n)=0.
504 mact(k,m,n)=0.
505 enddo
506 enddo
507 zn(k)=1./(cs(k)*dz(k))
508 if(k>kts)then
509 ekd(k)=kvh(i,k,j)
510 ekd(k)=max(ekd(k),zkmin)
511 ekd(k)=min(ekd(k),zkmax)
512 else
513 ekd(k)=0
514 endif
515 ! diagnose subgrid vertical velocity from diffusivity
516 if(k.eq.kts)then
517 wtke(k)=sq2pi*depvel_drop
518 ! wtke(k)=sq2pi*kvh(i,k,j)
519 ! wtke(k)=max(wtke(k),wmixmin)
520 else
521 wtke(k)=sq2pi*ekd(k)/dz(k)
522 endif
523 wtke(k)=max(wtke(k),wmixmin)
524 nsource(i,k,j)=0.
525 enddo
526
527 ! calculate surface rate and mass mixing ratio for aerosol
528
529 surfratemax = 0.0
530 nsav=1
531 nnew=2
532 surfrate_drop=depvel_drop/dz(kts)
533 surfratemax = max( surfratemax, surfrate_drop )
534 do n=1,ntype_aer
535 do m=1,nsize_aer(n)
536 lnum=numptr_aer(m,n,ai_phase)
537 lnumcw=numptr_aer(m,n,cw_phase)
538 if(lnum>0)then
539 surfrate(lnum)=ddvel(i,j,lnum)/dz(kts)
540 surfrate(lnumcw)=surfrate_drop
541 surfratemax = max( surfratemax, surfrate(lnum) )
542 ! scale = 1000./mwdry ! moles/kg
543 scale = 1.
544 raercol(kts:kte-1,lnumcw,nsav)=chem(i,kts:kte-1,j,lnumcw)*scale ! #/kg
545 raercol(kts:kte-1,lnum,nsav)=chem(i,kts:kte-1,j,lnum)*scale
546 endif
547 do l=1,ncomp(n)
548 lmass=massptr_aer(l,m,n,ai_phase)
549 lmasscw=massptr_aer(l,m,n,cw_phase)
550 ! scale = mw_aer(l,n)/mwdry
551 scale = 1.e-9 ! kg/ug
552 surfrate(lmass)=ddvel(i,j,lmass)/dz(kts)
553 surfrate(lmasscw)=surfrate_drop
554 surfratemax = max( surfratemax, surfrate(lmass) )
555 raercol(kts:kte-1,lmasscw,nsav)=chem(i,kts:kte-1,j,lmasscw)*scale ! kg/kg
556 raercol(kts:kte-1,lmass,nsav)=chem(i,kts:kte-1,j,lmass)*scale ! kg/kg
557 enddo
558 lwater=waterptr_aer(m,n)
559 if(lwater>0)then
560 surfrate(lwater)=ddvel(i,j,lwater)/dz(kts)
561 surfratemax = max( surfratemax, surfrate(lwater) )
562 raercol(kts:kte-1,lwater,nsav)=chem(i,kts:kte-1,j,lwater) ! don't bother to convert units,
563 ! because it doesn't contribute to aerosol mass
564 endif
565 enddo ! size
566 enddo ! type
567
568
569 ! droplet nucleation/aerosol activation
570
571 ! k-loop for growing/shrinking cloud calcs .............................
572
573 do k=kts,kte-1
574 km1=max0(k-1,1)
575 kp1=min0(k+1,kte-1)
576
577 if(lcldfra(k)-lcldfra_old(k).gt.0.01)then
578 ! go to 10
579
580 ! growing cloud
581
582 ! wmix=wtke(k)
583 wbar=w(i,k,j)+wtke(k)
584 wmix=0.
585 wmin=0.
586 ! 06-nov-2005 rce - increase wmax from 10 to 50 (deep convective clouds)
587 wmax=50.
588 wdiab=0
589
590 ! load aerosol properties, assuming external mixtures
591
592 do n=1,ntype_aer
593 do m=1,nsize_aer(n)
594 call loadaer(raercol(1,1,nsav),k,kms,kme,num_chem, &
595 cs(k), npv(m,n), dlo_sect(m,n),dhi_sect(m,n), &
596 maxd_acomp, ncomp(n), &
597 grid_id, ktau, i, j, m, n, &
598 numptr_aer(m,n,ai_phase),numptr_aer(m,n,cw_phase), &
599 dens_aer(1,n), &
600 massptr_aer(1,m,n,ai_phase), massptr_aer(1,m,n,cw_phase), &
601 maerosol(1,m,n), maerosolcw(1,m,n), &
602 maerosol_tot(m,n), maerosol_totcw(m,n), &
603 naerosol(m,n), naerosolcw(m,n), &
604 vaerosol(m,n), vaerosolcw(m,n) )
605
606 hygro_aer(m,n)=hygro(i,k,j,m,n)
607 enddo
608 enddo
609
610 ! 06-nov-2005 rce - grid_id & ktau added to arg list
611 call activate(wbar,wmix,wdiab,wmin,wmax,temp(i,k,j),cs(k), &
612 msectional, maxd_atype, ntype_aer, maxd_asize, nsize_aer, &
613 naerosol, vaerosol, &
614 dlo_sect,dhi_sect,sigmag_aer,hygro_aer, &
615 fn,fs,fm,fluxn,fluxs,fluxm, grid_id, ktau, i, j, k )
616
617 dumc=(lcldfra(k)-lcldfra_old(k))
618 do n=1,ntype_aer
619 do m=1,nsize_aer(n)
620 lnum=numptr_aer(m,n,ai_phase)
621 lnumcw=numptr_aer(m,n,cw_phase)
622 dact=dumc*fn(m,n)*(raercol(k,lnum,nsav)) ! interstitial only
623 ! print *,'fn=',fn(m,n),' for m,n=',m,n
624 ! print *,'growing cloud dumc=',dumc,' fn=',fn(m,n)
625 qndrop(k)=qndrop(k)+dact
626 nsource(i,k,j)=nsource(i,k,j)+dact*dtinv
627 if(lnum.gt.0)then
628 raercol(k,lnumcw,nsav) = raercol(k,lnumcw,nsav)+dact
629 raercol(k,lnum,nsav) = raercol(k,lnum,nsav)-dact
630 endif
631 do l=1,ncomp(n)
632 lmass=massptr_aer(l,m,n,ai_phase)
633 lmasscw=massptr_aer(l,m,n,cw_phase)
634 ! rce 07-jul-2005 - changed dact for mass to mimic that used for number
635 ! dact=dum*(raercol(k,lmass,nsav)) ! interstitial only
636 dact=dumc*fm(m,n)*(raercol(k,lmass,nsav)) ! interstitial only
637 raercol(k,lmasscw,nsav) = raercol(k,lmasscw,nsav)+dact
638 raercol(k,lmass,nsav) = raercol(k,lmass,nsav)-dact
639 enddo
640 enddo
641 enddo
642 ! 10 continue
643 endif
644
645 if(lcldfra(k) < lcldfra_old(k) .and. lcldfra_old(k) > 1.e-20)then
646 ! go to 20
647
648 ! shrinking cloud ......................................................
649
650 ! droplet loss in decaying cloud
651 nsource(i,k,j)=nsource(i,k,j)+qndrop(k)*(lcldfra(k)-lcldfra_old(k))*dtinv
652 qndrop(k)=qndrop(k)*(1.+lcldfra(k)-lcldfra_old(k))
653 ! convert activated aerosol to interstitial in decaying cloud
654
655 dumc=(lcldfra(k)-lcldfra_old(k))/lcldfra_old(k)
656 ! print *,'shrinking cloud dumc=',dumc
657 do n=1,ntype_aer
658 do m=1,nsize_aer(n)
659 lnum=numptr_aer(m,n,ai_phase)
660 lnumcw=numptr_aer(m,n,cw_phase)
661 if(lnum.gt.0)then
662 dact=raercol(k,lnumcw,nsav)*dumc
663 raercol(k,lnumcw,nsav)=raercol(k,lnumcw,nsav)+dact
664 raercol(k,lnum,nsav)=raercol(k,lnum,nsav)-dact
665 endif
666 do l=1,ncomp(n)
667 lmass=massptr_aer(l,m,n,ai_phase)
668 lmasscw=massptr_aer(l,m,n,cw_phase)
669 dact=raercol(k,lmasscw,nsav)*dumc
670 raercol(k,lmasscw,nsav)=raercol(k,lmasscw,nsav)+dact
671 raercol(k,lmass,nsav)=raercol(k,lmass,nsav)-dact
672 enddo
673 enddo
674 enddo
675 ! 20 continue
676 endif
677
678 enddo !k loop
679
680 ! end of k-loop for growing/shrinking cloud calcs ......................
681
682
683 ! ......................................................................
684 ! start of k-loop for calc of old cloud activation tendencies ..........
685
686 do k=kts,kte-1
687 km1=max0(k-1,kts)
688 kp1=min0(k+1,kte-1)
689 if(lcldfra(k).gt.0.01)then
690 if(lcldfra_old(k).gt.0.01)then
691 ! go to 30
692
693 ! old cloud
694
695 if(lcldfra_old(k)-lcldfra_old(km1).gt.0.01.or.k.eq.kts)then
696
697 ! interior cloud
698
699 ! cloud base
700
701 wdiab=0
702 wmix=wtke(k) ! spectrum of updrafts
703 wbar=w(i,k,j) ! spectrum of updrafts
704 ! wmix=0. ! single updraft
705 ! wbar=wtke(k) ! single updraft
706 ! 06-nov-2005 rce - increase wmax from 10 to 50 (deep convective clouds)
707 wmax=50.
708 top=.false.
709 ekd(k)=wtke(k)*dz(k)/sq2pi
710 alogarg=max(1.e-20,1/lcldfra_old(k)-1.)
711 wmin=wbar+wmix*0.25*sq2pi*alog(alogarg)
712
713 do n=1,ntype_aer
714 do m=1,nsize_aer(n)
715 call loadaer(raercol(1,1,nsav),km1,kms,kme,num_chem, &
716 cs(k), npv(m,n),dlo_sect(m,n),dhi_sect(m,n), &
717 maxd_acomp, ncomp(n), &
718 grid_id, ktau, i, j, m, n, &
719 numptr_aer(m,n,ai_phase),numptr_aer(m,n,cw_phase), &
720 dens_aer(1,n), &
721 massptr_aer(1,m,n,ai_phase), massptr_aer(1,m,n,cw_phase), &
722 maerosol(1,m,n), maerosolcw(1,m,n), &
723 maerosol_tot(m,n), maerosol_totcw(m,n), &
724 naerosol(m,n), naerosolcw(m,n), &
725 vaerosol(m,n), vaerosolcw(m,n) )
726
727 hygro_aer(m,n)=hygro(i,k,j,m,n)
728
729 enddo
730 enddo
731 ! print *,'old cloud wbar,wmix=',wbar,wmix
732
733 call activate(wbar,wmix,wdiab,wmin,wmax,temp(i,k,j),cs(k), &
734 msectional, maxd_atype, ntype_aer, maxd_asize, nsize_aer, &
735 naerosol, vaerosol, &
736 dlo_sect,dhi_sect, sigmag_aer,hygro_aer, &
737 fn,fs,fm,fluxn,fluxs,fluxm, grid_id, ktau, i, j, k )
738
739 if(k.gt.kts)then
740 dumc = lcldfra_old(k)-lcldfra_old(km1)
741 else
742 dumc=lcldfra_old(k)
743 endif
744 dum=1./(dz(k))
745 fluxntot=0.
746 do n=1,ntype_aer
747 do m=1,nsize_aer(n)
748 fluxn(m,n)=fluxn(m,n)*dumc
749 ! fluxs(m,n)=fluxs(m,n)*dumc
750 fluxm(m,n)=fluxm(m,n)*dumc
751 lnum=numptr_aer(m,n,ai_phase)
752 fluxntot=fluxntot+fluxn(m,n)*raercol(km1,lnum,nsav)
753 ! print *,'fn=',fn(m,n),' for m,n=',m,n
754 ! print *,'old cloud dumc=',dumc,' fn=',fn(m,n),' for m,n=',m,n
755 nact(k,m,n)=nact(k,m,n)+fluxn(m,n)*dum
756 mact(k,m,n)=mact(k,m,n)+fluxm(m,n)*dum
757 enddo
758 enddo
759 nsource(i,k,j)=nsource(i,k,j)+fluxntot*zs(k)
760 fluxntot=fluxntot*cs(k)
761 endif
762 ! 30 continue
763 endif
764 else
765 ! go to 40
766 ! no cloud
767 if(qndrop(k).gt.10000.e6)then
768 print *,'i,k,j,lcldfra,qndrop=',i,k,j,lcldfra(k),qndrop(k)
769 print *,'cldfra,ql,qi',cldfra(i,k,j),qc(i,k,j),qi(i,k,j)
770 endif
771 nsource(i,k,j)=nsource(i,k,j)-qndrop(k)*dtinv
772 qndrop(k)=0.
773 ! convert activated aerosol to interstitial in decaying cloud
774 do n=1,ntype_aer
775 do m=1,nsize_aer(n)
776 lnum=numptr_aer(m,n,ai_phase)
777 lnumcw=numptr_aer(m,n,cw_phase)
778 if(lnum.gt.0)then
779 raercol(k,lnum,nsav)=raercol(k,lnum,nsav)+raercol(k,lnumcw,nsav)
780 raercol(k,lnumcw,nsav)=0.
781 endif
782 do l=1,ncomp(n)
783 lmass=massptr_aer(l,m,n,ai_phase)
784 lmasscw=massptr_aer(l,m,n,cw_phase)
785 raercol(k,lmass,nsav)=raercol(k,lmass,nsav)+raercol(k,lmasscw,nsav)
786 raercol(k,lmasscw,nsav)=0.
787 enddo
788 enddo
789 enddo
790 ! 40 continue
791 endif
792 enddo
793 ! 50 continue
794
795 ! go to 100
796
797 ! switch nsav, nnew so that nnew is the updated aerosol
798
799 ntemp=nsav
800 nsav=nnew
801 nnew=ntemp
802
803 ! load new droplets in layers above, below clouds
804
805 dtmin=dtstep
806 ekk(kts)=0.0
807 do k=kts+1,kte-1
808 ekk(k)=ekd(k)*p_at_w(i,k,j)/(r_d*t_at_w(i,k,j))
809 enddo
810 ekk(kte)=0.0
811 do k=kts,kte-1
812 ekkp(k)=zn(k)*ekk(k+1)*zs(k+1)
813 ekkm(k)=zn(k)*ekk(k)*zs(k)
814 tinv=ekkp(k)+ekkm(k)
815 if(k.eq.kts)tinv=tinv+surfratemax
816 if(tinv.gt.1.e-6)then
817 dtt=1./tinv
818 dtmin=min(dtmin,dtt)
819 endif
820 enddo
821 dtmix=0.9*dtmin
822 nsubmix=dtstep/dtmix+1
823 if(nsubmix>100)then
824 nsubmix_bnd=100
825 else
826 nsubmix_bnd=nsubmix
827 endif
828 count_submix(nsubmix_bnd)=count_submix(nsubmix_bnd)+1
829 dtmix=dtstep/nsubmix
830 fac_srflx = -1.0/(zn(1)*nsubmix)
831
832 do k=kts,kte-1
833 kp1=min(k+1,kte-1)
834 km1=max(k-1,1)
835 if(lcldfra(kp1).gt.0)then
836 overlapp(k)=min(lcldfra(k)/lcldfra(kp1),1.)
837 else
838 overlapp(k)=1.
839 endif
840 if(lcldfra(km1).gt.0)then
841 overlapm(k)=min(lcldfra(k)/lcldfra(km1),1.)
842 else
843 overlapm(k)=1.
844 endif
845 enddo
846
847 do nsub=1,nsubmix
848 qndrop_new(kts:kte-1)=qndrop(kts:kte-1)
849 ! switch nsav, nnew so that nsav is the updated aerosol
850 ntemp=nsav
851 nsav=nnew
852 nnew=ntemp
853 srcn(:)=0.0
854 do n=1,ntype_aer
855 do m=1,nsize_aer(n)
856 lnum=numptr_aer(m,n,ai_phase)
857 ! update droplet source
858 srcn(kts:kte-1)=srcn(kts:kte-1)+nact(kts:kte-1,m,n)*(raercol(kts:kte-1,lnum,nsav))
859 enddo
860 enddo
861
862 call explmix(qndrop,srcn,ekkp,ekkm,overlapp,overlapm, &
863 qndrop_new,surfrate_drop,kts,kte-1,dtmix,.false.)
864 do n=1,ntype_aer
865 do m=1,nsize_aer(n)
866 lnum=numptr_aer(m,n,ai_phase)
867 lnumcw=numptr_aer(m,n,cw_phase)
868 if(lnum>0)then
869 source(kts:kte-1)= nact(kts:kte-1,m,n)*(raercol(kts:kte-1,lnum,nsav))
870 call explmix(raercol(1,lnumcw,nnew),source,ekkp,ekkm,overlapp,overlapm, &
871 raercol(1,lnumcw,nsav),surfrate(lnumcw),kts,kte-1,dtmix,&
872 .false.)
873 call explmix(raercol(1,lnum,nnew),source,ekkp,ekkm,overlapp,overlapm, &
874 raercol(1,lnum,nsav),surfrate(lnum),kts,kte-1,dtmix, &
875 .true.,raercol(1,lnumcw,nsav))
876 qsrflx(i,j,lnum) = qsrflx(i,j,lnum) + fac_srflx* &
877 raercol(kts,lnum,nsav)*surfrate(lnum)
878 qsrflx(i,j,lnumcw) = qsrflx(i,j,lnumcw) + fac_srflx* &
879 raercol(kts,lnumcw,nsav)*surfrate(lnumcw)
880 endif
881 do l=1,ncomp(n)
882 lmass=massptr_aer(l,m,n,ai_phase)
883 lmasscw=massptr_aer(l,m,n,cw_phase)
884 source(kts:kte-1)= mact(kts:kte-1,m,n)*(raercol(kts:kte-1,lmass,nsav))
885 call explmix(raercol(1,lmasscw,nnew),source,ekkp,ekkm,overlapp,overlapm, &
886 raercol(1,lmasscw,nsav),surfrate(lmasscw),kts,kte-1,dtmix, &
887 .false.)
888 call explmix(raercol(1,lmass,nnew),source,ekkp,ekkm,overlapp,overlapm, &
889 raercol(1,lmass,nsav),surfrate(lmass),kts,kte-1,dtmix, &
890 .true.,raercol(1,lmasscw,nsav))
891 qsrflx(i,j,lmass) = qsrflx(i,j,lmass) + fac_srflx* &
892 raercol(kts,lmass,nsav)*surfrate(lmass)
893 qsrflx(i,j,lmasscw) = qsrflx(i,j,lmasscw) + fac_srflx* &
894 raercol(kts,lmasscw,nsav)*surfrate(lmasscw)
895 enddo
896 lwater=waterptr_aer(m,n) ! aerosol water
897 if(lwater>0)then
898 source(:)=0.
899 call explmix( raercol(1,lwater,nnew),source,ekkp,ekkm,overlapp,overlapm, &
900 raercol(1,lwater,nsav),surfrate(lwater),kts,kte-1,dtmix, &
901 .true.,source)
902 endif
903 enddo ! size
904 enddo ! type
905
906 enddo !nsub
907
908 ! go to 100
909
910 ! evaporate particles again if no cloud
911
912 do k=kts,kte-1
913 if(lcldfra(k).eq.0.)then
914
915 ! no cloud
916
917 qndrop(k)=0.
918 ! convert activated aerosol to interstitial in decaying cloud
919 do n=1,ntype_aer
920 do m=1,nsize_aer(n)
921 lnum=numptr_aer(m,n,ai_phase)
922 lnumcw=numptr_aer(m,n,cw_phase)
923 if(lnum.gt.0)then
924 raercol(k,lnum,nnew)=raercol(k,lnum,nnew)+raercol(k,lnumcw,nnew)
925 raercol(k,lnumcw,nnew)=0.
926 endif
927 do l=1,ncomp(n)
928 lmass=massptr_aer(l,m,n,ai_phase)
929 lmasscw=massptr_aer(l,m,n,cw_phase)
930 raercol(k,lmass,nnew)=raercol(k,lmass,nnew)+raercol(k,lmasscw,nnew)
931 raercol(k,lmasscw,nnew)=0.
932 enddo
933 enddo
934 enddo
935 endif
936 enddo
937
938 ! go to 100
939 ! droplet number
940
941 do k=kts,kte-1
942 ! if(lcldfra(k).gt.0.1)then
943 ! write(6,'(a,3i5,f12.1)')'i,j,k,qndrop=',i,j,k,qndrop(k)
944 ! endif
945 if(qndrop(k).lt.-10.e6.or.qndrop(k).gt.1.e12)then
946 write(6,'(a,g12.2,a,3i5)')'after qndrop=',qndrop(k),' for i,k,j=',i,k,j
947 ! call exit(1)
948 endif
949
950 qndrop3d(i,k,j) = max(qndrop(k),1.e-6)
951
952 if(qndrop3d(i,k,j).lt.-10.e6.or.qndrop3d(i,k,j).gt.1.E20)then
953 write(6,'(a,g12.2,a,3i5)')'after qndrop=',qndrop3d(i,k,j),' for i,k,j=',i,k,j
954 ! call exit(1)
955 endif
956 if(qc(i,k,j).lt.-1..or.qc(i,k,j).gt.1.)then
957 write(6,'(a,g12.2,a,3i5)')'qc=',qc(i,k,j),' for i,k,j=',i,k,j
958 call exit(1)
959 endif
960 if(qi(i,k,j).lt.-1..or.qi(i,k,j).gt.1.)then
961 write(6,'(a,g12.2,a,3i5)')'qi=',qi(i,k,j),' for i,k,j=',i,k,j
962 call exit(1)
963 endif
964 if(qv(i,k,j).lt.-1..or.qv(i,k,j).gt.1.)then
965 write(6,'(a,g12.2,a,3i5)')'qv=',qv(i,k,j),' for i,k,j=',i,k,j
966 call exit(1)
967 endif
968 cldfra_old(i,k,j) = cldfra(i,k,j)
969 ! if(k.gt.6.and.k.lt.11)cldfra_old(i,k,j)=1.
970 enddo
971
972
973
974 ! go to 100
975 ! update chem and convert back to mole/mole
976
977 ccn(:,:) = 0.
978 do n=1,ntype_aer
979 do m=1,nsize_aer(n)
980 lnum=numptr_aer(m,n,ai_phase)
981 lnumcw=numptr_aer(m,n,cw_phase)
982 if(lnum.gt.0)then
983 ! scale=mwdry*0.001
984 scale = 1.
985 chem(i,kts:kte-1,j,lnumcw)= raercol(kts:kte-1,lnumcw,nnew)*scale
986 chem(i,kts:kte-1,j,lnum)= raercol(kts:kte-1,lnum,nnew)*scale
987 endif
988 do l=1,ncomp(n)
989 lmass=massptr_aer(l,m,n,ai_phase)
990 lmasscw=massptr_aer(l,m,n,cw_phase)
991 ! scale = mwdry/mw_aer(l,n)
992 scale = 1.e9
993 chem(i,kts:kte-1,j,lmasscw)=raercol(kts:kte-1,lmasscw,nnew)*scale ! ug/kg
994 chem(i,kts:kte-1,j,lmass)=raercol(kts:kte-1,lmass,nnew)*scale ! ug/kg
995 enddo
996 lwater=waterptr_aer(m,n)
997 if(lwater>0)chem(i,kts:kte-1,j,lwater)=raercol(kts:kte-1,lwater,nnew) ! don't convert units
998 do k=kts,kte-1
999 sm=2.*aten*sqrt(aten/(27.*hygro(i,k,j,m,n)*amcube(m,n)))
1000 do l=1,psat
1001 arg=argfactor(m,n)*log(sm/super(l))
1002 if(arg<2)then
1003 if(arg<-2)then
1004 ccnfact(l,m,n)=1.e-6 ! convert from #/m3 to #/cm3
1005 else
1006 ccnfact(l,m,n)=1.e-6*0.5*ERFC_NUM_RECIPES(arg)
1007 endif
1008 else
1009 ccnfact(l,m,n) = 0.
1010 endif
1011 ! ccn concentration as diagnostic
1012 ! assume same hygroscopicity and ccnfact for cloud-phase and aerosol phase particles
1013 ccn(k,l)=ccn(k,l)+(raercol(k,lnum,nnew)+raercol(k,lnumcw,nnew))*cs(k)*ccnfact(l,m,n)
1014 enddo
1015 enddo
1016 enddo
1017 enddo
1018 do l=1,psat
1019 !wig, 22-Nov-2006: added vertical bounds to prevent out-of-bounds at top
1020 if(l.eq.1)ccn1(i,kts:kte,j)=ccn(:,l)
1021 if(l.eq.2)ccn2(i,kts:kte,j)=ccn(:,l)
1022 if(l.eq.3)ccn3(i,kts:kte,j)=ccn(:,l)
1023 if(l.eq.4)ccn4(i,kts:kte,j)=ccn(:,l)
1024 if(l.eq.5)ccn5(i,kts:kte,j)=ccn(:,l)
1025 if(l.eq.6)ccn6(i,kts:kte,j)=ccn(:,l)
1026 end do
1027
1028 100 continue ! end of main loop over i
1029 120 continue ! end of main loop over j
1030
1031
1032 return
1033 end subroutine mixactivate
1034
1035
1036 !----------------------------------------------------------------------
1037 !----------------------------------------------------------------------
1038 subroutine explmix( q, src, ekkp, ekkm, overlapp, overlapm, &
1039 qold, surfrate, kts, kte, dt, is_unact, &
1040 qactold )
1041
1042 ! explicit integration of droplet/aerosol mixing
1043 ! with source due to activation/nucleation
1044
1045
1046 implicit none
1047 integer, intent(in) :: kts,kte ! number of levels
1048 real, intent(inout) :: q(kts:kte) ! mixing ratio to be updated
1049 real, intent(in) :: qold(kts:kte) ! mixing ratio from previous time step
1050 real, intent(in) :: src(kts:kte) ! source due to activation/nucleation (/s)
1051 real, intent(in) :: ekkp(kts:kte) ! zn*zs*density*diffusivity (kg/m3 m2/s) at interface
1052 ! below layer k (k,k+1 interface)
1053 real, intent(in) :: ekkm(kts:kte) ! zn*zs*density*diffusivity (kg/m3 m2/s) at interface
1054 ! above layer k (k,k+1 interface)
1055 real, intent(in) :: overlapp(kts:kte) ! cloud overlap below
1056 real, intent(in) :: overlapm(kts:kte) ! cloud overlap above
1057 real, intent(in) :: surfrate ! surface exchange rate (/s)
1058 real, intent(in) :: dt ! time step (s)
1059 logical, intent(in) :: is_unact ! true if this is an unactivated species
1060 real, intent(in),optional :: qactold(kts:kte)
1061 ! mixing ratio of ACTIVATED species from previous step
1062 ! *** this should only be present
1063 ! if the current species is unactivated number/sfc/mass
1064
1065 integer k,kp1,km1
1066
1067 if ( is_unact ) then
1068 ! the qactold*(1-overlap) terms are resuspension of activated material
1069 do k=kts,kte
1070 kp1=min(k+1,kte)
1071 km1=max(k-1,kts)
1072 q(k) = qold(k) + dt*( - src(k) + ekkp(k)*(qold(kp1) - qold(k) + &
1073 qactold(kp1)*(1.0-overlapp(k))) &
1074 + ekkm(k)*(qold(km1) - qold(k) + &
1075 qactold(km1)*(1.0-overlapm(k))) )
1076 ! if(q(k)<-1.e-30)then ! force to non-negative
1077 ! print *,'q=',q(k),' in explmix'
1078 q(k)=max(q(k),0.)
1079 ! endif
1080 end do
1081 else
1082 do k=kts,kte
1083 kp1=min(k+1,kte)
1084 km1=max(k-1,kts)
1085 q(k) = qold(k) + dt*(src(k) + ekkp(k)*(overlapp(k)*qold(kp1)-qold(k)) + &
1086 ekkm(k)*(overlapm(k)*qold(km1)-qold(k)) )
1087 ! if(q(k)<-1.e-30)then ! force to non-negative
1088 ! print *,'q=',q(k),' in explmix'
1089 q(k)=max(q(k),0.)
1090 ! endif
1091 end do
1092 end if
1093 ! diffusion loss at base of lowest layer
1094 q(kts)=q(kts)-surfrate*qold(kts)*dt
1095
1096 ! if(q(kts)<-1.e-30)then ! force to non-negative
1097 ! print *,'q=',q(kts),' in explmix'
1098 q(kts)=max(q(kts),0.)
1099 ! endif
1100
1101 return
1102 end subroutine explmix
1103
1104 !----------------------------------------------------------------------
1105 !----------------------------------------------------------------------
1106 ! 06-nov-2005 rce - grid_id & ktau added to arg list
1107 subroutine activate(wbar, sigw, wdiab, wminf, wmaxf, tair, rhoair, &
1108 msectional, maxd_atype, ntype_aer, maxd_asize, nsize_aer, &
1109 na, volc, dlo_sect,dhi_sect,sigman, hygro, &
1110 fn, fs, fm, fluxn, fluxs, fluxm, &
1111 grid_id, ktau, ii, jj, kk )
1112
1113 ! calculates number, surface, and mass fraction of aerosols activated as CCN
1114 ! calculates flux of cloud droplets, surface area, and aerosol mass into cloud
1115 ! assumes an internal mixture within each of up to pmaxd_atype X pmaxd_asize
1116 ! multiple aerosol modes.
1117 ! A sectional treatment within each type is assumed if ntype_aer >7.
1118 ! A gaussiam spectrum of updrafts can be treated.
1119
1120 ! mks units
1121
1122 ! Abdul-Razzak and Ghan, A parameterization of aerosol activation.
1123 ! 2. Multiple aerosol types. J. Geophys. Res., 105, 6837-6844.
1124
1125 USE module_model_constants, only: g,rhowater, xlv, cp, rvovrd, r_d, r_v, &
1126 mwdry,svp1,svp2,svp3,ep_2
1127
1128 implicit none
1129
1130
1131 ! input
1132
1133 integer,intent(in) :: maxd_atype ! dimension of types
1134 integer,intent(in) :: maxd_asize ! dimension of sizes
1135 integer,intent(in) :: ntype_aer ! number of types
1136 integer,intent(in) :: nsize_aer(maxd_atype) ! number of sizes for type
1137 integer,intent(in) :: msectional ! 1 for sectional, 0 for modal
1138 integer,intent(in) :: grid_id ! WRF grid%id
1139 integer,intent(in) :: ktau ! WRF time step count
1140 integer,intent(in) :: ii, jj, kk ! i,j,k of current grid cell
1141 real,intent(in) :: wbar ! grid cell mean vertical velocity (m/s)
1142 real,intent(in) :: sigw ! subgrid standard deviation of vertical vel (m/s)
1143 real,intent(in) :: wdiab ! diabatic vertical velocity (0 if adiabatic)
1144 real,intent(in) :: wminf ! minimum updraft velocity for integration (m/s)
1145 real,intent(in) :: wmaxf ! maximum updraft velocity for integration (m/s)
1146 real,intent(in) :: tair ! air temperature (K)
1147 real,intent(in) :: rhoair ! air density (kg/m3)
1148 real,intent(in) :: na(maxd_asize,maxd_atype) ! aerosol number concentration (/m3)
1149 real,intent(in) :: sigman(maxd_asize,maxd_atype) ! geometric standard deviation of aerosol size distribution
1150 real,intent(in) :: hygro(maxd_asize,maxd_atype) ! bulk hygroscopicity of aerosol mode
1151 real,intent(in) :: volc(maxd_asize,maxd_atype) ! total aerosol volume concentration (m3/m3)
1152 real,intent(in) :: dlo_sect( maxd_asize, maxd_atype ), & ! minimum size of section (cm)
1153 dhi_sect( maxd_asize, maxd_atype ) ! maximum size of section (cm)
1154
1155 ! output
1156
1157 real,intent(inout) :: fn(maxd_asize,maxd_atype) ! number fraction of aerosols activated
1158 real,intent(inout) :: fs(maxd_asize,maxd_atype) ! surface fraction of aerosols activated
1159 real,intent(inout) :: fm(maxd_asize,maxd_atype) ! mass fraction of aerosols activated
1160 real,intent(inout) :: fluxn(maxd_asize,maxd_atype) ! flux of activated aerosol number fraction into cloud (m/s)
1161 real,intent(inout) :: fluxs(maxd_asize,maxd_atype) ! flux of activated aerosol surface fraction (m/s)
1162 real,intent(inout) :: fluxm(maxd_asize,maxd_atype) ! flux of activated aerosol mass fraction into cloud (m/s)
1163
1164 ! local
1165
1166 !!$ external erf,erfc
1167 !!$ real erf,erfc
1168 ! external qsat_water
1169 integer, parameter:: nx=200
1170 integer iquasisect_option, isectional
1171 real integ,integf
1172 real surften ! surface tension of water w/respect to air (N/m)
1173 data surften/0.076/
1174 save surften
1175 real p0 ! reference pressure (Pa)
1176 real t0 ! reference temperature (K)
1177 data p0/1013.25e2/,t0/273.15/
1178 save p0,t0
1179 real ylo(maxd_asize,maxd_atype),yhi(maxd_asize,maxd_atype) ! 1-particle volume at section interfaces
1180 real ymean(maxd_asize,maxd_atype) ! 1-particle volume at r=rmean
1181 real ycut, lnycut, betayy, betayy2, gammayy, phiyy
1182 real surfc(maxd_asize,maxd_atype) ! surface concentration (m2/m3)
1183 real sign(maxd_asize,maxd_atype) ! geometric standard deviation of size distribution
1184 real alnsign(maxd_asize,maxd_atype) ! natl log of geometric standard dev of aerosol
1185 real am(maxd_asize,maxd_atype) ! number mode radius of dry aerosol (m)
1186 real lnhygro(maxd_asize,maxd_atype) ! ln(b)
1187 real pres ! pressure (Pa)
1188 real path ! mean free path (m)
1189 real diff ! diffusivity (m2/s)
1190 real conduct ! thermal conductivity (Joule/m/sec/deg)
1191 real diff0,conduct0
1192 real es ! saturation vapor pressure
1193 real qs ! water vapor saturation mixing ratio
1194 real dqsdt ! change in qs with temperature
1195 real dqsdp ! change in qs with pressure
1196 real gg ! thermodynamic function (m2/s)
1197 real sqrtg ! sqrt(gg)
1198 real sm(maxd_asize,maxd_atype) ! critical supersaturation for number mode radius
1199 real lnsm(maxd_asize,maxd_atype) ! ln( sm )
1200 real zeta, eta(maxd_asize,maxd_atype)
1201 real lnsmax ! ln(smax)
1202 real alpha
1203 real gamma
1204 real beta
1205 real gaus
1206 logical top ! true if cloud top, false if cloud base or new cloud
1207 data top/.false./
1208 save top
1209 real asub(maxd_asize,maxd_atype),bsub(maxd_asize,maxd_atype) ! coefficients of submode size distribution N=a+bx
1210 real totn(maxd_atype) ! total aerosol number concentration
1211 real aten ! surface tension parameter
1212 real gmrad(maxd_atype) ! geometric mean radius
1213 real gmradsq(maxd_atype) ! geometric mean of radius squared
1214 real gmlnsig(maxd_atype) ! geometric standard deviation
1215 real gmsm(maxd_atype) ! critical supersaturation at radius gmrad
1216 real sumflxn(maxd_asize,maxd_atype)
1217 real sumflxs(maxd_asize,maxd_atype)
1218 real sumflxm(maxd_asize,maxd_atype)
1219 real sumfn(maxd_asize,maxd_atype)
1220 real sumfs(maxd_asize,maxd_atype)
1221 real sumfm(maxd_asize,maxd_atype)
1222 real sumns(maxd_atype)
1223 real fnold(maxd_asize,maxd_atype) ! number fraction activated
1224 real fsold(maxd_asize,maxd_atype) ! surface fraction activated
1225 real fmold(maxd_asize,maxd_atype) ! mass fraction activated
1226 real wold,gold
1227 real alogten,alog2,alog3,alogaten
1228 real alogam
1229 real rlo(maxd_asize,maxd_atype), rhi(maxd_asize,maxd_atype)
1230 real rmean(maxd_asize,maxd_atype)
1231 ! mean radius (m) for the section (not used with modal)
1232 ! calculated from current volume & number
1233 real ccc
1234 real dumaa,dumbb
1235 real wmin,wmax,w,dw,dwmax,dwmin,wnuc,dwnew,wb
1236 real dfmin,dfmax,fnew,fold,fnmin,fnbar,fsbar,fmbar
1237 real alw,sqrtalw
1238 real smax
1239 real x,arg
1240 real xmincoeff,xcut
1241 real z,z1,z2,wf1,wf2,zf1,zf2,gf1,gf2,gf
1242 real etafactor1,etafactor2(maxd_asize,maxd_atype),etafactor2max
1243 integer m,n,nw,nwmax
1244
1245 ! numerical integration parameters
1246 real eps,fmax,sds
1247 data eps/0.3/,fmax/0.99/,sds/3./
1248
1249 ! mathematical constants
1250 real third, twothird, sixth, zero, one, two, three
1251 ! 04-nov-2005 rce - make this more precise
1252 ! data third/0.333333/, twothird/0.66666667/, sixth/0.166666667/,zero/0./,one/1./,two/2./,three/3./
1253 ! data third/0.33333333333/, twothird/0.66666666667/, sixth/0.16666666667/
1254 ! data zero/0./,one/1./,two/2./,three/3./
1255 ! save third, sixth,twothird,zero,one,two,three
1256
1257 real sq2, sqpi, pi
1258 ! 04-nov-2005 rce - make this more precise
1259 ! data sq2/1.4142136/, sqpi/1.7724539/,pi/3.14159/
1260 data sq2/1.4142135624/, sqpi/1.7724538509/,pi/3.1415926536/
1261 save sq2,sqpi,pi
1262
1263 integer ndist(nx) ! accumulates frequency distribution of integration bins required
1264 data ndist/nx*0/
1265 save eps,fmax,sds,ndist
1266
1267 ! for nsize_aer>7, a sectional approach is used and isectional = iquasisect_option
1268 ! activation fractions (fn,fs,fm) are computed as follows
1269 ! iquasisect_option = 1,3 - each section treated as a narrow lognormal
1270 ! iquasisect_option = 2,4 - within-section dn/dx = a + b*x, x = ln(r)
1271 ! smax is computed as follows (when explicit activation is OFF)
1272 ! iquasisect_option = 1,2 - razzak-ghan modal parameterization with
1273 ! single mode having same ntot, dgnum, sigmag as the combined sections
1274 ! iquasisect_option = 3,4 - razzak-ghan sectional parameterization
1275 ! for nsize_aer=<9, a modal approach is used and isectional = 0
1276
1277 ! rce 08-jul-2005
1278 ! if either (na(n,m) < nsmall) or (volc(n,m) < vsmall)
1279 ! then treat bin/mode (n,m) as being empty, and set its fn/fs/fm=0.0
1280 ! (for single precision, gradual underflow starts around 1.0e-38,
1281 ! and strange things can happen when in that region)
1282 real, parameter :: nsmall = 1.0e-20 ! aer number conc in #/m3
1283 real, parameter :: vsmall = 1.0e-37 ! aer volume conc in m3/m3
1284 logical bin_is_empty(maxd_asize,maxd_atype), all_bins_empty
1285 logical bin_is_narrow(maxd_asize,maxd_atype)
1286
1287 integer idiagaa, ipass_nwloop
1288 integer idiag_dndy_neg, idiag_fnsm_prob
1289
1290 !.......................................................................
1291 !
1292 ! start calc. of modal or sectional activation properties (start of section 1)
1293 !
1294 !.......................................................................
1295 idiag_dndy_neg = 1 ! set this to 0 to turn off
1296 ! warnings about dn/dy < 0
1297 idiag_fnsm_prob = 1 ! set this to 0 to turn off
1298 ! warnings about fn/fs/fm misbehavior
1299
1300 iquasisect_option = 2
1301 if(msectional.gt.0)then
1302 isectional = iquasisect_option
1303 else
1304 isectional = 0
1305 endif
1306
1307 do n=1,ntype_aer
1308 ! print *,'ntype_aer,n,nsize_aer(n)=',ntype_aer,n,nsize_aer(n)
1309
1310 if(ntype_aer.eq.1.and.nsize_aer(n).eq.1.and.na(1,1).lt.1.e-20)then
1311 fn(1,1)=0.
1312 fs(1,1)=0.
1313 fm(1,1)=0.
1314 fluxn(1,1)=0.
1315 fluxs(1,1)=0.
1316 fluxm(1,1)=0.
1317 return
1318 endif
1319 enddo
1320
1321 zero = 0.0
1322 one = 1.0
1323 two = 2.0
1324 three = 3.0
1325 third = 1.0/3.0
1326 twothird = 2.0/6.0
1327 sixth = 1.0/6.0
1328
1329 pres=r_d*rhoair*tair
1330 diff0=0.211e-4*(p0/pres)*(tair/t0)**1.94
1331 conduct0=(5.69+0.017*(tair-t0))*4.186e2*1.e-5 ! convert to J/m/s/deg
1332 es=1000.*svp1*exp( svp2*(tair-t0)/(tair-svp3) )
1333 qs=ep_2*es/(pres-es)
1334 dqsdt=xlv/(r_v*tair*tair)*qs
1335 alpha=g*(xlv/(cp*r_v*tair*tair)-1./(r_d*tair))
1336 gamma=(1+xlv/cp*dqsdt)/(rhoair*qs)
1337 gg=1./(rhowater/(diff0*rhoair*qs)+xlv*rhowater/(conduct0*tair)*(xlv/(r_v*tair)-1.))
1338 sqrtg=sqrt(gg)
1339 beta=4.*pi*rhowater*gg*gamma
1340 aten=2.*surften/(r_v*tair*rhowater)
1341 alogaten=log(aten)
1342 alog2=log(two)
1343 alog3=log(three)
1344 ccc=4.*pi*third
1345 etafactor2max=1.e10/(alpha*wmaxf)**1.5 ! this should make eta big if na is very small.
1346
1347 all_bins_empty = .true.
1348 do n=1,ntype_aer
1349 totn(n)=0.
1350 gmrad(n)=0.
1351 gmradsq(n)=0.
1352 sumns(n)=0.
1353 do m=1,nsize_aer(n)
1354 alnsign(m,n)=log(sigman(m,n))
1355 ! internal mixture of aerosols
1356
1357 bin_is_empty(m,n) = .true.
1358 if (volc(m,n).gt.vsmall .and. na(m,n).gt.nsmall) then
1359 bin_is_empty(m,n) = .false.
1360 all_bins_empty = .false.
1361 lnhygro(m,n)=log(hygro(m,n))
1362 ! number mode radius (m,n)
1363 ! write(6,*)'alnsign,volc,na=',alnsign(m,n),volc(m,n),na(m,n)
1364 am(m,n)=exp(-1.5*alnsign(m,n)*alnsign(m,n))* &
1365 (3.*volc(m,n)/(4.*pi*na(m,n)))**third
1366
1367 if (isectional .gt. 0) then
1368 ! sectional model.
1369 ! need to use bulk properties because parameterization doesn't
1370 ! work well for narrow bins.
1371 totn(n)=totn(n)+na(m,n)
1372 alogam=log(am(m,n))
1373 gmrad(n)=gmrad(n)+na(m,n)*alogam
1374 gmradsq(n)=gmradsq(n)+na(m,n)*alogam*alogam
1375 endif
1376 etafactor2(m,n)=1./(na(m,n)*beta*sqrtg)
1377
1378 if(hygro(m,n).gt.1.e-10)then
1379 sm(m,n)=2.*aten/(3.*am(m,n))*sqrt(aten/(3.*hygro(m,n)*am(m,n)))
1380 else
1381 sm(m,n)=100.
1382 endif
1383 ! write(6,*)'sm,hygro,am=',sm(m,n),hygro(m,n),am(m,n)
1384 else
1385 sm(m,n)=1.
1386 etafactor2(m,n)=etafactor2max ! this should make eta big if na is very small.
1387
1388 endif
1389 lnsm(m,n)=log(sm(m,n))
1390 if ((isectional .eq. 3) .or. (isectional .eq. 4)) then
1391 sumns(n)=sumns(n)+na(m,n)/sm(m,n)**twothird
1392 endif
1393 ! write(6,'(a,i4,6g12.2)')'m,na,am,hygro,lnhygro,sm,lnsm=',m,na(m,n),am(m,n),hygro(m,n),lnhygro(m,n),sm(m,n),lnsm(m,n)
1394 end do ! size
1395 end do ! type
1396
1397 ! if all bins are empty, set all activation fractions to zero and exit
1398 if ( all_bins_empty ) then
1399 do n=1,ntype_aer
1400 do m=1,nsize_aer(n)
1401 fluxn(m,n)=0.
1402 fn(m,n)=0.
1403 fluxs(m,n)=0.
1404 fs(m,n)=0.
1405 fluxm(m,n)=0.
1406 fm(m,n)=0.
1407 end do
1408 end do
1409 return
1410 endif
1411
1412
1413
1414 if (isectional .le. 0) goto 30000
1415
1416 do n=1,ntype_aer
1417 !wig 19-Oct-2006: Add zero trap based May 2006 e-mail from
1418 !Ghan. Transport can clear out a cell leading to
1419 !inconsistencies with the mass.
1420 gmrad(n)=gmrad(n)/max(totn(n),1e-20)
1421 gmlnsig=gmradsq(n)/totn(n)-gmrad(n)*gmrad(n) ! [ln(sigmag)]**2
1422 gmlnsig(n)=sqrt( max( 1.e-4, gmlnsig(n) ) )
1423 gmrad(n)=exp(gmrad(n))
1424 if ((isectional .eq. 3) .or. (isectional .eq. 4)) then
1425 gmsm(n)=totn(n)/sumns(n)
1426 gmsm(n)=gmsm(n)*gmsm(n)*gmsm(n)
1427 gmsm(n)=sqrt(gmsm(n))
1428 else
1429 ! gmsm(n)=2.*aten/(3.*gmrad(n))*sqrt(aten/(3.*hygro(1,n)*gmrad(n)))
1430 gmsm(n)=2.*aten/(3.*gmrad(n))*sqrt(aten/(3.*hygro(nsize_aer(n),n)*gmrad(n)))
1431 endif
1432 enddo
1433
1434 !.......................................................................
1435 ! calculate sectional "sub-bin" size distribution
1436 !
1437 ! dn/dy = nt*( a + b*y ) for ylo < y < yhi
1438 !
1439 ! nt = na(m,n) = number mixing ratio of the bin
1440 ! y = v/vhi
1441 ! v = (4pi/3)*r**3 = particle volume
1442 ! vhi = v at r=rhi (upper bin boundary)
1443 ! ylo = y at lower bin boundary = vlo/vhi = (rlo/rhi)**3
1444 ! yhi = y at upper bin boundary = 1.0
1445 !
1446 ! dv/dy = v * dn/dy = nt*vhi*( a*y + b*y*y )
1447 !
1448 !.......................................................................
1449 ! 02-may-2006 - this dn/dy replaces the previous
1450 ! dn/dx = a + b*x where l = ln(r)
1451 ! the old dn/dx was overly complicated for cases of rmean near rlo or rhi
1452 ! the new dn/dy is consistent with that used in the movesect routine,
1453 ! which does continuous growth by condensation and aqueous chemistry
1454 !.......................................................................
1455 do 25002 n = 1,ntype_aer
1456 do 25000 m = 1,nsize_aer(n)
1457
1458 ! convert from diameter in cm to radius in m
1459 rlo(m,n) = 0.5*0.01*dlo_sect(m,n)
1460 rhi(m,n) = 0.5*0.01*dhi_sect(m,n)
1461 ylo(m,n) = (rlo(m,n)/rhi(m,n))**3
1462 yhi(m,n) = 1.0
1463
1464 ! 04-nov-2005 - extremely narrow bins will be treated using 0/1 activation
1465 ! this is to avoid potential numberical problems
1466 bin_is_narrow(m,n) = .false.
1467 if ((rhi(m,n)/rlo(m,n)) .le. 1.01) bin_is_narrow(m,n) = .true.
1468
1469 ! rmean is mass mean radius for the bin; xmean = log(rmean)
1470 ! just use section midpoint if bin is empty
1471 if ( bin_is_empty(m,n) ) then
1472 rmean(m,n) = sqrt(rlo(m,n)*rhi(m,n))
1473 ymean(m,n) = (rmean(m,n)/rhi(m,n))**3
1474 goto 25000
1475 end if
1476
1477 rmean(m,n) = (volc(m,n)/(ccc*na(m,n)))**third
1478 rmean(m,n) = max( rlo(m,n), min( rhi(m,n), rmean(m,n) ) )
1479 ymean(m,n) = (rmean(m,n)/rhi(m,n))**3
1480 if ( bin_is_narrow(m,n) ) goto 25000
1481
1482 ! if rmean is extremely close to either rlo or rhi,
1483 ! treat the bin as extremely narrow
1484 if ((rhi(m,n)/rmean(m,n)) .le. 1.01) then
1485 bin_is_narrow(m,n) = .true.
1486 rlo(m,n) = min( rmean(m,n), (rhi(m,n)/1.01) )
1487 ylo(m,n) = (rlo(m,n)/rhi(m,n))**3
1488 goto 25000
1489 else if ((rmean(m,n)/rlo(m,n)) .le. 1.01) then
1490 bin_is_narrow(m,n) = .true.
1491 rhi(m,n) = max( rmean(m,n), (rlo(m,n)*1.01) )
1492 ylo(m,n) = (rlo(m,n)/rhi(m,n))**3
1493 ymean(m,n) = (rmean(m,n)/rhi(m,n))**3
1494 goto 25000
1495 endif
1496
1497 ! if rmean is somewhat close to either rlo or rhi, then dn/dy will be
1498 ! negative near the upper or lower bin boundary
1499 ! in these cases, assume that all the particles are in a subset of the full bin,
1500 ! and adjust rlo or rhi so that rmean will be near the center of this subset
1501 ! note that the bin is made narrower LOCALLY/TEMPORARILY,
1502 ! just for the purposes of the activation calculation
1503 gammayy = (ymean(m,n)-ylo(m,n)) / (yhi(m,n)-ylo(m,n))
1504 if (gammayy .lt. 0.34) then
1505 dumaa = ylo(m,n) + (yhi(m,n)-ylo(m,n))*(gammayy/0.34)
1506 rhi(m,n) = rhi(m,n)*(dumaa**third)
1507 ylo(m,n) = (rlo(m,n)/rhi(m,n))**3
1508 ymean(m,n) = (rmean(m,n)/rhi(m,n))**3
1509 else if (gammayy .ge. 0.66) then
1510 dumaa = ylo(m,n) + (yhi(m,n)-ylo(m,n))*((gammayy-0.66)/0.34)
1511 ylo(m,n) = dumaa
1512 rlo(m,n) = rhi(m,n)*(dumaa**third)
1513 end if
1514 if ((rhi(m,n)/rlo(m,n)) .le. 1.01) then
1515 bin_is_narrow(m,n) = .true.
1516 goto 25000
1517 end if
1518
1519 betayy = ylo(m,n)/yhi(m,n)
1520 betayy2 = betayy*betayy
1521 bsub(m,n) = (12.0*ymean(m,n) - 6.0*(1.0+betayy)) / &
1522 (4.0*(1.0-betayy2*betayy) - 3.0*(1.0-betayy2)*(1.0+betayy))
1523 asub(m,n) = (1.0 - bsub(m,n)*(1.0-betayy2)*0.5) / (1.0-betayy)
1524
1525 if ( asub(m,n)+bsub(m,n)*ylo(m,n) .lt. 0. ) then
1526 if (idiag_dndy_neg .gt. 0) then
1527 print *,'dndy<0 at lower boundary'
1528 print *,'n,m=',n,m
1529 print *,'na=',na(m,n),' volc=',volc(m,n)
1530 print *,'volc/(na*pi*4/3)=', (volc(m,n)/(na(m,n)*ccc))
1531 print *,'rlo(m,n),rhi(m,n)=',rlo(m,n),rhi(m,n)
1532 print *,'dlo_sect/2,dhi_sect/2=', &
1533 (0.005*dlo_sect(m,n)),(0.005*dhi_sect(m,n))
1534 print *,'asub,bsub,ylo,yhi=',asub(m,n),bsub(m,n),ylo(m,n),yhi(m,n)
1535 print *,'asub+bsub*ylo=', &
1536 (asub(m,n)+bsub(m,n)*ylo(m,n))
1537 print *,'subr activate error 11 - i,j,k =', ii, jj, kk
1538 ! 07-nov-2005 rce - don't stop for this, it's not fatal
1539 ! stop
1540 endif
1541 endif
1542 if ( asub(m,n)+bsub(m,n)*yhi(m,n) .lt. 0. ) then
1543 if (idiag_dndy_neg .gt. 0) then
1544 print *,'dndy<0 at upper boundary'
1545 print *,'n,m=',n,m
1546 print *,'na=',na(m,n),' volc=',volc(m,n)
1547 print *,'volc/(na*pi*4/3)=', (volc(m,n)/(na(m,n)*ccc))
1548 print *,'rlo(m,n),rhi(m,n)=',rlo(m,n),rhi(m,n)
1549 print *,'dlo_sect/2,dhi_sect/2=', &
1550 (0.005*dlo_sect(m,n)),(0.005*dhi_sect(m,n))
1551 print *,'asub,bsub,ylo,yhi=',asub(m,n),bsub(m,n),ylo(m,n),yhi(m,n)
1552 print *,'asub+bsub*yhi=', &
1553 (asub(m,n)+bsub(m,n)*yhi(m,n))
1554 print *,'subr activate error 12 - i,j,k =', ii, jj, kk
1555 ! stop
1556 endif
1557 endif
1558
1559 25000 continue ! m=1,nsize_aer(n)
1560 25002 continue ! n=1,ntype_aer
1561
1562
1563 30000 continue
1564 !.......................................................................
1565 !
1566 ! end calc. of modal or sectional activation properties (end of section 1)
1567 !
1568 !.......................................................................
1569
1570
1571
1572 ! sjg 7-16-98 upward
1573 ! print *,'wbar,sigw=',wbar,sigw
1574
1575 if(sigw.le.1.e-5) goto 50000
1576
1577 !.......................................................................
1578 !
1579 ! start calc. of activation fractions/fluxes
1580 ! for spectrum of updrafts (start of section 2)
1581 !
1582 !.......................................................................
1583 ipass_nwloop = 1
1584 idiagaa = 0
1585 ! 06-nov-2005 rce - set idiagaa=1 for testing/debugging
1586 ! if ((grid_id.eq.1) .and. (ktau.eq.167) .and. &
1587 ! (ii.eq.24) .and. (jj.eq. 1) .and. (kk.eq.14)) idiagaa = 1
1588
1589 40000 continue
1590 if(top)then
1591 wmax=0.
1592 wmin=min(zero,-wdiab)
1593 else
1594 wmax=min(wmaxf,wbar+sds*sigw)
1595 wmin=max(wminf,-wdiab)
1596 endif
1597 wmin=max(wmin,wbar-sds*sigw)
1598 w=wmin
1599 dwmax=eps*sigw
1600 dw=dwmax
1601 dfmax=0.2
1602 dfmin=0.1
1603 if(wmax.le.w)then
1604 do n=1,ntype_aer
1605 do m=1,nsize_aer(n)
1606 fluxn(m,n)=0.
1607 fn(m,n)=0.
1608 fluxs(m,n)=0.
1609 fs(m,n)=0.
1610 fluxm(m,n)=0.
1611 fm(m,n)=0.
1612 end do
1613 end do
1614 return
1615 endif
1616 do n=1,ntype_aer
1617 do m=1,nsize_aer(n)
1618 sumflxn(m,n)=0.
1619 sumfn(m,n)=0.
1620 fnold(m,n)=0.
1621 sumflxs(m,n)=0.
1622 sumfs(m,n)=0.
1623 fsold(m,n)=0.
1624 sumflxm(m,n)=0.
1625 sumfm(m,n)=0.
1626 fmold(m,n)=0.
1627 enddo
1628 enddo
1629
1630 fold=0
1631 gold=0
1632 ! 06-nov-2005 rce - set wold=w here
1633 ! wold=0
1634 wold=w
1635
1636
1637 ! 06-nov-2005 rce - define nwmax; calc dwmin from nwmax
1638 nwmax = 200
1639 ! dwmin = min( dwmax, 0.01 )
1640 dwmin = (wmax - wmin)/(nwmax-1)
1641 dwmin = min( dwmax, dwmin )
1642 dwmin = max( 0.01, dwmin )
1643
1644 !
1645 ! loop over updrafts, incrementing sums as you go
1646 ! the "200" is (arbitrary) upper limit for number of updrafts
1647 ! if integration finishes before this, OK; otherwise, ERROR
1648 !
1649 if (idiagaa.gt.0) then
1650 write(*,94700) ktau, grid_id, ii, jj, kk, nwmax
1651 write(*,94710) 'wbar,sigw,wdiab=', wbar, sigw, wdiab
1652 write(*,94710) 'wmin,wmax,dwmin,dwmax=', wmin, wmax, dwmin, dwmax
1653 write(*,94720) -1, w, wold, dw
1654 end if
1655 94700 format( / 'activate 47000 - ktau,id,ii,jj,kk,nwmax=', 6i5 )
1656 94710 format( 'activate 47000 - ', a, 6(1x,f11.5) )
1657 94720 format( 'activate 47000 - nw,w,wold,dw=', i5, 3(1x,f11.5) )
1658
1659 do 47000 nw = 1, nwmax
1660 41000 wnuc=w+wdiab
1661
1662 if (idiagaa.gt.0) write(*,94720) nw, w, wold, dw
1663
1664 ! write(6,*)'wnuc=',wnuc
1665 alw=alpha*wnuc
1666 sqrtalw=sqrt(alw)
1667 zeta=2.*sqrtalw*aten/(3.*sqrtg)
1668 etafactor1=2.*alw*sqrtalw
1669 if (isectional .gt. 0) then
1670 ! sectional model.
1671 ! use bulk properties
1672
1673 do n=1,ntype_aer
1674 if(totn(n).gt.1.e-10)then
1675 eta(1,n)=etafactor1/(totn(n)*beta*sqrtg)
1676 else
1677 eta(1,n)=1.e10
1678 endif
1679 enddo
1680 call maxsat(zeta,eta,maxd_atype,ntype_aer,maxd_asize,(/1/),gmsm,gmlnsig,smax)
1681 ! call maxsat(zeta,eta,maxd_atype,ntype_aer,maxd_asize,nsize_aer,gmsm,gmlnsig,smax)
1682 lnsmax=log(smax)
1683 x=2*(log(gmsm(1))-lnsmax)/(3*sq2*gmlnsig(1))
1684 fnew=0.5*(1.-ERF_ALT(x))
1685
1686 else
1687
1688 do n=1,ntype_aer
1689 do m=1,nsize_aer(n)
1690 eta(m,n)=etafactor1*etafactor2(m,n)
1691 enddo
1692 enddo
1693
1694 call maxsat(zeta,eta,maxd_atype,ntype_aer,maxd_asize,nsize_aer,sm,alnsign,smax)
1695 ! write(6,*)'w,smax=',w,smax
1696
1697 lnsmax=log(smax)
1698
1699 x=2*(lnsm(nsize_aer(1),1)-lnsmax)/(3*sq2*alnsign(nsize_aer(1),1))
1700 fnew=0.5*(1.-ERF_ALT(x))
1701
1702 endif
1703
1704 dwnew = dw
1705 ! 06-nov-2005 rce - "n" here should be "nw" (?)
1706 ! if(fnew-fold.gt.dfmax.and.n.gt.1)then
1707 if(fnew-fold.gt.dfmax.and.nw.gt.1)then
1708 ! reduce updraft increment for greater accuracy in integration
1709 if (dw .gt. 1.01*dwmin) then
1710 dw=0.7*dw
1711 dw=max(dw,dwmin)
1712 w=wold+dw
1713 go to 41000
1714 else
1715 dwnew = dwmin
1716 endif
1717 endif
1718
1719 if(fnew-fold.lt.dfmin)then
1720 ! increase updraft increment to accelerate integration
1721 dwnew=min(1.5*dw,dwmax)
1722 endif
1723 fold=fnew
1724
1725 z=(w-wbar)/(sigw*sq2)
1726 gaus=exp(-z*z)
1727 fnmin=1.
1728 xmincoeff=alogaten-2.*third*(lnsmax-alog2)-alog3
1729 ! write(6,*)'xmincoeff=',xmincoeff
1730
1731
1732 do 44002 n=1,ntype_aer
1733 do 44000 m=1,nsize_aer(n)
1734 if ( bin_is_empty(m,n) ) then
1735 fn(m,n)=0.
1736 fs(m,n)=0.
1737 fm(m,n)=0.
1738 else if ((isectional .eq. 2) .or. (isectional .eq. 4)) then
1739 ! sectional
1740 ! within-section dn/dx = a + b*x
1741 xcut=xmincoeff-third*lnhygro(m,n)
1742 ! ycut=(exp(xcut)/rhi(m,n))**3
1743 ! 07-jul-2006 rce - the above line gave a (rare) overflow when smax=1.0e-20
1744 ! if (ycut > yhi), then actual value of ycut is unimportant,
1745 ! so do the following to avoid overflow
1746 lnycut = 3.0 * ( xcut - log(rhi(m,n)) )
1747 lnycut = min( lnycut, log(yhi(m,n)*1.0e5) )
1748 ycut=exp(lnycut)
1749 ! write(6,*)'m,n,rcut,rlo,rhi=',m,n,exp(xcut),rlo(m,n),rhi(m,n)
1750 ! if(lnsmax.lt.lnsmn(m,n))then
1751 if(ycut.gt.yhi(m,n))then
1752 fn(m,n)=0.
1753 fs(m,n)=0.
1754 fm(m,n)=0.
1755 elseif(ycut.lt.ylo(m,n))then
1756 fn(m,n)=1.
1757 fs(m,n)=1.
1758 fm(m,n)=1.
1759 elseif ( bin_is_narrow(m,n) ) then
1760 ! 04-nov-2005 rce - for extremely narrow bins,
1761 ! do zero activation if xcut>xmean, 100% activation otherwise
1762 if (ycut.gt.ymean(m,n)) then
1763 fn(m,n)=0.
1764 fs(m,n)=0.
1765 fm(m,n)=0.
1766 else
1767 fn(m,n)=1.
1768 fs(m,n)=1.
1769 fm(m,n)=1.
1770 endif
1771 else
1772 phiyy=ycut/yhi(m,n)
1773 fn(m,n) = asub(m,n)*(1.0-phiyy) + 0.5*bsub(m,n)*(1.0-phiyy*phiyy)
1774 if (fn(m,n).lt.zero .or. fn(m,n).gt.one) then
1775 if (idiag_fnsm_prob .gt. 0) then
1776 print *,'fn(',m,n,')=',fn(m,n),' outside 0,1 - activate err21'
1777 print *,'na,volc =', na(m,n), volc(m,n)
1778 print *,'asub,bsub =', asub(m,n), bsub(m,n)
1779 print *,'yhi,ycut =', yhi(m,n), ycut
1780 endif
1781 endif
1782
1783 if (fn(m,n) .le. zero) then
1784 ! 10-nov-2005 rce - if fn=0, then fs & fm must be 0
1785 fn(m,n)=zero
1786 fs(m,n)=zero
1787 fm(m,n)=zero
1788 else if (fn(m,n) .ge. one) then
1789 ! 10-nov-2005 rce - if fn=1, then fs & fm must be 1
1790 fn(m,n)=one
1791 fs(m,n)=one
1792 fm(m,n)=one
1793 else
1794 ! 10-nov-2005 rce - otherwise, calc fm and check it
1795 fm(m,n) = (yhi(m,n)/ymean(m,n)) * (0.5*asub(m,n)*(1.0-phiyy*phiyy) + &
1796 third*bsub(m,n)*(1.0-phiyy*phiyy*phiyy))
1797 if (fm(m,n).lt.fn(m,n) .or. fm(m,n).gt.one) then
1798 if (idiag_fnsm_prob .gt. 0) then
1799 print *,'fm(',m,n,')=',fm(m,n),' outside fn,1 - activate err22'
1800 print *,'na,volc,fn =', na(m,n), volc(m,n), fn(m,n)
1801 print *,'asub,bsub =', asub(m,n), bsub(m,n)
1802 print *,'yhi,ycut =', yhi(m,n), ycut
1803 endif
1804 endif
1805 if (fm(m,n) .le. fn(m,n)) then
1806 ! 10-nov-2005 rce - if fm=fn, then fs must =fn
1807 fm(m,n)=fn(m,n)
1808 fs(m,n)=fn(m,n)
1809 else if (fm(m,n) .ge. one) then
1810 ! 10-nov-2005 rce - if fm=1, then fs & fn must be 1
1811 fm(m,n)=one
1812 fs(m,n)=one
1813 fn(m,n)=one
1814 else
1815 ! 10-nov-2005 rce - these two checks assure that the mean size
1816 ! of the activated & interstitial particles will be between rlo & rhi
1817 dumaa = fn(m,n)*(yhi(m,n)/ymean(m,n))
1818 fm(m,n) = min( fm(m,n), dumaa )
1819 dumaa = 1.0 + (fn(m,n)-1.0)*(ylo(m,n)/ymean(m,n))
1820 fm(m,n) = min( fm(m,n), dumaa )
1821 ! 10-nov-2005 rce - now calculate fs and bound it by fn, fm
1822 betayy = ylo(m,n)/yhi(m,n)
1823 dumaa = phiyy**twothird
1824 dumbb = betayy**twothird
1825 fs(m,n) = &
1826 (asub(m,n)*(1.0-phiyy*dumaa) + &
1827 0.625*bsub(m,n)*(1.0-phiyy*phiyy*dumaa)) / &
1828 (asub(m,n)*(1.0-betayy*dumbb) + &
1829 0.625*bsub(m,n)*(1.0-betayy*betayy*dumbb))
1830 fs(m,n)=max(fs(m,n),fn(m,n))
1831 fs(m,n)=min(fs(m,n),fm(m,n))
1832 endif
1833 endif
1834 endif
1835
1836 else
1837 ! modal
1838 x=2*(lnsm(m,n)-lnsmax)/(3*sq2*alnsign(m,n))
1839 fn(m,n)=0.5*(1.-ERF_ALT(x))
1840 arg=x-sq2*alnsign(m,n)
1841 fs(m,n)=0.5*(1.-ERF_ALT(arg))
1842 arg=x-1.5*sq2*alnsign(m,n)
1843 fm(m,n)=0.5*(1.-ERF_ALT(arg))
1844 ! print *,'w,x,fn,fs,fm=',w,x,fn(m,n),fs(m,n),fm(m,n)
1845 endif
1846
1847 ! fn(m,n)=1. !test
1848 ! fs(m,n)=1.
1849 ! fm(m,n)=1.
1850 fnmin=min(fn(m,n),fnmin)
1851 ! integration is second order accurate
1852 ! assumes linear variation of f*gaus with w
1853 wb=(w+wold)
1854 fnbar=(fn(m,n)*gaus+fnold(m,n)*gold)
1855 fsbar=(fs(m,n)*gaus+fsold(m,n)*gold)
1856 fmbar=(fm(m,n)*gaus+fmold(m,n)*gold)
1857 if((top.and.w.lt.0.).or.(.not.top.and.w.gt.0.))then
1858 sumflxn(m,n)=sumflxn(m,n)+sixth*(wb*fnbar &
1859 +(fn(m,n)*gaus*w+fnold(m,n)*gold*wold))*dw
1860 sumflxs(m,n)=sumflxs(m,n)+sixth*(wb*fsbar &
1861 +(fs(m,n)*gaus*w+fsold(m,n)*gold*wold))*dw
1862 sumflxm(m,n)=sumflxm(m,n)+sixth*(wb*fmbar &
1863 +(fm(m,n)*gaus*w+fmold(m,n)*gold*wold))*dw
1864 endif
1865 sumfn(m,n)=sumfn(m,n)+0.5*fnbar*dw
1866 ! write(6,'(a,9g10.2)')'lnsmax,lnsm(m,n),x,fn(m,n),fnold(m,n),g,gold,fnbar,dw=', &
1867 ! lnsmax,lnsm(m,n),x,fn(m,n),fnold(m,n),g,gold,fnbar,dw
1868 fnold(m,n)=fn(m,n)
1869 sumfs(m,n)=sumfs(m,n)+0.5*fsbar*dw
1870 fsold(m,n)=fs(m,n)
1871 sumfm(m,n)=sumfm(m,n)+0.5*fmbar*dw
1872 fmold(m,n)=fm(m,n)
1873
1874 44000 continue ! m=1,nsize_aer(n)
1875 44002 continue ! n=1,ntype_aer
1876
1877 ! sumg=sumg+0.5*(gaus+gold)*dw
1878 gold=gaus
1879 wold=w
1880 dw=dwnew
1881
1882 if(nw.gt.1.and.(w.gt.wmax.or.fnmin.gt.fmax))go to 48000
1883 w=w+dw
1884
1885 47000 continue ! nw = 1, nwmax
1886
1887
1888 print *,'do loop is too short in activate'
1889 print *,'wmin=',wmin,' w=',w,' wmax=',wmax,' dw=',dw
1890 print *,'wbar=',wbar,' sigw=',sigw,' wdiab=',wdiab
1891 print *,'wnuc=',wnuc
1892 do n=1,ntype_aer
1893 print *,'ntype=',n
1894 print *,'na=',(na(m,n),m=1,nsize_aer(n))
1895 print *,'fn=',(fn(m,n),m=1,nsize_aer(n))
1896 end do
1897 ! dump all subr parameters to allow testing with standalone code
1898 ! (build a driver that will read input and call activate)
1899 print *,'top,wbar,sigw,wdiab,tair,rhoair,ntype_aer='
1900 print *, top,wbar,sigw,wdiab,tair,rhoair,ntype_aer
1901 print *,'na='
1902 print *, na
1903 print *,'volc='
1904 print *, volc
1905 print *,'sigman='
1906 print *, sigman
1907 print *,'hygro='
1908 print *, hygro
1909
1910 print *,'subr activate error 31 - i,j,k =', ii, jj, kk
1911 ! 06-nov-2005 rce - if integration fails, repeat it once with additional diagnostics
1912 if (ipass_nwloop .eq. 1) then
1913 ipass_nwloop = 2
1914 idiagaa = 2
1915 goto 40000
1916 end if
1917 stop
1918
1919 48000 continue
1920
1921
1922 ndist(n)=ndist(n)+1
1923 if(.not.top.and.w.lt.wmaxf)then
1924
1925 ! contribution from all updrafts stronger than wmax
1926 ! assuming constant f (close to fmax)
1927 wnuc=w+wdiab
1928
1929 z1=(w-wbar)/(sigw*sq2)
1930 z2=(wmaxf-wbar)/(sigw*sq2)
1931 integ=sigw*0.5*sq2*sqpi*(ERFC_NUM_RECIPES(z1)-ERFC_NUM_RECIPES(z2))
1932 ! consider only upward flow into cloud base when estimating flux
1933 wf1=max(w,zero)
1934 zf1=(wf1-wbar)/(sigw*sq2)
1935 gf1=exp(-zf1*zf1)
1936 wf2=max(wmaxf,zero)
1937 zf2=(wf2-wbar)/(sigw*sq2)
1938 gf2=exp(-zf2*zf2)
1939 gf=(gf1-gf2)
1940 integf=wbar*sigw*0.5*sq2*sqpi*(ERFC_NUM_RECIPES(zf1)-ERFC_NUM_RECIPES(zf2))+sigw*sigw*gf
1941
1942 do n=1,ntype_aer
1943 do m=1,nsize_aer(n)
1944 sumflxn(m,n)=sumflxn(m,n)+integf*fn(m,n)
1945 sumfn(m,n)=sumfn(m,n)+fn(m,n)*integ
1946 sumflxs(m,n)=sumflxs(m,n)+integf*fs(m,n)
1947 sumfs(m,n)=sumfs(m,n)+fs(m,n)*integ
1948 sumflxm(m,n)=sumflxm(m,n)+integf*fm(m,n)
1949 sumfm(m,n)=sumfm(m,n)+fm(m,n)*integ
1950 end do
1951 end do
1952 ! sumg=sumg+integ
1953 endif
1954
1955
1956 do n=1,ntype_aer
1957 do m=1,nsize_aer(n)
1958
1959 ! fn(m,n)=sumfn(m,n)/(sumg)
1960 fn(m,n)=sumfn(m,n)/(sq2*sqpi*sigw)
1961 fluxn(m,n)=sumflxn(m,n)/(sq2*sqpi*sigw)
1962 if(fn(m,n).gt.1.01)then
1963 if (idiag_fnsm_prob .gt. 0) then
1964 print *,'fn=',fn(m,n),' > 1 - activate err41'
1965 print *,'w,m,n,na,am=',w,m,n,na(m,n),am(m,n)
1966 print *,'integ,sumfn,sigw=',integ,sumfn(m,n),sigw
1967 print *,'subr activate error - i,j,k =', ii, jj, kk
1968 ! call exit
1969 endif
1970 fluxn(m,n) = fluxn(m,n)/fn(m,n)
1971 endif
1972
1973 fs(m,n)=sumfs(m,n)/(sq2*sqpi*sigw)
1974 fluxs(m,n)=sumflxs(m,n)/(sq2*sqpi*sigw)
1975 if(fs(m,n).gt.1.01)then
1976 if (idiag_fnsm_prob .gt. 0) then
1977 print *,'fs=',fs(m,n),' > 1 - activate err42'
1978 print *,'m,n,isectional=',m,n,isectional
1979 print *,'alnsign(m,n)=',alnsign(m,n)
1980 print *,'rcut,rlo(m,n),rhi(m,n)',exp(xcut),rlo(m,n),rhi(m,n)
1981 print *,'w,m,na,am=',w,m,na(m,n),am(m,n)
1982 print *,'integ,sumfs,sigw=',integ,sumfs(m,n),sigw
1983 endif
1984 fluxs(m,n) = fluxs(m,n)/fs(m,n)
1985 endif
1986
1987 ! fm(m,n)=sumfm(m,n)/(sumg)
1988 fm(m,n)=sumfm(m,n)/(sq2*sqpi*sigw)
1989 fluxm(m,n)=sumflxm(m,n)/(sq2*sqpi*sigw)
1990 if(fm(m,n).gt.1.01)then
1991 if (idiag_fnsm_prob .gt. 0) then
1992 print *,'fm(',m,n,')=',fm(m,n),' > 1 - activate err43'
1993 endif
1994 fluxm(m,n) = fluxm(m,n)/fm(m,n)
1995 endif
1996
1997 end do
1998 end do
1999
2000 goto 60000
2001 !.......................................................................
2002 !
2003 ! end calc. of activation fractions/fluxes
2004 ! for spectrum of updrafts (end of section 2)
2005 !
2006 !.......................................................................
2007
2008
2009 !.......................................................................
2010 !
2011 ! start calc. of activation fractions/fluxes
2012 ! for (single) uniform updraft (start of section 3)
2013 !
2014 !.......................................................................
2015 50000 continue
2016 wnuc=wbar+wdiab
2017 ! write(6,*)'uniform updraft =',wnuc
2018
2019 ! 04-nov-2005 rce - moved the code for "wnuc.le.0" code to here
2020 if(wnuc.le.0.)then
2021 do n=1,ntype_aer
2022 do m=1,nsize_aer(n)
2023 fn(m,n)=0
2024 fluxn(m,n)=0
2025 fs(m,n)=0
2026 fluxs(m,n)=0
2027 fm(m,n)=0
2028 fluxm(m,n)=0
2029 end do
2030 end do
2031 return
2032 endif
2033
2034 w=wbar
2035 alw=alpha*wnuc
2036 sqrtalw=sqrt(alw)
2037 zeta=2.*sqrtalw*aten/(3.*sqrtg)
2038
2039 if (isectional .gt. 0) then
2040 ! sectional model.
2041 ! use bulk properties
2042 do n=1,ntype_aer
2043 if(totn(n).gt.1.e-10)then
2044 eta(1,n)=2*alw*sqrtalw/(totn(n)*beta*sqrtg)
2045 else
2046 eta(1,n)=1.e10
2047 endif
2048 end do
2049 call maxsat(zeta,eta,maxd_atype,ntype_aer,maxd_asize,(/1/),gmsm,gmlnsig,smax)
2050 ! call maxsat(zeta,eta,maxd_atype,ntype_aer,maxd_asize,nsize_aer,gmsm,gmlnsig,smax)
2051
2052 else
2053
2054 do n=1,ntype_aer
2055 do m=1,nsize_aer(n)
2056 if(na(m,n).gt.1.e-10)then
2057 eta(m,n)=2*alw*sqrtalw/(na(m,n)*beta*sqrtg)
2058 else
2059 eta(m,n)=1.e10
2060 endif
2061 end do
2062 end do
2063
2064 call maxsat(zeta,eta,maxd_atype,ntype_aer,maxd_asize,nsize_aer,sm,alnsign,smax)
2065
2066 endif
2067
2068 lnsmax=log(smax)
2069 xmincoeff=alogaten-2.*third*(lnsmax-alog2)-alog3
2070
2071 ! print *,'smax=',smax
2072
2073
2074 do 55002 n=1,ntype_aer
2075 do 55000 m=1,nsize_aer(n)
2076
2077 ! 04-nov-2005 rce - check for bin_is_empty here too, just like earlier
2078 if ( bin_is_empty(m,n) ) then
2079 fn(m,n)=0.
2080 fs(m,n)=0.
2081 fm(m,n)=0.
2082
2083 else if ((isectional .eq. 2) .or. (isectional .eq. 4)) then
2084 ! sectional
2085 ! within-section dn/dx = a + b*x
2086 xcut=xmincoeff-third*lnhygro(m,n)
2087 ! ycut=(exp(xcut)/rhi(m,n))**3
2088 ! 07-jul-2006 rce - the above line gave a (rare) overflow when smax=1.0e-20
2089 ! if (ycut > yhi), then actual value of ycut is unimportant,
2090 ! so do the following to avoid overflow
2091 lnycut = 3.0 * ( xcut - log(rhi(m,n)) )
2092 lnycut = min( lnycut, log(yhi(m,n)*1.0e5) )
2093 ycut=exp(lnycut)
2094 ! write(6,*)'m,n,rcut,rlo,rhi=',m,n,exp(xcut),rlo(m,n),rhi(m,n)
2095 ! if(lnsmax.lt.lnsmn(m,n))then
2096 if(ycut.gt.yhi(m,n))then
2097 fn(m,n)=0.
2098 fs(m,n)=0.
2099 fm(m,n)=0.
2100 ! elseif(lnsmax.gt.lnsmx(m,n))then
2101 elseif(ycut.lt.ylo(m,n))then
2102 fn(m,n)=1.
2103 fs(m,n)=1.
2104 fm(m,n)=1.
2105 elseif ( bin_is_narrow(m,n) ) then
2106 ! 04-nov-2005 rce - for extremely narrow bins,
2107 ! do zero activation if xcut>xmean, 100% activation otherwise
2108 if (ycut.gt.ymean(m,n)) then
2109 fn(m,n)=0.
2110 fs(m,n)=0.
2111 fm(m,n)=0.
2112 else
2113 fn(m,n)=1.
2114 fs(m,n)=1.
2115 fm(m,n)=1.
2116 endif
2117 else
2118 phiyy=ycut/yhi(m,n)
2119 fn(m,n) = asub(m,n)*(1.0-phiyy) + 0.5*bsub(m,n)*(1.0-phiyy*phiyy)
2120 if (fn(m,n).lt.zero .or. fn(m,n).gt.one) then
2121 if (idiag_fnsm_prob .gt. 0) then
2122 print *,'fn(',m,n,')=',fn(m,n),' outside 0,1 - activate err21'
2123 print *,'na,volc =', na(m,n), volc(m,n)
2124 print *,'asub,bsub =', asub(m,n), bsub(m,n)
2125 print *,'yhi,ycut =', yhi(m,n), ycut
2126 endif
2127 endif
2128
2129 if (fn(m,n) .le. zero) then
2130 ! 10-nov-2005 rce - if fn=0, then fs & fm must be 0
2131 fn(m,n)=zero
2132 fs(m,n)=zero
2133 fm(m,n)=zero
2134 else if (fn(m,n) .ge. one) then
2135 ! 10-nov-2005 rce - if fn=1, then fs & fm must be 1
2136 fn(m,n)=one
2137 fs(m,n)=one
2138 fm(m,n)=one
2139 else
2140 ! 10-nov-2005 rce - otherwise, calc fm and check it
2141 fm(m,n) = (yhi(m,n)/ymean(m,n)) * (0.5*asub(m,n)*(1.0-phiyy*phiyy) + &
2142 third*bsub(m,n)*(1.0-phiyy*phiyy*phiyy))
2143 if (fm(m,n).lt.fn(m,n) .or. fm(m,n).gt.one) then
2144 if (idiag_fnsm_prob .gt. 0) then
2145 print *,'fm(',m,n,')=',fm(m,n),' outside fn,1 - activate err22'
2146 print *,'na,volc,fn =', na(m,n), volc(m,n), fn(m,n)
2147 print *,'asub,bsub =', asub(m,n), bsub(m,n)
2148 print *,'yhi,ycut =', yhi(m,n), ycut
2149 endif
2150 endif
2151 if (fm(m,n) .le. fn(m,n)) then
2152 ! 10-nov-2005 rce - if fm=fn, then fs must =fn
2153 fm(m,n)=fn(m,n)
2154 fs(m,n)=fn(m,n)
2155 else if (fm(m,n) .ge. one) then
2156 ! 10-nov-2005 rce - if fm=1, then fs & fn must be 1
2157 fm(m,n)=one
2158 fs(m,n)=one
2159 fn(m,n)=one
2160 else
2161 ! 10-nov-2005 rce - these two checks assure that the mean size
2162 ! of the activated & interstitial particles will be between rlo & rhi
2163 dumaa = fn(m,n)*(yhi(m,n)/ymean(m,n))
2164 fm(m,n) = min( fm(m,n), dumaa )
2165 dumaa = 1.0 + (fn(m,n)-1.0)*(ylo(m,n)/ymean(m,n))
2166 fm(m,n) = min( fm(m,n), dumaa )
2167 ! 10-nov-2005 rce - now calculate fs and bound it by fn, fm
2168 betayy = ylo(m,n)/yhi(m,n)
2169 dumaa = phiyy**twothird
2170 dumbb = betayy**twothird
2171 fs(m,n) = &
2172 (asub(m,n)*(1.0-phiyy*dumaa) + &
2173 0.625*bsub(m,n)*(1.0-phiyy*phiyy*dumaa)) / &
2174 (asub(m,n)*(1.0-betayy*dumbb) + &
2175 0.625*bsub(m,n)*(1.0-betayy*betayy*dumbb))
2176 fs(m,n)=max(fs(m,n),fn(m,n))
2177 fs(m,n)=min(fs(m,n),fm(m,n))
2178 endif
2179 endif
2180
2181 endif
2182
2183 else
2184 ! modal
2185 x=2*(lnsm(m,n)-lnsmax)/(3*sq2*alnsign(m,n))
2186 fn(m,n)=0.5*(1.-ERF_ALT(x))
2187 arg=x-sq2*alnsign(m,n)
2188 fs(m,n)=0.5*(1.-ERF_ALT(arg))
2189 arg=x-1.5*sq2*alnsign(m,n)
2190 fm(m,n)=0.5*(1.-ERF_ALT(arg))
2191 endif
2192
2193 ! fn(m,n)=1. ! test
2194 ! fs(m,n)=1.
2195 ! fm(m,n)=1.
2196 if((top.and.wbar.lt.0.).or.(.not.top.and.wbar.gt.0.))then
2197 fluxn(m,n)=fn(m,n)*w
2198 fluxs(m,n)=fs(m,n)*w
2199 fluxm(m,n)=fm(m,n)*w
2200 else
2201 fluxn(m,n)=0
2202 fluxs(m,n)=0
2203 fluxm(m,n)=0
2204 endif
2205
2206 55000 continue ! m=1,nsize_aer(n)
2207 55002 continue ! n=1,ntype_aer
2208
2209 ! 04-nov-2005 rce - moved the code for "wnuc.le.0" from here
2210 ! to near the start the uniform undraft section
2211
2212 !.......................................................................
2213 !
2214 ! end calc. of activation fractions/fluxes
2215 ! for (single) uniform updraft (end of section 3)
2216 !
2217 !.......................................................................
2218
2219
2220
2221 60000 continue
2222
2223
2224 ! do n=1,ntype_aer
2225 ! do m=1,nsize_aer(n)
2226 ! write(6,'(a,2i3,5e10.1)')'n,m,na,wbar,sigw,fn,fm=',n,m,na(m,n),wbar,sigw,fn(m,n),fm(m,n)
2227 ! end do
2228 ! end do
2229
2230
2231 return
2232 end subroutine activate
2233
2234
2235
2236 !----------------------------------------------------------------------
2237 !----------------------------------------------------------------------
2238 subroutine maxsat(zeta,eta,maxd_atype,ntype_aer,maxd_asize,nsize_aer, &
2239 sm,alnsign,smax)
2240
2241 ! calculates maximum supersaturation for multiple
2242 ! competing aerosol modes.
2243
2244 ! Abdul-Razzak and Ghan, A parameterization of aerosol activation.
2245 ! 2. Multiple aerosol types. J. Geophys. Res., 105, 6837-6844.
2246
2247 integer maxd_atype
2248 integer ntype_aer
2249 integer maxd_asize
2250 integer nsize_aer(maxd_atype) ! number of size bins
2251 real sm(maxd_asize,maxd_atype) ! critical supersaturation for number mode radius
2252 real zeta, eta(maxd_asize,maxd_atype)
2253 real alnsign(maxd_asize,maxd_atype) ! ln(sigma)
2254 integer pmax
2255 parameter (pmax=100)
2256 real f1(pmax,pmax)
2257 real smax ! maximum supersaturation
2258 save f1
2259 logical first
2260 data first/.true./
2261 save first
2262 real twothird,sum
2263 ! 04-nov-2005 rce - make this more precise
2264 ! data twothird/0.666666666/
2265 data twothird/0.66666666667/
2266 save twothird
2267 integer m ! size index
2268 integer n ! type index
2269
2270 if(first)then
2271 ! calculate and save f1(sigma). assumes sigma is invariant.
2272 do n=1,ntype_aer
2273 do m=1,nsize_aer(n)
2274 if(ntype_aer>pmax)then
2275 print *,'pmax < ',ntype_aer,' in maxsat'
2276 call exit
2277 endif
2278 if(nsize_aer(n)>pmax)then
2279 print *,'pmax < ',nsize_aer(n),' in maxsat'
2280 call exit
2281 endif
2282 f1(m,n)=0.5*exp(2.5*alnsign(m,n)*alnsign(m,n))
2283 end do
2284 end do
2285 first=.false.
2286 endif
2287
2288 do n=1,ntype_aer
2289 do m=1,nsize_aer(n)
2290 if(zeta.gt.1.e5*eta(m,n).or.sm(m,n)*sm(m,n).gt.1.e5*eta(m,n))then
2291 ! weak forcing. essentially none activated
2292 smax=1.e-20
2293 else
2294 ! significant activation of this mode. calc activation all modes.
2295 go to 1
2296 endif
2297 end do
2298 end do
2299
2300 return
2301
2302 1 continue
2303
2304 sum=0
2305 do n=1,ntype_aer
2306 do m=1,nsize_aer(n)
2307 if(eta(m,n).gt.1.e-20)then
2308 g1=sqrt(zeta/eta(m,n))
2309 g1=g1*g1*g1
2310 g2=sm(m,n)/sqrt(eta(m,n)+3*zeta)
2311 g2=sqrt(g2)
2312 g2=g2*g2*g2
2313 sum=sum+(f1(m,n)*g1+(1.+0.25*alnsign(m,n))*g2)/(sm(m,n)*sm(m,n))
2314 else
2315 sum=1.e20
2316 endif
2317 end do
2318 end do
2319
2320 smax=1./sqrt(sum)
2321
2322 return
2323 end subroutine maxsat
2324
2325
2326
2327
2328 !----------------------------------------------------------------------
2329 !----------------------------------------------------------------------
2330 ! 25-apr-2006 rce - dens_aer is (g/cm3), NOT (kg/m3);
2331 ! grid_id, ktau, i, j, isize, itype added to arg list to assist debugging
2332 subroutine loadaer(chem,k,kmn,kmx,num_chem,cs,npv, &
2333 dlo_sect,dhi_sect,maxd_acomp, ncomp, &
2334 grid_id, ktau, i, j, isize, itype, &
2335 numptr_aer, numptrcw_aer, dens_aer, &
2336 massptr_aer, massptrcw_aer, &
2337 maerosol, maerosolcw, &
2338 maerosol_tot, maerosol_totcw, &
2339 naerosol, naerosolcw, &
2340 vaerosol, vaerosolcw)
2341
2342 implicit none
2343
2344 ! load aerosol number, surface, mass concentrations
2345
2346 ! input
2347
2348 integer num_chem ! maximum number of consituents
2349 integer k,kmn,kmx
2350 real chem(kmn:kmx,num_chem) ! aerosol mass, number mixing ratios
2351 real cs ! air density (kg/m3)
2352 real npv ! number per volume concentration (/m3)
2353 integer maxd_acomp,ncomp
2354 integer numptr_aer,numptrcw_aer
2355 integer massptr_aer(maxd_acomp), massptrcw_aer(maxd_acomp)
2356 real dens_aer(maxd_acomp) ! aerosol material density (g/cm3)
2357 real dlo_sect,dhi_sect ! minimum, maximum diameter of section (cm)
2358 integer grid_id, ktau, i, j, isize, itype
2359
2360 ! output
2361
2362 real naerosol ! interstitial number conc (/m3)
2363 real naerosolcw ! activated number conc (/m3)
2364 real maerosol(maxd_acomp) ! interstitial mass conc (kg/m3)
2365 real maerosolcw(maxd_acomp) ! activated mass conc (kg/m3)
2366 real maerosol_tot ! total-over-species interstitial mass conc (kg/m3)
2367 real maerosol_totcw ! total-over-species activated mass conc (kg/m3)
2368 real vaerosol ! interstitial volume conc (m3/m3)
2369 real vaerosolcw ! activated volume conc (m3/m3)
2370
2371 ! internal
2372
2373 integer lnum,lnumcw,l,ltype,lmass,lmasscw,lsfc,lsfccw
2374 real num_at_dhi, num_at_dlo
2375 real npv_at_dhi, npv_at_dlo
2376 real pi
2377 real specvol ! inverse aerosol material density (m3/kg)
2378 ! 04-nov-2005 rce - make this more precise
2379 ! data pi/3.14159/
2380 data pi/3.1415926526/
2381 save pi
2382
2383
2384 lnum=numptr_aer
2385 lnumcw=numptrcw_aer
2386 maerosol_tot=0.
2387 maerosol_totcw=0.
2388 vaerosol=0.
2389 vaerosolcw=0.
2390 do l=1,ncomp
2391 lmass=massptr_aer(l)
2392 lmasscw=massptrcw_aer(l)
2393 maerosol(l)=chem(k,lmass)*cs
2394 maerosol(l)=max(maerosol(l),0.)
2395 maerosolcw(l)=chem(k,lmasscw)*cs
2396 maerosolcw(l)=max(maerosolcw(l),0.)
2397 maerosol_tot=maerosol_tot+maerosol(l)
2398 maerosol_totcw=maerosol_totcw+maerosolcw(l)
2399 ! [ 1.e-3 factor because dens_aer is (g/cm3), specvol is (m3/kg) ]
2400 specvol=1.0e-3/dens_aer(l)
2401 vaerosol=vaerosol+maerosol(l)*specvol
2402 vaerosolcw=vaerosolcw+maerosolcw(l)*specvol
2403 ! write(6,'(a,3e12.2)')'maerosol,dens_aer,vaerosol=',maerosol(l),dens_aer(l),vaerosol
2404 enddo
2405
2406 if(lnum.gt.0)then
2407 ! aerosol number predicted
2408 ! [ 1.0e6 factor because because dhi_ & dlo_sect are (cm), vaerosol is (m3) ]
2409 npv_at_dhi = 6.0e6/(pi*dhi_sect*dhi_sect*dhi_sect)
2410 npv_at_dlo = 6.0e6/(pi*dlo_sect*dlo_sect*dlo_sect)
2411
2412 naerosol=chem(k,lnum)*cs
2413 naerosolcw=chem(k,lnumcw)*cs
2414 num_at_dhi = vaerosol*npv_at_dhi
2415 num_at_dlo = vaerosol*npv_at_dlo
2416 naerosol = max( num_at_dhi, min( num_at_dlo, naerosol ) )
2417 ! write(6,'(a,5e10.1)')'naerosol,num_at_dhi,num_at_dlo,dhi_sect,dlo_sect', &
2418 ! naerosol,num_at_dhi,num_at_dlo,dhi_sect,dlo_sect
2419 num_at_dhi = vaerosolcw*npv_at_dhi
2420 num_at_dlo = vaerosolcw*npv_at_dlo
2421 naerosolcw = max( num_at_dhi, min( num_at_dlo, naerosolcw ) )
2422 else
2423 ! aerosol number diagnosed from mass and prescribed size
2424 naerosol=vaerosol*npv
2425 naerosol=max(naerosol,0.)
2426 naerosolcw=vaerosolcw*npv
2427 naerosolcw=max(naerosolcw,0.)
2428 endif
2429
2430
2431 return
2432 end subroutine loadaer
2433
2434
2435
2436 !-----------------------------------------------------------------------
2437 real function erfc_num_recipes( x )
2438 !
2439 ! from press et al, numerical recipes, 1990, page 164
2440 !
2441 implicit none
2442 real x
2443 double precision erfc_dbl, dum, t, zz
2444
2445 zz = abs(x)
2446 t = 1.0/(1.0 + 0.5*zz)
2447
2448 ! erfc_num_recipes =
2449 ! & t*exp( -zz*zz - 1.26551223 + t*(1.00002368 + t*(0.37409196 +
2450 ! & t*(0.09678418 + t*(-0.18628806 + t*(0.27886807 +
2451 ! & t*(-1.13520398 +
2452 ! & t*(1.48851587 + t*(-0.82215223 + t*0.17087277 )))))))))
2453
2454 dum = ( -zz*zz - 1.26551223 + t*(1.00002368 + t*(0.37409196 + &
2455 t*(0.09678418 + t*(-0.18628806 + t*(0.27886807 + &
2456 t*(-1.13520398 + &
2457 t*(1.48851587 + t*(-0.82215223 + t*0.17087277 )))))))))
2458
2459 erfc_dbl = t * exp(dum)
2460 if (x .lt. 0.0) erfc_dbl = 2.0d0 - erfc_dbl
2461
2462 erfc_num_recipes = erfc_dbl
2463
2464 return
2465 end function erfc_num_recipes
2466
2467 !-----------------------------------------------------------------------
2468 real function erf_alt( x )
2469
2470 implicit none
2471
2472 real,intent(in) :: x
2473
2474 erf_alt = 1. - erfc_num_recipes(x)
2475
2476 end function erf_alt
2477
2478 END MODULE module_mixactivate