MODULE real16 INTEGER,PARAMETER :: REAL_16 = SELECTED_REAL_KIND (p=14,r=100) END MODULE real16 program binarytest use real16 implicit none real(real_16) :: t,rh,rhoa,j,x,ntot,r,rhotres logical :: kinetic t = 240.0 !kelvin rh = 0.5 !saturation ratio, unitless rhoa = 1e14 !m-3 call newbinapara(t,rh,rhoa,j,x,ntot,r,kinetic,rhotres) print *,'t ',t print *,'rh ',rh print *,'rhoa ',rhoa print *,'j ',j print *,'x ',x print *,'ntot ',ntot print *,'r ',r print *,'kinetic? ',kinetic print *,'[acid] for J=1 ',rhotres end program binarytest SUBROUTINE newbinapara(t,rh,rhoa,jnuc,x,ntot,rc,kinetic,rhotres) ! Fortran 90 subroutine binapara ! ! Calculates parametrized values of particle formation rate, ! mole fraction of sulphuric acid, ! total number of particles, and the radius of the critical cluster ! in H2O-H2SO4 system if temperature, saturatio ratio of water and ! sulfuric acid concentration are given. ! ! Calculates also the kinetic limit and the particle formation rate ! above this limit (in which case we set ntot=1 and x=1) ! ! ! Copyright (C)2016 Määttänen et al. 2016 ! ! ! anni.maattanen@latmos.ipsl.fr ! joonas.merikanto@fmi.fi ! hanna.vehkamaki@helsinki.fi use real16 IMPLICIT NONE real(REAL_16),intent(inout) :: t !temperature in K (165-400 K) real(REAL_16),intent(inout) :: rh !saturatio ratio of water (0.00001-1) real(REAL_16),intent(inout) :: rhoa !sulfuric acid concentration in 1/m3 (10^10 - 10^19) real(REAL_16),intent(out) :: jnuc !nucleation rate in 1/cm3s (J>10^-7 1/cm3s) real(REAL_16),intent(out) :: ntot !total number of molecules in the critical cluster real(REAL_16),intent(out) :: x ! mole fraction of H2SO4 in the critical cluster real(REAL_16),intent(out) :: rc !radius of the critical cluster in nm real(REAL_16),intent(out) :: rhotres ! treshold concentration of h2so4 (1/cm^3) !which produces nucleation rate 1/(cm^3 s) as a function of rh and t real(REAL_16),intent(out) :: kinrhotres ! treshold concentration of h2so4 (1/cm^3) !of the kinetic limit as a function of rh and t real(REAL_16) :: nac !local variable for number of acid molecules logical,intent(out) :: kinetic kinetic=.false. jnuc=0.0 x=0.0 ntot=0.0 rc=0.0 rhotres=0.0 kinrhotres=0.0 if (rh .ge. 1.e-2 .and. rh .lt. 1.) then kinrhotres=exp(7.8920778706888086e+1 + 7.3665492897447082*rh - 1.2420166571163805e+4/t & & + (-6.1831234251470971e+2*rh)/t - 2.4501159970109945e-2*t & & -1.3463066443605762e-2*rh*t + 8.3736373989909194e-06*t**2 & & -1.4673887785408892*Log(rh) + (-3.2141890006517094e+1*Log(rh))/t & & + 2.7137429081917556e-3*t*Log(rh)) if(kinrhotres*1.E6.lt.rhoa) kinetic=.true. endif if (rh .gt. 1.e-4 .and. rh .lt. 1.e-2) then kinrhotres=exp(7.9074383049843647e+1 - 2.8746005462158347e+1*rh - 1.2070272068458380e+4/t & & + (-5.9205040320056632e+3*rh)/t - 2.4800372593452726e-2*t & & -4.3983007681295948e-2*rh*t + 2.5943854791342071e-5*t**2 & & -2.3141363245211317*Log(rh) + (9.9186787997857735e+1*Log(rh))/t & & + 5.6819382556144681e-3*t*Log(rh)) if(kinrhotres*1.E6.lt.rhoa) kinetic=.true. endif if (rh .gt. 5.e-6 .and. rh .lt. 1.e-4) then kinrhotres=exp(8.5599712000361677e+1 + 2.7335119660796581e+3*rh - 1.1842350246291651e+4/t & & + (-1.2439843468881438e+6*rh)/t - 5.4536964974944230e-2*t & & + 5.0886987425326087*rh*t + 7.1964722655507067e-5*t**2 & & -2.4472627526306372*Log(rh) + (1.7561478001423779e+2*Log(rh))/t & & + 6.2640132818141811e-3*t*Log(rh)) if(kinrhotres*1.E6.lt.rhoa) kinetic=.true. endif if(kinetic) then jnuc=(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./1.E6 ntot=1. x=1. rc=0.3E-9 else if(t < 165.) then print *,'Warning: temperature < 165.0 K, using 165.0 K' t=165.0 end if if(t > 400.) then print *,'Warning: temperature > 400. K, using 400. K' t=400. end if if(rh < 1.E-5) then print *,'Warning: saturation ratio of water < 1.e-5, using 1.e-5' rh=1.E-5 end if if(rh > 1.0) then print *,'Warning: saturation ratio of water > 1 using 1' rh=1.0 end if if(rhoa < 1.e10) then print *,'Warning: sulfuric acid concentration < 1e4 1/cm3, using 1e4 1/cm3' rhoa=1.e10 end if if(rhoa > 1.e19) then print *,'Warning: sulfuric acid concentration > 1e13 1/cm3, using 1e13 1/cm3' rhoa=1.e19 end if x= 7.9036365428891719e-1 - 2.8414059650092153e-3*t + 1.4976802556584141e-2*LOG(rh) & & - 2.4511581740839115e-4*t*LOG(rh) + 3.4319869471066424e-3 *LOG(rh)**2 & & - 2.8799393617748428e-5*t*LOG(rh)**2 + 3.0174314126331765e-4*LOG(rh)**3 & & - 2.2673492408841294e-6*t*LOG(rh)**3 - 4.3948464567032377e-3*LOG(rhoa/1.e6)& & + 5.3305314722492146e-5*t*LOG(rhoa/1.e6) jnuc= 2.1361182605986115e-1 + 3.3827029855551838 *t -3.2423555796175563e-2*t**2 + & & 7.0120069477221989e-5*t**3 +8.0286874752695141/x + & & -2.6939840579762231e-1*LOG(rh) +1.6079879299099518*t*LOG(rh) + & & -1.9667486968141933e-2*t**2*LOG(rh) + & & 5.5244755979770844e-5*t**3*LOG(rh) + (7.8884704837892468*LOG(rh))/x + & & 4.6374659198909596*LOG(rh)**2 - 8.2002809894792153e-2*t*LOG(rh)**2 + & & 8.5077424451172196e-4*t**2*LOG(rh)**2 + & & -2.6518510168987462e-6*t**3*LOG(rh)**2 + & & (-1.4625482500575278*LOG(rh)**2)/x - 5.2413002989192037e-1*LOG(rh)**3 + & & 5.2755117653715865e-3*t*LOG(rh)**3 + & & -2.9491061332113830e-6*t**2*LOG(rh)**3 + & & -2.4815454194486752e-8*t**3*LOG(rh)**3 + & & (-5.2663760117394626e-2*LOG(rh)**3)/x + & & 1.6496664658266762*LOG(rhoa/1.e6) + & & -8.0809397859218401e-1*t*LOG(rhoa/1.e6) + & & 8.9302927091946642e-3*t**2*LOG(rhoa/1.e6) + & & -1.9583649496497497e-5*t**3*LOG(rhoa/1.e6) + & & (-8.9505572676891685*LOG(rhoa/1.e6))/x + & & -3.0025283601622881e+1*LOG(rh)*LOG(rhoa/1.e6) + & & 3.0783365644763633e-1*t*LOG(rh)*LOG(rhoa/1.e6) + & & -7.4521756337984706e-4*t**2*LOG(rh)*LOG(rhoa/1.e6) + & & -5.7651433870681853e-7*t**3*LOG(rh)*LOG(rhoa/1.e6) + & & (1.2872868529673207*LOG(rh)*LOG(rhoa/1.e6))/x + & & -6.1739867501526535e-1*LOG(rh)**2*LOG(rhoa/1.e6) + & & 7.2347385705333975e-3*t*LOG(rh)**2*LOG(rhoa/1.e6) + & & -3.0640494530822439e-5*t**2*LOG(rh)**2*LOG(rhoa/1.e6) + & & 6.5944609194346214e-8*t**3*LOG(rh)**2*LOG(rhoa/1.e6) + & & (-2.8681650332461055e-2*LOG(rh)**2*LOG(rhoa/1.e6))/x + & & 6.5213802375160306*LOG(rhoa/1.e6)**2 + & & -4.7907162004793016e-2*t*LOG(rhoa/1.e6)**2 + & & -1.0727890114215117e-4*t**2*LOG(rhoa/1.e6)**2 + & & 5.6401818280534507e-7*t**3*LOG(rhoa/1.e6)**2 + & & (5.4113070888923009e-1*LOG(rhoa/1.e6)**2)/x + & & 5.2062808476476330e-1*LOG(rh)*LOG(rhoa/1.e6)**2 + & & -6.0696882500824584e-3*t*LOG(rh)*LOG(rhoa/1.e6)**2 + & & 2.3851383302608477e-5*t**2*LOG(rh)*LOG(rhoa/1.e6)**2 + & & -1.5243837103067096e-8*t**3*LOG(rh)*LOG(rhoa/1.e6)**2 + & & (-5.6543192378015687e-2*LOG(rh)*LOG(rhoa/1.e6)**2)/x + & & -1.1630806410696815e-1*LOG(rhoa/1.e6)**3 + & & 1.3806404273119610e-3*t*LOG(rhoa/1.e6)**3 + & & -2.0199865087650833e-6*t**2*LOG(rhoa/1.e6)**3 + & & -3.0200284885763192e-9*t**3*LOG(rhoa/1.e6)**3 + & & (-6.9425267104126316e-3*LOG(rhoa/1.e6)**3)/x ntot =-3.5863435141979573e-3 - 1.0098670235841110e-1 *t + 8.9741268319259721e-4 *t**2 - 1.4855098605195757e-6*t**3 & & - 1.2080330016937095e-1/x + 1.1902674923928015e-3*LOG(rh) - 1.9211358507172177e-2*t*LOG(rh) + & & 2.4648094311204255e-4*t**2*LOG(rh) - 7.5641448594711666e-7*t**3*LOG(rh) + & & (-2.0668639384228818e-02*LOG(rh))/x - 3.7593072011595188e-2*LOG(rh)**2 + 9.0993182774415718e-4 *t*LOG(rh)**2 + & & -9.5698412164297149e-6*t**2*LOG(rh)**2 + 3.7163166416110421e-8*t**3*LOG(rh)**2 + & & (1.1026579525210847e-2*LOG(rh)**2)/x + 1.1530844115561925e-2 *LOG(rh)**3 + & & - 1.8083253906466668e-4 *t*LOG(rh)**3 + 8.0213604053330654e-7*t**2*LOG(rh)**3 + & & -8.5797885383051337e-10*t**3*LOG(rh)**3 + (1.0243693899717402e-3*LOG(rh)**3)/x + & & -1.7248695296299649e-2*LOG(rhoa/1.e6) + 1.1294004162437157e-2*t*LOG(rhoa/1.e6) + & & -1.2283640163189278e-4*t**2*LOG(rhoa/1.e6) + 2.7391732258259009e-7*t**3*LOG(rhoa/1.e6) + & & (6.8505583974029602e-2*LOG(rhoa/1.e6))/x +2.9750968179523635e-1*LOG(rh)*LOG(rhoa/1.e6) + & & -3.6681154503992296e-3 *t*LOG(rh)*LOG(rhoa/1.e6) + 1.0636473034653114e-5*t**2*LOG(rh)*LOG(rhoa/1.e6) + & & 5.8687098466515866e-9*t**3*LOG(rh)*LOG(rhoa/1.e6) + (-5.2028866094191509e-3*LOG(rh)*LOG(rhoa/1.e6))/x + & & 7.6971988880587231e-4*LOG(rh)**2*LOG(rhoa/1.e6) - 2.4605575820433763e-5*t*LOG(rh)**2*LOG(rhoa/1.e6) + & & 2.3818484400893008e-7*t**2*LOG(rh)**2*LOG(rhoa/1.e6) + & & -8.8474102392445200e-10*t**3*LOG(rh)**2*LOG(rhoa/1.e6) + & & (-1.6640566678168968e-4*LOG(rh)**2*LOG(rhoa/1.e6))/x - 7.7390093776705471e-2*LOG(rhoa/1.e6)**2 + & & 5.8220163188828482e-4*t*LOG(rhoa/1.e6)**2 + 1.2291679321523287e-6*t**2*LOG(rhoa/1.e6)**2 + & & -7.4690997508075749e-9*t**3*LOG(rhoa/1.e6)**2 + (-5.6357941220497648e-3*LOG(rhoa/1.e6)**2)/x + & & -4.7170109625089768e-3*LOG(rh)*LOG(rhoa/1.e6)**2 + 6.9828868534370193e-5*t*LOG(rh)*LOG(rhoa/1.e6)**2 + & & -3.1738912157036403e-7*t**2*LOG(rh)*LOG(rhoa/1.e6)**2 + & & 2.3975538706787416e-10*t**3*LOG(rh)*LOG(rhoa/1.e6)**2 + & & (4.2304213386288567e-4*LOG(rh)*LOG(rhoa/1.e6)**2)/x + 1.3696520973423231e-3*LOG(rhoa/1.e6)**3 + & & -1.6863387574788199e-5*t*LOG(rhoa/1.e6)**3 + 2.7959499278844516e-8*t**2*LOG(rhoa/1.e6)**3 + & & 3.9423927013227455e-11*t**3*LOG(rhoa/1.e6)**3 + (8.6136359966337272e-5*LOG(rhoa/1.e6)**3)/x ntot=EXP(ntot) rc=EXP(-22.378268374023630+0.44462953606125100*x+0.33499495707849131*LOG(ntot)) jnuc=EXP(jnuc) if(jnuc < 1.e-7) then print *,'Warning: nucleation rate < 1e-7/cm3s, using 0.0/cm3s,' jnuc=0.0 endif nac =x*ntot if (nac .lt. 1.) then print *, 'Warning: number of acid molecules < 1 in nucleation regime, setting x=1' x=1. endif endif if (t .ge. 310.) then rhotres=exp(-3.1820396091231999e+2 + 7.2451289153199676*rh + 2.6729355170089486e+4/t & & + (-7.1492506076423069e+2*rh)/t + 1.2617291148391978*t & & -1.6438112080468487e-2*rh*t -1.4185518234553220e-3*t**2 & & -9.2864597847386694*Log(rh) + (1.2607421852455602e+3*Log(rh))/t & & + 1.3324434472218746e-2*t*Log(rh)) endif if (t .gt. 190. .and. t .lt. 310.) then rhotres=exp(-3.1820396091231999e+2 + 7.2451289153199676*rh + 2.6729355170089486e+4/t & & + (-7.1492506076423069e+2*rh)/t + 1.2617291148391978*t & & - 1.6438112080468487e-2*rh*t -1.4185518234553220e-3*t**2 & & -9.2864597847386694*Log(rh) + (1.2607421852455602e+3*Log(rh))/t & & + 1.3324434472218746e-2*t*Log(rh)) endif if (t .lt. 185. .and. t .gt. 155.) then rhotres=1.1788859232398459e+5 - 1.0244255702550814e+4*rh + 4.6815029684321962e+3*rh**2 -1.6755952338499657e+2*t endif RETURN END SUBROUTINE newbinapara