source: CONFIG/publications/ICOLMDZORINCA_CO2_Transport_GMD_2023/INCA/src/INCA_SRC/soa.F90 @ 6610

Last change on this file since 6610 was 6610, checked in by acosce, 10 months ago

INCA used for ICOLMDZORINCA_CO2_Transport_GMD_2023
File size: 23.1 KB
Line 
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) 
52SUBROUTINE 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
124END SUBROUTINE soa_main
125
126
127SUBROUTINE 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
182END SUBROUTINE soa_init
183
184
185
186
187SUBROUTINE 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
36715  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
39911  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    !
45930  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    !
47120  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
50060  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
511END SUBROUTINE soa_aerosolphase
512
513
514REAL 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 
540END FUNCTION KpCALC
541
542
543REAL 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 
554END FUNCTION molec2ug
555
556
557#endif
Note: See TracBrowser for help on using the repository browser.