- Timestamp:
- Jul 24, 2024, 12:17:33 PM (2 months ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/Amaury_dev/libf/phylmd/StratAer/traccoag_mod.F90
r5110 r5111 1 1 MODULE traccoag_mod 2 2 3 ! This module calculates the concentration of aerosol particles in certain size bins4 ! considering coagulation and sedimentation.5 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 6 CONTAINS 7 7 8 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_mode15 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 16 USE dimphy 17 17 USE infotrac_phy, ONLY: nbtr_bin, nbtr_sulgas, nbtr, id_SO2_strat … … 19 19 USE geometry_mod, ONLY: cell_area, boundslat 20 20 USE lmdz_grid_phy 21 USE lmdz_phys_mpi_data, ONLY: 21 USE lmdz_phys_mpi_data, ONLY: is_mpi_root 22 22 USE lmdz_phys_para, only: gather, scatter 23 23 USE phys_cal_mod, ONLY: year_len, year_cur, mth_cur, day_cur, hour … … 27 27 USE print_control_mod, ONLY: lunout 28 28 USE strataer_local_var_mod 29 29 30 30 IMPLICIT NONE 31 31 32 ! Input argument33 !---------------34 REAL, INTENT(IN):: pdtphys ! Pas d'integration pour la physique (seconde)35 REAL, INTENT(IN):: gmtime ! Heure courante36 LOGICAL, INTENT(IN) :: debutphy ! le flag de l'initialisation de la physique37 INTEGER, INTENT(IN) :: julien ! Jour julien38 39 REAL, DIMENSION(klev),INTENT(IN):: presnivs! pressions approximat. des milieux couches (en PA)40 REAL, DIMENSION(klon),INTENT(IN) :: xlat ! latitudes pour chaque point41 REAL, DIMENSION(klon),INTENT(IN) :: xlon ! longitudes pour chaque point42 REAL, DIMENSION(klon),INTENT(IN):: pphis ! geopotentiel du sol43 REAL, DIMENSION(klon,klev),INTENT(IN):: pphi ! geopotentiel de chaque couche44 45 REAL, DIMENSION(klon,klev),INTENT(IN):: t_seri ! Temperature46 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 specifique49 REAL, DIMENSION(klon,klev),INTENT(IN) :: rh ! humidite relative50 51 ! Output argument52 !----------------53 REAL, DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA]54 55 ! Local variables56 !----------------57 REAL 58 REAL 59 INTEGER 60 LOGICAL, DIMENSION(klon,klev):: is_strato ! true = above tropopause, false = below61 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 m64 REAL, DIMENSION(klev):: f_lay_emiss ! fraction of emission for every vertical layer65 REAL 66 REAL 67 INTEGER, PARAMETER :: n_int_alt=10 ! number of subintervals for integration over Gaussian emission profile68 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 73 REAL 74 REAL, DIMENSION(klev):: zdm ! mass of atm. model layer in kg75 REAL, DIMENSION(klon,klev):: dens_aer ! density of aerosol particles [kg/m3 aerosol] with default H2SO4 mass fraction76 REAL 77 REAL 78 REAL 79 REAL :: latmin,latmax,lonmin,lonmax ! lat/lon min/max for injection80 REAL 81 REAL :: pdt,stretchlong ! physic timestep, stretch emission over one day82 83 INTEGER 84 INTEGER 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 85 86 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,'-',hour88 WRITE(lunout,*) 'IN traccoag flag_emit: ',flag_emit87 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 89 ENDIF 90 90 91 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) 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) 104 124 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)) 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.) 107 135 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) 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 152 SELECT CASE(flag_emit) 153 153 … … 158 158 !--only emit on day of eruption 159 159 ! stretch emission over one day of Pinatubo eruption 160 DO ieru=1, nErupt161 162 day_cur>=day_emit_vol(ieru).AND.day_cur<(day_emit_vol(ieru)+injdur)) THEN163 164 165 m_aer=m_aer_emiss_vol(ieru,1)/(REAL(injdur)*REAL(ponde_lonlat_vol(ieru)))166 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=',ieru172 WRITE(lunout,*) 'IN traccoag, dlon=',dlon173 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 179 180 pdt=pdtphys181 182 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 189 190 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 191 CASE(2) ! stratospheric aerosol injections (SAI) 192 192 193 ! Computing duration of SAI in days...194 ! ... starting from 0...195 injdur_sai = 0196 ! ... then adding whole years from first to (n-1)th...197 DO yr = year_emit_sai_start, year_emit_sai_end-1198 ! (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_bissext201 ENDDO202 ! ... 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 -31207 CASE(3)208 injdur_sai = injdur_sai -31-28-is_bissext209 CASE(4)210 injdur_sai = injdur_sai -31-28-is_bissext-31211 CASE(5)212 injdur_sai = injdur_sai -31-28-is_bissext-31-30213 CASE(6)214 injdur_sai = injdur_sai -31-28-is_bissext-31-30-31215 CASE(7)216 injdur_sai = injdur_sai -31-28-is_bissext-31-30-31-30217 CASE(8)218 injdur_sai = injdur_sai -31-28-is_bissext-31-30-31-30-31219 CASE(9)220 injdur_sai = injdur_sai -31-28-is_bissext-31-30-31-30-31-31221 CASE(10)222 injdur_sai = injdur_sai -31-28-is_bissext-31-30-31-30-31-31-30223 CASE(11)224 injdur_sai = injdur_sai -31-28-is_bissext-31-30-31-30-31-31-30-31225 CASE(12)226 injdur_sai = injdur_sai -31-28-is_bissext-31-30-31-30-31-31-30-31-30227 END SELECT228 injdur_sai = injdur_sai-day_emit_sai_start+1229 ! ... then adding part of the n-th year230 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 +31234 CASE(3)235 injdur_sai = injdur_sai +31+28+is_bissext236 CASE(4)237 injdur_sai = injdur_sai +31+28+is_bissext+31238 CASE(5)239 injdur_sai = injdur_sai +31+28+is_bissext+31+30240 CASE(6)241 injdur_sai = injdur_sai +31+28+is_bissext+31+30+31242 CASE(7)243 injdur_sai = injdur_sai +31+28+is_bissext+31+30+31+30244 CASE(8)245 injdur_sai = injdur_sai +31+28+is_bissext+31+30+31+30+31246 CASE(9)247 injdur_sai = injdur_sai +31+28+is_bissext+31+30+31+30+31+31248 CASE(10)249 injdur_sai = injdur_sai +31+28+is_bissext+31+30+31+30+31+31+30250 CASE(11)251 injdur_sai = injdur_sai +31+28+is_bissext+31+30+31+30+31+31+30+31252 CASE(12)253 injdur_sai = injdur_sai +31+28+is_bissext+31+30+31+30+31+31+30+31+30254 END SELECT255 injdur_sai = injdur_sai+day_emit_sai_end256 ! A security: are SAI dates of injection consistent?257 IF (injdur_sai <= 0) THEN258 CALL abort_physic('traccoag_mod', 'Pb in SAI dates of injection.', 1)259 ENDIF260 ! Injection in itself261 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)) THEN267 268 m_aer=m_aer_emiss_sai269 !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)270 m_aer=m_aer*(mSO2mol/mSatom)271 272 latmin=xlat_sai273 latmax=xlat_sai274 lonmin=xlon_sai275 lonmax=xlon_sai276 altemiss = altemiss_sai277 sigma_alt = sigma_alt_sai278 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/mSO2mol287 ENDIF ! Condition over injection dates193 ! 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 288 289 289 CASE(3) ! --- SAI injection over a single band of longitude and between 290 291 292 m_aer=m_aer_emiss_sai293 294 m_aer=m_aer*(mSO2mol/mSatom)295 296 latmin=xlat_min_sai297 latmax=xlat_max_sai298 lonmin=xlon_sai299 lonmax=xlon_sai300 301 302 pdt=0.303 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/mSO2mol311 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 312 END SELECT ! emission scenario (flag_emit) 313 313 314 !--read background concentrations of OCS and SO2 and lifetimes from input file315 !--update the variables defined in phys_local_var_mod316 CALL interp_sulf_input(debutphy, pdtphys,paprs,tr_seri)317 318 !--convert OCS to SO2 in the stratosphere319 CALL ocs_to_so2(pdtphys, tr_seri,t_seri,pplay,paprs,is_strato)320 321 !--convert SO2 to H2SO4322 CALL so2_to_h2so4(pdtphys, tr_seri,t_seri,pplay,paprs,is_strato)323 324 !--common routine for nucleation and condensation/evaporation with adaptive timestep325 CALL micphy_tstep(pdtphys, tr_seri,t_seri,pplay,paprs,rh,is_strato)326 327 !--CALL coagulation routine328 CALL coagulate(pdtphys, mdw,tr_seri,t_seri,pplay,dens_aer,is_strato)329 330 !--CALL sedimentation routine314 !--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 331 CALL aer_sedimnt(pdtphys, t_seri, pplay, paprs, tr_seri, dens_aer) 332 332 333 !--compute mass concentration of PM2.5 sulfate particles (wet diameter and mass) at the surface for health studies334 surf_PM25_sulf(:) =0.0335 DO i =1,klon336 DO it =1, nbtr_bin333 !--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 337 IF (mdw(it) < 2.5e-6) THEN 338 338 !surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas)*m_part(i,1,it) & 339 339 !assume that particles consist of ammonium sulfate at the surface (132g/mol) 340 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.e9341 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 344 ENDIF 345 345 ENDDO 346 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 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 391 390 ENDDO 392 ENDDO 393 391 ENDDO 392 ENDDO 393 394 394 END SUBROUTINE traccoag 395 395
Note: See TracChangeset
for help on using the changeset viewer.
