1 | ! |
---|
2 | ! $Id: coagulate.F90 3677 2020-05-06 15:18:32Z nfevrier $ |
---|
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 | 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 4*RPI/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 4*RPI/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/(zrho*thvelair) |
---|
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)=RKBOL*t_seri(ilon,ilev)/(6.*RPI*radius(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.*RPI*radius(i)**3.*DENSO4(ilon,ilev)*1000. |
---|
159 | thvelpar(i)=sqrt(8.*RKBOL*t_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)/(RPI*thvelpar(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/(RKBOL*t_seri(ilon,ilev))*(4.*radius(i)*radius(j))/(radius(i)+radius(j))**2. |
---|
191 | xvdW = LOG(1.+AvdWi) |
---|
192 | EvdW = 1. + avdW1*xvdW + avdW3*xvdW**3 |
---|
193 | ELSEIF (ok_vdw .EQ. 2) THEN !--E(infinity) case |
---|
194 | AvdWi = AvdW/(RKBOL*t_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.+bvdW0*SQRT(AvdWi)) + bvdW1*xvdW + bvdW3*xvdW**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)+pdtcoag*num)/(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 |
---|