C======================================================================= C C *** MADM CODE C *** SUBROUTINE EQUAER C *** THIS SUBROUTINE CALLS ISORROPIA+ FOR CALCULATING THE AEROSOL C EQUILIBRIUM CONCENTRATIONS C C ======================== ARGUMENTS / USAGE =========================== C C INPUT: c IPROB = 1=Reverse prob, 0=Foreward prob C 1. [WI] C DOUBLE PRECISION array of length [NCOMP=5]. [NCOMP] is defined in C include file 'isrpia.inc'. C Total aerosol concentrations, expressed in moles/m3. C WI(1) - sodium, expressed as [Na] C WI(2) - sulfate, expressed as [H2SO4] C WI(3) - ammonium, expressed as [NH3] C WI(4) - nitrate, expressed as [HNO3] C WI(5) - chloride, expressed as [HCl] C C 2. [RHI] C DOUBLE PRECISION variable. C Ambient relative humidity expressed on a (0,1) scale. C C 3. [TEMPI] C DOUBLE PRECISION variable. C Ambient temperature expressed in Kelvins. C C OUTPUT: C 1. [AEROSI] C DOUBLE PRECISION array of length [NINTI=21]. C Aerosol species concentrations, expressed in ug/m3. C c(K,01) - H2O c(K,12) - (NH4)2SO4(s) C c(K,02) - H+(aq) c(K,13) - NH4NO3(s) C c(K,03) - Na+(aq) c(K,14) - NH4Cl(s) C c(K,04) - NH4+(aq) c(K,15) - H2SO4(aq) (Undefined) C c(K,05) - Cl-(aq) c(K,16) - NH4HSO4(s) C c(K,06) - SO4--(aq) c(K,17) - NaHSO4(s) C c(K,07) - HSO4-(aq) c(K,18) - (NH4)4H(SO4)2(s) C c(K,08) - NO3-(aq) c(K,19( - OH-(aq) C c(K,09) - NaCl(s) C c(K,10) - Na2SO4(s) C c(K,11) - NaNO3(s) C C 2. [ps] C DOUBLE PRECISION array of length [NGAS]. C Gaseous species partial pressures, expressed in Pascals. C ps(k,1) - NH3 C ps(k,2) - HNO3 C ps(k,3) - H2SO4 C ps(k,4) - HCl C C 3. [DRYF] C LOGICAL variable. C Contains information about the physical state of the system. C .TRUE. - There is no aqueous phase present C .FALSE.- There is an aqueous phase present C C======================================================================= C SUBROUTINE EQUAER (k, q, iprob) C INCLUDE 'dynamic.inc' ! MADM declarations INCLUDE 'equaer.inc' ! ISORROPIA declarations C real*8 q(ntotal) ! local conc. array DIMENSION WI(NCOMP) dimension pgas(ngas) integer k ! master array section indices integer nsec0 ! last section of master array NOT in q integer nsec1 ! last section of master array occuring in q C C *** CONVERT INPUT CONCENTRATIONS TO moles/m3 ************************** C nn = (k-1)*nsp ! aerosol ng = nsp*nsec ! gases if(q(nn+KNa).le.tinys.and.q(nn+KSO4).le.tinys.and. # q(nn+KNH4).le.tinys.and.q(nn+KNO3).le.tinys.and. # q(nn+KCL).le.tinys) go to 999 100 WI(1) = q(nn+KNa) /emw(KNa) *1.D-6 WI(2) = q(nn+KSO4)/emw(KSO4)*1.D-6 WI(3) = q(nn+KNH4)/emw(KNH4)*1.D-6 WI(4) = q(nn+KNO3)/emw(KNO3)*1.D-6 WI(5) = q(nn+KCL) /emw(KCL) *1.D-6 c iprob is always 1; the codes below are never used - bkoo if(iprob.eq.0) then ! add gas phase prod=rgas*temp/(1.01325D5*pres) ! conversion from ppm to umoles/m3 ! pres (bkoo, 06/09/00) WI(1) = WI(1) + 0.0D0 WI(2) = WI(2) + Q(ng+ih2SO4)/PROD*1.0d-6 WI(3) = WI(3) + q(ng+iNH3) /PROD*1.0d-6 WI(4) = WI(4) + q(ng+ihNO3) /PROD*1.0d-6 WI(5) = WI(5) + q(ng+ihCL) /PROD*1.0d-6 endif C C *** CALL ISORROPIA+ *************************************************** C RHI=rh TEMPI=temp c IPROB = 1 ! 1=Reverse prob, 0=Foreward prob metstbl=ims(k) ! bkoo (03/05/02) CALL ISRPIA ( WI, RHI, TEMPI, IPROB ) metstbl=0 ! bkoo (03/05/02) C C *** CONVERT OUTPUT CONCENTRATIONS TO ug/m3 **************************** C C AEROSOL SPECIES C q(nn+KH2O)= WATER *1.D9 ! H2O C(k,01) = q(nn+KH2O) ! H2O C(k,02) = MOLAL(01) *1.D6 ! H+(aq) C(k,03) = MOLAL(02)*intmw(03)*1.D6 ! Na+(aq) C(k,04) = MOLAL(03)*intmw(04)*1.D6 ! NH4+(aq) C(k,05) = MOLAL(04)*intmw(05)*1.D6 ! Cl-(aq) C(k,06) = MOLAL(05)*intmw(06)*1.D6 ! SO4--(aq) C(k,07) = MOLAL(06)*intmw(07)*1.D6 ! HSO4-(aq) C(k,08) = MOLAL(07)*intmw(08)*1.D6 ! NO3-(aq) C(k,09) = CNACL *intmw(09)*1.D6 ! NaCl(s) C(k,10) = CNA2SO4 *intmw(10)*1.D6 ! Na2SO4(s) C(k,11) = CNANO3 *intmw(11)*1.D6 ! NaNO3(s) C(k,12) = CNH42S4 *intmw(12)*1.D6 ! (NH4)2SO4(s) C(k,13) = CNH4NO3 *intmw(13)*1.D6 ! NH4NO3(s) C(k,14) = CNH4CL *intmw(14)*1.D6 ! NH4Cl(s) C(k,15) = 0.0D0 ! H2SO4(aq) (not used by isorropia) C(k,16) = CNH4HS4 *intmw(16)*1.D6 ! NH4HSO4(s) C(k,17) = CNAHSO4 *intmw(17)*1.D6 ! NaHSO4(s) C(k,18) = CLC *intmw(18)*1.D6 ! (NH4)4H(SO4)2(s) C(k,19) = COH *intmw(19)*1.D6 ! bkoo (02/14/02) c c Since issorropia ignores no3 and cl concentrations if the particle is c found to be both acidic and dry, we need to evaporate these species c instantaneously. Since NH4Cl forms from H2SO4 condensation before c NH4NO3, we start transfer cloride before nitrate. c We then call issorropia again to calculate the new ci's. c c Since statements except [if] and [write] have been blocked, however, c nothing actually happens. Therefore the whole part was made blocked. c (if unblocked, it should be modified so that pressure terms are included) c - bkoo c c if(dryf) then c if(c(k,18).gt.tinys.or.c(k,17).gt.tinys.or.c(k,16).gt.tinys) then c irecalc=0 c acid=clc*1.d6+cnahso4*1.d6+cnh4hs4*1.d6 c if(q(nn+kcl).gt.tinys) then c if(q(nn+kcl).lt.acid*emw(kcl)) then c write(6,*) tcom,'SOLID ACID ',k,q(nn+kcl), c & 'ug/m3 cl- moved to gas1' c write(90,*) tcom,'SOLID ACID ',k,q(nn+kcl), c & 'ug/m3 cl- moved to gas1' c q(ng+ihcl)=q(ng+ihcl)+ c & q(nn+kcl)*rgas*temp/(1.015d5*gmw(ihcl)) c q(nn+kcl)=0.0d0 c acid=acid-q(nn+kcl)/float(emw(kcl)) c else c write(6,*) tcom,'SOLID ACID ',k,acid*emw(kcl), c & 'ug/m3 cl- moved to gas2' c write(90,*) tcom,'SOLID ACID ',k,acid*emw(kcl), c & 'ug/m3 cl- moved to gas2' c q(ng+ihcl)=q(ng+ihcl)+ c & acid*emw(kcl)*rgas*temp/(1.015d5*gmw(ihcl)) c q(nn+kcl)=q(nn+kcl)-acid*emw(kcl) c acid=0.0d0 c endif c irecalc=1 c endif c if(q(nn+kno3).gt.tinys.and.acid.gt.0.0d0) then c if(q(nn+kno3).lt.acid*emw(kno3)) then c write(6,*) tcom,'SOLID ACID ',k,q(nn+kno3), c & 'ug/m3 NO3- moved to gas1' c write(90,*) tcom,'SOLID ACID ',k,q(nn+kno3), c & 'ug/m3 NO3- moved to gas1' c q(ng+ihno3)=q(ng+ihno3)+ c & q(nn+kno3)*rgas*temp/(1.015d5*gmw(ihno3)) c q(nn+kno3)=0.0d0 c else c write(6,*) tcom,'SOLID ACID ',k,acid*emw(kno3), c & 'ug/m3 NO3- moved to gas2' c write(90,*) tcom,'SOLID ACID ',k,acid*emw(kno3), c & 'ug/m3 NO3- moved to gas2' c q(ng+ihno3)=q(ng+ihno3)+ c & acid*emw(kno3)*rgas*temp/(1.015d5*gmw(ihno3)) c q(nn+kno3)=q(nn+kno3)-acid*emw(kno3) c endif c irecalc=1 c endif c if(irecalc.eq.1) goto 100 c endif c endif C C GASEOUS SPECIES C c Originally Pilinis tried to slow things down by capping the "ps" values c at 1.0. I had to add the acidity constraint to achieve control in a more c robust way. Because the time scales of concern for atmospheric modeling c are large enough, artificially slowing things down to the degree in HYBR c seemed appropriate to improve stability. There is nothing special however c about a "ps" of 1.0 though. It just seemed large enough not to seriously c affect accuracy and small enough to allow the code to run stably in a c reasonable amount of time. For MADM I expect I took this out because c I was more concerned about accuracy and less concerned about runtime. c Also, MADM is somewhat more stable since there is no division between c "equilibrium" and "dynamic" particles. - Kevin, 5/30/00 c (bkoo, 06/06/00) c if(aerm.eq.'MADM') then ! bkoo (01/08/01) ps(k,INH3 ) = GNH3 *rgas*temp ! NH3(g) (Pa) ps(k,IHNO3 ) = GHNO3*rgas*temp ! HNO3(g) (Pa) ps(k,IH2SO4) = 0.0D0 ! H2SO4(g) (Pa) ps(k,IHCL ) = GHCL *rgas*temp ! HCL(g) (Pa) c else c ps(k,INH3 ) = min(1.0d0,GNH3 *rgas*temp) c ps(k,IHNO3 ) = min(1.0d0,GHNO3*rgas*temp) c ps(k,IH2SO4) = 0.0D0 c ps(k,IHCL ) = min(1.0d0,GHCL *rgas*temp) c endif C C 'DRY' FLAG C dry(k) = DRYF c RETURN 999 CONTINUE q(nn+KH2O)= 0.0 ! H2O C(k,01) = 0.0 ! H2O C(k,02) = 0.0D0 ! H+(aq) C(k,03) = 0.0D0 ! Na+(aq) C(k,04) = 0.0D0 ! NH4+(aq) C(k,05) = 0.0D0 ! Cl-(aq) C(k,06) = 0.0D0 ! SO4--(aq) C(k,07) = 0.0D0 ! HSO4-(aq) C(k,08) = 0.0D0 ! NO3-(aq) C(k,09) = 0.0D0 ! NaCl(s) C(k,10) = 0.0D0 ! Na2SO4(s) C(k,11) = 0.0D0 ! NaNO3(s) C(k,12) = 0.0D0 ! (NH4)2SO4(s) C(k,13) = 0.0D0 ! NH4NO3(s) C(k,14) = 0.0D0 ! NH4Cl(s) C(k,15) = 0.0D0 ! H2SO4(aq) (not used by isorropia) C(k,16) = 0.0D0 ! NH4HSO4(s) C(k,17) = 0.0D0 ! NaHSO4(s) C(k,18) = 0.0D0 ! (NH4)4H(SO4)2(s) C(k,19) = 0.0D0 ! bkoo (02/14/02) ps(k,INH3 ) = 0.0D0 ! NH3(g) (Pa) ps(k,IHNO3 ) = 0.0D0 ! HNO3(g) (Pa) ps(k,IH2SO4) = 0.0D0 ! H2SO4(g) (Pa) ps(k,IHCL ) = 0.0D0 ! HCL(g) (Pa) RETURN C C *** END OF SUBROUTINE EQUAER ****************************************** C END