Context Navigation

coagulate.f90 @ 5274

Last change on this file since 5274 was 5274, checked in by abarral, 31 hours ago
Replace yomcst.h by existing module
Property svn:keywords set to `Id`
File size: 12.8 KB

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

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

Download in other formats:

Original Format