Ignore:
Timestamp:
Jul 13, 2017, 9:33:32 PM (7 years ago)
Author:
oboucher
Message:

Adding enhancement factor to coagulation kernels due to van der Waals forces.
Two options are possible, default is no enhancement factor.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • LMDZ5/trunk/libf/phylmd/StratAer/coagulate.F90

    r2690 r2949  
    6363  REAL                                          :: denom
    6464
     65  ! Additional variables for coagulation enhancement factor due to van der Waals forces
     66  ! Taken from Chan and Mozurkewich, Measurement of the coagulation rate constant for sulfuric acid
     67  ! particles as a function of particle size using TDMA, Aerosol Science, 32, 321-339, 2001.
     68  !--ok_vdw is 0 for no vdW forces, 1 for E(0), 2 for E(infinity)
     69  INTEGER, PARAMETER                            :: ok_vdw = 0
     70  REAL, PARAMETER                               :: avdW1 = 0.0757
     71  REAL, PARAMETER                               :: avdW3 = 0.0015
     72  REAL, PARAMETER                               :: bvdW0 = 0.0151
     73  REAL, PARAMETER                               :: bvdW1 = -0.186
     74  REAL, PARAMETER                               :: bvdW3 = -0.0163
     75  REAL, PARAMETER                               :: AvdW = 6.4e-20 !Hamaker constant in J = 1e7 erg
     76  REAL                                          :: AvdWi
     77  REAL                                          :: xvdW
     78  REAL                                          :: EvdW
     79
    6580  DO i=1, nbtr_bin
    6681   radius(i)=mdw(i)/2.
     
    156171  ENDDO
    157172
    158 !--pre-compute the beta(i,j) from equation 17
     173!--pre-compute the beta(i,j) from equation 17 in Jacobson
    159174  num=0.0
    160175  DO j=1, nbtr_bin
    161176  DO i=1, nbtr_bin
     177!
    162178   num=4.*RPI*(radius(i)+radius(j))*(Di(i)+Di(j))
    163179   denom=(radius(i)+radius(j))/(radius(i)+radius(j)+sqrt(delta(i)**2.+delta(j)**2.))+ &
    164180        & 4.*(Di(i)+Di(j))/(sqrt(thvelpar(i)**2.+thvelpar(j)**2.)*(radius(i)+radius(j)))
    165181   beta(i,j)=num/denom
     182!
     183!--compute enhancement factor due to van der Waals forces
     184   IF (ok_vdw .EQ. 0) THEN      !--no enhancement factor
     185     Evdw=1.0
     186   ELSEIF (ok_vdw .EQ. 1) THEN  !--E(0) case
     187     AvdWi = AvdW/(k_B*t_seri(ilon,ilev))*(4.*radius(i)*radius(j))/(radius(i)+radius(j))**2.
     188     xvdW = LOG(1.+AvdWi)
     189     EvdW = 1. + avdW1*xvdW + avdW3*xvdW**3
     190   ELSEIF (ok_vdw .EQ. 2) THEN  !--E(infinity) case
     191     AvdWi = AvdW/(k_B*t_seri(ilon,ilev))*(4.*radius(i)*radius(j))/(radius(i)+radius(j))**2.
     192     xvdW = LOG(1.+AvdWi)
     193     EvdW = 1. + SQRT(AvdWi/3.)/(1.+bvdW0*SQRT(AvdWi)) + bvdW1*xvdW + bvdW3*xvdW**3.
     194   ENDIF
     195!
     196   beta(i,j)=beta(i,j)*EvdW
     197
    166198  ENDDO
    167199  ENDDO
Note: See TracChangeset for help on using the changeset viewer.