| 1 | #include <inca_define.h> |
|---|
| 2 | #ifdef STRAT |
|---|
| 3 | MODULE HETCHEM |
|---|
| 4 | |
|---|
| 5 | ! This module contains the heterogeneous chemistry routines for |
|---|
| 6 | ! the stratosphere from the REPROBUS CTM provided by F. Lefevre. |
|---|
| 7 | ! They were adapted to LMDz-INCA by L. Jourdain and reworked by |
|---|
| 8 | ! by D. Hauglustaine and F. Jegou. DH, LSCE, 06/2005 |
|---|
| 9 | ! |
|---|
| 10 | ! Updated version 2018, DH. |
|---|
| 11 | |
|---|
| 12 | USE INCA_DIM |
|---|
| 13 | IMPLICIT NONE |
|---|
| 14 | |
|---|
| 15 | INTEGER, PARAMETER :: nlat = 1 |
|---|
| 16 | INTEGER, PARAMETER :: nivtop = 32 |
|---|
| 17 | INTEGER, PARAMETER :: nivbot = 15 |
|---|
| 18 | INTEGER, PARAMETER :: nbcon = PCNST |
|---|
| 19 | |
|---|
| 20 | REAL, SAVE, ALLOCATABLE :: parthno3(:,:) |
|---|
| 21 | REAL, SAVE, ALLOCATABLE :: parthno3s(:,:) |
|---|
| 22 | REAL, SAVE, ALLOCATABLE :: parthno3l(:,:) |
|---|
| 23 | REAL, SAVE, ALLOCATABLE :: parthno3sl(:,:) |
|---|
| 24 | |
|---|
| 25 | REAL, SAVE, ALLOCATABLE :: parth2o(:,:) |
|---|
| 26 | REAL, SAVE, ALLOCATABLE :: parth2os(:,:) |
|---|
| 27 | REAL, SAVE, ALLOCATABLE :: parth2ol(:,:) |
|---|
| 28 | |
|---|
| 29 | REAL, SAVE, ALLOCATABLE :: parthcl(:,:) |
|---|
| 30 | REAL, SAVE, ALLOCATABLE :: parthcll(:,:) |
|---|
| 31 | REAL, SAVE, ALLOCATABLE :: parthcls(:,:) |
|---|
| 32 | |
|---|
| 33 | REAL, SAVE, ALLOCATABLE :: parthbr(:,:) |
|---|
| 34 | REAL, SAVE, ALLOCATABLE :: parthbrl(:,:) |
|---|
| 35 | REAL, SAVE, ALLOCATABLE :: parthbrs(:,:) |
|---|
| 36 | |
|---|
| 37 | REAL, SAVE, ALLOCATABLE :: vsed(:,:) |
|---|
| 38 | REAL, SAVE, ALLOCATABLE :: surfarealiq(:,:) |
|---|
| 39 | REAL, SAVE, ALLOCATABLE :: surfareanat(:,:) |
|---|
| 40 | REAL, SAVE, ALLOCATABLE :: surfareaice(:,:) |
|---|
| 41 | REAL, SAVE, ALLOCATABLE :: qh2so4(:,:) |
|---|
| 42 | |
|---|
| 43 | !$OMP THREADPRIVATE(parthno3,parthno3s,parthno3l,parthno3sl,parth2o, parth2os,parth2ol) |
|---|
| 44 | !$OMP THREADPRIVATE(parthcl, parthcll,parthcls, parthbr,parthbrl,parthbrs,vsed,qh2so4) |
|---|
| 45 | !$OMP THREADPRIVATE(surfarealiq,surfareanat,surfareaice) |
|---|
| 46 | |
|---|
| 47 | CONTAINS |
|---|
| 48 | |
|---|
| 49 | SUBROUTINE INIT_HETCHEM |
|---|
| 50 | USE INCA_DIM |
|---|
| 51 | IMPLICIT NONE |
|---|
| 52 | |
|---|
| 53 | ALLOCATE ( vsed (PLON,PLEV)) |
|---|
| 54 | ALLOCATE ( parthno3 (PLON,PLEV)) |
|---|
| 55 | ALLOCATE ( parthno3s (PLON,PLEV)) |
|---|
| 56 | ALLOCATE ( parthno3l (PLON,PLEV)) |
|---|
| 57 | ALLOCATE ( parthno3sl (PLON,PLEV)) |
|---|
| 58 | ALLOCATE ( parth2o (PLON,PLEV)) |
|---|
| 59 | ALLOCATE ( parth2os (PLON,PLEV)) |
|---|
| 60 | ALLOCATE ( parth2ol (PLON,PLEV)) |
|---|
| 61 | ALLOCATE ( parthcl (PLON,PLEV)) |
|---|
| 62 | ALLOCATE ( parthcll (PLON,PLEV)) |
|---|
| 63 | ALLOCATE ( parthcls (PLON,PLEV)) |
|---|
| 64 | ALLOCATE ( parthbr (PLON,PLEV)) |
|---|
| 65 | ALLOCATE ( parthbrl (PLON,PLEV)) |
|---|
| 66 | ALLOCATE ( parthbrs (PLON,PLEV)) |
|---|
| 67 | ALLOCATE ( qh2so4 (PLON,PLEV)) |
|---|
| 68 | ALLOCATE ( surfarealiq (PLON,PLEV)) |
|---|
| 69 | ALLOCATE ( surfareanat (PLON,PLEV)) |
|---|
| 70 | ALLOCATE ( surfareaice (PLON,PLEV)) |
|---|
| 71 | |
|---|
| 72 | vsed(:,:) = 0. |
|---|
| 73 | parthno3(:,:) = 1. |
|---|
| 74 | parthno3s(:,:) = 1. |
|---|
| 75 | parthno3l(:,:) = 1. |
|---|
| 76 | parthno3sl(:,:) = 1. |
|---|
| 77 | parth2o(:,:) = 1. |
|---|
| 78 | parth2os(:,:) = 1. |
|---|
| 79 | parth2ol(:,:) = 1. |
|---|
| 80 | parthcl(:,:) = 1. |
|---|
| 81 | parthcll(:,:) = 1. |
|---|
| 82 | parthcls(:,:) = 1. |
|---|
| 83 | parthbr (:,:) = 1. |
|---|
| 84 | parthbrl(:,:) = 1. |
|---|
| 85 | parthbrs(:,:) = 1. |
|---|
| 86 | qh2so4(:,:) = 0. |
|---|
| 87 | surfarealiq(:,:) = 0. |
|---|
| 88 | surfareanat(:,:) = 0. |
|---|
| 89 | surfareaice(:,:) = 0. |
|---|
| 90 | |
|---|
| 91 | END SUBROUTINE INIT_HETCHEM |
|---|
| 92 | |
|---|
| 93 | SUBROUTINE ANALYTIC( & |
|---|
| 94 | temp, & |
|---|
| 95 | pmid, & |
|---|
| 96 | hnm, & |
|---|
| 97 | qj, & |
|---|
| 98 | k1, & |
|---|
| 99 | k2, & |
|---|
| 100 | k3, & |
|---|
| 101 | k4, & |
|---|
| 102 | k5, & |
|---|
| 103 | k6, & |
|---|
| 104 | k7, & |
|---|
| 105 | k8, & |
|---|
| 106 | k9, & |
|---|
| 107 | h2o, & |
|---|
| 108 | mair ) |
|---|
| 109 | |
|---|
| 110 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 111 | !c F. LEFEVRE (LATMOS/CNRS) from REPROBUS |
|---|
| 112 | !c |
|---|
| 113 | !c Adapted to INCA: |
|---|
| 114 | !c |
|---|
| 115 | !c D. HAUGLUSTAINE (LSCE) |
|---|
| 116 | !c L. JOURDAIN (LATMOS) |
|---|
| 117 | !c F. JEGOU (LSCE) |
|---|
| 118 | !c R. VALORSO (LIVE) |
|---|
| 119 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 120 | |
|---|
| 121 | USE SPECIES_NAMES |
|---|
| 122 | USE SAD_COM, ONLY : MASS,VOL |
|---|
| 123 | USE INCA_DIM |
|---|
| 124 | USE AIRPLANE_SRC, ONLY : ptrop |
|---|
| 125 | |
|---|
| 126 | REAL, PARAMETER :: r=8.205e-5 |
|---|
| 127 | REAL, PARAMETER :: ctoa=7.336e21 |
|---|
| 128 | REAL, PARAMETER :: sigma=1.6 |
|---|
| 129 | |
|---|
| 130 | INTEGER :: i, ilon, iniv, j |
|---|
| 131 | REAL :: t, ph2o, zh2o, zh2oeq |
|---|
| 132 | REAL :: pi, amu |
|---|
| 133 | REAL :: qh2o, qhno3, qhcl |
|---|
| 134 | REAL :: zmt, zbt, zhno3eq |
|---|
| 135 | REAL :: pn0t, pn |
|---|
| 136 | REAL :: hno3s, h2os |
|---|
| 137 | REAL :: xsb, hsb, xnb, hnb |
|---|
| 138 | REAL :: tf, tt, tr, pr |
|---|
| 139 | REAL :: phi |
|---|
| 140 | REAL :: wsf, wnf, partf |
|---|
| 141 | REAL :: vliq, rmode |
|---|
| 142 | REAL :: ph2o0, phcl1 |
|---|
| 143 | REAL :: sclono2, shocl |
|---|
| 144 | REAL :: rpart, rhopart |
|---|
| 145 | REAL :: ynre |
|---|
| 146 | REAL :: vsedliq, shapek, vsedsol |
|---|
| 147 | REAL :: wts |
|---|
| 148 | |
|---|
| 149 | REAL temp(PLON,PLEV), ptot(PLON,PLEV), hnm(PLON,PLEV) |
|---|
| 150 | REAL ck(PLON,PLEV,nbcon),qj(PLON,PLEV,nbcon) |
|---|
| 151 | REAL k1(PLON,PLEV), k2(PLON,PLEV), k3(PLON,PLEV) |
|---|
| 152 | REAL k4(PLON,PLEV), k5(PLON,PLEV), k6(PLON,PLEV) |
|---|
| 153 | REAL k7(PLON,PLEV), k8(PLON,PLEV), k9(PLON,PLEV) |
|---|
| 154 | REAL h2o(PLON,PLEV),pmid(PLON,PLEV) |
|---|
| 155 | REAL a,b,c,det,msb,mnb,mcl,mnf,msf |
|---|
| 156 | REAL qn(10),qs(10),kn(7),ks(7) |
|---|
| 157 | SAVE qn,qs,kn,ks |
|---|
| 158 | !$OMP THREADPRIVATE(qn,qs,kn,ks) |
|---|
| 159 | |
|---|
| 160 | REAL, INTENT(in) :: mair(PLON,PLEV) ! total atm density |
|---|
| 161 | REAL :: VOL_S(PLON,PLEV) |
|---|
| 162 | |
|---|
| 163 | REAL dens(PLON,PLEV) |
|---|
| 164 | REAL ns(PLON,PLEV) |
|---|
| 165 | REAL pn0(PLON,PLEV), phcl0(PLON,PLEV) |
|---|
| 166 | REAL pw(PLON,PLEV), ms(PLON,PLEV), mn(PLON,PLEV) |
|---|
| 167 | REAL ws(PLON,PLEV), wn(PLON,PLEV) |
|---|
| 168 | REAL hhcl(PLON,PLEV), hhocl(PLON,PLEV) |
|---|
| 169 | REAL hhobr(PLON,PLEV), hclono2(PLON,PLEV) |
|---|
| 170 | REAL hhbr(PLON,PLEV) |
|---|
| 171 | REAL tice(PLON,PLEV), rice(PLON,PLEV) |
|---|
| 172 | REAL rnat(PLON,PLEV), aliq(PLON,PLEV) |
|---|
| 173 | REAL rmean(PLON,PLEV) |
|---|
| 174 | REAL mh2so4c, bh2so4, ph2so4, muh2so4 |
|---|
| 175 | |
|---|
| 176 | INTEGER condliq(PLON,PLEV) |
|---|
| 177 | INTEGER lnat(PLON,PLEV), lice(PLON,PLEV) |
|---|
| 178 | !c |
|---|
| 179 | !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 180 | !c ajouts jpl 2000 |
|---|
| 181 | !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 182 | !c |
|---|
| 183 | REAL aw(PLON,PLEV) |
|---|
| 184 | REAL wt(PLON,PLEV), mh2so4(PLON,PLEV) |
|---|
| 185 | REAL m, y1, y2, bigx |
|---|
| 186 | REAL a1(3), b1(3), c1(3), d1(3) |
|---|
| 187 | REAL a2(3), b2(3), c2(3), d2(3) |
|---|
| 188 | |
|---|
| 189 | DATA a1/12.37208932, 11.820654354, -180.06541028/ |
|---|
| 190 | DATA b1/-0.16125516114, -0.20786404244, -0.38601102592/ |
|---|
| 191 | DATA c1/-30.490657554, -4.807306373, -93.317846778/ |
|---|
| 192 | DATA d1/-2.1133114241, -5.1727540348, 273.88132245/ |
|---|
| 193 | DATA a2/13.455394705, 12.891938068, -176.95814097/ |
|---|
| 194 | DATA b2/-0.1921312255,-0.232338847708,-0.36257048154/ |
|---|
| 195 | DATA c2/-34.285174607, -6.4261237757, -90.469744201/ |
|---|
| 196 | DATA d2/-1.7620073078, -4.9005471319, 267.45509988/ |
|---|
| 197 | |
|---|
| 198 | REAL ndrop(PLON,PLEV) |
|---|
| 199 | REAL nice(PLON,PLEV), nnat(PLON,PLEV) |
|---|
| 200 | REAL nucleimin |
|---|
| 201 | |
|---|
| 202 | REAL nnat_small_max, vnat_small_max |
|---|
| 203 | REAL nnats(PLON,PLEV) |
|---|
| 204 | REAL nnatl(PLON,PLEV) |
|---|
| 205 | REAL rnats, rnatl |
|---|
| 206 | REAL anat(PLON,PLEV), vnat(PLON,PLEV), vnatl(PLON,PLEV) |
|---|
| 207 | REAL aice(PLON,PLEV), vice(PLON,PLEV) |
|---|
| 208 | |
|---|
| 209 | REAL vsed(PLON,PLEV) |
|---|
| 210 | REAL rverif(PLON*PLEV) |
|---|
| 211 | REAL rhoair, rhoice, rhonat, rholiq |
|---|
| 212 | REAL netha, lambda, lambda0, nre, us, cdnre2 |
|---|
| 213 | REAL bnre(0:6), tcel |
|---|
| 214 | !c |
|---|
| 215 | !c pruppacher and klett, 1988 |
|---|
| 216 | !c |
|---|
| 217 | DATA bnre /-3.18657,0.992696,-0.153193e-2, & |
|---|
| 218 | -0.987059e-3,-0.578878e-3,0.855176e-4,& |
|---|
| 219 | -0.327815e-5/ |
|---|
| 220 | !c |
|---|
| 221 | !c carslaw et al., 1995 |
|---|
| 222 | !c |
|---|
| 223 | DATA qn/14.5734,0.0615994,-1.14895,0.691693, -0.098863,& |
|---|
| 224 | 0.0051579,0.123472,-0.115574,0.0110113,0.0097914/ |
|---|
| 225 | DATA qs/14.4700,0.0638795,-3.29597,1.778224,-0.223244, & |
|---|
| 226 | 0.0086486,0.536695,-0.335164,0.0265153,0.0157550/ |
|---|
| 227 | DATA kn/-39.136,6358.4,83.29,-17650.0,198.53,-11948,-28.469/ |
|---|
| 228 | DATA ks/-21.661,2724.2,51.81,-15732.0,47.004,-6969.0,-4.6183/ |
|---|
| 229 | |
|---|
| 230 | LOGICAL pscliq, pscnat, pscice |
|---|
| 231 | !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 232 | ! choix du schema de psc: |
|---|
| 233 | ! |
|---|
| 234 | ! pscnat = true : autoriser les nat |
|---|
| 235 | ! pscnat + pscice = true : autoriser les nat et la glace |
|---|
| 236 | ! pscnat + pscice + pscliq = true : autoriser nat, glace et liquide |
|---|
| 237 | ! |
|---|
| 238 | !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 239 | pscnat = .true. |
|---|
| 240 | pscice = .true. |
|---|
| 241 | pscliq = .true. |
|---|
| 242 | |
|---|
| 243 | !c |
|---|
| 244 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 245 | !c initialisations |
|---|
| 246 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 247 | !c |
|---|
| 248 | |
|---|
| 249 | nnat_small_max = 1. ! densite maximale de nat (cm-3) |
|---|
| 250 | rnats = 0.5e-4 ! rayon du small mode de nat (cm) |
|---|
| 251 | rnatl = 6.5e-4 ! rayon du large mode de nat (cm) |
|---|
| 252 | rice = 10.e-4 ! rayon de la glace |
|---|
| 253 | nucleimin = 1.e-3 |
|---|
| 254 | pi = 2.*ASIN(1.0) |
|---|
| 255 | amu = 1.66e-24 |
|---|
| 256 | ptot(:,:)=pmid(:,:)/100. |
|---|
| 257 | |
|---|
| 258 | vsed = 0. |
|---|
| 259 | ws(:,:) = 0. |
|---|
| 260 | |
|---|
| 261 | ! Remove values below tropopause for sulfates |
|---|
| 262 | DO i=1,PLON |
|---|
| 263 | DO j=1,PLEV |
|---|
| 264 | IF (pmid(i,j).GT.ptrop(i)) THEN |
|---|
| 265 | VOL_S(i,j)=0 |
|---|
| 266 | ELSE |
|---|
| 267 | VOL_S(i,j)=VOL(i,j) |
|---|
| 268 | ENDIF |
|---|
| 269 | ENDDO |
|---|
| 270 | ENDDO |
|---|
| 271 | |
|---|
| 272 | DO ilon = 1, PLON |
|---|
| 273 | DO iniv = nivbot, nivtop |
|---|
| 274 | |
|---|
| 275 | qh2o = h2o(ilon,iniv) |
|---|
| 276 | pw(ilon,iniv) = qh2o*ptot(ilon,iniv)/1013.0 |
|---|
| 277 | pw(ilon,iniv) = MIN(pw(ilon,iniv),2.e-3/1013.) |
|---|
| 278 | |
|---|
| 279 | qhno3 = qj(ilon,iniv,id_HNO3) |
|---|
| 280 | pn0(ilon,iniv)= ptot(ilon,iniv)/1013.0*qhno3 |
|---|
| 281 | |
|---|
| 282 | qhcl = qj(ilon,iniv,id_HCl) |
|---|
| 283 | phcl0(ilon,iniv)=ptot(ilon,iniv)/1013.0*qhcl |
|---|
| 284 | |
|---|
| 285 | !c |
|---|
| 286 | !c all h2so4 is assumed to be in the droplets. you might |
|---|
| 287 | !c want to input directly the moles h2so4/m3, which is a more |
|---|
| 288 | !c normal unit for h2so4 amounts, but remember that this will |
|---|
| 289 | !c vary with pressure and temperature. pinatubo can be simulated |
|---|
| 290 | !c by multiplying qh2so4 by different factors. |
|---|
| 291 | !c |
|---|
| 292 | !c DH. Note: so4 is volume_aerosols/volume_air and is based on aerosol surface |
|---|
| 293 | !c density from Considine provided from SAGE and in um2/cm3. Then we transform |
|---|
| 294 | !c the volume into nb moles below. |
|---|
| 295 | |
|---|
| 296 | !qh2so4 = MAX((MASS(ilon,iniv)/mair(ilon,iniv)),1.e-20) |
|---|
| 297 | !qh2so4 = MAX(VOL_S(ilon,iniv)*1e-12,1.e-20) |
|---|
| 298 | !qh2so4 = MAX(so4(ilon,iniv),1.e-20) |
|---|
| 299 | !qh2so4 = 1.e-20 |
|---|
| 300 | |
|---|
| 301 | !c The line below replaced by the calculation below: |
|---|
| 302 | !c ns(ilon,iniv)=ptot(ilon,iniv)*100.0*qh2so4/8.314/t |
|---|
| 303 | |
|---|
| 304 | ndrop(ilon,iniv) = 10. |
|---|
| 305 | |
|---|
| 306 | END DO |
|---|
| 307 | END DO |
|---|
| 308 | |
|---|
| 309 | DO ilon = 1, PLON |
|---|
| 310 | DO iniv = nivbot, nivtop |
|---|
| 311 | lnat(ilon,iniv) = 0 |
|---|
| 312 | lice(ilon,iniv) = 0 |
|---|
| 313 | rmean(ilon,iniv) = 0. |
|---|
| 314 | aliq(ilon,iniv) = 0. |
|---|
| 315 | nnats(ilon,iniv) = 0. |
|---|
| 316 | nnatl(ilon,iniv) = 0. |
|---|
| 317 | nnat(ilon,iniv) = 0. |
|---|
| 318 | nice(ilon,iniv) = 0. |
|---|
| 319 | rnat(ilon,iniv) = 0. |
|---|
| 320 | vnat(ilon,iniv) = 0. |
|---|
| 321 | vnatl(ilon,iniv) = 0. |
|---|
| 322 | anat(ilon,iniv) = 0. |
|---|
| 323 | vice(ilon,iniv) = 0. |
|---|
| 324 | aice(ilon,iniv) = 0. |
|---|
| 325 | vsed(ilon,iniv) = 0. |
|---|
| 326 | parthno3s(ilon,iniv) = 1. |
|---|
| 327 | parthno3sl(ilon,iniv) = 1. |
|---|
| 328 | parthno3l(ilon,iniv) = 1. |
|---|
| 329 | parthno3(ilon,iniv) = 1. |
|---|
| 330 | parth2o(ilon,iniv) = 1. |
|---|
| 331 | parthcl(ilon,iniv) = 1. |
|---|
| 332 | END DO |
|---|
| 333 | END DO |
|---|
| 334 | |
|---|
| 335 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 336 | !c traitement des aerosols solides |
|---|
| 337 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 338 | |
|---|
| 339 | IF (pscnat) THEN |
|---|
| 340 | |
|---|
| 341 | !c |
|---|
| 342 | !cc traitement du nat |
|---|
| 343 | !c |
|---|
| 344 | |
|---|
| 345 | DO ilon = 1, PLON |
|---|
| 346 | DO iniv = nivbot, nivtop |
|---|
| 347 | |
|---|
| 348 | t = MIN(temp(ilon,iniv),273.) |
|---|
| 349 | |
|---|
| 350 | pn0t = pn0(ilon,iniv)*1013.*0.75 |
|---|
| 351 | zmt = -2.7836 - 0.00088*t |
|---|
| 352 | zbt = 38.9855 - 11397.0/t + 0.009179*t |
|---|
| 353 | zhno3eq = 10.0**(zmt*LOG10(pw(ilon,iniv)*1013.*0.75) + zbt) |
|---|
| 354 | |
|---|
| 355 | IF (pn0t/zhno3eq .GE. 1.) THEN |
|---|
| 356 | lnat(ilon,iniv) = 1 |
|---|
| 357 | pn = zhno3eq/(0.75*1013.) |
|---|
| 358 | hno3s = (pn0(ilon,iniv) - pn)*1013. & |
|---|
| 359 | *hnm(ilon,iniv)/ptot(ilon,iniv) |
|---|
| 360 | hno3s = MAX(hno3s,0.) |
|---|
| 361 | vnat(ilon,iniv) = hno3s*amu*117./1.35 |
|---|
| 362 | |
|---|
| 363 | ! volume maximal du small mode |
|---|
| 364 | vnat_small_max = nnat_small_max*(4./3.*pi*rnats**3.) |
|---|
| 365 | |
|---|
| 366 | ! densite de nat (cm-3) pour les deux modes |
|---|
| 367 | |
|---|
| 368 | if (vnat(ilon,iniv) .gt. vnat_small_max) then |
|---|
| 369 | nnats(ilon,iniv) = nnat_small_max |
|---|
| 370 | nnatl(ilon,iniv) = (vnat(ilon,iniv) - vnat_small_max)/(4./3.*pi*rnatl**3.) |
|---|
| 371 | vnatl(ilon,iniv) = vnat(ilon,iniv) - vnat_small_max |
|---|
| 372 | else |
|---|
| 373 | nnats(ilon,iniv) = vnat(ilon,iniv)/(4./3.*pi*rnats**3.) |
|---|
| 374 | nnatl(ilon,iniv) = 0. |
|---|
| 375 | vnatl(ilon,iniv) = 0. |
|---|
| 376 | end if |
|---|
| 377 | |
|---|
| 378 | ! surface du nat |
|---|
| 379 | |
|---|
| 380 | anat(ilon,iniv) = nnats(ilon,iniv)*4*pi*rnats**2& |
|---|
| 381 | + nnatl(ilon,iniv)*4*pi*rnatl**2 |
|---|
| 382 | |
|---|
| 383 | ! parthno3s: part en phase gazeuse |
|---|
| 384 | |
|---|
| 385 | parthno3s(ilon,iniv) = 1.0 - (pn0(ilon,iniv)-pn)& |
|---|
| 386 | /pn0(ilon,iniv) |
|---|
| 387 | parthno3s(ilon,iniv) = max(parthno3s(ilon,iniv),0.0) |
|---|
| 388 | parthno3s(ilon,iniv) = min(parthno3s(ilon,iniv),1.0) |
|---|
| 389 | |
|---|
| 390 | ! parthno3sl: part non sedimentee (gaz + petit mode de nat) |
|---|
| 391 | |
|---|
| 392 | parthno3sl(ilon,iniv) = 1.0& |
|---|
| 393 | - (pn0(ilon,iniv)-pn)*vnatl(ilon,iniv)/vnat(ilon,iniv)/pn0(ilon,iniv) |
|---|
| 394 | parthno3sl(ilon,iniv) = max(parthno3sl(ilon,iniv),0.0) |
|---|
| 395 | parthno3sl(ilon,iniv) = min(parthno3sl(ilon,iniv),1.0) |
|---|
| 396 | |
|---|
| 397 | ENDIF |
|---|
| 398 | |
|---|
| 399 | ENDDO |
|---|
| 400 | ENDDO |
|---|
| 401 | |
|---|
| 402 | ENDIF !pscnat |
|---|
| 403 | |
|---|
| 404 | IF ( pscice ) THEN |
|---|
| 405 | |
|---|
| 406 | DO ilon = 1, PLON |
|---|
| 407 | DO iniv = nivbot, nivtop |
|---|
| 408 | |
|---|
| 409 | !cc traitement de la glace |
|---|
| 410 | |
|---|
| 411 | t = MIN(temp(ilon,iniv),273.) |
|---|
| 412 | tice(ilon,iniv) = 2668.70/(10.4310-(log(pw(ilon,iniv))+log(760.0))/log(10.0)) |
|---|
| 413 | |
|---|
| 414 | IF (t .LE. tice(ilon,iniv) ) THEN |
|---|
| 415 | |
|---|
| 416 | lice(ilon,iniv) = 1 |
|---|
| 417 | lnat(ilon,iniv) = 0 |
|---|
| 418 | anat(ilon,iniv) = 0. |
|---|
| 419 | |
|---|
| 420 | zh2oeq = 10.0**(-2668.7/t + 10.4310) |
|---|
| 421 | ph2o = zh2oeq/760.0 |
|---|
| 422 | h2os = MAX(pw(ilon,iniv) - ph2o,0.) & |
|---|
| 423 | *1013.*hnm(ilon,iniv)/ptot(ilon,iniv) |
|---|
| 424 | |
|---|
| 425 | vice(ilon,iniv) = h2os*amu*18./0.928 |
|---|
| 426 | nice(ilon,iniv) = vice(ilon,iniv)/(4./3.*pi*rice(ilon,iniv)**3.) |
|---|
| 427 | aice(ilon,iniv) = nice(ilon,iniv)*4*pi*rice(ilon,iniv)**2 |
|---|
| 428 | |
|---|
| 429 | ! parth2o : part en phase gazeuse |
|---|
| 430 | parth2o(ilon,iniv) = (1.0-MAX(pw(ilon,iniv)-ph2o,0.)/pw(ilon,iniv)) |
|---|
| 431 | parth2o(ilon,iniv) = MAX(parth2o(ilon,iniv),0.) |
|---|
| 432 | parth2o(ilon,iniv) = MIN(parth2o(ilon,iniv),1.) |
|---|
| 433 | |
|---|
| 434 | ! On retire hno3 de la phase gazeuse sous forme de nat |
|---|
| 435 | zmt = -2.7836 - 0.00088*t |
|---|
| 436 | zbt = 38.9855 - 11397.0/t + 0.009179*t |
|---|
| 437 | zhno3eq = 10.0**(zmt*log10(pw(ilon,iniv)& |
|---|
| 438 | *parth2o(ilon,iniv)*1013.*0.75)+ zbt) |
|---|
| 439 | pn = zhno3eq/(0.75*1013.) |
|---|
| 440 | |
|---|
| 441 | ! parthno3s : part en phase gazeuse |
|---|
| 442 | parthno3s(ilon,iniv) = 1.0 - (pn0(ilon,iniv)-pn)& |
|---|
| 443 | /pn0(ilon,iniv) |
|---|
| 444 | parthno3s(ilon,iniv) = MAX(parthno3s(ilon,iniv),0.) |
|---|
| 445 | parthno3s(ilon,iniv) = MIN(parthno3s(ilon,iniv),1.) |
|---|
| 446 | parthno3sl(ilon,iniv) = parthno3s(ilon,iniv) |
|---|
| 447 | |
|---|
| 448 | ENDIF |
|---|
| 449 | |
|---|
| 450 | END DO |
|---|
| 451 | END DO |
|---|
| 452 | |
|---|
| 453 | ENDIF !pscice |
|---|
| 454 | |
|---|
| 455 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 456 | !c traitement des aerosols liquides |
|---|
| 457 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 458 | |
|---|
| 459 | IF (pscliq) THEN |
|---|
| 460 | |
|---|
| 461 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 462 | !c Calcul du nb de moles par unite de volume |
|---|
| 463 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 464 | |
|---|
| 465 | DO ilon = 1, PLON |
|---|
| 466 | DO iniv = nivbot, nivtop |
|---|
| 467 | |
|---|
| 468 | t = MAX(temp(ilon,iniv),MAX(tice(ilon,iniv)-3.0,185.)) |
|---|
| 469 | |
|---|
| 470 | qh2o = h2o(ilon,iniv)*parth2o(ilon,iniv) |
|---|
| 471 | pw(ilon,iniv) = qh2o*ptot(ilon,iniv)/1013.0 |
|---|
| 472 | pw(ilon,iniv) = MIN(pw(ilon,iniv),2.e-3/1013.) |
|---|
| 473 | pw(ilon,iniv) = MAX(pw(ilon,iniv),2.e-5/1013.) |
|---|
| 474 | |
|---|
| 475 | xsb = 1.0/(2.0*(ks(3)+ks(4)/t))*( -ks(1)-ks(2)/t- & |
|---|
| 476 | ((ks(1)+ks(2)/t)**2-4.0*(ks(3)+ks(4)/t) & |
|---|
| 477 | *(ks(5)+ks(6)/t+ks(7)*LOG(t)-LOG(pw(ilon,iniv))))& |
|---|
| 478 | **0.5) |
|---|
| 479 | |
|---|
| 480 | msb = 55.51*xsb/(1.0-xsb) |
|---|
| 481 | ms(ilon,iniv) = msb |
|---|
| 482 | mn(ilon,iniv) = 0.0 |
|---|
| 483 | |
|---|
| 484 | ws(ilon,iniv) = msb*0.098076/(1.0+msb*0.098076) |
|---|
| 485 | wn(ilon,iniv) = 0.0 |
|---|
| 486 | ws(ilon,iniv) = MAX(ws(ilon,iniv), 0.005) !dilution minimale 0.005% |
|---|
| 487 | |
|---|
| 488 | END DO |
|---|
| 489 | END DO |
|---|
| 490 | |
|---|
| 491 | CALL density(ws,wn,temp,dens) |
|---|
| 492 | |
|---|
| 493 | DO ilon = 1, PLON |
|---|
| 494 | DO iniv = nivbot, nivtop |
|---|
| 495 | qh2so4(ilon,iniv) = MAX(VOL_S(ilon,iniv)*1e-12,1.e-20) |
|---|
| 496 | ns(ilon,iniv)=qh2so4(ilon,iniv)*1.e6*ws(ilon,iniv)*dens(ilon,iniv)/98.076 |
|---|
| 497 | END DO |
|---|
| 498 | END DO |
|---|
| 499 | |
|---|
| 500 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 501 | !ccc aerosols liquides binaires h2so4/h2o |
|---|
| 502 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 503 | |
|---|
| 504 | DO ilon = 1, PLON |
|---|
| 505 | DO iniv = nivbot, nivtop |
|---|
| 506 | |
|---|
| 507 | condliq(ilon,iniv) = 1 |
|---|
| 508 | |
|---|
| 509 | t = MAX(temp(ilon,iniv),MAX(tice(ilon,iniv)-3.0,185.)) |
|---|
| 510 | |
|---|
| 511 | qh2o = h2o(ilon,iniv)*parth2o(ilon,iniv) |
|---|
| 512 | pw(ilon,iniv) = qh2o*ptot(ilon,iniv)/1013.0 |
|---|
| 513 | pw(ilon,iniv) = MIN(pw(ilon,iniv),2.e-3/1013.) |
|---|
| 514 | pw(ilon,iniv) = MAX(pw(ilon,iniv),2.e-5/1013.) |
|---|
| 515 | |
|---|
| 516 | qhno3 = qj(ilon,iniv,id_HNO3)*parthno3s(ilon,iniv) |
|---|
| 517 | pn0(ilon,iniv)= ptot(ilon,iniv)/1013.0*qhno3 |
|---|
| 518 | |
|---|
| 519 | a = ks(3) + ks(4)/t |
|---|
| 520 | b = ks(1) + ks(2)/t |
|---|
| 521 | c = ks(5) + ks(6)/t + ks(7)*log(t) - log(pw(ilon,iniv)) |
|---|
| 522 | det = b**2 - 4.*a*c |
|---|
| 523 | |
|---|
| 524 | if (det .gt. 0.) then |
|---|
| 525 | xsb = (-b - sqrt(det))/(2.*a) |
|---|
| 526 | msb = 55.51*xsb/(1.0 - xsb) |
|---|
| 527 | else |
|---|
| 528 | msb = 0. |
|---|
| 529 | condliq(ilon,iniv) = 0 |
|---|
| 530 | end if |
|---|
| 531 | |
|---|
| 532 | ms(ilon,iniv) = msb |
|---|
| 533 | ws(ilon,iniv) = msb*0.098076/(1.0 + msb*0.098076) |
|---|
| 534 | mn(ilon,iniv) = 0.0 |
|---|
| 535 | wn(ilon,iniv) = 0.0 |
|---|
| 536 | |
|---|
| 537 | !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 538 | ! pression saturante de h2so4 |
|---|
| 539 | ! d apres ayers et al., grl, 1980 |
|---|
| 540 | ! et giauque et al., j. am. chem. soc., 1960 |
|---|
| 541 | !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 542 | |
|---|
| 543 | wts = max(41., ws(ilon,iniv)*100.) |
|---|
| 544 | muh2so4 = 4.184*(1.514e4 - 286.*(wts - 40.)& |
|---|
| 545 | + 1.08*(wts - 40.)**2& |
|---|
| 546 | - 3941./(wts - 40.)**0.1) |
|---|
| 547 | |
|---|
| 548 | ! pression saturante en atmosphere |
|---|
| 549 | |
|---|
| 550 | ph2so4 = exp(-10156./t + 16.2590 - muh2so4/(8.314*t)) |
|---|
| 551 | |
|---|
| 552 | if (qh2so4(ilon,iniv)*ptot(ilon,iniv)/1013. .le. ph2so4) then |
|---|
| 553 | ms(ilon,iniv) = 0. |
|---|
| 554 | ws(ilon,iniv) = 0. |
|---|
| 555 | condliq(ilon,iniv) = 0 |
|---|
| 556 | end if |
|---|
| 557 | |
|---|
| 558 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 559 | !ccc aerosols liquides ternaires hno3/h2so4/h2o |
|---|
| 560 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 561 | |
|---|
| 562 | IF ( condliq(ilon,iniv) .eq. 1 .and. t .lt. 215.) THEN |
|---|
| 563 | |
|---|
| 564 | tf = t |
|---|
| 565 | tt = r*tf*ns(ilon,iniv) |
|---|
| 566 | tr = 1.0e4/tf-43.4782608 |
|---|
| 567 | pr = LOG(pw(ilon,iniv))+18.4 |
|---|
| 568 | |
|---|
| 569 | hsb = qs(1)+qs(2)*tr**2+(qs(3)+qs(4)*tr+qs(5)*tr**2 & |
|---|
| 570 | +qs(6)*tr**3)*pr + (qs(7)+qs(8)*tr+qs(9)*tr**2) & |
|---|
| 571 | *pr**2+qs(10)*tr*pr**3 |
|---|
| 572 | hsb = EXP(hsb) |
|---|
| 573 | xnb = 1.0/(2.0*(kn(3)+kn(4)/tf))*(-kn(1)-kn(2)/tf- & |
|---|
| 574 | ((kn(1)+kn(2)/tf)**2-4.0*(kn(3)+kn(4)/tf) & |
|---|
| 575 | *(kn(5)+kn(6)/tf +kn(7)*LOG(tf) & |
|---|
| 576 | -LOG(pw(ilon,iniv))))**0.5) |
|---|
| 577 | mnb = 55.51*xnb/(1.0-xnb) |
|---|
| 578 | hnb = qn(1)+qn(2)*tr**2+(qn(3)+qn(4)*tr+qn(5)*tr**2 & |
|---|
| 579 | +qn(6)*tr**3)*pr + (qn(7)+qn(8)*tr+qn(9)*tr**2) & |
|---|
| 580 | *pr**2+qn(10)*tr*pr**3 |
|---|
| 581 | hnb = EXP(hnb) |
|---|
| 582 | |
|---|
| 583 | a = (tt*hnb*mnb**2 - tt*hsb*mnb*msb - 2.0*mnb**2*msb & |
|---|
| 584 | + mnb*msb**2 + hnb*mnb*msb*pn0(ilon,iniv) & |
|---|
| 585 | - hsb*msb**2*pn0(ilon,iniv))/(mnb**2 - mnb*msb) |
|---|
| 586 | b = msb*(-2.0*tt*hnb*mnb+tt*hsb*msb+mnb*msb & |
|---|
| 587 | -hnb*msb*pn0(ilon,iniv))/(mnb-msb) |
|---|
| 588 | c = (tt*hnb*mnb*msb**2)/(mnb - msb) |
|---|
| 589 | |
|---|
| 590 | phi = ATAN(SQRT(MAX(4.0*(a**2-3.0*b)**3-(-2.0*a**3+9.0*a*b & |
|---|
| 591 | -27.0*c)**2,0.))/(-2.0*a**3+9.0*a*b-27.0*c) ) |
|---|
| 592 | IF (phi .LT. 0.) THEN |
|---|
| 593 | phi = phi + pi |
|---|
| 594 | ENDIF |
|---|
| 595 | |
|---|
| 596 | msf = -1.0/3.0*(a+2.0*SQRT(a**2-3.0*b)*COS((pi+phi)/3.0)) |
|---|
| 597 | msf = MAX(msf,0.) |
|---|
| 598 | |
|---|
| 599 | mnf = mnb*(1.0-msf/msb) |
|---|
| 600 | mnf = MAX(mnf,0.) |
|---|
| 601 | pn = mnf/(hnb*mnf/(mnf+msf)+hsb*msf/(mnf+msf)) |
|---|
| 602 | wsf = msf*0.098076/(1.0+msf*0.098076+mnf*0.063012) |
|---|
| 603 | wnf = mnf*0.063012/(1.0+msf*0.098076+mnf*0.063012) |
|---|
| 604 | |
|---|
| 605 | ms(ilon,iniv) = msf |
|---|
| 606 | mn(ilon,iniv) = mnf |
|---|
| 607 | ws(ilon,iniv) = wsf |
|---|
| 608 | wn(ilon,iniv) = wnf |
|---|
| 609 | |
|---|
| 610 | parthno3l(ilon,iniv) = (1.0-(pn0(ilon,iniv)-pn)/pn0(ilon,iniv)) |
|---|
| 611 | parthno3l(ilon,iniv) = MIN(parthno3l(ilon,iniv),1.) |
|---|
| 612 | |
|---|
| 613 | ENDIF |
|---|
| 614 | |
|---|
| 615 | END DO |
|---|
| 616 | END DO |
|---|
| 617 | |
|---|
| 618 | !c dilution minimale de h2so4: 0.5% |
|---|
| 619 | WHERE ( ws < .005 ) ws = 0.005 |
|---|
| 620 | |
|---|
| 621 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 622 | !ccc densite des aerosols ternaires (g/cm3) |
|---|
| 623 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 624 | CALL density(ws,wn,temp,dens) |
|---|
| 625 | |
|---|
| 626 | DO ilon = 1, PLON |
|---|
| 627 | DO iniv = nivbot, nivtop |
|---|
| 628 | !c |
|---|
| 629 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 630 | !ccc total specific liquid aerosol volume (dimensionless) |
|---|
| 631 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 632 | !c |
|---|
| 633 | IF ( condliq(ilon,iniv) .EQ. 1) THEN |
|---|
| 634 | vliq = ns(ilon,iniv)*98.076e-6/ws(ilon,iniv)/dens(ilon,iniv) |
|---|
| 635 | !c |
|---|
| 636 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 637 | !ccc liquid aerosol radius and surface area |
|---|
| 638 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 639 | !c |
|---|
| 640 | !c rmean=mean radius of liquid droplets (units = cm) |
|---|
| 641 | !c this varies with the total liquid volume per unit volume of air |
|---|
| 642 | !c (called vliq in main program). rmean can be calculated |
|---|
| 643 | !c from vliq using |
|---|
| 644 | !c |
|---|
| 645 | !c rmean=rmode*exp(0.5*(log(sigma))**2), where |
|---|
| 646 | !c |
|---|
| 647 | !c rmode=(3.0*vliq/(4.*pi*ndrop)* |
|---|
| 648 | !c exp(-9.0/2.0*(log(sigma)**2)))**(1./3.) |
|---|
| 649 | !c |
|---|
| 650 | !c ndrop=number of droplets per cm3 air |
|---|
| 651 | !c sigma=width of lognormal (use sigma=1.8). |
|---|
| 652 | !c |
|---|
| 653 | rmode=(3.0*vliq/(4.*pi*ndrop(ilon,iniv))*EXP(-9.0/2.0*(LOG(sigma)**2)))**(1./3.) |
|---|
| 654 | rmean(ilon,iniv)=rmode*exp(0.5*(log(sigma))**2) |
|---|
| 655 | |
|---|
| 656 | aliq(ilon,iniv)=4.0*pi*rmode**2*ndrop(ilon,iniv) & |
|---|
| 657 | *EXP(2.0*(LOG(sigma))**2) |
|---|
| 658 | |
|---|
| 659 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 660 | !ccc the solubility of hcl |
|---|
| 661 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 662 | !c |
|---|
| 663 | !c hhcl est calcule directement a partir de jpl 2000, |
|---|
| 664 | !c lui-meme utilisant carslaw et al., j. phys. chem., 1995. |
|---|
| 665 | !c l'ensemble des calculs suppose un aerosol binaire: |
|---|
| 666 | !c la molalite de h2so4 (wt) est donc recalculee dans |
|---|
| 667 | !c ces conditions. |
|---|
| 668 | !c |
|---|
| 669 | t = MAX(temp(ilon,iniv),MAX(tice(ilon,iniv)-3.0,185.)) |
|---|
| 670 | !c |
|---|
| 671 | !c ph2o0 : saturation water vapor pressure, hpa |
|---|
| 672 | !c |
|---|
| 673 | ph2o0 = EXP(18.452406985 - 3505.1578807/t & |
|---|
| 674 | - 330918.55082/t**2 + 12725068.262/t**3) |
|---|
| 675 | !c |
|---|
| 676 | !c aw : water activity |
|---|
| 677 | !c |
|---|
| 678 | aw(ilon,iniv) = pw(ilon,iniv)*1013./ph2o0 |
|---|
| 679 | aw(ilon,iniv) = MIN(aw(ilon,iniv), 1.) |
|---|
| 680 | |
|---|
| 681 | IF (aw(ilon,iniv) .LE. 0.05) THEN |
|---|
| 682 | y1 = a1(1)*aw(ilon,iniv)**b1(1) + c1(1)*aw(ilon,iniv) + d1(1) |
|---|
| 683 | y2 = a2(1)*aw(ilon,iniv)**b2(1) + c2(1)*aw(ilon,iniv) + d2(1) |
|---|
| 684 | else if (aw(ilon,iniv) .lt. 0.85) then |
|---|
| 685 | y1 = a1(2)*aw(ilon,iniv)**b1(2) + c1(2)*aw(ilon,iniv) + d1(2) |
|---|
| 686 | y2 = a2(2)*aw(ilon,iniv)**b2(2) + c2(2)*aw(ilon,iniv) + d2(2) |
|---|
| 687 | else |
|---|
| 688 | y1 = a1(3)*aw(ilon,iniv)**b1(3) + c1(3)*aw(ilon,iniv) + d1(3) |
|---|
| 689 | y2 = a2(3)*aw(ilon,iniv)**b2(3) + c2(3)*aw(ilon,iniv) + d2(3) |
|---|
| 690 | endif |
|---|
| 691 | !c |
|---|
| 692 | !c m: h2so4 molality, mol/kg |
|---|
| 693 | !c |
|---|
| 694 | m = y1 + (t - 190.)*(y2 - y1)/70. |
|---|
| 695 | !c |
|---|
| 696 | !c wt: h2so4 weight percentage, % |
|---|
| 697 | !c |
|---|
| 698 | wt(ilon,iniv) = 9800.*m/(98.*m + 1000.) |
|---|
| 699 | wt(ilon,iniv) = max(wt(ilon,iniv), 0.5) |
|---|
| 700 | |
|---|
| 701 | !c |
|---|
| 702 | !c mh2so4: h2so4 molarity, mol l-1 |
|---|
| 703 | !c |
|---|
| 704 | mh2so4(ilon,iniv) = dens(ilon,iniv)*wt(ilon,iniv)/9.8 |
|---|
| 705 | |
|---|
| 706 | bigx = wt(ilon,iniv)/(wt(ilon,iniv) & |
|---|
| 707 | + (100. - wt(ilon,iniv))*98./18.) |
|---|
| 708 | |
|---|
| 709 | hhcl(ilon,iniv) = (0.094 - 0.61*bigx + 1.2*bigx**2) & |
|---|
| 710 | *exp(-8.68 + (8515. - 10718.*(bigx**0.7))/t) |
|---|
| 711 | !c |
|---|
| 712 | !c phcl1 : pression de hcl restant en phase gazeuse |
|---|
| 713 | !c ns/ms : kilogrammes (litres) d'eau par metre cube |
|---|
| 714 | !c |
|---|
| 715 | phcl1 = (phcl0(ilon,iniv)/(r*t)) & |
|---|
| 716 | /(1./(r*t) + hhcl(ilon,iniv)*ns(ilon,iniv)/ms(ilon,iniv)) |
|---|
| 717 | |
|---|
| 718 | parthcl(ilon,iniv) = 1.0-(phcl0(ilon,iniv)-phcl1)/phcl0(ilon,iniv) |
|---|
| 719 | |
|---|
| 720 | !c |
|---|
| 721 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 722 | !ccc the solubility of clono2 |
|---|
| 723 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 724 | !c |
|---|
| 725 | !c on peut negliger la diminution de la concentration |
|---|
| 726 | !c en phase gazeuse. |
|---|
| 727 | !c |
|---|
| 728 | !c sclono2 : setchenow coefficient m-1 |
|---|
| 729 | !c |
|---|
| 730 | sclono2 = 0.306 + 24./t |
|---|
| 731 | !c |
|---|
| 732 | !c hclono2 : shi et al., m atm-1 |
|---|
| 733 | !c |
|---|
| 734 | hclono2(ilon,iniv) = 1.6e-6*exp(4710./t) & |
|---|
| 735 | *exp(-sclono2*mh2so4(ilon,iniv)) |
|---|
| 736 | !c |
|---|
| 737 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 738 | !ccc the solubility of hocl |
|---|
| 739 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 740 | !c |
|---|
| 741 | !c hhocl d'apres jpl 2000 (shi et al.) |
|---|
| 742 | !c on peut negliger la diminution de la concentration |
|---|
| 743 | !c en phase gazeuse. |
|---|
| 744 | !c |
|---|
| 745 | !c shocl : setchenow coefficient m-1 |
|---|
| 746 | !c |
|---|
| 747 | shocl = 0.0776 + 59.18/t |
|---|
| 748 | !c |
|---|
| 749 | !c hhocl : shi et al. |
|---|
| 750 | !c |
|---|
| 751 | hhocl(ilon,iniv) = 1.91e-6*exp(5862.4/t) & |
|---|
| 752 | *exp(-shocl*mh2so4(ilon,iniv)) |
|---|
| 753 | !c |
|---|
| 754 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 755 | !ccc the solubility of hobr |
|---|
| 756 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 757 | !c |
|---|
| 758 | !c hhobr d'apres waschewsky and abbatt, j. phys. chem. a, |
|---|
| 759 | !c 5312-5320, 1999. |
|---|
| 760 | !c on peut negliger la diminution de la concentration |
|---|
| 761 | !c en phase gazeuse. |
|---|
| 762 | !c |
|---|
| 763 | hhobr(ilon,iniv) = 4.6e-4*exp(4.5e3/t) |
|---|
| 764 | !c |
|---|
| 765 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 766 | !ccc the solubility of hbr |
|---|
| 767 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 768 | !c |
|---|
| 769 | !c hhbr d'apres kleffmann et al., j. phys. chem. a, |
|---|
| 770 | !c 8489-8495, 2000. attention a l'erreur typographique |
|---|
| 771 | !c contenue dans la version imprimee de l'article |
|---|
| 772 | !c (m3 = - 4.445 imprime au lieu de m3 = + 4.445) |
|---|
| 773 | !c |
|---|
| 774 | mh2so4c = -1.977e-4*wt(ilon,iniv)**2. & |
|---|
| 775 | -2.096e-2*wt(ilon,iniv) + 4.445 |
|---|
| 776 | |
|---|
| 777 | bh2so4 = -8.979e-5*wt(ilon,iniv)**2. & |
|---|
| 778 | +2.141e-2*wt(ilon,iniv) - 6.067 |
|---|
| 779 | |
|---|
| 780 | hhbr(ilon,iniv) = 10.**(mh2so4c*1000./t + bh2so4) |
|---|
| 781 | |
|---|
| 782 | ENDIF !condliq |
|---|
| 783 | |
|---|
| 784 | ENDDO |
|---|
| 785 | ENDDO |
|---|
| 786 | |
|---|
| 787 | ENDIF !pscliq |
|---|
| 788 | |
|---|
| 789 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 790 | !c surfaces des PSC |
|---|
| 791 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 792 | |
|---|
| 793 | do ilon = 1,PLON |
|---|
| 794 | do iniv = nivbot,nivtop |
|---|
| 795 | |
|---|
| 796 | surfarealiq(ilon,iniv) = aliq(ilon,iniv) |
|---|
| 797 | surfareanat(ilon,iniv) = anat(ilon,iniv)*lnat(ilon,iniv) |
|---|
| 798 | surfareaice(ilon,iniv) = aice(ilon,iniv)*lice(ilon,iniv) |
|---|
| 799 | |
|---|
| 800 | end do |
|---|
| 801 | end do |
|---|
| 802 | |
|---|
| 803 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 804 | !c partitions |
|---|
| 805 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 806 | |
|---|
| 807 | do ilon = 1, PLON |
|---|
| 808 | do iniv = nivbot,nivtop |
|---|
| 809 | |
|---|
| 810 | parthno3(ilon,iniv) = parthno3l(ilon,iniv) * parthno3s(ilon,iniv) |
|---|
| 811 | parthno3(ilon,iniv) = MIN ( parthno3(ilon,iniv) , 1. ) |
|---|
| 812 | |
|---|
| 813 | end do |
|---|
| 814 | end do |
|---|
| 815 | |
|---|
| 816 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 817 | !c sedimentation |
|---|
| 818 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 819 | |
|---|
| 820 | rhoice = 0.928*1000.0 |
|---|
| 821 | rhonat = 1.350*1000.0 |
|---|
| 822 | |
|---|
| 823 | do ilon = 1, PLON |
|---|
| 824 | do iniv = nivbot,nivtop |
|---|
| 825 | |
|---|
| 826 | IF ( lice(ilon,iniv) .EQ. 1 .OR. lnat(ilon,iniv) .EQ. 1 ) THEN |
|---|
| 827 | |
|---|
| 828 | t = temp(ilon,iniv) |
|---|
| 829 | |
|---|
| 830 | rpart=(lice(ilon,iniv)*rice(ilon,iniv) & |
|---|
| 831 | +lnat(ilon,iniv)*rnat(ilon,iniv)) *1.e-2 |
|---|
| 832 | rpart=min(rpart,535.e-6) |
|---|
| 833 | |
|---|
| 834 | if (rpart .ge. 5.e-7) then |
|---|
| 835 | |
|---|
| 836 | rhopart=lice(ilon,iniv)*rhoice +lnat(ilon,iniv)*rhonat |
|---|
| 837 | rhoair=(hnm(ilon,iniv)*28.97/6.022e23) *1000.0 |
|---|
| 838 | |
|---|
| 839 | tcel = t - 273.15 |
|---|
| 840 | netha = (1.718 + 0.0049*tcel -1.2e-5*tcel*tcel)*1.e-5 |
|---|
| 841 | |
|---|
| 842 | cdnre2 = 32.0*(rpart)**3*(rhopart-rhoair) & |
|---|
| 843 | *rhoair*9.80665 /(3.0*(netha)**2) |
|---|
| 844 | |
|---|
| 845 | !c pruppacher and klett 1998 |
|---|
| 846 | ynre = bnre(0)+bnre(1)*log(cdnre2) & |
|---|
| 847 | +bnre(2)*(log(cdnre2))**2 & |
|---|
| 848 | +bnre(3)*(log(cdnre2))**3 & |
|---|
| 849 | +bnre(4)*(log(cdnre2))**4 & |
|---|
| 850 | +bnre(5)*(log(cdnre2))**5 & |
|---|
| 851 | +bnre(6)*(log(cdnre2))**6 |
|---|
| 852 | nre = exp(ynre) |
|---|
| 853 | |
|---|
| 854 | if (nre .ge. 1.e-2) then |
|---|
| 855 | !c |
|---|
| 856 | !c big particles |
|---|
| 857 | !c (1.e-2 < reynolds number < 300) |
|---|
| 858 | !c (10.0 < radius < 535 microns) |
|---|
| 859 | !c |
|---|
| 860 | vsed(ilon,iniv) = netha*nre & |
|---|
| 861 | /(2*rhoair*rpart) |
|---|
| 862 | |
|---|
| 863 | else |
|---|
| 864 | !c |
|---|
| 865 | !c small particules |
|---|
| 866 | !c (1.e-6 < reynolds number < 1.e-2) |
|---|
| 867 | !c (0.5 < radius < 10 microns) |
|---|
| 868 | !c |
|---|
| 869 | !c lambda0 : libre parcours moyen (m) a 1013.25 hpa |
|---|
| 870 | !c et 293.15k. voir pruppacher and klett |
|---|
| 871 | !c page 416. |
|---|
| 872 | !c |
|---|
| 873 | lambda0 = 6.6e-8 |
|---|
| 874 | lambda = lambda0 & |
|---|
| 875 | *(1013.25/ptot(ilon,iniv)) & |
|---|
| 876 | *(t/293.15) |
|---|
| 877 | |
|---|
| 878 | us = 2.0*rpart**2*9.80665 & |
|---|
| 879 | *(rhopart-rhoair)/(9.0*netha) |
|---|
| 880 | |
|---|
| 881 | !c |
|---|
| 882 | !ccc solid particles: |
|---|
| 883 | !c |
|---|
| 884 | !c niels larsen - 1991 |
|---|
| 885 | !c columnar or elipsoid shape |
|---|
| 886 | !c with a shape factor of 1.2 (shapek) |
|---|
| 887 | !c |
|---|
| 888 | shapek = 1.2 |
|---|
| 889 | |
|---|
| 890 | vsedsol = us/shapek & |
|---|
| 891 | *(1.0+1.246*lambda/rpart & |
|---|
| 892 | +0.42*lambda/rpart & |
|---|
| 893 | *exp(-0.87*rpart/lambda)) |
|---|
| 894 | |
|---|
| 895 | vsed(ilon,iniv) = vsedsol |
|---|
| 896 | endif |
|---|
| 897 | endif |
|---|
| 898 | |
|---|
| 899 | ENDIF |
|---|
| 900 | |
|---|
| 901 | end do |
|---|
| 902 | end do |
|---|
| 903 | |
|---|
| 904 | !c |
|---|
| 905 | !c only the parameters that are needed in 'hetero' are |
|---|
| 906 | !c included in the variable list. output variables should |
|---|
| 907 | !c be included in the list as needed! |
|---|
| 908 | !c |
|---|
| 909 | call hetero(aw,temp,ptot,hnm,qj,ws, & |
|---|
| 910 | hhcl,hhocl,hhobr,hclono2,hhbr, & |
|---|
| 911 | rice,nice,rnat,nnat, & |
|---|
| 912 | aliq,rmean,lnat,lice,condliq, & |
|---|
| 913 | k1,k2,k3,k4,k5,k6,k7,k8,k9,wt,h2o) |
|---|
| 914 | |
|---|
| 915 | END SUBROUTINE analytic |
|---|
| 916 | |
|---|
| 917 | SUBROUTINE HETERO(aw,temp,ptot,hnm,qj,ws, & |
|---|
| 918 | hhcl,hhocl,hhobr,hclono2,hhbr, & |
|---|
| 919 | rice,nice,rnat,nnat, & |
|---|
| 920 | aliq,rmean,lnat,lice,condliq, & |
|---|
| 921 | k1,k2,k3,k4,k5,k6,k7,k8,k9,wt,h2o) |
|---|
| 922 | |
|---|
| 923 | !C ROUTINE TO CALCULATE UPTAKE COEFFICIENTS (GAMMA VALUES). |
|---|
| 924 | !C GAMMA VALUES ARE INDICATED BY VARIABLES WITH PREFIX 'G', FOR |
|---|
| 925 | !C EXAMPLE GHOCLHCL IS THE GAMMA VALUE OF HOCL DUE TO REACTION WITH |
|---|
| 926 | !C HCL IN THE DROPLETS. |
|---|
| 927 | !C FROM THE GAMMA VALUES, SECOND ORDER RATE CONSTANTS ARE CALCULATED. |
|---|
| 928 | !C THESE HAVE THE PREFIX 'R' AND HAVE UNITS CM3 MOLECULE-1 S-1. FOR |
|---|
| 929 | !C EXAMPLE, THE LOSS OF CLNO3 AND HCL DUE TO THE HETEROGENEOUS REACTION |
|---|
| 930 | !C CLNO3+HCL -> CL2+HNO3 IS D(CLNO3)/DT (UNITS MOLECULE CM-3 S-1) = |
|---|
| 931 | !C -RCLNO3HCL.[CLNO3].[HCL], WHERE [CLNO3] AND [HCL] ARE THE |
|---|
| 932 | !C ****TOTAL**** AMOUNTS OF THESE SPECIES IN UNITS MOLECULE CM-3. |
|---|
| 933 | |
|---|
| 934 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 935 | !c F. LEFEVRE (SA/CNRS) from REPROBUS |
|---|
| 936 | !c |
|---|
| 937 | !c Adapted to LMDz-INCA: |
|---|
| 938 | !c |
|---|
| 939 | !c L. JOURDAIN (SA) 2004 |
|---|
| 940 | !c D. HAUGLUSTAINE (LSCE) 06/2005 |
|---|
| 941 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 942 | |
|---|
| 943 | use SPECIES_NAMES |
|---|
| 944 | |
|---|
| 945 | INTEGER :: i, ilon, iniv |
|---|
| 946 | |
|---|
| 947 | REAL, PARAMETER :: PI=3.1415 |
|---|
| 948 | REAL, PARAMETER :: ZERO = 1.e-20 |
|---|
| 949 | |
|---|
| 950 | REAL :: t0, t, dt |
|---|
| 951 | REAL :: zh2o, zhcl, zhocl, zhobr |
|---|
| 952 | REAL :: biga, eta, ah |
|---|
| 953 | REAL :: dclono2, cclono2 |
|---|
| 954 | REAL :: gbh2o, grxn, gbhcl, fhcl, gs, gps, gbhclp, gb |
|---|
| 955 | REAL :: gclono2, gclono2hcl, gclono2h2o |
|---|
| 956 | REAL :: rclono2h2o, rclono2hcl |
|---|
| 957 | REAL :: chocl, dhocl, fhocl, ghoclhcl, rhoclhcl |
|---|
| 958 | REAL :: ratio |
|---|
| 959 | REAL :: dhobr, fhobr, chobr |
|---|
| 960 | REAL :: ghobrhcl, rhobrhcl |
|---|
| 961 | REAL :: dhcl, chcl |
|---|
| 962 | REAL :: cbar |
|---|
| 963 | REAL :: gxrn |
|---|
| 964 | REAL :: gn2o5h2o, rn2o5h2o |
|---|
| 965 | REAL :: gbrno3h2o, rbrno3h2o |
|---|
| 966 | REAL :: gclno3hclnat, gclno3h2onat |
|---|
| 967 | REAL :: gn2o5h2onat, gn2o5hclnat, ghoclhclnat, gbrno3h2onat, ghobrhclnat, ghobrhbrnat, ghoclhbrnat |
|---|
| 968 | REAL :: rclno3h2onat, rclno3hclnat, rn2o5h2onat, rbrno3h2onat, rn2o5hclnat, rhoclhclnat, rhobrhclnat |
|---|
| 969 | REAL :: gclno3hclice, gclno3h2oice, gn2o5h2oice, gn2o5hclice, ghoclhclice, gbrno3h2oice, ghobrhclice, ghobrhbrice, ghoclhbrice |
|---|
| 970 | REAL :: rclno3hclice, rclno3h2oice, rn2o5h2oice, rn2o5hclice, rhoclhclice, rbrno3h2oice, rhobrhclice, rhobrhbrice |
|---|
| 971 | |
|---|
| 972 | REAL TEMP(PLON,PLEV), PTOT(PLON,PLEV), HNM(PLON,PLEV) |
|---|
| 973 | REAL CK(PLON,PLEV,NBCON), qj(PLON,PLEV,NBCON) |
|---|
| 974 | |
|---|
| 975 | REAL, DIMENSION(PLON,PLEV) :: k1 |
|---|
| 976 | REAL, DIMENSION(PLON,PLEV) :: k2 |
|---|
| 977 | REAL, DIMENSION(PLON,PLEV) :: k3 |
|---|
| 978 | REAL, DIMENSION(PLON,PLEV) :: k4 |
|---|
| 979 | REAL, DIMENSION(PLON,PLEV) :: k5 |
|---|
| 980 | REAL, DIMENSION(PLON,PLEV) :: k6 |
|---|
| 981 | REAL, DIMENSION(PLON,PLEV) :: k7 |
|---|
| 982 | REAL, DIMENSION(PLON,PLEV) :: k8 |
|---|
| 983 | REAL, DIMENSION(PLON,PLEV) :: k9 |
|---|
| 984 | |
|---|
| 985 | REAL, DIMENSION(PLON,PLEV) :: k1l |
|---|
| 986 | REAL, DIMENSION(PLON,PLEV) :: k2l |
|---|
| 987 | REAL, DIMENSION(PLON,PLEV) :: k3l |
|---|
| 988 | REAL, DIMENSION(PLON,PLEV) :: k4l |
|---|
| 989 | REAL, DIMENSION(PLON,PLEV) :: k5l |
|---|
| 990 | REAL, DIMENSION(PLON,PLEV) :: k6l |
|---|
| 991 | REAL, DIMENSION(PLON,PLEV) :: k7l |
|---|
| 992 | REAL, DIMENSION(PLON,PLEV) :: k8l |
|---|
| 993 | REAL, DIMENSION(PLON,PLEV) :: k9l |
|---|
| 994 | |
|---|
| 995 | REAL, DIMENSION(PLON,PLEV) :: k1s |
|---|
| 996 | REAL, DIMENSION(PLON,PLEV) :: k2s |
|---|
| 997 | REAL, DIMENSION(PLON,PLEV) :: k3s |
|---|
| 998 | REAL, DIMENSION(PLON,PLEV) :: k4s |
|---|
| 999 | REAL, DIMENSION(PLON,PLEV) :: k5s |
|---|
| 1000 | REAL, DIMENSION(PLON,PLEV) :: k6s |
|---|
| 1001 | REAL, DIMENSION(PLON,PLEV) :: k7s |
|---|
| 1002 | REAL, DIMENSION(PLON,PLEV) :: k8s |
|---|
| 1003 | REAL, DIMENSION(PLON,PLEV) :: k9s |
|---|
| 1004 | |
|---|
| 1005 | REAL H2O(PLON,PLEV), H2OCONC(PLON,PLEV) |
|---|
| 1006 | |
|---|
| 1007 | REAL WS(PLON,PLEV) |
|---|
| 1008 | REAL HHCL(PLON,PLEV), HHOCL(PLON,PLEV) |
|---|
| 1009 | REAL HHOBR(PLON,PLEV), HCLONO2(PLON,PLEV) |
|---|
| 1010 | REAL HHBR(PLON,PLEV) |
|---|
| 1011 | REAL RICE(PLON,PLEV) |
|---|
| 1012 | REAL RNAT(PLON,PLEV), ALIQ(PLON,PLEV) |
|---|
| 1013 | REAL RMEAN(PLON,PLEV) |
|---|
| 1014 | REAL NNAT(PLON,PLEV), NICE(PLON,PLEV) |
|---|
| 1015 | |
|---|
| 1016 | INTEGER CONDLIQ(PLON,PLEV) |
|---|
| 1017 | INTEGER LNAT(PLON,PLEV), LICE(PLON,PLEV) |
|---|
| 1018 | !c |
|---|
| 1019 | !c additions pour jpl 2000 |
|---|
| 1020 | !c |
|---|
| 1021 | real aw(PLON,PLEV) |
|---|
| 1022 | real wt(PLON,PLEV) |
|---|
| 1023 | real kh, kh2o, khydr |
|---|
| 1024 | real phcl, mhcl, khcl |
|---|
| 1025 | real pclono2, lclono2, fclono2 |
|---|
| 1026 | real phobr, mhobr, k1hobrhcl, k2hobrhcl |
|---|
| 1027 | real phbr, mhbr, lhbr, k1hoclhbr, k2hoclhbr |
|---|
| 1028 | real khoclhcl, lhocl, phocl, mhocl |
|---|
| 1029 | real lhobr, lhcl |
|---|
| 1030 | |
|---|
| 1031 | |
|---|
| 1032 | ! Initialisation |
|---|
| 1033 | |
|---|
| 1034 | k1l(:,:) = 0. |
|---|
| 1035 | k2l(:,:) = 0. |
|---|
| 1036 | k3l(:,:) = 0. |
|---|
| 1037 | k4l(:,:) = 0. |
|---|
| 1038 | k5l(:,:) = 0. |
|---|
| 1039 | k6l(:,:) = 0. |
|---|
| 1040 | k7l(:,:) = 0. |
|---|
| 1041 | k8l(:,:) = 0. |
|---|
| 1042 | k9l(:,:) = 0. |
|---|
| 1043 | k1s(:,:) = 0. |
|---|
| 1044 | k2s(:,:) = 0. |
|---|
| 1045 | k3s(:,:) = 0. |
|---|
| 1046 | k4s(:,:) = 0. |
|---|
| 1047 | k5s(:,:) = 0. |
|---|
| 1048 | k6s(:,:) = 0. |
|---|
| 1049 | k7s(:,:) = 0. |
|---|
| 1050 | k8s(:,:) = 0. |
|---|
| 1051 | k9s(:,:) = 0. |
|---|
| 1052 | |
|---|
| 1053 | h2oconc(:,:) = h2o(:,:)*hnm(:,:) |
|---|
| 1054 | DO i = 1, nbcon |
|---|
| 1055 | ck(:,:,i) = qj(:,:,i)*hnm(:,:) |
|---|
| 1056 | ENDDO |
|---|
| 1057 | |
|---|
| 1058 | !C |
|---|
| 1059 | !C ****************************************************** |
|---|
| 1060 | !C THE FOLLOWING PARAMETERS ARE NEEDED IN THIS ROUTINE!! |
|---|
| 1061 | !C ****************************************************** |
|---|
| 1062 | !C |
|---|
| 1063 | !C RNAT=MEAN RADIUS OF NAT PARTICLES (CM) |
|---|
| 1064 | !C NNAT=NUMBER OF NAT PARTICLES PER CM3 AIR |
|---|
| 1065 | !C RICE=MEAN RADIUS OF ICE PARTICLES (CM) |
|---|
| 1066 | !C NICE=NUMBER OF ICE PARTICLES PER CM3 AIR |
|---|
| 1067 | !C RMEAN=MEAN RADIUS OF LIQUID AEROSOL DISTRIBUTION (CM) |
|---|
| 1068 | !C ALIQ=TOTAL AREA OF LIQUID AEROSOLS (CM2 CM-3 AIR) |
|---|
| 1069 | |
|---|
| 1070 | DO ILON = 1, PLON |
|---|
| 1071 | DO INIV = nivbot, nivtop |
|---|
| 1072 | |
|---|
| 1073 | IF ( condliq(ilon,iniv) .EQ. 1 ) THEN |
|---|
| 1074 | |
|---|
| 1075 | T = TEMP(ILON,INIV) |
|---|
| 1076 | !C |
|---|
| 1077 | !C Z VARIABLES ARE SET TO AVOID DIVISION BY ZERO |
|---|
| 1078 | !C |
|---|
| 1079 | ZH2O = MAX(h2oconc(ILON,INIV),1.E-08*HNM(ILON,INIV)) |
|---|
| 1080 | ZHCL = MAX(CK(ILON,INIV,id_HCl),1.E-13*HNM(ILON,INIV)) |
|---|
| 1081 | ZHOCL = MAX(CK(ILON,INIV,id_HOCl),1.E-13*HNM(ILON,INIV)) |
|---|
| 1082 | ZHOBR = MAX(CK(ILON,INIV,id_HOBr),1.E-13*HNM(ILON,INIV)) |
|---|
| 1083 | !c |
|---|
| 1084 | !c pressions de hcl, clono2, hobr, et hbr (atm) |
|---|
| 1085 | !c |
|---|
| 1086 | phcl = (ck(ilon,iniv,id_HCl)/hnm(ilon,iniv))*ptot(ilon,iniv)/1013. |
|---|
| 1087 | pclono2 = (ck(ilon,iniv,id_ClONO2)/hnm(ilon,iniv))*ptot(ilon,iniv)/1013. |
|---|
| 1088 | phocl = (ck(ilon,iniv,id_HOCl)/hnm(ilon,iniv))*ptot(ilon,iniv)/1013. |
|---|
| 1089 | phobr = (ck(ilon,iniv,id_HOBr)/hnm(ilon,iniv))*ptot(ilon,iniv)/1013. |
|---|
| 1090 | phbr = (ck(ilon,iniv,id_HBr)/hnm(ilon,iniv))*ptot(ilon,iniv)/1013. |
|---|
| 1091 | !c |
|---|
| 1092 | !C ===================================================================== |
|---|
| 1093 | !C ====================== CLONO2+HCL -> CL2+HNO3 ======================= |
|---|
| 1094 | !C ====================== CLONO2+H2O -> HOCL+HNO3 ====================== |
|---|
| 1095 | !C ========================= ON LIQUID AEROSOL ========================= |
|---|
| 1096 | !C ===================================================================== |
|---|
| 1097 | !c |
|---|
| 1098 | !c d'apres jpl 2000. On ne tient pas compte de la deposition |
|---|
| 1099 | !c de hno3 (solutions ternaires) sur les probabilites de reaction. |
|---|
| 1100 | !c |
|---|
| 1101 | t0 = 144.11 + 0.166*wt(ilon,iniv) - 0.015*wt(ilon,iniv)**2 & |
|---|
| 1102 | + 2.18e-4*wt(ilon,iniv)**3 |
|---|
| 1103 | biga = 169.5 + 5.18*wt(ilon,iniv) - 0.0825*wt(ilon,iniv)**2 & |
|---|
| 1104 | + 3.27e-3*wt(ilon,iniv)**3 |
|---|
| 1105 | !c |
|---|
| 1106 | !c eta : viscosity of h2so4 solution, cp |
|---|
| 1107 | !c |
|---|
| 1108 | ! Temperature check |
|---|
| 1109 | dt = t - t0 |
|---|
| 1110 | IF ( dt < 1.) THEN |
|---|
| 1111 | dt = 1. |
|---|
| 1112 | PRINT *, 'Warning : temperature problem in hetchem (ILON, T, T0) : ', ilon, t, t0 |
|---|
| 1113 | END IF |
|---|
| 1114 | eta = biga*t**(-1.43)*exp(448./dt) |
|---|
| 1115 | !c |
|---|
| 1116 | !c ah : acid activity in molarity |
|---|
| 1117 | !c |
|---|
| 1118 | ah = exp(60.51 - 0.095*wt(ilon,iniv) + 0.0077*wt(ilon,iniv)**2 & |
|---|
| 1119 | -1.61e-5*wt(ilon,iniv)**3 - (1.76 + & |
|---|
| 1120 | 2.52e-4*wt(ilon,iniv)**2)*sqrt(t) & |
|---|
| 1121 | + (-805.89 + 253.05*wt(ilon,iniv)**0.076)/sqrt(t)) |
|---|
| 1122 | !c |
|---|
| 1123 | !c dclono2 : clono2 diffusivity, klassen et al., cm2 cm-1 |
|---|
| 1124 | !c |
|---|
| 1125 | dclono2 = 5.e-8*t/eta |
|---|
| 1126 | !c |
|---|
| 1127 | !c cbar |
|---|
| 1128 | !c |
|---|
| 1129 | cclono2 = 1474.*sqrt(t) |
|---|
| 1130 | !c |
|---|
| 1131 | !c kh2o : shi et al., s-1 |
|---|
| 1132 | !c |
|---|
| 1133 | kh2o = 1.95e10*exp(-2800./t) |
|---|
| 1134 | !c |
|---|
| 1135 | !c kh : shi et al., s-1 |
|---|
| 1136 | !c |
|---|
| 1137 | kh = 1.22e12*exp(-6200./t) |
|---|
| 1138 | !c |
|---|
| 1139 | !c khydr : shi et al., s-1 |
|---|
| 1140 | !c |
|---|
| 1141 | khydr = kh2o*aw(ilon,iniv) + kh*ah*aw(ilon,iniv) |
|---|
| 1142 | |
|---|
| 1143 | gbh2o = 4*hclono2(ilon,iniv)*0.082*t*sqrt(dclono2*khydr) & |
|---|
| 1144 | / cclono2 |
|---|
| 1145 | |
|---|
| 1146 | mhcl = hhcl(ilon,iniv)*phcl |
|---|
| 1147 | !c |
|---|
| 1148 | !c khcl : shi et al., s-1 |
|---|
| 1149 | !c |
|---|
| 1150 | khcl = 7.9e11*ah*dclono2*mhcl |
|---|
| 1151 | !c |
|---|
| 1152 | !c lclono2 : reacto-diffusive length, cm |
|---|
| 1153 | !c |
|---|
| 1154 | lclono2 = sqrt(dclono2/(khydr + khcl)) |
|---|
| 1155 | !c |
|---|
| 1156 | !c fclono2 |
|---|
| 1157 | !c |
|---|
| 1158 | fclono2 = 1./tanh(rmean(ilon,iniv)/lclono2) & |
|---|
| 1159 | - lclono2/rmean(ilon,iniv) |
|---|
| 1160 | |
|---|
| 1161 | grxn = fclono2*gbh2o*sqrt(1. + khcl/khydr) |
|---|
| 1162 | gbhcl = grxn*khcl/(khcl + khydr) |
|---|
| 1163 | gs = 66.12*hclono2(ilon,iniv)*mhcl*exp(-1374./t) |
|---|
| 1164 | |
|---|
| 1165 | fhcl = 1. & |
|---|
| 1166 | /(1 + 0.612*(gs + gbhcl)*pclono2/phcl) |
|---|
| 1167 | |
|---|
| 1168 | gps = fhcl*gs |
|---|
| 1169 | gbhclp = fhcl*gbhcl |
|---|
| 1170 | gb = gbhclp + grxn*khydr/(khcl + khydr) |
|---|
| 1171 | |
|---|
| 1172 | gclono2 = 1./(1. + 1./(gps + gb)) |
|---|
| 1173 | |
|---|
| 1174 | gclono2hcl = gclono2*(gps + gbhclp)/(gps + gb) |
|---|
| 1175 | gclono2h2o = gclono2 - gclono2hcl |
|---|
| 1176 | |
|---|
| 1177 | rclono2h2o = 0.25*gclono2h2o*cclono2*aliq(ilon,iniv)/zh2o |
|---|
| 1178 | k1l(ilon,iniv) = rclono2h2o |
|---|
| 1179 | |
|---|
| 1180 | rclono2hcl = 0.25*gclono2hcl*cclono2*aliq(ilon,iniv)/zhcl |
|---|
| 1181 | k2l(ilon,iniv) = rclono2hcl |
|---|
| 1182 | !c |
|---|
| 1183 | !C ===================================================================== |
|---|
| 1184 | !C ========================== HOCL + HCL =============================== |
|---|
| 1185 | !C ======================= ON LIQUID AEROSOL =========================== |
|---|
| 1186 | !C ===================================================================== |
|---|
| 1187 | !c |
|---|
| 1188 | !c d'apres jpl 2000 |
|---|
| 1189 | !c |
|---|
| 1190 | !c cbar |
|---|
| 1191 | !c |
|---|
| 1192 | chocl = 2009.*SQRT(t) |
|---|
| 1193 | !c |
|---|
| 1194 | !c dhocl : hocl diffusivity, klassen et al., cm2 cm-1 |
|---|
| 1195 | !c |
|---|
| 1196 | dhocl = 6.4e-8*t/eta |
|---|
| 1197 | !c |
|---|
| 1198 | !c khoclhcl : shi et al., s-1 |
|---|
| 1199 | !c |
|---|
| 1200 | khoclhcl = 1.25e9*ah*dhocl*mhcl |
|---|
| 1201 | !c |
|---|
| 1202 | !c lhocl : reacto-diffusive length, cm |
|---|
| 1203 | !c |
|---|
| 1204 | lhocl = SQRT(dhocl/khoclhcl) |
|---|
| 1205 | !c |
|---|
| 1206 | !c fhocl |
|---|
| 1207 | !c |
|---|
| 1208 | fhocl = 1./TANH(rmean(ilon,iniv)/lhocl) & |
|---|
| 1209 | - lhocl/rmean(ilon,iniv) |
|---|
| 1210 | |
|---|
| 1211 | grxn = 4.*hhocl(ilon,iniv)*0.082*t*SQRT(dhocl*khoclhcl)/chocl |
|---|
| 1212 | !c |
|---|
| 1213 | !c modif line ghoclhcl commenter |
|---|
| 1214 | !c ghoclhcl = 1./(1. + 1./(fhocl*grxn*fhcl)) |
|---|
| 1215 | !cc a la place |
|---|
| 1216 | |
|---|
| 1217 | IF ( (fhocl.LE.zero) .OR. (fhcl.LE.zero) .OR. (grxn.LE.zero) ) THEN |
|---|
| 1218 | ghoclhcl = 0. |
|---|
| 1219 | ELSE |
|---|
| 1220 | ghoclhcl = 1./(1. + 1./(fhocl*grxn*fhcl)) |
|---|
| 1221 | ENDIF |
|---|
| 1222 | |
|---|
| 1223 | rhoclhcl = 0.25*ghoclhcl*chocl*aliq(ilon,iniv)/zhcl |
|---|
| 1224 | k5l(ilon,iniv) = rhoclhcl |
|---|
| 1225 | !C |
|---|
| 1226 | !C ==================================================================== |
|---|
| 1227 | !C =========================== HOBR + HCL ============================= |
|---|
| 1228 | !C ======================= ON LIQUID AEROSOL ========================== |
|---|
| 1229 | !C ==================================================================== |
|---|
| 1230 | !c |
|---|
| 1231 | !c d'apres waschewsky and abbatt, j. phys. chem. a, |
|---|
| 1232 | !c 5312-5320, 1999. |
|---|
| 1233 | !c |
|---|
| 1234 | !c mhobr : concentration de hobr en phase liquide |
|---|
| 1235 | !c |
|---|
| 1236 | mhobr = hhobr(ilon,iniv)*phobr |
|---|
| 1237 | |
|---|
| 1238 | k2hobrhcl = EXP(0.542*wt(ilon,iniv) - 6.44e3/t + 10.3) |
|---|
| 1239 | |
|---|
| 1240 | ratio = mhobr/mhcl |
|---|
| 1241 | |
|---|
| 1242 | IF (ratio .LT. 1.) THEN |
|---|
| 1243 | !c |
|---|
| 1244 | !c hcl est en exces dans l'aerosol. |
|---|
| 1245 | !c hobr est le reactif limitant. |
|---|
| 1246 | !c |
|---|
| 1247 | k1hobrhcl = k2hobrhcl*mhcl |
|---|
| 1248 | !c |
|---|
| 1249 | !c dhobr : hobr diffusivity, klassen et al., cm2 cm-1 |
|---|
| 1250 | !c |
|---|
| 1251 | dhobr = 6.2e-8*t/eta |
|---|
| 1252 | !c |
|---|
| 1253 | !c lhobr : reacto-diffusive length, cm |
|---|
| 1254 | !c |
|---|
| 1255 | lhobr = SQRT(dhobr/k1hobrhcl) |
|---|
| 1256 | |
|---|
| 1257 | fhobr = 1./TANH(rmean(ilon,iniv)/lhobr) & |
|---|
| 1258 | - lhobr/rmean(ilon,iniv) |
|---|
| 1259 | |
|---|
| 1260 | chobr = 1477.*SQRT(t) |
|---|
| 1261 | |
|---|
| 1262 | ghobrhcl = 4.*hhobr(ilon,iniv)*0.082*t & |
|---|
| 1263 | *SQRT(dhobr*k1hobrhcl)*fhobr/chobr |
|---|
| 1264 | |
|---|
| 1265 | rhobrhcl = 0.25*ghobrhcl*chobr*aliq(ilon,iniv)/zhcl |
|---|
| 1266 | ELSE |
|---|
| 1267 | !c |
|---|
| 1268 | !c hobr est en exces dans l'aerosol. |
|---|
| 1269 | !c hcl est le reactif limitant. |
|---|
| 1270 | !c |
|---|
| 1271 | k1hobrhcl = k2hobrhcl*mhobr |
|---|
| 1272 | !c |
|---|
| 1273 | !c dhcl : hcl diffusivity, klassen et al., cm2 cm-1 |
|---|
| 1274 | !c |
|---|
| 1275 | dhcl = 7.8e-8*t/eta |
|---|
| 1276 | !c |
|---|
| 1277 | !c lhcl : reacto-diffusive length, cm |
|---|
| 1278 | !c |
|---|
| 1279 | lhcl = SQRT(dhcl/k1hobrhcl) |
|---|
| 1280 | |
|---|
| 1281 | fhcl = 1./TANH(rmean(ilon,iniv)/lhcl) & |
|---|
| 1282 | - lhcl/rmean(ilon,iniv) |
|---|
| 1283 | |
|---|
| 1284 | chcl = 2408.*SQRT(t) |
|---|
| 1285 | |
|---|
| 1286 | ghobrhcl = 4.*hhcl(ilon,iniv)*0.082*t & |
|---|
| 1287 | *SQRT(dhcl*k1hobrhcl)*fhcl/chcl |
|---|
| 1288 | |
|---|
| 1289 | rhobrhcl = 0.25*ghobrhcl*chcl*aliq(ilon,iniv)/zhobr |
|---|
| 1290 | ENDIF |
|---|
| 1291 | |
|---|
| 1292 | k7l(ilon,iniv) = rhobrhcl |
|---|
| 1293 | !C |
|---|
| 1294 | !C ==================================================================== |
|---|
| 1295 | !C =========================== HOBR + HBR ============================= |
|---|
| 1296 | !C ======================= ON LIQUID AEROSOL ========================== |
|---|
| 1297 | !C ==================================================================== |
|---|
| 1298 | !c |
|---|
| 1299 | k8l(ilon,iniv) = 0. |
|---|
| 1300 | !C |
|---|
| 1301 | !C ==================================================================== |
|---|
| 1302 | !C =========================== HOCL + HBR ============================= |
|---|
| 1303 | !C ======================= ON LIQUID AEROSOL ========================== |
|---|
| 1304 | !C ==================================================================== |
|---|
| 1305 | !c |
|---|
| 1306 | !c mhbr : concentration de hbr en phase liquide |
|---|
| 1307 | !c mhocl : concentration de hocl en phase liquide |
|---|
| 1308 | !c |
|---|
| 1309 | !c mhbr = hhbr(ilon,iniv)*phbr |
|---|
| 1310 | !c mhocl = hhocl(ilon,iniv)*phocl |
|---|
| 1311 | !c |
|---|
| 1312 | !c k2hoclhbr: abbatt and nowak, j. phys. chem. a, 2131, 1997. |
|---|
| 1313 | !c valeur a 228k et 69.3% h2so4 |
|---|
| 1314 | !c |
|---|
| 1315 | !c k2hoclhbr = 2.e6 |
|---|
| 1316 | !c |
|---|
| 1317 | !c ratio = mhbr/mhocl |
|---|
| 1318 | !c |
|---|
| 1319 | !c if (ratio .lt. 1.) then |
|---|
| 1320 | !c |
|---|
| 1321 | !c hocl est en exces dans l'aerosol. |
|---|
| 1322 | !c hbr est le reactif limitant. |
|---|
| 1323 | !c |
|---|
| 1324 | !c k1hoclhbr = k2hoclhbr*mhocl |
|---|
| 1325 | !c |
|---|
| 1326 | !c dhbr : hbr diffusivity, klassen et al., cm2 cm-1 |
|---|
| 1327 | !c |
|---|
| 1328 | !c dhbr = 7.9e-8*t/eta |
|---|
| 1329 | !c chbr = 1616.*sqrt(t) |
|---|
| 1330 | !c |
|---|
| 1331 | !c lhbr : reacto-diffusive length, cm |
|---|
| 1332 | !c |
|---|
| 1333 | !c lhbr = sqrt(dhbr/k1hoclhbr) |
|---|
| 1334 | !c |
|---|
| 1335 | !c fhbr = 1./tanh(rmean(ilon,iniv)/lhbr) |
|---|
| 1336 | !c $ - lhbr/rmean(ilon,iniv) |
|---|
| 1337 | !c |
|---|
| 1338 | !c ghoclhbr = 4.*hhbr(ilon,iniv)*0.082*t |
|---|
| 1339 | !c $ *sqrt(dhbr*k1hoclhbr)*fhbr/chbr |
|---|
| 1340 | !c |
|---|
| 1341 | !c else |
|---|
| 1342 | !c |
|---|
| 1343 | !c hbr est en exces dans l'aerosol. |
|---|
| 1344 | !c hocl est le reactif limitant. |
|---|
| 1345 | !c |
|---|
| 1346 | !c k1hoclhbr = k2hoclhbr*mhbr |
|---|
| 1347 | !c |
|---|
| 1348 | !c lhocl : reacto-diffusive length, cm |
|---|
| 1349 | !c |
|---|
| 1350 | !c lhocl = sqrt(dhocl/k1hoclhbr) |
|---|
| 1351 | !c |
|---|
| 1352 | !c fhocl = 1./tanh(rmean(ilon,iniv)/lhocl) |
|---|
| 1353 | !c $ - lhocl/rmean(ilon,iniv) |
|---|
| 1354 | !c |
|---|
| 1355 | !c ghoclhbr = 4.*hhocl(ilon,iniv)*0.082*t |
|---|
| 1356 | !c $ *sqrt(dhocl*k1hoclhbr)*fhocl/chocl |
|---|
| 1357 | !c endif |
|---|
| 1358 | !c |
|---|
| 1359 | k9l(ilon,iniv) = 0. |
|---|
| 1360 | !C |
|---|
| 1361 | !C ===================================================================== |
|---|
| 1362 | !C ======================== N2O5 + H2O ================================= |
|---|
| 1363 | !C ====================== ON LIQUID AEROSOL ============================ |
|---|
| 1364 | !C ===================================================================== |
|---|
| 1365 | !C |
|---|
| 1366 | GN2O5H2O=0.1 |
|---|
| 1367 | CBAR=1400.1*SQRT(T) |
|---|
| 1368 | |
|---|
| 1369 | RN2O5H2O=0.25*GN2O5H2O*CBAR*ALIQ(ILON,INIV)/ZH2O |
|---|
| 1370 | K3L(ILON,INIV) = RN2O5H2O |
|---|
| 1371 | !c |
|---|
| 1372 | !C ===================================================================== |
|---|
| 1373 | !C ======================== N2O5 + HCl ================================= |
|---|
| 1374 | !C ====================== ON LIQUID AEROSOL ============================ |
|---|
| 1375 | !C ===================================================================== |
|---|
| 1376 | !C |
|---|
| 1377 | K4L(ILON,INIV) = 0. |
|---|
| 1378 | !C |
|---|
| 1379 | !C ===================================================================== |
|---|
| 1380 | !C ======================== BrONO2 + H2O =============================== |
|---|
| 1381 | !C ====================== ON LIQUID AEROSOL ============================ |
|---|
| 1382 | !C ===================================================================== |
|---|
| 1383 | !C |
|---|
| 1384 | !c parametrisation jpl 2000 (hanson, private comm.) |
|---|
| 1385 | !c |
|---|
| 1386 | gxrn = exp(29.24 - 0.396*100.*ws(ilon,iniv)) + 0.114 |
|---|
| 1387 | gbrno3h2o = 1./(1./0.805 + 1./gxrn) |
|---|
| 1388 | |
|---|
| 1389 | CBAR=1221.4*SQRT(T) |
|---|
| 1390 | |
|---|
| 1391 | RBRNO3H2O=0.25*GBRNO3H2O*CBAR*ALIQ(ILON,INIV)/ZH2O |
|---|
| 1392 | K6L(ILON,INIV) = RBRNO3H2O |
|---|
| 1393 | |
|---|
| 1394 | ELSE !condliq |
|---|
| 1395 | |
|---|
| 1396 | K1L(ILON,INIV) = 0. |
|---|
| 1397 | K2L(ILON,INIV) = 0. |
|---|
| 1398 | K3L(ILON,INIV) = 0. |
|---|
| 1399 | K4L(ILON,INIV) = 0. |
|---|
| 1400 | K5L(ILON,INIV) = 0. |
|---|
| 1401 | K6L(ILON,INIV) = 0. |
|---|
| 1402 | K7L(ILON,INIV) = 0. |
|---|
| 1403 | K8L(ILON,INIV) = 0. |
|---|
| 1404 | K9L(ILON,INIV) = 0. |
|---|
| 1405 | |
|---|
| 1406 | ENDIF !condliq |
|---|
| 1407 | |
|---|
| 1408 | END DO |
|---|
| 1409 | END DO |
|---|
| 1410 | |
|---|
| 1411 | DO ILON = 1, PLON |
|---|
| 1412 | DO iniv = nivbot, nivtop |
|---|
| 1413 | |
|---|
| 1414 | IF (LNAT(ILON,INIV) .EQ. 1) THEN |
|---|
| 1415 | |
|---|
| 1416 | T = TEMP(ILON,INIV) |
|---|
| 1417 | !C |
|---|
| 1418 | !C Z VARIABLES ARE SET TO AVOID DIVISION BY ZERO |
|---|
| 1419 | !C |
|---|
| 1420 | ZH2O = MAX(h2oconc(ILON,INIV),1.E-08*HNM(ILON,INIV)) |
|---|
| 1421 | ZHCL = MAX(CK(ILON,INIV,id_HCl),1.E-13*HNM(ILON,INIV)) |
|---|
| 1422 | ZHOCL = MAX(CK(ILON,INIV,id_HOCl),1.E-13*HNM(ILON,INIV)) |
|---|
| 1423 | ZHOBR = MAX(CK(ILON,INIV,id_HOBr),1.E-13*HNM(ILON,INIV)) |
|---|
| 1424 | !C |
|---|
| 1425 | !C ===================================================================== |
|---|
| 1426 | !C ====================== CLNO3 + HCL ================================== |
|---|
| 1427 | !C ====================== CLNO3 + H2O ================================== |
|---|
| 1428 | !C ====================== N2O5 + H2O ================================== |
|---|
| 1429 | !C ====================== N2O5 + HCL ================================== |
|---|
| 1430 | !C ====================== HOCL + HCL ================================== |
|---|
| 1431 | !C ====================== BRNO3 + H2O ================================= |
|---|
| 1432 | !C ====================== HOBR + HCL ================================= |
|---|
| 1433 | !C ====================== HOBR + HBR ================================= |
|---|
| 1434 | !C ====================== HOCL + HBR ================================= |
|---|
| 1435 | !C ======================= ON NAT ====================================== |
|---|
| 1436 | !C ===================================================================== |
|---|
| 1437 | !C |
|---|
| 1438 | !C GAMMA VALUES |
|---|
| 1439 | !C |
|---|
| 1440 | |
|---|
| 1441 | !DH 2018 |
|---|
| 1442 | !GCLNO3HCLNAT=0.1 |
|---|
| 1443 | GCLNO3HCLNAT=0.2 |
|---|
| 1444 | |
|---|
| 1445 | GCLNO3H2ONAT=0.004 |
|---|
| 1446 | GN2O5H2ONAT=4.0E-4 |
|---|
| 1447 | GN2O5HCLNAT=3.0E-3 |
|---|
| 1448 | GHOCLHCLNAT=0.1 |
|---|
| 1449 | GBRNO3H2ONAT=0. |
|---|
| 1450 | GHOBRHCLNAT=0. |
|---|
| 1451 | GHOBRHBRNAT=0. |
|---|
| 1452 | GHOCLHBRNAT=0. |
|---|
| 1453 | !C |
|---|
| 1454 | !C BECAUSE NAT PARTICLES CAN BE VERY LARGE, GAS DIFFUSION LIMITATION |
|---|
| 1455 | !C IS TAKEN INTO ACCOUNT. THE SECOND ORDER RATE CONSTANTS ARE |
|---|
| 1456 | !C CALCULATED FROM THE EQUATION R=G*PI*R**2*CBAR*N/(1+3G.R/(4.L)), WHERE |
|---|
| 1457 | !C G=GAMMA, R=PARTICLE RADIUS, CBAR=MEAN MOLECULAR SPEED, L=MEAN FREE |
|---|
| 1458 | !C PATH (TURCO ET AL, JGR, 94, 16493, 1989). |
|---|
| 1459 | !C |
|---|
| 1460 | RCLNO3H2ONAT = & |
|---|
| 1461 | GCLNO3H2ONAT*4.56E4*SQRT(T/97.45)*RNAT(ILON,INIV)**2 & |
|---|
| 1462 | *NNAT(ILON,INIV)/ & |
|---|
| 1463 | (1.0+3.3E4*GCLNO3H2ONAT*RNAT(ILON,INIV)*PTOT(ILON,INIV)/T) & |
|---|
| 1464 | /ZH2O |
|---|
| 1465 | K1S(ILON,INIV) = RCLNO3H2ONAT |
|---|
| 1466 | |
|---|
| 1467 | RCLNO3HCLNAT= & |
|---|
| 1468 | GCLNO3HCLNAT*4.56E4*SQRT(T/97.45)*RNAT(ILON,INIV)**2 & |
|---|
| 1469 | *NNAT(ILON,INIV)/ & |
|---|
| 1470 | (1.0+3.3E4*GCLNO3HCLNAT*RNAT(ILON,INIV)*PTOT(ILON,INIV)/T) & |
|---|
| 1471 | /ZHCL |
|---|
| 1472 | K2S(ILON,INIV) = RCLNO3HCLNAT |
|---|
| 1473 | |
|---|
| 1474 | RN2O5H2ONAT = & |
|---|
| 1475 | GN2O5H2ONAT*4.56E4*SQRT(T/108.0)*RNAT(ILON,INIV)**2 & |
|---|
| 1476 | *NNAT(ILON,INIV)/ & |
|---|
| 1477 | (1.0+3.3E4*GN2O5H2ONAT*RNAT(ILON,INIV)*PTOT(ILON,INIV)/T) & |
|---|
| 1478 | /ZH2O |
|---|
| 1479 | K3S(ILON,INIV) = RN2O5H2ONAT |
|---|
| 1480 | |
|---|
| 1481 | RBRNO3H2ONAT = & |
|---|
| 1482 | GBRNO3H2ONAT*4.56E4*SQRT(T/142.0)*RNAT(ILON,INIV)**2 & |
|---|
| 1483 | *NNAT(ILON,INIV)/ & |
|---|
| 1484 | (1.0+3.3E4*GBRNO3H2ONAT*RNAT(ILON,INIV)*PTOT(ILON,INIV)/T) & |
|---|
| 1485 | /ZH2O |
|---|
| 1486 | K6S(ILON,INIV) = RBRNO3H2ONAT |
|---|
| 1487 | |
|---|
| 1488 | RN2O5HCLNAT= & |
|---|
| 1489 | GN2O5HCLNAT*4.56E4*SQRT(T/108.0)*RNAT(ILON,INIV)**2 & |
|---|
| 1490 | *NNAT(ILON,INIV)/ & |
|---|
| 1491 | (1.0+3.3E4*GN2O5HCLNAT*RNAT(ILON,INIV)*PTOT(ILON,INIV)/T) & |
|---|
| 1492 | /ZHCL |
|---|
| 1493 | K4S(ILON,INIV) = RN2O5HCLNAT |
|---|
| 1494 | |
|---|
| 1495 | RHOCLHCLNAT= & |
|---|
| 1496 | GHOCLHCLNAT*4.56E4*SQRT(T/52.5)*RNAT(ILON,INIV)**2 & |
|---|
| 1497 | *NNAT(ILON,INIV)/ & |
|---|
| 1498 | (1.0+3.3E4*GHOCLHCLNAT*RNAT(ILON,INIV)*PTOT(ILON,INIV)/T) & |
|---|
| 1499 | /ZHCL |
|---|
| 1500 | K5S(ILON,INIV) = RHOCLHCLNAT |
|---|
| 1501 | |
|---|
| 1502 | RHOBRHCLNAT= & |
|---|
| 1503 | GHOBRHCLNAT*4.56E4*SQRT(T/96.9)*RNAT(ILON,INIV)**2 & |
|---|
| 1504 | *NNAT(ILON,INIV)/ & |
|---|
| 1505 | (1.0+3.3E4*GHOBRHCLNAT*RNAT(ILON,INIV)*PTOT(ILON,INIV)/T) & |
|---|
| 1506 | /ZHCL |
|---|
| 1507 | K7S(ILON,INIV) = RHOBRHCLNAT |
|---|
| 1508 | |
|---|
| 1509 | K8S(ILON,INIV) = 0. |
|---|
| 1510 | K9S(ILON,INIV) = 0. |
|---|
| 1511 | |
|---|
| 1512 | ENDIF |
|---|
| 1513 | |
|---|
| 1514 | END DO |
|---|
| 1515 | END DO |
|---|
| 1516 | |
|---|
| 1517 | DO ILON = 1, PLON |
|---|
| 1518 | DO iniv = nivbot, nivtop |
|---|
| 1519 | |
|---|
| 1520 | IF (LICE(ILON,INIV) .EQ. 1) THEN |
|---|
| 1521 | |
|---|
| 1522 | T = TEMP(ILON,INIV) |
|---|
| 1523 | !C |
|---|
| 1524 | !C Z VARIABLES ARE SET TO AVOID DIVISION BY ZERO |
|---|
| 1525 | !C |
|---|
| 1526 | ZH2O = MAX(h2oconc(ILON,INIV),1.E-08*HNM(ILON,INIV)) |
|---|
| 1527 | ZHCL = MAX(CK(ILON,INIV,id_HCl),1.E-13*HNM(ILON,INIV)) |
|---|
| 1528 | ZHOCL = MAX(CK(ILON,INIV,id_HOCl),1.E-13*HNM(ILON,INIV)) |
|---|
| 1529 | ZHOBR = MAX(CK(ILON,INIV,id_HOBr),1.E-13*HNM(ILON,INIV)) |
|---|
| 1530 | !C |
|---|
| 1531 | !C ===================================================================== |
|---|
| 1532 | !C ====================== CLNO3 + HCL ================================== |
|---|
| 1533 | !C ====================== CLNO3 + H2O ================================== |
|---|
| 1534 | !C ====================== N2O5 + H2O ================================== |
|---|
| 1535 | !C ====================== N2O5 + HCL ================================== |
|---|
| 1536 | !C ====================== HOCL + HCL ================================== |
|---|
| 1537 | !C ====================== BRNO3 + H2O ================================= |
|---|
| 1538 | !C ====================== HOBR + HCL ================================= |
|---|
| 1539 | !C ====================== HOBR + HBR ================================= |
|---|
| 1540 | !C ====================== HOCL + HBR ================================= |
|---|
| 1541 | !C ======================= ON ICE ====================================== |
|---|
| 1542 | !C ===================================================================== |
|---|
| 1543 | !C GAMMA VALUES |
|---|
| 1544 | !C |
|---|
| 1545 | GCLNO3HCLICE=0.3 |
|---|
| 1546 | GCLNO3H2OICE=0.3 |
|---|
| 1547 | GN2O5H2OICE=0.02 |
|---|
| 1548 | GN2O5HCLICE=0.03 |
|---|
| 1549 | GHOCLHCLICE=0.2 |
|---|
| 1550 | GBRNO3H2OICE=0.3 |
|---|
| 1551 | GHOBRHCLICE=0.3 |
|---|
| 1552 | GHOBRHBRICE=0.1 |
|---|
| 1553 | GHOCLHBRICE=0.1 !FJegou 13/10/06 JPL 2006 |
|---|
| 1554 | !C |
|---|
| 1555 | !C BECAUSE ICE PARTICLES CAN BE VERY LARGE, GAS DIFFUSION LIMITATION |
|---|
| 1556 | !C IS TAKEN INTO ACCOUNT. THE SECOND ORDER RATE CONSTANTS ARE |
|---|
| 1557 | !C CALCULATED FROM THE EQUATION R=G*PI*R**2*CBAR*N/(1+3G.R/(4.L)), WHERE |
|---|
| 1558 | !C G=GAMMA, R=PARTICLE RADIUS, CBAR=MEAN MOLECULAR SPEED, L=MEAN FREE |
|---|
| 1559 | !C PATH (TURCO ET AL, JGR, 94, 16493, 1989). |
|---|
| 1560 | !C |
|---|
| 1561 | RCLNO3HCLICE= & |
|---|
| 1562 | GCLNO3HCLICE*4.56E4*SQRT(T/97.45)*RICE(ILON,INIV)**2 & |
|---|
| 1563 | *NICE(ILON,INIV)/ & |
|---|
| 1564 | (1.0+3.3E4*GCLNO3HCLICE*RICE(ILON,INIV)*PTOT(ILON,INIV)/T) & |
|---|
| 1565 | /ZHCL |
|---|
| 1566 | K2S(ILON,INIV) = RCLNO3HCLICE |
|---|
| 1567 | |
|---|
| 1568 | RCLNO3H2OICE= & |
|---|
| 1569 | GCLNO3H2OICE*4.56E4*SQRT(T/97.45)*RICE(ILON,INIV)**2 & |
|---|
| 1570 | *NICE(ILON,INIV)/ & |
|---|
| 1571 | (1.0+3.3E4*GCLNO3H2OICE*RICE(ILON,INIV)*PTOT(ILON,INIV)/T) & |
|---|
| 1572 | /ZH2O |
|---|
| 1573 | K1S(ILON,INIV) = RCLNO3H2OICE |
|---|
| 1574 | |
|---|
| 1575 | RN2O5H2OICE= & |
|---|
| 1576 | GN2O5H2OICE*4.56E4*SQRT(T/108.0)*RICE(ILON,INIV)**2 & |
|---|
| 1577 | *NICE(ILON,INIV)/ & |
|---|
| 1578 | (1.0+3.3E4*GN2O5H2OICE*RICE(ILON,INIV)*PTOT(ILON,INIV)/T) & |
|---|
| 1579 | /ZH2O |
|---|
| 1580 | K3S(ILON,INIV) = RN2O5H2OICE |
|---|
| 1581 | |
|---|
| 1582 | RN2O5HCLICE= & |
|---|
| 1583 | GN2O5HCLICE*4.56E4*SQRT(T/108.0)*RICE(ILON,INIV)**2 & |
|---|
| 1584 | *NICE(ILON,INIV)/ & |
|---|
| 1585 | (1.0+3.3E4*GN2O5HCLICE*RICE(ILON,INIV)*PTOT(ILON,INIV)/T) & |
|---|
| 1586 | /ZHCL |
|---|
| 1587 | K4S(ILON,INIV) = RN2O5HCLICE |
|---|
| 1588 | |
|---|
| 1589 | RHOCLHCLICE= & |
|---|
| 1590 | GHOCLHCLICE*4.56E4*SQRT(T/52.5)*RICE(ILON,INIV)**2 & |
|---|
| 1591 | *NICE(ILON,INIV)/ & |
|---|
| 1592 | (1.0+3.3E4*GHOCLHCLICE*RICE(ILON,INIV)*PTOT(ILON,INIV)/T) & |
|---|
| 1593 | /ZHCL |
|---|
| 1594 | K5S(ILON,INIV) = RHOCLHCLICE |
|---|
| 1595 | |
|---|
| 1596 | RBRNO3H2OICE= & |
|---|
| 1597 | GBRNO3H2OICE*4.56E4*SQRT(T/142.0)*RICE(ILON,INIV)**2 & |
|---|
| 1598 | *NICE(ILON,INIV)/ & |
|---|
| 1599 | (1.0+3.3E4*GBRNO3H2OICE*RICE(ILON,INIV)*PTOT(ILON,INIV)/T) & |
|---|
| 1600 | /ZH2O |
|---|
| 1601 | K6S(ILON,INIV) = RBRNO3H2OICE |
|---|
| 1602 | |
|---|
| 1603 | RHOBRHCLICE= & |
|---|
| 1604 | GHOBRHCLICE*4.56E4*SQRT(T/96.9)*RICE(ILON,INIV)**2 & |
|---|
| 1605 | *NICE(ILON,INIV)/ & |
|---|
| 1606 | (1.0+3.3E4*GHOBRHCLICE*RICE(ILON,INIV)*PTOT(ILON,INIV)/T) & |
|---|
| 1607 | /ZHCL |
|---|
| 1608 | K7S(ILON,INIV) = RHOBRHCLICE |
|---|
| 1609 | |
|---|
| 1610 | RHOBRHBRICE= & |
|---|
| 1611 | GHOBRHBRICE*4.56E4*SQRT(T/96.9)*RICE(ILON,INIV)**2 & |
|---|
| 1612 | *NICE(ILON,INIV)/ & |
|---|
| 1613 | (1.0+3.3E4*GHOBRHBRICE*RICE(ILON,INIV)*PTOT(ILON,INIV)/T) & |
|---|
| 1614 | /ZHOBR |
|---|
| 1615 | K8S(ILON,INIV) = RHOBRHBRICE |
|---|
| 1616 | |
|---|
| 1617 | END IF |
|---|
| 1618 | |
|---|
| 1619 | END DO |
|---|
| 1620 | END DO |
|---|
| 1621 | |
|---|
| 1622 | DO ILON = 1, PLON |
|---|
| 1623 | DO iniv = nivbot, nivtop |
|---|
| 1624 | K1(ILON,INIV) = K1L(ILON,INIV)+ (LNAT(ILON,INIV) + LICE(ILON,INIV))*K1S(ILON,INIV) |
|---|
| 1625 | K2(ILON,INIV) = K2L(ILON,INIV)+ (LNAT(ILON,INIV) + LICE(ILON,INIV))*K2S(ILON,INIV) |
|---|
| 1626 | K3(ILON,INIV) = K3L(ILON,INIV)+ (LNAT(ILON,INIV) + LICE(ILON,INIV))*K3S(ILON,INIV) |
|---|
| 1627 | K4(ILON,INIV) = K4L(ILON,INIV)+ (LNAT(ILON,INIV) + LICE(ILON,INIV))*K4S(ILON,INIV) |
|---|
| 1628 | K5(ILON,INIV) = K5L(ILON,INIV)+ (LNAT(ILON,INIV) + LICE(ILON,INIV))*K5S(ILON,INIV) |
|---|
| 1629 | K6(ILON,INIV) = K6L(ILON,INIV)+ (LNAT(ILON,INIV) + LICE(ILON,INIV))*K6S(ILON,INIV) |
|---|
| 1630 | K7(ILON,INIV) = K7L(ILON,INIV)+ (LNAT(ILON,INIV) + LICE(ILON,INIV))*K7S(ILON,INIV) |
|---|
| 1631 | K8(ILON,INIV) = K8L(ILON,INIV)+ (LNAT(ILON,INIV) + LICE(ILON,INIV))*K8S(ILON,INIV) |
|---|
| 1632 | K9(ILON,INIV) = 0. |
|---|
| 1633 | END DO |
|---|
| 1634 | END DO |
|---|
| 1635 | |
|---|
| 1636 | DO ILON = 1, PLON |
|---|
| 1637 | DO iniv = nivbot, nivtop |
|---|
| 1638 | K1(ILON,INIV) = MAX(K1(ILON,INIV),0.) |
|---|
| 1639 | K2(ILON,INIV) = MAX(K2(ILON,INIV),0.) |
|---|
| 1640 | K3(ILON,INIV) = MAX(K3(ILON,INIV),0.) |
|---|
| 1641 | K4(ILON,INIV) = MAX(K4(ILON,INIV),0.) |
|---|
| 1642 | K5(ILON,INIV) = MAX(K5(ILON,INIV),0.) |
|---|
| 1643 | K6(ILON,INIV) = MAX(K6(ILON,INIV),0.) |
|---|
| 1644 | K7(ILON,INIV) = MAX(K7(ILON,INIV),0.) |
|---|
| 1645 | K8(ILON,INIV) = MAX(K8(ILON,INIV),0.) |
|---|
| 1646 | K9(ILON,INIV) = MAX(K9(ILON,INIV),0.) |
|---|
| 1647 | END DO |
|---|
| 1648 | END DO |
|---|
| 1649 | |
|---|
| 1650 | RETURN |
|---|
| 1651 | END SUBROUTINE hetero |
|---|
| 1652 | |
|---|
| 1653 | SUBROUTINE density(WS,WN,TEMP,DENS) |
|---|
| 1654 | !C |
|---|
| 1655 | !C DENSITY OF TERNARY SOLUTION IN G/CM3 |
|---|
| 1656 | !C WS ,WN ARE WT FRACTION, |
|---|
| 1657 | !C FITTED TO 0.05<WS+WN<0.70 WT FRACTION, BUT EXTRAPOLATES WELL |
|---|
| 1658 | !C 185 < T (K) |
|---|
| 1659 | |
|---|
| 1660 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 1661 | !c F. LEFEVRE (SA/CNRS) from REPROBUS |
|---|
| 1662 | !c |
|---|
| 1663 | !c Adapted to LMDz-INCA: |
|---|
| 1664 | !c |
|---|
| 1665 | !c L. JOURDAIN (SA) 2004 |
|---|
| 1666 | !c D. HAUGLUSTAINE (LSCE) 06/2005 |
|---|
| 1667 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 1668 | |
|---|
| 1669 | INTEGER :: ilon, iniv |
|---|
| 1670 | |
|---|
| 1671 | REAL :: t, w, wh, v1, vs, vn, vmcal |
|---|
| 1672 | |
|---|
| 1673 | REAL WS(PLON,PLEV),WN(PLON,PLEV),TEMP(PLON,PLEV) |
|---|
| 1674 | REAL DENS(PLON,PLEV) |
|---|
| 1675 | REAL X2(22) |
|---|
| 1676 | SAVE X2 |
|---|
| 1677 | !$OMP THREADPRIVATE(X2) |
|---|
| 1678 | DATA X2/ & |
|---|
| 1679 | 2.393284E-02,-4.359335E-05,7.961181E-08,0.0,-0.198716351, & |
|---|
| 1680 | 1.39564574E-03,-2.020633E-06,0.51684706,-3.0539E-03, & |
|---|
| 1681 | 4.505475E-06,-0.30119511,1.840408E-03,-2.7221253742E-06, & |
|---|
| 1682 | -0.11331674116,8.47763E-04,-1.22336185E-06,0.3455282, & |
|---|
| 1683 | -2.2111E-03,3.503768245E-06,-0.2315332,1.60074E-03, & |
|---|
| 1684 | -2.5827835E-06/ |
|---|
| 1685 | |
|---|
| 1686 | DO ILON = 1, PLON |
|---|
| 1687 | DO iniv = nivbot, nivtop |
|---|
| 1688 | |
|---|
| 1689 | T = TEMP(ILON,INIV) |
|---|
| 1690 | W = WS(ILON,INIV) + WN(ILON,INIV) |
|---|
| 1691 | WH= 1.0 - W |
|---|
| 1692 | V1=X2(1)+X2(2)*T+X2(3)*T**2+X2(4)*T**3 |
|---|
| 1693 | VS=X2(5)+X2(6)*T+X2(7)*T**2+(X2(8)+X2(9)*T+X2(10)*T**2)*W & |
|---|
| 1694 | +(X2(11)+X2(12)*T+X2(13)*T**2)*W*W |
|---|
| 1695 | VN=X2(14)+X2(15)*T+X2(16)*T**2+(X2(17)+X2(18)*T+X2(19)*T**2)*W & |
|---|
| 1696 | +(X2(20)+X2(21)*T+X2(22)*T**2)*W*W |
|---|
| 1697 | VMCAL=WH/18.0160*V1 + VS*WS(ILON,INIV)/98.080 & |
|---|
| 1698 | + VN*WN(ILON,INIV)/63.0160 |
|---|
| 1699 | DENS(ILON,INIV) = 1.0/VMCAL*0.001 |
|---|
| 1700 | |
|---|
| 1701 | END DO |
|---|
| 1702 | END DO |
|---|
| 1703 | |
|---|
| 1704 | RETURN |
|---|
| 1705 | END SUBROUTINE density |
|---|
| 1706 | |
|---|
| 1707 | SUBROUTINE STRSEDCALC (hnm, pdel, qj1, dt, h2o, h2oc) |
|---|
| 1708 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 1709 | !c F. LEFEVRE (SA/CNRS) from REPROBUS |
|---|
| 1710 | !c |
|---|
| 1711 | !c Adapted to LMDz-INCA: |
|---|
| 1712 | !c |
|---|
| 1713 | !c L. JOURDAIN (SA) 2004 |
|---|
| 1714 | !c D. HAUGLUSTAINE (LSCE) 06/2005 |
|---|
| 1715 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 1716 | |
|---|
| 1717 | USE SPECIES_NAMES |
|---|
| 1718 | USE AIRPLANE_SRC, ONLY : itrop |
|---|
| 1719 | |
|---|
| 1720 | INTEGER :: ilon, iniv |
|---|
| 1721 | !c |
|---|
| 1722 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 1723 | !c calcule les quantites de h2o et hno3 sedimentees |
|---|
| 1724 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 1725 | !c |
|---|
| 1726 | real pdel(PLON,PLEV) |
|---|
| 1727 | real qj1(PLON,PLEV,PCNST) |
|---|
| 1728 | real dt |
|---|
| 1729 | real hnm(PLON,PLEV) |
|---|
| 1730 | real h2oc(PLON,PLEV) |
|---|
| 1731 | real h2o(PLON,PLEV) |
|---|
| 1732 | real dz, fracup, fracdown |
|---|
| 1733 | !c |
|---|
| 1734 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 1735 | !c boucle sur les latitudes |
|---|
| 1736 | !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc |
|---|
| 1737 | !ccc calcul de la concentration |
|---|
| 1738 | !c |
|---|
| 1739 | !ccc calcul des fractions sedimentees |
|---|
| 1740 | !c |
|---|
| 1741 | DO ilon = 1,PLON |
|---|
| 1742 | DO iniv =itrop(ilon),nivtop-1,1 |
|---|
| 1743 | |
|---|
| 1744 | dz =pdel(ilon,iniv)/(hnm(ilon,iniv)*1.e06*28.9*1.e-03) & |
|---|
| 1745 | /9.8*6.02e23 |
|---|
| 1746 | |
|---|
| 1747 | fracup = vsed(ilon,iniv+1)*dt/dz |
|---|
| 1748 | fracdown = vsed(ilon,iniv )*dt/dz |
|---|
| 1749 | |
|---|
| 1750 | fracup = MIN(fracup,1.0) |
|---|
| 1751 | fracdown = MIN(fracdown,1.0) |
|---|
| 1752 | !c |
|---|
| 1753 | !ccc sedimentation de h2o |
|---|
| 1754 | |
|---|
| 1755 | h2oc(ilon,iniv) = & |
|---|
| 1756 | h2oc(ilon,iniv) & |
|---|
| 1757 | + h2oc(ilon,iniv+1) & |
|---|
| 1758 | *hnm(ilon,iniv+1)/hnm(ilon,iniv)*fracup & |
|---|
| 1759 | - h2oc(ilon,iniv)*fracdown |
|---|
| 1760 | |
|---|
| 1761 | !ccc update h2o |
|---|
| 1762 | |
|---|
| 1763 | h2o(ilon,iniv) = qj1(ilon,iniv,id_H2O) + h2oc(ilon,iniv) |
|---|
| 1764 | |
|---|
| 1765 | !ccc sedimentation de hno3 |
|---|
| 1766 | |
|---|
| 1767 | qj1(ilon,iniv,id_HNO3) = & |
|---|
| 1768 | qj1(ilon,iniv,id_HNO3) & |
|---|
| 1769 | + qj1(ilon,iniv+1,id_HNO3) & |
|---|
| 1770 | *(1. - parthno3sl(ilon,iniv+1)) & |
|---|
| 1771 | *hnm(ilon,iniv+1)/hnm(ilon,iniv)*fracup & |
|---|
| 1772 | - qj1(ilon,iniv,id_HNO3) & |
|---|
| 1773 | *(1. - parthno3sl(ilon,iniv))*fracdown |
|---|
| 1774 | |
|---|
| 1775 | END DO |
|---|
| 1776 | END DO |
|---|
| 1777 | |
|---|
| 1778 | RETURN |
|---|
| 1779 | END SUBROUTINE STRSEDCALC |
|---|
| 1780 | |
|---|
| 1781 | END MODULE hetchem |
|---|
| 1782 | #endif |
|---|
