source: LMDZ6/branches/Amaury_dev/libf/phylmd/StratAer/nucleation_tstep_mod.F90 @ 5409

Last change on this file since 5409 was 5105, checked in by abarral, 5 months ago

Replace 1DUTILS.h by module lmdz_1dutils.f90
Replace 1DConv.h by module lmdz_old_1dconv.f90 (it's only used by old_* files)
Convert *.F to *.f90
Fix gradsdef.h formatting
Remove unnecessary "RETURN" at the end of functions/subroutines
  • Property svn:keywords set to Id
File size: 47.8 KB
Line 
1
2! $Id: nucleation_tstep_mod.F90 5105 2024-07-23 17:14:34Z evignon $
3
4MODULE nucleation_tstep_mod
5
6CONTAINS
7
8SUBROUTINE nucleation_rate(rhoa,t_seri,pplay,rh,a_xm,b_xm,c_xm,nucl_rate,ntot_n,x_n)
9
10  USE aerophys
11  USE infotrac_phy
12  USE strataer_local_var_mod, ONLY: flag_new_nucl
13  USE lmdz_yomcst
14 
15  IMPLICIT NONE
16
17  ! input variables
18  REAL rhoa    ! H2SO4 number density [molecules/cm3]
19  REAL t_seri  ! temperature (K)
20  REAL pplay   ! pressure (Pa)
21  REAL rh      ! relative humidity (between 0 and 1)
22  REAL a_xm, b_xm, c_xm
23
24  ! output variables
25  REAL nucl_rate
26  REAL ntot_n ! total number of molecules in the critical cluster for neutral nucleation
27  REAL x_n    ! mole fraction of H2SO4 in the critical cluster for neutral nucleation
28
29  ! local variables
30  REAL jnuc_n ! nucleation rate in 1/cm3s (10^-7-10^10 1/cm3s) for neutral nucleation
31  REAL rc_n   ! radius of the critical cluster in nm for neutral nucleation
32  REAL na_n   ! sulfuric acid molecules in the neutral critical cluster (NOT IN USE)
33  LOGICAL kinetic_n ! true if kinetic neutral nucleation (NOT IN USE)
34  LOGICAL kinetic_i ! true if kinetic ion-induced nucleation (NOT IN USE)
35  REAL VH2SO4mol
36  REAL ntot_i, x_i, jnuc_i, rc_i, na_i, n_i ! quantities for charged nucleation (NOT IN USE)
37  REAL csi     ! Ion condensation sink (s-1) NOT IN USE
38  REAL airn    ! Air molecule concentration in (cm-3) NOT IN USE
39  REAL ipr     ! Ion pair production rate (cm-3 s-1) NOT IN USE
40
41  ! CALL nucleation routine
42  IF (.NOT.flag_new_nucl) THEN
43    ! Use older routine from Hanna Vehkamäki (FMI)
44    CALL binapara(t_seri,rh,rhoa,jnuc_n,x_n,ntot_n,rc_n)
45    ! when total number of molecules is too small
46    ! then set jnuc_n to collision rate of two H2SO4 molecules (following personal communication of Ulrike Niemeier and Hanna Vehkamäki)
47    IF (ntot_n < 4.0) THEN
48      VH2SO4mol=mH2SO4mol/(1.E-3*(a_xm+t_seri*(b_xm+t_seri*c_xm))) !cm3
49      jnuc_n = rhoa**2. *(3./4.*RPI)**(1./6.) *(12.*RKBOL*t_seri/mH2SO4mol)**0.5 &
50   *100.*(2.*VH2SO4mol**(1./3.))**2. !1/(cm3s)
51      ntot_n=2.0
52      x_n=1.0
53    ENDIF
54  ELSE
55    ! Use new routine from Anni Maattanen (LATMOS)
56    csi=0.0   ! no charged nucleation for now
57    ipr=-1.0  ! dummy value to make sure charged nucleation does not occur
58    airn=0.0  ! NOT IN USE
59!   airn=pplay/t_seri/RD/1.E3*RNAVO/RMD ! molec cm-3 (for future use, to be confirmed)
60    CALL newbinapara(t_seri,rh,rhoa,csi,airn,ipr,jnuc_n,ntot_n,jnuc_i,ntot_i, &
61   x_n,x_i,na_n,na_i,rc_n,rc_i,n_i,kinetic_n,kinetic_i)
62  ENDIF
63
64  ! convert jnuc_n from particles/cm3/s to kg(H2SO4)/kgA/s
65  nucl_rate=jnuc_n*ntot_n*x_n*mH2SO4mol/(pplay/t_seri/RD/1.E6)
66
67END SUBROUTINE nucleation_rate
68
69!--------------------------------------------------------------------------------------------------
70
71SUBROUTINE nucleation_part(nucl_rate,ntot,x,dt,tr_seri)
72
73  USE aerophys
74  USE infotrac_phy
75  USE strataer_local_var_mod, ONLY: Vbin
76 
77  IMPLICIT NONE
78
79  ! input variable
80  REAL nucl_rate
81  REAL ntot ! total number of molecules in the critical cluster
82  REAL x    ! mole raction of H2SO4 in the critical cluster
83  REAL dt
84 
85  ! output variable
86  REAL tr_seri(nbtr)
87
88  ! local variable
89  INTEGER k
90  REAL Vnew
91  REAL ff(nbtr_bin)
92
93  ! dry volume of nucleated particle (at reference temperature)
94  !  dens_aer_dry (density of dry aerosol []kg/m3)
95  !  H2SO4mol: mass of an H2SO4 molecule [kg]
96  !  Vnew [m3]
97  Vnew=mH2SO4mol*ntot*x/dens_aer_dry
98
99  ! compute distribution factor for particles of intermediate size (from Jacobson 1994, equation 13)
100  ff(:)=0.0
101  DO k=1, nbtr_bin
102  ! CK 20160531: bug fix for first bin
103    IF (k<=nbtr_bin-1) THEN
104      IF (Vbin(k)<=Vnew.AND.Vnew<Vbin(k+1)) THEN
105        ff(k)= Vbin(k)/Vnew*(Vbin(k+1)-Vnew)/(Vbin(k+1)-Vbin(k))
106      ENDIF
107    ENDIF
108    IF (k==1.AND.Vnew<=Vbin(k)) THEN
109      ff(k)= 1.
110    ENDIF
111    IF (k>1) THEN
112      IF (Vbin(k-1)<Vnew.AND.Vnew<Vbin(k)) THEN
113        ff(k)= 1.-ff(k-1)
114      ENDIF
115    ENDIF
116    IF (k==nbtr_bin.AND.Vnew>=Vbin(k)) THEN
117      ff(k)= 1.
118    ENDIF
119  ENDDO
120  ! correction of ff for volume conservation
121  DO k=1, nbtr_bin         
122    ff(k)=ff(k)*Vnew/Vbin(k)
123  ENDDO
124
125  ! add nucleated H2SO4 to corresponding aerosol bin
126  DO k=1, nbtr_bin
127    tr_seri(k+nbtr_sulgas)=tr_seri(k+nbtr_sulgas)+nucl_rate*dt/(mH2SO4mol*ntot*x)*ff(k)
128  ENDDO
129
130END SUBROUTINE nucleation_part
131
132!---------------------------------------------------------------------------------------------------
133
134SUBROUTINE binapara(pt,prh,rhoa_in,jnuc,x,ntot,rc)
135
136
137  !    Fortran 90 SUBROUTINE binapara
138
139  !    Calculates parametrized values of nucleation rate,
140  !    mole fraction of sulphuric acid
141  !    total number of particles, and the radius of the critical cluster
142  !    in H2O-H2SO4 system IF temperature, saturatio ratio of water and
143  !    sulfuric acid concentration  are given.
144
145  !    Copyright (C) 2002 Hanna Vehkamäki
146
147  !    Division of Atmospheric Sciences
148  !    Department of Physical Sciences
149  !    P.O. Box 64
150  !    FIN-00014 University of Helsinki
151  !    Finland
152
153  !    hanna.vehkamaki@helsinki.fi
154
155  IMPLICIT NONE
156
157  REAL :: pt,t     !temperature in K (190.15-300.15K)
158  REAL :: prh,rh    !saturation ratio of water (0.0001-1)
159  REAL,INTENT(IN) :: rhoa_in    !sulfuric acid concentration in 1/cm3 (10^4-10^11 1/cm3)
160  REAL,INTENT(OUT) :: jnuc    !nucleation rate in 1/cm3s (10^-7-10^10 1/cm3s)
161  REAL,INTENT(OUT) :: ntot !total number of molecules in the critical cluster (ntot>4)
162  REAL,INTENT(OUT) :: x    ! mole fraction of H2SO4 in the critical cluster       
163  REAL,INTENT(OUT) :: rc    !radius of the critical cluster in nm     
164  REAL :: rhotres    ! threshold concentration of h2so4 (1/cm^3)
165                     ! which produces nucleation rate   1/(cm^3 s) as a function of rh and t
166  REAL rhoa
167
168! CK: use intermediate variable to avoid overwriting
169  t=pt
170  rh=prh
171  rhoa=rhoa_in
172
173  IF (t < 190.15) THEN
174     t=190.15
175  ENDIF
176
177  IF (t > 300.15) THEN
178     t=300.15
179  ENDIF
180
181  IF (rh < 0.0001) THEN
182     rh=0.0001
183  ENDIF
184
185  IF (rh > 1.0) THEN
186     rh=1.0
187  ENDIF
188
189  IF (rhoa < 1.e4) THEN
190     rhoa=1.e4
191  ENDIF
192
193  IF (rhoa > 1.e11) THEN
194     rhoa=1.e11
195  ENDIF
196
197  x=  0.7409967177282139 - 0.002663785665140117*t + 0.002010478847383187*LOG(rh)  &
198   - 0.0001832894131464668*t*LOG(rh) + 0.001574072538464286*LOG(rh)**2      &
199   - 0.00001790589121766952*t*LOG(rh)**2 + 0.0001844027436573778*LOG(rh)**3 &
200   -  1.503452308794887E-6*t*LOG(rh)**3 - 0.003499978417957668*LOG(rhoa)    &
201   + 0.0000504021689382576*t*LOG(rhoa)
202
203  jnuc= 0.1430901615568665 + 2.219563673425199*t - 0.02739106114964264*t**2 + &
204    0.00007228107239317088*t**3 + 5.91822263375044/x +                  &
205    0.1174886643003278*LOG(rh) + 0.4625315047693772*t*LOG(rh) -          &
206    0.01180591129059253*t**2*LOG(rh) +                                  &
207    0.0000404196487152575*t**3*LOG(rh) + (15.79628615047088*LOG(rh))/x - &
208    0.215553951893509*LOG(rh)**2 - 0.0810269192332194*t*LOG(rh)**2 +    &
209    0.001435808434184642*t**2*LOG(rh)**2 -                               &
210    4.775796947178588E-6*t**3*LOG(rh)**2 -                               &
211    (2.912974063702185*LOG(rh)**2)/x - 3.588557942822751*LOG(rh)**3 +   &
212    0.04950795302831703*t*LOG(rh)**3 -                                   &
213    0.0002138195118737068*t**2*LOG(rh)**3 +                             &
214    3.108005107949533E-7*t**3*LOG(rh)**3 -                               &
215    (0.02933332747098296*LOG(rh)**3)/x +                                &
216    1.145983818561277*LOG(rhoa) -                                        &
217    0.6007956227856778*t*LOG(rhoa) +                                    &
218    0.00864244733283759*t**2*LOG(rhoa) -                                 &
219    0.00002289467254710888*t**3*LOG(rhoa) -                              &
220    (8.44984513869014*LOG(rhoa))/x +                                    &
221    2.158548369286559*LOG(rh)*LOG(rhoa) +                               &
222    0.0808121412840917*t*LOG(rh)*LOG(rhoa) -                             &
223    0.0004073815255395214*t**2*LOG(rh)*LOG(rhoa) -                       &
224    4.019572560156515E-7*t**3*LOG(rh)*LOG(rhoa) +                       &
225    (0.7213255852557236*LOG(rh)*LOG(rhoa))/x +                          &
226    1.62409850488771*LOG(rh)**2*LOG(rhoa) -                              &
227    0.01601062035325362*t*LOG(rh)**2*LOG(rhoa) +                        &
228    0.00003771238979714162*t**2*LOG(rh)**2*LOG(rhoa) +                  &
229    3.217942606371182E-8*t**3*LOG(rh)**2*LOG(rhoa) -                     &
230    (0.01132550810022116*LOG(rh)**2*LOG(rhoa))/x +                      &
231    9.71681713056504*LOG(rhoa)**2 -                                      &
232    0.1150478558347306*t*LOG(rhoa)**2 +                                 &
233    0.0001570982486038294*t**2*LOG(rhoa)**2 +                           &
234    4.009144680125015E-7*t**3*LOG(rhoa)**2 +                            &
235    (0.7118597859976135*LOG(rhoa)**2)/x -                                &
236    1.056105824379897*LOG(rh)*LOG(rhoa)**2 +                            &
237    0.00903377584628419*t*LOG(rh)*LOG(rhoa)**2 -                         &
238    0.00001984167387090606*t**2*LOG(rh)*LOG(rhoa)**2 +                  &
239    2.460478196482179E-8*t**3*LOG(rh)*LOG(rhoa)**2 -                     &
240    (0.05790872906645181*LOG(rh)*LOG(rhoa)**2)/x -                       &
241    0.1487119673397459*LOG(rhoa)**3 +                                   &
242    0.002835082097822667*t*LOG(rhoa)**3 -                                &
243    9.24618825471694E-6*t**2*LOG(rhoa)**3 +                             &
244    5.004267665960894E-9*t**3*LOG(rhoa)**3 -                             &
245    (0.01270805101481648*LOG(rhoa)**3)/x
246  jnuc=EXP(jnuc) !1/(cm3s)
247
248  ntot =-0.002954125078716302 - 0.0976834264241286*t + 0.001024847927067835*t**2 - 2.186459697726116E-6*t**3 -    &
249     0.1017165718716887/x - 0.002050640345231486*LOG(rh) - 0.007585041382707174*t*LOG(rh) +                &
250     0.0001926539658089536*t**2*LOG(rh) - 6.70429719683894E-7*t**3*LOG(rh) -                                &
251     (0.2557744774673163*LOG(rh))/x + 0.003223076552477191*LOG(rh)**2 + 0.000852636632240633*t*LOG(rh)**2 - &
252     0.00001547571354871789*t**2*LOG(rh)**2 + 5.666608424980593E-8*t**3*LOG(rh)**2 +                       &
253     (0.03384437400744206*LOG(rh)**2)/x + 0.04743226764572505*LOG(rh)**3 -                                  &
254     0.0006251042204583412*t*LOG(rh)**3 + 2.650663328519478E-6*t**2*LOG(rh)**3 -                            &
255     3.674710848763778E-9*t**3*LOG(rh)**3 - (0.0002672510825259393*LOG(rh)**3)/x -                          &
256     0.01252108546759328*LOG(rhoa) + 0.005806550506277202*t*LOG(rhoa) -                                     &
257     0.0001016735312443444*t**2*LOG(rhoa) + 2.881946187214505E-7*t**3*LOG(rhoa) +                          &
258     (0.0942243379396279*LOG(rhoa))/x - 0.0385459592773097*LOG(rh)*LOG(rhoa) -                              &
259     0.0006723156277391984*t*LOG(rh)*LOG(rhoa) + 2.602884877659698E-6*t**2*LOG(rh)*LOG(rhoa) +             &
260     1.194163699688297E-8*t**3*LOG(rh)*LOG(rhoa) - (0.00851515345806281*LOG(rh)*LOG(rhoa))/x -              &
261     0.01837488495738111*LOG(rh)**2*LOG(rhoa) + 0.0001720723574407498*t*LOG(rh)**2*LOG(rhoa) -              &
262     3.717657974086814E-7*t**2*LOG(rh)**2*LOG(rhoa) -                                                       &
263     5.148746022615196E-10*t**3*LOG(rh)**2*LOG(rhoa) +                                                     &
264     (0.0002686602132926594*LOG(rh)**2*LOG(rhoa))/x - 0.06199739728812199*LOG(rhoa)**2 +                   &
265     0.000906958053583576*t*LOG(rhoa)**2 - 9.11727926129757E-7*t**2*LOG(rhoa)**2 -                          &
266     5.367963396508457E-9*t**3*LOG(rhoa)**2 - (0.007742343393937707*LOG(rhoa)**2)/x +                      &
267     0.0121827103101659*LOG(rh)*LOG(rhoa)**2 - 0.0001066499571188091*t*LOG(rh)*LOG(rhoa)**2 +              &
268     2.534598655067518E-7*t**2*LOG(rh)*LOG(rhoa)**2 -                                                       &
269     3.635186504599571E-10*t**3*LOG(rh)*LOG(rhoa)**2 +                                                     &
270     (0.0006100650851863252*LOG(rh)*LOG(rhoa)**2)/x + 0.0003201836700403512*LOG(rhoa)**3 -                  &
271     0.0000174761713262546*t*LOG(rhoa)**3 + 6.065037668052182E-8*t**2*LOG(rhoa)**3 -                        &
272     1.421771723004557E-11*t**3*LOG(rhoa)**3 + (0.0001357509859501723*LOG(rhoa)**3)/x
273  ntot=EXP(ntot)
274
275  rc=EXP(-1.6524245+0.42316402*x+0.33466487*LOG(ntot)) !nm
276
277  IF (jnuc < 1.E-7) THEN
278     jnuc=0.0
279  ENDIF
280
281  IF (jnuc > 1.e10) THEN
282     jnuc=1.e10
283  ENDIF
284
285  rhotres=EXP( -279.2430007512709 + 11.73439886096903*rh + 22700.92970508331/t &
286   - (1088.644983466801*rh)/t + 1.144362942094912*t                      &
287   - 0.03023314602163684*rh*t - 0.001302541390154324*t**2                &
288   - 6.386965238433532*LOG(rh) + (854.980361026715*LOG(rh))/t            &
289   + 0.00879662256826497*t*LOG(rh)) !1/cm3
290
291
292
293END SUBROUTINE binapara
294
295!---------------------------------------------------------------------------------------------------
296
297SUBROUTINE newbinapara(t,satrat,rhoa,csi,airn,ipr,jnuc_n_real,ntot_n_real,jnuc_i_real,ntot_i_real,        &
298     x_n_real,x_i_real,na_n_real,na_i_real,rc_n_real,rc_i_real,n_i_real,                &
299     kinetic_n,kinetic_i)
300
301  !    Fortran 90 SUBROUTINE newbinapara
302
303  !    Calculates parametrized values for neutral and ion-induced sulfuric acid-water particle formation rate
304  !    of critical clusters,
305  !    number of particle in the critical clusters, the radii of the critical clusters
306  !    in H2O-H2SO4-ion system if temperature, saturation ratio of water, sulfuric acid concentration,
307  !    and, optionally, either condensation sink due to pre-existing particle and ion pair production rate,
308  !    or atmospheric concentration of negative ions are given.
309
310  !    The code calculates also the kinetic limit and the particle formation rate
311  !    above this limit (in which case we set ntot=1 and na=1)
312
313  !    Copyright (C)2018 Määttänen et al. 2018
314
315  !    anni.maattanen@latmos.ipsl.fr
316  !    joonas.merikanto@fmi.fi
317  !    hanna.vehkamaki@helsinki.fi
318
319  !    References
320  !    A. Määttänen, J. Merikanto, H. Henschel, J. Duplissy, R. Makkonen,
321  !    I. K. Ortega and H. Vehkamäki (2018), New parameterizations for
322  !    neutral and ion-induced sulfuric acid-water particle formation in
323  !    nucleation and kinetic regimes, J. Geophys. Res. Atmos., 122, doi:10.1002/2017JD027429.
324
325  !    Brasseur, G., and A.  Chatel (1983),  paper  presented  at  the  9th  Annual  Meeting  of  the 
326  !    European Geophysical Society, Leeds, Great Britain, August 1982.
327
328  !    Dunne, Eimear M., et al.(2016), Global atmospheric particle formation from CERN CLOUD measurements,
329  !    Science 354.6316, 1119-1124.   
330
331  USE aerophys
332  USE lmdz_yomcst
333
334  IMPLICIT NONE
335
336  !----------------------------------------------------
337 
338  !Global intent in
339  REAL,INTENT(IN) :: t         ! temperature in K
340  REAL,INTENT(IN) :: satrat    ! saturatio ratio of water (between zero and 1)
341  REAL,INTENT(IN) :: rhoa      ! sulfuric acid concentration in 1/cm3
342  REAL,INTENT(IN) :: csi       ! Ion condensation sink (s-1)
343  REAL,INTENT(IN) :: airn      ! Air molecule concentration in (cm-3)
344  REAL,INTENT(IN) :: ipr       ! Ion pair production rate (cm-3 s-1)
345  !Global intent out
346  REAL,INTENT(OUT) :: jnuc_n_real   ! Neutral nucleation rate in 1/cm3s (J>10^-7 1/cm3s)
347  REAL,INTENT(OUT) :: ntot_n_real   ! total number of molecules in the neutral critical cluster
348  REAL,INTENT(OUT) :: jnuc_i_real   ! Charged nucleation rate in 1/cm3s (J>10^-7 1/cm3s)
349  REAL,INTENT(OUT) :: ntot_i_real   ! total number of molecules in the charged critical cluster
350  REAL,INTENT(OUT) :: x_n_real      ! mole fraction of H2SO4 in the neutral critical cluster
351  REAL,INTENT(OUT) :: x_i_real      ! mole fraction of H2SO4 in the charged critical cluster
352                                           ! (note that x_n=x_i in nucleation regime)
353  REAL,INTENT(OUT) :: na_n_real     ! sulfuric acid molecules in the neutral critical cluster
354  REAL,INTENT(OUT) :: na_i_real     ! sulfuric molecules in the charged critical cluster
355  REAL,INTENT(OUT) :: rc_n_real     ! radius of the charged critical cluster in nm
356  REAL,INTENT(OUT) :: rc_i_real     ! radius of the charged critical cluster in nm
357  REAL,INTENT(OUT) :: n_i_real      ! number of ion pairs in air (cm-3)
358  LOGICAL,INTENT(OUT)  :: kinetic_n        ! true if kinetic neutral nucleation
359  LOGICAL,INTENT(OUT)  :: kinetic_i        ! true if kinetic ion-induced nucleation
360
361  ! Local
362  DOUBLE PRECISION :: jnuc_n      ! Neutral nucleation rate in 1/cm3s (J>10^-7 1/cm3s)
363  DOUBLE PRECISION :: ntot_n      ! total number of molecules in the neutral critical cluster
364  DOUBLE PRECISION :: jnuc_i      ! Charged nucleation rate in 1/cm3s (J>10^-7 1/cm3s)
365  DOUBLE PRECISION :: ntot_i      ! total number of molecules in the charged critical cluster
366  DOUBLE PRECISION :: x_n         ! mole fraction of H2SO4 in the neutral critical cluster
367  DOUBLE PRECISION :: x_i         ! mole fraction of H2SO4 in the charged critical cluster
368                                              ! (note that x_n=x_i in nucleation regime)
369  DOUBLE PRECISION :: na_n        ! sulfuric acid molecules in the neutral critical cluster
370  DOUBLE PRECISION :: na_i        ! sulfuric molecules in the charged critical cluster
371  DOUBLE PRECISION :: rc_n        ! radius of the charged critical cluster in nm
372  DOUBLE PRECISION :: rc_i        ! radius of the charged critical cluster in nm
373  DOUBLE PRECISION :: n_i         ! number of ion pairs in air (cm-3)
374  DOUBLE PRECISION :: x           ! mole fraction of H2SO4 in the critical cluster
375  DOUBLE PRECISION :: satratln    ! bounded water saturation ratio for neutral case (between 5.E-6 - 1.0)
376  DOUBLE PRECISION :: satratli    ! bounded water saturation ratio for ion-induced case (between 1.E-7 - 0.95)
377  DOUBLE PRECISION :: rhoaln      ! bounded concentration of h2so4 for neutral case (between 10^10 - 10^19 m-3)
378  DOUBLE PRECISION :: rhoali      ! bounded concentration of h2so4 for ion-induced case (between 10^10 - 10^22 m-3)
379  DOUBLE PRECISION :: tln         ! bounded temperature for neutral case (between 165-400 K)
380  DOUBLE PRECISION :: tli         ! bounded temperature for ion-induced case (195-400 K)
381  DOUBLE PRECISION :: kinrhotresn ! threshold sulfuric acid for neutral kinetic nucleation   
382  DOUBLE PRECISION :: kinrhotresi ! threshold sulfuric acid for ion-induced kinetic nucleation
383  DOUBLE PRECISION :: jnuc_i1     ! Ion-induced rate for n_i=1 cm-3
384  DOUBLE PRECISION :: xloss       ! Ion loss rate
385  DOUBLE PRECISION :: recomb      ! Ion-ion recombination rate
386
387  !--- 0) Initializations:
388
389  kinetic_n=.FALSE.
390  kinetic_i=.FALSE.
391  jnuc_n=0.0
392  jnuc_i=0.0
393  ntot_n=0.0
394  ntot_i=0.0
395  na_n=0.0
396  na_i=0.0
397  rc_n=0.0
398  rc_i=0.0
399  x=0.0
400  x_n=0.0
401  x_i=0.0
402  satratln=satrat
403  satratli=satrat
404  rhoaln=rhoa
405  rhoali=rhoa
406  tln=t
407  tli=t
408  n_i=0.0
409
410  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
411
412  !Boundary values according to parameterization limits   
413
414  !Temperature bounds
415  IF (t<=165.) THEN
416     tln=165.0
417  ENDIF
418  IF (t<=195.) THEN
419     tli=195.0
420  ENDIF
421  IF (t>=400.) THEN
422     tln=400.
423     tli=400.
424  ENDIF
425 
426  ! Saturation ratio bounds
427  IF (satrat<1.E-7) THEN
428     satratli=1.E-7
429  ENDIF
430  IF (satrat<1.E-5) THEN
431     satratln=1.E-5
432  ENDIF
433  IF (satrat>0.95) THEN
434     satratli=0.95
435  ENDIF
436  IF (satrat>1.0) THEN
437     satratln=1.0
438  ENDIF
439 
440  ! Sulfuric acid concentration bounds
441  IF (rhoa<=1.E4) THEN
442     rhoaln=1.E4
443     rhoali=1.E4
444  ENDIF
445  IF (rhoa>1.E13) THEN
446     rhoaln=1.E13
447  ENDIF
448  IF (rhoa>1.E16) THEN
449     rhoali=1.E16
450  ENDIF
451 
452  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
453 
454  !Critical cluster composition (valid for both cases, bounds not used here)
455  x_n=  7.9036365428891719E-1 - 2.8414059650092153E-3*tln + 1.4976802556584141E-2*LOG(satratln)  &
456   - 2.4511581740839115E-4*tln*LOG(satratln) + 3.4319869471066424E-3 *LOG(satratln)**2     &
457   - 2.8799393617748428E-5*tln*LOG(satratln)**2 + 3.0174314126331765E-4*LOG(satratln)**3   &
458   - 2.2673492408841294E-6*tln*LOG(satratln)**3 - 4.3948464567032377E-3*LOG(rhoaln)        &
459   + 5.3305314722492146E-5*tln*LOG(rhoaln)
460  x_i=  7.9036365428891719E-1 - 2.8414059650092153E-3*tli + 1.4976802556584141E-2*LOG(satratli)  &
461   - 2.4511581740839115E-4*tli*LOG(satratli) + 3.4319869471066424E-3 *LOG(satratli)**2     &
462   - 2.8799393617748428E-5*tli*LOG(satratli)**2 + 3.0174314126331765E-4*LOG(satratli)**3   &
463   - 2.2673492408841294E-6*tli*LOG(satratli)**3 - 4.3948464567032377E-3*LOG(rhoali)        &
464   + 5.3305314722492146E-5*tli*LOG(rhoali)
465       
466  x_n=MIN(MAX(x_n,1.E-30),1.)
467  x_i=MIN(MAX(x_i,1.E-30),1.)
468 
469  !Neutral nucleation
470 
471  !Kinetic limit check
472  IF (satratln >= 1.E-2 .AND. satratln <= 1.) THEN
473     kinrhotresn=EXP(7.8920778706888086E+1 + 7.3665492897447082*satratln - 1.2420166571163805E+4/tln &
474   + (-6.1831234251470971E+2*satratln)/tln - 2.4501159970109945E-2*tln                      &
475   -1.3463066443605762E-2*satratln*tln + 8.3736373989909194E-06*tln**2                      &
476   -1.4673887785408892*LOG(satratln) + (-3.2141890006517094E+1*LOG(satratln))/tln           &
477   + 2.7137429081917556E-3*tln*LOG(satratln)) !1/cm3
478     IF (kinrhotresn<rhoaln) kinetic_n=.TRUE.
479  ENDIF
480
481  IF (satratln >= 1.E-4  .AND. satratln < 1.E-2) THEN
482     kinrhotresn=EXP(7.9074383049843647E+1 - 2.8746005462158347E+1*satratln - 1.2070272068458380E+4/tln &
483   + (-5.9205040320056632E+3*satratln)/tln - 2.4800372593452726E-2*tln                         &
484   -4.3983007681295948E-2*satratln*tln + 2.5943854791342071E-5*tln**2                          &
485   -2.3141363245211317*LOG(satratln) + (9.9186787997857735E+1*LOG(satratln))/tln               &
486   + 5.6819382556144681E-3*tln*LOG(satratln)) !1/cm3
487     IF (kinrhotresn<rhoaln) kinetic_n=.TRUE.
488  ENDIF
489
490  IF (satratln >= 5.E-6  .AND. satratln < 1.E-4) THEN
491     kinrhotresn=EXP(8.5599712000361677E+1 + 2.7335119660796581E+3*satratln - 1.1842350246291651E+4/tln &
492   + (-1.2439843468881438E+6*satratln)/tln - 5.4536964974944230E-2*tln                         &
493   + 5.0886987425326087*satratln*tln + 7.1964722655507067E-5*tln**2                            &
494   -2.4472627526306372*LOG(satratln) + (1.7561478001423779E+2*LOG(satratln))/tln               &
495   + 6.2640132818141811E-3*tln*LOG(satratln)) !1/cm3
496     IF (kinrhotresn<rhoaln) kinetic_n=.TRUE.
497  ENDIF
498 
499  IF (kinetic_n) THEN   
500     ! Dimer formation rate
501     jnuc_n=1.E6*(2.*0.3E-9)**2.*SQRT(8.*RPI*RKBOL*(1./mH2SO4mol+1./mH2SO4mol))/2.*SQRT(t)*rhoa**2.
502     ntot_n=1. !set to 1
503     na_n=1.   ! The critical cluster contains one molecules but the produced cluster contains 2 molecules
504     x_n=na_n/ntot_n  ! so also set this to 1
505     rc_n=0.3E-9
506  ELSE
507     jnuc_n= 2.1361182605986115E-1 + 3.3827029855551838*tln -3.2423555796175563E-2*tln**2 +        &
508    7.0120069477221989E-5*tln**3 +8.0286874752695141/x_n +                                &
509    (-2.6939840579762231E-1)*LOG(satratln) +1.6079879299099518*tln*LOG(satratln) +        &
510    (-1.9667486968141933E-2)*tln**2*LOG(satratln) +                                       &
511    5.5244755979770844E-5*tln**3*LOG(satratln) + (7.8884704837892468*LOG(satratln))/x_n + &
512    4.6374659198909596*LOG(satratln)**2 - 8.2002809894792153E-2*tln*LOG(satratln)**2 +    &
513    8.5077424451172196E-4*tln**2*LOG(satratln)**2 +                                       &
514    (-2.6518510168987462E-6)*tln**3*LOG(satratln)**2 +                                    &
515    (-1.4625482500575278*LOG(satratln)**2)/x_n - 5.2413002989192037E-1*LOG(satratln)**3 + &
516    5.2755117653715865E-3*tln*LOG(satratln)**3 +                                          &
517    (-2.9491061332113830E-6)*tln**2*LOG(satratln)**3 +                                    &
518    (-2.4815454194486752E-8)*tln**3*LOG(satratln)**3 +                                    &
519    (-5.2663760117394626E-2*LOG(satratln)**3)/x_n +                                       &
520    1.6496664658266762*LOG(rhoaln) +                                                      &
521    (-8.0809397859218401E-1)*tln*LOG(rhoaln) +                                            &
522    8.9302927091946642E-3*tln**2*LOG(rhoaln) +                                            &
523    (-1.9583649496497497E-5)*tln**3*LOG(rhoaln) +                                         &
524    (-8.9505572676891685*LOG(rhoaln))/x_n +                                               &
525    (-3.0025283601622881E+1)*LOG(satratln)*LOG(rhoaln) +                                  &
526    3.0783365644763633E-1*tln*LOG(satratln)*LOG(rhoaln) +                                 &
527    (-7.4521756337984706E-4)*tln**2*LOG(satratln)*LOG(rhoaln) +                           &
528    (-5.7651433870681853E-7)*tln**3*LOG(satratln)*LOG(rhoaln) +                           &
529    (1.2872868529673207*LOG(satratln)*LOG(rhoaln))/x_n +                                  &
530    (-6.1739867501526535E-1)*LOG(satratln)**2*LOG(rhoaln) +                               &
531    7.2347385705333975E-3*tln*LOG(satratln)**2*LOG(rhoaln) +                              &
532    (-3.0640494530822439E-5)*tln**2*LOG(satratln)**2*LOG(rhoaln) +                        &
533    6.5944609194346214E-8*tln**3*LOG(satratln)**2*LOG(rhoaln) +                           &
534    (-2.8681650332461055E-2*LOG(satratln)**2*LOG(rhoaln))/x_n +                           &
535    6.5213802375160306*LOG(rhoaln)**2 +                                                   &
536    (-4.7907162004793016E-2)*tln*LOG(rhoaln)**2 +                                         &
537    (-1.0727890114215117E-4)*tln**2*LOG(rhoaln)**2 +                                      &
538    5.6401818280534507E-7*tln**3*LOG(rhoaln)**2 +                                         &
539    (5.4113070888923009E-1*LOG(rhoaln)**2)/x_n +                                          &
540    5.2062808476476330E-1*LOG(satratln)*LOG(rhoaln)**2 +                                  &
541    (-6.0696882500824584E-3)*tln*LOG(satratln)*LOG(rhoaln)**2 +                           &
542    2.3851383302608477E-5*tln**2*LOG(satratln)*LOG(rhoaln)**2 +                           &
543    (-1.5243837103067096E-8)*tln**3*LOG(satratln)*LOG(rhoaln)**2 +                        &
544    (-5.6543192378015687E-2*LOG(satratln)*LOG(rhoaln)**2)/x_n +                           &
545    (-1.1630806410696815E-1)*LOG(rhoaln)**3 +                                             &
546    1.3806404273119610E-3*tln*LOG(rhoaln)**3 +                                            &
547    (-2.0199865087650833E-6)*tln**2*LOG(rhoaln)**3 +                                      &
548    (-3.0200284885763192E-9)*tln**3*LOG(rhoaln)**3 +                                      &
549    (-6.9425267104126316E-3*LOG(rhoaln)**3)/x_n
550     jnuc_n=EXP(jnuc_n)
551     
552     ntot_n =-3.5863435141979573E-3 - 1.0098670235841110E-1*tln + 8.9741268319259721E-4*tln**2 - 1.4855098605195757E-6*tln**3  &
553     - 1.2080330016937095E-1/x_n + 1.1902674923928015E-3*LOG(satratln) - 1.9211358507172177E-2*tln*LOG(satratln) +    &
554     2.4648094311204255E-4*tln**2*LOG(satratln) - 7.5641448594711666E-7*tln**3*LOG(satratln) +                        &
555     (-2.0668639384228818E-02*LOG(satratln))/x_n - 3.7593072011595188E-2*LOG(satratln)**2 +                           &
556     9.0993182774415718E-4 *tln*LOG(satratln)**2 +                                                                    &
557     (-9.5698412164297149E-6)*tln**2*LOG(satratln)**2 + 3.7163166416110421E-8*tln**3*LOG(satratln)**2 +               &
558     (1.1026579525210847E-2*LOG(satratln)**2)/x_n + 1.1530844115561925E-2 *LOG(satratln)**3 +                         &
559     (-1.8083253906466668E-4)*tln*LOG(satratln)**3 + 8.0213604053330654E-7*tln**2*LOG(satratln)**3 +                  &
560     (-8.5797885383051337E-10)*tln**3*LOG(satratln)**3 + (1.0243693899717402E-3*LOG(satratln)**3)/x_n +               &
561     (-1.7248695296299649E-2)*LOG(rhoaln) + 1.1294004162437157E-2*tln*LOG(rhoaln) +                                   &
562     (-1.2283640163189278E-4)*tln**2*LOG(rhoaln) + 2.7391732258259009E-7*tln**3*LOG(rhoaln) +                         &
563     (6.8505583974029602E-2*LOG(rhoaln))/x_n +2.9750968179523635E-1*LOG(satratln)*LOG(rhoaln) +                       &
564     (-3.6681154503992296E-3)*tln*LOG(satratln)*LOG(rhoaln) + 1.0636473034653114E-5*tln**2*LOG(satratln)*LOG(rhoaln)+ &
565     5.8687098466515866E-9*tln**3*LOG(satratln)*LOG(rhoaln) + (-5.2028866094191509E-3*LOG(satratln)*LOG(rhoaln))/x_n+ &
566     7.6971988880587231E-4*LOG(satratln)**2*LOG(rhoaln) - 2.4605575820433763E-5*tln*LOG(satratln)**2*LOG(rhoaln) +    &
567     2.3818484400893008E-7*tln**2*LOG(satratln)**2*LOG(rhoaln) +                                                      &
568     (-8.8474102392445200E-10)*tln**3*LOG(satratln)**2*LOG(rhoaln) +                                                  &
569     (-1.6640566678168968E-4*LOG(satratln)**2*LOG(rhoaln))/x_n - 7.7390093776705471E-2*LOG(rhoaln)**2 +               &
570     5.8220163188828482E-4*tln*LOG(rhoaln)**2 + 1.2291679321523287E-6*tln**2*LOG(rhoaln)**2 +                         &
571     (-7.4690997508075749E-9)*tln**3*LOG(rhoaln)**2 + (-5.6357941220497648E-3*LOG(rhoaln)**2)/x_n +                   &
572     (-4.7170109625089768E-3)*LOG(satratln)*LOG(rhoaln)**2 + 6.9828868534370193E-5*tln*LOG(satratln)*LOG(rhoaln)**2 + &
573     (-3.1738912157036403E-7)*tln**2*LOG(satratln)*LOG(rhoaln)**2 +                                                   &
574     2.3975538706787416E-10*tln**3*LOG(satratln)*LOG(rhoaln)**2 +                                                     &
575     (4.2304213386288567E-4*LOG(satratln)*LOG(rhoaln)**2)/x_n + 1.3696520973423231E-3*LOG(rhoaln)**3 +                &
576     (-1.6863387574788199E-5)*tln*LOG(rhoaln)**3 + 2.7959499278844516E-8*tln**2*LOG(rhoaln)**3 +                      &
577     3.9423927013227455E-11*tln**3*LOG(rhoaln)**3 + (8.6136359966337272E-5*LOG(rhoaln)**3)/x_n
578     ntot_n=EXP(ntot_n)
579     
580     rc_n=EXP(-22.378268374023630+0.44462953606125100*x_n+0.33499495707849131*LOG(ntot_n)) !in meters
581     
582     na_n=x_n*ntot_n
583     IF (na_n < 1.) THEN
584        na_n=1.0
585     ENDIF
586  ENDIF
587 
588  ! Set the neutral nucleation rate to 0.0 if less than 1.0E-7     
589  IF (jnuc_n<1.E-7) THEN
590     jnuc_n=0.0
591  ENDIF
592 
593  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
594 
595  ! Ion-induced nucleation parameterization
596 
597  IF (ipr>0.0) THEN ! if the ion production rate is above zero
598     
599     ! Calculate the ion induced nucleation rate wrt. concentration of 1 ion/cm3
600     
601     kinrhotresi = 5.3742280876674478E1  - 6.6837931590012266E-3 *LOG(satratli)**(-2)                                     &
602   - 1.0142598385422842E-01 * LOG(satratli)**(-1) - 6.4170597272606873E+00 * LOG(satratli)                       &
603   - 6.4315798914824518E-01 * LOG(satratli)**2 - 2.4428391714772721E-02 * LOG(satratli)**3                       &
604   - 3.5356658734539019E-04 * LOG(satratli)**4 + 2.5400015099140506E-05 * tli * LOG(satratli)**(-2)              &
605   - 2.7928900816637790E-04 * tli * LOG(satratli)**(-1) + 4.4108573484923690E-02 * tli * LOG(satratli)           &
606   + 6.3943789012475532E-03 * tli * LOG(satratli)**(2) + 2.3164296174966580E-04 * tli * LOG(satratli)**(3)       &
607   + 3.0372070669934950E-06 * tli * LOG(satratli)**4 + 3.8255873977423475E-06 * tli**2 * LOG(satratli)**(-1)     &
608   - 1.2344793083561629E-04 * tli**2 * LOG(satratli) - 1.7959048869810192E-05 * tli**2 * LOG(satratli)**(2)      &
609   - 3.2165622558722767E-07 * tli**2 * LOG(satratli)**3 - 4.7136923780988659E-09 * tli**3 * LOG(satratli)**(-1)  &
610   + 1.1873317184482216E-07 * tli**3 * LOG(satratli) + 1.5685860354866621E-08 * tli**3 * LOG(satratli)**2        &
611   - 1.4329645891059557E+04 * tli**(-1) + 1.3842599842575321E-01 * tli                                           &
612   - 4.1376265912842938E-04 * tli**(2) + 3.9147639775826004E-07 * tli**3
613     
614     kinrhotresi=EXP(kinrhotresi) !1/cm3
615     
616     IF (kinrhotresi<rhoali) kinetic_i=.TRUE.
617     
618     IF (kinetic_i) THEN   
619        jnuc_i1=1.0E6*(0.3E-9 + 0.487E-9)**2.*SQRT(8.*RPI*RKBOL*(1./mH2SO4mol+1./mH2SO4mol))*  &
620    SQRT(tli)*rhoali !1/cm3s
621        ntot_i=1. !set to 1
622        na_i=1.
623        x_i=na_i/ntot_i  ! so also set this to 1
624        rc_i=0.487E-9
625     ELSE
626        jnuc_i1 = 3.0108954259038608E+01+tli*6.1176722090512577E+01+(tli**2)*8.7240333618891663E-01+(tli**3)*                &
627   (-4.6191788649375719E-03)+(tli**(-1))*8.3537059107024481E-01 +                                                &
628   (1.5028549216690628E+01+tli*(-1.9310989753720623E-01)+(tli**2)*8.0155514634860480E-04+(tli**3)*               &
629   (-1.0832730707799128E-06)+(tli**(-1))*1.7577660457989019)*(LOG(satratli)**(-2)) +                             &
630   (-2.0487870170216488E-01 + tli * 1.3263949252910405E-03 + (tli**2) *(-8.4195688402450274E-06) +               &
631   (tli**3)*1.6154895940993287E-08 + (tli**(-1))*3.8734212545203874E+01) * (LOG(satratli)**(-2)*LOG(rhoali)) +   &
632   (1.4955918863858371 + tli * 9.2290004245522454E+01 + (tli**2) *(-8.9006965195392618E-01) +                    &
633   (tli**3) * 2.2319123411013099E-03 + (tli**(-1)) * 4.0180079996840852E-03) *                                   &
634   (LOG(satratli)**(-1) * LOG(rhoali)**(-1)) +                                                                   &
635   (7.9018031228561085 + tli *(-1.1649433968658949E+01) + (tli**2) * 1.1400827854910951E-01 +                    &
636   (tli**3) *(-3.1941526492127755E-04) + (tli**(-1)) *(-3.7662115740271446E-01)) * (LOG(satratli)**(-1)) +       &
637   (1.5725237111225979E+02 + tli *(-1.0051649979836277) + (tli**2) * 1.1866484014507624E-03 +                    &
638   (tli**3) * 7.3557614998540389E-06 + (tli**(-1)) * 2.6270197023115189) * (LOG(satratli)**(-1) * LOG(rhoali)) + &
639   (-1.6973840122470968E+01 + tli * 1.1258423691432135E-01 + (tli**2) *(-2.9850139351463793E-04) + (tli**3) *    &
640   1.4301286324827064E-07 + (tli**(-1)) * 1.3163389235253725E+01) * (LOG(satratli)**(-1) * LOG(rhoali)**2) +     &
641   (-1.0399591631839757 + tli * 2.7022055588257691E-03 + (tli**2) *(-2.1507467231330936E-06) + (tli**3) *        &
642   3.8059489037584171E-10 + (tli**(-1)) * 1.5000492788553410E+02) * (LOG(satratli)**(-1) * LOG(rhoali)**3) +     &
643   (1.2250990965305315 + tli * 3.0495946490079444E+01 + (tli**2) * 2.1051563135187106E+01 + (tli**3) *           &
644   (-8.2200682916580878E-02) + (tli**(-1)) * 2.9965871386685029E-02) * (LOG(rhoali)**(-2)) +                     &
645   (4.8281605955680433 + tli * 1.7346551710836445E+02 + (tli**2) *(-1.0113602140796010E+01) + (tli**3) *         &
646   3.7482518458685089E-02 + (tli**(-1)) *(-1.4449998158558205E-01)) * (LOG(rhoali)**(-1)) +                      &
647   (2.3399230964451237E+02 + tli *(-2.3099267235261948E+01) + (tli**2) * 8.0122962140916354E-02 +                &
648   (tli**3) * 6.1542576994557088E-05 + (tli**(-1)) * 5.3718413254843007) * (LOG(rhoali)) +                       &
649   (1.0299715519499360E+02 + tli *(-6.4663357203364136E-02) + (tli**2) *(-2.0487150565050316E-03) +              &
650   (tli**3) * 8.7935289055530897E-07 + (tli**(-1)) * 3.6013204601215229E+01) * (LOG(rhoali)**2) +                &
651   (-3.5452115439584042 + tli * 1.7083445731159330E-02 + (tli**2) *(-1.2552625290862626E-05) + (tli**3) *        &
652   1.2968447449182847E-09 + (tli**(-1)) * 1.5748687512056560E+02) * (LOG(rhoali)**3) +                           &
653   (2.2338490119517975 + tli * 1.0229410216045540E+02 + (tli**2) *(-3.2103611955174052) + (tli**3) *             &
654   1.3397152304977591E-02 + (tli**(-1)) *(-2.4155187776460030E-02)) * (LOG(satratli)* LOG(rhoali)**(-2)) +       &
655   (3.7592282990713963 + tli *(-1.5257988769009816E+02) + (tli**2) * 2.6113805420558802 + (tli**3) *             &
656   (-9.0380721653694363E-03) + (tli**(-1)) *(-1.3974197138171082E-01)) * (LOG(satratli)* LOG(rhoali)**(-1)) +    &
657   (1.8293600730573988E+01 + tli * 1.8344728606002992E+01 + (tli**2) *(-4.0063363221106751E-01) + (tli**3)       &
658   * 1.4842749371258522E-03 + (tli**(-1)) * 1.1848846003282287) * (LOG(satratli)) +                              &
659   (-1.7634531623032314E+02 + tli * 4.9011762441271278 + (tli**2) *(-1.3195821562746339E-02) + (tli**3) *        &
660   (-2.8668619526430859E-05) + (tli**(-1)) *(-2.9823396976393551E-01)) * (LOG(satratli)* LOG(rhoali)) +          &
661   (-3.2944043694275727E+01 + tli * 1.2517571921051887E-01 + (tli**2) * 8.3239769771186714E-05 + (tli**3) *      &
662   2.8191859341519507E-07 + (tli**(-1)) *(-2.7352880736682319E+01)) * (LOG(satratli)* LOG(rhoali)**2) +          &
663   (-1.1451811137553243 + tli * 2.0625997485732494E-03 + (tli**2) *(-3.4225389469233624E-06) + (tli**3) *        &
664   4.4437613496984567E-10 + (tli**(-1)) * 1.8666644332606754E+02) * (LOG(satratli)* LOG(rhoali)**3) +            &
665   (3.2270897099493567E+01 + tli * 7.7898447327513687E-01 + (tli**2) *(-6.5662738484679626E-03) + (tli**3) *     &
666   3.7899330796456790E-06 + (tli**(-1)) * 7.1106427501756542E-01) * (LOG(satratli)**2 * LOG(rhoali)**(-1)) +     &
667   (-2.8901906781697811E+01 + tli *(-1.5356398793054860) + (tli**2) * 1.9267271774384788E-02 + (tli**3) *        &
668   (-5.3886270475516162E-05) + (tli**(-1)) * 5.0490415975693426E-01) * (LOG(satratli)**2) +                      &
669   (3.3365683645733924E+01 + tli *(-3.6114561564894537E-01) + (tli**2) * 9.2977354471929262E-04 + (tli**3) *     &
670   1.9549769069511355E-07 + (tli**(-1)) *(-8.8865930095112855)) * (LOG(satratli)**2 * LOG(rhoali)) +             &
671   (2.4592563042806375 + tli *(-8.3227071743101084E-03) + (tli**2) * 8.2563338043447783E-06 + (tli**3) *         &
672   (-8.4374976698593496E-09) + (tli**(-1)) *(-2.0938173949893473E+02)) * (LOG(satratli)**2 * LOG(rhoali)**2) +   &
673   (4.4099823444352317E+01 + tli * 2.5915665826835252 + (tli**2) *(-1.6449091819482634E-02) + (tli**3) *         &
674   2.6797249816144721E-05 + (tli**(-1)) * 5.5045672663909995E-01)* satratli
675        jnuc_i1=EXP(jnuc_i1)
676       
677        ntot_i = ABS((-4.8324296064013375E+04 + tli * 5.0469120697428906E+02 + (tli**2) *(-1.1528940488496042E+00) +         &
678   (tli**(-1)) *(-8.6892744676239192E+02) + (tli**(3)) * 4.0030302028120469E-04) +                               &
679   (-6.7259105232039847E+03 + tli * 1.9197488157452008E+02 + (tli**2) *(-1.3602976930126354E+00) +               &
680   (tli**(-1)) *(-1.1212637938360332E+02) + (tli**(3)) * 2.8515597265933207E-03) *                               &
681   LOG(satratli)**(-2) * LOG(rhoali)**(-2) +                                                                     &
682   (2.6216455217763342E+02 + tli *(-2.3687553252750821E+00) + (tli**2) * 7.4074554767517521E-03 +                &
683   (tli**(-1)) *(-1.9213956820114927E+03) + (tli**(3)) *(-9.3839114856129453E-06)) * LOG(satratli)**(-2) +       &
684   (3.9652478944137344E+00 + tli * 1.2469375098256536E-02 + (tli**2) *(-9.9837754694045633E-05) + (tli**(-1)) *  &
685   (-5.1919499210175138E+02) + (tli**(3)) * 1.6489001324583862E-07) * LOG(satratli)**(-2) * LOG(rhoali) +        &
686   (2.4975714429096206E+02 + tli * 1.7107594562445172E+02 + (tli**2) *(-7.8988711365135289E-01) + (tli**(-1)) *  &
687   (-2.2243599782483177E+01) + (tli**(3)) *(-1.6291523004095427E-04)) * LOG(satratli)**(-1) * LOG(rhoali)**(-2) +&
688   (-8.9270715592533611E+02 + tli * 1.2053538883338946E+02 + (tli**2) *(-1.5490408828541018E+00) + (tli**(-1)) * &
689   (-1.1243275579419826E+01) + (tli**(3)) * 4.8053105606904655E-03) * LOG(satratli)**(-1) * LOG(rhoali)**(-1) +  &
690   (7.6426441642091631E+03 + tli *(-7.1785462414656578E+01) + (tli**2) * 2.3851864923199523E-01 + (tli**(-1)) *  &
691   8.5591775688708395E+01 + (tli**(3)) *(-3.7000473243342858E-04)) * LOG(satratli)**(-1) +                       &
692   (-5.1516826398607911E+01 + tli * 9.1385720811460558E-01 + (tli**2) *(-3.5477100262158974E-03) +               &
693   (tli**(-1)) * 2.7545544507625586E+03 + (tli**(3)) * 5.4708262093640928E-06) * LOG(satratli)**(-1) * LOG(rhoali) + &
694   (-3.0386767129196176E+02 + tli *(-1.1033438883583569E+04) + (tli**2) * 8.1296859732896067E+01 + (tli**(-1)) * &
695   1.2625883141097162E+01 + (tli**(3)) *(-1.2728497822219101E-01)) * LOG(rhoali)**(-2) +                         &
696   (-3.3763494256461472E+03 + tli * 3.1916579136391006E+03 + (tli**2) *(-2.7234339474441143E+01) + (tli**(-1)) * &
697   (-2.1897653262707397E+01) + (tli**(3)) * 5.1788505812259071E-02) * LOG(rhoali)**(-1) +                        &
698   (-1.8817843873687068E+03 + tli * 4.3038072285882070E+00 + (tli**2) * 6.6244087689671860E-03 + (tli**(-1)) *   &
699   (-2.7133073605696295E+03) + (tli**(3)) *(-1.7951557394285043E-05)) * LOG(rhoali) +                            &
700   (-1.7668827539244447E+02 + tli * 4.8160932330629913E-01 + (tli**2) *(-6.3133007671100293E-04) + (tli**(-1)) * &
701   2.5631774669873157E+04 + (tli**(3)) * 4.1534484127873519E-07) * LOG(rhoali)**(2) +                            &
702   (-1.6661835889222382E+03 + tli * 1.3708900504682877E+03 + (tli**2) *(-1.7919060052198969E+01) + (tli**(-1)) * &
703   (-3.5145029804436405E+01) + (tli**(3)) * 5.1047240947371224E-02) * LOG(satratli)* LOG(rhoali)**(-2) +         &
704   (1.0843549363030939E+04 + tli *(-7.3557073636139577E+01) + (tli**2) * 1.2054625131778862E+00 + (tli**(-1)) *  &
705   1.9358737917864391E+02 + (tli**(3)) *(-4.2871620775911338E-03)) * LOG(satratli)* LOG(rhoali)**(-1) +          &
706   (-2.4269802549752835E+03 + tli * 1.1348265061941714E+01 + (tli**2) *(-5.0430423939495157E-02) + (tli**(-1)) * &
707   2.3709874548950634E+03 + (tli**(3)) * 1.4091851828620244E-04) * LOG(satratli) +                               &
708   (5.2745372575251588E+02 + tli *(-2.6080675912627314E+00) + (tli**2) * 5.6902218056670145E-03 + (tli**(-1)) *  &
709   (-3.2149319482897838E+04) + (tli**(3)) *(-5.4121996056745853E-06)) * LOG(satratli)* LOG(rhoali) +             &
710   (-1.6401959518360403E+01 + tli * 2.4322962162439640E-01 + (tli**2) * 1.1744366627725344E-03 + (tli**(-1)) *   &
711   (-8.2694427518413195E+03) + (tli**(3)) *(-5.0028379203873102E-06))* LOG(satratli)**(2) +                      &
712   (-2.7556572017167782E+03 + tli * 4.9293344495058264E+01 + (tli**2) *(-2.6503456520676050E-01) + (tli**(-1)) * &
713   1.2130698030982167E+03 + (tli**(3)) * 4.3530610668042957E-04)* LOG(satratli)**2 * LOG(rhoali)**(-1) +         &
714   (-6.3419182228959192E+00 + tli * 4.0636212834605827E-02 + (tli**2) *(-1.0450112687842742E-04) + (tli**(-1)) * &
715   3.1035882189759656E+02 + (tli**(3)) * 9.4328418657873500E-08)* LOG(satratli)**(-3) +                          &
716   (3.0189213304689042E+03 + tli *(-2.3804654203861684E+01) + (tli**2) * 6.8113013411972942E-02 + (tli**(-1)) *  &
717   6.3112071081188913E+02 + (tli**(3)) *(-9.4460854261685723E-05))* (satratli) * LOG(rhoali) +                   &
718   (1.1924791930673702E+04 + tli *(-1.1973824959206000E+02) + (tli**2) * 1.6888713097971020E-01 + (tli**(-1)) *  &
719   1.8735938211539585E+02 + (tli**(3)) * 5.0974564680442852E-04)* (satratli) +                                   &
720   (3.6409071302482083E+01 + tli * 1.7919859306449623E-01 + (tli**2) *(-1.0020116255895206E-03) + (tli**(-1)) *  &
721   (-8.3521083354432303E+03) + (tli**(3)) * 1.5879900546795635E-06)* satratli * LOG(rhoali)**(2))
722         
723        rc_i = (-3.6318550637865524E-08 + tli * 2.1740704135789128E-09   + (tli**2) *                            &
724   (-8.5521429066506161E-12) + (tli**3) *(-9.3538647454573390E-15)) +                                &
725   (2.1366936839394922E-08 + tli *(-2.4087168827395623E-10) + (tli**2) * 8.7969869277074319E-13 +    &
726   (tli**3) *(-1.0294466881303291E-15))* LOG(satratli)**(-2) * LOG(rhoali)**(-1) +                   &
727   (-7.7804007761164303E-10 + tli * 1.0327058173517932E-11 + (tli**2) *(-4.2557697639692428E-14) +   &
728   (tli**3) * 5.4082507061618662E-17)* LOG(satratli)**(-2) +                                         &
729   (3.2628927397420860E-12 + tli *(-7.6475692919751066E-14) + (tli**2) * 4.1985816845259788E-16 +    &
730   (tli**3) *(-6.2281395889592719E-19))* LOG(satratli)**(-2) * LOG(rhoali) +                         &
731   (2.0442205540818555E-09 + tli * 4.0441858911249830E-08 + (tli**2) *(-3.3423487629482825E-10) +    &
732   (tli**3) * 6.8000404742985678E-13)* LOG(satratli)**(-1) * LOG(rhoali)**(-2) +                     &
733   (1.8381489183824627E-08 + tli *(-8.9853322951518919E-09) + (tli**2) * 7.5888799566036185E-11 +    &
734   (tli**3) *(-1.5823457864755549E-13))* LOG(satratli)**(-1) * LOG(rhoali)**(-1) +                   &
735   (1.1795760639695057E-07 + tli *(-8.1046722896375875E-10) + (tli**2) * 9.1868604369041857E-14 +    &
736   (tli**3) * 4.7882428237444610E-15)* LOG(satratli)**(-1) +                                         &
737   (-4.4028846582545952E-09 + tli * 4.6541269232626618E-11 + (tli**2) *(-1.1939929984285194E-13) +   &
738   (tli**3) * 2.3602037016614437E-17)* LOG(satratli)**(-1) * LOG(rhoali) +                           &
739   (2.7885056884209128E-11 + tli *(-4.5167129624119121E-13) + (tli**2) * 1.6558404997394422E-15 +    &
740   (tli**3) *(-1.2037336621218054E-18))* LOG(satratli)**(-1) * LOG(rhoali)**2 +                      &
741   (-2.3719627171699983E-09 + tli *(-1.5260127909292053E-07) + (tli**2) * 1.7177017944754134E-09 +   &
742   (tli**3) *(-4.7031737537526395E-12))* LOG(rhoali)**(-2) +                                         &
743   (-5.6946433724699646E-09 + tli * 8.4629788237081735E-09 + (tli**2) *(-1.7674135187061521E-10) +   &
744   (tli**3) * 6.6236547903091862E-13)* LOG(rhoali)**(-1) +                                           &
745   (-2.2808617930606012E-08 + tli * 1.4773376696847775E-10 + (tli**2) *(-1.3076953119957355E-13) +   &
746   (tli**3) * 2.3625301497914000E-16)* LOG(rhoali) +                                                 &
747   (1.4014269939947841E-10 + tli *(-2.3675117757377632E-12) + (tli**2) * 5.1514033966707879E-15 +    &
748   (tli**3) *(-4.8864233454747856E-18))* LOG(rhoali)**2 +                                            &
749   (6.5464943868885886E-11 + tli * 1.6494354816942769E-08 + (tli**2) *(-1.7480097393483653E-10) +    &
750   (tli**3) * 4.7460075628523984E-13)* LOG(satratli)* LOG(rhoali)**(-2) +                            &
751   (8.4737893183927871E-09 + tli *(-6.0243327445597118E-09) + (tli**2) * 5.8766070529814883E-11 +    &
752   (tli**3) *(-1.4926748560042018E-13))* LOG(satratli)* LOG(rhoali)**(-1) +                          &
753   (1.0761964135701397E-07 + tli *(-1.0142496009071148E-09) + (tli**2) * 2.1337312466519190E-12 +    &
754   (tli**3) * 1.6376014957685404E-15)* LOG(satratli) +                                               &
755   (-3.5621571395968670E-09 + tli * 4.1175339587760905E-11 + (tli**2) *(-1.3535372357998504E-13) +   &
756   (tli**3) * 8.9334219536920720E-17)* LOG(satratli)* LOG(rhoali) +                                  &
757   (2.0700482083136289E-11 + tli *(-3.9238944562717421E-13) + (tli**2) * 1.5850961422040196E-15 +    &
758   (tli**3) *(-1.5336775610911665E-18))* LOG(satratli)* LOG(rhoali)**2 +                             &
759   (1.8524255464416206E-09 + tli *(-2.1959816152743264E-11) + (tli**2) *(-6.4478119501677012E-14) +  &
760   (tli**3) * 5.5135243833766056E-16)* LOG(satratli)**2 * LOG(rhoali)**(-1) +                        &
761   (1.9349488650922679E-09 + tli *(-2.2647295919976428E-11) + (tli**2) * 9.2917479748268751E-14 +    &
762   (tli**3) *(-1.2741959892173170E-16))* LOG(satratli)**2 +                                          &
763   (2.1484978031650972E-11 + tli *(-9.3976642475838013E-14) + (tli**2) *(-4.8892738002751923E-16) +  &
764   (tli**3) * 1.4676120441783832E-18)* LOG(satratli)**2 * LOG(rhoali) +                              &
765   (6.7565715216420310E-13 + tli *(-3.5421162549480807E-15) + (tli**2) *(-3.4201196868693569E-18) +  &
766   (tli**3) * 2.2260187650412392E-20)* LOG(satratli)**3 * LOG(rhoali)
767                   
768        na_i=x_i*ntot_i
769        IF (na_i < 1.) THEN
770           na_n=1.0
771        ENDIF
772     ENDIF
773   
774     jnuc_i=jnuc_i1
775     ! Ion loss rate (1/s)
776     xloss=csi+jnuc_i
777     
778     ! Recombination (here following Brasseur and Chatel, 1983)   
779     recomb=6.0E-8*SQRT(300./tli)+6.0E-26*airn*(300./tli)**4
780     
781     ! Small ion concentration in air (1/cm3) (following Dunne et al., 2016)
782     ! max function is to avoid n_i to go practically zero at very high J_ion
783     n_i=MAX(0.01,(SQRT(xloss**2.0+4.0*recomb*ipr)-xloss)/(2.0*recomb))
784     
785     ! Ion-induced nucleation rate
786     ! Min function is to ensure that max function above does not cause J_ion to overshoot
787     jnuc_i=MIN(ipr,n_i*jnuc_i1)
788     ! Set the ion-induced nucleation rate to 0.0 if less than 1.0E-7     
789     IF (jnuc_i<1.E-7) THEN
790        jnuc_i=0.0
791     ENDIF
792
793  ENDIF
794
795!--conversion from double precision to float in case the model is run in single precision
796  jnuc_n_real   = jnuc_n
797  ntot_n_real   = ntot_n
798  jnuc_i_real   = jnuc_i
799  ntot_i_real   = ntot_i
800  x_n_real      = x_n
801  x_i_real      = x_i
802  na_n_real     = na_n
803  na_i_real     = na_i
804  rc_n_real     = rc_n
805  rc_i_real     = rc_i
806  n_i_real      = n_i
807 
808END SUBROUTINE newbinapara
809
810END MODULE nucleation_tstep_mod
Note: See TracBrowser for help on using the repository browser.