source: LMDZ6/branches/Amaury_dev/libf/phylmd/StratAer/coagulate.F90 @ 5101

Last change on this file since 5101 was 5101, checked in by abarral, 2 months ago

Handle DEBUG_IO in lmdz_cppkeys_wrapper.F90
Transform some files .F -> .[fF]90
[ne compile pas à cause de writefield_u non défini - en attente de réponse Laurent]
  • Property svn:keywords set to Id
File size: 11.9 KB
RevLine 
[5099]1
[3526]2! $Id: coagulate.F90 5101 2024-07-23 06:22:55Z abarral $
[5099]3
[2690]4SUBROUTINE COAGULATE(pdtcoag,mdw,tr_seri,t_seri,pplay,dens_aer,is_strato)
5  !     -----------------------------------------------------------------------
[5099]6
[2690]7  !     Author : Christoph Kleinschmitt (with Olivier Boucher)
8  !     ------
[5099]9
[2690]10  !     purpose
11  !     -------
[5099]12
[2690]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]
[5099]19
[2690]20  !     method
21  !     ------
[5099]22
[2690]23  !     -----------------------------------------------------------------------
24
[5101]25  USE dimphy, ONLY: klon,klev
[2690]26  USE aerophys
[3677]27  USE infotrac_phy
[4950]28  USE phys_local_var_mod, ONLY: DENSO4, DENSO4B, f_r_wet, f_r_wetB
29  USE strataer_local_var_mod, ONLY: flag_new_strat_compo
[5098]30  USE lmdz_yomcst
[4950]31 
[2690]32  IMPLICIT NONE
33
34  !--------------------------------------------------------
35
36  ! transfer variables when calling this routine
37  REAL,INTENT(IN)                               :: pdtcoag ! Time step in coagulation routine [s]
38  REAL,DIMENSION(nbtr_bin),INTENT(IN)           :: mdw     ! aerosol particle diameter in each bin [m]
39  REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT)  :: tr_seri ! Concentration Traceur [U/KgA]
40  REAL,DIMENSION(klon,klev),INTENT(IN)          :: t_seri  ! Temperature
41  REAL,DIMENSION(klon,klev),INTENT(IN)          :: pplay   ! pression pour le mileu de chaque couche (en Pa)
42  REAL,DIMENSION(klon,klev)                     :: dens_aer! density of aerosol [kg/m3 aerosol] with default H2SO4 mass
43  LOGICAL,DIMENSION(klon,klev),INTENT(IN)       :: is_strato
44
45  ! local variables in coagulation routine
46  INTEGER                                       :: i,j,k,nb,ilon,ilev
[4950]47  REAL, DIMENSION(nbtr_bin)                     :: radiusdry ! dry aerosol particle radius in each bin [m]
48  REAL, DIMENSION(nbtr_bin)                     :: radiuswet ! wet aerosol particle radius in each bin [m]
[2690]49  REAL, DIMENSION(klon,klev,nbtr_bin)           :: tr_t ! Concentration Traceur at time t [U/KgA]
50  REAL, DIMENSION(klon,klev,nbtr_bin)           :: tr_tp1 ! Concentration Traceur at time t+1 [U/KgA]
51  REAL, DIMENSION(nbtr_bin,nbtr_bin,nbtr_bin)   :: ff   ! Volume fraction of intermediate particles
[4950]52  REAL, DIMENSION(nbtr_bin)                     :: Vdry ! Volume dry of bins
53  REAL, DIMENSION(nbtr_bin)                     :: Vwet ! Volume wet of bins
[2690]54  REAL, DIMENSION(nbtr_bin,nbtr_bin)            :: Vij  ! Volume sum of i and j
55  REAL                                          :: eta  ! Dynamic viscosity of air
56  REAL, PARAMETER                               :: mair=4.8097E-26 ! Average mass of an air molecule [kg]
57  REAL                                          :: zrho ! Density of air
58  REAL                                          :: mnfrpth ! Mean free path of air
59  REAL, DIMENSION(nbtr_bin)                     :: Kn   ! Knudsen number of particle i
60  REAL, DIMENSION(nbtr_bin)                     :: Di   ! Particle diffusion coefficient
61  REAL, DIMENSION(nbtr_bin)                     :: m_par   ! Mass of particle i
62  REAL, DIMENSION(nbtr_bin)                     :: thvelpar! Thermal velocity of particle i
63  REAL, DIMENSION(nbtr_bin)                     :: mfppar  ! Mean free path of particle i
64  REAL, DIMENSION(nbtr_bin)                     :: delta! delta of particle i (from equation 21)
65  REAL, DIMENSION(nbtr_bin,nbtr_bin)            :: beta ! Coagulation kernel from Brownian diffusion
66  REAL                                          :: beta_const ! Constant coagulation kernel (for comparison)
67  REAL                                          :: num
68  REAL                                          :: numi
69  REAL                                          :: denom
70
[2949]71  ! Additional variables for coagulation enhancement factor due to van der Waals forces
72  ! Taken from Chan and Mozurkewich, Measurement of the coagulation rate constant for sulfuric acid
73  ! particles as a function of particle size using TDMA, Aerosol Science, 32, 321-339, 2001.
74  !--ok_vdw is 0 for no vdW forces, 1 for E(0), 2 for E(infinity)
75  INTEGER, PARAMETER                            :: ok_vdw = 0
76  REAL, PARAMETER                               :: avdW1 = 0.0757
77  REAL, PARAMETER                               :: avdW3 = 0.0015
78  REAL, PARAMETER                               :: bvdW0 = 0.0151
79  REAL, PARAMETER                               :: bvdW1 = -0.186
80  REAL, PARAMETER                               :: bvdW3 = -0.0163
81  REAL, PARAMETER                               :: AvdW = 6.4e-20 !Hamaker constant in J = 1e7 erg
82  REAL                                          :: AvdWi
83  REAL                                          :: xvdW
84  REAL                                          :: EvdW
85
[4950]86! ff(i,j,k): Volume fraction of Vi,j that is partitioned to each model bin k
87! just need to be calculated in model initialization because mdw(:) size is fixed
88! no need to recalculate radius, Vdry, Vij, and ff every timestep because it is for 
89! dry aerosols
[2690]90  DO i=1, nbtr_bin
[4950]91     radiusdry(i)=mdw(i)/2.
92     Vdry(i)=radiusdry(i)**3.  !neglecting factor 4*RPI/3
93     Vwet(i)=0.0
[2690]94  ENDDO
95
96  DO j=1, nbtr_bin
[4950]97     DO i=1, nbtr_bin
98        Vij(i,j)= Vdry(i)+Vdry(j)
99     ENDDO
[2690]100  ENDDO
[4950]101 
[2690]102!--pre-compute the f(i,j,k) from Jacobson equation 13
103  ff=0.0
104  DO k=1, nbtr_bin
105  DO j=1, nbtr_bin
106  DO i=1, nbtr_bin
[5082]107    IF (k==1) THEN
[2690]108      ff(i,j,k)= 0.0
[5082]109    ELSEIF (k>1.AND.Vdry(k-1)<Vij(i,j).AND.Vij(i,j)<Vdry(k)) THEN
[2690]110      ff(i,j,k)= 1.-ff(i,j,k-1)
[5082]111    ELSEIF (k==nbtr_bin) THEN
112      IF (Vij(i,j)>=Vdry(k)) THEN
[2690]113        ff(i,j,k)= 1.
114      ELSE
115        ff(i,j,k)= 0.0
116      ENDIF
[5082]117    ELSEIF (k<=(nbtr_bin-1).AND.Vdry(k)<=Vij(i,j).AND.Vij(i,j)<Vdry(k+1)) THEN
[4950]118      ff(i,j,k)= Vdry(k)/Vij(i,j)*(Vdry(k+1)-Vij(i,j))/(Vdry(k+1)-Vdry(k))
[2690]119    ENDIF
120  ENDDO
121  ENDDO
122  ENDDO
[4950]123! End of just need to be calculated at initialization because mdw(:) size is fixed
124 
[2690]125  DO ilon=1, klon
126  DO ilev=1, klev
127  !only in the stratosphere
128  IF (is_strato(ilon,ilev)) THEN
129  !compute actual wet particle radius & volume for every grid box
[4950]130  IF(flag_new_strat_compo) THEN
131     DO i=1, nbtr_bin
132        radiuswet(i)=f_r_wetB(ilon,ilev,i)*mdw(i)/2.
133        Vwet(i)= radiuswet(i)**3.  !neglecting factor 4*RPI/3
134!!      Vwet(i)= Vdry(i)*(f_r_wetB(ilon,ilev,i)**3)
135     ENDDO
136  ELSE
137     DO i=1, nbtr_bin
138        radiuswet(i)=f_r_wet(ilon,ilev)*mdw(i)/2.
139        Vwet(i)= radiuswet(i)**3.  !neglecting factor 4*RPI/3
140!!      Vwet(i)= Vdry(i)*(f_r_wet(ilon,ilev)**3)
141     ENDDO
142  ENDIF
143 
[2690]144!--Calculations for the coagulation kernel---------------------------------------------------------
145
146  zrho=pplay(ilon,ilev)/t_seri(ilon,ilev)/RD
147
148!--initialize the tracer at time t and convert from [number/KgA] to [number/m3]
149  DO i=1, nbtr_bin
150  tr_t(ilon,ilev,i) = tr_seri(ilon,ilev,i+nbtr_sulgas) * zrho
151  ENDDO
152
153  ! mean free path of air (Pruppacher and Klett, 2010, p.417) [m]
154  mnfrpth=6.6E-8*(1.01325E+5/pplay(ilon,ilev))*(t_seri(ilon,ilev)/293.15)
155  ! mnfrpth=2.*eta/(zrho*thvelair)
156  ! mean free path of air (Prupp. Klett) in [10^-6 m]
157  ! ZLAIR = 0.066 *(1.01325E+5/PPLAY)*(T_SERI/293.15)*1.E-06
158
159  ! dynamic viscosity of air (Pruppacher and Klett, 2010, p.417) [kg/(m*s)]
[5082]160  IF (t_seri(ilon,ilev)>=273.15) THEN
[2690]161    eta=(1.718+0.0049*(t_seri(ilon,ilev)-273.15))*1.E-5
162  ELSE
163    eta=(1.718+0.0049*(t_seri(ilon,ilev)-273.15)-1.2E-5*(t_seri(ilon,ilev)-273.15)**2)*1.E-5
164  ENDIF
165
166!--pre-compute the particle diffusion coefficient Di(i) from equation 18
167  Di=0.0
168  DO i=1, nbtr_bin
[4950]169      Kn(i)=mnfrpth/radiuswet(i)
170      Di(i)=RKBOL*t_seri(ilon,ilev)/(6.*RPI*radiuswet(i)*eta)*(1.+Kn(i)*(1.249+0.42*exp(-0.87/Kn(i))))
[2690]171  ENDDO
172
173!--pre-compute the thermal velocity of a particle thvelpar(i) from equation 20
174  thvelpar=0.0
[4950]175  IF(flag_new_strat_compo) THEN
176     DO i=1, nbtr_bin
177        m_par(i)=4./3.*RPI*radiuswet(i)**3.*DENSO4B(ilon,ilev,i)*1000.
178        thvelpar(i)=sqrt(8.*RKBOL*t_seri(ilon,ilev)/(RPI*m_par(i)))
179     ENDDO
180  ELSE
181     DO i=1, nbtr_bin
182        m_par(i)=4./3.*RPI*radiuswet(i)**3.*DENSO4(ilon,ilev)*1000.
183        thvelpar(i)=sqrt(8.*RKBOL*t_seri(ilon,ilev)/(RPI*m_par(i)))
184     ENDDO
185  ENDIF
[2690]186
187!--pre-compute the particle mean free path mfppar(i) from equation 22
188  mfppar=0.0
189  DO i=1, nbtr_bin
190   mfppar(i)=8.*Di(i)/(RPI*thvelpar(i))
191  ENDDO
192
193!--pre-compute the mean distance delta(i) from the center of a sphere reached by particles
194!--leaving the surface of the sphere and traveling a distance of particle mfppar(i) from equation 21
195  delta=0.0
196  DO i=1, nbtr_bin
[4950]197      delta(i)=((2.*radiuswet(i)+mfppar(i))**3.-(4.*radiuswet(i)**2.+mfppar(i)**2.)**1.5)/ &
[5087]198   (6.*radiuswet(i)*mfppar(i))-2.*radiuswet(i)
[2690]199  ENDDO
200
[4950]201!   beta(i,j): coagulation kernel (rate coefficient) of 2 colliding particles i,j
[2949]202!--pre-compute the beta(i,j) from equation 17 in Jacobson
[2690]203  num=0.0
204  DO j=1, nbtr_bin
205  DO i=1, nbtr_bin
[5099]206
[4950]207     num=4.*RPI*(radiuswet(i)+radiuswet(j))*(Di(i)+Di(j))
208     denom=(radiuswet(i)+radiuswet(j))/(radiuswet(i)+radiuswet(j)+sqrt(delta(i)**2.+delta(j)**2.))+ &
[5087]209   4.*(Di(i)+Di(j))/(sqrt(thvelpar(i)**2.+thvelpar(j)**2.)*(radiuswet(i)+radiuswet(j)))
[4950]210     beta(i,j)=num/denom
[5099]211
[2949]212!--compute enhancement factor due to van der Waals forces
[5082]213   IF (ok_vdw == 0) THEN      !--no enhancement factor
[4950]214      Evdw=1.0
[5082]215   ELSEIF (ok_vdw == 1) THEN  !--E(0) case
[4950]216      AvdWi = AvdW/(RKBOL*t_seri(ilon,ilev))*(4.*radiuswet(i)*radiuswet(j))/(radiuswet(i)+radiuswet(j))**2.
217      xvdW = LOG(1.+AvdWi)
218      EvdW = 1. + avdW1*xvdW + avdW3*xvdW**3
[5082]219   ELSEIF (ok_vdw == 2) THEN  !--E(infinity) case
[4950]220      AvdWi = AvdW/(RKBOL*t_seri(ilon,ilev))*(4.*radiuswet(i)*radiuswet(j))/(radiuswet(i)+radiuswet(j))**2.
221      xvdW = LOG(1.+AvdWi)
222      EvdW = 1. + SQRT(AvdWi/3.)/(1.+bvdW0*SQRT(AvdWi)) + bvdW1*xvdW + bvdW3*xvdW**3.
[2949]223   ENDIF
[5099]224
[2949]225   beta(i,j)=beta(i,j)*EvdW
226
[2690]227  ENDDO
228  ENDDO
229
230!--external loop for equation 14
231  DO k=1, nbtr_bin
232
233!--calculating denominator sum
234  denom=0.0
235  DO j=1, nbtr_bin
[4950]236     !    fraction of coagulation of k and j that is not giving k
237     denom=denom+(1.-ff(k,j,k))*beta(k,j)*tr_t(ilon,ilev,j)
[2690]238  ENDDO
239
[5082]240  IF (k==1) THEN
[2690]241!--calculate new concentration of smallest bin
242    tr_tp1(ilon,ilev,k)=tr_t(ilon,ilev,k)/(1.+pdtcoag*denom)
243    ELSE
244!--calculating double sum terms in numerator of eq 14
245    num=0.0
246    DO j=1, k
[4950]247       numi=0.0
248       DO i=1, k-1
[5099]249
[4950]250!           see Jacobson: " In order to conserve volume and volume concentration (which
251!           coagulation physically does) while giving up some accuracy in number concentration"
[5099]252
[4950]253!           Coagulation of i and j giving k
254!           with V(i) and then V(j) because it considers i,j and j,i with the double loop
[5099]255
[4950]256!           BUT WHY WET VOLUME V(i) in old STRATAER? tracers are already dry aerosols and coagulation
257!           kernel beta(i,j) accounts for wet aerosols -> reply below
[5099]258
[4950]259!             numi=numi+ff(i,j,k)*beta(i,j)*V(i)*tr_tp1(ilon,ilev,i)*tr_t(ilon,ilev,j)
260            numi=numi+ff(i,j,k)*beta(i,j)*Vdry(i)*tr_tp1(ilon,ilev,i)*tr_t(ilon,ilev,j)
261       ENDDO
262       num=num+numi
[2690]263    ENDDO
264
265!--calculate new concentration of other bins
[4950]266!      tr_tp1(ilon,ilev,k)=(V(k)*tr_t(ilon,ilev,k)+pdtcoag*num)/( (1.+pdtcoag*denom)*V(k) )
267    tr_tp1(ilon,ilev,k)=(Vdry(k)*tr_t(ilon,ilev,k)+pdtcoag*num)/( (1.+pdtcoag*denom)*Vdry(k) )
[5099]268
[4950]269!       In constant composition (no dependency on aerosol size because no kelvin effect)
270!       V(l)= (f_r_wet(ilon,ilev)**3)*((mdw(l)/2.)**3) = (f_r_wet(ilon,ilev)**3)*Vdry(i)
271!       so numi and num are proportional (f_r_wet(ilon,ilev)**3)
272!       and so
273!        tr_tp1(ilon,ilev,k)=(V(k)*tr_t(ilon,ilev,k)+pdtcoag*num)/( (1.+pdtcoag*denom)*V(k) )
274!                     =(Vdry(k)*tr_t(ilon,ilev,k)+pdtcoag*num_dry)/( (1.+pdtcoag*denom)*Vdry(k) )
275!          with num_dry=...beta(i,j)*Vdry(i)*....
276!       so in old STRATAER (.not.flag_new_strat_compo), it was correct
[2690]277  ENDIF
278
279  ENDDO !--end of loop k
280
281!--convert tracer concentration back from [number/m3] to [number/KgA] and write into tr_seri
282  DO i=1, nbtr_bin
[4950]283     tr_seri(ilon,ilev,i+nbtr_sulgas) = tr_tp1(ilon,ilev,i) / zrho
[2690]284  ENDDO
285
286  ENDIF ! IF IN STRATOSPHERE
287  ENDDO !--end of loop klev
288  ENDDO !--end of loop klon
[4950]289! *********************************************
[2690]290
291END SUBROUTINE COAGULATE
Note: See TracBrowser for help on using the repository browser.