source: LMDZ6/trunk/libf/phylmd/StratAer/coagulate.F90 @ 4738

Last change on this file since 4738 was 3677, checked in by oboucher, 5 years ago

Changed the way to initialise nbtr_bin and other dimensions and indices
in the StratAer? module based on infotrac_phy rather than infotrac.

Also added a missing $OMP THREADPRIVATE(nqperes)
  • Property svn:keywords set to Id
File size: 9.2 KB
Line 
1!
2! $Id: coagulate.F90 3677 2020-05-06 15:18:32Z oboucher $
3!
4SUBROUTINE COAGULATE(pdtcoag,mdw,tr_seri,t_seri,pplay,dens_aer,is_strato)
5  !     -----------------------------------------------------------------------
6  !
7  !     Author : Christoph Kleinschmitt (with Olivier Boucher)
8  !     ------
9  !
10  !     purpose
11  !     -------
12  !
13  !     interface
14  !     ---------
15  !      input
16  !       pdtphys        time step duration                 [sec]
17  !       tr_seri        tracer mixing ratios               [kg/kg]
18  !       mdw             # or mass median diameter          [m]
19  !
20  !     method
21  !     ------
22  !
23  !     -----------------------------------------------------------------------
24
25  USE dimphy, ONLY : klon,klev
26  USE aerophys
27  USE infotrac_phy
28  USE phys_local_var_mod, ONLY: DENSO4, f_r_wet
29  USE YOMCST
30
31  IMPLICIT NONE
32
33  !--------------------------------------------------------
34
35  ! transfer variables when calling this routine
36  REAL,INTENT(IN)                               :: pdtcoag ! Time step in coagulation routine [s]
37  REAL,DIMENSION(nbtr_bin),INTENT(IN)           :: mdw     ! aerosol particle diameter in each bin [m]
38  REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT)  :: tr_seri ! Concentration Traceur [U/KgA]
39  REAL,DIMENSION(klon,klev),INTENT(IN)          :: t_seri  ! Temperature
40  REAL,DIMENSION(klon,klev),INTENT(IN)          :: pplay   ! pression pour le mileu de chaque couche (en Pa)
41  REAL,DIMENSION(klon,klev)                     :: dens_aer! density of aerosol [kg/m3 aerosol] with default H2SO4 mass
42  LOGICAL,DIMENSION(klon,klev),INTENT(IN)       :: is_strato
43
44  ! local variables in coagulation routine
45  INTEGER                                       :: i,j,k,nb,ilon,ilev
46  REAL, DIMENSION(nbtr_bin)                     :: radius ! aerosol particle radius in each bin [m]
47  REAL, DIMENSION(klon,klev,nbtr_bin)           :: tr_t ! Concentration Traceur at time t [U/KgA]
48  REAL, DIMENSION(klon,klev,nbtr_bin)           :: tr_tp1 ! Concentration Traceur at time t+1 [U/KgA]
49  REAL, DIMENSION(nbtr_bin,nbtr_bin,nbtr_bin)   :: ff   ! Volume fraction of intermediate particles
50  REAL, DIMENSION(nbtr_bin)                     :: V    ! Volume of bins
51  REAL, DIMENSION(nbtr_bin,nbtr_bin)            :: Vij  ! Volume sum of i and j
52  REAL                                          :: eta  ! Dynamic viscosity of air
53  REAL, PARAMETER                               :: mair=4.8097E-26 ! Average mass of an air molecule [kg]
54  REAL                                          :: zrho ! Density of air
55  REAL                                          :: mnfrpth ! Mean free path of air
56  REAL, DIMENSION(nbtr_bin)                     :: Kn   ! Knudsen number of particle i
57  REAL, DIMENSION(nbtr_bin)                     :: Di   ! Particle diffusion coefficient
58  REAL, DIMENSION(nbtr_bin)                     :: m_par   ! Mass of particle i
59  REAL, DIMENSION(nbtr_bin)                     :: thvelpar! Thermal velocity of particle i
60  REAL, DIMENSION(nbtr_bin)                     :: mfppar  ! Mean free path of particle i
61  REAL, DIMENSION(nbtr_bin)                     :: delta! delta of particle i (from equation 21)
62  REAL, DIMENSION(nbtr_bin,nbtr_bin)            :: beta ! Coagulation kernel from Brownian diffusion
63  REAL                                          :: beta_const ! Constant coagulation kernel (for comparison)
64  REAL                                          :: num
65  REAL                                          :: numi
66  REAL                                          :: denom
67
68  ! Additional variables for coagulation enhancement factor due to van der Waals forces
69  ! Taken from Chan and Mozurkewich, Measurement of the coagulation rate constant for sulfuric acid
70  ! particles as a function of particle size using TDMA, Aerosol Science, 32, 321-339, 2001.
71  !--ok_vdw is 0 for no vdW forces, 1 for E(0), 2 for E(infinity)
72  INTEGER, PARAMETER                            :: ok_vdw = 0
73  REAL, PARAMETER                               :: avdW1 = 0.0757
74  REAL, PARAMETER                               :: avdW3 = 0.0015
75  REAL, PARAMETER                               :: bvdW0 = 0.0151
76  REAL, PARAMETER                               :: bvdW1 = -0.186
77  REAL, PARAMETER                               :: bvdW3 = -0.0163
78  REAL, PARAMETER                               :: AvdW = 6.4e-20 !Hamaker constant in J = 1e7 erg
79  REAL                                          :: AvdWi
80  REAL                                          :: xvdW
81  REAL                                          :: EvdW
82
83  DO i=1, nbtr_bin
84   radius(i)=mdw(i)/2.
85   V(i)= radius(i)**3.  !neglecting factor 4*RPI/3
86  ENDDO
87
88  DO j=1, nbtr_bin
89  DO i=1, nbtr_bin
90   Vij(i,j)= V(i)+V(j)
91  ENDDO
92  ENDDO
93
94!--pre-compute the f(i,j,k) from Jacobson equation 13
95  ff=0.0
96  DO k=1, nbtr_bin
97  DO j=1, nbtr_bin
98  DO i=1, nbtr_bin
99    IF (k.EQ.1) THEN
100      ff(i,j,k)= 0.0
101    ELSEIF (k.GT.1.AND.V(k-1).LT.Vij(i,j).AND.Vij(i,j).LT.V(k)) THEN
102      ff(i,j,k)= 1.-ff(i,j,k-1)
103    ELSEIF (k.EQ.nbtr_bin) THEN
104      IF (Vij(i,j).GE.v(k)) THEN
105        ff(i,j,k)= 1.
106      ELSE
107        ff(i,j,k)= 0.0
108      ENDIF
109    ELSEIF (k.LE.(nbtr_bin-1).AND.V(k).LE.Vij(i,j).AND.Vij(i,j).LT.V(k+1)) THEN
110      ff(i,j,k)= V(k)/Vij(i,j)*(V(k+1)-Vij(i,j))/(V(k+1)-V(k))
111    ENDIF
112  ENDDO
113  ENDDO
114  ENDDO
115
116  DO ilon=1, klon
117  DO ilev=1, klev
118  !only in the stratosphere
119  IF (is_strato(ilon,ilev)) THEN
120  !compute actual wet particle radius & volume for every grid box
121  DO i=1, nbtr_bin
122   radius(i)=f_r_wet(ilon,ilev)*mdw(i)/2.
123   V(i)= radius(i)**3.  !neglecting factor 4*RPI/3
124  ENDDO
125
126!--Calculations for the coagulation kernel---------------------------------------------------------
127
128  zrho=pplay(ilon,ilev)/t_seri(ilon,ilev)/RD
129
130!--initialize the tracer at time t and convert from [number/KgA] to [number/m3]
131  DO i=1, nbtr_bin
132  tr_t(ilon,ilev,i) = tr_seri(ilon,ilev,i+nbtr_sulgas) * zrho
133  ENDDO
134
135  ! mean free path of air (Pruppacher and Klett, 2010, p.417) [m]
136  mnfrpth=6.6E-8*(1.01325E+5/pplay(ilon,ilev))*(t_seri(ilon,ilev)/293.15)
137  ! mnfrpth=2.*eta/(zrho*thvelair)
138  ! mean free path of air (Prupp. Klett) in [10^-6 m]
139  ! ZLAIR = 0.066 *(1.01325E+5/PPLAY)*(T_SERI/293.15)*1.E-06
140
141  ! dynamic viscosity of air (Pruppacher and Klett, 2010, p.417) [kg/(m*s)]
142  IF (t_seri(ilon,ilev).GE.273.15) THEN
143    eta=(1.718+0.0049*(t_seri(ilon,ilev)-273.15))*1.E-5
144  ELSE
145    eta=(1.718+0.0049*(t_seri(ilon,ilev)-273.15)-1.2E-5*(t_seri(ilon,ilev)-273.15)**2)*1.E-5
146  ENDIF
147
148!--pre-compute the particle diffusion coefficient Di(i) from equation 18
149  Di=0.0
150  DO i=1, nbtr_bin
151   Kn(i)=mnfrpth/radius(i)
152   Di(i)=RKBOL*t_seri(ilon,ilev)/(6.*RPI*radius(i)*eta)*(1.+Kn(i)*(1.249+0.42*exp(-0.87/Kn(i))))
153  ENDDO
154
155!--pre-compute the thermal velocity of a particle thvelpar(i) from equation 20
156  thvelpar=0.0
157  DO i=1, nbtr_bin
158   m_par(i)=4./3.*RPI*radius(i)**3.*DENSO4(ilon,ilev)*1000.
159   thvelpar(i)=sqrt(8.*RKBOL*t_seri(ilon,ilev)/(RPI*m_par(i)))
160  ENDDO
161
162!--pre-compute the particle mean free path mfppar(i) from equation 22
163  mfppar=0.0
164  DO i=1, nbtr_bin
165   mfppar(i)=8.*Di(i)/(RPI*thvelpar(i))
166  ENDDO
167
168!--pre-compute the mean distance delta(i) from the center of a sphere reached by particles
169!--leaving the surface of the sphere and traveling a distance of particle mfppar(i) from equation 21
170  delta=0.0
171  DO i=1, nbtr_bin
172   delta(i)=((2.*radius(i)+mfppar(i))**3.-(4.*radius(i)**2.+mfppar(i)**2.)**1.5)/ &
173           & (6.*radius(i)*mfppar(i))-2.*radius(i)
174  ENDDO
175
176!--pre-compute the beta(i,j) from equation 17 in Jacobson
177  num=0.0
178  DO j=1, nbtr_bin
179  DO i=1, nbtr_bin
180!
181   num=4.*RPI*(radius(i)+radius(j))*(Di(i)+Di(j))
182   denom=(radius(i)+radius(j))/(radius(i)+radius(j)+sqrt(delta(i)**2.+delta(j)**2.))+ &
183        & 4.*(Di(i)+Di(j))/(sqrt(thvelpar(i)**2.+thvelpar(j)**2.)*(radius(i)+radius(j)))
184   beta(i,j)=num/denom
185!
186!--compute enhancement factor due to van der Waals forces
187   IF (ok_vdw .EQ. 0) THEN      !--no enhancement factor
188     Evdw=1.0
189   ELSEIF (ok_vdw .EQ. 1) THEN  !--E(0) case
190     AvdWi = AvdW/(RKBOL*t_seri(ilon,ilev))*(4.*radius(i)*radius(j))/(radius(i)+radius(j))**2.
191     xvdW = LOG(1.+AvdWi)
192     EvdW = 1. + avdW1*xvdW + avdW3*xvdW**3
193   ELSEIF (ok_vdw .EQ. 2) THEN  !--E(infinity) case
194     AvdWi = AvdW/(RKBOL*t_seri(ilon,ilev))*(4.*radius(i)*radius(j))/(radius(i)+radius(j))**2.
195     xvdW = LOG(1.+AvdWi)
196     EvdW = 1. + SQRT(AvdWi/3.)/(1.+bvdW0*SQRT(AvdWi)) + bvdW1*xvdW + bvdW3*xvdW**3.
197   ENDIF
198!
199   beta(i,j)=beta(i,j)*EvdW
200
201  ENDDO
202  ENDDO
203
204!--external loop for equation 14
205  DO k=1, nbtr_bin
206
207!--calculating denominator sum
208  denom=0.0
209  DO j=1, nbtr_bin
210  denom=denom+(1.-ff(k,j,k))*beta(k,j)*tr_t(ilon,ilev,j)
211  ENDDO
212
213  IF (k.EQ.1) THEN
214!--calculate new concentration of smallest bin
215    tr_tp1(ilon,ilev,k)=tr_t(ilon,ilev,k)/(1.+pdtcoag*denom)
216    ELSE
217!--calculating double sum terms in numerator of eq 14
218    num=0.0
219    DO j=1, k
220    numi=0.0
221    DO i=1, k-1
222    numi=numi+ff(i,j,k)*beta(i,j)*V(i)*tr_tp1(ilon,ilev,i)*tr_t(ilon,ilev,j)
223    ENDDO
224    num=num+numi
225    ENDDO
226
227!--calculate new concentration of other bins
228    tr_tp1(ilon,ilev,k)=(V(k)*tr_t(ilon,ilev,k)+pdtcoag*num)/(1.+pdtcoag*denom)/V(k)
229  ENDIF
230
231  ENDDO !--end of loop k
232
233!--convert tracer concentration back from [number/m3] to [number/KgA] and write into tr_seri
234  DO i=1, nbtr_bin
235   tr_seri(ilon,ilev,i+nbtr_sulgas) = tr_tp1(ilon,ilev,i) / zrho
236  ENDDO
237
238  ENDIF ! IF IN STRATOSPHERE
239  ENDDO !--end of loop klev
240  ENDDO !--end of loop klon
241
242END SUBROUTINE COAGULATE
Note: See TracBrowser for help on using the repository browser.