source: trunk/LMDZ.VENUS/libf/phyvenus/cloudvenus/nucleations.F90 @ 3556

Last change on this file since 3556 was 1661, checked in by slebonnois, 8 years ago

SL: Cloud model for Venus. Not validated yet.
File size: 13.1 KB
Line 
1
2!  SUBROUTINE    HONUMO        HOmogeneous NUcleation with MOments
3!  SUBROUTINE    HETNUMO       HETerogeneous NUcleation with MOments
4!  SUBROUTINE    newnuklefit   Calculates jnuc and rc
5!  SUBROUTINE    heteronuk     Calculates the sum of activated aerosols
6
7
8!****************************************************************************
9! HOmogeneous NUcleation with MOments
10!
11! This subroutine uses the newnuclefit function.
12! Hanna Vehkamäki et al., 2002 : An improved parameterization for
13! sulfuric acid/water nucleation rates for tropospheric
14! and stratospheric conditions, () J. Geophys. Res., 107, pp. 4622-4631
15!
16! we determine here the tendency of M0
17! the created paticles will be in mode 1 (the smallest particles).
18!
19! The inputs mk0 and mk3 are not use in this function.
20!
21! An implicit method are used (easy)
22!
23! by S.Guilbon (LATMOS), 2016.
24!
25SUBROUTINE HONUMO(TAIR,dt,mk0,mk3,dM0_homo,dM3_homo)
26
27  use free_param
28  use donnees
29
30  IMPLICIT NONE
31
32  ! INputs
33  real, intent(in), dimension(2) :: mk0   ! first moment
34  real, intent(in), dimension(3) :: mk3   ! third moment
35  real, intent(in)  :: TAIR, dt           ! temperature, timestep
36  ! Outputs
37  real, intent(out) :: dM0_homo, dM3_homo ! Variations of moments
38
39  !Local variables
40  real :: jnuc, rc        ! Nucleation rate and critical radius
41  real :: m_H2SO4, m_H2O  ! Total condensed masses
42  real :: CWV, ntotal
43
44  ! Calculate jnuc and rc
45  CALL newnuklefit(TAIR,h2so4_m3,jnuc,rc,ntotal)     
46
47  ! initialization
48  dM0_homo = 0.D0
49  dM3_homo = 0.D0
50
51  ! Conditions to have nucleation process
52  IF (jnuc .ge. 1.D-1 .and. jnuc .le. 1.D16 .and. ntotal .gt. 4.0) then
53     CWV = MAIR/MWV
54     rc = 1.D-9  !m  Because of the calculation of WSA
55
56     ! Modal tendancies: Total nb density and vol. of created particles
57     ! Same thing with implicit scheme (cf thesis)
58     dM0_homo = dt * jnuc          !#.m(air)-3
59     dM3_homo = dt * jnuc * rc**3  !m3.m(air)-3
60  ENDIF
61
62  RETURN
63
64END SUBROUTINE HONUMO
65
66
67!****************************************************************************
68! HETerogeneous NUcleation with MOment
69!
70! Here, we use a simple representation of nucleation rate
71! and heterogeneous nucleation, like in VenLA (see Maattanen Anni).
72!
73! The number of activated ccn will become the number of created droplets
74! ra is the mean radius of the aerosol distribution.
75!
76! for mk0 calculation, method 1 and method 2 give same results :)
77!  mk0 = Jhetsum ! Number of activated CCN, in # - (method 1)
78! Here, we use implicit scheme (cf thesis)
79!
80SUBROUTINE HETNUMO(TAIR,dt,ra,mk0aer,mk3aer,mk0drop, mk3drop, &
81     & dM0_aer,dM3_aer,dM0_het,dM3_het)
82
83  use free_param
84  use donnees
85 
86  IMPLICIT NONE
87 
88  ! Inputs and Outputs
89  real, intent(in)  :: dt,TAIR,ra           ! time step, temperature, mean radius
90  real, intent(in)  :: mk0aer, mk3aer       ! original moments M(t)
91  real, intent(in), dimension(2) :: mk0drop ! original moment M(t)
92  real, intent(in), dimension(3) :: mk3drop ! original moment M(t)
93  real, intent(out) :: dM0_het, dM3_het     ! tendencies of liquid droplets
94  real, intent(out) :: dM0_aer, dM3_aer     ! tendencies of aerosols
95  ! Local variables
96  real :: Jhetsum, alpha_loc1, alpha_loc2
97  real :: mk0, mk3, alpha_k
98
99  call build_radius_grid(nbin,rmin,rmax,vratio)
100  call heteronuk(TAIR,ra,mk0aer,Jhetsum)
101
102  write(*,*)'mk0aer',mk0aer
103
104  alpha_loc1 = alpha_k(2,sig_aer)
105  mk0 = 1.D0 + (4.D0*PI*Jhetsum*r_aer*r_aer*alpha_loc1)*dt
106  mk0 = (1.D0/mk0)   
107  mk0 = mk0 * mk0aer ! - (method 2)
108  mk0 = mk0 - mk0aer ! Delta betwen 2 time steps = tendancy
109
110  alpha_loc2 = alpha_k(5,sig_aer)/alpha_k(3,sig_aer)
111  mk3 = 1.D0 + (4.D0*PI*Jhetsum*r_aer*r_aer*alpha_loc2)*dt
112  mk3 = 1.D0/mk3
113  mk3 = mk3 * mk3aer
114  mk3 = mk3 - mk3aer ! Delta betwen 2 time steps = tendancy
115  write(*,*)'HET - mk3 - methode 2', mk3
116
117  ! Modal tendancies
118  ! the lost quantity in aerosol will go to droplets.
119  dM0_aer = mk0
120  dM3_aer = mk3
121  dM0_het = mk0  * (-1)     
122  dM3_het = mk3  * (-1)
123
124END SUBROUTINE HETNUMO
125
126
127!*****************************************************************************
128SUBROUTINE newnuklefit(t,rhoa,jnuc,rc,ntotal)
129
130  !    Fortran 90 subroutine binapara
131  !
132  !    Calculates parametrized values of particle formation rate,
133  !    mole fraction of sulphuric acid,
134  !    total number of particles, and the radius of the critical cluster
135  !    in H2O-H2SO4 system if temperature, saturatio ratio of water and
136  !    sulfuric acid concentration  are given.
137  !
138  !    Calculates also the kinetic limit and the particle formation rate
139  !    above this limit (in which case we set ntotal=1 and x=1)
140  !
141  !    NB: produces nucleation rate   1/(cm^3 s) as a function of rh and t
142  !    of the kinetic limit as a function of rh and t
143  !
144  !    Copyright (C)2016 Määttänen et al. 2016
145  !
146  !    anni.maattanen@latmos.ipsl.fr
147  !    joonas.merikanto@fmi.fi
148  !    hanna.vehkamaki@helsinki.fi
149  !
150  !    modified by S.Guilbon for the Venus GCM (2016)
151  !
152  ! SG: change the unity of jnuc here
153  ! SG: change inputs and outputs
154
155  use free_param
156  use donnees
157
158  ! Inputs
159  real,intent(inout) :: t        ! temperature in K (165-400 K)
160  real,intent(inout) :: rhoa     ! sulfuric acid concentration in 1/m3 (10^10 - 10^19)
161
162  ! Outputs
163  real,intent(out) :: jnuc       ! nucleation rate in 1/m3s (J>10^-1 1/m3s)
164  real,intent(out) :: rc         ! radius of the critical cluster in m
165  real,intent(out) :: ntotal       ! total number of molecules in the critical cluster
166
167  ! Local variables
168  real :: x          ! mole fraction of H2SO4 in the critical cluster     
169  real :: nac        ! local variable for number of acid molecules
170
171  ! Initialization
172  jnuc = 0.D0
173  ntotal = 0.D0
174  rc = 0.D0
175  x = 0.D0   
176
177  if(t < 165.D0) then
178     print *,'Warning: temperature < 165.0 K, using 165.0 K'
179     t=165.D0
180  end if
181
182  if(t > 400.D0) then
183     print *,'Warning: temperature > 400. K, using 400. K'
184     t=400.D0
185  end if
186
187  if(rh < 1.D-5) then
188     print *,'Warning: saturation ratio of water < 1.e-5, using 1.e-5'
189     rh=1.D-5
190  end if
191
192  if(rh > 1.D0) then
193     print *,'Warning: saturation ratio of water > 1 using 1'
194     rh=1.D0
195  end if
196
197  if(rhoa < 1.D10) then
198     print *,'Warning: sulfuric acid concentration < 1e4 1/cm3, using 1e4 1/cm3'
199     rhoa=1.D10
200  end if
201
202  if(rhoa > 1.D19) then
203     print *,'Warning: sulfuric acid concentration > 1e13 1/cm3, using 1e13 1/cm3'
204     rhoa=1.D19
205  end if
206
207  x=  7.9036365428891719D-1 - 2.8414059650092153D-3*t + 1.4976802556584141D-2*LOG(rh) &
208       & - 2.4511581740839115D-4*t*LOG(rh) + 3.4319869471066424D-3 *LOG(rh)**2     & 
209       & - 2.8799393617748428D-5*t*LOG(rh)**2 + 3.0174314126331765D-4*LOG(rh)**3 &
210       & - 2.2673492408841294D-6*t*LOG(rh)**3 - 4.3948464567032377D-3*LOG(rhoa/1.D6)&
211       & + 5.3305314722492146D-5*t*LOG(rhoa/1.D6)
212
213
214  jnuc= 2.1361182605986115D-1 + 3.3827029855551838D0 *t -3.2423555796175563D-2*t**2 +  &
215       &  7.0120069477221989D-5*t**3 + 8.0286874752695141D0/x +  &
216       &  -2.6939840579762231D-1*LOG(rh) + 1.6079879299099518D0*t*LOG(rh) +  &
217       &  -1.9667486968141933D-2*t**2*LOG(rh) +  &
218       &  5.5244755979770844D-5*t**3*LOG(rh) + (7.8884704837892468D0*LOG(rh))/x +  &
219       &  4.6374659198909596D0*LOG(rh)**2 - 8.2002809894792153D-2*t*LOG(rh)**2 +  &
220       &  8.5077424451172196D-4*t**2*LOG(rh)**2 +  &
221       &  -2.6518510168987462D-6*t**3*LOG(rh)**2 +  &
222       &  (-1.4625482500575278D0*LOG(rh)**2)/x - 5.2413002989192037D-1*LOG(rh)**3 +  &
223       &  5.2755117653715865D-3*t*LOG(rh)**3 +  &
224       &  -2.9491061332113830D-6*t**2*LOG(rh)**3 +  &
225       &  -2.4815454194486752D-8*t**3*LOG(rh)**3 +  &
226       &  (-5.2663760117394626D-2*LOG(rh)**3)/x +  &
227       &  1.6496664658266762D0*LOG(rhoa/1.D6) +  &
228       &  -8.0809397859218401D-1*t*LOG(rhoa/1.D6) +  &
229       &  8.9302927091946642D-3*t**2*LOG(rhoa/1.D6) +  &
230       &  -1.9583649496497497D-5*t**3*LOG(rhoa/1.D6) +  &
231       &  (-8.9505572676891685D0*LOG(rhoa/1.D6))/x +  &
232       &  -3.0025283601622881D+1*LOG(rh)*LOG(rhoa/1.D6) +  &
233       &  3.0783365644763633D-1*t*LOG(rh)*LOG(rhoa/1.D6) +  &
234       &  -7.4521756337984706D-4*t**2*LOG(rh)*LOG(rhoa/1.D6) +  &
235       &  -5.7651433870681853D-7*t**3*LOG(rh)*LOG(rhoa/1.D6) +  &
236       &  (1.2872868529673207D0*LOG(rh)*LOG(rhoa/1.D6))/x +  &
237       &  -6.1739867501526535D-1*LOG(rh)**2*LOG(rhoa/1.D6) +  &
238       &  7.2347385705333975D-3*t*LOG(rh)**2*LOG(rhoa/1.D6) +  &
239       &  -3.0640494530822439D-5*t**2*LOG(rh)**2*LOG(rhoa/1.D6) +  &
240       &  6.5944609194346214D-8*t**3*LOG(rh)**2*LOG(rhoa/1.D6) +  &
241       &  (-2.8681650332461055D-2*LOG(rh)**2*LOG(rhoa/1.D6))/x +  &
242       &  6.5213802375160306D0*LOG(rhoa/1.D6)**2 +  &
243       &  -4.7907162004793016D-2*t*LOG(rhoa/1.D6)**2 +  &
244       &  -1.0727890114215117D-4*t**2*LOG(rhoa/1.D6)**2 +  &
245       &  5.6401818280534507D-7*t**3*LOG(rhoa/1.D6)**2 +  &
246       &  (5.4113070888923009D-1*LOG(rhoa/1.D6)**2)/x +  &
247       &  5.2062808476476330D-1*LOG(rh)*LOG(rhoa/1.D6)**2 +  &
248       &  -6.0696882500824584D-3*t*LOG(rh)*LOG(rhoa/1.D6)**2 +  &
249       &  2.3851383302608477D-5*t**2*LOG(rh)*LOG(rhoa/1.D6)**2 +  &
250       &  -1.5243837103067096D-8*t**3*LOG(rh)*LOG(rhoa/1.D6)**2 +  &
251       &  (-5.6543192378015687D-2*LOG(rh)*LOG(rhoa/1.D6)**2)/x +  &
252       &  -1.1630806410696815D-1*LOG(rhoa/1.D6)**3 +  &
253       &  1.3806404273119610D-3*t*LOG(rhoa/1.D6)**3 +  &
254       &  -2.0199865087650833D-6*t**2*LOG(rhoa/1.D6)**3 +  &
255       &  -3.0200284885763192D-9*t**3*LOG(rhoa/1.D6)**3 +  &
256       &  (-6.9425267104126316D-3*LOG(rhoa/1.D6)**3)/x
257
258  ntotal =-3.5863435141979573D-3 - 1.0098670235841110D-1 *t + 8.9741268319259721D-4 *t**2 - 1.4855098605195757D-6*t**3  &
259       &   - 1.2080330016937095D-1 /x + 1.1902674923928015D-3*LOG(rh) - 1.9211358507172177D-2*t*LOG(rh) +  &
260       &   2.4648094311204255D-4*t**2*LOG(rh) - 7.5641448594711666D-7*t**3*LOG(rh) +  &
261       &   (-2.0668639384228818D-02*LOG(rh))/x - 3.7593072011595188D-2*LOG(rh)**2 + 9.0993182774415718D-4 *t*LOG(rh)**2 +  &
262       &   -9.5698412164297149D-6*t**2*LOG(rh)**2 + 3.7163166416110421D-8*t**3*LOG(rh)**2 +  &
263       &   (1.1026579525210847D-2*LOG(rh)**2)/x + 1.1530844115561925D-2 *LOG(rh)**3 +  &
264       &   - 1.8083253906466668D-4 *t*LOG(rh)**3 + 8.0213604053330654D-7*t**2*LOG(rh)**3 +  &
265       &   -8.5797885383051337D-10*t**3*LOG(rh)**3 + (1.0243693899717402D-3*LOG(rh)**3)/x +  &
266       &   -1.7248695296299649D-2*LOG(rhoa/1.D6) + 1.1294004162437157D-2*t*LOG(rhoa/1.D6) +  &
267       &   -1.2283640163189278D-4*t**2*LOG(rhoa/1.D6) + 2.7391732258259009D-7*t**3*LOG(rhoa/1.D6) +  &
268       &   (6.8505583974029602D-2*LOG(rhoa/1.D6))/x +2.9750968179523635D-1*LOG(rh)*LOG(rhoa/1.D6) +  &
269       &   -3.6681154503992296D-3 *t*LOG(rh)*LOG(rhoa/1.D6) + 1.0636473034653114D-5*t**2*LOG(rh)*LOG(rhoa/1.D6) +  &
270       &   5.8687098466515866D-9*t**3*LOG(rh)*LOG(rhoa/1.D6) + (-5.2028866094191509D-3*LOG(rh)*LOG(rhoa/1.D6))/x +  &
271       &   7.6971988880587231D-4*LOG(rh)**2*LOG(rhoa/1.D6) - 2.4605575820433763D-5*t*LOG(rh)**2*LOG(rhoa/1.D6) +  &
272       &   2.3818484400893008D-7*t**2*LOG(rh)**2*LOG(rhoa/1.D6) +  &
273       &   -8.8474102392445200D-10*t**3*LOG(rh)**2*LOG(rhoa/1.D6) +  &
274       &   (-1.6640566678168968D-4*LOG(rh)**2*LOG(rhoa/1.D6))/x - 7.7390093776705471D-2*LOG(rhoa/1.D6)**2 +  &
275       &   5.8220163188828482D-4*t*LOG(rhoa/1.D6)**2 + 1.2291679321523287D-6*t**2*LOG(rhoa/1.D6)**2 +  &
276       &   -7.4690997508075749D-9*t**3*LOG(rhoa/1.D6)**2 + (-5.6357941220497648D-3*LOG(rhoa/1.D6)**2)/x +  &
277       &   -4.7170109625089768D-3*LOG(rh)*LOG(rhoa/1.D6)**2 + 6.9828868534370193D-5*t*LOG(rh)*LOG(rhoa/1.D6)**2 +  &
278       &   -3.1738912157036403D-7*t**2*LOG(rh)*LOG(rhoa/1.D6)**2 +  &
279       &   2.3975538706787416D-10*t**3*LOG(rh)*LOG(rhoa/1.D6)**2 +  &
280       &   (4.2304213386288567D-4*LOG(rh)*LOG(rhoa/1.D6)**2)/x + 1.3696520973423231D-3*LOG(rhoa/1.D6)**3 +  &
281       &   -1.6863387574788199D-5*t*LOG(rhoa/1.D6)**3 + 2.7959499278844516D-8*t**2*LOG(rhoa/1.D6)**3 +  &
282       &   3.9423927013227455D-11*t**3*LOG(rhoa/1.D6)**3 + (8.6136359966337272D-5*LOG(rhoa/1.D6)**3)/x
283
284  ntotal=EXP(ntotal)
285
286  rc=EXP(-22.378268374023630D0 + 0.44462953606125100D0*x + 0.33499495707849131D0*LOG(ntotal))
287
288  jnuc=EXP(jnuc)
289
290  if(jnuc < 1.D-7) then
291     print *,'Warning: nucleation rate < 1e-7/cm3s, using 0.0/cm3s,'
292     jnuc = 0.0D0
293  endif
294
295  jnuc = jnuc*1.D6 !1/(m^3s)
296
297  nac =x*ntotal
298
299  if (nac .lt. 1.D0) then
300     print *, 'Warning: number of acid molecules < 1 in nucleation regime, setting x=1'
301     x=1.D0
302  endif
303
304  RETURN
305
306END SUBROUTINE newnuklefit
307
308
309!*****************************************************************************
310SUBROUTINE heteronuk(TAIR,ra,mk0i,ACTOT)
311
312  use donnees
313  use free_param
314  IMPLICIT NONE
315
316  ! inputs/outputs
317  REAL, intent(in) :: TAIR, mk0i, ra
318  REAL, intent(out) :: ACTOT
319
320  ! Local variables
321  INTEGER :: i
322  INTEGER :: AER
323  REAL :: NDnuk(nbin,2)
324  REAL :: PSASAS
325
326  ! Initialization
327  AER = 2
328  ntot_aer = mk0i
329  ACTOT=0.D0
330
331  call logdist(sig_aer,ntot_aer,ra,rad_cld,NDnuk(:,AER)) ! NDnuk in #
332
333  do i=1,nbin
334     ! Saturation ratio with curvature effect
335     PSASAS = RHOsasat * EXP((2.D0*MSA*ST)/(RGAS*TAIR*RHOSA*rad_cld(i)))
336     SH2SO4 = PPSA / PSASAS
337   
338     if (SH2SO4 .ge. 1.D0 .and. NDnuk(i,AER) .gt. 0.) then
339        ! ACTOT = ACTOT + NDnuk(i,AER)/(vratio * rad_cld(i))
340        ACTOT = ACTOT + NDnuk(i,AER) ! Total AER number, #
341        NDnuk(i,AER) = 0.D0          ! Now, it's empty !
342     endif
343  enddo
344END SUBROUTINE heteronuk
Note: See TracBrowser for help on using the repository browser.