| 1 | ! |
|---|
| 2 | ! $Id: traccoag_mod.F90 4293 2022-10-04 19:42:41Z fhourdin $ |
|---|
| 3 | ! |
|---|
| 4 | MODULE traccoag_mod |
|---|
| 5 | ! |
|---|
| 6 | ! This module calculates the concentration of aerosol particles in certain size bins |
|---|
| 7 | ! considering coagulation and sedimentation. |
|---|
| 8 | ! |
|---|
| 9 | CONTAINS |
|---|
| 10 | |
|---|
| 11 | SUBROUTINE traccoag(pdtphys, gmtime, debutphy, julien, & |
|---|
| 12 | presnivs, xlat, xlon, pphis, pphi, & |
|---|
| 13 | t_seri, pplay, paprs, sh, rh, tr_seri) |
|---|
| 14 | |
|---|
| 15 | USE phys_local_var_mod, ONLY: mdw, R2SO4, DENSO4, f_r_wet, surf_PM25_sulf, & |
|---|
| 16 | & budg_emi_ocs, budg_emi_so2, budg_emi_h2so4, budg_emi_part |
|---|
| 17 | |
|---|
| 18 | USE dimphy |
|---|
| 19 | USE infotrac_phy, ONLY : nbtr_bin, nbtr_sulgas, nbtr, id_SO2_strat |
|---|
| 20 | USE aerophys |
|---|
| 21 | USE geometry_mod, ONLY : cell_area, boundslat |
|---|
| 22 | USE mod_grid_phy_lmdz |
|---|
| 23 | USE mod_phys_lmdz_mpi_data, ONLY : is_mpi_root |
|---|
| 24 | USE mod_phys_lmdz_para, only: gather, scatter |
|---|
| 25 | USE phys_cal_mod, ONLY : year_len, year_cur,mth_cur, day_cur, hour |
|---|
| 26 | USE sulfate_aer_mod |
|---|
| 27 | USE phys_local_var_mod, ONLY: stratomask |
|---|
| 28 | USE YOMCST |
|---|
| 29 | USE print_control_mod, ONLY: lunout |
|---|
| 30 | USE strataer_mod |
|---|
| 31 | |
|---|
| 32 | IMPLICIT NONE |
|---|
| 33 | |
|---|
| 34 | ! Input argument |
|---|
| 35 | !--------------- |
|---|
| 36 | REAL,INTENT(IN) :: pdtphys ! Pas d'integration pour la physique (seconde) |
|---|
| 37 | REAL,INTENT(IN) :: gmtime ! Heure courante |
|---|
| 38 | LOGICAL,INTENT(IN) :: debutphy ! le flag de l'initialisation de la physique |
|---|
| 39 | INTEGER,INTENT(IN) :: julien ! Jour julien |
|---|
| 40 | |
|---|
| 41 | REAL,DIMENSION(klev),INTENT(IN) :: presnivs! pressions approximat. des milieux couches (en PA) |
|---|
| 42 | REAL,DIMENSION(klon),INTENT(IN) :: xlat ! latitudes pour chaque point |
|---|
| 43 | REAL,DIMENSION(klon),INTENT(IN) :: xlon ! longitudes pour chaque point |
|---|
| 44 | REAL,DIMENSION(klon),INTENT(IN) :: pphis ! geopotentiel du sol |
|---|
| 45 | REAL,DIMENSION(klon,klev),INTENT(IN) :: pphi ! geopotentiel de chaque couche |
|---|
| 46 | |
|---|
| 47 | REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature |
|---|
| 48 | REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa) |
|---|
| 49 | REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression pour chaque inter-couche (en Pa) |
|---|
| 50 | REAL,DIMENSION(klon,klev),INTENT(IN) :: sh ! humidite specifique |
|---|
| 51 | REAL,DIMENSION(klon,klev),INTENT(IN) :: rh ! humidite relative |
|---|
| 52 | |
|---|
| 53 | ! Output argument |
|---|
| 54 | !---------------- |
|---|
| 55 | REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA] |
|---|
| 56 | |
|---|
| 57 | ! Local variables |
|---|
| 58 | !---------------- |
|---|
| 59 | REAL :: m_aer_emiss_vol_daily ! daily injection mass emission |
|---|
| 60 | INTEGER :: it, k, i, ilon, ilev, itime, i_int, ieru |
|---|
| 61 | LOGICAL,DIMENSION(klon,klev) :: is_strato ! true = above tropopause, false = below |
|---|
| 62 | REAL,DIMENSION(klon,klev) :: m_air_gridbox ! mass of air in every grid box [kg] |
|---|
| 63 | REAL,DIMENSION(klon_glo,klev,nbtr) :: tr_seri_glo ! Concentration Traceur [U/KgA] |
|---|
| 64 | REAL,DIMENSION(klev+1) :: altLMDz ! altitude of layer interfaces in m |
|---|
| 65 | REAL,DIMENSION(klev) :: f_lay_emiss ! fraction of emission for every vertical layer |
|---|
| 66 | REAL :: f_lay_sum ! sum of layer emission fractions |
|---|
| 67 | REAL :: alt ! altitude for integral calculation |
|---|
| 68 | INTEGER,PARAMETER :: n_int_alt=10 ! number of subintervals for integration over Gaussian emission profile |
|---|
| 69 | REAL,DIMENSION(nbtr_bin) :: r_bin ! particle radius in size bin [m] |
|---|
| 70 | REAL,DIMENSION(nbtr_bin) :: r_lower ! particle radius at lower bin boundary [m] |
|---|
| 71 | REAL,DIMENSION(nbtr_bin) :: r_upper ! particle radius at upper bin boundary [m] |
|---|
| 72 | REAL,DIMENSION(nbtr_bin) :: m_part_dry ! mass of one dry particle in size bin [kg] |
|---|
| 73 | REAL :: zrho ! Density of air [kg/m3] |
|---|
| 74 | REAL :: zdz ! thickness of atm. model layer in m |
|---|
| 75 | REAL,DIMENSION(klev) :: zdm ! mass of atm. model layer in kg |
|---|
| 76 | REAL,DIMENSION(klon,klev) :: dens_aer ! density of aerosol particles [kg/m3 aerosol] with default H2SO4 mass fraction |
|---|
| 77 | REAL :: emission ! emission |
|---|
| 78 | REAL :: theta_min, theta_max ! for SAI computation between two latitudes |
|---|
| 79 | REAL :: dlat_loc |
|---|
| 80 | |
|---|
| 81 | IF (is_mpi_root) THEN |
|---|
| 82 | WRITE(lunout,*) 'in traccoag: date from phys_cal_mod =',year_cur,'-',mth_cur,'-',day_cur,'-',hour |
|---|
| 83 | WRITE(lunout,*) 'IN traccoag flag_sulf_emit: ',flag_sulf_emit |
|---|
| 84 | ENDIF |
|---|
| 85 | |
|---|
| 86 | DO it=1, nbtr_bin |
|---|
| 87 | r_bin(it)=mdw(it)/2. |
|---|
| 88 | ENDDO |
|---|
| 89 | |
|---|
| 90 | !--set boundaries of size bins |
|---|
| 91 | DO it=1, nbtr_bin |
|---|
| 92 | IF (it.EQ.1) THEN |
|---|
| 93 | r_upper(it)=sqrt(r_bin(it+1)*r_bin(it)) |
|---|
| 94 | r_lower(it)=r_bin(it)**2./r_upper(it) |
|---|
| 95 | ELSEIF (it.EQ.nbtr_bin) THEN |
|---|
| 96 | r_lower(it)=sqrt(r_bin(it)*r_bin(it-1)) |
|---|
| 97 | r_upper(it)=r_bin(it)**2./r_lower(it) |
|---|
| 98 | ELSE |
|---|
| 99 | r_lower(it)=sqrt(r_bin(it)*r_bin(it-1)) |
|---|
| 100 | r_upper(it)=sqrt(r_bin(it+1)*r_bin(it)) |
|---|
| 101 | ENDIF |
|---|
| 102 | ENDDO |
|---|
| 103 | |
|---|
| 104 | IF (debutphy .and. is_mpi_root) THEN |
|---|
| 105 | DO it=1, nbtr_bin |
|---|
| 106 | WRITE(lunout,*) 'radius bin', it, ':', r_bin(it), '(from', r_lower(it), 'to', r_upper(it), ')' |
|---|
| 107 | ENDDO |
|---|
| 108 | ENDIF |
|---|
| 109 | |
|---|
| 110 | !--initialising logical is_strato from stratomask |
|---|
| 111 | is_strato(:,:)=.FALSE. |
|---|
| 112 | WHERE (stratomask.GT.0.5) is_strato=.TRUE. |
|---|
| 113 | |
|---|
| 114 | ! STRACOMP (H2O, P, t_seri -> aerosol composition (R2SO4)) |
|---|
| 115 | ! H2SO4 mass fraction in aerosol (%) |
|---|
| 116 | CALL stracomp(sh,t_seri,pplay) |
|---|
| 117 | |
|---|
| 118 | ! aerosol density (gr/cm3) |
|---|
| 119 | CALL denh2sa(t_seri) |
|---|
| 120 | |
|---|
| 121 | ! compute factor for converting dry to wet radius (for every grid box) |
|---|
| 122 | f_r_wet(:,:) = (dens_aer_dry/(DENSO4(:,:)*1000.)/(R2SO4(:,:)/100.))**(1./3.) |
|---|
| 123 | |
|---|
| 124 | !--calculate mass of air in every grid box |
|---|
| 125 | DO ilon=1, klon |
|---|
| 126 | DO ilev=1, klev |
|---|
| 127 | m_air_gridbox(ilon,ilev)=(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG*cell_area(ilon) |
|---|
| 128 | ENDDO |
|---|
| 129 | ENDDO |
|---|
| 130 | |
|---|
| 131 | ! IF (debutphy) THEN |
|---|
| 132 | ! CALL gather(tr_seri, tr_seri_glo) |
|---|
| 133 | ! IF (MAXVAL(tr_seri_glo).LT.1.e-30) THEN |
|---|
| 134 | !--initialising tracer concentrations to zero |
|---|
| 135 | ! DO it=1, nbtr |
|---|
| 136 | ! tr_seri(:,:,it)=0.0 |
|---|
| 137 | ! ENDDO |
|---|
| 138 | ! ENDIF |
|---|
| 139 | ! ENDIF |
|---|
| 140 | |
|---|
| 141 | !--initialise emission diagnostics |
|---|
| 142 | budg_emi_ocs(:)=0.0 |
|---|
| 143 | budg_emi_so2(:)=0.0 |
|---|
| 144 | budg_emi_h2so4(:)=0.0 |
|---|
| 145 | budg_emi_part(:)=0.0 |
|---|
| 146 | |
|---|
| 147 | !--sulfur emission, depending on chosen scenario (flag_sulf_emit) |
|---|
| 148 | SELECT CASE(flag_sulf_emit) |
|---|
| 149 | |
|---|
| 150 | CASE(0) ! background aerosol |
|---|
| 151 | ! do nothing (no emission) |
|---|
| 152 | |
|---|
| 153 | CASE(1) ! volcanic eruption |
|---|
| 154 | !--only emit on day of eruption |
|---|
| 155 | ! stretch emission over one day of Pinatubo eruption |
|---|
| 156 | DO ieru=1, nErupt |
|---|
| 157 | IF (year_cur==year_emit_vol(ieru).AND.mth_cur==mth_emit_vol(ieru).AND.& |
|---|
| 158 | day_cur>=day_emit_vol(ieru).AND.day_cur<(day_emit_vol(ieru)+injdur)) THEN |
|---|
| 159 | ! |
|---|
| 160 | ! daily injection mass emission - NL |
|---|
| 161 | m_aer_emiss_vol_daily = m_aer_emiss_vol(ieru)/(REAL(injdur)*REAL(ponde_lonlat_vol(ieru))) |
|---|
| 162 | WRITE(lunout,*) 'IN traccoag DD m_aer_emiss_vol(ieru)=',m_aer_emiss_vol(ieru), & |
|---|
| 163 | 'ponde_lonlat_vol(ieru)=',ponde_lonlat_vol(ieru),'(injdur*ponde_lonlat_vol(ieru))', & |
|---|
| 164 | (injdur*ponde_lonlat_vol(ieru)),'m_aer_emiss_vol_daily=',m_aer_emiss_vol_daily,'ieru=',ieru |
|---|
| 165 | WRITE(lunout,*) 'IN traccoag, dlon=',dlon |
|---|
| 166 | DO i=1,klon |
|---|
| 167 | !Pinatubo eruption at 15.14N, 120.35E |
|---|
| 168 | dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes |
|---|
| 169 | WRITE(lunout,*) 'IN traccoag, dlat=',dlat_loc |
|---|
| 170 | IF ( xlat(i).GE.xlat_min_vol(ieru)-dlat_loc .AND. xlat(i).LT.xlat_max_vol(ieru)+dlat_loc .AND. & |
|---|
| 171 | xlon(i).GE.xlon_min_vol(ieru)-dlon .AND. xlon(i).LT.xlon_max_vol(ieru)+dlon ) THEN |
|---|
| 172 | ! |
|---|
| 173 | WRITE(lunout,*) 'coordinates of volcanic injection point=',xlat(i),xlon(i),day_cur,mth_cur,year_cur |
|---|
| 174 | WRITE(lunout,*) 'DD m_aer_emiss_vol_daily=',m_aer_emiss_vol_daily |
|---|
| 175 | ! compute altLMDz |
|---|
| 176 | altLMDz(:)=0.0 |
|---|
| 177 | DO k=1, klev |
|---|
| 178 | zrho=pplay(i,k)/t_seri(i,k)/RD !air density in kg/m3 |
|---|
| 179 | zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG !mass of layer in kg |
|---|
| 180 | zdz=zdm(k)/zrho !thickness of layer in m |
|---|
| 181 | altLMDz(k+1)=altLMDz(k)+zdz !altitude of interface |
|---|
| 182 | ENDDO |
|---|
| 183 | |
|---|
| 184 | SELECT CASE(flag_sulf_emit_distrib) |
|---|
| 185 | |
|---|
| 186 | CASE(0) ! Gaussian distribution |
|---|
| 187 | !compute distribution of emission to vertical model layers (based on Gaussian peak in altitude) |
|---|
| 188 | f_lay_sum=0.0 |
|---|
| 189 | DO k=1, klev |
|---|
| 190 | f_lay_emiss(k)=0.0 |
|---|
| 191 | DO i_int=1, n_int_alt |
|---|
| 192 | alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt) |
|---|
| 193 | f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_vol(ieru))* & |
|---|
| 194 | & exp(-0.5*((alt-altemiss_vol(ieru))/sigma_alt_vol(ieru))**2.)* & |
|---|
| 195 | & (altLMDz(k+1)-altLMDz(k))/float(n_int_alt) |
|---|
| 196 | ENDDO |
|---|
| 197 | f_lay_sum=f_lay_sum+f_lay_emiss(k) |
|---|
| 198 | ENDDO |
|---|
| 199 | |
|---|
| 200 | CASE(1) ! Uniform distribution |
|---|
| 201 | ! In this case, parameter sigma_alt_vol(ieru) is considered to be half the |
|---|
| 202 | ! height of the injection, centered around altemiss_vol(ieru) |
|---|
| 203 | DO k=1, klev |
|---|
| 204 | f_lay_emiss(k)=max(min(altemiss_vol(ieru)+sigma_alt_vol(ieru),altLMDz(k+1))- & |
|---|
| 205 | & max(altemiss_vol(ieru)-sigma_alt_vol(ieru),altLMDz(k)),0.)/(2.*sigma_alt_vol(ieru)) |
|---|
| 206 | f_lay_sum=f_lay_sum+f_lay_emiss(k) |
|---|
| 207 | ENDDO |
|---|
| 208 | |
|---|
| 209 | END SELECT ! End CASE over flag_sulf_emit_distrib) |
|---|
| 210 | |
|---|
| 211 | WRITE(lunout,*) "IN traccoag m_aer_emiss_vol=",m_aer_emiss_vol(ieru) |
|---|
| 212 | WRITE(lunout,*) "IN traccoag f_lay_emiss=",f_lay_emiss |
|---|
| 213 | !correct for step integration error |
|---|
| 214 | f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum |
|---|
| 215 | !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss) |
|---|
| 216 | !vertically distributed emission |
|---|
| 217 | DO k=1, klev |
|---|
| 218 | ! stretch emission over one day of Pinatubo eruption |
|---|
| 219 | emission=m_aer_emiss_vol_daily*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/1./(86400.-pdtphys) |
|---|
| 220 | tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys |
|---|
| 221 | budg_emi_so2(i)=budg_emi_so2(i)+emission*zdm(k)*mSatom/mSO2mol |
|---|
| 222 | ENDDO |
|---|
| 223 | ENDIF ! emission grid cell |
|---|
| 224 | ENDDO ! klon loop |
|---|
| 225 | WRITE(lunout,*) "IN traccoag (ieru=",ieru,") m_aer_emiss_vol_daily=",m_aer_emiss_vol_daily |
|---|
| 226 | ENDIF ! emission period |
|---|
| 227 | ENDDO ! eruption number |
|---|
| 228 | |
|---|
| 229 | CASE(2) ! stratospheric aerosol injections (SAI) |
|---|
| 230 | ! |
|---|
| 231 | DO i=1,klon |
|---|
| 232 | ! SAI standard scenario with continuous emission from 1 grid point at the equator |
|---|
| 233 | ! SAI emission on single month |
|---|
| 234 | ! SAI continuous emission o |
|---|
| 235 | dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes |
|---|
| 236 | IF ( xlat(i).GE.xlat_sai-dlat_loc .AND. xlat(i).LT.xlat_sai+dlat_loc .AND. & |
|---|
| 237 | & xlon(i).GE.xlon_sai-dlon .AND. xlon(i).LT.xlon_sai+dlon ) THEN |
|---|
| 238 | ! |
|---|
| 239 | ! compute altLMDz |
|---|
| 240 | altLMDz(:)=0.0 |
|---|
| 241 | DO k=1, klev |
|---|
| 242 | zrho=pplay(i,k)/t_seri(i,k)/RD !air density in kg/m3 |
|---|
| 243 | zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG !mass of layer in kg |
|---|
| 244 | zdz=zdm(k)/zrho !thickness of layer in m |
|---|
| 245 | altLMDz(k+1)=altLMDz(k)+zdz !altitude of interface |
|---|
| 246 | ENDDO |
|---|
| 247 | |
|---|
| 248 | SELECT CASE(flag_sulf_emit_distrib) |
|---|
| 249 | |
|---|
| 250 | CASE(0) ! Gaussian distribution |
|---|
| 251 | !compute distribution of emission to vertical model layers (based on Gaussian peak in altitude) |
|---|
| 252 | f_lay_sum=0.0 |
|---|
| 253 | DO k=1, klev |
|---|
| 254 | f_lay_emiss(k)=0.0 |
|---|
| 255 | DO i_int=1, n_int_alt |
|---|
| 256 | alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt) |
|---|
| 257 | f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_sai)* & |
|---|
| 258 | & exp(-0.5*((alt-altemiss_sai)/sigma_alt_sai)**2.)* & |
|---|
| 259 | & (altLMDz(k+1)-altLMDz(k))/float(n_int_alt) |
|---|
| 260 | ENDDO |
|---|
| 261 | f_lay_sum=f_lay_sum+f_lay_emiss(k) |
|---|
| 262 | ENDDO |
|---|
| 263 | |
|---|
| 264 | CASE(1) ! Uniform distribution |
|---|
| 265 | f_lay_sum=0.0 |
|---|
| 266 | ! In this case, parameter sigma_alt_vol(ieru) is considered to be half |
|---|
| 267 | ! the height of the injection, centered around altemiss_sai |
|---|
| 268 | DO k=1, klev |
|---|
| 269 | f_lay_emiss(k)=max(min(altemiss_sai+sigma_alt_sai,altLMDz(k+1))- & |
|---|
| 270 | & max(altemiss_sai-sigma_alt_sai,altLMDz(k)),0.)/(2.*sigma_alt_sai) |
|---|
| 271 | f_lay_sum=f_lay_sum+f_lay_emiss(k) |
|---|
| 272 | ENDDO |
|---|
| 273 | |
|---|
| 274 | END SELECT ! Gaussian or uniform distribution |
|---|
| 275 | |
|---|
| 276 | !correct for step integration error |
|---|
| 277 | f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum |
|---|
| 278 | !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss) |
|---|
| 279 | !vertically distributed emission |
|---|
| 280 | DO k=1, klev |
|---|
| 281 | ! stretch emission over whole year (360d) |
|---|
| 282 | emission=m_aer_emiss_sai*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/FLOAT(year_len)/86400. |
|---|
| 283 | tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys |
|---|
| 284 | budg_emi_so2(i)=budg_emi_so2(i)+emission*zdm(k)*mSatom/mSO2mol |
|---|
| 285 | ENDDO |
|---|
| 286 | |
|---|
| 287 | ! !emission as monodisperse particles with 0.1um dry radius (BIN21) |
|---|
| 288 | ! !vertically distributed emission |
|---|
| 289 | ! DO k=1, klev |
|---|
| 290 | ! ! stretch emission over whole year (360d) |
|---|
| 291 | ! emission=m_aer_emiss*(mH2SO4mol/mSatom)/m_part_dry(21)/m_air_gridbox(i,k)*f_lay_emiss(k)/FLOAT(year_len)/86400. |
|---|
| 292 | ! tr_seri(i,k,id_BIN01_strat+20)=tr_seri(i,k,id_BIN01_strat+20)+emission*pdtphys |
|---|
| 293 | ! budg_emi_part(i)=budg_emi_part(i)+emission*zdm(k)*mSatom/mH2SO4mol |
|---|
| 294 | ! ENDDO |
|---|
| 295 | ENDIF ! emission grid cell |
|---|
| 296 | ENDDO ! klon loop |
|---|
| 297 | |
|---|
| 298 | CASE(3) ! --- SAI injection over a single band of longitude and between |
|---|
| 299 | ! lat_min and lat_max |
|---|
| 300 | |
|---|
| 301 | DO i=1,klon |
|---|
| 302 | ! SAI scenario with continuous emission |
|---|
| 303 | dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes |
|---|
| 304 | theta_min = max(xlat(i)-dlat_loc,xlat_min_sai) |
|---|
| 305 | theta_max = min(xlat(i)+dlat_loc,xlat_max_sai) |
|---|
| 306 | IF ( xlat(i).GE.xlat_min_sai-dlat_loc .AND. xlat(i).LT.xlat_max_sai+dlat_loc .AND. & |
|---|
| 307 | & xlon(i).GE.xlon_sai-dlon .AND. xlon(i).LT.xlon_sai+dlon ) THEN |
|---|
| 308 | ! |
|---|
| 309 | ! compute altLMDz |
|---|
| 310 | altLMDz(:)=0.0 |
|---|
| 311 | DO k=1, klev |
|---|
| 312 | zrho=pplay(i,k)/t_seri(i,k)/RD !air density in kg/m3 |
|---|
| 313 | zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG !mass of layer in kg |
|---|
| 314 | zdz=zdm(k)/zrho !thickness of layer in m |
|---|
| 315 | altLMDz(k+1)=altLMDz(k)+zdz !altitude of interface |
|---|
| 316 | ENDDO |
|---|
| 317 | |
|---|
| 318 | SELECT CASE(flag_sulf_emit_distrib) |
|---|
| 319 | |
|---|
| 320 | CASE(0) ! Gaussian distribution |
|---|
| 321 | !compute distribution of emission to vertical model layers (based on |
|---|
| 322 | !Gaussian peak in altitude) |
|---|
| 323 | f_lay_sum=0.0 |
|---|
| 324 | DO k=1, klev |
|---|
| 325 | f_lay_emiss(k)=0.0 |
|---|
| 326 | DO i_int=1, n_int_alt |
|---|
| 327 | alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt) |
|---|
| 328 | f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_sai)* & |
|---|
| 329 | & exp(-0.5*((alt-altemiss_sai)/sigma_alt_sai)**2.)* & |
|---|
| 330 | & (altLMDz(k+1)-altLMDz(k))/float(n_int_alt) |
|---|
| 331 | ENDDO |
|---|
| 332 | f_lay_sum=f_lay_sum+f_lay_emiss(k) |
|---|
| 333 | ENDDO |
|---|
| 334 | |
|---|
| 335 | CASE(1) ! Uniform distribution |
|---|
| 336 | f_lay_sum=0.0 |
|---|
| 337 | ! In this case, parameter sigma_alt_vol(ieru) is considered to be half |
|---|
| 338 | ! the height of the injection, centered around altemiss_sai |
|---|
| 339 | DO k=1, klev |
|---|
| 340 | f_lay_emiss(k)=max(min(altemiss_sai+sigma_alt_sai,altLMDz(k+1))- & |
|---|
| 341 | & max(altemiss_sai-sigma_alt_sai,altLMDz(k)),0.)/(2.*sigma_alt_sai) |
|---|
| 342 | f_lay_sum=f_lay_sum+f_lay_emiss(k) |
|---|
| 343 | ENDDO |
|---|
| 344 | |
|---|
| 345 | END SELECT ! Gaussian or uniform distribution |
|---|
| 346 | |
|---|
| 347 | !correct for step integration error |
|---|
| 348 | f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum |
|---|
| 349 | !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss) |
|---|
| 350 | !vertically distributed emission |
|---|
| 351 | DO k=1, klev |
|---|
| 352 | ! stretch emission over whole year (360d) |
|---|
| 353 | emission=m_aer_emiss_sai*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/ & |
|---|
| 354 | & FLOAT(year_len)/86400.*(sin(theta_max/180.*RPI)-sin(theta_min/180.*RPI))/ & |
|---|
| 355 | & (sin(xlat_max_sai/180.*RPI)-sin(xlat_min_sai/180.*RPI)) |
|---|
| 356 | tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys |
|---|
| 357 | budg_emi_so2(i)=budg_emi_so2(i)+emission*zdm(k)*mSatom/mSO2mol |
|---|
| 358 | ENDDO |
|---|
| 359 | |
|---|
| 360 | ! !emission as monodisperse particles with 0.1um dry radius (BIN21) |
|---|
| 361 | ! !vertically distributed emission |
|---|
| 362 | ! DO k=1, klev |
|---|
| 363 | ! ! stretch emission over whole year (360d) |
|---|
| 364 | ! emission=m_aer_emiss*(mH2SO4mol/mSatom)/m_part_dry(21)/m_air_gridbox(i,k)*f_lay_emiss(k)/year_len/86400 |
|---|
| 365 | ! tr_seri(i,k,id_BIN01_strat+20)=tr_seri(i,k,id_BIN01_strat+20)+emission*pdtphys |
|---|
| 366 | ! budg_emi_part(i)=budg_emi_part(i)+emission*zdm(k)*mSatom/mH2SO4mol |
|---|
| 367 | ! ENDDO |
|---|
| 368 | ENDIF ! emission grid cell |
|---|
| 369 | ENDDO ! klon loop |
|---|
| 370 | |
|---|
| 371 | END SELECT ! emission scenario (flag_sulf_emit) |
|---|
| 372 | |
|---|
| 373 | !--read background concentrations of OCS and SO2 and lifetimes from input file |
|---|
| 374 | !--update the variables defined in phys_local_var_mod |
|---|
| 375 | CALL interp_sulf_input(debutphy,pdtphys,paprs,tr_seri) |
|---|
| 376 | |
|---|
| 377 | !--convert OCS to SO2 in the stratosphere |
|---|
| 378 | CALL ocs_to_so2(pdtphys,tr_seri,t_seri,pplay,paprs,is_strato) |
|---|
| 379 | |
|---|
| 380 | !--convert SO2 to H2SO4 |
|---|
| 381 | CALL so2_to_h2so4(pdtphys,tr_seri,t_seri,pplay,paprs,is_strato) |
|---|
| 382 | |
|---|
| 383 | !--common routine for nucleation and condensation/evaporation with adaptive timestep |
|---|
| 384 | CALL micphy_tstep(pdtphys,tr_seri,t_seri,pplay,paprs,rh,is_strato) |
|---|
| 385 | |
|---|
| 386 | !--call coagulation routine |
|---|
| 387 | CALL coagulate(pdtphys,mdw,tr_seri,t_seri,pplay,dens_aer,is_strato) |
|---|
| 388 | |
|---|
| 389 | !--call sedimentation routine |
|---|
| 390 | CALL aer_sedimnt(pdtphys, t_seri, pplay, paprs, tr_seri, dens_aer) |
|---|
| 391 | |
|---|
| 392 | !--compute mass concentration of PM2.5 sulfate particles (wet diameter and mass) at the surface for health studies |
|---|
| 393 | surf_PM25_sulf(:)=0.0 |
|---|
| 394 | DO i=1,klon |
|---|
| 395 | DO it=1, nbtr_bin |
|---|
| 396 | IF (mdw(it) .LT. 2.5e-6) THEN |
|---|
| 397 | !surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas)*m_part(i,1,it) & |
|---|
| 398 | !assume that particles consist of ammonium sulfate at the surface (132g/mol) |
|---|
| 399 | !and are dry at T = 20 deg. C and 50 perc. humidity |
|---|
| 400 | surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas) & |
|---|
| 401 | & *132./98.*dens_aer_dry*4./3.*RPI*(mdw(it)/2.)**3 & |
|---|
| 402 | & *pplay(i,1)/t_seri(i,1)/RD*1.e9 |
|---|
| 403 | ENDIF |
|---|
| 404 | ENDDO |
|---|
| 405 | ENDDO |
|---|
| 406 | |
|---|
| 407 | ! CALL minmaxsimple(tr_seri(:,:,id_SO2_strat),0.0,0.0,'fin SO2') |
|---|
| 408 | ! DO it=1, nbtr_bin |
|---|
| 409 | ! CALL minmaxsimple(tr_seri(:,:,nbtr_sulgas+it),0.0,0.0,'fin bin ') |
|---|
| 410 | ! ENDDO |
|---|
| 411 | |
|---|
| 412 | END SUBROUTINE traccoag |
|---|
| 413 | |
|---|
| 414 | END MODULE traccoag_mod |
|---|
