Context Navigation

coagulate.F90 @ 3526

Last change on this file since 3526 was 3526, checked in by Laurent Fairhead, 5 years ago
Olivier's modifications + addition of Id keyword
Property svn:keywords set to `Id`
File size: 9.2 KB

Line
1	!
2	! $Id: coagulate.F90 3526 2019-05-28 13:00:44Z fairhead $
3	!
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
27	USE infotrac
28	USE phys_local_var_mod, ONLY: DENSO4, f_r_wet
29	USE YOMCST
30
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
46	REAL, DIMENSION(nbtr_bin) :: radius ! aerosol particle radius in each bin [m]
47	REAL, DIMENSION(klon,klev,nbtr_bin) :: tr_t ! Concentration Traceur at time t [U/KgA]
48	REAL, DIMENSION(klon,klev,nbtr_bin) :: tr_tp1 ! Concentration Traceur at time t+1 [U/KgA]
49	REAL, DIMENSION(nbtr_bin,nbtr_bin,nbtr_bin) :: ff ! Volume fraction of intermediate particles
50	REAL, DIMENSION(nbtr_bin) :: V ! Volume of bins
51	REAL, DIMENSION(nbtr_bin,nbtr_bin) :: Vij ! Volume sum of i and j
52	REAL :: eta ! Dynamic viscosity of air
53	REAL, PARAMETER :: mair=4.8097E-26 ! Average mass of an air molecule [kg]
54	REAL :: zrho ! Density of air
55	REAL :: mnfrpth ! Mean free path of air
56	REAL, DIMENSION(nbtr_bin) :: Kn ! Knudsen number of particle i
57	REAL, DIMENSION(nbtr_bin) :: Di ! Particle diffusion coefficient
58	REAL, DIMENSION(nbtr_bin) :: m_par ! Mass of particle i
59	REAL, DIMENSION(nbtr_bin) :: thvelpar! Thermal velocity of particle i
60	REAL, DIMENSION(nbtr_bin) :: mfppar ! Mean free path of particle i
61	REAL, DIMENSION(nbtr_bin) :: delta! delta of particle i (from equation 21)
62	REAL, DIMENSION(nbtr_bin,nbtr_bin) :: beta ! Coagulation kernel from Brownian diffusion
63	REAL :: beta_const ! Constant coagulation kernel (for comparison)
64	REAL :: num
65	REAL :: numi
66	REAL :: denom
67
68	! Additional variables for coagulation enhancement factor due to van der Waals forces
69	! Taken from Chan and Mozurkewich, Measurement of the coagulation rate constant for sulfuric acid
70	! particles as a function of particle size using TDMA, Aerosol Science, 32, 321-339, 2001.
71	!--ok_vdw is 0 for no vdW forces, 1 for E(0), 2 for E(infinity)
72	INTEGER, PARAMETER :: ok_vdw = 0
73	REAL, PARAMETER :: avdW1 = 0.0757
74	REAL, PARAMETER :: avdW3 = 0.0015
75	REAL, PARAMETER :: bvdW0 = 0.0151
76	REAL, PARAMETER :: bvdW1 = -0.186
77	REAL, PARAMETER :: bvdW3 = -0.0163
78	REAL, PARAMETER :: AvdW = 6.4e-20 !Hamaker constant in J = 1e7 erg
79	REAL :: AvdWi
80	REAL :: xvdW
81	REAL :: EvdW
82
83	DO i=1, nbtr_bin
84	radius(i)=mdw(i)/2.
85	V(i)= radius(i)*3. !neglecting factor 4RPI/3
86	ENDDO
87
88	DO j=1, nbtr_bin
89	DO i=1, nbtr_bin
90	Vij(i,j)= V(i)+V(j)
91	ENDDO
92	ENDDO
93
94	!--pre-compute the f(i,j,k) from Jacobson equation 13
95	ff=0.0
96	DO k=1, nbtr_bin
97	DO j=1, nbtr_bin
98	DO i=1, nbtr_bin
99	IF (k.EQ.1) THEN
100	ff(i,j,k)= 0.0
101	ELSEIF (k.GT.1.AND.V(k-1).LT.Vij(i,j).AND.Vij(i,j).LT.V(k)) THEN
102	ff(i,j,k)= 1.-ff(i,j,k-1)
103	ELSEIF (k.EQ.nbtr_bin) THEN
104	IF (Vij(i,j).GE.v(k)) THEN
105	ff(i,j,k)= 1.
106	ELSE
107	ff(i,j,k)= 0.0
108	ENDIF
109	ELSEIF (k.LE.(nbtr_bin-1).AND.V(k).LE.Vij(i,j).AND.Vij(i,j).LT.V(k+1)) THEN
110	ff(i,j,k)= V(k)/Vij(i,j)*(V(k+1)-Vij(i,j))/(V(k+1)-V(k))
111	ENDIF
112	ENDDO
113	ENDDO
114	ENDDO
115
116	DO ilon=1, klon
117	DO ilev=1, klev
118	!only in the stratosphere
119	IF (is_strato(ilon,ilev)) THEN
120	!compute actual wet particle radius & volume for every grid box
121	DO i=1, nbtr_bin
122	radius(i)=f_r_wet(ilon,ilev)*mdw(i)/2.
123	V(i)= radius(i)*3. !neglecting factor 4RPI/3
124	ENDDO
125
126	!--Calculations for the coagulation kernel---------------------------------------------------------
127
128	zrho=pplay(ilon,ilev)/t_seri(ilon,ilev)/RD
129
130	!--initialize the tracer at time t and convert from [number/KgA] to [number/m3]
131	DO i=1, nbtr_bin
132	tr_t(ilon,ilev,i) = tr_seri(ilon,ilev,i+nbtr_sulgas) * zrho
133	ENDDO
134
135	! mean free path of air (Pruppacher and Klett, 2010, p.417) [m]
136	mnfrpth=6.6E-8(1.01325E+5/pplay(ilon,ilev))(t_seri(ilon,ilev)/293.15)
137	! mnfrpth=2.eta/(zrhothvelair)
138	! mean free path of air (Prupp. Klett) in [10^-6 m]
139	! ZLAIR = 0.066 (1.01325E+5/PPLAY)(T_SERI/293.15)*1.E-06
140
141	! dynamic viscosity of air (Pruppacher and Klett, 2010, p.417) [kg/(m*s)]
142	IF (t_seri(ilon,ilev).GE.273.15) THEN
143	eta=(1.718+0.0049(t_seri(ilon,ilev)-273.15))1.E-5
144	ELSE
145	eta=(1.718+0.0049(t_seri(ilon,ilev)-273.15)-1.2E-5(t_seri(ilon,ilev)-273.15)*2)1.E-5
146	ENDIF
147
148	!--pre-compute the particle diffusion coefficient Di(i) from equation 18
149	Di=0.0
150	DO i=1, nbtr_bin
151	Kn(i)=mnfrpth/radius(i)
152	Di(i)=RKBOLt_seri(ilon,ilev)/(6.RPIradius(i)eta)(1.+Kn(i)(1.249+0.42*exp(-0.87/Kn(i))))
153	ENDDO
154
155	!--pre-compute the thermal velocity of a particle thvelpar(i) from equation 20
156	thvelpar=0.0
157	DO i=1, nbtr_bin
158	m_par(i)=4./3.RPIradius(i)*3.DENSO4(ilon,ilev)*1000.
159	thvelpar(i)=sqrt(8.RKBOLt_seri(ilon,ilev)/(RPI*m_par(i)))
160	ENDDO
161
162	!--pre-compute the particle mean free path mfppar(i) from equation 22
163	mfppar=0.0
164	DO i=1, nbtr_bin
165	mfppar(i)=8.Di(i)/(RPIthvelpar(i))
166	ENDDO
167
168	!--pre-compute the mean distance delta(i) from the center of a sphere reached by particles
169	!--leaving the surface of the sphere and traveling a distance of particle mfppar(i) from equation 21
170	delta=0.0
171	DO i=1, nbtr_bin
172	delta(i)=((2.radius(i)+mfppar(i))3.-(4.radius(i)2.+mfppar(i)2.)**1.5)/ &
173	& (6.radius(i)mfppar(i))-2.*radius(i)
174	ENDDO
175
176	!--pre-compute the beta(i,j) from equation 17 in Jacobson
177	num=0.0
178	DO j=1, nbtr_bin
179	DO i=1, nbtr_bin
180	!
181	num=4.RPI(radius(i)+radius(j))*(Di(i)+Di(j))
182	denom=(radius(i)+radius(j))/(radius(i)+radius(j)+sqrt(delta(i)2.+delta(j)2.))+ &
183	& 4.(Di(i)+Di(j))/(sqrt(thvelpar(i)2.+thvelpar(j)2.)(radius(i)+radius(j)))
184	beta(i,j)=num/denom
185	!
186	!--compute enhancement factor due to van der Waals forces
187	IF (ok_vdw .EQ. 0) THEN !--no enhancement factor
188	Evdw=1.0
189	ELSEIF (ok_vdw .EQ. 1) THEN !--E(0) case
190	AvdWi = AvdW/(RKBOLt_seri(ilon,ilev))(4.radius(i)radius(j))/(radius(i)+radius(j))**2.
191	xvdW = LOG(1.+AvdWi)
192	EvdW = 1. + avdW1xvdW + avdW3xvdW**3
193	ELSEIF (ok_vdw .EQ. 2) THEN !--E(infinity) case
194	AvdWi = AvdW/(RKBOLt_seri(ilon,ilev))(4.radius(i)radius(j))/(radius(i)+radius(j))**2.
195	xvdW = LOG(1.+AvdWi)
196	EvdW = 1. + SQRT(AvdWi/3.)/(1.+bvdW0SQRT(AvdWi)) + bvdW1xvdW + bvdW3xvdW*3.
197	ENDIF
198	!
199	beta(i,j)=beta(i,j)*EvdW
200
201	ENDDO
202	ENDDO
203
204	!--external loop for equation 14
205	DO k=1, nbtr_bin
206
207	!--calculating denominator sum
208	denom=0.0
209	DO j=1, nbtr_bin
210	denom=denom+(1.-ff(k,j,k))beta(k,j)tr_t(ilon,ilev,j)
211	ENDDO
212
213	IF (k.EQ.1) THEN
214	!--calculate new concentration of smallest bin
215	tr_tp1(ilon,ilev,k)=tr_t(ilon,ilev,k)/(1.+pdtcoag*denom)
216	ELSE
217	!--calculating double sum terms in numerator of eq 14
218	num=0.0
219	DO j=1, k
220	numi=0.0
221	DO i=1, k-1
222	numi=numi+ff(i,j,k)beta(i,j)V(i)tr_tp1(ilon,ilev,i)tr_t(ilon,ilev,j)
223	ENDDO
224	num=num+numi
225	ENDDO
226
227	!--calculate new concentration of other bins
228	tr_tp1(ilon,ilev,k)=(V(k)tr_t(ilon,ilev,k)+pdtcoagnum)/(1.+pdtcoag*denom)/V(k)
229	ENDIF
230
231	ENDDO !--end of loop k
232
233	!--convert tracer concentration back from [number/m3] to [number/KgA] and write into tr_seri
234	DO i=1, nbtr_bin
235	tr_seri(ilon,ilev,i+nbtr_sulgas) = tr_tp1(ilon,ilev,i) / zrho
236	ENDDO
237
238	ENDIF ! IF IN STRATOSPHERE
239	ENDDO !--end of loop klev
240	ENDDO !--end of loop klon
241
242	END SUBROUTINE COAGULATE

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

Download in other formats:

Original Format