!### macro's ##################################################### ! #define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr #define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if #define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if ! #include "tm5.inc" ! !################################################################# module eqsam_param ! ! use GO, only : gol, goPr, goErr ! ! 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,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) ! ! ! 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] ! ! ! 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. ! ! ! 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) ! ! 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] ! 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 ! enddo ! ! ! end subroutine eqsam_v03d ! ! end module eqsam_param