Changeset 4601 for LMDZ6/trunk/libf/phylmd/StratAer/traccoag_mod.F90
- Timestamp:
- Jul 1, 2023, 12:07:30 AM (11 months ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
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.
