Changeset 3534 for LMDZ6/trunk/libf


Ignore:
Timestamp:
Jun 10, 2019, 9:19:46 PM (6 years ago)
Author:
oboucher
Message:

Introducing new nucleation routine from Anni
(not quite finished but it's getting in)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/phylmd/StratAer/nucleation_tstep_mod.F90

    r3527 r3534  
    1010  USE aerophys
    1111  USE infotrac
    12   USE YOMCST, ONLY : RPI, RD
     12  USE YOMCST, ONLY : RPI, RD, RKBOL
    1313
    1414  IMPLICIT NONE
    1515
    16   ! input variables
     16  ! input variable
     17  LOGICAL, PARAMETER :: flag_new_nucl=.TRUE.
    1718  REAL rhoa !H2SO4 number density [molecules/cm3]
    1819  REAL t_seri
    … …  
    2122  REAL a_xm, b_xm, c_xm
    2223
    23   ! output variables
     24  ! output variable
    2425  REAL nucl_rate
    2526  REAL ntot ! total number of molecules in the critical cluster
    26   REAL x    ! molefraction of H2SO4 in the critical cluster
    27 
    28   ! local variables
    29   REAL, PARAMETER                               :: k_B=1.3806E-23  ! Boltzmann constant [J/K]
     27  REAL x    ! mole fraction of H2SO4 in the critical cluster
     28
     29  ! local variable
    3030  REAL                                          :: jnuc !nucleation rate in 1/cm3s (10^-7-10^10 1/cm3s)
    3131  REAL                                          :: rc   !radius of the critical cluster in nm
    … …  
    3333
    3434  ! call nucleation routine
    35   CALL binapara(t_seri,rh,rhoa,jnuc,x,ntot,rc)
     35!  IF (.NOT.flag_new_nucl) THEN
     36    CALL binapara(t_seri,rh,rhoa,jnuc,x,ntot,rc)
     37!  ELSE
     38!    CALL newbinapara(t_seri,rh,rhoa,jnuc,x,ntot,rc)
     39!  ENDIF
    3640
    3741  IF (ntot < 4.0) THEN
    3842    !set jnuc to collision rate of two H2SO4 molecules (following personal communication of Ulrike Niemeier and Hanna Vehkamäki)
    39     VH2SO4mol=mH2SO4mol/(1.e-3*(a_xm+t_seri*(b_xm+t_seri*c_xm))) !cm3
    40     jnuc = rhoa**2. *(3./4.*RPI)**(1./6.) *(12.*k_B*t_seri/mH2SO4mol)**0.5 &
     43    VH2SO4mol=mH2SO4mol/(1.E-3*(a_xm+t_seri*(b_xm+t_seri*c_xm))) !cm3
     44    jnuc = rhoa**2. *(3./4.*RPI)**(1./6.) *(12.*RKBOL*t_seri/mH2SO4mol)**0.5 &
    4145         & *100.*(2.*VH2SO4mol**(1./3.))**2. !1/(cm3s)
    4246    ntot=2.0
    … …  
    5862  IMPLICIT NONE
    5963
    60   ! input variables
     64  ! input variable
    6165  REAL nucl_rate
    6266  REAL ntot ! total number of molecules in the critical cluster
    63   REAL x    ! molefraction of H2SO4 in the critical cluster
     67  REAL x    ! mole raction of H2SO4 in the critical cluster
    6468  REAL dt
    6569  REAL Vbin(nbtr_bin)
    6670
    67   ! output variables
     71  ! output variable
    6872  REAL tr_seri(nbtr)
    6973
    70   ! local variables
     74  ! local variable
    7175  INTEGER k
    7276  REAL Vnew
    … …  
    8084  DO k=1, nbtr_bin
    8185  ! CK 20160531: bug fix for first bin
    82     IF (k.LE.(nbtr_bin-1)) THEN
     86    IF (k.LE.nbtr_bin-1) THEN
    8387      IF (Vbin(k).LE.Vnew.AND.Vnew.LT.Vbin(k+1)) THEN
    8488        ff(k)= Vbin(k)/Vnew*(Vbin(k+1)-Vnew)/(Vbin(k+1)-Vbin(k))
    … …  
    113117SUBROUTINE binapara(pt,prh,rhoa_in,jnuc,x,ntot,rc)
    114118
     119
    115120  !    Fortran 90 subroutine binapara
    116121  !
    … …  
    134139
    135140  REAL :: pt,t     !temperature in K (190.15-300.15K)
    136   REAL :: prh,rh    !saturatio ratio of water (0.0001-1)
    137   REAL,intent(in) :: rhoa_in    !sulfuric acid concentration in 1/cm3 (10^4-10^11 1/cm3)
    138   REAL,intent(out) :: jnuc    !nucleation rate in 1/cm3s (10^-7-10^10 1/cm3s)
    139   REAL,intent(out) :: ntot !total number of molecules in the critical cluster (ntot>4)
    140   REAL,intent(out) :: x    ! molefraction of H2SO4 in the critical cluster       
    141   REAL,intent(out) :: rc    !radius of the critical cluster in nm     
    142   REAL :: rhotres    ! treshold concentration of h2so4 (1/cm^3)
     141  REAL :: prh,rh    !saturation ratio of water (0.0001-1)
     142  REAL,INTENT(IN) :: rhoa_in    !sulfuric acid concentration in 1/cm3 (10^4-10^11 1/cm3)
     143  REAL,INTENT(OUT) :: jnuc    !nucleation rate in 1/cm3s (10^-7-10^10 1/cm3s)
     144  REAL,INTENT(OUT) :: ntot !total number of molecules in the critical cluster (ntot>4)
     145  REAL,INTENT(OUT) :: x    ! mole fraction of H2SO4 in the critical cluster       
     146  REAL,INTENT(OUT) :: rc    !radius of the critical cluster in nm     
     147  REAL :: rhotres    ! threshold concentration of h2so4 (1/cm^3)
    143148                     ! which produces nucleation rate   1/(cm^3 s) as a function of rh and t
    144149  REAL rhoa
    145150
    146 ! CK: use intermediate variables to avoid overwriting
     151! CK: use intermediate variable to avoid overwriting
    147152  t=pt
    148153  rh=prh
    … …  
    150155
    151156  IF (t < 190.15) THEN
    152 !     print *,'Warning (ilon=',ilon,'ilev=',ilev,'): temperature < 190.15 K, using 190.15 K (T=',t,'K)'
     157!     print *,'Warning (ilon=',ilon,'ilat=',ilat,'): temperature < 190.15 K, using 190.15 K (T=',t,'K)'
    153158     t=190.15
    154159  ENDIF
    155160
    156161  IF (t > 300.15) THEN
    157 !     print *,'Warning (ilon=',ilon,'ilev=',ilev,'): temperature > 300.15 K, using 300.15 K'
     162!     print *,'Warning (ilon=',ilon,'ilat=',ilat,'): temperature > 300.15 K, using 300.15 K'
    158163     t=300.15
    159164  ENDIF
    160165
    161166  IF (rh < 0.0001) THEN
    162 !     print *,'Warning (ilon=',ilon,'ilev=',ilev,'): saturation ratio of water < 0.0001, using 0.0001'
     167!     print *,'Warning (ilon=',ilon,'ilat=',ilat,'): saturation ratio of water < 0.0001, using 0.0001'
    163168     rh=0.0001
    164169  ENDIF
    165170
    166171  IF (rh > 1.0) THEN
    167 !     print *,'Warning (ilon=',ilon,'ilev=',ilev,'): saturation ratio of water > 1 using 1'
     172!     print *,'Warning (ilon=',ilon,'ilat=',ilat,'): saturation ratio of water > 1 using 1'
    168173!     print *, 'rh=',rh
    169174     rh=1.0
    … …  
    171176
    172177  IF (rhoa < 1.e4) THEN
    173 !     print *,'Warning (ilon=',ilon,'ilev=',ilev,'): sulfuric acid concentration < 1e4 1/cm3, using 1e4 1/cm3'
     178!     print *,'Warning (ilon=',ilon,'ilat=',ilat,'): sulfuric acid concentration < 1e4 1/cm3, using 1e4 1/cm3'
    174179     rhoa=1.e4
    175180  ENDIF
    176181
    177182  IF (rhoa > 1.e11) THEN
    178 !     print *,'Warning (ilon=',ilon,'ilev=',ilev,'): sulfuric acid concentration > 1e11 1/cm3, using 1e11 1/cm3'
     183!     print *,'Warning (ilon=',ilon,'ilat=',ilat,'): sulfuric acid concentration > 1e11 1/cm3, using 1e11 1/cm3'
    179184     rhoa=1.e11
    180185  ENDIF
    181186
    182   x=  0.7409967177282139 - 0.002663785665140117*t + 0.002010478847383187*Log(rh)  &
    183        & - 0.0001832894131464668*t*Log(rh) + 0.001574072538464286*Log(rh)**2      &
    184        & - 0.00001790589121766952*t*Log(rh)**2 + 0.0001844027436573778*Log(rh)**3 &
    185        & -  1.503452308794887e-6*t*Log(rh)**3 - 0.003499978417957668*Log(rhoa)    &
    186        & + 0.0000504021689382576*t*Log(rhoa)
     187  x=  0.7409967177282139 - 0.002663785665140117*t + 0.002010478847383187*LOG(rh)  &
     188       & - 0.0001832894131464668*t*LOG(rh) + 0.001574072538464286*LOG(rh)**2      &
     189       & - 0.00001790589121766952*t*LOG(rh)**2 + 0.0001844027436573778*LOG(rh)**3 &
     190       & -  1.503452308794887E-6*t*LOG(rh)**3 - 0.003499978417957668*LOG(rhoa)    &
     191       & + 0.0000504021689382576*t*LOG(rhoa)
    187192
    188193  jnuc= 0.1430901615568665 + 2.219563673425199*t - 0.02739106114964264*t**2 +  &
    189194       &  0.00007228107239317088*t**3 + 5.91822263375044/x +                   &
    190        &  0.1174886643003278*Log(rh) + 0.4625315047693772*t*Log(rh) -          &
    191        &  0.01180591129059253*t**2*Log(rh) +                                   &
    192        &  0.0000404196487152575*t**3*Log(rh) + (15.79628615047088*Log(rh))/x - &
    193        &  0.215553951893509*Log(rh)**2 - 0.0810269192332194*t*Log(rh)**2 +     &
    194        &  0.001435808434184642*t**2*Log(rh)**2 -                               &
    195        &  4.775796947178588e-6*t**3*Log(rh)**2 -                               &
    196        &  (2.912974063702185*Log(rh)**2)/x - 3.588557942822751*Log(rh)**3 +    &
    197        &  0.04950795302831703*t*Log(rh)**3 -                                   &
    198        &  0.0002138195118737068*t**2*Log(rh)**3 +                              &
    199        &  3.108005107949533e-7*t**3*Log(rh)**3 -                               &
    200        &  (0.02933332747098296*Log(rh)**3)/x +                                 &
    201        &  1.145983818561277*Log(rhoa) -                                        &
    202        &  0.6007956227856778*t*Log(rhoa) +                                     &
    203        &  0.00864244733283759*t**2*Log(rhoa) -                                 &
    204        &  0.00002289467254710888*t**3*Log(rhoa) -                              &
    205        &  (8.44984513869014*Log(rhoa))/x +                                     &
    206        &  2.158548369286559*Log(rh)*Log(rhoa) +                                &
    207        &  0.0808121412840917*t*Log(rh)*Log(rhoa) -                             &
    208        &  0.0004073815255395214*t**2*Log(rh)*Log(rhoa) -                       &
    209        &  4.019572560156515e-7*t**3*Log(rh)*Log(rhoa) +                        &
    210        &  (0.7213255852557236*Log(rh)*Log(rhoa))/x +                           &
    211        &  1.62409850488771*Log(rh)**2*Log(rhoa) -                              &
    212        &  0.01601062035325362*t*Log(rh)**2*Log(rhoa) +                         &
    213        &  0.00003771238979714162*t**2*Log(rh)**2*Log(rhoa) +                   &
    214        &  3.217942606371182e-8*t**3*Log(rh)**2*Log(rhoa) -                     &
    215        &  (0.01132550810022116*Log(rh)**2*Log(rhoa))/x +                       &
    216        &  9.71681713056504*Log(rhoa)**2 -                                      &
    217        &  0.1150478558347306*t*Log(rhoa)**2 +                                  &
    218        &  0.0001570982486038294*t**2*Log(rhoa)**2 +                            &
    219        &  4.009144680125015e-7*t**3*Log(rhoa)**2 +                             &
    220        &  (0.7118597859976135*Log(rhoa)**2)/x -                                &
    221        &  1.056105824379897*Log(rh)*Log(rhoa)**2 +                             &
    222        &  0.00903377584628419*t*Log(rh)*Log(rhoa)**2 -                         &
    223        &  0.00001984167387090606*t**2*Log(rh)*Log(rhoa)**2 +                   &
    224        &  2.460478196482179e-8*t**3*Log(rh)*Log(rhoa)**2 -                     &
    225        &  (0.05790872906645181*Log(rh)*Log(rhoa)**2)/x -                       &
    226        &  0.1487119673397459*Log(rhoa)**3 +                                    &
    227        &  0.002835082097822667*t*Log(rhoa)**3 -                                &
    228        &  9.24618825471694e-6*t**2*Log(rhoa)**3 +                              &
    229        &  5.004267665960894e-9*t**3*Log(rhoa)**3 -                             &
    230        &  (0.01270805101481648*Log(rhoa)**3)/x
    231   jnuc=exp(jnuc) !1/(cm3s)
    232 
    233   ntot =-0.002954125078716302 - 0.0976834264241286*t + 0.001024847927067835*t**2 - 2.186459697726116e-6*t**3 -    &
    234        &   0.1017165718716887/x - 0.002050640345231486*Log(rh) - 0.007585041382707174*t*Log(rh) +                 &
    235        &   0.0001926539658089536*t**2*Log(rh) - 6.70429719683894e-7*t**3*Log(rh) -                                &
    236        &   (0.2557744774673163*Log(rh))/x + 0.003223076552477191*Log(rh)**2 + 0.000852636632240633*t*Log(rh)**2 - &
    237        &   0.00001547571354871789*t**2*Log(rh)**2 + 5.666608424980593e-8*t**3*Log(rh)**2 +                        &
    238        &   (0.03384437400744206*Log(rh)**2)/x + 0.04743226764572505*Log(rh)**3 -                                  &
    239        &   0.0006251042204583412*t*Log(rh)**3 + 2.650663328519478e-6*t**2*Log(rh)**3 -                            &
    240        &   3.674710848763778e-9*t**3*Log(rh)**3 - (0.0002672510825259393*Log(rh)**3)/x -                          &
    241        &   0.01252108546759328*Log(rhoa) + 0.005806550506277202*t*Log(rhoa) -                                     &
    242        &   0.0001016735312443444*t**2*Log(rhoa) + 2.881946187214505e-7*t**3*Log(rhoa) +                           &
    243        &   (0.0942243379396279*Log(rhoa))/x - 0.0385459592773097*Log(rh)*Log(rhoa) -                              &
    244        &   0.0006723156277391984*t*Log(rh)*Log(rhoa) + 2.602884877659698e-6*t**2*Log(rh)*Log(rhoa) +              &
    245        &   1.194163699688297e-8*t**3*Log(rh)*Log(rhoa) - (0.00851515345806281*Log(rh)*Log(rhoa))/x -              &
    246        &   0.01837488495738111*Log(rh)**2*Log(rhoa) + 0.0001720723574407498*t*Log(rh)**2*Log(rhoa) -              &
    247        &   3.717657974086814e-7*t**2*Log(rh)**2*Log(rhoa) -                                                       &
    248        &   5.148746022615196e-10*t**3*Log(rh)**2*Log(rhoa) +                                                      &
    249        &   (0.0002686602132926594*Log(rh)**2*Log(rhoa))/x - 0.06199739728812199*Log(rhoa)**2 +                    &
    250        &   0.000906958053583576*t*Log(rhoa)**2 - 9.11727926129757e-7*t**2*Log(rhoa)**2 -                          &
    251        &   5.367963396508457e-9*t**3*Log(rhoa)**2 - (0.007742343393937707*Log(rhoa)**2)/x +                       &
    252        &   0.0121827103101659*Log(rh)*Log(rhoa)**2 - 0.0001066499571188091*t*Log(rh)*Log(rhoa)**2 +               &
    253        &   2.534598655067518e-7*t**2*Log(rh)*Log(rhoa)**2 -                                                       &
    254        &   3.635186504599571e-10*t**3*Log(rh)*Log(rhoa)**2 +                                                      &
    255        &   (0.0006100650851863252*Log(rh)*Log(rhoa)**2)/x + 0.0003201836700403512*Log(rhoa)**3 -                  &
    256        &   0.0000174761713262546*t*Log(rhoa)**3 + 6.065037668052182e-8*t**2*Log(rhoa)**3 -                        &
    257        &   1.421771723004557e-11*t**3*Log(rhoa)**3 + (0.0001357509859501723*Log(rhoa)**3)/x
    258   ntot=exp(ntot)
    259 
    260   rc=exp(-1.6524245+0.42316402*x+0.33466487*log(ntot)) !nm
    261 
    262   IF (jnuc < 1.e-7) THEN
    263 !     print *,'Warning (ilon=',ilon,'ilev=',ilev'): nucleation rate < 1e-7/cm3s, using 0.0/cm3s,'
     195       &  0.1174886643003278*LOG(rh) + 0.4625315047693772*t*LOG(rh) -          &
     196       &  0.01180591129059253*t**2*LOG(rh) +                                   &
     197       &  0.0000404196487152575*t**3*LOG(rh) + (15.79628615047088*LOG(rh))/x - &
     198       &  0.215553951893509*LOG(rh)**2 - 0.0810269192332194*t*LOG(rh)**2 +     &
     199       &  0.001435808434184642*t**2*LOG(rh)**2 -                               &
     200       &  4.775796947178588E-6*t**3*LOG(rh)**2 -                               &
     201       &  (2.912974063702185*LOG(rh)**2)/x - 3.588557942822751*LOG(rh)**3 +    &
     202       &  0.04950795302831703*t*LOG(rh)**3 -                                   &
     203       &  0.0002138195118737068*t**2*LOG(rh)**3 +                              &
     204       &  3.108005107949533E-7*t**3*LOG(rh)**3 -                               &
     205       &  (0.02933332747098296*LOG(rh)**3)/x +                                 &
     206       &  1.145983818561277*LOG(rhoa) -                                        &
     207       &  0.6007956227856778*t*LOG(rhoa) +                                     &
     208       &  0.00864244733283759*t**2*LOG(rhoa) -                                 &
     209       &  0.00002289467254710888*t**3*LOG(rhoa) -                              &
     210       &  (8.44984513869014*LOG(rhoa))/x +                                     &
     211       &  2.158548369286559*LOG(rh)*LOG(rhoa) +                                &
     212       &  0.0808121412840917*t*LOG(rh)*LOG(rhoa) -                             &
     213       &  0.0004073815255395214*t**2*LOG(rh)*LOG(rhoa) -                       &
     214       &  4.019572560156515E-7*t**3*LOG(rh)*LOG(rhoa) +                        &
     215       &  (0.7213255852557236*LOG(rh)*LOG(rhoa))/x +                           &
     216       &  1.62409850488771*LOG(rh)**2*LOG(rhoa) -                              &
     217       &  0.01601062035325362*t*LOG(rh)**2*LOG(rhoa) +                         &
     218       &  0.00003771238979714162*t**2*LOG(rh)**2*LOG(rhoa) +                   &
     219       &  3.217942606371182E-8*t**3*LOG(rh)**2*LOG(rhoa) -                     &
     220       &  (0.01132550810022116*LOG(rh)**2*LOG(rhoa))/x +                       &
     221       &  9.71681713056504*LOG(rhoa)**2 -                                      &
     222       &  0.1150478558347306*t*LOG(rhoa)**2 +                                  &
     223       &  0.0001570982486038294*t**2*LOG(rhoa)**2 +                            &
     224       &  4.009144680125015E-7*t**3*LOG(rhoa)**2 +                             &
     225       &  (0.7118597859976135*LOG(rhoa)**2)/x -                                &
     226       &  1.056105824379897*LOG(rh)*LOG(rhoa)**2 +                             &
     227       &  0.00903377584628419*t*LOG(rh)*LOG(rhoa)**2 -                         &
     228       &  0.00001984167387090606*t**2*LOG(rh)*LOG(rhoa)**2 +                   &
     229       &  2.460478196482179E-8*t**3*LOG(rh)*LOG(rhoa)**2 -                     &
     230       &  (0.05790872906645181*LOG(rh)*LOG(rhoa)**2)/x -                       &
     231       &  0.1487119673397459*LOG(rhoa)**3 +                                    &
     232       &  0.002835082097822667*t*LOG(rhoa)**3 -                                &
     233       &  9.24618825471694E-6*t**2*LOG(rhoa)**3 +                              &
     234       &  5.004267665960894E-9*t**3*LOG(rhoa)**3 -                             &
     235       &  (0.01270805101481648*LOG(rhoa)**3)/x
     236  jnuc=EXP(jnuc) !1/(cm3s)
     237
     238  ntot =-0.002954125078716302 - 0.0976834264241286*t + 0.001024847927067835*t**2 - 2.186459697726116E-6*t**3 -    &
     239       &   0.1017165718716887/x - 0.002050640345231486*LOG(rh) - 0.007585041382707174*t*LOG(rh) +                 &
     240       &   0.0001926539658089536*t**2*LOG(rh) - 6.70429719683894E-7*t**3*LOG(rh) -                                &
     241       &   (0.2557744774673163*LOG(rh))/x + 0.003223076552477191*LOG(rh)**2 + 0.000852636632240633*t*LOG(rh)**2 - &
     242       &   0.00001547571354871789*t**2*LOG(rh)**2 + 5.666608424980593E-8*t**3*LOG(rh)**2 +                        &
     243       &   (0.03384437400744206*LOG(rh)**2)/x + 0.04743226764572505*LOG(rh)**3 -                                  &
     244       &   0.0006251042204583412*t*LOG(rh)**3 + 2.650663328519478E-6*t**2*LOG(rh)**3 -                            &
     245       &   3.674710848763778E-9*t**3*LOG(rh)**3 - (0.0002672510825259393*LOG(rh)**3)/x -                          &
     246       &   0.01252108546759328*LOG(rhoa) + 0.005806550506277202*t*LOG(rhoa) -                                     &
     247       &   0.0001016735312443444*t**2*LOG(rhoa) + 2.881946187214505E-7*t**3*LOG(rhoa) +                           &
     248       &   (0.0942243379396279*LOG(rhoa))/x - 0.0385459592773097*LOG(rh)*LOG(rhoa) -                              &
     249       &   0.0006723156277391984*t*LOG(rh)*LOG(rhoa) + 2.602884877659698E-6*t**2*LOG(rh)*LOG(rhoa) +              &
     250       &   1.194163699688297E-8*t**3*LOG(rh)*LOG(rhoa) - (0.00851515345806281*LOG(rh)*LOG(rhoa))/x -              &
     251       &   0.01837488495738111*LOG(rh)**2*LOG(rhoa) + 0.0001720723574407498*t*LOG(rh)**2*LOG(rhoa) -              &
     252       &   3.717657974086814E-7*t**2*LOG(rh)**2*LOG(rhoa) -                                                       &
     253       &   5.148746022615196E-10*t**3*LOG(rh)**2*LOG(rhoa) +                                                      &
     254       &   (0.0002686602132926594*LOG(rh)**2*LOG(rhoa))/x - 0.06199739728812199*LOG(rhoa)**2 +                    &
     255       &   0.000906958053583576*t*LOG(rhoa)**2 - 9.11727926129757E-7*t**2*LOG(rhoa)**2 -                          &
     256       &   5.367963396508457E-9*t**3*LOG(rhoa)**2 - (0.007742343393937707*LOG(rhoa)**2)/x +                       &
     257       &   0.0121827103101659*LOG(rh)*LOG(rhoa)**2 - 0.0001066499571188091*t*LOG(rh)*LOG(rhoa)**2 +               &
     258       &   2.534598655067518E-7*t**2*LOG(rh)*LOG(rhoa)**2 -                                                       &
     259       &   3.635186504599571E-10*t**3*LOG(rh)*LOG(rhoa)**2 +                                                      &
     260       &   (0.0006100650851863252*LOG(rh)*LOG(rhoa)**2)/x + 0.0003201836700403512*LOG(rhoa)**3 -                  &
     261       &   0.0000174761713262546*t*LOG(rhoa)**3 + 6.065037668052182E-8*t**2*LOG(rhoa)**3 -                        &
     262       &   1.421771723004557E-11*t**3*LOG(rhoa)**3 + (0.0001357509859501723*LOG(rhoa)**3)/x
     263  ntot=EXP(ntot)
     264
     265  rc=EXP(-1.6524245+0.42316402*x+0.33466487*LOG(ntot)) !nm
     266
     267  IF (jnuc < 1.E-7) THEN
     268!     print *,'Warning (ilon=',ilon,'ilat=',ilat'): nucleation rate < 1E-7/cm3s, using 0.0/cm3s,'
    264269     jnuc=0.0
    265270  ENDIF
    266271
    267272  IF (jnuc > 1.e10) THEN
    268 !     print *,'Warning (ilon=',ilon,'ilev=',ilev, &
     273!     print *,'Warning (ilon=',ilon,'ilat=',ilat, &
    269274!        & '): nucleation rate > 1e10/cm3s, using 1e10/cm3s'
    270275     jnuc=1.e10
    271276  ENDIF
    272277
    273   rhotres=exp( -279.2430007512709 + 11.73439886096903*rh + 22700.92970508331/t &
     278  rhotres=EXP( -279.2430007512709 + 11.73439886096903*rh + 22700.92970508331/t &
    274279       & - (1088.644983466801*rh)/t + 1.144362942094912*t                      &
    275280       & - 0.03023314602163684*rh*t - 0.001302541390154324*t**2                &
    276        & - 6.386965238433532*Log(rh) + (854.980361026715*Log(rh))/t            &
    277        & + 0.00879662256826497*t*Log(rh)) !1/cm3
     281       & - 6.386965238433532*LOG(rh) + (854.980361026715*LOG(rh))/t            &
     282       & + 0.00879662256826497*t*LOG(rh)) !1/cm3
    278283
    279284  RETURN
    … …  
    281286END SUBROUTINE binapara
    282287
     288!---------------------------------------------------------------------------------------------------
     289
     290SUBROUTINE newbinapara(t,satrat,rhoa,csi,airn,ipr,jnuc_n,ntot_n,jnuc_i,ntot_i,        &
     291                   &   x_n,x_i,na_n,na_i,rc_n,rc_i,n_i,kinetic_n,kinetic_i,rhoatres)
     292
     293  !    Fortran 90 subroutine newbinapara
     294  !
     295  !    Calculates parametrized values for neutral and ion-induced sulfuric acid-water particle formation rate
     296  !    of critical clusters,
     297  !    number of particle in the critical clusters, the radii of the critical clusters
     298  !    in H2O-H2SO4-ion system if temperature, saturation ratio of water, sulfuric acid concentration,
     299  !    and, optionally, either condensation sink due to pre-existing particle and ion pair production rate,
     300  !    or atmospheric concentration of negative ions are given.
     301  !
     302  !    The code calculates also the kinetic limit and the particle formation rate
     303  !    above this limit (in which case we set ntot=1 and na=1)
     304  !
     305  !    Copyright (C)2018 Määttänen et al. 2018
     306  !   
     307  !    anni.maattanen@latmos.ipsl.fr
     308  !    joonas.merikanto@fmi.fi
     309  !    hanna.vehkamaki@helsinki.fi
     310  !
     311  !    References
     312  !    A. Määttänen, J. Merikanto, H. Henschel, J. Duplissy, R. Makkonen,
     313  !    I. K. Ortega and H. Vehkamäki (2018), New parameterizations for
     314  !    neutral and ion-induced sulfuric acid-water particle formation in
     315  !    nucleation and kinetic regimes, J. Geophys. Res. Atmos., 122, doi:10.1002/2017JD027429.
     316  !
     317  !    Brasseur, G., and A.  Chatel (1983),  paper  presented  at  the  9th  Annual  Meeting  of  the 
     318  !    European Geophysical Society, Leeds, Great Britain, August 1982.
     319  ! 
     320  !    Dunne, Eimear M., et al.(2016), Global atmospheric particle formation from CERN CLOUD measurements,
     321  !    Science 354.6316, 1119-1124.   
     322  !
     323
     324  USE aerophys
     325  USE YOMCST, ONLY : RPI, RKBOL
     326
     327  IMPLICIT NONE
     328
     329  !----------------------------------------------------
     330 
     331  !Global intent in
     332  DOUBLE PRECISION,INTENT(IN) :: t         ! temperature in K
     333  DOUBLE PRECISION,INTENT(IN) :: satrat    ! saturatio ratio of water (between zero and 1)
     334  DOUBLE PRECISION,INTENT(IN) :: rhoa      ! sulfuric acid concentration in 1/cm3
     335  DOUBLE PRECISION,INTENT(IN) :: csi       ! Ion condensation sink (s-1)
     336  DOUBLE PRECISION,INTENT(IN) :: airn      ! Air molecule concentration in (cm-3)
     337  DOUBLE PRECISION,INTENT(IN) :: ipr       ! Ion pair production rate (cm-3 s-1)
     338  !Global intent out
     339  DOUBLE PRECISION,INTENT(OUT) :: jnuc_n   ! Neutral nucleation rate in 1/cm3s (J>10^-7 1/cm3s)
     340  DOUBLE PRECISION,INTENT(OUT) :: ntot_n   ! total number of molecules in the neutral critical cluster
     341  DOUBLE PRECISION,INTENT(OUT) :: jnuc_i   ! Charged nucleation rate in 1/cm3s (J>10^-7 1/cm3s)
     342  DOUBLE PRECISION,INTENT(OUT) :: ntot_i   ! total number of molecules in the charged critical cluster
     343  DOUBLE PRECISION,INTENT(OUT) :: x_n      ! mole fraction of H2SO4 in the neutral critical cluster
     344  DOUBLE PRECISION,INTENT(OUT) :: x_i      ! mole fraction of H2SO4 in the charged critical cluster
     345                                           ! (note that x_n=x_i in nucleation regime)
     346  DOUBLE PRECISION,INTENT(OUT) :: na_n     ! sulfuric acid molecules in the neutral critical cluster
     347  DOUBLE PRECISION,INTENT(OUT) :: na_i     ! sulfuric molecules in the charged critical cluster
     348  DOUBLE PRECISION,INTENT(OUT) :: rc_n     ! radius of the charged critical cluster in nm
     349  DOUBLE PRECISION,INTENT(OUT) :: rc_i     ! radius of the charged critical cluster in nm
     350  DOUBLE PRECISION,INTENT(OUT) :: n_i      ! number of ion pairs in air (cm-3)
     351  LOGICAL,INTENT(OUT)  :: kinetic_n        ! true if kinetic neutral nucleation
     352  LOGICAL,INTENT(OUT)  :: kinetic_i        ! true if kinetic ion-induced nucleation
     353  DOUBLE PRECISION,INTENT(OUT) :: rhoatres ! threshold concentration of H2SO4 (1/cm^3) for neutral kinetic nucleation
     354
     355  ! Local
     356  DOUBLE PRECISION :: x           ! mole fraction of H2SO4 in the critical cluster
     357  DOUBLE PRECISION :: satratln    ! bounded water saturation ratio for neutral case (between 5.E-6 - 1.0)
     358  DOUBLE PRECISION :: satratli    ! bounded water saturation ratio for ion-induced case (between 1.E-7 - 0.95)
     359  DOUBLE PRECISION :: rhoaln      ! bounded concentration of h2so4 for neutral case (between 10^10 - 10^19 m-3)
     360  DOUBLE PRECISION :: rhoali      ! bounded concentration of h2so4 for ion-induced case (between 10^10 - 10^22 m-3)
     361  DOUBLE PRECISION :: tln         ! bounded temperature for neutral case (between 165-400 K)
     362  DOUBLE PRECISION :: tli         ! bounded temperature for ion-induced case (195-400 K)
     363  DOUBLE PRECISION :: kinrhotresn ! threshold sulfuric acid for neutral kinetic nucleation   
     364  DOUBLE PRECISION :: kinrhotresi ! threshold sulfuric acid for ion-induced kinetic nucleation
     365  DOUBLE PRECISION :: jnuc_i1     ! Ion-induced rate for n_i=1 cm-3
     366  DOUBLE PRECISION :: xloss       ! Ion loss rate
     367  DOUBLE PRECISION :: recomb      ! Ion-ion recombination rate
     368
     369  !--- 0) Initializations:
     370
     371  kinetic_n=.FALSE.
     372  kinetic_i=.FALSE.
     373  jnuc_n=0.0
     374  jnuc_i=0.0
     375  ntot_n=0.0
     376  ntot_i=0.0
     377  na_n=0.0
     378  na_i=0.0
     379  rc_n=0.0
     380  rc_i=0.0
     381  x=0.0
     382  x_n=0.0
     383  x_i=0.0
     384  satratln=satrat
     385  satratli=satrat
     386  rhoaln=rhoa
     387  rhoali=rhoa
     388  tln=t
     389  tli=t
     390  n_i=0.0
     391
     392  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     393
     394  !Boundary values according to parameterization limits   
     395
     396  !Temperature bounds
     397  IF (t.LE.165.) THEN
     398!     print *,'Warning: temperature < 165.0 K, using 165.0 K in neutral nucleation calculation'
     399     tln=165.0
     400  ENDIF
     401  IF (t.LE.195.) THEN
     402!     print *,'Warning: temperature < 195.0 K, using 195.0 K in ion-induced nucleation calculation'
     403     tli=195.0
     404  ENDIF
     405  IF (t.GE.400.) THEN
     406!     print *,'Warning: temperature > 400. K, using 400. K in nucleation calculations'
     407     tln=400.
     408     tli=400.
     409  ENDIF
     410 
     411  ! Saturation ratio bounds
     412  IF (satrat.LT.1.E-7) THEN
     413!     print *,'Warning: saturation ratio of water < 1.E-7, using 1.E-7 in ion-induced nucleation calculation'
     414     satratli=1.E-7
     415  ENDIF
     416  IF (satrat.LT.1.E-5) THEN
     417!     print *,'Warning: saturation ratio of water < 1.E-5, using 1.E-5 in neutral nucleation calculation'
     418     satratln=1.E-5
     419  ENDIF
     420  IF (satrat.GT.0.95) THEN
     421!     print *,'Warning: saturation ratio of water > 0.95, using 0.95 in ion-induced nucleation calculation'
     422     satratli=0.95
     423  ENDIF
     424  IF (satrat.GT.1.0) THEN
     425!     print *,'Warning: saturation ratio of water > 1 using 1 in neutral nucleation calculation'
     426     satratln=1.0
     427  ENDIF
     428 
     429  ! Sulfuric acid concentration bounds
     430  IF (rhoa.LE.1.E4) THEN
     431!     print *,'Warning: sulfuric acid < 1e4 1/cm3, using 1e4 1/cm3 in nucleation calculation'
     432     rhoaln=1.E4
     433     rhoali=1.E4
     434  ENDIF
     435  IF (rhoa.GT.1.E13) THEN
     436!     print *,'Warning: sulfuric acid > 1e13 1/cm3, using 1e13 1/cm3 in neutral nucleation calculation'
     437     rhoaln=1.E13
     438  ENDIF
     439  IF (rhoa.GT.1.E16) THEN
     440!     print *,'Warning: sulfuric acid concentration > 1e16 1/cm3, using 1e16 1/cm3 in ion-induced nucleation calculation'
     441     rhoali=1.E16
     442  ENDIF
     443 
     444  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     445 
     446  !Critical cluster composition (valid for both cases, bounds not used here)
     447  x_n=  7.9036365428891719E-1 - 2.8414059650092153E-3*tln + 1.4976802556584141E-2*LOG(satratln)  &
     448       & - 2.4511581740839115E-4*tln*LOG(satratln) + 3.4319869471066424E-3 *LOG(satratln)**2     &
     449       & - 2.8799393617748428E-5*tln*LOG(satratln)**2 + 3.0174314126331765E-4*LOG(satratln)**3   &
     450       & - 2.2673492408841294E-6*tln*LOG(satratln)**3 - 4.3948464567032377E-3*LOG(rhoaln)        &
     451       & + 5.3305314722492146E-5*tln*LOG(rhoaln)
     452  x_i=  7.9036365428891719E-1 - 2.8414059650092153E-3*tli + 1.4976802556584141E-2*LOG(satratli)  &
     453       & - 2.4511581740839115E-4*tli*LOG(satratli) + 3.4319869471066424E-3 *LOG(satratli)**2     &
     454       & - 2.8799393617748428E-5*tli*LOG(satratli)**2 + 3.0174314126331765E-4*LOG(satratli)**3   &
     455       & - 2.2673492408841294E-6*tli*LOG(satratli)**3 - 4.3948464567032377E-3*LOG(rhoali)        &
     456       & + 5.3305314722492146E-5*tli*LOG(rhoali)
     457       
     458  x_n=MIN(MAX(x_n,1.E-30),1.)
     459  x_i=MIN(MAX(x_i,1.E-30),1.)
     460 
     461  !Neutral nucleation
     462 
     463  !Kinetic limit check
     464  IF (satratln .GE. 1.E-2 .AND. satratln .LE. 1.) THEN
     465     kinrhotresn=EXP(7.8920778706888086E+1 + 7.3665492897447082*satratln - 1.2420166571163805E+4/tln &
     466          & + (-6.1831234251470971E+2*satratln)/tln - 2.4501159970109945E-2*tln                      &
     467          & -1.3463066443605762E-2*satratln*tln + 8.3736373989909194E-06*tln**2                      &
     468          & -1.4673887785408892*LOG(satratln) + (-3.2141890006517094E+1*LOG(satratln))/tln           &
     469          & + 2.7137429081917556E-3*tln*LOG(satratln)) !1/cm3     
     470     IF (kinrhotresn.LT.rhoaln) kinetic_n=.TRUE.
     471  ENDIF
     472
     473  IF (satratln .GE. 1.E-4  .AND. satratln .LT. 1.E-2) THEN     
     474     kinrhotresn=EXP(7.9074383049843647E+1 - 2.8746005462158347E+1*satratln - 1.2070272068458380E+4/tln &
     475          & + (-5.9205040320056632E+3*satratln)/tln - 2.4800372593452726E-2*tln                         &
     476          & -4.3983007681295948E-2*satratln*tln + 2.5943854791342071E-5*tln**2                          &
     477          & -2.3141363245211317*LOG(satratln) + (9.9186787997857735E+1*LOG(satratln))/tln               &
     478          & + 5.6819382556144681E-3*tln*LOG(satratln)) !1/cm3
     479     IF (kinrhotresn.LT.rhoaln) kinetic_n=.TRUE.
     480  ENDIF
     481
     482  IF (satratln .GE. 5.E-6  .AND. satratln .LT. 1.E-4) THEN
     483     kinrhotresn=EXP(8.5599712000361677E+1 + 2.7335119660796581E+3*satratln - 1.1842350246291651E+4/tln &
     484          & + (-1.2439843468881438E+6*satratln)/tln - 5.4536964974944230E-2*tln                         &
     485          & + 5.0886987425326087*satratln*tln + 7.1964722655507067E-5*tln**2                            &
     486          & -2.4472627526306372*LOG(satratln) + (1.7561478001423779E+2*LOG(satratln))/tln               &
     487          & + 6.2640132818141811E-3*tln*LOG(satratln)) !1/cm3
     488     IF (kinrhotresn.LT.rhoaln) kinetic_n=.TRUE.
     489  ENDIF
     490 
     491  IF (kinetic_n) THEN   
     492     ! Dimer formation rate
     493     jnuc_n=1.E6*(2.*0.3E-9)**2.*SQRT(8.*RPI*RKBOL*(1./mH2SO4mol+1./mH2SO4mol))/2.*SQRT(t)*rhoa**2.
     494!     jnuc_n=1.E6*(2.*0.3E-9)**2.*SQRT(8.*3.141593*1.38E-23*(1./(1.661E-27*98.07)+1./(1.661E-27*98.07)))/2.*SQRT(t)*rhoa**2.
     495     ntot_n=1. !set to 1
     496     na_n=1.   ! The critical cluster contains one molecules but the produced cluster contains 2 molecules
     497     x_n=na_n/ntot_n  ! so also set this to 1
     498     rc_n=0.3E-9
     499  ELSE
     500     jnuc_n= 2.1361182605986115E-1 + 3.3827029855551838*tln -3.2423555796175563E-2*tln**2 +         &
     501          &  7.0120069477221989E-5*tln**3 +8.0286874752695141/x_n +                                 &
     502          &  -2.6939840579762231E-1*LOG(satratln) +1.6079879299099518*tln*LOG(satratln) +           &
     503          &  -1.9667486968141933E-2*tln**2*LOG(satratln) +                                          &
     504          &  5.5244755979770844E-5*tln**3*LOG(satratln) + (7.8884704837892468*LOG(satratln))/x_n +  &
     505          &  4.6374659198909596*LOG(satratln)**2 - 8.2002809894792153E-2*tln*LOG(satratln)**2 +     &
     506          &  8.5077424451172196E-4*tln**2*LOG(satratln)**2 +                                        &
     507          &  -2.6518510168987462E-6*tln**3*LOG(satratln)**2 +                                       &
     508          &  (-1.4625482500575278*LOG(satratln)**2)/x_n - 5.2413002989192037E-1*LOG(satratln)**3 +  &
     509          &  5.2755117653715865E-3*tln*LOG(satratln)**3 +                                           &
     510          &  -2.9491061332113830E-6*tln**2*LOG(satratln)**3 +                                       &
     511          &  -2.4815454194486752E-8*tln**3*LOG(satratln)**3 +                                       &
     512          &  (-5.2663760117394626E-2*LOG(satratln)**3)/x_n +                                        &
     513          &  1.6496664658266762*LOG(rhoaln) +                                                       &
     514          &  -8.0809397859218401E-1*tln*LOG(rhoaln) +                                               &
     515          &  8.9302927091946642E-3*tln**2*LOG(rhoaln) +                                             &
     516          &  -1.9583649496497497E-5*tln**3*LOG(rhoaln) +                                            &
     517          &  (-8.9505572676891685*LOG(rhoaln))/x_n +                                                &
     518          &  -3.0025283601622881E+1*LOG(satratln)*LOG(rhoaln) +                                     &
     519          &  3.0783365644763633E-1*tln*LOG(satratln)*LOG(rhoaln) +                                  &
     520          &  -7.4521756337984706E-4*tln**2*LOG(satratln)*LOG(rhoaln) +                              &
     521          &  -5.7651433870681853E-7*tln**3*LOG(satratln)*LOG(rhoaln) +                              &
     522          &  (1.2872868529673207*LOG(satratln)*LOG(rhoaln))/x_n +                                   &
     523          &  -6.1739867501526535E-1*LOG(satratln)**2*LOG(rhoaln) +                                  &
     524          &  7.2347385705333975E-3*tln*LOG(satratln)**2*LOG(rhoaln) +                               &
     525          &  -3.0640494530822439E-5*tln**2*LOG(satratln)**2*LOG(rhoaln) +                           &
     526          &  6.5944609194346214E-8*tln**3*LOG(satratln)**2*LOG(rhoaln) +                            &
     527          &  (-2.8681650332461055E-2*LOG(satratln)**2*LOG(rhoaln))/x_n +                            &
     528          &  6.5213802375160306*LOG(rhoaln)**2 +                                                    &
     529          &  -4.7907162004793016E-2*tln*LOG(rhoaln)**2 +                                            &
     530          &  -1.0727890114215117E-4*tln**2*LOG(rhoaln)**2 +                                         &
     531          &  5.6401818280534507E-7*tln**3*LOG(rhoaln)**2 +                                          &
     532          &  (5.4113070888923009E-1*LOG(rhoaln)**2)/x_n +                                           &
     533          &  5.2062808476476330E-1*LOG(satratln)*LOG(rhoaln)**2 +                                   &
     534          &  -6.0696882500824584E-3*tln*LOG(satratln)*LOG(rhoaln)**2 +                              &
     535          &  2.3851383302608477E-5*tln**2*LOG(satratln)*LOG(rhoaln)**2 +                            &
     536          &  -1.5243837103067096E-8*tln**3*LOG(satratln)*LOG(rhoaln)**2 +                           &
     537          &  (-5.6543192378015687E-2*LOG(satratln)*LOG(rhoaln)**2)/x_n +                            &
     538          &  -1.1630806410696815E-1*LOG(rhoaln)**3 +                                                &
     539          &  1.3806404273119610E-3*tln*LOG(rhoaln)**3 +                                             &
     540          &  -2.0199865087650833E-6*tln**2*LOG(rhoaln)**3 +                                         &
     541          &  -3.0200284885763192E-9*tln**3*LOG(rhoaln)**3 +                                         &
     542          &  (-6.9425267104126316E-3*LOG(rhoaln)**3)/x_n
     543     jnuc_n=EXP(jnuc_n)
     544     
     545     ntot_n =-3.5863435141979573E-3 - 1.0098670235841110E-1*tln + 8.9741268319259721E-4*tln**2 - 1.4855098605195757E-6*tln**3 &
     546          &   - 1.2080330016937095E-1/x_n + 1.1902674923928015E-3*LOG(satratln) - 1.9211358507172177E-2*tln*LOG(satratln) +   &
     547          &   2.4648094311204255E-4*tln**2*LOG(satratln) - 7.5641448594711666E-7*tln**3*LOG(satratln) +                       &
     548          &   (-2.0668639384228818E-02*LOG(satratln))/x_n - 3.7593072011595188E-2*LOG(satratln)**2 +                          &
     549          &   9.0993182774415718E-4 *tln*LOG(satratln)**2 +                                                                   &
     550          &   -9.5698412164297149E-6*tln**2*LOG(satratln)**2 + 3.7163166416110421E-8*tln**3*LOG(satratln)**2 +                &
     551          &   (1.1026579525210847E-2*LOG(satratln)**2)/x_n + 1.1530844115561925E-2 *LOG(satratln)**3 +                        &
     552          &   - 1.8083253906466668E-4 *tln*LOG(satratln)**3 + 8.0213604053330654E-7*tln**2*LOG(satratln)**3 +                 &
     553          &   -8.5797885383051337E-10*tln**3*LOG(satratln)**3 + (1.0243693899717402E-3*LOG(satratln)**3)/x_n +                &
     554          &   -1.7248695296299649E-2*LOG(rhoaln) + 1.1294004162437157E-2*tln*LOG(rhoaln) +                                    &
     555          &   -1.2283640163189278E-4*tln**2*LOG(rhoaln) + 2.7391732258259009E-7*tln**3*LOG(rhoaln) +                          &
     556          &   (6.8505583974029602E-2*LOG(rhoaln))/x_n +2.9750968179523635E-1*LOG(satratln)*LOG(rhoaln) +                      &
     557          &   -3.6681154503992296E-3 *tln*LOG(satratln)*LOG(rhoaln) + 1.0636473034653114E-5*tln**2*LOG(satratln)*LOG(rhoaln)+ &
     558          &   5.8687098466515866E-9*tln**3*LOG(satratln)*LOG(rhoaln) + (-5.2028866094191509E-3*LOG(satratln)*LOG(rhoaln))/x_n+&
     559          &   7.6971988880587231E-4*LOG(satratln)**2*LOG(rhoaln) - 2.4605575820433763E-5*tln*LOG(satratln)**2*LOG(rhoaln) +   &
     560          &   2.3818484400893008E-7*tln**2*LOG(satratln)**2*LOG(rhoaln) +                                                     &
     561          &   -8.8474102392445200E-10*tln**3*LOG(satratln)**2*LOG(rhoaln) +                                                   &
     562          &   (-1.6640566678168968E-4*LOG(satratln)**2*LOG(rhoaln))/x_n - 7.7390093776705471E-2*LOG(rhoaln)**2 +              &
     563          &   5.8220163188828482E-4*tln*LOG(rhoaln)**2 + 1.2291679321523287E-6*tln**2*LOG(rhoaln)**2 +                        &
     564          &   -7.4690997508075749E-9*tln**3*LOG(rhoaln)**2 + (-5.6357941220497648E-3*LOG(rhoaln)**2)/x_n +                    &
     565          &   -4.7170109625089768E-3*LOG(satratln)*LOG(rhoaln)**2 + 6.9828868534370193E-5*tln*LOG(satratln)*LOG(rhoaln)**2 +  &
     566          &   -3.1738912157036403E-7*tln**2*LOG(satratln)*LOG(rhoaln)**2 +                                                    &
     567          &   2.3975538706787416E-10*tln**3*LOG(satratln)*LOG(rhoaln)**2 +                                                    &
     568          &   (4.2304213386288567E-4*LOG(satratln)*LOG(rhoaln)**2)/x_n + 1.3696520973423231E-3*LOG(rhoaln)**3 +               &
     569          &   -1.6863387574788199E-5*tln*LOG(rhoaln)**3 + 2.7959499278844516E-8*tln**2*LOG(rhoaln)**3 +                       &
     570          &   3.9423927013227455E-11*tln**3*LOG(rhoaln)**3 + (8.6136359966337272E-5*LOG(rhoaln)**3)/x_n
     571     ntot_n=EXP(ntot_n)
     572     
     573     rc_n=EXP(-22.378268374023630+0.44462953606125100*x_n+0.33499495707849131*LOG(ntot_n)) !in meters
     574     
     575     na_n=x_n*ntot_n
     576     IF (na_n .LT. 1.) THEN
     577        print *, 'Warning: number of acid molecules < 1 in nucleation regime, setting na_n=1'
     578        na_n=1.0
     579     ENDIF
     580  ENDIF
     581 
     582  ! Set the neutral nucleation rate to 0.0 if less than 1.0E-7     
     583  IF (jnuc_n.LT.1.E-7) THEN
     584     jnuc_n=0.0
     585  ENDIF
     586 
     587  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     588 
     589  ! Threshold neutral nucleation rate (j > 1/cm3s) parameterization (can be commented out if not needed)
     590  IF (tln .GE. 310.) THEN
     591     rhoatres=EXP(-2.8220714121794250 + 1.1492362322651116E+1*satratln -3.3034839106184218E+3/tln &
     592          & + (-7.1828571490168133E+2*satratln)/tln + 1.4649510835204091E-1*tln                   &
     593          & -3.0442736551916524E-2*satratln*tln -9.3258567137451497E-5*tln**2                     &
     594          & -1.1583992506895649E+1*LOG(satratln) + (1.5184848765906165E+3*LOG(satratln))/tln      &
     595          & + 1.8144983916747057E-2*tln*LOG(satratln)) !1/cm3
     596  ENDIF
     597
     598  IF (tln .GT. 190. .AND. tln .LT. 310.) THEN
     599     rhoatres=EXP(-3.1820396091231999E+2 + 7.2451289153199676*satratln + 2.6729355170089486E+4/tln &
     600          & + (-7.1492506076423069E+2*satratln)/tln + 1.2617291148391978*tln                       &
     601          & - 1.6438112080468487E-2*satratln*tln -1.4185518234553220E-3*tln**2                     &
     602          & -9.2864597847386694*LOG(satratln) + (1.2607421852455602E+3*LOG(satratln))/tln          &
     603          & + 1.3324434472218746E-2*tln*LOG(satratln)) !1/cm3
     604  ENDIF
     605
     606  IF (tln .LT. 185. .AND. tln .GT. 155.) THEN
     607     rhoatres=1.1788859232398459E+5 - 1.0244255702550814E+4*satratln +   &
     608          & 4.6815029684321962E+3*satratln**2 -1.6755952338499657E+2*tln
     609  ENDIF
     610 
     611  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     612 
     613  ! Ion-induced nucleation parameterization
     614 
     615  IF (ipr.GT.0.0) THEN ! if the ion production rate is above zero
     616     
     617     ! Calculate the ion induced nucleation rate wrt. concentration of 1 ion/cm3
     618     
     619     kinrhotresi = 5.3742280876674478E1  - 6.6837931590012266E-3 *LOG(satratli)**(-2)                                     &
     620          & - 1.0142598385422842E-01 * LOG(satratli)**(-1) - 6.4170597272606873E+00 * LOG(satratli)                       &
     621          & - 6.4315798914824518E-01 * LOG(satratli)**2 - 2.4428391714772721E-02 * LOG(satratli)**3                       &
     622          & - 3.5356658734539019E-04 * LOG(satratli)**4 + 2.5400015099140506E-05 * tli * LOG(satratli)**(-2)              &
     623          & - 2.7928900816637790E-04 * tli * LOG(satratli)**(-1) + 4.4108573484923690E-02 * tli * LOG(satratli)           &
     624          & + 6.3943789012475532E-03 * tli * LOG(satratli)**(2) + 2.3164296174966580E-04 * tli * LOG(satratli)**(3)       &
     625          & + 3.0372070669934950E-06 * tli * LOG(satratli)**4 + 3.8255873977423475E-06 * tli**2 * LOG(satratli)**(-1)     &
     626          & - 1.2344793083561629E-04 * tli**2 * LOG(satratli) - 1.7959048869810192E-05 * tli**2 * LOG(satratli)**(2)      &
     627          & - 3.2165622558722767E-07 * tli**2 * LOG(satratli)**3 - 4.7136923780988659E-09 * tli**3 * LOG(satratli)**(-1)  &
     628          & + 1.1873317184482216E-07 * tli**3 * LOG(satratli) + 1.5685860354866621E-08 * tli**3 * LOG(satratli)**2        &
     629          & - 1.4329645891059557E+04 * tli**(-1) + 1.3842599842575321E-01 * tli                                           &
     630          & - 4.1376265912842938E-04 * tli**(2) + 3.9147639775826004E-07 * tli**3
     631     
     632     kinrhotresi=EXP(kinrhotresi) !1/cm3
     633     
     634     IF (kinrhotresi.LT.rhoali) kinetic_i=.true.
     635     
     636     IF (kinetic_i) THEN   
     637!        jnuc_i1=1.0E6*(0.3E-9 + 0.487E-9)**2.*SQRT(8.*3.141593*1.38E-23*(1./(1.661E-27*98.07)+1./(1.661E-27*98.07)))*  &
     638        jnuc_i1=1.0E6*(0.3E-9 + 0.487E-9)**2.*SQRT(8.*RPI*RKBOL*(1./mH2SO4mol+1./mH2SO4mol))*  &
     639             &  SQRT(tli)*rhoali !1/cm3s 
     640        ntot_i=1. !set to 1
     641        na_i=1.
     642        x_i=na_i/ntot_i  ! so also set this to 1
     643        rc_i=0.487E-9
     644     ELSE
     645        jnuc_i1 = 3.0108954259038608E+01+tli*6.1176722090512577E+01+(tli**2)*8.7240333618891663E-01+(tli**3)*                 &
     646             & -4.6191788649375719E-03+(tli**(-1))*8.3537059107024481E-01 +                                                   &
     647             & (1.5028549216690628E+01+tli*-1.9310989753720623E-01+(tli**2)*8.0155514634860480E-04+(tli**3)*                  &
     648             & -1.0832730707799128E-06+(tli**(-1))*1.7577660457989019)*(LOG(satratli)**(-2)) +                                &
     649             & (-2.0487870170216488E-01 +  tli * 1.3263949252910405E-03 +  (tli**2) * -8.4195688402450274E-06 +               &
     650             & (tli**3)*1.6154895940993287E-08 + (tli**(-1))*3.8734212545203874E+01) * (LOG(satratli)**(-2)*LOG(rhoali)) +    &
     651             & (1.4955918863858371 +  tli * 9.2290004245522454E+01 +  (tli**2) * -8.9006965195392618E-01 +                    &
     652             & (tli**3) * 2.2319123411013099E-03 + (tli**(-1)) * 4.0180079996840852E-03) *                                    &
     653             & (LOG(satratli)**(-1) * LOG(rhoali)**(-1)) +                                                                    &
     654             & (7.9018031228561085 +  tli * -1.1649433968658949E+01 +  (tli**2) * 1.1400827854910951E-01 +                    &
     655             & (tli**3) * -3.1941526492127755E-04 + (tli**(-1)) * -3.7662115740271446E-01) * (LOG(satratli)**(-1)) +          &
     656             & (1.5725237111225979E+02 +  tli * -1.0051649979836277 +  (tli**2) * 1.1866484014507624E-03 +                    &
     657             & (tli**3) * 7.3557614998540389E-06 + (tli**(-1)) * 2.6270197023115189) * (LOG(satratli)**(-1) * LOG(rhoali)) +  &
     658             & (-1.6973840122470968E+01 +  tli * 1.1258423691432135E-01 +  (tli**2) * -2.9850139351463793E-04 + (tli**3) *    &
     659             & 1.4301286324827064E-07 + (tli**(-1)) * 1.3163389235253725E+01) * (LOG(satratli)**(-1) * LOG(rhoali)**2) +      &
     660             & (-1.0399591631839757 +  tli * 2.7022055588257691E-03 +  (tli**2) * -2.1507467231330936E-06 + (tli**3) *        &
     661             & 3.8059489037584171E-10 + (tli**(-1)) * 1.5000492788553410E+02) * (LOG(satratli)**(-1) * LOG(rhoali)**3) +      &
     662             & (1.2250990965305315 +  tli * 3.0495946490079444E+01 +  (tli**2) * 2.1051563135187106E+01 + (tli**3) *          &
     663             & -8.2200682916580878E-02 + (tli**(-1)) * 2.9965871386685029E-02) * (LOG(rhoali)**(-2)) +                        &
     664             & (4.8281605955680433 +  tli * 1.7346551710836445E+02 +  (tli**2) * -1.0113602140796010E+01 + (tli**3) *         &
     665             & 3.7482518458685089E-02 + (tli**(-1)) * -1.4449998158558205E-01) * (LOG(rhoali)**(-1)) +                        &
     666             & (2.3399230964451237E+02 +  tli * -2.3099267235261948E+01 +  (tli**2) * 8.0122962140916354E-02 +                &
     667             & (tli**3) * 6.1542576994557088E-05 + (tli**(-1)) * 5.3718413254843007) * (LOG(rhoali)) +                        &
     668             & (1.0299715519499360E+02 +  tli * -6.4663357203364136E-02 +  (tli**2) * -2.0487150565050316E-03 +               &
     669             & (tli**3) * 8.7935289055530897E-07 + (tli**(-1)) * 3.6013204601215229E+01) * (LOG(rhoali)**2) +                 &
     670             & (-3.5452115439584042 +  tli * 1.7083445731159330E-02 +  (tli**2) * -1.2552625290862626E-05 + (tli**3) *        &
     671             & 1.2968447449182847E-09 + (tli**(-1)) * 1.5748687512056560E+02) * (LOG(rhoali)**3) +                            &
     672             & (2.2338490119517975 +  tli * 1.0229410216045540E+02 +  (tli**2) * -3.2103611955174052 + (tli**3) *             &
     673             & 1.3397152304977591E-02 + (tli**(-1)) * -2.4155187776460030E-02) * (LOG(satratli)* LOG(rhoali)**(-2)) +         &
     674             & (3.7592282990713963 +  tli * -1.5257988769009816E+02 +  (tli**2) * 2.6113805420558802 + (tli**3) *             &
     675             & -9.0380721653694363E-03 + (tli**(-1)) * -1.3974197138171082E-01) * (LOG(satratli)* LOG(rhoali)**(-1)) +        &
     676             & (1.8293600730573988E+01 +  tli * 1.8344728606002992E+01 +  (tli**2) * -4.0063363221106751E-01 + (tli**3)       &
     677             & * 1.4842749371258522E-03 + (tli**(-1)) * 1.1848846003282287) * (LOG(satratli)) +                               &
     678             & (-1.7634531623032314E+02 +  tli * 4.9011762441271278 +  (tli**2) * -1.3195821562746339E-02 + (tli**3) *        &
     679             & -2.8668619526430859E-05 + (tli**(-1)) * -2.9823396976393551E-01) * (LOG(satratli)* LOG(rhoali)) +              &
     680             & (-3.2944043694275727E+01 +  tli * 1.2517571921051887E-01 +  (tli**2) * 8.3239769771186714E-05 + (tli**3) *     &
     681             & 2.8191859341519507E-07 + (tli**(-1)) * -2.7352880736682319E+01) * (LOG(satratli)* LOG(rhoali)**2) +            &
     682             & (-1.1451811137553243 +  tli * 2.0625997485732494E-03 +  (tli**2) * -3.4225389469233624E-06 + (tli**3) *        &
     683             & 4.4437613496984567E-10 + (tli**(-1)) * 1.8666644332606754E+02) * (LOG(satratli)* LOG(rhoali)**3) +             &
     684             & (3.2270897099493567E+01 +  tli * 7.7898447327513687E-01 +  (tli**2) * -6.5662738484679626E-03 + (tli**3) *     &
     685             & 3.7899330796456790E-06 + (tli**(-1)) * 7.1106427501756542E-01) * (LOG(satratli)**2 * LOG(rhoali)**(-1)) +      &
     686             & (-2.8901906781697811E+01 +  tli * -1.5356398793054860 +  (tli**2) * 1.9267271774384788E-02 + (tli**3) *        &
     687             & -5.3886270475516162E-05 + (tli**(-1)) * 5.0490415975693426E-01) * (LOG(satratli)**2) +                         &
     688             & (3.3365683645733924E+01 +  tli * -3.6114561564894537E-01 +  (tli**2) * 9.2977354471929262E-04 + (tli**3) *     &
     689             & 1.9549769069511355E-07 + (tli**(-1)) * -8.8865930095112855) * (LOG(satratli)**2 * LOG(rhoali)) +               &
     690             & (2.4592563042806375 +  tli * -8.3227071743101084E-03 +  (tli**2) * 8.2563338043447783E-06 + (tli**3) *         &
     691             & -8.4374976698593496E-09 + (tli**(-1)) * -2.0938173949893473E+02) * (LOG(satratli)**2 * LOG(rhoali)**2) +       &
     692             & (4.4099823444352317E+01 +  tli * 2.5915665826835252 +  (tli**2) * -1.6449091819482634E-02 + (tli**3) *         &
     693             & 2.6797249816144721E-05 + (tli**(-1)) * 5.5045672663909995E-01)* satratli
     694        jnuc_i1=EXP(jnuc_i1)
     695       
     696        ntot_i = ABS((-4.8324296064013375E+04 +  tli * 5.0469120697428906E+02 +  (tli**2) * -1.1528940488496042E+00 +         &
     697             & (tli**(-1)) * -8.6892744676239192E+02 + (tli**(3)) * 4.0030302028120469E-04) +                                 &
     698             & (-6.7259105232039847E+03 +  tli * 1.9197488157452008E+02 +  (tli**2) * -1.3602976930126354E+00 +               &
     699             & (tli**(-1)) * -1.1212637938360332E+02 + (tli**(3)) * 2.8515597265933207E-03) *                                 &
     700             & LOG(satratli)**(-2) * LOG(rhoali)**(-2) +                                                                      &
     701             & (2.6216455217763342E+02 +  tli * -2.3687553252750821E+00 +  (tli**2) * 7.4074554767517521E-03 +                &
     702             & (tli**(-1)) * -1.9213956820114927E+03 + (tli**(3)) * -9.3839114856129453E-06) * LOG(satratli)**(-2) +          &
     703             & (3.9652478944137344E+00 +  tli * 1.2469375098256536E-02 +  (tli**2) * -9.9837754694045633E-05 + (tli**(-1)) *  &
     704             & -5.1919499210175138E+02 + (tli**(3)) * 1.6489001324583862E-07) * LOG(satratli)**(-2) * LOG(rhoali) +           &
     705             & (2.4975714429096206E+02 +  tli * 1.7107594562445172E+02 +  (tli**2) * -7.8988711365135289E-01 + (tli**(-1)) *  &
     706             & -2.2243599782483177E+01 + (tli**(3)) * -1.6291523004095427E-04) * LOG(satratli)**(-1) * LOG(rhoali)**(-2) +    &
     707             & (-8.9270715592533611E+02 +  tli * 1.2053538883338946E+02 +  (tli**2) * -1.5490408828541018E+00 + (tli**(-1)) * &
     708             & -1.1243275579419826E+01 + (tli**(3)) * 4.8053105606904655E-03) * LOG(satratli)**(-1) * LOG(rhoali)**(-1) +     &
     709             & (7.6426441642091631E+03 +  tli * -7.1785462414656578E+01 +  (tli**2) * 2.3851864923199523E-01 + (tli**(-1)) *  &
     710             & 8.5591775688708395E+01 + (tli**(3)) * -3.7000473243342858E-04) * LOG(satratli)**(-1) +                         &
     711             & (-5.1516826398607911E+01 +  tli * 9.1385720811460558E-01 +  (tli**2) * -3.5477100262158974E-03 +               &
     712             & (tli**(-1)) * 2.7545544507625586E+03 + (tli**(3)) * 5.4708262093640928E-06) * LOG(satratli)**(-1) * LOG(rhoali) + &
     713             & (-3.0386767129196176E+02 +  tli * -1.1033438883583569E+04 +  (tli**2) * 8.1296859732896067E+01 + (tli**(-1)) * &
     714             & 1.2625883141097162E+01 + (tli**(3)) * -1.2728497822219101E-01) * LOG(rhoali)**(-2) +                           &
     715             & (-3.3763494256461472E+03 +  tli * 3.1916579136391006E+03 +  (tli**2) * -2.7234339474441143E+01 + (tli**(-1)) * &
     716             & -2.1897653262707397E+01 + (tli**(3)) * 5.1788505812259071E-02) * LOG(rhoali)**(-1) +                           &
     717             & (-1.8817843873687068E+03 +  tli * 4.3038072285882070E+00 +  (tli**2) * 6.6244087689671860E-03 + (tli**(-1)) *  &
     718             & -2.7133073605696295E+03 + (tli**(3)) * -1.7951557394285043E-05) * LOG(rhoali) +                                &
     719             & (-1.7668827539244447E+02 +  tli * 4.8160932330629913E-01 +  (tli**2) * -6.3133007671100293E-04 + (tli**(-1)) * &
     720             & 2.5631774669873157E+04 + (tli**(3)) * 4.1534484127873519E-07) * LOG(rhoali)**(2) +                             &
     721             & (-1.6661835889222382E+03 +  tli * 1.3708900504682877E+03 +  (tli**2) * -1.7919060052198969E+01 + (tli**(-1)) * &
     722             & -3.5145029804436405E+01 + (tli**(3)) * 5.1047240947371224E-02) * LOG(satratli)* LOG(rhoali)**(-2) +            &
     723             & (1.0843549363030939E+04 +  tli * -7.3557073636139577E+01 +  (tli**2) * 1.2054625131778862E+00 + (tli**(-1)) *  &
     724             & 1.9358737917864391E+02 + (tli**(3)) * -4.2871620775911338E-03) * LOG(satratli)* LOG(rhoali)**(-1) +            &
     725             & (-2.4269802549752835E+03 +  tli * 1.1348265061941714E+01 +  (tli**2) * -5.0430423939495157E-02 + (tli**(-1)) * &
     726             & 2.3709874548950634E+03 + (tli**(3)) * 1.4091851828620244E-04) * LOG(satratli) +                                &
     727             & (5.2745372575251588E+02 +  tli * -2.6080675912627314E+00 +  (tli**2) * 5.6902218056670145E-03 + (tli**(-1)) *  &
     728             & -3.2149319482897838E+04 + (tli**(3)) * -5.4121996056745853E-06) * LOG(satratli)* LOG(rhoali) +                 &
     729             & (-1.6401959518360403E+01 +  tli * 2.4322962162439640E-01 +  (tli**2) * 1.1744366627725344E-03 + (tli**(-1)) *  &
     730             & -8.2694427518413195E+03 + (tli**(3)) * -5.0028379203873102E-06)* LOG(satratli)**(2) +                          &
     731             & (-2.7556572017167782E+03 +  tli * 4.9293344495058264E+01 +  (tli**2) * -2.6503456520676050E-01 + (tli**(-1)) * &
     732             & 1.2130698030982167E+03 + (tli**(3)) * 4.3530610668042957E-04)* LOG(satratli)**2 * LOG(rhoali)**(-1) +          &
     733             & (-6.3419182228959192E+00 +  tli * 4.0636212834605827E-02 +  (tli**2) * -1.0450112687842742E-04 + (tli**(-1)) * &
     734             & 3.1035882189759656E+02 + (tli**(3)) * 9.4328418657873500E-08)* LOG(satratli)**(-3) +                           &
     735             & (3.0189213304689042E+03 +  tli * -2.3804654203861684E+01 +  (tli**2) * 6.8113013411972942E-02 + (tli**(-1)) *  &
     736             & 6.3112071081188913E+02 + (tli**(3)) * -9.4460854261685723E-05)* (satratli) * LOG(rhoali) +                     &
     737             & (1.1924791930673702E+04 +  tli * -1.1973824959206000E+02 +  (tli**2) * 1.6888713097971020E-01 + (tli**(-1)) *  &
     738             & 1.8735938211539585E+02 + (tli**(3)) * 5.0974564680442852E-04)* (satratli) +                                    &
     739             & (3.6409071302482083E+01 +  tli * 1.7919859306449623E-01 +  (tli**2) * -1.0020116255895206E-03 + (tli**(-1)) *  &
     740             & -8.3521083354432303E+03 + (tli**(3)) * 1.5879900546795635E-06)* satratli * LOG(rhoali)**(2))
     741         
     742        rc_i = (-3.6318550637865524E-08 +  tli * 2.1740704135789128E-09   +  (tli**2) *                          &
     743             & -8.5521429066506161E-12 + (tli**3) * -9.3538647454573390E-15) +                                   &
     744             & (2.1366936839394922E-08 +  tli * -2.4087168827395623E-10 +  (tli**2) * 8.7969869277074319E-13 +   &
     745             & (tli**3) * -1.0294466881303291E-15)* LOG(satratli)**(-2) * LOG(rhoali)**(-1) +                    &
     746             & (-7.7804007761164303E-10 +  tli * 1.0327058173517932E-11 +  (tli**2) * -4.2557697639692428E-14 +  &
     747             & (tli**3) * 5.4082507061618662E-17)* LOG(satratli)**(-2) +                                         &
     748             & (3.2628927397420860E-12 +  tli * -7.6475692919751066E-14 +  (tli**2) * 4.1985816845259788E-16 +   &
     749             & (tli**3) * -6.2281395889592719E-19)* LOG(satratli)**(-2) * LOG(rhoali) +                          &
     750             & (2.0442205540818555E-09 +  tli * 4.0441858911249830E-08 +  (tli**2) * -3.3423487629482825E-10 +   &
     751             & (tli**3) * 6.8000404742985678E-13)* LOG(satratli)**(-1) * LOG(rhoali)**(-2) +                     &
     752             & (1.8381489183824627E-08 +  tli * -8.9853322951518919E-09 +  (tli**2) * 7.5888799566036185E-11 +   &
     753             & (tli**3) * -1.5823457864755549E-13)* LOG(satratli)**(-1) * LOG(rhoali)**(-1) +                    &
     754             & (1.1795760639695057E-07 +  tli * -8.1046722896375875E-10 +  (tli**2) * 9.1868604369041857E-14 +   &
     755             & (tli**3) * 4.7882428237444610E-15)* LOG(satratli)**(-1) +                                         &
     756             & (-4.4028846582545952E-09 +  tli * 4.6541269232626618E-11 +  (tli**2) * -1.1939929984285194E-13 +  &
     757             & (tli**3) * 2.3602037016614437E-17)* LOG(satratli)**(-1) * LOG(rhoali) +                           &
     758             & (2.7885056884209128E-11 +  tli * -4.5167129624119121E-13 +  (tli**2) * 1.6558404997394422E-15 +   &
     759             & (tli**3) * -1.2037336621218054E-18)* LOG(satratli)**(-1) * LOG(rhoali)**2 +                       &
     760             & (-2.3719627171699983E-09 +  tli * -1.5260127909292053E-07 +  (tli**2) * 1.7177017944754134E-09 +  &
     761             & (tli**3) * -4.7031737537526395E-12)* LOG(rhoali)**(-2) +                                          &
     762             & (-5.6946433724699646E-09 +  tli * 8.4629788237081735E-09 +  (tli**2) * -1.7674135187061521E-10 +  &
     763             & (tli**3) * 6.6236547903091862E-13)* LOG(rhoali)**(-1) +                                           &
     764             & (-2.2808617930606012E-08 +  tli * 1.4773376696847775E-10 +  (tli**2) * -1.3076953119957355E-13 +  &
     765             & (tli**3) * 2.3625301497914000E-16)* LOG(rhoali) +                                                 &
     766             & (1.4014269939947841E-10 +  tli * -2.3675117757377632E-12 +  (tli**2) * 5.1514033966707879E-15 +   &
     767             & (tli**3) * -4.8864233454747856E-18)* LOG(rhoali)**2 +                                             &
     768             & (6.5464943868885886E-11 +  tli * 1.6494354816942769E-08 +  (tli**2) * -1.7480097393483653E-10 +   &
     769             & (tli**3) * 4.7460075628523984E-13)* LOG(satratli)* LOG(rhoali)**(-2) +                            &
     770             & (8.4737893183927871E-09 +  tli * -6.0243327445597118E-09 +  (tli**2) * 5.8766070529814883E-11 +   &
     771             & (tli**3) * -1.4926748560042018E-13)* LOG(satratli)* LOG(rhoali)**(-1) +                           &
     772             & (1.0761964135701397E-07 +  tli * -1.0142496009071148E-09 +  (tli**2) * 2.1337312466519190E-12 +   &
     773             & (tli**3) * 1.6376014957685404E-15)* LOG(satratli) +                                               &
     774             & (-3.5621571395968670E-09 +  tli * 4.1175339587760905E-11 +  (tli**2) * -1.3535372357998504E-13 +  &
     775             & (tli**3) * 8.9334219536920720E-17)* LOG(satratli)* LOG(rhoali) +                                  &
     776             & (2.0700482083136289E-11 +  tli * -3.9238944562717421E-13 +  (tli**2) * 1.5850961422040196E-15 +   &
     777             & (tli**3) * -1.5336775610911665E-18)* LOG(satratli)* LOG(rhoali)**2 +                              &
     778             & (1.8524255464416206E-09 +  tli * -2.1959816152743264E-11 +  (tli**2) * -6.4478119501677012E-14 +  &
     779             & (tli**3) * 5.5135243833766056E-16)* LOG(satratli)**2 * LOG(rhoali)**(-1) +                        &
     780             & (1.9349488650922679E-09 +  tli * -2.2647295919976428E-11 +  (tli**2) * 9.2917479748268751E-14 +   &
     781             & (tli**3) * -1.2741959892173170E-16)* LOG(satratli)**2 +                                           &
     782             & (2.1484978031650972E-11 +  tli * -9.3976642475838013E-14 +  (tli**2) * -4.8892738002751923E-16 +  &
     783             & (tli**3) * 1.4676120441783832E-18)* LOG(satratli)**2 * LOG(rhoali) +                              &
     784             & (6.7565715216420310E-13 +  tli * -3.5421162549480807E-15 +  (tli**2) * -3.4201196868693569E-18 +  &
     785             & (tli**3) * 2.2260187650412392E-20)* LOG(satratli)**3 * LOG(rhoali)
     786                   
     787        na_i=x_i*ntot_i
     788        IF (na_i .LT. 1.) THEN
     789!           print *, 'Warning: number of acid molecules < 1 in nucleation regime, setting na_n=1'
     790           na_n=1.0
     791        ENDIF
     792     ENDIF
     793   
     794     jnuc_i=jnuc_i1
     795     ! Ion loss rate (1/s)
     796     xloss=csi+jnuc_i
     797     
     798     ! Recombination (here following Brasseur and Chatel, 1983)   
     799     recomb=6.0E-8*SQRT(300./tli)+6.0E-26*airn*(300./tli)**4
     800     
     801     ! Small ion concentration in air (1/cm3) (following Dunne et al., 2016)
     802     ! max function is to avoid n_i to go practically zero at very high J_ion
     803     n_i=MAX(0.01,(SQRT(xloss**2.0+4.0*recomb*ipr)-xloss)/(2.0*recomb))
     804     
     805     ! Ion-induced nucleation rate
     806     ! Min function is to ensure that max function above does not cause J_ion to overshoot
     807     jnuc_i=MIN(ipr,n_i*jnuc_i1)
     808     ! Set the ion-induced nucleation rate to 0.0 if less than 1.0E-7     
     809     IF (jnuc_i.LT.1.E-7) THEN
     810        jnuc_i=0.0
     811     ENDIF
     812
     813  ENDIF
     814 
     815END SUBROUTINE newbinapara
     816
    283817END MODULE nucleation_tstep_mod
Note: See TracChangeset for help on using the changeset viewer.