c========================================================================== c 12/17/02: tmg c - switched order so equilibrium is first again c 10/22/02: tmg c - added parameter to subroutine step because not all sections c are always used c 01/25/02: tmg c - use dry and wet basis for section diameter as appropriate c 01/27/01: bkoo c - change the order of equilibrium and dynamic calculations c - change the iteration number c 01/08/01: bkoo c - HYBR-EQUI shift c 10/11/00: bkoo c - establish electroneutrality in the first dynamic section c 09/20/00: bkoo c - combine eqparth & eqpart c - change the order of step & diameter at the end of heqdyn c 06/08/00: bkoo c - add coagulation & nucleation c MAY 2000: bkoo c - correct negative gas conc. after EQPARTh c - EQPARTh : EQPART for HYBRID c cgy supress error messages to units 39 and 90 (6/12/00) c c=========================================================================+ c-------------------------------------------------------------------------- c HEQDYN - HYBRID EQUILIBRIUM AND DYNAMIC GAS/AERSOL PARTITIONING C ROUTINE FOR USE OF BOTH EQUILIBRIUM AND DYNAMIC GAS/AERSOL PARTITIONING C C THE FIRST 'NEQSEC' SECTIONS USE EQUILIBRIUM PARTITIONING AND THE C REMAINING AEROSOL SECTIONS USE DYNAMIC PARTITIONING C C CONSIDERATIONS FOR WHERE TO CUT THE SECTIONS C 1. THE EQUILIBRIUM SECTIONS SHOULD BE SMALL ENOUGH TO REACH EQUILIBRIUM C IN A SINGLE TIME STEP. DP OF ABOUT 1 UM IS PROBABLY SAFE FOR A 10 MIN C TIME STEP (SEE WEXLER AND SEINFELD, AE-V24, 1990). C 2. ONE SHOULD ALSO CONSIDER THE CHEMICAL HOMOGENEITY OF THE SECTIONS TO C BE PLACED IN EQUILIBRIUM. IF THE INITIAL CHEMICAL COMPOSITIONS OF C THE EQUILIBRIUM SIZE SECTIONS ARE TOO DIFFERENT ERRORS WILL BE INTRODUCED C SUBROUTINE HEQDYN(T0,T1,Q) include 'dynamic.inc' c master conc. array: aerosols (ug/m3) then gas (ppm) real*8 q(ntotal) c conc. array for equilibrium: aerosols (ug/m3) then gas (ppm) real*8 qe(ntotale) c conc. array for dynamics: aerosols (ug/m3) then gas (ppm) real*8 qd(ntotald) common / counters / icount, ifail ! bkoo (01/08/01) c gas phase "source" due to equilibrium with small sections (ppm/s) c real*8 sdgas(ngas),sgasrate(ngas),residual(ngas) c common / eqgas / sdgas,sgasrate,residual,t0c c Equation for transition regime (Dahneke 1983) f(yKnD, alfa)=(1.0d0+yKnD)/(1.0d0+2.0d0*yKnD*(1.0d0+yKnD)/alfa) c c STEP 1/3: CALCULATE NUCLEATION RATE FOR THE WHOLE STEP c call nucl(q) c c STEP 2/3: CALCULATE COAGULATION RATE FOR THE WHOLE STEP c call coagul(q) c call step(nsec,q) ! tmg (10/22/02) c c t0c=t0 c set useful indices ng =nsec*nsp nge=neqsec*nsp ngd=(nsec-neqsec)*nsp c c We iterate to minimize bias towards large sections caused by emmissions c only occuring in the dynamic portion of the code c beware however, that too many iterations will bias towards the small c sections because instantaneous mass transfer will be a bad assumption c iterations=2 tstart=t0 tend=tstart+(t1-t0)/float(iterations) do iterate=1,iterations c c transfer concentrations to equilibrium array do isec=1,neqsec indx=(isec-1)*nsp do isp=1,nsp qe(indx+isp)=q(indx+isp) enddo enddo do isp=1,ngas qe(nge+isp)=q(ng+isp) ! bkoo (01/27/01) enddo c call equilibrium partitioning code call eqpart(tstart,qe) c c transfer concentrations to dynamics array do isec=neqsec+1,nsec indx=(isec-1)*nsp indxd=(isec-neqsec-1)*nsp do isp=1,nsp qd(indxd+isp)=q(indx+isp) enddo enddo c this step dictates sequential changes to the gas phase - c first equilibrium then dynamic sections changes c Now, first dynamic then equilibrium: bkoo (01/27/01) c Switched back to equilibrium first: tmg (12/17/02) do isp=1,ngas qd(ngd+isp)=qe(nge+isp) ! bkoo (01/27/01) enddo c call dynamics partitioning code call madm(tstart,tend,qd) c if HYBR failed return (bkoo, 01/08/01) if(ifail.eq.1) return c set the gas "source" term parameters to represent the change in the gas c phase concentrations due to exchange with the equilibrium sections. c We assume that the change in the gas phase conc. decays exponentially c with a rate equal to the sum of the sectional rates - then NH4 is set c equal the rate of negative ions exchanging with the gas phase. The c residual delta Cgas after a time step of exponential decay will be c linearized and added to the exponential "source" to preserve mass balance. c The decay rate given here is a maximum rate. If a sepecies completely c evaporates in a section within the time step, the decay rate will be c slower. This is neglected here for simplicity. c c do k=1,ngas c sgasrate(k)=0.0d0 c sdgas(k)=qe(nge+k)-q(ng+k) c do i=1,neqsec c indx=(i-1)*nsp c if(sdgas(k).gt.0.0d0.and. ! evaporation but c & ((k.eq.ihcl.and.q(indx+kcl).lt.tinys).or. ! species can't evap. c & (k.eq.ihno3.and.q(indx+kno3).lt.tinys))) goto 100 c dp = 1.0d-6 * dsec(i) ! (moving sections) c KnD=2.0d0*lamda(k)/dp c sgasrate(k)=sgasrate(k)+2.0d0*pi*dp*diffus(k)* c & f(KnD, delta(k))*qn(i)*1.0d6 ! in sec-1 c 100 enddo c enddo c sgasrate(inh3)=2*sgasrate(ih2so4)+sgasrate(ihno3)+sgasrate(ihcl) c write(50,*)(sgasrate(k),k=1,ngas) c do k=1,ngas c if(sgasrate(k)*(t1-t0).gt.20) then c residual(k)=0.0d0 ! exp(-20)=~1e-9 c else c residual(k)=sdgas(k)*exp(-sgasrate(k)*(t1-t0))/(t1-t0) c endif c enddo c c EQPARTh can give negative gas concentrations ! done in negchk tmg ! (01/18/02) c do isp=1,ngas c if(qe(nge+isp).lt.0.0d0) then cgy write(90,*)'NEG_GAS: gas(',isp,')=',qe(nge+isp),' t=',tstart c qe(nge+isp)=0.0d0 c endif c enddo c c return new concentrations to the original array do isec=1,neqsec ! from equilibrium sections indx=(isec-1)*nsp do isp=1,nsp q(indx+isp)=qe(indx+isp) enddo enddo do isec=neqsec+1,nsec ! from dynamic sections indx=(isec-1)*nsp indxd=(isec-neqsec-1)*nsp do isp=1,nsp q(indx+isp)=qd(indxd+isp) enddo enddo do isp=1,ngas q(ng+isp)=qd(ngd+isp) ! bkoo (01/27/01) enddo call negchk(tend,q,nsec) c calculate water c call step(nsec,q) ! tmg (10/22/02) c calculate new diameters call ddiameter(q) tstart=tend tend=tend+(t1-t0)/float(iterations) enddo return end