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

Last change on this file since 4257 was 3677, checked in by oboucher, 5 years ago

Changed the way to initialise nbtr_bin and other dimensions and indices
in the StratAer? module based on infotrac_phy rather than infotrac.

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