subroutine aqoperator2(tbeg,deltat,gas,aerosol,lwc, & rh,temp,p,iaq) c c TOTAL GAS AQUEOUS c----------------------------------------------------------------------- c 1. FORMALDEHYDE TOTAL 1. SO2 1. S(IV) c 2. FORMIC ACID TOTAL 2. H2O2 2. H2O2(aq) c 3. NH3 3. NITRATE c 4. CHLORIDE c 5. AMMONIUM c 6. SULFATE c 7. HSO5- c 8. HMSA c c c THE AQUEOUS-PHASE EQUATIONS ARE SOLVED FOR EACH SPECIES c c c.....Y MATRIX STRUCTURE EVERYTHING IS (ug/m3) c ONLY 21 ODEs ARE SOLVED FOR THE WHOLE DISTRIBUTION c c YAQ(1) = TOTAL FORMALDEHYDE c c YAQ(2) = SO2 (g) c YAQ(3) = H2O2 (g) c YAQ(4) = NH3 (g) c YAQ(5) = HCOOH (total) c c YAQ(6) = SO2 (aq) SECTION 1 c YAQ(7) = H2O2 (aq) SECTION 1 c YAQ(8) = NITRATE (aq) SECTION 1 c YAQ(9) = CHLORIDE (aq) SECTION 1 c YAQ(10) = AMMONIUM (aq) SECTION 1 c YAQ(11) = SULFATE (aq) SECTION 1 c YAQ(12) = HSO5- (aq) SECTION 1 c YAQ(13) = HMSA (aq) SECTION 1 c c YAQ(14) = SO2 (aq) SECTION 2 c YAQ(15) = H2O2 (aq) SECTION 2 c YAQ(16) = NITRATE (aq) SECTION 2 c YAQ(17) = CHLORIDE (aq) SECTION 2 c YAQ(18) = AMMONIUM (aq) SECTION 2 c YAQ(19) = SULFATE (aq) SECTION 2 c YAQ(20) = HSO5- (aq) SECTION 2 c YAQ(21) = HMSA (aq) SECTION 2 c include 'aerpar.inc' include 'droppar.inc' include 'dropcom.inc' real*4 gas(ngas_aq), aerosol(nsect,naers) real*4 sodium(nsect) real*4 deltat, rh, temp real*4 gascon(ngas_aq) real*4 yaq(meqn2) real*4 watcont(2) real*4 salt(2), crustal(2), lwc(2) real*4 check(3,naers) c common / gascon / gascon common / epcomy / ymin, hmax common / integration / wat1,relh1,therm,dt,dtemp common / aqrate / sodium, pres common / count/ icount, iprint common/aqchem2/watcont,crustal,salt,watvap c c INITIALIZATION (ONLY ON FIRST CALL) c (CALCULATION OF SECTIONAL DIAMETERS, LOADING OF MOLEC. WEIGHTS) c if (iaq .eq. 1) then call dropinit iaq = 0 endif pres = p ! pressure in atm call constants(temp) c c ZERO THE YAQ MATRIX c do i=1, meqn2 yaq(i)=0.0 enddo c c TRANSFER VIA COMMON VARIABLES TO AQRATE c watcont(1) = lwc(1) watcont(2) = lwc(2) relh1 = rh therm = temp dt = deltat icount=0 iprint=20 c c TRANSFER OF ALL GAS-PHASE CONCENTRATIONS TO GASCON AND THEN c THROUGH COMMON BLOCK TO RATES c do i=1, ngas_aq gascon(i) = gas(i) enddo c c UNIT CONVERSION FACTORS c fso2=1000.*p*wso2/(8.314e-2*temp) !SO2 conver. factor from ppm to ug/m3 fh2o2=1000.*p*wh2o2/(8.314e-2*temp) !H2O2 conver. factor from ppm to ug/m3 fhcho=1000.*p*whcho/(8.314e-2*temp) !HCHO conver. factor from ppm to ug/m3 fhcooh=1000.*p*whcooh/(8.314e-2*temp) !HCOOH conver. factor from ppm to ug/m3 fammon=1000.*p*wnh3/(8.314e-2*temp) !NH3 conver. factor from ppm to ug/m3 c c c *** BEGINNING OF AQUEOUS-PHASE CALCULATION c LOADING OF THE CONCENTRATIONS IN THE YAQ VECTOR c c THE TOTAL FORMALDEHYDE/FORMIC ACID IS TRANSFERRED AS GAS-PHASE SPECIES c (IT IS NOT INCLUDED IN THE AQUEOUS OR AEROSOL MATRICES) yaq(1) = gas(nghcho)*fhcho ! total HCHO (ug/m3) c c GAS-PHASE SPECIES c yaq(2) = gas(ngso2)*fso2 ! total SO2 (ug/m3) yaq(3) = gas(ngh2o2)*fh2o2 ! total H2O2 (ug/m3) yaq(4) = gas(nga) *fammon ! NH3(g) in ug/m3 yaq(5) = gas(nghcooh)*fhcooh ! HCOOH (total) (ug/m3) c c *** CALCULATION FOR FIRST AQUEOUS SECTION *** c c c A FIXED SECTIONAL DIAMETER IS USED FOR THE ACTIVATION c CALCULATION. THIS IS COMPATIBLE WITH THE ASSUMPTION THAT c A CLOUD OR FOG EXISTS IF THE RH EXCEEDS A THRESHOLD VALUE c (FOR EXAMPLE RH>85%). c c ** THE DROPLETS ARE ASSUMED TO BY WITHOUT SO2 OR c H2O2 IN THE BEGINNING OF A TIMESTEP ** c (THIS IS DONE TO AVOID CARRYING THE S(IV)/H2O2 CONCENTRATIONS c IN THE 3D SIMULATION ** yaq(6) = 1.e-10 ! SO2 (aq) in ug/m3 yaq(7) = 1.e-10 ! H2O2 (aq) in ug/m3 c do isect=1, nsect if (daer(isect) .ge. dactiv .AND. daer(isect) .lt. dsep) then yaq(8) = yaq(8) + aerosol(isect,nan) ! NITRATE (aq) in ug/m3 yaq(9) = yaq(9) + aerosol(isect,nac) ! CHLORIDE (aq) in ug/m3 yaq(10) = yaq(10) + aerosol(isect,naa) ! AMMONIUM (aq) in ug/m3 yaq(11) = yaq(11) + aerosol(isect,na4) ! SULFATE (aq) in ug/m c yaq(12) = yaq(12) + aerosol(isect,nahso5) ! HSO5- in ug/m3 yaq(13) = yaq(13) + aerosol(isect,nahmsa) ! HMSA in ug/m3 endif enddo c *** CALCULATION FOR SECOND AQUEOUS SECTION *** c yaq(14) = 1.e-10 ! SO2 (aq) in ug/m3 yaq(15) = 1.e-10 ! H2O2 (aq) in ug/m3 c do isect=1, nsect if (daer(isect) .ge. dsep) then yaq(16) = yaq(16) + aerosol(isect,nan) ! NITRATE (aq) in ug/m3 yaq(17) = yaq(17) + aerosol(isect,nac) ! CHLORIDE (aq) in ug/m3 yaq(18) = yaq(18) + aerosol(isect,naa) ! AMMONIUM (aq) in ug/m3 yaq(19) = yaq(19) + aerosol(isect,na4) ! SULFATE (aq) in ug/m3 yaq(20) = yaq(20) + aerosol(isect,nahso5) ! HSO5- in ug/m3 yaq(21) = yaq(21) + aerosol(isect,nahmsa) ! HMSA in ug/m3 endif enddo c c CALCULATION OF THE DISSOLUTION OF THE AVAILABLE HNO3 AND HCl c TO THE DROPLETS/PARTICLES. WE ASSUME THAT THE TIMESCALE OF DISSOLUTION c IS SMALL AND THAT ALL HNO3 AND HCl ARE DISSOLVED (ZERO VAPOR PRESSURE) c fhno3=1000.*whno3*pres/(8.314e-2*temp) !HNO3 conver. factor from ppm to ug/m3 fhcl=1000.*whcl*pres/(8.314e-2*temp) !HCl conver. factor from ppm to ug/m3 c hno3 = gas(ngcn)*fhno3 ! HNO3(g) (in ug/m3) hcl = gas(ngcc)*fhcl ! HCl(g) (in ug/m3) c c SAVE THE OLD AQUEOUS-PHASE CONCENTRATIONS BEFORE THE INTEGRATION c tnitold1= yaq(8) chlorold1= yaq(9) ammonold1= yaq(10) sulfateold1= yaq(11) hso5old1 = yaq(12) hmsaold1 = yaq(13) tnitold2= yaq(16) chlorold2= yaq(17) ammonold2= yaq(18) sulfateold2= yaq(19) hso5old2 = yaq(20) hmsaold2 = yaq(21) c c ** ADD THE GAS-PHASE CONCENTRATIONS TO THE TOTAL FOR HCl AND HNO3 ** c c * REMOVE IT FROM THE GAS-PHASE * gas(ngcn) = 0.0 ! all HNO3 is dissolved gas(ngcc) = 0.0 ! all HCl is dissolved c * CALCULATE PARTITIONING FACTORS * z1 = lwc(1) / diamet1 z2 = lwc(2) / diamet2 z = z1+z2 c fr1 = z1/z ! fraction that goes to section 1 fr2 = z2/z ! fraction that goes to section 2 c yaq(8) = yaq(8)+hno3*fr1 ! HNO3 increase Section 1 yaq(9) = yaq(9)+hcl*fr1 ! HCl increase Section 1 yaq(16) = yaq(16)+hno3*fr2 ! HNO3 increase Section 2 yaq(17) = yaq(17)+hcl*fr2 ! HCl increase Section 2 c c c.....THE YAQ MATRIX HAS BEEN INITIALIZED c CALCULATION OF SODIUM VECTOR (NEEDED FOR THE INTEGRATION ROUTINE) c (TRANSFERRED VIA COMMON STATEMENT AQRATES1) c do isect=1, nsect sodium(isect) = aerosol(isect,nas) enddo c c VARIABLES FOR INTEGRATION c stime = tbeg*60. ! in seconds stout = (tbeg+deltat) * 60. ! in seconds c c c CALCULATE CRUSTAL SPECIES CONCENTRATION INSIDE DROPLETS c (IT IS USED IN THE AQUEOUS-PHASE CHEMISTRY CALCULATION TO c ESTIMATE Fe and Mn CONCENTRATIONS c crustal(1) = aerosol(5,kcru)+aerosol(6,kcru)+aerosol(7,kcru) crustal(2) = aerosol(8,kcru)+aerosol(9,kcru)+aerosol(10,kcru) salt(1) = aerosol(5,ksod)+aerosol(6,ksod)+aerosol(7,ksod) salt(2) = aerosol(8,ksod)+aerosol(9,ksod)+aerosol(10,ksod) c c INTEGRATE c call aqintegr2(meqn2, yaq, stime, stout, temp, p) c c CALCULATE THE DIFFERENCES c dnit1 = yaq(8) - tnitold1 dchlor1 = yaq(9) - chlorold1 dammon1 = yaq(10) - ammonold1 dsulf1 = yaq(11) - sulfateold1 dhso51 = yaq(12) - hso5old1 dhmsa1 = yaq(13) - hmsaold1 c dnit2 = yaq(16) - tnitold2 dchlor2 = yaq(17) - chlorold2 dammon2 = yaq(18) - ammonold2 dsulf2 = yaq(19) - sulfateold2 dhso52 = yaq(20) - hso5old2 dhmsa2 = yaq(21) - hmsaold2 c c ADD THE MASS CHANGE TO THE CORRESPONDING SECTION c c do isect=5, 7 ! for sections 0.1 to 10 um do isect=4, 6 ! for sections 0.04 to 40 um aerosol(isect,nan)=aerosol(isect,nan)+dnit1*fdist2(isect) ! NITRATE (aq) in ug/m3 aerosol(isect,nac)=aerosol(isect,nac)+dchlor1*fdist2(isect) ! CHLORIDE (aq) in ug/m3 aerosol(isect,naa)=aerosol(isect,naa)+dammon1*fdist2(isect) ! AMMONIUM (aq) in ug/m3 aerosol(isect,na4)=aerosol(isect,na4)+dsulf1*fdist2(isect) ! SULFATE (aq) in ug/m3 aerosol(isect,nahso5)=aerosol(isect,nahso5)+dhso51*fdist2(isect) aerosol(isect,nahmsa)=aerosol(isect,nahmsa)+dhmsa1*fdist2(isect) enddo c c do isect=8, 10 ! for sections 0.1 to 10 um do isect=7, 8 ! for sections 0.04 to 40 um aerosol(isect,nan)=aerosol(isect,nan)+dnit2*fdist2(isect) ! NITRATE (aq) in ug/m3 aerosol(isect,nac)=aerosol(isect,nac)+dchlor2*fdist2(isect) ! CHLORIDE (aq) in ug/m3 aerosol(isect,naa)=aerosol(isect,naa)+dammon2*fdist2(isect) ! AMMONIUM (aq) in ug/m3 aerosol(isect,na4)=aerosol(isect,na4)+dsulf2*fdist2(isect) ! SULFATE (aq) in ug/m3 aerosol(isect,nahso5)=aerosol(isect,nahso5)+dhso52*fdist2(isect) aerosol(isect,nahmsa)=aerosol(isect,nahmsa)+dhmsa2*fdist2(isect) enddo c CHECK IF THE ALGORITHM RESULTED IN NEGATIVE CONCENTRATIONS c REPORT THE CORRECTIONS AT THE WARNING FILE c c check(1,isp) = total aerosol concentration for species isp c check(2,isp) = concentration below zero for species isp c check(3,isp) = (-) total mass added to the system to get a positive conc. do j = 1,3 do isp = 1,naers check(j,isp) = 0.0 enddo enddo do isp = 1,naers do isect=1,nsect check(1,isp) = check(1,isp) + aerosol(isect,isp) enddo enddo c If negative in a section, change to a small positive value c Account for this mass change by subtracting from other c non-negative sections. do isp = 1,naers do isect = 1,nsect if (aerosol(isect,isp) .lt. 0.0) then check(2,isp) = check(2,isp)+aerosol(isect,isp) check(3,isp) = check(3,isp)+aerosol(isect,isp)-1.e-12 ! bkoo (10/07/03) aerosol(isect,isp) = 1.e-12 endif enddo enddo do isect = 1,nsect if (aerosol(isect,nan) .gt. 1.e-12) then aerosol(isect,nan)=aerosol(isect,nan)+(aerosol(isect,nan))* & check(3,nan)/(check(1,nan)-check(2,nan)) endif if (aerosol(isect,nac) .gt. 1.e-12) then aerosol(isect,nac)=aerosol(isect,nac)+(aerosol(isect,nac))* & check(3,nac)/(check(1,nac)-check(2,nac)) endif if (aerosol(isect,naa) .gt. 1.e-12) then aerosol(isect,naa)=aerosol(isect,naa)+(aerosol(isect,naa))* & check(3,naa)/(check(1,naa)-check(2,naa)) endif c For HSO5 and HMSA, if negative and there is not enough mass in the c remaining aerosol sections to account for the added mass, subtract c from sulfate. if (check(1,nahso5)-check(2,nahso5) .lt. 0.-check(3,nahso5)) then check(3,na4) = check(3,na4)+check(3,nahso5)*96./113. else if (aerosol(isect,nahso5) .gt. 1.e-12) then aerosol(isect,nahso5)=aerosol(isect,nahso5)+ & (aerosol(isect,nahso5))*check(3,nahso5)/ & (check(1,nahso5)-check(2,nahso5)) endif if (check(1,nahmsa)-check(2,nahmsa) .lt. 0.-check(3,nahmsa)) then check(3,na4) = check(3,na4)+check(3,nahmsa)*96./111. else if (aerosol(isect,nahmsa) .gt. 1.e-12) then aerosol(isect,nahmsa)=aerosol(isect,nahmsa)+ & (aerosol(isect,nahmsa))*check(3,nahmsa)/ & (check(1,nahmsa)-check(2,nahmsa)) endif if (aerosol(isect,na4) .gt. 1.e-12) then aerosol(isect,na4)=aerosol(isect,na4)+(aerosol(isect,na4))* & check(3,na4)/(check(1,na4)-check(2,na4)) endif enddo c c RETURN REST OF YAQ RESULTS TO THEIR ORIGINAL MATRICES (GAS,AEROSOL) c c ** GAS-PHASE SPECIES ** c IMPORTANT : THE DISSOLVED S(IV) and H2O2 ARE TRANSFERRED c BACK TO THE GAS PHASE c gas(nghcho) = yaq(1)/fhcho ! total HCHO (ppm) gas(ngso2)= (yaq(2)+yaq(6)*0.7901+yaq(14)*0.7901)/fso2 ! SO2(g) (ppm) gas(ngh2o2)=(yaq(3)+yaq(7)+yaq(15))/fh2o2 ! H2O2(g) (ppm) gas(nga)=yaq(4)/fammon ! NH3(g) (ppm) gas(nghcooh)= yaq(5)/fhcooh ! total HCOOH (ppm) c c THE CODE NEGLECTS THE CHANGE IN THE SIZE DISTRIBUTION SHAPE c OF THE NONVOLATILE AEROSOL COMPONENTS BECAUSE OF THE SULFATE c PRODUCTION. THE CHANGE IN THE SULFATE DISTRIBUTION IS CALCULATED c USING THE DISTRIBUTION FUNCTIONS WHILE THE CHANGE IN THE c DISTRIBUTION OF THE VOLATILE INORGANIC AEROSOL COMPONENTS WILL c CALCULATED BY THE AEROSOL MODULE AFTER CLOUD EVAPORATION c c c.....END OF CODE return end