Changeset 4727 for LMDZ6/branches/LMDZ_ECRad/libf/phylmd/StratAer
- Timestamp:
- Oct 19, 2023, 4:02:57 PM (15 months ago)
- Location:
- LMDZ6/branches/LMDZ_ECRad
- Files:
-
- 1 deleted
- 8 edited
- 6 copied
-
. (modified) (1 prop)
-
libf/phylmd/StratAer/aerophys.F90 (modified) (1 diff)
-
libf/phylmd/StratAer/interp_sulf_input.F90 (modified) (15 diffs)
-
libf/phylmd/StratAer/micphy_tstep.F90 (modified) (1 diff)
-
libf/phylmd/StratAer/nucleation_tstep_mod.F90 (modified) (10 diffs)
-
libf/phylmd/StratAer/ocs_to_so2.F90 (modified) (3 diffs)
-
libf/phylmd/StratAer/so2_to_h2so4.F90 (modified) (3 diffs)
-
libf/phylmd/StratAer/stratH2O_methox.F90 (copied) (copied from LMDZ6/trunk/libf/phylmd/StratAer/stratH2O_methox.F90)
-
libf/phylmd/StratAer/strataer_emiss_mod.F90 (copied) (copied from LMDZ6/trunk/libf/phylmd/StratAer/strataer_emiss_mod.F90)
-
libf/phylmd/StratAer/strataer_local_var_mod.F90 (copied) (copied from LMDZ6/trunk/libf/phylmd/StratAer/strataer_local_var_mod.F90)
-
libf/phylmd/StratAer/strataer_mod.F90 (deleted)
-
libf/phylmd/StratAer/strataer_nuc_mod.F90 (copied) (copied from LMDZ6/trunk/libf/phylmd/StratAer/strataer_nuc_mod.F90)
-
libf/phylmd/StratAer/stratdistrib.F90 (copied) (copied from LMDZ6/trunk/libf/phylmd/StratAer/stratdistrib.F90)
-
libf/phylmd/StratAer/stratemit.F90 (copied) (copied from LMDZ6/trunk/libf/phylmd/StratAer/stratemit.F90)
-
libf/phylmd/StratAer/traccoag_mod.F90 (modified) (9 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/LMDZ_ECRad
- Property svn:mergeinfo changed
-
LMDZ6/branches/LMDZ_ECRad/libf/phylmd/StratAer/aerophys.F90
r3526 r4727 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/branches/LMDZ_ECRad/libf/phylmd/StratAer/interp_sulf_input.F90
r3677 r4727 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 … … 42 43 REAL paprs_glo(klon_glo,klev+1) 43 44 44 REAL, POINTER:: latitude(:)45 REAL, ALLOCATABLE :: latitude(:) 45 46 ! (of input data sorted in strictly ascending order) 46 47 47 REAL, POINTER:: longitude(:)48 REAL, ALLOCATABLE :: longitude(:) 48 49 ! (of input data sorted in strictly ascending order) 49 50 50 REAL, POINTER:: time(:)51 REAL, ALLOCATABLE :: time(:) 51 52 ! (of input data sorted in strictly ascending order) 52 53 53 REAL, POINTER:: lev(:)54 REAL, ALLOCATABLE :: lev(:) 54 55 ! levels of input data 55 56 … … 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/branches/LMDZ_ECRad/libf/phylmd/StratAer/micphy_tstep.F90
r4482 r4727 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/branches/LMDZ_ECRad/libf/phylmd/StratAer/nucleation_tstep_mod.F90
r4482 r4727 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/branches/LMDZ_ECRad/libf/phylmd/StratAer/ocs_to_so2.F90
r3677 r4727 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/branches/LMDZ_ECRad/libf/phylmd/StratAer/so2_to_h2so4.F90
r3677 r4727 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/branches/LMDZ_ECRad/libf/phylmd/StratAer/traccoag_mod.F90
r4482 r4727 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, 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 80 81 IF (is_mpi_root) THEN 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 85 INTEGER :: injdur_sai ! injection duration for SAI case [days] 86 INTEGER :: yr, is_bissext 87 88 IF (is_mpi_root .AND. flag_verbose_strataer) THEN 82 89 WRITE(lunout,*) 'in traccoag: date from phys_cal_mod =',year_cur,'-',mth_cur,'-',day_cur,'-',hour 83 WRITE(lunout,*) 'IN traccoag flag_ sulf_emit: ',flag_sulf_emit90 WRITE(lunout,*) 'IN traccoag flag_emit: ',flag_emit 84 91 ENDIF 85 92 … … 124 131 !--calculate mass of air in every grid box 125 132 DO ilon=1, klon 126 DO ilev=1, klev 127 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 128 136 ENDDO 129 ENDDO 130 131 ! IF (debutphy) THEN 132 ! CALL gather(tr_seri, tr_seri_glo) 133 ! IF (MAXVAL(tr_seri_glo).LT.1.e-30) THEN 134 !--initialising tracer concentrations to zero 135 ! DO it=1, nbtr 136 ! tr_seri(:,:,it)=0.0 137 ! ENDDO 138 ! ENDIF 139 ! ENDIF 140 137 141 138 !--initialise emission diagnostics 139 budg_emi(:,1)=0.0 142 140 budg_emi_ocs(:)=0.0 143 141 budg_emi_so2(:)=0.0 … … 145 143 budg_emi_part(:)=0.0 146 144 147 !--sulfur emission, depending on chosen scenario (flag_ sulf_emit)148 SELECT CASE(flag_ sulf_emit)145 !--sulfur emission, depending on chosen scenario (flag_emit) 146 SELECT CASE(flag_emit) 149 147 150 148 CASE(0) ! background aerosol … … 157 155 IF (year_cur==year_emit_vol(ieru).AND.mth_cur==mth_emit_vol(ieru).AND.& 158 156 day_cur>=day_emit_vol(ieru).AND.day_cur<(day_emit_vol(ieru)+injdur)) THEN 159 ! 160 ! daily injection mass emission - NL 161 m_aer_emiss_vol_daily = m_aer_emiss_vol(ieru)/(REAL(injdur)*REAL(ponde_lonlat_vol(ieru))) 162 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), & 163 164 'ponde_lonlat_vol(ieru)=',ponde_lonlat_vol(ieru),'(injdur*ponde_lonlat_vol(ieru))', & 164 (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 165 166 WRITE(lunout,*) 'IN traccoag, dlon=',dlon 166 DO i=1,klon 167 !Pinatubo eruption at 15.14N, 120.35E 168 dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes 169 WRITE(lunout,*) 'IN traccoag, dlat=',dlat_loc 170 IF ( xlat(i).GE.xlat_min_vol(ieru)-dlat_loc .AND. xlat(i).LT.xlat_max_vol(ieru)+dlat_loc .AND. & 171 xlon(i).GE.xlon_min_vol(ieru)-dlon .AND. xlon(i).LT.xlon_max_vol(ieru)+dlon ) THEN 172 ! 173 WRITE(lunout,*) 'coordinates of volcanic injection point=',xlat(i),xlon(i),day_cur,mth_cur,year_cur 174 WRITE(lunout,*) 'DD m_aer_emiss_vol_daily=',m_aer_emiss_vol_daily 175 ! compute altLMDz 176 altLMDz(:)=0.0 177 DO k=1, klev 178 zrho=pplay(i,k)/t_seri(i,k)/RD !air density in kg/m3 179 zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG !mass of layer in kg 180 zdz=zdm(k)/zrho !thickness of layer in m 181 altLMDz(k+1)=altLMDz(k)+zdz !altitude of interface 182 ENDDO 183 184 SELECT CASE(flag_sulf_emit_distrib) 185 186 CASE(0) ! Gaussian distribution 187 !compute distribution of emission to vertical model layers (based on Gaussian peak in altitude) 188 f_lay_sum=0.0 189 DO k=1, klev 190 f_lay_emiss(k)=0.0 191 DO i_int=1, n_int_alt 192 alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt) 193 f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_vol(ieru))* & 194 & exp(-0.5*((alt-altemiss_vol(ieru))/sigma_alt_vol(ieru))**2.)* & 195 & (altLMDz(k+1)-altLMDz(k))/float(n_int_alt) 196 ENDDO 197 f_lay_sum=f_lay_sum+f_lay_emiss(k) 198 ENDDO 199 200 CASE(1) ! Uniform distribution 201 ! In this case, parameter sigma_alt_vol(ieru) is considered to be half the 202 ! height of the injection, centered around altemiss_vol(ieru) 203 DO k=1, klev 204 f_lay_emiss(k)=max(min(altemiss_vol(ieru)+sigma_alt_vol(ieru),altLMDz(k+1))- & 205 & max(altemiss_vol(ieru)-sigma_alt_vol(ieru),altLMDz(k)),0.)/(2.*sigma_alt_vol(ieru)) 206 f_lay_sum=f_lay_sum+f_lay_emiss(k) 207 ENDDO 208 209 END SELECT ! End CASE over flag_sulf_emit_distrib) 210 211 WRITE(lunout,*) "IN traccoag m_aer_emiss_vol=",m_aer_emiss_vol(ieru) 212 WRITE(lunout,*) "IN traccoag f_lay_emiss=",f_lay_emiss 213 !correct for step integration error 214 f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum 215 !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss) 216 !vertically distributed emission 217 DO k=1, klev 218 ! stretch emission over one day of Pinatubo eruption 219 emission=m_aer_emiss_vol_daily*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/1./(86400.-pdtphys) 220 tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys 221 budg_emi_so2(i)=budg_emi_so2(i)+emission*zdm(k)*mSatom/mSO2mol 222 ENDDO 223 ENDIF ! emission grid cell 224 ENDDO ! klon loop 225 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 226 182 ENDIF ! emission period 227 183 ENDDO ! eruption number … … 229 185 CASE(2) ! stratospheric aerosol injections (SAI) 230 186 ! 231 DO i=1,klon 232 ! SAI standard scenario with continuous emission from 1 grid point at the equator 233 ! SAI emission on single month 234 ! SAI continuous emission o 235 dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes 236 IF ( xlat(i).GE.xlat_sai-dlat_loc .AND. xlat(i).LT.xlat_sai+dlat_loc .AND. & 237 & xlon(i).GE.xlon_sai-dlon .AND. xlon(i).LT.xlon_sai+dlon ) THEN 238 ! 239 ! compute altLMDz 240 altLMDz(:)=0.0 241 DO k=1, klev 242 zrho=pplay(i,k)/t_seri(i,k)/RD !air density in kg/m3 243 zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG !mass of layer in kg 244 zdz=zdm(k)/zrho !thickness of layer in m 245 altLMDz(k+1)=altLMDz(k)+zdz !altitude of interface 246 ENDDO 247 248 SELECT CASE(flag_sulf_emit_distrib) 249 250 CASE(0) ! Gaussian distribution 251 !compute distribution of emission to vertical model layers (based on Gaussian peak in altitude) 252 f_lay_sum=0.0 253 DO k=1, klev 254 f_lay_emiss(k)=0.0 255 DO i_int=1, n_int_alt 256 alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt) 257 f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_sai)* & 258 & exp(-0.5*((alt-altemiss_sai)/sigma_alt_sai)**2.)* & 259 & (altLMDz(k+1)-altLMDz(k))/float(n_int_alt) 260 ENDDO 261 f_lay_sum=f_lay_sum+f_lay_emiss(k) 262 ENDDO 263 264 CASE(1) ! Uniform distribution 265 f_lay_sum=0.0 266 ! In this case, parameter sigma_alt_vol(ieru) is considered to be half 267 ! the height of the injection, centered around altemiss_sai 268 DO k=1, klev 269 f_lay_emiss(k)=max(min(altemiss_sai+sigma_alt_sai,altLMDz(k+1))- & 270 & max(altemiss_sai-sigma_alt_sai,altLMDz(k)),0.)/(2.*sigma_alt_sai) 271 f_lay_sum=f_lay_sum+f_lay_emiss(k) 272 ENDDO 273 274 END SELECT ! Gaussian or uniform distribution 275 276 !correct for step integration error 277 f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum 278 !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss) 279 !vertically distributed emission 280 DO k=1, klev 281 ! stretch emission over whole year (360d) 282 emission=m_aer_emiss_sai*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/FLOAT(year_len)/86400. 283 tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys 284 budg_emi_so2(i)=budg_emi_so2(i)+emission*zdm(k)*mSatom/mSO2mol 285 ENDDO 286 287 ! !emission as monodisperse particles with 0.1um dry radius (BIN21) 288 ! !vertically distributed emission 289 ! DO k=1, klev 290 ! ! stretch emission over whole year (360d) 291 ! emission=m_aer_emiss*(mH2SO4mol/mSatom)/m_part_dry(21)/m_air_gridbox(i,k)*f_lay_emiss(k)/FLOAT(year_len)/86400. 292 ! tr_seri(i,k,id_BIN01_strat+20)=tr_seri(i,k,id_BIN01_strat+20)+emission*pdtphys 293 ! budg_emi_part(i)=budg_emi_part(i)+emission*zdm(k)*mSatom/mH2SO4mol 294 ! ENDDO 295 ENDIF ! emission grid cell 296 ENDDO ! klon loop 187 ! Computing duration of SAI in days... 188 ! ... starting from 0... 189 injdur_sai = 0 190 ! ... then adding whole years from first to (n-1)th... 191 DO yr = year_emit_sai_start, year_emit_sai_end-1 192 ! (n % 4 == 0) and (n % 100 != 0 or n % 400 == 0) 193 is_bissext = (MOD(yr,4)==0) .AND. (MOD(yr,100) /= 0 .OR. MOD(yr,400) == 0) 194 injdur_sai = injdur_sai+365+is_bissext 195 ENDDO 196 ! ... then subtracting part of the first year where no injection yet... 197 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) 198 SELECT CASE(mth_emit_sai_start) 199 CASE(2) 200 injdur_sai = injdur_sai-31 201 CASE(3) 202 injdur_sai = injdur_sai-31-28-is_bissext 203 CASE(4) 204 injdur_sai = injdur_sai-31-28-is_bissext-31 205 CASE(5) 206 injdur_sai = injdur_sai-31-28-is_bissext-31-30 207 CASE(6) 208 injdur_sai = injdur_sai-31-28-is_bissext-31-30-31 209 CASE(7) 210 injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30 211 CASE(8) 212 injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31 213 CASE(9) 214 injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31-31 215 CASE(10) 216 injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31-31-30 217 CASE(11) 218 injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31-31-30-31 219 CASE(12) 220 injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31-31-30-31-30 221 END SELECT 222 injdur_sai = injdur_sai-day_emit_sai_start+1 223 ! ... then adding part of the n-th year 224 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) 225 SELECT CASE(mth_emit_sai_end) 226 CASE(2) 227 injdur_sai = injdur_sai+31 228 CASE(3) 229 injdur_sai = injdur_sai+31+28+is_bissext 230 CASE(4) 231 injdur_sai = injdur_sai+31+28+is_bissext+31 232 CASE(5) 233 injdur_sai = injdur_sai+31+28+is_bissext+31+30 234 CASE(6) 235 injdur_sai = injdur_sai+31+28+is_bissext+31+30+31 236 CASE(7) 237 injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30 238 CASE(8) 239 injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31 240 CASE(9) 241 injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31+31 242 CASE(10) 243 injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31+31+30 244 CASE(11) 245 injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31+31+30+31 246 CASE(12) 247 injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31+31+30+31+30 248 END SELECT 249 injdur_sai = injdur_sai+day_emit_sai_end 250 ! A security: are SAI dates of injection consistent? 251 IF (injdur_sai <= 0) THEN 252 CALL abort_physic('traccoag_mod', 'Pb in SAI dates of injection.',1) 253 ENDIF 254 ! Injection in itself 255 IF (( year_emit_sai_start <= year_cur ) & 256 .AND. ( year_cur <= year_emit_sai_end ) & 257 .AND. ( mth_emit_sai_start <= mth_cur .OR. year_emit_sai_start < year_cur ) & 258 .AND. ( mth_cur <= mth_emit_sai_end .OR. year_cur < year_emit_sai_end ) & 259 .AND. ( day_emit_sai_start <= day_cur .OR. mth_emit_sai_start < mth_cur .OR. year_emit_sai_start < year_cur ) & 260 .AND. ( day_cur <= day_emit_sai_end .OR. mth_cur < mth_emit_sai_end .OR. year_cur < year_emit_sai_end )) THEN 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 281 ENDIF ! Condition over injection dates 297 282 298 283 CASE(3) ! --- SAI injection over a single band of longitude and between 299 284 ! lat_min and lat_max 300 285 301 DO i=1,klon 302 ! SAI scenario with continuous emission 303 dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes 304 theta_min = max(xlat(i)-dlat_loc,xlat_min_sai) 305 theta_max = min(xlat(i)+dlat_loc,xlat_max_sai) 306 IF ( xlat(i).GE.xlat_min_sai-dlat_loc .AND. xlat(i).LT.xlat_max_sai+dlat_loc .AND. & 307 & xlon(i).GE.xlon_sai-dlon .AND. xlon(i).LT.xlon_sai+dlon ) THEN 308 ! 309 ! compute altLMDz 310 altLMDz(:)=0.0 311 DO k=1, klev 312 zrho=pplay(i,k)/t_seri(i,k)/RD !air density in kg/m3 313 zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG !mass of layer in kg 314 zdz=zdm(k)/zrho !thickness of layer in m 315 altLMDz(k+1)=altLMDz(k)+zdz !altitude of interface 316 ENDDO 317 318 SELECT CASE(flag_sulf_emit_distrib) 319 320 CASE(0) ! Gaussian distribution 321 !compute distribution of emission to vertical model layers (based on 322 !Gaussian peak in altitude) 323 f_lay_sum=0.0 324 DO k=1, klev 325 f_lay_emiss(k)=0.0 326 DO i_int=1, n_int_alt 327 alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt) 328 f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.*RPI)*sigma_alt_sai)* & 329 & exp(-0.5*((alt-altemiss_sai)/sigma_alt_sai)**2.)* & 330 & (altLMDz(k+1)-altLMDz(k))/float(n_int_alt) 331 ENDDO 332 f_lay_sum=f_lay_sum+f_lay_emiss(k) 333 ENDDO 334 335 CASE(1) ! Uniform distribution 336 f_lay_sum=0.0 337 ! In this case, parameter sigma_alt_vol(ieru) is considered to be half 338 ! the height of the injection, centered around altemiss_sai 339 DO k=1, klev 340 f_lay_emiss(k)=max(min(altemiss_sai+sigma_alt_sai,altLMDz(k+1))- & 341 & max(altemiss_sai-sigma_alt_sai,altLMDz(k)),0.)/(2.*sigma_alt_sai) 342 f_lay_sum=f_lay_sum+f_lay_emiss(k) 343 ENDDO 344 345 END SELECT ! Gaussian or uniform distribution 346 347 !correct for step integration error 348 f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum 349 !emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss) 350 !vertically distributed emission 351 DO k=1, klev 352 ! stretch emission over whole year (360d) 353 emission=m_aer_emiss_sai*(mSO2mol/mSatom)/m_air_gridbox(i,k)*f_lay_emiss(k)/ & 354 & FLOAT(year_len)/86400.*(sin(theta_max/180.*RPI)-sin(theta_min/180.*RPI))/ & 355 & (sin(xlat_max_sai/180.*RPI)-sin(xlat_min_sai/180.*RPI)) 356 tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys 357 budg_emi_so2(i)=budg_emi_so2(i)+emission*zdm(k)*mSatom/mSO2mol 358 ENDDO 359 360 ! !emission as monodisperse particles with 0.1um dry radius (BIN21) 361 ! !vertically distributed emission 362 ! DO k=1, klev 363 ! ! stretch emission over whole year (360d) 364 ! emission=m_aer_emiss*(mH2SO4mol/mSatom)/m_part_dry(21)/m_air_gridbox(i,k)*f_lay_emiss(k)/year_len/86400 365 ! tr_seri(i,k,id_BIN01_strat+20)=tr_seri(i,k,id_BIN01_strat+20)+emission*pdtphys 366 ! budg_emi_part(i)=budg_emi_part(i)+emission*zdm(k)*mSatom/mH2SO4mol 367 ! ENDDO 368 ENDIF ! emission grid cell 369 ENDDO ! klon loop 370 371 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) 372 307 373 308 !--read background concentrations of OCS and SO2 and lifetimes from input file … … 387 322 CALL coagulate(pdtphys,mdw,tr_seri,t_seri,pplay,dens_aer,is_strato) 388 323 389 !--call sedimentation routine 324 !--call sedimentation routine 390 325 CALL aer_sedimnt(pdtphys, t_seri, pplay, paprs, tr_seri, dens_aer) 391 326 … … 404 339 ENDDO 405 340 ENDDO 406 407 ! CALL minmaxsimple(tr_seri(:,:,id_SO2_strat),0.0,0.0,'fin SO2') 408 ! DO it=1, nbtr_bin 409 ! CALL minmaxsimple(tr_seri(:,:,nbtr_sulgas+it),0.0,0.0,'fin bin ') 410 ! ENDDO 411 341 412 342 END SUBROUTINE traccoag 413 343
Note: See TracChangeset
for help on using the changeset viewer.
