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

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