#include "tm5.inc" MODULE mo_aero_m7 ! *mo_aero_m7* Contains parameters, switches and initialization ! routines for the m7 aerosol scheme. ! ! Authors: ! -------- ! E. Vignati (JRC/IES) 2005 ! P. Stier (MPI) 2001/2002 ! E. Vignati and J. Wilson (JRC / EI) 2000 IMPLICIT NONE SAVE !--- 1) Define and pre-set switches for the processes of M7: ----------------------- !--- Physical: LOGICAL :: lsnucl = .TRUE., & ! nucleation lscoag = .TRUE., & ! coagulation lscond = .TRUE. ! condensation of H2SO4 INTEGER :: nnucl = 1 ! Choice of the nucleation scheme: ! ! nnucl = 1 Vehkamaeki (2002) ! = 2 Kulmala (1998) -NOT RECOMMENDED- !--- Technical: LOGICAL :: lmass_diag = .FALSE. ! mass balance check in m7_interface !--- 2) Numbers of compounds and modes of m7: -------------------------------------- INTEGER, PARAMETER :: naermod=23, & !number of all compounds nmod=7, & !number of modes nss=2, & !number of sea salt compounds nsol=4, & !number of soluble compounds ngas=3, & !number of gaseous compounds nsulf=4,& !number of sulfate compounds ncomp=5 !number of compounds !--- 3) List of indexes corresponding to the compound masses and mode numbers:------ !--- 3.1) Mass index (in array aerml and ttn): ! ! Attention: ! The mass of sulfate compounds is always given in [molec. cm-3] ! whilst the mass of other compounds is given in [ug cm-3]. ! ! Compounds: ! ! so4 = sulphate ! bc = black carbon ! oc = organic carbon, ! ss = sea salt ! du = dust ! ! Modes: ! ! n = nucleation mode ! k = aitken mode ! a = accumulation mode ! c = coarse mode ! ! Type: ! ! s = soluble mode ! i = insoluble mode ! COMPOUND: INTEGER, PARAMETER :: & iso4ns=1, iso4ks=2, iso4as=3, iso4cs=4, & !- Sulfate ibcks =5, ibcas =6, ibccs =7, ibcki =8, & !- Black Carbon iocks =9, iocas=10, ioccs=11, iocki=12, & !- Organic Carbon issas=13, isscs=14, & !- Sea Salt iduas=15, iducs=16, iduai=17, iduci=18,& !- Dust isoans=19,isoaks=20,isoaas=21,isoacs=22,isoaki=23 !- SOA ! MODE: | | | | | ! nucl. | aitk. | acc. | coar. | aitk. | acc. | coar. | ! soluble | soluble | soluble | soluble | insol. | insol. | insol. | !--- 3.2) Number index (in array aernl): ! INTEGER, PARAMETER :: & inucs=1, iaits=2, iaccs=3, icoas=4, iaiti=5, iacci=6, icoai=7 ! MODE: | | | | | ! nucl. | aitk. | acc. | coar. | aitk. | acc. | coar. | ! soluble | soluble | soluble | soluble | insol. | insol. | insol. | !--- 4) Definition of the modes of M7: ------------------------------------------------------ !--- 4.1) Threshold radii between the different modes [cm]: ! Used for the repartititioning in m7_dconc. ! crdiv(jmod) is the lower bound and crdiv(jmod+1) is ! the upper bound of the respective geometric mode ! Default value for nucleation mode is modified by the ! choice of the nuclation scheme. REAL :: crdiv(4)=(/ 0.0005E-4, 0.005E-4, 0.05E-4, 0.5E-4 /) ! | | | ! | | | ! nucleation -- aitken - accum -- coarse mode !--- 4.2) Standard deviation for the modes: REAL, PARAMETER :: sigma(nmod)=(/ 1.59, 1.59, 1.59, 2.00, 1.59, 1.59, 2.00 /) !--- Natural logarithm of the standard deviation of each mode: ! Calulated in m7_initialize. REAL :: sigmaln(nmod) !--- 5) Conversion factors for lognormal particle size distributions: ------------- ! Calulated in m7_initialize. REAL :: cmr2ras(nmod) ! Conversion factor: count median radius to radius of average surface REAL :: cmr2mmr(nmod) ! Conversion factor: count median radius to mass mean radius REAL :: cmedr2mmedr(nmod) ! Conversion factor: count median radius to mass median radius REAL :: cmr2ram(nmod) ! Conversion factor: count median radius to radius of average mass REAL :: ram2cmr(nmod) ! Conversion factor: radius of average mass to count median radius !--- 6) Assumed thresholds for occurence of specific quantities: ------------- !@@@ To be done! ! REAL, PARAMETER :: cmin_aerml = 1.E-15 , ! Aerosol mass ! cmin_aernl = 1.E-10 , ! Aerosol number ! !--- 7) Chemical constants: ---------------------------------------------------- ! !--- Accomodation coefficient of H2SO4 on aerosols: ! (reduced for insoluble modes) REAL, PARAMETER :: caccso4(nmod) = (/ 1.0, 1.0, 1.0, 1.0, 0.3, 0.3, 0.3 /) !--- Critical relative humidity: REAL, PARAMETER :: crh = 0.45 ! Assumed relative humidity for the ! Na2SO4 / NaCl system below which ! crystalization occurs. ! (estimated from Tang, I.N.; JGR 102, D2 1883-1893) !--- 8) Physical constants: ---------------------------------------------------- ! !--- 8.1) General physical constants: REAL, PARAMETER :: bk = 1.38e-16, & ! Bolzman constant [] avo = 6.02217E+23,& ! Avogadro number [mol-1] rerg = 8.314E+7, & ! Ideal gas constant [erg.K-1.mole-1] r_kcal = 1.986E-3 ! Ideal gas constant [kcal K-1.mole-1] !--- 8.2) Type specific physical constants: ! REAL, PARAMETER :: dh2so4 = 1.841, & ! Density H2SO4 [g cm-3] ddust = 2.650, & ! Density du [g cm-3] !>>> TvN ! The density of BC is in the range 1.7 to 1.9 g/cm3. ! (Bond and Bergstrom, Aerosol Sci. Technol., 2006). ! We therefore adopt a value of 1.8 g/cm3. ! Details can be found in Bond et al. (JGR, 2013), ! and references therein: ! Park et al. (J. Nanoparticle Research, 2004) measured ! 1.77 +- 0.07 g/cm3 for the non-volatile components of diesel soot, ! and give a range 1.7-1.8 g/cm3 in their conclusions. ! Kondo et al. (Aerosol Sci. Techn., 2011) measured ! 1.718 +- 0.004 g/cm3 for fullerene soot. ! Schmid et al. (Environ. Sci. Technol., 2009) ! derive a value 1.8 +- 0.2 g/cm3 ! for elemental carbon from biomass burning. ! For comparison, in GLOMAP a value of 1.5 g/cm3 is used. ! Note that these density estimates measure the mass per volume ! occupied by the spherules, as should be the case: ! "If the radiative forcing of BC particles is to be ! calculated from their mass concentrations, ! as it is usually the case, density should represent ! the material density of the spherules, and not that ! of their ramiform (branched) or aciniform (packed) aggregates." ! (A. Gelencser, Carbonaceous Aerosol, Springer, 2004, p. 228). ! This is explained in more detail by Bond and Bergstrom ! (Aerosol Sci. Technol., 2005). ! !dbc = 2., & ! Density bc [g cm-3] dbc = 1.8, & ! Density bc [g cm-3] ! The density of OA is highly variable, ! but in any case substantially lower than 2 g/cm3. ! We adopt an average value of 1.3 g/cm3, ! based on a number of studies: ! Turpin and Lim (Aerosol Sci. Technol., 2001) ! suggest that 1.2 g/cm3 is a reasonable estimate. ! Lee et al. (ACP, 2010) measured 1.26 +- 0.24 g/cm3 ! at during FAME-2009 at Finokalia after evaporation ! Cross et al. (Aerosol Sci. Technol., 2007) assume an ! average bulk density of 1.27 g/cm3 ! Nakao et al. (Atmos. Environ., 2013) measure average densities ! for SOA between 1.22 and 1.42 g/cm3, depending of species. ! This is in line with predictions from the OA density model ! by Kuwata et al. (Environ. Sci. Technol., 2012), ! who find a range between 1.23 and 1.46 g/cm3 for SOA. ! This model can be used to estimate OA density ! as function of O:C and H:C elemental ratios, ! with an accuracy of 12% or more. ! As a further simplication, it is often assumed that ! H:C = 2 - O:C (e.g. Murphy et al., ACP, 2011). ! The model, however, is restricted to particle components ! having negligible quantities of additional elements, ! most notably nitrogen. ! Schmid et al. (Environ. Sci. Technol., 2009) ! derive a value of 1.39 +- 0.13 for OA from biomass burning. ! !doc = 2., & ! Density oc [g cm-3] doc = 1.3, & ! Density POM [g cm-3] !<<< TvN dnacl = 2.165, & ! Density NaCl [g cm-3] dna2so4 = 2.68, & ! Density Na2SO4 [g cm-3] dnahso4 = 2.435, & ! Density NaHSO4 [g cm-3] dh2o = 1.0, & ! Density H2O [g cm-3] wh2so4 = 98.0734, & ! Molecular weight H2SO4 [g mol-1] wh2o = 18.0, & ! Molecular weight H2O [g mol-1] wso4 = 96.0576, & ! Molecular weight SO4 [g mol-1] wso2 = 64.0, & ! Molecular weight SO2 [g mol-1] wna = 22.99, & ! Atomic weight Na [g mol-1] wcl = 35.453, & ! Atomic weight Cl [g mol-1] wnacl = 58.443, & ! Molecular weight NaCl [g mol-1] wna2so4 = 142.0376, & ! Molecular weight Na2SO4 [g mol-1] wnahso4 = 120.0555 ! Molecular weight NaHSO4 [g mol-1] !--- 9) Assumed parameters: ------------------------------------------------------ REAL, PARAMETER :: critn=100., & ! Assumed mass of an nucleated sulfate ! particle for the Kulmala scheme [molecules] fmax=0.95, & ! Factor that limits the condensation ! of sulfate to fmax times the available ! sulfate in the gas phase [1]. ! (m7_dgas) cLayerThickness = 1.0 ! Assumed required layer thickness of ! sulfate to transfer an insoluble ! particle to a soluble mode. It is ! given in units of layers of ! monomolecular sulfate. Determines the ! transfer rate from insoluble to ! soluble modes. !--- 10) Computational constants: ------------------------------------------------ REAL, PARAMETER :: sqrt2=1.4142136, & pi=3.141592654 !--- 11) Data used for the calculation of the aerosol properties ----------------- ! under ambient conditions: ! (Included the conversion from Pa to hPa in the first parameter.) REAL, PARAMETER :: wvb(17)= & (/ 95.80188, -28.5257, -1.082153, 0.1466501, & -20627.51, 0.0461242, -0.003935, -3.36115, & -0.00024137, 0.067938345, 0.00000649899, 8616124.373, & 1.168155578, -0.021317481, 0.000270358, -1353332314.0, & -0.002403805 /) REAL, PARAMETER :: gmb(9)= & (/ 1.036391467, 0.00728531, -0.011013887, -0.068887407, & 0.001047842, 0.001049607, 0.000740534, -1.081202685, & -0.0000029113 /) !--- 4) Logical mask for coagulation kernel: ------------------------------------- ! (The coagulation kernel mask is symmetric and not all ! values are used for physical considerations. As its ! calculation is very expensive, a mask is used to ! calculate only the necessarey elements.) LOGICAL :: locoagmask(nmod,nmod) DATA locoagmask(1:nmod,1) / .TRUE., .TRUE., .TRUE., .TRUE., .TRUE., .TRUE., .TRUE. / DATA locoagmask(1:nmod,2) / .FALSE., .TRUE., .TRUE., .TRUE., .TRUE., .TRUE., .TRUE. / DATA locoagmask(1:nmod,3) / .FALSE., .FALSE., .TRUE., .FALSE., .TRUE., .FALSE., .FALSE. / DATA locoagmask(1:nmod,4) / .FALSE., .FALSE., .FALSE., .FALSE., .FALSE., .FALSE., .FALSE. / DATA locoagmask(1:nmod,5) / .FALSE., .FALSE., .FALSE., .FALSE., .TRUE., .FALSE., .FALSE. / DATA locoagmask(1:nmod,6) / .FALSE., .FALSE., .FALSE., .FALSE., .FALSE., .FALSE., .FALSE. / DATA locoagmask(1:nmod,7) / .FALSE., .FALSE., .FALSE., .FALSE., .FALSE., .FALSE., .FALSE. / !--- 12) Service routines for initialization and auxiliary computations ---------- CONTAINS SUBROUTINE m7_initialize ! Purpose: ! --------- ! Initializes constants and parameters ! used in the m7 aerosol model. ! ! Author: ! --------- ! Philip Stier, MPI may 2001 ! ! Interface: ! --------- ! *m7_initialize* is called from *call_init_submodels* in *call_submodels* ! IMPLICIT NONE INTEGER :: jmod DO jmod=1, nmod !--- 1) Calculate conversion factors for lognormal distributions:---- ! Radius of average mass (ram) to count median radius (cmr) and ! vice versa. Count median radius to radius of average ! mass (ram). ! These factors depend on the standard deviation (sigma) ! of the lognormal distribution. ! (Based on the Hatch-Choate Conversins Equations; ! see Hinds, Chapter 4.5, 4.6 for more details. ! In particular equation 4.53.) !--- Count Median Radius to Mass Median Radius: cmedr2mmedr(jmod) = EXP(3.0*(LOG(sigma(jmod)))**2) !--- Count Median Radius to Mass Mean Radius: cmr2mmr(jmod) = EXP(3.5*(LOG(sigma(jmod)))**2) !--- Count Median Radius to Radius of Average Mass: cmr2ram(jmod) = EXP(1.5*(LOG(sigma(jmod)))**2) !--- Radius of Average Mass to Count Median Radius: ram2cmr(jmod) = 1. / cmr2ram(jmod) !--- Count Median Radius to Radius of Average Surface: cmr2ras(jmod) = EXP(1.0*(LOG(sigma(jmod)))**2) !--- 2) Calculate the natural logarithm of the standard deviation: sigmaln(jmod) = LOG(sigma(jmod)) END DO !--- 3) Set the lower mode boundary for the nucleation mode: ! (Depends on the choice of the nucleation scheme.) SELECT CASE (nnucl) CASE(1) crdiv(1)=0.0005E-4 CASE(2) crdiv(1)=( critn * wh2so4/avo/dh2so4*0.75/pi )**(1./3.) END SELECT !-------------------------------------------------------------------- END SUBROUTINE m7_initialize SUBROUTINE m7_coat(pm6rp_lon_lev, pcrtcst) ! Purpose: ! --------- ! *m7_coat* calculates the number of sulfate ! molecules required to coat a particle ! with cLayerThickness of sulfate ! ! Author: ! --------- ! Philip Stier, MPI 2001 ! ! Interface: ! --------- ! *m7_coat* is called from *m7_concoag* ! IMPLICIT NONE INTEGER :: jmod REAL :: pm6rp_lon_lev(nmod) ! Ambient radii for current ! longitude and level [cm] REAL :: pcrtcst(nmod) ! Critical constant, i.e. number of ! sulfate to cover an average particle ! of the mode with a layer of the ! thickness determined by cLayerThickness. REAL :: zras(nmod) ! Radius of average surface ! for a single particle [cm] REAL :: zas(nmod) ! Average surface ! for single particle [cm+2] REAL, PARAMETER :: csurf_molec = 2.39E-15 ! Average cross-section ! of a single H2SO4 molecule [cm+2] !--- 1) Calculate the radii of average surface for modes 5-7: zras(5) = pm6rp_lon_lev(5) * cmr2ras(5) zras(6) = pm6rp_lon_lev(6) * cmr2ras(6) zras(7) = pm6rp_lon_lev(7) * cmr2ras(7) DO jmod=5, 7 !--- 2) Calculate the average surface of an particle for modes 5-7: zas(jmod) = 4 * zras(jmod)**2 * pi !--- 3) Determine the number of sulfate molecules needed to form ! a cLayerThickness thick layer of sulfate on the particles ! in modes 5-7: pcrtcst(jmod) = (zas(jmod) / csurf_molec) * cLayerThickness END DO END SUBROUTINE m7_coat SUBROUTINE condensation_sink(NMODE,N,MEAN_IN,LAMBDA_IN,OUTPUT) !this is very approximative tool to calculate the division factors for organic condensation ! Assumes sigma = 1.5 (will underestimate condensation to coarse mode by factor of ~1-5) ! AA (2006) USE mo_kind,ONLY: dp INTEGER,INTENT(IN) :: NMODE !number of modes considered REAL(dp), DIMENSION(NMODE),INTENT(OUT) :: OUTPUT !output coefficients REAL(dp), DIMENSION(NMODE),INTENT(IN) :: N !number concentrations [arbitrary] REAL(dp), DIMENSION(NMODE),INTENT(IN) :: MEAN_IN !current geom median diameters of the modes [m] REAL(dp), INTENT(IN) :: LAMBDA_IN REAL(dp), DIMENSION(NMODE) :: MEAN REAL(dp) :: LAMBDA !local variables. REAL(dp), DIMENSION(3) :: COEF !the actual fit of CS to log10("mean") !fit of fit parameters in respect to log10(lambda) REAL(dp), DIMENSION(3), PARAMETER :: C1=(/ 0.0292 , 0.4376 , 1.5011/) !2nd order REAL(dp), DIMENSION(3), PARAMETER :: C2=(/0.4289 , 6.6937 , 25.4629/) !1st order REAL(dp), DIMENSION(3), PARAMETER :: C3=(/1.4346 , 23.0498 , 87.5852/) !0th order !REAL, DIMENSION(3) :: COEF_COARSE = (/ !Note : CS is without the standard non-size dependent !factor of Diff*2*pi, which cancels out REAL(dp), DIMENSION(NMODE) :: CS,LMEAN REAL(dp) :: LLAM CONTINUE LAMBDA = LAMBDA_IN MEAN = MEAN_IN !Check parameter validity IF (LAMBDA>7.1220e-07) LAMBDA =7.1220e-07 IF (LAMBDA<6.5100e-08) LAMBDA =6.5100e-08 WHERE (MEAN<1e-10) MEAN = 1E-10 WHERE (MEAN>10E-6) MEAN = 10E-6 !logs LMEAN = LOG10(MEAN) LLAM = LOG10(LAMBDA) !coefficents COEF(1) = LLAM*(LLAM*C1(1) + C1(2)) + C1(3) COEF(2) = LLAM*(LLAM*C2(1) + C2(2)) + C2(3) COEF(3) = LLAM*(LLAM*C3(1) + C3(2)) + C3(3) !Cond-sink CS= 10**(LMEAN*(LMEAN*COEF(1)+COEF(2))+COEF(3)) IF (SUM(N*CS) .GT. 1E-20) THEN OUTPUT = N*CS/SUM(N*CS) ELSE OUTPUT = 1./NMODE END IF END SUBROUTINE condensation_sink !!$ SUBROUTINE setaeroM7 !!$ !!$ ! *setaeroM7* modifies pre-set switches of the aeroM7ctl !!$ ! namelist for the configuration of the !!$ ! M7 component of the ECHAM/HAM aerosol model !!$ ! !!$ ! Authors: !!$ ! -------- !!$ ! Philip Stier, MPI-MET 12/2002 !!$ ! !!$ ! *setaeroM7* is called from *call_init_submodels* in *call_submodels* !!$ ! !!$ !!$ USE mo_mpi, ONLY: p_parallel, p_parallel_io, p_bcast, p_io !!$ USE mo_namelist, ONLY: position_nml, nnml, POSITIONED !!$ USE mo_doctor, ONLY: nout !!$ USE mo_exception, ONLY: finish !!$ !!$ IMPLICIT NONE !!$ !!$ INCLUDE 'aerom7ctl.inc' !!$ !!$ !--- Local variables: !!$ !!$ INTEGER :: ierr !!$ !!$ !--- 1) Read namelist: !!$ !!$ IF (p_parallel_io) THEN !!$ CALL position_nml ('AEROM7CTL', status=ierr) !!$ SELECT CASE (ierr) !!$ CASE (POSITIONED) !!$ READ (nnml, aerom7ctl) !!$ END SELECT !!$ ENDIF !!$ !!$ !--- 2) Broadcast over processors: !!$ !!$ IF (p_parallel) THEN !!$ CALL p_bcast (lsnucl, p_io) !!$ CALL p_bcast (lscoag, p_io) !!$ CALL p_bcast (lscond, p_io) !!$ CALL p_bcast (nnucl, p_io) !!$ CALL p_bcast (lmass_diag, p_io) !!$ END IF !!$ !!$ IF (.NOT. p_parallel) THEN !!$ WRITE(nout,*) '' !!$ WRITE(nout,*) '' !!$ WRITE(nout,*) '----------------------------------------------------------' !!$ WRITE(nout,*) '----------------------------------------------------------' !!$ WRITE(nout,*) '--- Initialization of settings for aerosol module M7 ---' !!$ WRITE(nout,*) '---' !!$ WRITE(nout,*) '--- Default values of aeroctl modified by setaero:' !!$ WRITE(nout,*) '---' !!$ WRITE(nout,*) '--- New settings: lsnucl = ', lsnucl !!$ WRITE(nout,*) '--- lscoag = ', lscoag !!$ WRITE(nout,*) '--- lscond = ', lscond !!$ IF (nnucl==1) THEN !!$ WRITE(nout,*) '--- nnucl = ', nnucl !!$ WRITE(nout,*) '--- => Vehkamaeki et al., 2002' !!$ ELSE IF (nnucl==2) THEN !!$ WRITE(nout,*) '--- nnucl = ', nnucl !!$ WRITE(nout,*) '--- => Kulmala et al., 1998' !!$ ELSE IF (lsnucl .AND. (nnucl/=1 .OR. nnucl/=2)) THEN !!$ WRITE(nout,*) '--- Error:' !!$ CALL finish('setaerom7', 'nucleation requested but no scheme selected') !!$ END IF !!$ WRITE(nout,*) '---' !!$ IF (lmass_diag) THEN !!$ WRITE(nout,*) '--- Mass balance check in m7_interface activated' !!$ ELSE !!$ WRITE(nout,*) '--- Mass balance check in m7_interface deactivated' !!$ END IF !!$ WRITE(nout,*) '---' !!$ WRITE(nout,*) '----------------------------------------------------------' !!$ WRITE(nout,*) '----------------------------------------------------------' !!$ WRITE(nout,*) '' !!$ WRITE(nout,*) '' !!$ END IF !!$ !!$ END SUBROUTINE setaeroM7 END MODULE mo_aero_m7