source: lmdz_wrf/trunk/WRFV3/phys/module_mixactivate.F @ 1544

Last change on this file since 1544 was 1, checked in by lfita, 10 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

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