| 1 | MODULE traccoag_mod |
|---|
| 2 | |
|---|
| 3 | ! This module calculates the concentration of aerosol particles in certain size bins |
|---|
| 4 | ! considering coagulation and sedimentation. |
|---|
| 5 | USE lmdz_abort_physic, ONLY: abort_physic |
|---|
| 6 | CONTAINS |
|---|
| 7 | |
|---|
| 8 | SUBROUTINE traccoag(pdtphys, gmtime, debutphy, julien, & |
|---|
| 9 | presnivs, xlat, xlon, pphis, pphi, & |
|---|
| 10 | t_seri, pplay, paprs, sh, rh, tr_seri) |
|---|
| 11 | |
|---|
| 12 | USE phys_local_var_mod, ONLY: mdw, R2SO4, DENSO4, f_r_wet, surf_PM25_sulf, & |
|---|
| 13 | budg_emi_ocs, budg_emi_so2, budg_emi_h2so4, budg_emi_part, & |
|---|
| 14 | R2SO4B, DENSO4B, f_r_wetB, sulfmmr, SAD_sulfate, sulfmmr_mode, nd_mode |
|---|
| 15 | |
|---|
| 16 | USE dimphy |
|---|
| 17 | USE infotrac_phy, ONLY: nbtr_bin, nbtr_sulgas, nbtr, id_SO2_strat |
|---|
| 18 | USE aerophys |
|---|
| 19 | USE lmdz_geometry, ONLY: cell_area, boundslat |
|---|
| 20 | USE lmdz_grid_phy |
|---|
| 21 | USE lmdz_phys_mpi_data, ONLY: is_mpi_root |
|---|
| 22 | USE lmdz_phys_para, ONLY: gather, scatter |
|---|
| 23 | USE phys_cal_mod, ONLY: year_len, year_cur, mth_cur, day_cur, hour |
|---|
| 24 | USE sulfate_aer_mod |
|---|
| 25 | USE phys_local_var_mod, ONLY: stratomask |
|---|
| 26 | USE lmdz_yomcst |
|---|
| 27 | USE lmdz_print_control, ONLY: lunout |
|---|
| 28 | USE strataer_local_var_mod |
|---|
| 29 | |
|---|
| 30 | IMPLICIT NONE |
|---|
| 31 | |
|---|
| 32 | ! Input argument |
|---|
| 33 | !--------------- |
|---|
| 34 | REAL, INTENT(IN) :: pdtphys ! Pas d'integration pour la physique (seconde) |
|---|
| 35 | REAL, INTENT(IN) :: gmtime ! Heure courante |
|---|
| 36 | LOGICAL, INTENT(IN) :: debutphy ! le flag de l'initialisation de la physique |
|---|
| 37 | INTEGER, INTENT(IN) :: julien ! Jour julien |
|---|
| 38 | |
|---|
| 39 | REAL, DIMENSION(klev), INTENT(IN) :: presnivs! pressions approximat. des milieux couches (en PA) |
|---|
| 40 | REAL, DIMENSION(klon), INTENT(IN) :: xlat ! latitudes pour chaque point |
|---|
| 41 | REAL, DIMENSION(klon), INTENT(IN) :: xlon ! longitudes pour chaque point |
|---|
| 42 | REAL, DIMENSION(klon), INTENT(IN) :: pphis ! geopotentiel du sol |
|---|
| 43 | REAL, DIMENSION(klon, klev), INTENT(IN) :: pphi ! geopotentiel de chaque couche |
|---|
| 44 | |
|---|
| 45 | REAL, DIMENSION(klon, klev), INTENT(IN) :: t_seri ! Temperature |
|---|
| 46 | REAL, DIMENSION(klon, klev), INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa) |
|---|
| 47 | REAL, DIMENSION(klon, klev + 1), INTENT(IN) :: paprs ! pression pour chaque inter-couche (en Pa) |
|---|
| 48 | REAL, DIMENSION(klon, klev), INTENT(IN) :: sh ! humidite specifique |
|---|
| 49 | REAL, DIMENSION(klon, klev), INTENT(IN) :: rh ! humidite relative |
|---|
| 50 | |
|---|
| 51 | ! Output argument |
|---|
| 52 | !---------------- |
|---|
| 53 | REAL, DIMENSION(klon, klev, nbtr), INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA] |
|---|
| 54 | |
|---|
| 55 | ! Local variables |
|---|
| 56 | !---------------- |
|---|
| 57 | REAL :: m_aer_emiss_vol_daily ! daily injection mass emission |
|---|
| 58 | REAL :: m_aer ! aerosol mass |
|---|
| 59 | INTEGER :: it, k, i, j, ilon, ilev, itime, i_int, ieru |
|---|
| 60 | LOGICAL, DIMENSION(klon, klev) :: is_strato ! true = above tropopause, false = below |
|---|
| 61 | REAL, DIMENSION(klon, klev) :: m_air_gridbox ! mass of air in every grid box [kg] |
|---|
| 62 | REAL, DIMENSION(klon_glo, klev, nbtr) :: tr_seri_glo ! Concentration Traceur [U/KgA] |
|---|
| 63 | REAL, DIMENSION(klev + 1) :: altLMDz ! altitude of layer interfaces in m |
|---|
| 64 | REAL, DIMENSION(klev) :: f_lay_emiss ! fraction of emission for every vertical layer |
|---|
| 65 | REAL :: f_lay_sum ! sum of layer emission fractions |
|---|
| 66 | REAL :: alt ! altitude for integral calculation |
|---|
| 67 | INTEGER, PARAMETER :: n_int_alt = 10 ! number of subintervals for integration over Gaussian emission profile |
|---|
| 68 | REAL, DIMENSION(nbtr_bin) :: r_bin ! particle radius in size bin [m] |
|---|
| 69 | REAL, DIMENSION(nbtr_bin) :: r_lower ! particle radius at lower bin boundary [m] |
|---|
| 70 | REAL, DIMENSION(nbtr_bin) :: r_upper ! particle radius at upper bin boundary [m] |
|---|
| 71 | REAL, DIMENSION(nbtr_bin) :: m_part_dry ! mass of one dry particle in size bin [kg] |
|---|
| 72 | REAL :: zrho ! Density of air [kg/m3] |
|---|
| 73 | REAL :: zdz ! thickness of atm. model layer in m |
|---|
| 74 | REAL, DIMENSION(klev) :: zdm ! mass of atm. model layer in kg |
|---|
| 75 | REAL, DIMENSION(klon, klev) :: dens_aer ! density of aerosol particles [kg/m3 aerosol] with default H2SO4 mass fraction |
|---|
| 76 | REAL :: emission ! emission |
|---|
| 77 | REAL :: theta_min, theta_max ! for SAI computation between two latitudes |
|---|
| 78 | REAL :: dlat_loc |
|---|
| 79 | REAL :: latmin, latmax, lonmin, lonmax ! lat/lon min/max for injection |
|---|
| 80 | REAL :: sigma_alt, altemiss ! injection altitude + sigma for distrib |
|---|
| 81 | REAL :: pdt, stretchlong ! physic timestep, stretch emission over one day |
|---|
| 82 | |
|---|
| 83 | INTEGER :: injdur_sai ! injection duration for SAI case [days] |
|---|
| 84 | INTEGER :: yr, is_bissext |
|---|
| 85 | |
|---|
| 86 | IF (is_mpi_root .AND. flag_verbose_strataer) THEN |
|---|
| 87 | WRITE(lunout, *) 'in traccoag: date from phys_cal_mod =', year_cur, '-', mth_cur, '-', day_cur, '-', hour |
|---|
| 88 | WRITE(lunout, *) 'IN traccoag flag_emit: ', flag_emit |
|---|
| 89 | ENDIF |
|---|
| 90 | |
|---|
| 91 | ! radius [m] |
|---|
| 92 | DO it = 1, nbtr_bin |
|---|
| 93 | r_bin(it) = mdw(it) / 2. |
|---|
| 94 | ENDDO |
|---|
| 95 | |
|---|
| 96 | !--set boundaries of size bins |
|---|
| 97 | DO it = 1, nbtr_bin |
|---|
| 98 | IF (it==1) THEN |
|---|
| 99 | r_upper(it) = sqrt(r_bin(it + 1) * r_bin(it)) |
|---|
| 100 | r_lower(it) = r_bin(it)**2. / r_upper(it) |
|---|
| 101 | ELSEIF (it==nbtr_bin) THEN |
|---|
| 102 | r_lower(it) = sqrt(r_bin(it) * r_bin(it - 1)) |
|---|
| 103 | r_upper(it) = r_bin(it)**2. / r_lower(it) |
|---|
| 104 | ELSE |
|---|
| 105 | r_lower(it) = sqrt(r_bin(it) * r_bin(it - 1)) |
|---|
| 106 | r_upper(it) = sqrt(r_bin(it + 1) * r_bin(it)) |
|---|
| 107 | ENDIF |
|---|
| 108 | ENDDO |
|---|
| 109 | |
|---|
| 110 | IF (debutphy .AND. is_mpi_root) THEN |
|---|
| 111 | DO it = 1, nbtr_bin |
|---|
| 112 | WRITE(lunout, *) 'radius bin', it, ':', r_bin(it), '(from', r_lower(it), 'to', r_upper(it), ')' |
|---|
| 113 | ENDDO |
|---|
| 114 | ENDIF |
|---|
| 115 | |
|---|
| 116 | !--initialising LOGICAL is_strato from stratomask |
|---|
| 117 | is_strato(:, :) = .FALSE. |
|---|
| 118 | WHERE (stratomask>0.5) is_strato = .TRUE. |
|---|
| 119 | |
|---|
| 120 | IF(flag_new_strat_compo) THEN |
|---|
| 121 | IF(debutphy) WRITE(lunout, *) 'traccoag: COMPO/DENSITY (Tabazadeh 97) + H2O kelvin effect', flag_new_strat_compo |
|---|
| 122 | ! STRACOMP (H2O, P, t_seri, R -> R2SO4 + Kelvin effect) : Taba97, Socol, etc... |
|---|
| 123 | CALL stracomp_kelvin(sh, t_seri, pplay) |
|---|
| 124 | ELSE |
|---|
| 125 | IF(debutphy) WRITE(lunout, *) 'traccoag: COMPO from Bekki 2D model', flag_new_strat_compo |
|---|
| 126 | ! STRACOMP (H2O, P, t_seri -> aerosol composition (R2SO4)) |
|---|
| 127 | ! H2SO4 mass fraction in aerosol (%) |
|---|
| 128 | CALL stracomp(sh, t_seri, pplay) |
|---|
| 129 | |
|---|
| 130 | ! aerosol density (gr/cm3) |
|---|
| 131 | CALL denh2sa(t_seri) |
|---|
| 132 | |
|---|
| 133 | ! compute factor for converting dry to wet radius (for every grid box) |
|---|
| 134 | f_r_wet(:, :) = (dens_aer_dry / (DENSO4(:, :) * 1000.) / (R2SO4(:, :) / 100.))**(1. / 3.) |
|---|
| 135 | ENDIF |
|---|
| 136 | |
|---|
| 137 | !--calculate mass of air in every grid box |
|---|
| 138 | DO ilon = 1, klon |
|---|
| 139 | DO ilev = 1, klev |
|---|
| 140 | m_air_gridbox(ilon, ilev) = (paprs(ilon, ilev) - paprs(ilon, ilev + 1)) / RG * cell_area(ilon) |
|---|
| 141 | ENDDO |
|---|
| 142 | ENDDO |
|---|
| 143 | |
|---|
| 144 | !--initialise emission diagnostics |
|---|
| 145 | IF (nErupt > 0 .AND. (flag_emit == 1 .OR. flag_emit == 4)) budg_emi(:, 1) = 0.0 |
|---|
| 146 | budg_emi_ocs(:) = 0.0 |
|---|
| 147 | budg_emi_so2(:) = 0.0 |
|---|
| 148 | budg_emi_h2so4(:) = 0.0 |
|---|
| 149 | budg_emi_part(:) = 0.0 |
|---|
| 150 | |
|---|
| 151 | !--sulfur emission, depending on chosen scenario (flag_emit) |
|---|
| 152 | SELECT CASE(flag_emit) |
|---|
| 153 | |
|---|
| 154 | CASE(0) ! background aerosol |
|---|
| 155 | ! do nothing (no emission) |
|---|
| 156 | |
|---|
| 157 | CASE(1) ! volcanic eruption |
|---|
| 158 | !--only emit on day of eruption |
|---|
| 159 | ! stretch emission over one day of Pinatubo eruption |
|---|
| 160 | DO ieru = 1, nErupt |
|---|
| 161 | IF (year_cur==year_emit_vol(ieru).AND.mth_cur==mth_emit_vol(ieru).AND.& |
|---|
| 162 | day_cur>=day_emit_vol(ieru).AND.day_cur<(day_emit_vol(ieru) + injdur)) THEN |
|---|
| 163 | |
|---|
| 164 | ! daily injection mass emission |
|---|
| 165 | m_aer = m_aer_emiss_vol(ieru, 1) / (REAL(injdur) * REAL(ponde_lonlat_vol(ieru))) |
|---|
| 166 | !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss) |
|---|
| 167 | m_aer = m_aer * (mSO2mol / mSatom) |
|---|
| 168 | |
|---|
| 169 | WRITE(lunout, *) 'IN traccoag m_aer_emiss_vol(ieru)=', m_aer_emiss_vol(ieru, 1), & |
|---|
| 170 | 'ponde_lonlat_vol(ieru)=', ponde_lonlat_vol(ieru), '(injdur*ponde_lonlat_vol(ieru))', & |
|---|
| 171 | (injdur * ponde_lonlat_vol(ieru)), 'm_aer_emiss_vol_daily=', m_aer, 'ieru=', ieru |
|---|
| 172 | WRITE(lunout, *) 'IN traccoag, dlon=', dlon |
|---|
| 173 | |
|---|
| 174 | latmin = xlat_min_vol(ieru) |
|---|
| 175 | latmax = xlat_max_vol(ieru) |
|---|
| 176 | lonmin = xlon_min_vol(ieru) |
|---|
| 177 | lonmax = xlon_max_vol(ieru) |
|---|
| 178 | altemiss = altemiss_vol(ieru) |
|---|
| 179 | sigma_alt = sigma_alt_vol(ieru) |
|---|
| 180 | pdt = pdtphys |
|---|
| 181 | ! stretch emission over one day of eruption |
|---|
| 182 | stretchlong = 1. |
|---|
| 183 | |
|---|
| 184 | CALL STRATEMIT(pdtphys, pdt, xlat, xlon, t_seri, pplay, paprs, tr_seri, & |
|---|
| 185 | m_aer, latmin, latmax, lonmin, lonmax, altemiss, sigma_alt, id_SO2_strat, & |
|---|
| 186 | stretchlong, 1, 0) |
|---|
| 187 | |
|---|
| 188 | ENDIF ! emission period |
|---|
| 189 | ENDDO ! eruption number |
|---|
| 190 | |
|---|
| 191 | CASE(2) ! stratospheric aerosol injections (SAI) |
|---|
| 192 | |
|---|
| 193 | ! Computing duration of SAI in days... |
|---|
| 194 | ! ... starting from 0... |
|---|
| 195 | injdur_sai = 0 |
|---|
| 196 | ! ... then adding whole years from first to (n-1)th... |
|---|
| 197 | DO yr = year_emit_sai_start, year_emit_sai_end - 1 |
|---|
| 198 | ! (n % 4 == 0) and (n % 100 != 0 or n % 400 == 0) |
|---|
| 199 | is_bissext = (MOD(yr, 4)==0) .AND. (MOD(yr, 100) /= 0 .OR. MOD(yr, 400) == 0) |
|---|
| 200 | injdur_sai = injdur_sai + 365 + is_bissext |
|---|
| 201 | ENDDO |
|---|
| 202 | ! ... then subtracting part of the first year where no injection yet... |
|---|
| 203 | 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) |
|---|
| 204 | SELECT CASE(mth_emit_sai_start) |
|---|
| 205 | CASE(2) |
|---|
| 206 | injdur_sai = injdur_sai - 31 |
|---|
| 207 | CASE(3) |
|---|
| 208 | injdur_sai = injdur_sai - 31 - 28 - is_bissext |
|---|
| 209 | CASE(4) |
|---|
| 210 | injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 |
|---|
| 211 | CASE(5) |
|---|
| 212 | injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30 |
|---|
| 213 | CASE(6) |
|---|
| 214 | injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30 - 31 |
|---|
| 215 | CASE(7) |
|---|
| 216 | injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30 - 31 - 30 |
|---|
| 217 | CASE(8) |
|---|
| 218 | injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30 - 31 - 30 - 31 |
|---|
| 219 | CASE(9) |
|---|
| 220 | injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30 - 31 - 30 - 31 - 31 |
|---|
| 221 | CASE(10) |
|---|
| 222 | injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30 - 31 - 30 - 31 - 31 - 30 |
|---|
| 223 | CASE(11) |
|---|
| 224 | injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30 - 31 - 30 - 31 - 31 - 30 - 31 |
|---|
| 225 | CASE(12) |
|---|
| 226 | injdur_sai = injdur_sai - 31 - 28 - is_bissext - 31 - 30 - 31 - 30 - 31 - 31 - 30 - 31 - 30 |
|---|
| 227 | END SELECT |
|---|
| 228 | injdur_sai = injdur_sai - day_emit_sai_start + 1 |
|---|
| 229 | ! ... then adding part of the n-th year |
|---|
| 230 | 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) |
|---|
| 231 | SELECT CASE(mth_emit_sai_end) |
|---|
| 232 | CASE(2) |
|---|
| 233 | injdur_sai = injdur_sai + 31 |
|---|
| 234 | CASE(3) |
|---|
| 235 | injdur_sai = injdur_sai + 31 + 28 + is_bissext |
|---|
| 236 | CASE(4) |
|---|
| 237 | injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 |
|---|
| 238 | CASE(5) |
|---|
| 239 | injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30 |
|---|
| 240 | CASE(6) |
|---|
| 241 | injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30 + 31 |
|---|
| 242 | CASE(7) |
|---|
| 243 | injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30 + 31 + 30 |
|---|
| 244 | CASE(8) |
|---|
| 245 | injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30 + 31 + 30 + 31 |
|---|
| 246 | CASE(9) |
|---|
| 247 | injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30 + 31 + 30 + 31 + 31 |
|---|
| 248 | CASE(10) |
|---|
| 249 | injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30 + 31 + 30 + 31 + 31 + 30 |
|---|
| 250 | CASE(11) |
|---|
| 251 | injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30 + 31 + 30 + 31 + 31 + 30 + 31 |
|---|
| 252 | CASE(12) |
|---|
| 253 | injdur_sai = injdur_sai + 31 + 28 + is_bissext + 31 + 30 + 31 + 30 + 31 + 31 + 30 + 31 + 30 |
|---|
| 254 | END SELECT |
|---|
| 255 | injdur_sai = injdur_sai + day_emit_sai_end |
|---|
| 256 | ! A security: are SAI dates of injection consistent? |
|---|
| 257 | IF (injdur_sai <= 0) THEN |
|---|
| 258 | CALL abort_physic('traccoag_mod', 'Pb in SAI dates of injection.', 1) |
|---|
| 259 | ENDIF |
|---|
| 260 | ! Injection in itself |
|---|
| 261 | IF ((year_emit_sai_start <= year_cur) & |
|---|
| 262 | .AND. (year_cur <= year_emit_sai_end) & |
|---|
| 263 | .AND. (mth_emit_sai_start <= mth_cur .OR. year_emit_sai_start < year_cur) & |
|---|
| 264 | .AND. (mth_cur <= mth_emit_sai_end .OR. year_cur < year_emit_sai_end) & |
|---|
| 265 | .AND. (day_emit_sai_start <= day_cur .OR. mth_emit_sai_start < mth_cur .OR. year_emit_sai_start < year_cur) & |
|---|
| 266 | .AND. (day_cur <= day_emit_sai_end .OR. mth_cur < mth_emit_sai_end .OR. year_cur < year_emit_sai_end)) THEN |
|---|
| 267 | |
|---|
| 268 | m_aer = m_aer_emiss_sai |
|---|
| 269 | !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss) |
|---|
| 270 | m_aer = m_aer * (mSO2mol / mSatom) |
|---|
| 271 | |
|---|
| 272 | latmin = xlat_sai |
|---|
| 273 | latmax = xlat_sai |
|---|
| 274 | lonmin = xlon_sai |
|---|
| 275 | lonmax = xlon_sai |
|---|
| 276 | altemiss = altemiss_sai |
|---|
| 277 | sigma_alt = sigma_alt_sai |
|---|
| 278 | pdt = 0. |
|---|
| 279 | ! stretch emission over whole year (360d) |
|---|
| 280 | stretchlong = FLOAT(year_len) |
|---|
| 281 | |
|---|
| 282 | CALL STRATEMIT(pdtphys, pdt, xlat, xlon, t_seri, pplay, paprs, m_air_gridbox, tr_seri, & |
|---|
| 283 | m_aer, latmin, latmax, lonmin, lonmax, altemiss, sigma_alt, id_SO2_strat, & |
|---|
| 284 | stretchlong, 1, 0) |
|---|
| 285 | |
|---|
| 286 | budg_emi_so2(:) = budg_emi(:, 1) * mSatom / mSO2mol |
|---|
| 287 | ENDIF ! Condition over injection dates |
|---|
| 288 | |
|---|
| 289 | CASE(3) ! --- SAI injection over a single band of longitude and between |
|---|
| 290 | ! lat_min and lat_max |
|---|
| 291 | |
|---|
| 292 | m_aer = m_aer_emiss_sai |
|---|
| 293 | !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss) |
|---|
| 294 | m_aer = m_aer * (mSO2mol / mSatom) |
|---|
| 295 | |
|---|
| 296 | latmin = xlat_min_sai |
|---|
| 297 | latmax = xlat_max_sai |
|---|
| 298 | lonmin = xlon_sai |
|---|
| 299 | lonmax = xlon_sai |
|---|
| 300 | altemiss = altemiss_sai |
|---|
| 301 | sigma_alt = sigma_alt_sai |
|---|
| 302 | pdt = 0. |
|---|
| 303 | ! stretch emission over whole year (360d) |
|---|
| 304 | stretchlong = FLOAT(year_len) |
|---|
| 305 | |
|---|
| 306 | CALL STRATEMIT(pdtphys, pdt, xlat, xlon, t_seri, pplay, paprs, m_air_gridbox, tr_seri, & |
|---|
| 307 | m_aer, latmin, latmax, lonmin, lonmax, altemiss, sigma_alt, id_SO2_strat, & |
|---|
| 308 | stretchlong, 1, 0) |
|---|
| 309 | |
|---|
| 310 | budg_emi_so2(:) = budg_emi(:, 1) * mSatom / mSO2mol |
|---|
| 311 | |
|---|
| 312 | END SELECT ! emission scenario (flag_emit) |
|---|
| 313 | |
|---|
| 314 | !--read background concentrations of OCS and SO2 and lifetimes from input file |
|---|
| 315 | !--update the variables defined in phys_local_var_mod |
|---|
| 316 | CALL interp_sulf_input(debutphy, pdtphys, paprs, tr_seri) |
|---|
| 317 | |
|---|
| 318 | !--convert OCS to SO2 in the stratosphere |
|---|
| 319 | CALL ocs_to_so2(pdtphys, tr_seri, t_seri, pplay, paprs, is_strato) |
|---|
| 320 | |
|---|
| 321 | !--convert SO2 to H2SO4 |
|---|
| 322 | CALL so2_to_h2so4(pdtphys, tr_seri, t_seri, pplay, paprs, is_strato) |
|---|
| 323 | |
|---|
| 324 | !--common routine for nucleation and condensation/evaporation with adaptive timestep |
|---|
| 325 | CALL micphy_tstep(pdtphys, tr_seri, t_seri, pplay, paprs, rh, is_strato) |
|---|
| 326 | |
|---|
| 327 | !--CALL coagulation routine |
|---|
| 328 | CALL coagulate(pdtphys, mdw, tr_seri, t_seri, pplay, dens_aer, is_strato) |
|---|
| 329 | |
|---|
| 330 | !--CALL sedimentation routine |
|---|
| 331 | CALL aer_sedimnt(pdtphys, t_seri, pplay, paprs, tr_seri, dens_aer) |
|---|
| 332 | |
|---|
| 333 | !--compute mass concentration of PM2.5 sulfate particles (wet diameter and mass) at the surface for health studies |
|---|
| 334 | surf_PM25_sulf(:) = 0.0 |
|---|
| 335 | DO i = 1, klon |
|---|
| 336 | DO it = 1, nbtr_bin |
|---|
| 337 | IF (mdw(it) < 2.5e-6) THEN |
|---|
| 338 | !surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas)*m_part(i,1,it) & |
|---|
| 339 | !assume that particles consist of ammonium sulfate at the surface (132g/mol) |
|---|
| 340 | !and are dry at T = 20 deg. C and 50 perc. humidity |
|---|
| 341 | surf_PM25_sulf(i) = surf_PM25_sulf(i) + tr_seri(i, 1, it + nbtr_sulgas) & |
|---|
| 342 | * 132. / 98. * dens_aer_dry * 4. / 3. * RPI * (mdw(it) / 2.)**3 & |
|---|
| 343 | * pplay(i, 1) / t_seri(i, 1) / RD * 1.e9 |
|---|
| 344 | ENDIF |
|---|
| 345 | ENDDO |
|---|
| 346 | ENDDO |
|---|
| 347 | |
|---|
| 348 | !--compute |
|---|
| 349 | ! sulfmmr: Sulfate aerosol concentration (dry mixing ratio) (condensed H2SO4 mmr) |
|---|
| 350 | ! SAD_sulfate: SAD all aerosols (cm2/cm3) (must be WET) |
|---|
| 351 | ! sulfmmr_mode: sulfate(=H2SO4 if dry) MMR in different modes (ambiguous but based on sulfmmr, it mus be DRY(?) mmr) |
|---|
| 352 | ! nd_mode: DRY(?) particle concentration in different modes (part/m3) |
|---|
| 353 | sulfmmr(:, :) = 0.0 |
|---|
| 354 | SAD_sulfate(:, :) = 0.0 |
|---|
| 355 | sulfmmr_mode(:, :, :) = 0.0 |
|---|
| 356 | nd_mode(:, :, :) = 0.0 |
|---|
| 357 | |
|---|
| 358 | DO i = 1, klon |
|---|
| 359 | DO j = 1, klev |
|---|
| 360 | DO it = 1, nbtr_bin |
|---|
| 361 | !surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas)*m_part(i,1,it) & |
|---|
| 362 | !assume that particles consist of ammonium sulfate at the surface (132g/mol) |
|---|
| 363 | !and are dry at T = 20 deg. C and 50 perc. humidity |
|---|
| 364 | |
|---|
| 365 | ! sulfmmr_mode: sulfate(=H2SO4 if dry) MMR in different modes (based on sulfmmr, it must be DRY mmr) |
|---|
| 366 | ! equivalent to condensed H2SO4 mmr= H2SO4 kg / kgA in bin it |
|---|
| 367 | sulfmmr_mode(i, j, it) = tr_seri(i, j, it + nbtr_sulgas) & ! [DRY part/kgA in bin it] |
|---|
| 368 | * (4. / 3.) * RPI * (mdw(it) / 2.)**3. & ! [mdw: dry diameter in m] |
|---|
| 369 | * dens_aer_dry ! [dry aerosol mass density in kg/m3] |
|---|
| 370 | |
|---|
| 371 | ! sulfmmr: Sulfate aerosol concentration (dry mass mixing ratio) |
|---|
| 372 | ! equivalent to total condensed H2SO4 mmr (H2SO4 kg / kgA |
|---|
| 373 | sulfmmr(i, j) = sulfmmr(i, j) + sulfmmr_mode(i, j, it) |
|---|
| 374 | |
|---|
| 375 | ! nd_mode: particle concentration in different modes (DRY part/m3) |
|---|
| 376 | nd_mode(i, j, it) = tr_seri(i, j, it + nbtr_sulgas) & ! [DRY part/kgA in bin it] |
|---|
| 377 | * pplay(i, j) / t_seri(i, j) / RD ! [air mass concentration in kg air /m3A] |
|---|
| 378 | |
|---|
| 379 | IF(flag_new_strat_compo) THEN |
|---|
| 380 | ! SAD_sulfate: SAD WET sulfate aerosols (cm2/cm3) |
|---|
| 381 | SAD_sulfate(i, j) = SAD_sulfate(i, j) + nd_mode(i, j, it) & ! [DRY part/m3A (in bin it)] |
|---|
| 382 | * 4. * RPI * (mdw(it) * f_r_wetB(i, j, it) / 2.)**2. & ! [WET SA of part it in m2] |
|---|
| 383 | * 1.e-2 ! conversion from m2/m3 to cm2/cm3A |
|---|
| 384 | ELSE |
|---|
| 385 | ! SAD_sulfate: SAD WET sulfate aerosols (cm2/cm3) |
|---|
| 386 | SAD_sulfate(i, j) = SAD_sulfate(i, j) + nd_mode(i, j, it) & ! [DRY part/m3A (in bin it)] |
|---|
| 387 | * 4. * RPI * (mdw(it) * f_r_wet(i, j) / 2.)**2. & ! [WET SA of part it in m2] |
|---|
| 388 | * 1.e-2 ! conversion from m2/m3 to cm2/cm3A |
|---|
| 389 | ENDIF |
|---|
| 390 | ENDDO |
|---|
| 391 | ENDDO |
|---|
| 392 | ENDDO |
|---|
| 393 | |
|---|
| 394 | END SUBROUTINE traccoag |
|---|
| 395 | |
|---|
| 396 | END MODULE traccoag_mod |
|---|
