Changeset 3526 for LMDZ6/trunk/libf/phylmd/StratAer
- Timestamp:
- May 28, 2019, 3:00:44 PM (6 years ago)
- Location:
- LMDZ6/trunk/libf/phylmd/StratAer
- Files:
-
- 1 added
- 14 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/phylmd/StratAer/aer_sedimnt.F90
-
Property
svn:keywords
set to
Id
r2752 r3526 1 ! 2 ! $Id$ 3 ! 1 4 SUBROUTINE AER_SEDIMNT(pdtphys, t_seri, pplay, paprs, tr_seri, dens_aer) 2 5 -
Property
svn:keywords
set to
-
LMDZ6/trunk/libf/phylmd/StratAer/aerophys.F90
-
Property
svn:keywords
set to
Id
-
Property
svn:keywords
set to
-
LMDZ6/trunk/libf/phylmd/StratAer/calcaerosolstrato_rrtm.F90
-
Property
svn:keywords
set to
Id
r2715 r3526 1 ! 2 ! $Id$ 3 ! 1 4 SUBROUTINE calcaerosolstrato_rrtm(pplay,t_seri,paprs,debut) 2 5 -
Property
svn:keywords
set to
-
LMDZ6/trunk/libf/phylmd/StratAer/coagulate.F90
-
Property
svn:keywords
set to
Id
r2950 r3526 1 ! 2 ! $Id$ 3 ! 1 4 SUBROUTINE COAGULATE(pdtcoag,mdw,tr_seri,t_seri,pplay,dens_aer,is_strato) 2 5 ! ----------------------------------------------------------------------- -
Property
svn:keywords
set to
-
LMDZ6/trunk/libf/phylmd/StratAer/cond_evap_tstep_mod.F90
-
Property
svn:keywords
set to
Id
r2695 r3526 1 ! 2 ! $Id$ 3 ! 1 4 MODULE cond_evap_tstep_mod 2 5 -
Property
svn:keywords
set to
-
LMDZ6/trunk/libf/phylmd/StratAer/interp_sulf_input.F90
-
Property
svn:keywords
set to
Id
r3097 r3526 1 ! 2 ! $Id$ 3 ! 1 4 SUBROUTINE interp_sulf_input(debutphy,pdtphys,paprs,tr_seri) 2 5 -
Property
svn:keywords
set to
-
LMDZ6/trunk/libf/phylmd/StratAer/micphy_tstep.F90
-
Property
svn:keywords
set to
Id
r3098 r3526 1 ! 2 ! $Id$ 3 ! 1 4 SUBROUTINE micphy_tstep(pdtphys,tr_seri,t_seri,pplay,paprs,rh,is_strato) 2 5 6 USE geometry_mod, ONLY : latitude_deg !NL- latitude corr. to local domain 3 7 USE dimphy, ONLY : klon,klev 4 8 USE aerophys … … 9 13 USE sulfate_aer_mod, ONLY : STRAACT 10 14 USE YOMCST, ONLY : RPI, RD, RG 11 15 USE print_control_mod, ONLY: lunout 16 USE strataer_mod 17 12 18 IMPLICIT NONE 13 19 … … 89 95 ! compute nucleation rate in kg(H2SO4)/kgA/s 90 96 CALL nucleation_rate(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev),rh(ilon,ilev), & 91 & a_xm,b_xm,c_xm,nucl_rate,ntot,x) 97 & a_xm,b_xm,c_xm,nucl_rate,ntot,x) 98 !NL - add nucleation box (if flag on) 99 IF (flag_nuc_rate_box) THEN 100 IF (latitude_deg(ilon).LE.(nuclat_min) .OR. latitude_deg(ilon).GE.(nuclat_max) & 101 .OR. pplay(ilon,ilev).GE.nucpres_max .AND. pplay(ilon,ilev) .LE. nucpres_min ) THEN 102 nucl_rate=0.0 103 ENDIF 104 ENDIF 92 105 ! compute cond/evap rate in kg(H2SO4)/kgA/s 93 106 CALL condens_evapor_rate(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev), & … … 160 173 DO it=1, nbtr 161 174 IF (tr_seri(ilon,ilev,it).LT.0.0) THEN 162 PRINT *,'micphy_tstep: negative concentration', tr_seri(ilon,ilev,it), ilon, ilev, it175 WRITE(lunout,*) 'micphy_tstep: negative concentration', tr_seri(ilon,ilev,it), ilon, ilev, it 163 176 ENDIF 164 177 ENDDO -
Property
svn:keywords
set to
-
LMDZ6/trunk/libf/phylmd/StratAer/miecalc_aer.F90
-
Property
svn:keywords
set to
Id
r2948 r3526 1 ! 2 ! $Id$ 3 ! 1 4 SUBROUTINE MIECALC_AER(tau_strat, piz_strat, cg_strat, tau_strat_wave, tau_lw_abs_rrtm, paprs, debut) 2 5 -
Property
svn:keywords
set to
-
LMDZ6/trunk/libf/phylmd/StratAer/minmaxsimple.F90
-
Property
svn:keywords
set to
Id
r2690 r3526 1 1 ! 2 ! $Id : minmaxsimple.F90 1910 2013-11-29 08:40:25Z fairhead$2 ! $Id$ 3 3 ! 4 4 SUBROUTINE minmaxsimple(zq,qmin,qmax,comment) -
Property
svn:keywords
set to
-
LMDZ6/trunk/libf/phylmd/StratAer/nucleation_tstep_mod.F90
-
Property
svn:keywords
set to
Id
r2695 r3526 1 ! 2 ! $Id$ 3 ! 1 4 MODULE nucleation_tstep_mod 2 5 -
Property
svn:keywords
set to
-
LMDZ6/trunk/libf/phylmd/StratAer/ocs_to_so2.F90
-
Property
svn:keywords
set to
Id
r2752 r3526 1 ! 2 ! $Id$ 3 ! 1 4 SUBROUTINE ocs_to_so2(pdtphys,tr_seri,t_seri,pplay,paprs,is_strato) 2 5 -
Property
svn:keywords
set to
-
LMDZ6/trunk/libf/phylmd/StratAer/so2_to_h2so4.F90
-
Property
svn:keywords
set to
Id
r2752 r3526 1 ! 2 ! $Id$ 3 ! 1 4 SUBROUTINE SO2_TO_H2SO4(pdtphys,tr_seri,t_seri,pplay,paprs,is_strato) 2 5 -
Property
svn:keywords
set to
-
LMDZ6/trunk/libf/phylmd/StratAer/sulfate_aer_mod.F90
-
Property
svn:keywords
set to
Id
r2690 r3526 1 ! 2 ! $Id$ 3 ! 1 4 MODULE sulfate_aer_mod 2 5 -
Property
svn:keywords
set to
-
LMDZ6/trunk/libf/phylmd/StratAer/traccoag_mod.F90
-
Property
svn:keywords
set to
Id
r3114 r3526 1 ! 2 ! $Id$ 3 ! 1 4 MODULE traccoag_mod 2 5 ! … … 16 19 USE infotrac 17 20 USE aerophys 18 USE geometry_mod, ONLY : cell_area 21 USE geometry_mod, ONLY : cell_area, boundslat 19 22 USE mod_grid_phy_lmdz 20 23 USE mod_phys_lmdz_mpi_data, ONLY : is_mpi_root … … 24 27 USE phys_local_var_mod, ONLY: stratomask 25 28 USE YOMCST 29 USE print_control_mod, ONLY: lunout 30 USE strataer_mod 31 USE phys_cal_mod, ONLY : year_len 26 32 27 33 IMPLICIT NONE … … 52 58 ! Local variables 53 59 !---------------- 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 71 REAL,PARAMETER :: xlat_sai=0.01 ! latitude of SAI in degree 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 60 REAL :: m_aer_emiss_vol_daily ! daily injection mass emission 61 REAL :: sum_emi_so2 ! Test sum of all LON for budg_emi_so2 62 INTEGER :: it, k, i, ilon, ilev, itime, i_int, ieru 76 63 LOGICAL,DIMENSION(klon,klev) :: is_strato ! true = above tropopause, false = below 77 64 REAL,DIMENSION(klon,klev) :: m_air_gridbox ! mass of air in every grid box [kg] … … 90 77 REAL,DIMENSION(klev) :: zdm ! mass of atm. model layer in kg 91 78 REAL,DIMENSION(klon,klev) :: dens_aer ! density of aerosol particles [kg/m3 aerosol] with default H2SO4 mass fraction 92 REAL :: dlat, dlon ! d latitude and d longitude of grid in degree93 79 REAL :: emission ! emission 80 REAL :: theta_min, theta_max ! for SAI computation between two latitudes 81 REAL :: dlat_loc 94 82 95 83 IF (is_mpi_root) THEN 96 PRINT *,'in traccoag: date from phys_cal_mod =',year_cur,'-',mth_cur,'-',day_cur,'-',hour 84 WRITE(lunout,*) 'in traccoag: date from phys_cal_mod =',year_cur,'-',mth_cur,'-',day_cur,'-',hour 85 WRITE(lunout,*) 'IN traccoag flag_sulf_emit: ',flag_sulf_emit 86 IF (flag_sulf_emit == 1) THEN 87 WRITE(lunout,*) 'IN traccoag nErupt: ',nErupt 88 WRITE(lunout,*) 'IN traccoag injdur: ',injdur 89 WRITE(lunout,*) 'IN traccoag : year_emit_vol',year_emit_vol 90 WRITE(lunout,*) 'IN traccoag : mth_emit_vol',mth_emit_vol 91 WRITE(lunout,*) 'IN traccoag : day_emit_vol',day_emit_vol 92 WRITE(lunout,*) 'IN traccoag : m_aer_emiss_vol',m_aer_emiss_vol 93 WRITE(lunout,*) 'IN traccoag : altemiss_vol',altemiss_vol 94 WRITE(lunout,*) 'IN traccoag : sigma_alt_vol',sigma_alt_vol 95 WRITE(lunout,*) 'IN traccoag : ponde_lonlat_vol',ponde_lonlat_vol 96 WRITE(lunout,*) 'IN traccoag : xlat_min_vol',xlat_min_vol 97 WRITE(lunout,*) 'IN traccoag : xlat_max_vol',xlat_max_vol 98 WRITE(lunout,*) 'IN traccoag : xlon_min_vol',xlon_min_vol 99 WRITE(lunout,*) 'IN traccoag : xlon_max_vol',xlon_max_vol 100 ELSEIF (flag_sulf_emit == 2) THEN 101 WRITE(lunout,*) 'IN traccoag : m_aer_emiss_sai',m_aer_emiss_sai 102 WRITE(lunout,*) 'IN traccoag : altemiss_sai',altemiss_sai 103 WRITE(lunout,*) 'IN traccoag : sigma_alt_sai',sigma_alt_sai 104 WRITE(lunout,*) 'IN traccoag : xlat_sai',xlat_sai 105 WRITE(lunout,*) 'IN traccoag : xlon_sai',xlon_sai 106 ELSEIF (flag_sulf_emit == 3) THEN 107 WRITE(lunout,*) 'IN traccoag : m_aer_emiss_sai',m_aer_emiss_sai 108 WRITE(lunout,*) 'IN traccoag : altemiss_sai',altemiss_sai 109 WRITE(lunout,*) 'IN traccoag : sigma_alt_sai',sigma_alt_sai 110 WRITE(lunout,*) 'IN traccoag : xlat_min_sai',xlat_min_sai 111 WRITE(lunout,*) 'IN traccoag : xlat_max_sai',xlat_max_sai 112 WRITE(lunout,*) 'IN traccoag : xlon_sai',xlon_sai 113 ENDIF 114 WRITE(lunout,*) 'IN traccoag : flag_nuc_rate_box = ',flag_nuc_rate_box 115 IF (flag_nuc_rate_box) THEN 116 WRITE(lunout,*) 'IN traccoag : nuclat_min = ',nuclat_min,', nuclat_max = ',nuclat_max 117 WRITE(lunout,*) 'IN traccoag : nucpres_min = ',nucpres_min,', nucpres_max = ',nucpres_max 118 ENDIF 97 119 ENDIF 98 99 dlat=180./2./FLOAT(nbp_lat) ! d latitude in degree 100 dlon=360./2./FLOAT(nbp_lon) ! d longitude in degree 101 120 102 121 DO it=1, nbtr_bin 103 122 r_bin(it)=mdw(it)/2. … … 120 139 IF (debutphy .and. is_mpi_root) THEN 121 140 DO it=1, nbtr_bin 122 PRINT *,'radius bin', it, ':', r_bin(it), '(from', r_lower(it), 'to', r_upper(it), ')'141 WRITE(lunout,*) 'radius bin', it, ':', r_bin(it), '(from', r_lower(it), 'to', r_upper(it), ')' 123 142 ENDDO 124 143 ENDIF … … 170 189 !--only emit on day of eruption 171 190 ! 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 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 178 ! 179 PRINT *,'coordinates of volcanic injection point=',xlat(i), xlon(i), day_cur, mth_cur, year_cur 180 ! compute altLMDz 181 altLMDz(:)=0.0 182 DO k=1, klev 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 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 193 alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt) 194 f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_vol)* & 195 & exp(-0.5*((alt-altemiss_vol)/sigma_alt_vol)**2.)* & 196 & (altLMDz(k+1)-altLMDz(k))/float(n_int_alt) 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 205 ! stretch emission over one day (minus one timestep) of Pinatubo eruption 206 emission=m_aer_emiss_vol*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/1./(86400.-pdtphys) 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 209 ENDDO 210 ENDIF ! emission grid cell 211 ENDDO ! klon loop 212 ENDIF ! emission period 213 191 DO ieru=1, nErupt 192 IF (is_mpi_root) THEN 193 sum_emi_so2 = 0.0 ! Init sum 194 ENDIF 195 IF (year_cur==year_emit_vol(ieru).AND.mth_cur==mth_emit_vol(ieru).AND.& 196 day_cur>=day_emit_vol(ieru).AND.day_cur<(day_emit_vol(ieru)+injdur)) THEN 197 ! 198 ! daily injection mass emission - NL 199 m_aer_emiss_vol_daily = m_aer_emiss_vol(ieru)/(REAL(injdur)*REAL(ponde_lonlat_vol(ieru))) 200 WRITE(lunout,*) 'IN traccoag DD m_aer_emiss_vol(ieru)=',m_aer_emiss_vol(ieru), & 201 'ponde_lonlat_vol(ieru)=',ponde_lonlat_vol(ieru),'(injdur*ponde_lonlat_vol(ieru))', & 202 (injdur*ponde_lonlat_vol(ieru)),'m_aer_emiss_vol_daily=',m_aer_emiss_vol_daily,'ieru=',ieru 203 WRITE(lunout,*) 'IN traccoag, dlon=',dlon 204 DO i=1,klon 205 !Pinatubo eruption at 15.14N, 120.35E 206 dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes 207 WRITE(lunout,*) 'IN traccoag, dlat=',dlat_loc 208 IF ( xlat(i).GE.xlat_min_vol(ieru)-dlat_loc .AND. xlat(i).LT.xlat_max_vol(ieru)+dlat_loc .AND. & 209 xlon(i).GE.xlon_min_vol(ieru)-dlon .AND. xlon(i).LT.xlon_max_vol(ieru)+dlon ) THEN 210 ! 211 WRITE(lunout,*) 'coordinates of volcanic injection point=',xlat(i),xlon(i),day_cur,mth_cur,year_cur 212 WRITE(lunout,*) 'DD m_aer_emiss_vol_daily=',m_aer_emiss_vol_daily 213 ! compute altLMDz 214 altLMDz(:)=0.0 215 DO k=1, klev 216 zrho=pplay(i,k)/t_seri(i,k)/RD !air density in kg/m3 217 zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG !mass of layer in kg 218 zdz=zdm(k)/zrho !thickness of layer in m 219 altLMDz(k+1)=altLMDz(k)+zdz !altitude of interface 220 ENDDO 221 222 SELECT CASE(flag_sulf_emit_distrib) 223 224 CASE(0) ! Gaussian distribution 225 !compute distribution of emission to vertical model layers (based on Gaussian peak in altitude) 226 f_lay_sum=0.0 227 DO k=1, klev 228 f_lay_emiss(k)=0.0 229 DO i_int=1, n_int_alt 230 alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt) 231 f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_vol(ieru))* & 232 & exp(-0.5*((alt-altemiss_vol(ieru))/sigma_alt_vol(ieru))**2.)* & 233 & (altLMDz(k+1)-altLMDz(k))/float(n_int_alt) 234 ENDDO 235 f_lay_sum=f_lay_sum+f_lay_emiss(k) 236 ENDDO 237 238 CASE(1) ! Uniform distribution 239 ! In this case, parameter sigma_alt_vol(ieru) is considered to be half the 240 ! height of the injection, centered around altemiss_vol(ieru) 241 DO k=1, klev 242 f_lay_emiss(k)=max(min(altemiss_vol(ieru)+sigma_alt_vol(ieru),altLMDz(k+1))- & 243 & max(altemiss_vol(ieru)-sigma_alt_vol(ieru),altLMDz(k)),0.)/(2.*sigma_alt_vol(ieru)) 244 f_lay_sum=f_lay_sum+f_lay_emiss(k) 245 ENDDO 246 247 END SELECT ! End CASE over flag_sulf_emit_distrib) 248 249 WRITE(lunout,*) "IN traccoag m_aer_emiss_vol=",m_aer_emiss_vol(ieru) 250 WRITE(lunout,*) "IN traccoag f_lay_emiss=",f_lay_emiss 251 !correct for step integration error 252 f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum 253 !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss) 254 !vertically distributed emission 255 DO k=1, klev 256 ! stretch emission over one day of Pinatubo eruption 257 emission=m_aer_emiss_vol_daily*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/1./(86400.-pdtphys) 258 tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys 259 budg_emi_so2(i)=budg_emi_so2(i)+emission*zdm(k)*mSatom/mSO2mol 260 ENDDO 261 sum_emi_so2 = sum_emi_so2 + budg_emi_so2(i) ! Sum all LON 262 ENDIF ! emission grid cell 263 ENDDO ! klon loop 264 WRITE(lunout,*) "IN traccoag (ieru=",ieru,") global sum_emi_so2=",sum_emi_so2 265 WRITE(lunout,*) "IN traccoag (ieru=",ieru,") m_aer_emiss_vol_daily=",m_aer_emiss_vol_daily 266 ENDIF ! emission period 267 ENDDO ! eruption number 268 214 269 CASE(2) ! stratospheric aerosol injections (SAI) 215 270 ! … … 217 272 ! SAI standard scenario with continuous emission from 1 grid point at the equator 218 273 ! SAI emission on single month 219 ! IF ((mth_cur==4 .AND. &220 274 ! SAI continuous emission o 221 IF ( xlat(i).GE.xlat_sai-dlat .AND. xlat(i).LT.xlat_sai+dlat .AND. & 275 dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes 276 WRITE(lunout,*) "IN traccoag, dlon=",dlon 277 WRITE(lunout,*) "IN traccoag, dlat=",dlat_loc 278 IF ( xlat(i).GE.xlat_sai-dlat_loc .AND. xlat(i).LT.xlat_sai+dlat_loc .AND. & 222 279 & xlon(i).GE.xlon_sai-dlon .AND. xlon(i).LT.xlon_sai+dlon ) THEN 223 280 ! 224 PRINT *,'coordinates of SAI point=',xlat(i), xlon(i), day_cur, mth_cur, year_cur281 WRITE(lunout,*) 'coordinates of SAI point=',xlat(i), xlon(i), day_cur, mth_cur, year_cur 225 282 ! compute altLMDz 226 283 altLMDz(:)=0.0 … … 231 288 altLMDz(k+1)=altLMDz(k)+zdz !altitude of interface 232 289 ENDDO 290 291 SELECT CASE(flag_sulf_emit_distrib) 292 293 CASE(0) ! Gaussian distribution 233 294 !compute distribution of emission to vertical model layers (based on Gaussian peak in altitude) 234 295 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 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) 242 ENDDO 243 f_lay_sum=f_lay_sum+f_lay_emiss(k) 244 ENDDO 296 DO k=1, klev 297 f_lay_emiss(k)=0.0 298 DO i_int=1, n_int_alt 299 alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt) 300 f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_sai)* & 301 & exp(-0.5*((alt-altemiss_sai)/sigma_alt_sai)**2.)* & 302 & (altLMDz(k+1)-altLMDz(k))/float(n_int_alt) 303 ENDDO 304 f_lay_sum=f_lay_sum+f_lay_emiss(k) 305 ENDDO 306 307 CASE(1) ! Uniform distribution 308 f_lay_sum=0.0 309 ! In this case, parameter sigma_alt_vol(ieru) is considered to be half 310 ! the height of the injection, centered around altemiss_sai 311 DO k=1, klev 312 f_lay_emiss(k)=max(min(altemiss_sai+sigma_alt_sai,altLMDz(k+1))- & 313 & max(altemiss_sai-sigma_alt_sai,altLMDz(k)),0.)/(2.*sigma_alt_sai) 314 f_lay_sum=f_lay_sum+f_lay_emiss(k) 315 ENDDO 316 317 END SELECT ! Gaussian or uniform distribution 318 245 319 !correct for step integration error 246 320 f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum … … 249 323 DO k=1, klev 250 324 ! 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.325 emission=m_aer_emiss_sai*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/year_len/86400. 252 326 tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys 253 327 budg_emi_so2(i)=budg_emi_so2(i)+emission*zdm(k)*mSatom/mSO2mol 254 328 ENDDO 329 255 330 ! !emission as monodisperse particles with 0.1um dry radius (BIN21) 256 331 ! !vertically distributed emission 257 332 ! DO k=1, klev 258 333 ! ! 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 334 ! emission=m_aer_emiss*(mH2SO4mol/mSatom)/m_part_dry(21)/m_air_gridbox(i,k)*f_lay_emiss(k)/year_len/86400 335 ! tr_seri(i,k,id_BIN01_strat+20)=tr_seri(i,k,id_BIN01_strat+20)+emission*pdtphys 336 ! budg_emi_part(i)=budg_emi_part(i)+emission*zdm(k)*mSatom/mH2SO4mol 337 ! ENDDO 338 ENDIF ! emission grid cell 339 ENDDO ! klon loop 340 341 CASE(3) ! --- SAI injection over a single band of longitude and between 342 ! lat_min and lat_max 343 344 WRITE(lunout,*) 'IN traccoag, dlon=',dlon 345 DO i=1,klon 346 ! SAI scenario with continuous emission 347 dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes 348 WRITE(lunout,*) 'IN traccoag, dlat = ',dlat_loc 349 theta_min = max(xlat(i)-dlat_loc,xlat_min_sai) 350 theta_max = min(xlat(i)+dlat_loc,xlat_max_sai) 351 IF ( xlat(i).GE.xlat_min_sai-dlat_loc .AND. xlat(i).LT.xlat_max_sai+dlat_loc .AND. & 352 & xlon(i).GE.xlon_sai-dlon .AND. xlon(i).LT.xlon_sai+dlon ) THEN 353 ! 354 ! compute altLMDz 355 altLMDz(:)=0.0 356 DO k=1, klev 357 zrho=pplay(i,k)/t_seri(i,k)/RD !air density in kg/m3 358 zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG !mass of layer in kg 359 zdz=zdm(k)/zrho !thickness of layer in m 360 altLMDz(k+1)=altLMDz(k)+zdz !altitude of interface 361 ENDDO 362 363 SELECT CASE(flag_sulf_emit_distrib) 364 365 CASE(0) ! Gaussian distribution 366 !compute distribution of emission to vertical model layers (based on 367 !Gaussian peak in altitude) 368 f_lay_sum=0.0 369 DO k=1, klev 370 f_lay_emiss(k)=0.0 371 DO i_int=1, n_int_alt 372 alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt) 373 f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_sai)* & 374 & exp(-0.5*((alt-altemiss_sai)/sigma_alt_sai)**2.)* & 375 & (altLMDz(k+1)-altLMDz(k))/float(n_int_alt) 376 ENDDO 377 f_lay_sum=f_lay_sum+f_lay_emiss(k) 378 ENDDO 379 380 CASE(1) ! Uniform distribution 381 f_lay_sum=0.0 382 ! In this case, parameter sigma_alt_vol(ieru) is considered to be half 383 ! the height of the injection, centered around altemiss_sai 384 DO k=1, klev 385 f_lay_emiss(k)=max(min(altemiss_sai+sigma_alt_sai,altLMDz(k+1))- & 386 & max(altemiss_sai-sigma_alt_sai,altLMDz(k)),0.)/(2.*sigma_alt_sai) 387 f_lay_sum=f_lay_sum+f_lay_emiss(k) 388 ENDDO 389 390 END SELECT ! Gaussian or uniform distribution 391 392 !correct for step integration error 393 f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum 394 !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss) 395 !vertically distributed emission 396 DO k=1, klev 397 ! stretch emission over whole year (360d) 398 emission=m_aer_emiss_sai*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/ & 399 & year_len/86400.*(sin(theta_max/180.*RPI)-sin(theta_min/180.*RPI))/ & 400 & (sin(xlat_max_sai/180.*RPI)-sin(xlat_min_sai/180.*RPI)) 401 tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys 402 budg_emi_so2(i)=budg_emi_so2(i)+emission*zdm(k)*mSatom/mSO2mol 403 ENDDO 404 405 ! !emission as monodisperse particles with 0.1um dry radius (BIN21) 406 ! !vertically distributed emission 407 ! DO k=1, klev 408 ! ! stretch emission over whole year (360d) 409 ! emission=m_aer_emiss*(mH2SO4mol/mSatom)/m_part_dry(21)/m_air_gridbox(i,k)*f_lay_emiss(k)/year_len/86400 260 410 ! tr_seri(i,k,id_BIN01_strat+20)=tr_seri(i,k,id_BIN01_strat+20)+emission*pdtphys 261 411 ! budg_emi_part(i)=budg_emi_part(i)+emission*zdm(k)*mSatom/mH2SO4mol … … 291 441 IF (mdw(it) .LT. 2.5e-6) THEN 292 442 !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 443 !assume that particles consist of ammonium sulfate at the surface (132g/mol) 444 !and are dry at T = 20 deg. C and 50 perc. humidity 294 445 surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas) & 295 446 & *132./98.*dens_aer_dry*4./3.*RPI*(mdw(it)/2.)**3 & -
Property
svn:keywords
set to
Note: See TracChangeset
for help on using the changeset viewer.