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

Last change on this file since 4768 was 4762, checked in by lguez, 11 months ago

Use YOMCST from LMDZ, not from RRTM or ECRad

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