!$Id: !! ========================================================================= !! INCA - INteraction with Chemistry and Aerosols !! !! Copyright Laboratoire des Sciences du Climat et de l'Environnement (LSCE) !! Unite mixte CEA-CNRS-UVSQ !! !! Contributors to this INCA subroutine: !! !! Didier Hauglustaine, LSCE, hauglustaine@cea.fr !! !! This software is a computer program whose purpose is to simulate the !! atmospheric gas phase and aerosol composition. The model is designed to be !! used within a transport model or a general circulation model. This version !! of INCA was designed to be coupled to the LMDz GCM. LMDz-INCA accounts !! for emissions, transport (resolved and sub-grid scale), photochemical !! transformations, and scavenging (dry deposition and washout) of chemical !! species and aerosols interactively in the GCM. Several versions of the INCA !! model are currently used depending on the envisaged applications with the !! chemistry-climate model. !! !! This software is governed by the CeCILL license under French law and !! abiding by the rules of distribution of free software. You can use, !! modify and/ or redistribute the software under the terms of the CeCILL !! license as circulated by CEA, CNRS and INRIA at the following URL !! "http://www.cecill.info". !! !! As a counterpart to the access to the source code and rights to copy, !! modify and redistribute granted by the license, users are provided only !! with a limited warranty and the software's author, the holder of the !! economic rights, and the successive licensors have only limited !! liability. !! !! In this respect, the user's attention is drawn to the risks associated !! with loading, using, modifying and/or developing or reproducing the !! software by the user in light of its specific status of free software, !! that may mean that it is complicated to manipulate, and that also !! therefore means that it is reserved for developers and experienced !! professionals having in-depth computer knowledge. Users are therefore !! encouraged to load and test the software's suitability as regards their !! requirements in conditions enabling the security of their systems and/or !! data to be ensured and, more generally, to use and operate it in the !! same conditions as regards security. !! !! The fact that you are presently reading this means that you have had !! knowledge of the CeCILL license and that you accept its terms. !! ========================================================================= #include #if defined(NMHC) && defined(AER) SUBROUTINE soa_main(mmr,& mbar,& pmid,& tfld,& ASAPp1a_p,& ASAPp2a_p,& ASARp1a_p,& ASARp2a_p,& delt) USE CHEM_MODS, ONLY: nadv_mass,invariants USE AEROSOL_METEO, ONLY: zdens USE SOA_MOD USE INCA_DIM IMPLICIT NONE REAL, INTENT(in) :: delt REAL, DIMENSION(PLON,PLEV), INTENT(in) :: tfld ! temperature (K) REAL, DIMENSION(PLON,PLEV), INTENT(in) :: pmid ! air pressure at the center of the box (Pa) REAL, DIMENSION(PLON,PLEV), INTENT(in) :: mbar ! mean wet atmospheric mass (amu) REAL, DIMENSION(PLON,PLEV,PCNST), INTENT(inout) :: mmr ! mass mixing ratio of species REAL, DIMENSION(PLON,PLEV), INTENT(inout) :: ASAPp1a_p REAL, DIMENSION(PLON,PLEV), INTENT(inout) :: ASAPp2a_p REAL, DIMENSION(PLON,PLEV), INTENT(inout) :: ASARp1a_p REAL, DIMENSION(PLON,PLEV), INTENT(inout) :: ASARp2a_p REAL, DIMENSION(PLON) :: temp ! temperature (K) REAL, DIMENSION(PLON) :: fac1 ! mmr to molec/cm3 conversion factor INTEGER :: l,i IF (ifirst) CALL soa_init DO l=1,PLEV ! check that temperature range is within allowed limits DO i=1,PLON IF (ifix(tfld(i,l)) < tmin) PRINT*,'ERROR!!! Temperature in SOA is below the allowed value of ',tmin IF (ifix(tfld(i,l)) > tmax) PRINT*,'ERROR!!! Temperature in SOA is above the allowed value of ',tmax ENDDO ! change units from mmr to molecules/cm3 using fac1 fac1(:)=7.24e16*pmid(:,l)/tfld(:,l)*mbar(:,l) ! 7.24e16 is the inverse of Boltzmann constant *1.e-6 for twice g->kg DO i=1,PCNST if (ra(i) .eq. 0) cycle kty0(:,i)=mmr(:,l,i)/ra(i)*fac1(:) ENDDO ! initialize other needed parameters ! rr(:,:)=reaction_rates(:,l,:)/invariants(1,1,1) ! convert in order to use with molecules cm-3 concentrations kty=kty0 temp(:)=tfld(:,l) ! call main soa subroutines ! call soa_gasphase(dt) CALL soa_aerosolphase(temp) ! change units from molecules/cm3 to mmr DO i=1,PCNST if (ra(i) .eq. 0) cycle mmr(:,l,i)=kty(:,i)*ra(i)/fac1(:) ENDDO !Store changes for diagnostics (molec/cm3/s) ASAPp1a_p(:,l) = (kty(:,id_ASAPp1a)-kty0(:,id_ASAPp1a)) /delt ASAPp2a_p(:,l) = (kty(:,id_ASAPp2a)-kty0(:,id_ASAPp2a)) /delt ASARp1a_p(:,l) = (kty(:,id_ASARp1a)-kty0(:,id_ASARp1a)) /delt ASARp2a_p(:,l) = (kty(:,id_ASARp2a)-kty0(:,id_ASARp2a)) /delt ENDDO END SUBROUTINE soa_main SUBROUTINE soa_init USE CHEM_MODS, ONLY: adv_mass,nadv_mass USE SOA_MOD USE INCA_DIM IMPLICIT NONE REAL, PARAMETER :: insbc2om=1.1 ! conversion of insoluble black carbon to insoluble organic mass REAL, PARAMETER :: solbc2om=insbc2om*1.26 ! conversion of soluble black carbon to soluble organic mass ! (increased, based on Katrib et al., ACP, 2005) REAL, PARAMETER :: insoc2om=1.4 ! conversion of insoluble organic carbon to insoluble organic mass REAL, PARAMETER :: soloc2om=insoc2om*1.26 ! conversion of soluble organic carbon to soluble organic mass ! (increased, based on Katrib et al., ACP, 2005) INTEGER :: i,l,imf,jmf REAL :: KpCALC,molec2ug ! Initialize molecular weights ra=adv_mass ! These are the real molecular weights, needed for the mole fraction calculations mw=ra ! mw(id_AIBCM)=170. ! mw(id_ASBCM)=170.*1.26 !kt Katrib et al., ACP, 2005 ! mw(id_AIPOMM)=170. ! mw(id_ASPOMM)=170.*1.26 !kt Katrib et al., ACP, 2005 mw(id_AIBCM)=mw(id_AIBCM)*insbc2om mw(id_ASBCM)=mw(id_ASBCM)*solbc2om mw(id_AIPOMM)=mw(id_AIPOMM)*insoc2om mw(id_ASPOMM)=mw(id_ASPOMM)*soloc2om ! Initialize soaindex and tabulate m2ug soaindex=0 DO i=1,PCNST DO l=1,nsoa IF (issoa(l) == i) soaindex(i)=l ENDDO m2ug(i)=molec2ug(ra(i)) ENDDO DO imf=1,soacomp DO jmf=1,soacomp IF (Lamda(imf,jmf).EQ.-1.0) Lamda(imf,jmf)=Lamda(jmf,imf) ! symmetric interactions ENDDO ENDDO ! Tabulate kpart DO l=tmin,tmax kpart(soaindex(id_ASAPp1a),l)=KpCALC(72.9,1./15.7,float(l),295.) ! Presto et al., EST, 2005 kpart(soaindex(id_ASAPp2a),l)=KpCALC(72.9,1./385.,float(l),295.) ! Presto et al., EST, 2005 kpart(soaindex(id_ASARp1a),l)=KpCALC(72.9,0.301,float(l),300.) ! Song et al., EST, 2005 kpart(soaindex(id_ASARp2a),l)=KpCALC(72.9,0.008,float(l),300.) ! Song et al., EST, 2005 ENDDO ifirst=.FALSE. END SUBROUTINE soa_init SUBROUTINE soa_aerosolphase(temp) USE SOA_MOD USE INCA_DIM IMPLICIT NONE LOGICAL, PARAMETER :: Spart=.TRUE. ! partitioning on sulfur-containing particles (true) or not (false) LOGICAL, PARAMETER :: Npart=.TRUE. ! partitioning on nitrogen-containing particles (true) or not (false) LOGICAL, PARAMETER :: SOAevap=.TRUE. ! possibility of evaporation (true) or not (false) REAL, DIMENSION(PCNST) :: y0_ug,y_ug ! real,intent(in) :: dt REAL, DIMENSION(PLON), INTENT(in) :: temp INTEGER :: jl ! ! For correction of Kp due to composition ! REAL :: AEROtot,sigma_ij,sigma_ik,sigma_kj,zc,meanmw INTEGER :: imf,jmf,kmf REAL, DIMENSION(soacomp) :: xmf,zcoef REAL, DIMENSION(nsoa) :: kp ! partitioning coefficient ! Other SOA parameters REAL :: M0,M0old,PCP,M0a,M0b,M0err_curr,soamasstot,kty0tot!,soarat!,M0temp REAL, PARAMETER :: M0err_min=1.e-6,M0err_max=TINY(M0)*1.e10!1.e-15 REAL, DIMENSION(nsoa) :: soamass,partfact!,soaeva,soacond,soachemp INTEGER :: i,iternum ! character(len=11), parameter :: itermethod="chord" ! valid values are "bisectional" and "chord" REAL :: x1,x2,y1,y2,a,b DO jl=1,PLON ! ! Initialization ! y0_ug=0. y_ug=0. M0old=0. ! ! Change concentrations to ug/m3 ! DO i=1,nsoa y0_ug(issoa(i)-soaoffset)=kty0(jl,issoa(i)-soaoffset)*m2ug(issoa(i)-soaoffset) y0_ug(issoa(i))=kty0(jl,issoa(i))*m2ug(issoa(i)) y_ug(issoa(i)-soaoffset)=kty(jl,issoa(i)-soaoffset)*m2ug(issoa(i)-soaoffset) y_ug(issoa(i))=kty(jl,issoa(i))*m2ug(issoa(i)) ENDDO ! ! Calculate the primary aerosol able to absorb semivolatile compounds !kt 09/06/04 Isoprene SOA added, _without_ evaporation, like POC and BC ! Partitioning in both OC and BC ! Partitioning in (SO4+MSA) and/or NH4/NO3 is also an option ! ! PCP=kty(jl,iisopsoa)*molec2ug(ra(iisopsoa))+& PCP=kty0(jl,id_AIPOMM)*m2ug(id_AIPOMM) & +kty0(jl,id_ASPOMM)*m2ug(id_ASPOMM) & +kty0(jl,id_AIBCM)*m2ug(id_AIBCM) & +kty0(jl,id_ASBCM)*m2ug(id_ASBCM) IF(Spart) PCP=PCP+kty0(jl,id_CSSO4M)*m2ug(id_CSSO4M) & +kty0(jl,id_ASSO4M)*m2ug(id_ASSO4M)& +kty0(jl,id_CSMSAM)*m2ug(id_CSMSAM)& +kty0(jl,id_ASMSAM)*m2ug(id_ASMSAM) IF(Npart) PCP=PCP+kty0(jl,id_CSNH4M)*m2ug(id_CSNH4M) & +kty0(jl,id_ASNH4M)*m2ug(id_ASNH4M) & +kty0(jl,id_CINO3M)*m2ug(id_CINO3M) & +kty0(jl,id_CSNO3M)*m2ug(id_CSNO3M) & +kty0(jl,id_ASNO3M)*m2ug(id_ASNO3M) ! ! Correct the Kp to take into account the change of activity coefficient due to change in composition ! using the Wilson equation. Method described in Bowman and Karamalegos, EST, 2002, 36, 2701-2707. ! ! AEROtot=kty(jl,iisopsoa)+kty(jl,id_AIPOMM)+kty(jl,id_ASPOMM)+kty(jl,id_AIBCM)+kty(jl,id_ASBCM) AEROtot=kty0(jl,id_AIPOMM)/ra(id_AIPOMM)*mw(id_AIPOMM)+kty0(jl,id_ASPOMM)/ra(id_ASPOMM)*mw(id_ASPOMM)+ & kty0(jl,id_AIBCM)/ra(id_AIBCM)*mw(id_AIBCM)+kty0(jl,id_ASBCM)/ra(id_ASBCM)*mw(id_ASBCM) IF(Spart) AEROtot=AEROtot+kty0(jl,id_CSSO4M)/ra(id_CSSO4M)*mw(id_CSSO4M)+kty0(jl,id_ASSO4M)/ra(id_ASSO4M)*mw(id_ASSO4M)+ & kty0(jl,id_CSMSAM)/ra(id_CSMSAM)*mw(id_CSMSAM)+kty0(jl,id_ASMSAM)/ra(id_ASMSAM)*mw(id_ASMSAM) IF(Npart) AEROtot=AEROtot+kty0(jl,id_CSNH4M)/ra(id_CSNH4M)*mw(id_CSNH4M)+kty0(jl,id_ASNH4M)/ra(id_ASNH4M)*mw(id_ASNH4M)+ & kty0(jl,id_CINO3M)/ra(id_CINO3M)*mw(id_CINO3M)+kty0(jl,id_CSNO3M)/ra(id_CSNO3M)*mw(id_CSNO3M)+ & kty0(jl,id_ASNO3M)/ra(id_ASNO3M)*mw(id_ASNO3M) DO i=1,nsoa AEROtot=AEROtot+kty0(jl,issoa(i))/ra(issoa(i))*mw(issoa(i)) ENDDO IF (AEROtot.LE.0.) GOTO 15 ! ! Calculate mole fraction of individual species. ! xmf=0. xmf(imfter)=(kty0(jl,id_ASAPp1a)/ra(id_ASAPp1a)*mw(id_ASAPp1a)+kty0(jl,id_ASAPp2a)/ra(id_ASAPp2a)*mw(id_ASAPp2a))/AEROtot ! kty(jl,ibpinp1a_nox)+kty(jl,ibpinp2a_nox)+kty(jl,ibpinp1a_hox)+kty(jl,ibpinp2a_hox))/AEROtot xmf(imfaro)=(kty0(jl,id_ASARp1a)/ra(id_ASARp1a)*mw(id_ASARp1a)+kty0(jl,id_ASARp2a)/ra(id_ASARp2a)*mw(id_ASARp2a))/AEROtot ! kty(jl,ixylp1a_nox)+kty(jl,ixylp2a_nox)+kty(jl,ixylp1a_hox)+kty(jl,ixylp2a_hox))/AEROtot ! xmf(imfisopsoa)=kty(jl,iisopsoa)/AEROtot xmf(imfpocins)=kty0(jl,id_AIPOMM)/ra(id_AIPOMM)*mw(id_AIPOMM)/AEROtot xmf(imfpocsol)=kty0(jl,id_ASPOMM)/ra(id_ASPOMM)*mw(id_ASPOMM)/AEROtot xmf(imfbcins)=kty0(jl,id_AIBCM)/ra(id_AIBCM)*mw(id_AIBCM)/AEROtot xmf(imfbcsol)=kty0(jl,id_ASBCM)/ra(id_ASBCM)*mw(id_ASBCM)/AEROtot ! ! If partitioning does not occur on some aerosol species, it's mole fraction equals to zero, no matter the real ! concentration is, in order not to affect the calculation of the activity coefficient. ! IF(Spart) THEN xmf(imfinorg)=(kty0(jl,id_CSSO4M)/ra(id_CSSO4M)*mw(id_CSSO4M)+kty0(jl,id_ASSO4M)/ra(id_ASSO4M)*mw(id_ASSO4M))/AEROtot 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 ENDIF IF(Npart) THEN 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 xmf(imfinorg)=xmf(imfinorg)+(kty0(jl,id_CINO3M)/ra(id_CINO3M)*mw(id_CINO3M)+kty0(jl,id_CSNO3M)/ra(id_CSNO3M)*mw(id_CSNO3M)+ & kty0(jl,id_ASNO3M)/ra(id_ASNO3M)*mw(id_ASNO3M))/AEROtot ENDIF ! ! Calculate activity coefficient zcoef ! zcoef=1. ! do kmf=1,soacomp DO kmf=imfter,imfaro sigma_kj=0. DO jmf=1,soacomp sigma_kj=sigma_kj+xmf(jmf)*Lamda(kmf,jmf) ENDDO sigma_ik=0. DO imf=1,soacomp sigma_ij=0. DO jmf=1,soacomp sigma_ij=sigma_ij+xmf(jmf)*Lamda(imf,jmf) ENDDO sigma_ik=sigma_ik+xmf(imf)*Lamda(imf,kmf)/sigma_ij ENDDO zc=EXP(1.-LOG(sigma_kj)-sigma_ik) zcoef(kmf)=MAX(0.3,MIN(5.0,zc)) IF((zc.LT.0.3).OR.(zc.GT.5.0)) THEN PRINT*,'WARNING: zcoef set to',zcoef(kmf),' (was',zc,')' ENDIF ENDDO ! ! Calculate mean molecular weight ! meanmw=0. ! meanmw=kty0(jl,iisopsoa)*ra(iisopsoa)/AEROtot meanmw=meanmw+kty0(jl,id_AIPOMM)/ra(id_AIPOMM)*mw(id_AIPOMM)*mw(id_AIPOMM)/AEROtot meanmw=meanmw+kty0(jl,id_ASPOMM)/ra(id_ASPOMM)*mw(id_ASPOMM)*mw(id_ASPOMM)/AEROtot meanmw=meanmw+kty0(jl,id_AIBCM)/ra(id_AIBCM)*mw(id_AIBCM)*mw(id_AIBCM)/AEROtot meanmw=meanmw+kty0(jl,id_ASBCM)/ra(id_ASBCM)*mw(id_ASBCM)*mw(id_ASBCM)/AEROtot IF(Spart) THEN meanmw=meanmw+kty0(jl,id_CSSO4M)/ra(id_CSSO4M)*mw(id_CSSO4M)*mw(id_CSSO4M)/AEROtot meanmw=meanmw+kty0(jl,id_ASSO4M)/ra(id_ASSO4M)*mw(id_ASSO4M)*mw(id_ASSO4M)/AEROtot meanmw=meanmw+kty0(jl,id_CSMSAM)/ra(id_CSMSAM)*mw(id_CSMSAM)*mw(id_CSMSAM)/AEROtot meanmw=meanmw+kty0(jl,id_ASMSAM)/ra(id_ASMSAM)*mw(id_ASMSAM)*mw(id_ASMSAM)/AEROtot ENDIF IF(Npart) THEN meanmw=meanmw+kty0(jl,id_CSNH4M)/ra(id_CSNH4M)*mw(id_CSNH4M)*mw(id_CSNH4M)/AEROtot meanmw=meanmw+kty0(jl,id_ASNH4M)/ra(id_ASNH4M)*mw(id_ASNH4M)*mw(id_ASNH4M)/AEROtot meanmw=meanmw+kty0(jl,id_CINO3M)/ra(id_CINO3M)*mw(id_CINO3M)*mw(id_CINO3M)/AEROtot meanmw=meanmw+kty0(jl,id_CSNO3M)/ra(id_CSNO3M)*mw(id_CSNO3M)*mw(id_CSNO3M)/AEROtot meanmw=meanmw+kty0(jl,id_ASNO3M)/ra(id_ASNO3M)*mw(id_ASNO3M)*mw(id_ASNO3M)/AEROtot ENDIF DO i=1,nsoa meanmw=meanmw+kty0(jl,issoa(i))/ra(issoa(i))*mw(issoa(i))*mw(issoa(i))/AEROtot ENDDO ! ! Calculate final value of partitioning coefficient ! kp(soaindex(id_ASAPp1a))=kpart(soaindex(id_ASAPp1a),ifix(temp(jl)))/zcoef(imfter)*mw(id_ASAPp1a)/meanmw kp(soaindex(id_ASAPp2a))=kpart(soaindex(id_ASAPp2a),ifix(temp(jl)))/zcoef(imfter)*mw(id_ASAPp2a)/meanmw kp(soaindex(id_ASARp1a))=kpart(soaindex(id_ASARp1a),ifix(temp(jl)))/zcoef(imfaro)*mw(id_ASARp1a)/meanmw kp(soaindex(id_ASARp2a))=kpart(soaindex(id_ASARp2a),ifix(temp(jl)))/zcoef(imfaro)*mw(id_ASARp2a)/meanmw ! kp(ibpinp2a_nox)=kpart(jl,ibpinp1a_nox)/zcoef(imfter)*ra(ibpinp1a_nox)/meanmw ! kp(ibpinp2a_nox)=kpart(jl,ibpinp2a_nox)/zcoef(imfter)*ra(ibpinp2a_nox)/meanmw ! kp(ibpinp1a_hox)=kpart(jl,ibpinp1a_hox)/zcoef(imfter)*ra(ibpinp1a_hox)/meanmw ! kp(ibpinp2a_hox)=kpart(jl,ibpinp2a_hox)/zcoef(imfter)*ra(ibpinp2a_hox)/meanmw ! kp(itolp1a_nox)=kpart(jl,itolp1a_nox)/zcoef(imfaro)*ra(itolp1a_nox)/meanmw ! kp(itolp2a_nox)=kpart(jl,itolp2a_nox)/zcoef(imfaro)*ra(itolp2a_nox)/meanmw ! kp(itolp1a_hox)=kpart(jl,itolp1a_hox)/zcoef(imfaro)*ra(itolp1a_hox)/meanmw ! kp(itolp2a_hox)=kpart(jl,itolp2a_hox)/zcoef(imfaro)*ra(itolp2a_hox)/meanmw ! kp(ixylp1a_nox)=kpart(jl,ixylp1a_nox)/zcoef(imfaro)*ra(ixylp1a_nox)/meanmw ! kp(ixylp2a_nox)=kpart(jl,ixylp2a_nox)/zcoef(imfaro)*ra(ixylp2a_nox)/meanmw ! kp(ixylp1a_hox)=kpart(jl,ixylp1a_hox)/zcoef(imfaro)*ra(ixylp1a_hox)/meanmw ! kp(ixylp2a_hox)=kpart(jl,ixylp2a_hox)/zcoef(imfaro)*ra(ixylp2a_hox)/meanmw 15 CONTINUE !AEROtot=0. ! ! Add previous gas phase (soagas - after destruction), previous particulate phase and soachem in variable soamass (ug m-3) ! If no evaporation, add only soagas and soachem. soamass contains ONLY the mass that can be in the gas-phase ! DO i=1,nsoa ! soamass(i)=soagas(jl,i)+soachem(jl,i) soamass(i)=y0_ug(issoa(i)-soaoffset) IF(SOAevap) soamass(i)=soamass(i)+y0_ug(issoa(i)) ENDDO soamasstot=SUM(soamass) kty0tot=SUM(y0_ug(issoa(:))) IF (soamasstot == 0.) GOTO 60 ! no possible gas-phase species, nothing to do ! if (soamasstot == 0.) then ! no possible gas-phase species, nothing to do ! M0=0. ! goto 20 ! endif M0a=PCP M0b=PCP+soamasstot ! do i=1,nsoa ! M0b=M0b+soamass(i) IF(.NOT.SOAevap) THEN M0a=M0a+kty0tot M0b=M0b+kty0tot ! M0a=M0a+y0_ug(issoa(i)) ! M0b=M0b+y0_ug(issoa(i)) ENDIF ! enddo iternum=0 M0err_curr=MIN(M0err_min,MAX(M0err_max,soamasstot*M0err_min)) ! M0=M0b 11 CONTINUE ! Try with another M0 IF (iternum.EQ.0) THEN x2=M0 x1=x2-M0err_curr ELSE x1=M0 x2=x1+M0err_curr ENDIF ! if (x2 == x1) then ! PCP is by far larger than SOA (based on model precision), aerosol mass increase neglected IF ((x2 == x1) .OR. (soamasstot < M0err_min)) THEN ! PCP is by far larger than SOA (based on model precision), aerosol mass increase neglected ! if ((x2 == x1) .or. (soamasstot < M0err_curr)) then ! PCP is by far larger than SOA (based on model precision), aerosol mass increase neglected ! M0=x1 partfact(:)=kp(:)*M0/(1.+kp(:)*M0) GOTO 20 ENDIF y1=PCP-x1 y2=PCP-x2 IF(.NOT.SOAevap) THEN ! do i=1,nsoa y1=y1+kty0tot y2=y2+kty0tot ! y1=y1+y0_ug(issoa(i)) ! y2=y2+y0_ug(issoa(i)) ! enddo ENDIF DO i=1,nsoa partfact(i)=kp(i)*x2/(1.+kp(i)*x2) y2=y2+soamass(i)*partfact(i) partfact(i)=kp(i)*x1/(1.+kp(i)*x1) y1=y1+soamass(i)*partfact(i) ENDDO iternum=iternum+1 IF (ABS(y1).LT.M0err_curr/2.) GOTO 20 ! ! Stop iteration after 100 iterations (no solution found) ! IF(iternum.GE.100) GOTO 30 ! ! If -M0err/2 ug/m3. 1/molec2ug converts ug/m3 --> molec/cm3 !------------------------------------------------------------------ IMPLICIT NONE REAL, PARAMETER :: avo=6.0221367e23 ! Avogadro number REAL, PARAMETER :: AVOug=1.e12/avo REAL :: mw molec2ug=mw*AVOug END FUNCTION molec2ug #endif