[6610] | 1 | !$Id: |
---|
| 2 | !! ========================================================================= |
---|
| 3 | !! INCA - INteraction with Chemistry and Aerosols |
---|
| 4 | !! |
---|
| 5 | !! Copyright Laboratoire des Sciences du Climat et de l'Environnement (LSCE) |
---|
| 6 | !! Unite mixte CEA-CNRS-UVSQ |
---|
| 7 | !! |
---|
| 8 | !! Contributors to this INCA subroutine: |
---|
| 9 | !! |
---|
| 10 | !! Didier Hauglustaine, LSCE, hauglustaine@cea.fr |
---|
| 11 | !! |
---|
| 12 | !! This software is a computer program whose purpose is to simulate the |
---|
| 13 | !! atmospheric gas phase and aerosol composition. The model is designed to be |
---|
| 14 | !! used within a transport model or a general circulation model. This version |
---|
| 15 | !! of INCA was designed to be coupled to the LMDz GCM. LMDz-INCA accounts |
---|
| 16 | !! for emissions, transport (resolved and sub-grid scale), photochemical |
---|
| 17 | !! transformations, and scavenging (dry deposition and washout) of chemical |
---|
| 18 | !! species and aerosols interactively in the GCM. Several versions of the INCA |
---|
| 19 | !! model are currently used depending on the envisaged applications with the |
---|
| 20 | !! chemistry-climate model. |
---|
| 21 | !! |
---|
| 22 | !! This software is governed by the CeCILL license under French law and |
---|
| 23 | !! abiding by the rules of distribution of free software. You can use, |
---|
| 24 | !! modify and/ or redistribute the software under the terms of the CeCILL |
---|
| 25 | !! license as circulated by CEA, CNRS and INRIA at the following URL |
---|
| 26 | !! "http://www.cecill.info". |
---|
| 27 | !! |
---|
| 28 | !! As a counterpart to the access to the source code and rights to copy, |
---|
| 29 | !! modify and redistribute granted by the license, users are provided only |
---|
| 30 | !! with a limited warranty and the software's author, the holder of the |
---|
| 31 | !! economic rights, and the successive licensors have only limited |
---|
| 32 | !! liability. |
---|
| 33 | !! |
---|
| 34 | !! In this respect, the user's attention is drawn to the risks associated |
---|
| 35 | !! with loading, using, modifying and/or developing or reproducing the |
---|
| 36 | !! software by the user in light of its specific status of free software, |
---|
| 37 | !! that may mean that it is complicated to manipulate, and that also |
---|
| 38 | !! therefore means that it is reserved for developers and experienced |
---|
| 39 | !! professionals having in-depth computer knowledge. Users are therefore |
---|
| 40 | !! encouraged to load and test the software's suitability as regards their |
---|
| 41 | !! requirements in conditions enabling the security of their systems and/or |
---|
| 42 | !! data to be ensured and, more generally, to use and operate it in the |
---|
| 43 | !! same conditions as regards security. |
---|
| 44 | !! |
---|
| 45 | !! The fact that you are presently reading this means that you have had |
---|
| 46 | !! knowledge of the CeCILL license and that you accept its terms. |
---|
| 47 | !! ========================================================================= |
---|
| 48 | |
---|
| 49 | #include <inca_define.h> |
---|
| 50 | |
---|
| 51 | #if defined(NMHC) && defined(AER) |
---|
| 52 | SUBROUTINE soa_main(mmr,& |
---|
| 53 | mbar,& |
---|
| 54 | pmid,& |
---|
| 55 | tfld,& |
---|
| 56 | ASAPp1a_p,& |
---|
| 57 | ASAPp2a_p,& |
---|
| 58 | ASARp1a_p,& |
---|
| 59 | ASARp2a_p,& |
---|
| 60 | delt) |
---|
| 61 | |
---|
| 62 | USE CHEM_MODS, ONLY: nadv_mass,invariants |
---|
| 63 | USE AEROSOL_METEO, ONLY: zdens |
---|
| 64 | USE SOA_MOD |
---|
| 65 | USE INCA_DIM |
---|
| 66 | |
---|
| 67 | IMPLICIT NONE |
---|
| 68 | |
---|
| 69 | REAL, INTENT(in) :: delt |
---|
| 70 | REAL, DIMENSION(PLON,PLEV), INTENT(in) :: tfld ! temperature (K) |
---|
| 71 | REAL, DIMENSION(PLON,PLEV), INTENT(in) :: pmid ! air pressure at the center of the box (Pa) |
---|
| 72 | REAL, DIMENSION(PLON,PLEV), INTENT(in) :: mbar ! mean wet atmospheric mass (amu) |
---|
| 73 | REAL, DIMENSION(PLON,PLEV,PCNST), INTENT(inout) :: mmr ! mass mixing ratio of species |
---|
| 74 | REAL, DIMENSION(PLON,PLEV), INTENT(inout) :: ASAPp1a_p |
---|
| 75 | REAL, DIMENSION(PLON,PLEV), INTENT(inout) :: ASAPp2a_p |
---|
| 76 | REAL, DIMENSION(PLON,PLEV), INTENT(inout) :: ASARp1a_p |
---|
| 77 | REAL, DIMENSION(PLON,PLEV), INTENT(inout) :: ASARp2a_p |
---|
| 78 | REAL, DIMENSION(PLON) :: temp ! temperature (K) |
---|
| 79 | REAL, DIMENSION(PLON) :: fac1 ! mmr to molec/cm3 conversion factor |
---|
| 80 | |
---|
| 81 | INTEGER :: l,i |
---|
| 82 | |
---|
| 83 | |
---|
| 84 | IF (ifirst) CALL soa_init |
---|
| 85 | |
---|
| 86 | |
---|
| 87 | DO l=1,PLEV |
---|
| 88 | ! check that temperature range is within allowed limits |
---|
| 89 | DO i=1,PLON |
---|
| 90 | IF (ifix(tfld(i,l)) < tmin) PRINT*,'ERROR!!! Temperature in SOA is below the allowed value of ',tmin |
---|
| 91 | IF (ifix(tfld(i,l)) > tmax) PRINT*,'ERROR!!! Temperature in SOA is above the allowed value of ',tmax |
---|
| 92 | ENDDO |
---|
| 93 | ! change units from mmr to molecules/cm3 using fac1 |
---|
| 94 | fac1(:)=7.24e16*pmid(:,l)/tfld(:,l)*mbar(:,l) ! 7.24e16 is the inverse of Boltzmann constant *1.e-6 for twice g->kg |
---|
| 95 | DO i=1,PCNST |
---|
| 96 | if (ra(i) .eq. 0) cycle |
---|
| 97 | kty0(:,i)=mmr(:,l,i)/ra(i)*fac1(:) |
---|
| 98 | ENDDO |
---|
| 99 | |
---|
| 100 | ! initialize other needed parameters |
---|
| 101 | ! rr(:,:)=reaction_rates(:,l,:)/invariants(1,1,1) ! convert in order to use with molecules cm-3 concentrations |
---|
| 102 | kty=kty0 |
---|
| 103 | temp(:)=tfld(:,l) |
---|
| 104 | |
---|
| 105 | ! call main soa subroutines |
---|
| 106 | ! call soa_gasphase(dt) |
---|
| 107 | |
---|
| 108 | CALL soa_aerosolphase(temp) |
---|
| 109 | |
---|
| 110 | ! change units from molecules/cm3 to mmr |
---|
| 111 | DO i=1,PCNST |
---|
| 112 | if (ra(i) .eq. 0) cycle |
---|
| 113 | mmr(:,l,i)=kty(:,i)*ra(i)/fac1(:) |
---|
| 114 | ENDDO |
---|
| 115 | |
---|
| 116 | !Store changes for diagnostics (molec/cm3/s) |
---|
| 117 | ASAPp1a_p(:,l) = (kty(:,id_ASAPp1a)-kty0(:,id_ASAPp1a)) /delt |
---|
| 118 | ASAPp2a_p(:,l) = (kty(:,id_ASAPp2a)-kty0(:,id_ASAPp2a)) /delt |
---|
| 119 | ASARp1a_p(:,l) = (kty(:,id_ASARp1a)-kty0(:,id_ASARp1a)) /delt |
---|
| 120 | ASARp2a_p(:,l) = (kty(:,id_ASARp2a)-kty0(:,id_ASARp2a)) /delt |
---|
| 121 | |
---|
| 122 | ENDDO |
---|
| 123 | |
---|
| 124 | END SUBROUTINE soa_main |
---|
| 125 | |
---|
| 126 | |
---|
| 127 | SUBROUTINE soa_init |
---|
| 128 | |
---|
| 129 | USE CHEM_MODS, ONLY: adv_mass,nadv_mass |
---|
| 130 | USE SOA_MOD |
---|
| 131 | USE INCA_DIM |
---|
| 132 | |
---|
| 133 | IMPLICIT NONE |
---|
| 134 | |
---|
| 135 | REAL, PARAMETER :: insbc2om=1.1 ! conversion of insoluble black carbon to insoluble organic mass |
---|
| 136 | REAL, PARAMETER :: solbc2om=insbc2om*1.26 ! conversion of soluble black carbon to soluble organic mass |
---|
| 137 | ! (increased, based on Katrib et al., ACP, 2005) |
---|
| 138 | REAL, PARAMETER :: insoc2om=1.4 ! conversion of insoluble organic carbon to insoluble organic mass |
---|
| 139 | REAL, PARAMETER :: soloc2om=insoc2om*1.26 ! conversion of soluble organic carbon to soluble organic mass |
---|
| 140 | ! (increased, based on Katrib et al., ACP, 2005) |
---|
| 141 | INTEGER :: i,l,imf,jmf |
---|
| 142 | REAL :: KpCALC,molec2ug |
---|
| 143 | |
---|
| 144 | ! Initialize molecular weights |
---|
| 145 | ra=adv_mass |
---|
| 146 | ! These are the real molecular weights, needed for the mole fraction calculations |
---|
| 147 | mw=ra |
---|
| 148 | ! mw(id_AIBCM)=170. |
---|
| 149 | ! mw(id_ASBCM)=170.*1.26 !kt Katrib et al., ACP, 2005 |
---|
| 150 | ! mw(id_AIPOMM)=170. |
---|
| 151 | ! mw(id_ASPOMM)=170.*1.26 !kt Katrib et al., ACP, 2005 |
---|
| 152 | mw(id_AIBCM)=mw(id_AIBCM)*insbc2om |
---|
| 153 | mw(id_ASBCM)=mw(id_ASBCM)*solbc2om |
---|
| 154 | mw(id_AIPOMM)=mw(id_AIPOMM)*insoc2om |
---|
| 155 | mw(id_ASPOMM)=mw(id_ASPOMM)*soloc2om |
---|
| 156 | |
---|
| 157 | ! Initialize soaindex and tabulate m2ug |
---|
| 158 | soaindex=0 |
---|
| 159 | DO i=1,PCNST |
---|
| 160 | DO l=1,nsoa |
---|
| 161 | IF (issoa(l) == i) soaindex(i)=l |
---|
| 162 | ENDDO |
---|
| 163 | m2ug(i)=molec2ug(ra(i)) |
---|
| 164 | ENDDO |
---|
| 165 | |
---|
| 166 | DO imf=1,soacomp |
---|
| 167 | DO jmf=1,soacomp |
---|
| 168 | IF (Lamda(imf,jmf).EQ.-1.0) Lamda(imf,jmf)=Lamda(jmf,imf) ! symmetric interactions |
---|
| 169 | ENDDO |
---|
| 170 | ENDDO |
---|
| 171 | |
---|
| 172 | ! Tabulate kpart |
---|
| 173 | DO l=tmin,tmax |
---|
| 174 | kpart(soaindex(id_ASAPp1a),l)=KpCALC(72.9,1./15.7,float(l),295.) ! Presto et al., EST, 2005 |
---|
| 175 | kpart(soaindex(id_ASAPp2a),l)=KpCALC(72.9,1./385.,float(l),295.) ! Presto et al., EST, 2005 |
---|
| 176 | kpart(soaindex(id_ASARp1a),l)=KpCALC(72.9,0.301,float(l),300.) ! Song et al., EST, 2005 |
---|
| 177 | kpart(soaindex(id_ASARp2a),l)=KpCALC(72.9,0.008,float(l),300.) ! Song et al., EST, 2005 |
---|
| 178 | ENDDO |
---|
| 179 | |
---|
| 180 | ifirst=.FALSE. |
---|
| 181 | |
---|
| 182 | END SUBROUTINE soa_init |
---|
| 183 | |
---|
| 184 | |
---|
| 185 | |
---|
| 186 | |
---|
| 187 | SUBROUTINE soa_aerosolphase(temp) |
---|
| 188 | |
---|
| 189 | USE SOA_MOD |
---|
| 190 | USE INCA_DIM |
---|
| 191 | |
---|
| 192 | |
---|
| 193 | IMPLICIT NONE |
---|
| 194 | |
---|
| 195 | LOGICAL, PARAMETER :: Spart=.TRUE. ! partitioning on sulfur-containing particles (true) or not (false) |
---|
| 196 | LOGICAL, PARAMETER :: Npart=.TRUE. ! partitioning on nitrogen-containing particles (true) or not (false) |
---|
| 197 | LOGICAL, PARAMETER :: SOAevap=.TRUE. ! possibility of evaporation (true) or not (false) |
---|
| 198 | |
---|
| 199 | REAL, DIMENSION(PCNST) :: y0_ug,y_ug |
---|
| 200 | ! real,intent(in) :: dt |
---|
| 201 | REAL, DIMENSION(PLON), INTENT(in) :: temp |
---|
| 202 | INTEGER :: jl |
---|
| 203 | ! |
---|
| 204 | ! For correction of Kp due to composition |
---|
| 205 | ! |
---|
| 206 | REAL :: AEROtot,sigma_ij,sigma_ik,sigma_kj,zc,meanmw |
---|
| 207 | INTEGER :: imf,jmf,kmf |
---|
| 208 | REAL, DIMENSION(soacomp) :: xmf,zcoef |
---|
| 209 | |
---|
| 210 | REAL, DIMENSION(nsoa) :: kp ! partitioning coefficient |
---|
| 211 | ! Other SOA parameters |
---|
| 212 | REAL :: M0,M0old,PCP,M0a,M0b,M0err_curr,soamasstot,kty0tot!,soarat!,M0temp |
---|
| 213 | REAL, PARAMETER :: M0err_min=1.e-6,M0err_max=TINY(M0)*1.e10!1.e-15 |
---|
| 214 | REAL, DIMENSION(nsoa) :: soamass,partfact!,soaeva,soacond,soachemp |
---|
| 215 | INTEGER :: i,iternum |
---|
| 216 | ! character(len=11), parameter :: itermethod="chord" ! valid values are "bisectional" and "chord" |
---|
| 217 | REAL :: x1,x2,y1,y2,a,b |
---|
| 218 | |
---|
| 219 | DO jl=1,PLON |
---|
| 220 | |
---|
| 221 | ! |
---|
| 222 | ! Initialization |
---|
| 223 | ! |
---|
| 224 | y0_ug=0. |
---|
| 225 | y_ug=0. |
---|
| 226 | M0old=0. |
---|
| 227 | ! |
---|
| 228 | ! Change concentrations to ug/m3 |
---|
| 229 | ! |
---|
| 230 | DO i=1,nsoa |
---|
| 231 | y0_ug(issoa(i)-soaoffset)=kty0(jl,issoa(i)-soaoffset)*m2ug(issoa(i)-soaoffset) |
---|
| 232 | y0_ug(issoa(i))=kty0(jl,issoa(i))*m2ug(issoa(i)) |
---|
| 233 | y_ug(issoa(i)-soaoffset)=kty(jl,issoa(i)-soaoffset)*m2ug(issoa(i)-soaoffset) |
---|
| 234 | y_ug(issoa(i))=kty(jl,issoa(i))*m2ug(issoa(i)) |
---|
| 235 | ENDDO |
---|
| 236 | ! |
---|
| 237 | ! Calculate the primary aerosol able to absorb semivolatile compounds |
---|
| 238 | !kt 09/06/04 Isoprene SOA added, _without_ evaporation, like POC and BC |
---|
| 239 | ! Partitioning in both OC and BC |
---|
| 240 | ! Partitioning in (SO4+MSA) and/or NH4/NO3 is also an option |
---|
| 241 | ! |
---|
| 242 | ! PCP=kty(jl,iisopsoa)*molec2ug(ra(iisopsoa))+& |
---|
| 243 | PCP=kty0(jl,id_AIPOMM)*m2ug(id_AIPOMM) & |
---|
| 244 | +kty0(jl,id_ASPOMM)*m2ug(id_ASPOMM) & |
---|
| 245 | +kty0(jl,id_AIBCM)*m2ug(id_AIBCM) & |
---|
| 246 | +kty0(jl,id_ASBCM)*m2ug(id_ASBCM) |
---|
| 247 | |
---|
| 248 | IF(Spart) PCP=PCP+kty0(jl,id_CSSO4M)*m2ug(id_CSSO4M) & |
---|
| 249 | +kty0(jl,id_ASSO4M)*m2ug(id_ASSO4M)& |
---|
| 250 | +kty0(jl,id_CSMSAM)*m2ug(id_CSMSAM)& |
---|
| 251 | +kty0(jl,id_ASMSAM)*m2ug(id_ASMSAM) |
---|
| 252 | IF(Npart) PCP=PCP+kty0(jl,id_CSNH4M)*m2ug(id_CSNH4M) & |
---|
| 253 | +kty0(jl,id_ASNH4M)*m2ug(id_ASNH4M) & |
---|
| 254 | +kty0(jl,id_CINO3M)*m2ug(id_CINO3M) & |
---|
| 255 | +kty0(jl,id_CSNO3M)*m2ug(id_CSNO3M) & |
---|
| 256 | +kty0(jl,id_ASNO3M)*m2ug(id_ASNO3M) |
---|
| 257 | ! |
---|
| 258 | ! Correct the Kp to take into account the change of activity coefficient due to change in composition |
---|
| 259 | ! using the Wilson equation. Method described in Bowman and Karamalegos, EST, 2002, 36, 2701-2707. |
---|
| 260 | ! |
---|
| 261 | ! AEROtot=kty(jl,iisopsoa)+kty(jl,id_AIPOMM)+kty(jl,id_ASPOMM)+kty(jl,id_AIBCM)+kty(jl,id_ASBCM) |
---|
| 262 | AEROtot=kty0(jl,id_AIPOMM)/ra(id_AIPOMM)*mw(id_AIPOMM)+kty0(jl,id_ASPOMM)/ra(id_ASPOMM)*mw(id_ASPOMM)+ & |
---|
| 263 | kty0(jl,id_AIBCM)/ra(id_AIBCM)*mw(id_AIBCM)+kty0(jl,id_ASBCM)/ra(id_ASBCM)*mw(id_ASBCM) |
---|
| 264 | IF(Spart) AEROtot=AEROtot+kty0(jl,id_CSSO4M)/ra(id_CSSO4M)*mw(id_CSSO4M)+kty0(jl,id_ASSO4M)/ra(id_ASSO4M)*mw(id_ASSO4M)+ & |
---|
| 265 | kty0(jl,id_CSMSAM)/ra(id_CSMSAM)*mw(id_CSMSAM)+kty0(jl,id_ASMSAM)/ra(id_ASMSAM)*mw(id_ASMSAM) |
---|
| 266 | IF(Npart) AEROtot=AEROtot+kty0(jl,id_CSNH4M)/ra(id_CSNH4M)*mw(id_CSNH4M)+kty0(jl,id_ASNH4M)/ra(id_ASNH4M)*mw(id_ASNH4M)+ & |
---|
| 267 | kty0(jl,id_CINO3M)/ra(id_CINO3M)*mw(id_CINO3M)+kty0(jl,id_CSNO3M)/ra(id_CSNO3M)*mw(id_CSNO3M)+ & |
---|
| 268 | kty0(jl,id_ASNO3M)/ra(id_ASNO3M)*mw(id_ASNO3M) |
---|
| 269 | DO i=1,nsoa |
---|
| 270 | AEROtot=AEROtot+kty0(jl,issoa(i))/ra(issoa(i))*mw(issoa(i)) |
---|
| 271 | ENDDO |
---|
| 272 | IF (AEROtot.LE.0.) GOTO 15 |
---|
| 273 | ! |
---|
| 274 | ! Calculate mole fraction of individual species. |
---|
| 275 | ! |
---|
| 276 | xmf=0. |
---|
| 277 | xmf(imfter)=(kty0(jl,id_ASAPp1a)/ra(id_ASAPp1a)*mw(id_ASAPp1a)+kty0(jl,id_ASAPp2a)/ra(id_ASAPp2a)*mw(id_ASAPp2a))/AEROtot |
---|
| 278 | ! kty(jl,ibpinp1a_nox)+kty(jl,ibpinp2a_nox)+kty(jl,ibpinp1a_hox)+kty(jl,ibpinp2a_hox))/AEROtot |
---|
| 279 | xmf(imfaro)=(kty0(jl,id_ASARp1a)/ra(id_ASARp1a)*mw(id_ASARp1a)+kty0(jl,id_ASARp2a)/ra(id_ASARp2a)*mw(id_ASARp2a))/AEROtot |
---|
| 280 | ! kty(jl,ixylp1a_nox)+kty(jl,ixylp2a_nox)+kty(jl,ixylp1a_hox)+kty(jl,ixylp2a_hox))/AEROtot |
---|
| 281 | ! xmf(imfisopsoa)=kty(jl,iisopsoa)/AEROtot |
---|
| 282 | xmf(imfpocins)=kty0(jl,id_AIPOMM)/ra(id_AIPOMM)*mw(id_AIPOMM)/AEROtot |
---|
| 283 | xmf(imfpocsol)=kty0(jl,id_ASPOMM)/ra(id_ASPOMM)*mw(id_ASPOMM)/AEROtot |
---|
| 284 | xmf(imfbcins)=kty0(jl,id_AIBCM)/ra(id_AIBCM)*mw(id_AIBCM)/AEROtot |
---|
| 285 | xmf(imfbcsol)=kty0(jl,id_ASBCM)/ra(id_ASBCM)*mw(id_ASBCM)/AEROtot |
---|
| 286 | ! |
---|
| 287 | ! If partitioning does not occur on some aerosol species, it's mole fraction equals to zero, no matter the real |
---|
| 288 | ! concentration is, in order not to affect the calculation of the activity coefficient. |
---|
| 289 | ! |
---|
| 290 | IF(Spart) THEN |
---|
| 291 | xmf(imfinorg)=(kty0(jl,id_CSSO4M)/ra(id_CSSO4M)*mw(id_CSSO4M)+kty0(jl,id_ASSO4M)/ra(id_ASSO4M)*mw(id_ASSO4M))/AEROtot |
---|
| 292 | xmf(imfinorg)=xmf(imfinorg)+(kty0(jl,id_CSMSAM)/ra(id_CSMSAM)*mw(id_CSMSAM)+kty0(jl,id_ASMSAM)/ra(id_ASMSAM)*mw(id_ASMSAM))/AEROtot |
---|
| 293 | ENDIF |
---|
| 294 | IF(Npart) THEN |
---|
| 295 | xmf(imfinorg)=xmf(imfinorg)+(kty0(jl,id_CSNH4M)/ra(id_CSNH4M)*mw(id_CSNH4M)+kty0(jl,id_ASNH4M)/ra(id_ASNH4M)*mw(id_ASNH4M))/AEROtot |
---|
| 296 | xmf(imfinorg)=xmf(imfinorg)+(kty0(jl,id_CINO3M)/ra(id_CINO3M)*mw(id_CINO3M)+kty0(jl,id_CSNO3M)/ra(id_CSNO3M)*mw(id_CSNO3M)+ & |
---|
| 297 | kty0(jl,id_ASNO3M)/ra(id_ASNO3M)*mw(id_ASNO3M))/AEROtot |
---|
| 298 | ENDIF |
---|
| 299 | ! |
---|
| 300 | ! Calculate activity coefficient zcoef |
---|
| 301 | ! |
---|
| 302 | zcoef=1. |
---|
| 303 | ! do kmf=1,soacomp |
---|
| 304 | DO kmf=imfter,imfaro |
---|
| 305 | sigma_kj=0. |
---|
| 306 | DO jmf=1,soacomp |
---|
| 307 | sigma_kj=sigma_kj+xmf(jmf)*Lamda(kmf,jmf) |
---|
| 308 | ENDDO |
---|
| 309 | sigma_ik=0. |
---|
| 310 | DO imf=1,soacomp |
---|
| 311 | sigma_ij=0. |
---|
| 312 | DO jmf=1,soacomp |
---|
| 313 | sigma_ij=sigma_ij+xmf(jmf)*Lamda(imf,jmf) |
---|
| 314 | ENDDO |
---|
| 315 | sigma_ik=sigma_ik+xmf(imf)*Lamda(imf,kmf)/sigma_ij |
---|
| 316 | ENDDO |
---|
| 317 | zc=EXP(1.-LOG(sigma_kj)-sigma_ik) |
---|
| 318 | zcoef(kmf)=MAX(0.3,MIN(5.0,zc)) |
---|
| 319 | IF((zc.LT.0.3).OR.(zc.GT.5.0)) THEN |
---|
| 320 | PRINT*,'WARNING: zcoef set to',zcoef(kmf),' (was',zc,')' |
---|
| 321 | ENDIF |
---|
| 322 | ENDDO |
---|
| 323 | ! |
---|
| 324 | ! Calculate mean molecular weight |
---|
| 325 | ! |
---|
| 326 | meanmw=0. |
---|
| 327 | ! meanmw=kty0(jl,iisopsoa)*ra(iisopsoa)/AEROtot |
---|
| 328 | meanmw=meanmw+kty0(jl,id_AIPOMM)/ra(id_AIPOMM)*mw(id_AIPOMM)*mw(id_AIPOMM)/AEROtot |
---|
| 329 | meanmw=meanmw+kty0(jl,id_ASPOMM)/ra(id_ASPOMM)*mw(id_ASPOMM)*mw(id_ASPOMM)/AEROtot |
---|
| 330 | meanmw=meanmw+kty0(jl,id_AIBCM)/ra(id_AIBCM)*mw(id_AIBCM)*mw(id_AIBCM)/AEROtot |
---|
| 331 | meanmw=meanmw+kty0(jl,id_ASBCM)/ra(id_ASBCM)*mw(id_ASBCM)*mw(id_ASBCM)/AEROtot |
---|
| 332 | IF(Spart) THEN |
---|
| 333 | meanmw=meanmw+kty0(jl,id_CSSO4M)/ra(id_CSSO4M)*mw(id_CSSO4M)*mw(id_CSSO4M)/AEROtot |
---|
| 334 | meanmw=meanmw+kty0(jl,id_ASSO4M)/ra(id_ASSO4M)*mw(id_ASSO4M)*mw(id_ASSO4M)/AEROtot |
---|
| 335 | meanmw=meanmw+kty0(jl,id_CSMSAM)/ra(id_CSMSAM)*mw(id_CSMSAM)*mw(id_CSMSAM)/AEROtot |
---|
| 336 | meanmw=meanmw+kty0(jl,id_ASMSAM)/ra(id_ASMSAM)*mw(id_ASMSAM)*mw(id_ASMSAM)/AEROtot |
---|
| 337 | ENDIF |
---|
| 338 | IF(Npart) THEN |
---|
| 339 | meanmw=meanmw+kty0(jl,id_CSNH4M)/ra(id_CSNH4M)*mw(id_CSNH4M)*mw(id_CSNH4M)/AEROtot |
---|
| 340 | meanmw=meanmw+kty0(jl,id_ASNH4M)/ra(id_ASNH4M)*mw(id_ASNH4M)*mw(id_ASNH4M)/AEROtot |
---|
| 341 | meanmw=meanmw+kty0(jl,id_CINO3M)/ra(id_CINO3M)*mw(id_CINO3M)*mw(id_CINO3M)/AEROtot |
---|
| 342 | meanmw=meanmw+kty0(jl,id_CSNO3M)/ra(id_CSNO3M)*mw(id_CSNO3M)*mw(id_CSNO3M)/AEROtot |
---|
| 343 | meanmw=meanmw+kty0(jl,id_ASNO3M)/ra(id_ASNO3M)*mw(id_ASNO3M)*mw(id_ASNO3M)/AEROtot |
---|
| 344 | ENDIF |
---|
| 345 | DO i=1,nsoa |
---|
| 346 | meanmw=meanmw+kty0(jl,issoa(i))/ra(issoa(i))*mw(issoa(i))*mw(issoa(i))/AEROtot |
---|
| 347 | ENDDO |
---|
| 348 | ! |
---|
| 349 | ! Calculate final value of partitioning coefficient |
---|
| 350 | ! |
---|
| 351 | kp(soaindex(id_ASAPp1a))=kpart(soaindex(id_ASAPp1a),ifix(temp(jl)))/zcoef(imfter)*mw(id_ASAPp1a)/meanmw |
---|
| 352 | kp(soaindex(id_ASAPp2a))=kpart(soaindex(id_ASAPp2a),ifix(temp(jl)))/zcoef(imfter)*mw(id_ASAPp2a)/meanmw |
---|
| 353 | kp(soaindex(id_ASARp1a))=kpart(soaindex(id_ASARp1a),ifix(temp(jl)))/zcoef(imfaro)*mw(id_ASARp1a)/meanmw |
---|
| 354 | kp(soaindex(id_ASARp2a))=kpart(soaindex(id_ASARp2a),ifix(temp(jl)))/zcoef(imfaro)*mw(id_ASARp2a)/meanmw |
---|
| 355 | ! kp(ibpinp2a_nox)=kpart(jl,ibpinp1a_nox)/zcoef(imfter)*ra(ibpinp1a_nox)/meanmw |
---|
| 356 | ! kp(ibpinp2a_nox)=kpart(jl,ibpinp2a_nox)/zcoef(imfter)*ra(ibpinp2a_nox)/meanmw |
---|
| 357 | ! kp(ibpinp1a_hox)=kpart(jl,ibpinp1a_hox)/zcoef(imfter)*ra(ibpinp1a_hox)/meanmw |
---|
| 358 | ! kp(ibpinp2a_hox)=kpart(jl,ibpinp2a_hox)/zcoef(imfter)*ra(ibpinp2a_hox)/meanmw |
---|
| 359 | ! kp(itolp1a_nox)=kpart(jl,itolp1a_nox)/zcoef(imfaro)*ra(itolp1a_nox)/meanmw |
---|
| 360 | ! kp(itolp2a_nox)=kpart(jl,itolp2a_nox)/zcoef(imfaro)*ra(itolp2a_nox)/meanmw |
---|
| 361 | ! kp(itolp1a_hox)=kpart(jl,itolp1a_hox)/zcoef(imfaro)*ra(itolp1a_hox)/meanmw |
---|
| 362 | ! kp(itolp2a_hox)=kpart(jl,itolp2a_hox)/zcoef(imfaro)*ra(itolp2a_hox)/meanmw |
---|
| 363 | ! kp(ixylp1a_nox)=kpart(jl,ixylp1a_nox)/zcoef(imfaro)*ra(ixylp1a_nox)/meanmw |
---|
| 364 | ! kp(ixylp2a_nox)=kpart(jl,ixylp2a_nox)/zcoef(imfaro)*ra(ixylp2a_nox)/meanmw |
---|
| 365 | ! kp(ixylp1a_hox)=kpart(jl,ixylp1a_hox)/zcoef(imfaro)*ra(ixylp1a_hox)/meanmw |
---|
| 366 | ! kp(ixylp2a_hox)=kpart(jl,ixylp2a_hox)/zcoef(imfaro)*ra(ixylp2a_hox)/meanmw |
---|
| 367 | 15 CONTINUE !AEROtot=0. |
---|
| 368 | ! |
---|
| 369 | ! Add previous gas phase (soagas - after destruction), previous particulate phase and soachem in variable soamass (ug m-3) |
---|
| 370 | ! If no evaporation, add only soagas and soachem. soamass contains ONLY the mass that can be in the gas-phase |
---|
| 371 | ! |
---|
| 372 | DO i=1,nsoa |
---|
| 373 | ! soamass(i)=soagas(jl,i)+soachem(jl,i) |
---|
| 374 | soamass(i)=y0_ug(issoa(i)-soaoffset) |
---|
| 375 | IF(SOAevap) soamass(i)=soamass(i)+y0_ug(issoa(i)) |
---|
| 376 | ENDDO |
---|
| 377 | soamasstot=SUM(soamass) |
---|
| 378 | kty0tot=SUM(y0_ug(issoa(:))) |
---|
| 379 | IF (soamasstot == 0.) GOTO 60 ! no possible gas-phase species, nothing to do |
---|
| 380 | ! if (soamasstot == 0.) then ! no possible gas-phase species, nothing to do |
---|
| 381 | ! M0=0. |
---|
| 382 | ! goto 20 |
---|
| 383 | ! endif |
---|
| 384 | M0a=PCP |
---|
| 385 | M0b=PCP+soamasstot |
---|
| 386 | ! do i=1,nsoa |
---|
| 387 | ! M0b=M0b+soamass(i) |
---|
| 388 | IF(.NOT.SOAevap) THEN |
---|
| 389 | M0a=M0a+kty0tot |
---|
| 390 | M0b=M0b+kty0tot |
---|
| 391 | ! M0a=M0a+y0_ug(issoa(i)) |
---|
| 392 | ! M0b=M0b+y0_ug(issoa(i)) |
---|
| 393 | ENDIF |
---|
| 394 | ! enddo |
---|
| 395 | iternum=0 |
---|
| 396 | M0err_curr=MIN(M0err_min,MAX(M0err_max,soamasstot*M0err_min)) |
---|
| 397 | ! |
---|
| 398 | M0=M0b |
---|
| 399 | 11 CONTINUE ! Try with another M0 |
---|
| 400 | IF (iternum.EQ.0) THEN |
---|
| 401 | x2=M0 |
---|
| 402 | x1=x2-M0err_curr |
---|
| 403 | ELSE |
---|
| 404 | x1=M0 |
---|
| 405 | x2=x1+M0err_curr |
---|
| 406 | ENDIF |
---|
| 407 | ! if (x2 == x1) then ! PCP is by far larger than SOA (based on model precision), aerosol mass increase neglected |
---|
| 408 | IF ((x2 == x1) .OR. (soamasstot < M0err_min)) THEN ! PCP is by far larger than SOA (based on model precision), aerosol mass increase neglected |
---|
| 409 | ! if ((x2 == x1) .or. (soamasstot < M0err_curr)) then ! PCP is by far larger than SOA (based on model precision), aerosol mass increase neglected |
---|
| 410 | ! M0=x1 |
---|
| 411 | partfact(:)=kp(:)*M0/(1.+kp(:)*M0) |
---|
| 412 | GOTO 20 |
---|
| 413 | ENDIF |
---|
| 414 | y1=PCP-x1 |
---|
| 415 | y2=PCP-x2 |
---|
| 416 | IF(.NOT.SOAevap) THEN |
---|
| 417 | ! do i=1,nsoa |
---|
| 418 | y1=y1+kty0tot |
---|
| 419 | y2=y2+kty0tot |
---|
| 420 | ! y1=y1+y0_ug(issoa(i)) |
---|
| 421 | ! y2=y2+y0_ug(issoa(i)) |
---|
| 422 | ! enddo |
---|
| 423 | ENDIF |
---|
| 424 | DO i=1,nsoa |
---|
| 425 | partfact(i)=kp(i)*x2/(1.+kp(i)*x2) |
---|
| 426 | y2=y2+soamass(i)*partfact(i) |
---|
| 427 | partfact(i)=kp(i)*x1/(1.+kp(i)*x1) |
---|
| 428 | y1=y1+soamass(i)*partfact(i) |
---|
| 429 | ENDDO |
---|
| 430 | iternum=iternum+1 |
---|
| 431 | IF (ABS(y1).LT.M0err_curr/2.) GOTO 20 |
---|
| 432 | ! |
---|
| 433 | ! Stop iteration after 100 iterations (no solution found) |
---|
| 434 | ! |
---|
| 435 | IF(iternum.GE.100) GOTO 30 |
---|
| 436 | ! |
---|
| 437 | ! If -M0err/2<M0temp<M0err/2, M0 is the solution (goto 20) |
---|
| 438 | ! Otherwise change the limits of M0a and M0b and restart (goto 11) |
---|
| 439 | ! |
---|
| 440 | a=(y2-y1)/(x2-x1) |
---|
| 441 | b=y1-a*x1 |
---|
| 442 | M0old=M0 |
---|
| 443 | M0=-b/a |
---|
| 444 | IF (M0 == M0old) THEN ! solution will not be able to be found due to very small mass |
---|
| 445 | partfact(i)=kp(i)*M0/(1.+kp(i)*M0) |
---|
| 446 | GOTO 20 |
---|
| 447 | ENDIF |
---|
| 448 | IF (M0.LT.0.) THEN |
---|
| 449 | PRINT*,'WARNING: M0 set to zero (was ',M0,')' |
---|
| 450 | M0=0. |
---|
| 451 | ENDIF |
---|
| 452 | GOTO 11 |
---|
| 453 | ! else |
---|
| 454 | ! stop "Invalid iteration method" |
---|
| 455 | ! endif |
---|
| 456 | ! |
---|
| 457 | ! No solution found, aerosol phase set to previous values, chemistry still applies to gas-phase species |
---|
| 458 | ! |
---|
| 459 | 30 CONTINUE |
---|
| 460 | WRITE(*,"('WARNING: Too many iterations. SOA forced to the previous values. M0-PCP=',e, & |
---|
| 461 | ' soamass=',e,' M0err_curr=',e,' iternum=',i)") & |
---|
| 462 | M0-PCP,soamasstot,M0err_curr,iternum |
---|
| 463 | DO i=1,nsoa |
---|
| 464 | y_ug(issoa(i))=y0_ug(issoa(i)) |
---|
| 465 | y_ug(issoa(i)-soaoffset)=y0_ug(issoa(i)-soaoffset) |
---|
| 466 | ENDDO |
---|
| 467 | GOTO 60 |
---|
| 468 | ! |
---|
| 469 | ! Found solution for M0 |
---|
| 470 | ! |
---|
| 471 | 20 CONTINUE |
---|
| 472 | IF (M0.NE.0.) THEN |
---|
| 473 | DO i=1,nsoa |
---|
| 474 | partfact(i)=MAX(0.,partfact(i)) |
---|
| 475 | IF (SOAevap) THEN |
---|
| 476 | y_ug(issoa(i))=soamass(i)*partfact(i) ! Calculate new aerosol-phase concentration |
---|
| 477 | y_ug(issoa(i)-soaoffset)=y_ug(issoa(i))/(kp(i)*M0) ! Calculate new gas-phase concentration |
---|
| 478 | ELSE |
---|
| 479 | y_ug(issoa(i))=y0_ug(issoa(i))+soamass(i)*partfact(i) |
---|
| 480 | y_ug(issoa(i)-soaoffset)=MAX(0.,y_ug(issoa(i)-soaoffset)-soamass(i)*partfact(i)) |
---|
| 481 | ENDIF |
---|
| 482 | ! Variables needed for budget calculations (stored in molecules cm-3) |
---|
| 483 | ! soaeva(i)=kty0(jl,issoa(i))*(partfact(i)-1.) |
---|
| 484 | ! soacond(i)=soagas(jl,i)*partfact(i)/molec2ug(ra(issoa(i)-soaoffset)) |
---|
| 485 | ! soachemp(i)=soachem(jl,i)*partfact(i)/molec2ug(ra(issoa(i)-soaoffset)) |
---|
| 486 | ENDDO |
---|
| 487 | ELSE |
---|
| 488 | PRINT*,'WARNING: SOA forced to previous values (M0=',M0,', PCP=',PCP,')' |
---|
| 489 | DO i=1,nsoa |
---|
| 490 | y_ug(issoa(i))=y0_ug(issoa(i)) |
---|
| 491 | y_ug(issoa(i)-soaoffset)=y0_ug(issoa(i)-soaoffset) |
---|
| 492 | ! soaeva(i)=0. |
---|
| 493 | ! soacond(i)=0. |
---|
| 494 | ! soachemp(i)=0. |
---|
| 495 | ENDDO |
---|
| 496 | ENDIF |
---|
| 497 | ! |
---|
| 498 | ! Check mass conservation |
---|
| 499 | !kt not needed, only for test purposes. Generally soarat~1e-7-1e-8 |
---|
| 500 | 60 CONTINUE |
---|
| 501 | ! |
---|
| 502 | ! Change concentrations to molec/cm3 |
---|
| 503 | ! |
---|
| 504 | DO i=1,nsoa |
---|
| 505 | kty(jl,issoa(i)-soaoffset)=y_ug(issoa(i)-soaoffset)/m2ug(issoa(i)-soaoffset) |
---|
| 506 | kty(jl,issoa(i))=y_ug(issoa(i))/m2ug(issoa(i)) |
---|
| 507 | ENDDO |
---|
| 508 | |
---|
| 509 | ENDDO !jl |
---|
| 510 | |
---|
| 511 | END SUBROUTINE soa_aerosolphase |
---|
| 512 | |
---|
| 513 | |
---|
| 514 | REAL FUNCTION KpCALC(dH,Ksc,temp,Tsc) |
---|
| 515 | !---------------------------------------------------------------------- |
---|
| 516 | !**** KpCALC calculates the temperature dependence of the partitioning coefficients |
---|
| 517 | ! |
---|
| 518 | ! interface |
---|
| 519 | ! --------- |
---|
| 520 | ! call KpCALC(dH,Ksc,temp) |
---|
| 521 | ! where: dH is the enthalpy of vaporization of the selected species, in cal/mol |
---|
| 522 | ! Ksc is the partitioning coefficient of the selected species at temperature Tsc |
---|
| 523 | ! temp is the temperature in current box |
---|
| 524 | ! |
---|
| 525 | ! reference |
---|
| 526 | ! --------- |
---|
| 527 | ! Chung and Seinfeld, JGR, 2002 |
---|
| 528 | !------------------------------------------------------------------ |
---|
| 529 | IMPLICIT NONE |
---|
| 530 | ! real, parameter :: dH=0. ! No temperature dependence on Kp due to vapor pressure |
---|
| 531 | ! real, parameter :: dH=4.2e4 ! Chung and Seinfeld, 2002, in cal/mol |
---|
| 532 | ! real, parameter :: dH=7.29e4 ! Pun et al., June 2002, submitted (Pandis group), in cal/mol |
---|
| 533 | ! real, parameter :: dH=7.9e4 ! Andersson-Skold and Simpson, 2001, in cal/mol |
---|
| 534 | REAL, INTENT(in):: dH ! Kcal/mol |
---|
| 535 | REAL, INTENT(in):: Ksc,temp,Tsc |
---|
| 536 | REAL, PARAMETER :: rgas=8.3144 ! universal gas constant (J mol-1 K-1) |
---|
| 537 | |
---|
| 538 | KpCALC=Ksc*(temp/Tsc)*EXP(dH*1000./rgas*(1./temp-1./Tsc)) |
---|
| 539 | |
---|
| 540 | END FUNCTION KpCALC |
---|
| 541 | |
---|
| 542 | |
---|
| 543 | REAL FUNCTION molec2ug(mw) |
---|
| 544 | !---------------------------------------------------------------------- |
---|
| 545 | !**** molec2ug converts molec/cm3 --> ug/m3. 1/molec2ug converts ug/m3 --> molec/cm3 |
---|
| 546 | !------------------------------------------------------------------ |
---|
| 547 | IMPLICIT NONE |
---|
| 548 | REAL, PARAMETER :: avo=6.0221367e23 ! Avogadro number |
---|
| 549 | REAL, PARAMETER :: AVOug=1.e12/avo |
---|
| 550 | REAL :: mw |
---|
| 551 | |
---|
| 552 | molec2ug=mw*AVOug |
---|
| 553 | |
---|
| 554 | END FUNCTION molec2ug |
---|
| 555 | |
---|
| 556 | |
---|
| 557 | #endif |
---|