MODULE traccoag_mod ! ! This module calculates the concentration of aerosol particles in certain size bins ! considering coagulation and sedimentation. ! IMPLICIT NONE; PRIVATE PUBLIC traccoag 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, reff_sulfate USE dimphy USE infotrac_phy, ONLY : nbtr_bin, nbtr_sulgas, nbtr, id_SO2_strat USE aerophys USE geometry_mod, ONLY : cell_area, boundslat USE mod_grid_phy_lmdz USE mod_phys_lmdz_mpi_data, ONLY : is_mpi_root USE mod_phys_lmdz_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 yomcst_mod_h USE print_control_mod, 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 REAL :: samoment2, samoment3! 2nd and 3rd order moments of size distribution 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.EQ.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.EQ.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.GT.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) .LT. 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 reff_sulfate(:,:)=0.0 DO i=1,klon DO j=1,klev samoment2=0.0 samoment3=0.0 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 ! samoment2 : 2nd order moment of WET sulfate aerosols (m2/m3) samoment2 = samoment2 + nd_mode(i,j,it) & ! [DRY part/m3A (in bin it)] & *( mdw(it)*f_r_wetB(i,j,it)/2. )**2. ! [WET SA of part it in m2] ! samoment3 : 3nd order moment of WET sulfate aerosols (cm2/cm3) samoment3 = samoment3 + nd_mode(i,j,it) & ! [DRY part/m3A (in bin it)] & *( mdw(it)*f_r_wetB(i,j,it)/2. )**3. ! [WET SA of part it in m2] 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 ! samoment2 : 2nd order moment of WET sulfate aerosols (m2/m3) samoment2 = samoment2 + nd_mode(i,j,it) & ! [DRY part/m3A (in bin it)] & *( mdw(it)*f_r_wet(i,j)/2. )**2. ! [WET SA of part it in m2] ! samoment3 : 3nd order moment of WET sulfate aerosols (cm2/cm3) samoment3 = samoment3 + nd_mode(i,j,it) & ! [DRY part/m3A (in bin it)] & *( mdw(it)*f_r_wet(i,j)/2. )**3. ! [WET SA of part it in m2] ENDIF ENDDO ! reff_sulfate: effective radius of WET sulfate aerosols (cm) IF(samoment2 > 1.e-30) THEN reff_sulfate(i,j) = (samoment3 / samoment2) & & *1.e2 ! conversion from m to cm ! Sanity check IF(reff_sulfate(i,j) > 5.e-4) reff_sulfate(i,j) = 5.e-4 ! reff_sulfate max = 5 micron IF(reff_sulfate(i,j) < 1.e-6) reff_sulfate(i,j) = 1.e-6 ! reff_sulfate min = 10 nm ELSE reff_sulfate(i,j) = 1.e-5 ! reff_sulfate N part = nul (ref = 100 nm) ENDIF ENDDO ENDDO END SUBROUTINE traccoag END MODULE traccoag_mod