Context Navigation

traccoag_mod.f90 @ 5410

Last change on this file since 5410 was 5367, checked in by abarral, 2 weeks ago
(WIP) Turn implicit into explicit declarations
Property svn:keywords set to `Id`
File size: 19.9 KB

Rev	Line
[2690]	1	MODULE traccoag_mod
	2	!
	3	! This module calculates the concentration of aerosol particles in certain size bins
	4	! considering coagulation and sedimentation.
	5	!
[5367]	6	IMPLICIT NONE; PRIVATE
	7	PUBLIC traccoag
[2690]	8	CONTAINS
	9
	10	SUBROUTINE traccoag(pdtphys, gmtime, debutphy, julien, &
	11	presnivs, xlat, xlon, pphis, pphi, &
[2752]	12	t_seri, pplay, paprs, sh, rh, tr_seri)
[4950]	13
[2752]	14	USE phys_local_var_mod, ONLY: mdw, R2SO4, DENSO4, f_r_wet, surf_PM25_sulf, &
[4950]	15	& budg_emi_ocs, budg_emi_so2, budg_emi_h2so4, budg_emi_part, &
[5150]	16	& R2SO4B, DENSO4B, f_r_wetB, sulfmmr, SAD_sulfate, sulfmmr_mode, nd_mode, reff_sulfate
[4950]	17
[2690]	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
[4601]	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
[5264]	28	USE yomcst_mod_h
[3526]	29	USE print_control_mod, ONLY: lunout
[4601]	30	USE strataer_local_var_mod
	31
[2690]	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
[4601]	60	REAL :: m_aer ! aerosol mass
[4998]	61	INTEGER :: it, k, i, j, 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
[4601]	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
[4513]	85	INTEGER :: injdur_sai ! injection duration for SAI case [days]
	86	INTEGER :: yr, is_bissext
[5150]	87	REAL :: samoment2, samoment3! 2nd and 3rd order moments of size distribution
[2690]	88
[4601]	89	IF (is_mpi_root .AND. flag_verbose_strataer) THEN
[3526]	90	WRITE(lunout,*) 'in traccoag: date from phys_cal_mod =',year_cur,'-',mth_cur,'-',day_cur,'-',hour
[4601]	91	WRITE(lunout,*) 'IN traccoag flag_emit: ',flag_emit
[2690]	92	ENDIF
[3526]	93
[4950]	94	! radius [m]
[2690]	95	DO it=1, nbtr_bin
	96	r_bin(it)=mdw(it)/2.
	97	ENDDO
	98
	99	!--set boundaries of size bins
	100	DO it=1, nbtr_bin
	101	IF (it.EQ.1) THEN
	102	r_upper(it)=sqrt(r_bin(it+1)*r_bin(it))
	103	r_lower(it)=r_bin(it)**2./r_upper(it)
	104	ELSEIF (it.EQ.nbtr_bin) THEN
	105	r_lower(it)=sqrt(r_bin(it)*r_bin(it-1))
	106	r_upper(it)=r_bin(it)**2./r_lower(it)
	107	ELSE
	108	r_lower(it)=sqrt(r_bin(it)*r_bin(it-1))
	109	r_upper(it)=sqrt(r_bin(it+1)*r_bin(it))
	110	ENDIF
	111	ENDDO
	112
	113	IF (debutphy .and. is_mpi_root) THEN
	114	DO it=1, nbtr_bin
[3526]	115	WRITE(lunout,*) 'radius bin', it, ':', r_bin(it), '(from', r_lower(it), 'to', r_upper(it), ')'
[2690]	116	ENDDO
	117	ENDIF
	118
	119	!--initialising logical is_strato from stratomask
	120	is_strato(:,:)=.FALSE.
[2695]	121	WHERE (stratomask.GT.0.5) is_strato=.TRUE.
[2690]	122
[4750]	123	IF(flag_new_strat_compo) THEN
[4950]	124	IF(debutphy) WRITE(lunout,*) 'traccoag: COMPO/DENSITY (Tabazadeh 97) + H2O kelvin effect', flag_new_strat_compo
	125	! STRACOMP (H2O, P, t_seri, R -> R2SO4 + Kelvin effect) : Taba97, Socol, etc...
	126	CALL stracomp_kelvin(sh,t_seri,pplay)
[4750]	127	ELSE
[4950]	128	IF(debutphy) WRITE(lunout,*) 'traccoag: COMPO from Bekki 2D model', flag_new_strat_compo
[4750]	129	! STRACOMP (H2O, P, t_seri -> aerosol composition (R2SO4))
	130	! H2SO4 mass fraction in aerosol (%)
	131	CALL stracomp(sh,t_seri,pplay)
	132
	133	! aerosol density (gr/cm3)
	134	CALL denh2sa(t_seri)
[4950]	135
	136	! compute factor for converting dry to wet radius (for every grid box)
	137	f_r_wet(:,:) = (dens_aer_dry/(DENSO4(:,:)1000.)/(R2SO4(:,:)/100.))*(1./3.)
[4750]	138	ENDIF
	139
[2690]	140	!--calculate mass of air in every grid box
	141	DO ilon=1, klon
[4601]	142	DO ilev=1, klev
	143	m_air_gridbox(ilon,ilev)=(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG*cell_area(ilon)
	144	ENDDO
[2690]	145	ENDDO
[4601]	146
[2752]	147	!--initialise emission diagnostics
[4769]	148	if (nErupt > 0 .and. (flag_emit == 1 .or. flag_emit == 4)) budg_emi(:,1)=0.0
[2752]	149	budg_emi_ocs(:)=0.0
	150	budg_emi_so2(:)=0.0
	151	budg_emi_h2so4(:)=0.0
	152	budg_emi_part(:)=0.0
	153
[4601]	154	!--sulfur emission, depending on chosen scenario (flag_emit)
	155	SELECT CASE(flag_emit)
[2690]	156
	157	CASE(0) ! background aerosol
	158	! do nothing (no emission)
	159
	160	CASE(1) ! volcanic eruption
	161	!--only emit on day of eruption
	162	! stretch emission over one day of Pinatubo eruption
[3526]	163	DO ieru=1, nErupt
	164	IF (year_cur==year_emit_vol(ieru).AND.mth_cur==mth_emit_vol(ieru).AND.&
	165	day_cur>=day_emit_vol(ieru).AND.day_cur<(day_emit_vol(ieru)+injdur)) THEN
[4601]	166
	167	! daily injection mass emission
	168	m_aer=m_aer_emiss_vol(ieru,1)/(REAL(injdur)*REAL(ponde_lonlat_vol(ieru)))
	169	!emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
	170	m_aer=m_aer*(mSO2mol/mSatom)
	171
	172	WRITE(lunout,*) 'IN traccoag m_aer_emiss_vol(ieru)=',m_aer_emiss_vol(ieru,1), &
[3526]	173	'ponde_lonlat_vol(ieru)=',ponde_lonlat_vol(ieru),'(injdur*ponde_lonlat_vol(ieru))', &
[4601]	174	(injdur*ponde_lonlat_vol(ieru)),'m_aer_emiss_vol_daily=',m_aer,'ieru=',ieru
[3526]	175	WRITE(lunout,*) 'IN traccoag, dlon=',dlon
[4601]	176
	177	latmin=xlat_min_vol(ieru)
	178	latmax=xlat_max_vol(ieru)
	179	lonmin=xlon_min_vol(ieru)
	180	lonmax=xlon_max_vol(ieru)
	181	altemiss = altemiss_vol(ieru)
	182	sigma_alt = sigma_alt_vol(ieru)
	183	pdt=pdtphys
	184	! stretch emission over one day of eruption
	185	stretchlong = 1.
	186
	187	CALL STRATEMIT(pdtphys,pdt,xlat,xlon,t_seri,pplay,paprs,tr_seri,&
	188	m_aer,latmin,latmax,lonmin,lonmax,altemiss,sigma_alt,id_SO2_strat, &
	189	stretchlong,1,0)
	190
[3526]	191	ENDIF ! emission period
	192	ENDDO ! eruption number
	193
[2690]	194	CASE(2) ! stratospheric aerosol injections (SAI)
	195	!
[4513]	196	! Computing duration of SAI in days...
	197	! ... starting from 0...
	198	injdur_sai = 0
	199	! ... then adding whole years from first to (n-1)th...
	200	DO yr = year_emit_sai_start, year_emit_sai_end-1
	201	! (n % 4 == 0) and (n % 100 != 0 or n % 400 == 0)
	202	is_bissext = (MOD(yr,4)==0) .AND. (MOD(yr,100) /= 0 .OR. MOD(yr,400) == 0)
	203	injdur_sai = injdur_sai+365+is_bissext
	204	ENDDO
	205	! ... then subtracting part of the first year where no injection yet...
	206	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)
	207	SELECT CASE(mth_emit_sai_start)
	208	CASE(2)
	209	injdur_sai = injdur_sai-31
	210	CASE(3)
	211	injdur_sai = injdur_sai-31-28-is_bissext
	212	CASE(4)
	213	injdur_sai = injdur_sai-31-28-is_bissext-31
	214	CASE(5)
	215	injdur_sai = injdur_sai-31-28-is_bissext-31-30
	216	CASE(6)
	217	injdur_sai = injdur_sai-31-28-is_bissext-31-30-31
	218	CASE(7)
	219	injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30
	220	CASE(8)
	221	injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31
	222	CASE(9)
	223	injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31-31
	224	CASE(10)
	225	injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31-31-30
	226	CASE(11)
	227	injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31-31-30-31
	228	CASE(12)
	229	injdur_sai = injdur_sai-31-28-is_bissext-31-30-31-30-31-31-30-31-30
	230	END SELECT
	231	injdur_sai = injdur_sai-day_emit_sai_start+1
	232	! ... then adding part of the n-th year
	233	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)
	234	SELECT CASE(mth_emit_sai_end)
	235	CASE(2)
	236	injdur_sai = injdur_sai+31
	237	CASE(3)
	238	injdur_sai = injdur_sai+31+28+is_bissext
	239	CASE(4)
	240	injdur_sai = injdur_sai+31+28+is_bissext+31
	241	CASE(5)
	242	injdur_sai = injdur_sai+31+28+is_bissext+31+30
	243	CASE(6)
	244	injdur_sai = injdur_sai+31+28+is_bissext+31+30+31
	245	CASE(7)
	246	injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30
	247	CASE(8)
	248	injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31
	249	CASE(9)
	250	injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31+31
	251	CASE(10)
	252	injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31+31+30
	253	CASE(11)
	254	injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31+31+30+31
	255	CASE(12)
	256	injdur_sai = injdur_sai+31+28+is_bissext+31+30+31+30+31+31+30+31+30
	257	END SELECT
	258	injdur_sai = injdur_sai+day_emit_sai_end
	259	! A security: are SAI dates of injection consistent?
	260	IF (injdur_sai <= 0) THEN
	261	CALL abort_physic('traccoag_mod', 'Pb in SAI dates of injection.',1)
	262	ENDIF
	263	! Injection in itself
	264	IF (( year_emit_sai_start <= year_cur ) &
	265	.AND. ( year_cur <= year_emit_sai_end ) &
	266	.AND. ( mth_emit_sai_start <= mth_cur .OR. year_emit_sai_start < year_cur ) &
	267	.AND. ( mth_cur <= mth_emit_sai_end .OR. year_cur < year_emit_sai_end ) &
	268	.AND. ( day_emit_sai_start <= day_cur .OR. mth_emit_sai_start < mth_cur .OR. year_emit_sai_start < year_cur ) &
	269	.AND. ( day_cur <= day_emit_sai_end .OR. mth_cur < mth_emit_sai_end .OR. year_cur < year_emit_sai_end )) THEN
[4601]	270
	271	m_aer=m_aer_emiss_sai
	272	!emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
	273	m_aer=m_aer*(mSO2mol/mSatom)
	274
	275	latmin=xlat_sai
	276	latmax=xlat_sai
	277	lonmin=xlon_sai
	278	lonmax=xlon_sai
	279	altemiss = altemiss_sai
	280	sigma_alt = sigma_alt_sai
	281	pdt=0.
	282	! stretch emission over whole year (360d)
	283	stretchlong=FLOAT(year_len)
	284
	285	CALL STRATEMIT(pdtphys,pdt,xlat,xlon,t_seri,pplay,paprs,m_air_gridbox,tr_seri,&
	286	m_aer,latmin,latmax,lonmin,lonmax,altemiss,sigma_alt,id_SO2_strat, &
	287	stretchlong,1,0)
	288
	289	budg_emi_so2(:) = budg_emi(:,1)*mSatom/mSO2mol
[4513]	290	ENDIF ! Condition over injection dates
	291
[3526]	292	CASE(3) ! --- SAI injection over a single band of longitude and between
	293	! lat_min and lat_max
	294
[4601]	295	m_aer=m_aer_emiss_sai
	296	!emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
	297	m_aer=m_aer*(mSO2mol/mSatom)
[3526]	298
[4601]	299	latmin=xlat_min_sai
	300	latmax=xlat_max_sai
	301	lonmin=xlon_sai
	302	lonmax=xlon_sai
	303	altemiss = altemiss_sai
	304	sigma_alt = sigma_alt_sai
	305	pdt=0.
	306	! stretch emission over whole year (360d)
	307	stretchlong=FLOAT(year_len)
[3526]	308
[4601]	309	CALL STRATEMIT(pdtphys,pdt,xlat,xlon,t_seri,pplay,paprs,m_air_gridbox,tr_seri,&
	310	m_aer,latmin,latmax,lonmin,lonmax,altemiss,sigma_alt,id_SO2_strat, &
	311	stretchlong,1,0)
[3526]	312
[4601]	313	budg_emi_so2(:) = budg_emi(:,1)*mSatom/mSO2mol
	314
	315	END SELECT ! emission scenario (flag_emit)
[3526]	316
[2690]	317	!--read background concentrations of OCS and SO2 and lifetimes from input file
[2695]	318	!--update the variables defined in phys_local_var_mod
	319	CALL interp_sulf_input(debutphy,pdtphys,paprs,tr_seri)
[2690]	320
	321	!--convert OCS to SO2 in the stratosphere
[2752]	322	CALL ocs_to_so2(pdtphys,tr_seri,t_seri,pplay,paprs,is_strato)
[2690]	323
	324	!--convert SO2 to H2SO4
[2752]	325	CALL so2_to_h2so4(pdtphys,tr_seri,t_seri,pplay,paprs,is_strato)
[2690]	326
	327	!--common routine for nucleation and condensation/evaporation with adaptive timestep
	328	CALL micphy_tstep(pdtphys,tr_seri,t_seri,pplay,paprs,rh,is_strato)
	329
	330	!--call coagulation routine
	331	CALL coagulate(pdtphys,mdw,tr_seri,t_seri,pplay,dens_aer,is_strato)
	332
[4601]	333	!--call sedimentation routine
[2690]	334	CALL aer_sedimnt(pdtphys, t_seri, pplay, paprs, tr_seri, dens_aer)
	335
	336	!--compute mass concentration of PM2.5 sulfate particles (wet diameter and mass) at the surface for health studies
	337	surf_PM25_sulf(:)=0.0
	338	DO i=1,klon
	339	DO it=1, nbtr_bin
	340	IF (mdw(it) .LT. 2.5e-6) THEN
	341	!surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas)*m_part(i,1,it) &
[3526]	342	!assume that particles consist of ammonium sulfate at the surface (132g/mol)
	343	!and are dry at T = 20 deg. C and 50 perc. humidity
[2690]	344	surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas) &
	345	& 132./98.dens_aer_dry4./3.RPI(mdw(it)/2.)*3 &
[2752]	346	& pplay(i,1)/t_seri(i,1)/RD1.e9
[2690]	347	ENDIF
	348	ENDDO
	349	ENDDO
[4601]	350
[4998]	351	!--compute
	352	! sulfmmr: Sulfate aerosol concentration (dry mixing ratio) (condensed H2SO4 mmr)
	353	! SAD_sulfate: SAD all aerosols (cm2/cm3) (must be WET)
	354	! sulfmmr_mode: sulfate(=H2SO4 if dry) MMR in different modes (ambiguous but based on sulfmmr, it mus be DRY(?) mmr)
	355	! nd_mode: DRY(?) particle concentration in different modes (part/m3)
	356	sulfmmr(:,:)=0.0
	357	SAD_sulfate(:,:)=0.0
	358	sulfmmr_mode(:,:,:)=0.0
	359	nd_mode(:,:,:)=0.0
[5150]	360	reff_sulfate(:,:)=0.0
[4998]	361
	362	DO i=1,klon
	363	DO j=1,klev
[5150]	364	samoment2=0.0
	365	samoment3=0.0
[4998]	366	DO it=1, nbtr_bin
	367	!surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas)*m_part(i,1,it) &
	368	!assume that particles consist of ammonium sulfate at the surface (132g/mol)
	369	!and are dry at T = 20 deg. C and 50 perc. humidity
	370
	371	! sulfmmr_mode: sulfate(=H2SO4 if dry) MMR in different modes (based on sulfmmr, it must be DRY mmr)
	372	! equivalent to condensed H2SO4 mmr= H2SO4 kg / kgA in bin it
	373	sulfmmr_mode(i,j,it) = tr_seri(i,j,it+nbtr_sulgas) & ! [DRY part/kgA in bin it]
	374	& (4./3.)RPI(mdw(it)/2.)*3. & ! [mdw: dry diameter in m]
	375	& *dens_aer_dry ! [dry aerosol mass density in kg/m3]
	376
	377	! sulfmmr: Sulfate aerosol concentration (dry mass mixing ratio)
	378	! equivalent to total condensed H2SO4 mmr (H2SO4 kg / kgA
	379	sulfmmr(i,j) = sulfmmr(i,j) + sulfmmr_mode(i,j,it)
	380
	381	! nd_mode: particle concentration in different modes (DRY part/m3)
	382	nd_mode(i,j,it) = tr_seri(i,j,it+nbtr_sulgas) & ! [DRY part/kgA in bin it]
	383	& *pplay(i,j)/t_seri(i,j)/RD ! [air mass concentration in kg air /m3A]
	384
	385	IF(flag_new_strat_compo) THEN
	386	! SAD_sulfate: SAD WET sulfate aerosols (cm2/cm3)
	387	SAD_sulfate(i,j) = SAD_sulfate(i,j) + nd_mode(i,j,it) & ! [DRY part/m3A (in bin it)]
	388	& 4.RPI( mdw(it)f_r_wetB(i,j,it)/2. )**2. & ! [WET SA of part it in m2]
	389	& *1.e-2 ! conversion from m2/m3 to cm2/cm3A
[5150]	390	! samoment2 : 2nd order moment of WET sulfate aerosols (m2/m3)
	391	samoment2 = samoment2 + nd_mode(i,j,it) & ! [DRY part/m3A (in bin it)]
	392	& ( mdw(it)f_r_wetB(i,j,it)/2. )**2. ! [WET SA of part it in m2]
	393	! samoment3 : 3nd order moment of WET sulfate aerosols (cm2/cm3)
	394	samoment3 = samoment3 + nd_mode(i,j,it) & ! [DRY part/m3A (in bin it)]
	395	& ( mdw(it)f_r_wetB(i,j,it)/2. )**3. ! [WET SA of part it in m2]
[4998]	396	ELSE
[5150]	397	! SAD_sulfate: SAD WET sulfate aerosols (cm2/cm3)
[4998]	398	SAD_sulfate(i,j) = SAD_sulfate(i,j) + nd_mode(i,j,it) & ! [DRY part/m3A (in bin it)]
	399	& 4.RPI( mdw(it)f_r_wet(i,j)/2. )**2. & ! [WET SA of part it in m2]
	400	& *1.e-2 ! conversion from m2/m3 to cm2/cm3A
[5150]	401	! samoment2 : 2nd order moment of WET sulfate aerosols (m2/m3)
	402	samoment2 = samoment2 + nd_mode(i,j,it) & ! [DRY part/m3A (in bin it)]
	403	& ( mdw(it)f_r_wet(i,j)/2. )**2. ! [WET SA of part it in m2]
	404	! samoment3 : 3nd order moment of WET sulfate aerosols (cm2/cm3)
	405	samoment3 = samoment3 + nd_mode(i,j,it) & ! [DRY part/m3A (in bin it)]
	406	& ( mdw(it)f_r_wet(i,j)/2. )**3. ! [WET SA of part it in m2]
[4998]	407	ENDIF
	408	ENDDO
[5150]	409	! reff_sulfate: effective radius of WET sulfate aerosols (cm)
[5239]	410	IF(samoment2 > 1.e-30) THEN
	411	reff_sulfate(i,j) = (samoment3 / samoment2) &
	412	& *1.e2 ! conversion from m to cm
	413
	414	! Sanity check
	415	IF(reff_sulfate(i,j) > 5.e-4) reff_sulfate(i,j) = 5.e-4 ! reff_sulfate max = 5 micron
	416	IF(reff_sulfate(i,j) < 1.e-6) reff_sulfate(i,j) = 1.e-6 ! reff_sulfate min = 10 nm
	417	ELSE
	418	reff_sulfate(i,j) = 1.e-5 ! reff_sulfate N part = nul (ref = 100 nm)
	419	ENDIF
[4998]	420	ENDDO
	421	ENDDO
	422
[2690]	423	END SUBROUTINE traccoag
	424
	425	END MODULE traccoag_mod

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

Download in other formats:

Original Format