| 1 | !$Id: chimieaq.F90 163 2010-02-22 15:41:45Z acosce $ |
|---|
| 2 | !! ========================================================================= |
|---|
| 3 | !! INCA - INteraction with Chemistry and Aerosols |
|---|
| 4 | !! |
|---|
| 5 | !! Copyright Laboratoire des Sciences du Climat et de l'Environnement (LSCE) |
|---|
| 6 | !! Unite mixte CEA-CNRS-UVSQ |
|---|
| 7 | !! |
|---|
| 8 | !! Contributors to this INCA subroutine: |
|---|
| 9 | !! |
|---|
| 10 | !! Olivier Boucher |
|---|
| 11 | !! Christiane Textor |
|---|
| 12 | !! Michael Schulz, LSCE, Michael.Schulz@cea.fr |
|---|
| 13 | !! |
|---|
| 14 | !! Anne Cozic, LSCE, anne.cozic@cea.fr |
|---|
| 15 | !! Yann Meurdesoif, LSCE, yann.meurdesoif@cea.fr |
|---|
| 16 | !! |
|---|
| 17 | !! This software is a computer program whose purpose is to simulate the |
|---|
| 18 | !! atmospheric gas phase and aerosol composition. The model is designed to be |
|---|
| 19 | !! used within a transport model or a general circulation model. This version |
|---|
| 20 | !! of INCA was designed to be coupled to the LMDz GCM. LMDz-INCA accounts |
|---|
| 21 | !! for emissions, transport (resolved and sub-grid scale), photochemical |
|---|
| 22 | !! transformations, and scavenging (dry deposition and washout) of chemical |
|---|
| 23 | !! species and aerosols interactively in the GCM. Several versions of the INCA |
|---|
| 24 | !! model are currently used depending on the envisaged applications with the |
|---|
| 25 | !! chemistry-climate model. |
|---|
| 26 | !! |
|---|
| 27 | !! This software is governed by the CeCILL license under French law and |
|---|
| 28 | !! abiding by the rules of distribution of free software. You can use, |
|---|
| 29 | !! modify and/ or redistribute the software under the terms of the CeCILL |
|---|
| 30 | !! license as circulated by CEA, CNRS and INRIA at the following URL |
|---|
| 31 | !! "http://www.cecill.info". |
|---|
| 32 | !! |
|---|
| 33 | !! As a counterpart to the access to the source code and rights to copy, |
|---|
| 34 | !! modify and redistribute granted by the license, users are provided only |
|---|
| 35 | !! with a limited warranty and the software's author, the holder of the |
|---|
| 36 | !! economic rights, and the successive licensors have only limited |
|---|
| 37 | !! liability. |
|---|
| 38 | !! |
|---|
| 39 | !! In this respect, the user's attention is drawn to the risks associated |
|---|
| 40 | !! with loading, using, modifying and/or developing or reproducing the |
|---|
| 41 | !! software by the user in light of its specific status of free software, |
|---|
| 42 | !! that may mean that it is complicated to manipulate, and that also |
|---|
| 43 | !! therefore means that it is reserved for developers and experienced |
|---|
| 44 | !! professionals having in-depth computer knowledge. Users are therefore |
|---|
| 45 | !! encouraged to load and test the software's suitability as regards their |
|---|
| 46 | !! requirements in conditions enabling the security of their systems and/or |
|---|
| 47 | !! data to be ensured and, more generally, to use and operate it in the |
|---|
| 48 | !! same conditions as regards security. |
|---|
| 49 | !! |
|---|
| 50 | !! The fact that you are presently reading this means that you have had |
|---|
| 51 | !! knowledge of the CeCILL license and that you accept its terms. |
|---|
| 52 | !! ========================================================================= |
|---|
| 53 | |
|---|
| 54 | #include <inca_define.h> |
|---|
| 55 | |
|---|
| 56 | #ifndef DUSS |
|---|
| 57 | #ifdef AER |
|---|
| 58 | SUBROUTINE chimieaq(& |
|---|
| 59 | pdtphys & ! time step of physics |
|---|
| 60 | ,rneb & ! cloud volume fraction in grid box |
|---|
| 61 | ,qliq & ! kg kg-1 prescribed |
|---|
| 62 | ,pplay & ! =pmid pressure in grid box center |
|---|
| 63 | ,t & ! =t_seri temperautre |
|---|
| 64 | ,tr_seri & ! tracer |
|---|
| 65 | #ifdef AERONLY |
|---|
| 66 | ,invariants & |
|---|
| 67 | #endif |
|---|
| 68 | ) |
|---|
| 69 | |
|---|
| 70 | ! calculation of aqueous phase chemistry of sulfur |
|---|
| 71 | ! written by Olivier Boucher |
|---|
| 72 | ! units are [molecules/cm^3] |
|---|
| 73 | ! modified by Christiane Textor for use in INCA |
|---|
| 74 | ! |
|---|
| 75 | |
|---|
| 76 | USE SPECIES_NAMES |
|---|
| 77 | USE CHEM_MODS, ONLY : adv_mass, nadv_mass, ph_hist,asso4m_p_so2h2o2,asso4m_p_so2o3 |
|---|
| 78 | USE RATE_INDEX_MOD |
|---|
| 79 | USE INCA_DIM |
|---|
| 80 | USE PRINT_INCA |
|---|
| 81 | USE RATE_INDEX_MOD |
|---|
| 82 | |
|---|
| 83 | IMPLICIT NONE |
|---|
| 84 | |
|---|
| 85 | INCLUDE "YOMCST_I.h" |
|---|
| 86 | |
|---|
| 87 | !--Parameters of the scheme |
|---|
| 88 | REAL:: pdtaq !--timestep for aqueous chemistry [s] |
|---|
| 89 | PARAMETER (pdtaq=120.) |
|---|
| 90 | REAL:: t_melt, t0, t1 !--temperature range for aquous-phase chemistry |
|---|
| 91 | PARAMETER (t_melt=273.15) |
|---|
| 92 | PARAMETER (t0=t_melt-30.) |
|---|
| 93 | PARAMETER (t1=t_melt-10.) |
|---|
| 94 | REAL:: xx_clean, xx_polluted, xx !--NH4+/SO4= ratio in droplet |
|---|
| 95 | PARAMETER (xx_clean=1.0, xx_polluted=1.5) |
|---|
| 96 | INTEGER:: niter !--No of iteration for calculating pH |
|---|
| 97 | PARAMETER (niter=5) |
|---|
| 98 | |
|---|
| 99 | !--Input and output parameters |
|---|
| 100 | REAL, INTENT(in) :: pdtphys ! timestep in seconds of physics |
|---|
| 101 | REAL, INTENT(in) :: t(PLON,PLEV) ! temperature |
|---|
| 102 | REAL, INTENT(in) :: pplay(PLON,PLEV) ! pression |
|---|
| 103 | REAL, INTENT(inout) :: tr_seri(PLON,PLEV,PCNST) ! traceurs ?? |
|---|
| 104 | REAL, INTENT(in) :: qliq ! kg kg-1 prescribed |
|---|
| 105 | REAL, INTENT(in) :: rneb(PLON,PLEV) ! cloud volume fraction in grid box |
|---|
| 106 | #ifdef AERONLY |
|---|
| 107 | REAL, INTENT(in) :: invariants(PLON,PLEV,NFS) ! invariant array |
|---|
| 108 | #endif |
|---|
| 109 | |
|---|
| 110 | !--Local variables |
|---|
| 111 | REAL zh2o2(PLON,PLEV), zo3(PLON,PLEV) |
|---|
| 112 | REAL zso2(PLON,PLEV), zso4(PLON,PLEV) |
|---|
| 113 | REAL tend, prod, beta |
|---|
| 114 | REAL k1(PLON,PLEV), k2(PLON,PLEV), KE |
|---|
| 115 | PARAMETER (KE=1.e-14) |
|---|
| 116 | REAL qliq2(PLON,PLEV), prod1, prod2, tend1, tend2 |
|---|
| 117 | REAL henry_h2o2, kk_h2o2 |
|---|
| 118 | REAL henry_o3, kk_o3 |
|---|
| 119 | REAL henry_so2, kk_so2 |
|---|
| 120 | PARAMETER (henry_h2o2=8.33e4,kk_h2o2=7379.) |
|---|
| 121 | PARAMETER (henry_o3=1.15e-2, kk_o3=2560.) |
|---|
| 122 | PARAMETER (henry_so2=1.2, kk_so2=3200.) |
|---|
| 123 | REAL scavh2o2(PLON,PLEV), scavo3(PLON,PLEV) |
|---|
| 124 | REAL so2_aq(PLON), h2o2_aq, o3_aq |
|---|
| 125 | REAL so4_aq(PLON), hso3m_aq(PLON) |
|---|
| 126 | REAL so3mm_aq(PLON), no3m_aq(PLON) |
|---|
| 127 | REAL bb, cc, hplus(PLON), ph, correc |
|---|
| 128 | INTEGER i, k, pas, nbpas, iter |
|---|
| 129 | REAL zrho, henry_t, f_a(PLON,PLEV) |
|---|
| 130 | REAL mmr2moleqcm, moleqcm2mmr |
|---|
| 131 | REAL tfactor, f_a_o3, f_a_h2o2 |
|---|
| 132 | REAL rho_water_inv,tfac,fac,fac1 |
|---|
| 133 | PARAMETER (rho_water_inv=1.e-3) !--kg m-3 |
|---|
| 134 | |
|---|
| 135 | nbpas=pdtphys/pdtaq |
|---|
| 136 | IF ((FLOAT(nbpas)*pdtaq-pdtphys).GT.0.01) THEN |
|---|
| 137 | CALL print_err(3, "CHIMIEAQ", 'pdtphys doit etre un multiple de pdtaq', '', '') |
|---|
| 138 | ENDIF |
|---|
| 139 | |
|---|
| 140 | !--initialisations and conversion of units from mmr to molecules/cm^3 |
|---|
| 141 | |
|---|
| 142 | fac1=1.e-3*RNAVO |
|---|
| 143 | DO k = 1, PLEV |
|---|
| 144 | DO i = 1, PLON |
|---|
| 145 | ph_hist(i,k)=0.0 |
|---|
| 146 | asso4m_p_so2o3(i,k) =0.0 |
|---|
| 147 | asso4m_p_so2h2o2(i,k)=0.0 |
|---|
| 148 | |
|---|
| 149 | zrho=pplay(i,k)/(t(i,k)*RD) ! density of air in [kg/m^3] |
|---|
| 150 | mmr2moleqcm=fac1*zrho |
|---|
| 151 | |
|---|
| 152 | ! convert from mmr to molec cm-3 |
|---|
| 153 | #ifndef AERONLY |
|---|
| 154 | tr_seri(i,k,id_h2o2) =tr_seri(i,k,id_h2o2) *mmr2moleqcm/adv_mass(id_h2o2) |
|---|
| 155 | tr_seri(i,k,id_o3) =tr_seri(i,k,id_o3) *mmr2moleqcm/adv_mass(id_o3) |
|---|
| 156 | #endif |
|---|
| 157 | tr_seri(i,k,id_so2) =tr_seri(i,k,id_so2) *mmr2moleqcm/adv_mass(id_so2) |
|---|
| 158 | tr_seri(i,k,id_asso4m)=tr_seri(i,k,id_asso4m)*mmr2moleqcm/adv_mass(id_asso4m) |
|---|
| 159 | # if GRPCNT != 0 |
|---|
| 160 | tr_seri(i,k,id_ox) =tr_seri(i,k,id_ox) *mmr2moleqcm/nadv_mass(id_o3) |
|---|
| 161 | # endif |
|---|
| 162 | #ifndef AERONLY |
|---|
| 163 | zh2o2(i,k) =tr_seri(i,k,id_h2o2) |
|---|
| 164 | #else |
|---|
| 165 | zh2o2(i,k) =invariants(i,k,inv_H2O2) |
|---|
| 166 | #endif |
|---|
| 167 | zso2(i,k) =tr_seri(i,k,id_so2) |
|---|
| 168 | zso4(i,k) =tr_seri(i,k,id_asso4m) |
|---|
| 169 | #ifndef AERONLY |
|---|
| 170 | zo3(i,k) =tr_seri(i,k,id_o3) |
|---|
| 171 | #else |
|---|
| 172 | zo3(i,k) = invariants(i,k,inv_O3) |
|---|
| 173 | #endif |
|---|
| 174 | qliq2(i,k) =qliq*zrho*rho_water_inv*1.e-3 !--l cm-3 |
|---|
| 175 | |
|---|
| 176 | !--Henry's law constants |
|---|
| 177 | fac=1./101.325*R*t(i,k)*qliq*zrho*rho_water_inv |
|---|
| 178 | tfac=1./298.-1./t(i,k) |
|---|
| 179 | henry_t=henry_h2o2*EXP(-kk_h2o2*tfac) !--mol/l/atm |
|---|
| 180 | f_a_h2o2=henry_t*fac |
|---|
| 181 | scavh2o2(i,k)=f_a_h2o2/(1.+f_a_h2o2) |
|---|
| 182 | |
|---|
| 183 | henry_t=henry_o3*EXP(-kk_o3*tfac) !--mol/l/atm |
|---|
| 184 | f_a_o3=henry_t*fac |
|---|
| 185 | scavo3(i,k)=f_a_o3/(1.+f_a_o3) |
|---|
| 186 | |
|---|
| 187 | henry_t=henry_so2*EXP(-kk_so2*tfac) !--mol/l/atm |
|---|
| 188 | f_a(i,k)=henry_t*fac |
|---|
| 189 | |
|---|
| 190 | K1(i,k)=1.7e-2*EXP(-2090.*tfac) !-SO2/HSO3- |
|---|
| 191 | K2(i,k)=6.0e-8*EXP(-1120.*tfac) !-HSO3-/SO3= |
|---|
| 192 | |
|---|
| 193 | ENDDO |
|---|
| 194 | ENDDO |
|---|
| 195 | |
|---|
| 196 | DO pas=1, nbpas |
|---|
| 197 | DO k = 1, PLEV |
|---|
| 198 | DO i=1, PLON |
|---|
| 199 | so4_aq(i)=zso4(i,k)/qliq2(i,k)/RNAVO !--mol l-1 |
|---|
| 200 | hso3m_aq(i)=0.0 |
|---|
| 201 | so3mm_aq(i)=0.0 |
|---|
| 202 | no3m_aq(i)=0.0 |
|---|
| 203 | ENDDO |
|---|
| 204 | |
|---|
| 205 | DO iter=1, niter !--iterative calculation of pH |
|---|
| 206 | DO i = 1, PLON |
|---|
| 207 | |
|---|
| 208 | !--calculation of hplus -- see document aq.txt |
|---|
| 209 | xx=xx_clean |
|---|
| 210 | IF (so4_aq(i).GT.15.e-6) xx=xx_polluted |
|---|
| 211 | |
|---|
| 212 | bb=(2.-xx)*so4_aq(i)+2.*so3mm_aq(i)+hso3m_aq(i)+no3m_aq(i) |
|---|
| 213 | cc=KE |
|---|
| 214 | hplus(i)=0.5*(bb + SQRT(bb*bb+4*cc) ) |
|---|
| 215 | |
|---|
| 216 | correc=1.+K1(i,k)/hplus(i)*(1.+K2(i,k)/hplus(i)) |
|---|
| 217 | so2_aq(i)=f_a(i,k)*zso2(i,k)/RNAVO/qliq2(i,k)/ & |
|---|
| 218 | (1.+f_a(i,k)*correc) !--mol l-1 |
|---|
| 219 | hso3m_aq(i)=so2_aq(i)*K1(i,k)/hplus(i) !--mol l-1 |
|---|
| 220 | so3mm_aq(i)=hso3m_aq(i)*K2(i,k)/hplus(i) !--mol l-1 |
|---|
| 221 | |
|---|
| 222 | ENDDO |
|---|
| 223 | ENDDO |
|---|
| 224 | |
|---|
| 225 | DO i=1, PLON |
|---|
| 226 | !--calculation of pH |
|---|
| 227 | ph=-LOG(hplus(i))/LOG(10.) |
|---|
| 228 | ph_hist(i,k)=ph_hist(i,k)+ph/FLOAT(nbpas) |
|---|
| 229 | |
|---|
| 230 | !--calculation of other aqueous concentrations |
|---|
| 231 | h2o2_aq=scavh2o2(i,k)*zh2o2(i,k)/qliq2(i,k)/RNAVO !--mol l-1 |
|---|
| 232 | o3_aq=scavo3(i,k)*zo3(i,k)/qliq2(i,k)/RNAVO !--mol l-1 |
|---|
| 233 | ! |
|---|
| 234 | !--aqueous-phase chemistry |
|---|
| 235 | !--Reaction rates by Seinfeld and Pandis, Atmospheric Chemistry and Physics |
|---|
| 236 | !--page 363 and ff. |
|---|
| 237 | ! |
|---|
| 238 | tfac=1./298.-1./t(i,k) |
|---|
| 239 | prod1=7.5e7*EXP(4430.*tfac) !--M-2 s-1 |
|---|
| 240 | prod1=prod1*hplus(i)*h2o2_aq*hso3m_aq(i)/(1.+13.*hplus(i)) !--M s-1 |
|---|
| 241 | tend1=prod1*pdtaq !--M |
|---|
| 242 | tend1=tend1*qliq2(i,k)*RNAVO !--molec cm-3 |
|---|
| 243 | |
|---|
| 244 | prod2=2.4e4*so2_aq(i) !--s-1 |
|---|
| 245 | prod2=prod2+ 3.7e5*EXP(5530.*tfac)*hso3m_aq(i) !--s-1 |
|---|
| 246 | prod2=prod2+ 1.5e9*EXP(5280.*tfac)*so3mm_aq(i) !--s-1 |
|---|
| 247 | tend2=prod2*o3_aq*pdtaq !--M |
|---|
| 248 | tend2=tend2*qliq2(i,k)*RNAVO !--molec cm-3 |
|---|
| 249 | |
|---|
| 250 | !--limitation by temperature |
|---|
| 251 | tfactor=MIN(1.,MAX(0.,(t(i,k)-t0)/(t1-t0))) |
|---|
| 252 | tend1=tfactor*tend1 |
|---|
| 253 | tend2=tfactor*tend2 |
|---|
| 254 | |
|---|
| 255 | !--limitation by availability of ambient H2O2, O3, and SO2 |
|---|
| 256 | tend1=MIN(MIN(tend1,zh2o2(i,k)),zso2(i,k)) |
|---|
| 257 | tend2=MIN(MIN(tend2,zo3(i,k)),zso2(i,k)) |
|---|
| 258 | |
|---|
| 259 | !---do not oxidize too much |
|---|
| 260 | IF (tend1+tend2.GT.zso2(i,k)) THEN |
|---|
| 261 | tend1=tend1*zso2(i,k)/(tend1+tend2) |
|---|
| 262 | tend2=zso2(i,k)-tend1 |
|---|
| 263 | ENDIF |
|---|
| 264 | |
|---|
| 265 | !--S(IV)+H2O2 for history [molec/cm3/pdtphys] |
|---|
| 266 | asso4m_p_so2h2o2(i,k) = asso4m_p_so2h2o2(i,k) + tend1*rneb(i,k) |
|---|
| 267 | |
|---|
| 268 | !--S(IV)+O3 for history [molec/cm3/pdtphys] |
|---|
| 269 | asso4m_p_so2o3(i,k) = asso4m_p_so2o3(i,k) + tend2*rneb(i,k) |
|---|
| 270 | |
|---|
| 271 | !--update of concentration in the cloudy part |
|---|
| 272 | zh2o2(i,k)=zh2o2(i,k)-tend1 |
|---|
| 273 | zo3(i,k) =zo3(i,k) -tend2 |
|---|
| 274 | zso4(i,k) =zso4(i,k) +(tend1+tend2) |
|---|
| 275 | zso2(i,k) =zso2(i,k) -(tend1+tend2) |
|---|
| 276 | zso2(i,k) =MAX(zso2(i,k),0.0) !--au cas ou zso2 ~ 0 |
|---|
| 277 | |
|---|
| 278 | ENDDO !--fin i |
|---|
| 279 | ENDDO !--fin k |
|---|
| 280 | ENDDO !--fin pas |
|---|
| 281 | |
|---|
| 282 | !---update of total concentration (cloudy and clear) |
|---|
| 283 | ! and conversion of units from molecules/cm^3 to mmr |
|---|
| 284 | |
|---|
| 285 | DO k = 1, PLEV |
|---|
| 286 | DO i = 1, PLON |
|---|
| 287 | tr_seri(i,k,id_so2) =rneb(i,k)*zso2(i,k) + (1.-rneb(i,k))*tr_seri(i,k,id_so2) |
|---|
| 288 | tr_seri(i,k,id_asso4m) =rneb(i,k)*zso4(i,k) + (1.-rneb(i,k))*tr_seri(i,k,id_asso4m) |
|---|
| 289 | |
|---|
| 290 | zrho=pplay(i,k)/(t(i,k)*RD) |
|---|
| 291 | moleqcm2mmr=1./(zrho*fac1) |
|---|
| 292 | |
|---|
| 293 | ! convert from molec cm-3 to mmr |
|---|
| 294 | #ifndef AERONLY |
|---|
| 295 | tr_seri(i,k,id_h2o2) =tr_seri(i,k,id_h2o2) *moleqcm2mmr*adv_mass(id_h2o2) |
|---|
| 296 | tr_seri(i,k,id_o3) =tr_seri(i,k,id_o3) *moleqcm2mmr*adv_mass(id_o3) |
|---|
| 297 | #endif |
|---|
| 298 | tr_seri(i,k,id_so2) =tr_seri(i,k,id_so2) *moleqcm2mmr*adv_mass(id_so2) |
|---|
| 299 | tr_seri(i,k,id_asso4m)=tr_seri(i,k,id_asso4m)*moleqcm2mmr*adv_mass(id_asso4m) |
|---|
| 300 | # if GRPCNT != 0 |
|---|
| 301 | tr_seri(i,k,id_ox) =tr_seri(i,k,id_ox) *moleqcm2mmr*nadv_mass(id_o3) |
|---|
| 302 | # endif |
|---|
| 303 | |
|---|
| 304 | ! calculate chemical transfers per second |
|---|
| 305 | !--S(IV)+H2O2 for history [molec/cm3/s] |
|---|
| 306 | asso4m_p_so2h2o2(i,k)=asso4m_p_so2h2o2(i,k)/pdtphys |
|---|
| 307 | asso4m_p_so2o3(i,k) =asso4m_p_so2o3(i,k)/pdtphys |
|---|
| 308 | ENDDO |
|---|
| 309 | ENDDO |
|---|
| 310 | |
|---|
| 311 | RETURN |
|---|
| 312 | END SUBROUTINE chimieaq |
|---|
| 313 | #endif |
|---|
| 314 | #endif |
|---|
