source: trunk/LMDZ.MARS/libf/phymars/dust_coagulation_mod.F90

Last change on this file was 4064, checked in by jmauxion, 5 weeks ago

Mars PCM:
Fix some typos in callkeys_mod.F90 and dust_coagulation_mod.F90.
JM
File size: 24.3 KB
Line 
1module dust_coagulation_mod
2
3! ------------------------
4! Dust coagulation scheme
5! Based on Bertrand et al., 2022
6! Author: T. Bertrand
7! ------------------------
8
9use comcstfi_h, only: pi,g,r,mugaz
10use tracer_mod, only: igcm_dust_mass,igcm_dust_number, &
11                      igcm_stormdust_mass,igcm_stormdust_number, &
12                      igcm_topdust_mass,igcm_topdust_number
13use microphys_h, only: kbz,m0co2
14implicit none
15
16public ::  dust_coagulation_init
17
18logical,save :: full_coag_equations ! coagulation of dust particles using full equations instead of lookup tables
19logical,save :: dust_coag_kernel_b ! flag to activate Brownian coagulation
20logical,save :: dust_coag_kernel_g !                  gravitational
21logical,save :: dust_coag_kernel_de !                 diffusion enhancement
22logical,save :: dust_coag_kernel_ti ! Not implemented yet.
23!$OMP THREADPRIVATE(full_coag_equations,dust_coag_kernel_b,dust_coag_kernel_g,dust_coag_kernel_de,dust_coag_kernel_ti)
24integer,save :: coal_coeff_mode ! Flag to choose how to calculate coalescence coefficient
25!$OMP THREADPRIVATE(coal_coeff_mode)
26
27! -------------------------------
28! Coagulation Input Parameters
29! -------------------------------
30real    ::  effb=1.  !  Sticking efficiency for B coagulation
31real    ::  coal_fac=1.  ! Overall coalescence factor
32! Bin discretization of particle sizes
33real*8, save :: r1_coag = 0.01e-6 ! min radius
34real*8, save :: rn_coag = 40.e-6  ! max radius
35integer,parameter :: nres_coag = 10  ! nb bins
36real*8, save :: vrat_coag   ! volume ratio
37real*8, save, allocatable :: rads_coag(:),vols_coag(:),deltar_coag(:)
38
39!$OMP THREADPRIVATE(r1_coag,rn_coag,vrat_coag)
40!$OMP THREADPRIVATE(rads_coag,vols_coag,deltar_coag)
41
42! parameters of the look up tables of coagulation coefficients
43integer,parameter :: table_numt = 15   ! Nb of temperature bins
44integer,parameter :: table_nump = 25   ! Nb of pressure bins
45integer,parameter :: table_numm = 25   ! Nb of mass bins
46real, save :: table_pt(table_numt)
47real, save :: table_pres(table_nump)
48real, save :: table_mfp(table_numm)
49real, save :: table_b(table_numt,table_nump,nres_coag,nres_coag) ! Brownian
50real, save :: table_g(table_numt,table_nump,nres_coag,nres_coag) ! Gravitational
51real, save :: table_de(table_numt,table_nump,nres_coag,nres_coag) ! Brownian diffusion enhancement
52real, save :: table_ti(table_numt,table_nump,nres_coag,nres_coag) ! Turbulent
53real*8, public :: dev_dt = 0.63676    ! standard deviation of dust distribution
54real, public ::rho_dust=2500.  ! Mars dust density (kg.m-3)
55
56!$OMP THREADPRIVATE(table_pt,table_pres,table_mfp)
57!$OMP THREADPRIVATE(table_b,table_g,table_de,table_ti)
58
59contains
60
61!#######################################################################
62
63subroutine dust_coagulation_init()
64
65! Initialisation of the particle bin discretization
66
67integer  :: i
68allocate(rads_coag(nres_coag))
69allocate(vols_coag(nres_coag))
70allocate(deltar_coag(nres_coag))
71
72vrat_coag=(rn_coag/r1_coag)**(3./(nres_coag-1))
73
74do i = 1, nres_coag
75        rads_coag(i)  = r1_coag*vrat_coag**((i-1)/3.)
76        vols_coag(i)  = 4./3.*pi*r1_coag**3*vrat_coag**(i-1)
77enddo
78! diameter width
79deltar_coag(:)=2.*rads_coag(:)*2.**(1/3.)*(vrat_coag**(1./3.)-1)/(1+vrat_coag)**(1./3.)
80
81! Making Lookup Tables
82if (.not.full_coag_equations) then
83 call make_tables()
84endif
85!deallocate(rads_coag)
86!deallocate(vols_coag)
87!deallocate(deltar_coag)
88
89end subroutine dust_coagulation_init
90
91!#######################################################################
92
93SUBROUTINE dust_coagulation_main( &
94   ngrid, nlayer, nq, ptime, ptimestep, pq, pdqfi, pt, pdtfi, &
95   pplay, pplev, pdqcoag)
96
97!--------------------------------------------------------
98! Input Variables
99!--------------------------------------------------------
100INTEGER, INTENT(IN) :: ngrid ! number of horizontal grid points
101INTEGER, INTENT(IN) :: nlayer ! number of vertical grid points
102INTEGER, INTENT(IN) :: nq ! number of tracer species
103REAL, INTENT(IN) :: ptime
104REAL, INTENT(IN) :: ptimestep
105
106REAL, INTENT(IN) :: pq(ngrid,nlayer,nq) ! advected field nq
107REAL, INTENT(IN) :: pdqfi(ngrid,nlayer,nq)! tendancy field pq
108REAL, INTENT(IN) :: pt(ngrid,nlayer)      ! temperature at mid-layer (K)
109REAL, INTENT(IN) :: pdtfi(ngrid,nlayer)   ! tendancy temperature at mid-layer (K)
110REAL, INTENT(IN) :: pplay(ngrid,nlayer)     ! pressure at middle of the layers
111REAL, INTENT(IN) :: pplev(ngrid,nlayer+1) ! pressure at intermediate levels
112
113!--------------------------------------------------------
114! Output Variables
115!--------------------------------------------------------
116REAL, INTENT(OUT) :: pdqcoag(ngrid,nlayer,nq) ! tendancy field for dust coagulation
117
118!-----------------------------------------------------------------------
119!   Local variables
120!-----------------------------------------------------------------------
121integer  :: ie, je, id, jd, iq, i, j, k, l, ii, jj
122real dens,dev,cst,sig0
123real mfp0,pt0,rad0,rad10,rad20,rho0,pres0
124real term1,term2,kernel,norm
125real, dimension(ngrid, nlayer) :: r0, rho, rn, ntot, mtot, mtotnew
126real, dimension(ngrid, nlayer, nq) :: numb_new, numb_ini, numb_ini_dis, mass_ini, mass_new
127real, dimension(ngrid, nlayer, nres_coag) :: ndis, ndis_new
128real, dimension(ngrid,nlayer,nres_coag,nq) :: ndis_tab,rat_tab,ndis_tab_new
129real, dimension(ngrid, nlayer) :: zt
130real, dimension(ngrid, nlayer, nq) :: zq
131
132logical, save :: firstcall = .true.
133!$OMP THREADPRIVATE(firstcall)
134!for boxinterp function:
135integer :: t1,t2,p1,p2,m1,m2
136
137!======================================================================
138!   Main routine
139
140!-----------------------------------------------------------------------
141!*** 1) Initialisations
142!-----------------------------------------------------------------------
143
144!!! Particles properties
145dens    = rho_dust ! density dust
146sig0     = dev_dt    ! Standard deviation of the dust distribution
147cst     = .75 / (pi*dens) * exp( -4.5*sig0**2. )   ! fixed distribution
148
149!!! Radius and Volumes defined in coagul_init
150! rads_coag, vols_coag, deltar_coag, vrat_coag
151
152!!! Getting tracers and temperature fields
153zq(:,:,:)=pq(:,:,:)+pdqfi(:,:,:)*ptimestep
154zt(:,:)=pt(:,:)+pdtfi(:,:)*ptimestep
155
156!!! Atmospheric state
157DO l=1,nlayer
158   rho(:,l)=pplay(:,l)/(r*pt(:,l))
159ENDDO
160
161!!! Initializing tracers and other fields
162pdqcoag(:,:,:)=0.d0
163ndis_tab(:,:,:,:)=0.d0
164ndis_new(:,:,:)=0.d0
165mtot(:,:)=0.d0
166rat_tab(:,:,:,:) = 0.d0
167numb_ini_dis(:,:,:) = 0.d0
168numb_ini(:,:,:) = 0.d0
169numb_new(:,:,:) = 0.d0
170mass_ini(:,:,:) = 0.d0
171mass_new(:,:,:) = 0.d0
172
173!!! Preparing the initial log normal distribution
174do iq=1,nq
175     if ((iq.eq.igcm_dust_mass).or.(iq.eq.igcm_stormdust_mass).or.(iq.eq.igcm_topdust_mass)) then
176       mass_ini(:,:,iq)=max(zq(:,:,iq),1e-15)
177       numb_ini(:,:,iq)=max(zq(:,:,iq-1),1e-15)
178
179       !!! Get total mass to ensure mass conservation
180       mtot(:,:)=mtot(:,:)+mass_ini(:,:,iq)
181
182       !!! Initial Lognormal distribution
183       ! sig0 remains constant
184       r0(:,:)=(3/4.*mass_ini(:,:,iq)/(pi*numb_ini(:,:,iq)*dens))**(1/3.)*exp(-1.5*sig0**2)
185
186       do i=1,nres_coag ! Sum of the ndis for all dust modes
187          ndis_tab(:,:,i,iq)=numb_ini(:,:,iq)*rho(:,:)*deltar_coag(i)/(2.*rads_coag(i)*(2*pi)**0.5*sig0)*exp(-0.5*(log(rads_coag(i)/r0(:,:)))**2/sig0**2)
188       enddo
189       !!! Initial number mixing ratio from bin discretization nb/kg
190       numb_ini_dis(:,:,iq)=sum(ndis_tab(:,:,:,iq),3)/rho(:,:)
191       numb_new(:,:,iq)=numb_ini_dis(:,:,iq)
192     endif
193enddo
194
195
196!!! Defining initial number ratio for each mode
197do iq=1,nq
198     if ((iq.eq.igcm_dust_mass).or.(iq.eq.igcm_stormdust_mass).or.(iq.eq.igcm_topdust_mass)) then
199       rat_tab(:,:,:,iq)=ndis_tab(:,:,:,iq)/sum(ndis_tab(:,:,:,:),4)
200     endif
201enddo
202
203!!! Normalizing the ratio (make sure sum=1)
204do iq=1,nq
205     if ((iq.eq.igcm_dust_mass).or.(iq.eq.igcm_stormdust_mass).or.(iq.eq.igcm_topdust_mass)) then
206       rat_tab(:,:,:,iq)=rat_tab(:,:,:,iq)/sum(rat_tab(:,:,:,:),4)
207   endif
208enddo
209
210!!! Total initial distribution nb/m3 per bin
211ndis(:,:,:)=sum(ndis_tab(:,:,:,:),4)
212!!! Total initial distribution nb/kg
213ntot(:,:)=sum(ndis(:,:,:),3)/rho(:,:)
214
215!-----------------------------------------------------------------------
216!*** 2) Coagulation
217!-----------------------------------------------------------------------
218! Spatial loop => 1D
219
220if (full_coag_equations) then ! computing full coagulation equations
221 do i=1,ngrid
222  do j=1,nlayer
223      if (ntot(i,j).gt.1000.) then ! activating coagulation above this threshold
224       pt0=pt(i,j)
225       rho0=rho(i,j)
226       mfp0=mfp(pt0,rho0)
227
228       ! Bin loop
229       do k=1,nres_coag
230
231         ! Term 1
232         term1=0.
233         do jj=1,k
234         do ii=1,k-1
235             kernel=0.
236             if (dust_coag_kernel_b) kernel=kernel+betab(pt0,mfp0,rads_coag(ii),rads_coag(jj))
237             if (dust_coag_kernel_g) kernel=kernel+betag(pt0,mfp0,rads_coag(ii),rads_coag(jj))
238             if (dust_coag_kernel_de) kernel=kernel+betade(pt0,mfp0,rho0,rads_coag(ii),rads_coag(jj))
239             if (dust_coag_kernel_ti) kernel=kernel+betati(pt0,mfp0,rho0,rads_coag(ii),rads_coag(jj))
240             term1=term1+frac(ii,jj,k)*vols_coag(ii)*kernel*ndis_new(i,j,ii)*ndis(i,j,jj)
241          enddo
242         enddo
243         term1=term1*ptimestep/vols_coag(k)
244
245         ! Term 2
246         term2=0.
247         do jj=1,nres_coag
248            kernel=0.
249            if (dust_coag_kernel_b) kernel=kernel+betab(pt0,mfp0,rads_coag(k),rads_coag(jj))
250            if (dust_coag_kernel_g) kernel=kernel+betag(pt0,mfp0,rads_coag(k),rads_coag(jj))
251            if (dust_coag_kernel_de) kernel=kernel+betade(pt0,mfp0,rho0,rads_coag(k),rads_coag(jj))
252            if (dust_coag_kernel_ti) kernel=kernel+betati(pt0,mfp0,rho0,rads_coag(k),rads_coag(jj))
253            term2=term2+(1-frac(k,jj,k))*kernel*ndis(i,j,jj)
254         enddo
255         term2=term2*ptimestep
256
257         ndis_new(i,j,k)=(ndis(i,j,k)+coal_fac*term1)/(1.+coal_fac*term2)
258
259       enddo
260
261       norm=sum(ndis(i,j,:)*vols_coag(:))/sum(ndis_new(i,j,:)*vols_coag(:))
262       ndis_new(i,j,:)=ndis_new(i,j,:)*norm
263
264      endif
265  enddo
266 enddo
267
268else ! full_coag_equations
269
270 ! Using Lookup tables
271 do i=1,ngrid
272  do j=1,nlayer
273      if (ntot(i,j).gt.1000.) then ! activating coagulation above this threshold
274       pt0=pt(i,j)
275       rho0=rho(i,j)
276       mfp0=mfp(pt0,rho0)
277       pres0=pplay(i,j)
278       !find the corners for interpolation       
279       t2=findval(table_pt,pt0)
280       if (t2==0) t2=table_numt
281       if (t2==1) t2=2
282       t1=t2-1
283       m2=findval(table_mfp,mfp0)
284       if (m2==0) m2=table_numm
285       if (m2==1) m2=2
286       m1=m2-1
287       p2=findval(table_pres,pres0)
288       if (p2==0) p2=table_numt
289       if (p2==1) p2=2
290       p1=p2-1
291
292       ! Bin loop
293       do k=1,nres_coag
294
295         ! Term 1
296         term1=0.
297         do jj=1,k
298          do ii=1,k-1
299             kernel=0.
300             if (dust_coag_kernel_b) kernel=kernel+boxinterp(table_b(t1,m1,ii,jj),       &   
301                                                   table_b(t1,m2,ii,jj),       &
302                                                   table_b(t2,m1,ii,jj),       &
303                                                   table_b(t2,m2,ii,jj),       &
304                                                   table_pt(t1),table_pt(t2),log10(table_mfp(m1)),log10(table_mfp(m2)),pt0,log10(mfp0))
305             if (dust_coag_kernel_g) kernel=kernel+boxinterp(table_g(t1,m1,ii,jj),       &
306                                                   table_g(t1,m2,ii,jj),       &
307                                                   table_g(t2,m1,ii,jj),       &
308                                                   table_g(t2,m2,ii,jj),       &
309                                                   table_pt(t1),table_pt(t2),log10(table_mfp(m1)),log10(table_mfp(m2)),pt0,log10(mfp0))
310             if (dust_coag_kernel_de) kernel=kernel+boxinterp(table_de(t1,p1,ii,jj),    &
311                                                   table_de(t1,p2,ii,jj),     &
312                                                   table_de(t2,p1,ii,jj),     &
313                                                   table_de(t2,p2,ii,jj),     &
314                                                   table_pt(t1),table_pt(t2),log10(table_pres(p1)),log10(table_pres(p2)),pt0,log10(pres0))
315             ! if (dust_coag_kernel_ti) kernel=kernel+betati(pt0,mfp0,rho(i,j,l),ii,jj)  !not implemented yet
316             term1=term1+frac(ii,jj,k)*vols_coag(ii)*kernel*ndis_new(i,j,ii)*ndis(i,j,jj)
317          enddo
318         enddo
319         term1=term1*ptimestep/vols_coag(k)
320
321
322         ! Term 2
323         term2=0.
324         do jj=1,nres_coag
325            kernel=0.
326            if (dust_coag_kernel_b) kernel=kernel+boxinterp(table_b(t1,m1,k,jj),       &
327                                                  table_b(t1,m2,k,jj),       &
328                                                  table_b(t2,m1,k,jj),       &
329                                                  table_b(t2,m2,k,jj),       &
330                                                  table_pt(t1),table_pt(t2),log10(table_mfp(m1)),log10(table_mfp(m2)),pt0,log10(mfp0))
331             if (dust_coag_kernel_g) kernel=kernel+boxinterp(table_g(t1,m1,k,jj),       &
332                                                   table_g(t1,m2,k,jj),       &
333                                                   table_g(t2,m1,k,jj),       &
334                                                   table_g(t2,m2,k,jj),       &
335                                                   table_pt(t1),table_pt(t2),log10(table_mfp(m1)),log10(table_mfp(m2)),pt0,log10(mfp0))
336             if (dust_coag_kernel_de) kernel=kernel+boxinterp(table_de(t1,p1,k,jj),       &
337                                                   table_de(t1,p2,k,jj),       &
338                                                   table_de(t2,p1,k,jj),       &
339                                                   table_de(t2,p2,k,jj),       &
340                                                   table_pt(t1),table_pt(t2),log10(table_pres(p1)),log10(table_pres(p2)),pt0,log10(pres0))
341            ! if (dust_coag_kernel_ti) kernel=kernel+betati(pt0,mfp0,rho(i,j,l),k,jj)  !not implemented yet
342            term2=term2+(1-frac(k,jj,k))*kernel*ndis(i,j,jj)
343         enddo
344         term2=term2*ptimestep
345
346         ndis_new(i,j,k)=(ndis(i,j,k)+coal_fac*term1)/(1.+coal_fac*term2)
347       enddo
348     
349      !! Volume conservation
350      norm=sum(ndis(i,j,:)*vols_coag(:))/sum(ndis_new(i,j,:)*vols_coag(:))
351      ndis_new(i,j,:)=ndis_new(i,j,:)*norm
352
353      endif
354
355  enddo
356 enddo
357endif ! full_coag_equations
358
359!-----------------------------------------------------------------------
360!*** 3) new moments mass and number
361!-----------------------------------------------------------------------
362mtotnew(:,:)=1.d-15
363
364do iq=1,nq
365     if ((iq.eq.igcm_dust_mass).or.(iq.eq.igcm_stormdust_mass).or.(iq.eq.igcm_topdust_mass)) then
366      !! Estimate new distribution in each dust mode using the initial number ratios
367      ndis_tab_new(:,:,:,iq)=ndis_new(:,:,:)*rat_tab(:,:,:,iq)     
368
369      do i=1,ngrid
370       do j=1,nlayer
371         !! New number and mass mixing ratio in each mode
372         numb_new(i,j,iq)=sum(ndis_tab_new(i,j,:,iq))/rho(i,j)
373         mass_new(i,j,iq)=sum(ndis_tab_new(i,j,:,iq)*vols_coag(:))*dens/rho(i,j)
374       enddo
375      enddo
376
377      !! New total mass
378      mtotnew(:,:)=mtotnew(:,:)+mass_new(:,:,iq)
379   endif
380enddo
381
382do iq=1,nq
383     if ((iq.eq.igcm_dust_mass).or.(iq.eq.igcm_stormdust_mass).or.(iq.eq.igcm_topdust_mass)) then
384      !! mass conservation
385      mass_new(:,:,iq)=mass_new(:,:,iq)/mtotnew(:,:)*mtot(:,:)
386      !! New tendancy
387      pdqcoag(:,:,iq)=(mass_new(:,:,iq)-mass_ini(:,:,iq))/ptimestep
388      !pdqcoag(:,:,iq-1)=(numb_new(:,:,iq)-numb_ini(:,:,iq))/ptimestep
389      pdqcoag(:,:,iq-1)=(numb_new(:,:,iq)-numb_ini_dis(:,:,iq))/ptimestep
390      !! check
391      where (zq(:,:,iq)+pdqcoag(:,:,iq)*ptimestep.le.0.)
392        pdqcoag(:,:,iq)=-zq(:,:,iq)/ptimestep+1.e-15
393      end where
394      where (zq(:,:,iq-1)+pdqcoag(:,:,iq-1)*ptimestep.le.0.)
395        pdqcoag(:,:,iq-1)=-zq(:,:,iq-1)/ptimestep+1.e-14
396      end where
397   endif
398enddo
399
400end subroutine dust_coagulation_main
401
402
403!======================================================================
404!======================================================================
405real function dynvis(pt)
406  implicit none
407  ! dynamic viscosity following sutherland's formula for CO2
408  ! Pa.s = kg m-1 s-1
409  real, intent(in) :: pt
410  dynvis=1.37e-5*(273.15 + 222.)/(pt + 222.)*(pt/273.15)**(3./2.)
411end function dynvis
412 
413real function kinvis(pt,rho)
414  implicit none
415  real, intent(in) :: pt,rho
416  kinvis=dynvis(pt)/rho
417end function kinvis
418
419real function thervel(pt,massm)
420  implicit none
421  ! thermal velocity of an air molecule or a dust particle
422  ! m s-1
423  real, intent(in) :: pt, massm
424  thervel=(8*kbz*pt/(pi*massm))**0.5
425end function thervel
426
427real function mfp(pt,rho)
428  implicit none
429  ! thermal velocity of an air molecule
430  ! m s-1
431  real, intent(in) :: pt, rho
432  mfp=2.*kinvis(pt,rho)/thervel(pt,m0co2) ! mean free path (m)
433end function mfp
434
435real function cun(pt,mfp,rad)
436  implicit none
437  !! Cunningham slip flow correction
438  real, intent(in) :: pt, mfp, rad
439  cun=1.+mfp/rad*(1.246+0.42*exp(-0.87/(mfp/rad)))
440end function cun
441
442real function diff(pt,mfp,rad)
443  implicit none
444  !! diffusion coefficient
445  real, intent(in) :: pt, mfp, rad
446  diff=kbz*pt/(6.*pi*rad*dynvis(pt))*cun(pt,mfp,rad)
447end function diff
448
449real function lambda(pt, mfp, rad)
450  implicit none
451  !! particle mean free path
452  real, intent(in) :: pt,mfp,rad
453  real massm,dens
454  dens    = rho_dust  ! density dust
455  massm=dens*4./3.*pi*rad**3
456  lambda=2.*diff(pt, mfp, rad)/(pi*thervel(pt,massm))
457end function lambda
458
459real function delta(pt,mfp,rad)
460  implicit none
461  !! mean distance
462  real, intent(in) :: pt,mfp,rad
463  real :: ltmp
464
465  ltmp=lambda(pt,mfp,rad)
466  delta=( (2.*rad+ltmp)**3 - (4.*rad**2+ltmp**2)**1.5 ) /(6*rad*ltmp) - 2*rad
467end function delta
468
469real function fallv(pt,mfp,rad)
470  implicit none
471  !! fall velocity
472  real, intent(in) :: pt,mfp,rad
473  real dens
474  dens    = rho_dust  ! density dust
475  fallv=2./9.*rad**2*dens*g/dynvis(pt)*cun(pt,mfp,rad)
476end function fallv
477
478real function reyn(pt,mfp,rad,rho)
479  implicit none
480  !! Reynold number
481  real, intent(in) :: pt,rad,mfp,rho
482  reyn=2.*rad*fallv(pt,mfp,rad)/kinvis(pt,rho)
483end function reyn
484
485real function stokes(pt,mfp,rad1,rad2)
486  implicit none
487  !! Stokes number
488  real, intent(in) :: pt,mfp,rad1,rad2
489  real :: f1
490
491  f1=fallv(pt,mfp,rad1)
492  stokes=f1*abs(fallv(pt,mfp,rad2)-f1)/(rad2*g)
493end function stokes
494
495real function schmidt(pt,mfp,rho,rad)
496  implicit none
497  !! Schmidt number
498  real, intent(in) :: pt,mfp,rho,rad
499  schmidt=kinvis(pt,rho)/diff(pt,mfp,rad)
500end function schmidt
501
502real function coal(pt,mfp,rad1,rad2,rho)
503  implicit none
504  !! Coalescence efficiency
505  real, intent(in) :: pt,mfp,rad1,rad2,rho
506  real :: coalea,coalev
507  real :: stmp,rtmp
508
509  stmp=stokes(pt,mfp,rad1,rad2)
510  rtmp=reyn(pt,mfp,rad2,rho)
511  coalea=stmp**2/(stmp+0.5)**2
512  if (stmp>1.214) then
513    coalev=(1+0.75*log(2.*stmp)/(stmp-1.214))**(-2)
514  else
515    coalev=0.
516  endif
517  coal=(60.*coalev+coalea*rtmp)/(60.+rtmp)
518end function coal
519
520real function betab(pt,mfp,rad1,rad2)
521  implicit none
522  real, intent(in) :: pt,mfp,rad1,rad2
523  real :: m1,m2,num,den1,den2
524  real :: dens
525  real :: d1,d2
526
527  dens    = rho_dust  ! density dust
528  d1=diff(pt,mfp,rad1)
529  d2=diff(pt,mfp,rad2)
530  !! kernel transition regime
531  ! Mass of 1 particle 1 and 2
532  m1=dens*4./3.*pi*rad1**3
533  m2=dens*4./3.*pi*rad2**3
534  num=4.*pi*(rad1+rad2)*(d1+d2)
535  den1=(rad1+rad2)/(rad1+rad2+(delta(pt,mfp,rad1)**2+delta(pt,mfp,rad2)**2)**0.5)
536  den2=4./effb*(d1+d2)/((thervel(pt,m1)**2+thervel(pt,m2)**2)**0.5*(rad1+rad2))
537  betab=num/(den1+den2)
538end function betab
539
540real function betag(pt,mfp,rad1,rad2)
541  implicit none
542  real, intent(in) :: pt,mfp,rad1,rad2
543  real w1,w2,ee
544  real :: radt
545  !! kernel gravitation
546  ! fall vel
547  w1=fallv(pt,mfp,rad1)
548  w2=fallv(pt,mfp,rad2)
549  radt=(rad1+rad2)**2
550  if (coal_coeff_mode.eq.0) then ! E=1
551    ee=1.
552  elseif (coal_coeff_mode.eq.1) then
553    ee=1.5*(min(rad1,rad2))**2/radt
554  elseif (coal_coeff_mode.eq.2) then
555    ee=0.25*(min(rad1,rad2))**2/radt
556  endif
557  !b=coal(pt,r1,r2)*pi*(r1+r2)**2*abs(w1-w2)
558  betag=ee*pi*radt*abs(w1-w2)
559end function betag
560
561real function betade(pt,mfp,rho,rad1,rad2)
562  implicit none
563  real, intent(in) :: pt,mfp,rad1,rad2,rho
564  real rd,sc
565  !! kernel diffusion enhancement
566  rd=reyn(pt,mfp,max(rad1,rad2),rho)
567  sc=schmidt(pt,mfp,rho,min(rad1,rad2))
568  if (rd.le.1.) then
569    betade=betab(pt,mfp,rad1,rad2)*0.45*rd**(1/3.)*sc**(1/3.)
570  else
571    betade=betab(pt,mfp,rad1,rad2)*0.45*rd**(1/2.)*sc**(1/3.)
572  endif
573end function betade
574
575real function betade_tabfunc(pt,pres,rad1,rad2)
576  implicit none
577  real, intent(in) :: pt,pres,rad1,rad2
578  real rd,sc,rho,mf
579  !! kernel diffusion enhancement
580  rho=pres/(pt*r)
581  mf=mfp(pt,rho)
582  rd=reyn(pt,mf,max(rad1,rad2),rho)
583  sc=schmidt(pt,mf,rho,min(rad1,rad2))
584  if (rd.le.1.) then
585    betade_tabfunc=betab(pt,mf,rad1,rad2)*0.45*rd**(1/3.)*sc**(1/3.)
586  else
587    betade_tabfunc=betab(pt,mf,rad1,rad2)*0.45*rd**(1/2.)*sc**(1/3.)
588  endif
589end function betade_tabfunc
590
591real function betati(pt,mfp,rho,rad1,rad2)
592  implicit none
593  real, intent(in) :: pt,mfp,rad1,rad2,rho
594  real w1,w2,eps
595  ! fall vel
596  w1=fallv(pt,mfp,rad1)
597  w2=fallv(pt,mfp,rad2)
598  betati=eps**(0.75)*pi/(g*kinvis(pt,rho)**(0.25))*(rad1+rad2)**2*abs(w1-w2)
599end function betati
600
601real function betati_tabfunc(pt,pres,rad1,rad2)
602  implicit none
603  real, intent(in) :: pt,rad1,rad2
604  real w1,w2,eps,mf,rho,pres
605  rho=pres/(pt*r)
606  mf=mfp(pt,rho)
607  ! fall vel
608  w1=fallv(pt,mf,rad1)
609  w2=fallv(pt,mf,rad2)
610  betati_tabfunc=eps**(0.75)*pi/(g*kinvis(pt,rho)**(0.25))*(rad1+rad2)**2*abs(w1-w2)
611end function betati_tabfunc
612
613real function frac(i,j,k)
614  implicit none
615  integer, intent(in) :: i,j,k
616  real vint
617  ! intermediate volume
618  vint=vols_coag(i)+vols_coag(j)
619  if (k.lt.nres_coag) then
620     if ( (vint.ge.vols_coag(k)) .and. (vint.lt.vols_coag(k+1))) then
621         frac=(vols_coag(k+1)-vint)/(vols_coag(k+1)-vols_coag(k))*vols_coag(k)/vint
622         return
623     endif
624  endif
625  if ((k.eq.nres_coag) .and. (vint.ge.vols_coag(k))) then
626     frac=1.
627     return
628  endif
629  if (k.gt.1) then
630     if ( (vint.lt.vols_coag(k)) .and. (vint.gt.vols_coag(k-1))) then
631        frac=(-vols_coag(k-1)+vint)/(-vols_coag(k-1)+vols_coag(k))*vols_coag(k)/vint
632        return
633     endif
634  endif
635
636  frac=0.
637  return
638end function frac
639
640! *********************************************************
641! Building lookup tables : called in the firstcall
642! *********************************************************
643
644subroutine make_tables()
645
646implicit none
647integer :: i,j,k,l
648real :: t1,t2,p1,p2
649
650do i = 1,table_numt
651    table_pt(i)=50.+25.*real(i)
652end do
653
654do j = 1,table_nump
655    table_pres(j)=(1.e-8)*10.**(real(j)/2.)
656end do
657
658do j = 1,table_numm
659    table_mfp(j)=(1.e-8)*10.**(real(j)/2.)
660end do
661
662do i = 1,table_numt
663    do j = 1,table_numm
664        do k = 1,nres_coag
665            do l = 1,nres_coag
666                table_b(i,j,k,l)=betab(table_pt(i),table_mfp(j),rads_coag(k),rads_coag(l))
667                table_g(i,j,k,l)=betag(table_pt(i),table_mfp(j),rads_coag(k),rads_coag(l))
668            end do
669        end do
670    end do
671end do
672
673do i = 1,table_numt
674    do j = 1,table_nump
675        do k = 1,nres_coag
676            do l = 1,nres_coag
677                table_de(i,j,k,l)=betade_tabfunc(table_pt(i),table_pres(j),rads_coag(k),rads_coag(l))
678            end do
679        end do
680    end do
681end do
682
683end subroutine make_tables
684
685!======================================================================C
686!======================================================================C
687
688real function boxinterp(ft1p1,ft2p1,ft1p2,ft2p2,tref1,tref2,pref1,  &
689         pref2,tmid,pmid)
690!
691!   Calculate 2d interpolation from pressure and temperature
692!----------------------------------------------------------------------!
693!                   T2 .FT2P1                    .FT2P2
694!
695!
696!                   T1 .FT1P1                    .FT1P2
697!                      P1                        P2
698!----------------------------------------------------------------------!
699
700implicit none
701
702real :: ft1p1, ft2p1, ft1p2, ft2p2, tref1, tref2
703real :: pref1, pref2, tmid, pmid, ans1, ans2, ans
704
705!======================================================================C
706
707ans1 = ft1p1 + (ft2p1 - ft1p1)*(tmid - tref1)/(tref2 - tref1)
708ans2 = ft1p2 + (ft2p2 - ft1p2)*(tmid - tref1)/(tref2 - tref1)
709boxinterp  = ans1 + (ans2 - ans1)*(pmid - pref1)/(pref2 - pref1)
710
711end function boxinterp
712
713!======================================================================C
714!======================================================================C
715
716integer function findval(array,value)
717implicit none
718real :: array(:)
719real :: value
720
721findval=minloc(array,dim=1,mask = (array > value))
722
723end function findval
724
725end module dust_coagulation_mod
726
Note: See TracBrowser for help on using the repository browser.