Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

coagulate.F90

Last change on this file was 4950, checked in by lebasn, 4 months ago
StratAer?: New model version (microphysic, composition routine, code cleaning, new params...)
Property svn:keywords set to `Id`
File size: 12.0 KB

Rev	Line
[3526]	1	!
	2	! $Id: coagulate.F90 4950 2024-05-22 13:16:36Z evignon $
	3	!
[2690]	4	SUBROUTINE COAGULATE(pdtcoag,mdw,tr_seri,t_seri,pplay,dens_aer,is_strato)
	5	! -----------------------------------------------------------------------
	6	!
	7	! Author : Christoph Kleinschmitt (with Olivier Boucher)
	8	! ------
	9	!
	10	! purpose
	11	! -------
	12	!
	13	! interface
	14	! ---------
	15	! input
	16	! pdtphys time step duration [sec]
	17	! tr_seri tracer mixing ratios [kg/kg]
	18	! mdw # or mass median diameter [m]
	19	!
	20	! method
	21	! ------
	22	!
	23	! -----------------------------------------------------------------------
	24
	25	USE dimphy, ONLY : klon,klev
	26	USE aerophys
[3677]	27	USE infotrac_phy
[4950]	28	USE phys_local_var_mod, ONLY: DENSO4, DENSO4B, f_r_wet, f_r_wetB
	29	USE strataer_local_var_mod, ONLY: flag_new_strat_compo
	30
[2690]	31	IMPLICIT NONE
	32
	33	!--------------------------------------------------------
	34
	35	! transfer variables when calling this routine
	36	REAL,INTENT(IN) :: pdtcoag ! Time step in coagulation routine [s]
	37	REAL,DIMENSION(nbtr_bin),INTENT(IN) :: mdw ! aerosol particle diameter in each bin [m]
	38	REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA]
	39	REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature
	40	REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa)
	41	REAL,DIMENSION(klon,klev) :: dens_aer! density of aerosol [kg/m3 aerosol] with default H2SO4 mass
	42	LOGICAL,DIMENSION(klon,klev),INTENT(IN) :: is_strato
	43
	44	! local variables in coagulation routine
	45	INTEGER :: i,j,k,nb,ilon,ilev
[4950]	46	REAL, DIMENSION(nbtr_bin) :: radiusdry ! dry aerosol particle radius in each bin [m]
	47	REAL, DIMENSION(nbtr_bin) :: radiuswet ! wet aerosol particle radius in each bin [m]
[2690]	48	REAL, DIMENSION(klon,klev,nbtr_bin) :: tr_t ! Concentration Traceur at time t [U/KgA]
	49	REAL, DIMENSION(klon,klev,nbtr_bin) :: tr_tp1 ! Concentration Traceur at time t+1 [U/KgA]
	50	REAL, DIMENSION(nbtr_bin,nbtr_bin,nbtr_bin) :: ff ! Volume fraction of intermediate particles
[4950]	51	REAL, DIMENSION(nbtr_bin) :: Vdry ! Volume dry of bins
	52	REAL, DIMENSION(nbtr_bin) :: Vwet ! Volume wet of bins
[2690]	53	REAL, DIMENSION(nbtr_bin,nbtr_bin) :: Vij ! Volume sum of i and j
	54	REAL :: eta ! Dynamic viscosity of air
	55	REAL, PARAMETER :: mair=4.8097E-26 ! Average mass of an air molecule [kg]
	56	REAL :: zrho ! Density of air
	57	REAL :: mnfrpth ! Mean free path of air
	58	REAL, DIMENSION(nbtr_bin) :: Kn ! Knudsen number of particle i
	59	REAL, DIMENSION(nbtr_bin) :: Di ! Particle diffusion coefficient
	60	REAL, DIMENSION(nbtr_bin) :: m_par ! Mass of particle i
	61	REAL, DIMENSION(nbtr_bin) :: thvelpar! Thermal velocity of particle i
	62	REAL, DIMENSION(nbtr_bin) :: mfppar ! Mean free path of particle i
	63	REAL, DIMENSION(nbtr_bin) :: delta! delta of particle i (from equation 21)
	64	REAL, DIMENSION(nbtr_bin,nbtr_bin) :: beta ! Coagulation kernel from Brownian diffusion
	65	REAL :: beta_const ! Constant coagulation kernel (for comparison)
	66	REAL :: num
	67	REAL :: numi
	68	REAL :: denom
	69
[2949]	70	! Additional variables for coagulation enhancement factor due to van der Waals forces
	71	! Taken from Chan and Mozurkewich, Measurement of the coagulation rate constant for sulfuric acid
	72	! particles as a function of particle size using TDMA, Aerosol Science, 32, 321-339, 2001.
	73	!--ok_vdw is 0 for no vdW forces, 1 for E(0), 2 for E(infinity)
	74	INTEGER, PARAMETER :: ok_vdw = 0
	75	REAL, PARAMETER :: avdW1 = 0.0757
	76	REAL, PARAMETER :: avdW3 = 0.0015
	77	REAL, PARAMETER :: bvdW0 = 0.0151
	78	REAL, PARAMETER :: bvdW1 = -0.186
	79	REAL, PARAMETER :: bvdW3 = -0.0163
	80	REAL, PARAMETER :: AvdW = 6.4e-20 !Hamaker constant in J = 1e7 erg
	81	REAL :: AvdWi
	82	REAL :: xvdW
	83	REAL :: EvdW
	84
[4762]	85	include "YOMCST.h"
	86
[4950]	87	! ff(i,j,k): Volume fraction of Vi,j that is partitioned to each model bin k
	88	! just need to be calculated in model initialization because mdw(:) size is fixed
	89	! no need to recalculate radius, Vdry, Vij, and ff every timestep because it is for
	90	! dry aerosols
[2690]	91	DO i=1, nbtr_bin
[4950]	92	radiusdry(i)=mdw(i)/2.
	93	Vdry(i)=radiusdry(i)*3. !neglecting factor 4RPI/3
	94	Vwet(i)=0.0
[2690]	95	ENDDO
	96
	97	DO j=1, nbtr_bin
[4950]	98	DO i=1, nbtr_bin
	99	Vij(i,j)= Vdry(i)+Vdry(j)
	100	ENDDO
[2690]	101	ENDDO
[4950]	102
[2690]	103	!--pre-compute the f(i,j,k) from Jacobson equation 13
	104	ff=0.0
	105	DO k=1, nbtr_bin
	106	DO j=1, nbtr_bin
	107	DO i=1, nbtr_bin
	108	IF (k.EQ.1) THEN
	109	ff(i,j,k)= 0.0
[4950]	110	ELSEIF (k.GT.1.AND.Vdry(k-1).LT.Vij(i,j).AND.Vij(i,j).LT.Vdry(k)) THEN
[2690]	111	ff(i,j,k)= 1.-ff(i,j,k-1)
	112	ELSEIF (k.EQ.nbtr_bin) THEN
[4950]	113	IF (Vij(i,j).GE.Vdry(k)) THEN
[2690]	114	ff(i,j,k)= 1.
	115	ELSE
	116	ff(i,j,k)= 0.0
	117	ENDIF
[4950]	118	ELSEIF (k.LE.(nbtr_bin-1).AND.Vdry(k).LE.Vij(i,j).AND.Vij(i,j).LT.Vdry(k+1)) THEN
	119	ff(i,j,k)= Vdry(k)/Vij(i,j)*(Vdry(k+1)-Vij(i,j))/(Vdry(k+1)-Vdry(k))
[2690]	120	ENDIF
	121	ENDDO
	122	ENDDO
	123	ENDDO
[4950]	124	! End of just need to be calculated at initialization because mdw(:) size is fixed
	125
[2690]	126	DO ilon=1, klon
	127	DO ilev=1, klev
	128	!only in the stratosphere
	129	IF (is_strato(ilon,ilev)) THEN
	130	!compute actual wet particle radius & volume for every grid box
[4950]	131	IF(flag_new_strat_compo) THEN
	132	DO i=1, nbtr_bin
	133	radiuswet(i)=f_r_wetB(ilon,ilev,i)*mdw(i)/2.
	134	Vwet(i)= radiuswet(i)*3. !neglecting factor 4RPI/3
	135	!! Vwet(i)= Vdry(i)(f_r_wetB(ilon,ilev,i)*3)
	136	ENDDO
	137	ELSE
	138	DO i=1, nbtr_bin
	139	radiuswet(i)=f_r_wet(ilon,ilev)*mdw(i)/2.
	140	Vwet(i)= radiuswet(i)*3. !neglecting factor 4RPI/3
	141	!! Vwet(i)= Vdry(i)(f_r_wet(ilon,ilev)*3)
	142	ENDDO
	143	ENDIF
	144
[2690]	145	!--Calculations for the coagulation kernel---------------------------------------------------------
	146
	147	zrho=pplay(ilon,ilev)/t_seri(ilon,ilev)/RD
	148
	149	!--initialize the tracer at time t and convert from [number/KgA] to [number/m3]
	150	DO i=1, nbtr_bin
	151	tr_t(ilon,ilev,i) = tr_seri(ilon,ilev,i+nbtr_sulgas) * zrho
	152	ENDDO
	153
	154	! mean free path of air (Pruppacher and Klett, 2010, p.417) [m]
	155	mnfrpth=6.6E-8(1.01325E+5/pplay(ilon,ilev))(t_seri(ilon,ilev)/293.15)
	156	! mnfrpth=2.eta/(zrhothvelair)
	157	! mean free path of air (Prupp. Klett) in [10^-6 m]
	158	! ZLAIR = 0.066 (1.01325E+5/PPLAY)(T_SERI/293.15)*1.E-06
	159
	160	! dynamic viscosity of air (Pruppacher and Klett, 2010, p.417) [kg/(m*s)]
	161	IF (t_seri(ilon,ilev).GE.273.15) THEN
	162	eta=(1.718+0.0049(t_seri(ilon,ilev)-273.15))1.E-5
	163	ELSE
	164	eta=(1.718+0.0049(t_seri(ilon,ilev)-273.15)-1.2E-5(t_seri(ilon,ilev)-273.15)*2)1.E-5
	165	ENDIF
	166
	167	!--pre-compute the particle diffusion coefficient Di(i) from equation 18
	168	Di=0.0
	169	DO i=1, nbtr_bin
[4950]	170	Kn(i)=mnfrpth/radiuswet(i)
	171	Di(i)=RKBOLt_seri(ilon,ilev)/(6.RPIradiuswet(i)eta)(1.+Kn(i)(1.249+0.42*exp(-0.87/Kn(i))))
[2690]	172	ENDDO
	173
	174	!--pre-compute the thermal velocity of a particle thvelpar(i) from equation 20
	175	thvelpar=0.0
[4950]	176	IF(flag_new_strat_compo) THEN
	177	DO i=1, nbtr_bin
	178	m_par(i)=4./3.RPIradiuswet(i)*3.DENSO4B(ilon,ilev,i)*1000.
	179	thvelpar(i)=sqrt(8.RKBOLt_seri(ilon,ilev)/(RPI*m_par(i)))
	180	ENDDO
	181	ELSE
	182	DO i=1, nbtr_bin
	183	m_par(i)=4./3.RPIradiuswet(i)*3.DENSO4(ilon,ilev)*1000.
	184	thvelpar(i)=sqrt(8.RKBOLt_seri(ilon,ilev)/(RPI*m_par(i)))
	185	ENDDO
	186	ENDIF
[2690]	187
	188	!--pre-compute the particle mean free path mfppar(i) from equation 22
	189	mfppar=0.0
	190	DO i=1, nbtr_bin
	191	mfppar(i)=8.Di(i)/(RPIthvelpar(i))
	192	ENDDO
	193
	194	!--pre-compute the mean distance delta(i) from the center of a sphere reached by particles
	195	!--leaving the surface of the sphere and traveling a distance of particle mfppar(i) from equation 21
	196	delta=0.0
	197	DO i=1, nbtr_bin
[4950]	198	delta(i)=((2.radiuswet(i)+mfppar(i))3.-(4.radiuswet(i)2.+mfppar(i)2.)**1.5)/ &
	199	& (6.radiuswet(i)mfppar(i))-2.*radiuswet(i)
[2690]	200	ENDDO
	201
[4950]	202	! beta(i,j): coagulation kernel (rate coefficient) of 2 colliding particles i,j
[2949]	203	!--pre-compute the beta(i,j) from equation 17 in Jacobson
[2690]	204	num=0.0
	205	DO j=1, nbtr_bin
	206	DO i=1, nbtr_bin
[2949]	207	!
[4950]	208	num=4.RPI(radiuswet(i)+radiuswet(j))*(Di(i)+Di(j))
	209	denom=(radiuswet(i)+radiuswet(j))/(radiuswet(i)+radiuswet(j)+sqrt(delta(i)2.+delta(j)2.))+ &
	210	& 4.(Di(i)+Di(j))/(sqrt(thvelpar(i)2.+thvelpar(j)2.)(radiuswet(i)+radiuswet(j)))
	211	beta(i,j)=num/denom
[2949]	212	!
	213	!--compute enhancement factor due to van der Waals forces
	214	IF (ok_vdw .EQ. 0) THEN !--no enhancement factor
[4950]	215	Evdw=1.0
[2949]	216	ELSEIF (ok_vdw .EQ. 1) THEN !--E(0) case
[4950]	217	AvdWi = AvdW/(RKBOLt_seri(ilon,ilev))(4.radiuswet(i)radiuswet(j))/(radiuswet(i)+radiuswet(j))**2.
	218	xvdW = LOG(1.+AvdWi)
	219	EvdW = 1. + avdW1xvdW + avdW3xvdW**3
[2949]	220	ELSEIF (ok_vdw .EQ. 2) THEN !--E(infinity) case
[4950]	221	AvdWi = AvdW/(RKBOLt_seri(ilon,ilev))(4.radiuswet(i)radiuswet(j))/(radiuswet(i)+radiuswet(j))**2.
	222	xvdW = LOG(1.+AvdWi)
	223	EvdW = 1. + SQRT(AvdWi/3.)/(1.+bvdW0SQRT(AvdWi)) + bvdW1xvdW + bvdW3xvdW*3.
[2949]	224	ENDIF
	225	!
	226	beta(i,j)=beta(i,j)*EvdW
	227
[2690]	228	ENDDO
	229	ENDDO
	230
	231	!--external loop for equation 14
	232	DO k=1, nbtr_bin
	233
	234	!--calculating denominator sum
	235	denom=0.0
	236	DO j=1, nbtr_bin
[4950]	237	! fraction of coagulation of k and j that is not giving k
	238	denom=denom+(1.-ff(k,j,k))beta(k,j)tr_t(ilon,ilev,j)
[2690]	239	ENDDO
	240
	241	IF (k.EQ.1) THEN
	242	!--calculate new concentration of smallest bin
	243	tr_tp1(ilon,ilev,k)=tr_t(ilon,ilev,k)/(1.+pdtcoag*denom)
	244	ELSE
	245	!--calculating double sum terms in numerator of eq 14
	246	num=0.0
	247	DO j=1, k
[4950]	248	numi=0.0
	249	DO i=1, k-1
	250	!
	251	! see Jacobson: " In order to conserve volume and volume concentration (which
	252	! coagulation physically does) while giving up some accuracy in number concentration"
	253	!
	254	! Coagulation of i and j giving k
	255	! with V(i) and then V(j) because it considers i,j and j,i with the double loop
	256	!
	257	! BUT WHY WET VOLUME V(i) in old STRATAER? tracers are already dry aerosols and coagulation
	258	! kernel beta(i,j) accounts for wet aerosols -> reply below
	259	!
	260	! numi=numi+ff(i,j,k)beta(i,j)V(i)tr_tp1(ilon,ilev,i)tr_t(ilon,ilev,j)
	261	numi=numi+ff(i,j,k)beta(i,j)Vdry(i)tr_tp1(ilon,ilev,i)tr_t(ilon,ilev,j)
	262	ENDDO
	263	num=num+numi
[2690]	264	ENDDO
	265
	266	!--calculate new concentration of other bins
[4950]	267	! tr_tp1(ilon,ilev,k)=(V(k)tr_t(ilon,ilev,k)+pdtcoagnum)/( (1.+pdtcoagdenom)V(k) )
	268	tr_tp1(ilon,ilev,k)=(Vdry(k)tr_t(ilon,ilev,k)+pdtcoagnum)/( (1.+pdtcoagdenom)Vdry(k) )
	269	!
	270	! In constant composition (no dependency on aerosol size because no kelvin effect)
	271	! V(l)= (f_r_wet(ilon,ilev)*3)((mdw(l)/2.)3) = (f_r_wet(ilon,ilev)3)*Vdry(i)
	272	! so numi and num are proportional (f_r_wet(ilon,ilev)**3)
	273	! and so
	274	! tr_tp1(ilon,ilev,k)=(V(k)tr_t(ilon,ilev,k)+pdtcoagnum)/( (1.+pdtcoagdenom)V(k) )
	275	! =(Vdry(k)tr_t(ilon,ilev,k)+pdtcoagnum_dry)/( (1.+pdtcoagdenom)Vdry(k) )
	276	! with num_dry=...beta(i,j)Vdry(i)....
	277	! so in old STRATAER (.not.flag_new_strat_compo), it was correct
[2690]	278	ENDIF
	279
	280	ENDDO !--end of loop k
	281
	282	!--convert tracer concentration back from [number/m3] to [number/KgA] and write into tr_seri
	283	DO i=1, nbtr_bin
[4950]	284	tr_seri(ilon,ilev,i+nbtr_sulgas) = tr_tp1(ilon,ilev,i) / zrho
[2690]	285	ENDDO
	286
	287	ENDIF ! IF IN STRATOSPHERE
	288	ENDDO !--end of loop klev
	289	ENDDO !--end of loop klon
[4950]	290	! *********************************************
[2690]	291
	292	END SUBROUTINE COAGULATE

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

Download in other formats:

Original Format