Changeset 2949 for LMDZ5/trunk/libf/phylmd

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:

• LMDZ5/trunk/libf/phylmd/StratAer/coagulate.F90 (modified) (2 diffs)

LMDZ5/trunk/libf/phylmd/StratAer/coagulate.F90

@@ -63 +63 @@
   REAL                                          :: denom
 
+  ! Additional variables for coagulation enhancement factor due to van der Waals forces
+  ! Taken from Chan and Mozurkewich, Measurement of the coagulation rate constant for sulfuric acid 
+  ! particles as a function of particle size using TDMA, Aerosol Science, 32, 321-339, 2001.
+  !--ok_vdw is 0 for no vdW forces, 1 for E(0), 2 for E(infinity)
+  INTEGER, PARAMETER                            :: ok_vdw = 0
+  REAL, PARAMETER                               :: avdW1 = 0.0757
+  REAL, PARAMETER                               :: avdW3 = 0.0015
+  REAL, PARAMETER                               :: bvdW0 = 0.0151
+  REAL, PARAMETER                               :: bvdW1 = -0.186
+  REAL, PARAMETER                               :: bvdW3 = -0.0163
+  REAL, PARAMETER                               :: AvdW = 6.4e-20 !Hamaker constant in J = 1e7 erg
+  REAL                                          :: AvdWi
+  REAL                                          :: xvdW
+  REAL                                          :: EvdW
+
   DO i=1, nbtr_bin
    radius(i)=mdw(i)/2.
@@ -156 +171 @@
   ENDDO
 
-!--pre-compute the beta(i,j) from equation 17
+!--pre-compute the beta(i,j) from equation 17 in Jacobson
   num=0.0
   DO j=1, nbtr_bin
   DO i=1, nbtr_bin
+!
    num=4.*RPI*(radius(i)+radius(j))*(Di(i)+Di(j))
    denom=(radius(i)+radius(j))/(radius(i)+radius(j)+sqrt(delta(i)**2.+delta(j)**2.))+ &
         & 4.*(Di(i)+Di(j))/(sqrt(thvelpar(i)**2.+thvelpar(j)**2.)*(radius(i)+radius(j)))
    beta(i,j)=num/denom
+!
+!--compute enhancement factor due to van der Waals forces
+   IF (ok_vdw .EQ. 0) THEN      !--no enhancement factor
+     Evdw=1.0
+   ELSEIF (ok_vdw .EQ. 1) THEN  !--E(0) case
+     AvdWi = AvdW/(k_B*t_seri(ilon,ilev))*(4.*radius(i)*radius(j))/(radius(i)+radius(j))**2.
+     xvdW = LOG(1.+AvdWi)
+     EvdW = 1. + avdW1*xvdW + avdW3*xvdW**3
+   ELSEIF (ok_vdw .EQ. 2) THEN  !--E(infinity) case
+     AvdWi = AvdW/(k_B*t_seri(ilon,ilev))*(4.*radius(i)*radius(j))/(radius(i)+radius(j))**2.
+     xvdW = LOG(1.+AvdWi)
+     EvdW = 1. + SQRT(AvdWi/3.)/(1.+bvdW0*SQRT(AvdWi)) + bvdW1*xvdW + bvdW3*xvdW**3.
+   ENDIF
+!
+   beta(i,j)=beta(i,j)*EvdW
+
   ENDDO
   ENDDO