84   subroutine ShashiLowerBound(an_r)  argument
85     PetscScalar an_r(26)
89       an_r(i) = 1000.0/6.023D+23
93   subroutine ShashiInitialGuess(an_r)  argument
105     PetscScalar an_r(26)
123     an_r(1) = 1.66000D-24
124     an_r(2) = 1.66030D-22
125     an_r(3) = 5.00000D-01
126     an_r(4) = 1.66030D-22
127     an_r(5) = 1.66030D-22
128     an_r(6) = 1.88000D+00
129     an_r(7) = 1.66030D-22
130     an_r(8) = 1.66030D-22
131     an_r(9) = 1.66030D-22
132     an_r(10) = 1.66030D-22
133     an_r(11) = 1.66030D-22
134     an_r(12) = 1.66030D-22
135     an_r(13) = 1.66030D-22
136     an_r(14) = 1.00000D+00
137     an_r(15) = 1.66030D-22
138     an_r(16) = 1.66030D-22
139     an_r(17) = 1.66000D-24
140     an_r(18) = 1.66030D-24
141     an_r(19) = 1.66030D-24
142     an_r(20) = 1.66030D-24
143     an_r(21) = 1.66030D-24
144     an_r(22) = 1.66030D-24
145     an_r(23) = 1.66030D-24
146     an_r(24) = 1.66030D-24
147     an_r(25) = 1.66030D-24
148     an_r(26) = 1.66030D-24
150     an_r = 0
151     an_r(3) = .5
152     an_r(6) = 1.88000
153     an_r(14) = 1.
156     an_r(1) = 3.802208D-33
157     an_r(2) = 1.298287D-29
158     an_r(3) = 2.533067D-04
159     an_r(4) = 6.865078D-22
160     an_r(5) = 9.993125D-01
161     an_r(6) = 1.879964D+00
162     an_r(7) = 4.449489D-13
163     an_r(8) = 3.428687D-07
164     an_r(9) = 7.105138D-05
165     an_r(10) = 1.094368D-04
166     an_r(11) = 2.362305D-06
167     an_r(12) = 1.107145D-09
168     an_r(13) = 1.276162D-24
169     an_r(14) = 6.315538D-04
170     an_r(15) = 2.356540D-09
171     an_r(16) = 2.048248D-09
172     an_r(17) = 1.966187D-22
173     an_r(18) = 7.856497D-29
174     an_r(19) = 1.987840D-36
175     an_r(20) = 8.182441D-22
176     an_r(21) = 2.684880D-16
177     an_r(22) = 2.680473D-16
178     an_r(23) = 6.594967D-18
179     an_r(24) = 2.509714D-21
180     an_r(25) = 3.096459D-21
181     an_r(26) = 6.149551D-18
186   subroutine ShashiFormFunction(an_r, f_eq)  argument
198     PetscScalar an_r(26), k_eq(23), f_eq(26)
276       an_t = an_t + an_r(i)
280               - (an_h(1)*an_r(1) + an_h_additive*an_r(2) &
281                  + 2*an_r(5) + an_r(10) + an_r(11) + 2*an_r(14) &
282                  + an_r(16) + 2*an_r(17) + an_r(19) &
283                  + an_r(20) + 3*an_r(22) + an_r(26))
286               - (an_o_additive*an_r(2) + 2*an_r(3) &
287                  + 2*an_r(4) + an_r(5) &
288                  + an_r(8) + an_r(9) + an_r(10) + an_r(12) + an_r(13) &
289                  + 2*an_r(15) + 2*an_r(16) + an_r(20) + an_r(22) &
290                  + an_r(23) + 2*an_r(24) + 1*an_r(25) + an_r(26))
292     f_eq(3) = an_r(2) - 1.0d-150
295               - (an_c(1)*an_r(1) + an_c_additive*an_r(2) &
296                  + an_r(4) + an_r(13) + 2*an_r(17) + an_r(18) &
297                  + an_r(19) + an_r(20))
300       part_p(ip) = (an_r(ip)/an_t)*pt
309     f_eq(5) = k_eq(11)*an_r(1)*an_r(3)**a_io2 &
310               - (an_r(4)**i_cc)*(an_r(5)**i_h2o)*((pt/an_t)**idiff)
314               - (2*an_r(6) + an_r(7) + an_r(9) + 2*an_r(12) &
315                  + an_r(15) &
316                  + an_r(23))
373     f_eq(26) = -(an_r(20) + an_r(22) + an_r(23)) &
374                + (an_r(21) + an_r(24) + an_r(25) + an_r(26))
382   subroutine ShashiFormJacobian(an_r, d_eq)  argument
404     PetscScalar an_r(26), k_eq(23)
491       an_t = an_t + an_r(id)
576       -(an_r(4)**i_cc)*(an_r(5)**i_h2o) &
583     d_eq(5, 1) = d_eq(5, 1) + k_eq(11)*(an_r(3)**ai_o2)
585     d_eq(5, 3) = d_eq(5, 3) + k_eq(11)*(ai_o2*an_r(3)**(ai_o2 - 1)) &
586                 &                           *an_r(1)
588     d_eq(5, 4) = d_eq(5, 4) - (i_cc*an_r(4)**(i_cc - 1))* &
589                  (an_r(5)**i_h2o)*((pt/an_t)**idiff)
591                  - (i_h2o*(an_r(5)**(i_h2o - 1))) &
592                  *(an_r(4)**i_cc)*((pt/an_t)**idiff)
594     d_eq(3, 1) = -(an_r(4)**2)*(an_r(5)**3)*(pt/an_t)*(-1.0/an_t)
599     d_eq(3, 2) = d_eq(3, 2) + k_eq(12)*(an_r(3)**3)
601     d_eq(3, 3) = d_eq(3, 3) + k_eq(12)*(3*an_r(3)**2)*an_r(2)
603     d_eq(3, 4) = d_eq(3, 4) - 2*an_r(4)*(an_r(5)**3)*(pt/an_t)
605     d_eq(3, 5) = d_eq(3, 5) - 3*(an_r(5)**2)*(an_r(4)**2)*(pt/an_t)
613     d_eq(7, 1) = pt*an_r(11)*(-1.0)/const2 &
614                  - k_eq(1)*sqrt(pt)*sqrt(an_r(14) + 1d-50)*(-0.5/const3)
622                   - k_eq(1)*sqrt(pt)*(0.5/(sqrt((an_r(14) + 1d-50)*an_t)))
624     d_eq(8, 1) = pt*an_r(8)*(-1.0)/const2 &
625                  - k_eq(2)*sqrt(pt)*sqrt(an_r(3) + 1.0d-50)*(-0.5/const3)
632                  - k_eq(2)*sqrt(pt)*(0.5/(sqrt((an_r(3) + 1.0d-50)*an_t)))
635     d_eq(9, 1) = pt*an_r(7)*(-1.0)/const2 &
636                  - k_eq(3)*sqrt(pt)*sqrt(an_r(6))*(-0.5/const3)
644                  - k_eq(3)*sqrt(pt)*(0.5/(sqrt(an_r(6)*an_t)))
646     d_eq(10, 1) = pt*an_r(10)*(-1.0)/const2 &
647                   - k_eq(4)*(pt)*sqrt((an_r(3) + 1.0d-50) &
648                                       *an_r(14))*(-1.0/const2)
654                   - k_eq(4)*(pt)*sqrt(an_r(14)) &
655                   *(0.5/(sqrt(an_r(3) + 1.0d-50)*an_t))
658                    - k_eq(4)*(pt)*sqrt(an_r(3) + 1.0d-50) &
659                    *(0.5/(sqrt(an_r(14) + 1.0d-50)*an_t))
661     d_eq(11, 1) = pt*an_r(9)*(-1.0)/const2 &
662                   - k_eq(5)*(pt)*sqrt((an_r(3) + 1.0d-50)*an_r(6)) &
670                   - k_eq(5)*(pt)*sqrt(an_r(6))*(0.5/ &
671                                                 (sqrt(an_r(3) + 1.0d-50)*an_t))
673                   - k_eq(5)*(pt)*sqrt(an_r(3) + 1.0d-50) &
674                   *(0.5/(sqrt(an_r(6))*an_t))
677     d_eq(12, 1) = pt*an_r(5)*(-1.0)/const2 &
678                   - k_eq(6)*(pt**1.5)*sqrt(an_r(3) + 1.0d-50) &
679                   *(an_r(14))*(-1.5/const5)
686                   - k_eq(6)*(pt**1.5)*((an_r(14) + 1.0d-50)/const3) &
687                   *(0.5/sqrt(an_r(3) + 1.0d-50))
691                    - k_eq(6)*(pt**1.5)*(sqrt(an_r(3) + 1.0d-50)/const3)
693     d_eq(13, 1) = pt*an_r(4)*(-1.0)/const2 &
694                   - k_eq(7)*(pt**1.5)*sqrt(an_r(3) + 1.0d-50) &
695                   *(an_r(13))*(-1.5/const5)
702                   - k_eq(7)*(pt**1.5)*(an_r(13)/const3) &
703                   *(0.5/sqrt(an_r(3) + 1.0d-50))
707                    - k_eq(7)*(pt**1.5)*(sqrt(an_r(3) + 1.0d-50)/const3)
709     d_eq(14, 1) = pt*an_r(15)*(-1.0)/const2 &
710                   - k_eq(8)*(pt**1.5)*sqrt(an_r(3) + 1.0d-50) &
711                   *(an_r(9))*(-1.5/const5)
718                   - k_eq(8)*(pt**1.5)*(an_r(9)/const3) &
719                   *(0.5/sqrt(an_r(3) + 1.0d-50))
721                   - k_eq(8)*(pt**1.5)*(sqrt(an_r(3) + 1.0d-50)/const3)
724     d_eq(15, 1) = pt*an_r(16)*(-1.0)/const2 &
725                   - k_eq(9)*(pt**1.5)*sqrt(an_r(14) + 1.0d-50) &
726                   *(an_r(3))*(-1.5/const5)
733                   - k_eq(9)*(pt**1.5)*(sqrt(an_r(14) + 1.0d-50)/const3)
735                    - k_eq(9)*(pt**1.5)*(an_r(3)/const3) &
736                    *(0.5/sqrt(an_r(14) + 1.0d-50))
739     d_eq(16, 1) = pt*an_r(12)*(-1.0)/const2 &
740                   - k_eq(10)*(pt**1.5)*sqrt(an_r(3) + 1.0d-50) &
741                   *(an_r(6))*(-1.5/const5)
748                   - k_eq(10)*(pt**1.5)*(an_r(6)/const3) &
749                   *(0.5/sqrt(an_r(3) + 1.0d-50))
752                   - k_eq(10)*(pt**1.5)*(sqrt(an_r(3) + 1.0d-50)/const3)
758     d_eq(17, 1) = an_r(14)*an_r(18)*an_r(18)*(pt**3)*(-3/const_four) &
759                   - k_eq(15)*an_r(17)*pt*(-1/const2)
763     d_eq(17, 14) = d_eq(17, 14) + an_r(18)*an_r(18)*(pt**3)/const_cube
765     d_eq(17, 18) = d_eq(17, 18) + 2*an_r(18)*an_r(14)                 &
768     d_eq(18, 1) = an_r(13)*an_r(13)*(pt**2)*(-2/const_cube) &
769                   - k_eq(16)*an_r(3)*an_r(18)*an_r(18) &
775                   - k_eq(16)*an_r(18)*an_r(18)*pt*pt*pt/const_cube
777                    + 2*an_r(13)*pt*pt/const2
778     d_eq(18, 18) = d_eq(18, 18) - k_eq(16)*an_r(3)                     &
779                   &              *2*an_r(18)*pt*pt*pt/const_cube
783     d_eq(19, 1) = an_r(3)*an_r(19)*(pt**2)*(-2/const_cube) &
784                   - k_eq(17)*an_r(13)*an_r(10)*pt*pt*(-2/const_cube)
789                    - k_eq(17)*an_r(10)*pt*pt/const2
791                    - k_eq(17)*an_r(13)*pt*pt/const2
792     d_eq(19, 3) = d_eq(19, 3) + an_r(19)*pt*pt/const2
793     d_eq(19, 19) = d_eq(19, 19) + an_r(3)*pt*pt/const2
796     d_eq(20, 1) = an_r(21)*an_r(20)*(pt**2)*(-2/const_cube) &
797                   - k_eq(18)*an_r(19)*an_r(8)*pt*pt*(-2/const_cube)
802                   - k_eq(18)*an_r(19)*pt*pt/const2
804                    - k_eq(18)*an_r(8)*pt*pt/const2
805     d_eq(20, 20) = d_eq(20, 20) + an_r(21)*pt*pt/const2
806     d_eq(20, 21) = d_eq(20, 21) + an_r(20)*pt*pt/const2
811     d_eq(21, 1) = an_r(21)*an_r(23)*(pt**2)*(-2/const_cube) &
812                   - k_eq(19)*an_r(7)*an_r(8)*pt*pt*(-2/const_cube)
817                   - k_eq(19)*an_r(8)*pt*pt/const2
819                   - k_eq(19)*an_r(7)*pt*pt/const2
820     d_eq(21, 21) = d_eq(21, 21) + an_r(23)*pt*pt/const2
821     d_eq(21, 23) = d_eq(21, 23) + an_r(21)*pt*pt/const2
825     d_eq(22, 1) = an_r(5)*an_r(11)*(pt**2)*(-2/const_cube) &
826                   - k_eq(20)*an_r(21)*an_r(22)*pt*pt*(-2/const_cube)
831                    - k_eq(20)*an_r(22)*pt*pt/const2
833                    - k_eq(20)*an_r(21)*pt*pt/const2
834     d_eq(22, 11) = d_eq(22, 11) + an_r(5)*pt*pt/(const2)
835     d_eq(22, 5) = d_eq(22, 5) + an_r(11)*pt*pt/(const2)
840     d_eq(23, 1) = an_r(24)*(pt)*(-1/const2) &
841                   - k_eq(21)*an_r(21)*an_r(3)*pt*pt*(-2/const_cube)
846                   - k_eq(21)*an_r(21)*pt*pt/const2
848                    - k_eq(21)*an_r(3)*pt*pt/const2
853     d_eq(24, 1) = an_r(3)*an_r(25)*(pt**2)*(-2/const_cube) &
854                   - k_eq(22)*an_r(24)*an_r(8)*pt*pt*(-2/const_cube)
859                   - k_eq(22)*an_r(24)*pt*pt/const2
861                    - k_eq(22)*an_r(8)*pt*pt/const2
862     d_eq(24, 3) = d_eq(24, 3) + an_r(25)*pt*pt/const2
863     d_eq(24, 25) = d_eq(24, 25) + an_r(3)*pt*pt/const2
868     d_eq(25, 1) = an_r(26)*(pt)*(-1/const2) &
869                   - k_eq(23)*an_r(21)*an_r(10)*pt*pt*(-2/const_cube)
874                    - k_eq(23)*an_r(21)*pt*pt/const2
876                    - k_eq(23)*an_r(10)*pt*pt/const2