source: LMDZ5/branches/testing/libf/phylmd/StratAer/coagulate.F90 @ 5215

Last change on this file since 5215 was 2720, checked in by Laurent Fairhead, 8 years ago

Merged trunk changes r2664:2719 into testing branch
File size: 7.6 KB
RevLine 
[2690]1SUBROUTINE COAGULATE(pdtcoag,mdw,tr_seri,t_seri,pplay,dens_aer,is_strato)
2  !     -----------------------------------------------------------------------
3  !
4  !     Author : Christoph Kleinschmitt (with Olivier Boucher)
5  !     ------
6  !
7  !     purpose
8  !     -------
9  !
10  !     interface
11  !     ---------
12  !      input
13  !       pdtphys        time step duration                 [sec]
14  !       tr_seri        tracer mixing ratios               [kg/kg]
15  !       mdw             # or mass median diameter          [m]
16  !
17  !     method
18  !     ------
19  !
20  !     -----------------------------------------------------------------------
21
22  USE dimphy, ONLY : klon,klev
23  USE aerophys
24  USE infotrac
25  USE phys_local_var_mod, ONLY: DENSO4, f_r_wet
26  USE YOMCST
27
28  IMPLICIT NONE
29
30  !--------------------------------------------------------
31
32  ! transfer variables when calling this routine
33  REAL,INTENT(IN)                               :: pdtcoag ! Time step in coagulation routine [s]
34  REAL,DIMENSION(nbtr_bin),INTENT(IN)           :: mdw     ! aerosol particle diameter in each bin [m]
35  REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT)  :: tr_seri ! Concentration Traceur [U/KgA]
36  REAL,DIMENSION(klon,klev),INTENT(IN)          :: t_seri  ! Temperature
37  REAL,DIMENSION(klon,klev),INTENT(IN)          :: pplay   ! pression pour le mileu de chaque couche (en Pa)
38  REAL,DIMENSION(klon,klev)                     :: dens_aer! density of aerosol [kg/m3 aerosol] with default H2SO4 mass
39  LOGICAL,DIMENSION(klon,klev),INTENT(IN)       :: is_strato
40
41  ! local variables in coagulation routine
42  INTEGER                                       :: i,j,k,nb,ilon,ilev
43  REAL, DIMENSION(nbtr_bin)                     :: radius ! aerosol particle radius in each bin [m]
44  REAL, DIMENSION(klon,klev,nbtr_bin)           :: tr_t ! Concentration Traceur at time t [U/KgA]
45  REAL, DIMENSION(klon,klev,nbtr_bin)           :: tr_tp1 ! Concentration Traceur at time t+1 [U/KgA]
46  REAL, DIMENSION(nbtr_bin,nbtr_bin,nbtr_bin)   :: ff   ! Volume fraction of intermediate particles
47  REAL, DIMENSION(nbtr_bin)                     :: V    ! Volume of bins
48  REAL, DIMENSION(nbtr_bin,nbtr_bin)            :: Vij  ! Volume sum of i and j
49  REAL                                          :: eta  ! Dynamic viscosity of air
50  REAL, PARAMETER                               :: mair=4.8097E-26 ! Average mass of an air molecule [kg]
51  REAL                                          :: zrho ! Density of air
52  REAL                                          :: mnfrpth ! Mean free path of air
53  REAL, DIMENSION(nbtr_bin)                     :: Kn   ! Knudsen number of particle i
54  REAL, DIMENSION(nbtr_bin)                     :: Di   ! Particle diffusion coefficient
55  REAL, DIMENSION(nbtr_bin)                     :: m_par   ! Mass of particle i
56  REAL, DIMENSION(nbtr_bin)                     :: thvelpar! Thermal velocity of particle i
57  REAL, DIMENSION(nbtr_bin)                     :: mfppar  ! Mean free path of particle i
58  REAL, DIMENSION(nbtr_bin)                     :: delta! delta of particle i (from equation 21)
59  REAL, DIMENSION(nbtr_bin,nbtr_bin)            :: beta ! Coagulation kernel from Brownian diffusion
60  REAL                                          :: beta_const ! Constant coagulation kernel (for comparison)
61  REAL                                          :: num
62  REAL                                          :: numi
63  REAL                                          :: denom
64
65  DO i=1, nbtr_bin
66   radius(i)=mdw(i)/2.
67   V(i)= radius(i)**3.  !neglecting factor 4*RPI/3
68  ENDDO
69
70  DO j=1, nbtr_bin
71  DO i=1, nbtr_bin
72   Vij(i,j)= V(i)+V(j)
73  ENDDO
74  ENDDO
75
76!--pre-compute the f(i,j,k) from Jacobson equation 13
77  ff=0.0
78  DO k=1, nbtr_bin
79  DO j=1, nbtr_bin
80  DO i=1, nbtr_bin
81    IF (k.EQ.1) THEN
82      ff(i,j,k)= 0.0
83    ELSEIF (k.GT.1.AND.V(k-1).LT.Vij(i,j).AND.Vij(i,j).LT.V(k)) THEN
84      ff(i,j,k)= 1.-ff(i,j,k-1)
85    ELSEIF (k.EQ.nbtr_bin) THEN
86      IF (Vij(i,j).GE.v(k)) THEN
87        ff(i,j,k)= 1.
88      ELSE
89        ff(i,j,k)= 0.0
90      ENDIF
91    ELSEIF (k.LE.(nbtr_bin-1).AND.V(k).LE.Vij(i,j).AND.Vij(i,j).LT.V(k+1)) THEN
92      ff(i,j,k)= V(k)/Vij(i,j)*(V(k+1)-Vij(i,j))/(V(k+1)-V(k))
93    ENDIF
94  ENDDO
95  ENDDO
96  ENDDO
97
98  DO ilon=1, klon
99  DO ilev=1, klev
100  !only in the stratosphere
101  IF (is_strato(ilon,ilev)) THEN
102  !compute actual wet particle radius & volume for every grid box
103  DO i=1, nbtr_bin
104   radius(i)=f_r_wet(ilon,ilev)*mdw(i)/2.
105   V(i)= radius(i)**3.  !neglecting factor 4*RPI/3
106  ENDDO
107
108!--Calculations for the coagulation kernel---------------------------------------------------------
109
110  zrho=pplay(ilon,ilev)/t_seri(ilon,ilev)/RD
111
112!--initialize the tracer at time t and convert from [number/KgA] to [number/m3]
113  DO i=1, nbtr_bin
114  tr_t(ilon,ilev,i) = tr_seri(ilon,ilev,i+nbtr_sulgas) * zrho
115  ENDDO
116
117  ! mean free path of air (Pruppacher and Klett, 2010, p.417) [m]
118  mnfrpth=6.6E-8*(1.01325E+5/pplay(ilon,ilev))*(t_seri(ilon,ilev)/293.15)
119  ! mnfrpth=2.*eta/(zrho*thvelair)
120  ! mean free path of air (Prupp. Klett) in [10^-6 m]
121  ! ZLAIR = 0.066 *(1.01325E+5/PPLAY)*(T_SERI/293.15)*1.E-06
122
123  ! dynamic viscosity of air (Pruppacher and Klett, 2010, p.417) [kg/(m*s)]
124  IF (t_seri(ilon,ilev).GE.273.15) THEN
125    eta=(1.718+0.0049*(t_seri(ilon,ilev)-273.15))*1.E-5
126  ELSE
127    eta=(1.718+0.0049*(t_seri(ilon,ilev)-273.15)-1.2E-5*(t_seri(ilon,ilev)-273.15)**2)*1.E-5
128  ENDIF
129
130!--pre-compute the particle diffusion coefficient Di(i) from equation 18
131  Di=0.0
132  DO i=1, nbtr_bin
133   Kn(i)=mnfrpth/radius(i)
134   Di(i)=RKBOL*t_seri(ilon,ilev)/(6.*RPI*radius(i)*eta)*(1.+Kn(i)*(1.249+0.42*exp(-0.87/Kn(i))))
135  ENDDO
136
137!--pre-compute the thermal velocity of a particle thvelpar(i) from equation 20
138  thvelpar=0.0
139  DO i=1, nbtr_bin
140   m_par(i)=4./3.*RPI*radius(i)**3.*DENSO4(ilon,ilev)*1000.
141   thvelpar(i)=sqrt(8.*RKBOL*t_seri(ilon,ilev)/(RPI*m_par(i)))
142  ENDDO
143
144!--pre-compute the particle mean free path mfppar(i) from equation 22
145  mfppar=0.0
146  DO i=1, nbtr_bin
147   mfppar(i)=8.*Di(i)/(RPI*thvelpar(i))
148  ENDDO
149
150!--pre-compute the mean distance delta(i) from the center of a sphere reached by particles
151!--leaving the surface of the sphere and traveling a distance of particle mfppar(i) from equation 21
152  delta=0.0
153  DO i=1, nbtr_bin
154   delta(i)=((2.*radius(i)+mfppar(i))**3.-(4.*radius(i)**2.+mfppar(i)**2.)**1.5)/ &
155           & (6.*radius(i)*mfppar(i))-2.*radius(i)
156  ENDDO
157
158!--pre-compute the beta(i,j) from equation 17
159  num=0.0
160  DO j=1, nbtr_bin
161  DO i=1, nbtr_bin
162   num=4.*RPI*(radius(i)+radius(j))*(Di(i)+Di(j))
163   denom=(radius(i)+radius(j))/(radius(i)+radius(j)+sqrt(delta(i)**2.+delta(j)**2.))+ &
164        & 4.*(Di(i)+Di(j))/(sqrt(thvelpar(i)**2.+thvelpar(j)**2.)*(radius(i)+radius(j)))
165   beta(i,j)=num/denom
166  ENDDO
167  ENDDO
168
169!--external loop for equation 14
170  DO k=1, nbtr_bin
171
172!--calculating denominator sum
173  denom=0.0
174  DO j=1, nbtr_bin
175  denom=denom+(1.-ff(k,j,k))*beta(k,j)*tr_t(ilon,ilev,j)
176  ENDDO
177
178  IF (k.EQ.1) THEN
179!--calculate new concentration of smallest bin
180    tr_tp1(ilon,ilev,k)=tr_t(ilon,ilev,k)/(1.+pdtcoag*denom)
181    ELSE
182!--calculating double sum terms in numerator of eq 14
183    num=0.0
184    DO j=1, k
185    numi=0.0
186    DO i=1, k-1
187    numi=numi+ff(i,j,k)*beta(i,j)*V(i)*tr_tp1(ilon,ilev,i)*tr_t(ilon,ilev,j)
188    ENDDO
189    num=num+numi
190    ENDDO
191
192!--calculate new concentration of other bins
193    tr_tp1(ilon,ilev,k)=(V(k)*tr_t(ilon,ilev,k)+pdtcoag*num)/(1.+pdtcoag*denom)/V(k)
194  ENDIF
195
196  ENDDO !--end of loop k
197
198!--convert tracer concentration back from [number/m3] to [number/KgA] and write into tr_seri
199  DO i=1, nbtr_bin
200   tr_seri(ilon,ilev,i+nbtr_sulgas) = tr_tp1(ilon,ilev,i) / zrho
201  ENDDO
202
203  ENDIF ! IF IN STRATOSPHERE
204  ENDDO !--end of loop klev
205  ENDDO !--end of loop klon
206
207END SUBROUTINE COAGULATE
Note: See TracBrowser for help on using the repository browser.