module ndrop 2,3
#if (defined MODAL_AERO)
use shr_kind_mod
, only: r8 => shr_kind_r8
use abortutils, only: endrun
use modal_aero_data, only: maxd_amode
use cam_logfile, only: iulog
implicit none
private
save
public activate_init, dropmixnuc
real(r8) :: npv(maxd_amode) ! number per volume concentration
real(r8) :: alogsig(maxd_amode) ! natl log of geometric standard dev of aerosol
real(r8) :: exp45logsig(maxd_amode)
real(r8) :: argfactor(maxd_amode)
real(r8) :: f1(maxd_amode),f2(maxd_amode) ! abdul-razzak functions of width
real(r8) :: t0 ! reference temperature
real(r8) :: aten
real(r8) :: surften ! surface tension of water w/respect to air (N/m)
real(r8) :: alogten,alog2,alog3,alogaten
real(r8) :: third, twothird, sixth, zero
real(r8) :: sq2, sqpi, pi
contains
subroutine dropmixnuc(lchnk, ncol, ncldwtr,tendnd, temp,omega, & 1,30
pmid,pint,pdel,rpdel,zm,kvh,wsub,cldn,cldo, &
raer, cflx, raertend, qqcw, dtmicro)
! vertical diffusion and nucleation of cloud droplets
! assume cloud presence controlled by cloud fraction
! doesn't distinguish between warm, cold clouds
use ppgrid, only: pcols, pver, pverp
use physconst, only: gravit, rhoh2o, latvap, cpair, epsilo, rair, mwh2o, r_universal
use time_manager, only: get_nstep
use constituents, only: pcnst, qmin, cnst_get_ind, cnst_name
use error_function, only: erf
! use aerosol_radiation_interface, only: collect_sw_aer_masses, aer_mass
use modal_aero_data
use cam_history, only: fieldname_len
use phys_grid, only: get_lat_all_p, get_lon_all_p
use physics_types, only: physics_state
use cam_history, only: outfld
implicit none
! input
integer, intent(in) :: lchnk ! chunk identifier
integer, intent(in) :: ncol ! number of columns
! type(physics_state), intent(in) :: state ! Physics state variables
real(r8), intent(in) :: dtmicro ! time step for microphysics (s)
real(r8), intent(in) :: temp(pcols,pver) ! temperature (K)
real(r8), intent(inout) :: ncldwtr(pcols,pver) ! droplet number concentration (#/kg)
real(r8), intent(inout) :: tendnd(pcols,pver) ! change in droplet number concentration (#/kg/s)
real(r8), intent(in) :: omega(pcols,pver) ! vertical velocity (Pa/s)
real(r8), intent(in) :: pmid(pcols,pver) ! mid-level pressure (Pa)
real(r8), intent(in) :: pint(pcols,pverp) ! pressure at layer interfaces (Pa)
real(r8), intent(in) :: pdel(pcols,pver) ! pressure thickess of layer (Pa)
real(r8), intent(in) :: rpdel(pcols,pver) ! inverse of pressure thickess of layer (/Pa)
real(r8), intent(in) :: zm(pcols,pver) ! geopotential height of level (m)
real(r8), intent(in) :: kvh(pcols,pverp) ! vertical diffusivity (m2/s)
real(r8), intent(in) :: wsub(pcols,pver) ! subgrid vertical velocity
real(r8), intent(in) :: cldo(pcols,pver) ! cloud fraction on previous time step
real(r8), intent(in) :: cldn(pcols,pver) ! cloud fraction
! output
real(r8), intent(in) :: raer(pcols,pver,pcnst) ! aerosol mass, number mixing ratios
real(r8), intent(in) :: cflx(pcols,pcnst) ! constituent flux from surface
real(r8), intent(out) :: raertend(pcols,pver,pcnst) ! tendency of aerosol mass, number mixing ratios
real(r8), intent(inout) :: qqcw(pcols,pver,pcnst) ! cloud-borne aerosol mass, number mixing ratios
!--------------------Local storage-------------------------------------
!
real(r8) depvel(pcols,pcnst)! deposition velocity for droplets (m/s)
real(r8) qcld(pver) ! cloud droplet number mixing ratio (#/kg)
real(r8) wtke(pcols,pver) ! turbulent vertical velocity at base of layer k (m/s)
real(r8) wtke_cen(pcols,pver) ! turbulent vertical velocity at center of layer k (m/s)
real(r8) zn(pver) ! g/pdel (m2/g) for layer
real(r8) zs(pver) ! inverse of distance between levels (m)
real(r8), parameter :: zkmin=0.01_r8,zkmax=100._r8
real(r8) w(pcols,pver) ! vertical velocity (m/s)
real(r8) cs(pcols,pver) ! air density (kg/m3)
real(r8) dz(pcols,pver) ! geometric thickness of layers (m)
real(r8) zero,flxconv ! convergence of flux into lowest layer
real(r8) wdiab ! diabatic vertical velocity
real(r8), parameter :: wmixmin = 0.1 ! minimum turbulence vertical velocity (m/s)
! real(r8), parameter :: wmixmin = 0.2 ! minimum turbulence vertical velocity (m/s)
! real(r8), parameter :: wmixmin = 1.0 ! minimum turbulence vertical velocity (m/s)
real(r8) qncld(pver) ! droplet number nucleated on cloud boundaries
real(r8) ekd(pver) ! diffusivity for droplets (m2/s)
real(r8) ekk(0:pver) ! density*diffusivity for droplets (kg/m3 m2/s)
real(r8) srcn(pver) ! droplet source rate (/s)
real(r8), parameter :: sq2pi=2.5066283_r8
real(r8) dtinv
logical top ! true if cloud top, false if cloud base or new cloud
integer km1,kp1
real(r8) wbar,wmix,wmin,wmax
real(r8) dum,dumc
real(r8) dact
real(r8) fluxntot ! (#/cm2/s)
real(r8) fac_srflx
real(r8) surfrate(pcnst) ! surface exchange rate (/s)
real(r8) surfratemax ! max surfrate for all species treated here
real(r8) dtmin,tinv,dtt
integer nsubmix,nsubmix_bnd
integer i,k,m,n
real(r8) tempr4 ! temperature
real(r8) dtmix
real(r8) alogarg
integer nstep
real(r8) pi
integer nnew,nsav,ntemp
real(r8) overlapp(pver),overlapm(pver) ! cloud overlap
real(r8) ekkp(pver),ekkm(pver) ! zn*zs*density*diffusivity
integer count_submix(100)
save count_submix
real(r8) kvhmax
save kvhmax
real(r8) nsource(pcols,pver) ! droplet number source (#/kg/s)
real(r8) ndropmix(pcols,pver) ! droplet number mixing (#/kg/s)
real(r8) ndropcol(pcols) ! column droplet number (#/m2)
integer lnum, lnumcw, lmass, lmasscw, lphase, &
lsfc, lsfccw, lsig, lspec, ltype, lwater
integer ntype(maxd_amode)
real(r8) na(pcols),va(pcols),hy(pcols)
real(r8) naermod(maxd_amode) ! (/m3)
real(r8) sigmag(maxd_amode) ! geometric standard deviation of aerosol size dist
real(r8) hygro(maxd_amode) ! hygroscopicity of aerosol mode
real(r8) vaerosol(maxd_amode) ! interstit+activated aerosol volume conc (cm3/cm3)
real(r8) raercol(pver,pcnst,2) ! aerosol mass, number mixing ratios
real(r8) raercol_cw(pver,pcnst,2) ! cloud-borne aerosol mass, number mixing ratios
real(r8) surfrate_cw(pcnst) ! surface exchange rate for cloud-borne aerosol (/s)
real(r8) source(pver) !
real(r8) fn(maxd_amode) ! activation fraction for aerosol number
real(r8) fm(maxd_amode) ! activation fraction for aerosol mass
real(r8) fluxn(maxd_amode) ! number activation fraction flux (cm/s)
real(r8) fluxm(maxd_amode) ! mass activation fraction flux (cm/s)
real(r8) sum,sumcw,flux,fluxcw
! note: activation fraction fluxes are defined as
! fluxn = [flux of activated aero. number into cloud (#/cm2/s)]
! / [aero. number conc. in updraft, just below cloudbase (#/cm3)]
real(r8) nact(pver,maxd_amode) ! fractional aero. number activation rate (/s)
real(r8) mact(pver,maxd_amode) ! fractional aero. mass activation rate (/s)
real(r8) sigmag_amode_cur(maxd_amode) ! sigmag for current grid cell
save sigmag_amode_cur
real(r8) :: qqcwtend(pcols,pver,pcnst) ! tendency of cloudborne aerosol mass, number mixing ratios
real(r8) :: coltend(pcols) ! column tendency
real(r8) :: tmpa
integer :: lc
character(len=fieldname_len) :: tmpname
character(len=fieldname_len+3) :: fieldname
real(r8) :: csbot(pver) ! air density at bottom (interface) of layer (kg/m3)
real(r8) :: csbot_cscen(pver) ! csbot(i)/cs(i,k)
real(r8) :: flux_fullact(pver) ! 100% activation fraction flux (cm/s)
real(r8) :: taumix_internal_pver_inv ! 1/(internal mixing time scale for k=pver) (1/s)
real(r8) dactn,taunuc,damp
real(r8) :: cldo_tmp, cldn_tmp
real(r8) :: tau_cld_regenerate
logical cldinterior
integer ixndrop, l
real(r8) naer_tot,naer(pcols)
integer, parameter :: psat=6 ! number of supersaturations to calc ccn concentration
real(r8) :: supersat(psat)= & ! supersaturation (%) to determine ccn concentration
(/0.02,0.05,0.1,0.2,0.5,1.0/)
real(r8) ccn(pcols,pver,psat) ! number conc of aerosols activated at supersat
character(len=8), dimension(psat) :: ccn_name(psat)= &
(/'CCN1','CCN2','CCN3','CCN4','CCN5','CCN6'/)
real(r8) arg
integer phase ! phase of aerosol
arg = 1.0_r8
if (abs(0.8427_r8-erf(arg))/0.8427_r8>0.001_r8) then
write (iulog,*) 'erf(1.0) = ',ERF(arg)
call endrun('dropmixnuc: Error function error')
endif
arg = 0.0
if (erf(arg) /= 0.0) then
write (iulog,*) 'erf(0.0) = ',erf(arg)
write (iulog,*) 'dropmixnuc: Error function error'
call endrun('dropmixnuc: Error function error')
endif
zero=0._r8
pi = 4._r8*atan(1.0_r8)
dtinv=1./dtmicro
nstep = get_nstep()
! call t_startf('load_aerosol')
call cnst_get_ind('NUMLIQ', ixndrop)
depvel(:,:) = 0.0_r8 ! droplet number is done in pkg_cld_sediment, aerosols in mz_aerosols_intr
! depvel(:,ixndrop) = 0.1_r8 ! prescribed here rather than getting it from MIRAGE depvel_part
! call t_stopf('load_aerosol')
! start loop over columns
overall_main_i_loop: &
do i=1,ncol
! load number nucleated into qcld on cloud boundaries
! initialization for current i .........................................
! call t_startf('calc_wtke')
do k=1,pver-1
zs(k)=1._r8/(zm(i,k)-zm(i,k+1))
enddo
zs(pver)=zs(pver-1)
do k=1,pver
qcld(k)=ncldwtr(i,k)
qncld(k)=0._r8
srcn(k)=0._r8
cs(i,k)=pmid(i,k)/(rair*temp(i,k)) ! air density (kg/m3)
dz(i,k)=1._r8/(cs(i,k)*gravit*rpdel(i,k)) ! layer thickness in m
w(i,k)=-1._r8*omega(i,k)/(cs(i,k)*gravit) ! large-scale vertical velocity m/s
do m=1,ntot_amode
nact(k,m)=0._r8
mact(k,m)=0._r8
enddo
zn(k)=gravit*rpdel(i,k)
kvhmax=max(kvh(i,k),kvhmax)
if(k<pver)then
ekd(k)=kvh(i,k+1)
ekd(k)=max(ekd(k),zkmin)
ekd(k)=min(ekd(k),zkmax)
csbot(k)=2.0_r8*pint(i,k+1)/(rair*(temp(i,k)+temp(i,k+1)))
csbot_cscen(k) = csbot(k)/cs(i,k)
else
ekd(k)=0._r8
csbot(k)=cs(i,k)
csbot_cscen(k) = 1.0_r8
endif
! diagnose subgrid vertical velocity from diffusivity
if(k.eq.pver)then
wtke(i,k)=sq2pi*depvel(i,ixndrop)
! wtke(i,k)=sq2pi*kvh(i,k+1)
wtke(i,k)=max(wtke(i,k),wmixmin)
else
wtke(i,k)=sq2pi*ekd(k)/dz(i,k)
endif
! get sub-grid vertical velocity from diff coef.
! following morrison et al. 2005, JAS
! assume mixing length of 30 m
! rce-comment - define wtke at layer centers for new-cloud activation
! and at layer boundaries for old-cloud activation
!++ag
! wtke_cen(i,k)=0.5_r8*(kvh(i,k)+kvh(i,k+1))/30._r8
! wtke(i,k)=kvh(i,k+1)/30._r8
wtke_cen(i,k)=wsub(i,k)
wtke(i,k)=wsub(i,k)
!--ag
wtke_cen(i,k)=max(wtke_cen(i,k),wmixmin)
wtke(i,k)=max(wtke(i,k),wmixmin)
nsource(i,k)=0._r8
enddo
surfratemax = 0.0_r8
nsav=1
nnew=2
surfrate(ixndrop)=depvel(i,ixndrop)/dz(i,pver)
surfratemax = max( surfratemax, surfrate(ixndrop) )
do m=1,ntot_amode
lnum=numptr_amode(m)
lnumcw=numptrcw_amode(m)
if(lnum>0)then
surfrate(lnum)=depvel(i,lnum)/dz(i,pver)
! surfrate(lnumcw)=surfrate(ixndrop)
surfrate_cw(lnumcw)=surfrate(ixndrop)
surfratemax = max( surfratemax, surfrate(lnum) )
! raercol(:,lnumcw,nsav)=raer(i,:,lnumcw)
raercol_cw(:,lnumcw,nsav)=qqcw(i,:,lnumcw) ! use cloud-borne aerosol array
raercol(:,lnum,nsav)=raer(i,:,lnum)
endif
do l=1,nspec_amode(m)
lmass=lmassptr_amode(l,m)
lmasscw=lmassptrcw_amode(l,m)
surfrate(lmass)=depvel(i,lmass)/dz(i,pver)
! surfrate(lmasscw)=surfrate(ixndrop)
surfrate_cw(lmasscw)=surfrate(ixndrop)
surfratemax = max( surfratemax, surfrate(lmass) )
! raercol(:,lmasscw,nsav)=raer(i,:,lmasscw)
raercol_cw(:,lmasscw,nsav)=qqcw(i,:,lmasscw) ! use cloud-borne aerosol array
raercol(:,lmass,nsav)=raer(i,:,lmass)
enddo
! lwater=lwaterptr_amode(m)
! surfrate(lwater)=depvel(i,lwater)/dz(i,pver)
! surfratemax = max( surfratemax, surfrate(lwater) )
! raercol(:,lwater,nsav)=raer(i,:,lwater)
enddo
! call t_stopf('calc_wtke')
! droplet nucleation/aerosol activation
! tau_cld_regenerate = time scale for regeneration of cloudy air
! by (horizontal) exchange with clear air
tau_cld_regenerate = 3600.0_r8 * 3.0_r8
! k-loop for growing/shrinking cloud calcs .............................
grow_shrink_main_k_loop: &
do k=1,pver
km1=max0(k-1,1)
kp1=min0(k+1,pver)
! shrinking cloud ......................................................
! treat the reduction of cloud fraction from when cldn(i,k) < cldo(i,k)
! and also dissipate the portion of the cloud that will be regenerated
cldo_tmp = cldo(i,k)
cldn_tmp = cldn(i,k) * exp( -dtmicro/tau_cld_regenerate )
! alternate formulation
! cldn_tmp = cldn(i,k) * max( 0.0_r8, (1.0_r8-dtmicro/tau_cld_regenerate) )
if(cldn_tmp.lt.cldo_tmp)then
! droplet loss in decaying cloud
!++ sungsup
nsource(i,k)=nsource(i,k)+qcld(k)*(cldn_tmp-cldo_tmp)/cldo_tmp*dtinv
qcld(k)=qcld(k)*(1.+(cldn_tmp-cldo_tmp)/cldo_tmp)
!-- sungsup
! convert activated aerosol to interstitial in decaying cloud
dumc=(cldn_tmp-cldo_tmp)/cldo_tmp
do m=1,ntot_amode
lnum=numptr_amode(m)
lnumcw=numptrcw_amode(m)
if(lnum.gt.0)then
dact=raercol_cw(k,lnumcw,nsav)*dumc
! raercol(k,lnumcw,nsav)=raercol(k,lnumcw,nsav)+dact
raercol_cw(k,lnumcw,nsav)=raercol_cw(k,lnumcw,nsav)+dact ! cloud-borne aerosol
raercol(k,lnum,nsav)=raercol(k,lnum,nsav)-dact
endif
do l=1,nspec_amode(m)
lmass=lmassptr_amode(l,m)
lmasscw=lmassptrcw_amode(l,m)
dact=raercol_cw(k,lmasscw,nsav)*dumc
! raercol(k,lmasscw,nsav)=raercol(k,lmasscw,nsav)+dact
raercol_cw(k,lmasscw,nsav)=raercol_cw(k,lmasscw,nsav)+dact ! cloud-borne aerosol
raercol(k,lmass,nsav)=raercol(k,lmass,nsav)-dact
enddo
enddo
endif
! growing cloud ......................................................
! treat the increase of cloud fraction from when cldn(i,k) > cldo(i,k)
! and also regenerate part of the cloud
cldo_tmp = cldn_tmp
cldn_tmp = cldn(i,k)
if(cldn_tmp-cldo_tmp.gt.0.01)then
! wmix=wtke(i,k)
! wbar=wtke(i,k)
! rce-comment - use wtke at layer centers for new-cloud activation
wbar=wtke_cen(i,k)
wmix=0._r8
wmin=0._r8
wmax=10._r8
wdiab=0
! load aerosol properties, assuming external mixtures
phase=1 ! interstitial
do m=1,ntot_amode
call loadaer(raer,qqcw,i,i,k,m,cs,npv(m),phase, &
na, va, hy )
naermod(m)=na(i)
vaerosol(m)=va(i)
hygro(m)=hy(i)
end do
! m=1
! naermod(m)=1000.e6
! vaerosol(m)=naermod(m)/(density*num_to_mass_aer(m))
! call t_startf ('activate')
call activate_modal(wbar,wmix,wdiab,wmin,wmax,temp(i,k),cs(i,k), &
naermod, ntot_amode,ntot_amode, vaerosol, sigmag,hygro, &
fn,fm,fluxn,fluxm,flux_fullact(k))
! call t_stopf ('activate')
! write(iulog,'(a,15f8.2)')'wtke,naer,fn=',wtke(i,k),naer_tot*1.e-6,(fn(m),m=1,ntot_amode)
dumc=(cldn_tmp-cldo_tmp)
do m=1,ntot_amode
lnum=numptr_amode(m)
lnumcw=numptrcw_amode(m)
dact=dumc*fn(m)*raer(i,k,lnum) ! interstitial only
qcld(k)=qcld(k)+dact
nsource(i,k)=nsource(i,k)+dact*dtinv
if(lnum.gt.0)then
raercol_cw(k,lnumcw,nsav)=raercol_cw(k,lnumcw,nsav)+dact ! cloud-borne aerosol
raercol(k,lnum,nsav)=raercol(k,lnum,nsav)-dact
endif
dum=dumc*fm(m)
do l=1,nspec_amode(m)
lmass=lmassptr_amode(l,m)
lmasscw=lmassptrcw_amode(l,m)
dact=dum*(raer(i,k,lmass)) ! interstitial only
raercol_cw(k,lmasscw,nsav)=raercol_cw(k,lmasscw,nsav)+dact ! cloud-borne aerosol
raercol(k,lmass,nsav)=raercol(k,lmass,nsav)-dact
enddo
enddo
endif
enddo grow_shrink_main_k_loop
! end of k-loop for growing/shrinking cloud calcs ......................
! ......................................................................
! start of k-loop for calc of old cloud activation tendencies ..........
!
! rce-comment
! changed this part of code to use current cloud fraction (cldn) exclusively
! consider case of cldo(:)=0, cldn(k)=1, cldn(k+1)=0
! previous code (which used cldo below here) would have no cloud-base activation
! into layer k. however, activated particles in k mix out to k+1,
! so they are incorrectly depleted with no replacement
!
old_cloud_main_k_loop: &
do k=1,pver
km1=max0(k-1,1)
kp1=min0(k+1,pver)
taumix_internal_pver_inv = 0.0
if(cldn(i,k).gt.0.01)then
! rce-comment
! if(cldo(i,k).gt.0.01)then ! with cldo changed to cldn this is redundant
wdiab=0
! wmix=wtke(i,k) ! spectrum of updrafts
! wbar=w(i,k) ! spectrum of updrafts
wmix=0._r8 ! single updraft
wbar=wtke(i,k) ! single updraft
! rce-comment - different treatment for k=pver
if (k == pver) wbar=wtke_cen(i,k) ! single updraft
wmax=10._r8
wmin=0._r8
! old cloud
if(cldn(i,k)-cldn(i,kp1).gt.0.01.or.k.eq.pver)then
! cloud base
ekd(k)=wtke(i,k)*dz(i,k)/sq2pi
! rce-comments
! first, should probably have 1/zs(k) here rather than dz(i,k) because
! the turbulent flux is proportional to ekd(k)*zs(k),
! while the dz(i,k) is used to get flux divergences
! and mixing ratio tendency/change
! second and more importantly, using a single updraft velocity here
! means having monodisperse turbulent updraft and downdrafts.
! The sq2pi factor assumes a normal draft spectrum.
! The fluxn/fluxm from activate must be consistent with the
! fluxes calculated in explmix.
ekd(k)=wbar/zs(k)
alogarg=max(1.e-20_r8,1/cldn(i,k)-1._r8)
wmin=wbar+wmix*0.25*sq2pi*log(alogarg)
phase=1 ! interstitial
do m=1,ntot_amode
! rce-comment - use kp1 here as old-cloud activation involves
! aerosol from layer below
! call loadaer(raer,qqcw,i,i,k,m,cs, npv(m),phase, &
call loadaer(raer,qqcw,i,i,kp1,m,cs, npv(m),phase, &
na, va, hy )
naermod(m)=na(i)
vaerosol(m)=va(i)
hygro(m)=hy(i)
end do
! m=1
! naermod(m)=1000.e6
! vaerosol(1,m)=naermod(m)/(density*num_to_mass_aer(m))
! call t_startf ('activate')
call activate_modal(wbar,wmix,wdiab,wmin,wmax,temp(i,k),cs(i,k), &
naermod, ntot_amode,ntot_amode, vaerosol, sigmag, hygro, &
fn,fm,fluxn,fluxm,flux_fullact(k))
! call t_stopf ('activate')
! write(iulog,'(a,15f8.2)')'wtke,naer,fn=',wtke(i,k),naer_tot*1.e-6,(fn(m),m=1,ntot_amode)
if(k.lt.pver)then
dumc = cldn(i,k)-cldn(i,kp1)
else
dumc=cldn(i,k)
endif
fluxntot=0
! rce-comment 1
! flux of activated mass into layer k (in kg/m2/s)
! = "actmassflux" = dumc*fluxm*raercol(kp1,lmass)*csbot(k)
! source of activated mass (in kg/kg/s) = flux divergence
! = actmassflux/(cs(i,k)*dz(i,k))
! so need factor of csbot_cscen = csbot(k)/cs(i,k)
! dum=1./(dz(i,k))
dum=csbot_cscen(k)/(dz(i,k))
! rce-comment 2
! code for k=pver was changed to use the following conceptual model
! in k=pver, there can be no cloud-base activation unless one considers
! a scenario such as the layer being partially cloudy,
! with clear air at bottom and cloudy air at top
! assume this scenario, and that the clear/cloudy portions mix with
! a timescale taumix_internal = dz(i,pver)/wtke_cen(i,pver)
! in the absence of other sources/sinks, qact (the activated particle
! mixratio) attains a steady state value given by
! qact_ss = fcloud*fact*qtot
! where fcloud is cloud fraction, fact is activation fraction,
! qtot=qact+qint, qint is interstitial particle mixratio
! the activation rate (from mixing within the layer) can now be
! written as
! d(qact)/dt = (qact_ss - qact)/taumix_internal
! = qtot*(fcloud*fact*wtke/dz) - qact*(wtke/dz)
! note that (fcloud*fact*wtke/dz) is equal to the nact/mact
! also, d(qact)/dt can be negative. in the code below
! it is forced to be >= 0
!
! steve --
! you will likely want to change this. i did not really understand
! what was previously being done in k=pver
! in the cam3_5_3 code, wtke(i,pver) appears to be equal to the
! droplet deposition velocity which is quite small
! in the cam3_5_37 version, wtke is done differently and is much
! larger in k=pver, so the activation is stronger there
!
if (k == pver) then
taumix_internal_pver_inv = flux_fullact(k)/dz(i,k)
end if
do m=1,ntot_amode
fluxn(m)=fluxn(m)*dumc
fluxm(m)=fluxm(m)*dumc
lnum=numptr_amode(m)
nact(k,m)=nact(k,m)+fluxn(m)*dum
mact(k,m)=mact(k,m)+fluxm(m)*dum
if (k > pver) then
! note that kp1 is used here
fluxntot = fluxntot &
+fluxn(m)*raercol(kp1,lnum,nsav)*cs(i,k)
else
lnumcw=numptrcw_amode(m)
tmpa = raercol(kp1,lnum,nsav)*fluxn(m) &
+ raercol(kp1,lnumcw,nsav)*(fluxn(m) &
-taumix_internal_pver_inv)
fluxntot = fluxntot + max(0.0_r8,tmpa)*cs(i,k)
end if
enddo
srcn(k)=srcn(k)+fluxntot/(cs(i,k)*dz(i,k))
nsource(i,k)=nsource(i,k)+fluxntot/(cs(i,k)*dz(i,k))
endif ! (cldn(i,k)-cldn(i,kp1).gt.0.01 .or. k.eq.pver)
! rce-comment
! endif ! (cldo(i,k).gt.0.01) ! with cldo changed to cldn this is redundant
else
! no cloud
nsource(i,k)=nsource(i,k)-qcld(k)*dtinv
qcld(k)=0
! convert activated aerosol to interstitial in decaying cloud
do m=1,ntot_amode
lnum=numptr_amode(m)
lnumcw=numptrcw_amode(m)
if(lnum.gt.0)then
raercol(k,lnum,nsav)=raercol(k,lnum,nsav)+raercol_cw(k,lnumcw,nsav) ! cloud-borne aerosol
raercol_cw(k,lnumcw,nsav)=0.
endif
do l=1,nspec_amode(m)
lmass=lmassptr_amode(l,m)
lmasscw=lmassptrcw_amode(l,m)
raercol(k,lmass,nsav)=raercol(k,lmass,nsav)+raercol_cw(k,lmasscw,nsav) ! cloud-borne aerosol
raercol_cw(k,lmasscw,nsav)=0.
enddo
enddo
endif
enddo old_cloud_main_k_loop
! switch nsav, nnew so that nnew is the updated aerosol
ntemp=nsav
nsav=nnew
nnew=ntemp
! call t_startf ('nsubmix')
! load new droplets in layers above, below clouds
dtmin=dtmicro
ekk(0)=0.0
ekk(pver)=0.0
do k=1,pver-1
! rce-comment -- ekd(k) is eddy-diffusivity at k/k+1 interface
! want ekk(k) = ekd(k) * (density at k/k+1 interface)
! so use pint(i,k+1) as pint is 1:pverp
! ekk(k)=ekd(k)*2.*pint(i,k)/(rair*(temp(i,k)+temp(i,k+1)))
! ekk(k)=ekd(k)*2.*pint(i,k+1)/(rair*(temp(i,k)+temp(i,k+1)))
ekk(k)=ekd(k)*csbot(k)
enddo
do k=1,pver
km1=max0(k-1,1)
ekkp(k)=zn(k)*ekk(k)*zs(k)
ekkm(k)=zn(k)*ekk(k-1)*zs(km1)
tinv=ekkp(k)+ekkm(k)
if(k.eq.pver)tinv=tinv+surfratemax
! rce-comment -- tinv is the sum of all first-order-loss-rates
! for the layer. for most layers, the activation loss rate
! (for interstitial particles) is accounted for by the loss by
! turb-transfer to the layer above.
! k=pver is special, and the loss rate for activation within
! the layer must be added to tinv. if not, the time step
! can be too big, and explmix can produce negative values.
! the negative values are reset to zero, resulting in an
! artificial source.
if(k.eq.pver)tinv=tinv+taumix_internal_pver_inv
if(tinv.gt.1.e-6)then
dtt=1./tinv
dtmin=min(dtmin,dtt)
endif
enddo
dtmix=0.9*dtmin
nsubmix=dtmicro/dtmix+1
if(nsubmix>100)then
nsubmix_bnd=100
else
nsubmix_bnd=nsubmix
endif
count_submix(nsubmix_bnd)=count_submix(nsubmix_bnd)+1
dtmix=dtmicro/nsubmix
fac_srflx = -1.0/(zn(pver)*nsubmix)
do k=1,pver
kp1=min(k+1,pver)
km1=max(k-1,1)
! maximum overlap assumption
if(cldn(i,kp1).gt.0)then
overlapp(k)=min(cldn(i,k)/cldn(i,kp1),1._r8)
else
overlapp(k)=1.
endif
if(cldn(i,km1).gt.0)then
overlapm(k)=min(cldn(i,k)/cldn(i,km1),1._r8)
else
overlapm(k)=1.
endif
enddo
! rce-comment
! the activation source(k) = mact(k,m)*raercol(kp1,lmass)
! should not exceed the rate of transfer of unactivated particles
! from kp1 to k which = ekkp(k)*raercol(kp1,lmass)
! however it might if things are not "just right" in subr activate
! the following is a safety measure to avoid negatives in explmix
do k = 1, pver-1
do m = 1, ntot_amode
nact(k,m) = min( nact(k,m), ekkp(k) )
mact(k,m) = min( mact(k,m), ekkp(k) )
end do
end do
old_cloud_nsubmix_loop: &
do n=1,nsubmix
qncld(:)=qcld(:)
! switch nsav, nnew so that nsav is the updated aerosol
ntemp=nsav
nsav=nnew
nnew=ntemp
srcn(:)=0.0
do m=1,ntot_amode
lnum=numptr_amode(m)
lnumcw=numptrcw_amode(m)
! update droplet source
! rce-comment- activation source in layer k involves particles from k+1
! srcn(:)=srcn(:)+nact(:,m)*(raercol(:,lnum,nsav))
srcn(1:pver-1)=srcn(1:pver-1)+nact(1:pver-1,m)*(raercol(2:pver,lnum,nsav))
! rce-comment- new formulation for k=pver
! srcn( pver )=srcn( pver )+nact( pver ,m)*(raercol( pver,lnum,nsav))
tmpa = raercol(pver,lnum,nsav)*nact(pver,m) &
+ raercol(pver,lnumcw,nsav)*(nact(pver,m) - taumix_internal_pver_inv)
srcn(pver)=srcn(pver) + max(0.0_r8,tmpa)
enddo
call explmix(qcld,srcn,ekkp,ekkm,overlapp,overlapm, &
qncld,surfrate(ixndrop),zero,pver,dtmix,.false.)
! rce-comment
! the interstitial particle mixratio is different in clear/cloudy portions
! of a layer, and generally higher in the clear portion. (we have/had
! a method for diagnosing the the clear/cloudy mixratios.) the activation
! source terms involve clear air (from below) moving into cloudy air (above).
! in theory, the clear-portion mixratio should be used when calculating
! source terms
do m=1,ntot_amode
lnum=numptr_amode(m)
lnumcw=numptrcw_amode(m)
if(lnum>0)then
! rce-comment - activation source in layer k involves particles from k+1
! source(:)= nact(:,m)*(raercol(:,lnum,nsav))
source(1:pver-1)= nact(1:pver-1,m)*(raercol(2:pver,lnum,nsav))
! rce-comment- new formulation for k=pver
! source( pver )= nact( pver, m)*(raercol( pver,lnum,nsav))
tmpa = raercol(pver,lnum,nsav)*nact(pver,m) &
+ raercol(pver,lnumcw,nsav)*(nact(pver,m) - taumix_internal_pver_inv)
source(pver) = max(0.0_r8,tmpa)
! flxconv=cflx(i,lnum)*zn(pver)
flxconv=0.
call explmix(raercol_cw(1,lnumcw,nnew),source,ekkp,ekkm,overlapp,overlapm, &
raercol_cw(1,lnumcw,nsav),surfrate_cw(lnumcw),zero,pver,dtmix,&
.false.)
call explmix(raercol(1,lnum,nnew),source,ekkp,ekkm,overlapp,overlapm, &
raercol(1,lnum,nsav),surfrate(lnum),flxconv,pver,dtmix, &
.true.,raercol_cw(1,lnumcw,nsav))
endif
do l=1,nspec_amode(m)
lmass=lmassptr_amode(l,m)
lmasscw=lmassptrcw_amode(l,m)
! rce-comment - activation source in layer k involves particles from k+1
! source(:)= mact(:,m)*(raercol(:,lmass,nsav))
source(1:pver-1)= mact(1:pver-1,m)*(raercol(2:pver,lmass,nsav))
! rce-comment- new formulation for k=pver
! source( pver )= mact( pver ,m)*(raercol( pver,lmass,nsav))
tmpa = raercol(pver,lmass,nsav)*mact(pver,m) &
+ raercol(pver,lmasscw,nsav)*(mact(pver,m) - taumix_internal_pver_inv)
source(pver) = max(0.0_r8,tmpa)
! flxconv=cflx(i,lmass)*zn(pver)
flxconv=0.
call explmix(raercol_cw(1,lmasscw,nnew),source,ekkp,ekkm,overlapp,overlapm, &
raercol_cw(1,lmasscw,nsav),surfrate_cw(lmasscw),zero,pver,dtmix,&
.false.)
call explmix(raercol(1,lmass,nnew),source,ekkp,ekkm,overlapp,overlapm, &
raercol(1,lmass,nsav),surfrate(lmass),flxconv, pver,dtmix,&
.true.,raercol_cw(1,lmasscw,nsav))
enddo ! l
! lwater=lwaterptr_amode(m) ! aerosol water
! source(:)=0.
! call explmix( raercol(1,lwater,nnew),source,ekkp,ekkm,overlapp,overlapm, &
! raercol(1,lwater,nsav),surfrate(lwater),zero,pver,dtmix,&
! .true.,source)
enddo ! m
enddo old_cloud_nsubmix_loop
! call t_stopf ('nsubmix')
! evaporate particles again if no cloud
do k=1,pver
if(cldn(i,k).eq.0.)then
! no cloud
qcld(k)=0.
! convert activated aerosol to interstitial in decaying cloud
do m=1,ntot_amode
lnum=numptr_amode(m)
lnumcw=numptrcw_amode(m)
if(lnum.gt.0)then
raercol(k,lnum,nnew)=raercol(k,lnum,nnew)+raercol_cw(k,lnumcw,nnew)
raercol_cw(k,lnumcw,nnew)=0.
endif
!
do l=1,nspec_amode(m)
lmass=lmassptr_amode(l,m)
lmasscw=lmassptrcw_amode(l,m)
! raercol(k,lmass,nnew)=raercol(k,lmass,nnew)+raercol(k,lmasscw,nnew)
! raercol(k,lmasscw,nnew)=0.
raercol(k,lmass,nnew)=raercol(k,lmass,nnew)+raercol_cw(k,lmasscw,nnew)
raercol_cw(k,lmasscw,nnew)=0.
enddo
enddo
endif
enddo
! droplet number
ndropcol(i)=0.
do k=1,pver
ndropmix(i,k)=(qcld(k)-ncldwtr(i,k))*dtinv - nsource(i,k)
tendnd(i,k) = (max(qcld(k),1.e-6_r8)-ncldwtr(i,k))*dtinv
ndropcol(i) = ndropcol(i) + ncldwtr(i,k)*pdel(i,k)
enddo
ndropcol(i) = ndropcol(i)/gravit
! aerosol tendency
do m=1,ntot_amode
lnum=numptr_amode(m)
lnumcw=numptrcw_amode(m)
if(lnum.gt.0)then
qqcwtend(i,:,lnumcw)=(raercol_cw(:,lnumcw,nnew)-qqcw(i,:,lnumcw))*dtinv
qqcw(i,:,lnumcw)=raercol_cw(:,lnumcw,nnew) ! update cloud-borne aerosol
raertend(i,:,lnum)= (raercol(:,lnum,nnew)-raer(i,:,lnum))*dtinv
endif
do l=1,nspec_amode(m)
lmass=lmassptr_amode(l,m)
lmasscw=lmassptrcw_amode(l,m)
qqcwtend(i,:,lmasscw)=(raercol_cw(:,lmasscw,nnew)-qqcw(i,:,lmasscw))*dtinv
qqcw(i,:,lmasscw)=raercol_cw(:,lmasscw,nnew) ! update cloud-borne aerosol
raertend(i,:,lmass)=(raercol(:,lmass,nnew)-raer(i,:,lmass))*dtinv
enddo
! lwater=lwaterptr_amode(m)
! raertend(i,:,lwater)=(raercol(:,lwater,nnew)-raer(i,:,lwater))*dtinv
enddo
enddo overall_main_i_loop
! end of main loop over i/longitude ....................................
call outfld('NDROPCOL', ndropcol , pcols, lchnk )
call outfld('NDROPSRC', nsource , pcols, lchnk )
call outfld('NDROPMIX', ndropmix , pcols, lchnk )
call outfld('LCLOUD ', cldn , pcols, lchnk )
call outfld('WTKE ', wtke , pcols, lchnk )
call ccncalc(lchnk,ncol,temp,cs,raer,qqcw,ccn,psat,supersat,alogsig,npv)
do l=1,psat
call outfld(ccn_name(l), ccn(1,1,l) , pcols, lchnk )
enddo
! do column tendencies
do m = 1, ntot_amode
do lphase = 1, 2
do lspec = 0, nspec_amode(m)+1 ! loop over number + chem constituents + water
if (lspec == 0) then ! number
if (lphase == 1) then
l = numptr_amode(m)
else
l = numptrcw_amode(m)
endif
else if (lspec <= nspec_amode(m)) then ! non-water mass
if (lphase == 1) then
l = lmassptr_amode(lspec,m)
else
l = lmassptrcw_amode(lspec,m)
endif
else ! water mass
! if (lphase == 1) then
! l = lwaterptr_amode(m)
! else
cycle
! end if
end if
if (l <= 0) cycle
! mixnuc1 tendency for interstitial = (total tendency) - cflx
! mixnuc1 tendency for activated = (total tendency)
coltend(:) = 0.0
if (lphase == 1) then
tmpname = cnst_name(l)
do i = 1, ncol
coltend(i) = sum( pdel(i,:)*raertend(i,:,l) )/gravit
! coltend(i) = coltend(i) - cflx(i,l)
end do
else
tmpname = cnst_name_cw(l)
do i = 1, ncol
coltend(i) = sum( pdel(i,:)*qqcwtend(i,:,l) )/gravit
end do
end if
fieldname = trim(tmpname) // '_mixnuc1'
call outfld( fieldname, coltend, pcols, lchnk)
end do ! lspec
end do ! lphase
end do ! m
return
end subroutine dropmixnuc
subroutine explmix( q, src, ekkp, ekkm, overlapp, overlapm, & 5
qold, surfrate, flxconv, pver, dt, is_unact, qactold )
! explicit integration of droplet/aerosol mixing
! with source due to activation/nucleation
implicit none
integer, intent(in) :: pver ! number of levels
real(r8), intent(out) :: q(pver) ! mixing ratio to be updated
real(r8), intent(in) :: qold(pver) ! mixing ratio from previous time step
real(r8), intent(in) :: src(pver) ! source due to activation/nucleation (/s)
real(r8), intent(in) :: ekkp(pver) ! zn*zs*density*diffusivity (kg/m3 m2/s) at interface
! below layer k (k,k+1 interface)
real(r8), intent(in) :: ekkm(pver) ! zn*zs*density*diffusivity (kg/m3 m2/s) at interface
! above layer k (k,k+1 interface)
real(r8), intent(in) :: overlapp(pver) ! cloud overlap below
real(r8), intent(in) :: overlapm(pver) ! cloud overlap above
real(r8), intent(in) :: surfrate ! surface exchange rate (/s)
real(r8), intent(in) :: flxconv ! convergence of flux from surface
real(r8), intent(in) :: dt ! time step (s)
logical, intent(in) :: is_unact ! true if this is an unactivated species
real(r8), intent(in),optional :: qactold(pver)
! mixing ratio of ACTIVATED species from previous step
! *** this should only be present
! if the current species is unactivated number/sfc/mass
integer k,kp1,km1
if ( is_unact ) then
! the qactold*(1-overlap) terms are resuspension of activated material
do k=1,pver
kp1=min(k+1,pver)
km1=max(k-1,1)
q(k) = qold(k) + dt*( - src(k) + ekkp(k)*(qold(kp1) - qold(k) + &
qactold(kp1)*(1.0-overlapp(k))) &
+ ekkm(k)*(qold(km1) - qold(k) + &
qactold(km1)*(1.0-overlapm(k))) )
! force to non-negative
! if(q(k)<-1.e-30)then
! write(iulog,*)'q=',q(k),' in explmix'
q(k)=max(q(k),0._r8)
! endif
end do
! diffusion loss at base of lowest layer
q(pver)=q(pver)-surfrate*qold(pver)*dt+flxconv*dt
! force to non-negative
! if(q(pver)<-1.e-30)then
! write(iulog,*)'q=',q(pver),' in explmix'
q(pver)=max(q(pver),0._r8)
! endif
else
do k=1,pver
kp1=min(k+1,pver)
km1=max(k-1,1)
q(k) = qold(k) + dt*(src(k) + ekkp(k)*(overlapp(k)*qold(kp1)-qold(k)) + &
ekkm(k)*(overlapm(k)*qold(km1)-qold(k)) )
! force to non-negative
! if(q(k)<-1.e-30)then
! write(iulog,*)'q=',q(k),' in explmix'
q(k)=max(q(k),0._r8)
! endif
end do
! diffusion loss at base of lowest layer
q(pver)=q(pver)-surfrate*qold(pver)*dt+flxconv*dt
! force to non-negative
! if(q(pver)<-1.e-30)then
! write(iulog,*)'q=',q(pver),' in explmix'
q(pver)=max(q(pver),0._r8)
! endif
end if
return
end subroutine explmix
subroutine activate_init 1,17
use ppgrid, only: pver
use modal_aero_data
use cam_history, only: addfld, add_default, phys_decomp
use physconst, only: rhoh2o, mwh2o, r_universal
use error_function, only: erf
implicit none
integer l,m
real(r8) arg
integer lnum, lmass
character*16 ccn_longname
! mathematical constants
zero=0._r8
third=1./3._r8
twothird=2.*third
sixth=1./6._r8
sq2=sqrt(2._r8)
pi=4._r8*atan(1.0_r8)
sqpi=sqrt(pi)
t0=273.
surften=0.076_r8
aten=2.*mwh2o*surften/(r_universal*t0*rhoh2o)
alogaten=log(aten)
alog2=log(2._r8)
alog3=log(3._r8)
do m=1,ntot_amode
! use only if width of size distribution is prescribed
alogsig(m)=log(sigmag_amode(m))
exp45logsig(m)=exp(4.5*alogsig(m)*alogsig(m))
argfactor(m)=2./(3.*sqrt(2.)*alogsig(m))
f1(m)=0.5*exp(2.5*alogsig(m)*alogsig(m))
f2(m)=1.+0.25*alogsig(m)
lnum=numptr_amode(m)
do l=1,nspec_amode(m)
lmass=lmassptr_amode(l,m)
if(lmass.le.0)then
write(iulog,*)'lmassptr_amode(',l,m,')=',lmass
call endrun
endif
enddo
npv(m)=6./(pi*dgnum_amode(m)**3*exp45logsig(m))
enddo
call addfld ('WTKE ', 'm/s ', pver, 'A', 'Standard deviation of updraft velocity' ,phys_decomp)
call addfld ('LCLOUD ', ' ', pver, 'A', 'Liquid cloud fraction' ,phys_decomp)
call addfld ('NDROPMIX ', '#/kg/s ', pver, 'A', 'Droplet number mixing' ,phys_decomp)
call addfld ('NDROPSRC ', '#/kg/s ', pver, 'A', 'Droplet number source' ,phys_decomp)
call addfld ('NDROPSNK ', '#/kg/s ', pver, 'A', 'Droplet number loss by microphysics' ,phys_decomp)
call addfld ('NDROPCOL ', '#/m2 ', 1, 'A', 'Column droplet number' ,phys_decomp)
call add_default ('WTKE ', 1, ' ')
call add_default ('LCLOUD ', 1, ' ')
call add_default ('NDROPSRC', 1, ' ')
call add_default ('NDROPMIX', 1, ' ')
call add_default ('NDROPSNK', 1, ' ')
call add_default ('NDROPCOL', 1, ' ')
return
end subroutine activate_init
subroutine activate_modal(wbar, sigw, wdiab, wminf, wmaxf, tair, rhoair, & 2,9
na, pmode, nmode, volume, sigman, hygro, &
fn, fm, fluxn, fluxm, flux_fullact )
! calculates number, surface, and mass fraction of aerosols activated as CCN
! calculates flux of cloud droplets, surface area, and aerosol mass into cloud
! assumes an internal mixture within each of up to pmode multiple aerosol modes
! a gaussiam spectrum of updrafts can be treated.
! mks units
! Abdul-Razzak and Ghan, A parameterization of aerosol activation.
! 2. Multiple aerosol types. J. Geophys. Res., 105, 6837-6844.
use physconst, only: rair, epsilo, cpair, rh2o, latvap, gravit, &
rhoh2o, mwh2o, r_universal
use wv_saturation, only: estblf, epsqs
use error_function, only: erf
implicit none
! input
integer pmode ! dimension of modes
real(r8) wbar ! grid cell mean vertical velocity (m/s)
real(r8) sigw ! subgrid standard deviation of vertical vel (m/s)
real(r8) wdiab ! diabatic vertical velocity (0 if adiabatic)
real(r8) wminf ! minimum updraft velocity for integration (m/s)
real(r8) wmaxf ! maximum updraft velocity for integration (m/s)
real(r8) tair ! air temperature (K)
real(r8) rhoair ! air density (kg/m3)
real(r8) na(pmode) ! aerosol number concentration (/m3)
integer nmode ! number of aerosol modes
real(r8) volume(pmode) ! aerosol volume concentration (m3/m3)
real(r8) sigman(pmode) ! geometric standard deviation of aerosol size distribution
real(r8) hygro(pmode) ! hygroscopicity of aerosol mode
! output
real(r8) fn(pmode) ! number fraction of aerosols activated
real(r8) fm(pmode) ! mass fraction of aerosols activated
real(r8) fluxn(pmode) ! flux of activated aerosol number fraction into cloud (cm/s)
real(r8) fluxm(pmode) ! flux of activated aerosol mass fraction into cloud (cm/s)
real(r8) flux_fullact ! flux of activated aerosol fraction assuming 100% activation (cm/s)
! rce-comment
! used for consistency check -- this should match (ekd(k)*zs(k))
! also, fluxm/flux_fullact gives fraction of aerosol mass flux
! that is activated
! local
integer, parameter:: nx=200
integer iquasisect_option, isectional
real(r8) integ,integf
real(r8), parameter :: p0 = 1013.25e2_r8 ! reference pressure (Pa)
real(r8) xmin(maxd_amode),xmax(maxd_amode) ! ln(r) at section interfaces
real(r8) volmin(maxd_amode),volmax(maxd_amode) ! volume at interfaces
real(r8) tmass ! total aerosol mass concentration (g/cm3)
real(r8) sign(maxd_amode) ! geometric standard deviation of size distribution
real(r8) rm ! number mode radius of aerosol at max supersat (cm)
real(r8) pres ! pressure (Pa)
real(r8) path ! mean free path (m)
real(r8) diff ! diffusivity (m2/s)
real(r8) conduct ! thermal conductivity (Joule/m/sec/deg)
real(r8) diff0,conduct0
real(r8) es ! saturation vapor pressure
real(r8) qs ! water vapor saturation mixing ratio
real(r8) dqsdt ! change in qs with temperature
real(r8) dqsdp ! change in qs with pressure
real(r8) g ! thermodynamic function (m2/s)
real(r8) zeta(maxd_amode), eta(maxd_amode)
real(r8) lnsmax ! ln(smax)
real(r8) alpha
real(r8) gamma
real(r8) beta
real(r8) sqrtg(maxd_amode)
real(r8) :: amcube(maxd_amode) ! cube of dry mode radius (m)
real(r8) :: smcrit(maxd_amode) ! critical supersatuation for activation
real(r8) :: lnsm(maxd_amode) ! ln(smcrit)
real(r8) smc(maxd_amode) ! critical supersaturation for number mode radius
real(r8) sumflx_fullact
real(r8) sumflxn(maxd_amode)
real(r8) sumflxm(maxd_amode)
real(r8) sumfn(maxd_amode)
real(r8) sumfm(maxd_amode)
real(r8) fnold(maxd_amode) ! number fraction activated
real(r8) fmold(maxd_amode) ! mass fraction activated
real(r8) wold,gold
real(r8) alogam
real(r8) rlo,rhi,xint1,xint2,xint3,xint4
real(r8) wmin,wmax,w,dw,dwmax,dwmin,wnuc,dwnew,wb
real(r8) dfmin,dfmax,fnew,fold,fnmin,fnbar,fsbar,fmbar
real(r8) alw,sqrtalw
real(r8) smax
real(r8) x,arg
real(r8) xmincoeff,xcut,volcut,surfcut
real(r8) z,z1,z2,wf1,wf2,zf1,zf2,gf1,gf2,gf
real(r8) etafactor1,etafactor2(maxd_amode),etafactor2max
integer m,n
! numerical integration parameters
real(r8), parameter :: eps=0.3_r8,fmax=0.99_r8,sds=3._r8
real(r8), parameter :: namin=1.e6_r8 ! minimum aerosol number concentration (/m3)
integer ndist(nx) ! accumulates frequency distribution of integration bins required
data ndist/nx*0/
save ndist
if(maxd_amode<pmode)then
write(iulog,*)'maxd_amode,pmode in activate =',maxd_amode,pmode
call endrun('activate')
endif
fn(:)=0._r8
fm(:)=0._r8
fluxn(:)=0._r8
fluxm(:)=0._r8
flux_fullact=0._r8
if(nmode.eq.1.and.na(1).lt.1.e-20_r8)return
if(sigw.le.1.e-5_r8.and.wbar.le.0.)return
pres=rair*rhoair*tair
diff0=0.211e-4_r8*(p0/pres)*(tair/t0)**1.94
conduct0=(5.69_r8+0.017_r8*(tair-t0))*4.186e2_r8*1.e-5_r8 ! convert to J/m/s/deg
es = estblf(tair)
qs = epsilo*es/(pres-(1.0_r8 - epsqs)*es)
dqsdt=latvap/(rh2o*tair*tair)*qs
alpha=gravit*(latvap/(cpair*rh2o*tair*tair)-1./(rair*tair))
gamma=(1+latvap/cpair*dqsdt)/(rhoair*qs)
etafactor2max=1.e10/(alpha*wmaxf)**1.5 ! this should make eta big if na is very small.
do m=1,nmode
if(volume(m).gt.1.e-39_r8.and.na(m).gt.1.e-39_r8)then
! number mode radius (m)
! write(iulog,*)'alogsig,volc,na=',alogsig(m),volc(m),na(m)
amcube(m)=(3.*volume(m)/(4.*pi*exp45logsig(m)*na(m))) ! only if variable size dist
! growth coefficent Abdul-Razzak & Ghan 1998 eqn 16
! should depend on mean radius of mode to account for gas kinetic effects
! see Fountoukis and Nenes, JGR2005 and Meskhidze et al., JGR2006
! for approriate size to use for effective diffusivity.
g=1._r8/(rhoh2o/(diff0*rhoair*qs) &
+latvap*rhoh2o/(conduct0*tair)*(latvap/(rh2o*tair)-1._r8))
sqrtg(m)=sqrt(g)
beta=2._r8*pi*rhoh2o*g*gamma
etafactor2(m)=1._r8/(na(m)*beta*sqrtg(m))
if(hygro(m).gt.1.e-10)then
smc(m)=2.*aten*sqrt(aten/(27.*hygro(m)*amcube(m))) ! only if variable size dist
else
smc(m)=100.
endif
! write(iulog,*)'sm,hygro,amcube=',smcrit(m),hygro(m),amcube(m)
else
g=1._r8/(rhoh2o/(diff0*rhoair*qs) &
+latvap*rhoh2o/(conduct0*tair)*(latvap/(rh2o*tair)-1._r8))
sqrtg(m)=sqrt(g)
smc(m)=1._r8
etafactor2(m)=etafactor2max ! this should make eta big if na is very small.
endif
lnsm(m)=log(smc(m)) ! only if variable size dist
! write(iulog,'(a,i4,4g12.2)')'m,na,amcube,hygro,sm,lnsm=', &
! m,na(m),amcube(m),hygro(m),sm(m),lnsm(m)
enddo
if(sigw.gt.1.e-5_r8)then ! spectrum of updrafts
wmax=min(wmaxf,wbar+sds*sigw)
wmin=max(wminf,-wdiab)
wmin=max(wmin,wbar-sds*sigw)
w=wmin
dwmax=eps*sigw
dw=dwmax
dfmax=0.2_r8
dfmin=0.1_r8
if(wmax.le.w)then
do m=1,nmode
fluxn(m)=0.
fn(m)=0.
fluxm(m)=0.
fm(m)=0.
enddo
flux_fullact=0._r8
return
endif
do m=1,nmode
sumflxn(m)=0._r8
sumfn(m)=0._r8
fnold(m)=0._r8
sumflxm(m)=0._r8
sumfm(m)=0._r8
fmold(m)=0._r8
enddo
sumflx_fullact=0._r8
fold=0._r8
wold=0._r8
gold=0._r8
dwmin = min( dwmax, 0.01_r8 )
do n=1,200
100 wnuc=w+wdiab
! write(iulog,*)'wnuc=',wnuc
alw=alpha*wnuc
sqrtalw=sqrt(alw)
etafactor1=alw*sqrtalw
do m=1,nmode
eta(m)=etafactor1*etafactor2(m)
zeta(m)=twothird*sqrtalw*aten/sqrtg(m)
enddo
call maxsat(zeta,eta,nmode,smc,smax)
! write(iulog,*)'w,smax=',w,smax
lnsmax=log(smax)
x=twothird*(lnsm(nmode)-lnsmax)/(sq2*alogsig(nmode))
fnew=0.5_r8*(1._r8-erf(x))
dwnew = dw
if(fnew-fold.gt.dfmax.and.n.gt.1)then
! reduce updraft increment for greater accuracy in integration
if (dw .gt. 1.01*dwmin) then
dw=0.7_r8*dw
dw=max(dw,dwmin)
w=wold+dw
go to 100
else
dwnew = dwmin
endif
endif
if(fnew-fold.lt.dfmin)then
! increase updraft increment to accelerate integration
dwnew=min(1.5_r8*dw,dwmax)
endif
fold=fnew
z=(w-wbar)/(sigw*sq2)
g=exp(-z*z)
fnmin=1._r8
xmincoeff=alogaten-twothird*(lnsmax-alog2)-alog3
do m=1,nmode
! modal
x=twothird*(lnsm(m)-lnsmax)/(sq2*alogsig(m))
fn(m)=0.5_r8*(1.-erf(x))
fnmin=min(fn(m),fnmin)
! integration is second order accurate
! assumes linear variation of f*g with w
fnbar=(fn(m)*g+fnold(m)*gold)
arg=x-1.5_r8*sq2*alogsig(m)
fm(m)=0.5_r8*(1._r8-erf(arg))
fmbar=(fm(m)*g+fmold(m)*gold)
wb=(w+wold)
if(w.gt.0.)then
sumflxn(m)=sumflxn(m)+sixth*(wb*fnbar &
+(fn(m)*g*w+fnold(m)*gold*wold))*dw
sumflxm(m)=sumflxm(m)+sixth*(wb*fmbar &
+(fm(m)*g*w+fmold(m)*gold*wold))*dw
endif
sumfn(m)=sumfn(m)+0.5_r8*fnbar*dw
! write(iulog,'(a,9g10.2)')'lnsmax,lnsm(m),x,fn(m),fnold(m),g,gold,fnbar,dw=',lnsmax,lnsm(m),x,fn(m),fnold(m),g,gold,fnbar,dw
fnold(m)=fn(m)
sumfm(m)=sumfm(m)+0.5_r8*fmbar*dw
fmold(m)=fm(m)
enddo
! same form as sumflxm but replace the fm with 1.0
sumflx_fullact = sumflx_fullact &
+ sixth*(wb*(g+gold) + (g*w+gold*wold))*dw
! sumg=sumg+0.5_r8*(g+gold)*dw
gold=g
wold=w
dw=dwnew
if(n.gt.1.and.(w.gt.wmax.or.fnmin.gt.fmax))go to 20
w=w+dw
enddo
write(iulog,*)'do loop is too short in activate'
write(iulog,*)'wmin=',wmin,' w=',w,' wmax=',wmax,' dw=',dw
write(iulog,*)'wbar=',wbar,' sigw=',sigw,' wdiab=',wdiab
write(iulog,*)'wnuc=',wnuc
write(iulog,*)'na=',(na(m),m=1,nmode)
write(iulog,*)'fn=',(fn(m),m=1,nmode)
! dump all subr parameters to allow testing with standalone code
! (build a driver that will read input and call activate)
write(iulog,*)'wbar,sigw,wdiab,tair,rhoair,nmode='
write(iulog,*) wbar,sigw,wdiab,tair,rhoair,nmode
write(iulog,*)'na=',na
write(iulog,*)'volume=', (volume(m),m=1,nmode)
write(iulog,*)'sigman=',sigman
write(iulog,*)'hydro='
write(iulog,*) hygro
call endrun
20 continue
ndist(n)=ndist(n)+1
if(w.lt.wmaxf)then
! contribution from all updrafts stronger than wmax
! assuming constant f (close to fmax)
wnuc=w+wdiab
z1=(w-wbar)/(sigw*sq2)
z2=(wmaxf-wbar)/(sigw*sq2)
g=exp(-z1*z1)
integ=sigw*0.5*sq2*sqpi*(erf(z2)-erf(z1))
! consider only upward flow into cloud base when estimating flux
wf1=max(w,zero)
zf1=(wf1-wbar)/(sigw*sq2)
gf1=exp(-zf1*zf1)
wf2=max(wmaxf,zero)
zf2=(wf2-wbar)/(sigw*sq2)
gf2=exp(-zf2*zf2)
gf=(gf1-gf2)
integf=wbar*sigw*0.5*sq2*sqpi*(erf(zf2)-erf(zf1))+sigw*sigw*gf
do m=1,nmode
sumflxn(m)=sumflxn(m)+integf*fn(m)
sumfn(m)=sumfn(m)+fn(m)*integ
sumflxm(m)=sumflxm(m)+integf*fm(m)
sumfm(m)=sumfm(m)+fm(m)*integ
enddo
! same form as sumflxm but replace the fm with 1.0
sumflx_fullact = sumflx_fullact + integf
! sumg=sumg+integ
endif
do m=1,nmode
fn(m)=sumfn(m)/(sq2*sqpi*sigw)
! fn(m)=sumfn(m)/(sumg)
if(fn(m).gt.1.01)then
write(iulog,*)'fn=',fn(m),' > 1 in activate'
write(iulog,*)'w,m,na,amcube=',w,m,na(m),amcube(m)
write(iulog,*)'integ,sumfn,sigw=',integ,sumfn(m),sigw
call endrun('activate')
endif
fluxn(m)=sumflxn(m)/(sq2*sqpi*sigw)
fm(m)=sumfm(m)/(sq2*sqpi*sigw)
! fm(m)=sumfm(m)/(sumg)
if(fm(m).gt.1.01)then
write(iulog,*)'fm=',fm(m),' > 1 in activate'
endif
fluxm(m)=sumflxm(m)/(sq2*sqpi*sigw)
enddo
! same form as fluxm
flux_fullact = sumflx_fullact/(sq2*sqpi*sigw)
else
! single updraft
wnuc=wbar+wdiab
if(wnuc.gt.0._r8)then
w=wbar
alw=alpha*wnuc
sqrtalw=sqrt(alw)
etafactor1=alw*sqrtalw
do m=1,nmode
eta(m)=etafactor1*etafactor2(m)
zeta(m)=twothird*sqrtalw*aten/sqrtg(m)
enddo
call maxsat(zeta,eta,nmode,smc,smax)
lnsmax=log(smax)
xmincoeff=alogaten-twothird*(lnsmax-alog2)-alog3
do m=1,nmode
! modal
x=twothird*(lnsm(m)-lnsmax)/(sq2*alogsig(m))
fn(m)=0.5_r8*(1._r8-erf(x))
arg=x-1.5_r8*sq2*alogsig(m)
fm(m)=0.5_r8*(1._r8-erf(arg))
if(wbar.gt.0._r8)then
fluxn(m)=fn(m)*w
fluxm(m)=fm(m)*w
endif
enddo
flux_fullact = w
endif
endif
return
end subroutine activate_modal
subroutine maxsat(zeta,eta,nmode,smc,smax) 3
! calculates maximum supersaturation for multiple
! competing aerosol modes.
! Abdul-Razzak and Ghan, A parameterization of aerosol activation.
! 2. Multiple aerosol types. J. Geophys. Res., 105, 6837-6844.
implicit none
integer nmode ! number of modes
real(r8) smc(maxd_amode) ! critical supersaturation for number mode radius
real(r8) zeta(maxd_amode), eta(maxd_amode)
real(r8) smax ! maximum supersaturation
integer m ! mode index
real(r8) sum, g1, g2, g1sqrt, g2sqrt
do m=1,nmode
if(zeta(m).gt.1.e5_r8*eta(m).or.smc(m)*smc(m).gt.1.e5_r8*eta(m))then
! weak forcing. essentially none activated
smax=1.e-20_r8
else
! significant activation of this mode. calc activation all modes.
go to 1
endif
enddo
return
1 continue
sum=0
do m=1,nmode
if(eta(m).gt.1.e-20_r8)then
g1=zeta(m)/eta(m)
g1sqrt=sqrt(g1)
g1=g1sqrt*g1
g1=g1sqrt*g1
g2=smc(m)/sqrt(eta(m)+3._r8*zeta(m))
g2sqrt=sqrt(g2)
g2=g2sqrt*g2
sum=sum+(f1(m)*g1+f2(m)*g2)/(smc(m)*smc(m))
else
sum=1.e20_r8
endif
enddo
smax=1._r8/sqrt(sum)
return
end subroutine maxsat
subroutine ccncalc(lchnk,ncol,tair,cs,raer,qqcw,ccn,psat,supersat,alogsig,npv) 1,6
! calculates number concentration of aerosols activated as CCN at
! supersaturation supersat.
! assumes an internal mixture of a multiple externally-mixed aerosol modes
! cgs units
! Ghan et al., Atmos. Res., 1993, 198-221.
use shr_kind_mod, only: r8 => shr_kind_r8
use ppgrid, only: pcols, pver, pverp
use constituents, only: pcnst
use physconst, only: rhoh2o, mwh2o, r_universal
use modal_aero_data
use error_function, only: erf
implicit none
! input
integer lchnk ! chunk index
integer, intent(in) :: ncol ! number of columns
integer, intent(in) :: psat ! number of supesaturations
real(r8), intent(in) :: raer(pcols,pver,pcnst) ! aerosol mass, number mixing ratios
real(r8), intent(in) :: qqcw(pcols,pver,pcnst) ! cloud-borne aerosol mass, number mixing ratios
real(r8), intent(in) :: tair(pcols,pver) ! air temperature (K)
real(r8), intent(in) :: cs(pcols,pver) ! air density (kg/m3)
real(r8), intent(in) :: supersat(psat)
real(r8), intent(in) :: npv(maxd_amode) ! number per volume concentration
real(r8), intent(in) :: alogsig(maxd_amode) ! natl log of geometric standard dev of aerosol
real(r8), intent(out) :: ccn(pcols,pver,psat) ! number conc of aerosols activated at supersat (#/m3)
! local
real(r8) naerosol(pcols) ! interstit+activated aerosol number conc (/m3)
real(r8) vaerosol(pcols) ! interstit+activated aerosol volume conc (m3/m3)
real(r8) exp45logsig(maxd_amode) ! exp(4.5*alogsig**2)
real(r8) amcube(pcols)
real(r8) super(psat) ! supersaturation
real(r8) amcubecoef(maxd_amode)
real(r8) :: surften ! surface tension of water w/respect to air (N/m)
real(r8) surften_coef
real(r8) a(pcols) ! surface tension parameter
real(r8) hygro(pcols) ! aerosol hygroscopicity
real(r8) sm(pcols) ! critical supersaturation at mode radius
real(r8) argfactor(maxd_amode),arg(pcols)
! mathematical constants
real(r8) pi
real(r8) twothird,sq2
integer l,m,n,i,k
real(r8) log,cc
real(r8) smcoefcoef,smcoef(pcols)
integer phase ! phase of aerosol
super(:)=supersat(:)*0.01
pi = 4.*atan(1.0)
sq2=sqrt(2._r8)
twothird=2._r8/3._r8
surften=0.076_r8
surften_coef=2._r8*mwh2o*surften/(r_universal*rhoh2o)
smcoefcoef=2._r8/sqrt(27._r8)
do m=1,ntot_amode
exp45logsig(m)=exp(4.5*alogsig(m)*alogsig(m))
amcubecoef(m)=3./(4.*pi*exp45logsig(m))
argfactor(m)=twothird/(sq2*alogsig(m))
end do
do k=1,pver
do i=1,ncol
ccn(i,k,:)=0.
a(i)=surften_coef/tair(i,k)
smcoef(i)=smcoefcoef*a(i)*sqrt(a(i))
end do
do m=1,ntot_amode
phase=3 ! interstitial+cloudborne
call loadaer(raer,qqcw,1,ncol,k,m,cs,npv(m), phase, &
naerosol, vaerosol, hygro )
where(naerosol(:ncol)>1.e-3)
amcube(:ncol)=amcubecoef(m)*vaerosol(:ncol)/naerosol(:ncol)
sm(:ncol)=smcoef(:ncol)/sqrt(hygro(:ncol)*amcube(:ncol)) ! critical supersaturation
elsewhere
sm(:ncol)=1. ! value shouldn't matter much since naerosol is small
endwhere
do l=1,psat
do i=1,ncol
arg(i)=argfactor(m)*log(sm(i)/super(l))
ccn(i,k,l)=ccn(i,k,l)+naerosol(i)*0.5*(1._r8-erf(arg(i)))
enddo
enddo
enddo
enddo
ccn(:ncol,:,:)=ccn(:ncol,:,:)*1.e-6 ! convert from #/m3 to #/cm3
return
end subroutine ccncalc
subroutine loadaer(raer,qqcw,istart,istop,k,m,cs,npv1, phase, & 3,4
naerosol, vaerosol, hygro )
use shr_kind_mod, only: r8 => shr_kind_r8
use abortutils, only: endrun
use ppgrid, only: pcols, pver
use modal_aero_data
implicit none
! load aerosol number, volume concentrations, and bulk hygroscopicity
real(r8), intent(in) :: raer(pcols,pver,pcnst) ! aerosol mass, number mixing ratios
real(r8), intent(in) :: qqcw(pcols,pver,pcnst) ! cloud-borne aerosol mass, number mixing ratios
integer, intent(in) :: istart, istop ! start, stop column index
integer, intent(in) :: m ! m=mode index
integer, intent(in) :: k ! level index
real(r8), intent(in) :: cs(pcols,pver) ! air density (kg/m3)
integer, intent(in) :: phase ! phase of aerosol: 1 for interstitial, 2 for cloud-borne, 3 for sum
real(r8), intent(in) :: npv1 ! number per volume
real(r8), intent(out) :: naerosol(pcols) ! interstitial number conc (/m3)
real(r8), intent(out) :: vaerosol(pcols) ! interstitial+activated volume conc (m3/m3)
real(r8), intent(out) :: hygro(pcols) ! bulk hygroscopicity of mode
! internal
real(r8) vol(pcols) ! aerosol volume mixing ratio
integer i,lnum,lnumcw,l,ltype,lmass,lmasscw
do i=istart,istop
vaerosol(i)=0._r8
hygro(i)=0._r8
end do
do l=1,nspec_amode(m)
lmass=lmassptr_amode(l,m) ! interstitial
lmasscw=lmassptrcw_amode(l,m) ! cloud-borne
ltype=lspectype_amode(l,m)
if(phase.eq.3)then
do i=istart,istop
! vol(i)=max(raer(i,k,lmass)+raer(i,k,lmasscw),0._r8)/specdens_amode(ltype)
vol(i)=max(raer(i,k,lmass)+qqcw(i,k,lmasscw),0._r8)/specdens_amode(ltype)
end do
elseif(phase.eq.2)then
do i=istart,istop
vol(i)=max(qqcw(i,k,lmasscw),0._r8)/specdens_amode(ltype)
end do
elseif(phase.eq.1)then
do i=istart,istop
vol(i)=max(raer(i,k,lmass),0._r8)/specdens_amode(ltype)
end do
else
write(iulog,*)'phase=',phase,' in loadaer'
call endrun('phase error in loadaer')
endif
do i=istart,istop
vaerosol(i)=vaerosol(i)+vol(i)
hygro(i)=hygro(i)+vol(i)*spechygro(ltype)
end do
enddo
do i=istart,istop
if (vaerosol(i) > 1.0e-30_r8) then ! +++xl add 8/2/2007
hygro(i)=hygro(i)/(vaerosol(i))
vaerosol(i)=vaerosol(i)*cs(i,k)
else
hygro(i)=0.0_r8
vaerosol(i)=0.0_r8
endif
end do
lnum=numptr_amode(m)
lnumcw=numptrcw_amode(m)
if(lnum.gt.0)then
! aerosol number predicted
if(phase.eq.3)then
do i=istart,istop
! naerosol(i)=(raer(i,k,lnum)+raer(i,k,lnumcw))*cs(i,k)
naerosol(i)=(raer(i,k,lnum)+qqcw(i,k,lnumcw))*cs(i,k)
end do
elseif(phase.eq.2)then
do i=istart,istop
naerosol(i)=qqcw(i,k,lnumcw)*cs(i,k)
end do
else
do i=istart,istop
naerosol(i)=raer(i,k,lnum)*cs(i,k)
end do
endif
! adjust number so that dgnumlo < dgnum < dgnumhi
do i=istart,istop
naerosol(i) = max( naerosol(i), vaerosol(i)*voltonumbhi_amode(m) )
naerosol(i) = min( naerosol(i), vaerosol(i)*voltonumblo_amode(m) )
end do
else
! aerosol number diagnosed from mass and prescribed size
do i=istart,istop
naerosol(i)=vaerosol(i)*npv1
naerosol(i)=max(naerosol(i),0._r8)
end do
endif
return
end subroutine loadaer
#endif
end module ndrop