c========================================================================= c 03/09/03: bkoo c - added call to orgps2 to calculate PS for organics c 03/07/03: bkoo c - removed redundant qtemp since step has now parameter for nsec c - removed redundant IF statement c - commented out icount-initializing code c (note that HYBR calls MADM multiple times) c 11/20/02: tmg c - modified timestep check to include factor of qt (and changed c error tolerances c 10/22/02: tmg c - added parameter to subroutine step because not all sections c are always in use c 07/02/01: tmg c - changed to trajectory grid method for solving dynamic equations c 01/27/01: bkoo c - change the STEP position c - water change in HYBR c 01/26/01: bkoo c - c 01/07/01: bkoo c - try a new scheme to estimate xx c 10/08/00: bkoo c - move declaration statements before inline function declaration c 09/19/00: bkoo c - combine diffundm & diffundh, dryinm & dryinh c MAY 2000: bkoo c - combine MADM with HYBRID c - diffundm & dryinm: diffund & dryin for MADM c - diffundh & dryinh: diffund & dryin for HYBRID cgy write error messages only to unit 6, only before a stop c c========================================================================+ subroutine diffund(neq, t, q, tout) c c This subroutine calculates the derivatives of the concetrations of all c the aerosol species, for all sections, as well as the concentrations of c the gases for the new generation c multicomponent aerosol dynamics model (MADM) c ************************************************************************* c WRITTEN BY DR. CHRISTODOULOS PILINIS c January-February 1998 c ************************************************************************* c include 'dynamic.inc' real*8 dq(ntotal) real*8 q(ntotal), q2(ntotal), qlast(ntotal) real*8 dqdt(ntotal) ! Aerosol module master array for c concentration derivatives of aerosol species in ug/m3/sec and c gases in ppm/sec. double precision KnD dimension dqsum(nsp) ! dqdt mass balance for each species dimension pgas(ngas) ! partial pressures of condensible gases (Pa) c add ht array back in (tmg, 07/02/01) dimension ht(nsec) ! total rate of change real*8 merr(nsec), htavg(nsec) c parameters for gas "source" term by the equilibrium sections c real*8 sdgas(ngas),sgasrate(ngas),residual(ngas) integer icase(nsec) c new variables (tmg, 07/02/01) dimension qtlast(nsec), hplast(nsec), hptlast(nsec) dimension htlast(nsec), hilast(nsec,ngas) real*8 masserr, errtol, ltol common / counters / icount, ifail ! bkoo (01/08/01) c common / eqgas / sdgas,sgasrate,residual,t0c common /debug/ icnt common /drycase/ icase c Equation for transition regime (Dahneke 1983) f(yKnD, alfa)=(1.0d0+yKnD)/(1.0d0+2.0d0*yKnD*(1.0d0+yKnD)/alfa) masserr = 0.0d0 tmincr=10.0d0 cbk icount=0 ! bkoo (03/07/03) tpres=t errtol = 1.0d-3 ltol = 1.0d-6 c tinys=1.0d-9 ! minimum mass in ug/m3 tiny2= tinys * govra ! 02/14/02 ng = neq - ngas 20 do i=1,neq dqdt(i)=0.0d0 ! Initialize derivatives enddo c if (icount.gt.200005) return ! tmg (12/17/02) do i=1,nsecx qt(i)=0.0d0 ! Initialize total masses ht(i)=0.0d0 ! Initialize total rates enddo c STEP subroutine updates c, ps, dry, q(H2O), keeping dqdt intact cbk if(aerm.eq.'MADM') call step(nsecx,q) ! (10/22/02) cbk if(aerm.eq.'HYBR') then cbk do i=1,ntotal cbk if(i.le.ng) then cbk qtemp(i)=q(i) cbk else cbk qtemp(i)=0.0d0 cbk endif cbk enddo call step(nsecx,q) ! bkoo (03/07/03) cbk call step(nsecx,qtemp) ! this will "front end" common variables c ps(i,j) and c(i,j). So call STEP again before c using these variables outside of the dynamics c integrator (for equilibrium code, or printing, etc.) cbk do i=1,nsecx cbk q((i-1)*nsp+KH2O)=qtemp((i-1)*nsp+KH2O) ! bkoo (01/27/01) cbk enddo cbk endif cbk calculate the saturation concentrations for organics - bkoo (04/18/02) (03/09/03) call orgps2(q) c If an aerosol oscillates between wet and dry states the numerical c integration can hardly converge because the PS calculations for c the two states are much different from each other. In that case, c we consider the aerosol as metastable liquid state to avoid the c oscillation. c - bkoo (03/05/02) if(isfirst) then isfirst=.FALSE. do i=1,nsecx icdry(i)=0 dold(i)=dry(i) enddo else do i=1,nsecx if(ims(i).eq.0) then ! if not metastable if(dold(i) .and. .not.dry(i)) then if(icdry(i).eq.100) then ims(i)=1 c call errprt('METASTABLE ',t,0,i,1.d0,0) endif icdry(i)=icdry(i)+1 endif dold(i)=dry(i) endif enddo endif do i=1,nsecx ! dbg do k=1,ngas ! dbg if(ps(i,k).lt.0. .and. ps(i,k).ge.-1.e-15) ps(i,k)=0. ! dbg if(ps(i,k).lt.0.) then ! dbg c call errprt('NEG-PS ',t,k,i,ps(i,k),1) ! dbg c stop ! dbg ps(i,k)=0. ! dbg endif ! dbg enddo ! dbg do k=1,ninti ! dbg if(c(i,k).lt.0.d0) c(i,k)=0.d0 ! bkoo (02/14/02) dbg enddo ! dbg enddo ! dbg c do i=1,nsecp1 c do k=1,nsp c intr(i,k)=0.0d0 c end do c end do do i=1,nsecx ! Calculate total mass in each section do k=1,nsp qt(i)=qt(i)+q((i-1)*nsp+k) enddo if(qn(i+nsecx2).gt.0.and.qt(i).gt.tinys) then !avoid empty sections dsec(i+nsecx2)=(qt(i)/qn(i+nsecx2))**(0.33333) cgy if(dsec(i+nsecx2).lt.0.005) write(90,*) 'at t= ',t, cgy & ' Diameter of section ',i+neqsec, cgy & ' is smaller than 0.005 um at Dp = ',dsec(i+neqsec) cgy if(dsec(i+nsecx2).gt.20.0) write(90,*) 'at t= ',t, cgy & ' Diameter of section ',i+neqsec, cgy & ' is larger than 20.0 um at Dp = ',dsec(i+neqsec) endif enddo c calculate the partial pressures of the consensible gases c Note that since for gases, q is in ppm, to convert to Pa c we multiply by (1.0d-6 * 1.015d5) c we multiply by (1.01325d-1 * pres) instead of 1.015d-1 (bkoo, 06/09/00) c pgas(INH3) =q(ng+INH3) * 1.01325d-1 * pres pgas(IHNO3) =q(ng+IHNO3) * 1.01325d-1 * pres pgas(IH2SO4)=q(ng+IH2SO4) * 1.01325d-1 * pres pgas(IHCL) =q(ng+IHCL) * 1.01325d-1 * pres c c Do now the organic gases (organics follow HCl) c do i=1,norg-1 kk=IHCL+i pgas(kk)=q(ng+kk) * 1.01325d-1 * pres enddo call flush(6) do i=1,nsecx if(qt(i).le.tinys) then ! no aerosol in this section--> do k=1,ngas hi(i,k)=0.0d0 ! no condensation/evaporation enddo else indx=(i-1)*nsp dp = 1.0d-6 * dsec(i+nsecx2) ! representative particle diameter ( m ) if(.not.dry(i)) then ! WET AEROSOL icase(i)=0 do k=1,IHCL ! inorganic call flush(6) KnD=2.0d0*lamda(k)/dp hi(i,k)=12.0d0*diffus(k)*(pgas(k)-ps(i,k))* & f(KnD, delta(k))*(gmw(k)*1.0d-3) / & (1.0d3*dp**2*rgas*temp) enddo c **** restrict NH3 transfer to force electroneutrality *********** c **** modified by bkoo (02/14/02) c GG=2.d0*c(i,6)/intmw(6)+c(i,7)/intmw(7) ! [umol/m3] & +c(i,8)/intmw(8)+c(i,5)/intmw(5)+c(i,19)/intmw(19) & -c(i,3)/intmw(3)-c(i,4)/intmw(4) c alpha=0.1*ABS(GG) ! in umoles/m3/s 10% change/s in H+ conc. c c First check for electroneutrality of the fluxes elec=qt(i)*(2*hi(i,ih2so4)/gmw(ih2so4)+hi(i,ihno3)/gmw(ihno3)+ & hi(i,ihcl)/gmw(ihcl)-hi(i,inh3)/gmw(inh3)) c c If not within alpha find xx. Iteration is occationally required c due to precission issues and the sensitivity of the new elec to xx c if(elec.lt.-alpha) then ! flux over basic limit if(hi(i,ihno3).lt.0.d0.and.hi(i,ihcl).lt.0.d0) then ! both NO3 & Cl evaporate x1=-alpha/qt(i)-2.d0*hi(i,ih2so4)/gmw(ih2so4) & +hi(i,inh3)/gmw(inh3) if(x1.eq.0.d0) then hi(i,ihno3)=0.d0 hi(i,ihcl)=0.d0 else xx=(hi(i,ihno3)/gmw(ihno3)+hi(i,ihcl)/gmw(ihcl))/x1 if(xx.ge.1.d0) then hi(i,ihno3)=hi(i,ihno3)/xx hi(i,ihcl)=hi(i,ihcl)/xx else hi(i,ihno3)=0.d0 hi(i,ihcl)=0.d0 endif endif elseif(hi(i,ihno3).lt.0.d0) then ! Cl condenses x1=-alpha/qt(i)-2.d0*hi(i,ih2so4)/gmw(ih2so4) & -hi(i,ihcl)/gmw(ihcl)+hi(i,inh3)/gmw(inh3) if(x1.eq.0.d0) then hi(i,ihno3)=0.d0 else xx=hi(i,ihno3)/gmw(ihno3)/x1 if(xx.ge.1.d0) then hi(i,ihno3)=hi(i,ihno3)/xx else hi(i,ihno3)=0.d0 endif endif elseif(hi(i,ihcl).lt.0.d0) then ! NO3 condenses x1=-alpha/qt(i)-2.d0*hi(i,ih2so4)/gmw(ih2so4) & -hi(i,ihno3)/gmw(ihno3)+hi(i,inh3)/gmw(inh3) if(x1.eq.0.d0) then hi(i,ihcl)=0.d0 else xx=hi(i,ihcl)/gmw(ihcl)/x1 if(xx.ge.1.d0) then hi(i,ihcl)=hi(i,ihcl)/xx else hi(i,ihcl)=0.d0 endif endif endif ! both NO3 & Cl evaporate elseif(elec.gt.alpha .and. hi(i,inh3).lt.0.d0) then ! flux over acidic limit x1=2.d0*hi(i,ih2so4)/gmw(ih2so4)+ ! & NH4 evaporates & hi(i,ihno3)/gmw(ihno3)+hi(i,ihcl)/gmw(ihcl)-alpha/qt(i) if(x1.eq.0.d0) then hi(i,inh3)=0.d0 else xx=hi(i,inh3)/gmw(inh3)/x1 if(xx.ge.1.d0) then hi(i,inh3)=hi(i,inh3)/xx else hi(i,inh3)=0.d0 endif endif endif ! flux over basic limit c elec2=qt(i)*(2.d0*hi(i,ih2so4)/gmw(ih2so4)+hi(i,ihno3)/ & gmw(ihno3)+hi(i,ihcl)/gmw(ihcl)-hi(i,inh3)/gmw(inh3)) if(elec2.lt.-alpha) then hmax=12.d0*diffus(inh3)*pgas(inh3) & *f(2.d0*lamda(inh3)/dp,delta(inh3)) & *(gmw(inh3)*1.d-3)/(1.d3*dp**2*rgas*temp) hi(i,inh3)=MIN((2.d0*hi(i,ih2so4)/gmw(ih2so4) & +hi(i,ihno3)/gmw(ihno3)+hi(i,ihcl)/gmw(ihcl) & +alpha/qt(i))*gmw(inh3),hmax) elseif(elec2.gt.alpha) then hmax1=12.d0*diffus(ihno3)*pgas(ihno3) & *f(2.d0*lamda(ihno3)/dp,delta(ihno3)) & *(gmw(ihno3)*1.d-3)/(1.d3*dp**2*rgas*temp) hmax2=12.d0*diffus(ihcl)*pgas(ihcl) & *f(2.d0*lamda(ihcl)/dp,delta(ihcl)) & *(gmw(ihcl)*1.d-3)/(1.d3*dp**2*rgas*temp) if(hi(i,ihno3).gt.0.d0 .and. hi(i,ihcl).gt.0.d0) then wf=hi(i,ihcl)/gmw(ihcl) / & (hi(i,ihcl)/gmw(ihcl)+hi(i,ihno3)/gmw(ihno3)) hi(i,ihno3)=MIN((alpha/qt(i)+hi(i,inh3)/gmw(inh3) & -2.d0*hi(i,ih2so4)/gmw(ih2so4))*gmw(ihno3) & *(1-wf),hmax1) hi(i,ihcl)=MIN((alpha/qt(i)+hi(i,inh3)/gmw(inh3) & -2.d0*hi(i,ih2so4)/gmw(ih2so4))*gmw(ihcl) & *wf,hmax2) elseif(hi(i,ihno3).gt.0.d0) then hi(i,ihno3)=MIN((alpha/qt(i)+hi(i,inh3)/gmw(inh3) & -2.d0*hi(i,ih2so4)/gmw(ih2so4) & -hi(i,ihcl)/gmw(ihcl))*gmw(ihno3),hmax1) elseif(hi(i,ihcl).gt.0.d0) then hi(i,ihcl)=MIN((alpha/qt(i)+hi(i,inh3)/gmw(inh3) & -2.d0*hi(i,ih2so4)/gmw(ih2so4) & -hi(i,ihno3)/gmw(ihno3))*gmw(ihcl),hmax2) endif endif ! dbg c ******************************************************************* else ! DRY AEROSOL call dryin(i,dp,ht,pgas) ! dry inorganics endif ! if WET or DRY c c if very small quantities are evaporating - inorganics (02/14/02) c if(q(indx+KNO3).lt.tinys .and. hi(i,ihno3).lt.0.d0) then ! KNO3 tmp=hi(i,inh3)-hi(i,ihno3)/gmw(ihno3)*gmw(inh3) if((q(indx+KNH4).ge.tinys .and. tmp.le.0.d0) .or. & (tmp.ge.0.d0 .and. q(ng+inh3).ge.tiny2)) then hi(i,inh3)=tmp else hi(i,inh3)=0.d0 if(q(indx+KCL).ge.tinys) then hi(i,ihcl)=-2.d0*hi(i,ih2so4)/gmw(ih2so4)*gmw(ihcl) else hi(i,ihcl)=0.d0 endif endif hi(i,ihno3)=0.d0 endif ! KNO3 if(q(indx+KCL) .lt.tinys .and. hi(i,ihcl) .lt.0.d0) then ! KCL tmp=hi(i,inh3)-hi(i,ihcl)/gmw(ihcl)*gmw(inh3) if((q(indx+KNH4).ge.tinys .and. tmp.le.0.d0) .or. & (tmp.ge.0.d0 .and. q(ng+inh3).ge.tiny2)) then hi(i,inh3)=tmp else hi(i,inh3)=0.d0 if(q(indx+KNO3).ge.tinys) then hi(i,ihno3)=-2.d0*hi(i,ih2so4)/gmw(ih2so4)*gmw(ihno3) else hi(i,ihno3)=0.d0 endif endif hi(i,ihcl)=0.d0 endif ! KCL ! now, "hi < 0" means "q(aer) >= tinys" for HNO3 & HCl if(q(indx+KNH4).lt.tinys .and. hi(i,inh3) .lt.0.d0) then ! KNH4 if(hi(i,ihno3).lt.0.d0 .and. hi(i,ihcl).lt.0.d0) then wf=hi(i,ihcl)/gmw(ihcl) / & (hi(i,ihcl)/gmw(ihcl)+hi(i,ihno3)/gmw(ihno3)) hi(i,ihcl)=hi(i,ihcl)-hi(i,inh3)/gmw(inh3)*gmw(ihcl)*wf hi(i,ihno3)=hi(i,ihno3)-hi(i,inh3)/gmw(inh3)*gmw(ihno3) & *(1-wf) if(hi(i,ihcl).gt.0.d0 .and. q(ng+ihcl).lt.tiny2) then hi(i,ihcl)=0.d0 hi(i,ihno3)=-2.d0*hi(i,ih2so4)/gmw(ih2so4)*gmw(ihno3) endif if(hi(i,ihno3).gt.0.d0 .and. q(ng+ihno3).lt.tiny2) then hi(i,ihno3)=0.d0 hi(i,ihcl)=-2.d0*hi(i,ih2so4)/gmw(ih2so4)*gmw(ihcl) endif elseif(hi(i,ihno3).lt.0.d0) then hi(i,ihno3)=hi(i,ihno3)-hi(i,inh3)/gmw(inh3)*gmw(ihno3) if(hi(i,ihno3).gt.0.d0 .and. q(ng+ihno3).lt.tiny2) then hi(i,ihno3)=0.d0 hi(i,ihcl)=-2.d0*hi(i,ih2so4)/gmw(ih2so4)*gmw(ihcl) endif elseif(hi(i,ihcl).lt.0.d0) then hi(i,ihcl)=hi(i,ihcl)-hi(i,inh3)/gmw(inh3)*gmw(ihcl) if(hi(i,ihcl).gt.0.d0 .and. q(ng+ihcl).lt.tiny2) then hi(i,ihcl)=0.d0 hi(i,ihno3)=-2.d0*hi(i,ih2so4)/gmw(ih2so4)*gmw(ihno3) endif endif hi(i,inh3)=0.d0 endif ! KNH4 c Special case for if a dry particle with NaNO3 and case 1.1.3 runs out of c NH4NO3 to evaporate. In this case NO3 from NaNO3 will be lost from the aerosol c and NaCl will form which will drive towards a different equilibrium. c To prevent the inappropriate formation of NaCl when the NH4NO3 runs out c we will set the NO3 flux to zero if(icase(i).eq.113) then if(hi(i,ihno3).lt.0.d0 .and. c(i,13).le.tinys) then hi(i,inh3)=hi(i,inh3)-hi(i,ihno3)/gmw(ihno3)*gmw(inh3) hi(i,ihno3)=0.d0 if(hi(i,inh3).gt.0.d0 .and. q(ng+inh3).lt.tiny2) then hi(i,inh3)=0.d0 if(q(indx+KCL).ge.tinys) then hi(i,ihcl)=-2.d0*hi(i,ih2so4)/gmw(ih2so4)*gmw(ihcl) else hi(i,ihcl)=0.d0 endif endif endif endif c Special case: A dry particle contains NH4Cl and uses case 1.1.2 but runs out of c NaCl to evaporate. In this case NH4NO3 will form from NaNO3 regardless of the c potential of NH4NO3 to evaporate To prevent the c inappropriate formation of NH4NO3 when the NaCl runs out we will set the Cl flux c to equal the ammonia flux if(icase(i).eq.112) then if(hi(i,ihcl).lt.0.d0 .and. c(i,9).le.tinys) then hi(i,ihcl)=hi(i,inh3)/gmw(inh3)*gmw(ihcl) if(q(indx+KNO3).ge.tinys) then hi(i,ihno3)=-2.d0*hi(i,ih2so4)/gmw(ih2so4)*gmw(ihno3) else hi(i,ihno3)=0.d0 endif if(hi(i,ihcl).gt.0.d0 .and. q(ng+ihcl).lt.tiny2) then hi(i,ihcl)=0.d0 hi(i,inh3)=0.d0 endif endif endif c do k=1,norg-1 ! do organics kk=IHCL+k KnD=2.0d0*lamda(kk)/dp hi(i,kk)=12.0d0*diffus(kk)*(pgas(kk)-ps(i,kk))* & f(KnD, delta(kk))*(gmw(kk)*1.0d-3) / & (1.0d3*dp**2*rgas*temp) c c if very small quantities are evaporating - organics (02/14/02) c ii=nexti+k if(q(indx+ii).lt.tinys .and. hi(i,kk).lt.0.d0) hi(i,kk)=0.d0 enddo endif ! no aerosol in this section enddo ! do nsecx c ************************************************************** c MODIFIED this section to calculate as in Chock and Winkler c tmg(07/02/01) c Take the derivative of the total flux for each particle c Calculate the derivatives using Equation 11 of Pilinis et al. (1987) c c do i=1,nsecx ht(i)=0.0d0 do k=1,ngas ht(i)=ht(i)+hi(i,k) enddo enddo c makes sure previous step wasn't too large (tmg, 07/02/01) if (tpres.gt.t) then do i=1,nsecx htavg(i)=(ht(i)+htlast(i))/2. merr(i) = qt(i)*(htavg(i)-htlast(i))*tmincr ! tmg (11/20/02) masserr = max(masserr,abs(merr(i))) enddo c if error small, increase tmincr if (masserr.lt.ltol/100) then tmincr= tmincr*2.0d0 else if (masserr.lt.errtol/10.) then tmincr = tmincr*1.1d0 endif c if error large cut timestep in half and redo last iteration if (masserr.gt.errtol.and.tmincr.gt.1.0d-4) then 30 tmincr = tmincr/2.0d0 do k=1,ntotal q(k) = qlast(k) enddo do i=1,nsecx ht(i)=htlast(i) qt(i)=qtlast(i) do k=1,ngas hi(i,k)=hilast(i,k) enddo enddo tpres=tlast else if (tpres.eq.(tout)) then 40 t=tpres return endif endif c make sure tout will not be exceeded if ((tpres+tmincr).gt.(tout)) then tmincr= tout-tpres endif masserr = 0.0d0 c ************************************************************** c MODIFIED this section to calculate as in Chock and Winkler c tmg(07/02/01) c Take the derivative of the total flux for each particle c Calculate the derivatives using Equation 11 of Pilinis et al. (1987) c do i=1,nsecx indx=(i-1)*nsp c dqdt(indx+KH2O)= 0.0d0 c dqdt(indx+KNa) = 0.0d0 dqdt(indx+KSO4)= hi(i,IH2SO4)*qt(i) dqdt(indx+KNO3)= hi(i,IHNO3) *qt(i) dqdt(indx+KNH4)= hi(i,INH3) *qt(i) dqdt(indx+KCL) = hi(i,IHCL) *qt(i) do k=1,norg-1 ii=nexti+k kk=IHCL+k dqdt(indx+ii) = hi(i,kk) *qt(i) enddo enddo c c Done with aerosol rates. Start gases. c c changed according to T-G model (tmg, 07/02/01) c initialize do i = 1,ngas dq(ng+i) = 0.0d0 enddo do i=1,nsecx dqdt(ng+IH2SO4)=-hi(i,IH2SO4)*qt(i) dqdt(ng+IHNO3) =-hi(i,IHNO3) *qt(i) dqdt(ng+INH3) =-hi(i,INH3) *qt(i) dqdt(ng+IHCL) =-hi(i,IHCL) *qt(i) if (ht(i) .ne. 0.0d0) then dq(ng+IH2SO4)= dq(ng+IH2SO4) + & dqdt(ng+IH2SO4)/ht(i)*(exp(ht(i)*tmincr)-1.0d0) dq(ng+IHNO3) = dq(ng+IHNO3) & +dqdt(ng+IHNO3)/ht(i)*(exp(ht(i)*tmincr)-1.0d0) dq(ng+INH3) = dq(ng+INH3) & +dqdt(ng+INH3)/ht(i)*(exp(ht(i)*tmincr)-1.0d0) dq(ng+IHCL) = dq(ng+IHCL) & +dqdt(ng+IHCL)/ht(i)*(exp(ht(i)*tmincr)-1.0d0) else dq(ng+IH2SO4) = dq(ng+IH2SO4)+dqdt(ng+IH2SO4)*tmincr dq(ng+IHNO3) = dq(ng+IHNO3)+dqdt(ng+IHNO3)*tmincr dq(ng+INH3) = dq(ng+INH3)+dqdt(ng+INH3)*tmincr dq(ng+IHCL) = dq(ng+IHCL)+dqdt(ng+IHCL)*tmincr endif do k=1,norg-1 ii=ng+IHCL+k km=IHCL+k dqdt(ii) =-hi(i,km)*qt(i) if (ht(i) .ne. 0.0d0) then dq(ii) = dq(ii) & +dqdt(ii)/ht(i)*(exp(ht(i)*tmincr)-1.0d0) else dq(ii) = dq(ii) + dqdt(ii)*tmincr endif enddo enddo c c Check derivative balance to ensure no change in mass will occur c Ignore for T-G method (tmg, 07/02/01) c To change the order of the summation later... (bkoo) dqdtsum=0.0d0 c do i=1,neq c dqdtsum=dqdtsum+dqdt(i) c enddo c c -- an error of 1e-5 ug/m3/s will produce 1ug/m3 in 24 hours -- if(abs(dqdtsum).ge.1e-5) then cgy write(6,*)'mass error of ',dqdtsum,' ug/m3/s occured at ',t cgy write(90,*)'mass error of ',dqdtsum,' ug/m3/s occured at ',t do ispec=1,nsp dqsum(ispec)=0.0d0 do isec=1,nsecx dqsum(ispec)=dqsum(ispec)+dqdt((isec-1)*nsp+ispec) enddo if(ispec.eq.kso4) dqsum(ispec)=dqsum(ispec)+dqdt(ng+ih2so4) if(ispec.eq.kno3) dqsum(ispec)=dqsum(ispec)+dqdt(ng+ihno3) if(ispec.eq.knh4) dqsum(ispec)=dqsum(ispec)+dqdt(ng+inh3) if(ispec.eq.kcl) dqsum(ispec)=dqsum(ispec)+dqdt(ng+ihcl) if(ispec.gt.kcl) then dqsum(ispec)=dqsum(ispec)+dqdt(ng+ihcl+ispec-kcl) endif ctmg do isec=1,nsecx cgy write(90,21)t,ispec,isec+neqsec,dqsum(ispec)/emw(ispec), cgy & dqdt((isec-1)*nsp+ispec)/emw(ispec) ctmg enddo enddo endif 21 format('dqdtsum error: ',e14.5,2i5,2e14.5) c c Add gas phase source term in ug/m3/sec c Modified according to T-G method (tmg, 07/02/01) do igas=1,ngas dqdt(ng+igas)=gsource(igas) enddo c c Change the units of the gases from ug/m3/sec to ppm/sec c add pres (bkoo, 06/09/00) do i=1,ngas dqdt(ng+i)=dqdt(ng+i)*rgas*temp/(1.01325d5*pres*gmw(i)) dq(ng+i)=dq(ng+i)*rgas*temp/(1.01325d5*pres*gmw(i)) enddo c additional code for T-G method (tmg, 07/02/01) do k=1,ntotal qlast(k) = q(k) enddo c update concentrations according to new method do i = 1,ngas q(ng+i) = q(ng+i) + dq(ng+i)+ dqdt(ng+i)*tmincr enddo do i=1,nsecx htlast(i) = ht(i) qtlast(i) = qt(i) do k=1,ngas hilast(i,k)=hi(i,k) enddo enddo do i=1,nsecx do k=2,nsp if (ht(i).ne.0.0d0) then q((i-1)*nsp+k)=q((i-1)*nsp+k)+dqdt((i-1)*nsp+k)/ht(i)* & (exp(ht(i)*tmincr)-1.0d0) else q((i-1)*nsp+k)=q((i-1)*nsp+k)+dqdt((i-1)*nsp+k)*tmincr endif enddo enddo cbk if(aerm.eq.'HYBR') then cbk call negchk(tpres,q,nsecd) cbk else cbk call negchk(tpres,q,nsec) cbk endif call negchk(tpres,q,nsecx) ! bkoo (03/07/03) icount = icount + 1 ! tmg (12/17/02) tlast = tpres tpres = tpres+tmincr go to 20 end subroutine dryin(i,dp,ht,pgas) c This subroutine calculates the intrasectional c coefficients for the 4 inorganic condensible gases c in the subcase of dry aerosols, for the new c multicomponent aerosol dynamics model (MADM) c ************************************************************************* c WRITTEN BY DR. CHRISTODOULOS PILINIS c January 1998 c ************************************************************************* c MODIFIED TO INCLUDE SUPERSATURATION FLAGS BY C KEVIN P. CAPALDO C JUNE 1998 c ************************************************************************* c For the Ammonia rich case 9 subcases are used which utilize at most 2 c of the following three reversible reactions c (1) NH4NO3(s) <=> NH3(g) + HNO3(g) c (2) NH4Cl(s) <=> NH3(g) + HCl(g) c (3) NaCl(s) + HNO3(g) <=> NaNO3(s) + HCl(g) c Thermodynamically all three equations cannot be in equilibrium simultaneously c The choice of which reactions to use is determined by the compounds present c in the aerosol phase and the gas phase. Because all 4 of the salts NaCl, NaNO3, c NH4Cl, and NH4NO3 cannot exist at the same time in the aerosol phase, c the assumption that NH4NO3 and NaCl will never coexist (these compounds will c react to form NaNO3 and NH4Cl) is used. Because of this assumption, c supersaturations in the gas phase do not always invoke the use of c the corresponding reaction. c include 'dynamic.inc' include 'equaer.inc' dimension pgas(ngas) ! partial pressures of condensible gase (Pa) dimension ht(nsecx) ! total rate of change integer icase(nsec) common /drycase/ icase c Equation for transition regime (Dahneke 1983) f(yKnD, alfa)=(1.0d0+yKnD)/(1.0d0+2.0d0*yKnD*(1.0d0+yKnD)/alfa) c tinys =1.0d-9 ! minimum mass in ug/m3 if(qt(i).le.tinys) then ! no aerosol in this section----> hi(i,IH2SO4)=0.0d0 ! no condensation/evaporation hi(i,INH3) =0.0d0 hi(i,IHNO3) =0.0d0 hi(i,IHCL) =0.0d0 goto 200 endif c change the units for R to use it for equilibrium constants RGAS1=RGAS/(1.01325d5*pres) ! add pres (bkoo, 06/09/00) hno3kn=2.0d0*lamda(ihno3)/dp fhno3=f(hno3kn, delta(ihno3)) dhno3=diffus(ihno3)*fhno3 ammokn=2.0d0*lamda(inh3)/dp fammo=f(ammokn, delta(inh3)) dammo=diffus(inh3)*fammo hclkn=2.0d0*lamda(ihcl)/dp fhcl=f(hclkn, delta(ihcl)) dhcl=diffus(ihcl)*fhcl sulkn=2.0d0*lamda(IH2SO4)/dp fsul=f(sulkn, delta(IH2SO4)) dsul=diffus(IH2SO4)*fsul c c Calculate the rate constants and flags for condensation c Make the reaction consistent with the reactions used c by ISORROPIA. The equilibrium constants used by MADM c are caclulated as a combination of the equilibrium c constants used by ISORROPIA. c RKNH4NO3=xk10/(rgas1*temp)**2 RKNA=XK4*XK8/(XK3*XK9) RKNH4CL=xk6/(rgas1*temp)**2 satnh4no3=pgas(ihno3)*pgas(inh3)/(rgas*temp)**2 satnh4cl=pgas(ihcl)*pgas(inh3)/(rgas*temp)**2 if(satnh4no3.gt.rknh4no3) then ifgnh4no3 = 1 ! NH4NO3 forms else ifgnh4no3 = 0 endif if(satnh4cl.gt.rknh4cl) then ifgnh4cl = 1 ! NH4CL forms else ifgnh4cl = 0 endif if(pgas(ihcl).gt.rkna*pgas(ihno3)) then ifgna = 1 ! NaCl forms if NaNO3 exists else ifgna = 0 ! NaNO3 forms if NaCl exists endif c if(ifgnh4no3.eq.1.and.ifgnh4cl.eq.0.and.ifgna.eq.1) then cgy write(90,*)'error in gas phase conc. (or solid eq rates)' write(6,*)'error in gas phase conc. (or solid eq rates)' stop elseif(c(i,9).gt.tinys.and.c(i,13).gt.tinys) then cgy write(90,*)'error: NaCl and NH4NO3 should not coexist' write(6,*)'error: NaCl and NH4NO3 should not coexist' stop elseif(c(i,12).gt.tinys.and. & (c(i,9).gt.tinys.or.c(i,11).gt.tinys)) then cgy write(90,*)'(NH4)2SO4 and NaCl or NaNO3 should not coexist' write(6,*)'(NH4)2SO4 and NaCl or NaNO3 should not coexist' stop endif c k=IH2SO4 ! do H2SO4 hi(i,k)=12.0d0*dsul*pgas(k)* & (gmw(k)*1.0d-3)/(1.0d3*dp**2*rgas*temp) c c Now find the subcases to calculate the fluxes c for HNO3, NH3 and HCl c c Ammonia poor case if(c(i,16).gt.tinys.or.c(i,17).gt.tinys.or.c(i,18).gt.tinys) & goto 120 ! Acidic aerosol c c Ammonia Rich cases c c No Sodium (except Na2SO4) if(c(i,9).le.tinys.and.c(i,11).le.tinys) then ! ~NaCl and ~NaNO3 if(c(i,13).gt.tinys.or.ifgnh4no3.eq.1) then ! NH4NO3 if(c(i,14).gt.tinys.or.ifgnh4cl.eq.1) goto 113 ! NH4Cl goto 116 ! ~NH4Cl elseif(c(i,14).gt.tinys.or.ifgnh4cl.eq.1) then goto 117 ! ~NH4NO3 and NH4Cl else goto 114 ! ~NH4NO3 and ~NH4Cl endif c NaCl present elseif(c(i,9).gt.tinys) then ! NaCl if(c(i,11).le.tinys.and.ifgna.eq.1) goto 117 ! ~NaNO3 goto 112 ! NaNO3 c NaNO3 present else ! NaNO3 if(c(i,13).gt.tinys) then ! NH4NO3 and NaCl wont form if(c(i,14).le.tinys.and.ifgnh4cl.eq.0.and.ifgna.eq.0) & goto 116 ! ~NH4Cl goto 113 ! NH4Cl elseif(ifgna.eq.1) then ! NaCl may form if(ifgnh4no3.eq.1) goto 111 ! NH4NO3 may form goto 112 ! NH4NO3 may form elseif(c(i,14).gt.tinys.or.ifgnh4cl.eq.1) then ! NH4Cl goto 113 else ! ~NH4Cl goto 116 endif endif c c NH3 poor case, i.e. acidic aerosol (case 1.2) c 120 icase(i)=120 cgy if(c(i,9).gt.tinys.or.c(i,11).gt.tinys.or.c(i,13).gt.tinys.or. cgy # c(i,14).gt.tinys)write(90,*)'ERROR: SALTS AND ACID PRESENT ', cgy # c(i,9),c(i,11),c(i,13),c(i,14) c k=INH3 ! do NH3 hi(i,k)=12.0d0*dammo*pgas(k)*(gmw(k)*1.0d-3)/ # (1.0d3*dp**2*rgas*temp) hi(i,IHCL)=0.0d0 hi(i,IHNO3)=0.0d0 goto 200 c c case 1.1.1 *************************************************** c Reactions 1 and 3 c No NH4Cl but NaNO3, NaCl, and NH4NO3 111 icase(i)=111 c this case is not applied because NH4NO3 and NaCl cannot coexist c c However, for the specific cases where NaCl and NH4NO3 are absent from c the aerosol phase but are suggested to form by supersaterations c of the gas phase components of reactions 1 and 3 (ifgnh4no3=1 and c ifgna=1), we must determine which compound (NaCl or NH4NO3) will form. c This is done by determining whether NO3 will condense or evaporate c if both reactions 1 and 3 are employed. If case 1.1.1 indicates net c condensation of NO3 then any NaCl that forms will be imeadiately c converted to NaNO3 and NH4Cl, if net evaporation of NO3 is indicated c then reation 3 dominates and the NaCl formed is sufficient to react c with any NH4NO3 that condenses and keep the NH4NO3 concentration zero. c que=(dammo*pgas(inh3)-dhno3*pgas(ihno3)-dhcl*pgas(ihcl) # -2.0d0*dsul*pgas(ih2so4))/(rgas*temp) c CHNO3=(-que+sqrt(que**2+4.0*dammo*RKNH4NO3*(dhno3+RKNA*dhcl # )))/(2.0d0*(dhno3+RKNA*dhcl)) CNH3=RKNH4NO3/CHNO3 CHCL=RKNA*CHNO3 k=IHNO3 ! do HNO3 hi(i,k)=12.0d0*dhno3*(pgas(k)-CHNO3*rgas*temp)* & (gmw(k)*1.0d-3)/(1.0d3*dp**2*rgas*temp) k=INH3 ! do NH3 hi(i,k)=12.0d0*dammo*(pgas(k)-CNH3*rgas*temp)* & (gmw(k)*1.0d-3)/(1.0d3*dp**2*rgas*temp) k=IHCL ! do HCL hi(i,k)=12.0d0*dhcl*(pgas(k)-CHCL*rgas*temp)* & (gmw(k)*1.0d-3)/(1.0d3*dp**2*rgas*temp) c if NO3 evaporates then NaCl forms and NH4NO3 does not form if(hi(i,ihno3).lt.0.0d0) goto 112 c else NH4NO3 forms and NaCl does not form goto 113 c case 1.1.2 *************************************************** c Reactions 2 and 3 c No NH4NO3 but NaCl, NaNO3, and NH4Cl can all occur 112 icase(i)=112 c que=(dammo*pgas(inh3)-dhno3*pgas(ihno3)-dhcl*pgas(ihcl) # -2.0d0*dsul*pgas(ih2so4))/(rgas*temp) CHNO3=(-que+sqrt(que**2+4.0*dammo*RKNH4CL*(dhno3+RKNA*dhcl # )/RKNA))/(2.0d0*(dhno3+RKNA*dhcl)) CNH3=RKNH4CL/RKNA/CHNO3 CHCL=RKNA*CHNO3 k=IHNO3 ! do HNO3 hi(i,k)=12.0d0*dhno3*(pgas(k)-CHNO3*rgas*temp)* & (gmw(k)*1.0d-3)/(1.0d3*dp**2*rgas*temp) k=INH3 ! do NH3 hi(i,k)=12.0d0*dammo*(pgas(k)-CNH3*rgas*temp)* & (gmw(k)*1.0d-3)/(1.0d3*dp**2*rgas*temp) k=IHCL ! do HCL hi(i,k)=12.0d0*dhcl*(pgas(k)-CHCL*rgas*temp)* & (gmw(k)*1.0d-3)/(1.0d3*dp**2*rgas*temp) c c if NaCl doesn't exist and case 1.1.2 has been employed then c both NaNO3 must exist and reaction 3 indicates formation of c NaCl (ifgna=1). If H2SO4 condenses then we must determine whether c there is sufficient NaCl formed to react with the condensing C H2SO4 (to form Na2SO4 and NH4Cl) or if NaNO3 reacts with H2SO4 c (and NH4NO3 is formed). If the later is the case then 1.1.3 c is required. To determine this the HCl flux is compared to the c NH3 flux. If JNH3 > JHCl (in molar units) then NaCl does not form. if(c(i,9).le.tinys.and. & (hi(i,inh3)/gmw(inh3)).gt.(hi(i,ihcl)/gmw(ihcl))) goto 113 c c if there is no NH4 to evaporate then if(c(i,14).le.tinys.and.hi(i,inh3).lt.0.0d0) goto 115 c if there is no NO3 to evaporate then if(c(i,11).le.tinys.and.hi(i,ihno3).lt.0.0d0) goto 117 c if there is no Cl to evaporate then if(c(i,9).le.tinys.and.c(i,14).le.tinys.and. & hi(i,ihcl).lt.0.0d0) goto 116 goto 200 c c case 1.1.3 ********************************************* c Reactios 1 and 2 c No NaCl but NH4Cl, NH4NO3, and NaNO3 can occur 113 icase(i)=113 que=(dammo*pgas(inh3)-dhno3*pgas(ihno3)-dhcl*pgas(ihcl) # -2.0d0*dsul*pgas(ih2so4))/(rgas*temp) CHNO3=(-que+sqrt(que**2+4.0*dammo*RKNH4NO3*(dhno3+dhcl # *RKNH4CL/RKNH4NO3)))/(2.0d0*(dhno3+dhcl*RKNH4CL/RKNH4NO3)) CNH3=RKNH4NO3/CHNO3 CHCL=CHNO3*RKNH4CL/RKNH4NO3 k=IHNO3 ! do HNO3 hi(i,k)=12.0d0*dhno3*(pgas(k)-CHNO3*rgas*temp)* & (gmw(k)*1.0d-3)/(1.0d3*dp**2*rgas*temp) k=INH3 ! do NH3 hi(i,k)=12.0d0*dammo*(pgas(k)-CNH3*rgas*temp)* & (gmw(k)*1.0d-3)/(1.0d3*dp**2*rgas*temp) k=IHCL ! do HCL hi(i,k)=12.0d0*dhcl*(pgas(k)-CHCL*rgas*temp)* & (gmw(k)*1.0d-3)/(1.0d3*dp**2*rgas*temp) c c if there is no NH4 to evaporate then if(c(i,13).le.tinys.and.c(i,14).le.tinys.and. & hi(i,inh3).lt.0.0d0) then if(c(i,11).le.tinys) goto 114 goto 119 endif c if there is no Cl to evaporate then if(c(i,14).le.tinys.and.hi(i,ihcl).lt.0.0d0) goto 116 c if there is no NO3 to evaporate then if(c(i,13).le.tinys.and.c(i,11).le.tinys.and. & hi(i,ihno3).lt.0.0d0) goto 117 goto 200 c c case 1.1.4 ********************************************* c None of reactions 1, 2, or 3 c Sulfates only (no nitrates or chlorides) 114 icase(i)=114 k=INH3 ! do NH3 hi(i,k)=2.0d0*hi(i,IH2SO4)*(gmw(k)/gmw(IH2SO4)) hi(i,k)=min(hi(i,k),12.0d0*dammo*pgas(k)*(gmw(k)*1.0d-3)/ # (1.0d3*dp**2*rgas*temp)) hi(i,IHCL)=0.0d0 hi(i,IHNO3)=0.0d0 goto 200 c c case 1.1.5 ********************************************* c Reaction 3 only with the flux of NH3 set to zero c Sodium system (No ammonium salts) 115 icase(i)=115 ss=(dhno3*pgas(ihno3)+dhcl*pgas(ihcl) # +2.0d0*dsul*pgas(ih2so4))/(rgas*temp) CHNO3=ss/(dhno3+dhcl*RKNA) CHCL=CHNO3*RKNA k=IHNO3 ! do HNO3 hi(i,k)=12.0d0*dhno3*(pgas(k)-CHNO3*rgas*temp)* & (gmw(k)*1.0d-3)/(1.0d3*dp**2*rgas*temp) k=INH3 ! do NH3 hi(i,k)=0.0d0 k=IHCL ! do HCL hi(i,k)=12.0d0*dhcl*(pgas(k)-CHCL*rgas*temp)* & (gmw(k)*1.0d-3)/(1.0d3*dp**2*rgas*temp) c We assume here that case 1.1.5 is only called through case 1.1.2 c if there is no Cl to evaporate then if(c(i,9).le.tinys.and.hi(i,ihcl).lt.0.0d0) goto 119 c if there is no NO3 to evaporate then if(c(i,11).le.tinys.and.hi(i,ihno3).lt.0.0d0) goto 118 goto 200 c c case 1.1.6 ********************************************* c Only reaction 1 c Nitrate system (no Chloride) 116 icase(i)=116 p=(dammo*pgas(inh3)-dhno3*pgas(ihno3) # -2.0d0*dsul*pgas(ih2so4))/(rgas*temp) CHNO3=(-p+sqrt(p**2+4.0*dammo*RKNH4NO3*dhno3)) # /(2.0d0*dhno3) CNH3=RKNH4NO3/CHNO3 k=IHNO3 ! do HNO3 hi(i,k)=12.0d0*dhno3*(pgas(k)-CHNO3*rgas*temp)* & (gmw(k)*1.0d-3)/(1.0d3*dp**2*rgas*temp) k=INH3 ! do NH3 hi(i,k)=12.0d0*dammo*(pgas(k)-CNH3*rgas*temp)* & (gmw(k)*1.0d-3)/(1.0d3*dp**2*rgas*temp) k=IHCL ! do HCL hi(i,k)=0.0d0 c c if there is no NH4 to evaporate then if(c(i,13).le.tinys.and.hi(i,inh3).lt.0.0d0) goto 119 c if there is no NO3 to evaporate then if(c(i,13).le.tinys.and.c(i,11).le.tinys.and. & hi(i,ihno3).lt.0.0d0) goto 114 goto 200 c c case 1.1.7 ********************************************* c Only reaction 2 c Chloride system (no Nitrates) 117 icase(i)=117 p=(-dammo*pgas(inh3)+dhcl*pgas(ihcl) # +2.0d0*dsul*pgas(ih2so4))/(rgas*temp) RKNH4CL=xk6/(rgas1*temp)**2 CHNO3=0.0d0 CNH3=(-p+sqrt(p**2+4.0d0*dammo*dhcl*RKNH4CL))/(2.0*dammo) CHCL=RKNH4CL/CNH3 k=IHNO3 ! do HNO3 hi(i,k)=0.0d0 k=INH3 ! do NH3 hi(i,k)=12.0d0*dammo*(pgas(k)-CNH3*rgas*temp)* & (gmw(k)*1.0d-3)/(1.0d3*dp**2*rgas*temp) k=IHCL ! do HCL hi(i,k)=12.0d0*dhcl*(pgas(k)-CHCL*rgas*temp)* & (gmw(k)*1.0d-3)/(1.0d3*dp**2*rgas*temp) c if there is no NH4 to evaporate then if(c(i,14).le.tinys.and.hi(i,inh3).lt.0.0d0) goto 118 c if there is no Cl to evaporate then if(c(i,14).le.tinys.and.c(i,9).le.tinys.and. & hi(i,ihcl).lt.0.0d0) goto 114 goto 200 c c case 1.1.8 ********************************************* c none of reactions 1, 2, or 3 but JNH3 is set to zero c NaCl system 118 icase(i)=118 hi(i,INH3)=0.0d0 ! do NH3 hi(i,IHCL)=-2.0d0*hi(i,IH2SO4)/gmw(iH2SO4)*gmw(iHCL) ! HCL hi(i,IHNO3)=0.0d0 ! do HNO3 goto 200 c c case 1.1.9 ********************************************* c none of reactions 1, 2, or 3 but JNH3 is set to zero c NaNO3 system 119 icase(i)=119 hi(i,INH3)=0.0D0 ! do NH3 hi(i,IHNO3)=-2.0d0*hi(i,IH2SO4)/gmw(iH2SO4)*gmw(iHNO3) ! HNO3 hi(i,IHCL)=0.0d0 ! do HCL c 200 ht(i)=0.0d0 200 ht(i)=ht(i)+hi(i,IH2SO4)+hi(i,IHCL)+hi(i,IHNO3)+hi(i,INH3) c write(6,*)'section ',i,'case ',icase(i) return end