source: LMDZ6/trunk/libf/phylmd/StratAer/cond_evap_tstep_mod.F90

Last change on this file was 4950, checked in by lebasn, 5 weeks ago

StratAer?: New model version (microphysic, composition routine, code cleaning, new params...)
  • Property svn:keywords set to Id
File size: 12.3 KB
Line 
1!
2! $Id: cond_evap_tstep_mod.F90 4950 2024-05-22 13:16:36Z abarral $
3!
4MODULE cond_evap_tstep_mod
5
6! based on UPMC aerosol model by Slimane Bekki
7! adapted for stratospheric sulfate aerosol in LMDZ by Christoph Kleinschmitt
8
9CONTAINS
10
11      SUBROUTINE condens_evapor_rate_kelvin(R2SO4G,t_seri,pplay,R2SO4, &
12          & DENSO4,f_r_wet,R2SO4ik,DENSO4ik,f_r_wetik,FL,ASO4,DNDR)
13!
14!     INPUT:
15!     R2SO4G: number density of gaseous H2SO4 [molecules/cm3]
16!     t_seri: temperature (K)
17!     pplay: pressure (Pa)
18!     R2SO4: aerosol H2SO4 weight fraction (percent) - flat surface (does not depend on aerosol size)
19!     DENSO4: aerosol density (gr/cm3)
20!     f_r_wet: factor for converting dry to wet radius
21!        assuming 'flat surface' composition (does not depend on aerosol size)
22!     variables that depends on aerosol size because of Kelvin effect
23!     R2SO4Gik: number density of gaseous H2SO4 [molecules/cm3] - depends on aerosol size
24!     DENSO4ik: aerosol density (gr/cm3) - depends on aerosol size
25!     f_r_wetik: factor for converting dry to wet radius - depends on aerosol size
26!     RRSI: radius [cm]
27
28      USE aerophys
29      USE infotrac_phy
30      USE YOMCST, ONLY : RPI
31      USE sulfate_aer_mod, ONLY : wph2so4, surftension, solh2so4, rpmvh2so4
32      USE strataer_local_var_mod, ONLY : ALPH2SO4, RRSI
33     
34      IMPLICIT NONE
35     
36      REAL, PARAMETER :: third=1./3.
37     
38      ! input variables
39      REAL            :: R2SO4G !H2SO4 number density [molecules/cm3]
40      REAL            :: t_seri
41      REAL            :: pplay
42      REAL            :: R2SO4
43      REAL            :: DENSO4
44      REAL            :: f_r_wet
45      REAL            :: R2SO4ik(nbtr_bin), DENSO4ik(nbtr_bin), f_r_wetik(nbtr_bin)
46     
47      ! output variables
48      REAL            :: FL(nbtr_bin)
49      REAL            :: ASO4(nbtr_bin)
50      REAL            :: DNDR(nbtr_bin)
51     
52      ! local variables
53      INTEGER            :: IK
54      REAL            :: ALPHA,CST
55      REAL            :: WH2(nbtr_bin)
56      REAL            :: RP,VTK,AA,FL1,RKNUD
57      REAL            :: DND
58      REAL            :: ATOT,AH2O
59      REAL            :: RRSI_wet(nbtr_bin)
60      REAL            :: FPATH, WPP, XA, FKELVIN
61      REAL            :: surtens, mvh2so4, temp
62     
63! ///    MOLEC CONDENSATION GROWTH (DUE TO CHANGES IN H2SO4 AND SO H2O)
64!    ------------------------------------------------------------------
65!                                  EXCEPT CN
66!       RK:H2SO4 WEIGHT PERCENT DOESN'T CHANGE
67!     BE CAREFUL,H2SO4 WEIGHT PERCENTAGE
68
69!                   MOLECULAR ACCOMODATION OF H2SO4
70!     H2SO4 accommodation  coefficient [condensation/evaporation]
71      ALPHA = ALPH2SO4
72!      FPLAIR=(2.281238E-5)*TAIR/PAIR
73!     1.E2 (m to cm),
74      CST=1.E2*2.281238E-5
75!     same expression as in coagulate
76!     in coagulate: mean free path of air (Pruppacher and Klett, 2010, p.417) [m]
77!     mnfrpth=6.6E-8*(1.01325E+5/pplay(ilon,ilev))*(t_seri(ilon,ilev)/293.15)
78!     mnfrpth=2.28E-5*t_seri/pplay
79     
80      temp = min( max(t_seri, 190.), 300.) ! 190K <= temp <= 300K
81     
82      RRSI_wet(:)=RRSI(:)*f_r_wetik(:)
83
84!     Pruppa and Klett
85      FPATH=CST*t_seri/pplay
86   
87!     H2SO4 mass fraction in aerosol
88      WH2(:)=R2SO4ik(:)*1.0E-2
89
90!                               ACTIVITY COEFFICIENT(SEE GIAUQUE,1951)
91!                               AYERS ET AL (1980)
92!                                  (MU-MU0)
93!      RP=-10156.0/t_seri +16.259-(ACTSO4*4.184)/(8.31441*t_seri)
94!                                  DROPLET H2SO4 PRESSURE IN DYN.CM-2
95!      RP=EXP(RP)*1.01325E6/0.086
96!!      RP=EXP(RP)*1.01325E6
97!                                  H2SO4 NUMBER DENSITY NEAR DROPLET
98
99!      DND=RP*6.02E23/(8.31E7*t_seri)
100
101!                                 KELVIN EFFECT FACTOR
102!CK 20160613: bug fix, removed factor 250 (from original code by S. Bekki)
103!!      AA =2.0*MH2O*72.0/(DENSO4*BOLZ*t_seri*250.0)
104!      AA =2.0*MH2O*72.0/(DENSO4*BOLZ*t_seri)
105
106!                                  MEAN KINETIC VELOCITY
107!     DYN*CM*K/(K*GR)=(CM/SEC2)*CM
108!                                  IN CM/SEC
109      VTK=SQRT(8.0*BOLZ*t_seri/(RPI*MH2SO4))
110!                                 KELVIN EFFECT FACTOR
111
112!     Loop on bin radius (RRSI in cm)
113      DO IK=1,nbtr_bin
114
115      IF(R2SO4ik(IK) > 0.0) THEN
116
117!       h2so4 mass fraction (0<wpp<1)
118        wpp=R2SO4ik(IK)*1.e-2   
119        xa=18.*wpp/(18.*wpp+98.*(1.-wpp))
120!       equilibrium h2so4 number density over H2SO4/H2O solution (molec/cm3)
121        DND=solh2so4(t_seri,xa)
122!          KELVIN EFFECT: 
123!       surface tension (mN/m=1.e-3.kg/s2) = f(T,h2so4 mole fraction)
124        surtens=surftension(temp,xa)
125!       partial molar volume of h2so4 (cm3.mol-1 =1.e-6.m3.mol-1)
126        mvh2so4= rpmvh2so4(temp,R2SO4ik(IK))
127!       Kelvin factor (MKS)
128        fkelvin=exp( 2.*1.e-3*surtens*1.e-6*mvh2so4/ (1.e-2*RRSI_wet(IK)*rgas*temp) )
129!                             
130        DNDR(IK) =DND*fkelvin
131
132        FL1=RPI*ALPHA*VTK*(R2SO4G-DNDR(IK))
133
134!       TURCO(1979) FOR HNO3:ALH2SO4 CONDENSATION= ALH2SO4 EVAPORATION
135!       RPI*R2*VTK IS EQUIVALENT TO DIFFUSION COEFFICIENT
136!       EXTENSION OF THE RELATION FOR DIFFUSION KINETICS
137!       KNUDSEN NUMBER FPATH/RRSI
138!       NEW VERSION (SEE NOTES)
139        RKNUD=FPATH/RRSI_wet(IK)
140!       SENFELD
141        FL(IK)=FL1*RRSI_wet(IK)**2*( 1.0 +RKNUD ) &
142     &     /( 1.0 +ALPHA/(2.0*RKNUD) +RKNUD )
143!       TURCO
144!        RL= (4.0/3.0 +0.71/RKNUD)/(1.0+1.0/RKNUD)
145!     *         +4.0*(1.0-ALPHA)/(3.0*ALPHA)
146!        FL=FL1*RRSI(IK)*RRSI(IK)
147!     *         /( (3.0*ALPHA/4.0)*(1.0/RKNUD+RL*ALPHA) )
148
149!                         INITIAL NUMBER OF H2SO4 MOLEC OF 1 DROPLET
150        ATOT=4.0*RPI*DENSO4ik(IK)*(RRSI_wet(IK)**3)/3.0 !attention: g and cm
151        ASO4(IK)=WH2(IK)*ATOT/MH2SO4 !attention: g
152!        ATOT=4.0*RPI*dens_aer(I,J)/1000.*(RRSI(IK)**3)/3.0
153!        ASO4=mfrac_H2SO4*ATOT/MH2SO4
154!                        INITIAL NUMBER OF H2O MOLEC OF 1 DROPLET
155        AH2O=(1.0-WH2(IK))*ATOT/MH2O !attention: g
156
157!       CHANGE OF THE NUMBER OF H2SO4 MOLEC OF 1 DROPLET DURING DT
158!       IT IS FOR KEM BUT THERE ARE OTHER WAYS
159
160      ENDIF 
161
162      ENDDO !loop over bins
163
164      END SUBROUTINE condens_evapor_rate_kelvin
165     
166!********************************************************************
167      SUBROUTINE condens_evapor_rate(R2SO4G,t_seri,pplay,ACTSO4,R2SO4, &
168                   & DENSO4,f_r_wet,FL,ASO4,DNDR)
169!
170!     INPUT:
171!     R2SO4: aerosol H2SO4 weight fraction (percent)
172!     ACTSO4: H2SO4 activity
173!     R2SO4G: number density of gaseous H2SO4 [molecules/cm3]
174!     t_seri: temperature (K)
175!     DENSO4: aerosol density (gr/cm3)
176
177      USE aerophys
178      USE infotrac_phy
179      USE YOMCST, ONLY : RPI
180      USE strataer_local_var_mod, ONLY : ALPH2SO4, RRSI
181
182      IMPLICIT NONE
183
184      ! input variables
185      REAL R2SO4G !H2SO4 number density [molecules/cm3]
186      REAL t_seri
187      REAL pplay
188      REAL ACTSO4
189      REAL R2SO4
190      REAL DENSO4
191      REAL f_r_wet
192     
193      ! output variables
194      REAL FL(nbtr_bin)
195      REAL ASO4(nbtr_bin)
196      REAL DNDR(nbtr_bin)
197
198      ! local variables
199      INTEGER IK
200      REAL ALPHA,CST
201      REAL WH2,RP,VTK,AA,FL1,RKNUD
202      REAL DND
203      REAL ATOT,AH2O
204      REAL RRSI_wet(nbtr_bin)
205      REAL FPATH
206
207! ///    MOLEC CONDENSATION GROWTH (DUE TO CHANGES IN H2SO4 AND SO H2O)
208!    ------------------------------------------------------------------
209!                                  EXCEPT CN
210!       RK:H2SO4 WEIGHT PERCENT DOESN'T CHANGE
211!     BE CAREFUL,H2SO4 WEIGHT PERCENTAGE
212
213!                   MOLECULAR ACCOMODATION OF H2SO4
214!     H2SO4 accommodation coefficient [condensation/evaporation]
215      ALPHA = ALPH2SO4
216!      FPLAIR=(2.281238E-5)*TAIR/PAIR
217!     1.E2 (m to cm),
218      CST=1.E2*2.281238E-5
219
220      ! compute local wet particle radius [cm]
221      RRSI_wet(:)=RRSI(:)*f_r_wet
222     
223!     Pruppa and Klett
224      FPATH=CST*t_seri/pplay
225
226     
227!     H2SO4 mass fraction in aerosol
228      WH2=R2SO4*1.0E-2
229      IF(WH2.EQ.0.0) RETURN
230!                               ACTIVITY COEFFICIENT(SEE GIAUQUE,1951)
231!                               AYERS ET AL (1980)
232!                                  (MU-MU0)
233      RP=-10156.0/t_seri +16.259-(ACTSO4*4.184)/(8.31441*t_seri)
234!                                  DROPLET H2SO4 PRESSURE IN DYN.CM-2
235      RP=EXP(RP)*1.01325E6/0.086
236!      RP=EXP(RP)*1.01325E6
237!                                  H2SO4 NUMBER DENSITY NEAR DROPLET
238!     R=8.31E7 DYN.CM.MOL-1*K-1
239!                               R/AVOGADRO NUMBER=DYN.CM.MOLEC-1*K-1
240      DND=RP*6.02E23/(8.31E7*t_seri)
241!                                  MEAN KINETIC VELOCITY
242!     DYN*CM*K/(K*GR)=(CM/SEC2)*CM
243!                                  IN CM/SEC
244      VTK=SQRT(8.0*BOLZ*t_seri/(RPI*MH2SO4))
245!                                 KELVIN EFFECT FACTOR
246!CK 20160613: bug fix, removed factor 250 (from original code by S. Bekki)
247!      AA =2.0*MH2O*72.0/(DENSO4*BOLZ*t_seri*250.0)
248      AA =2.0*MH2O*72.0/(DENSO4*BOLZ*t_seri)
249
250!     Loop on bin radius (RRSI in cm)
251      DO IK=1,nbtr_bin
252!                                 KELVIN EFFECT
253        DNDR(IK) =DND*EXP(AA/RRSI_wet(IK))
254
255        FL1=RPI*ALPHA*VTK*(R2SO4G-DNDR(IK))
256
257!       TURCO(1979) FOR HNO3:ALH2SO4 CONDENSATION= ALH2SO4 EVAPORATION
258!       RPI*R2*VTK IS EQUIVALENT TO DIFFUSION COEFFICIENT
259!       EXTENSION OF THE RELATION FOR DIFFUSION KINETICS
260!       KNUDSEN NUMBER FPATH/RRSI
261!       NEW VERSION (SEE NOTES)
262        RKNUD=FPATH/RRSI_wet(IK)
263!       SENFELD
264        FL(IK)=FL1*RRSI_wet(IK)**2*( 1.0 +RKNUD ) &
265     &     /( 1.0 +ALPHA/(2.0*RKNUD) +RKNUD )
266!       TURCO
267!        RL= (4.0/3.0 +0.71/RKNUD)/(1.0+1.0/RKNUD)
268!     *         +4.0*(1.0-ALPHA)/(3.0*ALPHA)
269!        FL=FL1*RRSI(IK)*RRSI(IK)
270!     *         /( (3.0*ALPHA/4.0)*(1.0/RKNUD+RL*ALPHA) )
271
272!                         INITIAL NUMBER OF H2SO4 MOLEC OF 1 DROPLET
273        ATOT=4.0*RPI*DENSO4*(RRSI_wet(IK)**3)/3.0 !attention: g and cm
274        ASO4(IK)=WH2*ATOT/MH2SO4 !attention: g
275!        ATOT=4.0*RPI*dens_aer(I,J)/1000.*(RRSI(IK)**3)/3.0
276!        ASO4=mfrac_H2SO4*ATOT/MH2SO4
277!                        INITIAL NUMBER OF H2O MOLEC OF 1 DROPLET
278        AH2O=(1.0-WH2)*ATOT/MH2O !attention: g
279
280!       CHANGE OF THE NUMBER OF H2SO4 MOLEC OF 1 DROPLET DURING DT
281!       IT IS FOR KEM BUT THERE ARE OTHER WAYS
282
283      ENDDO !loop over bins
284
285      END SUBROUTINE condens_evapor_rate
286
287!********************************************************************
288      SUBROUTINE condens_evapor_part(dt,FL,ASO4,f_r_wet,tr_seri)
289
290      USE aerophys
291      USE infotrac_phy
292      USE YOMCST, ONLY : RPI
293      USE strataer_local_var_mod, ONLY : RRSI,Vbin
294     
295      IMPLICIT NONE
296
297      ! input variables
298      REAL dt
299      REAL FL(nbtr_bin)
300      REAL ASO4(nbtr_bin)
301      REAL f_r_wet
302     
303      ! output variables
304      REAL tr_seri(nbtr)
305     
306      ! local variables
307      REAL tr_seri_new(nbtr)
308      INTEGER IK,JK,k
309      REAL Vnew
310      REAL RRSI_wet(nbtr_bin)
311      REAL Vbin_wet(nbtr_bin)
312      REAL sum_IK(nbtr_bin)
313      REAL ff(nbtr_bin,nbtr_bin)
314
315      tr_seri_new(:)=tr_seri(:)
316
317      ! compute local wet particle radius and volume
318      RRSI_wet(:)=RRSI(:)*f_r_wet
319      Vbin_wet(:)=Vbin(:)*f_r_wet**3 *1.e6 !Vbin_wet in cm3 (as Vnew)
320
321      ! compute distribution factor for particles of intermediate size (from Jacobson 1994, equation 13)
322      DO IK=1,nbtr_bin
323        Vnew=4.0*RPI*(RRSI_wet(IK)**3)/3.0*(1.+FL(IK)*dt/ASO4(IK))
324        ff(IK,:)=0.0
325        DO k=1, nbtr_bin
326          IF (k.LE.(nbtr_bin-1)) THEN
327            IF (Vbin_wet(k).LE.Vnew.AND.Vnew.LT.Vbin_wet(k+1)) THEN
328              ff(IK,k)= Vbin_wet(k)/Vnew*(Vbin_wet(k+1)-Vnew)/(Vbin_wet(k+1)-Vbin_wet(k))
329            ENDIF
330          ENDIF
331          IF (k.EQ.1.AND.Vnew.LE.Vbin_wet(k)) THEN
332            ff(IK,k)= 1.
333          ENDIF
334          IF (k.GT.1) THEN
335            IF (Vbin_wet(k-1).LT.Vnew.AND.Vnew.LT.Vbin_wet(k)) THEN
336              ff(IK,k)= 1.-ff(IK,k-1)
337            ENDIF
338          ENDIF
339          IF (k.EQ.nbtr_bin.AND.Vnew.GE.Vbin_wet(k)) THEN
340            ff(IK,k)= 1.
341          ENDIF
342        ENDDO
343        ! correction of ff for volume conservation
344        DO k=1, nbtr_bin         
345          ff(IK,k)=ff(IK,k)*Vnew/Vbin_wet(k)
346        ENDDO
347      ENDDO !loop over bins
348
349      DO IK=1, nbtr_bin
350        sum_IK(IK)=0.0
351        DO JK=1, nbtr_bin
352          sum_IK(IK)=sum_IK(IK)+tr_seri(JK+nbtr_sulgas)*ff(JK,IK)
353        ENDDO
354        ! compute new particle concentrations
355        tr_seri_new(IK+nbtr_sulgas)=sum_IK(IK)
356      ENDDO
357
358      tr_seri(:)=tr_seri_new(:)
359
360      END SUBROUTINE condens_evapor_part
361
362END MODULE cond_evap_tstep_mod
Note: See TracBrowser for help on using the repository browser.