source: LMDZ5/trunk/libf/phylmd/StratAer/nucleation_tstep_mod.F90 @ 2902

Last change on this file since 2902 was 2695, checked in by oboucher, 8 years ago

Removing useless global variables
Interpolation only done once a month
Changing test orders in implicit coagulation routine as some compilers didn't like it
File size: 12.7 KB
Line 
1MODULE nucleation_tstep_mod
2
3CONTAINS
4
5SUBROUTINE nucleation_rate(rhoa,t_seri,pplay,rh,a_xm,b_xm,c_xm,nucl_rate,ntot,x)
6
7  USE aerophys
8  USE infotrac
9  USE YOMCST, ONLY : RPI, RD
10
11  IMPLICIT NONE
12
13  ! input variables
14  REAL rhoa !H2SO4 number density [molecules/cm3]
15  REAL t_seri
16  REAL pplay
17  REAL rh
18  REAL a_xm, b_xm, c_xm
19
20  ! output variables
21  REAL nucl_rate
22  REAL ntot ! total number of molecules in the critical cluster
23  REAL x    ! molefraction of H2SO4 in the critical cluster
24
25  ! local variables
26  REAL, PARAMETER                               :: k_B=1.3806E-23  ! Boltzmann constant [J/K]
27  REAL                                          :: jnuc !nucleation rate in 1/cm3s (10^-7-10^10 1/cm3s)
28  REAL                                          :: rc   !radius of the critical cluster in nm
29  REAL VH2SO4mol
30
31  ! call nucleation routine
32  CALL binapara(t_seri,rh,rhoa,jnuc,x,ntot,rc)
33
34  IF (ntot < 4.0) THEN
35    !set jnuc to collision rate of two H2SO4 molecules (following personal communication of Ulrike Niemeier and Hanna Vehkamäki)
36    VH2SO4mol=mH2SO4mol/(1.e-3*(a_xm+t_seri*(b_xm+t_seri*c_xm))) !cm3
37    jnuc = rhoa**2. *(3./4.*RPI)**(1./6.) *(12.*k_B*t_seri/mH2SO4mol)**0.5 &
38         & *100.*(2.*VH2SO4mol**(1./3.))**2. !1/(cm3s)
39    ntot=2.0
40    x=1.0
41  ENDIF
42
43  ! convert jnuc from particles/cm3/s to kg(H2SO4)/kgA/s
44  nucl_rate=jnuc*ntot*x*mH2SO4mol/(pplay/t_seri/RD/1.E6)
45
46END SUBROUTINE nucleation_rate
47
48!--------------------------------------------------------------------------------------------------
49
50SUBROUTINE nucleation_part(nucl_rate,ntot,x,dt,Vbin,tr_seri)
51
52  USE aerophys
53  USE infotrac
54
55  IMPLICIT NONE
56
57  ! input variables
58  REAL nucl_rate
59  REAL ntot ! total number of molecules in the critical cluster
60  REAL x    ! molefraction of H2SO4 in the critical cluster
61  REAL dt
62  REAL Vbin(nbtr_bin)
63
64  ! output variables
65  REAL tr_seri(nbtr)
66
67  ! local variables
68  INTEGER k
69  REAL Vnew
70  REAL ff(nbtr_bin)
71
72  ! dry volume of nucleated particle (at reference temperature)
73  Vnew=mH2SO4mol*ntot*x/dens_aer_dry
74
75  ! compute distribution factor for particles of intermediate size (from Jacobson 1994, equation 13)
76  ff(:)=0.0
77  DO k=1, nbtr_bin
78  ! CK 20160531: bug fix for first bin
79    IF (k.LE.(nbtr_bin-1)) THEN
80      IF (Vbin(k).LE.Vnew.AND.Vnew.LT.Vbin(k+1)) THEN
81        ff(k)= Vbin(k)/Vnew*(Vbin(k+1)-Vnew)/(Vbin(k+1)-Vbin(k))
82      ENDIF
83    ENDIF
84    IF (k.EQ.1.AND.Vnew.LE.Vbin(k)) THEN
85      ff(k)= 1.
86    ENDIF
87    IF (k.GT.1) THEN
88      IF (Vbin(k-1).LT.Vnew.AND.Vnew.LT.Vbin(k)) THEN
89        ff(k)= 1.-ff(k-1)
90      ENDIF
91    ENDIF
92    IF (k.EQ.nbtr_bin.AND.Vnew.GE.Vbin(k)) THEN
93      ff(k)= 1.
94    ENDIF
95  ENDDO
96  ! correction of ff for volume conservation
97  DO k=1, nbtr_bin         
98    ff(k)=ff(k)*Vnew/Vbin(k)
99  ENDDO
100
101  ! add nucleated H2SO4 to corresponding aerosol bin
102  DO k=1, nbtr_bin
103    tr_seri(k+nbtr_sulgas)=tr_seri(k+nbtr_sulgas)+nucl_rate*dt/(mH2SO4mol*ntot*x)*ff(k)
104  ENDDO
105
106END SUBROUTINE nucleation_part
107
108!---------------------------------------------------------------------------------------------------
109
110SUBROUTINE binapara(pt,prh,rhoa_in,jnuc,x,ntot,rc)
111
112
113  !    Fortran 90 subroutine binapara
114  !
115  !    Calculates parametrized values of nucleation rate,
116  !    mole fraction of sulphuric acid
117  !    total number of particles, and the radius of the critical cluster
118  !    in H2O-H2SO4 system IF temperature, saturatio ratio of water and
119  !    sulfuric acid concentration  are given.
120  !
121  !    Copyright (C) 2002 Hanna Vehkamäki
122  !
123  !    Division of Atmospheric Sciences
124  !    Department of Physical Sciences
125  !    P.O. Box 64
126  !    FIN-00014 University of Helsinki
127  !    Finland
128  !   
129  !    hanna.vehkamaki@helsinki.fi
130
131  IMPLICIT NONE
132
133  REAL :: pt,t     !temperature in K (190.15-300.15K)
134  REAL :: prh,rh    !saturatio ratio of water (0.0001-1)
135  REAL,intent(in) :: rhoa_in    !sulfuric acid concentration in 1/cm3 (10^4-10^11 1/cm3)
136  REAL,intent(out) :: jnuc    !nucleation rate in 1/cm3s (10^-7-10^10 1/cm3s)
137  REAL,intent(out) :: ntot !total number of molecules in the critical cluster (ntot>4)
138  REAL,intent(out) :: x    ! molefraction of H2SO4 in the critical cluster       
139  REAL,intent(out) :: rc    !radius of the critical cluster in nm     
140  REAL :: rhotres    ! treshold concentration of h2so4 (1/cm^3)
141                     ! which produces nucleation rate   1/(cm^3 s) as a function of rh and t
142  REAL rhoa
143
144! CK: use intermediate variables to avoid overwriting
145  t=pt
146  rh=prh
147  rhoa=rhoa_in
148
149  IF (t < 190.15) THEN
150!     print *,'Warning (ilon=',ilon,'ilev=',ilev,'): temperature < 190.15 K, using 190.15 K (T=',t,'K)'
151     t=190.15
152  ENDIF
153
154  IF (t > 300.15) THEN
155!     print *,'Warning (ilon=',ilon,'ilev=',ilev,'): temperature > 300.15 K, using 300.15 K'
156     t=300.15
157  ENDIF
158
159  IF (rh < 0.0001) THEN
160!     print *,'Warning (ilon=',ilon,'ilev=',ilev,'): saturation ratio of water < 0.0001, using 0.0001'
161     rh=0.0001
162  ENDIF
163
164  IF (rh > 1.0) THEN
165!     print *,'Warning (ilon=',ilon,'ilev=',ilev,'): saturation ratio of water > 1 using 1'
166!     print *, 'rh=',rh
167     rh=1.0
168  ENDIF
169
170  IF (rhoa < 1.e4) THEN
171!     print *,'Warning (ilon=',ilon,'ilev=',ilev,'): sulfuric acid concentration < 1e4 1/cm3, using 1e4 1/cm3'
172     rhoa=1.e4
173  ENDIF
174
175  IF (rhoa > 1.e11) THEN
176!     print *,'Warning (ilon=',ilon,'ilev=',ilev,'): sulfuric acid concentration > 1e11 1/cm3, using 1e11 1/cm3'
177     rhoa=1.e11
178  ENDIF
179
180  x=  0.7409967177282139 - 0.002663785665140117*t + 0.002010478847383187*Log(rh)  &
181       & - 0.0001832894131464668*t*Log(rh) + 0.001574072538464286*Log(rh)**2      &
182       & - 0.00001790589121766952*t*Log(rh)**2 + 0.0001844027436573778*Log(rh)**3 &
183       & -  1.503452308794887e-6*t*Log(rh)**3 - 0.003499978417957668*Log(rhoa)    &
184       & + 0.0000504021689382576*t*Log(rhoa)
185
186  jnuc= 0.1430901615568665 + 2.219563673425199*t - 0.02739106114964264*t**2 +  &
187       &  0.00007228107239317088*t**3 + 5.91822263375044/x +                   &
188       &  0.1174886643003278*Log(rh) + 0.4625315047693772*t*Log(rh) -          &
189       &  0.01180591129059253*t**2*Log(rh) +                                   &
190       &  0.0000404196487152575*t**3*Log(rh) + (15.79628615047088*Log(rh))/x - &
191       &  0.215553951893509*Log(rh)**2 - 0.0810269192332194*t*Log(rh)**2 +     &
192       &  0.001435808434184642*t**2*Log(rh)**2 -                               &
193       &  4.775796947178588e-6*t**3*Log(rh)**2 -                               &
194       &  (2.912974063702185*Log(rh)**2)/x - 3.588557942822751*Log(rh)**3 +    &
195       &  0.04950795302831703*t*Log(rh)**3 -                                   &
196       &  0.0002138195118737068*t**2*Log(rh)**3 +                              &
197       &  3.108005107949533e-7*t**3*Log(rh)**3 -                               &
198       &  (0.02933332747098296*Log(rh)**3)/x +                                 &
199       &  1.145983818561277*Log(rhoa) -                                        &
200       &  0.6007956227856778*t*Log(rhoa) +                                     &
201       &  0.00864244733283759*t**2*Log(rhoa) -                                 &
202       &  0.00002289467254710888*t**3*Log(rhoa) -                              &
203       &  (8.44984513869014*Log(rhoa))/x +                                     &
204       &  2.158548369286559*Log(rh)*Log(rhoa) +                                &
205       &  0.0808121412840917*t*Log(rh)*Log(rhoa) -                             &
206       &  0.0004073815255395214*t**2*Log(rh)*Log(rhoa) -                       &
207       &  4.019572560156515e-7*t**3*Log(rh)*Log(rhoa) +                        &
208       &  (0.7213255852557236*Log(rh)*Log(rhoa))/x +                           &
209       &  1.62409850488771*Log(rh)**2*Log(rhoa) -                              &
210       &  0.01601062035325362*t*Log(rh)**2*Log(rhoa) +                         &
211       &  0.00003771238979714162*t**2*Log(rh)**2*Log(rhoa) +                   &
212       &  3.217942606371182e-8*t**3*Log(rh)**2*Log(rhoa) -                     &
213       &  (0.01132550810022116*Log(rh)**2*Log(rhoa))/x +                       &
214       &  9.71681713056504*Log(rhoa)**2 -                                      &
215       &  0.1150478558347306*t*Log(rhoa)**2 +                                  &
216       &  0.0001570982486038294*t**2*Log(rhoa)**2 +                            &
217       &  4.009144680125015e-7*t**3*Log(rhoa)**2 +                             &
218       &  (0.7118597859976135*Log(rhoa)**2)/x -                                &
219       &  1.056105824379897*Log(rh)*Log(rhoa)**2 +                             &
220       &  0.00903377584628419*t*Log(rh)*Log(rhoa)**2 -                         &
221       &  0.00001984167387090606*t**2*Log(rh)*Log(rhoa)**2 +                   &
222       &  2.460478196482179e-8*t**3*Log(rh)*Log(rhoa)**2 -                     &
223       &  (0.05790872906645181*Log(rh)*Log(rhoa)**2)/x -                       &
224       &  0.1487119673397459*Log(rhoa)**3 +                                    &
225       &  0.002835082097822667*t*Log(rhoa)**3 -                                &
226       &  9.24618825471694e-6*t**2*Log(rhoa)**3 +                              &
227       &  5.004267665960894e-9*t**3*Log(rhoa)**3 -                             &
228       &  (0.01270805101481648*Log(rhoa)**3)/x
229  jnuc=exp(jnuc) !1/(cm3s)
230
231  ntot =-0.002954125078716302 - 0.0976834264241286*t + 0.001024847927067835*t**2 - 2.186459697726116e-6*t**3 -    &
232       &   0.1017165718716887/x - 0.002050640345231486*Log(rh) - 0.007585041382707174*t*Log(rh) +                 &
233       &   0.0001926539658089536*t**2*Log(rh) - 6.70429719683894e-7*t**3*Log(rh) -                                &
234       &   (0.2557744774673163*Log(rh))/x + 0.003223076552477191*Log(rh)**2 + 0.000852636632240633*t*Log(rh)**2 - &
235       &   0.00001547571354871789*t**2*Log(rh)**2 + 5.666608424980593e-8*t**3*Log(rh)**2 +                        &
236       &   (0.03384437400744206*Log(rh)**2)/x + 0.04743226764572505*Log(rh)**3 -                                  &
237       &   0.0006251042204583412*t*Log(rh)**3 + 2.650663328519478e-6*t**2*Log(rh)**3 -                            &
238       &   3.674710848763778e-9*t**3*Log(rh)**3 - (0.0002672510825259393*Log(rh)**3)/x -                          &
239       &   0.01252108546759328*Log(rhoa) + 0.005806550506277202*t*Log(rhoa) -                                     &
240       &   0.0001016735312443444*t**2*Log(rhoa) + 2.881946187214505e-7*t**3*Log(rhoa) +                           &
241       &   (0.0942243379396279*Log(rhoa))/x - 0.0385459592773097*Log(rh)*Log(rhoa) -                              &
242       &   0.0006723156277391984*t*Log(rh)*Log(rhoa) + 2.602884877659698e-6*t**2*Log(rh)*Log(rhoa) +              &
243       &   1.194163699688297e-8*t**3*Log(rh)*Log(rhoa) - (0.00851515345806281*Log(rh)*Log(rhoa))/x -              &
244       &   0.01837488495738111*Log(rh)**2*Log(rhoa) + 0.0001720723574407498*t*Log(rh)**2*Log(rhoa) -              &
245       &   3.717657974086814e-7*t**2*Log(rh)**2*Log(rhoa) -                                                       &
246       &   5.148746022615196e-10*t**3*Log(rh)**2*Log(rhoa) +                                                      &
247       &   (0.0002686602132926594*Log(rh)**2*Log(rhoa))/x - 0.06199739728812199*Log(rhoa)**2 +                    &
248       &   0.000906958053583576*t*Log(rhoa)**2 - 9.11727926129757e-7*t**2*Log(rhoa)**2 -                          &
249       &   5.367963396508457e-9*t**3*Log(rhoa)**2 - (0.007742343393937707*Log(rhoa)**2)/x +                       &
250       &   0.0121827103101659*Log(rh)*Log(rhoa)**2 - 0.0001066499571188091*t*Log(rh)*Log(rhoa)**2 +               &
251       &   2.534598655067518e-7*t**2*Log(rh)*Log(rhoa)**2 -                                                       &
252       &   3.635186504599571e-10*t**3*Log(rh)*Log(rhoa)**2 +                                                      &
253       &   (0.0006100650851863252*Log(rh)*Log(rhoa)**2)/x + 0.0003201836700403512*Log(rhoa)**3 -                  &
254       &   0.0000174761713262546*t*Log(rhoa)**3 + 6.065037668052182e-8*t**2*Log(rhoa)**3 -                        &
255       &   1.421771723004557e-11*t**3*Log(rhoa)**3 + (0.0001357509859501723*Log(rhoa)**3)/x
256  ntot=exp(ntot)
257
258  rc=exp(-1.6524245+0.42316402*x+0.33466487*log(ntot)) !nm
259
260  IF (jnuc < 1.e-7) THEN
261!     print *,'Warning (ilon=',ilon,'ilev=',ilev'): nucleation rate < 1e-7/cm3s, using 0.0/cm3s,'
262     jnuc=0.0
263  ENDIF
264
265  IF (jnuc > 1.e10) THEN
266!     print *,'Warning (ilon=',ilon,'ilev=',ilev, &
267!        & '): nucleation rate > 1e10/cm3s, using 1e10/cm3s'
268     jnuc=1.e10
269  ENDIF
270
271  rhotres=exp( -279.2430007512709 + 11.73439886096903*rh + 22700.92970508331/t &
272       & - (1088.644983466801*rh)/t + 1.144362942094912*t                      &
273       & - 0.03023314602163684*rh*t - 0.001302541390154324*t**2                &
274       & - 6.386965238433532*Log(rh) + (854.980361026715*Log(rh))/t            &
275       & + 0.00879662256826497*t*Log(rh)) !1/cm3
276
277  RETURN
278
279END SUBROUTINE binapara
280
281END MODULE nucleation_tstep_mod
Note: See TracBrowser for help on using the repository browser.