source: trunk/WRF.COMMON/WRFV3/phys/module_mixactivate.F @ 3553

Last change on this file since 3553 was 2759, checked in by aslmd, 2 years ago

adding unmodified code from WRFV3.0.1.1, expurged from useless data +1M size

File size: 96.2 KB
RevLine 
[2759]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
10MODULE module_mixactivate
11PRIVATE
12PUBLIC prescribe_aerosol_mixactivate, mixactivate
13CONTAINS
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)
205subroutine 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
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
1014100  continue ! end of main loop over i
1015120  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
155325000        continue      ! m=1,nsize_aer(n)
155425002        continue      ! n=1,ntype_aer
1555
1556
155730000    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
158340000    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
164994700    format( / 'activate 47000 - ktau,id,ii,jj,kk,nwmax=', 6i5 )
165094710    format( 'activate 47000 - ', a, 6(1x,f11.5) )
165194720    format( 'activate 47000 - nw,w,wold,dw=', i5, 3(1x,f11.5) )
1652
1653         do 47000 nw = 1, nwmax
165441000       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
186944000       continue      ! m=1,nsize_aer(n)
187044002       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
188047000    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
191448000    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!.......................................................................
200950000 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
219955000       continue      ! m=1,nsize_aer(n)
220055002       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
221460000 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
2483END MODULE module_mixactivate
Note: See TracBrowser for help on using the repository browser.