MODULE traccoag_mod ! This module calculates the concentration of aerosol particles in certain size bins ! considering coagulation and sedimentation. USE lmdz_abort_physic, ONLY: abort_physic CONTAINS SUBROUTINE traccoag(pdtphys, gmtime, debutphy, julien, & presnivs, xlat, xlon, pphis, pphi, & t_seri, pplay, paprs, sh, rh, tr_seri) USE phys_local_var_mod, ONLY: mdw, R2SO4, DENSO4, f_r_wet, surf_PM25_sulf, & budg_emi_ocs, budg_emi_so2, budg_emi_h2so4, budg_emi_part, & R2SO4B, DENSO4B, f_r_wetB, sulfmmr, SAD_sulfate, sulfmmr_mode, nd_mode USE dimphy USE infotrac_phy, ONLY: nbtr_bin, nbtr_sulgas, nbtr, id_SO2_strat USE aerophys USE lmdz_geometry, ONLY: cell_area, boundslat USE lmdz_grid_phy USE lmdz_phys_mpi_data, ONLY: is_mpi_root USE lmdz_phys_para, ONLY: gather, scatter USE phys_cal_mod, ONLY: year_len, year_cur, mth_cur, day_cur, hour USE sulfate_aer_mod USE phys_local_var_mod, ONLY: stratomask USE lmdz_yomcst USE lmdz_print_control, ONLY: lunout USE strataer_local_var_mod IMPLICIT NONE ! Input argument !--------------- REAL, INTENT(IN) :: pdtphys ! Pas d'integration pour la physique (seconde) REAL, INTENT(IN) :: gmtime ! Heure courante LOGICAL, INTENT(IN) :: debutphy ! le flag de l'initialisation de la physique INTEGER, INTENT(IN) :: julien ! Jour julien REAL, DIMENSION(klev), INTENT(IN) :: presnivs! pressions approximat. des milieux couches (en PA) REAL, DIMENSION(klon), INTENT(IN) :: xlat ! latitudes pour chaque point REAL, DIMENSION(klon), INTENT(IN) :: xlon ! longitudes pour chaque point REAL, DIMENSION(klon), INTENT(IN) :: pphis ! geopotentiel du sol REAL, DIMENSION(klon, klev), INTENT(IN) :: pphi ! geopotentiel de chaque couche REAL, DIMENSION(klon, klev), INTENT(IN) :: t_seri ! Temperature REAL, DIMENSION(klon, klev), INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa) REAL, DIMENSION(klon, klev + 1), INTENT(IN) :: paprs ! pression pour chaque inter-couche (en Pa) REAL, DIMENSION(klon, klev), INTENT(IN) :: sh ! humidite specifique REAL, DIMENSION(klon, klev), INTENT(IN) :: rh ! humidite relative ! Output argument !---------------- REAL, DIMENSION(klon, klev, nbtr), INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA] ! Local variables !---------------- REAL :: m_aer_emiss_vol_daily ! daily injection mass emission REAL :: m_aer ! aerosol mass INTEGER :: it, k, i, j, ilon, ilev, itime, i_int, ieru LOGICAL, DIMENSION(klon, klev) :: is_strato ! true = above tropopause, false = below REAL, DIMENSION(klon, klev) :: m_air_gridbox ! mass of air in every grid box [kg] REAL, DIMENSION(klon_glo, klev, nbtr) :: tr_seri_glo ! Concentration Traceur [U/KgA] REAL, DIMENSION(klev + 1) :: altLMDz ! altitude of layer interfaces in m REAL, DIMENSION(klev) :: f_lay_emiss ! fraction of emission for every vertical layer REAL :: f_lay_sum ! sum of layer emission fractions REAL :: alt ! altitude for integral calculation INTEGER, PARAMETER :: n_int_alt = 10 ! number of subintervals for integration over Gaussian emission profile REAL, DIMENSION(nbtr_bin) :: r_bin ! particle radius in size bin [m] REAL, DIMENSION(nbtr_bin) :: r_lower ! particle radius at lower bin boundary [m] REAL, DIMENSION(nbtr_bin) :: r_upper ! particle radius at upper bin boundary [m] REAL, DIMENSION(nbtr_bin) :: m_part_dry ! mass of one dry particle in size bin [kg] REAL :: zrho ! Density of air [kg/m3] REAL :: zdz ! thickness of atm. model layer in m REAL, DIMENSION(klev) :: zdm ! mass of atm. model layer in kg REAL, DIMENSION(klon, klev) :: dens_aer ! density of aerosol particles [kg/m3 aerosol] with default H2SO4 mass fraction REAL :: emission ! emission REAL :: theta_min, theta_max ! for SAI computation between two latitudes REAL :: dlat_loc REAL :: latmin, latmax, lonmin, lonmax ! lat/lon min/max for injection REAL :: sigma_alt, altemiss ! injection altitude + sigma for distrib REAL :: pdt, stretchlong ! physic timestep, stretch emission over one day INTEGER :: injdur_sai ! injection duration for SAI case [days] INTEGER :: yr, is_bissext IF (is_mpi_root .AND. flag_verbose_strataer) THEN WRITE(lunout, *) 'in traccoag: date from phys_cal_mod =', year_cur, '-', mth_cur, '-', day_cur, '-', hour WRITE(lunout, *) 'IN traccoag flag_emit: ', flag_emit ENDIF ! radius [m] DO it = 1, nbtr_bin r_bin(it) = mdw(it) / 2. ENDDO !--set boundaries of size bins DO it = 1, nbtr_bin IF (it==1) THEN r_upper(it) = sqrt(r_bin(it + 1) * r_bin(it)) r_lower(it) = r_bin(it)**2. / r_upper(it) ELSEIF (it==nbtr_bin) THEN r_lower(it) = sqrt(r_bin(it) * r_bin(it - 1)) r_upper(it) = r_bin(it)**2. / r_lower(it) ELSE r_lower(it) = sqrt(r_bin(it) * r_bin(it - 1)) r_upper(it) = sqrt(r_bin(it + 1) * r_bin(it)) ENDIF ENDDO IF (debutphy .AND. is_mpi_root) THEN DO it = 1, nbtr_bin WRITE(lunout, *) 'radius bin', it, ':', r_bin(it), '(from', r_lower(it), 'to', r_upper(it), ')' ENDDO ENDIF !--initialising LOGICAL is_strato from stratomask is_strato(:, :) = .FALSE. WHERE (stratomask>0.5) is_strato = .TRUE. IF(flag_new_strat_compo) THEN IF(debutphy) WRITE(lunout, *) 'traccoag: COMPO/DENSITY (Tabazadeh 97) + H2O kelvin effect', flag_new_strat_compo ! STRACOMP (H2O, P, t_seri, R -> R2SO4 + Kelvin effect) : Taba97, Socol, etc... CALL stracomp_kelvin(sh, t_seri, pplay) ELSE IF(debutphy) WRITE(lunout, *) 'traccoag: COMPO from Bekki 2D model', flag_new_strat_compo ! STRACOMP (H2O, P, t_seri -> aerosol composition (R2SO4)) ! H2SO4 mass fraction in aerosol (%) CALL stracomp(sh, t_seri, pplay) ! aerosol density (gr/cm3) CALL denh2sa(t_seri) ! compute factor for converting dry to wet radius (for every grid box) f_r_wet(:, :) = (dens_aer_dry / (DENSO4(:, :) * 1000.) / (R2SO4(:, :) / 100.))**(1. / 3.) ENDIF !--calculate mass of air in every grid box DO ilon = 1, klon DO ilev = 1, klev m_air_gridbox(ilon, ilev) = (paprs(ilon, ilev) - paprs(ilon, ilev + 1)) / RG * cell_area(ilon) ENDDO ENDDO !--initialise emission diagnostics IF (nErupt > 0 .AND. (flag_emit == 1 .OR. flag_emit == 4)) budg_emi(:, 1) = 0.0 budg_emi_ocs(:) = 0.0 budg_emi_so2(:) = 0.0 budg_emi_h2so4(:) = 0.0 budg_emi_part(:) = 0.0 !--sulfur emission, depending on chosen scenario (flag_emit) SELECT CASE(flag_emit) CASE(0) ! background aerosol ! do nothing (no emission) CASE(1) ! volcanic eruption !--only emit on day of eruption ! stretch emission over one day of Pinatubo eruption DO ieru = 1, nErupt IF (year_cur==year_emit_vol(ieru).AND.mth_cur==mth_emit_vol(ieru).AND.& day_cur>=day_emit_vol(ieru).AND.day_cur<(day_emit_vol(ieru) + injdur)) THEN ! daily injection mass emission m_aer = m_aer_emiss_vol(ieru, 1) / (REAL(injdur) * REAL(ponde_lonlat_vol(ieru))) !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss) m_aer = m_aer * (mSO2mol / mSatom) WRITE(lunout, *) 'IN traccoag m_aer_emiss_vol(ieru)=', m_aer_emiss_vol(ieru, 1), & 'ponde_lonlat_vol(ieru)=', ponde_lonlat_vol(ieru), '(injdur*ponde_lonlat_vol(ieru))', & (injdur * ponde_lonlat_vol(ieru)), 'm_aer_emiss_vol_daily=', m_aer, 'ieru=', ieru WRITE(lunout, *) 'IN traccoag, dlon=', dlon latmin = xlat_min_vol(ieru) latmax = xlat_max_vol(ieru) lonmin = xlon_min_vol(ieru) lonmax = xlon_max_vol(ieru) altemiss = altemiss_vol(ieru) sigma_alt = sigma_alt_vol(ieru) pdt = pdtphys ! stretch emission over one day of eruption stretchlong = 1. CALL STRATEMIT(pdtphys, pdt, xlat, xlon, t_seri, pplay, paprs, tr_seri, & m_aer, latmin, latmax, lonmin, lonmax, altemiss, sigma_alt, id_SO2_strat, & stretchlong, 1, 0) ENDIF ! emission period ENDDO ! eruption number CASE(2) ! stratospheric aerosol injections (SAI) ! Computing duration of SAI in days... ! ... starting from 0... injdur_sai = 0 ! ... then adding whole years from first to (n-1)th... DO yr = year_emit_sai_start, year_emit_sai_end - 1 ! (n % 4 == 0) and (n % 100 != 0 or n % 400 == 0) is_bissext = (MOD(yr, 4)==0) .AND. (MOD(yr, 100) /= 0 .OR. MOD(yr, 400) == 0) injdur_sai = injdur_sai + 365 + is_bissext ENDDO ! ... then subtracting part of the first year where no injection yet... is_bissext = (MOD(year_emit_sai_start, 4)==0) .AND. (MOD(year_emit_sai_start, 100) /= 0 .OR. MOD(year_emit_sai_start, 400) == 0) SELECT CASE(mth_emit_sai_start) CASE(2) injdur_sai = injdur_sai - 31 CASE(3) injdur_sai = injdur_sai - 31 - 28 - is_bissext CASE(4) injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 CASE(5) injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30 CASE(6) injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30 - 31 CASE(7) injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30 - 31 - 30 CASE(8) injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30 - 31 - 30 - 31 CASE(9) injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30 - 31 - 30 - 31 - 31 CASE(10) injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30 - 31 - 30 - 31 - 31 - 30 CASE(11) injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30 - 31 - 30 - 31 - 31 - 30 - 31 CASE(12) injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30 - 31 - 30 - 31 - 31 - 30 - 31 - 30 END SELECT injdur_sai = injdur_sai - day_emit_sai_start + 1 ! ... then adding part of the n-th year is_bissext = (MOD(year_emit_sai_end, 4)==0) .AND. (MOD(year_emit_sai_end, 100) /= 0 .OR. MOD(year_emit_sai_end, 400) == 0) SELECT CASE(mth_emit_sai_end) CASE(2) injdur_sai = injdur_sai + 31 CASE(3) injdur_sai = injdur_sai + 31 + 28 + is_bissext CASE(4) injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 CASE(5) injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30 CASE(6) injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30 + 31 CASE(7) injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30 + 31 + 30 CASE(8) injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30 + 31 + 30 + 31 CASE(9) injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30 + 31 + 30 + 31 + 31 CASE(10) injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30 + 31 + 30 + 31 + 31 + 30 CASE(11) injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30 + 31 + 30 + 31 + 31 + 30 + 31 CASE(12) injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30 + 31 + 30 + 31 + 31 + 30 + 31 + 30 END SELECT injdur_sai = injdur_sai + day_emit_sai_end ! A security: are SAI dates of injection consistent? IF (injdur_sai <= 0) THEN CALL abort_physic('traccoag_mod', 'Pb in SAI dates of injection.', 1) ENDIF ! Injection in itself IF ((year_emit_sai_start <= year_cur) & .AND. (year_cur <= year_emit_sai_end) & .AND. (mth_emit_sai_start <= mth_cur .OR. year_emit_sai_start < year_cur) & .AND. (mth_cur <= mth_emit_sai_end .OR. year_cur < year_emit_sai_end) & .AND. (day_emit_sai_start <= day_cur .OR. mth_emit_sai_start < mth_cur .OR. year_emit_sai_start < year_cur) & .AND. (day_cur <= day_emit_sai_end .OR. mth_cur < mth_emit_sai_end .OR. year_cur < year_emit_sai_end)) THEN m_aer = m_aer_emiss_sai !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss) m_aer = m_aer * (mSO2mol / mSatom) latmin = xlat_sai latmax = xlat_sai lonmin = xlon_sai lonmax = xlon_sai altemiss = altemiss_sai sigma_alt = sigma_alt_sai pdt = 0. ! stretch emission over whole year (360d) stretchlong = FLOAT(year_len) CALL STRATEMIT(pdtphys, pdt, xlat, xlon, t_seri, pplay, paprs, m_air_gridbox, tr_seri, & m_aer, latmin, latmax, lonmin, lonmax, altemiss, sigma_alt, id_SO2_strat, & stretchlong, 1, 0) budg_emi_so2(:) = budg_emi(:, 1) * mSatom / mSO2mol ENDIF ! Condition over injection dates CASE(3) ! --- SAI injection over a single band of longitude and between ! lat_min and lat_max m_aer = m_aer_emiss_sai !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss) m_aer = m_aer * (mSO2mol / mSatom) latmin = xlat_min_sai latmax = xlat_max_sai lonmin = xlon_sai lonmax = xlon_sai altemiss = altemiss_sai sigma_alt = sigma_alt_sai pdt = 0. ! stretch emission over whole year (360d) stretchlong = FLOAT(year_len) CALL STRATEMIT(pdtphys, pdt, xlat, xlon, t_seri, pplay, paprs, m_air_gridbox, tr_seri, & m_aer, latmin, latmax, lonmin, lonmax, altemiss, sigma_alt, id_SO2_strat, & stretchlong, 1, 0) budg_emi_so2(:) = budg_emi(:, 1) * mSatom / mSO2mol END SELECT ! emission scenario (flag_emit) !--read background concentrations of OCS and SO2 and lifetimes from input file !--update the variables defined in phys_local_var_mod CALL interp_sulf_input(debutphy, pdtphys, paprs, tr_seri) !--convert OCS to SO2 in the stratosphere CALL ocs_to_so2(pdtphys, tr_seri, t_seri, pplay, paprs, is_strato) !--convert SO2 to H2SO4 CALL so2_to_h2so4(pdtphys, tr_seri, t_seri, pplay, paprs, is_strato) !--common routine for nucleation and condensation/evaporation with adaptive timestep CALL micphy_tstep(pdtphys, tr_seri, t_seri, pplay, paprs, rh, is_strato) !--CALL coagulation routine CALL coagulate(pdtphys, mdw, tr_seri, t_seri, pplay, dens_aer, is_strato) !--CALL sedimentation routine CALL aer_sedimnt(pdtphys, t_seri, pplay, paprs, tr_seri, dens_aer) !--compute mass concentration of PM2.5 sulfate particles (wet diameter and mass) at the surface for health studies surf_PM25_sulf(:) = 0.0 DO i = 1, klon DO it = 1, nbtr_bin IF (mdw(it) < 2.5e-6) THEN !surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas)*m_part(i,1,it) & !assume that particles consist of ammonium sulfate at the surface (132g/mol) !and are dry at T = 20 deg. C and 50 perc. humidity surf_PM25_sulf(i) = surf_PM25_sulf(i) + tr_seri(i, 1, it + nbtr_sulgas) & * 132. / 98. * dens_aer_dry * 4. / 3. * RPI * (mdw(it) / 2.)**3 & * pplay(i, 1) / t_seri(i, 1) / RD * 1.e9 ENDIF ENDDO ENDDO !--compute ! sulfmmr: Sulfate aerosol concentration (dry mixing ratio) (condensed H2SO4 mmr) ! SAD_sulfate: SAD all aerosols (cm2/cm3) (must be WET) ! sulfmmr_mode: sulfate(=H2SO4 if dry) MMR in different modes (ambiguous but based on sulfmmr, it mus be DRY(?) mmr) ! nd_mode: DRY(?) particle concentration in different modes (part/m3) sulfmmr(:, :) = 0.0 SAD_sulfate(:, :) = 0.0 sulfmmr_mode(:, :, :) = 0.0 nd_mode(:, :, :) = 0.0 DO i = 1, klon DO j = 1, klev DO it = 1, nbtr_bin !surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas)*m_part(i,1,it) & !assume that particles consist of ammonium sulfate at the surface (132g/mol) !and are dry at T = 20 deg. C and 50 perc. humidity ! sulfmmr_mode: sulfate(=H2SO4 if dry) MMR in different modes (based on sulfmmr, it must be DRY mmr) ! equivalent to condensed H2SO4 mmr= H2SO4 kg / kgA in bin it sulfmmr_mode(i, j, it) = tr_seri(i, j, it + nbtr_sulgas) & ! [DRY part/kgA in bin it] * (4. / 3.) * RPI * (mdw(it) / 2.)**3. & ! [mdw: dry diameter in m] * dens_aer_dry ! [dry aerosol mass density in kg/m3] ! sulfmmr: Sulfate aerosol concentration (dry mass mixing ratio) ! equivalent to total condensed H2SO4 mmr (H2SO4 kg / kgA sulfmmr(i, j) = sulfmmr(i, j) + sulfmmr_mode(i, j, it) ! nd_mode: particle concentration in different modes (DRY part/m3) nd_mode(i, j, it) = tr_seri(i, j, it + nbtr_sulgas) & ! [DRY part/kgA in bin it] * pplay(i, j) / t_seri(i, j) / RD ! [air mass concentration in kg air /m3A] IF(flag_new_strat_compo) THEN ! SAD_sulfate: SAD WET sulfate aerosols (cm2/cm3) SAD_sulfate(i, j) = SAD_sulfate(i, j) + nd_mode(i, j, it) & ! [DRY part/m3A (in bin it)] * 4. * RPI * (mdw(it) * f_r_wetB(i, j, it) / 2.)**2. & ! [WET SA of part it in m2] * 1.e-2 ! conversion from m2/m3 to cm2/cm3A ELSE ! SAD_sulfate: SAD WET sulfate aerosols (cm2/cm3) SAD_sulfate(i, j) = SAD_sulfate(i, j) + nd_mode(i, j, it) & ! [DRY part/m3A (in bin it)] * 4. * RPI * (mdw(it) * f_r_wet(i, j) / 2.)**2. & ! [WET SA of part it in m2] * 1.e-2 ! conversion from m2/m3 to cm2/cm3A ENDIF ENDDO ENDDO ENDDO END SUBROUTINE traccoag END MODULE traccoag_mod