Changeset 4601 for LMDZ6/trunk/libf/phylmd/StratAer
- Timestamp:
- Jul 1, 2023, 12:07:30 AM (20 months ago)
- Location:
- LMDZ6/trunk/libf/phylmd/StratAer
- Files:
-
- 5 added
- 1 deleted
- 7 edited
-
aerophys.F90 (modified) (1 diff)
-
interp_sulf_input.F90 (modified) (14 diffs)
-
micphy_tstep.F90 (modified) (1 diff)
-
nucleation_tstep_mod.F90 (modified) (10 diffs)
-
ocs_to_so2.F90 (modified) (3 diffs)
-
so2_to_h2so4.F90 (modified) (3 diffs)
-
strataer_emiss_mod.F90 (added)
-
strataer_local_var_mod.F90 (added)
-
strataer_mod.F90 (deleted)
-
strataer_nuc_mod.F90 (added)
-
stratdistrib.F90 (added)
-
stratemit.F90 (added)
-
traccoag_mod.F90 (modified) (10 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/phylmd/StratAer/aerophys.F90
r3526 r4601 17 17 REAL, PARAMETER :: mSatom=32.06*1.66E-27 ! Mass of a S atom [kg] 18 18 REAL, PARAMETER :: mOCSmol=60.07*1.66E-27 ! Mass of an OCS molecule [kg] 19 REAL, PARAMETER :: mClatom=35.45*1.66E-27 ! Mass of an Cl atom [kg] 20 REAL, PARAMETER :: mHClmol=36.46*1.66E-27 ! Mass of an HCl molecule [kg] 21 REAL, PARAMETER :: mBratom=79.90*1.66E-27 ! Mass of an Br atom [kg] 22 REAL, PARAMETER :: mHBrmol=80.92*1.66E-27 ! Mass of an HBr molecule [kg] 23 REAL, PARAMETER :: mNOmol=30.01*1.66E-27 ! Mass of an NO molecule [kg] 24 REAL, PARAMETER :: mNO2mol=46.01*1.66E-27 ! Mass of an NO2 molecule [kg] 25 REAL, PARAMETER :: mNatome=14.0067*1.66E-27 ! Mass of an N atome [kg] 19 26 ! 20 27 END MODULE aerophys -
LMDZ6/trunk/libf/phylmd/StratAer/interp_sulf_input.F90
r3677 r4601 12 12 USE mod_phys_lmdz_omp_data, ONLY : is_omp_root 13 13 USE phys_local_var_mod, ONLY : budg_3D_backgr_ocs, budg_3D_backgr_so2 14 USE phys_local_var_mod, ONLY : OCS_lifetime, SO2_lifetime 14 USE phys_local_var_mod, ONLY : OCS_lifetime, SO2_lifetime, H2SO4_lifetime, O3_clim 15 15 USE mod_phys_lmdz_para 16 16 USE dimphy … … 19 19 USE aerophys 20 20 USE YOMCST 21 USE strataer_local_var_mod, ONLY : flag_newclim_file,flag_verbose_strataer 21 22 22 23 IMPLICIT NONE … … 62 63 REAL OCS_clim_glo(klon_glo,klev) 63 64 REAL SO2_clim_glo(klon_glo,klev) 65 66 REAL, ALLOCATABLE :: O3_clim_in(:, :, :, :) 67 REAL, ALLOCATABLE :: O3_clim_mth(:, :, :) 68 REAL, ALLOCATABLE :: O3_clim_tmp(:, :) 69 REAL O3_clim_glo(klon_glo,klev) 70 71 ! fixed climatos 64 72 REAL, ALLOCATABLE :: OCS_lifetime_in(:, :, :, :) 65 73 REAL, ALLOCATABLE :: SO2_lifetime_in(:, :, :, :) … … 70 78 REAL OCS_lifetime_glo(klon_glo,klev) 71 79 REAL SO2_lifetime_glo(klon_glo,klev) 80 81 REAL, ALLOCATABLE :: H2SO4_lifetime_in(:, :, :, :) 82 REAL, ALLOCATABLE :: H2SO4_lifetime_mth(:, :, :) 83 REAL, ALLOCATABLE :: H2SO4_lifetime_tmp(:, :) 84 85 REAL H2SO4_lifetime_glo(klon_glo,klev) 72 86 ! 73 87 REAL, ALLOCATABLE, SAVE :: OCS_clim(:,:) … … 79 93 80 94 ! For NetCDF: 81 INTEGER ncid_in ! IDs for input files 82 INTEGER varid, ncerr 95 INTEGER :: ncid_in ! IDs for input files 96 INTEGER :: varid, ncerr 97 CHARACTER (LEN=3) :: nc_lat, nc_lon 98 CHARACTER (LEN=256) :: nc_fname 83 99 84 100 INTEGER, PARAMETER :: lev_input=17 … … 103 119 IF (is_mpi_root.AND.is_omp_root) THEN 104 120 105 !--reading emission files 106 CALL nf95_open("ocs_so2_annual_lmdz.nc", nf90_nowrite, ncid_in) 107 121 !--init ncdf variables 122 IF(flag_newclim_file) THEN 123 nc_fname = "ocs_so2_h2so4_annual_lmdz.nc" 124 nc_lat = "LAT" 125 nc_lon = "LON" 126 ELSE 127 ! old file for retro compatibility 128 nc_fname = "ocs_so2_annual_lmdz.nc" 129 nc_lat = "lat" 130 nc_lon = "lon" 131 ENDIF 132 133 !--reading emission files 134 CALL nf95_open(nc_fname, nf90_nowrite, ncid_in) 135 IF(flag_verbose_strataer) print *, "Open file ", nc_fname 136 108 137 CALL nf95_inq_varid(ncid_in, "LEV", varid) 109 138 CALL nf95_gw_var(ncid_in, varid, lev) 110 139 n_lev = size(lev) 111 112 CALL nf95_inq_varid(ncid_in, "lat", varid)140 141 CALL nf95_inq_varid(ncid_in, nc_lat, varid) 113 142 CALL nf95_gw_var(ncid_in, varid, latitude) 114 143 n_lat = size(latitude) 115 116 CALL nf95_inq_varid(ncid_in, "lon", varid)144 145 CALL nf95_inq_varid(ncid_in, nc_lon, varid) 117 146 CALL nf95_gw_var(ncid_in, varid, longitude) 118 147 n_lon = size(longitude) 119 148 120 149 CALL nf95_inq_varid(ncid_in, "TIME", varid) 121 150 CALL nf95_gw_var(ncid_in, varid, time) 122 151 n_mth = size(time) 123 152 124 153 IF (.NOT.ALLOCATED(OCS_clim_in)) ALLOCATE(OCS_clim_in(n_lon, n_lat, n_lev, n_mth)) 125 154 IF (.NOT.ALLOCATED(SO2_clim_in)) ALLOCATE(SO2_clim_in(n_lon, n_lat, n_lev, n_mth)) 155 IF (.NOT.ALLOCATED(O3_clim_in)) ALLOCATE(O3_clim_in(n_lon, n_lat, n_lev, n_mth)) 156 126 157 IF (.NOT.ALLOCATED(OCS_lifetime_in)) ALLOCATE(OCS_lifetime_in(n_lon, n_lat, n_lev, n_mth)) 127 158 IF (.NOT.ALLOCATED(SO2_lifetime_in)) ALLOCATE(SO2_lifetime_in(n_lon, n_lat, n_lev, n_mth)) 128 159 IF (.NOT.ALLOCATED(H2SO4_lifetime_in)) ALLOCATE(H2SO4_lifetime_in(n_lon, n_lat, n_lev, n_mth)) 160 129 161 CALL nf95_inq_varid(ncid_in, "OCS", varid) 130 162 ncerr = nf90_get_var(ncid_in, varid, OCS_clim_in) 131 print *,'code erreur OCS=', ncerr, varid163 IF(flag_verbose_strataer) print *,'code erreur OCS=', ncerr, varid 132 164 133 165 CALL nf95_inq_varid(ncid_in, "SO2", varid) 134 166 ncerr = nf90_get_var(ncid_in, varid, SO2_clim_in) 135 print *,'code erreur SO2=', ncerr, varid167 IF(flag_verbose_strataer) print *,'code erreur SO2=', ncerr, varid 136 168 137 169 CALL nf95_inq_varid(ncid_in, "OCS_LIFET", varid) 138 170 ncerr = nf90_get_var(ncid_in, varid, OCS_lifetime_in) 139 print *,'code erreur OCS_lifetime_in=', ncerr, varid171 IF(flag_verbose_strataer) print *,'code erreur OCS_lifetime_in=', ncerr, varid 140 172 141 173 CALL nf95_inq_varid(ncid_in, "SO2_LIFET", varid) 142 174 ncerr = nf90_get_var(ncid_in, varid, SO2_lifetime_in) 143 print *,'code erreur SO2_lifetime_in=', ncerr, varid 144 175 IF(flag_verbose_strataer) print *,'code erreur SO2_lifetime_in=', ncerr, varid 176 177 IF(flag_newclim_file) THEN 178 CALL nf95_inq_varid(ncid_in, "O3", varid) 179 ncerr = nf90_get_var(ncid_in, varid, O3_clim_in) 180 IF(flag_verbose_strataer) print *,'code erreur O3=', ncerr, varid 181 182 CALL nf95_inq_varid(ncid_in, "H2SO4_LIFET", varid) 183 ncerr = nf90_get_var(ncid_in, varid, H2SO4_lifetime_in) 184 IF(flag_verbose_strataer) print *,'code erreur H2SO4_lifetime_in=', ncerr, varid 185 ENDIF 186 145 187 CALL nf95_close(ncid_in) 146 188 … … 149 191 IF (.NOT.ALLOCATED(OCS_clim_tmp)) ALLOCATE(OCS_clim_tmp(klon_glo, n_lev)) 150 192 IF (.NOT.ALLOCATED(SO2_clim_tmp)) ALLOCATE(SO2_clim_tmp(klon_glo, n_lev)) 193 IF (.NOT.ALLOCATED(O3_clim_mth)) ALLOCATE(O3_clim_mth(n_lon, n_lat, n_lev)) 194 IF (.NOT.ALLOCATED(O3_clim_tmp)) ALLOCATE(O3_clim_tmp(klon_glo, n_lev)) 195 151 196 IF (.NOT.ALLOCATED(OCS_lifetime_mth)) ALLOCATE(OCS_lifetime_mth(n_lon, n_lat, n_lev)) 152 197 IF (.NOT.ALLOCATED(SO2_lifetime_mth)) ALLOCATE(SO2_lifetime_mth(n_lon, n_lat, n_lev)) 153 198 IF (.NOT.ALLOCATED(OCS_lifetime_tmp)) ALLOCATE(OCS_lifetime_tmp(klon_glo, n_lev)) 154 199 IF (.NOT.ALLOCATED(SO2_lifetime_tmp)) ALLOCATE(SO2_lifetime_tmp(klon_glo, n_lev)) 155 200 IF (.NOT.ALLOCATED(H2SO4_lifetime_mth)) ALLOCATE(H2SO4_lifetime_mth(n_lon, n_lat, n_lev)) 201 IF (.NOT.ALLOCATED(H2SO4_lifetime_tmp)) ALLOCATE(H2SO4_lifetime_tmp(klon_glo, n_lev)) 202 156 203 !---select the correct month, undo multiplication with 1.e12 (precision reasons) 157 204 !---correct latitudinal order and convert input from volume mixing ratio to mass mixing ratio … … 161 208 SO2_lifetime_mth(:,j,:) = SO2_lifetime_in(:,n_lat+1-j,:,mth_cur) 162 209 OCS_lifetime_mth(:,j,:) = OCS_lifetime_in(:,n_lat+1-j,:,mth_cur) 210 211 ! O3 from 2d model is not tracer, in VMR 212 IF(flag_newclim_file) THEN 213 H2SO4_lifetime_mth(:,j,:) = H2SO4_lifetime_in(:,n_lat+1-j,:,mth_cur) 214 ! new input files 215 O3_clim_mth(:,j,:) = 1.e-6*O3_clim_in(:,n_lat+1-j,:,mth_cur) 216 ELSE 217 H2SO4_lifetime_mth(:,j,:) = 1.e-6 218 O3_clim_mth(:,j,:) = 1.e-6 219 ENDIF 163 220 ENDDO 164 221 … … 166 223 CALL grid2dTo1d_glo(OCS_clim_mth,OCS_clim_tmp) 167 224 CALL grid2dTo1d_glo(SO2_clim_mth,SO2_clim_tmp) 225 CALL grid2dTo1d_glo(O3_clim_mth,O3_clim_tmp) 226 168 227 CALL grid2dTo1d_glo(OCS_lifetime_mth,OCS_lifetime_tmp) 169 228 CALL grid2dTo1d_glo(SO2_lifetime_mth,SO2_lifetime_tmp) 229 CALL grid2dTo1d_glo(H2SO4_lifetime_mth,H2SO4_lifetime_tmp) 170 230 171 231 !--set lifetime to very high value in uninsolated areas 172 232 DO i=1, klon_glo 173 233 DO kk=1, n_lev 234 !--added tests on VMR of species for realism 235 IF (OCS_clim_tmp(i,kk).LT.1.e-20) THEN 236 OCS_clim_tmp(i,kk)=1.0e-20 237 ENDIF 238 IF (SO2_clim_tmp(i,kk).LT.1.e-20) THEN 239 SO2_clim_tmp(i,kk)=1.0e-20 240 ENDIF 241 ! 50 ppbv min for O3 242 IF (O3_clim_tmp(i,kk).LT.50.0e-9) THEN 243 O3_clim_tmp(i,kk)=50.0e-9 244 ENDIF 245 174 246 IF (OCS_lifetime_tmp(i,kk)==0.0) THEN 175 247 OCS_lifetime_tmp(i,kk)=1.0e12 … … 178 250 SO2_lifetime_tmp(i,kk)=1.0e12 179 251 ENDIF 252 IF (H2SO4_lifetime_tmp(i,kk)==0.0) THEN 253 H2SO4_lifetime_tmp(i,kk)=1.0e12 254 ENDIF 180 255 ENDDO 181 256 ENDDO … … 186 261 OCS_lifetime_glo(i,k)=0.0 187 262 SO2_lifetime_glo(i,k)=0.0 263 O3_clim_glo(i,k)=0.0 264 188 265 OCS_clim_glo(i,k)=0.0 189 266 SO2_clim_glo(i,k)=0.0 267 H2SO4_lifetime_glo(i,k)=0.0 268 190 269 DO kk=1, n_lev 191 270 OCS_lifetime_glo(i,k)=OCS_lifetime_glo(i,k)+ & … … 195 274 MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk))-MAX(paprs_glo(i,k+1),paprs_input(kk+1))) & 196 275 *SO2_lifetime_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1)) 276 IF(flag_newclim_file) THEN 277 H2SO4_lifetime_glo(i,k)=H2SO4_lifetime_glo(i,k)+ & 278 MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk)) & 279 -MAX(paprs_glo(i,k+1),paprs_input(kk+1))) & 280 *H2SO4_lifetime_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1)) 281 ENDIF 282 197 283 OCS_clim_glo(i,k)=OCS_clim_glo(i,k)+ & 198 284 MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk))-MAX(paprs_glo(i,k+1),paprs_input(kk+1))) & … … 201 287 MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk))-MAX(paprs_glo(i,k+1),paprs_input(kk+1))) & 202 288 *SO2_clim_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1)) 289 IF(flag_newclim_file) THEN 290 O3_clim_glo(i,k)=O3_clim_glo(i,k)+ & 291 MAX(0.0,MIN(paprs_glo(i,k),paprs_input(kk)) & 292 -MAX(paprs_glo(i,k+1),paprs_input(kk+1))) & 293 *O3_clim_tmp(i,kk)/(paprs_glo(i,k)-paprs_glo(i,k+1)) 294 ENDIF 203 295 ENDDO 204 296 ENDDO … … 213 305 CALL scatter(OCS_clim_glo, OCS_clim) 214 306 CALL scatter(SO2_clim_glo, SO2_clim) 307 CALL scatter(O3_clim_glo, O3_clim) 215 308 CALL scatter(OCS_lifetime_glo, OCS_lifetime) 216 309 CALL scatter(SO2_lifetime_glo, SO2_lifetime) 217 310 CALL scatter(H2SO4_lifetime_glo, H2SO4_lifetime) 311 218 312 IF (is_mpi_root.AND.is_omp_root) THEN 219 313 ! 220 DEALLOCATE(OCS_clim_in,SO2_clim_in )221 DEALLOCATE(OCS_clim_mth,SO2_clim_mth )222 DEALLOCATE(OCS_clim_tmp,SO2_clim_tmp )223 DEALLOCATE(OCS_lifetime_in, SO2_lifetime_in)224 DEALLOCATE(OCS_lifetime_mth, SO2_lifetime_mth)225 DEALLOCATE(OCS_lifetime_tmp, SO2_lifetime_tmp)314 DEALLOCATE(OCS_clim_in,SO2_clim_in,O3_clim_in) 315 DEALLOCATE(OCS_clim_mth,SO2_clim_mth,O3_clim_mth) 316 DEALLOCATE(OCS_clim_tmp,SO2_clim_tmp,O3_clim_tmp) 317 DEALLOCATE(OCS_lifetime_in, SO2_lifetime_in, H2SO4_lifetime_in) 318 DEALLOCATE(OCS_lifetime_mth, SO2_lifetime_mth, H2SO4_lifetime_mth) 319 DEALLOCATE(OCS_lifetime_tmp, SO2_lifetime_tmp, H2SO4_lifetime_tmp) 226 320 ! 227 321 ENDIF !--is_mpi_root and is_omp_root -
LMDZ6/trunk/libf/phylmd/StratAer/micphy_tstep.F90
r4293 r4601 14 14 USE YOMCST, ONLY : RPI, RD, RG 15 15 USE print_control_mod, ONLY: lunout 16 USE strataer_ mod16 USE strataer_local_var_mod 17 17 18 18 IMPLICIT NONE -
LMDZ6/trunk/libf/phylmd/StratAer/nucleation_tstep_mod.F90
r4293 r4601 11 11 USE infotrac_phy 12 12 USE YOMCST, ONLY : RPI, RD, RMD, RKBOL, RNAVO 13 13 USE strataer_local_var_mod, ONLY : flag_new_nucl 14 14 15 IMPLICIT NONE 15 16 16 17 ! input variables 17 LOGICAL, PARAMETER :: flag_new_nucl=.TRUE.18 18 REAL rhoa ! H2SO4 number density [molecules/cm3] 19 19 REAL t_seri ! temperature (K) … … 170 170 171 171 IF (t < 190.15) THEN 172 ! print *,'Warning (ilon=',ilon,'ilev=',ilev,'): temperature < 190.15 K, using 190.15 K (T=',t,'K)'173 172 t=190.15 174 173 ENDIF 175 174 176 175 IF (t > 300.15) THEN 177 ! print *,'Warning (ilon=',ilon,'ilev=',ilev,'): temperature > 300.15 K, using 300.15 K'178 176 t=300.15 179 177 ENDIF 180 178 181 179 IF (rh < 0.0001) THEN 182 ! print *,'Warning (ilon=',ilon,'ilev=',ilev,'): saturation ratio of water < 0.0001, using 0.0001'183 180 rh=0.0001 184 181 ENDIF 185 182 186 183 IF (rh > 1.0) THEN 187 ! print *,'Warning (ilon=',ilon,'ilev=',ilev,'): saturation ratio of water > 1 using 1'188 ! print *, 'rh=',rh189 184 rh=1.0 190 185 ENDIF 191 186 192 187 IF (rhoa < 1.e4) THEN 193 ! print *,'Warning (ilon=',ilon,'ilev=',ilev,'): sulfuric acid concentration < 1e4 1/cm3, using 1e4 1/cm3'194 188 rhoa=1.e4 195 189 ENDIF 196 190 197 191 IF (rhoa > 1.e11) THEN 198 ! print *,'Warning (ilon=',ilon,'ilev=',ilev,'): sulfuric acid concentration > 1e11 1/cm3, using 1e11 1/cm3'199 192 rhoa=1.e11 200 193 ENDIF … … 281 274 282 275 IF (jnuc < 1.E-7) THEN 283 ! print *,'Warning (ilon=',ilon,'ilev=',ilev'): nucleation rate < 1E-7/cm3s, using 0.0/cm3s,'284 276 jnuc=0.0 285 277 ENDIF 286 278 287 279 IF (jnuc > 1.e10) THEN 288 ! print *,'Warning (ilon=',ilon,'ilev=',ilev, &289 ! & '): nucleation rate > 1e10/cm3s, using 1e10/cm3s'290 280 jnuc=1.e10 291 281 ENDIF … … 425 415 !Temperature bounds 426 416 IF (t.LE.165.) THEN 427 ! print *,'Warning: temperature < 165.0 K, using 165.0 K in neutral nucleation calculation'428 417 tln=165.0 429 418 ENDIF 430 419 IF (t.LE.195.) THEN 431 ! print *,'Warning: temperature < 195.0 K, using 195.0 K in ion-induced nucleation calculation'432 420 tli=195.0 433 421 ENDIF 434 422 IF (t.GE.400.) THEN 435 ! print *,'Warning: temperature > 400. K, using 400. K in nucleation calculations'436 423 tln=400. 437 424 tli=400. … … 440 427 ! Saturation ratio bounds 441 428 IF (satrat.LT.1.E-7) THEN 442 ! print *,'Warning: saturation ratio of water < 1.E-7, using 1.E-7 in ion-induced nucleation calculation'443 429 satratli=1.E-7 444 430 ENDIF 445 431 IF (satrat.LT.1.E-5) THEN 446 ! print *,'Warning: saturation ratio of water < 1.E-5, using 1.E-5 in neutral nucleation calculation'447 432 satratln=1.E-5 448 433 ENDIF 449 434 IF (satrat.GT.0.95) THEN 450 ! print *,'Warning: saturation ratio of water > 0.95, using 0.95 in ion-induced nucleation calculation'451 435 satratli=0.95 452 436 ENDIF 453 437 IF (satrat.GT.1.0) THEN 454 ! print *,'Warning: saturation ratio of water > 1 using 1 in neutral nucleation calculation'455 438 satratln=1.0 456 439 ENDIF … … 458 441 ! Sulfuric acid concentration bounds 459 442 IF (rhoa.LE.1.E4) THEN 460 ! print *,'Warning: sulfuric acid < 1e4 1/cm3, using 1e4 1/cm3 in nucleation calculation'461 443 rhoaln=1.E4 462 444 rhoali=1.E4 463 445 ENDIF 464 446 IF (rhoa.GT.1.E13) THEN 465 ! print *,'Warning: sulfuric acid > 1e13 1/cm3, using 1e13 1/cm3 in neutral nucleation calculation'466 447 rhoaln=1.E13 467 448 ENDIF 468 449 IF (rhoa.GT.1.E16) THEN 469 ! print *,'Warning: sulfuric acid concentration > 1e16 1/cm3, using 1e16 1/cm3 in ion-induced nucleation calculation'470 450 rhoali=1.E16 471 451 ENDIF … … 521 501 ! Dimer formation rate 522 502 jnuc_n=1.E6*(2.*0.3E-9)**2.*SQRT(8.*RPI*RKBOL*(1./mH2SO4mol+1./mH2SO4mol))/2.*SQRT(t)*rhoa**2. 523 ! jnuc_n=1.E6*(2.*0.3E-9)**2.*SQRT(8.*3.141593*1.38E-23*(1./(1.661E-27*98.07)+1./(1.661E-27*98.07)))/2.*SQRT(t)*rhoa**2.524 503 ntot_n=1. !set to 1 525 504 na_n=1. ! The critical cluster contains one molecules but the produced cluster contains 2 molecules … … 604 583 na_n=x_n*ntot_n 605 584 IF (na_n .LT. 1.) THEN 606 ! print *, 'Warning: number of acid molecules < 1 in nucleation regime, setting na_n=1'607 585 na_n=1.0 608 586 ENDIF … … 664 642 665 643 IF (kinetic_i) THEN 666 ! jnuc_i1=1.0E6*(0.3E-9 + 0.487E-9)**2.*SQRT(8.*3.141593*1.38E-23*(1./(1.661E-27*98.07)+1./(1.661E-27*98.07)))* &667 644 jnuc_i1=1.0E6*(0.3E-9 + 0.487E-9)**2.*SQRT(8.*RPI*RKBOL*(1./mH2SO4mol+1./mH2SO4mol))* & 668 645 & SQRT(tli)*rhoali !1/cm3s … … 816 793 na_i=x_i*ntot_i 817 794 IF (na_i .LT. 1.) THEN 818 ! print *, 'Warning: number of acid molecules < 1 in nucleation regime, setting na_n=1'819 795 na_n=1.0 820 796 ENDIF -
LMDZ6/trunk/libf/phylmd/StratAer/ocs_to_so2.F90
r3677 r4601 8 8 USE infotrac_phy 9 9 USE YOMCST, ONLY : RG 10 USE phys_local_var_mod, ONLY : OCS_lifetime, budg_3D_ocs_to_so2, budg_ocs_to_so2 10 USE phys_local_var_mod, ONLY : OCS_lifetime, budg_3D_ocs_to_so2, budg_ocs_to_so2 11 USE strataer_local_var_mod, ONLY : flag_min_rreduce 11 12 12 13 IMPLICIT NONE … … 23 24 ! local variables 24 25 INTEGER :: i,j,k,nb,ilon,ilev 25 26 REAL :: rreduce 27 26 28 !--convert OCS to SO2 27 29 budg_3D_ocs_to_so2(:,:)=0.0 … … 29 31 30 32 DO ilon=1, klon 31 DO ilev=1, klev 32 !only in the stratosphere 33 IF (is_strato(ilon,ilev)) THEN 34 IF (OCS_lifetime(ilon,ilev).GT.0.0) THEN 35 budg_3D_ocs_to_so2(ilon,ilev)=tr_seri(ilon,ilev,id_OCS_strat)*(1.0-exp(-pdtphys/OCS_lifetime(ilon,ilev))) 36 ENDIF 37 tr_seri(ilon,ilev,id_OCS_strat)=tr_seri(ilon,ilev,id_OCS_strat) - budg_3D_ocs_to_so2(ilon,ilev) 38 tr_seri(ilon,ilev,id_SO2_strat)=tr_seri(ilon,ilev,id_SO2_strat) + mSO2mol/mOCSmol*budg_3D_ocs_to_so2(ilon,ilev) 39 !convert budget from kg(OCS)/kgA to kg(S)/m2/layer/s for saving as diagnostic 40 budg_3D_ocs_to_so2(ilon,ilev)=budg_3D_ocs_to_so2(ilon,ilev)*mSatom/mOCSmol*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG/pdtphys 41 budg_ocs_to_so2(ilon)=budg_ocs_to_so2(ilon)+budg_3D_ocs_to_so2(ilon,ilev) 42 ENDIF 43 ENDDO 33 DO ilev=1, klev 34 !only in the stratosphere 35 IF (is_strato(ilon,ilev)) THEN 36 IF (OCS_lifetime(ilon,ilev).GT.0.0.AND.OCS_lifetime(ilon,ilev).LT.1.E10) THEN 37 rreduce = OCS_lifetime(ilon,ilev) 38 ! Check lifetime rreduce < timestep*3 (such as H2SO4 loss > 0.28*H2SO4) with exp(-1/3)=0.72 39 IF(flag_min_rreduce) THEN 40 IF (rreduce .LT. (3.*pdtphys)) rreduce = 3.*pdtphys 41 ENDIF 42 budg_3D_ocs_to_so2(ilon,ilev)=tr_seri(ilon,ilev,id_OCS_strat)*(1.0-exp(-pdtphys/rreduce)) 43 tr_seri(ilon,ilev,id_OCS_strat)=tr_seri(ilon,ilev,id_OCS_strat) - budg_3D_ocs_to_so2(ilon,ilev) 44 tr_seri(ilon,ilev,id_SO2_strat)=tr_seri(ilon,ilev,id_SO2_strat) + mSO2mol/mOCSmol*budg_3D_ocs_to_so2(ilon,ilev) 45 !convert budget from kg(OCS)/kgA to kg(S)/m2/layer/s for saving as diagnostic 46 budg_3D_ocs_to_so2(ilon,ilev)=budg_3D_ocs_to_so2(ilon,ilev)*mSatom/mOCSmol*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG/pdtphys 47 budg_ocs_to_so2(ilon)=budg_ocs_to_so2(ilon)+budg_3D_ocs_to_so2(ilon,ilev) 48 ENDIF 49 ! END IF(OCS_lifetime(ilon,ilev).GT.0.0.AND.OCS_lifetime(ilon,ilev).LT.1.E10) 50 ENDIF 51 ! END IF(is_strato(ilon,ilev)) 52 ENDDO 44 53 ENDDO 45 54 -
LMDZ6/trunk/libf/phylmd/StratAer/so2_to_h2so4.F90
r3677 r4601 7 7 USE aerophys 8 8 USE infotrac_phy 9 USE YOMCST, ONLY : RG 10 USE phys_local_var_mod, ONLY : SO2_lifetime, budg_3D_so2_to_h2so4, budg_so2_to_h2so4 11 9 USE YOMCST, ONLY : RG, RD 10 ! lifetime (sec) et O3_clim (VMR) 11 USE phys_local_var_mod, ONLY : SO2_lifetime, H2SO4_lifetime, O3_clim, budg_3D_so2_to_h2so4, budg_so2_to_h2so4 12 USE strataer_local_var_mod, ONLY : flag_OH_reduced, flag_H2SO4_photolysis, flag_min_rreduce 13 12 14 IMPLICIT NONE 13 15 … … 24 26 ! local variables in coagulation routine 25 27 INTEGER :: i,j,k,nb,ilon,ilev 26 28 REAL :: dummyso4toso2,rkho2o3,rkoho3,rkso2oh 29 REAL :: rmairdens,rrak0,rrak1,rrate,rreduce 30 27 31 !--convert SO2 to H2SO4 28 32 budg_3D_so2_to_h2so4(:,:)=0.0 … … 30 34 31 35 DO ilon=1, klon 32 DO ilev=1, klev 33 !only in the stratosphere 34 IF (is_strato(ilon,ilev)) THEN 35 IF (SO2_lifetime(ilon,ilev).GT.0.0) THEN 36 budg_3D_so2_to_h2so4(ilon,ilev)=tr_seri(ilon,ilev,id_SO2_strat)*(1.0-exp(-pdtphys/SO2_lifetime(ilon,ilev))) 37 ENDIF 38 tr_seri(ilon,ilev,id_SO2_strat)=tr_seri(ilon,ilev,id_SO2_strat) - budg_3D_so2_to_h2so4(ilon,ilev) 39 tr_seri(ilon,ilev,id_H2SO4_strat)=tr_seri(ilon,ilev,id_H2SO4_strat) + mH2SO4mol/mSO2mol*budg_3D_so2_to_h2so4(ilon,ilev) 40 !convert budget from kg(SO2)/kgA to kg(S)/m2/layer/s for saving as diagnostic 41 budg_3D_so2_to_h2so4(ilon,ilev)=budg_3D_so2_to_h2so4(ilon,ilev)*mSatom/mSO2mol*(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG/pdtphys 42 budg_so2_to_h2so4(ilon)=budg_so2_to_h2so4(ilon)+budg_3D_so2_to_h2so4(ilon,ilev) 43 ENDIF 44 ENDDO 45 ENDDO 46 36 DO ilev=1, klev 37 !only in the stratosphere 38 IF (is_strato(ilon,ilev)) THEN 39 IF (SO2_lifetime(ilon,ilev).GT.0.0 .AND. SO2_lifetime(ilon,ilev).LT.1.E10) THEN 40 IF (flag_OH_reduced) THEN 41 !--convert SO2 to H2SO4 (slimane) 42 ! Reduce OH because of SO2 43 ! d[OH]/dt = k1[HO2][O] +k2[HO2][O3] +k3[HO2][NO] 44 ! − (k4[O] +k5[O3] +k6[CO] +ks[SO2])[OH] 45 ! In steady-state: d[OH]/dt = 0, so the ratio [OH]/[HO2] = 46 ! (k1[O] +k2[O3] +k3[NO]) / (k4[O] +k5[O3] +k6[CO] +ks[SO2]) 47 ! In the lower and middle stratosphere, 48 ! [OH]/[HO2] ~(k2[O3] +k3[NO]) / (k5[O3] +k6[CO] +ks[SO2]) 49 ! CO is only important near the tropopause. NO should not dominate over O3. 50 ! So when consideringthe big terms: [OH]/[HO2] ~ ~ k2[O3] / (k5[O3] +ks[SO2]) 51 ! For c: control run c (no effect of SO2), [OHc]/[HO2c] ~ ~ k2[O3] / k5[O3] = 1/ Rc 52 ! For p: perturbed run (high SO2), [OHp]/[HO2p] ~ ~ k2[O3] / (k5[O3] +ks[SO2]) =1/Rp 53 ! If HOx = OH + HO2 = constante, 54 ! OHc + HO2c = OHp + HO2p 55 ! OHc (1+Rc) = OHp (1+Rp) 56 ! OHp = OHc (1+Rc) / (1+Rp) 57 ! 1+Rc = (k2[O3] +k5[O3] ) / k2[O3] 58 ! 1+Rp = (k2[O3] +k5[O3] +ks[SO2]) / k2[O3] 59 ! (1+Rc) / (1+Rp)= (k2[O3] +k5[O3] )/ (k2[O3] +k5[O3] +ks[SO2]) 60 ! OHp = OHc * (k(HO2+OH)*[O3] +k(OH+O3)[O3] ) / 61 ! (k(HO2+OH)*[O3] +k(OH+O3)[O3] +k(SO2+OH)[SO2]) 62 ! Final: OHp = OHc * (k(HO2+O3)*[O3] +k(OH+O3)[O3] ) / 63 ! (k(HO2+O3)*[O3] +k(OH+O3)[O3] +k(SO2+OH)[SO2]) 64 65 ! latest jpl 2020 66 rkho2o3 = 1.0e-14*exp( -490./t_seri(ilon,ilev) ) 67 rkoho3 = 1.7e-12*exp( -940./t_seri(ilon,ilev) ) 68 69 ! rkso2oh= so2 + oh + m ->hso3 + m-> h2so4 + m 70 ! air density M(molecule/cm3) = Pa/(K*1.38e-19) 71 rmairdens = pplay(ilon,ilev)/(t_seri(ilon,ilev)*1.38e-19) 72 ! JPL 2020 73 rrak0 = 3.3e-31*( (300./t_seri(ilon,ilev))**4.3 ) 74 rrak1 = 1.6e-12 75 rrate = (rrak0*rmairdens) / ( 1. + rrak0*rmairdens/rrak1 ) 76 rreduce = 1./( 1. + log10((rrak0*rmairdens)/rrak1)**2 ) 77 rkso2oh = rrate*(0.6**rreduce) 78 79 ! O3 (molec/cm3): O3_clim in vmr 80 rrak0 = O3_clim(ilon,ilev)*rmairdens 81 ! SO2 (molec/cm3): convert from kg/kgA 82 rrak1 = tr_seri(ilon,ilev,id_SO2_strat) & 83 & *pplay(ilon,ilev)/t_seri(ilon,ilev)/RD/1.E6/mSO2mol 84 85 IF (rrak1 .GE. 0.0) THEN 86 rreduce =( (rkho2o3+rkoho3)*rrak0 + rkso2oh*rrak1 ) / & 87 ( (rkho2o3+rkoho3)*rrak0 ) 88 ! reduce OH, so extend SO2 lifetime 89 rreduce = rreduce*SO2_lifetime(ilon,ilev) 90 ELSE 91 rreduce = SO2_lifetime(ilon,ilev) 92 ENDIF 93 ! end of IF (rrak1) THEN 94 ELSE 95 rreduce = SO2_lifetime(ilon,ilev) 96 ENDIF 97 ! end of IF (flag_OH_reduced) THEN 98 99 ! Check lifetime rreduce < timestep*1.5 (such as SO2 loss > 0.5*SO2) with exp(-1/1.5)=0.52 100 ! Check lifetime rreduce < timestep*3 (such as SO2 loss > 0.28*SO2) with exp(-1/3)=0.72 101 IF(flag_min_rreduce) THEN 102 IF (rreduce .LT. (3.*pdtphys)) rreduce = 3.*pdtphys 103 ENDIF 104 budg_3D_so2_to_h2so4(ilon,ilev)=tr_seri(ilon,ilev,id_SO2_strat)*(1.0-exp(-pdtphys/rreduce)) 105 tr_seri(ilon,ilev,id_SO2_strat)=tr_seri(ilon,ilev,id_SO2_strat) - budg_3D_so2_to_h2so4(ilon,ilev) 106 tr_seri(ilon,ilev,id_H2SO4_strat)=tr_seri(ilon,ilev,id_H2SO4_strat) + mH2SO4mol/mSO2mol*budg_3D_so2_to_h2so4(ilon,ilev) 107 ENDIF 108 ! IF (SO2_lifetime(ilon,ilev).GT.0.0 .AND. SO2_lifetime(ilon,ilev).LT.1.E10) THEN 109 110 111 IF (flag_H2SO4_photolysis) THEN 112 !--convert H2SO4 to SO2 using HCl photolysis (Rinsland et al, 1995) 113 ! 4.6e-8 sec-1: J strato taken from Feierabend, et al., Chem. Phys. Lett., 2006. 114 ! error before: ...exp(-pdtphys/4.6e-8) so negligible photolysis, 115 ! should have been exp(-pdtphys*4.6e-8) 116 ! Also: Lane and Kjaergaard: ,J. Phys. Chem. A, 2008 117 ! We find that the dominant photodissociation mechanism of H2SO4 118 ! below 70 km is absorption in the visible region by OH stretching overtone 119 ! transitions, whereas above 70 km, absorption of Lyman-R ... 120 ! In 2D model, 0.3 Rinsland et al, 1995 but Burkholder, Mills and McKeen say 0.016 121 ! JH2SO4=JHCL*0.3 (*0.016 would be too slow) 122 ! test: quick sequential integration for SO2 to H2SO4 and reverse 123 124 IF (H2SO4_lifetime(ilon,ilev).GT.0.0 .AND. H2SO4_lifetime(ilon,ilev).LT.1.E10) THEN 125 126 rreduce = H2SO4_lifetime(ilon,ilev) 127 dummyso4toso2 = 0.0 128 129 ! Check lifetime rreduce < timestep*1.5 (such as H2SO4 loss > 0.5*H2SO4) with exp(-1/1.5)=0.52 130 ! Check lifetime rreduce < timestep*3 (such as H2SO4 loss > 0.28*H2SO4) with exp(-1/3)=0.72 131 IF(flag_min_rreduce) THEN 132 IF (rreduce .LT. (3.*pdtphys)) rreduce = 3.*pdtphys 133 ENDIF 134 dummyso4toso2 = (mSO2mol/mH2SO4mol)*tr_seri(ilon,ilev,id_H2SO4_strat)*(1.0-exp(-pdtphys/rreduce)) 135 budg_3D_so2_to_h2so4(ilon,ilev) = budg_3D_so2_to_h2so4(ilon,ilev) + dummyso4toso2 136 tr_seri(ilon,ilev,id_SO2_strat)=tr_seri(ilon,ilev,id_SO2_strat) + dummyso4toso2 137 tr_seri(ilon,ilev,id_H2SO4_strat)=tr_seri(ilon,ilev,id_H2SO4_strat) - (mH2SO4mol/mSO2mol)*dummyso4toso2 138 ENDIF 139 ! IF (H2SO4_lifetime(ilon,ilev).GT.0.0 .AND. H2SO4_lifetime(ilon,ilev).LT.1.E10) THEN 140 ENDIF 141 ! end of IF (flag_H2SO4_photolysis) THEN 142 143 !convert budget from kg(SO2)/kgA to kg(S)/m2/layer/s for saving as diagnostic 144 budg_3D_so2_to_h2so4(ilon,ilev)=budg_3D_so2_to_h2so4(ilon,ilev)*mSatom/mSO2mol* & 145 (paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG/pdtphys 146 budg_so2_to_h2so4(ilon)=budg_so2_to_h2so4(ilon)+budg_3D_so2_to_h2so4(ilon,ilev) 147 ENDIF 148 ! IF (is_strato(ilon,ilev)) THEN 149 ENDDO ! DO (klev) 150 ENDDO ! DO (klon) 151 47 152 END SUBROUTINE SO2_TO_H2SO4 -
LMDZ6/trunk/libf/phylmd/StratAer/traccoag_mod.F90
r4513 r4601 23 23 USE mod_phys_lmdz_mpi_data, ONLY : is_mpi_root 24 24 USE mod_phys_lmdz_para, only: gather, scatter 25 USE phys_cal_mod, ONLY : year_len, mth_len,year_cur, mth_cur, day_cur, hour25 USE phys_cal_mod, ONLY : year_len, year_cur, mth_cur, day_cur, hour 26 26 USE sulfate_aer_mod 27 27 USE phys_local_var_mod, ONLY: stratomask 28 28 USE YOMCST 29 29 USE print_control_mod, ONLY: lunout 30 USE strataer_ mod31 30 USE strataer_local_var_mod 31 32 32 IMPLICIT NONE 33 33 … … 58 58 !---------------- 59 59 REAL :: m_aer_emiss_vol_daily ! daily injection mass emission 60 REAL :: m_aer ! aerosol mass 60 61 INTEGER :: it, k, i, ilon, ilev, itime, i_int, ieru 61 62 LOGICAL,DIMENSION(klon,klev) :: is_strato ! true = above tropopause, false = below … … 78 79 REAL :: theta_min, theta_max ! for SAI computation between two latitudes 79 80 REAL :: dlat_loc 81 REAL :: latmin,latmax,lonmin,lonmax ! lat/lon min/max for injection 82 REAL :: sigma_alt, altemiss ! injection altitude + sigma for distrib 83 REAL :: pdt,stretchlong ! physic timestep, stretch emission over one day 84 80 85 INTEGER :: injdur_sai ! injection duration for SAI case [days] 81 86 INTEGER :: yr, is_bissext 82 87 83 IF (is_mpi_root ) THEN88 IF (is_mpi_root .AND. flag_verbose_strataer) THEN 84 89 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_emit90 WRITE(lunout,*) 'IN traccoag flag_emit: ',flag_emit 86 91 ENDIF 87 92 … … 126 131 !--calculate mass of air in every grid box 127 132 DO ilon=1, klon 128 DO ilev=1, klev 129 m_air_gridbox(ilon,ilev)=(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG*cell_area(ilon) 133 DO ilev=1, klev 134 m_air_gridbox(ilon,ilev)=(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG*cell_area(ilon) 135 ENDDO 130 136 ENDDO 131 ENDDO 132 133 ! IF (debutphy) THEN 134 ! CALL gather(tr_seri, tr_seri_glo) 135 ! IF (MAXVAL(tr_seri_glo).LT.1.e-30) THEN 136 !--initialising tracer concentrations to zero 137 ! DO it=1, nbtr 138 ! tr_seri(:,:,it)=0.0 139 ! ENDDO 140 ! ENDIF 141 ! ENDIF 142 137 143 138 !--initialise emission diagnostics 139 budg_emi(:,1)=0.0 144 140 budg_emi_ocs(:)=0.0 145 141 budg_emi_so2(:)=0.0 … … 147 143 budg_emi_part(:)=0.0 148 144 149 !--sulfur emission, depending on chosen scenario (flag_ sulf_emit)150 SELECT CASE(flag_ sulf_emit)145 !--sulfur emission, depending on chosen scenario (flag_emit) 146 SELECT CASE(flag_emit) 151 147 152 148 CASE(0) ! background aerosol … … 159 155 IF (year_cur==year_emit_vol(ieru).AND.mth_cur==mth_emit_vol(ieru).AND.& 160 156 day_cur>=day_emit_vol(ieru).AND.day_cur<(day_emit_vol(ieru)+injdur)) THEN 161 ! 162 ! daily injection mass emission - NL 163 m_aer_emiss_vol_daily = m_aer_emiss_vol(ieru)/(REAL(injdur)*REAL(ponde_lonlat_vol(ieru))) 164 WRITE(lunout,*) 'IN traccoag DD m_aer_emiss_vol(ieru)=',m_aer_emiss_vol(ieru), & 157 158 ! daily injection mass emission 159 m_aer=m_aer_emiss_vol(ieru,1)/(REAL(injdur)*REAL(ponde_lonlat_vol(ieru))) 160 !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss) 161 m_aer=m_aer*(mSO2mol/mSatom) 162 163 WRITE(lunout,*) 'IN traccoag m_aer_emiss_vol(ieru)=',m_aer_emiss_vol(ieru,1), & 165 164 'ponde_lonlat_vol(ieru)=',ponde_lonlat_vol(ieru),'(injdur*ponde_lonlat_vol(ieru))', & 166 (injdur*ponde_lonlat_vol(ieru)),'m_aer_emiss_vol_daily=',m_aer _emiss_vol_daily,'ieru=',ieru165 (injdur*ponde_lonlat_vol(ieru)),'m_aer_emiss_vol_daily=',m_aer,'ieru=',ieru 167 166 WRITE(lunout,*) 'IN traccoag, dlon=',dlon 168 DO i=1,klon 169 !Pinatubo eruption at 15.14N, 120.35E 170 dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes 171 WRITE(lunout,*) 'IN traccoag, dlat=',dlat_loc 172 IF ( xlat(i).GE.xlat_min_vol(ieru)-dlat_loc .AND. xlat(i).LT.xlat_max_vol(ieru)+dlat_loc .AND. & 173 xlon(i).GE.xlon_min_vol(ieru)-dlon .AND. xlon(i).LT.xlon_max_vol(ieru)+dlon ) THEN 174 ! 175 WRITE(lunout,*) 'coordinates of volcanic injection point=',xlat(i),xlon(i),day_cur,mth_cur,year_cur 176 WRITE(lunout,*) 'DD m_aer_emiss_vol_daily=',m_aer_emiss_vol_daily 177 ! compute altLMDz 178 altLMDz(:)=0.0 179 DO k=1, klev 180 zrho=pplay(i,k)/t_seri(i,k)/RD !air density in kg/m3 181 zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG !mass of layer in kg 182 zdz=zdm(k)/zrho !thickness of layer in m 183 altLMDz(k+1)=altLMDz(k)+zdz !altitude of interface 184 ENDDO 185 186 SELECT CASE(flag_sulf_emit_distrib) 187 188 CASE(0) ! Gaussian distribution 189 !compute distribution of emission to vertical model layers (based on Gaussian peak in altitude) 190 f_lay_sum=0.0 191 DO k=1, klev 192 f_lay_emiss(k)=0.0 193 DO i_int=1, n_int_alt 194 alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt) 195 f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_vol(ieru))* & 196 & exp(-0.5*((alt-altemiss_vol(ieru))/sigma_alt_vol(ieru))**2.)* & 197 & (altLMDz(k+1)-altLMDz(k))/float(n_int_alt) 198 ENDDO 199 f_lay_sum=f_lay_sum+f_lay_emiss(k) 200 ENDDO 201 202 CASE(1) ! Uniform distribution 203 ! In this case, parameter sigma_alt_vol(ieru) is considered to be half the 204 ! height of the injection, centered around altemiss_vol(ieru) 205 DO k=1, klev 206 f_lay_emiss(k)=max(min(altemiss_vol(ieru)+sigma_alt_vol(ieru),altLMDz(k+1))- & 207 & max(altemiss_vol(ieru)-sigma_alt_vol(ieru),altLMDz(k)),0.)/(2.*sigma_alt_vol(ieru)) 208 f_lay_sum=f_lay_sum+f_lay_emiss(k) 209 ENDDO 210 211 END SELECT ! End CASE over flag_sulf_emit_distrib) 212 213 WRITE(lunout,*) "IN traccoag m_aer_emiss_vol=",m_aer_emiss_vol(ieru) 214 WRITE(lunout,*) "IN traccoag f_lay_emiss=",f_lay_emiss 215 !correct for step integration error 216 f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum 217 !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss) 218 !vertically distributed emission 219 DO k=1, klev 220 ! stretch emission over one day of Pinatubo eruption 221 emission=m_aer_emiss_vol_daily*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/1./(86400.-pdtphys) 222 tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys 223 budg_emi_so2(i)=budg_emi_so2(i)+emission*zdm(k)*mSatom/mSO2mol 224 ENDDO 225 ENDIF ! emission grid cell 226 ENDDO ! klon loop 227 WRITE(lunout,*) "IN traccoag (ieru=",ieru,") m_aer_emiss_vol_daily=",m_aer_emiss_vol_daily 167 168 latmin=xlat_min_vol(ieru) 169 latmax=xlat_max_vol(ieru) 170 lonmin=xlon_min_vol(ieru) 171 lonmax=xlon_max_vol(ieru) 172 altemiss = altemiss_vol(ieru) 173 sigma_alt = sigma_alt_vol(ieru) 174 pdt=pdtphys 175 ! stretch emission over one day of eruption 176 stretchlong = 1. 177 178 CALL STRATEMIT(pdtphys,pdt,xlat,xlon,t_seri,pplay,paprs,tr_seri,& 179 m_aer,latmin,latmax,lonmin,lonmax,altemiss,sigma_alt,id_SO2_strat, & 180 stretchlong,1,0) 181 228 182 ENDIF ! emission period 229 183 ENDDO ! eruption number … … 305 259 .AND. ( day_emit_sai_start <= day_cur .OR. mth_emit_sai_start < mth_cur .OR. year_emit_sai_start < year_cur ) & 306 260 .AND. ( day_cur <= day_emit_sai_end .OR. mth_cur < mth_emit_sai_end .OR. year_cur < year_emit_sai_end )) THEN 307 308 DO i=1,klon 309 dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes 310 IF ( xlat(i).GE.xlat_sai-dlat_loc .AND. xlat(i).LT.xlat_sai+dlat_loc .AND. & 311 & xlon(i).GE.xlon_sai-dlon .AND. xlon(i).LT.xlon_sai+dlon ) THEN 312 ! 313 ! compute altLMDz 314 altLMDz(:)=0.0 315 DO k=1, klev 316 zrho=pplay(i,k)/t_seri(i,k)/RD !air density in kg/m3 317 zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG !mass of layer in kg 318 zdz=zdm(k)/zrho !thickness of layer in m 319 altLMDz(k+1)=altLMDz(k)+zdz !altitude of interface 320 ENDDO 321 322 SELECT CASE(flag_sulf_emit_distrib) 323 324 CASE(0) ! Gaussian distribution 325 !compute distribution of emission to vertical model layers (based on Gaussian peak in altitude) 326 f_lay_sum=0.0 327 DO k=1, klev 328 f_lay_emiss(k)=0.0 329 DO i_int=1, n_int_alt 330 alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt) 331 f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_sai)* & 332 & exp(-0.5*((alt-altemiss_sai)/sigma_alt_sai)**2.)* & 333 & (altLMDz(k+1)-altLMDz(k))/float(n_int_alt) 334 ENDDO 335 f_lay_sum=f_lay_sum+f_lay_emiss(k) 336 ENDDO 337 338 CASE(1) ! Uniform distribution 339 f_lay_sum=0.0 340 ! In this case, parameter sigma_alt_vol(ieru) is considered to be half 341 ! the height of the injection, centered around altemiss_sai 342 DO k=1, klev 343 f_lay_emiss(k)=max(min(altemiss_sai+sigma_alt_sai,altLMDz(k+1))- & 344 & max(altemiss_sai-sigma_alt_sai,altLMDz(k)),0.)/(2.*sigma_alt_sai) 345 f_lay_sum=f_lay_sum+f_lay_emiss(k) 346 ENDDO 347 348 END SELECT ! Gaussian or uniform distribution 349 350 !correct for step integration error 351 f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum 352 !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss) 353 !vertically distributed emission 354 DO k=1, klev 355 ! stretch emission over whole year (360d) 356 emission=m_aer_emiss_sai*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/FLOAT(injdur_sai)/86400. 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)/FLOAT(year_len)/86400. 366 ! tr_seri(i,k,id_BIN01_strat+20)=tr_seri(i,k,id_BIN01_strat+20)+emission*pdtphys 367 ! budg_emi_part(i)=budg_emi_part(i)+emission*zdm(k)*mSatom/mH2SO4mol 368 ! ENDDO 369 ENDIF ! emission grid cell 370 ENDDO ! klon loop 371 261 262 m_aer=m_aer_emiss_sai 263 !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss) 264 m_aer=m_aer*(mSO2mol/mSatom) 265 266 latmin=xlat_sai 267 latmax=xlat_sai 268 lonmin=xlon_sai 269 lonmax=xlon_sai 270 altemiss = altemiss_sai 271 sigma_alt = sigma_alt_sai 272 pdt=0. 273 ! stretch emission over whole year (360d) 274 stretchlong=FLOAT(year_len) 275 276 CALL STRATEMIT(pdtphys,pdt,xlat,xlon,t_seri,pplay,paprs,m_air_gridbox,tr_seri,& 277 m_aer,latmin,latmax,lonmin,lonmax,altemiss,sigma_alt,id_SO2_strat, & 278 stretchlong,1,0) 279 280 budg_emi_so2(:) = budg_emi(:,1)*mSatom/mSO2mol 372 281 ENDIF ! Condition over injection dates 373 282 … … 375 284 ! lat_min and lat_max 376 285 377 DO i=1,klon 378 ! SAI scenario with continuous emission 379 dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes 380 theta_min = max(xlat(i)-dlat_loc,xlat_min_sai) 381 theta_max = min(xlat(i)+dlat_loc,xlat_max_sai) 382 IF ( xlat(i).GE.xlat_min_sai-dlat_loc .AND. xlat(i).LT.xlat_max_sai+dlat_loc .AND. & 383 & xlon(i).GE.xlon_sai-dlon .AND. xlon(i).LT.xlon_sai+dlon ) THEN 384 ! 385 ! compute altLMDz 386 altLMDz(:)=0.0 387 DO k=1, klev 388 zrho=pplay(i,k)/t_seri(i,k)/RD !air density in kg/m3 389 zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG !mass of layer in kg 390 zdz=zdm(k)/zrho !thickness of layer in m 391 altLMDz(k+1)=altLMDz(k)+zdz !altitude of interface 392 ENDDO 393 394 SELECT CASE(flag_sulf_emit_distrib) 395 396 CASE(0) ! Gaussian distribution 397 !compute distribution of emission to vertical model layers (based on 398 !Gaussian peak in altitude) 399 f_lay_sum=0.0 400 DO k=1, klev 401 f_lay_emiss(k)=0.0 402 DO i_int=1, n_int_alt 403 alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt) 404 f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_sai)* & 405 & exp(-0.5*((alt-altemiss_sai)/sigma_alt_sai)**2.)* & 406 & (altLMDz(k+1)-altLMDz(k))/float(n_int_alt) 407 ENDDO 408 f_lay_sum=f_lay_sum+f_lay_emiss(k) 409 ENDDO 410 411 CASE(1) ! Uniform distribution 412 f_lay_sum=0.0 413 ! In this case, parameter sigma_alt_vol(ieru) is considered to be half 414 ! the height of the injection, centered around altemiss_sai 415 DO k=1, klev 416 f_lay_emiss(k)=max(min(altemiss_sai+sigma_alt_sai,altLMDz(k+1))- & 417 & max(altemiss_sai-sigma_alt_sai,altLMDz(k)),0.)/(2.*sigma_alt_sai) 418 f_lay_sum=f_lay_sum+f_lay_emiss(k) 419 ENDDO 420 421 END SELECT ! Gaussian or uniform distribution 422 423 !correct for step integration error 424 f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum 425 !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss) 426 !vertically distributed emission 427 DO k=1, klev 428 ! stretch emission over whole year (360d) 429 emission=m_aer_emiss_sai*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/ & 430 & FLOAT(year_len)/86400.*(sin(theta_max/180.*RPI)-sin(theta_min/180.*RPI))/ & 431 & (sin(xlat_max_sai/180.*RPI)-sin(xlat_min_sai/180.*RPI)) 432 tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys 433 budg_emi_so2(i)=budg_emi_so2(i)+emission*zdm(k)*mSatom/mSO2mol 434 ENDDO 435 436 ! !emission as monodisperse particles with 0.1um dry radius (BIN21) 437 ! !vertically distributed emission 438 ! DO k=1, klev 439 ! ! stretch emission over whole year (360d) 440 ! emission=m_aer_emiss*(mH2SO4mol/mSatom)/m_part_dry(21)/m_air_gridbox(i,k)*f_lay_emiss(k)/year_len/86400 441 ! tr_seri(i,k,id_BIN01_strat+20)=tr_seri(i,k,id_BIN01_strat+20)+emission*pdtphys 442 ! budg_emi_part(i)=budg_emi_part(i)+emission*zdm(k)*mSatom/mH2SO4mol 443 ! ENDDO 444 ENDIF ! emission grid cell 445 ENDDO ! klon loop 446 447 END SELECT ! emission scenario (flag_sulf_emit) 286 m_aer=m_aer_emiss_sai 287 !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss) 288 m_aer=m_aer*(mSO2mol/mSatom) 289 290 latmin=xlat_min_sai 291 latmax=xlat_max_sai 292 lonmin=xlon_sai 293 lonmax=xlon_sai 294 altemiss = altemiss_sai 295 sigma_alt = sigma_alt_sai 296 pdt=0. 297 ! stretch emission over whole year (360d) 298 stretchlong=FLOAT(year_len) 299 300 CALL STRATEMIT(pdtphys,pdt,xlat,xlon,t_seri,pplay,paprs,m_air_gridbox,tr_seri,& 301 m_aer,latmin,latmax,lonmin,lonmax,altemiss,sigma_alt,id_SO2_strat, & 302 stretchlong,1,0) 303 304 budg_emi_so2(:) = budg_emi(:,1)*mSatom/mSO2mol 305 306 END SELECT ! emission scenario (flag_emit) 448 307 449 308 !--read background concentrations of OCS and SO2 and lifetimes from input file … … 463 322 CALL coagulate(pdtphys,mdw,tr_seri,t_seri,pplay,dens_aer,is_strato) 464 323 465 !--call sedimentation routine 324 !--call sedimentation routine 466 325 CALL aer_sedimnt(pdtphys, t_seri, pplay, paprs, tr_seri, dens_aer) 467 326 … … 480 339 ENDDO 481 340 ENDDO 482 483 ! CALL minmaxsimple(tr_seri(:,:,id_SO2_strat),0.0,0.0,'fin SO2') 484 ! DO it=1, nbtr_bin 485 ! CALL minmaxsimple(tr_seri(:,:,nbtr_sulgas+it),0.0,0.0,'fin bin ') 486 ! ENDDO 487 341 488 342 END SUBROUTINE traccoag 489 343
Note: See TracChangeset
for help on using the changeset viewer.
