Context Navigation

traccoag_mod.F90 @ 3570

Last change on this file since 3570 was 3528, checked in by oboucher, 5 years ago
A little more cleaning with WRITEs
Property svn:keywords set to `Id`
File size: 19.4 KB

Rev	Line
[3526]	1	!
	2	! $Id: traccoag_mod.F90 3528 2019-05-30 14:13:37Z fhourdin $
	3	!
[2690]	4	MODULE traccoag_mod
	5	!
	6	! This module calculates the concentration of aerosol particles in certain size bins
	7	! considering coagulation and sedimentation.
	8	!
	9	CONTAINS
	10
	11	SUBROUTINE traccoag(pdtphys, gmtime, debutphy, julien, &
	12	presnivs, xlat, xlon, pphis, pphi, &
[2752]	13	t_seri, pplay, paprs, sh, rh, tr_seri)
[2690]	14
[2752]	15	USE phys_local_var_mod, ONLY: mdw, R2SO4, DENSO4, f_r_wet, surf_PM25_sulf, &
	16	& budg_emi_ocs, budg_emi_so2, budg_emi_h2so4, budg_emi_part
[2690]	17
	18	USE dimphy
	19	USE infotrac
	20	USE aerophys
[3526]	21	USE geometry_mod, ONLY : cell_area, boundslat
[2690]	22	USE mod_grid_phy_lmdz
	23	USE mod_phys_lmdz_mpi_data, ONLY : is_mpi_root
	24	USE mod_phys_lmdz_para, only: gather, scatter
	25	USE phys_cal_mod
	26	USE sulfate_aer_mod
	27	USE phys_local_var_mod, ONLY: stratomask
	28	USE YOMCST
[3526]	29	USE print_control_mod, ONLY: lunout
	30	USE strataer_mod
	31	USE phys_cal_mod, ONLY : year_len
[2690]	32
	33	IMPLICIT NONE
	34
	35	! Input argument
	36	!---------------
	37	REAL,INTENT(IN) :: pdtphys ! Pas d'integration pour la physique (seconde)
	38	REAL,INTENT(IN) :: gmtime ! Heure courante
	39	LOGICAL,INTENT(IN) :: debutphy ! le flag de l'initialisation de la physique
	40	INTEGER,INTENT(IN) :: julien ! Jour julien
	41
	42	REAL,DIMENSION(klev),INTENT(IN) :: presnivs! pressions approximat. des milieux couches (en PA)
	43	REAL,DIMENSION(klon),INTENT(IN) :: xlat ! latitudes pour chaque point
	44	REAL,DIMENSION(klon),INTENT(IN) :: xlon ! longitudes pour chaque point
	45	REAL,DIMENSION(klon),INTENT(IN) :: pphis ! geopotentiel du sol
	46	REAL,DIMENSION(klon,klev),INTENT(IN) :: pphi ! geopotentiel de chaque couche
	47
	48	REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature
	49	REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa)
	50	REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression pour chaque inter-couche (en Pa)
	51	REAL,DIMENSION(klon,klev),INTENT(IN) :: sh ! humidite specifique
	52	REAL,DIMENSION(klon,klev),INTENT(IN) :: rh ! humidite relative
	53
	54	! Output argument
	55	!----------------
	56	REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA]
	57
	58	! Local variables
	59	!----------------
[3526]	60	REAL :: m_aer_emiss_vol_daily ! daily injection mass emission
	61	INTEGER :: it, k, i, ilon, ilev, itime, i_int, ieru
[2690]	62	LOGICAL,DIMENSION(klon,klev) :: is_strato ! true = above tropopause, false = below
	63	REAL,DIMENSION(klon,klev) :: m_air_gridbox ! mass of air in every grid box [kg]
	64	REAL,DIMENSION(klon_glo,klev,nbtr) :: tr_seri_glo ! Concentration Traceur [U/KgA]
	65	REAL,DIMENSION(klev+1) :: altLMDz ! altitude of layer interfaces in m
	66	REAL,DIMENSION(klev) :: f_lay_emiss ! fraction of emission for every vertical layer
	67	REAL :: f_lay_sum ! sum of layer emission fractions
[2699]	68	REAL :: alt ! altitude for integral calculation
[2690]	69	INTEGER,PARAMETER :: n_int_alt=10 ! number of subintervals for integration over Gaussian emission profile
	70	REAL,DIMENSION(nbtr_bin) :: r_bin ! particle radius in size bin [m]
	71	REAL,DIMENSION(nbtr_bin) :: r_lower ! particle radius at lower bin boundary [m]
	72	REAL,DIMENSION(nbtr_bin) :: r_upper ! particle radius at upper bin boundary [m]
	73	REAL,DIMENSION(nbtr_bin) :: m_part_dry ! mass of one dry particle in size bin [kg]
	74	REAL :: zrho ! Density of air [kg/m3]
	75	REAL :: zdz ! thickness of atm. model layer in m
[2752]	76	REAL,DIMENSION(klev) :: zdm ! mass of atm. model layer in kg
[2690]	77	REAL,DIMENSION(klon,klev) :: dens_aer ! density of aerosol particles [kg/m3 aerosol] with default H2SO4 mass fraction
[2752]	78	REAL :: emission ! emission
[3526]	79	REAL :: theta_min, theta_max ! for SAI computation between two latitudes
	80	REAL :: dlat_loc
[2690]	81
	82	IF (is_mpi_root) THEN
[3526]	83	WRITE(lunout,*) 'in traccoag: date from phys_cal_mod =',year_cur,'-',mth_cur,'-',day_cur,'-',hour
	84	WRITE(lunout,*) 'IN traccoag flag_sulf_emit: ',flag_sulf_emit
[2690]	85	ENDIF
[3526]	86
[2690]	87	DO it=1, nbtr_bin
	88	r_bin(it)=mdw(it)/2.
	89	ENDDO
	90
	91	!--set boundaries of size bins
	92	DO it=1, nbtr_bin
	93	IF (it.EQ.1) THEN
	94	r_upper(it)=sqrt(r_bin(it+1)*r_bin(it))
	95	r_lower(it)=r_bin(it)**2./r_upper(it)
	96	ELSEIF (it.EQ.nbtr_bin) THEN
	97	r_lower(it)=sqrt(r_bin(it)*r_bin(it-1))
	98	r_upper(it)=r_bin(it)**2./r_lower(it)
	99	ELSE
	100	r_lower(it)=sqrt(r_bin(it)*r_bin(it-1))
	101	r_upper(it)=sqrt(r_bin(it+1)*r_bin(it))
	102	ENDIF
	103	ENDDO
	104
	105	IF (debutphy .and. is_mpi_root) THEN
	106	DO it=1, nbtr_bin
[3526]	107	WRITE(lunout,*) 'radius bin', it, ':', r_bin(it), '(from', r_lower(it), 'to', r_upper(it), ')'
[2690]	108	ENDDO
	109	ENDIF
	110
	111	!--initialising logical is_strato from stratomask
	112	is_strato(:,:)=.FALSE.
[2695]	113	WHERE (stratomask.GT.0.5) is_strato=.TRUE.
[2690]	114
	115	! STRACOMP (H2O, P, t_seri -> aerosol composition (R2SO4))
	116	! H2SO4 mass fraction in aerosol (%)
	117	CALL stracomp(sh,t_seri,pplay)
	118
	119	! aerosol density (gr/cm3)
	120	CALL denh2sa(t_seri)
	121
	122	! compute factor for converting dry to wet radius (for every grid box)
	123	f_r_wet(:,:) = (dens_aer_dry/(DENSO4(:,:)1000.)/(R2SO4(:,:)/100.))*(1./3.)
	124
	125	!--calculate mass of air in every grid box
	126	DO ilon=1, klon
	127	DO ilev=1, klev
[2752]	128	m_air_gridbox(ilon,ilev)=(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG*cell_area(ilon)
[2690]	129	ENDDO
	130	ENDDO
	131
[2752]	132	! IF (debutphy) THEN
	133	! CALL gather(tr_seri, tr_seri_glo)
	134	! IF (MAXVAL(tr_seri_glo).LT.1.e-30) THEN
[2690]	135	!--initialising tracer concentrations to zero
[2752]	136	! DO it=1, nbtr
	137	! tr_seri(:,:,it)=0.0
	138	! ENDDO
	139	! ENDIF
	140	! ENDIF
[2690]	141
[2752]	142	!--initialise emission diagnostics
	143	budg_emi_ocs(:)=0.0
	144	budg_emi_so2(:)=0.0
	145	budg_emi_h2so4(:)=0.0
	146	budg_emi_part(:)=0.0
	147
[2690]	148	!--sulfur emission, depending on chosen scenario (flag_sulf_emit)
	149	SELECT CASE(flag_sulf_emit)
	150
	151	CASE(0) ! background aerosol
	152	! do nothing (no emission)
	153
	154	CASE(1) ! volcanic eruption
	155	!--only emit on day of eruption
	156	! stretch emission over one day of Pinatubo eruption
[3526]	157	DO ieru=1, nErupt
	158	IF (year_cur==year_emit_vol(ieru).AND.mth_cur==mth_emit_vol(ieru).AND.&
	159	day_cur>=day_emit_vol(ieru).AND.day_cur<(day_emit_vol(ieru)+injdur)) THEN
	160	!
	161	! daily injection mass emission - NL
	162	m_aer_emiss_vol_daily = m_aer_emiss_vol(ieru)/(REAL(injdur)*REAL(ponde_lonlat_vol(ieru)))
	163	WRITE(lunout,*) 'IN traccoag DD m_aer_emiss_vol(ieru)=',m_aer_emiss_vol(ieru), &
	164	'ponde_lonlat_vol(ieru)=',ponde_lonlat_vol(ieru),'(injdur*ponde_lonlat_vol(ieru))', &
	165	(injdur*ponde_lonlat_vol(ieru)),'m_aer_emiss_vol_daily=',m_aer_emiss_vol_daily,'ieru=',ieru
	166	WRITE(lunout,*) 'IN traccoag, dlon=',dlon
	167	DO i=1,klon
	168	!Pinatubo eruption at 15.14N, 120.35E
	169	dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes
	170	WRITE(lunout,*) 'IN traccoag, dlat=',dlat_loc
	171	IF ( xlat(i).GE.xlat_min_vol(ieru)-dlat_loc .AND. xlat(i).LT.xlat_max_vol(ieru)+dlat_loc .AND. &
	172	xlon(i).GE.xlon_min_vol(ieru)-dlon .AND. xlon(i).LT.xlon_max_vol(ieru)+dlon ) THEN
	173	!
	174	WRITE(lunout,*) 'coordinates of volcanic injection point=',xlat(i),xlon(i),day_cur,mth_cur,year_cur
	175	WRITE(lunout,*) 'DD m_aer_emiss_vol_daily=',m_aer_emiss_vol_daily
	176	! compute altLMDz
	177	altLMDz(:)=0.0
	178	DO k=1, klev
	179	zrho=pplay(i,k)/t_seri(i,k)/RD !air density in kg/m3
	180	zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG !mass of layer in kg
	181	zdz=zdm(k)/zrho !thickness of layer in m
	182	altLMDz(k+1)=altLMDz(k)+zdz !altitude of interface
	183	ENDDO
[2690]	184
[3526]	185	SELECT CASE(flag_sulf_emit_distrib)
	186
	187	CASE(0) ! Gaussian distribution
	188	!compute distribution of emission to vertical model layers (based on Gaussian peak in altitude)
	189	f_lay_sum=0.0
	190	DO k=1, klev
	191	f_lay_emiss(k)=0.0
	192	DO i_int=1, n_int_alt
	193	alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
	194	f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.RPI)sigma_alt_vol(ieru))* &
	195	& exp(-0.5((alt-altemiss_vol(ieru))/sigma_alt_vol(ieru))2.) &
	196	& (altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
	197	ENDDO
	198	f_lay_sum=f_lay_sum+f_lay_emiss(k)
	199	ENDDO
	200
	201	CASE(1) ! Uniform distribution
	202	! In this case, parameter sigma_alt_vol(ieru) is considered to be half the
	203	! height of the injection, centered around altemiss_vol(ieru)
	204	DO k=1, klev
	205	f_lay_emiss(k)=max(min(altemiss_vol(ieru)+sigma_alt_vol(ieru),altLMDz(k+1))- &
	206	& max(altemiss_vol(ieru)-sigma_alt_vol(ieru),altLMDz(k)),0.)/(2.*sigma_alt_vol(ieru))
	207	f_lay_sum=f_lay_sum+f_lay_emiss(k)
	208	ENDDO
	209
	210	END SELECT ! End CASE over flag_sulf_emit_distrib)
	211
	212	WRITE(lunout,*) "IN traccoag m_aer_emiss_vol=",m_aer_emiss_vol(ieru)
	213	WRITE(lunout,*) "IN traccoag f_lay_emiss=",f_lay_emiss
	214	!correct for step integration error
	215	f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum
	216	!emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
	217	!vertically distributed emission
	218	DO k=1, klev
	219	! stretch emission over one day of Pinatubo eruption
	220	emission=m_aer_emiss_vol_daily(mSO2mol/mSatom)/m_air_gridbox(i,k)f_lay_emiss(k)/1./(86400.-pdtphys)
	221	tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys
	222	budg_emi_so2(i)=budg_emi_so2(i)+emissionzdm(k)mSatom/mSO2mol
	223	ENDDO
	224	ENDIF ! emission grid cell
	225	ENDDO ! klon loop
	226	WRITE(lunout,*) "IN traccoag (ieru=",ieru,") m_aer_emiss_vol_daily=",m_aer_emiss_vol_daily
	227	ENDIF ! emission period
	228	ENDDO ! eruption number
	229
[2690]	230	CASE(2) ! stratospheric aerosol injections (SAI)
	231	!
	232	DO i=1,klon
	233	! SAI standard scenario with continuous emission from 1 grid point at the equator
	234	! SAI emission on single month
	235	! SAI continuous emission o
[3526]	236	dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes
	237	IF ( xlat(i).GE.xlat_sai-dlat_loc .AND. xlat(i).LT.xlat_sai+dlat_loc .AND. &
[2704]	238	& xlon(i).GE.xlon_sai-dlon .AND. xlon(i).LT.xlon_sai+dlon ) THEN
[2752]	239	!
[2690]	240	! compute altLMDz
	241	altLMDz(:)=0.0
	242	DO k=1, klev
[2752]	243	zrho=pplay(i,k)/t_seri(i,k)/RD !air density in kg/m3
	244	zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG !mass of layer in kg
	245	zdz=zdm(k)/zrho !thickness of layer in m
	246	altLMDz(k+1)=altLMDz(k)+zdz !altitude of interface
[2690]	247	ENDDO
[3526]	248
	249	SELECT CASE(flag_sulf_emit_distrib)
	250
	251	CASE(0) ! Gaussian distribution
[2690]	252	!compute distribution of emission to vertical model layers (based on Gaussian peak in altitude)
	253	f_lay_sum=0.0
[3526]	254	DO k=1, klev
	255	f_lay_emiss(k)=0.0
	256	DO i_int=1, n_int_alt
	257	alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
	258	f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.RPI)sigma_alt_sai)* &
	259	& exp(-0.5((alt-altemiss_sai)/sigma_alt_sai)2.) &
	260	& (altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
	261	ENDDO
	262	f_lay_sum=f_lay_sum+f_lay_emiss(k)
	263	ENDDO
	264
	265	CASE(1) ! Uniform distribution
	266	f_lay_sum=0.0
	267	! In this case, parameter sigma_alt_vol(ieru) is considered to be half
	268	! the height of the injection, centered around altemiss_sai
	269	DO k=1, klev
	270	f_lay_emiss(k)=max(min(altemiss_sai+sigma_alt_sai,altLMDz(k+1))- &
	271	& max(altemiss_sai-sigma_alt_sai,altLMDz(k)),0.)/(2.*sigma_alt_sai)
	272	f_lay_sum=f_lay_sum+f_lay_emiss(k)
	273	ENDDO
	274
	275	END SELECT ! Gaussian or uniform distribution
	276
	277	!correct for step integration error
	278	f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum
	279	!emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
	280	!vertically distributed emission
[2690]	281	DO k=1, klev
[3526]	282	! stretch emission over whole year (360d)
[3527]	283	emission=m_aer_emiss_sai(mSO2mol/mSatom)/m_air_gridbox(i,k)f_lay_emiss(k)/FLOAT(year_len)/86400.
[3526]	284	tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys
	285	budg_emi_so2(i)=budg_emi_so2(i)+emissionzdm(k)mSatom/mSO2mol
[2690]	286	ENDDO
[3526]	287
	288	! !emission as monodisperse particles with 0.1um dry radius (BIN21)
	289	! !vertically distributed emission
	290	! DO k=1, klev
	291	! ! stretch emission over whole year (360d)
[3527]	292	! emission=m_aer_emiss(mH2SO4mol/mSatom)/m_part_dry(21)/m_air_gridbox(i,k)f_lay_emiss(k)/FLOAT(year_len)/86400.
[3526]	293	! tr_seri(i,k,id_BIN01_strat+20)=tr_seri(i,k,id_BIN01_strat+20)+emission*pdtphys
	294	! budg_emi_part(i)=budg_emi_part(i)+emissionzdm(k)mSatom/mH2SO4mol
	295	! ENDDO
	296	ENDIF ! emission grid cell
	297	ENDDO ! klon loop
	298
	299	CASE(3) ! --- SAI injection over a single band of longitude and between
	300	! lat_min and lat_max
	301
	302	DO i=1,klon
	303	! SAI scenario with continuous emission
	304	dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) ! dlat = half difference of boundary latitudes
	305	theta_min = max(xlat(i)-dlat_loc,xlat_min_sai)
	306	theta_max = min(xlat(i)+dlat_loc,xlat_max_sai)
	307	IF ( xlat(i).GE.xlat_min_sai-dlat_loc .AND. xlat(i).LT.xlat_max_sai+dlat_loc .AND. &
	308	& xlon(i).GE.xlon_sai-dlon .AND. xlon(i).LT.xlon_sai+dlon ) THEN
	309	!
	310	! compute altLMDz
	311	altLMDz(:)=0.0
	312	DO k=1, klev
	313	zrho=pplay(i,k)/t_seri(i,k)/RD !air density in kg/m3
	314	zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG !mass of layer in kg
	315	zdz=zdm(k)/zrho !thickness of layer in m
	316	altLMDz(k+1)=altLMDz(k)+zdz !altitude of interface
	317	ENDDO
	318
	319	SELECT CASE(flag_sulf_emit_distrib)
	320
	321	CASE(0) ! Gaussian distribution
	322	!compute distribution of emission to vertical model layers (based on
	323	!Gaussian peak in altitude)
	324	f_lay_sum=0.0
	325	DO k=1, klev
	326	f_lay_emiss(k)=0.0
	327	DO i_int=1, n_int_alt
	328	alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
	329	f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.RPI)sigma_alt_sai)* &
	330	& exp(-0.5((alt-altemiss_sai)/sigma_alt_sai)2.) &
	331	& (altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
	332	ENDDO
	333	f_lay_sum=f_lay_sum+f_lay_emiss(k)
	334	ENDDO
	335
	336	CASE(1) ! Uniform distribution
	337	f_lay_sum=0.0
	338	! In this case, parameter sigma_alt_vol(ieru) is considered to be half
	339	! the height of the injection, centered around altemiss_sai
	340	DO k=1, klev
	341	f_lay_emiss(k)=max(min(altemiss_sai+sigma_alt_sai,altLMDz(k+1))- &
	342	& max(altemiss_sai-sigma_alt_sai,altLMDz(k)),0.)/(2.*sigma_alt_sai)
	343	f_lay_sum=f_lay_sum+f_lay_emiss(k)
	344	ENDDO
	345
	346	END SELECT ! Gaussian or uniform distribution
	347
[2690]	348	!correct for step integration error
	349	f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum
	350	!emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
	351	!vertically distributed emission
	352	DO k=1, klev
[2752]	353	! stretch emission over whole year (360d)
[3526]	354	emission=m_aer_emiss_sai(mSO2mol/mSatom)/m_air_gridbox(i,k)f_lay_emiss(k)/ &
[3527]	355	& FLOAT(year_len)/86400.(sin(theta_max/180.RPI)-sin(theta_min/180.*RPI))/ &
[3526]	356	& (sin(xlat_max_sai/180.RPI)-sin(xlat_min_sai/180.RPI))
[2752]	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)+emissionzdm(k)mSatom/mSO2mol
[2690]	359	ENDDO
[3526]	360
[2690]	361	! !emission as monodisperse particles with 0.1um dry radius (BIN21)
	362	! !vertically distributed emission
	363	! DO k=1, klev
[2752]	364	! ! stretch emission over whole year (360d)
[3526]	365	! emission=m_aer_emiss(mH2SO4mol/mSatom)/m_part_dry(21)/m_air_gridbox(i,k)f_lay_emiss(k)/year_len/86400
[2752]	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)+emissionzdm(k)mSatom/mH2SO4mol
[2690]	368	! ENDDO
	369	ENDIF ! emission grid cell
	370	ENDDO ! klon loop
	371
	372	END SELECT ! emission scenario (flag_sulf_emit)
	373
	374	!--read background concentrations of OCS and SO2 and lifetimes from input file
[2695]	375	!--update the variables defined in phys_local_var_mod
	376	CALL interp_sulf_input(debutphy,pdtphys,paprs,tr_seri)
[2690]	377
	378	!--convert OCS to SO2 in the stratosphere
[2752]	379	CALL ocs_to_so2(pdtphys,tr_seri,t_seri,pplay,paprs,is_strato)
[2690]	380
	381	!--convert SO2 to H2SO4
[2752]	382	CALL so2_to_h2so4(pdtphys,tr_seri,t_seri,pplay,paprs,is_strato)
[2690]	383
	384	!--common routine for nucleation and condensation/evaporation with adaptive timestep
	385	CALL micphy_tstep(pdtphys,tr_seri,t_seri,pplay,paprs,rh,is_strato)
	386
	387	!--call coagulation routine
	388	CALL coagulate(pdtphys,mdw,tr_seri,t_seri,pplay,dens_aer,is_strato)
	389
	390	!--call sedimentation routine
	391	CALL aer_sedimnt(pdtphys, t_seri, pplay, paprs, tr_seri, dens_aer)
	392
	393	!--compute mass concentration of PM2.5 sulfate particles (wet diameter and mass) at the surface for health studies
	394	surf_PM25_sulf(:)=0.0
	395	DO i=1,klon
	396	DO it=1, nbtr_bin
	397	IF (mdw(it) .LT. 2.5e-6) THEN
	398	!surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas)*m_part(i,1,it) &
[3526]	399	!assume that particles consist of ammonium sulfate at the surface (132g/mol)
	400	!and are dry at T = 20 deg. C and 50 perc. humidity
[2690]	401	surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas) &
	402	& 132./98.dens_aer_dry4./3.RPI(mdw(it)/2.)*3 &
[2752]	403	& pplay(i,1)/t_seri(i,1)/RD1.e9
[2690]	404	ENDIF
	405	ENDDO
	406	ENDDO
	407
[2695]	408	! CALL minmaxsimple(tr_seri(:,:,id_SO2_strat),0.0,0.0,'fin SO2')
	409	! DO it=1, nbtr_bin
	410	! CALL minmaxsimple(tr_seri(:,:,nbtr_sulgas+it),0.0,0.0,'fin bin ')
	411	! ENDDO
	412
[2690]	413	END SUBROUTINE traccoag
	414
	415	END MODULE traccoag_mod

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format