Context Navigation

coagulate.F90 @ 2690

Last change on this file since 2690 was 2690, checked in by oboucher, 8 years ago
Adding a module for stratospheric aerosols with a bin scheme. The module gets activated with -strataer true compiling option. May not quite work yet, more testing needed, but should not affect the rest of LMDz as everything is under a CPP_StratAer key.
File size: 7.6 KB

Line
1	SUBROUTINE COAGULATE(pdtcoag,mdw,tr_seri,t_seri,pplay,dens_aer,is_strato)
2	! -----------------------------------------------------------------------
3	!
4	! Author : Christoph Kleinschmitt (with Olivier Boucher)
5	! ------
6	!
7	! purpose
8	! -------
9	!
10	! interface
11	! ---------
12	! input
13	! pdtphys time step duration [sec]
14	! tr_seri tracer mixing ratios [kg/kg]
15	! mdw # or mass median diameter [m]
16	!
17	! method
18	! ------
19	!
20	! -----------------------------------------------------------------------
21
22	USE dimphy, ONLY : klon,klev
23	USE aerophys
24	USE infotrac
25	USE phys_local_var_mod, ONLY: DENSO4, f_r_wet
26	USE YOMCST
27
28	IMPLICIT NONE
29
30	!--------------------------------------------------------
31
32	! transfer variables when calling this routine
33	REAL,INTENT(IN) :: pdtcoag ! Time step in coagulation routine [s]
34	REAL,DIMENSION(nbtr_bin),INTENT(IN) :: mdw ! aerosol particle diameter in each bin [m]
35	REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA]
36	REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature
37	REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa)
38	REAL,DIMENSION(klon,klev) :: dens_aer! density of aerosol [kg/m3 aerosol] with default H2SO4 mass
39	LOGICAL,DIMENSION(klon,klev),INTENT(IN) :: is_strato
40
41	! local variables in coagulation routine
42	INTEGER :: i,j,k,nb,ilon,ilev
43	REAL, DIMENSION(nbtr_bin) :: radius ! aerosol particle radius in each bin [m]
44	REAL, DIMENSION(klon,klev,nbtr_bin) :: tr_t ! Concentration Traceur at time t [U/KgA]
45	REAL, DIMENSION(klon,klev,nbtr_bin) :: tr_tp1 ! Concentration Traceur at time t+1 [U/KgA]
46	REAL, DIMENSION(nbtr_bin,nbtr_bin,nbtr_bin) :: ff ! Volume fraction of intermediate particles
47	REAL, DIMENSION(nbtr_bin) :: V ! Volume of bins
48	REAL, DIMENSION(nbtr_bin,nbtr_bin) :: Vij ! Volume sum of i and j
49	REAL :: eta ! Dynamic viscosity of air
50	REAL, PARAMETER :: mair=4.8097E-26 ! Average mass of an air molecule [kg]
51	REAL :: zrho ! Density of air
52	REAL :: mnfrpth ! Mean free path of air
53	REAL, DIMENSION(nbtr_bin) :: Kn ! Knudsen number of particle i
54	REAL, DIMENSION(nbtr_bin) :: Di ! Particle diffusion coefficient
55	REAL, DIMENSION(nbtr_bin) :: m_par ! Mass of particle i
56	REAL, DIMENSION(nbtr_bin) :: thvelpar! Thermal velocity of particle i
57	REAL, DIMENSION(nbtr_bin) :: mfppar ! Mean free path of particle i
58	REAL, DIMENSION(nbtr_bin) :: delta! delta of particle i (from equation 21)
59	REAL, DIMENSION(nbtr_bin,nbtr_bin) :: beta ! Coagulation kernel from Brownian diffusion
60	REAL :: beta_const ! Constant coagulation kernel (for comparison)
61	REAL :: num
62	REAL :: numi
63	REAL :: denom
64
65	DO i=1, nbtr_bin
66	radius(i)=mdw(i)/2.
67	V(i)= radius(i)*3. !neglecting factor 4RPI/3
68	ENDDO
69
70	DO j=1, nbtr_bin
71	DO i=1, nbtr_bin
72	Vij(i,j)= V(i)+V(j)
73	ENDDO
74	ENDDO
75
76	!--pre-compute the f(i,j,k) from Jacobson equation 13
77	ff=0.0
78	DO k=1, nbtr_bin
79	DO j=1, nbtr_bin
80	DO i=1, nbtr_bin
81	IF (k.EQ.1) THEN
82	ff(i,j,k)= 0.0
83	ELSEIF (k.GT.1.AND.V(k-1).LT.Vij(i,j).AND.Vij(i,j).LT.V(k)) THEN
84	ff(i,j,k)= 1.-ff(i,j,k-1)
85	ELSEIF (k.EQ.nbtr_bin) THEN
86	IF (Vij(i,j).GE.v(k)) THEN
87	ff(i,j,k)= 1.
88	ELSE
89	ff(i,j,k)= 0.0
90	ENDIF
91	ELSEIF (k.LE.(nbtr_bin-1).AND.V(k).LE.Vij(i,j).AND.Vij(i,j).LT.V(k+1)) THEN
92	ff(i,j,k)= V(k)/Vij(i,j)*(V(k+1)-Vij(i,j))/(V(k+1)-V(k))
93	ENDIF
94	ENDDO
95	ENDDO
96	ENDDO
97
98	DO ilon=1, klon
99	DO ilev=1, klev
100	!only in the stratosphere
101	IF (is_strato(ilon,ilev)) THEN
102	!compute actual wet particle radius & volume for every grid box
103	DO i=1, nbtr_bin
104	radius(i)=f_r_wet(ilon,ilev)*mdw(i)/2.
105	V(i)= radius(i)*3. !neglecting factor 4RPI/3
106	ENDDO
107
108	!--Calculations for the coagulation kernel---------------------------------------------------------
109
110	zrho=pplay(ilon,ilev)/t_seri(ilon,ilev)/RD
111
112	!--initialize the tracer at time t and convert from [number/KgA] to [number/m3]
113	DO i=1, nbtr_bin
114	tr_t(ilon,ilev,i) = tr_seri(ilon,ilev,i+nbtr_sulgas) * zrho
115	ENDDO
116
117	! mean free path of air (Pruppacher and Klett, 2010, p.417) [m]
118	mnfrpth=6.6E-8(1.01325E+5/pplay(ilon,ilev))(t_seri(ilon,ilev)/293.15)
119	! mnfrpth=2.eta/(zrhothvelair)
120	! mean free path of air (Prupp. Klett) in [10^-6 m]
121	! ZLAIR = 0.066 (1.01325E+5/PPLAY)(T_SERI/293.15)*1.E-06
122
123	! dynamic viscosity of air (Pruppacher and Klett, 2010, p.417) [kg/(m*s)]
124	IF (t_seri(ilon,ilev).GE.273.15) THEN
125	eta=(1.718+0.0049(t_seri(ilon,ilev)-273.15))1.E-5
126	ELSE
127	eta=(1.718+0.0049(t_seri(ilon,ilev)-273.15)-1.2E-5(t_seri(ilon,ilev)-273.15)*2)1.E-5
128	ENDIF
129
130	!--pre-compute the particle diffusion coefficient Di(i) from equation 18
131	Di=0.0
132	DO i=1, nbtr_bin
133	Kn(i)=mnfrpth/radius(i)
134	Di(i)=RKBOLt_seri(ilon,ilev)/(6.RPIradius(i)eta)(1.+Kn(i)(1.249+0.42*exp(-0.87/Kn(i))))
135	ENDDO
136
137	!--pre-compute the thermal velocity of a particle thvelpar(i) from equation 20
138	thvelpar=0.0
139	DO i=1, nbtr_bin
140	m_par(i)=4./3.RPIradius(i)*3.DENSO4(ilon,ilev)*1000.
141	thvelpar(i)=sqrt(8.RKBOLt_seri(ilon,ilev)/(RPI*m_par(i)))
142	ENDDO
143
144	!--pre-compute the particle mean free path mfppar(i) from equation 22
145	mfppar=0.0
146	DO i=1, nbtr_bin
147	mfppar(i)=8.Di(i)/(RPIthvelpar(i))
148	ENDDO
149
150	!--pre-compute the mean distance delta(i) from the center of a sphere reached by particles
151	!--leaving the surface of the sphere and traveling a distance of particle mfppar(i) from equation 21
152	delta=0.0
153	DO i=1, nbtr_bin
154	delta(i)=((2.radius(i)+mfppar(i))3.-(4.radius(i)2.+mfppar(i)2.)**1.5)/ &
155	& (6.radius(i)mfppar(i))-2.*radius(i)
156	ENDDO
157
158	!--pre-compute the beta(i,j) from equation 17
159	num=0.0
160	DO j=1, nbtr_bin
161	DO i=1, nbtr_bin
162	num=4.RPI(radius(i)+radius(j))*(Di(i)+Di(j))
163	denom=(radius(i)+radius(j))/(radius(i)+radius(j)+sqrt(delta(i)2.+delta(j)2.))+ &
164	& 4.(Di(i)+Di(j))/(sqrt(thvelpar(i)2.+thvelpar(j)2.)(radius(i)+radius(j)))
165	beta(i,j)=num/denom
166	ENDDO
167	ENDDO
168
169	!--external loop for equation 14
170	DO k=1, nbtr_bin
171
172	!--calculating denominator sum
173	denom=0.0
174	DO j=1, nbtr_bin
175	denom=denom+(1.-ff(k,j,k))beta(k,j)tr_t(ilon,ilev,j)
176	ENDDO
177
178	IF (k.EQ.1) THEN
179	!--calculate new concentration of smallest bin
180	tr_tp1(ilon,ilev,k)=tr_t(ilon,ilev,k)/(1.+pdtcoag*denom)
181	ELSE
182	!--calculating double sum terms in numerator of eq 14
183	num=0.0
184	DO j=1, k
185	numi=0.0
186	DO i=1, k-1
187	numi=numi+ff(i,j,k)beta(i,j)V(i)tr_tp1(ilon,ilev,i)tr_t(ilon,ilev,j)
188	ENDDO
189	num=num+numi
190	ENDDO
191
192	!--calculate new concentration of other bins
193	tr_tp1(ilon,ilev,k)=(V(k)tr_t(ilon,ilev,k)+pdtcoagnum)/(1.+pdtcoag*denom)/V(k)
194	ENDIF
195
196	ENDDO !--end of loop k
197
198	!--convert tracer concentration back from [number/m3] to [number/KgA] and write into tr_seri
199	DO i=1, nbtr_bin
200	tr_seri(ilon,ilev,i+nbtr_sulgas) = tr_tp1(ilon,ilev,i) / zrho
201	ENDDO
202
203	ENDIF ! IF IN STRATOSPHERE
204	ENDDO !--end of loop klev
205	ENDDO !--end of loop klon
206
207	END SUBROUTINE COAGULATE

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

Download in other formats:

Original Format