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