! #include "tm5.inc" ! module eqsam_param implicit none private public :: eqsam_v03d contains subroutine eqsam_v03d(yi,yo,nca,nco,iopt,loop,imax,ipunit,in) ! implicit none !___________________________________________________________________________________________________________________________________ ! Written by Swen Metzger 3/11/99. Modified 2002, 2003. ! ! Department of Atmospheric Chemistry, Max-Planck-Institute for Chemistry. ! email: metzger@mpch-mainz.mpg.de ! http://www.mpch-mainz.mpg.de/~metzger ! ! COPYRIGHT 1999-2003 ! ! purpose ! ------- ! EQSAM (EQuilibrium Simplified Aerosol Model) is a new and computationally efficient thermodynamic ! aerosol composition model that allows to calculate the gas/aerosol equilibrium partitioning, ! including aerosol water, sufficiently fast and accurate for global (or even regional) modeling. ! EQSAM is based on a number of parameterizations, including single solute molalities and activity ! coefficients (AC). The thermodynamic framework (domains and subdomains, internally mixed aerosols) ! is the same as of more sophisticated thermodynamic equilibrium models (EQMs), e.g. of ISORROPIA ! (Nenes et al., 1998). Details are given in the references below (and the references therein). ! ! The main assumption on which EQSAM/EQMs are based is thermodynamical and chemical equilibrium. ! From this assumption it directly follows that the aerosol water activity (aw) equals the ambient ! relative humidity (RH), if the water vapor pressure is sufficiently larger than the partial vapor ! pressure of the aerosol compounds. This is approximately true for tropospheric aerosols. Given the ! large amount of water vapor present, water vapor and aerosol water equilibrate relatively faster ! compared to all other aerosol compounds. This is subsequently also true for single aerosol compounds. ! The water activity of single solutes must also equal RH under this assumption. Therefore, the so ! called ZSR-relation is (and can be) used to calculate the aerosol associated water mass (simply ! from the sum of all water mass fractions that are derived from measured single solute molalities). ! ! In contrast to other EQMs, EQSAM utilizes the fact that the RH fixes the water activity ! (under the above assumptions) and the consequence that any changes in RH also causes changes in ! the aerosol water mass and, hence, aerosol activity (including activity coefficients). Thus, an decrease ! (increase) in RH decrease (increases) the aerosol water mass (and water activity). This can change the ! aerosol composition, e.g. due to condensation (evaporation/crystallization), because the vapor pressure ! above the aerosol reduces (increases). In turn, a vapor pressure reduction (increase) due to changes ! in the aerosol composition is compensated by an associated condensation (evaporation) of water vapor ! to maintain the aerosol molality to remain constant (because aw=RH). Furthermore, the aerosol water ! mainly depends on the aerosol mass and the type of solute, so that parameterizations of single solute ! molalities and activity coefficients can be defined, only depending on the type of solute and RH. ! The advantage of using such parameterizations is that the entire aerosol equilibrium composition ! can be solved analytically, i.e. non-iteratively, which considerably reduces the amount of CPU time ! that is usually need for aerosol thermodynamic calculations (especially if an EQM is incorporated in ! an aerosol dynamical model that is in turn embedded in a high resolution regional or global model). ! ! However, EQSAM should still be regarded as a starting point for further developments. There is still ! room for improvements. For instance, this code is not yet numerically optimized (vectorized) and a ! number of improvements with respect to an explicit treatment of additional equilibrium reactions, ! missing (or only implicit) dissociation, and a basic parameterization of the water uptake. ! ! Note that EQSAM was originally developed to calculate the gas/aerosol equilibrium partitioning of the ! ammonium-sulfate-nitrate-water system for climate models, excluding solid compounds. ! This version (eqsam_v03d.f90) is extended with respect to sea salt. Solids/hysteresis are treated in a ! simplified manner. Results of a box model comparison with ISORROPIA will be available from the web page. ! Please also note that the water uptake is based on additional (unpublished) parameterizations for single ! solute molalities, which are derived from tabulated measurements used in ISORROPIA. Note further that ! this extended version (eqsam_v03d.f90) is not yet published. A publication is in progress. ! ! ToDo: ! Split ion-pairs into ions for water parameterizations (since info is actually available) ! Include uptake/dissociation of NH3, HNO3, HCl (mainly to get pH right at near neutral conditions) ! Extension to K+,Ca++,Mg++, CO2/(CO3)2--/HCO3-,SOA,etc.. (maybe not) ! Vectorization. Translation of hardcoded formulas in array syntax. ! I/O Interface and program structure clean up. ! EQSAM info webpage. ! ! Version History: ! ! eqsam_v03d.f90 (MPI-CH, June 2003): ! - gama parameterizations now according to Metzger 2002 (JGR Appendix) ! - improved pH calculations (still restricted to strong acids) ! - removed bug that lead to too high nitrate formation at dry and cold regions (UT/LS) ! - removed bug in solid/hysteresis calculations ! (both bugs introduced in eqsam_v03b.f90 by cleaning up eqsam_v02a.f90) ! ! eqsam_v03c.f90 (MPI-CH, April 2003): ! - more accurate paramterizations of single solute molalities (Na, Cl species) ! - cleanded up RHD subdomain structure ! - improved water uptake (Na, Cl species) ! ! eqsam_v03b.f90 (MPI-CH, March 2003): ! System extended to HCl,Cl-/Na+. ! Parameterization (fit) of additional HNO3 uptake removed. ! Instead, complete analytical solution of equilibrium reactions, based on the AC-RH relationship. ! eqsam_v03.f90 (IMAU, October 1999): ! Test version (included in TM3). ! eqsam_v02a.f90 (IMAU, April 2000): ! Box model version. ! eqsam_v02.f90 (IMAU, October 1999): ! TM3 version. ! Version including solids and additional HNO3 uptake on acidic aerosols (parameterized). ! eqsam_v01b.f90 (MPI-CH, January 2003): ! Same as eqsam_v01a.f90 (additional lines though uncommented for test purposes only). ! eqsam_v01a.f90 (IMAU, April 2000): ! Box model version. ! eqsam_v01.f90 (IMAU, October 1999): ! TM3 version. ! First and most basic version (without solids) for better vectorization (for global modeling). ! System: NH3,NH4+/H2SO4+,HSO4-,SO4--/HNO3,NO3-, H2O ! based on equilibrium / internal mixture assumption / aw=rh / ZSR-relation ! parameterization of activcity coefficients (AC), i.e. an AC-RH relationship ! ! ! interface ! --------- ! call eqsam_v03d(yi,yo,nca,nco,iopt,loop,imax,ipunit,in) ! ! yi = input array (imax, nca) ! yo = output array (imax, nco) ! imax = max loop (e.g. time steps) ! nca >= 11 ! nco >= 35 ! iopt = 1 metastable ! iopt = 2 solids ! iopt = 3 hysteresis (metastable/solids) for online calculations ! iopt = 31 hysteresis lower branch ! iopt = 32 hysteresis upper branch ! ipunit = I/O unit (can be skipped) ! in = array (can be skipped) ! ! method ! ------ ! equilibrium / internal mixture assumption / aw=rh ! System: NH3,NH4+/H2SO4+,HSO4-,SO4--/HNO3,NO3-, HCl,Cl-/Na+, H2O ! (K+,Ca++,Mg++) ! external ! -------- ! program eqmd.f90 (driver only needed for the box model version) ! subroutine gribio.f90 (provides diagnostics output in grib/binary/ascii format) ! ! references ! --------- ! Swen Metzger Ph.D Thesis, University Utrecht, 2000. ! http://www.library.uu.nl/digiarchief/dip/diss/1930853/inhoud.htm ! ! Metzger, S. M., F. J. Dentener, J. Lelieveld, and S. N. Pandis, ! GAS/AEROSOL PARTITIONING I: A COMPUTATIONALLY EFFICIENT MODEL, ! J Geophys. Res., 107, D16, 10.1029/2001JD001102, 2002 ! http://www.agu.org/journals/jd/jd0216/2001JD001102/index.html ! Metzger, S. M., F. J. Dentener, A. Jeuken, and M. Krol, J. Lelieveld, ! GAS/AEROSOL PARTITIONING II: GLOBAL MODELING RESULTS, ! J Geophys. Res., 107, D16, 10.1029/2001JD001103, 2002. ! http://www.agu.org/journals/jd/jd0216/2001JD001103/index.html !___________________________________________________________________________________________________________________________________ real,parameter :: RH_HIST_DW=1.50 ! mean value for mixture of wet (2) and dry (1) gridboxes (needed for HYSTERESIS) real,parameter :: T0=298.15,T1=298.0,AVO=6.03e23,R=82.0567e-6, & ! in cu.m*atm/deg/mole r_kcal = 1.986E-3 ! Ideal gas constant [kcal K-1.mole-1] real,parameter :: RHMAX=0.99,RHMIN=0.0001 ! restrict to max / min RH real,parameter :: MWNH4=18.,MWSO4=96.,MWNO3=62.,MWCl=35.5 ! mole mass of species considered real,parameter :: MWNa=23.0,MWCa=40.1,MWN=14.0, MWS=32.1 real,parameter :: MWH20=55.51*18.01,ZERO=0.0 real,parameter :: GF1=0.25,GF2=0.50,GF3=0.40,GF4=1.00,K=2. ! exponents of AC-RH functions !______________________________________________ integer,parameter :: NPAIR=10 ! integer :: ii,il,IHYST integer,intent(in) :: nca,nco,imax,loop,ipunit integer,intent(inout) :: iopt !______________________________________________ integer,dimension(6),intent(in) :: in !______________________________________________ real :: T0T,TT,RH,PX,RHD,KAN,KAC,ZIONIC,RH_HIST,GAMA,GG,GF,GFN real :: X00,X01,X02,X03,X04,X05,X08,X09,X10,X11 real :: X0,X1,X2,X3,X4,X5,X6,XK10,XK6 real :: ZFLAG,ZKAN,ZKAC,PH,COEF,HPLUS,AKW,XKW,MOLAL real :: TNH4,TSO4,TNO3,TNa,TCl,TPo,TCa,TMg real :: PNH4,PSO4,PNO3,PCl,PNa,GNO3,GNH3,GSO4,GHCl real :: ASO4,ANO3,ANH4,ACl,ANa,SNH4,SSO4,SNO3,SCl,SNa real :: WH2O,PM,PMs,PMt,RINC,DON,RATIONS,GR,NO3P,NH4P Real :: H2O_NO3 !_______________________________________________ real,dimension(imax,nca),intent(in) :: yi real,dimension(imax,nco),intent(out) :: yo real,dimension(8) :: w1,w2 real,dimension(8) :: RHDA,RHDE,RHDX,RHDZ ! RHD / MRHD arrays for different aerosol types real,dimension(NPAIR) :: M0,MW,NW,ZW ! arrays of ion pairs ! ! salt solutes: ! 1 = NACl, 2 = (NA)2SO4, 3 = NANO3, 4 = (NH4)2SO4, 5 = NH4NO3, 6 = NH4CL, 7 = 2H-SO4 ! 8 = NH4HSO4, 9 = NAHSO4, 10 = (NH4)3H(SO4)2 ! DATA MW(1:NPAIR)/ 58.5, 142.0, 88.0, 132.0, 80.0, 53.5, 98.0, 115.0, 120.0, 247.0/ ! mole mass of the salt solute DATA NW(1:NPAIR)/ 2.0, 2.5, 2.5, 2.5, 3.5, 1.0, 4.5, 2.0, 2.0, 2.5/ ! square of max. dissocation number (not consistent) DATA ZW(1:NPAIR)/ 0.67, 1.0, 1.0, 1.0, 1.0, 1.0, 0.5, 1.0, 1.0, 1.0/ ! exponents of water activity functions ! DATA RHDA(1:8)/0.32840, 0.4906, 0.6183, 0.7997, 0.67500, 0.5000, 0.4000, 0.0000/ ! RHD / MRHD values as of ISORROPIA / SCAPE (T=298.15K) DATA RHDE(1:8)/-1860.0, -431.0, 852.00, 80.000, 262.000, 3951.0, 384.00, 0.0000/ ! Temp. coeff. !___________________________________________________________________________________________________________________________________ IHYST=2 IF(IOPT.EQ.31) THEN ! SOLID HYSTORY IHYST=1 IOPT=3 ELSEIF(IOPT.EQ.32) THEN ! WET HISTORY IHYST=2 IOPT=3 ENDIF !kt write(ipunit,*)'eqsam_v03d ...' !kt print*,' ' !kt print*,' EQuilibrium Simplified Aerosol Model (EQSAM)' !kt print*,' for global modeling ' !kt print*,' by ' !kt print*,' Swen Metzger, MPI-CH ' !kt print*,' Copyright 1999-2003 ' !kt print*,' >> metzger@mpch-mainz.mpg.de << ' !kt print*,' last change: 04. June, 2003 ' !kt print*,' (version 3.0d) ' !kt print*,' gas/aerosol calculations assuming ' !kt print*,' System: NH3,NH4+/H2SO4+,HSO4-,SO4-- ' !kt print*,' HNO3,NO3-, HCl,Cl-/Na+, H2O ' !kt if(iopt.eq.1) then !kt print*,' metastable aeorsols ' !kt elseif(iopt.eq.2) then !kt print*,' solid aeorsols ' !kt elseif(iopt.eq.3) then !kt print*,' hysteresis ' !kt print*,' (metastable/solids) ' !kt if(IHYST.eq.1) then !kt print*,' solid hystory ' !kt elseif(IHYST.eq.2) then !kt print*,' wet hystory ' !kt endif !kt endif !kt print*,' ' !kt print*,'loop over ',loop,' data sets' !kt print*,' ' !___________________________________________________________________________________________________________________________________ yo=0.;w1=0.;w2=0. ! init/reset !___________________________________________________________________________________________________________________________________ do il=1,loop ! get old relative humidity to calculate aerosol hysteresis (online only) RH_HIST = 2. ! WET HISTORY (DEFAULT) IF(IHYST.EQ.1.OR.IOPT.EQ.2) RH_HIST = 1. ! SET TO SOLIDS ! meteorology TT = yi(il,1) ! T [K] RH = yi(il,2) ! RH [0-1] PX = yi(il,11) ! p [hPa] ! ! gas+aerosol: w1(1) = yi(il,6) ! Na+ (ss + xsod) (a) [umol/m^3 air] w1(2) = yi(il,4) ! H2SO4 + SO4-- (p) [umol/m^3 air] w1(3) = yi(il,3) ! NH3 (g) + NH4+ (p) [umol/m^3 air] w1(4) = yi(il,5) ! HNO3 (g) + NO3- (p) [umol/m^3 air] w1(5) = yi(il,7) ! HCl (g) + Cl- (p) [umol/m^3 air] w1(6) = yi(il, 8) ! K+ (p) from Dust [umol/m^3 air] w1(7) = yi(il, 9) ! Ca++ (p) from Dust [umol/m^3 air] w1(8) = yi(il,10) ! Mg++ (p) from Dust [umol/m^3 air] !______________________________________________ zflag=1. w1=w1*1.0e-6 ! [mol/m^3 air] TNa = w1(1) ! total input sodium (g+p) TSO4 = w1(2) ! total input sulfate (g+p) TNH4 = w1(3) ! total input ammonium (g+p) TNO3 = w1(4) ! total input nitrate (g+p) TCl = w1(5) ! total input chloride (g+p) TPo = w1(6) ! total input potasium (g+p) TCa = w1(7) ! total input calcium (g+p) TMg = w1(8) ! total input magnesium(g+p) ! SULFATE RICH if((w1(1)+w1(3)+w1(6)+2.*(w1(7)+w1(8))).le.(2.*w1(2))) then zflag=3. endif ! SULFATE VERY RICH CASE if (NH4+Na+K+2(Ca+Mg))/SO4 < 1 if((w1(1)+w1(3)+w1(6)+2.*(w1(7)+w1(8))).le.w1(2)) then zflag=4. endif ! SULFATE NEUTRAL CASE if((w1(1)+w1(3)+w1(6)+2.*(w1(7)+w1(8))).gt.(2.*w1(2))) then zflag=2. endif ! SULFATE POOR AND CATION POOR CASE if((w1(1)+w1(6)+2.*(w1(7)+w1(8))).gt.(2.*w1(2))) then zflag=1. endif IF ( RH .LT. RHMIN ) RH=RHMIN IF ( RH .GT. RHMAX ) RH=RHMAX ! CALCULATE TEMPERATURE DEPENDENCY FOR SOME RHDs RHDX(:)=RHDA(:)*exp(RHDE(:)*(1./TT-1./T0)) RHDZ(:)=RHDX(:) ! ACCOUNT FOR VARIOUS AMMOMIUM/SODIUM SULFATE SALTS ACCORDING TO MEAN VALUE AS OF ISORROPIA GG=2.0 ! (Na)2SO4 / (NH4)2SO4 IS THE PREFFERED SPECIES FOR SULFATE DEFICIENT CASES IF(ZFLAG.EQ.3.) THEN IF(RH.LE.RHDZ(7)) THEN ! ACCOUNT FOR MIXTURE OF (NH4)2SO4(s) & NH4HSO4(s) & (NH4)3H(SO4)2(s) GG=1.677 ! (Na)2SO4 & NaHSO4 ! GG=1.5 ELSEIF(RH.GT.RHDZ(7).AND.RH.LE.RHDZ(5)) THEN ! MAINLY (Na)2SO4 / (NH4)2SO4(s) & (NH4)3H(SO4)2(s) GG=1.75 ! GG=1.5 ELSEIF(RH.GE.RHDZ(5)) THEN ! (NH4)2SO4(S) & NH4HSO4(S) & SO4-- & HSO4- GG=1.5 ! (Na)2SO4 & NaHSO4 ENDIF ENDIF IF(ZFLAG.EQ.4.) GG=1.0 ! IF SO4 NEUTRALIZED, THEN ONLY AS NaHSO4 / NH4HSO4(S) OR HSO4- / H2SO4 RHD=RH IF(IOPT.EQ.2.OR.RH_HIST.LT.RH_HIST_DW) THEN ! GET RHD FOR SOLIDS / HYSTERESIS ! ! GET LOWEST DELIQUESCENCE RELATIVE HUMIDITIES ACCORDING TO THE CONCENTRATION DOMAIN (APROXIMATION) ! BASED ON RHD / MRHD ISORROPIA/SCAPE ! w2(:)=1. do ii=1,8 if(w1(ii).le.1.e-12) w2(ii)=0. ! skip compound in RHD calculation if value is concentration is zero or rather small enddo ! GET LOWEST RHD ACCORDING TO THE CONCENTRATION DOMAIN ! zflag=1. (cation rich) ... ! 1. sea salt aerosol : RHDX(1)=MgCl2 ! 2. mineral dust aerosol : RHDX(2)=Ca(NO3)2 ! ! zflag=2. (sulfate neutral) ... ! 3. ammonium + nitrate : RHDX(3)= NH4NO3 ! 4. ammonium + sulfate : RHDX(4)=(NH4)2SO4 ! 5. ammonium + sulfate mixed salt : RHDX(5)=(NH4)3H(SO4)2, (NH4)2SO4 ! 6. ammonium + nitrate + sulfate : RHDX(6)=(NH4)2SO4, NH4NO3, NA2SO4, NH4CL ! ! zflag=3. (sulfate poor) ... ! 7. ammonium + sulfate (1:1,1.5) : RHDX(7)= NH4HSO4 ! ! zflag=4. (sulfate very poor) ... ! 8. sulfuric acid : RHDX(8)= H2SO4 IF(ZFLAG.EQ.1.)THEN RHD=W2(1)+W2(5) ! Na+ dependency IF(RHD.EQ.0.) RHDX(1)=1. RHD=W2(6)+W2(7)+W2(8) ! K+/Ca++/Mg++ dependency (incl. ss) IF(RHD.EQ.0.) RHDX(2)=1. RHD=MINVAL(RHDX(1:2)) ELSEIF(ZFLAG.EQ.2.)THEN RHD=W2(3)*W2(4) ! NH4+ & NO3- dependency IF(RHD.EQ.0.) RHDX(3)=1. RHD=W2(2)+W2(3) ! NH4+ & SO4-- dependency IF(GG.NE.2.) RHD=0. ! account only for pure (NH4)2SO4 IF(RHD.EQ.0.) RHDX(4)=1. RHD=W2(2)+W2(3) ! NH4+ & SO4-- dependency IF(RHD.EQ.0.) RHDX(5)=1. RHD=W2(2)+W2(3)+W2(4)+W2(5) ! (NH4)2SO4, NH4NO3, NA2SO4, NH4CL dependency IF(RHD.EQ.0.) RHDX(6)=1. ! RHD=MINVAL(RHDX(3:4)) RHD=MINVAL(RHDX(3:6)) ELSEIF(ZFLAG.EQ.3.)THEN RHD=W2(2)+W2(3) ! NH4+ & SO4-- dependency IF(RHD.EQ.0.) RHDX(7)=1. RHD=RHDX(7) ELSEIF(ZFLAG.EQ.4.)THEN RHD=W2(2) ! H2SO4 dependency (assume no dry aerosol) IF(RHD.EQ.0.) RHDX(8)=1. RHD=RHDX(8) ENDIF ! ZFLAG ENDIF ! SOLIDS ! GET WATER ACTIVITIES ACCORDING TO METZGER, 2000. ! FUNCTION DERIVED FROM ZSR RELATIONSHIP DATA (AS USED IN ISORROPIA) M0(:) = ((NW(:)*MWH20/MW(:)*(1./RH-1.)))**ZW(:) ! CALCULATE TEMPERATURE DEPENDENT EQUILIBRIUM CONSTANTS T0T=T0/TT COEF=1.0+LOG(T0T)-T0T ! EQUILIBRIUM CONSTANT NH4NO3(s) <==> NH3(g) + HNO3(g) [atm^2] (ISORROPIA) XK10 = 5.746e-17 XK10= XK10 * EXP(-74.38*(T0T-1.0) + 6.120*COEF) KAN = XK10/(R*TT)/(R*TT) ! EQUILIBRIUM CONSTANT NH4CL(s) <==> NH3(g) + HCL(g) [atm^2] (ISORROPIA) XK6 = 1.086e-16 XK6 = XK6 * EXP(-71.00*(T0T-1.0) + 2.400*COEF) KAC = XK6/(R*TT)/(R*TT) ! ! CALCULATE AUTODISSOCIATION CONSTANT (KW) FOR WATER H2O <==> H(aq) + OH(aq) [mol^2/kg^2] (ISORROPIA) XKW = 1.010e-14 XKW = XKW *EXP(-22.52*(T0T-1.0) + 26.920*COEF) ! GET MEAN MOLAL IONIC ACTIVITY COEFF ACCORDING TO METZGER, 2002. GAMA=0.0 IF(RH.GE.RHD) GAMA=(RH**ZFLAG/(1000./ZFLAG*(1.-RH)+ZFLAG)) GAMA = GAMA**GF1 ! ONLY GAMA TYPE OF NH4NO3, NaCl, etc. NEEDED SO FAR GAMA=0.0 GFN=K*K ! K=2, i.e. condensation of 2 water molecules per 1 mole ion pair GF=GFN*GF1 ! = GFN[=Nw=4] * GF1[=(1*1^1+1*1^1)/2/Nw=1/4] = 1 ! ONLY GAMA TYPE OF NH4NO3, NH4Cl, etc. needed so far IF(RH.GE.RHD) GAMA=RH**GF/((GFN*MWH20*(1./RH-1.)))**GF1 GAMA = min(GAMA,1.0) ! FOCUS ON 0-1 SCALE GAMA = max(GAMA,0.0) GAMA = (1.-GAMA)**K ! transplate into aqueous phase equillibrium and account for ! enhanced uptake of aerosol precursor gases with increasing RH ! (to match the results of ISORROPIA) ! CALCULATE RHD DEPENDENT EQ: IF RH < RHD => NH4NO3(s) <==> NH3 (g) + HNO3(g) (ISORROPIA) ! IF RH >> RHD => HNO3 (g) -> NO3 (aq) X00 = MAX(ZERO,MIN(TNa,GG*TSO4)) ! MAX SODIUM SULFATE X0 = MAX(ZERO,MIN(TNH4,GG*TSO4-X00)) ! MAX AMMOMIUM SULFATE X01 = MAX(ZERO,MIN(TNa-X00, TNO3)) ! MAX SODIUM NITRATE X1 = MAX(ZERO,MIN(TNH4-X0,TNO3-X01)) ! MAX AMMOMIUM NITRATE ! X02 = MAX(ZERO,MIN(TNa-X01-X00,TCl)) ! MAX SODIUM CHLORIDE X03 = MAX(ZERO,MIN(TNH4-X0-X1,TCl-X02))! MAX AMMOMIUM CHLORIDE X2 = MAX(TNH4-X1-X0-X03,ZERO) ! INTERIM RESIDUAL NH3 X3 = MAX(TNO3-X1-X01,ZERO) ! INTERIM RESIDUAL HNO3 X04 = MAX(TSO4-(X0+X00)/GG,ZERO) ! INTERIM RESIDUAL H2SO4 X05 = MAX(TCl-X03-X02,ZERO) ! INTERIM RESIDUAL HCl ! X06 = MAX(TNa-X02-X01-X00,ZERO) ! INTERIM RESIDUAL Na (should be zero for electro-neutrality in input data) ! ZKAN=2. IF(RH.GE.RHD) ZKAN=ZKAN*GAMA X4 = X2 + X3 X5 = SQRT(X4*X4+KAN*ZKAN*ZKAN) X6 = 0.5*(-X4+X5) X6 = MIN(X1,X6) GHCl = X05 ! INTERIM RESIDUAl HCl GNH3 = X2 + X6 ! INTERIM RESIDUAl NH3 GNO3 = X3 + X6 ! RESIDUAl HNO3 GSO4 = X04 ! RESIDUAl H2SO4 PNa = X02 + X01 + X00 ! RESIDUAl Na (neutralized) ZKAC=2. IF(RH.GE.RHD) ZKAC=ZKAC*GAMA X08 = GNH3 + GHCl X09 = SQRT(X08*X08+KAC*ZKAC*ZKAC) X10 = 0.5*(-X08+X09) X11 = MIN(X03,X10) GHCl = GHCl + X11 ! RESIDUAL HCl GNH3 = GNH3 + X11 ! RESIDUAL NH3 ! GO SAVE ... IF(GHCl.LT.0.) GHCl=0. IF(GSO4.LT.0.) GSO4=0. IF(GNH3.LT.0.) GNH3=0. IF(GNO3.LT.0.) GNO3=0. IF(PNa.LT.0.) PNa=0. IF(GSO4.GT.TSO4) GSO4=TSO4 IF(GNH3.GT.TNH4) GNH3=TNH4 IF(GNO3.GT.TNO3) GNO3=TNO3 IF(GHCl.GT.TCl) GHCl=TCl IF(PNa.GT.TNa) PNa=TNa ! IF(PNa.LT.TNa) print*,il,' PNa.LT.TNa => no electro-neutrality in input data! ',PNa,TNa ! DEFINE AQUEOUSE PHASE (NO SOLID NH4NO3 IF NO3/SO4>1, TEN BRINK, ET AL., 1996, ATMOS ENV, 24, 4251-4261) ! IF(TSO4.EQ.ZERO.AND.TNO3.GT.ZERO.OR.TNO3/TSO4.GE.1.) RHD=RH ! IF(IOPT.EQ.2.AND.RH.LT.RHD.OR.IOPT.EQ.2.AND.RH_HIST.LT.RH_HIST_DW) THEN ! SOLIDS / HYSTERESIS IF(RH_HIST.EQ.1.AND.RH.LT.RHD) THEN ! SOLIDS / HYSTERESIS ! EVERYTHING DRY, ONLY H2SO4 (GSO4) REMAINS IN THE AQUEOUSE PHASE ANH4 = 0. ASO4 = 0. ANO3 = 0. ACl = 0. ANa = 0. ELSE ! SUPERSATURATED SOLUTIONS NO SOLID FORMATION ASO4 = TSO4 - GSO4 ANH4 = TNH4 - GNH3 ANO3 = TNO3 - GNO3 ACl = TCl - GHCl ANa = PNa ENDIF ! SOLIDS / HYSTERESIS ! CALCULATE AEROSOL WATER [kg/m^3(air)] ! ! salt solutes: ! 1 = NACl, 2 = (NA)2SO4, 3 = NANO3, 4 = (NH4)2SO4, 5 = NH4NO3, 6 = NH4CL, 7 = 2H-SO4 ! 8 = NH4HSO4, 9 = NAHSO4, 10 = (NH4)3H(SO4)2 ! IF(ZFLAG.EQ.1.) WH2O = ASO4/M0( 2) + ANO3/M0(3) + ACl/M0(6) IF(ZFLAG.EQ.2.) WH2O = ASO4/M0( 9) + ANO3/M0(5) + ACl/M0(6) IF(ZFLAG.EQ.3.) WH2O = ASO4/M0( 8) + ANO3/M0(5) + ACl/M0(6) IF(ZFLAG.EQ.4.) WH2O = ASO4/M0( 8) + GSO4/M0(7) ! H2O_NO3 IF(ZFLAG.EQ.1.) H2O_NO3 = ANO3/M0(3) IF(ZFLAG.EQ.2.) H2O_NO3 = ANO3/M0(5) IF(ZFLAG.EQ.3.) H2O_NO3 = ANO3/M0(5) IF(ZFLAG.EQ.4.) H2O_NO3 = 0.0 ! CALCULATE AQUEOUS PHASE PROPERTIES ! PH = 9999. PH = 7. MOLAL = 0. HPLUS = 0. ZIONIC= 0. IF(WH2O.GT.0.) THEN ! CALCULATE AUTODISSOCIATION CONSTANT (KW) FOR WATER AKW=XKW*RH*WH2O*WH2O ! H2O <==> H+ + OH- with kw [mol^2/kg^2] AKW=AKW**0.5 ! [OH-] = [H+] [mol] ! trap zero [OH-] due to round-off: if ( AKW > 0.0 ) then ! Calculate hydrogen molality [mol/kg], i.e. H+ of the ions: ! Na+, NH4+, NO3-, Cl-, SO4--, HH-SO4- [mol/kg(water)] ! with [OH-] = kw/[H+] HPLUS = (-ANa/WH2O-ANH4/WH2O+ANO3/WH2O+ACl/WH2O+GG*ASO4/WH2O+GG*GSO4/WH2O+ & SQRT(( ANa/WH2O+ANH4/WH2O-ANO3/WH2O-ACl/WH2O-GG*ASO4/WH2O-GG*GSO4/WH2O)**2+XKW/AKW*WH2O))/2. ! trip zero [H+] due to round-off: if ( HPLUS > 0.0 ) then ! Calculate pH PH = -ALOG10(HPLUS) ! aerosol pH ! Calculate ionic strength [mol/kg] ZIONIC=0.5*(ANa+ANH4+ANO3+ACl+ASO4*GG*GG+GSO4*GG*GG+XKW/AKW*WH2O*WH2O) ZIONIC=ZIONIC/WH2O ! ionic strength [mol/kg] ! ZIONIC=min(ZIONIC,200.0) ! limit for output ! ZIONIC=max(ZIONIC,0.0) end if end if ! [OH-] > 0.0 ENDIF ! AQUEOUS PHASE ! !------------------------------------------------------- ! calculate diagnostic output consistent with other EQMs ... ! ASO4 = ASO4 + GSO4 ! assuming H2SO4 remains aqueous TNa = TNa * 1.e6 ! total input sodium (g+p) [umol/m^3] TSO4 = TSO4 * 1.e6 ! total input sulfate (g+p) [umol/m^3] TNH4 = TNH4 * 1.e6 ! total input ammonium (g+p) [umol/m^3] TNO3 = TNO3 * 1.e6 ! total input nitrate (g+p) [umol/m^3] TCl = TCl * 1.e6 ! total input chloride (g+p) [umol/m^3] TPo = TPo * 1.e6 ! total input potasium (g+p) [umol/m^3] TCa = TCa * 1.e6 ! total input calcium (g+p) [umol/m^3] TMg = TMg * 1.e6 ! total input magnesium(g+p) [umol/m^3] ! ! residual gas: GNH3 = GNH3 * 1.e6 ! residual NH3 GSO4 = GSO4 * 1.e6 ! residual H2SO4 GNO3 = GNO3 * 1.e6 ! residual HNO3 GHCl = GHCl * 1.e6 ! residual HCl ! total particulate matter (neutralized) PNH4=TNH4-GNH3 ! particulate ammonium [umol/m^3] !kt PNO3=TNO3-GNO3 ! particulate nitrate [umol/m^3] PNO3=max(0.,TNO3-GNO3) ! particulate nitrate [umol/m^3] PCl =TCl -GHCl ! particulate chloride [umol/m^3] PNa =TNa ! particulate sodium [umol/m^3] PSO4=TSO4 ! particulate sulfate [umol/m^3] ! liquid matter ASO4 = ASO4 * 1.e6 ! aqueous phase sulfate [umol/m^3] ANH4 = ANH4 * 1.e6 ! aqueous phase ammonium [umol/m^3] ANO3 = ANO3 * 1.e6 ! aqueous phase nitrate [umol/m^3] ACl = ACl * 1.e6 ! aqueous phase chloride [umol/m^3] ANa = ANa * 1.e6 ! aqueous phase sodium [umol/m^3] ! solid matter SNH4=PNH4-ANH4 ! solid phase ammonium [umol/m^3] SSO4=PSO4-ASO4 ! solid phase sulfate [umol/m^3] SNO3=PNO3-ANO3 ! solid phase nitrate [umol/m^3] SCl =PCl -ACl ! solid phase chloride [umol/m^3] SNa =PNa -ANa ! solid phase sodium [umol/m^3] ! GO SAVE ... IF(SNH4.LT.0.) SNH4=0. IF(SSO4.LT.0.) SSO4=0. IF(SNO3.LT.0.) SNO3=0. IF(SCl.LT.0.) SCl=0. IF(SNa.LT.0.) SNa=0. PM=SNH4+SSO4+SNO3+SNH4+SCl+SNa+ANH4+ASO4+ANO3+ACl+ANa ! total PM [umol/m^3] PMs=SNH4*MWNH4+SSO4*MWSO4+SNO3*MWNO3+SCl*MWCl+SNa*MWNa ! dry particulate matter (PM) [ug/m^3] PMt=PMs+ANH4*MWNH4+ASO4*MWSO4+ANO3*MWNO3+ACl*MWCl+ & ANa*MWNa ! total (dry + wet) PM, excl. H20 [ug/m^3] WH2O = WH2O * 1.e9 ! convert aerosol water from [kg/m^3] to [ug/m^3] H2O_NO3 = H2O_NO3 * 1.e9 ! convert NO3 water from [kg/m^3] to [ug/m^3] IF(WH2O.LT.1.e-3) WH2O=0. ! UPDATE HISTORY RH FOR HYSTERESIS (ONLINE CALCULATIONS ONLY) RH_HIST=2. ! wet IF(WH2O.EQ.0.) RH_HIST=1. ! dry RINC = 1. IF(PMt.GT.0.) RINC = (WH2O/PMt+1)**(1./3.) ! approx. radius increase due to water uptake IF(RINC.EQ.0.) RINC = 1. RATIONS = 0. IF(PSO4.GT.0.) RATIONS = PNO3/PSO4 ! nitrate / sulfate mol ratio GR = 0. IF(GNO3.GT.0.) GR = GNH3/GNO3 ! gas ratio = residual NH3 / residual HNO3 [-] DON = 0. IF((PNO3+2.*PSO4).GT.0.) DON = 100.*PNH4/(PNO3+2.*PSO4)! degree of neutralization by ammonia : ammonium / total nitrate + sulfate [%] NO3P = 0. IF(TNO3.GT.0.) NO3P = 100.*PNO3/TNO3 ! nitrate partitioning = nitrate / total nitrate [%] NH4P = 0. IF(TNH4.GT.0.) NH4P = 100.*PNH4/TNH4 ! ammonium partitioning = ammonium / total ammonium [%] ! ! store aerosol species for diagnostic output: !______________________________________________________________ ! Input values: yo(il, 1) = TT - 273.15 ! T [degC] yo(il, 2) = RH * 100.00 ! RH [%] yo(il, 3) = TNH4 ! total input ammonium (g+p) [umol/m^3] yo(il, 4) = TSO4 ! total input sulfate (g+p) [umol/m^3] yo(il, 5) = TNO3 ! total input nitrate (g+p) [umol/m^3] yo(il, 6) = TNa ! total input sodium (p) [umol/m^3] yo(il,33) = TCl ! total input chloride (g+p) [umol/m^3] yo(il, 7) = TPo ! total input potasium (p) [umol/m^3] yo(il,34) = TCa ! total input calcium (p) [umol/m^3] yo(il,35) = TMg ! total input magnesium(p) [umol/m^3] yo(il,25) = PX ! atmospheric pressure [hPa] ! Output values: yo(il, 8) = GHCL ! residual HCl (g) [umol/m^3] yo(il, 9) = GNO3 ! residual HNO3 (g) [umol/m^3] yo(il,10) = GNH3 ! residual NH3 (g) [umol/m^3] yo(il,11) = GSO4 ! residual H2SO4 (aq) [umol/m^3] yo(il,12) = WH2O ! aerosol Water (aq) [ug/m^3] yo(il,13) = PH ! aerosol pH [log] yo(il,14) = ZFLAG ! concnetration domain [1=SP,2=SN,3=SR,4=SVR] yo(il,15) = PM ! total particulate matter [umol/m^3] yo(il,16) = SNH4 ! solid ammonium (s) [umol/m^3] yo(il,17) = SNO3 ! solid nitrate (s) [umol/m^3] yo(il,18) = SSO4 ! solid sulfate (s) [umol/m^3] yo(il,19) = PNH4 ! particulate ammonium (p=a+s) [umol/m^3] yo(il,20) = PNO3 ! particulate nitrate (p=a+s) [umol/m^3] yo(il,21) = PSO4 ! particulate sulfate (p=a+s) [umol/m^3] yo(il,22) = RATIONS ! mol ratio Nitrate/Sulfate (p) [-] yo(il,23) = GAMA ! activity coefficient (e.g. NH4NO3) [-] yo(il,24) = ZIONIC ! ionic strength (aq) [mol/kg] yo(il,26) = PMt ! total PM (liquids & solids) [ug/m^3] yo(il,27) = PMs ! total PM (solid) [ug/m^3] yo(il,28) = RINC ! radius increase (H2O/PMt+1)**(1/3) [-] yo(il,29) = SCl ! solid chloride (s) [umol/m^3] yo(il,30) = SNa ! solid sodium (s) [umol/m^3] yo(il,31) = PCl ! particulate chloride (p=a+s) [umol/m^3] yo(il,32) = PNa ! particulate sodium (p=a+s) [umol/m^3] yo(il,36) = GG yo(il,37) = H2O_NO3 ! Water associated with nitrate [ug/m^3] enddo ! end subroutine eqsam_v03d end module eqsam_param