1!------------------------------------------------------------------------------- 2 3! This file is part of Code_Saturne, a general-purpose CFD tool. 4! 5! Copyright (C) 1998-2021 EDF S.A. 6! 7! This program is free software; you can redistribute it and/or modify it under 8! the terms of the GNU General Public License as published by the Free Software 9! Foundation; either version 2 of the License, or (at your option) any later 10! version. 11! 12! This program is distributed in the hope that it will be useful, but WITHOUT 13! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 14! FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 15! details. 16! 17! You should have received a copy of the GNU General Public License along with 18! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin 19! Street, Fifth Floor, Boston, MA 02110-1301, USA. 20 21!------------------------------------------------------------------------------- 22 23!> \file chem_luscheme2.f90 24!> \brief Routines for atmospheric chemical scheme 2 25!> \remarks 26!> These routines are automatically generated by SPACK 27!> See CEREA: http://cerea.enpc.fr/polyphemus 28 29!> kinetic_2 30!> \brief Computation of kinetic rates for atmospheric chemistry 31!------------------------------------------------------------------------------ 32 33!------------------------------------------------------------------------------ 34! Arguments 35!------------------------------------------------------------------------------ 36! mode name role 37!------------------------------------------------------------------------------ 38!> \param[in] nr total number of chemical reactions 39!> \param[in] option_photolysis flag to activate or not photolysis reactions 40!> \param[in] azi solar zenith angle 41!> \param[in] att atmospheric attenuation variable 42!> \param[in] temp temperature 43!> \param[in] press pressure 44!> \param[in] xlw water massic fraction 45!> \param[out] rk(nr) kinetic rates 46!______________________________________________________________________________ 47 48subroutine kinetic_2(nr,rk,temp,xlw,press,azi,att, & 49 option_photolysis) 50 51implicit none 52 53! Arguments 54 55integer nr 56double precision rk(nr),temp,xlw,press 57double precision azi, att 58integer option_photolysis 59 60! Local variables 61 62double precision effko,rapk,summ 63double precision ylh2o 64 65! Compute third body. 66! Conversion = Avogadro*1d-6/Perfect gas constant. 67! PRESS in Pascal, SUMM in molecules/cm3, TEMP in Kelvin 68 69summ = press * 7.243d16 / temp 70 71! Number of water molecules computed from the massic fraction 72! (absolute humidity) 73 74ylh2o = 29.d0*summ*xlw/(18.d0+11.d0*xlw) 75 76! For the zenithal angle at tropics 77 78azi=abs(azi) 79 80rk( 1) = dexp(-0.8860689615829534d+02 & 81 - ( -0.5300000000000000d+03 )/temp) 82rk( 1) = rk( 1) * summ * 0.2d0 83rk( 2) = dexp(-0.2653240882726044d+02 & 84 - ( 0.1500000000000000d+04 )/temp) 85 86if(option_photolysis.eq.2) then 87 rk( 3)= 0.d0 88elseif(option_photolysis.eq.1) then 89if(azi.lt.0.d0)then 90 stop 91elseif(azi.ge. 0.00d+00 .and. azi.lt. 0.10d+02) then 92 rk( 3)=-0.1302720567168795d-07 93 rk( 3)=-0.7822279432831311d-06+(azi- 0.00d+00) * rk( 3) 94 rk( 3)= 0.0000000000000000d+00+(azi- 0.00d+00) * rk( 3) 95 rk( 3)= 0.9310260000000001d-02+(azi- 0.00d+00) * rk( 3) 96elseif(azi.ge. 0.10d+02 .and. azi.lt. 0.20d+02) then 97 rk( 3)= 0.3771617015067078d-08 98 rk( 3)=-0.1173044113433769d-05+(azi- 0.10d+02) * rk( 3) 99 rk( 3)=-0.1955272056716901d-04+(azi- 0.10d+02) * rk( 3) 100 rk( 3)= 0.9219010000000000d-02+(azi- 0.10d+02) * rk( 3) 101elseif(azi.ge. 0.20d+02 .and. azi.lt. 0.30d+02) then 102 rk( 3)=-0.5859262388581815d-08 103 rk( 3)=-0.1059895602981758d-05+(azi- 0.20d+02) * rk( 3) 104 rk( 3)=-0.4188211773132428d-04+(azi- 0.20d+02) * rk( 3) 105 rk( 3)= 0.8909950000000000d-02+(azi- 0.20d+02) * rk( 3) 106elseif(azi.ge. 0.30d+02 .and. azi.lt. 0.40d+02) then 107 rk( 3)=-0.7024567460738029d-08 108 rk( 3)=-0.1235673474639213d-05+(azi- 0.30d+02) * rk( 3) 109 rk( 3)=-0.6483780850753392d-04+(azi- 0.30d+02) * rk( 3) 110 rk( 3)= 0.8379279999999999d-02+(azi- 0.30d+02) * rk( 3) 111elseif(azi.ge. 0.40d+02 .and. azi.lt. 0.50d+02) then 112 rk( 3)=-0.9202467768466835d-08 113 rk( 3)=-0.1446410498461367d-05+(azi- 0.40d+02) * rk( 3) 114 rk( 3)=-0.9165864823853972d-04+(azi- 0.40d+02) * rk( 3) 115 rk( 3)= 0.7600310000000000d-02+(azi- 0.40d+02) * rk( 3) 116elseif(azi.ge. 0.50d+02 .and. azi.lt. 0.60d+02) then 117 rk( 3)=-0.1612556146540100d-07 118 rk( 3)=-0.1722484531515342d-05+(azi- 0.50d+02) * rk( 3) 119 rk( 3)=-0.1233475985383066d-03+(azi- 0.50d+02) * rk( 3) 120 rk( 3)= 0.6529880000000000d-02+(azi- 0.50d+02) * rk( 3) 121elseif(azi.ge. 0.60d+02 .and. azi.lt. 0.70d+02) then 122 rk( 3)= 0.3226471363007382d-07 123 rk( 3)=-0.2206251375477548d-05+(azi- 0.60d+02) * rk( 3) 124 rk( 3)=-0.1626349576082332d-03+(azi- 0.60d+02) * rk( 3) 125 rk( 3)= 0.5108030000000000d-02+(azi- 0.60d+02) * rk( 3) 126elseif(azi.ge. 0.70d+02 .and. azi.lt. 0.78d+02) then 127 rk( 3)= 0.2027078243961372d-06 128 rk( 3)=-0.1238309966574737d-05+(azi- 0.70d+02) * rk( 3) 129 rk( 3)=-0.1970805710287543d-03+(azi- 0.70d+02) * rk( 3) 130 rk( 3)= 0.3293320000000000d-02+(azi- 0.70d+02) * rk( 3) 131elseif(azi.ge. 0.78d+02 .and. azi.lt. 0.86d+02) then 132 rk( 3)=-0.7448311471194499d-07 133 rk( 3)= 0.3626677818932555d-05+(azi- 0.78d+02) * rk( 3) 134 rk( 3)=-0.1779736282099126d-03+(azi- 0.78d+02) * rk( 3) 135 rk( 3)= 0.1741210000000000d-02+(azi- 0.78d+02) * rk( 3) 136elseif(azi.ge. 0.86d+02 .and. azi.lt. 0.90d+02) then 137 rk( 3)= 0.2490309929270573d-05 138 rk( 3)= 0.1839083065842406d-05+(azi- 0.86d+02) * rk( 3) 139 rk( 3)=-0.1342475411316713d-03+(azi- 0.86d+02) * rk( 3) 140 rk( 3)= 0.5113930000000000d-03+(azi- 0.86d+02) * rk( 3) 141elseif(azi.ge.90.d0)then 142 rk( 3)= 0.1632080000000000d-03 143endif 144if(att.lt.0.99999) rk( 3) = rk( 3) * att 145endif 146 147rk( 4) = summ * 6.0d-34 * (temp/3.d2) ** (-2.4d0) 148rk( 4) = rk( 4) * summ * 0.2d0 149 150if(option_photolysis.eq.2) then 151 rk( 5)= 0.d0 152elseif(option_photolysis.eq.1) then 153if(azi.lt.0.d0)then 154 stop 155elseif(azi.ge. 0.00d+00 .and. azi.lt. 0.10d+02) then 156 rk( 5)=-0.5928286648807996d-09 157 rk( 5)=-0.3096171335119280d-07+(azi- 0.00d+00) * rk( 5) 158 rk( 5)= 0.0000000000000000d+00+(azi- 0.00d+00) * rk( 5) 159 rk( 5)= 0.4927580000000000d-03+(azi- 0.00d+00) * rk( 5) 160elseif(azi.ge. 0.10d+02 .and. azi.lt. 0.20d+02) then 161 rk( 5)= 0.1444859946426876d-09 162 rk( 5)=-0.4874657329761677d-07+(azi- 0.10d+02) * rk( 5) 163 rk( 5)=-0.7970828664880956d-06+(azi- 0.10d+02) * rk( 5) 164 rk( 5)= 0.4890690000000000d-03+(azi- 0.10d+02) * rk( 5) 165elseif(azi.ge. 0.20d+02 .and. azi.lt. 0.30d+02) then 166 rk( 5)=-0.5511531369012520d-10 167 rk( 5)=-0.4441199345833616d-07+(azi- 0.20d+02) * rk( 5) 168 rk( 5)=-0.1728668534047625d-05+(azi- 0.20d+02) * rk( 5) 169 rk( 5)= 0.4763680000000000d-03+(azi- 0.20d+02) * rk( 5) 170elseif(azi.ge. 0.30d+02 .and. azi.lt. 0.40d+02) then 171 rk( 5)=-0.3000247398821449d-09 172 rk( 5)=-0.4606545286904014d-07+(azi- 0.30d+02) * rk( 5) 173 rk( 5)=-0.2633442997321385d-05+(azi- 0.30d+02) * rk( 5) 174 rk( 5)= 0.4545850000000000d-03+(azi- 0.30d+02) * rk( 5) 175elseif(azi.ge. 0.40d+02 .and. azi.lt. 0.50d+02) then 176 rk( 5)=-0.2397857267812366d-09 177 rk( 5)=-0.5506619506550444d-07+(azi- 0.40d+02) * rk( 5) 178 rk( 5)=-0.3644759476666826d-05+(azi- 0.40d+02) * rk( 5) 179 rk( 5)= 0.4233440000000000d-03+(azi- 0.40d+02) * rk( 5) 180elseif(azi.ge. 0.50d+02 .and. azi.lt. 0.60d+02) then 181 rk( 5)=-0.1844832352993067d-08 182 rk( 5)=-0.6225976686893995d-07+(azi- 0.50d+02) * rk( 5) 183 rk( 5)=-0.4818019096011278d-05+(azi- 0.50d+02) * rk( 5) 184 rk( 5)= 0.3811500000000000d-03+(azi- 0.50d+02) * rk( 5) 185elseif(azi.ge. 0.60d+02 .and. azi.lt. 0.70d+02) then 186 rk( 5)= 0.1101151387530619d-09 187 rk( 5)=-0.1176047374587370d-06+(azi- 0.60d+02) * rk( 5) 188 rk( 5)=-0.6616664139287950d-05+(azi- 0.60d+02) * rk( 5) 189 rk( 5)= 0.3248990000000000d-03+(azi- 0.60d+02) * rk( 5) 190elseif(azi.ge. 0.70d+02 .and. azi.lt. 0.78d+02) then 191 rk( 5)=-0.1557211541866474d-07 192 rk( 5)=-0.1143012832961418d-06+(azi- 0.70d+02) * rk( 5) 193 rk( 5)=-0.8935724346837023d-05+(azi- 0.70d+02) * rk( 5) 194 rk( 5)= 0.2470820000000000d-03+(azi- 0.70d+02) * rk( 5) 195elseif(azi.ge. 0.78d+02 .and. azi.lt. 0.86d+02) then 196 rk( 5)= 0.4048472604232427d-07 197 rk( 5)=-0.4880320533439059d-06+(azi- 0.78d+02) * rk( 5) 198 rk( 5)=-0.1375439103995816d-04+(azi- 0.78d+02) * rk( 5) 199 rk( 5)= 0.1603080000000000d-03+(azi- 0.78d+02) * rk( 5) 200elseif(azi.ge. 0.86d+02 .and. azi.lt. 0.90d+02) then 201 rk( 5)= 0.2066880316654862d-06 202 rk( 5)= 0.4836013716715513d-06+(azi- 0.86d+02) * rk( 5) 203 rk( 5)=-0.1378983649333310d-04+(azi- 0.86d+02) * rk( 5) 204 rk( 5)= 0.3976700000000000d-04+(azi- 0.86d+02) * rk( 5) 205elseif(azi.ge.90.d0)then 206 rk( 5)= 0.5573310000000000d-05 207endif 208if(att.lt.0.99999) rk( 5) = rk( 5) * att 209endif 210 211 212if(option_photolysis.eq.2) then 213 rk( 6)= 0.d0 214elseif(option_photolysis.eq.1) then 215if(azi.lt.0.d0)then 216 stop 217elseif(azi.ge. 0.00d+00 .and. azi.lt. 0.10d+02) then 218 rk( 6)=-0.8776629099833464d-10 219 rk( 6)=-0.1165033709001661d-07+(azi- 0.00d+00) * rk( 6) 220 rk( 6)= 0.0000000000000000d+00+(azi- 0.00d+00) * rk( 6) 221 rk( 6)= 0.3523480000000000d-04+(azi- 0.00d+00) * rk( 6) 222elseif(azi.ge. 0.10d+02 .and. azi.lt. 0.20d+02) then 223 rk( 6)= 0.1474988729949909d-09 224 rk( 6)=-0.1428332581996665d-07+(azi- 0.10d+02) * rk( 6) 225 rk( 6)=-0.2593366290998327d-06+(azi- 0.10d+02) * rk( 6) 226 rk( 6)= 0.3398200000000000d-04+(azi- 0.10d+02) * rk( 6) 227elseif(azi.ge. 0.20d+02 .and. azi.lt. 0.30d+02) then 228 rk( 6)= 0.1300707990183827d-09 229 rk( 6)=-0.9858359630116941d-08+(azi- 0.20d+02) * rk( 6) 230 rk( 6)=-0.5007534836006686d-06+(azi- 0.20d+02) * rk( 6) 231 rk( 6)= 0.3010780000000000d-04+(azi- 0.20d+02) * rk( 6) 232elseif(azi.ge. 0.30d+02 .and. azi.lt. 0.40d+02) then 233 rk( 6)= 0.1988179309314732d-09 234 rk( 6)=-0.5956235659565481d-08+(azi- 0.30d+02) * rk( 6) 235 rk( 6)=-0.6588994364974921d-06+(azi- 0.30d+02) * rk( 6) 236 rk( 6)= 0.2424450000000000d-04+(azi- 0.30d+02) * rk( 6) 237elseif(azi.ge. 0.40d+02 .and. azi.lt. 0.50d+02) then 238 rk( 6)= 0.2219574772557277d-09 239 rk( 6)= 0.8302268378835231d-11+(azi- 0.40d+02) * rk( 6) 240 rk( 6)=-0.7183787704093613d-06+(azi- 0.40d+02) * rk( 6) 241 rk( 6)= 0.1725870000000000d-04+(azi- 0.40d+02) * rk( 6) 242elseif(azi.ge. 0.50d+02 .and. azi.lt. 0.60d+02) then 243 rk( 6)= 0.1913521600455895d-09 244 rk( 6)= 0.6667026586050136d-08+(azi- 0.50d+02) * rk( 6) 245 rk( 6)=-0.6516254818650674d-06+(azi- 0.50d+02) * rk( 6) 246 rk( 6)= 0.1029770000000000d-04+(azi- 0.50d+02) * rk( 6) 247elseif(azi.ge. 0.60d+02 .and. azi.lt. 0.70d+02) then 248 rk( 6)= 0.1602388256152816d-10 249 rk( 6)= 0.1240759138741867d-07+(azi- 0.60d+02) * rk( 6) 250 rk( 6)=-0.4608793021303539d-06+(azi- 0.60d+02) * rk( 6) 251 rk( 6)= 0.4639500000000000d-05+(azi- 0.60d+02) * rk( 6) 252elseif(azi.ge. 0.70d+02 .and. azi.lt. 0.78d+02) then 253 rk( 6)=-0.3089359890783960d-09 254 rk( 6)= 0.1288830786428400d-07+(azi- 0.70d+02) * rk( 6) 255 rk( 6)=-0.2079203096133002d-06+(azi- 0.70d+02) * rk( 6) 256 rk( 6)= 0.1287490000000000d-05+(azi- 0.70d+02) * rk( 6) 257elseif(azi.ge. 0.78d+02 .and. azi.lt. 0.86d+02) then 258 rk( 6)=-0.2034952628632162d-09 259 rk( 6)= 0.5473844126395715d-08+(azi- 0.78d+02) * rk( 6) 260 rk( 6)=-0.6102309368797090d-07+(azi- 0.78d+02) * rk( 6) 261 rk( 6)= 0.2908040000000000d-06+(azi- 0.78d+02) * rk( 6) 262elseif(azi.ge. 0.86d+02 .and. azi.lt. 0.90d+02) then 263 rk( 6)= 0.1623544916128788d-09 264 rk( 6)= 0.5899578175158973d-09+(azi- 0.86d+02) * rk( 6) 265 rk( 6)=-0.1251267813581064d-07+(azi- 0.86d+02) * rk( 6) 266 rk( 6)= 0.4875570000000000d-07+(azi- 0.86d+02) * rk( 6) 267elseif(azi.ge.90.d0)then 268 rk( 6)= 0.1853500000000000d-07 269endif 270if(att.lt.0.99999) rk( 6) = rk( 6) * att 271endif 272 273rk( 7) = dexp(-0.2458649867820512d+02 & 274 - ( -0.1020000000000000d+03 )/temp) 275rk( 7) = rk( 7) * summ 276rk( 8) = 0.2200000000000000d-09 277rk( 8) = rk( 8) * ylh2o 278 279if(option_photolysis.eq.2) then 280 rk( 9)= 0.d0 281elseif(option_photolysis.eq.1) then 282if(azi.lt.0.d0)then 283 stop 284elseif(azi.ge. 0.00d+00 .and. azi.lt. 0.10d+02) then 285 rk( 9)=-0.2887225450832658d-08 286 rk( 9)=-0.1810277454916759d-06+(azi- 0.00d+00) * rk( 9) 287 rk( 9)= 0.0000000000000000d+00+(azi- 0.00d+00) * rk( 9) 288 rk( 9)= 0.2046710000000000d-02+(azi- 0.00d+00) * rk( 9) 289elseif(azi.ge. 0.10d+02 .and. azi.lt. 0.20d+02) then 290 rk( 9)= 0.8216763524985595d-09 291 rk( 9)=-0.2676445090166556d-06+(azi- 0.10d+02) * rk( 9) 292 rk( 9)=-0.4486722545083316d-05+(azi- 0.10d+02) * rk( 9) 293 rk( 9)= 0.2025720000000000d-02+(azi- 0.10d+02) * rk( 9) 294elseif(azi.ge. 0.20d+02 .and. azi.lt. 0.30d+02) then 295 rk( 9)=-0.1309479959161308d-08 296 rk( 9)=-0.2429942184416991d-06+(azi- 0.20d+02) * rk( 9) 297 rk( 9)=-0.9593109819666860d-05+(azi- 0.20d+02) * rk( 9) 298 rk( 9)= 0.1954910000000000d-02+(azi- 0.20d+02) * rk( 9) 299elseif(azi.ge. 0.30d+02 .and. azi.lt. 0.40d+02) then 300 rk( 9)=-0.1523756515853649d-08 301 rk( 9)=-0.2822786172165394d-06+(azi- 0.30d+02) * rk( 9) 302 rk( 9)=-0.1484583817624924d-04+(azi- 0.30d+02) * rk( 9) 303 rk( 9)= 0.1833370000000000d-02+(azi- 0.30d+02) * rk( 9) 304elseif(azi.ge. 0.40d+02 .and. azi.lt. 0.50d+02) then 305 rk( 9)=-0.1745493977424016d-08 306 rk( 9)=-0.3279913126921461d-06+(azi- 0.40d+02) * rk( 9) 307 rk( 9)=-0.2094853747533609d-04+(azi- 0.40d+02) * rk( 9) 308 rk( 9)= 0.1655160000000000d-02+(azi- 0.40d+02) * rk( 9) 309elseif(azi.ge. 0.50d+02 .and. azi.lt. 0.60d+02) then 310 rk( 9)=-0.2754267574450560d-08 311 rk( 9)=-0.3803561320148624d-06+(azi- 0.50d+02) * rk( 9) 312 rk( 9)=-0.2803201192240625d-04+(azi- 0.50d+02) * rk( 9) 313 rk( 9)= 0.1411130000000000d-02+(azi- 0.50d+02) * rk( 9) 314elseif(azi.ge. 0.60d+02 .and. azi.lt. 0.70d+02) then 315 rk( 9)= 0.1037656427523324d-07 316 rk( 9)=-0.4629841592484096d-06+(azi- 0.60d+02) * rk( 9) 317 rk( 9)=-0.3646541483503878d-04+(azi- 0.60d+02) * rk( 9) 318 rk( 9)= 0.1090020000000000d-02+(azi- 0.60d+02) * rk( 9) 319elseif(azi.ge. 0.70d+02 .and. azi.lt. 0.78d+02) then 320 rk( 9)= 0.4335158727139456d-07 321 rk( 9)=-0.1516872309913989d-06+(azi- 0.70d+02) * rk( 9) 322 rk( 9)=-0.4261212873743828d-04+(azi- 0.70d+02) * rk( 9) 323 rk( 9)= 0.6894440000000000d-03+(azi- 0.70d+02) * rk( 9) 324elseif(azi.ge. 0.78d+02 .and. azi.lt. 0.86d+02) then 325 rk( 9)=-0.2976581610911796d-07 326 rk( 9)= 0.8887508635220705d-06+(azi- 0.78d+02) * rk( 9) 327 rk( 9)=-0.3671561967719464d-04+(azi- 0.78d+02) * rk( 9) 328 rk( 9)= 0.3610350000000000d-03+(azi- 0.78d+02) * rk( 9) 329elseif(azi.ge. 0.86d+02 .and. azi.lt. 0.90d+02) then 330 rk( 9)= 0.5586598403857848d-06 331 rk( 9)= 0.1743712769036732d-06+(azi- 0.86d+02) * rk( 9) 332 rk( 9)=-0.2821064255377481d-04+(azi- 0.86d+02) * rk( 9) 333 rk( 9)= 0.1089500000000000d-03+(azi- 0.86d+02) * rk( 9) 334elseif(azi.ge.90.d0)then 335 rk( 9)= 0.3465160000000000d-04 336endif 337if(att.lt.0.99999) rk( 9) = rk( 9) * att 338endif 339 340 341if(option_photolysis.eq.2) then 342 rk( 10)= 0.d0 343elseif(option_photolysis.eq.1) then 344if(azi.lt.0.d0)then 345 stop 346elseif(azi.ge. 0.00d+00 .and. azi.lt. 0.10d+02) then 347 rk( 10)=-0.1441479345432039d-10 348 rk( 10)=-0.1242452065456794d-08+(azi- 0.00d+00) * rk( 10) 349 rk( 10)= 0.0000000000000000d+00+(azi- 0.00d+00) * rk( 10) 350 rk( 10)= 0.8394580000000000d-05+(azi- 0.00d+00) * rk( 10) 351elseif(azi.ge. 0.10d+02 .and. azi.lt. 0.20d+02) then 352 rk( 10)= 0.8244380362960478d-11 353 rk( 10)=-0.1674895869086405d-08+(azi- 0.10d+02) * rk( 10) 354 rk( 10)=-0.2917347934543199d-07+(azi- 0.10d+02) * rk( 10) 355 rk( 10)= 0.8255920000000000d-05+(azi- 0.10d+02) * rk( 10) 356elseif(azi.ge. 0.20d+02 .and. azi.lt. 0.30d+02) then 357 rk( 10)= 0.1172720024787734d-12 358 rk( 10)=-0.1427564458197592d-08+(azi- 0.20d+02) * rk( 10) 359 rk( 10)=-0.6019808261827194d-07+(azi- 0.20d+02) * rk( 10) 360 rk( 10)= 0.7804940000000000d-05+(azi- 0.20d+02) * rk( 10) 361elseif(azi.ge. 0.30d+02 .and. azi.lt. 0.40d+02) then 362 rk( 10)= 0.4506531627127201d-11 363 rk( 10)=-0.1424046298123240d-08+(azi- 0.30d+02) * rk( 10) 364 rk( 10)=-0.8871419018148013d-07+(azi- 0.30d+02) * rk( 10) 365 rk( 10)= 0.7060320000000000d-05+(azi- 0.30d+02) * rk( 10) 366elseif(azi.ge. 0.40d+02 .and. azi.lt. 0.50d+02) then 367 rk( 10)= 0.1086660148900755d-10 368 rk( 10)=-0.1288850349309390d-08+(azi- 0.40d+02) * rk( 10) 369 rk( 10)=-0.1158431566558070d-06+(azi- 0.40d+02) * rk( 10) 370 rk( 10)= 0.6035280000000000d-05+(azi- 0.40d+02) * rk( 10) 371elseif(azi.ge. 0.50d+02 .and. azi.lt. 0.60d+02) then 372 rk( 10)= 0.1803706241686108d-10 373 rk( 10)=-0.9628523046392104d-09+(azi- 0.50d+02) * rk( 10) 374 rk( 10)=-0.1383601831952927d-06+(azi- 0.50d+02) * rk( 10) 375 rk( 10)= 0.4758830000000000d-05+(azi- 0.50d+02) * rk( 10) 376elseif(azi.ge. 0.60d+02 .and. azi.lt. 0.70d+02) then 377 rk( 10)= 0.7329514884355585d-10 378 rk( 10)=-0.4217404321336693d-09+(azi- 0.60d+02) * rk( 10) 379 rk( 10)=-0.1522061105630206d-06+(azi- 0.60d+02) * rk( 10) 380 rk( 10)= 0.3296980000000000d-05+(azi- 0.60d+02) * rk( 10) 381elseif(azi.ge. 0.70d+02 .and. azi.lt. 0.78d+02) then 382 rk( 10)= 0.6172011386267615d-10 383 rk( 10)= 0.1777114033174383d-08+(azi- 0.70d+02) * rk( 10) 384 rk( 10)=-0.1386523745526241d-06+(azi- 0.70d+02) * rk( 10) 385 rk( 10)= 0.1806040000000000d-05+(azi- 0.70d+02) * rk( 10) 386elseif(azi.ge. 0.78d+02 .and. azi.lt. 0.86d+02) then 387 rk( 10)=-0.8279704635596216d-10 388 rk( 10)= 0.3258396765877763d-08+(azi- 0.78d+02) * rk( 10) 389 rk( 10)=-0.9836828816022045d-07+(azi- 0.78d+02) * rk( 10) 390 rk( 10)= 0.8421570000000000d-06+(azi- 0.78d+02) * rk( 10) 391elseif(azi.ge. 0.86d+02 .and. azi.lt. 0.90d+02) then 392 rk( 10)= 0.1082517324572734d-08 393 rk( 10)= 0.1271267653334671d-08+(azi- 0.86d+02) * rk( 10) 394 rk( 10)=-0.6213097280649386d-07+(azi- 0.86d+02) * rk( 10) 395 rk( 10)= 0.2213560000000000d-06+(azi- 0.86d+02) * rk( 10) 396elseif(azi.ge.90.d0)then 397 rk( 10)= 0.6245350000000000d-07 398endif 399if(att.lt.0.99999) rk( 10) = rk( 10) * att 400endif 401 402 403if(option_photolysis.eq.2) then 404 rk( 11)= 0.d0 405elseif(option_photolysis.eq.1) then 406if(azi.lt.0.d0)then 407 stop 408elseif(azi.ge. 0.00d+00 .and. azi.lt. 0.10d+02) then 409 rk( 11)=-0.5816837387219519d-10 410 rk( 11)=-0.5184316261278087d-08+(azi- 0.00d+00) * rk( 11) 411 rk( 11)= 0.0000000000000000d+00+(azi- 0.00d+00) * rk( 11) 412 rk( 11)= 0.3220560000000000d-04+(azi- 0.00d+00) * rk( 11) 413elseif(azi.ge. 0.10d+02 .and. azi.lt. 0.20d+02) then 414 rk( 11)= 0.3450512161659949d-10 415 rk( 11)=-0.6929367477443941d-08+(azi- 0.10d+02) * rk( 11) 416 rk( 11)=-0.1211368373872203d-06+(azi- 0.10d+02) * rk( 11) 417 rk( 11)= 0.3162900000000000d-04+(azi- 0.10d+02) * rk( 11) 418elseif(azi.ge. 0.20d+02 .and. azi.lt. 0.30d+02) then 419 rk( 11)= 0.3647887405787883d-11 420 rk( 11)=-0.5894213828945960d-08+(azi- 0.20d+02) * rk( 11) 421 rk( 11)=-0.2493726504511192d-06+(azi- 0.20d+02) * rk( 11) 422 rk( 11)= 0.2975920000000000d-04+(azi- 0.20d+02) * rk( 11) 423elseif(azi.ge. 0.30d+02 .and. azi.lt. 0.40d+02) then 424 rk( 11)= 0.2530332876025029d-10 425 rk( 11)=-0.5784777206772343d-08+(azi- 0.30d+02) * rk( 11) 426 rk( 11)=-0.3661625608083018d-06+(azi- 0.30d+02) * rk( 11) 427 rk( 11)= 0.2667970000000000d-04+(azi- 0.30d+02) * rk( 11) 428elseif(azi.ge. 0.40d+02 .and. azi.lt. 0.50d+02) then 429 rk( 11)= 0.6013879755321277d-10 430 rk( 11)=-0.5025677343964861d-08+(azi- 0.40d+02) * rk( 11) 431 rk( 11)=-0.4742671063156733d-06+(azi- 0.40d+02) * rk( 11) 432 rk( 11)= 0.2246490000000000d-04+(azi- 0.40d+02) * rk( 11) 433elseif(azi.ge. 0.50d+02 .and. azi.lt. 0.60d+02) then 434 rk( 11)= 0.1004414810269593d-09 435 rk( 11)=-0.3221513417368319d-08+(azi- 0.50d+02) * rk( 11) 436 rk( 11)=-0.5567390139290089d-06+(azi- 0.50d+02) * rk( 11) 437 rk( 11)= 0.1727980000000000d-04+(azi- 0.50d+02) * rk( 11) 438elseif(azi.ge. 0.60d+02 .and. azi.lt. 0.70d+02) then 439 rk( 11)= 0.3107352783389442d-09 440 rk( 11)=-0.2082689865616575d-09+(azi- 0.60d+02) * rk( 11) 441 rk( 11)=-0.5910368379682989d-06+(azi- 0.60d+02) * rk( 11) 442 rk( 11)= 0.1149070000000000d-04+(azi- 0.60d+02) * rk( 11) 443elseif(azi.ge. 0.70d+02 .and. azi.lt. 0.78d+02) then 444 rk( 11)= 0.1345128013872392d-09 445 rk( 11)= 0.9113789363612173d-08+(azi- 0.70d+02) * rk( 11) 446 rk( 11)=-0.5019816341977565d-06+(azi- 0.70d+02) * rk( 11) 447 rk( 11)= 0.5870240000000000d-05+(azi- 0.70d+02) * rk( 11) 448elseif(azi.ge. 0.78d+02 .and. azi.lt. 0.86d+02) then 449 rk( 11)=-0.1960445509611151d-09 450 rk( 11)= 0.1234209659691269d-07+(azi- 0.78d+02) * rk( 11) 451 rk( 11)=-0.3303345465135304d-06+(azi- 0.78d+02) * rk( 11) 452 rk( 11)= 0.2506540000000000d-05+(azi- 0.78d+02) * rk( 11) 453elseif(azi.ge. 0.86d+02 .and. azi.lt. 0.90d+02) then 454 rk( 11)= 0.2279277828298341d-08 455 rk( 11)= 0.7637027373778166d-08+(azi- 0.86d+02) * rk( 11) 456 rk( 11)=-0.1705015547476783d-06+(azi- 0.86d+02) * rk( 11) 457 rk( 11)= 0.5533830000000000d-06+(azi- 0.86d+02) * rk( 11) 458elseif(azi.ge.90.d0)then 459 rk( 11)= 0.1394430000000000d-06 460endif 461if(att.lt.0.99999) rk( 11) = rk( 11) * att 462endif 463 464 465if(option_photolysis.eq.2) then 466 rk( 12)= 0.d0 467elseif(option_photolysis.eq.1) then 468if(azi.lt.0.d0)then 469 stop 470elseif(azi.ge. 0.00d+00 .and. azi.lt. 0.10d+02) then 471 rk( 12)=-0.7261045436751545d-10 472 rk( 12)=-0.5576895456324842d-08+(azi- 0.00d+00) * rk( 12) 473 rk( 12)= 0.0000000000000000d+00+(azi- 0.00d+00) * rk( 12) 474 rk( 12)= 0.4669460000000000d-04+(azi- 0.00d+00) * rk( 12) 475elseif(azi.ge. 0.10d+02 .and. azi.lt. 0.20d+02) then 476 rk( 12)= 0.2853136310254372d-10 477 rk( 12)=-0.7755209087350302d-08+(azi- 0.10d+02) * rk( 12) 478 rk( 12)=-0.1333210454367515d-06+(azi- 0.10d+02) * rk( 12) 479 rk( 12)= 0.4606430000000000d-04+(azi- 0.10d+02) * rk( 12) 480elseif(azi.ge. 0.20d+02 .and. azi.lt. 0.30d+02) then 481 rk( 12)=-0.1901499804265726d-10 482 rk( 12)=-0.6899268194273995d-08+(azi- 0.20d+02) * rk( 12) 483 rk( 12)=-0.2798658182529944d-06+(azi- 0.20d+02) * rk( 12) 484 rk( 12)= 0.4398410000000000d-04+(azi- 0.20d+02) * rk( 12) 485elseif(azi.ge. 0.30d+02 .and. azi.lt. 0.40d+02) then 486 rk( 12)=-0.9471370931908783d-11 487 rk( 12)=-0.7469718135553710d-08+(azi- 0.30d+02) * rk( 12) 488 rk( 12)=-0.4235556815512711d-06+(azi- 0.30d+02) * rk( 12) 489 rk( 12)= 0.4047650000000000d-04+(azi- 0.30d+02) * rk( 12) 490elseif(azi.ge. 0.40d+02 .and. azi.lt. 0.50d+02) then 491 rk( 12)= 0.1440048177028887d-10 492 rk( 12)=-0.7753859263510966d-08+(azi- 0.40d+02) * rk( 12) 493 rk( 12)=-0.5757914555419174d-06+(azi- 0.40d+02) * rk( 12) 494 rk( 12)= 0.3548450000000000d-04+(azi- 0.40d+02) * rk( 12) 495elseif(azi.ge. 0.50d+02 .and. azi.lt. 0.60d+02) then 496 rk( 12)= 0.3856944385081257d-10 497 rk( 12)=-0.7321844810402538d-08+(azi- 0.50d+02) * rk( 12) 498 rk( 12)=-0.7265484962810543d-06+(azi- 0.50d+02) * rk( 12) 499 rk( 12)= 0.2896560000000000d-04+(azi- 0.50d+02) * rk( 12) 500elseif(azi.ge. 0.60d+02 .and. azi.lt. 0.70d+02) then 501 rk( 12)= 0.4136217428262900d-09 502 rk( 12)=-0.6164761494878585d-08+(azi- 0.60d+02) * rk( 12) 503 rk( 12)=-0.8614145593338596d-06+(azi- 0.60d+02) * rk( 12) 504 rk( 12)= 0.2100650000000000d-04+(azi- 0.60d+02) * rk( 12) 505elseif(azi.ge. 0.70d+02 .and. azi.lt. 0.78d+02) then 506 rk( 12)= 0.5376701572532502d-09 507 rk( 12)= 0.6243890789914774d-08+(azi- 0.70d+02) * rk( 12) 508 rk( 12)=-0.8606232663835604d-06+(azi- 0.70d+02) * rk( 12) 509 rk( 12)= 0.1218950000000000d-04+(azi- 0.70d+02) * rk( 12) 510elseif(azi.ge. 0.78d+02 .and. azi.lt. 0.86d+02) then 511 rk( 12)=-0.5508273899935273d-09 512 rk( 12)= 0.1914797456399955d-07+(azi- 0.78d+02) * rk( 12) 513 rk( 12)=-0.6574883435524898d-06+(azi- 0.78d+02) * rk( 12) 514 rk( 12)= 0.5979410000000000d-05+(azi- 0.78d+02) * rk( 12) 515elseif(azi.ge. 0.86d+02 .and. azi.lt. 0.90d+02) then 516 rk( 12)= 0.8530305661922505d-08 517 rk( 12)= 0.5928117204100688d-08+(azi- 0.86d+02) * rk( 12) 518 rk( 12)=-0.4568796094072541d-06+(azi- 0.86d+02) * rk( 12) 519 rk( 12)= 0.1662950000000000d-05+(azi- 0.86d+02) * rk( 12) 520elseif(azi.ge.90.d0)then 521 rk( 12)= 0.4762210000000000d-06 522endif 523if(att.lt.0.99999) rk( 12) = rk( 12) * att 524endif 525 526 527if(option_photolysis.eq.2) then 528 rk( 13)= 0.d0 529elseif(option_photolysis.eq.1) then 530if(azi.lt.0.d0)then 531 stop 532elseif(azi.ge. 0.00d+00 .and. azi.lt. 0.10d+02) then 533 rk( 13)=-0.1292568102250478d-10 534 rk( 13)=-0.1349143189774952d-08+(azi- 0.00d+00) * rk( 13) 535 rk( 13)= 0.0000000000000000d+00+(azi- 0.00d+00) * rk( 13) 536 rk( 13)= 0.6105070000000000d-05+(azi- 0.00d+00) * rk( 13) 537elseif(azi.ge. 0.10d+02 .and. azi.lt. 0.20d+02) then 538 rk( 13)= 0.1221704306751377d-10 539 rk( 13)=-0.1736913620450095d-08+(azi- 0.10d+02) * rk( 13) 540 rk( 13)=-0.3086056810225046d-07+(azi- 0.10d+02) * rk( 13) 541 rk( 13)= 0.5957230000000000d-05+(azi- 0.10d+02) * rk( 13) 542elseif(azi.ge. 0.20d+02 .and. azi.lt. 0.30d+02) then 543 rk( 13)= 0.7227508752451222d-11 544 rk( 13)=-0.1370402328424685d-08+(azi- 0.20d+02) * rk( 13) 545 rk( 13)=-0.6193372759099823d-07+(azi- 0.20d+02) * rk( 13) 546 rk( 13)= 0.5487150000000000d-05+(azi- 0.20d+02) * rk( 13) 547elseif(azi.ge. 0.30d+02 .and. azi.lt. 0.40d+02) then 548 rk( 13)= 0.1521292192268060d-10 549 rk( 13)=-0.1153577065851153d-08+(azi- 0.30d+02) * rk( 13) 550 rk( 13)=-0.8717352153375648d-07+(azi- 0.30d+02) * rk( 13) 551 rk( 13)= 0.4738000000000000d-05+(azi- 0.30d+02) * rk( 13) 552elseif(azi.ge. 0.40d+02 .and. azi.lt. 0.50d+02) then 553 rk( 13)= 0.2301080355683006d-10 554 rk( 13)=-0.6971894081707449d-09+(azi- 0.40d+02) * rk( 13) 555 rk( 13)=-0.1056811862739754d-06+(azi- 0.40d+02) * rk( 13) 556 rk( 13)= 0.3766120000000000d-05+(azi- 0.40d+02) * rk( 13) 557elseif(azi.ge. 0.50d+02 .and. azi.lt. 0.60d+02) then 558 rk( 13)= 0.3106386384999131d-10 559 rk( 13)=-0.6865301465823343d-11+(azi- 0.50d+02) * rk( 13) 560 rk( 13)=-0.1127217333703412d-06+(azi- 0.50d+02) * rk( 13) 561 rk( 13)= 0.2662600000000000d-05+(azi- 0.50d+02) * rk( 13) 562elseif(azi.ge. 0.60d+02 .and. azi.lt. 0.70d+02) then 563 rk( 13)= 0.4531874104319860d-10 564 rk( 13)= 0.9250506140339690d-09+(azi- 0.60d+02) * rk( 13) 565 rk( 13)=-0.1035398802446585d-06+(azi- 0.60d+02) * rk( 13) 566 rk( 13)= 0.1565760000000000d-05+(azi- 0.60d+02) * rk( 13) 567elseif(azi.ge. 0.70d+02 .and. azi.lt. 0.78d+02) then 568 rk( 13)=-0.1303761111962380d-10 569 rk( 13)= 0.2284612845330880d-08+(azi- 0.70d+02) * rk( 13) 570 rk( 13)=-0.7144324565100488d-07+(azi- 0.70d+02) * rk( 13) 571 rk( 13)= 0.6681850000000000d-06+(azi- 0.70d+02) * rk( 13) 572elseif(azi.ge. 0.78d+02 .and. azi.lt. 0.86d+02) then 573 rk( 13)=-0.4077238229824927d-10 574 rk( 13)= 0.1971710178462450d-08+(azi- 0.78d+02) * rk( 13) 575 rk( 13)=-0.3739266146067857d-07+(azi- 0.78d+02) * rk( 13) 576 rk( 13)= 0.2361790000000000d-06+(azi- 0.78d+02) * rk( 13) 577elseif(azi.ge. 0.86d+02 .and. azi.lt. 0.90d+02) then 578 rk( 13)= 0.1193377495821847d-09 579 rk( 13)= 0.9931730033112436d-09+(azi- 0.86d+02) * rk( 13) 580 rk( 13)=-0.1367359600643481d-07+(azi- 0.86d+02) * rk( 13) 581 rk( 13)= 0.4235170000000000d-07+(azi- 0.86d+02) * rk( 13) 582elseif(azi.ge.90.d0)then 583 rk( 13)= 0.1118570000000000d-07 584endif 585if(att.lt.0.99999) rk( 13) = rk( 13) * att 586endif 587 588effko = 0.4900000000000000d-02* (temp / 3.d2) & 589 **(- ( 0.0000000000000000d+00))* & 590 dexp(- 0.1210000000000000d+05/temp) 591rapk = 0.5400000000000000d+17* (temp / 3.d2) & 592 **(- ( 0.0000000000000000d+00))* & 593 dexp(- 0.1383000000000000d+05/temp) 594rk( 14) = (effko*summ / ( 1.0d0 + effko * summ / & 595 rapk)) * 0.3000d+00** (1.0d0 / (1.0d0 + & 596 (dlog10(effko * summ / rapk))**2)) 597effko = 0.2700000000000000d-27* (temp / 3.d2) & 598 **(- ( 0.7100000000000000d+01)) 599rapk = 0.1200000000000000d-10* (temp / 3.d2) & 600 **(- ( 0.9000000000000000d+00)) 601rk( 15) = (effko * summ / & 602 ( 1.0d0 + effko * summ / rapk)) * & 603 0.3000d+00** (1.0d0 / (1.0d0 + & 604 (dlog10(effko * summ / rapk))**2)) 605rk( 16) = dexp(-0.2553915705425015d+02 & 606 - ( -0.2700000000000000d+03 )/temp) 607rk( 17) = dexp(-0.2637825814743318d+02 & 608 - ( -0.2500000000000000d+03 )/temp) 609effko = 0.7000000000000000d-30* (temp / 3.d2) & 610 **(- ( 0.2600000000000000d+01)) 611rapk = 0.3600000000000000d-10* (temp / 3.d2) & 612 **(- ( 0.1000000000000000d+00)) 613rk( 18) = (effko * summ / & 614 ( 1.0d0 + effko * summ / rapk)) * & 615 0.6000d+00** (1.0d0 / (1.0d0 + & 616 (dlog10(effko * summ / rapk))**2)) 617rk( 19) = dexp(-0.8322449114623726d+01 & 618 + (-0.3000000000000000d+01 * dlog(temp)) ) 619effko = 0.3300000000000000d-30* (temp / 3.d2) & 620 **(- ( 0.3300000000000000d+01)) 621rapk = 0.1500000000000000d-11* (temp / 3.d2) & 622 **(- ( 0.0000000000000000d+00)) 623rk( 20) = (effko * summ / & 624 ( 1.0d0 + effko * summ / rapk)) * & 625 0.6000d+00** (1.0d0 / (1.0d0 + & 626 (dlog10(effko * summ / rapk))**2)) 627rk( 21) = dexp(-0.2975128465212864d+02 & 628 - ( 0.2450000000000000d+04 )/temp) 629rk( 22) = dexp(-0.2492297091482634d+02 & 630 - ( -0.1700000000000000d+03 )/temp) 631rk( 23) = dexp(-0.3073211390514037d+02 & 632 - ( 0.1260000000000000d+04 )/temp) 633effko = 0.2000000000000000d-29* (temp / 3.d2) & 634 **(- ( 0.4400000000000000d+01)) 635rapk = 0.1400000000000000d-11* (temp / 3.d2) & 636 **(- ( 0.7000000000000000d+00)) 637rk( 24) = (effko * summ / & 638 ( 1.0d0 + effko * summ / rapk)) * & 639 0.6000d+00** (1.0d0 / (1.0d0 + & 640 (dlog10(effko * summ / rapk))**2)) 641effko = 0.1300000000000000d-02* (temp / 3.d2) & 642 **(- ( 0.3500000000000000d+01))* & 643 dexp(- 0.1100000000000000d+05/temp) 644rapk = 0.9700000000000000d+15* (temp / 3.d2) & 645 **(- (-0.1000000000000000d+00))* & 646 dexp(- 0.1108000000000000d+05/temp) 647rk( 25) = (effko*summ / ( 1.0d0 + effko * summ / & 648 rapk)) * 0.4500d+00** (1.0d0 / (1.0d0 + & 649 (dlog10(effko * summ / rapk))**2)) 650rk( 26) = 0.1000000000000000d-21 651rk( 26) = rk( 26) * ylh2o 652rk( 27) = dexp(-0.2667550967090111d+02 & 653 - ( -0.3650000000000000d+03 )/temp) 654rk( 28) = 0.6800000000000000d-13 655rk( 29) = 2.3d-13 * dexp(600.0d0 / temp) & 656 + 1.7d-33* summ * dexp(1000.0d0 / temp) 657rk( 30) = 3.22d-34 * dexp(2800.0d0 / temp) + & 658 2.38d-54 * summ * dexp(3200.0d0 / temp) 659rk( 30) = rk( 30) * ylh2o 660 661if(option_photolysis.eq.2) then 662 rk( 31)= 0.d0 663elseif(option_photolysis.eq.1) then 664if(azi.lt.0.d0)then 665 stop 666elseif(azi.ge. 0.00d+00 .and. azi.lt. 0.10d+02) then 667 rk( 31)=-0.2293082537254139d-06 668 rk( 31)=-0.1232691746274577d-04+(azi- 0.00d+00) * rk( 31) 669 rk( 31)= 0.0000000000000000d+00+(azi- 0.00d+00) * rk( 31) 670 rk( 31)= 0.2395610000000000d+00+(azi- 0.00d+00) * rk( 31) 671elseif(azi.ge. 0.10d+02 .and. azi.lt. 0.20d+02) then 672 rk( 31)=-0.2107523882380334d-07 673 rk( 31)=-0.1920616507450818d-04+(azi- 0.10d+02) * rk( 31) 674 rk( 31)=-0.3153308253725395d-03+(azi- 0.10d+02) * rk( 31) 675 rk( 31)= 0.2380990000000000d+00+(azi- 0.10d+02) * rk( 31) 676elseif(azi.ge. 0.20d+02 .and. azi.lt. 0.30d+02) then 677 rk( 31)= 0.2060920902067111d-07 678 rk( 31)=-0.1983842223922230d-04+(azi- 0.20d+02) * rk( 31) 679 rk( 31)=-0.7057766985098441d-03+(azi- 0.20d+02) * rk( 31) 680 rk( 31)= 0.2330040000000000d+00+(azi- 0.20d+02) * rk( 31) 681elseif(azi.ge. 0.30d+02 .and. azi.lt. 0.40d+02) then 682 rk( 31)=-0.2323615972588770d-06 683 rk( 31)=-0.1922014596860219d-04+(azi- 0.30d+02) * rk( 31) 684 rk( 31)=-0.1096362380588088d-02+(azi- 0.30d+02) * rk( 31) 685 rk( 31)= 0.2239830000000000d+00+(azi- 0.30d+02) * rk( 31) 686elseif(azi.ge. 0.40d+02 .and. azi.lt. 0.50d+02) then 687 rk( 31)=-0.2281628199851812d-06 688 rk( 31)=-0.2619099388636854d-04+(azi- 0.40d+02) * rk( 31) 689 rk( 31)=-0.1550473779137796d-02+(azi- 0.40d+02) * rk( 31) 690 rk( 31)= 0.2108650000000000d+00+(azi- 0.40d+02) * rk( 31) 691elseif(azi.ge. 0.50d+02 .and. azi.lt. 0.60d+02) then 692 rk( 31)=-0.1187987122800317d-05 693 rk( 31)=-0.3303587848592447d-04+(azi- 0.50d+02) * rk( 31) 694 rk( 31)=-0.2142742502860728d-02+(azi- 0.50d+02) * rk( 31) 695 rk( 31)= 0.1925130000000000d+00+(azi- 0.50d+02) * rk( 31) 696elseif(azi.ge. 0.60d+02 .and. azi.lt. 0.70d+02) then 697 rk( 31)=-0.3438886888129482d-06 698 rk( 31)=-0.6867549216993743d-04+(azi- 0.60d+02) * rk( 31) 699 rk( 31)=-0.3159856209419318d-02+(azi- 0.60d+02) * rk( 31) 700 rk( 31)= 0.1665940000000000d+00+(azi- 0.60d+02) * rk( 31) 701elseif(azi.ge. 0.70d+02 .and. azi.lt. 0.78d+02) then 702 rk( 31)=-0.7936212779117297d-05 703 rk( 31)=-0.7899215283432154d-04+(azi- 0.70d+02) * rk( 31) 704 rk( 31)=-0.4636532659461956d-02+(azi- 0.70d+02) * rk( 31) 705 rk( 31)= 0.1277840000000000d+00+(azi- 0.70d+02) * rk( 31) 706elseif(azi.ge. 0.78d+02 .and. azi.lt. 0.86d+02) then 707 rk( 31)= 0.2654257866666759d-04 708 rk( 31)=-0.2694612595331436d-03+(azi- 0.78d+02) * rk( 31) 709 rk( 31)=-0.7424159958401733d-02+(azi- 0.78d+02) * rk( 31) 710 rk( 31)= 0.8157290000000000d-01+(azi- 0.78d+02) * rk( 31) 711elseif(azi.ge. 0.86d+02 .and. azi.lt. 0.90d+02) then 712 rk( 31)= 0.7705999956653109d-04 713 rk( 31)= 0.3675606284668786d-03+(azi- 0.86d+02) * rk( 31) 714 rk( 31)=-0.6639365006930298d-02+(azi- 0.86d+02) * rk( 31) 715 rk( 31)= 0.1852390000000000d-01+(azi- 0.86d+02) * rk( 31) 716elseif(azi.ge.90.d0)then 717 rk( 31)= 0.2779250000000000d-02 718endif 719if(att.lt.0.99999) rk( 31) = rk( 31) * att 720endif 721 722rk( 32) = dexp(-0.2590825451818744d+02 & 723 - ( -0.1800000000000000d+03 )/temp) 724effko = 0.2500000000000000d-30* (temp / 3.d2) & 725 **(- ( 0.1800000000000000d+01)) 726rapk = 0.2200000000000000d-10* (temp / 3.d2) & 727 **(- ( 0.7000000000000000d+00)) 728rk( 33) = (effko * summ / & 729 ( 1.0d0 + effko * summ / rapk)) * & 730 0.6000d+00** (1.0d0 / (1.0d0 + & 731 (dlog10(effko * summ / rapk))**2)) 732 733if(option_photolysis.eq.2) then 734 rk( 34)= 0.d0 735elseif(option_photolysis.eq.1) then 736if(azi.lt.0.d0)then 737 stop 738elseif(azi.ge. 0.00d+00 .and. azi.lt. 0.10d+02) then 739 rk( 34)=-0.1726633761595217d-06 740 rk( 34)=-0.3763226238404795d-05+(azi- 0.00d+00) * rk( 34) 741 rk( 34)= 0.0000000000000000d+00+(azi- 0.00d+00) * rk( 34) 742 rk( 34)= 0.2758968400000000d-01+(azi- 0.00d+00) * rk( 34) 743elseif(azi.ge. 0.10d+02 .and. azi.lt. 0.20d+02) then 744 rk( 34)= 0.8277361284785716d-06 745 rk( 34)=-0.8943127523190445d-05+(azi- 0.10d+02) * rk( 34) 746 rk( 34)=-0.1270635376159524d-03+(azi- 0.10d+02) * rk( 34) 747 rk( 34)= 0.2704069800000000d-01+(azi- 0.10d+02) * rk( 34) 748elseif(azi.ge. 0.20d+02 .and. azi.lt. 0.30d+02) then 749 rk( 34)=-0.1363361137754772d-05 750 rk( 34)= 0.1588895633116670d-04+(azi- 0.20d+02) * rk( 34) 751 rk( 34)=-0.5760524953618990d-04+(azi- 0.20d+02) * rk( 34) 752 rk( 34)= 0.2570348600000000d-01+(azi- 0.20d+02) * rk( 34) 753elseif(azi.ge. 0.30d+02 .and. azi.lt. 0.40d+02) then 754 rk( 34)= 0.9462274225405240d-06 755 rk( 34)=-0.2501187780147645d-04+(azi- 0.30d+02) * rk( 34) 756 rk( 34)=-0.1488344642392872d-03+(azi- 0.30d+02) * rk( 34) 757 rk( 34)= 0.2535296800000000d-01+(azi- 0.30d+02) * rk( 34) 758elseif(azi.ge. 0.40d+02 .and. azi.lt. 0.50d+02) then 759 rk( 34)=-0.2183585524073524d-06 760 rk( 34)= 0.3374944874739267d-05+(azi- 0.40d+02) * rk( 34) 761 rk( 34)=-0.3652037935066590d-03+(azi- 0.40d+02) * rk( 34) 762 rk( 34)= 0.2230966300000000d-01+(azi- 0.40d+02) * rk( 34) 763elseif(azi.ge. 0.50d+02 .and. azi.lt. 0.60d+02) then 764 rk( 34)= 0.1378647870888786d-06 765 rk( 34)=-0.3175811697481157d-05+(azi- 0.50d+02) * rk( 34) 766 rk( 34)=-0.3632124617340762d-03+(azi- 0.50d+02) * rk( 34) 767 rk( 34)= 0.1877676100000000d-01+(azi- 0.50d+02) * rk( 34) 768elseif(azi.ge. 0.60d+02 .and. azi.lt. 0.70d+02) then 769 rk( 34)=-0.3865635959481289d-06 770 rk( 34)= 0.9601319151840079d-06+(azi- 0.60d+02) * rk( 34) 771 rk( 34)=-0.3853692595570321d-03+(azi- 0.60d+02) * rk( 34) 772 rk( 34)= 0.1496492000000000d-01+(azi- 0.60d+02) * rk( 34) 773elseif(azi.ge. 0.70d+02 .and. azi.lt. 0.78d+02) then 774 rk( 34)=-0.1916508555648358d-06 775 rk( 34)=-0.1063677596325682d-04+(azi- 0.70d+02) * rk( 34) 776 rk( 34)=-0.4821357000377863d-03+(azi- 0.70d+02) * rk( 34) 777 rk( 34)= 0.1082067700000000d-01+(azi- 0.70d+02) * rk( 34) 778elseif(azi.ge. 0.78d+02 .and. azi.lt. 0.86d+02) then 779 rk( 34)= 0.2540478756918318d-05 780 rk( 34)=-0.1523639649680941d-04+(azi- 0.78d+02) * rk( 34) 781 rk( 34)=-0.6891210797183578d-03+(azi- 0.78d+02) * rk( 34) 782 rk( 34)= 0.6184712500000000d-02+(azi- 0.78d+02) * rk( 34) 783elseif(azi.ge. 0.86d+02 .and. azi.lt. 0.90d+02) then 784 rk( 34)= 0.1651057353835306d-05 785 rk( 34)= 0.4573509366928574d-04+(azi- 0.86d+02) * rk( 34) 786 rk( 34)=-0.4451315023386027d-03+(azi- 0.86d+02) * rk( 34) 787 rk( 34)= 0.9973396100000000d-03+(azi- 0.86d+02) * rk( 34) 788elseif(azi.ge.90.d0)then 789 rk( 34)= 0.5424277000000000d-04 790endif 791if(att.lt.0.99999) rk( 34) = rk( 34) * att 792endif 793 794 795return 796end subroutine kinetic_2 797 798 799!=============================================================================== 800 801!> fexchem_2 802!> \brief Computation of the chemical production terms 803!------------------------------------------------------------------------------ 804 805!------------------------------------------------------------------------------ 806! Arguments 807!------------------------------------------------------------------------------ 808! mode name role 809!------------------------------------------------------------------------------ 810!> \param[in] ns total number of chemical species 811!> \param[in] nr total number of chemical reactions 812!> \param[in] y concentrations vector 813!> \param[in] rk kinetic rates 814!> \param[in] zcsourc source term 815!> \param[in] convers_factor conversion factors 816!> \param[out] chem chemical production terms for every species 817!______________________________________________________________________________ 818 819subroutine fexchem_2(ns,nr,y,rk,zcsourc,convers_factor,chem) 820 821implicit none 822 823! Arguments 824 825integer nr,ns 826double precision rk(nr),y(ns),chem(ns),zcsourc(ns) 827double precision convers_factor(ns) 828 829! Local variables 830 831integer i 832double precision w(nr) 833double precision conc(ns) 834 835do i=1,ns 836 chem(i)=0.d0 837enddo 838 839! Conversion mg/kg to molecules/cm3. 840 841do i = 1, ns 842 conc(i) = y(i) * convers_factor(i) 843enddo 844 845! Compute reaction rates. 846 847call rates_2(ns,nr,rk,conc,w) 848 849! Chemical production terms. 850 851chem( 19) = chem( 19) + 0.2000000000000000d+01 * w( 1) 852chem( 20) = chem( 20) - 0.2000000000000000d+01 * w( 1) 853chem( 16) = chem( 16) - w( 2) 854chem( 19) = chem( 19) + w( 2) 855chem( 20) = chem( 20) - w( 2) 856chem( 17) = chem( 17) + w( 3) 857chem( 19) = chem( 19) - w( 3) 858chem( 20) = chem( 20) + w( 3) 859chem( 16) = chem( 16) + w( 4) 860chem( 17) = chem( 17) - w( 4) 861chem( 16) = chem( 16) - w( 5) 862chem( 17) = chem( 17) + w( 5) 863chem( 2) = chem( 2) + w( 6) 864chem( 16) = chem( 16) - w( 6) 865chem( 2) = chem( 2) - w( 7) 866chem( 17) = chem( 17) + w( 7) 867chem( 2) = chem( 2) - w( 8) 868chem( 15) = chem( 15) + 0.2000000000000000d+01 * w( 8) 869chem( 8) = chem( 8) - w( 9) 870chem( 15) = chem( 15) + w( 9) 871chem( 20) = chem( 20) + w( 9) 872chem( 3) = chem( 3) - w( 10) 873chem( 15) = chem( 15) + 0.2000000000000000d+01 * w( 10) 874chem( 4) = chem( 4) + w( 11) 875chem( 10) = chem( 10) - w( 11) 876chem( 14) = chem( 14) + 0.2000000000000000d+01 * w( 11) 877chem( 4) = chem( 4) + w( 12) 878chem( 10) = chem( 10) - w( 12) 879chem( 1) = chem( 1) - w( 13) 880chem( 4) = chem( 4) + w( 13) 881chem( 10) = chem( 10) + w( 13) 882chem( 13) = chem( 13) + w( 13) 883chem( 14) = chem( 14) + 0.2000000000000000d+01 * w( 13) 884chem( 11) = chem( 11) - w( 14) 885chem( 12) = chem( 12) + w( 14) 886chem( 19) = chem( 19) + w( 14) 887chem( 11) = chem( 11) + w( 15) 888chem( 12) = chem( 12) - w( 15) 889chem( 19) = chem( 19) - w( 15) 890chem( 10) = chem( 10) + w( 16) 891chem( 12) = chem( 12) - w( 16) 892chem( 13) = chem( 13) + w( 16) 893chem( 14) = chem( 14) + w( 16) 894chem( 19) = chem( 19) + w( 16) 895chem( 20) = chem( 20) - w( 16) 896chem( 14) = chem( 14) - w( 17) 897chem( 15) = chem( 15) + w( 17) 898chem( 19) = chem( 19) + w( 17) 899chem( 20) = chem( 20) - w( 17) 900chem( 8) = chem( 8) + w( 18) 901chem( 15) = chem( 15) - w( 18) 902chem( 20) = chem( 20) - w( 18) 903chem( 9) = chem( 9) + w( 19) 904chem( 15) = chem( 15) - w( 19) 905chem( 19) = chem( 19) - w( 19) 906chem( 5) = chem( 5) - w( 20) 907chem( 6) = chem( 6) + w( 20) 908chem( 14) = chem( 14) + w( 20) 909chem( 15) = chem( 15) - w( 20) 910chem( 16) = chem( 16) - w( 21) 911chem( 18) = chem( 18) + w( 21) 912chem( 19) = chem( 19) - w( 21) 913chem( 18) = chem( 18) - w( 22) 914chem( 19) = chem( 19) + 0.2000000000000000d+01 * w( 22) 915chem( 20) = chem( 20) - w( 22) 916chem( 18) = chem( 18) - w( 23) 917chem( 20) = chem( 20) + w( 23) 918chem( 7) = chem( 7) + w( 24) 919chem( 18) = chem( 18) - w( 24) 920chem( 19) = chem( 19) - w( 24) 921chem( 7) = chem( 7) - w( 25) 922chem( 18) = chem( 18) + w( 25) 923chem( 19) = chem( 19) + w( 25) 924chem( 7) = chem( 7) - w( 26) 925chem( 9) = chem( 9) + 0.2000000000000000d+01 * w( 26) 926chem( 13) = chem( 13) - w( 27) 927chem( 19) = chem( 19) + w( 27) 928chem( 20) = chem( 20) - w( 27) 929chem( 13) = chem( 13) - 0.2000000000000000d+01 * w( 28) 930chem( 3) = chem( 3) + w( 29) 931chem( 14) = chem( 14) - 0.2000000000000000d+01 * w( 29) 932chem( 3) = chem( 3) + w( 30) 933chem( 14) = chem( 14) - 0.2000000000000000d+01 * w( 30) 934chem( 17) = chem( 17) + 0.8900000000000000d+00 * w( 31) 935chem( 18) = chem( 18) - w( 31) 936chem( 19) = chem( 19) + 0.8900000000000000d+00 * w( 31) 937chem( 20) = chem( 20) + 0.1100000000000000d+00 * w( 31) 938chem( 17) = chem( 17) - w( 32) 939chem( 19) = chem( 19) - w( 32) 940chem( 20) = chem( 20) + w( 32) 941chem( 17) = chem( 17) - w( 33) 942chem( 18) = chem( 18) + w( 33) 943chem( 19) = chem( 19) - w( 33) 944chem( 7) = chem( 7) - w( 34) 945chem( 18) = chem( 18) + w( 34) 946chem( 19) = chem( 19) + w( 34) 947 948! Conversion molecules/cm3/sec to mg/kg/sec. 949 950do i = 1, ns 951 chem(i) = chem(i) / convers_factor(i) 952enddo 953 954! Volumic source terms. 955 956do i=1,ns 957chem(i)=chem(i)+zcsourc(i) 958enddo 959 960return 961end subroutine fexchem_2 962 963 964!=============================================================================== 965 966!> jacdchemdc_2 967!> \brief Computation of the Jacobian matrix for atmospheric chemistry 968!------------------------------------------------------------------------------ 969 970!------------------------------------------------------------------------------ 971! Arguments 972!------------------------------------------------------------------------------ 973! mode name role 974!------------------------------------------------------------------------------ 975!> \param[in] ns total number of chemical species 976!> \param[in] nr total number of chemical reactions 977!> \param[in] y concentrations vector 978!> \param[in] convers_factor conversion factors of mug/m3 to 979!> molecules/cm3 980!> \param[in] convers_factor_jac conversion factors for the Jacobian matrix 981!> (Wmol(i)/Wmol(j)) 982!> \param[in] rk kinetic rates 983!> \param[out] jacc Jacobian matrix 984!______________________________________________________________________________ 985 986subroutine jacdchemdc_2(ns,nr,y,convers_factor, & 987 convers_factor_jac,rk,jacc) 988 989implicit none 990 991! Arguments 992 993integer nr,ns 994double precision rk(nr),y(ns),jacc(ns,ns) 995double precision convers_factor(ns) 996double precision convers_factor_jac(ns,ns) 997 998! Local variables 999 1000integer i,j 1001double precision dw(nr,ns) 1002double precision conc(ns) 1003 1004do j=1,ns 1005 do i=1,ns 1006 jacc(i,j)=0.d0 1007 enddo 1008enddo 1009 1010! Conversion mg/kg to molecules/cm3. 1011 1012do i = 1, ns 1013 conc(i) = y(i) * convers_factor(i) 1014enddo 1015 1016call dratedc_2(ns,nr,rk,conc,dw) 1017 1018jacc( 19, 20) = jacc( 19, 20)+ 0.2000000000000000d+01*dw( 1, 20) 1019jacc( 19, 20) = jacc( 19, 20)+ 0.2000000000000000d+01*dw( 1, 20) 1020jacc( 20, 20) = jacc( 20, 20)- 0.2000000000000000d+01*dw( 1, 20) 1021jacc( 20, 20) = jacc( 20, 20)- 0.2000000000000000d+01*dw( 1, 20) 1022jacc( 16, 16) = jacc( 16, 16) - dw( 2, 16) 1023jacc( 16, 20) = jacc( 16, 20) - dw( 2, 20) 1024jacc( 19, 16) = jacc( 19, 16) + dw( 2, 16) 1025jacc( 19, 20) = jacc( 19, 20) + dw( 2, 20) 1026jacc( 20, 16) = jacc( 20, 16) - dw( 2, 16) 1027jacc( 20, 20) = jacc( 20, 20) - dw( 2, 20) 1028jacc( 17, 19) = jacc( 17, 19) + dw( 3, 19) 1029jacc( 19, 19) = jacc( 19, 19) - dw( 3, 19) 1030jacc( 20, 19) = jacc( 20, 19) + dw( 3, 19) 1031jacc( 16, 17) = jacc( 16, 17) + dw( 4, 17) 1032jacc( 17, 17) = jacc( 17, 17) - dw( 4, 17) 1033jacc( 16, 16) = jacc( 16, 16) - dw( 5, 16) 1034jacc( 17, 16) = jacc( 17, 16) + dw( 5, 16) 1035jacc( 2, 16) = jacc( 2, 16) + dw( 6, 16) 1036jacc( 16, 16) = jacc( 16, 16) - dw( 6, 16) 1037jacc( 2, 2) = jacc( 2, 2) - dw( 7, 2) 1038jacc( 17, 2) = jacc( 17, 2) + dw( 7, 2) 1039jacc( 2, 2) = jacc( 2, 2) - dw( 8, 2) 1040jacc( 15, 2) = jacc( 15, 2)+ 0.2000000000000000d+01*dw( 8, 2) 1041jacc( 8, 8) = jacc( 8, 8) - dw( 9, 8) 1042jacc( 15, 8) = jacc( 15, 8) + dw( 9, 8) 1043jacc( 20, 8) = jacc( 20, 8) + dw( 9, 8) 1044jacc( 3, 3) = jacc( 3, 3) - dw( 10, 3) 1045jacc( 15, 3) = jacc( 15, 3)+ 0.2000000000000000d+01*dw( 10, 3) 1046jacc( 4, 10) = jacc( 4, 10) + dw( 11, 10) 1047jacc( 10, 10) = jacc( 10, 10) - dw( 11, 10) 1048jacc( 14, 10) = jacc( 14, 10)+ 0.2000000000000000d+01*dw( 11, 10) 1049jacc( 4, 10) = jacc( 4, 10) + dw( 12, 10) 1050jacc( 10, 10) = jacc( 10, 10) - dw( 12, 10) 1051jacc( 1, 1) = jacc( 1, 1) - dw( 13, 1) 1052jacc( 4, 1) = jacc( 4, 1) + dw( 13, 1) 1053jacc( 10, 1) = jacc( 10, 1) + dw( 13, 1) 1054jacc( 13, 1) = jacc( 13, 1) + dw( 13, 1) 1055jacc( 14, 1) = jacc( 14, 1)+ 0.2000000000000000d+01*dw( 13, 1) 1056jacc( 11, 11) = jacc( 11, 11) - dw( 14, 11) 1057jacc( 12, 11) = jacc( 12, 11) + dw( 14, 11) 1058jacc( 19, 11) = jacc( 19, 11) + dw( 14, 11) 1059jacc( 11, 12) = jacc( 11, 12) + dw( 15, 12) 1060jacc( 11, 19) = jacc( 11, 19) + dw( 15, 19) 1061jacc( 12, 12) = jacc( 12, 12) - dw( 15, 12) 1062jacc( 12, 19) = jacc( 12, 19) - dw( 15, 19) 1063jacc( 19, 12) = jacc( 19, 12) - dw( 15, 12) 1064jacc( 19, 19) = jacc( 19, 19) - dw( 15, 19) 1065jacc( 10, 12) = jacc( 10, 12) + dw( 16, 12) 1066jacc( 10, 20) = jacc( 10, 20) + dw( 16, 20) 1067jacc( 12, 12) = jacc( 12, 12) - dw( 16, 12) 1068jacc( 12, 20) = jacc( 12, 20) - dw( 16, 20) 1069jacc( 13, 12) = jacc( 13, 12) + dw( 16, 12) 1070jacc( 13, 20) = jacc( 13, 20) + dw( 16, 20) 1071jacc( 14, 12) = jacc( 14, 12) + dw( 16, 12) 1072jacc( 14, 20) = jacc( 14, 20) + dw( 16, 20) 1073jacc( 19, 12) = jacc( 19, 12) + dw( 16, 12) 1074jacc( 19, 20) = jacc( 19, 20) + dw( 16, 20) 1075jacc( 20, 12) = jacc( 20, 12) - dw( 16, 12) 1076jacc( 20, 20) = jacc( 20, 20) - dw( 16, 20) 1077jacc( 14, 14) = jacc( 14, 14) - dw( 17, 14) 1078jacc( 14, 20) = jacc( 14, 20) - dw( 17, 20) 1079jacc( 15, 14) = jacc( 15, 14) + dw( 17, 14) 1080jacc( 15, 20) = jacc( 15, 20) + dw( 17, 20) 1081jacc( 19, 14) = jacc( 19, 14) + dw( 17, 14) 1082jacc( 19, 20) = jacc( 19, 20) + dw( 17, 20) 1083jacc( 20, 14) = jacc( 20, 14) - dw( 17, 14) 1084jacc( 20, 20) = jacc( 20, 20) - dw( 17, 20) 1085jacc( 8, 15) = jacc( 8, 15) + dw( 18, 15) 1086jacc( 8, 20) = jacc( 8, 20) + dw( 18, 20) 1087jacc( 15, 15) = jacc( 15, 15) - dw( 18, 15) 1088jacc( 15, 20) = jacc( 15, 20) - dw( 18, 20) 1089jacc( 20, 15) = jacc( 20, 15) - dw( 18, 15) 1090jacc( 20, 20) = jacc( 20, 20) - dw( 18, 20) 1091jacc( 9, 19) = jacc( 9, 19) + dw( 19, 19) 1092jacc( 9, 15) = jacc( 9, 15) + dw( 19, 15) 1093jacc( 15, 19) = jacc( 15, 19) - dw( 19, 19) 1094jacc( 15, 15) = jacc( 15, 15) - dw( 19, 15) 1095jacc( 19, 19) = jacc( 19, 19) - dw( 19, 19) 1096jacc( 19, 15) = jacc( 19, 15) - dw( 19, 15) 1097jacc( 5, 5) = jacc( 5, 5) - dw( 20, 5) 1098jacc( 5, 15) = jacc( 5, 15) - dw( 20, 15) 1099jacc( 6, 5) = jacc( 6, 5) + dw( 20, 5) 1100jacc( 6, 15) = jacc( 6, 15) + dw( 20, 15) 1101jacc( 14, 5) = jacc( 14, 5) + dw( 20, 5) 1102jacc( 14, 15) = jacc( 14, 15) + dw( 20, 15) 1103jacc( 15, 5) = jacc( 15, 5) - dw( 20, 5) 1104jacc( 15, 15) = jacc( 15, 15) - dw( 20, 15) 1105jacc( 16, 19) = jacc( 16, 19) - dw( 21, 19) 1106jacc( 16, 16) = jacc( 16, 16) - dw( 21, 16) 1107jacc( 18, 19) = jacc( 18, 19) + dw( 21, 19) 1108jacc( 18, 16) = jacc( 18, 16) + dw( 21, 16) 1109jacc( 19, 19) = jacc( 19, 19) - dw( 21, 19) 1110jacc( 19, 16) = jacc( 19, 16) - dw( 21, 16) 1111jacc( 18, 18) = jacc( 18, 18) - dw( 22, 18) 1112jacc( 18, 20) = jacc( 18, 20) - dw( 22, 20) 1113jacc( 19, 18) = jacc( 19, 18)+ 0.2000000000000000d+01*dw( 22, 18) 1114jacc( 19, 20) = jacc( 19, 20)+ 0.2000000000000000d+01*dw( 22, 20) 1115jacc( 20, 18) = jacc( 20, 18) - dw( 22, 18) 1116jacc( 20, 20) = jacc( 20, 20) - dw( 22, 20) 1117jacc( 18, 18) = jacc( 18, 18) - dw( 23, 18) 1118jacc( 18, 19) = jacc( 18, 19) - dw( 23, 19) 1119jacc( 20, 18) = jacc( 20, 18) + dw( 23, 18) 1120jacc( 20, 19) = jacc( 20, 19) + dw( 23, 19) 1121jacc( 7, 18) = jacc( 7, 18) + dw( 24, 18) 1122jacc( 7, 19) = jacc( 7, 19) + dw( 24, 19) 1123jacc( 18, 18) = jacc( 18, 18) - dw( 24, 18) 1124jacc( 18, 19) = jacc( 18, 19) - dw( 24, 19) 1125jacc( 19, 18) = jacc( 19, 18) - dw( 24, 18) 1126jacc( 19, 19) = jacc( 19, 19) - dw( 24, 19) 1127jacc( 7, 7) = jacc( 7, 7) - dw( 25, 7) 1128jacc( 18, 7) = jacc( 18, 7) + dw( 25, 7) 1129jacc( 19, 7) = jacc( 19, 7) + dw( 25, 7) 1130jacc( 7, 7) = jacc( 7, 7) - dw( 26, 7) 1131jacc( 9, 7) = jacc( 9, 7)+ 0.2000000000000000d+01*dw( 26, 7) 1132jacc( 13, 13) = jacc( 13, 13) - dw( 27, 13) 1133jacc( 13, 20) = jacc( 13, 20) - dw( 27, 20) 1134jacc( 19, 13) = jacc( 19, 13) + dw( 27, 13) 1135jacc( 19, 20) = jacc( 19, 20) + dw( 27, 20) 1136jacc( 20, 13) = jacc( 20, 13) - dw( 27, 13) 1137jacc( 20, 20) = jacc( 20, 20) - dw( 27, 20) 1138jacc( 13, 13) = jacc( 13, 13)- 0.2000000000000000d+01*dw( 28, 13) 1139jacc( 13, 13) = jacc( 13, 13)- 0.2000000000000000d+01*dw( 28, 13) 1140jacc( 3, 14) = jacc( 3, 14) + dw( 29, 14) 1141jacc( 3, 14) = jacc( 3, 14) + dw( 29, 14) 1142jacc( 14, 14) = jacc( 14, 14)- 0.2000000000000000d+01*dw( 29, 14) 1143jacc( 14, 14) = jacc( 14, 14)- 0.2000000000000000d+01*dw( 29, 14) 1144jacc( 3, 14) = jacc( 3, 14) + dw( 30, 14) 1145jacc( 3, 14) = jacc( 3, 14) + dw( 30, 14) 1146jacc( 14, 14) = jacc( 14, 14)- 0.2000000000000000d+01*dw( 30, 14) 1147jacc( 14, 14) = jacc( 14, 14)- 0.2000000000000000d+01*dw( 30, 14) 1148jacc( 17, 18) = jacc( 17, 18)+ 0.8900000000000000d+00*dw( 31, 18) 1149jacc( 18, 18) = jacc( 18, 18) - dw( 31, 18) 1150jacc( 19, 18) = jacc( 19, 18)+ 0.8900000000000000d+00*dw( 31, 18) 1151jacc( 20, 18) = jacc( 20, 18)+ 0.1100000000000000d+00*dw( 31, 18) 1152jacc( 17, 17) = jacc( 17, 17) - dw( 32, 17) 1153jacc( 17, 19) = jacc( 17, 19) - dw( 32, 19) 1154jacc( 19, 17) = jacc( 19, 17) - dw( 32, 17) 1155jacc( 19, 19) = jacc( 19, 19) - dw( 32, 19) 1156jacc( 20, 17) = jacc( 20, 17) + dw( 32, 17) 1157jacc( 20, 19) = jacc( 20, 19) + dw( 32, 19) 1158jacc( 17, 17) = jacc( 17, 17) - dw( 33, 17) 1159jacc( 17, 19) = jacc( 17, 19) - dw( 33, 19) 1160jacc( 18, 17) = jacc( 18, 17) + dw( 33, 17) 1161jacc( 18, 19) = jacc( 18, 19) + dw( 33, 19) 1162jacc( 19, 17) = jacc( 19, 17) - dw( 33, 17) 1163jacc( 19, 19) = jacc( 19, 19) - dw( 33, 19) 1164jacc( 7, 7) = jacc( 7, 7) - dw( 34, 7) 1165jacc( 18, 7) = jacc( 18, 7) + dw( 34, 7) 1166jacc( 19, 7) = jacc( 19, 7) + dw( 34, 7) 1167 1168do j = 1, ns 1169 do i = 1, ns 1170 jacc(i,j) = jacc(i,j) * convers_factor_jac(i,j) 1171 enddo 1172enddo 1173 1174 1175return 1176end subroutine jacdchemdc_2 1177 1178 1179!=============================================================================== 1180 1181!> rates_2 1182!> \brief Computation of reaction rates 1183!------------------------------------------------------------------------------ 1184 1185!------------------------------------------------------------------------------ 1186! Arguments 1187!------------------------------------------------------------------------------ 1188! mode name role 1189!------------------------------------------------------------------------------ 1190!> \param[in] ns total number of chemical species 1191!> \param[in] nr total number of chemical reactions 1192!> \param[in] rk kinetic rates 1193!> \param[in] y concentrations vector 1194!> \param[out] w reaction rates 1195!______________________________________________________________________________ 1196 1197subroutine rates_2(ns,nr,rk,y,w) 1198 1199implicit none 1200 1201! Arguments 1202 1203integer nr,ns 1204double precision rk(nr),y(ns) 1205 1206! Local variables 1207 1208double precision w(nr) 1209 1210w( 1) = rk( 1) * y( 20) * y( 20) 1211w( 2) = rk( 2) * y( 16) * y( 20) 1212w( 3) = rk( 3) * y( 19) 1213w( 4) = rk( 4) * y( 17) 1214w( 5) = rk( 5) * y( 16) 1215w( 6) = rk( 6) * y( 16) 1216w( 7) = rk( 7) * y( 2) 1217w( 8) = rk( 8) * y( 2) 1218w( 9) = rk( 9) * y( 8) 1219w( 10) = rk( 10) * y( 3) 1220w( 11) = rk( 11) * y( 10) 1221w( 12) = rk( 12) * y( 10) 1222w( 13) = rk( 13) * y( 1) 1223w( 14) = rk( 14) * y( 11) 1224w( 15) = rk( 15) * y( 12) * y( 19) 1225w( 16) = rk( 16) * y( 12) * y( 20) 1226w( 17) = rk( 17) * y( 14) * y( 20) 1227w( 18) = rk( 18) * y( 15) * y( 20) 1228w( 19) = rk( 19) * y( 19) * y( 15) 1229w( 20) = rk( 20) * y( 5) * y( 15) 1230w( 21) = rk( 21) * y( 19) * y( 16) 1231w( 22) = rk( 22) * y( 18) * y( 20) 1232w( 23) = rk( 23) * y( 18) * y( 19) 1233w( 24) = rk( 24) * y( 18) * y( 19) 1234w( 25) = rk( 25) * y( 7) 1235w( 26) = rk( 26) * y( 7) 1236w( 27) = rk( 27) * y( 13) * y( 20) 1237w( 28) = rk( 28) * y( 13) * y( 13) 1238w( 29) = rk( 29) * y( 14) * y( 14) 1239w( 30) = rk( 30) * y( 14) * y( 14) 1240w( 31) = rk( 31) * y( 18) 1241w( 32) = rk( 32) * y( 17) * y( 19) 1242w( 33) = rk( 33) * y( 17) * y( 19) 1243w( 34) = rk( 34) * y( 7) 1244 1245return 1246end subroutine rates_2 1247 1248 1249!=============================================================================== 1250 1251!> dratedc_2 1252!> \brief Computation of derivatives of reaction rates 1253!------------------------------------------------------------------------------ 1254 1255!------------------------------------------------------------------------------ 1256! Arguments 1257!------------------------------------------------------------------------------ 1258! mode name role 1259!------------------------------------------------------------------------------ 1260!> \param[in] nr total number of chemical reactions 1261!> \param[in] ns total number of chemical species 1262!> \param[in] rk kinetic rates 1263!> \param[in] y concentrations vector 1264!> \param[out] dw derivatives of reaction rates 1265!______________________________________________________________________________ 1266 1267subroutine dratedc_2(ns,nr,rk,y,dw) 1268 1269implicit none 1270 1271! Arguments 1272 1273integer nr,ns 1274double precision rk(nr),y(ns) 1275 1276! Local variables 1277 1278double precision dw(nr,ns) 1279 1280dw( 1, 20) = rk( 1) * y( 20) 1281dw( 1, 20) = rk( 1) * y( 20) 1282dw( 2, 16) = rk( 2) * y( 20) 1283dw( 2, 20) = rk( 2) * y( 16) 1284dw( 3, 19) = rk( 3) 1285dw( 4, 17) = rk( 4) 1286dw( 5, 16) = rk( 5) 1287dw( 6, 16) = rk( 6) 1288dw( 7, 2) = rk( 7) 1289dw( 8, 2) = rk( 8) 1290dw( 9, 8) = rk( 9) 1291dw( 10, 3) = rk( 10) 1292dw( 11, 10) = rk( 11) 1293dw( 12, 10) = rk( 12) 1294dw( 13, 1) = rk( 13) 1295dw( 14, 11) = rk( 14) 1296dw( 15, 12) = rk( 15) * y( 19) 1297dw( 15, 19) = rk( 15) * y( 12) 1298dw( 16, 12) = rk( 16) * y( 20) 1299dw( 16, 20) = rk( 16) * y( 12) 1300dw( 17, 14) = rk( 17) * y( 20) 1301dw( 17, 20) = rk( 17) * y( 14) 1302dw( 18, 15) = rk( 18) * y( 20) 1303dw( 18, 20) = rk( 18) * y( 15) 1304dw( 19, 19) = rk( 19) * y( 15) 1305dw( 19, 15) = rk( 19) * y( 19) 1306dw( 20, 5) = rk( 20) * y( 15) 1307dw( 20, 15) = rk( 20) * y( 5) 1308dw( 21, 19) = rk( 21) * y( 16) 1309dw( 21, 16) = rk( 21) * y( 19) 1310dw( 22, 18) = rk( 22) * y( 20) 1311dw( 22, 20) = rk( 22) * y( 18) 1312dw( 23, 18) = rk( 23) * y( 19) 1313dw( 23, 19) = rk( 23) * y( 18) 1314dw( 24, 18) = rk( 24) * y( 19) 1315dw( 24, 19) = rk( 24) * y( 18) 1316dw( 25, 7) = rk( 25) 1317dw( 26, 7) = rk( 26) 1318dw( 27, 13) = rk( 27) * y( 20) 1319dw( 27, 20) = rk( 27) * y( 13) 1320dw( 28, 13) = rk( 28) * y( 13) 1321dw( 28, 13) = rk( 28) * y( 13) 1322dw( 29, 14) = rk( 29) * y( 14) 1323dw( 29, 14) = rk( 29) * y( 14) 1324dw( 30, 14) = rk( 30) * y( 14) 1325dw( 30, 14) = rk( 30) * y( 14) 1326dw( 31, 18) = rk( 31) 1327dw( 32, 17) = rk( 32) * y( 19) 1328dw( 32, 19) = rk( 32) * y( 17) 1329dw( 33, 17) = rk( 33) * y( 19) 1330dw( 33, 19) = rk( 33) * y( 17) 1331dw( 34, 7) = rk( 34) 1332 1333return 1334end subroutine dratedc_2 1335 1336 1337!=============================================================================== 1338 1339!> lu_decompose_2 1340!> \brief Computation of LU factorization of matrix m 1341!------------------------------------------------------------------------------ 1342 1343!------------------------------------------------------------------------------ 1344! Arguments 1345!------------------------------------------------------------------------------ 1346! mode name role 1347!------------------------------------------------------------------------------ 1348!> \param[in] ns matrix row number from the chemical species 1349!> number 1350!> \param[in,out] m on entry, an invertible matrix. 1351!> On exit, an LU factorization of m 1352!______________________________________________________________________________ 1353 1354subroutine lu_decompose_2 (ns,m) 1355 1356implicit none 1357 1358! Arguments 1359 1360integer ns 1361double precision m(ns,ns) 1362 1363! Local variables 1364 1365double precision temp 1366 1367! Upper part. 1368m(2, 16) = m(2, 16) / m(2, 2) 1369 1370! Upper part. 1371m(3, 14) = m(3, 14) / m(3, 3) 1372 1373! Upper part. 1374m(4, 10) = m(4, 10) / m(4, 4) 1375 1376! Upper part. 1377m(5, 15) = m(5, 15) / m(5, 5) 1378 1379! Upper part. 1380temp = m(6, 5) * m(5, 15) 1381m(6, 15) = ( m(6, 15) - temp ) / m(6, 6) 1382 1383! Upper part. 1384m(7, 18) = m(7, 18) / m(7, 7) 1385! Upper part. 1386m(7, 19) = m(7, 19) / m(7, 7) 1387 1388! Upper part. 1389m(8, 15) = m(8, 15) / m(8, 8) 1390! Upper part. 1391m(8, 20) = m(8, 20) / m(8, 8) 1392 1393! Upper part. 1394m(9, 15) = m(9, 15) / m(9, 9) 1395! Upper part. 1396temp = m(9, 7) * m(7, 18) 1397m(9, 18) = ( m(9, 18) - temp ) / m(9, 9) 1398! Upper part. 1399temp = m(9, 7) * m(7, 19) 1400m(9, 19) = ( m(9, 19) - temp ) / m(9, 9) 1401 1402! Upper part. 1403m(10, 12) = m(10, 12) / m(10, 10) 1404! Upper part. 1405m(10, 20) = m(10, 20) / m(10, 10) 1406 1407! Upper part. 1408m(11, 12) = m(11, 12) / m(11, 11) 1409! Upper part. 1410m(11, 19) = m(11, 19) / m(11, 11) 1411 1412! Lower part. 1413temp = m(12, 11) * m(11, 12) 1414m(12, 12) = m(12, 12) - temp 1415! Lower part. 1416temp = m(14, 10) * m(10, 12) 1417m(14, 12) = m(14, 12) - temp 1418! Lower part. 1419temp = m(19, 11) * m(11, 12) 1420m(19, 12) = m(19, 12) - temp 1421! Upper part. 1422temp = m(12, 11) * m(11, 19) 1423m(12, 19) = ( m(12, 19) - temp ) / m(12, 12) 1424! Upper part. 1425m(12, 20) = m(12, 20) / m(12, 12) 1426 1427! Upper part. 1428temp = m(13, 12) * m(12, 19) 1429m(13, 19) = ( m(13, 19) - temp ) / m(13, 13) 1430! Upper part. 1431temp = m(13, 12) * m(12, 20) 1432m(13, 20) = ( m(13, 20) - temp ) / m(13, 13) 1433 1434! Lower part. 1435temp = m(15, 3) * m(3, 14) 1436m(15, 14) = m(15, 14) - temp 1437! Upper part. 1438temp = m(14, 5) * m(5, 15) 1439m(14, 15) = ( m(14, 15) - temp ) / m(14, 14) 1440! Upper part. 1441temp = m(14, 12) * m(12, 19) 1442m(14, 19) = ( m(14, 19) - temp ) / m(14, 14) 1443! Upper part. 1444temp = m(14, 10) * m(10, 20) 1445temp = temp + m(14, 12) * m(12, 20) 1446m(14, 20) = ( m(14, 20) - temp ) / m(14, 14) 1447 1448! Lower part. 1449temp = m(15, 5) * m(5, 15) 1450temp = temp + m(15, 8) * m(8, 15) 1451temp = temp + m(15, 14) * m(14, 15) 1452m(15, 15) = m(15, 15) - temp 1453! Lower part. 1454temp = m(19, 14) * m(14, 15) 1455m(19, 15) = m(19, 15) - temp 1456! Lower part. 1457temp = m(20, 8) * m(8, 15) 1458temp = temp + m(20, 14) * m(14, 15) 1459m(20, 15) = m(20, 15) - temp 1460! Upper part. 1461temp = m(15, 2) * m(2, 16) 1462m(15, 16) = ( m(15, 16) - temp ) / m(15, 15) 1463! Upper part. 1464temp = m(15, 14) * m(14, 19) 1465m(15, 19) = ( m(15, 19) - temp ) / m(15, 15) 1466! Upper part. 1467temp = m(15, 8) * m(8, 20) 1468temp = temp + m(15, 14) * m(14, 20) 1469m(15, 20) = ( m(15, 20) - temp ) / m(15, 15) 1470 1471! Lower part. 1472temp = m(17, 2) * m(2, 16) 1473m(17, 16) = m(17, 16) - temp 1474! Lower part. 1475temp = m(19, 15) * m(15, 16) 1476m(19, 16) = m(19, 16) - temp 1477! Lower part. 1478temp = m(20, 15) * m(15, 16) 1479m(20, 16) = m(20, 16) - temp 1480! Upper part. 1481m(16, 17) = m(16, 17) / m(16, 16) 1482! Upper part. 1483m(16, 19) = m(16, 19) / m(16, 16) 1484! Upper part. 1485m(16, 20) = m(16, 20) / m(16, 16) 1486 1487! Lower part. 1488temp = m(17, 16) * m(16, 17) 1489m(17, 17) = m(17, 17) - temp 1490! Lower part. 1491temp = m(18, 16) * m(16, 17) 1492m(18, 17) = m(18, 17) - temp 1493! Lower part. 1494temp = m(19, 16) * m(16, 17) 1495m(19, 17) = m(19, 17) - temp 1496! Lower part. 1497temp = m(20, 16) * m(16, 17) 1498m(20, 17) = m(20, 17) - temp 1499! Upper part. 1500m(17, 18) = m(17, 18) / m(17, 17) 1501! Upper part. 1502temp = m(17, 16) * m(16, 19) 1503m(17, 19) = ( m(17, 19) - temp ) / m(17, 17) 1504! Upper part. 1505temp = m(17, 16) * m(16, 20) 1506m(17, 20) = ( m(17, 20) - temp ) / m(17, 17) 1507 1508! Lower part. 1509temp = m(18, 7) * m(7, 18) 1510temp = temp + m(18, 17) * m(17, 18) 1511m(18, 18) = m(18, 18) - temp 1512! Lower part. 1513temp = m(19, 7) * m(7, 18) 1514temp = temp + m(19, 17) * m(17, 18) 1515m(19, 18) = m(19, 18) - temp 1516! Lower part. 1517temp = m(20, 17) * m(17, 18) 1518m(20, 18) = m(20, 18) - temp 1519! Upper part. 1520temp = m(18, 7) * m(7, 19) 1521temp = temp + m(18, 16) * m(16, 19) 1522temp = temp + m(18, 17) * m(17, 19) 1523m(18, 19) = ( m(18, 19) - temp ) / m(18, 18) 1524! Upper part. 1525temp = m(18, 16) * m(16, 20) 1526temp = temp + m(18, 17) * m(17, 20) 1527m(18, 20) = ( m(18, 20) - temp ) / m(18, 18) 1528 1529! Lower part. 1530temp = m(19, 7) * m(7, 19) 1531temp = temp + m(19, 11) * m(11, 19) 1532temp = temp + m(19, 12) * m(12, 19) 1533temp = temp + m(19, 13) * m(13, 19) 1534temp = temp + m(19, 14) * m(14, 19) 1535temp = temp + m(19, 15) * m(15, 19) 1536temp = temp + m(19, 16) * m(16, 19) 1537temp = temp + m(19, 17) * m(17, 19) 1538temp = temp + m(19, 18) * m(18, 19) 1539m(19, 19) = m(19, 19) - temp 1540! Lower part. 1541temp = m(20, 12) * m(12, 19) 1542temp = temp + m(20, 13) * m(13, 19) 1543temp = temp + m(20, 14) * m(14, 19) 1544temp = temp + m(20, 15) * m(15, 19) 1545temp = temp + m(20, 16) * m(16, 19) 1546temp = temp + m(20, 17) * m(17, 19) 1547temp = temp + m(20, 18) * m(18, 19) 1548m(20, 19) = m(20, 19) - temp 1549! Upper part. 1550temp = m(19, 12) * m(12, 20) 1551temp = temp + m(19, 13) * m(13, 20) 1552temp = temp + m(19, 14) * m(14, 20) 1553temp = temp + m(19, 15) * m(15, 20) 1554temp = temp + m(19, 16) * m(16, 20) 1555temp = temp + m(19, 17) * m(17, 20) 1556temp = temp + m(19, 18) * m(18, 20) 1557m(19, 20) = ( m(19, 20) - temp ) / m(19, 19) 1558 1559! Lower part. 1560temp = m(20, 8) * m(8, 20) 1561temp = temp + m(20, 12) * m(12, 20) 1562temp = temp + m(20, 13) * m(13, 20) 1563temp = temp + m(20, 14) * m(14, 20) 1564temp = temp + m(20, 15) * m(15, 20) 1565temp = temp + m(20, 16) * m(16, 20) 1566temp = temp + m(20, 17) * m(17, 20) 1567temp = temp + m(20, 18) * m(18, 20) 1568temp = temp + m(20, 19) * m(19, 20) 1569m(20, 20) = m(20, 20) - temp 1570 1571return 1572end subroutine lu_decompose_2 1573 1574 1575!=============================================================================== 1576 1577!> lu_solve_2 1578!> \brief Resolution of MY=X where M is an LU factorization computed by 1579!> lu_decompose_2 1580!------------------------------------------------------------------------------ 1581 1582!------------------------------------------------------------------------------ 1583! Arguments 1584!------------------------------------------------------------------------------ 1585! mode name role 1586!------------------------------------------------------------------------------ 1587!> \param[in] ns matrix row number from the chemical species number 1588!> \param[in] m an LU factorization computed by lu_decompose_2 1589!> \param[in,out] x on entry, the right-hand side of the equation. 1590! on exit, the solution of the equation 1591!______________________________________________________________________________ 1592 1593subroutine lu_solve_2 (ns, m, x) 1594 1595implicit none 1596 1597! Arguments 1598 1599integer ns 1600double precision m(ns,ns) 1601double precision x(ns) 1602 1603! Local variables 1604 1605double precision temp 1606 1607! Forward substitution. 1608 1609x(1) = x(1) / m(1, 1) 1610 1611x(2) = x(2) / m(2, 2) 1612 1613x(3) = x(3) / m(3, 3) 1614 1615temp = m(4, 1) * x(1) 1616x(4) = ( x(4) - temp ) / m(4, 4) 1617 1618x(5) = x(5) / m(5, 5) 1619 1620temp = m(6, 5) * x(5) 1621x(6) = ( x(6) - temp ) / m(6, 6) 1622 1623x(7) = x(7) / m(7, 7) 1624 1625x(8) = x(8) / m(8, 8) 1626 1627temp = m(9, 7) * x(7) 1628x(9) = ( x(9) - temp ) / m(9, 9) 1629 1630temp = m(10, 1) * x(1) 1631x(10) = ( x(10) - temp ) / m(10, 10) 1632 1633x(11) = x(11) / m(11, 11) 1634 1635temp = m(12, 11) * x(11) 1636x(12) = ( x(12) - temp ) / m(12, 12) 1637 1638temp = m(13, 1) * x(1) 1639temp = temp + m(13, 12) * x(12) 1640x(13) = ( x(13) - temp ) / m(13, 13) 1641 1642temp = m(14, 1) * x(1) 1643temp = temp + m(14, 5) * x(5) 1644temp = temp + m(14, 10) * x(10) 1645temp = temp + m(14, 12) * x(12) 1646x(14) = ( x(14) - temp ) / m(14, 14) 1647 1648temp = m(15, 2) * x(2) 1649temp = temp + m(15, 3) * x(3) 1650temp = temp + m(15, 5) * x(5) 1651temp = temp + m(15, 8) * x(8) 1652temp = temp + m(15, 14) * x(14) 1653x(15) = ( x(15) - temp ) / m(15, 15) 1654 1655x(16) = x(16) / m(16, 16) 1656 1657temp = m(17, 2) * x(2) 1658temp = temp + m(17, 16) * x(16) 1659x(17) = ( x(17) - temp ) / m(17, 17) 1660 1661temp = m(18, 7) * x(7) 1662temp = temp + m(18, 16) * x(16) 1663temp = temp + m(18, 17) * x(17) 1664x(18) = ( x(18) - temp ) / m(18, 18) 1665 1666temp = m(19, 7) * x(7) 1667temp = temp + m(19, 11) * x(11) 1668temp = temp + m(19, 12) * x(12) 1669temp = temp + m(19, 13) * x(13) 1670temp = temp + m(19, 14) * x(14) 1671temp = temp + m(19, 15) * x(15) 1672temp = temp + m(19, 16) * x(16) 1673temp = temp + m(19, 17) * x(17) 1674temp = temp + m(19, 18) * x(18) 1675x(19) = ( x(19) - temp ) / m(19, 19) 1676 1677temp = m(20, 8) * x(8) 1678temp = temp + m(20, 12) * x(12) 1679temp = temp + m(20, 13) * x(13) 1680temp = temp + m(20, 14) * x(14) 1681temp = temp + m(20, 15) * x(15) 1682temp = temp + m(20, 16) * x(16) 1683temp = temp + m(20, 17) * x(17) 1684temp = temp + m(20, 18) * x(18) 1685temp = temp + m(20, 19) * x(19) 1686x(20) = ( x(20) - temp ) / m(20, 20) 1687 1688 1689! Backward substitution. 1690 1691temp = m(19, 20) * x(20) 1692x(19) = x(19) - temp 1693 1694temp = m(18, 19) * x(19) 1695temp = temp + m(18, 20) * x(20) 1696x(18) = x(18) - temp 1697 1698temp = m(17, 18) * x(18) 1699temp = temp + m(17, 19) * x(19) 1700temp = temp + m(17, 20) * x(20) 1701x(17) = x(17) - temp 1702 1703temp = m(16, 17) * x(17) 1704temp = temp + m(16, 19) * x(19) 1705temp = temp + m(16, 20) * x(20) 1706x(16) = x(16) - temp 1707 1708temp = m(15, 16) * x(16) 1709temp = temp + m(15, 19) * x(19) 1710temp = temp + m(15, 20) * x(20) 1711x(15) = x(15) - temp 1712 1713temp = m(14, 15) * x(15) 1714temp = temp + m(14, 19) * x(19) 1715temp = temp + m(14, 20) * x(20) 1716x(14) = x(14) - temp 1717 1718temp = m(13, 19) * x(19) 1719temp = temp + m(13, 20) * x(20) 1720x(13) = x(13) - temp 1721 1722temp = m(12, 19) * x(19) 1723temp = temp + m(12, 20) * x(20) 1724x(12) = x(12) - temp 1725 1726temp = m(11, 12) * x(12) 1727temp = temp + m(11, 19) * x(19) 1728x(11) = x(11) - temp 1729 1730temp = m(10, 12) * x(12) 1731temp = temp + m(10, 20) * x(20) 1732x(10) = x(10) - temp 1733 1734temp = m(9, 15) * x(15) 1735temp = temp + m(9, 18) * x(18) 1736temp = temp + m(9, 19) * x(19) 1737x(9) = x(9) - temp 1738 1739temp = m(8, 15) * x(15) 1740temp = temp + m(8, 20) * x(20) 1741x(8) = x(8) - temp 1742 1743temp = m(7, 18) * x(18) 1744temp = temp + m(7, 19) * x(19) 1745x(7) = x(7) - temp 1746 1747temp = m(6, 15) * x(15) 1748x(6) = x(6) - temp 1749 1750temp = m(5, 15) * x(15) 1751x(5) = x(5) - temp 1752 1753temp = m(4, 10) * x(10) 1754x(4) = x(4) - temp 1755 1756temp = m(3, 14) * x(14) 1757x(3) = x(3) - temp 1758 1759temp = m(2, 16) * x(16) 1760x(2) = x(2) - temp 1761 1762return 1763end subroutine lu_solve_2 1764 1765