source: trunk/LMDZ.MARS/libf/phymars/coagulation_mod.F90 @ 4057

Last change on this file since 4057 was 4057, checked in by emillour, 8 days ago

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