source: LMDZ5/branches/testing/libf/phylmd/StratAer/cond_evap_tstep_mod.F90 @ 2787

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

Merged trunk changes r2664:2719 into testing branch
File size: 6.6 KB
Line 
1MODULE cond_evap_tstep_mod
2
3! based on UPMC aerosol model by Slimane Bekki
4! adapted for stratospheric sulfate aerosol in LMDZ by Christoph Kleinschmitt
5
6CONTAINS
7
8      SUBROUTINE condens_evapor_rate(R2SO4G,t_seri,pplay,ACTSO4,R2SO4, &
9                   & DENSO4,f_r_wet,RRSI,Vbin,FL,ASO4,DNDR)
10!
11!     INPUT:
12!     R2SO4: aerosol H2SO4 weight fraction (percent)
13!     ACTSO4: H2SO4 activity
14!     R2SO4G: number density of gaseous H2SO4 [molecules/cm3]
15!     t_seri: temperature (K)
16!     DENSO4: aerosol density (gr/cm3)
17
18      USE aerophys
19      USE infotrac
20      USE YOMCST, ONLY : RPI
21
22      IMPLICIT NONE
23
24      ! input variables
25      REAL R2SO4G !H2SO4 number density [molecules/cm3]
26      REAL t_seri
27      REAL pplay
28      REAL ACTSO4
29      REAL R2SO4
30      REAL DENSO4
31      REAL f_r_wet
32      REAL RRSI(nbtr_bin)
33      REAL Vbin(nbtr_bin)
34
35      ! output variables
36      REAL FL(nbtr_bin)
37      REAL ASO4(nbtr_bin)
38      REAL DNDR(nbtr_bin)
39
40      ! local variables
41      INTEGER IK
42      REAL ALPHA,CST
43      REAL WH2,RP,VTK,AA,FL1,RKNUD
44      REAL DND
45      REAL ATOT,AH2O
46      REAL RRSI_wet(nbtr_bin)
47      REAL Vbin_wet(nbtr_bin)
48      REAL MH2SO4,MH2O,BOLZ,FPATH
49
50! ///    MOLEC CONDENSATION GROWTH (DUE TO CHANGES IN H2SO4 AND SO H2O)
51!    ------------------------------------------------------------------
52!                                  EXCEPT CN
53!       RK:H2SO4 WEIGHT PERCENT DOESN'T CHANGE
54!     BE CAREFUL,H2SO4 WEIGHT PERCENTAGE
55
56!                   WEIGHT OF 1 MOLEC IN G
57      MH2O  =1000.*mH2Omol !18.016*1.66E-24
58      MH2SO4=1000.*mH2SO4mol !98.082*1.66E-24
59!                   BOLTZMANN CONSTANTE IN DYN.CM/K
60      BOLZ  =1.381E-16
61!                   MOLECULAR ACCOMODATION OF H2SO4
62!     raes and van dingen
63      ALPHA =0.1   
64!      FPLAIR=(2.281238E-5)*TAIR/PAIR
65!     1.E2 (m to cm),
66      CST=1.E2*2.281238E-5
67
68      ! compute local wet particle radius and volume
69      RRSI_wet(:)=RRSI(:)*f_r_wet
70      Vbin_wet(:)=Vbin(:)*f_r_wet**3
71
72!     Pruppa and Klett
73      FPATH=CST*t_seri/pplay
74
75     
76!     H2SO4 mass fraction in aerosol
77      WH2=R2SO4*1.0E-2
78      IF(WH2.EQ.0.0) RETURN
79!                               ACTIVITY COEFFICIENT(SEE GIAUQUE,1951)
80!                               AYERS ET AL (1980)
81!                                  (MU-MU0)
82      RP=-10156.0/t_seri +16.259-(ACTSO4*4.184)/(8.31441*t_seri)
83!                                  DROPLET H2SO4 PRESSURE IN DYN.CM-2
84      RP=EXP(RP)*1.01325E6/0.086
85!      RP=EXP(RP)*1.01325E6
86!                                  H2SO4 NUMBER DENSITY NEAR DROPLET
87!     R=8.31E7 DYN.CM.MOL-1*K-1
88!                               R/AVOGADRO NUMBER=DYN.CM.MOLEC-1*K-1
89      DND=RP*6.02E23/(8.31E7*t_seri)
90!                                  MEAN KINETIC VELOCITY
91!     DYN*CM*K/(K*GR)=(CM/SEC2)*CM
92!                                  IN CM/SEC
93      VTK=SQRT(8.0*BOLZ*t_seri/(RPI*MH2SO4))
94!                                 KELVIN EFFECT FACTOR
95!CK 20160613: bug fix, removed factor 250 (from original code by S. Bekki)
96!      AA =2.0*MH2O*72.0/(DENSO4*BOLZ*t_seri*250.0)
97      AA =2.0*MH2O*72.0/(DENSO4*BOLZ*t_seri)
98
99!     Loop on bin radius (RRSI in cm)
100      DO IK=1,nbtr_bin
101!                                 KELVIN EFFECT
102        DNDR(IK) =DND*EXP(AA/RRSI_wet(IK))
103
104        FL1=RPI*ALPHA*VTK*(R2SO4G-DNDR(IK))
105
106!       TURCO(1979) FOR HNO3:ALH2SO4 CONDENSATION= ALH2SO4 EVAPORATION
107!       RPI*R2*VTK IS EQUIVALENT TO DIFFUSION COEFFICIENT
108!       EXTENSION OF THE RELATION FOR DIFFUSION KINETICS
109!       KNUDSEN NUMBER FPATH/RRSI
110!       NEW VERSION (SEE NOTES)
111        RKNUD=FPATH/RRSI_wet(IK)
112!       SENFELD
113        FL(IK)=FL1*RRSI_wet(IK)**2*( 1.0 +RKNUD ) &
114     &     /( 1.0 +ALPHA/(2.0*RKNUD) +RKNUD )
115!       TURCO
116!        RL= (4.0/3.0 +0.71/RKNUD)/(1.0+1.0/RKNUD)
117!     *         +4.0*(1.0-ALPHA)/(3.0*ALPHA)
118!        FL=FL1*RRSI(IK)*RRSI(IK)
119!     *         /( (3.0*ALPHA/4.0)*(1.0/RKNUD+RL*ALPHA) )
120
121!                         INITIAL NUMBER OF H2SO4 MOLEC OF 1 DROPLET
122        ATOT=4.0*RPI*DENSO4*(RRSI_wet(IK)**3)/3.0 !attention: g and cm
123        ASO4(IK)=WH2*ATOT/MH2SO4 !attention: g
124!        ATOT=4.0*RPI*dens_aer(I,J)/1000.*(RRSI(IK)**3)/3.0
125!        ASO4=mfrac_H2SO4*ATOT/MH2SO4
126!                        INITIAL NUMBER OF H2O MOLEC OF 1 DROPLET
127        AH2O=(1.0-WH2)*ATOT/MH2O !attention: g
128
129!       CHANGE OF THE NUMBER OF H2SO4 MOLEC OF 1 DROPLET DURING DT
130!       IT IS FOR KEM BUT THERE ARE OTHER WAYS
131
132      ENDDO !loop over bins
133
134      END SUBROUTINE condens_evapor_rate
135
136!********************************************************************
137      SUBROUTINE cond_evap_part(dt,FL,ASO4,f_r_wet,RRSI,Vbin,tr_seri)
138
139      USE aerophys
140      USE infotrac
141      USE YOMCST, ONLY : RPI
142
143      IMPLICIT NONE
144
145      ! input variables
146      REAL dt
147      REAL FL(nbtr_bin)
148      REAL ASO4(nbtr_bin)
149      REAL f_r_wet
150      REAL RRSI(nbtr_bin)
151      REAL Vbin(nbtr_bin)
152
153      ! output variables
154      REAL tr_seri(nbtr)
155
156      ! local variables
157      REAL tr_seri_new(nbtr)
158      INTEGER IK,JK,k
159      REAL Vnew
160      REAL RRSI_wet(nbtr_bin)
161      REAL Vbin_wet(nbtr_bin)
162      REAL sum_IK(nbtr_bin)
163      REAL ff(nbtr_bin,nbtr_bin)
164
165      tr_seri_new(:)=tr_seri(:)
166
167      ! compute local wet particle radius and volume
168      RRSI_wet(:)=RRSI(:)*f_r_wet
169      Vbin_wet(:)=Vbin(:)*f_r_wet**3 *1.e6 !Vbin_wet in cm3 (as Vnew)
170
171      ! compute distribution factor for particles of intermediate size (from Jacobson 1994, equation 13)
172      DO IK=1,nbtr_bin
173        Vnew=4.0*RPI*(RRSI_wet(IK)**3)/3.0*(1.+FL(IK)*dt/ASO4(IK))
174        ff(IK,:)=0.0
175        DO k=1, nbtr_bin
176          IF (k.LE.(nbtr_bin-1)) THEN
177            IF (Vbin_wet(k).LE.Vnew.AND.Vnew.LT.Vbin_wet(k+1)) THEN
178              ff(IK,k)= Vbin_wet(k)/Vnew*(Vbin_wet(k+1)-Vnew)/(Vbin_wet(k+1)-Vbin_wet(k))
179            ENDIF
180          ENDIF
181          IF (k.EQ.1.AND.Vnew.LE.Vbin_wet(k)) THEN
182            ff(IK,k)= 1.
183          ENDIF
184          IF (k.GT.1) THEN
185            IF (Vbin_wet(k-1).LT.Vnew.AND.Vnew.LT.Vbin_wet(k)) THEN
186              ff(IK,k)= 1.-ff(IK,k-1)
187            ENDIF
188          ENDIF
189          IF (k.EQ.nbtr_bin.AND.Vnew.GE.Vbin_wet(k)) THEN
190            ff(IK,k)= 1.
191          ENDIF
192        ENDDO
193        ! correction of ff for volume conservation
194        DO k=1, nbtr_bin         
195          ff(IK,k)=ff(IK,k)*Vnew/Vbin_wet(k)
196        ENDDO
197      ENDDO !loop over bins
198
199      DO IK=1, nbtr_bin
200        sum_IK(IK)=0.0
201        DO JK=1, nbtr_bin
202          sum_IK(IK)=sum_IK(IK)+tr_seri(JK+nbtr_sulgas)*ff(JK,IK)
203        ENDDO
204        ! compute new particle concentrations
205        tr_seri_new(IK+nbtr_sulgas)=sum_IK(IK)
206      ENDDO
207
208      tr_seri(:)=tr_seri_new(:)
209
210      END SUBROUTINE cond_evap_part
211
212END MODULE cond_evap_tstep_mod
Note: See TracBrowser for help on using the repository browser.