[2690] | 1 | MODULE traccoag_mod |
---|
| 2 | ! |
---|
| 3 | ! This module calculates the concentration of aerosol particles in certain size bins |
---|
| 4 | ! considering coagulation and sedimentation. |
---|
| 5 | ! |
---|
| 6 | CONTAINS |
---|
| 7 | |
---|
| 8 | SUBROUTINE traccoag(pdtphys, gmtime, debutphy, julien, & |
---|
| 9 | presnivs, xlat, xlon, pphis, pphi, & |
---|
[2752] | 10 | t_seri, pplay, paprs, sh, rh, tr_seri) |
---|
[2690] | 11 | |
---|
[2752] | 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 |
---|
[2690] | 14 | |
---|
| 15 | USE dimphy |
---|
| 16 | USE infotrac |
---|
| 17 | USE aerophys |
---|
| 18 | USE geometry_mod, ONLY : cell_area |
---|
| 19 | USE mod_grid_phy_lmdz |
---|
| 20 | USE mod_phys_lmdz_mpi_data, ONLY : is_mpi_root |
---|
| 21 | USE mod_phys_lmdz_para, only: gather, scatter |
---|
| 22 | USE phys_cal_mod |
---|
| 23 | USE sulfate_aer_mod |
---|
| 24 | USE phys_local_var_mod, ONLY: stratomask |
---|
| 25 | USE YOMCST |
---|
| 26 | |
---|
| 27 | IMPLICIT NONE |
---|
| 28 | |
---|
| 29 | ! Input argument |
---|
| 30 | !--------------- |
---|
| 31 | REAL,INTENT(IN) :: pdtphys ! Pas d'integration pour la physique (seconde) |
---|
| 32 | REAL,INTENT(IN) :: gmtime ! Heure courante |
---|
| 33 | LOGICAL,INTENT(IN) :: debutphy ! le flag de l'initialisation de la physique |
---|
| 34 | INTEGER,INTENT(IN) :: julien ! Jour julien |
---|
| 35 | |
---|
| 36 | REAL,DIMENSION(klev),INTENT(IN) :: presnivs! pressions approximat. des milieux couches (en PA) |
---|
| 37 | REAL,DIMENSION(klon),INTENT(IN) :: xlat ! latitudes pour chaque point |
---|
| 38 | REAL,DIMENSION(klon),INTENT(IN) :: xlon ! longitudes pour chaque point |
---|
| 39 | REAL,DIMENSION(klon),INTENT(IN) :: pphis ! geopotentiel du sol |
---|
| 40 | REAL,DIMENSION(klon,klev),INTENT(IN) :: pphi ! geopotentiel de chaque couche |
---|
| 41 | |
---|
| 42 | REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature |
---|
| 43 | REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa) |
---|
| 44 | REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression pour chaque inter-couche (en Pa) |
---|
| 45 | REAL,DIMENSION(klon,klev),INTENT(IN) :: sh ! humidite specifique |
---|
| 46 | REAL,DIMENSION(klon,klev),INTENT(IN) :: rh ! humidite relative |
---|
| 47 | |
---|
| 48 | ! Output argument |
---|
| 49 | !---------------- |
---|
| 50 | REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA] |
---|
| 51 | |
---|
| 52 | ! Local variables |
---|
| 53 | !---------------- |
---|
| 54 | ! flag for sulfur emission scenario: (0) background aerosol ; (1) volcanic eruption ; (2) stratospheric aerosol injections (SAI) |
---|
| 55 | INTEGER,PARAMETER :: flag_sulf_emit=2 |
---|
| 56 | ! |
---|
| 57 | !--flag_sulf_emit=1 --example Pinatubo |
---|
| 58 | INTEGER,PARAMETER :: year_emit_vol=1991 ! year of emission date |
---|
| 59 | INTEGER,PARAMETER :: mth_emit_vol=6 ! month of emission date |
---|
| 60 | INTEGER,PARAMETER :: day_emit_vol=15 ! day of emission date |
---|
| 61 | REAL,PARAMETER :: m_aer_emiss_vol=7.e9 ! emitted sulfur mass in kgS, e.g. 7Tg(S)=14Tg(SO2) |
---|
| 62 | REAL,PARAMETER :: altemiss_vol=17.e3 ! emission altitude in m |
---|
| 63 | REAL,PARAMETER :: sigma_alt_vol=1.e3 ! standard deviation of emission altitude in m |
---|
| 64 | REAL,PARAMETER :: xlat_vol=15.14 ! latitude of volcano in degree |
---|
| 65 | REAL,PARAMETER :: xlon_vol=120.35 ! longitude of volcano in degree |
---|
| 66 | |
---|
| 67 | !--flag_sulf_emit=2 --SAI |
---|
| 68 | REAL,PARAMETER :: m_aer_emiss_sai=1.e10 ! emitted sulfur mass in kgS, eg 1e9=1TgS, 1e10=10TgS |
---|
| 69 | REAL,PARAMETER :: altemiss_sai=17.e3 ! emission altitude in m |
---|
| 70 | REAL,PARAMETER :: sigma_alt_sai=1.e3 ! standard deviation of emission altitude in m |
---|
[2752] | 71 | REAL,PARAMETER :: xlat_sai=0.01 ! latitude of SAI in degree |
---|
[2690] | 72 | REAL,PARAMETER :: xlon_sai=120.35 ! longitude of SAI in degree |
---|
| 73 | |
---|
| 74 | !--other local variables |
---|
| 75 | INTEGER :: it, k, i, ilon, ilev, itime, i_int |
---|
| 76 | LOGICAL,DIMENSION(klon,klev) :: is_strato ! true = above tropopause, false = below |
---|
| 77 | REAL,DIMENSION(klon,klev) :: m_air_gridbox ! mass of air in every grid box [kg] |
---|
| 78 | REAL,DIMENSION(klon_glo,klev,nbtr) :: tr_seri_glo ! Concentration Traceur [U/KgA] |
---|
| 79 | REAL,DIMENSION(klev+1) :: altLMDz ! altitude of layer interfaces in m |
---|
| 80 | REAL,DIMENSION(klev) :: f_lay_emiss ! fraction of emission for every vertical layer |
---|
| 81 | REAL :: f_lay_sum ! sum of layer emission fractions |
---|
[2699] | 82 | REAL :: alt ! altitude for integral calculation |
---|
[2690] | 83 | INTEGER,PARAMETER :: n_int_alt=10 ! number of subintervals for integration over Gaussian emission profile |
---|
| 84 | REAL,DIMENSION(nbtr_bin) :: r_bin ! particle radius in size bin [m] |
---|
| 85 | REAL,DIMENSION(nbtr_bin) :: r_lower ! particle radius at lower bin boundary [m] |
---|
| 86 | REAL,DIMENSION(nbtr_bin) :: r_upper ! particle radius at upper bin boundary [m] |
---|
| 87 | REAL,DIMENSION(nbtr_bin) :: m_part_dry ! mass of one dry particle in size bin [kg] |
---|
| 88 | REAL :: zrho ! Density of air [kg/m3] |
---|
| 89 | REAL :: zdz ! thickness of atm. model layer in m |
---|
[2752] | 90 | REAL,DIMENSION(klev) :: zdm ! mass of atm. model layer in kg |
---|
[2690] | 91 | REAL,DIMENSION(klon,klev) :: dens_aer ! density of aerosol particles [kg/m3 aerosol] with default H2SO4 mass fraction |
---|
[2704] | 92 | REAL :: dlat, dlon ! d latitude and d longitude of grid in degree |
---|
[2752] | 93 | REAL :: emission ! emission |
---|
[2690] | 94 | |
---|
| 95 | IF (is_mpi_root) THEN |
---|
| 96 | PRINT *,'in traccoag: date from phys_cal_mod =',year_cur,'-',mth_cur,'-',day_cur,'-',hour |
---|
| 97 | ENDIF |
---|
| 98 | |
---|
[2704] | 99 | dlat=180./2./FLOAT(nbp_lat) ! d latitude in degree |
---|
| 100 | dlon=360./2./FLOAT(nbp_lon) ! d longitude in degree |
---|
| 101 | |
---|
[2690] | 102 | DO it=1, nbtr_bin |
---|
| 103 | r_bin(it)=mdw(it)/2. |
---|
| 104 | ENDDO |
---|
| 105 | |
---|
| 106 | !--set boundaries of size bins |
---|
| 107 | DO it=1, nbtr_bin |
---|
| 108 | IF (it.EQ.1) THEN |
---|
| 109 | r_upper(it)=sqrt(r_bin(it+1)*r_bin(it)) |
---|
| 110 | r_lower(it)=r_bin(it)**2./r_upper(it) |
---|
| 111 | ELSEIF (it.EQ.nbtr_bin) THEN |
---|
| 112 | r_lower(it)=sqrt(r_bin(it)*r_bin(it-1)) |
---|
| 113 | r_upper(it)=r_bin(it)**2./r_lower(it) |
---|
| 114 | ELSE |
---|
| 115 | r_lower(it)=sqrt(r_bin(it)*r_bin(it-1)) |
---|
| 116 | r_upper(it)=sqrt(r_bin(it+1)*r_bin(it)) |
---|
| 117 | ENDIF |
---|
| 118 | ENDDO |
---|
| 119 | |
---|
| 120 | IF (debutphy .and. is_mpi_root) THEN |
---|
| 121 | DO it=1, nbtr_bin |
---|
| 122 | PRINT *,'radius bin', it, ':', r_bin(it), '(from', r_lower(it), 'to', r_upper(it), ')' |
---|
| 123 | ENDDO |
---|
| 124 | ENDIF |
---|
| 125 | |
---|
| 126 | !--initialising logical is_strato from stratomask |
---|
| 127 | is_strato(:,:)=.FALSE. |
---|
[2695] | 128 | WHERE (stratomask.GT.0.5) is_strato=.TRUE. |
---|
[2690] | 129 | |
---|
| 130 | ! STRACOMP (H2O, P, t_seri -> aerosol composition (R2SO4)) |
---|
| 131 | ! H2SO4 mass fraction in aerosol (%) |
---|
| 132 | CALL stracomp(sh,t_seri,pplay) |
---|
| 133 | |
---|
| 134 | ! aerosol density (gr/cm3) |
---|
| 135 | CALL denh2sa(t_seri) |
---|
| 136 | |
---|
| 137 | ! compute factor for converting dry to wet radius (for every grid box) |
---|
| 138 | f_r_wet(:,:) = (dens_aer_dry/(DENSO4(:,:)*1000.)/(R2SO4(:,:)/100.))**(1./3.) |
---|
| 139 | |
---|
| 140 | !--calculate mass of air in every grid box |
---|
| 141 | DO ilon=1, klon |
---|
| 142 | DO ilev=1, klev |
---|
[2752] | 143 | m_air_gridbox(ilon,ilev)=(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG*cell_area(ilon) |
---|
[2690] | 144 | ENDDO |
---|
| 145 | ENDDO |
---|
| 146 | |
---|
[2752] | 147 | ! IF (debutphy) THEN |
---|
| 148 | ! CALL gather(tr_seri, tr_seri_glo) |
---|
| 149 | ! IF (MAXVAL(tr_seri_glo).LT.1.e-30) THEN |
---|
[2690] | 150 | !--initialising tracer concentrations to zero |
---|
[2752] | 151 | ! DO it=1, nbtr |
---|
| 152 | ! tr_seri(:,:,it)=0.0 |
---|
| 153 | ! ENDDO |
---|
| 154 | ! ENDIF |
---|
| 155 | ! ENDIF |
---|
[2690] | 156 | |
---|
[2752] | 157 | !--initialise emission diagnostics |
---|
| 158 | budg_emi_ocs(:)=0.0 |
---|
| 159 | budg_emi_so2(:)=0.0 |
---|
| 160 | budg_emi_h2so4(:)=0.0 |
---|
| 161 | budg_emi_part(:)=0.0 |
---|
| 162 | |
---|
[2690] | 163 | !--sulfur emission, depending on chosen scenario (flag_sulf_emit) |
---|
| 164 | SELECT CASE(flag_sulf_emit) |
---|
| 165 | |
---|
| 166 | CASE(0) ! background aerosol |
---|
| 167 | ! do nothing (no emission) |
---|
| 168 | |
---|
| 169 | CASE(1) ! volcanic eruption |
---|
| 170 | !--only emit on day of eruption |
---|
| 171 | ! stretch emission over one day of Pinatubo eruption |
---|
| 172 | IF (year_cur==year_emit_vol.AND.mth_cur==mth_emit_vol.AND.day_cur==day_emit_vol) THEN |
---|
| 173 | ! |
---|
| 174 | DO i=1,klon |
---|
| 175 | !Pinatubo eruption at 15.14N, 120.35E |
---|
[2704] | 176 | IF ( xlat(i).GE.xlat_vol-dlat .AND. xlat(i).LT.xlat_vol+dlat .AND. & |
---|
| 177 | xlon(i).GE.xlon_vol-dlon .AND. xlon(i).LT.xlon_vol+dlon ) THEN |
---|
[2752] | 178 | ! |
---|
| 179 | PRINT *,'coordinates of volcanic injection point=',xlat(i), xlon(i), day_cur, mth_cur, year_cur |
---|
[2690] | 180 | ! compute altLMDz |
---|
| 181 | altLMDz(:)=0.0 |
---|
| 182 | DO k=1, klev |
---|
[2752] | 183 | zrho=pplay(i,k)/t_seri(i,k)/RD !air density in kg/m3 |
---|
| 184 | zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG !mass of layer in kg |
---|
| 185 | zdz=zdm(k)/zrho !thickness of layer in m |
---|
| 186 | altLMDz(k+1)=altLMDz(k)+zdz !altitude of interface |
---|
[2690] | 187 | ENDDO |
---|
| 188 | !compute distribution of emission to vertical model layers (based on Gaussian peak in altitude) |
---|
| 189 | f_lay_sum=0.0 |
---|
| 190 | DO k=1, klev |
---|
| 191 | f_lay_emiss(k)=0.0 |
---|
| 192 | DO i_int=1, n_int_alt |
---|
[2699] | 193 | alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt) |
---|
[2690] | 194 | f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_vol)* & |
---|
[2699] | 195 | & exp(-0.5*((alt-altemiss_vol)/sigma_alt_vol)**2.)* & |
---|
| 196 | & (altLMDz(k+1)-altLMDz(k))/float(n_int_alt) |
---|
[2690] | 197 | ENDDO |
---|
| 198 | f_lay_sum=f_lay_sum+f_lay_emiss(k) |
---|
| 199 | ENDDO |
---|
| 200 | !correct for step integration error |
---|
| 201 | f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum |
---|
| 202 | !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss) |
---|
| 203 | !vertically distributed emission |
---|
| 204 | DO k=1, klev |
---|
[3411] | 205 | ! stretch emission over one day of Pinatubo eruption |
---|
| 206 | emission=m_aer_emiss_vol*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/1./86400. |
---|
[2752] | 207 | tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys |
---|
| 208 | budg_emi_so2(i)=budg_emi_so2(i)+emission*zdm(k)*mSatom/mSO2mol |
---|
[2690] | 209 | ENDDO |
---|
| 210 | ENDIF ! emission grid cell |
---|
| 211 | ENDDO ! klon loop |
---|
| 212 | ENDIF ! emission period |
---|
| 213 | |
---|
| 214 | CASE(2) ! stratospheric aerosol injections (SAI) |
---|
| 215 | ! |
---|
| 216 | DO i=1,klon |
---|
| 217 | ! SAI standard scenario with continuous emission from 1 grid point at the equator |
---|
| 218 | ! SAI emission on single month |
---|
| 219 | ! IF ((mth_cur==4 .AND. & |
---|
| 220 | ! SAI continuous emission o |
---|
[2704] | 221 | IF ( xlat(i).GE.xlat_sai-dlat .AND. xlat(i).LT.xlat_sai+dlat .AND. & |
---|
| 222 | & xlon(i).GE.xlon_sai-dlon .AND. xlon(i).LT.xlon_sai+dlon ) THEN |
---|
[2752] | 223 | ! |
---|
| 224 | PRINT *,'coordinates of SAI point=',xlat(i), xlon(i), day_cur, mth_cur, year_cur |
---|
[2690] | 225 | ! compute altLMDz |
---|
| 226 | altLMDz(:)=0.0 |
---|
| 227 | DO k=1, klev |
---|
[2752] | 228 | zrho=pplay(i,k)/t_seri(i,k)/RD !air density in kg/m3 |
---|
| 229 | zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG !mass of layer in kg |
---|
| 230 | zdz=zdm(k)/zrho !thickness of layer in m |
---|
| 231 | altLMDz(k+1)=altLMDz(k)+zdz !altitude of interface |
---|
[2690] | 232 | ENDDO |
---|
| 233 | !compute distribution of emission to vertical model layers (based on Gaussian peak in altitude) |
---|
| 234 | f_lay_sum=0.0 |
---|
| 235 | DO k=1, klev |
---|
| 236 | f_lay_emiss(k)=0.0 |
---|
| 237 | DO i_int=1, n_int_alt |
---|
[2699] | 238 | alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt) |
---|
| 239 | f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_sai)* & |
---|
| 240 | & exp(-0.5*((alt-altemiss_sai)/sigma_alt_sai)**2.)* & |
---|
| 241 | & (altLMDz(k+1)-altLMDz(k))/float(n_int_alt) |
---|
[2690] | 242 | ENDDO |
---|
| 243 | f_lay_sum=f_lay_sum+f_lay_emiss(k) |
---|
| 244 | ENDDO |
---|
| 245 | !correct for step integration error |
---|
| 246 | f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum |
---|
| 247 | !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss) |
---|
| 248 | !vertically distributed emission |
---|
| 249 | DO k=1, klev |
---|
[2752] | 250 | ! stretch emission over whole year (360d) |
---|
| 251 | emission=m_aer_emiss_sai*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/360./86400. |
---|
| 252 | tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys |
---|
| 253 | budg_emi_so2(i)=budg_emi_so2(i)+emission*zdm(k)*mSatom/mSO2mol |
---|
[2690] | 254 | ENDDO |
---|
| 255 | ! !emission as monodisperse particles with 0.1um dry radius (BIN21) |
---|
| 256 | ! !vertically distributed emission |
---|
| 257 | ! DO k=1, klev |
---|
[2752] | 258 | ! ! stretch emission over whole year (360d) |
---|
| 259 | ! emission=m_aer_emiss*(mH2SO4mol/mSatom)/m_part_dry(21)/m_air_gridbox(i,k)*f_lay_emiss(k)/360./86400 |
---|
| 260 | ! tr_seri(i,k,id_BIN01_strat+20)=tr_seri(i,k,id_BIN01_strat+20)+emission*pdtphys |
---|
| 261 | ! budg_emi_part(i)=budg_emi_part(i)+emission*zdm(k)*mSatom/mH2SO4mol |
---|
[2690] | 262 | ! ENDDO |
---|
| 263 | ENDIF ! emission grid cell |
---|
| 264 | ENDDO ! klon loop |
---|
| 265 | |
---|
| 266 | END SELECT ! emission scenario (flag_sulf_emit) |
---|
| 267 | |
---|
| 268 | !--read background concentrations of OCS and SO2 and lifetimes from input file |
---|
[2695] | 269 | !--update the variables defined in phys_local_var_mod |
---|
| 270 | CALL interp_sulf_input(debutphy,pdtphys,paprs,tr_seri) |
---|
[2690] | 271 | |
---|
| 272 | !--convert OCS to SO2 in the stratosphere |
---|
[2752] | 273 | CALL ocs_to_so2(pdtphys,tr_seri,t_seri,pplay,paprs,is_strato) |
---|
[2690] | 274 | |
---|
| 275 | !--convert SO2 to H2SO4 |
---|
[2752] | 276 | CALL so2_to_h2so4(pdtphys,tr_seri,t_seri,pplay,paprs,is_strato) |
---|
[2690] | 277 | |
---|
| 278 | !--common routine for nucleation and condensation/evaporation with adaptive timestep |
---|
| 279 | CALL micphy_tstep(pdtphys,tr_seri,t_seri,pplay,paprs,rh,is_strato) |
---|
| 280 | |
---|
| 281 | !--call coagulation routine |
---|
| 282 | CALL coagulate(pdtphys,mdw,tr_seri,t_seri,pplay,dens_aer,is_strato) |
---|
| 283 | |
---|
| 284 | !--call sedimentation routine |
---|
| 285 | CALL aer_sedimnt(pdtphys, t_seri, pplay, paprs, tr_seri, dens_aer) |
---|
| 286 | |
---|
| 287 | !--compute mass concentration of PM2.5 sulfate particles (wet diameter and mass) at the surface for health studies |
---|
| 288 | surf_PM25_sulf(:)=0.0 |
---|
| 289 | DO i=1,klon |
---|
| 290 | DO it=1, nbtr_bin |
---|
| 291 | IF (mdw(it) .LT. 2.5e-6) THEN |
---|
| 292 | !surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas)*m_part(i,1,it) & |
---|
| 293 | !assume that particles consist of ammonium sulfate at the surface (132g/mol) and are dry at T = 20 deg. C and 50 perc. humidity |
---|
| 294 | surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas) & |
---|
| 295 | & *132./98.*dens_aer_dry*4./3.*RPI*(mdw(it)/2.)**3 & |
---|
[2752] | 296 | & *pplay(i,1)/t_seri(i,1)/RD*1.e9 |
---|
[2690] | 297 | ENDIF |
---|
| 298 | ENDDO |
---|
| 299 | ENDDO |
---|
| 300 | |
---|
[2695] | 301 | ! CALL minmaxsimple(tr_seri(:,:,id_SO2_strat),0.0,0.0,'fin SO2') |
---|
| 302 | ! DO it=1, nbtr_bin |
---|
| 303 | ! CALL minmaxsimple(tr_seri(:,:,nbtr_sulgas+it),0.0,0.0,'fin bin ') |
---|
| 304 | ! ENDDO |
---|
| 305 | |
---|
[2690] | 306 | END SUBROUTINE traccoag |
---|
| 307 | |
---|
| 308 | END MODULE traccoag_mod |
---|