!------------------------------------------------------------------------------- ! This file is part of Code_Saturne, a general-purpose CFD tool. ! ! Copyright (C) 1998-2021 EDF S.A. ! ! This program is free software; you can redistribute it and/or modify it under ! the terms of the GNU General Public License as published by the Free Software ! Foundation; either version 2 of the License, or (at your option) any later ! version. ! ! This program is distributed in the hope that it will be useful, but WITHOUT ! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS ! FOR A PARTICULAR PURPOSE. See the GNU General Public License for more ! details. ! ! You should have received a copy of the GNU General Public License along with ! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin ! Street, Fifth Floor, Boston, MA 02110-1301, USA. !------------------------------------------------------------------------------- !> \file chem_luscheme2.f90 !> \brief Routines for atmospheric chemical scheme 2 !> \remarks !> These routines are automatically generated by SPACK !> See CEREA: http://cerea.enpc.fr/polyphemus !> kinetic_2 !> \brief Computation of kinetic rates for atmospheric chemistry !------------------------------------------------------------------------------ !------------------------------------------------------------------------------ ! Arguments !------------------------------------------------------------------------------ ! mode name role !------------------------------------------------------------------------------ !> \param[in] nr total number of chemical reactions !> \param[in] option_photolysis flag to activate or not photolysis reactions !> \param[in] azi solar zenith angle !> \param[in] att atmospheric attenuation variable !> \param[in] temp temperature !> \param[in] press pressure !> \param[in] xlw water massic fraction !> \param[out] rk(nr) kinetic rates !______________________________________________________________________________ subroutine kinetic_2(nr,rk,temp,xlw,press,azi,att, & option_photolysis) implicit none ! Arguments integer nr double precision rk(nr),temp,xlw,press double precision azi, att integer option_photolysis ! Local variables double precision effko,rapk,summ double precision ylh2o ! Compute third body. ! Conversion = Avogadro*1d-6/Perfect gas constant. ! PRESS in Pascal, SUMM in molecules/cm3, TEMP in Kelvin summ = press * 7.243d16 / temp ! Number of water molecules computed from the massic fraction ! (absolute humidity) ylh2o = 29.d0*summ*xlw/(18.d0+11.d0*xlw) ! For the zenithal angle at tropics azi=abs(azi) rk( 1) = dexp(-0.8860689615829534d+02 & - ( -0.5300000000000000d+03 )/temp) rk( 1) = rk( 1) * summ * 0.2d0 rk( 2) = dexp(-0.2653240882726044d+02 & - ( 0.1500000000000000d+04 )/temp) if(option_photolysis.eq.2) then rk( 3)= 0.d0 elseif(option_photolysis.eq.1) then if(azi.lt.0.d0)then stop elseif(azi.ge. 0.00d+00 .and. azi.lt. 0.10d+02) then rk( 3)=-0.1302720567168795d-07 rk( 3)=-0.7822279432831311d-06+(azi- 0.00d+00) * rk( 3) rk( 3)= 0.0000000000000000d+00+(azi- 0.00d+00) * rk( 3) rk( 3)= 0.9310260000000001d-02+(azi- 0.00d+00) * rk( 3) elseif(azi.ge. 0.10d+02 .and. azi.lt. 0.20d+02) then rk( 3)= 0.3771617015067078d-08 rk( 3)=-0.1173044113433769d-05+(azi- 0.10d+02) * rk( 3) rk( 3)=-0.1955272056716901d-04+(azi- 0.10d+02) * rk( 3) rk( 3)= 0.9219010000000000d-02+(azi- 0.10d+02) * rk( 3) elseif(azi.ge. 0.20d+02 .and. azi.lt. 0.30d+02) then rk( 3)=-0.5859262388581815d-08 rk( 3)=-0.1059895602981758d-05+(azi- 0.20d+02) * rk( 3) rk( 3)=-0.4188211773132428d-04+(azi- 0.20d+02) * rk( 3) rk( 3)= 0.8909950000000000d-02+(azi- 0.20d+02) * rk( 3) elseif(azi.ge. 0.30d+02 .and. azi.lt. 0.40d+02) then rk( 3)=-0.7024567460738029d-08 rk( 3)=-0.1235673474639213d-05+(azi- 0.30d+02) * rk( 3) rk( 3)=-0.6483780850753392d-04+(azi- 0.30d+02) * rk( 3) rk( 3)= 0.8379279999999999d-02+(azi- 0.30d+02) * rk( 3) elseif(azi.ge. 0.40d+02 .and. azi.lt. 0.50d+02) then rk( 3)=-0.9202467768466835d-08 rk( 3)=-0.1446410498461367d-05+(azi- 0.40d+02) * rk( 3) rk( 3)=-0.9165864823853972d-04+(azi- 0.40d+02) * rk( 3) rk( 3)= 0.7600310000000000d-02+(azi- 0.40d+02) * rk( 3) elseif(azi.ge. 0.50d+02 .and. azi.lt. 0.60d+02) then rk( 3)=-0.1612556146540100d-07 rk( 3)=-0.1722484531515342d-05+(azi- 0.50d+02) * rk( 3) rk( 3)=-0.1233475985383066d-03+(azi- 0.50d+02) * rk( 3) rk( 3)= 0.6529880000000000d-02+(azi- 0.50d+02) * rk( 3) elseif(azi.ge. 0.60d+02 .and. azi.lt. 0.70d+02) then rk( 3)= 0.3226471363007382d-07 rk( 3)=-0.2206251375477548d-05+(azi- 0.60d+02) * rk( 3) rk( 3)=-0.1626349576082332d-03+(azi- 0.60d+02) * rk( 3) rk( 3)= 0.5108030000000000d-02+(azi- 0.60d+02) * rk( 3) elseif(azi.ge. 0.70d+02 .and. azi.lt. 0.78d+02) then rk( 3)= 0.2027078243961372d-06 rk( 3)=-0.1238309966574737d-05+(azi- 0.70d+02) * rk( 3) rk( 3)=-0.1970805710287543d-03+(azi- 0.70d+02) * rk( 3) rk( 3)= 0.3293320000000000d-02+(azi- 0.70d+02) * rk( 3) elseif(azi.ge. 0.78d+02 .and. azi.lt. 0.86d+02) then rk( 3)=-0.7448311471194499d-07 rk( 3)= 0.3626677818932555d-05+(azi- 0.78d+02) * rk( 3) rk( 3)=-0.1779736282099126d-03+(azi- 0.78d+02) * rk( 3) rk( 3)= 0.1741210000000000d-02+(azi- 0.78d+02) * rk( 3) elseif(azi.ge. 0.86d+02 .and. azi.lt. 0.90d+02) then rk( 3)= 0.2490309929270573d-05 rk( 3)= 0.1839083065842406d-05+(azi- 0.86d+02) * rk( 3) rk( 3)=-0.1342475411316713d-03+(azi- 0.86d+02) * rk( 3) rk( 3)= 0.5113930000000000d-03+(azi- 0.86d+02) * rk( 3) elseif(azi.ge.90.d0)then rk( 3)= 0.1632080000000000d-03 endif if(att.lt.0.99999) rk( 3) = rk( 3) * att endif rk( 4) = summ * 6.0d-34 * (temp/3.d2) ** (-2.4d0) rk( 4) = rk( 4) * summ * 0.2d0 if(option_photolysis.eq.2) then rk( 5)= 0.d0 elseif(option_photolysis.eq.1) then if(azi.lt.0.d0)then stop elseif(azi.ge. 0.00d+00 .and. azi.lt. 0.10d+02) then rk( 5)=-0.5928286648807996d-09 rk( 5)=-0.3096171335119280d-07+(azi- 0.00d+00) * rk( 5) rk( 5)= 0.0000000000000000d+00+(azi- 0.00d+00) * rk( 5) rk( 5)= 0.4927580000000000d-03+(azi- 0.00d+00) * rk( 5) elseif(azi.ge. 0.10d+02 .and. azi.lt. 0.20d+02) then rk( 5)= 0.1444859946426876d-09 rk( 5)=-0.4874657329761677d-07+(azi- 0.10d+02) * rk( 5) rk( 5)=-0.7970828664880956d-06+(azi- 0.10d+02) * rk( 5) rk( 5)= 0.4890690000000000d-03+(azi- 0.10d+02) * rk( 5) elseif(azi.ge. 0.20d+02 .and. azi.lt. 0.30d+02) then rk( 5)=-0.5511531369012520d-10 rk( 5)=-0.4441199345833616d-07+(azi- 0.20d+02) * rk( 5) rk( 5)=-0.1728668534047625d-05+(azi- 0.20d+02) * rk( 5) rk( 5)= 0.4763680000000000d-03+(azi- 0.20d+02) * rk( 5) elseif(azi.ge. 0.30d+02 .and. azi.lt. 0.40d+02) then rk( 5)=-0.3000247398821449d-09 rk( 5)=-0.4606545286904014d-07+(azi- 0.30d+02) * rk( 5) rk( 5)=-0.2633442997321385d-05+(azi- 0.30d+02) * rk( 5) rk( 5)= 0.4545850000000000d-03+(azi- 0.30d+02) * rk( 5) elseif(azi.ge. 0.40d+02 .and. azi.lt. 0.50d+02) then rk( 5)=-0.2397857267812366d-09 rk( 5)=-0.5506619506550444d-07+(azi- 0.40d+02) * rk( 5) rk( 5)=-0.3644759476666826d-05+(azi- 0.40d+02) * rk( 5) rk( 5)= 0.4233440000000000d-03+(azi- 0.40d+02) * rk( 5) elseif(azi.ge. 0.50d+02 .and. azi.lt. 0.60d+02) then rk( 5)=-0.1844832352993067d-08 rk( 5)=-0.6225976686893995d-07+(azi- 0.50d+02) * rk( 5) rk( 5)=-0.4818019096011278d-05+(azi- 0.50d+02) * rk( 5) rk( 5)= 0.3811500000000000d-03+(azi- 0.50d+02) * rk( 5) elseif(azi.ge. 0.60d+02 .and. azi.lt. 0.70d+02) then rk( 5)= 0.1101151387530619d-09 rk( 5)=-0.1176047374587370d-06+(azi- 0.60d+02) * rk( 5) rk( 5)=-0.6616664139287950d-05+(azi- 0.60d+02) * rk( 5) rk( 5)= 0.3248990000000000d-03+(azi- 0.60d+02) * rk( 5) elseif(azi.ge. 0.70d+02 .and. azi.lt. 0.78d+02) then rk( 5)=-0.1557211541866474d-07 rk( 5)=-0.1143012832961418d-06+(azi- 0.70d+02) * rk( 5) rk( 5)=-0.8935724346837023d-05+(azi- 0.70d+02) * rk( 5) rk( 5)= 0.2470820000000000d-03+(azi- 0.70d+02) * rk( 5) elseif(azi.ge. 0.78d+02 .and. azi.lt. 0.86d+02) then rk( 5)= 0.4048472604232427d-07 rk( 5)=-0.4880320533439059d-06+(azi- 0.78d+02) * rk( 5) rk( 5)=-0.1375439103995816d-04+(azi- 0.78d+02) * rk( 5) rk( 5)= 0.1603080000000000d-03+(azi- 0.78d+02) * rk( 5) elseif(azi.ge. 0.86d+02 .and. azi.lt. 0.90d+02) then rk( 5)= 0.2066880316654862d-06 rk( 5)= 0.4836013716715513d-06+(azi- 0.86d+02) * rk( 5) rk( 5)=-0.1378983649333310d-04+(azi- 0.86d+02) * rk( 5) rk( 5)= 0.3976700000000000d-04+(azi- 0.86d+02) * rk( 5) elseif(azi.ge.90.d0)then rk( 5)= 0.5573310000000000d-05 endif if(att.lt.0.99999) rk( 5) = rk( 5) * att endif if(option_photolysis.eq.2) then rk( 6)= 0.d0 elseif(option_photolysis.eq.1) then if(azi.lt.0.d0)then stop elseif(azi.ge. 0.00d+00 .and. azi.lt. 0.10d+02) then rk( 6)=-0.8776629099833464d-10 rk( 6)=-0.1165033709001661d-07+(azi- 0.00d+00) * rk( 6) rk( 6)= 0.0000000000000000d+00+(azi- 0.00d+00) * rk( 6) rk( 6)= 0.3523480000000000d-04+(azi- 0.00d+00) * rk( 6) elseif(azi.ge. 0.10d+02 .and. azi.lt. 0.20d+02) then rk( 6)= 0.1474988729949909d-09 rk( 6)=-0.1428332581996665d-07+(azi- 0.10d+02) * rk( 6) rk( 6)=-0.2593366290998327d-06+(azi- 0.10d+02) * rk( 6) rk( 6)= 0.3398200000000000d-04+(azi- 0.10d+02) * rk( 6) elseif(azi.ge. 0.20d+02 .and. azi.lt. 0.30d+02) then rk( 6)= 0.1300707990183827d-09 rk( 6)=-0.9858359630116941d-08+(azi- 0.20d+02) * rk( 6) rk( 6)=-0.5007534836006686d-06+(azi- 0.20d+02) * rk( 6) rk( 6)= 0.3010780000000000d-04+(azi- 0.20d+02) * rk( 6) elseif(azi.ge. 0.30d+02 .and. azi.lt. 0.40d+02) then rk( 6)= 0.1988179309314732d-09 rk( 6)=-0.5956235659565481d-08+(azi- 0.30d+02) * rk( 6) rk( 6)=-0.6588994364974921d-06+(azi- 0.30d+02) * rk( 6) rk( 6)= 0.2424450000000000d-04+(azi- 0.30d+02) * rk( 6) elseif(azi.ge. 0.40d+02 .and. azi.lt. 0.50d+02) then rk( 6)= 0.2219574772557277d-09 rk( 6)= 0.8302268378835231d-11+(azi- 0.40d+02) * rk( 6) rk( 6)=-0.7183787704093613d-06+(azi- 0.40d+02) * rk( 6) rk( 6)= 0.1725870000000000d-04+(azi- 0.40d+02) * rk( 6) elseif(azi.ge. 0.50d+02 .and. azi.lt. 0.60d+02) then rk( 6)= 0.1913521600455895d-09 rk( 6)= 0.6667026586050136d-08+(azi- 0.50d+02) * rk( 6) rk( 6)=-0.6516254818650674d-06+(azi- 0.50d+02) * rk( 6) rk( 6)= 0.1029770000000000d-04+(azi- 0.50d+02) * rk( 6) elseif(azi.ge. 0.60d+02 .and. azi.lt. 0.70d+02) then rk( 6)= 0.1602388256152816d-10 rk( 6)= 0.1240759138741867d-07+(azi- 0.60d+02) * rk( 6) rk( 6)=-0.4608793021303539d-06+(azi- 0.60d+02) * rk( 6) rk( 6)= 0.4639500000000000d-05+(azi- 0.60d+02) * rk( 6) elseif(azi.ge. 0.70d+02 .and. azi.lt. 0.78d+02) then rk( 6)=-0.3089359890783960d-09 rk( 6)= 0.1288830786428400d-07+(azi- 0.70d+02) * rk( 6) rk( 6)=-0.2079203096133002d-06+(azi- 0.70d+02) * rk( 6) rk( 6)= 0.1287490000000000d-05+(azi- 0.70d+02) * rk( 6) elseif(azi.ge. 0.78d+02 .and. azi.lt. 0.86d+02) then rk( 6)=-0.2034952628632162d-09 rk( 6)= 0.5473844126395715d-08+(azi- 0.78d+02) * rk( 6) rk( 6)=-0.6102309368797090d-07+(azi- 0.78d+02) * rk( 6) rk( 6)= 0.2908040000000000d-06+(azi- 0.78d+02) * rk( 6) elseif(azi.ge. 0.86d+02 .and. azi.lt. 0.90d+02) then rk( 6)= 0.1623544916128788d-09 rk( 6)= 0.5899578175158973d-09+(azi- 0.86d+02) * rk( 6) rk( 6)=-0.1251267813581064d-07+(azi- 0.86d+02) * rk( 6) rk( 6)= 0.4875570000000000d-07+(azi- 0.86d+02) * rk( 6) elseif(azi.ge.90.d0)then rk( 6)= 0.1853500000000000d-07 endif if(att.lt.0.99999) rk( 6) = rk( 6) * att endif rk( 7) = dexp(-0.2458649867820512d+02 & - ( -0.1020000000000000d+03 )/temp) rk( 7) = rk( 7) * summ rk( 8) = 0.2200000000000000d-09 rk( 8) = rk( 8) * ylh2o if(option_photolysis.eq.2) then rk( 9)= 0.d0 elseif(option_photolysis.eq.1) then if(azi.lt.0.d0)then stop elseif(azi.ge. 0.00d+00 .and. azi.lt. 0.10d+02) then rk( 9)=-0.2887225450832658d-08 rk( 9)=-0.1810277454916759d-06+(azi- 0.00d+00) * rk( 9) rk( 9)= 0.0000000000000000d+00+(azi- 0.00d+00) * rk( 9) rk( 9)= 0.2046710000000000d-02+(azi- 0.00d+00) * rk( 9) elseif(azi.ge. 0.10d+02 .and. azi.lt. 0.20d+02) then rk( 9)= 0.8216763524985595d-09 rk( 9)=-0.2676445090166556d-06+(azi- 0.10d+02) * rk( 9) rk( 9)=-0.4486722545083316d-05+(azi- 0.10d+02) * rk( 9) rk( 9)= 0.2025720000000000d-02+(azi- 0.10d+02) * rk( 9) elseif(azi.ge. 0.20d+02 .and. azi.lt. 0.30d+02) then rk( 9)=-0.1309479959161308d-08 rk( 9)=-0.2429942184416991d-06+(azi- 0.20d+02) * rk( 9) rk( 9)=-0.9593109819666860d-05+(azi- 0.20d+02) * rk( 9) rk( 9)= 0.1954910000000000d-02+(azi- 0.20d+02) * rk( 9) elseif(azi.ge. 0.30d+02 .and. azi.lt. 0.40d+02) then rk( 9)=-0.1523756515853649d-08 rk( 9)=-0.2822786172165394d-06+(azi- 0.30d+02) * rk( 9) rk( 9)=-0.1484583817624924d-04+(azi- 0.30d+02) * rk( 9) rk( 9)= 0.1833370000000000d-02+(azi- 0.30d+02) * rk( 9) elseif(azi.ge. 0.40d+02 .and. azi.lt. 0.50d+02) then rk( 9)=-0.1745493977424016d-08 rk( 9)=-0.3279913126921461d-06+(azi- 0.40d+02) * rk( 9) rk( 9)=-0.2094853747533609d-04+(azi- 0.40d+02) * rk( 9) rk( 9)= 0.1655160000000000d-02+(azi- 0.40d+02) * rk( 9) elseif(azi.ge. 0.50d+02 .and. azi.lt. 0.60d+02) then rk( 9)=-0.2754267574450560d-08 rk( 9)=-0.3803561320148624d-06+(azi- 0.50d+02) * rk( 9) rk( 9)=-0.2803201192240625d-04+(azi- 0.50d+02) * rk( 9) rk( 9)= 0.1411130000000000d-02+(azi- 0.50d+02) * rk( 9) elseif(azi.ge. 0.60d+02 .and. azi.lt. 0.70d+02) then rk( 9)= 0.1037656427523324d-07 rk( 9)=-0.4629841592484096d-06+(azi- 0.60d+02) * rk( 9) rk( 9)=-0.3646541483503878d-04+(azi- 0.60d+02) * rk( 9) rk( 9)= 0.1090020000000000d-02+(azi- 0.60d+02) * rk( 9) elseif(azi.ge. 0.70d+02 .and. azi.lt. 0.78d+02) then rk( 9)= 0.4335158727139456d-07 rk( 9)=-0.1516872309913989d-06+(azi- 0.70d+02) * rk( 9) rk( 9)=-0.4261212873743828d-04+(azi- 0.70d+02) * rk( 9) rk( 9)= 0.6894440000000000d-03+(azi- 0.70d+02) * rk( 9) elseif(azi.ge. 0.78d+02 .and. azi.lt. 0.86d+02) then rk( 9)=-0.2976581610911796d-07 rk( 9)= 0.8887508635220705d-06+(azi- 0.78d+02) * rk( 9) rk( 9)=-0.3671561967719464d-04+(azi- 0.78d+02) * rk( 9) rk( 9)= 0.3610350000000000d-03+(azi- 0.78d+02) * rk( 9) elseif(azi.ge. 0.86d+02 .and. azi.lt. 0.90d+02) then rk( 9)= 0.5586598403857848d-06 rk( 9)= 0.1743712769036732d-06+(azi- 0.86d+02) * rk( 9) rk( 9)=-0.2821064255377481d-04+(azi- 0.86d+02) * rk( 9) rk( 9)= 0.1089500000000000d-03+(azi- 0.86d+02) * rk( 9) elseif(azi.ge.90.d0)then rk( 9)= 0.3465160000000000d-04 endif if(att.lt.0.99999) rk( 9) = rk( 9) * att endif if(option_photolysis.eq.2) then rk( 10)= 0.d0 elseif(option_photolysis.eq.1) then if(azi.lt.0.d0)then stop elseif(azi.ge. 0.00d+00 .and. azi.lt. 0.10d+02) then rk( 10)=-0.1441479345432039d-10 rk( 10)=-0.1242452065456794d-08+(azi- 0.00d+00) * rk( 10) rk( 10)= 0.0000000000000000d+00+(azi- 0.00d+00) * rk( 10) rk( 10)= 0.8394580000000000d-05+(azi- 0.00d+00) * rk( 10) elseif(azi.ge. 0.10d+02 .and. azi.lt. 0.20d+02) then rk( 10)= 0.8244380362960478d-11 rk( 10)=-0.1674895869086405d-08+(azi- 0.10d+02) * rk( 10) rk( 10)=-0.2917347934543199d-07+(azi- 0.10d+02) * rk( 10) rk( 10)= 0.8255920000000000d-05+(azi- 0.10d+02) * rk( 10) elseif(azi.ge. 0.20d+02 .and. azi.lt. 0.30d+02) then rk( 10)= 0.1172720024787734d-12 rk( 10)=-0.1427564458197592d-08+(azi- 0.20d+02) * rk( 10) rk( 10)=-0.6019808261827194d-07+(azi- 0.20d+02) * rk( 10) rk( 10)= 0.7804940000000000d-05+(azi- 0.20d+02) * rk( 10) elseif(azi.ge. 0.30d+02 .and. azi.lt. 0.40d+02) then rk( 10)= 0.4506531627127201d-11 rk( 10)=-0.1424046298123240d-08+(azi- 0.30d+02) * rk( 10) rk( 10)=-0.8871419018148013d-07+(azi- 0.30d+02) * rk( 10) rk( 10)= 0.7060320000000000d-05+(azi- 0.30d+02) * rk( 10) elseif(azi.ge. 0.40d+02 .and. azi.lt. 0.50d+02) then rk( 10)= 0.1086660148900755d-10 rk( 10)=-0.1288850349309390d-08+(azi- 0.40d+02) * rk( 10) rk( 10)=-0.1158431566558070d-06+(azi- 0.40d+02) * rk( 10) rk( 10)= 0.6035280000000000d-05+(azi- 0.40d+02) * rk( 10) elseif(azi.ge. 0.50d+02 .and. azi.lt. 0.60d+02) then rk( 10)= 0.1803706241686108d-10 rk( 10)=-0.9628523046392104d-09+(azi- 0.50d+02) * rk( 10) rk( 10)=-0.1383601831952927d-06+(azi- 0.50d+02) * rk( 10) rk( 10)= 0.4758830000000000d-05+(azi- 0.50d+02) * rk( 10) elseif(azi.ge. 0.60d+02 .and. azi.lt. 0.70d+02) then rk( 10)= 0.7329514884355585d-10 rk( 10)=-0.4217404321336693d-09+(azi- 0.60d+02) * rk( 10) rk( 10)=-0.1522061105630206d-06+(azi- 0.60d+02) * rk( 10) rk( 10)= 0.3296980000000000d-05+(azi- 0.60d+02) * rk( 10) elseif(azi.ge. 0.70d+02 .and. azi.lt. 0.78d+02) then rk( 10)= 0.6172011386267615d-10 rk( 10)= 0.1777114033174383d-08+(azi- 0.70d+02) * rk( 10) rk( 10)=-0.1386523745526241d-06+(azi- 0.70d+02) * rk( 10) rk( 10)= 0.1806040000000000d-05+(azi- 0.70d+02) * rk( 10) elseif(azi.ge. 0.78d+02 .and. azi.lt. 0.86d+02) then rk( 10)=-0.8279704635596216d-10 rk( 10)= 0.3258396765877763d-08+(azi- 0.78d+02) * rk( 10) rk( 10)=-0.9836828816022045d-07+(azi- 0.78d+02) * rk( 10) rk( 10)= 0.8421570000000000d-06+(azi- 0.78d+02) * rk( 10) elseif(azi.ge. 0.86d+02 .and. azi.lt. 0.90d+02) then rk( 10)= 0.1082517324572734d-08 rk( 10)= 0.1271267653334671d-08+(azi- 0.86d+02) * rk( 10) rk( 10)=-0.6213097280649386d-07+(azi- 0.86d+02) * rk( 10) rk( 10)= 0.2213560000000000d-06+(azi- 0.86d+02) * rk( 10) elseif(azi.ge.90.d0)then rk( 10)= 0.6245350000000000d-07 endif if(att.lt.0.99999) rk( 10) = rk( 10) * att endif if(option_photolysis.eq.2) then rk( 11)= 0.d0 elseif(option_photolysis.eq.1) then if(azi.lt.0.d0)then stop elseif(azi.ge. 0.00d+00 .and. azi.lt. 0.10d+02) then rk( 11)=-0.5816837387219519d-10 rk( 11)=-0.5184316261278087d-08+(azi- 0.00d+00) * rk( 11) rk( 11)= 0.0000000000000000d+00+(azi- 0.00d+00) * rk( 11) rk( 11)= 0.3220560000000000d-04+(azi- 0.00d+00) * rk( 11) elseif(azi.ge. 0.10d+02 .and. azi.lt. 0.20d+02) then rk( 11)= 0.3450512161659949d-10 rk( 11)=-0.6929367477443941d-08+(azi- 0.10d+02) * rk( 11) rk( 11)=-0.1211368373872203d-06+(azi- 0.10d+02) * rk( 11) rk( 11)= 0.3162900000000000d-04+(azi- 0.10d+02) * rk( 11) elseif(azi.ge. 0.20d+02 .and. azi.lt. 0.30d+02) then rk( 11)= 0.3647887405787883d-11 rk( 11)=-0.5894213828945960d-08+(azi- 0.20d+02) * rk( 11) rk( 11)=-0.2493726504511192d-06+(azi- 0.20d+02) * rk( 11) rk( 11)= 0.2975920000000000d-04+(azi- 0.20d+02) * rk( 11) elseif(azi.ge. 0.30d+02 .and. azi.lt. 0.40d+02) then rk( 11)= 0.2530332876025029d-10 rk( 11)=-0.5784777206772343d-08+(azi- 0.30d+02) * rk( 11) rk( 11)=-0.3661625608083018d-06+(azi- 0.30d+02) * rk( 11) rk( 11)= 0.2667970000000000d-04+(azi- 0.30d+02) * rk( 11) elseif(azi.ge. 0.40d+02 .and. azi.lt. 0.50d+02) then rk( 11)= 0.6013879755321277d-10 rk( 11)=-0.5025677343964861d-08+(azi- 0.40d+02) * rk( 11) rk( 11)=-0.4742671063156733d-06+(azi- 0.40d+02) * rk( 11) rk( 11)= 0.2246490000000000d-04+(azi- 0.40d+02) * rk( 11) elseif(azi.ge. 0.50d+02 .and. azi.lt. 0.60d+02) then rk( 11)= 0.1004414810269593d-09 rk( 11)=-0.3221513417368319d-08+(azi- 0.50d+02) * rk( 11) rk( 11)=-0.5567390139290089d-06+(azi- 0.50d+02) * rk( 11) rk( 11)= 0.1727980000000000d-04+(azi- 0.50d+02) * rk( 11) elseif(azi.ge. 0.60d+02 .and. azi.lt. 0.70d+02) then rk( 11)= 0.3107352783389442d-09 rk( 11)=-0.2082689865616575d-09+(azi- 0.60d+02) * rk( 11) rk( 11)=-0.5910368379682989d-06+(azi- 0.60d+02) * rk( 11) rk( 11)= 0.1149070000000000d-04+(azi- 0.60d+02) * rk( 11) elseif(azi.ge. 0.70d+02 .and. azi.lt. 0.78d+02) then rk( 11)= 0.1345128013872392d-09 rk( 11)= 0.9113789363612173d-08+(azi- 0.70d+02) * rk( 11) rk( 11)=-0.5019816341977565d-06+(azi- 0.70d+02) * rk( 11) rk( 11)= 0.5870240000000000d-05+(azi- 0.70d+02) * rk( 11) elseif(azi.ge. 0.78d+02 .and. azi.lt. 0.86d+02) then rk( 11)=-0.1960445509611151d-09 rk( 11)= 0.1234209659691269d-07+(azi- 0.78d+02) * rk( 11) rk( 11)=-0.3303345465135304d-06+(azi- 0.78d+02) * rk( 11) rk( 11)= 0.2506540000000000d-05+(azi- 0.78d+02) * rk( 11) elseif(azi.ge. 0.86d+02 .and. azi.lt. 0.90d+02) then rk( 11)= 0.2279277828298341d-08 rk( 11)= 0.7637027373778166d-08+(azi- 0.86d+02) * rk( 11) rk( 11)=-0.1705015547476783d-06+(azi- 0.86d+02) * rk( 11) rk( 11)= 0.5533830000000000d-06+(azi- 0.86d+02) * rk( 11) elseif(azi.ge.90.d0)then rk( 11)= 0.1394430000000000d-06 endif if(att.lt.0.99999) rk( 11) = rk( 11) * att endif if(option_photolysis.eq.2) then rk( 12)= 0.d0 elseif(option_photolysis.eq.1) then if(azi.lt.0.d0)then stop elseif(azi.ge. 0.00d+00 .and. azi.lt. 0.10d+02) then rk( 12)=-0.7261045436751545d-10 rk( 12)=-0.5576895456324842d-08+(azi- 0.00d+00) * rk( 12) rk( 12)= 0.0000000000000000d+00+(azi- 0.00d+00) * rk( 12) rk( 12)= 0.4669460000000000d-04+(azi- 0.00d+00) * rk( 12) elseif(azi.ge. 0.10d+02 .and. azi.lt. 0.20d+02) then rk( 12)= 0.2853136310254372d-10 rk( 12)=-0.7755209087350302d-08+(azi- 0.10d+02) * rk( 12) rk( 12)=-0.1333210454367515d-06+(azi- 0.10d+02) * rk( 12) rk( 12)= 0.4606430000000000d-04+(azi- 0.10d+02) * rk( 12) elseif(azi.ge. 0.20d+02 .and. azi.lt. 0.30d+02) then rk( 12)=-0.1901499804265726d-10 rk( 12)=-0.6899268194273995d-08+(azi- 0.20d+02) * rk( 12) rk( 12)=-0.2798658182529944d-06+(azi- 0.20d+02) * rk( 12) rk( 12)= 0.4398410000000000d-04+(azi- 0.20d+02) * rk( 12) elseif(azi.ge. 0.30d+02 .and. azi.lt. 0.40d+02) then rk( 12)=-0.9471370931908783d-11 rk( 12)=-0.7469718135553710d-08+(azi- 0.30d+02) * rk( 12) rk( 12)=-0.4235556815512711d-06+(azi- 0.30d+02) * rk( 12) rk( 12)= 0.4047650000000000d-04+(azi- 0.30d+02) * rk( 12) elseif(azi.ge. 0.40d+02 .and. azi.lt. 0.50d+02) then rk( 12)= 0.1440048177028887d-10 rk( 12)=-0.7753859263510966d-08+(azi- 0.40d+02) * rk( 12) rk( 12)=-0.5757914555419174d-06+(azi- 0.40d+02) * rk( 12) rk( 12)= 0.3548450000000000d-04+(azi- 0.40d+02) * rk( 12) elseif(azi.ge. 0.50d+02 .and. azi.lt. 0.60d+02) then rk( 12)= 0.3856944385081257d-10 rk( 12)=-0.7321844810402538d-08+(azi- 0.50d+02) * rk( 12) rk( 12)=-0.7265484962810543d-06+(azi- 0.50d+02) * rk( 12) rk( 12)= 0.2896560000000000d-04+(azi- 0.50d+02) * rk( 12) elseif(azi.ge. 0.60d+02 .and. azi.lt. 0.70d+02) then rk( 12)= 0.4136217428262900d-09 rk( 12)=-0.6164761494878585d-08+(azi- 0.60d+02) * rk( 12) rk( 12)=-0.8614145593338596d-06+(azi- 0.60d+02) * rk( 12) rk( 12)= 0.2100650000000000d-04+(azi- 0.60d+02) * rk( 12) elseif(azi.ge. 0.70d+02 .and. azi.lt. 0.78d+02) then rk( 12)= 0.5376701572532502d-09 rk( 12)= 0.6243890789914774d-08+(azi- 0.70d+02) * rk( 12) rk( 12)=-0.8606232663835604d-06+(azi- 0.70d+02) * rk( 12) rk( 12)= 0.1218950000000000d-04+(azi- 0.70d+02) * rk( 12) elseif(azi.ge. 0.78d+02 .and. azi.lt. 0.86d+02) then rk( 12)=-0.5508273899935273d-09 rk( 12)= 0.1914797456399955d-07+(azi- 0.78d+02) * rk( 12) rk( 12)=-0.6574883435524898d-06+(azi- 0.78d+02) * rk( 12) rk( 12)= 0.5979410000000000d-05+(azi- 0.78d+02) * rk( 12) elseif(azi.ge. 0.86d+02 .and. azi.lt. 0.90d+02) then rk( 12)= 0.8530305661922505d-08 rk( 12)= 0.5928117204100688d-08+(azi- 0.86d+02) * rk( 12) rk( 12)=-0.4568796094072541d-06+(azi- 0.86d+02) * rk( 12) rk( 12)= 0.1662950000000000d-05+(azi- 0.86d+02) * rk( 12) elseif(azi.ge.90.d0)then rk( 12)= 0.4762210000000000d-06 endif if(att.lt.0.99999) rk( 12) = rk( 12) * att endif if(option_photolysis.eq.2) then rk( 13)= 0.d0 elseif(option_photolysis.eq.1) then if(azi.lt.0.d0)then stop elseif(azi.ge. 0.00d+00 .and. azi.lt. 0.10d+02) then rk( 13)=-0.1292568102250478d-10 rk( 13)=-0.1349143189774952d-08+(azi- 0.00d+00) * rk( 13) rk( 13)= 0.0000000000000000d+00+(azi- 0.00d+00) * rk( 13) rk( 13)= 0.6105070000000000d-05+(azi- 0.00d+00) * rk( 13) elseif(azi.ge. 0.10d+02 .and. azi.lt. 0.20d+02) then rk( 13)= 0.1221704306751377d-10 rk( 13)=-0.1736913620450095d-08+(azi- 0.10d+02) * rk( 13) rk( 13)=-0.3086056810225046d-07+(azi- 0.10d+02) * rk( 13) rk( 13)= 0.5957230000000000d-05+(azi- 0.10d+02) * rk( 13) elseif(azi.ge. 0.20d+02 .and. azi.lt. 0.30d+02) then rk( 13)= 0.7227508752451222d-11 rk( 13)=-0.1370402328424685d-08+(azi- 0.20d+02) * rk( 13) rk( 13)=-0.6193372759099823d-07+(azi- 0.20d+02) * rk( 13) rk( 13)= 0.5487150000000000d-05+(azi- 0.20d+02) * rk( 13) elseif(azi.ge. 0.30d+02 .and. azi.lt. 0.40d+02) then rk( 13)= 0.1521292192268060d-10 rk( 13)=-0.1153577065851153d-08+(azi- 0.30d+02) * rk( 13) rk( 13)=-0.8717352153375648d-07+(azi- 0.30d+02) * rk( 13) rk( 13)= 0.4738000000000000d-05+(azi- 0.30d+02) * rk( 13) elseif(azi.ge. 0.40d+02 .and. azi.lt. 0.50d+02) then rk( 13)= 0.2301080355683006d-10 rk( 13)=-0.6971894081707449d-09+(azi- 0.40d+02) * rk( 13) rk( 13)=-0.1056811862739754d-06+(azi- 0.40d+02) * rk( 13) rk( 13)= 0.3766120000000000d-05+(azi- 0.40d+02) * rk( 13) elseif(azi.ge. 0.50d+02 .and. azi.lt. 0.60d+02) then rk( 13)= 0.3106386384999131d-10 rk( 13)=-0.6865301465823343d-11+(azi- 0.50d+02) * rk( 13) rk( 13)=-0.1127217333703412d-06+(azi- 0.50d+02) * rk( 13) rk( 13)= 0.2662600000000000d-05+(azi- 0.50d+02) * rk( 13) elseif(azi.ge. 0.60d+02 .and. azi.lt. 0.70d+02) then rk( 13)= 0.4531874104319860d-10 rk( 13)= 0.9250506140339690d-09+(azi- 0.60d+02) * rk( 13) rk( 13)=-0.1035398802446585d-06+(azi- 0.60d+02) * rk( 13) rk( 13)= 0.1565760000000000d-05+(azi- 0.60d+02) * rk( 13) elseif(azi.ge. 0.70d+02 .and. azi.lt. 0.78d+02) then rk( 13)=-0.1303761111962380d-10 rk( 13)= 0.2284612845330880d-08+(azi- 0.70d+02) * rk( 13) rk( 13)=-0.7144324565100488d-07+(azi- 0.70d+02) * rk( 13) rk( 13)= 0.6681850000000000d-06+(azi- 0.70d+02) * rk( 13) elseif(azi.ge. 0.78d+02 .and. azi.lt. 0.86d+02) then rk( 13)=-0.4077238229824927d-10 rk( 13)= 0.1971710178462450d-08+(azi- 0.78d+02) * rk( 13) rk( 13)=-0.3739266146067857d-07+(azi- 0.78d+02) * rk( 13) rk( 13)= 0.2361790000000000d-06+(azi- 0.78d+02) * rk( 13) elseif(azi.ge. 0.86d+02 .and. azi.lt. 0.90d+02) then rk( 13)= 0.1193377495821847d-09 rk( 13)= 0.9931730033112436d-09+(azi- 0.86d+02) * rk( 13) rk( 13)=-0.1367359600643481d-07+(azi- 0.86d+02) * rk( 13) rk( 13)= 0.4235170000000000d-07+(azi- 0.86d+02) * rk( 13) elseif(azi.ge.90.d0)then rk( 13)= 0.1118570000000000d-07 endif if(att.lt.0.99999) rk( 13) = rk( 13) * att endif effko = 0.4900000000000000d-02* (temp / 3.d2) & **(- ( 0.0000000000000000d+00))* & dexp(- 0.1210000000000000d+05/temp) rapk = 0.5400000000000000d+17* (temp / 3.d2) & **(- ( 0.0000000000000000d+00))* & dexp(- 0.1383000000000000d+05/temp) rk( 14) = (effko*summ / ( 1.0d0 + effko * summ / & rapk)) * 0.3000d+00** (1.0d0 / (1.0d0 + & (dlog10(effko * summ / rapk))**2)) effko = 0.2700000000000000d-27* (temp / 3.d2) & **(- ( 0.7100000000000000d+01)) rapk = 0.1200000000000000d-10* (temp / 3.d2) & **(- ( 0.9000000000000000d+00)) rk( 15) = (effko * summ / & ( 1.0d0 + effko * summ / rapk)) * & 0.3000d+00** (1.0d0 / (1.0d0 + & (dlog10(effko * summ / rapk))**2)) rk( 16) = dexp(-0.2553915705425015d+02 & - ( -0.2700000000000000d+03 )/temp) rk( 17) = dexp(-0.2637825814743318d+02 & - ( -0.2500000000000000d+03 )/temp) effko = 0.7000000000000000d-30* (temp / 3.d2) & **(- ( 0.2600000000000000d+01)) rapk = 0.3600000000000000d-10* (temp / 3.d2) & **(- ( 0.1000000000000000d+00)) rk( 18) = (effko * summ / & ( 1.0d0 + effko * summ / rapk)) * & 0.6000d+00** (1.0d0 / (1.0d0 + & (dlog10(effko * summ / rapk))**2)) rk( 19) = dexp(-0.8322449114623726d+01 & + (-0.3000000000000000d+01 * dlog(temp)) ) effko = 0.3300000000000000d-30* (temp / 3.d2) & **(- ( 0.3300000000000000d+01)) rapk = 0.1500000000000000d-11* (temp / 3.d2) & **(- ( 0.0000000000000000d+00)) rk( 20) = (effko * summ / & ( 1.0d0 + effko * summ / rapk)) * & 0.6000d+00** (1.0d0 / (1.0d0 + & (dlog10(effko * summ / rapk))**2)) rk( 21) = dexp(-0.2975128465212864d+02 & - ( 0.2450000000000000d+04 )/temp) rk( 22) = dexp(-0.2492297091482634d+02 & - ( -0.1700000000000000d+03 )/temp) rk( 23) = dexp(-0.3073211390514037d+02 & - ( 0.1260000000000000d+04 )/temp) effko = 0.2000000000000000d-29* (temp / 3.d2) & **(- ( 0.4400000000000000d+01)) rapk = 0.1400000000000000d-11* (temp / 3.d2) & **(- ( 0.7000000000000000d+00)) rk( 24) = (effko * summ / & ( 1.0d0 + effko * summ / rapk)) * & 0.6000d+00** (1.0d0 / (1.0d0 + & (dlog10(effko * summ / rapk))**2)) effko = 0.1300000000000000d-02* (temp / 3.d2) & **(- ( 0.3500000000000000d+01))* & dexp(- 0.1100000000000000d+05/temp) rapk = 0.9700000000000000d+15* (temp / 3.d2) & **(- (-0.1000000000000000d+00))* & dexp(- 0.1108000000000000d+05/temp) rk( 25) = (effko*summ / ( 1.0d0 + effko * summ / & rapk)) * 0.4500d+00** (1.0d0 / (1.0d0 + & (dlog10(effko * summ / rapk))**2)) rk( 26) = 0.1000000000000000d-21 rk( 26) = rk( 26) * ylh2o rk( 27) = dexp(-0.2667550967090111d+02 & - ( -0.3650000000000000d+03 )/temp) rk( 28) = 0.6800000000000000d-13 rk( 29) = 2.3d-13 * dexp(600.0d0 / temp) & + 1.7d-33* summ * dexp(1000.0d0 / temp) rk( 30) = 3.22d-34 * dexp(2800.0d0 / temp) + & 2.38d-54 * summ * dexp(3200.0d0 / temp) rk( 30) = rk( 30) * ylh2o if(option_photolysis.eq.2) then rk( 31)= 0.d0 elseif(option_photolysis.eq.1) then if(azi.lt.0.d0)then stop elseif(azi.ge. 0.00d+00 .and. azi.lt. 0.10d+02) then rk( 31)=-0.2293082537254139d-06 rk( 31)=-0.1232691746274577d-04+(azi- 0.00d+00) * rk( 31) rk( 31)= 0.0000000000000000d+00+(azi- 0.00d+00) * rk( 31) rk( 31)= 0.2395610000000000d+00+(azi- 0.00d+00) * rk( 31) elseif(azi.ge. 0.10d+02 .and. azi.lt. 0.20d+02) then rk( 31)=-0.2107523882380334d-07 rk( 31)=-0.1920616507450818d-04+(azi- 0.10d+02) * rk( 31) rk( 31)=-0.3153308253725395d-03+(azi- 0.10d+02) * rk( 31) rk( 31)= 0.2380990000000000d+00+(azi- 0.10d+02) * rk( 31) elseif(azi.ge. 0.20d+02 .and. azi.lt. 0.30d+02) then rk( 31)= 0.2060920902067111d-07 rk( 31)=-0.1983842223922230d-04+(azi- 0.20d+02) * rk( 31) rk( 31)=-0.7057766985098441d-03+(azi- 0.20d+02) * rk( 31) rk( 31)= 0.2330040000000000d+00+(azi- 0.20d+02) * rk( 31) elseif(azi.ge. 0.30d+02 .and. azi.lt. 0.40d+02) then rk( 31)=-0.2323615972588770d-06 rk( 31)=-0.1922014596860219d-04+(azi- 0.30d+02) * rk( 31) rk( 31)=-0.1096362380588088d-02+(azi- 0.30d+02) * rk( 31) rk( 31)= 0.2239830000000000d+00+(azi- 0.30d+02) * rk( 31) elseif(azi.ge. 0.40d+02 .and. azi.lt. 0.50d+02) then rk( 31)=-0.2281628199851812d-06 rk( 31)=-0.2619099388636854d-04+(azi- 0.40d+02) * rk( 31) rk( 31)=-0.1550473779137796d-02+(azi- 0.40d+02) * rk( 31) rk( 31)= 0.2108650000000000d+00+(azi- 0.40d+02) * rk( 31) elseif(azi.ge. 0.50d+02 .and. azi.lt. 0.60d+02) then rk( 31)=-0.1187987122800317d-05 rk( 31)=-0.3303587848592447d-04+(azi- 0.50d+02) * rk( 31) rk( 31)=-0.2142742502860728d-02+(azi- 0.50d+02) * rk( 31) rk( 31)= 0.1925130000000000d+00+(azi- 0.50d+02) * rk( 31) elseif(azi.ge. 0.60d+02 .and. azi.lt. 0.70d+02) then rk( 31)=-0.3438886888129482d-06 rk( 31)=-0.6867549216993743d-04+(azi- 0.60d+02) * rk( 31) rk( 31)=-0.3159856209419318d-02+(azi- 0.60d+02) * rk( 31) rk( 31)= 0.1665940000000000d+00+(azi- 0.60d+02) * rk( 31) elseif(azi.ge. 0.70d+02 .and. azi.lt. 0.78d+02) then rk( 31)=-0.7936212779117297d-05 rk( 31)=-0.7899215283432154d-04+(azi- 0.70d+02) * rk( 31) rk( 31)=-0.4636532659461956d-02+(azi- 0.70d+02) * rk( 31) rk( 31)= 0.1277840000000000d+00+(azi- 0.70d+02) * rk( 31) elseif(azi.ge. 0.78d+02 .and. azi.lt. 0.86d+02) then rk( 31)= 0.2654257866666759d-04 rk( 31)=-0.2694612595331436d-03+(azi- 0.78d+02) * rk( 31) rk( 31)=-0.7424159958401733d-02+(azi- 0.78d+02) * rk( 31) rk( 31)= 0.8157290000000000d-01+(azi- 0.78d+02) * rk( 31) elseif(azi.ge. 0.86d+02 .and. azi.lt. 0.90d+02) then rk( 31)= 0.7705999956653109d-04 rk( 31)= 0.3675606284668786d-03+(azi- 0.86d+02) * rk( 31) rk( 31)=-0.6639365006930298d-02+(azi- 0.86d+02) * rk( 31) rk( 31)= 0.1852390000000000d-01+(azi- 0.86d+02) * rk( 31) elseif(azi.ge.90.d0)then rk( 31)= 0.2779250000000000d-02 endif if(att.lt.0.99999) rk( 31) = rk( 31) * att endif rk( 32) = dexp(-0.2590825451818744d+02 & - ( -0.1800000000000000d+03 )/temp) effko = 0.2500000000000000d-30* (temp / 3.d2) & **(- ( 0.1800000000000000d+01)) rapk = 0.2200000000000000d-10* (temp / 3.d2) & **(- ( 0.7000000000000000d+00)) rk( 33) = (effko * summ / & ( 1.0d0 + effko * summ / rapk)) * & 0.6000d+00** (1.0d0 / (1.0d0 + & (dlog10(effko * summ / rapk))**2)) if(option_photolysis.eq.2) then rk( 34)= 0.d0 elseif(option_photolysis.eq.1) then if(azi.lt.0.d0)then stop elseif(azi.ge. 0.00d+00 .and. azi.lt. 0.10d+02) then rk( 34)=-0.1726633761595217d-06 rk( 34)=-0.3763226238404795d-05+(azi- 0.00d+00) * rk( 34) rk( 34)= 0.0000000000000000d+00+(azi- 0.00d+00) * rk( 34) rk( 34)= 0.2758968400000000d-01+(azi- 0.00d+00) * rk( 34) elseif(azi.ge. 0.10d+02 .and. azi.lt. 0.20d+02) then rk( 34)= 0.8277361284785716d-06 rk( 34)=-0.8943127523190445d-05+(azi- 0.10d+02) * rk( 34) rk( 34)=-0.1270635376159524d-03+(azi- 0.10d+02) * rk( 34) rk( 34)= 0.2704069800000000d-01+(azi- 0.10d+02) * rk( 34) elseif(azi.ge. 0.20d+02 .and. azi.lt. 0.30d+02) then rk( 34)=-0.1363361137754772d-05 rk( 34)= 0.1588895633116670d-04+(azi- 0.20d+02) * rk( 34) rk( 34)=-0.5760524953618990d-04+(azi- 0.20d+02) * rk( 34) rk( 34)= 0.2570348600000000d-01+(azi- 0.20d+02) * rk( 34) elseif(azi.ge. 0.30d+02 .and. azi.lt. 0.40d+02) then rk( 34)= 0.9462274225405240d-06 rk( 34)=-0.2501187780147645d-04+(azi- 0.30d+02) * rk( 34) rk( 34)=-0.1488344642392872d-03+(azi- 0.30d+02) * rk( 34) rk( 34)= 0.2535296800000000d-01+(azi- 0.30d+02) * rk( 34) elseif(azi.ge. 0.40d+02 .and. azi.lt. 0.50d+02) then rk( 34)=-0.2183585524073524d-06 rk( 34)= 0.3374944874739267d-05+(azi- 0.40d+02) * rk( 34) rk( 34)=-0.3652037935066590d-03+(azi- 0.40d+02) * rk( 34) rk( 34)= 0.2230966300000000d-01+(azi- 0.40d+02) * rk( 34) elseif(azi.ge. 0.50d+02 .and. azi.lt. 0.60d+02) then rk( 34)= 0.1378647870888786d-06 rk( 34)=-0.3175811697481157d-05+(azi- 0.50d+02) * rk( 34) rk( 34)=-0.3632124617340762d-03+(azi- 0.50d+02) * rk( 34) rk( 34)= 0.1877676100000000d-01+(azi- 0.50d+02) * rk( 34) elseif(azi.ge. 0.60d+02 .and. azi.lt. 0.70d+02) then rk( 34)=-0.3865635959481289d-06 rk( 34)= 0.9601319151840079d-06+(azi- 0.60d+02) * rk( 34) rk( 34)=-0.3853692595570321d-03+(azi- 0.60d+02) * rk( 34) rk( 34)= 0.1496492000000000d-01+(azi- 0.60d+02) * rk( 34) elseif(azi.ge. 0.70d+02 .and. azi.lt. 0.78d+02) then rk( 34)=-0.1916508555648358d-06 rk( 34)=-0.1063677596325682d-04+(azi- 0.70d+02) * rk( 34) rk( 34)=-0.4821357000377863d-03+(azi- 0.70d+02) * rk( 34) rk( 34)= 0.1082067700000000d-01+(azi- 0.70d+02) * rk( 34) elseif(azi.ge. 0.78d+02 .and. azi.lt. 0.86d+02) then rk( 34)= 0.2540478756918318d-05 rk( 34)=-0.1523639649680941d-04+(azi- 0.78d+02) * rk( 34) rk( 34)=-0.6891210797183578d-03+(azi- 0.78d+02) * rk( 34) rk( 34)= 0.6184712500000000d-02+(azi- 0.78d+02) * rk( 34) elseif(azi.ge. 0.86d+02 .and. azi.lt. 0.90d+02) then rk( 34)= 0.1651057353835306d-05 rk( 34)= 0.4573509366928574d-04+(azi- 0.86d+02) * rk( 34) rk( 34)=-0.4451315023386027d-03+(azi- 0.86d+02) * rk( 34) rk( 34)= 0.9973396100000000d-03+(azi- 0.86d+02) * rk( 34) elseif(azi.ge.90.d0)then rk( 34)= 0.5424277000000000d-04 endif if(att.lt.0.99999) rk( 34) = rk( 34) * att endif return end subroutine kinetic_2 !=============================================================================== !> fexchem_2 !> \brief Computation of the chemical production terms !------------------------------------------------------------------------------ !------------------------------------------------------------------------------ ! Arguments !------------------------------------------------------------------------------ ! mode name role !------------------------------------------------------------------------------ !> \param[in] ns total number of chemical species !> \param[in] nr total number of chemical reactions !> \param[in] y concentrations vector !> \param[in] rk kinetic rates !> \param[in] zcsourc source term !> \param[in] convers_factor conversion factors !> \param[out] chem chemical production terms for every species !______________________________________________________________________________ subroutine fexchem_2(ns,nr,y,rk,zcsourc,convers_factor,chem) implicit none ! Arguments integer nr,ns double precision rk(nr),y(ns),chem(ns),zcsourc(ns) double precision convers_factor(ns) ! Local variables integer i double precision w(nr) double precision conc(ns) do i=1,ns chem(i)=0.d0 enddo ! Conversion mg/kg to molecules/cm3. do i = 1, ns conc(i) = y(i) * convers_factor(i) enddo ! Compute reaction rates. call rates_2(ns,nr,rk,conc,w) ! Chemical production terms. chem( 19) = chem( 19) + 0.2000000000000000d+01 * w( 1) chem( 20) = chem( 20) - 0.2000000000000000d+01 * w( 1) chem( 16) = chem( 16) - w( 2) chem( 19) = chem( 19) + w( 2) chem( 20) = chem( 20) - w( 2) chem( 17) = chem( 17) + w( 3) chem( 19) = chem( 19) - w( 3) chem( 20) = chem( 20) + w( 3) chem( 16) = chem( 16) + w( 4) chem( 17) = chem( 17) - w( 4) chem( 16) = chem( 16) - w( 5) chem( 17) = chem( 17) + w( 5) chem( 2) = chem( 2) + w( 6) chem( 16) = chem( 16) - w( 6) chem( 2) = chem( 2) - w( 7) chem( 17) = chem( 17) + w( 7) chem( 2) = chem( 2) - w( 8) chem( 15) = chem( 15) + 0.2000000000000000d+01 * w( 8) chem( 8) = chem( 8) - w( 9) chem( 15) = chem( 15) + w( 9) chem( 20) = chem( 20) + w( 9) chem( 3) = chem( 3) - w( 10) chem( 15) = chem( 15) + 0.2000000000000000d+01 * w( 10) chem( 4) = chem( 4) + w( 11) chem( 10) = chem( 10) - w( 11) chem( 14) = chem( 14) + 0.2000000000000000d+01 * w( 11) chem( 4) = chem( 4) + w( 12) chem( 10) = chem( 10) - w( 12) chem( 1) = chem( 1) - w( 13) chem( 4) = chem( 4) + w( 13) chem( 10) = chem( 10) + w( 13) chem( 13) = chem( 13) + w( 13) chem( 14) = chem( 14) + 0.2000000000000000d+01 * w( 13) chem( 11) = chem( 11) - w( 14) chem( 12) = chem( 12) + w( 14) chem( 19) = chem( 19) + w( 14) chem( 11) = chem( 11) + w( 15) chem( 12) = chem( 12) - w( 15) chem( 19) = chem( 19) - w( 15) chem( 10) = chem( 10) + w( 16) chem( 12) = chem( 12) - w( 16) chem( 13) = chem( 13) + w( 16) chem( 14) = chem( 14) + w( 16) chem( 19) = chem( 19) + w( 16) chem( 20) = chem( 20) - w( 16) chem( 14) = chem( 14) - w( 17) chem( 15) = chem( 15) + w( 17) chem( 19) = chem( 19) + w( 17) chem( 20) = chem( 20) - w( 17) chem( 8) = chem( 8) + w( 18) chem( 15) = chem( 15) - w( 18) chem( 20) = chem( 20) - w( 18) chem( 9) = chem( 9) + w( 19) chem( 15) = chem( 15) - w( 19) chem( 19) = chem( 19) - w( 19) chem( 5) = chem( 5) - w( 20) chem( 6) = chem( 6) + w( 20) chem( 14) = chem( 14) + w( 20) chem( 15) = chem( 15) - w( 20) chem( 16) = chem( 16) - w( 21) chem( 18) = chem( 18) + w( 21) chem( 19) = chem( 19) - w( 21) chem( 18) = chem( 18) - w( 22) chem( 19) = chem( 19) + 0.2000000000000000d+01 * w( 22) chem( 20) = chem( 20) - w( 22) chem( 18) = chem( 18) - w( 23) chem( 20) = chem( 20) + w( 23) chem( 7) = chem( 7) + w( 24) chem( 18) = chem( 18) - w( 24) chem( 19) = chem( 19) - w( 24) chem( 7) = chem( 7) - w( 25) chem( 18) = chem( 18) + w( 25) chem( 19) = chem( 19) + w( 25) chem( 7) = chem( 7) - w( 26) chem( 9) = chem( 9) + 0.2000000000000000d+01 * w( 26) chem( 13) = chem( 13) - w( 27) chem( 19) = chem( 19) + w( 27) chem( 20) = chem( 20) - w( 27) chem( 13) = chem( 13) - 0.2000000000000000d+01 * w( 28) chem( 3) = chem( 3) + w( 29) chem( 14) = chem( 14) - 0.2000000000000000d+01 * w( 29) chem( 3) = chem( 3) + w( 30) chem( 14) = chem( 14) - 0.2000000000000000d+01 * w( 30) chem( 17) = chem( 17) + 0.8900000000000000d+00 * w( 31) chem( 18) = chem( 18) - w( 31) chem( 19) = chem( 19) + 0.8900000000000000d+00 * w( 31) chem( 20) = chem( 20) + 0.1100000000000000d+00 * w( 31) chem( 17) = chem( 17) - w( 32) chem( 19) = chem( 19) - w( 32) chem( 20) = chem( 20) + w( 32) chem( 17) = chem( 17) - w( 33) chem( 18) = chem( 18) + w( 33) chem( 19) = chem( 19) - w( 33) chem( 7) = chem( 7) - w( 34) chem( 18) = chem( 18) + w( 34) chem( 19) = chem( 19) + w( 34) ! Conversion molecules/cm3/sec to mg/kg/sec. do i = 1, ns chem(i) = chem(i) / convers_factor(i) enddo ! Volumic source terms. do i=1,ns chem(i)=chem(i)+zcsourc(i) enddo return end subroutine fexchem_2 !=============================================================================== !> jacdchemdc_2 !> \brief Computation of the Jacobian matrix for atmospheric chemistry !------------------------------------------------------------------------------ !------------------------------------------------------------------------------ ! Arguments !------------------------------------------------------------------------------ ! mode name role !------------------------------------------------------------------------------ !> \param[in] ns total number of chemical species !> \param[in] nr total number of chemical reactions !> \param[in] y concentrations vector !> \param[in] convers_factor conversion factors of mug/m3 to !> molecules/cm3 !> \param[in] convers_factor_jac conversion factors for the Jacobian matrix !> (Wmol(i)/Wmol(j)) !> \param[in] rk kinetic rates !> \param[out] jacc Jacobian matrix !______________________________________________________________________________ subroutine jacdchemdc_2(ns,nr,y,convers_factor, & convers_factor_jac,rk,jacc) implicit none ! Arguments integer nr,ns double precision rk(nr),y(ns),jacc(ns,ns) double precision convers_factor(ns) double precision convers_factor_jac(ns,ns) ! Local variables integer i,j double precision dw(nr,ns) double precision conc(ns) do j=1,ns do i=1,ns jacc(i,j)=0.d0 enddo enddo ! Conversion mg/kg to molecules/cm3. do i = 1, ns conc(i) = y(i) * convers_factor(i) enddo call dratedc_2(ns,nr,rk,conc,dw) jacc( 19, 20) = jacc( 19, 20)+ 0.2000000000000000d+01*dw( 1, 20) jacc( 19, 20) = jacc( 19, 20)+ 0.2000000000000000d+01*dw( 1, 20) jacc( 20, 20) = jacc( 20, 20)- 0.2000000000000000d+01*dw( 1, 20) jacc( 20, 20) = jacc( 20, 20)- 0.2000000000000000d+01*dw( 1, 20) jacc( 16, 16) = jacc( 16, 16) - dw( 2, 16) jacc( 16, 20) = jacc( 16, 20) - dw( 2, 20) jacc( 19, 16) = jacc( 19, 16) + dw( 2, 16) jacc( 19, 20) = jacc( 19, 20) + dw( 2, 20) jacc( 20, 16) = jacc( 20, 16) - dw( 2, 16) jacc( 20, 20) = jacc( 20, 20) - dw( 2, 20) jacc( 17, 19) = jacc( 17, 19) + dw( 3, 19) jacc( 19, 19) = jacc( 19, 19) - dw( 3, 19) jacc( 20, 19) = jacc( 20, 19) + dw( 3, 19) jacc( 16, 17) = jacc( 16, 17) + dw( 4, 17) jacc( 17, 17) = jacc( 17, 17) - dw( 4, 17) jacc( 16, 16) = jacc( 16, 16) - dw( 5, 16) jacc( 17, 16) = jacc( 17, 16) + dw( 5, 16) jacc( 2, 16) = jacc( 2, 16) + dw( 6, 16) jacc( 16, 16) = jacc( 16, 16) - dw( 6, 16) jacc( 2, 2) = jacc( 2, 2) - dw( 7, 2) jacc( 17, 2) = jacc( 17, 2) + dw( 7, 2) jacc( 2, 2) = jacc( 2, 2) - dw( 8, 2) jacc( 15, 2) = jacc( 15, 2)+ 0.2000000000000000d+01*dw( 8, 2) jacc( 8, 8) = jacc( 8, 8) - dw( 9, 8) jacc( 15, 8) = jacc( 15, 8) + dw( 9, 8) jacc( 20, 8) = jacc( 20, 8) + dw( 9, 8) jacc( 3, 3) = jacc( 3, 3) - dw( 10, 3) jacc( 15, 3) = jacc( 15, 3)+ 0.2000000000000000d+01*dw( 10, 3) jacc( 4, 10) = jacc( 4, 10) + dw( 11, 10) jacc( 10, 10) = jacc( 10, 10) - dw( 11, 10) jacc( 14, 10) = jacc( 14, 10)+ 0.2000000000000000d+01*dw( 11, 10) jacc( 4, 10) = jacc( 4, 10) + dw( 12, 10) jacc( 10, 10) = jacc( 10, 10) - dw( 12, 10) jacc( 1, 1) = jacc( 1, 1) - dw( 13, 1) jacc( 4, 1) = jacc( 4, 1) + dw( 13, 1) jacc( 10, 1) = jacc( 10, 1) + dw( 13, 1) jacc( 13, 1) = jacc( 13, 1) + dw( 13, 1) jacc( 14, 1) = jacc( 14, 1)+ 0.2000000000000000d+01*dw( 13, 1) jacc( 11, 11) = jacc( 11, 11) - dw( 14, 11) jacc( 12, 11) = jacc( 12, 11) + dw( 14, 11) jacc( 19, 11) = jacc( 19, 11) + dw( 14, 11) jacc( 11, 12) = jacc( 11, 12) + dw( 15, 12) jacc( 11, 19) = jacc( 11, 19) + dw( 15, 19) jacc( 12, 12) = jacc( 12, 12) - dw( 15, 12) jacc( 12, 19) = jacc( 12, 19) - dw( 15, 19) jacc( 19, 12) = jacc( 19, 12) - dw( 15, 12) jacc( 19, 19) = jacc( 19, 19) - dw( 15, 19) jacc( 10, 12) = jacc( 10, 12) + dw( 16, 12) jacc( 10, 20) = jacc( 10, 20) + dw( 16, 20) jacc( 12, 12) = jacc( 12, 12) - dw( 16, 12) jacc( 12, 20) = jacc( 12, 20) - dw( 16, 20) jacc( 13, 12) = jacc( 13, 12) + dw( 16, 12) jacc( 13, 20) = jacc( 13, 20) + dw( 16, 20) jacc( 14, 12) = jacc( 14, 12) + dw( 16, 12) jacc( 14, 20) = jacc( 14, 20) + dw( 16, 20) jacc( 19, 12) = jacc( 19, 12) + dw( 16, 12) jacc( 19, 20) = jacc( 19, 20) + dw( 16, 20) jacc( 20, 12) = jacc( 20, 12) - dw( 16, 12) jacc( 20, 20) = jacc( 20, 20) - dw( 16, 20) jacc( 14, 14) = jacc( 14, 14) - dw( 17, 14) jacc( 14, 20) = jacc( 14, 20) - dw( 17, 20) jacc( 15, 14) = jacc( 15, 14) + dw( 17, 14) jacc( 15, 20) = jacc( 15, 20) + dw( 17, 20) jacc( 19, 14) = jacc( 19, 14) + dw( 17, 14) jacc( 19, 20) = jacc( 19, 20) + dw( 17, 20) jacc( 20, 14) = jacc( 20, 14) - dw( 17, 14) jacc( 20, 20) = jacc( 20, 20) - dw( 17, 20) jacc( 8, 15) = jacc( 8, 15) + dw( 18, 15) jacc( 8, 20) = jacc( 8, 20) + dw( 18, 20) jacc( 15, 15) = jacc( 15, 15) - dw( 18, 15) jacc( 15, 20) = jacc( 15, 20) - dw( 18, 20) jacc( 20, 15) = jacc( 20, 15) - dw( 18, 15) jacc( 20, 20) = jacc( 20, 20) - dw( 18, 20) jacc( 9, 19) = jacc( 9, 19) + dw( 19, 19) jacc( 9, 15) = jacc( 9, 15) + dw( 19, 15) jacc( 15, 19) = jacc( 15, 19) - dw( 19, 19) jacc( 15, 15) = jacc( 15, 15) - dw( 19, 15) jacc( 19, 19) = jacc( 19, 19) - dw( 19, 19) jacc( 19, 15) = jacc( 19, 15) - dw( 19, 15) jacc( 5, 5) = jacc( 5, 5) - dw( 20, 5) jacc( 5, 15) = jacc( 5, 15) - dw( 20, 15) jacc( 6, 5) = jacc( 6, 5) + dw( 20, 5) jacc( 6, 15) = jacc( 6, 15) + dw( 20, 15) jacc( 14, 5) = jacc( 14, 5) + dw( 20, 5) jacc( 14, 15) = jacc( 14, 15) + dw( 20, 15) jacc( 15, 5) = jacc( 15, 5) - dw( 20, 5) jacc( 15, 15) = jacc( 15, 15) - dw( 20, 15) jacc( 16, 19) = jacc( 16, 19) - dw( 21, 19) jacc( 16, 16) = jacc( 16, 16) - dw( 21, 16) jacc( 18, 19) = jacc( 18, 19) + dw( 21, 19) jacc( 18, 16) = jacc( 18, 16) + dw( 21, 16) jacc( 19, 19) = jacc( 19, 19) - dw( 21, 19) jacc( 19, 16) = jacc( 19, 16) - dw( 21, 16) jacc( 18, 18) = jacc( 18, 18) - dw( 22, 18) jacc( 18, 20) = jacc( 18, 20) - dw( 22, 20) jacc( 19, 18) = jacc( 19, 18)+ 0.2000000000000000d+01*dw( 22, 18) jacc( 19, 20) = jacc( 19, 20)+ 0.2000000000000000d+01*dw( 22, 20) jacc( 20, 18) = jacc( 20, 18) - dw( 22, 18) jacc( 20, 20) = jacc( 20, 20) - dw( 22, 20) jacc( 18, 18) = jacc( 18, 18) - dw( 23, 18) jacc( 18, 19) = jacc( 18, 19) - dw( 23, 19) jacc( 20, 18) = jacc( 20, 18) + dw( 23, 18) jacc( 20, 19) = jacc( 20, 19) + dw( 23, 19) jacc( 7, 18) = jacc( 7, 18) + dw( 24, 18) jacc( 7, 19) = jacc( 7, 19) + dw( 24, 19) jacc( 18, 18) = jacc( 18, 18) - dw( 24, 18) jacc( 18, 19) = jacc( 18, 19) - dw( 24, 19) jacc( 19, 18) = jacc( 19, 18) - dw( 24, 18) jacc( 19, 19) = jacc( 19, 19) - dw( 24, 19) jacc( 7, 7) = jacc( 7, 7) - dw( 25, 7) jacc( 18, 7) = jacc( 18, 7) + dw( 25, 7) jacc( 19, 7) = jacc( 19, 7) + dw( 25, 7) jacc( 7, 7) = jacc( 7, 7) - dw( 26, 7) jacc( 9, 7) = jacc( 9, 7)+ 0.2000000000000000d+01*dw( 26, 7) jacc( 13, 13) = jacc( 13, 13) - dw( 27, 13) jacc( 13, 20) = jacc( 13, 20) - dw( 27, 20) jacc( 19, 13) = jacc( 19, 13) + dw( 27, 13) jacc( 19, 20) = jacc( 19, 20) + dw( 27, 20) jacc( 20, 13) = jacc( 20, 13) - dw( 27, 13) jacc( 20, 20) = jacc( 20, 20) - dw( 27, 20) jacc( 13, 13) = jacc( 13, 13)- 0.2000000000000000d+01*dw( 28, 13) jacc( 13, 13) = jacc( 13, 13)- 0.2000000000000000d+01*dw( 28, 13) jacc( 3, 14) = jacc( 3, 14) + dw( 29, 14) jacc( 3, 14) = jacc( 3, 14) + dw( 29, 14) jacc( 14, 14) = jacc( 14, 14)- 0.2000000000000000d+01*dw( 29, 14) jacc( 14, 14) = jacc( 14, 14)- 0.2000000000000000d+01*dw( 29, 14) jacc( 3, 14) = jacc( 3, 14) + dw( 30, 14) jacc( 3, 14) = jacc( 3, 14) + dw( 30, 14) jacc( 14, 14) = jacc( 14, 14)- 0.2000000000000000d+01*dw( 30, 14) jacc( 14, 14) = jacc( 14, 14)- 0.2000000000000000d+01*dw( 30, 14) jacc( 17, 18) = jacc( 17, 18)+ 0.8900000000000000d+00*dw( 31, 18) jacc( 18, 18) = jacc( 18, 18) - dw( 31, 18) jacc( 19, 18) = jacc( 19, 18)+ 0.8900000000000000d+00*dw( 31, 18) jacc( 20, 18) = jacc( 20, 18)+ 0.1100000000000000d+00*dw( 31, 18) jacc( 17, 17) = jacc( 17, 17) - dw( 32, 17) jacc( 17, 19) = jacc( 17, 19) - dw( 32, 19) jacc( 19, 17) = jacc( 19, 17) - dw( 32, 17) jacc( 19, 19) = jacc( 19, 19) - dw( 32, 19) jacc( 20, 17) = jacc( 20, 17) + dw( 32, 17) jacc( 20, 19) = jacc( 20, 19) + dw( 32, 19) jacc( 17, 17) = jacc( 17, 17) - dw( 33, 17) jacc( 17, 19) = jacc( 17, 19) - dw( 33, 19) jacc( 18, 17) = jacc( 18, 17) + dw( 33, 17) jacc( 18, 19) = jacc( 18, 19) + dw( 33, 19) jacc( 19, 17) = jacc( 19, 17) - dw( 33, 17) jacc( 19, 19) = jacc( 19, 19) - dw( 33, 19) jacc( 7, 7) = jacc( 7, 7) - dw( 34, 7) jacc( 18, 7) = jacc( 18, 7) + dw( 34, 7) jacc( 19, 7) = jacc( 19, 7) + dw( 34, 7) do j = 1, ns do i = 1, ns jacc(i,j) = jacc(i,j) * convers_factor_jac(i,j) enddo enddo return end subroutine jacdchemdc_2 !=============================================================================== !> rates_2 !> \brief Computation of reaction rates !------------------------------------------------------------------------------ !------------------------------------------------------------------------------ ! Arguments !------------------------------------------------------------------------------ ! mode name role !------------------------------------------------------------------------------ !> \param[in] ns total number of chemical species !> \param[in] nr total number of chemical reactions !> \param[in] rk kinetic rates !> \param[in] y concentrations vector !> \param[out] w reaction rates !______________________________________________________________________________ subroutine rates_2(ns,nr,rk,y,w) implicit none ! Arguments integer nr,ns double precision rk(nr),y(ns) ! Local variables double precision w(nr) w( 1) = rk( 1) * y( 20) * y( 20) w( 2) = rk( 2) * y( 16) * y( 20) w( 3) = rk( 3) * y( 19) w( 4) = rk( 4) * y( 17) w( 5) = rk( 5) * y( 16) w( 6) = rk( 6) * y( 16) w( 7) = rk( 7) * y( 2) w( 8) = rk( 8) * y( 2) w( 9) = rk( 9) * y( 8) w( 10) = rk( 10) * y( 3) w( 11) = rk( 11) * y( 10) w( 12) = rk( 12) * y( 10) w( 13) = rk( 13) * y( 1) w( 14) = rk( 14) * y( 11) w( 15) = rk( 15) * y( 12) * y( 19) w( 16) = rk( 16) * y( 12) * y( 20) w( 17) = rk( 17) * y( 14) * y( 20) w( 18) = rk( 18) * y( 15) * y( 20) w( 19) = rk( 19) * y( 19) * y( 15) w( 20) = rk( 20) * y( 5) * y( 15) w( 21) = rk( 21) * y( 19) * y( 16) w( 22) = rk( 22) * y( 18) * y( 20) w( 23) = rk( 23) * y( 18) * y( 19) w( 24) = rk( 24) * y( 18) * y( 19) w( 25) = rk( 25) * y( 7) w( 26) = rk( 26) * y( 7) w( 27) = rk( 27) * y( 13) * y( 20) w( 28) = rk( 28) * y( 13) * y( 13) w( 29) = rk( 29) * y( 14) * y( 14) w( 30) = rk( 30) * y( 14) * y( 14) w( 31) = rk( 31) * y( 18) w( 32) = rk( 32) * y( 17) * y( 19) w( 33) = rk( 33) * y( 17) * y( 19) w( 34) = rk( 34) * y( 7) return end subroutine rates_2 !=============================================================================== !> dratedc_2 !> \brief Computation of derivatives of reaction rates !------------------------------------------------------------------------------ !------------------------------------------------------------------------------ ! Arguments !------------------------------------------------------------------------------ ! mode name role !------------------------------------------------------------------------------ !> \param[in] nr total number of chemical reactions !> \param[in] ns total number of chemical species !> \param[in] rk kinetic rates !> \param[in] y concentrations vector !> \param[out] dw derivatives of reaction rates !______________________________________________________________________________ subroutine dratedc_2(ns,nr,rk,y,dw) implicit none ! Arguments integer nr,ns double precision rk(nr),y(ns) ! Local variables double precision dw(nr,ns) dw( 1, 20) = rk( 1) * y( 20) dw( 1, 20) = rk( 1) * y( 20) dw( 2, 16) = rk( 2) * y( 20) dw( 2, 20) = rk( 2) * y( 16) dw( 3, 19) = rk( 3) dw( 4, 17) = rk( 4) dw( 5, 16) = rk( 5) dw( 6, 16) = rk( 6) dw( 7, 2) = rk( 7) dw( 8, 2) = rk( 8) dw( 9, 8) = rk( 9) dw( 10, 3) = rk( 10) dw( 11, 10) = rk( 11) dw( 12, 10) = rk( 12) dw( 13, 1) = rk( 13) dw( 14, 11) = rk( 14) dw( 15, 12) = rk( 15) * y( 19) dw( 15, 19) = rk( 15) * y( 12) dw( 16, 12) = rk( 16) * y( 20) dw( 16, 20) = rk( 16) * y( 12) dw( 17, 14) = rk( 17) * y( 20) dw( 17, 20) = rk( 17) * y( 14) dw( 18, 15) = rk( 18) * y( 20) dw( 18, 20) = rk( 18) * y( 15) dw( 19, 19) = rk( 19) * y( 15) dw( 19, 15) = rk( 19) * y( 19) dw( 20, 5) = rk( 20) * y( 15) dw( 20, 15) = rk( 20) * y( 5) dw( 21, 19) = rk( 21) * y( 16) dw( 21, 16) = rk( 21) * y( 19) dw( 22, 18) = rk( 22) * y( 20) dw( 22, 20) = rk( 22) * y( 18) dw( 23, 18) = rk( 23) * y( 19) dw( 23, 19) = rk( 23) * y( 18) dw( 24, 18) = rk( 24) * y( 19) dw( 24, 19) = rk( 24) * y( 18) dw( 25, 7) = rk( 25) dw( 26, 7) = rk( 26) dw( 27, 13) = rk( 27) * y( 20) dw( 27, 20) = rk( 27) * y( 13) dw( 28, 13) = rk( 28) * y( 13) dw( 28, 13) = rk( 28) * y( 13) dw( 29, 14) = rk( 29) * y( 14) dw( 29, 14) = rk( 29) * y( 14) dw( 30, 14) = rk( 30) * y( 14) dw( 30, 14) = rk( 30) * y( 14) dw( 31, 18) = rk( 31) dw( 32, 17) = rk( 32) * y( 19) dw( 32, 19) = rk( 32) * y( 17) dw( 33, 17) = rk( 33) * y( 19) dw( 33, 19) = rk( 33) * y( 17) dw( 34, 7) = rk( 34) return end subroutine dratedc_2 !=============================================================================== !> lu_decompose_2 !> \brief Computation of LU factorization of matrix m !------------------------------------------------------------------------------ !------------------------------------------------------------------------------ ! Arguments !------------------------------------------------------------------------------ ! mode name role !------------------------------------------------------------------------------ !> \param[in] ns matrix row number from the chemical species !> number !> \param[in,out] m on entry, an invertible matrix. !> On exit, an LU factorization of m !______________________________________________________________________________ subroutine lu_decompose_2 (ns,m) implicit none ! Arguments integer ns double precision m(ns,ns) ! Local variables double precision temp ! Upper part. m(2, 16) = m(2, 16) / m(2, 2) ! Upper part. m(3, 14) = m(3, 14) / m(3, 3) ! Upper part. m(4, 10) = m(4, 10) / m(4, 4) ! Upper part. m(5, 15) = m(5, 15) / m(5, 5) ! Upper part. temp = m(6, 5) * m(5, 15) m(6, 15) = ( m(6, 15) - temp ) / m(6, 6) ! Upper part. m(7, 18) = m(7, 18) / m(7, 7) ! Upper part. m(7, 19) = m(7, 19) / m(7, 7) ! Upper part. m(8, 15) = m(8, 15) / m(8, 8) ! Upper part. m(8, 20) = m(8, 20) / m(8, 8) ! Upper part. m(9, 15) = m(9, 15) / m(9, 9) ! Upper part. temp = m(9, 7) * m(7, 18) m(9, 18) = ( m(9, 18) - temp ) / m(9, 9) ! Upper part. temp = m(9, 7) * m(7, 19) m(9, 19) = ( m(9, 19) - temp ) / m(9, 9) ! Upper part. m(10, 12) = m(10, 12) / m(10, 10) ! Upper part. m(10, 20) = m(10, 20) / m(10, 10) ! Upper part. m(11, 12) = m(11, 12) / m(11, 11) ! Upper part. m(11, 19) = m(11, 19) / m(11, 11) ! Lower part. temp = m(12, 11) * m(11, 12) m(12, 12) = m(12, 12) - temp ! Lower part. temp = m(14, 10) * m(10, 12) m(14, 12) = m(14, 12) - temp ! Lower part. temp = m(19, 11) * m(11, 12) m(19, 12) = m(19, 12) - temp ! Upper part. temp = m(12, 11) * m(11, 19) m(12, 19) = ( m(12, 19) - temp ) / m(12, 12) ! Upper part. m(12, 20) = m(12, 20) / m(12, 12) ! Upper part. temp = m(13, 12) * m(12, 19) m(13, 19) = ( m(13, 19) - temp ) / m(13, 13) ! Upper part. temp = m(13, 12) * m(12, 20) m(13, 20) = ( m(13, 20) - temp ) / m(13, 13) ! Lower part. temp = m(15, 3) * m(3, 14) m(15, 14) = m(15, 14) - temp ! Upper part. temp = m(14, 5) * m(5, 15) m(14, 15) = ( m(14, 15) - temp ) / m(14, 14) ! Upper part. temp = m(14, 12) * m(12, 19) m(14, 19) = ( m(14, 19) - temp ) / m(14, 14) ! Upper part. temp = m(14, 10) * m(10, 20) temp = temp + m(14, 12) * m(12, 20) m(14, 20) = ( m(14, 20) - temp ) / m(14, 14) ! Lower part. temp = m(15, 5) * m(5, 15) temp = temp + m(15, 8) * m(8, 15) temp = temp + m(15, 14) * m(14, 15) m(15, 15) = m(15, 15) - temp ! Lower part. temp = m(19, 14) * m(14, 15) m(19, 15) = m(19, 15) - temp ! Lower part. temp = m(20, 8) * m(8, 15) temp = temp + m(20, 14) * m(14, 15) m(20, 15) = m(20, 15) - temp ! Upper part. temp = m(15, 2) * m(2, 16) m(15, 16) = ( m(15, 16) - temp ) / m(15, 15) ! Upper part. temp = m(15, 14) * m(14, 19) m(15, 19) = ( m(15, 19) - temp ) / m(15, 15) ! Upper part. temp = m(15, 8) * m(8, 20) temp = temp + m(15, 14) * m(14, 20) m(15, 20) = ( m(15, 20) - temp ) / m(15, 15) ! Lower part. temp = m(17, 2) * m(2, 16) m(17, 16) = m(17, 16) - temp ! Lower part. temp = m(19, 15) * m(15, 16) m(19, 16) = m(19, 16) - temp ! Lower part. temp = m(20, 15) * m(15, 16) m(20, 16) = m(20, 16) - temp ! Upper part. m(16, 17) = m(16, 17) / m(16, 16) ! Upper part. m(16, 19) = m(16, 19) / m(16, 16) ! Upper part. m(16, 20) = m(16, 20) / m(16, 16) ! Lower part. temp = m(17, 16) * m(16, 17) m(17, 17) = m(17, 17) - temp ! Lower part. temp = m(18, 16) * m(16, 17) m(18, 17) = m(18, 17) - temp ! Lower part. temp = m(19, 16) * m(16, 17) m(19, 17) = m(19, 17) - temp ! Lower part. temp = m(20, 16) * m(16, 17) m(20, 17) = m(20, 17) - temp ! Upper part. m(17, 18) = m(17, 18) / m(17, 17) ! Upper part. temp = m(17, 16) * m(16, 19) m(17, 19) = ( m(17, 19) - temp ) / m(17, 17) ! Upper part. temp = m(17, 16) * m(16, 20) m(17, 20) = ( m(17, 20) - temp ) / m(17, 17) ! Lower part. temp = m(18, 7) * m(7, 18) temp = temp + m(18, 17) * m(17, 18) m(18, 18) = m(18, 18) - temp ! Lower part. temp = m(19, 7) * m(7, 18) temp = temp + m(19, 17) * m(17, 18) m(19, 18) = m(19, 18) - temp ! Lower part. temp = m(20, 17) * m(17, 18) m(20, 18) = m(20, 18) - temp ! Upper part. temp = m(18, 7) * m(7, 19) temp = temp + m(18, 16) * m(16, 19) temp = temp + m(18, 17) * m(17, 19) m(18, 19) = ( m(18, 19) - temp ) / m(18, 18) ! Upper part. temp = m(18, 16) * m(16, 20) temp = temp + m(18, 17) * m(17, 20) m(18, 20) = ( m(18, 20) - temp ) / m(18, 18) ! Lower part. temp = m(19, 7) * m(7, 19) temp = temp + m(19, 11) * m(11, 19) temp = temp + m(19, 12) * m(12, 19) temp = temp + m(19, 13) * m(13, 19) temp = temp + m(19, 14) * m(14, 19) temp = temp + m(19, 15) * m(15, 19) temp = temp + m(19, 16) * m(16, 19) temp = temp + m(19, 17) * m(17, 19) temp = temp + m(19, 18) * m(18, 19) m(19, 19) = m(19, 19) - temp ! Lower part. temp = m(20, 12) * m(12, 19) temp = temp + m(20, 13) * m(13, 19) temp = temp + m(20, 14) * m(14, 19) temp = temp + m(20, 15) * m(15, 19) temp = temp + m(20, 16) * m(16, 19) temp = temp + m(20, 17) * m(17, 19) temp = temp + m(20, 18) * m(18, 19) m(20, 19) = m(20, 19) - temp ! Upper part. temp = m(19, 12) * m(12, 20) temp = temp + m(19, 13) * m(13, 20) temp = temp + m(19, 14) * m(14, 20) temp = temp + m(19, 15) * m(15, 20) temp = temp + m(19, 16) * m(16, 20) temp = temp + m(19, 17) * m(17, 20) temp = temp + m(19, 18) * m(18, 20) m(19, 20) = ( m(19, 20) - temp ) / m(19, 19) ! Lower part. temp = m(20, 8) * m(8, 20) temp = temp + m(20, 12) * m(12, 20) temp = temp + m(20, 13) * m(13, 20) temp = temp + m(20, 14) * m(14, 20) temp = temp + m(20, 15) * m(15, 20) temp = temp + m(20, 16) * m(16, 20) temp = temp + m(20, 17) * m(17, 20) temp = temp + m(20, 18) * m(18, 20) temp = temp + m(20, 19) * m(19, 20) m(20, 20) = m(20, 20) - temp return end subroutine lu_decompose_2 !=============================================================================== !> lu_solve_2 !> \brief Resolution of MY=X where M is an LU factorization computed by !> lu_decompose_2 !------------------------------------------------------------------------------ !------------------------------------------------------------------------------ ! Arguments !------------------------------------------------------------------------------ ! mode name role !------------------------------------------------------------------------------ !> \param[in] ns matrix row number from the chemical species number !> \param[in] m an LU factorization computed by lu_decompose_2 !> \param[in,out] x on entry, the right-hand side of the equation. ! on exit, the solution of the equation !______________________________________________________________________________ subroutine lu_solve_2 (ns, m, x) implicit none ! Arguments integer ns double precision m(ns,ns) double precision x(ns) ! Local variables double precision temp ! Forward substitution. x(1) = x(1) / m(1, 1) x(2) = x(2) / m(2, 2) x(3) = x(3) / m(3, 3) temp = m(4, 1) * x(1) x(4) = ( x(4) - temp ) / m(4, 4) x(5) = x(5) / m(5, 5) temp = m(6, 5) * x(5) x(6) = ( x(6) - temp ) / m(6, 6) x(7) = x(7) / m(7, 7) x(8) = x(8) / m(8, 8) temp = m(9, 7) * x(7) x(9) = ( x(9) - temp ) / m(9, 9) temp = m(10, 1) * x(1) x(10) = ( x(10) - temp ) / m(10, 10) x(11) = x(11) / m(11, 11) temp = m(12, 11) * x(11) x(12) = ( x(12) - temp ) / m(12, 12) temp = m(13, 1) * x(1) temp = temp + m(13, 12) * x(12) x(13) = ( x(13) - temp ) / m(13, 13) temp = m(14, 1) * x(1) temp = temp + m(14, 5) * x(5) temp = temp + m(14, 10) * x(10) temp = temp + m(14, 12) * x(12) x(14) = ( x(14) - temp ) / m(14, 14) temp = m(15, 2) * x(2) temp = temp + m(15, 3) * x(3) temp = temp + m(15, 5) * x(5) temp = temp + m(15, 8) * x(8) temp = temp + m(15, 14) * x(14) x(15) = ( x(15) - temp ) / m(15, 15) x(16) = x(16) / m(16, 16) temp = m(17, 2) * x(2) temp = temp + m(17, 16) * x(16) x(17) = ( x(17) - temp ) / m(17, 17) temp = m(18, 7) * x(7) temp = temp + m(18, 16) * x(16) temp = temp + m(18, 17) * x(17) x(18) = ( x(18) - temp ) / m(18, 18) temp = m(19, 7) * x(7) temp = temp + m(19, 11) * x(11) temp = temp + m(19, 12) * x(12) temp = temp + m(19, 13) * x(13) temp = temp + m(19, 14) * x(14) temp = temp + m(19, 15) * x(15) temp = temp + m(19, 16) * x(16) temp = temp + m(19, 17) * x(17) temp = temp + m(19, 18) * x(18) x(19) = ( x(19) - temp ) / m(19, 19) temp = m(20, 8) * x(8) temp = temp + m(20, 12) * x(12) temp = temp + m(20, 13) * x(13) temp = temp + m(20, 14) * x(14) temp = temp + m(20, 15) * x(15) temp = temp + m(20, 16) * x(16) temp = temp + m(20, 17) * x(17) temp = temp + m(20, 18) * x(18) temp = temp + m(20, 19) * x(19) x(20) = ( x(20) - temp ) / m(20, 20) ! Backward substitution. temp = m(19, 20) * x(20) x(19) = x(19) - temp temp = m(18, 19) * x(19) temp = temp + m(18, 20) * x(20) x(18) = x(18) - temp temp = m(17, 18) * x(18) temp = temp + m(17, 19) * x(19) temp = temp + m(17, 20) * x(20) x(17) = x(17) - temp temp = m(16, 17) * x(17) temp = temp + m(16, 19) * x(19) temp = temp + m(16, 20) * x(20) x(16) = x(16) - temp temp = m(15, 16) * x(16) temp = temp + m(15, 19) * x(19) temp = temp + m(15, 20) * x(20) x(15) = x(15) - temp temp = m(14, 15) * x(15) temp = temp + m(14, 19) * x(19) temp = temp + m(14, 20) * x(20) x(14) = x(14) - temp temp = m(13, 19) * x(19) temp = temp + m(13, 20) * x(20) x(13) = x(13) - temp temp = m(12, 19) * x(19) temp = temp + m(12, 20) * x(20) x(12) = x(12) - temp temp = m(11, 12) * x(12) temp = temp + m(11, 19) * x(19) x(11) = x(11) - temp temp = m(10, 12) * x(12) temp = temp + m(10, 20) * x(20) x(10) = x(10) - temp temp = m(9, 15) * x(15) temp = temp + m(9, 18) * x(18) temp = temp + m(9, 19) * x(19) x(9) = x(9) - temp temp = m(8, 15) * x(15) temp = temp + m(8, 20) * x(20) x(8) = x(8) - temp temp = m(7, 18) * x(18) temp = temp + m(7, 19) * x(19) x(7) = x(7) - temp temp = m(6, 15) * x(15) x(6) = x(6) - temp temp = m(5, 15) * x(15) x(5) = x(5) - temp temp = m(4, 10) * x(10) x(4) = x(4) - temp temp = m(3, 14) * x(14) x(3) = x(3) - temp temp = m(2, 16) * x(16) x(2) = x(2) - temp return end subroutine lu_solve_2