| 1 | MODULE sulfate_aer_mod |
|---|
| 2 | |
|---|
| 3 | ! microphysical routines based on UPMC aerosol model by Slimane Bekki |
|---|
| 4 | ! adapted for stratospheric sulfate aerosol in LMDZ by Christoph Kleinschmitt |
|---|
| 5 | |
|---|
| 6 | CONTAINS |
|---|
| 7 | |
|---|
| 8 | !******************************************************************* |
|---|
| 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) |
|---|
| 26 | |
|---|
| 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 |
|---|
| 34 | |
|---|
| 35 | IMPLICIT NONE |
|---|
| 36 | |
|---|
| 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 |
|---|
| 52 | |
|---|
| 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) > 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 | |
|---|
| 148 | |
|---|
| 149 | END SUBROUTINE STRACOMP_KELVIN |
|---|
| 150 | !******************************************************************** |
|---|
| 151 | SUBROUTINE STRACOMP(sh,t_seri,pplay) |
|---|
| 152 | |
|---|
| 153 | ! AEROSOL H2SO4 WEIGHT FRACTION AS A FUNCTION OF PH2O AND TEMPERATURE |
|---|
| 154 | ! ---------------------------------------------------------------- |
|---|
| 155 | ! INPUT: |
|---|
| 156 | ! H2O: VMR of H2O |
|---|
| 157 | ! t_seri: temperature (K) |
|---|
| 158 | ! PMB: pressure (mb) |
|---|
| 159 | ! klon: number of latitude bands in the model domain |
|---|
| 160 | ! klev: number of altitude bands in the model domain |
|---|
| 161 | ! for IFS: perhaps add another dimension for longitude |
|---|
| 162 | |
|---|
| 163 | ! OUTPUT: |
|---|
| 164 | ! R2SO4: aerosol H2SO4 weight fraction (percent) |
|---|
| 165 | |
|---|
| 166 | USE dimphy, ONLY: klon,klev |
|---|
| 167 | USE aerophys |
|---|
| 168 | USE phys_local_var_mod, ONLY: R2SO4 |
|---|
| 169 | |
|---|
| 170 | IMPLICIT NONE |
|---|
| 171 | |
|---|
| 172 | REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature |
|---|
| 173 | REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa) |
|---|
| 174 | REAL,DIMENSION(klon,klev),INTENT(IN) :: sh ! humidite specifique |
|---|
| 175 | |
|---|
| 176 | REAL PMB(klon,klev), H2O(klon,klev) |
|---|
| 177 | |
|---|
| 178 | ! working variables |
|---|
| 179 | INTEGER I,J,K |
|---|
| 180 | REAL TP, PH2O, VAL, A, B |
|---|
| 181 | ! local variables to be saved on exit |
|---|
| 182 | INTEGER INSTEP |
|---|
| 183 | INTEGER, PARAMETER :: N=16, M=28 |
|---|
| 184 | DATA INSTEP/0/ |
|---|
| 185 | REAL F(N,M) |
|---|
| 186 | REAL XC(N) |
|---|
| 187 | REAL YC(M) |
|---|
| 188 | REAL XC1, XC16, YC1, YC28 |
|---|
| 189 | |
|---|
| 190 | SAVE INSTEP,F,XC,YC,XC1,XC16,YC1,YC28 |
|---|
| 191 | !$OMP THREADPRIVATE(INSTEP,F,XC,YC,XC1,XC16,YC1,YC28) |
|---|
| 192 | |
|---|
| 193 | ! convert pplay (in Pa) to PMB (in mb) |
|---|
| 194 | PMB(:,:)=pplay(:,:)/100.0 |
|---|
| 195 | |
|---|
| 196 | ! convert specific humidity sh (in kg/kg) to VMR H2O |
|---|
| 197 | H2O(:,:)=sh(:,:)*mAIRmol/mH2Omol |
|---|
| 198 | |
|---|
| 199 | IF(INSTEP==0) THEN |
|---|
| 200 | |
|---|
| 201 | INSTEP=1 |
|---|
| 202 | XC(1)=0.01 |
|---|
| 203 | XC(2)=0.1 |
|---|
| 204 | XC(3)=0.5 |
|---|
| 205 | XC(4)=1.0 |
|---|
| 206 | XC(5)=1.5 |
|---|
| 207 | XC(6)=2.0 |
|---|
| 208 | XC(7)=3.0 |
|---|
| 209 | XC(8)=5.0 |
|---|
| 210 | XC(9)=6.0 |
|---|
| 211 | XC(10)=8.0 |
|---|
| 212 | XC(11)=10.0 |
|---|
| 213 | XC(12)=12.0 |
|---|
| 214 | XC(13)=15.0 |
|---|
| 215 | XC(14)=20.0 |
|---|
| 216 | XC(15)=30.0 |
|---|
| 217 | XC(16)=100.0 |
|---|
| 218 | |
|---|
| 219 | YC(1)=175.0 |
|---|
| 220 | DO I=2,28 |
|---|
| 221 | YC(I)=YC(I-1)+5.0 |
|---|
| 222 | ENDDO |
|---|
| 223 | |
|---|
| 224 | ! CONVERSION mb IN 1.0E-4mB |
|---|
| 225 | DO I=1,16 |
|---|
| 226 | XC(I)=XC(I)*1.0E-4 |
|---|
| 227 | ENDDO |
|---|
| 228 | |
|---|
| 229 | XC1=XC(1)+1.E-10 |
|---|
| 230 | XC16=XC(16)-1.E-8 |
|---|
| 231 | YC1=YC(1)+1.E-5 |
|---|
| 232 | YC28=YC(28)-1.E-5 |
|---|
| 233 | |
|---|
| 234 | F(6,4)=43.45 |
|---|
| 235 | F(6,5)=53.96 |
|---|
| 236 | F(6,6)=60.62 |
|---|
| 237 | F(6,7)=65.57 |
|---|
| 238 | F(6,8)=69.42 |
|---|
| 239 | F(6,9)=72.56 |
|---|
| 240 | F(6,10)=75.17 |
|---|
| 241 | F(6,11)=77.38 |
|---|
| 242 | F(6,12)=79.3 |
|---|
| 243 | F(6,13)=80.99 |
|---|
| 244 | F(6,14)=82.5 |
|---|
| 245 | F(6,15)=83.92 |
|---|
| 246 | F(6,16)=85.32 |
|---|
| 247 | F(6,17)=86.79 |
|---|
| 248 | F(6,18)=88.32 |
|---|
| 249 | |
|---|
| 250 | ! ADD FACTOR BECAUSE THE SLOP IS TOO IMPORTANT |
|---|
| 251 | ! NOT FOR THIS ONE BUT THE REST |
|---|
| 252 | ! LOG DOESN'T WORK |
|---|
| 253 | A=(F(6,5)-F(6,4))/( (YC(5)-YC(4))*2.0) |
|---|
| 254 | B=-A*YC(4) + F(6,4) |
|---|
| 255 | F(6,1)=A*YC(1) + B |
|---|
| 256 | F(6,2)=A*YC(2) + B |
|---|
| 257 | F(6,3)=A*YC(3) + B |
|---|
| 258 | |
|---|
| 259 | F(7,4)=37.02 |
|---|
| 260 | F(7,5)=49.46 |
|---|
| 261 | F(7,6)=57.51 |
|---|
| 262 | F(7,7)=63.12 |
|---|
| 263 | F(7,8)=67.42 |
|---|
| 264 | F(7,9)=70.85 |
|---|
| 265 | F(7,10)=73.70 |
|---|
| 266 | F(7,11)=76.09 |
|---|
| 267 | F(7,12)=78.15 |
|---|
| 268 | F(7,13)=79.96 |
|---|
| 269 | F(7,14)=81.56 |
|---|
| 270 | F(7,15)=83.02 |
|---|
| 271 | F(7,16)=84.43 |
|---|
| 272 | F(7,17)=85.85 |
|---|
| 273 | F(7,18)=87.33 |
|---|
| 274 | |
|---|
| 275 | A=(F(7,5)-F(7,4))/( (YC(5)-YC(4))*2.0) |
|---|
| 276 | B=-A*YC(4) + F(7,4) |
|---|
| 277 | F(7,1)=A*YC(1) + B |
|---|
| 278 | F(7,2)=A*YC(2) + B |
|---|
| 279 | F(7,3)=A*YC(3) + B |
|---|
| 280 | |
|---|
| 281 | F(8,4)=25.85 |
|---|
| 282 | F(8,5)=42.26 |
|---|
| 283 | F(8,6)=52.78 |
|---|
| 284 | F(8,7)=59.55 |
|---|
| 285 | F(8,8)=64.55 |
|---|
| 286 | F(8,9)=68.45 |
|---|
| 287 | F(8,10)=71.63 |
|---|
| 288 | F(8,11)=74.29 |
|---|
| 289 | F(8,12)=76.56 |
|---|
| 290 | F(8,13)=78.53 |
|---|
| 291 | F(8,14)=80.27 |
|---|
| 292 | F(8,15)=81.83 |
|---|
| 293 | F(8,16)=83.27 |
|---|
| 294 | F(8,17)=84.67 |
|---|
| 295 | F(8,18)=86.10 |
|---|
| 296 | |
|---|
| 297 | A=(F(8,5)-F(8,4))/( (YC(5)-YC(4))*2.5 ) |
|---|
| 298 | B=-A*YC(4) + F(8,4) |
|---|
| 299 | F(8,1)=A*YC(1) + B |
|---|
| 300 | F(8,2)=A*YC(2) + B |
|---|
| 301 | F(8,3)=A*YC(3) + B |
|---|
| 302 | |
|---|
| 303 | F(9,4)=15.38 |
|---|
| 304 | F(9,5)=39.35 |
|---|
| 305 | F(9,6)=50.73 |
|---|
| 306 | F(9,7)=58.11 |
|---|
| 307 | F(9,8)=63.41 |
|---|
| 308 | F(9,9)=67.52 |
|---|
| 309 | F(9,10)=70.83 |
|---|
| 310 | F(9,11)=73.6 |
|---|
| 311 | F(9,12)=75.95 |
|---|
| 312 | F(9,13)=77.98 |
|---|
| 313 | F(9,14)=79.77 |
|---|
| 314 | F(9,15)=81.38 |
|---|
| 315 | F(9,16)=82.84 |
|---|
| 316 | F(9,17)=84.25 |
|---|
| 317 | F(9,18)=85.66 |
|---|
| 318 | |
|---|
| 319 | A=(F(9,5)-F(9,4))/( (YC(5)-YC(4))*7.0) |
|---|
| 320 | B=-A*YC(4) + F(9,4) |
|---|
| 321 | F(9,1)=A*YC(1) + B |
|---|
| 322 | F(9,2)=A*YC(2) + B |
|---|
| 323 | F(9,3)=A*YC(3) + B |
|---|
| 324 | |
|---|
| 325 | F(10,4)=0.0 |
|---|
| 326 | F(10,5)=34.02 |
|---|
| 327 | F(10,6)=46.93 |
|---|
| 328 | F(10,7)=55.61 |
|---|
| 329 | F(10,8)=61.47 |
|---|
| 330 | F(10,9)=65.94 |
|---|
| 331 | F(10,10)=69.49 |
|---|
| 332 | F(10,11)=72.44 |
|---|
| 333 | F(10,12)=74.93 |
|---|
| 334 | F(10,13)=77.08 |
|---|
| 335 | F(10,14)=78.96 |
|---|
| 336 | F(10,15)=80.63 |
|---|
| 337 | F(10,16)=82.15 |
|---|
| 338 | F(10,17)=83.57 |
|---|
| 339 | F(10,18)=84.97 |
|---|
| 340 | |
|---|
| 341 | A=(F(10,6)-F(10,5))/( (YC(6)-YC(5))*1.5) |
|---|
| 342 | B=-A*YC(5) + F(10,5) |
|---|
| 343 | F(10,1)=A*YC(1) + B |
|---|
| 344 | F(10,2)=A*YC(2) + B |
|---|
| 345 | F(10,3)=A*YC(3) + B |
|---|
| 346 | F(10,4)=A*YC(4) + B |
|---|
| 347 | |
|---|
| 348 | F(11,4)=0.0 |
|---|
| 349 | F(11,5)=29.02 |
|---|
| 350 | F(11,6)=43.69 |
|---|
| 351 | F(11,7)=53.44 |
|---|
| 352 | F(11,8)=59.83 |
|---|
| 353 | F(11,9)=64.62 |
|---|
| 354 | F(11,10)=68.39 |
|---|
| 355 | F(11,11)=71.48 |
|---|
| 356 | F(11,12)=74.10 |
|---|
| 357 | F(11,13)=76.33 |
|---|
| 358 | F(11,14)=78.29 |
|---|
| 359 | F(11,15)=80.02 |
|---|
| 360 | F(11,16)=81.58 |
|---|
| 361 | F(11,17)=83.03 |
|---|
| 362 | F(11,18)=84.44 |
|---|
| 363 | |
|---|
| 364 | A=(F(11,6)-F(11,5))/( (YC(6)-YC(5))*2.5 ) |
|---|
| 365 | B=-A*YC(5) + F(11,5) |
|---|
| 366 | F(11,1)=A*YC(1) + B |
|---|
| 367 | F(11,2)=A*YC(2) + B |
|---|
| 368 | F(11,3)=A*YC(3) + B |
|---|
| 369 | F(11,4)=A*YC(4) + B |
|---|
| 370 | |
|---|
| 371 | F(12,4)=0.0 |
|---|
| 372 | F(12,5)=23.13 |
|---|
| 373 | F(12,6)=40.86 |
|---|
| 374 | F(12,7)=51.44 |
|---|
| 375 | F(12,8)=58.38 |
|---|
| 376 | F(12,9)=63.47 |
|---|
| 377 | F(12,10)=67.43 |
|---|
| 378 | F(12,11)=70.66 |
|---|
| 379 | F(12,12)=73.38 |
|---|
| 380 | F(12,13)=75.70 |
|---|
| 381 | F(12,14)=77.72 |
|---|
| 382 | F(12,15)=79.51 |
|---|
| 383 | F(12,16)=81.11 |
|---|
| 384 | F(12,17)=82.58 |
|---|
| 385 | F(12,18)=83.99 |
|---|
| 386 | |
|---|
| 387 | A=(F(12,6)-F(12,5))/( (YC(6)-YC(5))*3.5 ) |
|---|
| 388 | B=-A*YC(5) + F(12,5) |
|---|
| 389 | F(12,1)=A*YC(1) + B |
|---|
| 390 | F(12,2)=A*YC(2) + B |
|---|
| 391 | F(12,3)=A*YC(3) + B |
|---|
| 392 | F(12,4)=A*YC(4) + B |
|---|
| 393 | |
|---|
| 394 | F(13,4)=0.0 |
|---|
| 395 | F(13,5)=0.0 |
|---|
| 396 | F(13,6)=36.89 |
|---|
| 397 | F(13,7)=48.63 |
|---|
| 398 | F(13,8)=56.46 |
|---|
| 399 | F(13,9)=61.96 |
|---|
| 400 | F(13,10)=66.19 |
|---|
| 401 | F(13,11)=69.6 |
|---|
| 402 | F(13,12)=72.45 |
|---|
| 403 | F(13,13)=74.89 |
|---|
| 404 | F(13,14)=76.99 |
|---|
| 405 | F(13,15)=78.85 |
|---|
| 406 | F(13,16)=80.50 |
|---|
| 407 | F(13,17)=82.02 |
|---|
| 408 | F(13,18)=83.44 |
|---|
| 409 | |
|---|
| 410 | A=(F(13,7)-F(13,6))/( (YC(7)-YC(6))*2.0) |
|---|
| 411 | B=-A*YC(6) + F(13,6) |
|---|
| 412 | F(13,1)=A*YC(1) + B |
|---|
| 413 | F(13,2)=A*YC(2) + B |
|---|
| 414 | F(13,3)=A*YC(3) + B |
|---|
| 415 | F(13,4)=A*YC(4) + B |
|---|
| 416 | F(13,5)=A*YC(5) + B |
|---|
| 417 | |
|---|
| 418 | F(14,4)=0.0 |
|---|
| 419 | F(14,5)=0.0 |
|---|
| 420 | F(14,6)=30.82 |
|---|
| 421 | F(14,7)=44.49 |
|---|
| 422 | F(14,8)=53.69 |
|---|
| 423 | F(14,9)=59.83 |
|---|
| 424 | F(14,10)=64.47 |
|---|
| 425 | F(14,11)=68.15 |
|---|
| 426 | F(14,12)=71.19 |
|---|
| 427 | F(14,13)=73.77 |
|---|
| 428 | F(14,14)=76.0 |
|---|
| 429 | F(14,15)=77.95 |
|---|
| 430 | F(14,16)=79.69 |
|---|
| 431 | F(14,17)=81.26 |
|---|
| 432 | F(14,18)=82.72 |
|---|
| 433 | |
|---|
| 434 | A=(F(14,7)-F(14,6))/( (YC(7)-YC(6))*2.5 ) |
|---|
| 435 | B=-A*YC(6) + F(14,6) |
|---|
| 436 | F(14,1)=A*YC(1) + B |
|---|
| 437 | F(14,2)=A*YC(2) + B |
|---|
| 438 | F(14,3)=A*YC(3) + B |
|---|
| 439 | F(14,4)=A*YC(4) + B |
|---|
| 440 | F(14,5)=A*YC(5) + B |
|---|
| 441 | |
|---|
| 442 | F(15,4)=0.0 |
|---|
| 443 | F(15,5)=0.0 |
|---|
| 444 | F(15,6)=0.0 |
|---|
| 445 | F(15,7)=37.71 |
|---|
| 446 | F(15,8)=48.49 |
|---|
| 447 | F(15,9)=56.40 |
|---|
| 448 | F(15,10)=61.75 |
|---|
| 449 | F(15,11)=65.89 |
|---|
| 450 | F(15,12)=69.25 |
|---|
| 451 | F(15,13)=72.07 |
|---|
| 452 | F(15,14)=74.49 |
|---|
| 453 | F(15,15)=76.59 |
|---|
| 454 | F(15,16)=78.45 |
|---|
| 455 | F(15,17)=80.12 |
|---|
| 456 | F(15,18)=81.64 |
|---|
| 457 | |
|---|
| 458 | A=(F(15,8)-F(15,7))/( (YC(8)-YC(7))*1.5) |
|---|
| 459 | B=-A*YC(7) + F(15,7) |
|---|
| 460 | F(15,1)=A*YC(1) + B |
|---|
| 461 | F(15,2)=A*YC(2) + B |
|---|
| 462 | F(15,3)=A*YC(3) + B |
|---|
| 463 | F(15,4)=A*YC(4) + B |
|---|
| 464 | F(15,5)=A*YC(5) + B |
|---|
| 465 | F(15,6)=A*YC(6) + B |
|---|
| 466 | |
|---|
| 467 | ! SUPPOSE THAT AT GIVEN AND PH2O<2mB, |
|---|
| 468 | ! %H2SO4 = A *LOG(PH2O) +B |
|---|
| 469 | ! XC(1-5) :EXTENSION LEFT (LOW H2O) |
|---|
| 470 | DO J=1,18 |
|---|
| 471 | A=(F(6,J)-F(7,J))/(LOG(XC(6))-LOG(XC(7))) |
|---|
| 472 | B=-A*LOG(XC(6)) + F(6,J) |
|---|
| 473 | DO K=1,5 |
|---|
| 474 | F(K,J)=A*LOG(XC(K)) + B |
|---|
| 475 | ENDDO |
|---|
| 476 | ENDDO |
|---|
| 477 | |
|---|
| 478 | ! XC(16) :EXTENSION RIGHT (HIGH H2O) |
|---|
| 479 | DO J=1,18 |
|---|
| 480 | A=(F(15,J)-F(14,J))/(XC(15)-XC(14)) |
|---|
| 481 | B=-A*XC(15) + F(15,J) |
|---|
| 482 | F(16,J)=A*XC(16) + B |
|---|
| 483 | ! F(16,2)=1.0 |
|---|
| 484 | ENDDO |
|---|
| 485 | |
|---|
| 486 | ! YC(16-25) :EXTENSION DOWN (HIGH T) |
|---|
| 487 | DO I=1,16 |
|---|
| 488 | A=(F(I,18)-F(I,17))/(YC(18)-YC(17)) |
|---|
| 489 | B=-A*YC(18) + F(I,18) |
|---|
| 490 | DO K=19,28 |
|---|
| 491 | F(I,K)=A*YC(K) + B |
|---|
| 492 | ENDDO |
|---|
| 493 | ENDDO |
|---|
| 494 | |
|---|
| 495 | ! MANUAL CORRECTIONS |
|---|
| 496 | DO J=1,10 |
|---|
| 497 | F(1,J)=94.0 |
|---|
| 498 | ENDDO |
|---|
| 499 | |
|---|
| 500 | DO J=1,6 |
|---|
| 501 | F(2,J)=77.0 +REAL(J) |
|---|
| 502 | ENDDO |
|---|
| 503 | |
|---|
| 504 | DO J=1,7 |
|---|
| 505 | F(16,J)=9.0 |
|---|
| 506 | ENDDO |
|---|
| 507 | |
|---|
| 508 | DO I=1,16 |
|---|
| 509 | DO J=1,28 |
|---|
| 510 | IF (F(I,J)<9.0) F(I,J)=30.0 |
|---|
| 511 | IF (F(I,J)>99.99) F(I,J)=99.99 |
|---|
| 512 | ENDDO |
|---|
| 513 | ENDDO |
|---|
| 514 | |
|---|
| 515 | ENDIF |
|---|
| 516 | |
|---|
| 517 | DO I=1,klon |
|---|
| 518 | DO J=1,klev |
|---|
| 519 | TP=t_seri(I,J) |
|---|
| 520 | IF (TP<175.1) TP=175.1 |
|---|
| 521 | ! Partial pressure of H2O (mb) |
|---|
| 522 | PH2O =PMB(I,J)*H2O(I,J) |
|---|
| 523 | IF (PH2O<XC1) THEN |
|---|
| 524 | R2SO4(I,J)=99.99 |
|---|
| 525 | ! PH2O=XC(1)+1.0E-10 |
|---|
| 526 | ELSE |
|---|
| 527 | IF (PH2O>XC16) PH2O=XC16 |
|---|
| 528 | ! SIMPLE LINEAR INTERPOLATIONS |
|---|
| 529 | CALL FIND(PH2O,TP,XC,YC,F,VAL,N,M) |
|---|
| 530 | IF (PMB(I,J)>=10.0.AND.VAL<60.0) VAL=60.0 |
|---|
| 531 | R2SO4(I,J)=VAL |
|---|
| 532 | ENDIF |
|---|
| 533 | ENDDO |
|---|
| 534 | ENDDO |
|---|
| 535 | |
|---|
| 536 | END SUBROUTINE |
|---|
| 537 | |
|---|
| 538 | !**************************************************************** |
|---|
| 539 | SUBROUTINE STRAACT(ACTSO4) |
|---|
| 540 | |
|---|
| 541 | ! H2SO4 ACTIVITY (GIAUQUE) AS A FUNCTION OF H2SO4 WP |
|---|
| 542 | ! ---------------------------------------- |
|---|
| 543 | ! INPUT: |
|---|
| 544 | ! H2SO4: VMR of H2SO4 |
|---|
| 545 | ! klon: number of latitude bands in the model domain |
|---|
| 546 | ! klev: number of altitude bands in the model domain |
|---|
| 547 | ! for IFS: perhaps add another dimension for longitude |
|---|
| 548 | |
|---|
| 549 | ! OUTPUT: |
|---|
| 550 | ! ACTSO4: H2SO4 activity (percent) |
|---|
| 551 | |
|---|
| 552 | USE dimphy, ONLY: klon,klev |
|---|
| 553 | USE phys_local_var_mod, ONLY: R2SO4 |
|---|
| 554 | |
|---|
| 555 | IMPLICIT NONE |
|---|
| 556 | |
|---|
| 557 | REAL ACTSO4(klon,klev) |
|---|
| 558 | |
|---|
| 559 | ! Working variables |
|---|
| 560 | INTEGER NN,I,J,JX,JX1 |
|---|
| 561 | REAL TC,TB,TA,XT |
|---|
| 562 | PARAMETER (NN=109) |
|---|
| 563 | REAL XC(NN), X(NN) |
|---|
| 564 | |
|---|
| 565 | ! H2SO4 activity |
|---|
| 566 | DATA X/ & |
|---|
| 567 | 0.0,0.25,0.78,1.437,2.19,3.07,4.03,5.04,6.08 & |
|---|
| 568 | ,7.13,8.18,14.33,18.59,28.59,39.17,49.49 & |
|---|
| 569 | ,102.4,157.8,215.7,276.9,341.6,409.8,481.5,556.6 & |
|---|
| 570 | ,635.5,719.,808.,902.,1000.,1103.,1211.,1322.,1437.,1555. & |
|---|
| 571 | ,1677.,1800.,1926.,2054.,2183.,2312.,2442.,2572.,2701.,2829. & |
|---|
| 572 | ,2955.,3080.,3203.,3325.,3446.,3564.,3681.,3796.,3910.,4022. & |
|---|
| 573 | ,4134.,4351.,4564.,4771.,4974.,5171.,5364.,5551.,5732.,5908. & |
|---|
| 574 | ,6079.,6244.,6404.,6559.,6709.,6854.,6994.,7131.,7264.,7393. & |
|---|
| 575 | ,7520.,7821.,8105.,8373.,8627.,8867.,9093.,9308.,9511.,9703. & |
|---|
| 576 | ,9885.,10060.,10225.,10535.,10819.,11079.,11318.,11537. & |
|---|
| 577 | ,11740.,12097.,12407.,12676.,12915.,13126.,13564.,13910. & |
|---|
| 578 | ,14191.,14423.,14617.,14786.,10568.,15299.,15491.,15654. & |
|---|
| 579 | ,15811./ |
|---|
| 580 | ! H2SO4 weight fraction (percent) |
|---|
| 581 | DATA XC/ & |
|---|
| 582 | 100.0,99.982,99.963,99.945,99.927,99.908,99.890,99.872 & |
|---|
| 583 | ,99.853,99.835,99.817,99.725,99.634,99.452,99.270 & |
|---|
| 584 | ,99.090,98.196,97.319,96.457,95.610,94.777,93.959,93.156 & |
|---|
| 585 | ,92.365,91.588,90.824,90.073,89.334,88.607,87.892,87.188 & |
|---|
| 586 | ,86.495,85.814,85.143,84.482,83.832,83.191,82.560,81.939 & |
|---|
| 587 | ,81.327,80.724,80.130,79.545,78.968,78.399,77.839,77.286 & |
|---|
| 588 | ,76.741,76.204,75.675,75.152,74.637,74.129,73.628,73.133 & |
|---|
| 589 | ,72.164,71.220,70.300,69.404,68.530,67.678,66.847,66.037 & |
|---|
| 590 | ,65.245,64.472,63.718,62.981,62.261,61.557,60.868,60.195 & |
|---|
| 591 | ,59.537,58.893,58.263,57.646,56.159,54.747,53.405,52.126 & |
|---|
| 592 | ,50.908,49.745,48.634,47.572,46.555,45.580,44.646,43.749 & |
|---|
| 593 | ,42.059,40.495,39.043,37.691,36.430,35.251,33.107,31.209 & |
|---|
| 594 | ,29.517,27.999,26.629,23.728,21.397,19.482,17.882,16.525 & |
|---|
| 595 | ,15.360,13.461,11.980,10.792,9.819,8.932/ |
|---|
| 596 | |
|---|
| 597 | DO I=1,klon |
|---|
| 598 | DO J=1,klev |
|---|
| 599 | ! HERE LINEAR INTERPOLATIONS |
|---|
| 600 | XT=R2SO4(I,J) |
|---|
| 601 | CALL POSACT(XT,XC,NN,JX) |
|---|
| 602 | JX1=JX+1 |
|---|
| 603 | IF(JX==0) THEN |
|---|
| 604 | ACTSO4(I,J)=0.0 |
|---|
| 605 | ELSE IF(JX>=NN) THEN |
|---|
| 606 | ACTSO4(I,J)=15811.0 |
|---|
| 607 | ELSE |
|---|
| 608 | TC=XT -XC(JX) |
|---|
| 609 | TB=X(JX1) -X(JX) |
|---|
| 610 | TA=XC(JX1) -XC(JX) |
|---|
| 611 | TA=TB/TA |
|---|
| 612 | ACTSO4(I,J)=X(JX) + TA*TC |
|---|
| 613 | ENDIF |
|---|
| 614 | ENDDO |
|---|
| 615 | ENDDO |
|---|
| 616 | |
|---|
| 617 | END SUBROUTINE |
|---|
| 618 | |
|---|
| 619 | !**************************************************************** |
|---|
| 620 | SUBROUTINE DENH2SA(t_seri) |
|---|
| 621 | |
|---|
| 622 | ! AERSOL DENSITY AS A FUNCTION OF H2SO4 WEIGHT PERCENT AND T |
|---|
| 623 | ! --------------------------------------------- |
|---|
| 624 | ! VERY ROUGH APPROXIMATION (SEE FOR WATER IN HANDBOOK |
|---|
| 625 | ! LINEAR 2% FOR 30 DEGREES with RESPECT TO WATER) |
|---|
| 626 | |
|---|
| 627 | ! INPUT: |
|---|
| 628 | ! R2SO4: aerosol H2SO4 weight fraction (percent) |
|---|
| 629 | ! t_seri: temperature (K) |
|---|
| 630 | ! klon: number of latitude bands in the model domain |
|---|
| 631 | ! klev: number of altitude bands in the model domain |
|---|
| 632 | ! for IFS: perhaps add another dimension for longitude |
|---|
| 633 | |
|---|
| 634 | ! OUTPUT: |
|---|
| 635 | ! DENSO4: aerosol mass density (gr/cm3 = aerosol mass/aerosol volume) |
|---|
| 636 | |
|---|
| 637 | USE dimphy, ONLY: klon,klev |
|---|
| 638 | USE phys_local_var_mod, ONLY: R2SO4, DENSO4 |
|---|
| 639 | |
|---|
| 640 | IMPLICIT NONE |
|---|
| 641 | |
|---|
| 642 | REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature |
|---|
| 643 | |
|---|
| 644 | INTEGER I,J |
|---|
| 645 | |
|---|
| 646 | ! Loop on model domain (2 dimension for UPMC model; 3 for IFS) |
|---|
| 647 | DO I=1,klon |
|---|
| 648 | DO J=1,klev |
|---|
| 649 | ! RO AT 20C |
|---|
| 650 | DENSO4(I,J)=0.78681252E-5*R2SO4(I,J)*R2SO4(I,J)+ 0.82185978E-2*R2SO4(I,J)+0.97968381 |
|---|
| 651 | DENSO4(I,J)=DENSO4(I,J)* ( 1.0 - (t_seri(I,J)-293.0)*0.02/30.0 ) |
|---|
| 652 | ENDDO |
|---|
| 653 | ENDDO |
|---|
| 654 | |
|---|
| 655 | END SUBROUTINE |
|---|
| 656 | |
|---|
| 657 | !*********************************************************** |
|---|
| 658 | SUBROUTINE FIND(X,Y,XC,YC,F,VAL,N,M) |
|---|
| 659 | |
|---|
| 660 | ! BI-LINEAR INTERPOLATION |
|---|
| 661 | |
|---|
| 662 | ! INPUT: |
|---|
| 663 | ! X: Partial pressure of H2O (mb) |
|---|
| 664 | ! Y: temperature (K) |
|---|
| 665 | ! XC: Table partial pressure of H2O (mb) |
|---|
| 666 | ! YC: Table temperature (K) |
|---|
| 667 | ! F: Table aerosol H2SO4 weight fraction=f(XC,YC) (percent) |
|---|
| 668 | |
|---|
| 669 | ! OUTPUT: |
|---|
| 670 | ! VAL: aerosol H2SO4 weight fraction (percent) |
|---|
| 671 | |
|---|
| 672 | IMPLICIT NONE |
|---|
| 673 | |
|---|
| 674 | INTEGER N,M |
|---|
| 675 | REAL X,Y,XC(N),YC(M),F(N,M),VAL |
|---|
| 676 | |
|---|
| 677 | ! working variables |
|---|
| 678 | INTEGER IERX,IERY,JX,JY,JXP1,JYP1 |
|---|
| 679 | REAL SXY,SX1Y,SX1Y1,SXY1,TA,TB,T,UA,UB,U |
|---|
| 680 | |
|---|
| 681 | IERX=0 |
|---|
| 682 | IERY=0 |
|---|
| 683 | CALL POSITION(XC,X,N,JX,IERX) |
|---|
| 684 | CALL POSITION(YC,Y,M,JY,IERY) |
|---|
| 685 | |
|---|
| 686 | IF(JX==0.OR.IERY==1) THEN |
|---|
| 687 | VAL=99.99 |
|---|
| 688 | RETURN |
|---|
| 689 | ENDIF |
|---|
| 690 | |
|---|
| 691 | IF(JY==0.OR.IERX==1) THEN |
|---|
| 692 | VAL=9.0 |
|---|
| 693 | RETURN |
|---|
| 694 | ENDIF |
|---|
| 695 | |
|---|
| 696 | JXP1=JX+1 |
|---|
| 697 | JYP1=JY+1 |
|---|
| 698 | SXY=F(JX, JY ) |
|---|
| 699 | SX1Y=F(JXP1,JY ) |
|---|
| 700 | SX1Y1=F(JXP1,JYP1) |
|---|
| 701 | SXY1=F(JX, JYP1) |
|---|
| 702 | |
|---|
| 703 | ! x-slope. |
|---|
| 704 | TA=X -XC(JX) |
|---|
| 705 | TB=XC(JXP1)-XC(JX) |
|---|
| 706 | T=TA/TB |
|---|
| 707 | |
|---|
| 708 | ! y-slope. |
|---|
| 709 | UA=Y -YC(JY) |
|---|
| 710 | UB=YC(JYP1)-YC(JY) |
|---|
| 711 | U=UA/UB |
|---|
| 712 | |
|---|
| 713 | ! Use bilinear interpolation to determine function at point X,Y. |
|---|
| 714 | VAL=(1.-T)*(1.-U)*SXY + T*(1.0-U)*SX1Y + T*U*SX1Y1 + (1.0-T)*U*SXY1 |
|---|
| 715 | |
|---|
| 716 | IF(VAL<9.0) VAL=9.0 |
|---|
| 717 | IF(VAL>99.99) VAL=99.99 |
|---|
| 718 | |
|---|
| 719 | |
|---|
| 720 | END SUBROUTINE |
|---|
| 721 | !**************************************************************** |
|---|
| 722 | SUBROUTINE POSITION(XC,X,N,JX,IER) |
|---|
| 723 | |
|---|
| 724 | IMPLICIT NONE |
|---|
| 725 | |
|---|
| 726 | INTEGER N,JX,IER,I |
|---|
| 727 | REAL X,XC(N) |
|---|
| 728 | |
|---|
| 729 | IER=0 |
|---|
| 730 | IF(X<XC(1)) THEN |
|---|
| 731 | JX=0 |
|---|
| 732 | ELSE |
|---|
| 733 | DO I=1,N |
|---|
| 734 | IF (X<XC(I)) GO TO 20 |
|---|
| 735 | END DO |
|---|
| 736 | IER=1 |
|---|
| 737 | 20 JX=I-1 |
|---|
| 738 | ENDIF |
|---|
| 739 | |
|---|
| 740 | |
|---|
| 741 | END SUBROUTINE |
|---|
| 742 | !******************************************************************** |
|---|
| 743 | SUBROUTINE POSACT(XT,X,N,JX) |
|---|
| 744 | |
|---|
| 745 | ! POSITION OF XT IN THE ARRAY X |
|---|
| 746 | ! ----------------------------------------------- |
|---|
| 747 | |
|---|
| 748 | IMPLICIT NONE |
|---|
| 749 | |
|---|
| 750 | INTEGER N |
|---|
| 751 | REAL XT,X(N) |
|---|
| 752 | ! Working variables |
|---|
| 753 | INTEGER JX,I |
|---|
| 754 | |
|---|
| 755 | IF(XT>X(1)) THEN |
|---|
| 756 | JX=0 |
|---|
| 757 | ELSE |
|---|
| 758 | DO I=1,N |
|---|
| 759 | IF (XT>X(I)) GO TO 20 |
|---|
| 760 | END DO |
|---|
| 761 | 20 JX=I |
|---|
| 762 | ENDIF |
|---|
| 763 | |
|---|
| 764 | |
|---|
| 765 | END SUBROUTINE |
|---|
| 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 | |
|---|
| 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 | |
|---|
| 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>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 | |
|---|
| 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 | |
|---|
| 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 | |
|---|
| 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 <= 0.05 .AND. aw > 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 <= 0.85 .AND. aw > 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 | |
|---|
| 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 | |
|---|
| 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 | |
|---|
| 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 | |
|---|
| 1060 | END FUNCTION rmvh2o |
|---|
| 1061 | |
|---|
| 1062 | END MODULE sulfate_aer_mod |
|---|
