Context Navigation

coagulate.F90 @ 4766

Last change on this file since 4766 was 4762, checked in by lguez, 7 months ago

Use YOMCST from LMDZ, not from RRTM or ECRad

Because the modules yomcst in RRTM and ECRad do not contain the same
objects. In particular, RKBOL is in rrtm/yomcst.F90 but not in
ecrad/yomcst.F90. So compilation with -rad ecrad -strataer true
fails. And we want to avoid modifying files from ECRad.

Property svn:keywords set to Id

File size: 9.2 KB

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

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

Download in other formats:

Original Format