Context Navigation

traccoag_mod.F90 @ 4441

Last change on this file since 4441 was 4293, checked in by dcugnet, 2 years ago
Commit for Nicolas: fixes for StratAer?.
Property svn:keywords set to `Id`
File size: 19.5 KB

Rev	Line
[3526]	1	!
	2	! $Id: traccoag_mod.F90 4293 2022-10-04 19:42:41Z 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
[4293]	19	USE infotrac_phy, ONLY : nbtr_bin, nbtr_sulgas, nbtr, id_SO2_strat
[2690]	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
[4293]	25	USE phys_cal_mod, ONLY : year_len, year_cur,mth_cur, day_cur, hour
[2690]	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
[2690]	31
	32	IMPLICIT NONE
	33
	34	! Input argument
	35	!---------------
	36	REAL,INTENT(IN) :: pdtphys ! Pas d'integration pour la physique (seconde)
	37	REAL,INTENT(IN) :: gmtime ! Heure courante
	38	LOGICAL,INTENT(IN) :: debutphy ! le flag de l'initialisation de la physique
	39	INTEGER,INTENT(IN) :: julien ! Jour julien
	40
	41	REAL,DIMENSION(klev),INTENT(IN) :: presnivs! pressions approximat. des milieux couches (en PA)
	42	REAL,DIMENSION(klon),INTENT(IN) :: xlat ! latitudes pour chaque point
	43	REAL,DIMENSION(klon),INTENT(IN) :: xlon ! longitudes pour chaque point
	44	REAL,DIMENSION(klon),INTENT(IN) :: pphis ! geopotentiel du sol
	45	REAL,DIMENSION(klon,klev),INTENT(IN) :: pphi ! geopotentiel de chaque couche
	46
	47	REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature
	48	REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa)
	49	REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression pour chaque inter-couche (en Pa)
	50	REAL,DIMENSION(klon,klev),INTENT(IN) :: sh ! humidite specifique
	51	REAL,DIMENSION(klon,klev),INTENT(IN) :: rh ! humidite relative
	52
	53	! Output argument
	54	!----------------
	55	REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA]
	56
	57	! Local variables
	58	!----------------
[3526]	59	REAL :: m_aer_emiss_vol_daily ! daily injection mass emission
	60	INTEGER :: it, k, i, ilon, ilev, itime, i_int, ieru
[2690]	61	LOGICAL,DIMENSION(klon,klev) :: is_strato ! true = above tropopause, false = below
	62	REAL,DIMENSION(klon,klev) :: m_air_gridbox ! mass of air in every grid box [kg]
	63	REAL,DIMENSION(klon_glo,klev,nbtr) :: tr_seri_glo ! Concentration Traceur [U/KgA]
	64	REAL,DIMENSION(klev+1) :: altLMDz ! altitude of layer interfaces in m
	65	REAL,DIMENSION(klev) :: f_lay_emiss ! fraction of emission for every vertical layer
	66	REAL :: f_lay_sum ! sum of layer emission fractions
[2699]	67	REAL :: alt ! altitude for integral calculation
[2690]	68	INTEGER,PARAMETER :: n_int_alt=10 ! number of subintervals for integration over Gaussian emission profile
	69	REAL,DIMENSION(nbtr_bin) :: r_bin ! particle radius in size bin [m]
	70	REAL,DIMENSION(nbtr_bin) :: r_lower ! particle radius at lower bin boundary [m]
	71	REAL,DIMENSION(nbtr_bin) :: r_upper ! particle radius at upper bin boundary [m]
	72	REAL,DIMENSION(nbtr_bin) :: m_part_dry ! mass of one dry particle in size bin [kg]
	73	REAL :: zrho ! Density of air [kg/m3]
	74	REAL :: zdz ! thickness of atm. model layer in m
[2752]	75	REAL,DIMENSION(klev) :: zdm ! mass of atm. model layer in kg
[2690]	76	REAL,DIMENSION(klon,klev) :: dens_aer ! density of aerosol particles [kg/m3 aerosol] with default H2SO4 mass fraction
[2752]	77	REAL :: emission ! emission
[3526]	78	REAL :: theta_min, theta_max ! for SAI computation between two latitudes
	79	REAL :: dlat_loc
[2690]	80
	81	IF (is_mpi_root) THEN
[3526]	82	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_emit
[2690]	84	ENDIF
[3526]	85
[2690]	86	DO it=1, nbtr_bin
	87	r_bin(it)=mdw(it)/2.
	88	ENDDO
	89
	90	!--set boundaries of size bins
	91	DO it=1, nbtr_bin
	92	IF (it.EQ.1) THEN
	93	r_upper(it)=sqrt(r_bin(it+1)*r_bin(it))
	94	r_lower(it)=r_bin(it)**2./r_upper(it)
	95	ELSEIF (it.EQ.nbtr_bin) THEN
	96	r_lower(it)=sqrt(r_bin(it)*r_bin(it-1))
	97	r_upper(it)=r_bin(it)**2./r_lower(it)
	98	ELSE
	99	r_lower(it)=sqrt(r_bin(it)*r_bin(it-1))
	100	r_upper(it)=sqrt(r_bin(it+1)*r_bin(it))
	101	ENDIF
	102	ENDDO
	103
	104	IF (debutphy .and. is_mpi_root) THEN
	105	DO it=1, nbtr_bin
[3526]	106	WRITE(lunout,*) 'radius bin', it, ':', r_bin(it), '(from', r_lower(it), 'to', r_upper(it), ')'
[2690]	107	ENDDO
	108	ENDIF
	109
	110	!--initialising logical is_strato from stratomask
	111	is_strato(:,:)=.FALSE.
[2695]	112	WHERE (stratomask.GT.0.5) is_strato=.TRUE.
[2690]	113
	114	! STRACOMP (H2O, P, t_seri -> aerosol composition (R2SO4))
	115	! H2SO4 mass fraction in aerosol (%)
	116	CALL stracomp(sh,t_seri,pplay)
	117
	118	! aerosol density (gr/cm3)
	119	CALL denh2sa(t_seri)
	120
	121	! compute factor for converting dry to wet radius (for every grid box)
	122	f_r_wet(:,:) = (dens_aer_dry/(DENSO4(:,:)1000.)/(R2SO4(:,:)/100.))*(1./3.)
	123
	124	!--calculate mass of air in every grid box
	125	DO ilon=1, klon
	126	DO ilev=1, klev
[2752]	127	m_air_gridbox(ilon,ilev)=(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG*cell_area(ilon)
[2690]	128	ENDDO
	129	ENDDO
	130
[2752]	131	! IF (debutphy) THEN
	132	! CALL gather(tr_seri, tr_seri_glo)
	133	! IF (MAXVAL(tr_seri_glo).LT.1.e-30) THEN
[2690]	134	!--initialising tracer concentrations to zero
[2752]	135	! DO it=1, nbtr
	136	! tr_seri(:,:,it)=0.0
	137	! ENDDO
	138	! ENDIF
	139	! ENDIF
[2690]	140
[2752]	141	!--initialise emission diagnostics
	142	budg_emi_ocs(:)=0.0
	143	budg_emi_so2(:)=0.0
	144	budg_emi_h2so4(:)=0.0
	145	budg_emi_part(:)=0.0
	146
[2690]	147	!--sulfur emission, depending on chosen scenario (flag_sulf_emit)
	148	SELECT CASE(flag_sulf_emit)
	149
	150	CASE(0) ! background aerosol
	151	! do nothing (no emission)
	152
	153	CASE(1) ! volcanic eruption
	154	!--only emit on day of eruption
	155	! stretch emission over one day of Pinatubo eruption
[3526]	156	DO ieru=1, nErupt
	157	IF (year_cur==year_emit_vol(ieru).AND.mth_cur==mth_emit_vol(ieru).AND.&
	158	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), &
	163	'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=',ieru
	165	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
[2690]	183
[3526]	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)+emissionzdm(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
	226	ENDIF ! emission period
	227	ENDDO ! eruption number
	228
[2690]	229	CASE(2) ! stratospheric aerosol injections (SAI)
	230	!
	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
[3526]	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. &
[2704]	237	& xlon(i).GE.xlon_sai-dlon .AND. xlon(i).LT.xlon_sai+dlon ) THEN
[2752]	238	!
[2690]	239	! compute altLMDz
	240	altLMDz(:)=0.0
	241	DO k=1, klev
[2752]	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
[2690]	246	ENDDO
[3526]	247
	248	SELECT CASE(flag_sulf_emit_distrib)
	249
	250	CASE(0) ! Gaussian distribution
[2690]	251	!compute distribution of emission to vertical model layers (based on Gaussian peak in altitude)
	252	f_lay_sum=0.0
[3526]	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
[2690]	280	DO k=1, klev
[3526]	281	! stretch emission over whole year (360d)
[3527]	282	emission=m_aer_emiss_sai(mSO2mol/mSatom)/m_air_gridbox(i,k)f_lay_emiss(k)/FLOAT(year_len)/86400.
[3526]	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)+emissionzdm(k)mSatom/mSO2mol
[2690]	285	ENDDO
[3526]	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)
[3527]	291	! emission=m_aer_emiss(mH2SO4mol/mSatom)/m_part_dry(21)/m_air_gridbox(i,k)f_lay_emiss(k)/FLOAT(year_len)/86400.
[3526]	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)+emissionzdm(k)mSatom/mH2SO4mol
	294	! ENDDO
	295	ENDIF ! emission grid cell
	296	ENDDO ! klon loop
	297
	298	CASE(3) ! --- SAI injection over a single band of longitude and between
	299	! lat_min and lat_max
	300
	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
[2690]	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
[2752]	352	! stretch emission over whole year (360d)
[3526]	353	emission=m_aer_emiss_sai(mSO2mol/mSatom)/m_air_gridbox(i,k)f_lay_emiss(k)/ &
[3527]	354	& FLOAT(year_len)/86400.(sin(theta_max/180.RPI)-sin(theta_min/180.*RPI))/ &
[3526]	355	& (sin(xlat_max_sai/180.RPI)-sin(xlat_min_sai/180.RPI))
[2752]	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)+emissionzdm(k)mSatom/mSO2mol
[2690]	358	ENDDO
[3526]	359
[2690]	360	! !emission as monodisperse particles with 0.1um dry radius (BIN21)
	361	! !vertically distributed emission
	362	! DO k=1, klev
[2752]	363	! ! stretch emission over whole year (360d)
[3526]	364	! emission=m_aer_emiss(mH2SO4mol/mSatom)/m_part_dry(21)/m_air_gridbox(i,k)f_lay_emiss(k)/year_len/86400
[2752]	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)+emissionzdm(k)mSatom/mH2SO4mol
[2690]	367	! ENDDO
	368	ENDIF ! emission grid cell
	369	ENDDO ! klon loop
	370
	371	END SELECT ! emission scenario (flag_sulf_emit)
	372
	373	!--read background concentrations of OCS and SO2 and lifetimes from input file
[2695]	374	!--update the variables defined in phys_local_var_mod
	375	CALL interp_sulf_input(debutphy,pdtphys,paprs,tr_seri)
[2690]	376
	377	!--convert OCS to SO2 in the stratosphere
[2752]	378	CALL ocs_to_so2(pdtphys,tr_seri,t_seri,pplay,paprs,is_strato)
[2690]	379
	380	!--convert SO2 to H2SO4
[2752]	381	CALL so2_to_h2so4(pdtphys,tr_seri,t_seri,pplay,paprs,is_strato)
[2690]	382
	383	!--common routine for nucleation and condensation/evaporation with adaptive timestep
	384	CALL micphy_tstep(pdtphys,tr_seri,t_seri,pplay,paprs,rh,is_strato)
	385
	386	!--call coagulation routine
	387	CALL coagulate(pdtphys,mdw,tr_seri,t_seri,pplay,dens_aer,is_strato)
	388
	389	!--call sedimentation routine
	390	CALL aer_sedimnt(pdtphys, t_seri, pplay, paprs, tr_seri, dens_aer)
	391
	392	!--compute mass concentration of PM2.5 sulfate particles (wet diameter and mass) at the surface for health studies
	393	surf_PM25_sulf(:)=0.0
	394	DO i=1,klon
	395	DO it=1, nbtr_bin
	396	IF (mdw(it) .LT. 2.5e-6) THEN
	397	!surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas)*m_part(i,1,it) &
[3526]	398	!assume that particles consist of ammonium sulfate at the surface (132g/mol)
	399	!and are dry at T = 20 deg. C and 50 perc. humidity
[2690]	400	surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas) &
	401	& 132./98.dens_aer_dry4./3.RPI(mdw(it)/2.)*3 &
[2752]	402	& pplay(i,1)/t_seri(i,1)/RD1.e9
[2690]	403	ENDIF
	404	ENDDO
	405	ENDDO
	406
[2695]	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
[2690]	412	END SUBROUTINE traccoag
	413
	414	END MODULE traccoag_mod

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

Download in other formats:

Original Format