Changeset 4950 for LMDZ6/trunk/libf/phylmd/StratAer
- Timestamp:
- May 22, 2024, 3:16:36 PM (5 weeks ago)
- Location:
- LMDZ6/trunk/libf/phylmd/StratAer
- Files:
-
- 11 edited
-
aer_sedimnt.F90 (modified) (2 diffs)
-
aerophys.F90 (modified) (1 diff)
-
coagulate.F90 (modified) (12 diffs)
-
cond_evap_tstep_mod.F90 (modified) (8 diffs)
-
micphy_tstep.F90 (modified) (9 diffs)
-
miecalc_aer.F90 (modified) (2 diffs)
-
nucleation_tstep_mod.F90 (modified) (2 diffs)
-
strataer_local_var_mod.F90 (modified) (6 diffs)
-
strataer_nuc_mod.F90 (modified) (2 diffs)
-
sulfate_aer_mod.F90 (modified) (3 diffs)
-
traccoag_mod.F90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/phylmd/StratAer/aer_sedimnt.F90
r3677 r4950 17 17 !----------------------------------------------------------------------- 18 18 19 USE phys_local_var_mod, ONLY: mdw, budg_sed_part, DENSO4, f_r_wet, vsed_aer 19 USE phys_local_var_mod, ONLY: mdw, budg_sed_part, DENSO4, DENSO4B, f_r_wet, f_r_wetB, vsed_aer 20 USE strataer_local_var_mod, ONLY: flag_new_strat_compo 20 21 USE dimphy, ONLY : klon,klev 21 22 USE infotrac_phy … … 89 90 90 91 ! stokes-velocity with cunnigham slip- flow correction 91 ZVAER(JL,JK,nb) = 2./9.*(DENSO4(JL,JK)*1000.-ZRHO)*RG/zvis(JL,JK)*(f_r_wet(JL,JK)*mdw(nb)/2.)**2.* & 92 (1.+ 2.*zlair(JL,JK)/(f_r_wet(JL,JK)*mdw(nb))*(1.257+0.4*EXP(-0.55*f_r_wet(JL,JK)*mdw(nb)/zlair(JL,JK)))) 93 92 IF(flag_new_strat_compo) THEN 93 ! stokes-velocity with cunnigham slip- flow correction 94 ZVAER(JL,JK,nb) = 2./9.*(DENSO4B(JL,JK,nb)*1000.-ZRHO)*RG/zvis(JL,JK)*(f_r_wetB(JL,JK,nb)*mdw(nb)/2.)**2.* & 95 (1.+ 2.*zlair(JL,JK)/(f_r_wetB(JL,JK,nb)*mdw(nb))*(1.257+0.4*EXP(-0.55*f_r_wetB(JL,JK,nb)*mdw(nb)/zlair(JL,JK)))) 96 ELSE 97 ZVAER(JL,JK,nb) = 2./9.*(DENSO4(JL,JK)*1000.-ZRHO)*RG/zvis(JL,JK)*(f_r_wet(JL,JK)*mdw(nb)/2.)**2.* & 98 (1.+ 2.*zlair(JL,JK)/(f_r_wet(JL,JK)*mdw(nb))*(1.257+0.4*EXP(-0.55*f_r_wet(JL,JK)*mdw(nb)/zlair(JL,JK)))) 99 ENDIF 100 94 101 ZSEDFLX(JL,nb)=ZVAER(JL,JK,nb)*ZRHO 95 102 ZSOLAERB(nb)=ZSOLAERB(nb)+ZDTGDP*ZSEDFLX(JL,nb) -
LMDZ6/trunk/libf/phylmd/StratAer/aerophys.F90
r4601 r4950 5 5 IMPLICIT NONE 6 6 ! 7 REAL,PARAMETER :: ropx=1500.0 ! default aerosol particle mass density [kg/m3] 8 REAL,PARAMETER :: dens_aer_dry=1848.682308 ! dry aerosol particle mass density at T_0=293K[kg/m3] 9 REAL,PARAMETER :: dens_aer_ref=1483.905336 ! aerosol particle mass density with 75% H2SO4 at T_0=293K[kg/m3] 10 REAL,PARAMETER :: mdwmin=0.002e-6 ! dry diameter of smallest aerosol particles [m] 11 REAL,PARAMETER :: V_rat=2.0 ! volume ratio of neighboring size bins 12 REAL,PARAMETER :: mfrac_H2SO4=0.75 ! default mass fraction of H2SO4 in the aerosol 13 REAL, PARAMETER :: mAIRmol=28.949*1.66E-27 ! Average mass of an air molecule [kg] 14 REAL, PARAMETER :: mH2Omol=18.016*1.66E-27 ! Mass of an H2O molecule [kg] 15 REAL, PARAMETER :: mH2SO4mol=98.082*1.66E-27! Mass of an H2SO4 molecule [kg] 16 REAL, PARAMETER :: mSO2mol=64.06*1.66E-27 ! Mass of an SO2 molecule [kg] 17 REAL, PARAMETER :: mSatom=32.06*1.66E-27 ! Mass of a S atom [kg] 18 REAL, PARAMETER :: mOCSmol=60.07*1.66E-27 ! Mass of an OCS molecule [kg] 19 REAL, PARAMETER :: mClatom=35.45*1.66E-27 ! Mass of an Cl atom [kg] 20 REAL, PARAMETER :: mHClmol=36.46*1.66E-27 ! Mass of an HCl molecule [kg] 21 REAL, PARAMETER :: mBratom=79.90*1.66E-27 ! Mass of an Br atom [kg] 22 REAL, PARAMETER :: mHBrmol=80.92*1.66E-27 ! Mass of an HBr molecule [kg] 23 REAL, PARAMETER :: mNOmol=30.01*1.66E-27 ! Mass of an NO molecule [kg] 24 REAL, PARAMETER :: mNO2mol=46.01*1.66E-27 ! Mass of an NO2 molecule [kg] 25 REAL, PARAMETER :: mNatome=14.0067*1.66E-27 ! Mass of an N atome [kg] 7 REAL,PARAMETER :: ropx=1500.0 ! default aerosol particle mass density [kg/m3] 8 REAL,PARAMETER :: dens_aer_dry=1848.682308 ! dry aerosol particle mass density at T_0=293K[kg/m3] 9 REAL,PARAMETER :: dens_aer_ref=1483.905336 ! aerosol particle mass density with 75% H2SO4 at T_0=293K[kg/m3] 10 REAL,PARAMETER :: mdwmin=0.002e-6 ! dry diameter of smallest aerosol particles [m] 11 REAL,PARAMETER :: V_rat=2.0 ! volume ratio of neighboring size bins 12 REAL,PARAMETER :: mfrac_H2SO4=0.75 ! default mass fraction of H2SO4 in the aerosol 13 REAL, PARAMETER :: mAIRmol=28.949*1.66E-27 ! Average mass of an air molecule [kg] 14 REAL, PARAMETER :: mH2Omol=18.016*1.66E-27 ! Mass of an H2O molecule [kg] 15 REAL, PARAMETER :: mH2SO4mol=98.082*1.66E-27! Mass of an H2SO4 molecule [kg] 16 REAL, PARAMETER :: mSO2mol=64.06*1.66E-27 ! Mass of an SO2 molecule [kg] 17 REAL, PARAMETER :: mSatom=32.06*1.66E-27 ! Mass of a S atom [kg] 18 REAL, PARAMETER :: mOCSmol=60.07*1.66E-27 ! Mass of an OCS molecule [kg] 19 REAL, PARAMETER :: mClatom=35.45*1.66E-27 ! Mass of an Cl atom [kg] 20 REAL, PARAMETER :: mHClmol=36.46*1.66E-27 ! Mass of an HCl molecule [kg] 21 REAL, PARAMETER :: mBratom=79.90*1.66E-27 ! Mass of an Br atom [kg] 22 REAL, PARAMETER :: mHBrmol=80.92*1.66E-27 ! Mass of an HBr molecule [kg] 23 REAL, PARAMETER :: mNOmol=30.01*1.66E-27 ! Mass of an NO molecule [kg] 24 REAL, PARAMETER :: mNO2mol=46.01*1.66E-27 ! Mass of an NO2 molecule [kg] 25 REAL, PARAMETER :: mNatome=14.0067*1.66E-27 ! Mass of an N atome [kg] 26 REAL, PARAMETER :: rgas=8.3145 ! molar gas cste (J⋅K−1⋅mol−1=m3⋅Pa⋅K−1⋅mol−1=kg⋅m2⋅s−2⋅K−1⋅mol−1) 27 ! 28 REAL, PARAMETER :: MH2O =1000.*mH2Omol ! Mass of 1 molec [g] (18.016*1.66E-24) 29 REAL, PARAMETER :: MH2SO4=1000.*mH2SO4mol ! Mass of 1 molec [g] (98.082*1.66E-24) 30 REAL, PARAMETER :: BOLZ =1.381E-16 ! Boltzmann constant [dyn.cm/K] 26 31 ! 27 32 END MODULE aerophys -
LMDZ6/trunk/libf/phylmd/StratAer/coagulate.F90
r4762 r4950 26 26 USE aerophys 27 27 USE infotrac_phy 28 USE phys_local_var_mod, ONLY: DENSO4, f_r_wet 29 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 30 30 31 IMPLICIT NONE 31 32 … … 43 44 ! local variables in coagulation routine 44 45 INTEGER :: i,j,k,nb,ilon,ilev 45 REAL, DIMENSION(nbtr_bin) :: radius ! aerosol particle radius in each bin [m] 46 REAL, DIMENSION(nbtr_bin) :: radiusdry ! dry aerosol particle radius in each bin [m] 47 REAL, DIMENSION(nbtr_bin) :: radiuswet ! wet aerosol particle radius in each bin [m] 46 48 REAL, DIMENSION(klon,klev,nbtr_bin) :: tr_t ! Concentration Traceur at time t [U/KgA] 47 49 REAL, DIMENSION(klon,klev,nbtr_bin) :: tr_tp1 ! Concentration Traceur at time t+1 [U/KgA] 48 50 REAL, DIMENSION(nbtr_bin,nbtr_bin,nbtr_bin) :: ff ! Volume fraction of intermediate particles 49 REAL, DIMENSION(nbtr_bin) :: V ! Volume of bins 51 REAL, DIMENSION(nbtr_bin) :: Vdry ! Volume dry of bins 52 REAL, DIMENSION(nbtr_bin) :: Vwet ! Volume wet of bins 50 53 REAL, DIMENSION(nbtr_bin,nbtr_bin) :: Vij ! Volume sum of i and j 51 54 REAL :: eta ! Dynamic viscosity of air … … 82 85 include "YOMCST.h" 83 86 84 DO i=1, nbtr_bin 85 radius(i)=mdw(i)/2. 86 V(i)= radius(i)**3. !neglecting factor 4*RPI/3 87 ENDDO 88 89 DO j=1, nbtr_bin 90 DO i=1, nbtr_bin 91 Vij(i,j)= V(i)+V(j) 92 ENDDO 93 ENDDO 94 87 ! ff(i,j,k): Volume fraction of Vi,j that is partitioned to each model bin k 88 ! just need to be calculated in model initialization because mdw(:) size is fixed 89 ! no need to recalculate radius, Vdry, Vij, and ff every timestep because it is for 90 ! dry aerosols 91 DO i=1, nbtr_bin 92 radiusdry(i)=mdw(i)/2. 93 Vdry(i)=radiusdry(i)**3. !neglecting factor 4*RPI/3 94 Vwet(i)=0.0 95 ENDDO 96 97 DO j=1, nbtr_bin 98 DO i=1, nbtr_bin 99 Vij(i,j)= Vdry(i)+Vdry(j) 100 ENDDO 101 ENDDO 102 95 103 !--pre-compute the f(i,j,k) from Jacobson equation 13 96 104 ff=0.0 … … 100 108 IF (k.EQ.1) THEN 101 109 ff(i,j,k)= 0.0 102 ELSEIF (k.GT.1.AND.V (k-1).LT.Vij(i,j).AND.Vij(i,j).LT.V(k)) THEN110 ELSEIF (k.GT.1.AND.Vdry(k-1).LT.Vij(i,j).AND.Vij(i,j).LT.Vdry(k)) THEN 103 111 ff(i,j,k)= 1.-ff(i,j,k-1) 104 112 ELSEIF (k.EQ.nbtr_bin) THEN 105 IF (Vij(i,j).GE. v(k)) THEN113 IF (Vij(i,j).GE.Vdry(k)) THEN 106 114 ff(i,j,k)= 1. 107 115 ELSE 108 116 ff(i,j,k)= 0.0 109 117 ENDIF 110 ELSEIF (k.LE.(nbtr_bin-1).AND.V (k).LE.Vij(i,j).AND.Vij(i,j).LT.V(k+1)) THEN111 ff(i,j,k)= V (k)/Vij(i,j)*(V(k+1)-Vij(i,j))/(V(k+1)-V(k))118 ELSEIF (k.LE.(nbtr_bin-1).AND.Vdry(k).LE.Vij(i,j).AND.Vij(i,j).LT.Vdry(k+1)) THEN 119 ff(i,j,k)= Vdry(k)/Vij(i,j)*(Vdry(k+1)-Vij(i,j))/(Vdry(k+1)-Vdry(k)) 112 120 ENDIF 113 121 ENDDO 114 122 ENDDO 115 123 ENDDO 116 124 ! End of just need to be calculated at initialization because mdw(:) size is fixed 125 117 126 DO ilon=1, klon 118 127 DO ilev=1, klev … … 120 129 IF (is_strato(ilon,ilev)) THEN 121 130 !compute actual wet particle radius & volume for every grid box 122 DO i=1, nbtr_bin 123 radius(i)=f_r_wet(ilon,ilev)*mdw(i)/2. 124 V(i)= radius(i)**3. !neglecting factor 4*RPI/3 125 ENDDO 126 131 IF(flag_new_strat_compo) THEN 132 DO i=1, nbtr_bin 133 radiuswet(i)=f_r_wetB(ilon,ilev,i)*mdw(i)/2. 134 Vwet(i)= radiuswet(i)**3. !neglecting factor 4*RPI/3 135 !! Vwet(i)= Vdry(i)*(f_r_wetB(ilon,ilev,i)**3) 136 ENDDO 137 ELSE 138 DO i=1, nbtr_bin 139 radiuswet(i)=f_r_wet(ilon,ilev)*mdw(i)/2. 140 Vwet(i)= radiuswet(i)**3. !neglecting factor 4*RPI/3 141 !! Vwet(i)= Vdry(i)*(f_r_wet(ilon,ilev)**3) 142 ENDDO 143 ENDIF 144 127 145 !--Calculations for the coagulation kernel--------------------------------------------------------- 128 146 … … 150 168 Di=0.0 151 169 DO i=1, nbtr_bin 152 Kn(i)=mnfrpth/radius(i)153 Di(i)=RKBOL*t_seri(ilon,ilev)/(6.*RPI*radius(i)*eta)*(1.+Kn(i)*(1.249+0.42*exp(-0.87/Kn(i))))170 Kn(i)=mnfrpth/radiuswet(i) 171 Di(i)=RKBOL*t_seri(ilon,ilev)/(6.*RPI*radiuswet(i)*eta)*(1.+Kn(i)*(1.249+0.42*exp(-0.87/Kn(i)))) 154 172 ENDDO 155 173 156 174 !--pre-compute the thermal velocity of a particle thvelpar(i) from equation 20 157 175 thvelpar=0.0 158 DO i=1, nbtr_bin 159 m_par(i)=4./3.*RPI*radius(i)**3.*DENSO4(ilon,ilev)*1000. 160 thvelpar(i)=sqrt(8.*RKBOL*t_seri(ilon,ilev)/(RPI*m_par(i))) 161 ENDDO 176 IF(flag_new_strat_compo) THEN 177 DO i=1, nbtr_bin 178 m_par(i)=4./3.*RPI*radiuswet(i)**3.*DENSO4B(ilon,ilev,i)*1000. 179 thvelpar(i)=sqrt(8.*RKBOL*t_seri(ilon,ilev)/(RPI*m_par(i))) 180 ENDDO 181 ELSE 182 DO i=1, nbtr_bin 183 m_par(i)=4./3.*RPI*radiuswet(i)**3.*DENSO4(ilon,ilev)*1000. 184 thvelpar(i)=sqrt(8.*RKBOL*t_seri(ilon,ilev)/(RPI*m_par(i))) 185 ENDDO 186 ENDIF 162 187 163 188 !--pre-compute the particle mean free path mfppar(i) from equation 22 … … 171 196 delta=0.0 172 197 DO i=1, nbtr_bin 173 delta(i)=((2.*radius(i)+mfppar(i))**3.-(4.*radius(i)**2.+mfppar(i)**2.)**1.5)/ & 174 & (6.*radius(i)*mfppar(i))-2.*radius(i) 175 ENDDO 176 198 delta(i)=((2.*radiuswet(i)+mfppar(i))**3.-(4.*radiuswet(i)**2.+mfppar(i)**2.)**1.5)/ & 199 & (6.*radiuswet(i)*mfppar(i))-2.*radiuswet(i) 200 ENDDO 201 202 ! beta(i,j): coagulation kernel (rate coefficient) of 2 colliding particles i,j 177 203 !--pre-compute the beta(i,j) from equation 17 in Jacobson 178 204 num=0.0 … … 180 206 DO i=1, nbtr_bin 181 207 ! 182 num=4.*RPI*(radius(i)+radius(j))*(Di(i)+Di(j))183 denom=(radius(i)+radius(j))/(radius(i)+radius(j)+sqrt(delta(i)**2.+delta(j)**2.))+ &184 & 4.*(Di(i)+Di(j))/(sqrt(thvelpar(i)**2.+thvelpar(j)**2.)*(radius(i)+radius(j)))185 beta(i,j)=num/denom208 num=4.*RPI*(radiuswet(i)+radiuswet(j))*(Di(i)+Di(j)) 209 denom=(radiuswet(i)+radiuswet(j))/(radiuswet(i)+radiuswet(j)+sqrt(delta(i)**2.+delta(j)**2.))+ & 210 & 4.*(Di(i)+Di(j))/(sqrt(thvelpar(i)**2.+thvelpar(j)**2.)*(radiuswet(i)+radiuswet(j))) 211 beta(i,j)=num/denom 186 212 ! 187 213 !--compute enhancement factor due to van der Waals forces 188 214 IF (ok_vdw .EQ. 0) THEN !--no enhancement factor 189 Evdw=1.0215 Evdw=1.0 190 216 ELSEIF (ok_vdw .EQ. 1) THEN !--E(0) case 191 AvdWi = AvdW/(RKBOL*t_seri(ilon,ilev))*(4.*radius(i)*radius(j))/(radius(i)+radius(j))**2.192 xvdW = LOG(1.+AvdWi)193 EvdW = 1. + avdW1*xvdW + avdW3*xvdW**3217 AvdWi = AvdW/(RKBOL*t_seri(ilon,ilev))*(4.*radiuswet(i)*radiuswet(j))/(radiuswet(i)+radiuswet(j))**2. 218 xvdW = LOG(1.+AvdWi) 219 EvdW = 1. + avdW1*xvdW + avdW3*xvdW**3 194 220 ELSEIF (ok_vdw .EQ. 2) THEN !--E(infinity) case 195 AvdWi = AvdW/(RKBOL*t_seri(ilon,ilev))*(4.*radius(i)*radius(j))/(radius(i)+radius(j))**2.196 xvdW = LOG(1.+AvdWi)197 EvdW = 1. + SQRT(AvdWi/3.)/(1.+bvdW0*SQRT(AvdWi)) + bvdW1*xvdW + bvdW3*xvdW**3.221 AvdWi = AvdW/(RKBOL*t_seri(ilon,ilev))*(4.*radiuswet(i)*radiuswet(j))/(radiuswet(i)+radiuswet(j))**2. 222 xvdW = LOG(1.+AvdWi) 223 EvdW = 1. + SQRT(AvdWi/3.)/(1.+bvdW0*SQRT(AvdWi)) + bvdW1*xvdW + bvdW3*xvdW**3. 198 224 ENDIF 199 225 ! … … 209 235 denom=0.0 210 236 DO j=1, nbtr_bin 211 denom=denom+(1.-ff(k,j,k))*beta(k,j)*tr_t(ilon,ilev,j) 237 ! fraction of coagulation of k and j that is not giving k 238 denom=denom+(1.-ff(k,j,k))*beta(k,j)*tr_t(ilon,ilev,j) 212 239 ENDDO 213 240 … … 219 246 num=0.0 220 247 DO j=1, k 221 numi=0.0 222 DO i=1, k-1 223 numi=numi+ff(i,j,k)*beta(i,j)*V(i)*tr_tp1(ilon,ilev,i)*tr_t(ilon,ilev,j) 248 numi=0.0 249 DO i=1, k-1 250 ! 251 ! see Jacobson: " In order to conserve volume and volume concentration (which 252 ! coagulation physically does) while giving up some accuracy in number concentration" 253 ! 254 ! Coagulation of i and j giving k 255 ! with V(i) and then V(j) because it considers i,j and j,i with the double loop 256 ! 257 ! BUT WHY WET VOLUME V(i) in old STRATAER? tracers are already dry aerosols and coagulation 258 ! kernel beta(i,j) accounts for wet aerosols -> reply below 259 ! 260 ! numi=numi+ff(i,j,k)*beta(i,j)*V(i)*tr_tp1(ilon,ilev,i)*tr_t(ilon,ilev,j) 261 numi=numi+ff(i,j,k)*beta(i,j)*Vdry(i)*tr_tp1(ilon,ilev,i)*tr_t(ilon,ilev,j) 262 ENDDO 263 num=num+numi 224 264 ENDDO 225 num=num+numi226 ENDDO227 265 228 266 !--calculate new concentration of other bins 229 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)=(V(k)*tr_t(ilon,ilev,k)+pdtcoag*num)/( (1.+pdtcoag*denom)*V(k) ) 268 tr_tp1(ilon,ilev,k)=(Vdry(k)*tr_t(ilon,ilev,k)+pdtcoag*num)/( (1.+pdtcoag*denom)*Vdry(k) ) 269 ! 270 ! In constant composition (no dependency on aerosol size because no kelvin effect) 271 ! V(l)= (f_r_wet(ilon,ilev)**3)*((mdw(l)/2.)**3) = (f_r_wet(ilon,ilev)**3)*Vdry(i) 272 ! so numi and num are proportional (f_r_wet(ilon,ilev)**3) 273 ! and so 274 ! tr_tp1(ilon,ilev,k)=(V(k)*tr_t(ilon,ilev,k)+pdtcoag*num)/( (1.+pdtcoag*denom)*V(k) ) 275 ! =(Vdry(k)*tr_t(ilon,ilev,k)+pdtcoag*num_dry)/( (1.+pdtcoag*denom)*Vdry(k) ) 276 ! with num_dry=...beta(i,j)*Vdry(i)*.... 277 ! so in old STRATAER (.not.flag_new_strat_compo), it was correct 230 278 ENDIF 231 279 … … 234 282 !--convert tracer concentration back from [number/m3] to [number/KgA] and write into tr_seri 235 283 DO i=1, nbtr_bin 236 tr_seri(ilon,ilev,i+nbtr_sulgas) = tr_tp1(ilon,ilev,i) / zrho284 tr_seri(ilon,ilev,i+nbtr_sulgas) = tr_tp1(ilon,ilev,i) / zrho 237 285 ENDDO 238 286 … … 240 288 ENDDO !--end of loop klev 241 289 ENDDO !--end of loop klon 290 ! ********************************************* 242 291 243 292 END SUBROUTINE COAGULATE -
LMDZ6/trunk/libf/phylmd/StratAer/cond_evap_tstep_mod.F90
r3677 r4950 9 9 CONTAINS 10 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 !******************************************************************** 11 167 SUBROUTINE condens_evapor_rate(R2SO4G,t_seri,pplay,ACTSO4,R2SO4, & 12 & DENSO4,f_r_wet, RRSI,Vbin,FL,ASO4,DNDR)168 & DENSO4,f_r_wet,FL,ASO4,DNDR) 13 169 ! 14 170 ! INPUT: … … 22 178 USE infotrac_phy 23 179 USE YOMCST, ONLY : RPI 180 USE strataer_local_var_mod, ONLY : ALPH2SO4, RRSI 24 181 25 182 IMPLICIT NONE … … 33 190 REAL DENSO4 34 191 REAL f_r_wet 35 REAL RRSI(nbtr_bin) 36 REAL Vbin(nbtr_bin) 37 192 38 193 ! output variables 39 194 REAL FL(nbtr_bin) … … 48 203 REAL ATOT,AH2O 49 204 REAL RRSI_wet(nbtr_bin) 50 REAL Vbin_wet(nbtr_bin) 51 REAL MH2SO4,MH2O,BOLZ,FPATH 205 REAL FPATH 52 206 53 207 ! /// MOLEC CONDENSATION GROWTH (DUE TO CHANGES IN H2SO4 AND SO H2O) … … 57 211 ! BE CAREFUL,H2SO4 WEIGHT PERCENTAGE 58 212 59 ! WEIGHT OF 1 MOLEC IN G60 MH2O =1000.*mH2Omol !18.016*1.66E-2461 MH2SO4=1000.*mH2SO4mol !98.082*1.66E-2462 ! BOLTZMANN CONSTANTE IN DYN.CM/K63 BOLZ =1.381E-1664 213 ! MOLECULAR ACCOMODATION OF H2SO4 65 ! raes and van dingen66 ALPHA = 0.1214 ! H2SO4 accommodation coefficient [condensation/evaporation] 215 ALPHA = ALPH2SO4 67 216 ! FPLAIR=(2.281238E-5)*TAIR/PAIR 68 217 ! 1.E2 (m to cm), 69 218 CST=1.E2*2.281238E-5 70 219 71 ! compute local wet particle radius and volume220 ! compute local wet particle radius [cm] 72 221 RRSI_wet(:)=RRSI(:)*f_r_wet 73 Vbin_wet(:)=Vbin(:)*f_r_wet**3 74 222 75 223 ! Pruppa and Klett 76 224 FPATH=CST*t_seri/pplay … … 138 286 139 287 !******************************************************************** 140 SUBROUTINE cond _evap_part(dt,FL,ASO4,f_r_wet,RRSI,Vbin,tr_seri)288 SUBROUTINE condens_evapor_part(dt,FL,ASO4,f_r_wet,tr_seri) 141 289 142 290 USE aerophys 143 291 USE infotrac_phy 144 292 USE YOMCST, ONLY : RPI 145 293 USE strataer_local_var_mod, ONLY : RRSI,Vbin 294 146 295 IMPLICIT NONE 147 296 … … 151 300 REAL ASO4(nbtr_bin) 152 301 REAL f_r_wet 153 REAL RRSI(nbtr_bin) 154 REAL Vbin(nbtr_bin) 155 302 156 303 ! output variables 157 304 REAL tr_seri(nbtr) 158 305 159 306 ! local variables 160 307 REAL tr_seri_new(nbtr) … … 211 358 tr_seri(:)=tr_seri_new(:) 212 359 213 END SUBROUTINE cond _evap_part360 END SUBROUTINE condens_evapor_part 214 361 215 362 END MODULE cond_evap_tstep_mod -
LMDZ6/trunk/libf/phylmd/StratAer/micphy_tstep.F90
r4601 r4950 8 8 USE aerophys 9 9 USE infotrac_phy, ONLY : nbtr_bin, nbtr_sulgas, nbtr, id_H2SO4_strat 10 USE phys_local_var_mod, ONLY: mdw, budg_3D_nucl, budg_3D_cond_evap, budg_h2so4_to_part, R2SO4, DENSO4, f_r_wet 10 USE phys_local_var_mod, ONLY: mdw, budg_3D_nucl, budg_3D_cond_evap, budg_h2so4_to_part, R2SO4, DENSO4, & 11 f_r_wet, R2SO4B, DENSO4B, f_r_wetB 11 12 USE nucleation_tstep_mod 12 13 USE cond_evap_tstep_mod … … 14 15 USE YOMCST, ONLY : RPI, RD, RG 15 16 USE print_control_mod, ONLY: lunout 16 USE strataer_local_var_mod 17 USE strataer_local_var_mod ! contains also RRSI and Vbin 17 18 18 19 IMPLICIT NONE … … 35 36 REAL :: ntot !total number of molecules in the critical cluster (ntot>4) 36 37 REAL :: x ! molefraction of H2SO4 in the critical cluster 37 REAL Vbin(nbtr_bin)38 38 REAL a_xm, b_xm, c_xm 39 39 REAL PDT, dt 40 40 REAL H2SO4_init 41 41 REAL ACTSO4(klon,klev) 42 REAL RRSI(nbtr_bin)43 42 REAL nucl_rate 44 43 REAL cond_evap_rate … … 48 47 REAL DNDR(nbtr_bin) 49 48 REAL H2SO4_sat 50 51 DO it=1,nbtr_bin 52 Vbin(it)=4.0*RPI*((mdw(it)/2.)**3)/3.0 53 ENDDO 54 49 REAL R2SO4ik(nbtr_bin), DENSO4ik(nbtr_bin), f_r_wetik(nbtr_bin) 50 55 51 !coefficients for H2SO4 density parametrization used for nucleation if ntot<4 56 52 a_xm = 0.7681724 + 1.*(2.1847140 + 1.*(7.1630022 + 1.*(-44.31447 + & … … 61 57 & 1.*(7.990811e-4 + 1.*(-7.458060e-4 + 1.*2.58139e-4 ))))) 62 58 63 ! STRAACT (R2SO4, t_seri -> H2SO4 activity coefficient (ACTSO4)) for cond/evap 64 CALL STRAACT(ACTSO4) 65 66 ! compute particle radius in cm RRSI from diameter in m 67 DO it=1,nbtr_bin 68 RRSI(it)=mdw(it)/2.*100. 69 ENDDO 70 59 IF(.not.flag_new_strat_compo) THEN 60 ! STRAACT (R2SO4, t_seri -> H2SO4 activity coefficient (ACTSO4)) for cond/evap 61 CALL STRAACT(ACTSO4) 62 ENDIF 63 71 64 DO ilon=1, klon 72 65 ! … … 104 97 ENDIF 105 98 ! compute cond/evap rate in kg(H2SO4)/kgA/s 106 CALL condens_evapor_rate(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev), & 107 & ACTSO4(ilon,ilev),R2SO4(ilon,ilev),DENSO4(ilon,ilev),f_r_wet(ilon,ilev), & 108 & RRSI,Vbin,FL,ASO4,DNDR) 99 IF(flag_new_strat_compo) THEN 100 R2SO4ik(:) = R2SO4B(ilon,ilev,:) 101 DENSO4ik(:) = DENSO4B(ilon,ilev,:) 102 f_r_wetik(:) = f_r_wetB(ilon,ilev,:) 103 CALL condens_evapor_rate_kelvin(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev), & 104 & R2SO4(ilon,ilev),DENSO4(ilon,ilev),f_r_wet(ilon,ilev), & 105 & R2SO4ik,DENSO4ik,f_r_wetik,FL,ASO4,DNDR) 106 ELSE 107 CALL condens_evapor_rate(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev), & 108 & ACTSO4(ilon,ilev),R2SO4(ilon,ilev),DENSO4(ilon,ilev),f_r_wet(ilon,ilev), & 109 & FL,ASO4,DNDR) 110 ENDIF 109 111 ! Compute H2SO4 saturate vapor for big particules 110 112 H2SO4_sat = DNDR(nbtr_bin)/(pplay(ilon,ilev)/t_seri(ilon,ilev)/RD/1.E6/mH2SO4mol) … … 127 129 tr_seri(ilon,ilev,id_H2SO4_strat)=MAX(0.,tr_seri(ilon,ilev,id_H2SO4_strat)-(nucl_rate+cond_evap_rate)*dt) 128 130 ! apply cond to bins 129 CALL cond _evap_part(dt,FL,ASO4,f_r_wet(ilon,ilev),RRSI,Vbin,tr_seri(ilon,ilev,:))131 CALL condens_evapor_part(dt,FL,ASO4,f_r_wet(ilon,ilev),tr_seri(ilon,ilev,:)) 130 132 ! apply nucl. to bins 131 CALL nucleation_part(nucl_rate,ntot,x,dt, Vbin,tr_seri(ilon,ilev,:))133 CALL nucleation_part(nucl_rate,ntot,x,dt,tr_seri(ilon,ilev,:)) 132 134 ! compute fluxes as diagnostic in [kg(S)/m2/layer/s] (now - for evap and + for cond) 133 135 budg_3D_cond_evap(ilon,ilev)=budg_3D_cond_evap(ilon,ilev)+mSatom/mH2SO4mol & … … 142 144 & *pplay(ilon,ilev)/t_seri(ilon,ilev)/RD/1.E6/mH2SO4mol 143 145 ! compute cond/evap rate in kg(H2SO4)/kgA/s (now only evap for pdtphys) 144 CALL condens_evapor_rate(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev), & 145 & ACTSO4(ilon,ilev),R2SO4(ilon,ilev),DENSO4(ilon,ilev),f_r_wet(ilon,ilev), & 146 & RRSI,Vbin,FL,ASO4,DNDR) 146 IF(flag_new_strat_compo) THEN 147 CALL condens_evapor_rate_kelvin(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev), & 148 & R2SO4(ilon,ilev),DENSO4(ilon,ilev),f_r_wet(ilon,ilev), & 149 & R2SO4ik,DENSO4ik,f_r_wetik,FL,ASO4,DNDR) 150 ELSE 151 CALL condens_evapor_rate(rhoa,t_seri(ilon,ilev),pplay(ilon,ilev), & 152 & ACTSO4(ilon,ilev),R2SO4(ilon,ilev),DENSO4(ilon,ilev),f_r_wet(ilon,ilev), & 153 & FL,ASO4,DNDR) 154 ENDIF 147 155 ! limit evaporation (negative FL) over one physics time step to H2SO4 content of the droplet 148 156 DO it=1,nbtr_bin … … 159 167 tr_seri(ilon,ilev,id_H2SO4_strat)=MAX(0.,tr_seri(ilon,ilev,id_H2SO4_strat)-evap_rate*pdtphys) 160 168 ! apply evap to bins 161 CALL cond _evap_part(pdtphys,FL,ASO4,f_r_wet(ilon,ilev),RRSI,Vbin,tr_seri(ilon,ilev,:))169 CALL condens_evapor_part(pdtphys,FL,ASO4,f_r_wet(ilon,ilev),tr_seri(ilon,ilev,:)) 162 170 ! compute fluxes as diagnostic in [kg(S)/m2/layer/s] (now - for evap and + for cond) 163 171 budg_3D_cond_evap(ilon,ilev)=budg_3D_cond_evap(ilon,ilev)+mSatom/mH2SO4mol & -
LMDZ6/trunk/libf/phylmd/StratAer/miecalc_aer.F90
r3677 r4950 16 16 17 17 USE phys_local_var_mod, ONLY: tr_seri, mdw, alpha_bin, piz_bin, cg_bin 18 USE aerophys 18 USE aerophys, ONLY: dens_aer_dry, dens_aer_ref, V_rat 19 19 USE aero_mod 20 20 USE infotrac_phy, ONLY : nbtr, nbtr_bin, nbtr_sulgas, id_SO2_strat … … 226 226 40000.000, 0.2500, 1.48400, 1.0000E-08, & 227 227 50000.000, 0.2000, 1.49800, 1.0000E-08 /), (/nb_lambda_h2so4,4/), order=(/2,1/) ) 228 229 !--initialising dry diameters to geometrically spaced mass/volume (see Jacobson 1994) 230 mdw(1)=mdwmin 231 IF (V_rat.LT.1.62) THEN ! compensate for dip in second bin for lower volume ratio 232 mdw(2)=mdw(1)*2.**(1./3.) 233 DO it=3, nbtr_bin 234 mdw(it)=mdw(it-1)*V_rat**(1./3.) 235 ENDDO 236 ELSE 237 DO it=2, nbtr_bin 238 mdw(it)=mdw(it-1)*V_rat**(1./3.) 239 ENDDO 240 ENDIF 241 WRITE(lunout,*) 'init mdw=', mdw 242 228 243 229 !--compute particle radius for a composition of 75% H2SO4 / 25% H2O at T=293K 244 230 DO bin_number=1, nbtr_bin -
LMDZ6/trunk/libf/phylmd/StratAer/nucleation_tstep_mod.F90
r4912 r4950 70 70 !-------------------------------------------------------------------------------------------------- 71 71 72 SUBROUTINE nucleation_part(nucl_rate,ntot,x,dt, Vbin,tr_seri)72 SUBROUTINE nucleation_part(nucl_rate,ntot,x,dt,tr_seri) 73 73 74 74 USE aerophys 75 75 USE infotrac_phy 76 76 USE strataer_local_var_mod, ONLY : Vbin 77 77 78 IMPLICIT NONE 78 79 … … 82 83 REAL x ! mole raction of H2SO4 in the critical cluster 83 84 REAL dt 84 REAL Vbin(nbtr_bin) 85 85 86 86 ! output variable 87 87 REAL tr_seri(nbtr) -
LMDZ6/trunk/libf/phylmd/StratAer/strataer_local_var_mod.F90
r4767 r4950 51 51 52 52 !============= NUCLEATION VARS ============= 53 ! MOLECULAR ACCOMODATION OF H2SO4 (Raes and Van Dingen) 54 REAL,SAVE :: ALPH2SO4 ! H2SO4 accommodation coefficient [condensation/evaporation] 55 !$OMP THREADPRIVATE(ALPH2SO4) 56 53 57 ! flag to constraint nucleation rate in a lat/pres box 54 58 LOGICAL,SAVE :: flag_nuc_rate_box ! Nucleation rate limit or not to a lat/pres … … 64 68 INTEGER,SAVE :: flh2o ! ds stratemit : flh2o =0 (tr_seri), flh2o=1 (dq) 65 69 !$OMP THREADPRIVATE(flh2o) 66 ! REAL,ALLOCATABLE,SAVE :: d_q_emiss(:,:)67 ! !$OMP THREADPRIVATE(d_q_emiss)68 70 69 71 REAL,ALLOCATABLE,SAVE :: budg_emi(:,:) !DIMENSION(klon,n) … … 144 146 !$OMP THREADPRIVATE(day_emit_roc) 145 147 148 REAL,ALLOCATABLE,SAVE :: RRSI(:) ! radius [cm] for each aerosol size 149 REAL,ALLOCATABLE,SAVE :: Vbin(:) ! volume [m3] for each aerosol size 150 !$OMP THREADPRIVATE(RRSI, Vbin) 146 151 REAL,SAVE :: dlat, dlon ! delta latitude and d longitude of grid in degree 147 152 !$OMP THREADPRIVATE(dlat, dlon) … … 153 158 USE print_control_mod, ONLY : lunout 154 159 USE mod_phys_lmdz_para, ONLY : is_master 155 USE infotrac_phy, ONLY: id_OCS_strat,id_SO2_strat,id_H2SO4_strat,nbtr_sulgas 160 USE infotrac_phy, ONLY: id_OCS_strat,id_SO2_strat,id_H2SO4_strat,nbtr_sulgas,nbtr_bin 161 USE phys_local_var_mod, ONLY : mdw 162 USE aerophys, ONLY: mdwmin, V_rat 163 USE YOMCST , ONLY : RPI 164 165 INTEGER :: it 156 166 157 167 WRITE(lunout,*) 'IN STRATAER_LOCAL_VAR INIT WELCOME!' … … 185 195 186 196 ! nuc init 197 ALPH2SO4 = 0.1 187 198 flag_nuc_rate_box = .FALSE. 188 199 nuclat_min=0 ; nuclat_max=0 … … 238 249 ENDIF ! if master 239 250 251 !--initialising dry diameters to geometrically spaced mass/volume (see Jacobson 1994) 252 mdw(1)=mdwmin 253 IF (V_rat.LT.1.62) THEN ! compensate for dip in second bin for lower volume ratio 254 mdw(2)=mdw(1)*2.**(1./3.) 255 DO it=3, nbtr_bin 256 mdw(it)=mdw(it-1)*V_rat**(1./3.) 257 ENDDO 258 ELSE 259 DO it=2, nbtr_bin 260 mdw(it)=mdw(it-1)*V_rat**(1./3.) 261 ENDDO 262 ENDIF 263 IF (is_master) WRITE(lunout,*) 'init mdw=', mdw 264 265 ! compute particle radius RRSI [cm] and volume Vbin [m3] from diameter mdw [m] 266 ALLOCATE(RRSI(nbtr_bin), Vbin(nbtr_bin)) 267 268 DO it=1,nbtr_bin 269 ! [cm] 270 RRSI(it)=mdw(it)/2.*100. 271 ! [m3] 272 Vbin(it)=4.0*RPI*((mdw(it)/2.)**3)/3.0 273 ENDDO 274 275 IF (is_master) THEN 276 WRITE(lunout,*) 'init RRSI=', RRSI 277 WRITE(lunout,*) 'init Vbin=', Vbin 278 ENDIF 279 240 280 WRITE(lunout,*) 'IN STRATAER INIT END' 241 281 -
LMDZ6/trunk/libf/phylmd/StratAer/strataer_nuc_mod.F90
r4601 r4950 13 13 USE print_control_mod, ONLY : lunout 14 14 USE mod_phys_lmdz_para, ONLY : is_master 15 USE strataer_local_var_mod, ONLY: flag_nuc_rate_box,nuclat_min,nuclat_max,nucpres_min,nucpres_max 15 USE strataer_local_var_mod, ONLY: ALPH2SO4,flag_nuc_rate_box,nuclat_min,nuclat_max, & 16 nucpres_min,nucpres_max 16 17 17 18 !Config Key = flag_nuc_rate_box … … 30 31 CALL getin_p('nucpres_max',nucpres_max) 31 32 33 ! Read argument H2SO4 accommodation coefficient [condensation/evaporation] 34 CALL getin_p('alph2so4',ALPH2SO4) 35 32 36 !============= Print params ============= 33 37 IF (is_master) THEN 38 WRITE(lunout,*) 'IN STRATAER_NUC : ALPH2SO4 = ',alph2so4 34 39 WRITE(lunout,*) 'IN STRATAER_NUC : flag_nuc_rate_box = ',flag_nuc_rate_box 35 40 IF (flag_nuc_rate_box) THEN -
LMDZ6/trunk/libf/phylmd/StratAer/sulfate_aer_mod.F90
r4750 r4950 7 7 8 8 !******************************************************************* 9 SUBROUTINE STRACOMP_BIN(sh,t_seri,pplay) 10 ! 11 ! Aerosol H2SO4 weight fraction as a function of PH2O and temperature 12 ! INPUT: 13 ! sh: VMR of H2O 14 ! t_seri: temperature (K) 15 ! pplay: middle layer pression (Pa) 16 ! 17 ! OUTPUT: 18 ! R2SO4: aerosol H2SO4 weight fraction (percent) 9 SUBROUTINE STRACOMP_KELVIN(sh,t_seri,pplay) 10 ! 11 ! Aerosol H2SO4 weight fraction as a function of PH2O and temperature 12 ! INPUT: 13 ! sh: MMR of H2O 14 ! t_seri: temperature (K) 15 ! pplay: middle layer pression (Pa) 16 ! 17 ! Modified in modules: 18 ! R2SO4: aerosol H2SO4 weight fraction (percent) 19 ! R2SO4B: aerosol H2SO4 weight fraction (percent) for each aerosol bin 20 ! DENSO4: aerosol density (gr/cm3) 21 ! DENSO4B: aerosol density (gr/cm3)for each aerosol bin 22 ! f_r_wet: factor for converting dry to wet radius 23 ! assuming 'flat surface' composition (does not depend on aerosol size) 24 ! f_r_wetB: factor for converting dry to wet radius 25 ! assuming 'curved surface' composition (depends on aerosol size) 19 26 20 USE dimphy, ONLY : klon,klev ! nb of longitude and altitude bands 21 USE aerophys 22 USE phys_local_var_mod, ONLY: R2SO4 27 USE dimphy, ONLY : klon,klev ! nb of longitude and altitude bands 28 USE infotrac_phy, ONLY : nbtr_bin 29 USE aerophys 30 USE phys_local_var_mod, ONLY: R2SO4, R2SO4B, DENSO4, DENSO4B, f_r_wet, f_r_wetB 31 USE strataer_local_var_mod, ONLY: RRSI 32 ! WARNING: in phys_local_var_mod R2SO4B, DENSO4B, f_r_wetB (klon,klev,nbtr_bin) 33 ! and dens_aer_dry must be declared somewhere 23 34 24 IMPLICIT NONE35 IMPLICIT NONE 25 36 26 REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature 27 REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression in the middle of each layer (Pa) 28 REAL,DIMENSION(klon,klev),INTENT(IN) :: sh ! specific humidity 37 REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature 38 REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression in the middle of each layer (Pa) 39 REAL,DIMENSION(klon,klev),INTENT(IN) :: sh ! specific humidity (kg h2o/kg air) 40 41 ! local variables 42 integer :: ilon,ilev,ik 43 real, parameter :: rath2oair = mAIRmol/mH2Omol 44 real, parameter :: third = 1./3. 45 real :: pph2ogas(klon,klev) 46 real :: temp, wpp, xa, surtens, mvh2o, radwet, fkelvin, pph2okel, r2so4ik, denso4ik 47 !---------------------------------------- 48 49 ! gas-phase h2o partial pressure (Pa) 50 ! vmr=sh*rath2oair 51 pph2ogas(:,:) = pplay(:,:)*sh(:,:)*rath2oair 29 52 30 REAL ks(7) 31 REAL t,qh2o,ptot,pw 32 REAL a,b,c,det 33 REAL xsb,msb 53 DO ilon=1,klon 54 DO ilev=1,klev 55 56 temp = max(t_seri(ilon,ilev),190.) 57 temp = min(temp,300.) 58 59 ! *** H2SO4-H2O flat surface *** 60 !! equilibrium H2O pressure over pure flat liquid water (Pa) 61 !! pflath2o=psh2o(temp) 62 ! h2so4 weight percent(%) = f(P_h2o(Pa),T) 63 R2SO4(ilon,ilev)=wph2so4(pph2ogas(ilon,ilev),temp) 64 ! h2so4 mass fraction (0<wpp<1) 65 wpp=R2SO4(ilon,ilev)*1.e-2 66 ! mole fraction 67 xa=18.*wpp/(18.*wpp+98.*(1.-wpp)) 68 69 ! CHECK:compare h2so4 sat/ pressure (see Marti et al., 97 & reef. therein) 70 ! R2SO4(ilon,ilev)=70. temp=298.15 71 ! equilibrium h2so4 number density over H2SO4/H2O solution (molec/cm3) 72 ! include conversion from molec/cm3 to Pa 73 ! ph2so4=solh2so4(temp,xa)*(1.38065e-16*temp)/10. 74 ! print*,' ph2so4=',ph2so4,temp,R2SO4(ilon,ilev) 75 ! good match with Martin, et Ayers, not with Gmitro (the famous 0.086) 76 77 ! surface tension (mN/m=1.e-3.kg/s2) = f(T,h2so4 mole fraction) 78 surtens=surftension(temp,xa) 79 ! molar volume of pure h2o (cm3.mol-1 =1.e-6.m3.mol-1) 80 mvh2o= rmvh2o(temp) 81 ! aerosol density (gr/cm3) = f(T,h2so4 mass fraction) 82 DENSO4(ilon,ilev)=density(temp,wpp) 83 ! ->x1000., to have it in kg/m3 84 ! factor for converting dry to wet radius 85 f_r_wet(ilon,ilev) = (dens_aer_dry/(DENSO4(ilon,ilev)*1.e3)/ & 86 & (R2SO4(ilon,ilev)*1.e-2))**third 87 ! *** End of H2SO4-H2O flat surface *** 88 89 90 ! Loop on bin radius (RRSI in cm) 91 DO IK=1,nbtr_bin 92 93 ! *** H2SO4-H2O curved surface - Kelvin effect factor *** 94 ! wet radius (m) (RRSI(IK) in [cm]) 95 if (f_r_wetB(ilon,ilev,IK) .gt. 1.0) then 96 radwet = 1.e-2*RRSI(IK)*f_r_wetB(ilon,ilev,IK) 97 else 98 ! H2SO4-H2O flat surface, only on the first timestep 99 radwet = 1.e-2*RRSI(IK)*f_r_wet(ilon,ilev) 100 endif 101 ! Kelvin factor: 102 ! surface tension (mN/m=1.e-3.kg/s2) 103 ! molar volume of pure h2o (cm3.mol-1 =1.e-6.m3.mol-1) 104 fkelvin=exp( 2.*1.e-3*surtens*1.e-6*mvh2o/ (radwet*rgas*temp) ) 105 ! equilibrium: pph2o(gas) = pph2o(liq) = pph2o(liq_flat) * fkelvin 106 ! equilibrium: pph2o(liq_flat) = pph2o(gas) / fkelvin 107 ! h2o liquid partial pressure before Kelvin effect (Pa) 108 pph2okel = pph2ogas(ilon,ilev) / fkelvin 109 ! h2so4 weight percent(%) = f(P_h2o(Pa),temp) 110 r2so4ik=wph2so4(pph2okel,temp) 111 ! h2so4 mass fraction (0<wpp<1) 112 wpp=r2so4ik*1.e-2 113 ! mole fraction 114 xa=18.*wpp/(18.*wpp+98.*(1.-wpp)) 115 ! aerosol density (gr/cm3) = f(T,h2so4 mass fraction) 116 denso4ik=density(temp,wpp) 117 ! 118 ! recalculate Kelvin factor with surface tension and radwet 119 ! with new R2SO4B and DENSO4B 120 surtens=surftension(temp,xa) 121 ! wet radius (m) 122 radwet = 1.e-2*RRSI(IK)*(dens_aer_dry/(denso4ik*1.e3)/ & 123 & (r2so4ik*1.e-2))**third 124 fkelvin=exp( 2.*1.e-3*surtens*1.e-6*mvh2o / (radwet*rgas*temp) ) 125 pph2okel=pph2ogas(ilon,ilev) / fkelvin 126 ! h2so4 weight percent(%) = f(P_h2o(Pa),temp) 127 R2SO4B(ilon,ilev,IK)=wph2so4(pph2okel,temp) 128 ! h2so4 mass fraction (0<wpp<1) 129 wpp=R2SO4B(ilon,ilev,IK)*1.e-2 130 xa=18.*wpp/(18.*wpp+98.*(1.-wpp)) 131 ! aerosol density (gr/cm3) = f(T,h2so4 mass fraction) 132 DENSO4B(ilon,ilev,IK)=density(temp,wpp) 133 ! factor for converting dry to wet radius 134 f_r_wetB(ilon,ilev,IK) = (dens_aer_dry/(DENSO4B(ilon,ilev,IK)*1.e3)/ & 135 & (R2SO4B(ilon,ilev,IK)*1.e-2))**third 136 ! 137 ! print*,'R,Rwet(m),kelvin,h2so4(%),ro=',RRSI(ik),radwet,fkelvin, & 138 ! & R2SO4B(ilon,ilev,IK),DENSO4B(ilon,ilev,IK) 139 ! print*,' equil.h2so4(molec/cm3), & 140 ! & sigma',solh2so4(temp,xa),surftension(temp,xa) 141 142 ENDDO 143 144 ENDDO 145 ENDDO 146 147 RETURN 34 148 35 INTEGER ilon,ilev 36 DATA ks/-21.661,2724.2,51.81,-15732.0,47.004,-6969.0,-4.6183/ 37 38 !******************************************************************* 39 !*** liquid aerosols process 40 !******************************************************************* 41 ! BINARIES LIQUID AEROROLS: 42 43 DO ilon=1,klon 44 DO ilev=1,klev 45 46 t = max(t_seri(ilon,ilev),185.) 47 qh2o=sh(ilon,ilev)/18.*28.9 48 ptot=pplay(ilon,ilev)/100. 49 pw = qh2o*ptot/1013.0 50 pw = min(pw,2.e-3/1013.) 51 pw = max(pw,2.e-5/1013.) 52 53 !******************************************************************* 54 !*** binaries aerosols h2so4/h2o 55 !******************************************************************* 56 a = ks(3) + ks(4)/t 57 b = ks(1) + ks(2)/t 58 c = ks(5) + ks(6)/t + ks(7)*log(t) - log(pw) 59 60 det = b**2 - 4.*a*c 61 62 IF (det > 0.) THEN 63 xsb = (-b - sqrt(det))/(2.*a) 64 msb = 55.51*xsb/(1.0 - xsb) 65 ELSE 66 msb = 0. 67 ENDIF 68 R2SO4(ilon,ilev) = 100*msb*0.098076/(1.0 + msb*0.098076) 69 70 ! H2SO4 min dilution: 0.5% 71 R2SO4(ilon,ilev) = max( R2SO4(ilon,ilev), 0.005 ) 72 ENDDO 73 ENDDO 74 100 RETURN 75 76 END SUBROUTINE STRACOMP_BIN 77 149 END SUBROUTINE STRACOMP_KELVIN 78 150 !******************************************************************** 79 151 SUBROUTINE STRACOMP(sh,t_seri,pplay) … … 544 616 545 617 END SUBROUTINE 546 547 !****************************************************************548 SUBROUTINE DENH2SA_TABA(t_seri)549 550 ! AERSOL DENSITY AS A FUNCTION OF H2SO4 WEIGHT PERCENT AND T551 ! from Tabazadeh et al. (1994) abaques552 ! ---------------------------------------------553 554 !555 ! INPUT:556 ! R2SO4: aerosol H2SO4 weight fraction (percent)557 ! t_seri: temperature (K)558 ! klon: number of latitude bands in the model domain559 ! klev: number of altitude bands in the model domain560 ! for IFS: perhaps add another dimension for longitude561 !562 ! OUTPUT:563 ! DENSO4: aerosol mass density (gr/cm3 = aerosol mass/aerosol volume)564 !565 USE dimphy, ONLY : klon,klev566 USE phys_local_var_mod, ONLY: R2SO4, DENSO4567 568 IMPLICIT NONE569 570 REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature571 572 INTEGER i,j573 574 !----------------------------------------------------------------------575 ! ... Local variables576 !----------------------------------------------------------------------577 real, parameter :: a9 = -268.2616e4, a10 = 576.4288e3578 579 real :: a0, a1, a2, a3, a4, a5, a6, a7 ,a8580 real :: c1, c2, c3, c4, w581 582 583 ! Loop on model domain (2 dimension for UPMC model; 3 for IFS)584 DO i=1,klon585 DO j=1,klev586 !----------------------------------------------------------------------587 ! ... Temperature variables588 !----------------------------------------------------------------------589 c1 = t_seri(I,J)- 273.15590 c2 = c1**2591 c3 = c1*c2592 c4 = c1*c3593 !----------------------------------------------------------------------594 ! Polynomial Coefficients595 !----------------------------------------------------------------------596 a0 = 999.8426 + 334.5402e-4*c1 - 569.1304e-5*c2597 a1 = 547.2659 - 530.0445e-2*c1 + 118.7671e-4*c2 + 599.0008e-6*c3598 a2 = 526.295e1 + 372.0445e-1*c1 + 120.1909e-3*c2 - 414.8594e-5*c3 + 119.7973e-7*c4599 a3 = -621.3958e2 - 287.7670*c1 - 406.4638e-3*c2 + 111.9488e-4*c3 + 360.7768e-7*c4600 a4 = 409.0293e3 + 127.0854e1*c1 + 326.9710e-3*c2 - 137.7435e-4*c3 - 263.3585e-7*c4601 a5 = -159.6989e4 - 306.2836e1*c1 + 136.6499e-3*c2 + 637.3031e-5*c3602 a6 = 385.7411e4 + 408.3717e1*c1 - 192.7785e-3*c2603 a7 = -580.8064e4 - 284.4401e1*c1604 a8 = 530.1976e4 + 809.1053*c1605 !----------------------------------------------------------------------606 ! ... Summation607 !----------------------------------------------------------------------608 ! w : H2SO4 Weight fraction609 w=r2SO4(i,j)*0.01610 DENSO4(i,j) = 0.001*(a0 + w*(a1 + w*(a2 + w*(a3 + w*(a4 + &611 w*(a5 + w*(a6 + w*(a7 + w*(a8 + w*(a9 + w*a10))))))))))612 DENSO4(i,j) = max (0.0, DENSO4(i,j) )613 614 ENDDO615 ENDDO616 617 END SUBROUTINE DENH2SA_TABA618 618 619 619 !**************************************************************** … … 764 764 RETURN 765 765 END SUBROUTINE 766 766 !******************************************************************** 767 !----------------------------------------------------------------------- 768 real function psh2so4(T) result(psh2so4_out) 769 ! equilibrium H2SO4 pressure over pure H2SO4 solution (Pa) 770 ! 771 !---->Ayers et.al. (1980), GRL (7) pp 433-436 772 ! plus corrections for lower temperatures by Kulmala and Laaksonen (1990) 773 ! and Noppel et al. (1990) 774 775 implicit none 776 real, intent(in) :: T 777 real, parameter :: & 778 & b1=1.01325e5, & 779 & b2=11.5, & 780 & b3=1.0156e4, & 781 & b4=0.38/545., & 782 & tref=360.15 783 784 ! saturation vapor pressure ( N/m2 = Pa = kg/(m.s2) ) 785 psh2so4_out=b1*exp( -b2 +b3*( 1./tref-1./T & 786 & +b4*(1.+log(tref/T)-tref/T) ) ) 787 788 return 789 end function psh2so4 790 !----------------------------------------------------------------------- 791 real function ndsh2so4(T) result(ndsh2so4_out) 792 ! equilibrium H2SO4 number density over pure H2SO4 (molec/cm3) 793 794 implicit none 795 real, intent(in) :: T 796 real :: presat 797 798 ! Boltzmann constant ( 1.38065e-23 J/K = m2⋅kg/(s2⋅K) ) 799 ! akb idem in cm2⋅g/(s2⋅K) 800 real, parameter :: akb=1.38065e-16 801 802 ! pure h2so4 saturation vapor pressure (Pa) 803 presat=psh2so4(T) 804 ! saturation number density (1/cm3) - (molec/cm3) 805 ndsh2so4_out=presat*10./(akb*T) 806 807 return 808 end function ndsh2so4 809 !----------------------------------------------------------------------- 810 real function psh2o(T) result(psh2o_out) 811 ! equilibrium H2O pressure over pure liquid water (Pa) 812 ! 813 implicit none 814 real, intent(in) :: T 815 816 if(T.gt.229.) then 817 ! Preining et al., 1981 (from Kulmala et al., 1998) 818 ! saturation vapor pressure (N/m2 = 1 Pa = 1 kg/(m·s2)) 819 psh2o_out=exp( 77.34491296 -7235.424651/T & 820 & -8.2*log(T) + 5.7133e-3*T ) 821 else 822 ! Tabazadeh et al., 1997, parameterization for 185<T<260 823 ! saturation water vapor partial pressure (mb = hPa =1.E2 kg/(m·s2)) 824 ! or from Clegg and Brimblecombe , J. Chem. Eng., p43, 1995. 825 ; 826 psh2o_out=18.452406985 -3505.1578807/T & 827 & -330918.55082/(T*T) & 828 & +12725068.262/(T*T*T) 829 ! in Pa 830 psh2o_out=100.*exp(psh2o_out) 831 end if 832 ! print*,psh2o_out 833 834 return 835 end function psh2o 836 !----------------------------------------------------------------------- 837 real function density(T,so4mfrac) result(density_out) 838 ! calculation of particle density (gr/cm3) 839 840 ! requires Temperature (T) and acid mass fraction (so4mfrac) 841 !---->Vehkamaeki et al. (2002) 842 843 implicit none 844 real, intent(in) :: T, so4mfrac 845 real, parameter :: & 846 & a1= 0.7681724,& 847 & a2= 2.184714, & 848 & a3= 7.163002, & 849 & a4=-44.31447, & 850 & a5= 88.74606, & 851 & a6=-75.73729, & 852 & a7= 23.43228 853 real, parameter :: & 854 & b1= 1.808225e-3, & 855 & b2=-9.294656e-3, & 856 & b3=-3.742148e-2, & 857 & b4= 2.565321e-1, & 858 & b5=-5.362872e-1, & 859 & b6= 4.857736e-1, & 860 & b7=-1.629592e-1 861 real, parameter :: & 862 & c1=-3.478524e-6, & 863 & c2= 1.335867e-5, & 864 & c3= 5.195706e-5, & 865 & c4=-3.717636e-4, & 866 & c5= 7.990811e-4, & 867 & c6=-7.458060e-4, & 868 & c7= 2.581390e-4 869 real :: a,b,c,so4m2,so4m3,so4m4,so4m5,so4m6 870 871 so4m2=so4mfrac*so4mfrac 872 so4m3=so4mfrac*so4m2 873 so4m4=so4mfrac*so4m3 874 so4m5=so4mfrac*so4m4 875 so4m6=so4mfrac*so4m5 876 877 a=+a1+a2*so4mfrac+a3*so4m2+a4*so4m3 & 878 & +a5*so4m4+a6*so4m5+a7*so4m6 879 b=+b1+b2*so4mfrac+b3*so4m2+b4*so4m3 & 880 & +b5*so4m4+b6*so4m5+b7*so4m6 881 c=+c1+c2*so4mfrac+c3*so4m2+c4*so4m3 & 882 & +c5*so4m4+c6*so4m5+c7*so4m6 883 density_out=(a+b*T+c*T*T) ! units are gm/cm**3 884 885 return 886 end function density 887 !----------------------------------------------------------------------- 888 real function surftension(T,so4frac) result(surftension_out) 889 ! calculation of surface tension (mN/meter) 890 ! requires Temperature (T) and acid mole fraction (so4frac) 891 !---->Vehkamaeki et al. (2002) 892 893 implicit none 894 real,intent(in) :: T, so4frac 895 real :: a,b,so4mfrac,so4m2,so4m3,so4m4,so4m5,so4sig 896 real, parameter :: & 897 & a1= 0.11864, & 898 & a2=-0.11651, & 899 & a3= 0.76852, & 900 & a4=-2.40909, & 901 & a5= 2.95434, & 902 & a6=-1.25852 903 real, parameter :: & 904 & b1=-1.5709e-4, & 905 & b2= 4.0102e-4, & 906 & b3=-2.3995e-3, & 907 & b4= 7.611235e-3, & 908 & b5=-9.37386e-3, & 909 & b6= 3.89722e-3 910 real, parameter :: convfac=1.e3 ! convert from newton/m to dyne/cm 911 real, parameter :: Mw=18.01528, Ma=98.079 912 913 ! so4 mass fraction 914 so4mfrac=Ma*so4frac/( Ma*so4frac+Mw*(1.-so4frac) ) 915 so4m2=so4mfrac*so4mfrac 916 so4m3=so4mfrac*so4m2 917 so4m4=so4mfrac*so4m3 918 so4m5=so4mfrac*so4m4 919 920 a=+a1+a2*so4mfrac+a3*so4m2+a4*so4m3+a5*so4m4+a6*so4m5 921 b=+b1+b2*so4mfrac+b3*so4m2+b4*so4m3+b5*so4m4+b6*so4m5 922 so4sig=a+b*T 923 surftension_out=so4sig*convfac 924 925 return 926 end function surftension 927 !----------------------------------------------------------------------- 928 real function wph2so4(pph2o,T) result(wph2so4_out) 929 ! Calculates the equilibrium composition of h2so4 aerosols 930 ! as a function of temperature and H2O pressure, using 931 ! the parameterization of Tabazadeh et al., GRL, p1931, 1997. 932 ! 933 ! Parameters 934 ! 935 ! input: 936 ! T.....temperature (K) 937 ! pph2o..... amhbiant 2o pressure (Pa) 938 ! 939 ! output: 940 ! wph2so4......sulfuric acid composition (weight percent wt % h2so4) 941 ! = h2so4 mass fraction*100. 942 ! 943 implicit none 944 real, intent(in) :: pph2o, T 945 946 real :: aw, rh, y1, y2, sulfmolal 947 948 ! psh2o(T): equilibrium H2O pressure over pure liquid water (Pa) 949 ! relative humidity 950 rh=pph2o/psh2o(T) 951 ! water activity 952 ! aw=min( 0.999,max(1.e-3,rh) ) 953 aw=min( 0.999999999,max(1.e-8,rh) ) 954 955 ! composition 956 ! calculation of h2so4 molality 957 if(aw .le. 0.05 .and. aw .gt. 0.) then 958 y1=12.372089320*aw**(-0.16125516114) & 959 & -30.490657554*aw -2.1133114241 960 y2=13.455394705*aw**(-0.19213122550) & 961 & -34.285174607*aw -1.7620073078 962 else if(aw .le. 0.85 .and. aw .gt. 0.05) then 963 y1=11.820654354*aw**(-0.20786404244) & 964 & -4.8073063730*aw -5.1727540348 965 y2=12.891938068*aw**(-0.23233847708) & 966 & -6.4261237757*aw -4.9005471319 967 else 968 y1=-180.06541028*aw**(-0.38601102592) & 969 & -93.317846778*aw +273.88132245 970 y2=-176.95814097*aw**(-0.36257048154) & 971 & -90.469744201*aw +267.45509988 972 end if 973 ! h2so4 molality (m=moles of h2so4 (solute)/ kg of h2o(solvent)) 974 sulfmolal = y1+((T-190.)*(y2-y1)/70.) 975 976 ! for a solution containing mh2so4 and mh2o: 977 ! sulfmolal = (mh2so4(gr)/h2so4_molar_mass(gr/mole)) / (mh2o(gr)*1.e-3) 978 ! mh2o=1.e3*(mh2so4/Mh2so4)/sulfmolal=1.e3*mh2so4/(Mh2so4*sulfmolal) 979 ! h2so4_mass_fraction = mfh2so4 = mh2so4/(mh2o + mh2so4) 980 ! mh2o=mh2so4*(1-mfh2so4)/mfh2so4 981 ! combining the 2 equations 982 ! 1.e3*mh2so4/(Mh2so4*sulfmolal) = mh2so4*(1-mfh2so4)/mfh2so4 983 ! 1.e3/(Mh2so4*sulfmolal) = (1-mfh2so4)/mfh2so4 984 ! 1000*mfh2so4 = (1-mfh2so4)*Mh2so4*sulfmolal 985 ! mfh2so4*(1000.+Mh2so4*sulfmolal) = Mh2so4*sulfmolal 986 ! mfh2so4 = Mh2so4*sulfmolal / (1000.+Mh2so4*sulfmolal) 987 ! wph2so4 (% mass fraction)= 100.*Mh2so4*sulfmolal / (1000.+Mh2so4*sulfmolal) 988 ! recall activity of i = a_i = P_i/P_pure_i and 989 ! activity coefficient of i = gamma_i = a_i/X_i (X_i: mole fraction of i) 990 ! so P_i = gamma_i*X_i*P_pure_i 991 ! if ideal solution, gamma_i=1, P_i = X_i*P_pure_i 992 993 ! h2so4 weight precent 994 wph2so4_out = 9800.*sulfmolal/(98.*sulfmolal+1000.) 995 ! print*,rh,pph2o,psh2o(T),vpice(T) 996 ! print*,T,aw,sulfmolal,wph2so4_out 997 wph2so4_out = max(wph2so4_out,15.) 998 wph2so4_out = min(wph2so4_out,99.999) 999 1000 return 1001 end function wph2so4 1002 !----------------------------------------------------------------------- 1003 real function solh2so4(T,xa) result(solh2so4_out) 1004 ! equilibrium h2so4 number density over H2SO4/H2O solution (molec/cm3) 1005 1006 implicit none 1007 real, intent(in) :: T, xa ! T(K) xa(H2SO4 mass fraction) 1008 1009 real :: xw, a12,b12, cacta, presat 1010 1011 xw=1.0-xa 1012 1013 ! pure h2so4 saturation number density (molec/cm3) 1014 presat=ndsh2so4(T) 1015 ! compute activity of acid 1016 a12=5.672E3 -4.074E6/T +4.421E8/(T*T) 1017 b12=1./0.527 1018 cacta=10.**(a12*xw*xw/(xw+b12*xa)**2/T) 1019 ! h2so4 saturation number density over H2SO4/H2O solution (molec/cm3) 1020 solh2so4_out=cacta*xa*presat 1021 1022 return 1023 end function solh2so4 1024 !----------------------------------------------------------------------- 1025 real function rpmvh2so4(T,ws) result(rpmvh2so4_out) 1026 ! partial molar volume of h2so4 in h2so4/h2o solution (cm3/mole) 1027 1028 implicit none 1029 real, intent(in) :: T, ws 1030 real, dimension(22),parameter :: x=(/ & 1031 & 2.393284E-02,-4.359335E-05,7.961181E-08,0.0,-0.198716351, & 1032 & 1.39564574E-03,-2.020633E-06,0.51684706,-3.0539E-03,4.505475E-06, & 1033 & -0.30119511,1.840408E-03,-2.7221253742E-06,-0.11331674116, & 1034 & 8.47763E-04,-1.22336185E-06,0.3455282,-2.2111E-03,3.503768245E-06, & 1035 & -0.2315332,1.60074E-03,-2.5827835E-06/) 1036 1037 real :: w 1038 1039 w=ws*0.01 1040 rpmvh2so4_out=x(5)+x(6)*T+x(7)*T*T+(x(8)+x(9)*T+x(10)*T*T)*w & 1041 +(x(11)+x(12)*T+x(13)*T*T)*w*w 1042 ! h2so4 partial molar volume in h2so4/h2o solution (cm3/mole) 1043 rpmvh2so4_out=rpmvh2so4_out*1000. 1044 1045 return 1046 end function rpmvh2so4 1047 !----------------------------------------------------------------------- 1048 real function rmvh2o(T) result(rmvh2o_out) 1049 ! molar volume of pure h2o (cm3/mole) 1050 1051 implicit none 1052 real, intent(in) :: T 1053 real, parameter :: x1=2.393284E-02,x2=-4.359335E-05,x3=7.961181E-08 1054 1055 ! 1000: L/mole -> cm3/mole 1056 ! pure h2o molar volume (cm3/mole) 1057 rmvh2o_out=(x1+x2*T+x3*T*T)*1000. 1058 1059 return 1060 end function rmvh2o 1061 ! 767 1062 END MODULE sulfate_aer_mod -
LMDZ6/trunk/libf/phylmd/StratAer/traccoag_mod.F90
r4769 r4950 9 9 presnivs, xlat, xlon, pphis, pphi, & 10 10 t_seri, pplay, paprs, sh, rh, tr_seri) 11 11 12 12 USE phys_local_var_mod, ONLY: mdw, R2SO4, DENSO4, f_r_wet, surf_PM25_sulf, & 13 & budg_emi_ocs, budg_emi_so2, budg_emi_h2so4, budg_emi_part 14 13 & budg_emi_ocs, budg_emi_so2, budg_emi_h2so4, budg_emi_part, & 14 & R2SO4B, DENSO4B, f_r_wetB 15 15 16 USE dimphy 16 17 USE infotrac_phy, ONLY : nbtr_bin, nbtr_sulgas, nbtr, id_SO2_strat … … 88 89 ENDIF 89 90 91 ! radius [m] 90 92 DO it=1, nbtr_bin 91 93 r_bin(it)=mdw(it)/2. … … 117 119 118 120 IF(flag_new_strat_compo) THEN 119 IF(debutphy) WRITE(lunout,*) 'traccoag: USE STRAT COMPO from Tabazadeh 1994', flag_new_strat_compo 120 ! STRACOMP (H2O, P, t_seri -> aerosol composition (R2SO4)) : binary routine (from reprobus) 121 ! H2SO4 mass fraction in aerosol (%) from Tabazadeh et al. (1994). 122 CALL stracomp_bin(sh,t_seri,pplay) 123 124 ! aerosol density (gr/cm3) - from Tabazadeh 125 CALL denh2sa_taba(t_seri) 121 IF(debutphy) WRITE(lunout,*) 'traccoag: COMPO/DENSITY (Tabazadeh 97) + H2O kelvin effect', flag_new_strat_compo 122 ! STRACOMP (H2O, P, t_seri, R -> R2SO4 + Kelvin effect) : Taba97, Socol, etc... 123 CALL stracomp_kelvin(sh,t_seri,pplay) 126 124 ELSE 127 IF(debutphy) WRITE(lunout,*) 'traccoag: USE STRATCOMPO from Bekki 2D model', flag_new_strat_compo125 IF(debutphy) WRITE(lunout,*) 'traccoag: COMPO from Bekki 2D model', flag_new_strat_compo 128 126 ! STRACOMP (H2O, P, t_seri -> aerosol composition (R2SO4)) 129 127 ! H2SO4 mass fraction in aerosol (%) … … 132 130 ! aerosol density (gr/cm3) 133 131 CALL denh2sa(t_seri) 132 133 ! compute factor for converting dry to wet radius (for every grid box) 134 f_r_wet(:,:) = (dens_aer_dry/(DENSO4(:,:)*1000.)/(R2SO4(:,:)/100.))**(1./3.) 134 135 ENDIF 135 136 136 ! compute factor for converting dry to wet radius (for every grid box)137 f_r_wet(:,:) = (dens_aer_dry/(DENSO4(:,:)*1000.)/(R2SO4(:,:)/100.))**(1./3.)138 139 137 !--calculate mass of air in every grid box 140 138 DO ilon=1, klon
Note: See TracChangeset
for help on using the changeset viewer.
