- Timestamp:
- Nov 21, 2019, 4:43:45 PM (5 years ago)
- Location:
- LMDZ6/branches/Ocean_skin
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
libf/phylmd/StratAer/traccoag_mod.F90 (modified) (11 diffs, 1 prop)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/Ocean_skin
-
LMDZ6/branches/Ocean_skin/libf/phylmd/StratAer/traccoag_mod.F90
-
Property
svn:keywords
set to
Id
r3114 r3605 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 INTEGER :: it, k, i, ilon, ilev, itime, i_int, ieru 76 62 LOGICAL,DIMENSION(klon,klev) :: is_strato ! true = above tropopause, false = below 77 63 REAL,DIMENSION(klon,klev) :: m_air_gridbox ! mass of air in every grid box [kg] … … 90 76 REAL,DIMENSION(klev) :: zdm ! mass of atm. model layer in kg 91 77 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 78 REAL :: emission ! emission 79 REAL :: theta_min, theta_max ! for SAI computation between two latitudes 80 REAL :: dlat_loc 94 81 95 82 IF (is_mpi_root) THEN 96 PRINT *,'in traccoag: date from phys_cal_mod =',year_cur,'-',mth_cur,'-',day_cur,'-',hour 83 WRITE(lunout,*) 'in traccoag: date from phys_cal_mod =',year_cur,'-',mth_cur,'-',day_cur,'-',hour 84 WRITE(lunout,*) 'IN traccoag flag_sulf_emit: ',flag_sulf_emit 97 85 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 86 102 87 DO it=1, nbtr_bin 103 88 r_bin(it)=mdw(it)/2. … … 120 105 IF (debutphy .and. is_mpi_root) THEN 121 106 DO it=1, nbtr_bin 122 PRINT *,'radius bin', it, ':', r_bin(it), '(from', r_lower(it), 'to', r_upper(it), ')'107 WRITE(lunout,*) 'radius bin', it, ':', r_bin(it), '(from', r_lower(it), 'to', r_upper(it), ')' 123 108 ENDDO 124 109 ENDIF … … 170 155 !--only emit on day of eruption 171 156 ! 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 157 DO ieru=1, nErupt 158 IF (year_cur==year_emit_vol(ieru).AND.mth_cur==mth_emit_vol(ieru).AND.& 159 day_cur>=day_emit_vol(ieru).AND.day_cur<(day_emit_vol(ieru)+injdur)) THEN 160 ! 161 ! daily injection mass emission - NL 162 m_aer_emiss_vol_daily = m_aer_emiss_vol(ieru)/(REAL(injdur)*REAL(ponde_lonlat_vol(ieru))) 163 WRITE(lunout,*) 'IN traccoag DD m_aer_emiss_vol(ieru)=',m_aer_emiss_vol(ieru), & 164 'ponde_lonlat_vol(ieru)=',ponde_lonlat_vol(ieru),'(injdur*ponde_lonlat_vol(ieru))', & 165 (injdur*ponde_lonlat_vol(ieru)),'m_aer_emiss_vol_daily=',m_aer_emiss_vol_daily,'ieru=',ieru 166 WRITE(lunout,*) 'IN traccoag, dlon=',dlon 167 DO i=1,klon 168 !Pinatubo eruption at 15.14N, 120.35E 169 dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes 170 WRITE(lunout,*) 'IN traccoag, dlat=',dlat_loc 171 IF ( xlat(i).GE.xlat_min_vol(ieru)-dlat_loc .AND. xlat(i).LT.xlat_max_vol(ieru)+dlat_loc .AND. & 172 xlon(i).GE.xlon_min_vol(ieru)-dlon .AND. xlon(i).LT.xlon_max_vol(ieru)+dlon ) THEN 173 ! 174 WRITE(lunout,*) 'coordinates of volcanic injection point=',xlat(i),xlon(i),day_cur,mth_cur,year_cur 175 WRITE(lunout,*) 'DD m_aer_emiss_vol_daily=',m_aer_emiss_vol_daily 176 ! compute altLMDz 177 altLMDz(:)=0.0 178 DO k=1, klev 179 zrho=pplay(i,k)/t_seri(i,k)/RD !air density in kg/m3 180 zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG !mass of layer in kg 181 zdz=zdm(k)/zrho !thickness of layer in m 182 altLMDz(k+1)=altLMDz(k)+zdz !altitude of interface 183 ENDDO 184 185 SELECT CASE(flag_sulf_emit_distrib) 186 187 CASE(0) ! Gaussian distribution 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(ieru))* & 195 & exp(-0.5*((alt-altemiss_vol(ieru))/sigma_alt_vol(ieru))**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 201 CASE(1) ! Uniform distribution 202 ! In this case, parameter sigma_alt_vol(ieru) is considered to be half the 203 ! height of the injection, centered around altemiss_vol(ieru) 204 DO k=1, klev 205 f_lay_emiss(k)=max(min(altemiss_vol(ieru)+sigma_alt_vol(ieru),altLMDz(k+1))- & 206 & max(altemiss_vol(ieru)-sigma_alt_vol(ieru),altLMDz(k)),0.)/(2.*sigma_alt_vol(ieru)) 207 f_lay_sum=f_lay_sum+f_lay_emiss(k) 208 ENDDO 209 210 END SELECT ! End CASE over flag_sulf_emit_distrib) 211 212 WRITE(lunout,*) "IN traccoag m_aer_emiss_vol=",m_aer_emiss_vol(ieru) 213 WRITE(lunout,*) "IN traccoag f_lay_emiss=",f_lay_emiss 214 !correct for step integration error 215 f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum 216 !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss) 217 !vertically distributed emission 218 DO k=1, klev 219 ! stretch emission over one day of Pinatubo eruption 220 emission=m_aer_emiss_vol_daily*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/1./(86400.-pdtphys) 221 tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys 222 budg_emi_so2(i)=budg_emi_so2(i)+emission*zdm(k)*mSatom/mSO2mol 223 ENDDO 224 ENDIF ! emission grid cell 225 ENDDO ! klon loop 226 WRITE(lunout,*) "IN traccoag (ieru=",ieru,") m_aer_emiss_vol_daily=",m_aer_emiss_vol_daily 227 ENDIF ! emission period 228 ENDDO ! eruption number 229 214 230 CASE(2) ! stratospheric aerosol injections (SAI) 215 231 ! … … 217 233 ! SAI standard scenario with continuous emission from 1 grid point at the equator 218 234 ! SAI emission on single month 219 ! IF ((mth_cur==4 .AND. &220 235 ! SAI continuous emission o 221 IF ( xlat(i).GE.xlat_sai-dlat .AND. xlat(i).LT.xlat_sai+dlat .AND. & 236 dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes 237 IF ( xlat(i).GE.xlat_sai-dlat_loc .AND. xlat(i).LT.xlat_sai+dlat_loc .AND. & 222 238 & xlon(i).GE.xlon_sai-dlon .AND. xlon(i).LT.xlon_sai+dlon ) THEN 223 239 ! 224 PRINT *,'coordinates of SAI point=',xlat(i), xlon(i), day_cur, mth_cur, year_cur225 240 ! compute altLMDz 226 241 altLMDz(:)=0.0 … … 231 246 altLMDz(k+1)=altLMDz(k)+zdz !altitude of interface 232 247 ENDDO 248 249 SELECT CASE(flag_sulf_emit_distrib) 250 251 CASE(0) ! Gaussian distribution 233 252 !compute distribution of emission to vertical model layers (based on Gaussian peak in altitude) 234 253 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 254 DO k=1, klev 255 f_lay_emiss(k)=0.0 256 DO i_int=1, n_int_alt 257 alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt) 258 f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_sai)* & 259 & exp(-0.5*((alt-altemiss_sai)/sigma_alt_sai)**2.)* & 260 & (altLMDz(k+1)-altLMDz(k))/float(n_int_alt) 261 ENDDO 262 f_lay_sum=f_lay_sum+f_lay_emiss(k) 263 ENDDO 264 265 CASE(1) ! Uniform distribution 266 f_lay_sum=0.0 267 ! In this case, parameter sigma_alt_vol(ieru) is considered to be half 268 ! the height of the injection, centered around altemiss_sai 269 DO k=1, klev 270 f_lay_emiss(k)=max(min(altemiss_sai+sigma_alt_sai,altLMDz(k+1))- & 271 & max(altemiss_sai-sigma_alt_sai,altLMDz(k)),0.)/(2.*sigma_alt_sai) 272 f_lay_sum=f_lay_sum+f_lay_emiss(k) 273 ENDDO 274 275 END SELECT ! Gaussian or uniform distribution 276 245 277 !correct for step integration error 246 278 f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum … … 249 281 DO k=1, klev 250 282 ! 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.283 emission=m_aer_emiss_sai*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/FLOAT(year_len)/86400. 252 284 tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys 253 285 budg_emi_so2(i)=budg_emi_so2(i)+emission*zdm(k)*mSatom/mSO2mol 254 286 ENDDO 287 255 288 ! !emission as monodisperse particles with 0.1um dry radius (BIN21) 256 289 ! !vertically distributed emission 257 290 ! DO k=1, klev 258 291 ! ! 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 292 ! emission=m_aer_emiss*(mH2SO4mol/mSatom)/m_part_dry(21)/m_air_gridbox(i,k)*f_lay_emiss(k)/FLOAT(year_len)/86400. 293 ! tr_seri(i,k,id_BIN01_strat+20)=tr_seri(i,k,id_BIN01_strat+20)+emission*pdtphys 294 ! budg_emi_part(i)=budg_emi_part(i)+emission*zdm(k)*mSatom/mH2SO4mol 295 ! ENDDO 296 ENDIF ! emission grid cell 297 ENDDO ! klon loop 298 299 CASE(3) ! --- SAI injection over a single band of longitude and between 300 ! lat_min and lat_max 301 302 DO i=1,klon 303 ! SAI scenario with continuous emission 304 dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes 305 theta_min = max(xlat(i)-dlat_loc,xlat_min_sai) 306 theta_max = min(xlat(i)+dlat_loc,xlat_max_sai) 307 IF ( xlat(i).GE.xlat_min_sai-dlat_loc .AND. xlat(i).LT.xlat_max_sai+dlat_loc .AND. & 308 & xlon(i).GE.xlon_sai-dlon .AND. xlon(i).LT.xlon_sai+dlon ) THEN 309 ! 310 ! compute altLMDz 311 altLMDz(:)=0.0 312 DO k=1, klev 313 zrho=pplay(i,k)/t_seri(i,k)/RD !air density in kg/m3 314 zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG !mass of layer in kg 315 zdz=zdm(k)/zrho !thickness of layer in m 316 altLMDz(k+1)=altLMDz(k)+zdz !altitude of interface 317 ENDDO 318 319 SELECT CASE(flag_sulf_emit_distrib) 320 321 CASE(0) ! Gaussian distribution 322 !compute distribution of emission to vertical model layers (based on 323 !Gaussian peak in altitude) 324 f_lay_sum=0.0 325 DO k=1, klev 326 f_lay_emiss(k)=0.0 327 DO i_int=1, n_int_alt 328 alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt) 329 f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_sai)* & 330 & exp(-0.5*((alt-altemiss_sai)/sigma_alt_sai)**2.)* & 331 & (altLMDz(k+1)-altLMDz(k))/float(n_int_alt) 332 ENDDO 333 f_lay_sum=f_lay_sum+f_lay_emiss(k) 334 ENDDO 335 336 CASE(1) ! Uniform distribution 337 f_lay_sum=0.0 338 ! In this case, parameter sigma_alt_vol(ieru) is considered to be half 339 ! the height of the injection, centered around altemiss_sai 340 DO k=1, klev 341 f_lay_emiss(k)=max(min(altemiss_sai+sigma_alt_sai,altLMDz(k+1))- & 342 & max(altemiss_sai-sigma_alt_sai,altLMDz(k)),0.)/(2.*sigma_alt_sai) 343 f_lay_sum=f_lay_sum+f_lay_emiss(k) 344 ENDDO 345 346 END SELECT ! Gaussian or uniform distribution 347 348 !correct for step integration error 349 f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum 350 !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss) 351 !vertically distributed emission 352 DO k=1, klev 353 ! stretch emission over whole year (360d) 354 emission=m_aer_emiss_sai*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/ & 355 & FLOAT(year_len)/86400.*(sin(theta_max/180.*RPI)-sin(theta_min/180.*RPI))/ & 356 & (sin(xlat_max_sai/180.*RPI)-sin(xlat_min_sai/180.*RPI)) 357 tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys 358 budg_emi_so2(i)=budg_emi_so2(i)+emission*zdm(k)*mSatom/mSO2mol 359 ENDDO 360 361 ! !emission as monodisperse particles with 0.1um dry radius (BIN21) 362 ! !vertically distributed emission 363 ! DO k=1, klev 364 ! ! stretch emission over whole year (360d) 365 ! emission=m_aer_emiss*(mH2SO4mol/mSatom)/m_part_dry(21)/m_air_gridbox(i,k)*f_lay_emiss(k)/year_len/86400 260 366 ! tr_seri(i,k,id_BIN01_strat+20)=tr_seri(i,k,id_BIN01_strat+20)+emission*pdtphys 261 367 ! budg_emi_part(i)=budg_emi_part(i)+emission*zdm(k)*mSatom/mH2SO4mol … … 291 397 IF (mdw(it) .LT. 2.5e-6) THEN 292 398 !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 399 !assume that particles consist of ammonium sulfate at the surface (132g/mol) 400 !and are dry at T = 20 deg. C and 50 perc. humidity 294 401 surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas) & 295 402 & *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.
