Context Navigation

traccoag_mod.F90 @ 4400

Last change on this file since 4400 was 2787, checked in by Laurent Fairhead, 8 years ago
Merged trunk changes r2727:2785 into testing branch
File size: 13.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	!
	6	CONTAINS
	7
	8	SUBROUTINE traccoag(pdtphys, gmtime, debutphy, julien, &
	9	presnivs, xlat, xlon, pphis, pphi, &
[2787]	10	t_seri, pplay, paprs, sh, rh, tr_seri)
[2690]	11
[2787]	12	USE phys_local_var_mod, ONLY: mdw, R2SO4, DENSO4, f_r_wet, surf_PM25_sulf, &
	13	& budg_emi_ocs, budg_emi_so2, budg_emi_h2so4, budg_emi_part
[2690]	14
	15	USE dimphy
	16	USE infotrac
	17	USE aerophys
	18	USE geometry_mod, ONLY : cell_area
	19	USE mod_grid_phy_lmdz
	20	USE mod_phys_lmdz_mpi_data, ONLY : is_mpi_root
	21	USE mod_phys_lmdz_para, only: gather, scatter
	22	USE phys_cal_mod
	23	USE sulfate_aer_mod
	24	USE phys_local_var_mod, ONLY: stratomask
	25	USE YOMCST
	26
	27	IMPLICIT NONE
	28
	29	! Input argument
	30	!---------------
	31	REAL,INTENT(IN) :: pdtphys ! Pas d'integration pour la physique (seconde)
	32	REAL,INTENT(IN) :: gmtime ! Heure courante
	33	LOGICAL,INTENT(IN) :: debutphy ! le flag de l'initialisation de la physique
	34	INTEGER,INTENT(IN) :: julien ! Jour julien
	35
	36	REAL,DIMENSION(klev),INTENT(IN) :: presnivs! pressions approximat. des milieux couches (en PA)
	37	REAL,DIMENSION(klon),INTENT(IN) :: xlat ! latitudes pour chaque point
	38	REAL,DIMENSION(klon),INTENT(IN) :: xlon ! longitudes pour chaque point
	39	REAL,DIMENSION(klon),INTENT(IN) :: pphis ! geopotentiel du sol
	40	REAL,DIMENSION(klon,klev),INTENT(IN) :: pphi ! geopotentiel de chaque couche
	41
	42	REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature
	43	REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa)
	44	REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression pour chaque inter-couche (en Pa)
	45	REAL,DIMENSION(klon,klev),INTENT(IN) :: sh ! humidite specifique
	46	REAL,DIMENSION(klon,klev),INTENT(IN) :: rh ! humidite relative
	47
	48	! Output argument
	49	!----------------
	50	REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA]
	51
	52	! Local variables
	53	!----------------
	54	! flag for sulfur emission scenario: (0) background aerosol ; (1) volcanic eruption ; (2) stratospheric aerosol injections (SAI)
	55	INTEGER,PARAMETER :: flag_sulf_emit=2
	56	!
	57	!--flag_sulf_emit=1 --example Pinatubo
	58	INTEGER,PARAMETER :: year_emit_vol=1991 ! year of emission date
	59	INTEGER,PARAMETER :: mth_emit_vol=6 ! month of emission date
	60	INTEGER,PARAMETER :: day_emit_vol=15 ! day of emission date
	61	REAL,PARAMETER :: m_aer_emiss_vol=7.e9 ! emitted sulfur mass in kgS, e.g. 7Tg(S)=14Tg(SO2)
	62	REAL,PARAMETER :: altemiss_vol=17.e3 ! emission altitude in m
	63	REAL,PARAMETER :: sigma_alt_vol=1.e3 ! standard deviation of emission altitude in m
	64	REAL,PARAMETER :: xlat_vol=15.14 ! latitude of volcano in degree
	65	REAL,PARAMETER :: xlon_vol=120.35 ! longitude of volcano in degree
	66
	67	!--flag_sulf_emit=2 --SAI
	68	REAL,PARAMETER :: m_aer_emiss_sai=1.e10 ! emitted sulfur mass in kgS, eg 1e9=1TgS, 1e10=10TgS
	69	REAL,PARAMETER :: altemiss_sai=17.e3 ! emission altitude in m
	70	REAL,PARAMETER :: sigma_alt_sai=1.e3 ! standard deviation of emission altitude in m
[2787]	71	REAL,PARAMETER :: xlat_sai=0.01 ! latitude of SAI in degree
[2690]	72	REAL,PARAMETER :: xlon_sai=120.35 ! longitude of SAI in degree
	73
	74	!--other local variables
	75	INTEGER :: it, k, i, ilon, ilev, itime, i_int
	76	LOGICAL,DIMENSION(klon,klev) :: is_strato ! true = above tropopause, false = below
	77	REAL,DIMENSION(klon,klev) :: m_air_gridbox ! mass of air in every grid box [kg]
	78	REAL,DIMENSION(klon_glo,klev,nbtr) :: tr_seri_glo ! Concentration Traceur [U/KgA]
	79	REAL,DIMENSION(klev+1) :: altLMDz ! altitude of layer interfaces in m
	80	REAL,DIMENSION(klev) :: f_lay_emiss ! fraction of emission for every vertical layer
	81	REAL :: f_lay_sum ! sum of layer emission fractions
[2699]	82	REAL :: alt ! altitude for integral calculation
[2690]	83	INTEGER,PARAMETER :: n_int_alt=10 ! number of subintervals for integration over Gaussian emission profile
	84	REAL,DIMENSION(nbtr_bin) :: r_bin ! particle radius in size bin [m]
	85	REAL,DIMENSION(nbtr_bin) :: r_lower ! particle radius at lower bin boundary [m]
	86	REAL,DIMENSION(nbtr_bin) :: r_upper ! particle radius at upper bin boundary [m]
	87	REAL,DIMENSION(nbtr_bin) :: m_part_dry ! mass of one dry particle in size bin [kg]
	88	REAL :: zrho ! Density of air [kg/m3]
	89	REAL :: zdz ! thickness of atm. model layer in m
[2787]	90	REAL,DIMENSION(klev) :: zdm ! mass of atm. model layer in kg
[2690]	91	REAL,DIMENSION(klon,klev) :: dens_aer ! density of aerosol particles [kg/m3 aerosol] with default H2SO4 mass fraction
[2704]	92	REAL :: dlat, dlon ! d latitude and d longitude of grid in degree
[2787]	93	REAL :: emission ! emission
[2690]	94
	95	IF (is_mpi_root) THEN
	96	PRINT *,'in traccoag: date from phys_cal_mod =',year_cur,'-',mth_cur,'-',day_cur,'-',hour
	97	ENDIF
	98
[2704]	99	dlat=180./2./FLOAT(nbp_lat) ! d latitude in degree
	100	dlon=360./2./FLOAT(nbp_lon) ! d longitude in degree
	101
[2690]	102	DO it=1, nbtr_bin
	103	r_bin(it)=mdw(it)/2.
	104	ENDDO
	105
	106	!--set boundaries of size bins
	107	DO it=1, nbtr_bin
	108	IF (it.EQ.1) THEN
	109	r_upper(it)=sqrt(r_bin(it+1)*r_bin(it))
	110	r_lower(it)=r_bin(it)**2./r_upper(it)
	111	ELSEIF (it.EQ.nbtr_bin) THEN
	112	r_lower(it)=sqrt(r_bin(it)*r_bin(it-1))
	113	r_upper(it)=r_bin(it)**2./r_lower(it)
	114	ELSE
	115	r_lower(it)=sqrt(r_bin(it)*r_bin(it-1))
	116	r_upper(it)=sqrt(r_bin(it+1)*r_bin(it))
	117	ENDIF
	118	ENDDO
	119
	120	IF (debutphy .and. is_mpi_root) THEN
	121	DO it=1, nbtr_bin
	122	PRINT *,'radius bin', it, ':', r_bin(it), '(from', r_lower(it), 'to', r_upper(it), ')'
	123	ENDDO
	124	ENDIF
	125
	126	!--initialising logical is_strato from stratomask
	127	is_strato(:,:)=.FALSE.
[2695]	128	WHERE (stratomask.GT.0.5) is_strato=.TRUE.
[2690]	129
	130	! STRACOMP (H2O, P, t_seri -> aerosol composition (R2SO4))
	131	! H2SO4 mass fraction in aerosol (%)
	132	CALL stracomp(sh,t_seri,pplay)
	133
	134	! aerosol density (gr/cm3)
	135	CALL denh2sa(t_seri)
	136
	137	! compute factor for converting dry to wet radius (for every grid box)
	138	f_r_wet(:,:) = (dens_aer_dry/(DENSO4(:,:)1000.)/(R2SO4(:,:)/100.))*(1./3.)
	139
	140	!--calculate mass of air in every grid box
	141	DO ilon=1, klon
	142	DO ilev=1, klev
[2787]	143	m_air_gridbox(ilon,ilev)=(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG*cell_area(ilon)
[2690]	144	ENDDO
	145	ENDDO
	146
[2787]	147	! IF (debutphy) THEN
	148	! CALL gather(tr_seri, tr_seri_glo)
	149	! IF (MAXVAL(tr_seri_glo).LT.1.e-30) THEN
[2690]	150	!--initialising tracer concentrations to zero
[2787]	151	! DO it=1, nbtr
	152	! tr_seri(:,:,it)=0.0
	153	! ENDDO
	154	! ENDIF
	155	! ENDIF
[2690]	156
[2787]	157	!--initialise emission diagnostics
	158	budg_emi_ocs(:)=0.0
	159	budg_emi_so2(:)=0.0
	160	budg_emi_h2so4(:)=0.0
	161	budg_emi_part(:)=0.0
	162
[2690]	163	!--sulfur emission, depending on chosen scenario (flag_sulf_emit)
	164	SELECT CASE(flag_sulf_emit)
	165
	166	CASE(0) ! background aerosol
	167	! do nothing (no emission)
	168
	169	CASE(1) ! volcanic eruption
	170	!--only emit on day of eruption
	171	! stretch emission over one day of Pinatubo eruption
	172	IF (year_cur==year_emit_vol.AND.mth_cur==mth_emit_vol.AND.day_cur==day_emit_vol) THEN
	173	!
	174	DO i=1,klon
	175	!Pinatubo eruption at 15.14N, 120.35E
[2704]	176	IF ( xlat(i).GE.xlat_vol-dlat .AND. xlat(i).LT.xlat_vol+dlat .AND. &
	177	xlon(i).GE.xlon_vol-dlon .AND. xlon(i).LT.xlon_vol+dlon ) THEN
[2787]	178	!
	179	PRINT *,'coordinates of volcanic injection point=',xlat(i), xlon(i), day_cur, mth_cur, year_cur
[2690]	180	! compute altLMDz
	181	altLMDz(:)=0.0
	182	DO k=1, klev
[2787]	183	zrho=pplay(i,k)/t_seri(i,k)/RD !air density in kg/m3
	184	zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG !mass of layer in kg
	185	zdz=zdm(k)/zrho !thickness of layer in m
	186	altLMDz(k+1)=altLMDz(k)+zdz !altitude of interface
[2690]	187	ENDDO
	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
[2699]	193	alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
[2690]	194	f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.RPI)sigma_alt_vol)* &
[2699]	195	& exp(-0.5((alt-altemiss_vol)/sigma_alt_vol)2.) &
	196	& (altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
[2690]	197	ENDDO
	198	f_lay_sum=f_lay_sum+f_lay_emiss(k)
	199	ENDDO
	200	!correct for step integration error
	201	f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum
	202	!emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
	203	!vertically distributed emission
	204	DO k=1, klev
[2787]	205	! stretch emission over one day of Pinatubo eruption
	206	emission=m_aer_emiss_vol(mSO2mol/mSatom)/m_air_gridbox(i,k)f_lay_emiss(k)/1./86400.
	207	tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys
	208	budg_emi_so2(i)=budg_emi_so2(i)+emissionzdm(k)mSatom/mSO2mol
[2690]	209	ENDDO
	210	ENDIF ! emission grid cell
	211	ENDDO ! klon loop
	212	ENDIF ! emission period
	213
	214	CASE(2) ! stratospheric aerosol injections (SAI)
	215	!
	216	DO i=1,klon
	217	! SAI standard scenario with continuous emission from 1 grid point at the equator
	218	! SAI emission on single month
	219	! IF ((mth_cur==4 .AND. &
	220	! SAI continuous emission o
[2704]	221	IF ( xlat(i).GE.xlat_sai-dlat .AND. xlat(i).LT.xlat_sai+dlat .AND. &
	222	& xlon(i).GE.xlon_sai-dlon .AND. xlon(i).LT.xlon_sai+dlon ) THEN
[2787]	223	!
	224	PRINT *,'coordinates of SAI point=',xlat(i), xlon(i), day_cur, mth_cur, year_cur
[2690]	225	! compute altLMDz
	226	altLMDz(:)=0.0
	227	DO k=1, klev
[2787]	228	zrho=pplay(i,k)/t_seri(i,k)/RD !air density in kg/m3
	229	zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG !mass of layer in kg
	230	zdz=zdm(k)/zrho !thickness of layer in m
	231	altLMDz(k+1)=altLMDz(k)+zdz !altitude of interface
[2690]	232	ENDDO
	233	!compute distribution of emission to vertical model layers (based on Gaussian peak in altitude)
	234	f_lay_sum=0.0
	235	DO k=1, klev
	236	f_lay_emiss(k)=0.0
	237	DO i_int=1, n_int_alt
[2699]	238	alt=altLMDz(k)+float(i_int)*(altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
	239	f_lay_emiss(k)=f_lay_emiss(k)+1./(sqrt(2.RPI)sigma_alt_sai)* &
	240	& exp(-0.5((alt-altemiss_sai)/sigma_alt_sai)2.) &
	241	& (altLMDz(k+1)-altLMDz(k))/float(n_int_alt)
[2690]	242	ENDDO
	243	f_lay_sum=f_lay_sum+f_lay_emiss(k)
	244	ENDDO
	245	!correct for step integration error
	246	f_lay_emiss(:)=f_lay_emiss(:)/f_lay_sum
	247	!emission as SO2 gas (with m(SO2)=64/32*m_aer_emiss)
	248	!vertically distributed emission
	249	DO k=1, klev
[2787]	250	! stretch emission over whole year (360d)
	251	emission=m_aer_emiss_sai(mSO2mol/mSatom)/m_air_gridbox(i,k)f_lay_emiss(k)/360./86400.
	252	tr_seri(i,k,id_SO2_strat)=tr_seri(i,k,id_SO2_strat)+emission*pdtphys
	253	budg_emi_so2(i)=budg_emi_so2(i)+emissionzdm(k)mSatom/mSO2mol
[2690]	254	ENDDO
	255	! !emission as monodisperse particles with 0.1um dry radius (BIN21)
	256	! !vertically distributed emission
	257	! DO k=1, klev
[2787]	258	! ! stretch emission over whole year (360d)
	259	! emission=m_aer_emiss(mH2SO4mol/mSatom)/m_part_dry(21)/m_air_gridbox(i,k)f_lay_emiss(k)/360./86400
	260	! tr_seri(i,k,id_BIN01_strat+20)=tr_seri(i,k,id_BIN01_strat+20)+emission*pdtphys
	261	! budg_emi_part(i)=budg_emi_part(i)+emissionzdm(k)mSatom/mH2SO4mol
[2690]	262	! ENDDO
	263	ENDIF ! emission grid cell
	264	ENDDO ! klon loop
	265
	266	END SELECT ! emission scenario (flag_sulf_emit)
	267
	268	!--read background concentrations of OCS and SO2 and lifetimes from input file
[2695]	269	!--update the variables defined in phys_local_var_mod
	270	CALL interp_sulf_input(debutphy,pdtphys,paprs,tr_seri)
[2690]	271
	272	!--convert OCS to SO2 in the stratosphere
[2787]	273	CALL ocs_to_so2(pdtphys,tr_seri,t_seri,pplay,paprs,is_strato)
[2690]	274
	275	!--convert SO2 to H2SO4
[2787]	276	CALL so2_to_h2so4(pdtphys,tr_seri,t_seri,pplay,paprs,is_strato)
[2690]	277
	278	!--common routine for nucleation and condensation/evaporation with adaptive timestep
	279	CALL micphy_tstep(pdtphys,tr_seri,t_seri,pplay,paprs,rh,is_strato)
	280
	281	!--call coagulation routine
	282	CALL coagulate(pdtphys,mdw,tr_seri,t_seri,pplay,dens_aer,is_strato)
	283
	284	!--call sedimentation routine
	285	CALL aer_sedimnt(pdtphys, t_seri, pplay, paprs, tr_seri, dens_aer)
	286
	287	!--compute mass concentration of PM2.5 sulfate particles (wet diameter and mass) at the surface for health studies
	288	surf_PM25_sulf(:)=0.0
	289	DO i=1,klon
	290	DO it=1, nbtr_bin
	291	IF (mdw(it) .LT. 2.5e-6) THEN
	292	!surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas)*m_part(i,1,it) &
	293	!assume that particles consist of ammonium sulfate at the surface (132g/mol) and are dry at T = 20 deg. C and 50 perc. humidity
	294	surf_PM25_sulf(i)=surf_PM25_sulf(i)+tr_seri(i,1,it+nbtr_sulgas) &
	295	& 132./98.dens_aer_dry4./3.RPI(mdw(it)/2.)*3 &
[2787]	296	& pplay(i,1)/t_seri(i,1)/RD1.e9
[2690]	297	ENDIF
	298	ENDDO
	299	ENDDO
	300
[2695]	301	! CALL minmaxsimple(tr_seri(:,:,id_SO2_strat),0.0,0.0,'fin SO2')
	302	! DO it=1, nbtr_bin
	303	! CALL minmaxsimple(tr_seri(:,:,nbtr_sulgas+it),0.0,0.0,'fin bin ')
	304	! ENDDO
	305
[2690]	306	END SUBROUTINE traccoag
	307
	308	END MODULE traccoag_mod

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

Download in other formats:

Original Format