chem/KPP/kpp/kpp-2.1/int/atm_radau5.f

   1       SUBROUTINE INTEGRATE( TIN, TOUT )
   2
   3       IMPLICIT KPP_REAL (A-H,O-Z)
   4       INCLUDE 'KPP_ROOT_Parameters.h'
   5       INCLUDE 'KPP_ROOT_Global.h'
   6
   7 C TIN - Start Time
   8       KPP_REAL TIN
   9 C TOUT - End Time
  10       KPP_REAL TOUT, H
  11       SAVE H
  12       DATA H /0.0d0/
  13       INTEGER Nfun, Njac, Nstp, Nacc, Nrej, Ndec, Nsol
  14       SAVE Nstp, Nacc, Nrej
  15       DATA Nstp /0/
  16       DATA Nacc /0/
  17       DATA Nrej /0/
  18       INTEGER i
  19
  20       PARAMETER (LWORK=5*NVAR*NVAR+12*NVAR+20,LIWORK=3*NVAR+20)
  21       PARAMETER (LRCONT=4*NVAR+4)
  22
  23       KPP_REAL WORK(LWORK), RPAR(1)
  24       INTEGER IWORK(LIWORK), IPAR(1)
  25       COMMON /CONT/ ICONT(4),RCONT(LRCONT)
  26       EXTERNAL FUNC_CHEM,JAC_CHEM,SOLOUT
  27
  28
  29       ITOL=1     ! --- VECTOR TOLERANCES
  30       IJAC=1     ! --- COMPUTE THE JACOBIAN ANALYTICALLY
  31       IMAS=0     ! --- DIFFERENTIAL EQUATION IS IN EXPLICIT FORM
  32       IOUT=0     ! --- OUTPUT ROUTINE IS NOT USED DURING INTEGRATION
  33       MLJAC=NVAR ! --- JACOBIAN IS A FULL MATRIX
  34       MUJAC=NVAR ! --- JACOBIAN IS A FULL MATRIX
  35       MLMAS=NVAR ! --- JACOBIAN IS A FULL MATRIX
  36       MUMAS=NVAR ! --- JACOBIAN IS A FULL MATRIX
  37
  38       DO i=1,20
  39         IWORK(i) = 0
  40         WORK(i) = 0.D0
  41       ENDDO
  42
  43       IWORK(3) = 8    !  Max no. of Newton iterations
  44       IWORK(4) = 0    !  Starting values for Newton are interpolated (0) or zero (1)
  45       IWORK(8) = 2    ! Gustaffson (1) or classic(2) controller
  46       WORK(2)  = 0.9  ! Safety factor
  47
  48       CALL RADAU5(NVAR,FUNC_CHEM,TIN,VAR,TOUT,H,
  49      &                  RTOL,ATOL,ITOL,
  50      &                  JAC_CHEM ,IJAC,MLJAC,MUJAC,
  51      &                  FUNC_CHEM ,IMAS,MLMAS,MUMAS,
  52      &                  SOLOUT,IOUT,
  53      &                  WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID)
  54
  55
  56        Nfun = Nfun + IWORK(14)
  57        Njac = Njac + IWORK(15)
  58        Nstp = Nstp + IWORK(16)
  59        Nacc = Nacc + IWORK(17)
  60        Nrej = Nrej + IWORK(18)
  61        Ndec = Ndec + IWORK(19)
  62        Nsol = Nsol + IWORK(20)
  63
  64        print("('Nstep=',I5,' Nacc=',I6,' Nrej=',I6)"),Nstp, Nacc, Nrej
  65
  66       IF (IDID.LT.0) THEN
  67         print *,'RADAU: Unsucessfull exit at T=',
  68      &          TIN,' (IDID=',IDID,')'
  69       ENDIF
  70
  71       RETURN
  72       END
  73
  74
  75       SUBROUTINE RADAU5(N,FCN,X,Y,XEND,H,
  76      &                  RelTol,AbsTol,ITOL,
  77      &                  JAC ,IJAC,MLJAC,MUJAC,
  78      &                  MAS ,IMAS,MLMAS,MUMAS,
  79      &                  SOLOUT,IOUT,
  80      &                  WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID)
  81 C ----------------------------------------------------------
  82 C     NUMERICAL SOLUTION OF A STIFF (OR DIFFERENTIAL ALGEBRAIC)
  83 C     SYSTEM OF FIRST 0RDER ORDINARY DIFFERENTIAL EQUATIONS
  84 C                     M*Y'=F(X,Y).
  85 C     THE SYSTEM CAN BE (LINEARLY) IMPLICIT (MASS-MATRIX M .NE. I)
  86 C     OR EXPLICIT (M=I).
  87 C     THE METHOD USED IS AN IMPLICIT RUNGE-KUTTA METHOD (RADAU IIA)
  88 C     OF ORDER 5 WITH STEP SIZE CONTROL AND CONTINUOUS OUTPUT.
  89 C     C.F. SECTION IV.8
  90 C
  91 C     AUTHORS: E. HAIRER AND G. WANNER
  92 C              UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES
  93 C              CH-1211 GENEVE 24, SWITZERLAND
  94 C              E-MAIL:  HAIRER@DIVSUN.UNIGE.CH,  WANNER@DIVSUN.UNIGE.CH
  95 C
  96 C     THIS CODE IS PART OF THE BOOK:
  97 C         E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
  98 C         EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
  99 C         SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS 14,
 100 C         SPRINGER-VERLAG (1991)
 101 C
 102 C     VERSION OF SEPTEMBER 30, 1995
 103 C
 104 C     INPUT PARAMETERS
 105 C     ----------------
 106 C     N           DIMENSION OF THE SYSTEM
 107 C
 108 C     FCN         NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE
 109 C                 VALUE OF F(X,Y):
 110 C                    SUBROUTINE FCN(N,X,Y,F,RPAR,IPAR)
 111 C                    KPP_REAL X,Y(N),F(N)
 112 C                    F(1)=...   ETC.
 113 C                 RPAR, IPAR (SEE BELOW)
 114 C
 115 C     X           INITIAL X-VALUE
 116 C
 117 C     Y(N)        INITIAL VALUES FOR Y
 118 C
 119 C     XEND        FINAL X-VALUE (XEND-X MAY BE POSITIVE OR NEGATIVE)
 120 C
 121 C     H           INITIAL STEP SIZE GUESS;
 122 C                 FOR STIFF EQUATIONS WITH INITIAL TRANSIENT,
 123 C                 H=1.D0/(NORM OF F'), USUALLY 1.D-3 OR 1.D-5, IS GOOD.
 124 C                 THIS CHOICE IS NOT VERY IMPORTANT, THE STEP SIZE IS
 125 C                 QUICKLY ADAPTED. (IF H=0.D0, THE CODE PUTS H=1.D-6).
 126 C
 127 C     RelTol,AbsTol   RELATIVE AND ABSOLUTE ERROR TOLERANCES. THEY
 128 C                 CAN BE BOTH SCALARS OR ELSE BOTH VECTORS OF LENGTH N.
 129 C
 130 C     ITOL        SWITCH FOR RelTol AND AbsTol:
 131 C                   ITOL=0: BOTH RelTol AND AbsTol ARE SCALARS.
 132 C                     THE CODE KEEPS, ROUGHLY, THE LOCAL ERROR OF
 133 C                     Y(I) BELOW RelTol*ABS(Y(I))+AbsTol
 134 C                   ITOL=1: BOTH RelTol AND AbsTol ARE VECTORS.
 135 C                     THE CODE KEEPS THE LOCAL ERROR OF Y(I) BELOW
 136 C                     RelTol(I)*ABS(Y(I))+AbsTol(I).
 137 C
 138 C     JAC         NAME (EXTERNAL) OF THE SUBROUTINE WHICH COMPUTES
 139 C                 THE PARTIAL DERIVATIVES OF F(X,Y) WITH RESPECT TO Y
 140 C                 (THIS ROUTINE IS ONLY CALLED IF IJAC=1; SUPPLY
 141 C                 A DUMMY SUBROUTINE IN THE CASE IJAC=0).
 142 C                 FOR IJAC=1, THIS SUBROUTINE MUST HAVE THE FORM
 143 C                    SUBROUTINE JAC(N,X,Y,DFY,LDFY,RPAR,IPAR)
 144 C                    KPP_REAL X,Y(N),DFY(LDFY,N)
 145 C                    DFY(1,1)= ...
 146 C                 LDFY, THE COLUMN-LENGTH OF THE ARRAY, IS
 147 C                 FURNISHED BY THE CALLING PROGRAM.
 148 C                 IF (MLJAC.EQ.N) THE JACOBIAN IS SUPPOSED TO
 149 C                    BE FULL AND THE PARTIAL DERIVATIVES ARE
 150 C                    STORED IN DFY AS
 151 C                       DFY(I,J) = PARTIAL F(I) / PARTIAL Y(J)
 152 C                 ELSE, THE JACOBIAN IS TAKEN AS BANDED AND
 153 C                    THE PARTIAL DERIVATIVES ARE STORED
 154 C                    DIAGONAL-WISE AS
 155 C                       DFY(I-J+MUJAC+1,J) = PARTIAL F(I) / PARTIAL Y(J).
 156 C
 157 C     IJAC        SWITCH FOR THE COMPUTATION OF THE JACOBIAN:
 158 C                    IJAC=0: JACOBIAN IS COMPUTED INTERNALLY BY FINITE
 159 C                       DIFFERENCES, SUBROUTINE "JAC" IS NEVER CALLED.
 160 C                    IJAC=1: JACOBIAN IS SUPPLIED BY SUBROUTINE JAC.
 161 C
 162 C     MLJAC       SWITCH FOR THE BANDED STRUCTURE OF THE JACOBIAN:
 163 C                    MLJAC=N: JACOBIAN IS A FULL MATRIX. THE LINEAR
 164 C                       ALGEBRA IS DONE BY FULL-MATRIX GAUSS-ELIMINATION.
 165 C                    0<=MLJAC<N: MLJAC IS THE LOWER BANDWITH OF JACOBIAN
 166 C                       MATRIX (>= NUMBER OF NON-ZERO DIAGONALS BELOW
 167 C                       THE MAIN DIAGONAL).
 168 C
 169 C     MUJAC       UPPER BANDWITH OF JACOBIAN  MATRIX (>= NUMBER OF NON-
 170 C                 ZERO DIAGONALS ABOVE THE MAIN DIAGONAL).
 171 C                 NEED NOT BE DEFINED IF MLJAC=N.
 172 C
 173 C     ----   MAS,IMAS,MLMAS, AND MUMAS HAVE ANALOG MEANINGS      -----
 174 C     ----   FOR THE "MASS MATRIX" (THE MATRIX "M" OF SECTION IV.8): -
 175 C
 176 C     MAS         NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE MASS-
 177 C                 MATRIX M.
 178 C                 IF IMAS=0, THIS MATRIX IS ASSUMED TO BE THE IDENTITY
 179 C                 MATRIX AND NEEDS NOT TO BE DEFINED;
 180 C                 SUPPLY A DUMMY SUBROUTINE IN THIS CASE.
 181 C                 IF IMAS=1, THE SUBROUTINE MAS IS OF THE FORM
 182 C                    SUBROUTINE MAS(N,AM,LMAS,RPAR,IPAR)
 183 C                    KPP_REAL AM(LMAS,N)
 184 C                    AM(1,1)= ....
 185 C                    IF (MLMAS.EQ.N) THE MASS-MATRIX IS STORED
 186 C                    AS FULL MATRIX LIKE
 187 C                         AM(I,J) = M(I,J)
 188 C                    ELSE, THE MATRIX IS TAKEN AS BANDED AND STORED
 189 C                    DIAGONAL-WISE AS
 190 C                         AM(I-J+MUMAS+1,J) = M(I,J).
 191 C
 192 C     IMAS       GIVES INFORMATION ON THE MASS-MATRIX:
 193 C                    IMAS=0: M IS SUPPOSED TO BE THE IDENTITY
 194 C                       MATRIX, MAS IS NEVER CALLED.
 195 C                    IMAS=1: MASS-MATRIX  IS SUPPLIED.
 196 C
 197 C     MLMAS       SWITCH FOR THE BANDED STRUCTURE OF THE MASS-MATRIX:
 198 C                    MLMAS=N: THE FULL MATRIX CASE. THE LINEAR
 199 C                       ALGEBRA IS DONE BY FULL-MATRIX GAUSS-ELIMINATION.
 200 C                    0<=MLMAS<N: MLMAS IS THE LOWER BANDWITH OF THE
 201 C                       MATRIX (>= NUMBER OF NON-ZERO DIAGONALS BELOW
 202 C                       THE MAIN DIAGONAL).
 203 C                 MLMAS IS SUPPOSED TO BE .LE. MLJAC.
 204 C
 205 C     MUMAS       UPPER BANDWITH OF MASS-MATRIX (>= NUMBER OF NON-
 206 C                 ZERO DIAGONALS ABOVE THE MAIN DIAGONAL).
 207 C                 NEED NOT BE DEFINED IF MLMAS=N.
 208 C                 MUMAS IS SUPPOSED TO BE .LE. MUJAC.
 209 C
 210 C     SOLOUT      NAME (EXTERNAL) OF SUBROUTINE PROVIDING THE
 211 C                 NUMERICAL SOLUTION DURING INTEGRATION.
 212 C                 IF IOUT=1, IT IS CALLED AFTER EVERY SUCCESSFUL STEP.
 213 C                 SUPPLY A DUMMY SUBROUTINE IF IOUT=0.
 214 C                 IT MUST HAVE THE FORM
 215 C                    SUBROUTINE SOLOUT (NR,XOLD,X,Y,CONT,LRC,N,
 216 C                                       RPAR,IPAR,IRTRN)
 217 C                    KPP_REAL X,Y(N),CONT(LRC)
 218 C                    ....
 219 C                 SOLOUT FURNISHES THE SOLUTION "Y" AT THE NR-TH
 220 C                    GRID-POINT "X" (THEREBY THE INITIAL VALUE IS
 221 C                    THE FIRST GRID-POINT).
 222 C                 "XOLD" IS THE PRECEEDING GRID-POINT.
 223 C                 "IRTRN" SERVES TO INTERRUPT THE INTEGRATION. IF IRTRN
 224 C                    IS SET <0, RADAU5 RETURNS TO THE CALLING PROGRAM.
 225 C
 226 C          -----  CONTINUOUS OUTPUT: -----
 227 C                 DURING CALLS TO "SOLOUT", A CONTINUOUS SOLUTION
 228 C                 FOR THE INTERVAL [XOLD,X] IS AVAILABLE THROUGH
 229 C                 THE FUNCTION
 230 C                        >>>   CONTR5(I,S,CONT,LRC)   <<<
 231 C                 WHICH PROVIDES AN APPROXIMATION TO THE I-TH
 232 C                 COMPONENT OF THE SOLUTION AT THE POINT S. THE VALUE
 233 C                 S SHOULD LIE IN THE INTERVAL [XOLD,X].
 234 C                 DO NOT CHANGE THE ENTRIES OF CONT(LRC), IF THE
 235 C                 DENSE OUTPUT FUNCTION IS USED.
 236 C
 237 C     IOUT        SWITCH FOR CALLING THE SUBROUTINE SOLOUT:
 238 C                    IOUT=0: SUBROUTINE IS NEVER CALLED
 239 C                    IOUT=1: SUBROUTINE IS AVAILABLE FOR OUTPUT.
 240 C
 241 C     WORK        ARRAY OF WORKING SPACE OF LENGTH "LWORK".
 242 C                 WORK(1), WORK(2),.., WORK(20) SERVE AS PARAMETERS
 243 C                 FOR THE CODE. FOR STANDARD USE OF THE CODE
 244 C                 WORK(1),..,WORK(20) MUST BE SET TO ZERO BEFORE
 245 C                 CALLING. SEE BELOW FOR A MORE SOPHISTICATED USE.
 246 C                 WORK(21),..,WORK(LWORK) SERVE AS WORKING SPACE
 247 C                 FOR ALL VECTORS AND MATRICES.
 248 C                 "LWORK" MUST BE AT LEAST
 249 C                             N*(LJAC+LMAS+3*LE+12)+20
 250 C                 WHERE
 251 C                    LJAC=N              IF MLJAC=N (FULL JACOBIAN)
 252 C                    LJAC=MLJAC+MUJAC+1  IF MLJAC<N (BANDED JAC.)
 253 C                 AND
 254 C                    LMAS=0              IF IMAS=0
 255 C                    LMAS=N              IF IMAS=1 AND MLMAS=N (FULL)
 256 C                    LMAS=MLMAS+MUMAS+1  IF MLMAS<N (BANDED MASS-M.)
 257 C                 AND
 258 C                    LE=N               IF MLJAC=N (FULL JACOBIAN)
 259 C                    LE=2*MLJAC+MUJAC+1 IF MLJAC<N (BANDED JAC.)
 260 C
 261 C                 IN THE USUAL CASE WHERE THE JACOBIAN IS FULL AND THE
 262 C                 MASS-MATRIX IS THE INDENTITY (IMAS=0), THE MINIMUM
 263 C                 STORAGE REQUIREMENT IS
 264 C                             LWORK = 4*N*N+12*N+20.
 265 C                 IF IWORK(9)=M1>0 THEN "LWORK" MUST BE AT LEAST
 266 C                          N*(LJAC+12)+(N-M1)*(LMAS+3*LE)+20
 267 C                 WHERE IN THE DEFINITIONS OF LJAC, LMAS AND LE THE
 268 C                 NUMBER N CAN BE REPLACED BY N-M1.
 269 C
 270 C     LWORK       DECLARED LENGTH OF ARRAY "WORK".
 271 C
 272 C     IWORK       INTEGER WORKING SPACE OF LENGTH "LIWORK".
 273 C                 IWORK(1),IWORK(2),...,IWORK(20) SERVE AS PARAMETERS
 274 C                 FOR THE CODE. FOR STANDARD USE, SET IWORK(1),..,
 275 C                 IWORK(20) TO ZERO BEFORE CALLING.
 276 C                 IWORK(21),...,IWORK(LIWORK) SERVE AS WORKING AREA.
 277 C                 "LIWORK" MUST BE AT LEAST 3*N+20.
 278 C
 279 C     LIWORK      DECLARED LENGTH OF ARRAY "IWORK".
 280 C
 281 C     RPAR, IPAR  REAL AND INTEGER PARAMETERS (OR PARAMETER ARRAYS) WHICH
 282 C                 CAN BE USED FOR COMMUNICATION BETWEEN YOUR CALLING
 283 C                 PROGRAM AND THE FCN, JAC, MAS, SOLOUT SUBROUTINES.
 284 C
 285 C ----------------------------------------------------------------------
 286 C
 287 C     SOPHISTICATED SETTING OF PARAMETERS
 288 C     -----------------------------------
 289 C              SEVERAL PARAMETERS OF THE CODE ARE TUNED TO MAKE IT WORK
 290 C              WELL. THEY MAY BE DEFINED BY SETTING WORK(1),...
 291 C              AS WELL AS IWORK(1),... DIFFERENT FROM ZERO.
 292 C              FOR ZERO INPUT, THE CODE CHOOSES DEFAULT VALUES:
 293 C
 294 C    IWORK(1)  IF IWORK(1).NE.0, THE CODE TRANSFORMS THE JACOBIAN
 295 C              MATRIX TO HESSENBERG FORM. THIS IS PARTICULARLY
 296 C              ADVANTAGEOUS FOR LARGE SYSTEMS WITH FULL JACOBIAN.
 297 C              IT DOES NOT WORK FOR BANDED JACOBIAN (MLJAC<N)
 298 C              AND NOT FOR IMPLICIT SYSTEMS (IMAS=1).
 299 C
 300 C    IWORK(2)  THIS IS THE MAXIMAL NUMBER OF ALLOWED STEPS.
 301 C              THE DEFAULT VALUE (FOR IWORK(2)=0) IS 100000.
 302 C
 303 C    IWORK(3)  THE MAXIMUM NUMBER OF NEWTON ITERATIONS FOR THE
 304 C              SOLUTION OF THE IMPLICIT SYSTEM IN EACH STEP.
 305 C              THE DEFAULT VALUE (FOR IWORK(3)=0) IS 7.
 306 C
 307 C    IWORK(4)  IF IWORK(4).EQ.0 THE EXTRAPOLATED COLLOCATION SOLUTION
 308 C              IS TAKEN AS STARTING VALUE FOR NEWTON'S METHOD.
 309 C              IF IWORK(4).NE.0 ZERO STARTING VALUES ARE USED.
 310 C              THE LATTER IS RECOMMENDED IF NEWTON'S METHOD HAS
 311 C              DIFFICULTIES WITH CONVERGENCE (THIS IS THE CASE WHEN
 312 C              NSTEP IS LARGER THAN NACCPT + NREJCT; SEE OUTPUT PARAM.).
 313 C              DEFAULT IS IWORK(4)=0.
 314 C
 315 C       THE FOLLOWING 3 PARAMETERS ARE IMPORTANT FOR
 316 C       DIFFERENTIAL-ALGEBRAIC SYSTEMS OF INDEX > 1.
 317 C       THE FUNCTION-SUBROUTINE SHOULD BE WRITTEN SUCH THAT
 318 C       THE INDEX 1,2,3 VARIABLES APPEAR IN THIS ORDER.
 319 C       IN ESTIMATING THE ERROR THE INDEX 2 VARIABLES ARE
 320 C       MULTIPLIED BY H, THE INDEX 3 VARIABLES BY H**2.
 321 C
 322 C    IWORK(5)  DIMENSION OF THE INDEX 1 VARIABLES (MUST BE > 0). FOR
 323 C              ODE'S THIS EQUALS THE DIMENSION OF THE SYSTEM.
 324 C              DEFAULT IWORK(5)=N.
 325 C
 326 C    IWORK(6)  DIMENSION OF THE INDEX 2 VARIABLES. DEFAULT IWORK(6)=0.
 327 C
 328 C    IWORK(7)  DIMENSION OF THE INDEX 3 VARIABLES. DEFAULT IWORK(7)=0.
 329 C
 330 C    IWORK(8)  SWITCH FOR STEP SIZE STRATEGY
 331 C              IF IWORK(8).EQ.1  MOD. PREDICTIVE CONTROLLER (GUSTAFSSON)
 332 C              IF IWORK(8).EQ.2  CLASSICAL STEP SIZE CONTROL
 333 C              THE DEFAULT VALUE (FOR IWORK(8)=0) IS IWORK(8)=1.
 334 C              THE CHOICE IWORK(8).EQ.1 SEEMS TO PRODUCE SAFER RESULTS;
 335 C              FOR SIMPLE PROBLEMS, THE CHOICE IWORK(8).EQ.2 PRODUCES
 336 C              OFTEN SLIGHTLY FASTER RUNS
 337 C
 338 C       IF THE DIFFERENTIAL SYSTEM HAS THE SPECIAL STRUCTURE THAT
 339 C            Y(I)' = Y(I+M2)   FOR  I=1,...,M1,
 340 C       WITH M1 A MULTIPLE OF M2, A SUBSTANTIAL GAIN IN COMPUTERTIME
 341 C       CAN BE ACHIEVED BY SETTING THE FOLLOWING TWO PARAMETERS. E.G.,
 342 C       FOR SECOND ORDER SYSTEMS P'=V, V'=G(P,V), WHERE P AND V ARE
 343 C       VECTORS OF DIMENSION N/2, ONE HAS TO PUT M1=M2=N/2.
 344 C       FOR M1>0 SOME OF THE INPUT PARAMETERS HAVE DIFFERENT MEANINGS:
 345 C       - JAC: ONLY THE ELEMENTS OF THE NON-TRIVIAL PART OF THE
 346 C              JACOBIAN HAVE TO BE STORED
 347 C              IF (MLJAC.EQ.N-M1) THE JACOBIAN IS SUPPOSED TO BE FULL
 348 C                 DFY(I,J) = PARTIAL F(I+M1) / PARTIAL Y(J)
 349 C                FOR I=1,N-M1 AND J=1,N.
 350 C              ELSE, THE JACOBIAN IS BANDED ( M1 = M2 * MM )
 351 C                 DFY(I-J+MUJAC+1,J+K*M2) = PARTIAL F(I+M1) / PARTIAL Y(J+K*M2)
 352 C                FOR I=1,MLJAC+MUJAC+1 AND J=1,M2 AND K=0,MM.
 353 C       - MLJAC: MLJAC=N-M1: IF THE NON-TRIVIAL PART OF THE JACOBIAN IS FULL
 354 C                0<=MLJAC<N-M1: IF THE (MM+1) SUBMATRICES (FOR K=0,MM)
 355 C                     PARTIAL F(I+M1) / PARTIAL Y(J+K*M2),  I,J=1,M2
 356 C                    ARE BANDED, MLJAC IS THE MAXIMAL LOWER BANDWIDTH
 357 C                    OF THESE MM+1 SUBMATRICES
 358 C       - MUJAC: MAXIMAL UPPER BANDWIDTH OF THESE MM+1 SUBMATRICES
 359 C                NEED NOT BE DEFINED IF MLJAC=N-M1
 360 C       - MAS: IF IMAS=0 THIS MATRIX IS ASSUMED TO BE THE IDENTITY AND
 361 C              NEED NOT BE DEFINED. SUPPLY A DUMMY SUBROUTINE IN THIS CASE.
 362 C              IT IS ASSUMED THAT ONLY THE ELEMENTS OF RIGHT LOWER BLOCK OF
 363 C              DIMENSION N-M1 DIFFER FROM THAT OF THE IDENTITY MATRIX.
 364 C              IF (MLMAS.EQ.N-M1) THIS SUBMATRIX IS SUPPOSED TO BE FULL
 365 C                 AM(I,J) = M(I+M1,J+M1)     FOR I=1,N-M1 AND J=1,N-M1.
 366 C              ELSE, THE MASS MATRIX IS BANDED
 367 C                 AM(I-J+MUMAS+1,J) = M(I+M1,J+M1)
 368 C       - MLMAS: MLMAS=N-M1: IF THE NON-TRIVIAL PART OF M IS FULL
 369 C                0<=MLMAS<N-M1: LOWER BANDWIDTH OF THE MASS MATRIX
 370 C       - MUMAS: UPPER BANDWIDTH OF THE MASS MATRIX
 371 C                NEED NOT BE DEFINED IF MLMAS=N-M1
 372 C
 373 C    IWORK(9)  THE VALUE OF M1.  DEFAULT M1=0.
 374 C
 375 C    IWORK(10) THE VALUE OF M2.  DEFAULT M2=M1.
 376 C
 377 C ----------
 378 C
 379 C    WORK(1)   UROUND, THE ROUNDING UNIT, DEFAULT 1.D-16.
 380 C
 381 C    WORK(2)   THE SAFETY FACTOR IN STEP SIZE PREDICTION,
 382 C              DEFAULT 0.9D0.
 383 C
 384 C    WORK(3)   DECIDES WHETHER THE JACOBIAN SHOULD BE RECOMPUTED;
 385 C              INCREASE WORK(3), TO 0.1 SAY, WHEN JACOBIAN EVALUATIONS
 386 C              ARE COSTLY. FOR SMALL SYSTEMS WORK(3) SHOULD BE SMALLER
 387 C              (0.001D0, SAY). NEGATIV WORK(3) FORCES THE CODE TO
 388 C              COMPUTE THE JACOBIAN AFTER EVERY ACCEPTED STEP.
 389 C              DEFAULT 0.001D0.
 390 C
 391 C    WORK(4)   STOPPING CRITERION FOR NEWTON'S METHOD, USUALLY CHOSEN <1.
 392 C              SMALLER VALUES OF WORK(4) MAKE THE CODE SLOWER, BUT SAFER.
 393 C              DEFAULT 0.03D0.
 394 C
 395 C    WORK(5) AND WORK(6) : IF WORK(5) < HNEW/HOLD < WORK(6), THEN THE
 396 C              STEP SIZE IS NOT CHANGED. THIS SAVES, TOGETHER WITH A
 397 C              LARGE WORK(3), LU-DECOMPOSITIONS AND COMPUTING TIME FOR
 398 C              LARGE SYSTEMS. FOR SMALL SYSTEMS ONE MAY HAVE
 399 C              WORK(5)=1.D0, WORK(6)=1.2D0, FOR LARGE FULL SYSTEMS
 400 C              WORK(5)=0.99D0, WORK(6)=2.D0 MIGHT BE GOOD.
 401 C              DEFAULTS WORK(5)=1.D0, WORK(6)=1.2D0 .
 402 C
 403 C    WORK(7)   MAXIMAL STEP SIZE, DEFAULT XEND-X.
 404 C
 405 C    WORK(8), WORK(9)   PARAMETERS FOR STEP SIZE SELECTION
 406 C              THE NEW STEP SIZE IS CHOSEN SUBJECT TO THE RESTRICTION
 407 C                 WORK(8) <= HNEW/HOLD <= WORK(9)
 408 C              DEFAULT VALUES: WORK(8)=0.2D0, WORK(9)=8.D0
 409 C
 410 C-----------------------------------------------------------------------
 411 C
 412 C     OUTPUT PARAMETERS
 413 C     -----------------
 414 C     X           X-VALUE FOR WHICH THE SOLUTION HAS BEEN COMPUTED
 415 C                 (AFTER SUCCESSFUL RETURN X=XEND).
 416 C
 417 C     Y(N)        NUMERICAL SOLUTION AT X
 418 C
 419 C     H           PREDICTED STEP SIZE OF THE LAST ACCEPTED STEP
 420 C
 421 C     IDID        REPORTS ON SUCCESSFULNESS UPON RETURN:
 422 C                   IDID= 1  COMPUTATION SUCCESSFUL,
 423 C                   IDID= 2  COMPUT. SUCCESSFUL (INTERRUPTED BY SOLOUT)
 424 C                   IDID=-1  INPUT IS NOT CONSISTENT,
 425 C                   IDID=-2  LARGER NMAX IS NEEDED,
 426 C                   IDID=-3  STEP SIZE BECOMES TOO SMALL,
 427 C                   IDID=-4  MATRIX IS REPEATEDLY SINGULAR.
 428 C
 429 C   IWORK(14)  NFCN    NUMBER OF FUNCTION EVALUATIONS (THOSE FOR NUMERICAL
 430 C                      EVALUATION OF THE JACOBIAN ARE NOT COUNTED)
 431 C   IWORK(15)  NJAC    NUMBER OF JACOBIAN EVALUATIONS (EITHER ANALYTICALLY
 432 C                      OR NUMERICALLY)
 433 C   IWORK(16)  NSTEP   NUMBER OF COMPUTED STEPS
 434 C   IWORK(17)  NACCPT  NUMBER OF ACCEPTED STEPS
 435 C   IWORK(18)  NREJCT  NUMBER OF REJECTED STEPS (DUE TO ERROR TEST),
 436 C                      (STEP REJECTIONS IN THE FIRST STEP ARE NOT COUNTED)
 437 C   IWORK(19)  NDEC    NUMBER OF LU-DECOMPOSITIONS OF BOTH MATRICES
 438 C   IWORK(20)  NSOL    NUMBER OF FORWARD-BACKWARD SUBSTITUTIONS, OF BOTH
 439 C                      SYSTEMS; THE NSTEP FORWARD-BACKWARD SUBSTITUTIONS,
 440 C                      NEEDED FOR STEP SIZE SELECTION, ARE NOT COUNTED
 441 C-----------------------------------------------------------------------
 442 C *** *** *** *** *** *** *** *** *** *** *** *** ***
 443 C          DECLARATIONS
 444 C *** *** *** *** *** *** *** *** *** *** *** *** ***
 445       IMPLICIT KPP_REAL (A-H,O-Z)
 446       DIMENSION Y(N),AbsTol(*),RelTol(*),WORK(LWORK),IWORK(LIWORK)
 447       DIMENSION RPAR(*),IPAR(*)
 448       LOGICAL IMPLCT,JBAND,ARRET,STARTN,PRED
 449       EXTERNAL FCN,JAC,MAS,SOLOUT
 450 C *** *** *** *** *** *** ***
 451 C        SETTING THE PARAMETERS
 452 C *** *** *** *** *** *** ***
 453        NFCN=0
 454        NJAC=0
 455        NSTEP=0
 456        NACCPT=0
 457        NREJCT=0
 458        NDEC=0
 459        NSOL=0
 460        ARRET=.FALSE.
 461 C -------- NMAX , THE MAXIMAL NUMBER OF STEPS -----
 462       IF (IWORK(2).EQ.0) THEN
 463          NMAX=100000
 464       ELSE
 465          NMAX=IWORK(2)
 466          IF (NMAX.LE.0) THEN
 467             WRITE(6,*)' WRONG INPUT IWORK(2)=',IWORK(2)
 468             ARRET=.TRUE.
 469          END IF
 470       END IF
 471 C -------- NIT    MAXIMAL NUMBER OF NEWTON ITERATIONS
 472       IF (IWORK(3).EQ.0) THEN
 473          NIT=7
 474       ELSE
 475          NIT=IWORK(3)
 476          IF (NIT.LE.0) THEN
 477             WRITE(6,*)' CURIOUS INPUT IWORK(3)=',IWORK(3)
 478             ARRET=.TRUE.
 479          END IF
 480       END IF
 481 C -------- STARTN  SWITCH FOR STARTING VALUES OF NEWTON ITERATIONS
 482       IF(IWORK(4).EQ.0)THEN
 483          STARTN=.FALSE.
 484       ELSE
 485          STARTN=.TRUE.
 486       END IF
 487 C -------- PARAMETER FOR DIFFERENTIAL-ALGEBRAIC COMPONENTS
 488       NIND1=IWORK(5)
 489       NIND2=IWORK(6)
 490       NIND3=IWORK(7)
 491       IF (NIND1.EQ.0) NIND1=N
 492       IF (NIND1+NIND2+NIND3.NE.N) THEN
 493        WRITE(6,*)' CURIOUS INPUT FOR IWORK(5,6,7)=',NIND1,NIND2,NIND3
 494        ARRET=.TRUE.
 495       END IF
 496 C -------- PRED   STEP SIZE CONTROL
 497       IF(IWORK(8).LE.1)THEN
 498          PRED=.TRUE.
 499       ELSE
 500          PRED=.FALSE.
 501       END IF
 502 C -------- PARAMETER FOR SECOND ORDER EQUATIONS
 503       M1=IWORK(9)
 504       M2=IWORK(10)
 505       NM1=N-M1
 506       IF (M1.EQ.0) M2=N
 507       IF (M2.EQ.0) M2=M1
 508       IF (M1.LT.0.OR.M2.LT.0.OR.M1+M2.GT.N) THEN
 509        WRITE(6,*)' CURIOUS INPUT FOR IWORK(9,10)=',M1,M2
 510        ARRET=.TRUE.
 511       END IF
 512 C -------- UROUND   SMALLEST NUMBER SATISFYING 1.0D0+UROUND>1.0D0
 513       IF (WORK(1).EQ.0.0D0) THEN
 514          UROUND=1.0D-16
 515       ELSE
 516          UROUND=WORK(1)
 517          IF (UROUND.LE.1.0D-19.OR.UROUND.GE.1.0D0) THEN
 518             WRITE(6,*)' COEFFICIENTS HAVE 20 DIGITS, UROUND=',WORK(1)
 519             ARRET=.TRUE.
 520          END IF
 521       END IF
 522 C --------- SAFE     SAFETY FACTOR IN STEP SIZE PREDICTION
 523       IF (WORK(2).EQ.0.0D0) THEN
 524          SAFE=0.9D0
 525       ELSE
 526          SAFE=WORK(2)
 527          IF (SAFE.LE.0.001D0.OR.SAFE.GE.1.0D0) THEN
 528             WRITE(6,*)' CURIOUS INPUT FOR WORK(2)=',WORK(2)
 529             ARRET=.TRUE.
 530          END IF
 531       END IF
 532 C ------ THET     DECIDES WHETHER THE JACOBIAN SHOULD BE RECOMPUTED;
 533       IF (WORK(3).EQ.0.D0) THEN
 534          THET=0.001D0
 535       ELSE
 536          THET=WORK(3)
 537          IF (THET.LE.0.0D0.OR.THET.GE.1.0D0) THEN
 538             WRITE(6,*)' CURIOUS INPUT FOR WORK(3)=',WORK(3)
 539             ARRET=.TRUE.
 540          END IF
 541       END IF
 542 C --- FNEWT   STOPPING CRIERION FOR NEWTON'S METHOD, USUALLY CHOSEN <1.
 543       IF (WORK(4).EQ.0.D0) THEN
 544          FNEWT=0.03D0
 545       ELSE
 546          FNEWT=WORK(4)
 547          IF (FNEWT.LE.UROUND) THEN
 548             WRITE(6,*)' CURIOUS INPUT FOR WORK(4)=',WORK(4)
 549             ARRET=.TRUE.
 550          END IF
 551       END IF
 552 C --- QUOT1 AND QUOT2: IF QUOT1 < HNEW/HOLD < QUOT2, STEP SIZE = CONST.
 553       IF (WORK(5).EQ.0.D0) THEN
 554          QUOT1=1.D0
 555       ELSE
 556          QUOT1=WORK(5)
 557       END IF
 558       IF (WORK(6).EQ.0.D0) THEN
 559          QUOT2=1.2D0
 560       ELSE
 561          QUOT2=WORK(6)
 562       END IF
 563       IF (QUOT1.GT.1.0D0.OR.QUOT2.LT.1.0D0) THEN
 564          WRITE(6,*)' CURIOUS INPUT FOR WORK(5,6)=',QUOT1,QUOT2
 565          ARRET=.TRUE.
 566       END IF
 567 C -------- MAXIMAL STEP SIZE
 568       IF (WORK(7).EQ.0.D0) THEN
 569          HMAX=XEND-X
 570       ELSE
 571          HMAX=WORK(7)
 572       END IF
 573 C -------  FACL,FACR     PARAMETERS FOR STEP SIZE SELECTION
 574       IF(WORK(8).EQ.0.D0)THEN
 575          FACL=5.D0
 576       ELSE
 577          FACL=1.D0/WORK(8)
 578       END IF
 579       IF(WORK(9).EQ.0.D0)THEN
 580          FACR=1.D0/8.0D0
 581       ELSE
 582          FACR=1.D0/WORK(9)
 583       END IF
 584       IF (FACL.LT.1.0D0.OR.FACR.GT.1.0D0) THEN
 585             WRITE(6,*)' CURIOUS INPUT WORK(8,9)=',WORK(8),WORK(9)
 586             ARRET=.TRUE.
 587          END IF
 588 C --------- CHECK IF TOLERANCES ARE O.K.
 589       IF (ITOL.EQ.0) THEN
 590           IF (AbsTol(1).LE.0.D0.OR.RelTol(1).LE.10.D0*UROUND) THEN
 591               WRITE (6,*) ' TOLERANCES ARE TOO SMALL'
 592               ARRET=.TRUE.
 593           END IF
 594       ELSE
 595           DO I=1,N
 596           IF (AbsTol(I).LE.0.D0.OR.RelTol(I).LE.10.D0*UROUND) THEN
 597               WRITE (6,*) ' TOLERANCES(',I,') ARE TOO SMALL'
 598               ARRET=.TRUE.
 599           END IF
 600           END DO
 601       END IF
 602 C *** *** *** *** *** *** *** *** *** *** *** *** ***
 603 C         COMPUTATION OF ARRAY ENTRIES
 604 C *** *** *** *** *** *** *** *** *** *** *** *** ***
 605 C ---- IMPLICIT, BANDED OR NOT ?
 606       IMPLCT=IMAS.NE.0
 607       JBAND=MLJAC.LT.NM1
 608 C -------- COMPUTATION OF THE ROW-DIMENSIONS OF THE 2-ARRAYS ---
 609 C -- JACOBIAN  AND  MATRICES E1, E2
 610       IF (JBAND) THEN
 611          LDJAC=MLJAC+MUJAC+1
 612          LDE1=MLJAC+LDJAC
 613       ELSE
 614          MLJAC=NM1
 615          MUJAC=NM1
 616          LDJAC=NM1
 617          LDE1=NM1
 618       END IF
 619 C -- MASS MATRIX
 620       IF (IMPLCT) THEN
 621           IF (MLMAS.NE.NM1) THEN
 622               LDMAS=MLMAS+MUMAS+1
 623               IF (JBAND) THEN
 624                  IJOB=4
 625               ELSE
 626                  IJOB=3
 627               END IF
 628           ELSE
 629               LDMAS=NM1
 630               IJOB=5
 631           END IF
 632 C ------ BANDWITH OF "MAS" NOT SMALLER THAN BANDWITH OF "JAC"
 633           IF (MLMAS.GT.MLJAC.OR.MUMAS.GT.MUJAC) THEN
 634              WRITE (6,*) 'BANDWITH OF "MAS" NOT SMALLER THAN BANDWITH OF
 635      & "JAC"'
 636             ARRET=.TRUE.
 637           END IF
 638       ELSE
 639           LDMAS=0
 640           IF (JBAND) THEN
 641              IJOB=2
 642           ELSE
 643              IJOB=1
 644              IF (N.GT.2.AND.IWORK(1).NE.0) IJOB=7
 645           END IF
 646       END IF
 647       LDMAS2=MAX(1,LDMAS)
 648 C ------ HESSENBERG OPTION ONLY FOR EXPLICIT EQU. WITH FULL JACOBIAN
 649       IF ((IMPLCT.OR.JBAND).AND.IJOB.EQ.7) THEN
 650          WRITE(6,*)' HESSENBERG OPTION ONLY FOR EXPLICIT EQUATIONS WITH
 651      &FULL JACOBIAN'
 652          ARRET=.TRUE.
 653       END IF
 654 C ------- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN WORK -----
 655       IEZ1=21
 656       IEZ2=IEZ1+N
 657       IEZ3=IEZ2+N
 658       IEY0=IEZ3+N
 659       IESCAL=IEY0+N
 660       IEF1=IESCAL+N
 661       IEF2=IEF1+N
 662       IEF3=IEF2+N
 663       IECON=IEF3+N
 664       IEJAC=IECON+4*N
 665       IEMAS=IEJAC+N*LDJAC
 666       IEE1=IEMAS+NM1*LDMAS
 667       IEE2R=IEE1+NM1*LDE1
 668       IEE2I=IEE2R+NM1*LDE1
 669 C ------ TOTAL STORAGE REQUIREMENT -----------
 670       ISTORE=IEE2I+NM1*LDE1-1
 671       IF(ISTORE.GT.LWORK)THEN
 672          WRITE(6,*)' INSUFFICIENT STORAGE FOR WORK, MIN. LWORK=',ISTORE
 673          ARRET=.TRUE.
 674       END IF
 675 C ------- ENTRY POINTS FOR INTEGER WORKSPACE -----
 676       IEIP1=21
 677       IEIP2=IEIP1+NM1
 678       IEIPH=IEIP2+NM1
 679 C --------- TOTAL REQUIREMENT ---------------
 680       ISTORE=IEIPH+NM1-1
 681       IF (ISTORE.GT.LIWORK) THEN
 682          WRITE(6,*)' INSUFF. STORAGE FOR IWORK, MIN. LIWORK=',ISTORE
 683          ARRET=.TRUE.
 684       END IF
 685 C ------ WHEN A FAIL HAS OCCURED, WE RETURN WITH IDID=-1
 686       IF (ARRET) THEN
 687          IDID=-1
 688          RETURN
 689       END IF
 690 C -------- CALL TO CORE INTEGRATOR ------------
 691       CALL RADCOR(N,FCN,X,Y,XEND,HMAX,H,RelTol,AbsTol,ITOL,
 692      &   JAC,IJAC,MLJAC,MUJAC,MAS,MLMAS,MUMAS,SOLOUT,IOUT,IDID,
 693      &   NMAX,UROUND,SAFE,THET,FNEWT,QUOT1,QUOT2,NIT,IJOB,STARTN,
 694      &   NIND1,NIND2,NIND3,PRED,FACL,FACR,M1,M2,NM1,
 695      &   IMPLCT,JBAND,LDJAC,LDE1,LDMAS2,WORK(IEZ1),WORK(IEZ2),
 696      &   WORK(IEZ3),WORK(IEY0),WORK(IESCAL),WORK(IEF1),WORK(IEF2),
 697      &   WORK(IEF3),WORK(IEJAC),WORK(IEE1),WORK(IEE2R),WORK(IEE2I),
 698      &   WORK(IEMAS),IWORK(IEIP1),IWORK(IEIP2),IWORK(IEIPH),
 699      &   WORK(IECON),NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL,RPAR,IPAR)
 700       IWORK(14)=NFCN
 701       IWORK(15)=NJAC
 702       IWORK(16)=NSTEP
 703       IWORK(17)=NACCPT
 704       IWORK(18)=NREJCT
 705       IWORK(19)=NDEC
 706       IWORK(20)=NSOL
 707 C ----------- RETURN -----------
 708       RETURN
 709       END
 710 C
 711 C     END OF SUBROUTINE RADAU5
 712 C
 713 C ***********************************************************
 714 C
 715       SUBROUTINE RADCOR(N,FCN,X,Y,XEND,HMAX,H,RelTol,AbsTol,ITOL,
 716      &   JAC,IJAC,MLJAC,MUJAC,MAS,MLMAS,MUMAS,SOLOUT,IOUT,IDID,
 717      &   NMAX,UROUND,SAFE,THET,FNEWT,QUOT1,QUOT2,NIT,IJOB,STARTN,
 718      &   NIND1,NIND2,NIND3,PRED,FACL,FACR,M1,M2,NM1,
 719      &   IMPLCT,BANDED,LDJAC,LDE1,LDMAS,Z1,Z2,Z3,
 720      &   Y0,SCAL,F1,F2,F3,FJAC,E1,E2R,E2I,FMAS,IP1,IP2,IPHES,
 721      &   CONT,NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL,RPAR,IPAR)
 722 C ----------------------------------------------------------
 723 C     CORE INTEGRATOR FOR RADAU5
 724 C     PARAMETERS SAME AS IN RADAU5 WITH WORKSPACE ADDED
 725 C ----------------------------------------------------------
 726 C         DECLARATIONS
 727 C ----------------------------------------------------------
 728       IMPLICIT KPP_REAL (A-H,O-Z)
 729       DIMENSION Y(N),Z1(N),Z2(N),Z3(N),Y0(N),SCAL(N),F1(N),F2(N),F3(N)
 730       DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),CONT(4*N)
 731       DIMENSION E1(LDE1,NM1),E2R(LDE1,NM1),E2I(LDE1,NM1)
 732       DIMENSION AbsTol(*),RelTol(*),RPAR(*),IPAR(*)
 733       INTEGER IP1(NM1),IP2(NM1),IPHES(NM1)
 734       COMMON /CONRA5/NN,NN2,NN3,NN4,XSOL,HSOL,C2M1,C1M1
 735       COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
 736       LOGICAL REJECT,FIRST,IMPLCT,BANDED,CALJAC,STARTN,CALHES
 737       LOGICAL INDEX1,INDEX2,INDEX3,LAST,PRED
 738       EXTERNAL FCN, JAC, MAS, SOLOUT
 739 C *** *** *** *** *** *** ***
 740 C  INITIALISATIONS
 741 C *** *** *** *** *** *** ***
 742 C --------- DUPLIFY N FOR COMMON BLOCK CONT -----
 743       NN=N
 744       NN2=2*N
 745       NN3=3*N
 746       LRC=4*N
 747 C -------- CHECK THE INDEX OF THE PROBLEM -----
 748       INDEX1=NIND1.NE.0
 749       INDEX2=NIND2.NE.0
 750       INDEX3=NIND3.NE.0
 751 C ------- COMPUTE MASS MATRIX FOR IMPLICIT CASE ----------
 752       IF (IMPLCT) CALL MAS(NM1,FMAS,LDMAS,RPAR,IPAR)
 753 C ---------- CONSTANTS ---------
 754       SQ6=DSQRT(6.D0)
 755       C1=(4.D0-SQ6)/10.D0
 756       C2=(4.D0+SQ6)/10.D0
 757       C1M1=C1-1.D0
 758       C2M1=C2-1.D0
 759       C1MC2=C1-C2
 760       DD1=-(13.D0+7.D0*SQ6)/3.D0
 761       DD2=(-13.D0+7.D0*SQ6)/3.D0
 762       DD3=-1.D0/3.D0
 763       U1=(6.D0+81.D0**(1.D0/3.D0)-9.D0**(1.D0/3.D0))/30.D0
 764       ALPH=(12.D0-81.D0**(1.D0/3.D0)+9.D0**(1.D0/3.D0))/60.D0
 765       BETA=(81.D0**(1.D0/3.D0)+9.D0**(1.D0/3.D0))*DSQRT(3.D0)/60.D0
 766       CNO=ALPH**2+BETA**2
 767       U1=1.0D0/U1
 768       ALPH=ALPH/CNO
 769       BETA=BETA/CNO
 770       T11=9.1232394870892942792D-02
 771       T12=-0.14125529502095420843D0
 772       T13=-3.0029194105147424492D-02
 773       T21=0.24171793270710701896D0
 774       T22=0.20412935229379993199D0
 775       T23=0.38294211275726193779D0
 776       T31=0.96604818261509293619D0
 777       TI11=4.3255798900631553510D0
 778       TI12=0.33919925181580986954D0
 779       TI13=0.54177053993587487119D0
 780       TI21=-4.1787185915519047273D0
 781       TI22=-0.32768282076106238708D0
 782       TI23=0.47662355450055045196D0
 783       TI31=-0.50287263494578687595D0
 784       TI32=2.5719269498556054292D0
 785       TI33=-0.59603920482822492497D0
 786       IF (M1.GT.0) IJOB=IJOB+10
 787       POSNEG=SIGN(1.D0,XEND-X)
 788       HMAXN=MIN(ABS(HMAX),ABS(XEND-X))
 789       IF (ABS(H).LE.10.D0*UROUND) H=1.0D-6
 790       H=MIN(ABS(H),HMAXN)
 791       H=SIGN(H,POSNEG)
 792       HOLD=H
 793       REJECT=.FALSE.
 794       FIRST=.TRUE.
 795       LAST=.FALSE.
 796       IF ((X+H*1.0001D0-XEND)*POSNEG.GE.0.D0) THEN
 797          H=XEND-X
 798          LAST=.TRUE.
 799       END IF
 800       FACCON=1.D0
 801       CFAC=SAFE*(1+2*NIT)
 802       NSING=0
 803       XOLD=X
 804       IF (IOUT.NE.0) THEN
 805           IRTRN=1
 806           NRSOL=1
 807           XOSOL=XOLD
 808           XSOL=X
 809           DO I=1,N
 810              CONT(I)=Y(I)
 811           END DO
 812           NSOLU=N
 813           HSOL=HOLD
 814           CALL SOLOUT(NRSOL,XOSOL,XSOL,Y,CONT,LRC,NSOLU,
 815      &                RPAR,IPAR,IRTRN)
 816           IF (IRTRN.LT.0) GOTO 179
 817       END IF
 818       MLE=MLJAC
 819       MUE=MUJAC
 820       MBJAC=MLJAC+MUJAC+1
 821       MBB=MLMAS+MUMAS+1
 822       MDIAG=MLE+MUE+1
 823       MDIFF=MLE+MUE-MUMAS
 824       MBDIAG=MUMAS+1
 825       N2=2*N
 826       N3=3*N
 827       IF (ITOL.EQ.0) THEN
 828           DO I=1,N
 829              SCAL(I)=AbsTol(1)+RelTol(1)*ABS(Y(I))
 830           END DO
 831       ELSE
 832           DO I=1,N
 833              SCAL(I)=AbsTol(I)+RelTol(I)*ABS(Y(I))
 834           END DO
 835       END IF
 836       HHFAC=H
 837       CALL FCN(N,X,Y,Y0,RPAR,IPAR)
 838       NFCN=NFCN+1
 839 C --- BASIC INTEGRATION STEP
 840   10  CONTINUE
 841 C *** *** *** *** *** *** ***
 842 C  COMPUTATION OF THE JACOBIAN
 843 C *** *** *** *** *** *** ***
 844       NJAC=NJAC+1
 845       IF (IJAC.EQ.0) THEN
 846 C --- COMPUTE JACOBIAN MATRIX NUMERICALLY
 847          IF (BANDED) THEN
 848 C --- JACOBIAN IS BANDED
 849             MUJACP=MUJAC+1
 850             MD=MIN(MBJAC,M2)
 851             DO MM=1,M1/M2+1
 852                DO K=1,MD
 853                   J=K+(MM-1)*M2
 854  12               F1(J)=Y(J)
 855                   F2(J)=DSQRT(UROUND*MAX(1.D-5,ABS(Y(J))))
 856                   Y(J)=Y(J)+F2(J)
 857                   J=J+MD
 858                   IF (J.LE.MM*M2) GOTO 12
 859                   CALL FCN(N,X,Y,CONT,RPAR,IPAR)
 860                   J=K+(MM-1)*M2
 861                   J1=K
 862                   LBEG=MAX(1,J1-MUJAC)+M1
 863  14               LEND=MIN(M2,J1+MLJAC)+M1
 864                   Y(J)=F1(J)
 865                   MUJACJ=MUJACP-J1-M1
 866                   DO L=LBEG,LEND
 867                      FJAC(L+MUJACJ,J)=(CONT(L)-Y0(L))/F2(J)
 868                   END DO
 869                   J=J+MD
 870                   J1=J1+MD
 871                   LBEG=LEND+1
 872                   IF (J.LE.MM*M2) GOTO 14
 873                END DO
 874             END DO
 875          ELSE
 876 C --- JACOBIAN IS FULL
 877             DO I=1,N
 878                YSAFE=Y(I)
 879                DELT=DSQRT(UROUND*MAX(1.D-5,ABS(YSAFE)))
 880                Y(I)=YSAFE+DELT
 881                CALL FCN(N,X,Y,CONT,RPAR,IPAR)
 882                DO J=M1+1,N
 883                  FJAC(J-M1,I)=(CONT(J)-Y0(J))/DELT
 884                END DO
 885                Y(I)=YSAFE
 886             END DO
 887          END IF
 888       ELSE
 889 C --- COMPUTE JACOBIAN MATRIX ANALYTICALLY
 890       CALL JAC_CHEM(N,X,Y,FJAC,LDJAC,RPAR,IPAR)
 891       END IF
 892       CALJAC=.TRUE.
 893       CALHES=.TRUE.
 894   20  CONTINUE
 895 C --- COMPUTE THE MATRICES E1 AND E2 AND THEIR DECOMPOSITIONS
 896       FAC1=U1/H
 897       ALPHN=ALPH/H
 898       BETAN=BETA/H
 899       CALL DECOMR(N,FJAC,LDJAC,FMAS,LDMAS,MLMAS,MUMAS,
 900      &            M1,M2,NM1,FAC1,E1,LDE1,IP1,IER,IJOB,CALHES,IPHES)
 901       IF (IER.NE.0) GOTO 78
 902       CALL DECOMC(N,FJAC,LDJAC,FMAS,LDMAS,MLMAS,MUMAS,
 903      &            M1,M2,NM1,ALPHN,BETAN,E2R,E2I,LDE1,IP2,IER,IJOB)
 904       IF (IER.NE.0) GOTO 78
 905       NDEC=NDEC+1
 906   30  CONTINUE
 907       NSTEP=NSTEP+1
 908       IF (NSTEP.GT.NMAX) GOTO 178
 909       IF (0.1D0*ABS(H).LE.ABS(X)*UROUND) GOTO 177
 910           IF (INDEX2) THEN
 911              DO I=NIND1+1,NIND1+NIND2
 912                 SCAL(I)=SCAL(I)/HHFAC
 913              END DO
 914           END IF
 915           IF (INDEX3) THEN
 916              DO I=NIND1+NIND2+1,NIND1+NIND2+NIND3
 917                 SCAL(I)=SCAL(I)/(HHFAC*HHFAC)
 918              END DO
 919           END IF
 920       XPH=X+H
 921 C *** *** *** *** *** *** ***
 922 C  STARTING VALUES FOR NEWTON ITERATION
 923 C *** *** *** *** *** *** ***
 924       IF (FIRST.OR.STARTN) THEN
 925          DO I=1,N
 926             Z1(I)=0.D0
 927             Z2(I)=0.D0
 928             Z3(I)=0.D0
 929             F1(I)=0.D0
 930             F2(I)=0.D0
 931             F3(I)=0.D0
 932          END DO
 933       ELSE
 934          C3Q=H/HOLD
 935          C1Q=C1*C3Q
 936          C2Q=C2*C3Q
 937          DO I=1,N
 938             AK1=CONT(I+N)
 939             AK2=CONT(I+N2)
 940             AK3=CONT(I+N3)
 941             Z1I=C1Q*(AK1+(C1Q-C2M1)*(AK2+(C1Q-C1M1)*AK3))
 942             Z2I=C2Q*(AK1+(C2Q-C2M1)*(AK2+(C2Q-C1M1)*AK3))
 943             Z3I=C3Q*(AK1+(C3Q-C2M1)*(AK2+(C3Q-C1M1)*AK3))
 944             Z1(I)=Z1I
 945             Z2(I)=Z2I
 946             Z3(I)=Z3I
 947             F1(I)=TI11*Z1I+TI12*Z2I+TI13*Z3I
 948             F2(I)=TI21*Z1I+TI22*Z2I+TI23*Z3I
 949             F3(I)=TI31*Z1I+TI32*Z2I+TI33*Z3I
 950          END DO
 951       END IF
 952 C *** *** *** *** *** *** ***
 953 C  LOOP FOR THE SIMPLIFIED NEWTON ITERATION
 954 C *** *** *** *** *** *** ***
 955             NEWT=0
 956             FACCON=MAX(FACCON,UROUND)**0.8D0
 957             THETA=ABS(THET)
 958   40        CONTINUE
 959             IF (NEWT.GE.NIT) GOTO 78
 960 C ---     COMPUTE THE RIGHT-HAND SIDE
 961             DO I=1,N
 962                CONT(I)=Y(I)+Z1(I)
 963             END DO
 964             CALL FCN(N,X+C1*H,CONT,Z1,RPAR,IPAR)
 965             DO I=1,N
 966                CONT(I)=Y(I)+Z2(I)
 967             END DO
 968             CALL FCN(N,X+C2*H,CONT,Z2,RPAR,IPAR)
 969             DO I=1,N
 970                CONT(I)=Y(I)+Z3(I)
 971             END DO
 972             CALL FCN(N,XPH,CONT,Z3,RPAR,IPAR)
 973             NFCN=NFCN+3
 974 C ---     KppSolve THE LINEAR SYSTEMS
 975            DO I=1,N
 976               A1=Z1(I)
 977               A2=Z2(I)
 978               A3=Z3(I)
 979               Z1(I)=TI11*A1+TI12*A2+TI13*A3
 980               Z2(I)=TI21*A1+TI22*A2+TI23*A3
 981               Z3(I)=TI31*A1+TI32*A2+TI33*A3
 982            END DO
 983         CALL SLVRAD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
 984      &          M1,M2,NM1,FAC1,ALPHN,BETAN,E1,E2R,E2I,LDE1,Z1,Z2,Z3,
 985      &          F1,F2,F3,CONT,IP1,IP2,IPHES,IER,IJOB)
 986             NSOL=NSOL+1
 987             NEWT=NEWT+1
 988             DYNO=0.D0
 989             DO I=1,N
 990                DENOM=SCAL(I)
 991                DYNO=DYNO+(Z1(I)/DENOM)**2+(Z2(I)/DENOM)**2
 992      &          +(Z3(I)/DENOM)**2
 993             END DO
 994             DYNO=DSQRT(DYNO/N3)
 995 C ---     BAD CONVERGENCE OR NUMBER OF ITERATIONS TO LARGE
 996             IF (NEWT.GT.1.AND.NEWT.LT.NIT) THEN
 997                 THQ=DYNO/DYNOLD
 998                 IF (NEWT.EQ.2) THEN
 999                    THETA=THQ
1000                 ELSE
1001                    THETA=SQRT(THQ*THQOLD)
1002                 END IF
1003                 THQOLD=THQ
1004                 IF (THETA.LT.0.99D0) THEN
1005                     FACCON=THETA/(1.0D0-THETA)
1006                     DYTH=FACCON*DYNO*THETA**(NIT-1-NEWT)/FNEWT
1007                     IF (DYTH.GE.1.0D0) THEN
1008                          QNEWT=DMAX1(1.0D-4,DMIN1(20.0D0,DYTH))
1009                          HHFAC=.8D0*QNEWT**(-1.0D0/(4.0D0+NIT-1-NEWT))
1010                          H=HHFAC*H
1011                          REJECT=.TRUE.
1012                          LAST=.FALSE.
1013                          IF (CALJAC) GOTO 20
1014                          GOTO 10
1015                     END IF
1016                 ELSE
1017                     GOTO 78
1018                 END IF
1019             END IF
1020             DYNOLD=MAX(DYNO,UROUND)
1021             DO I=1,N
1022                F1I=F1(I)+Z1(I)
1023                F2I=F2(I)+Z2(I)
1024                F3I=F3(I)+Z3(I)
1025                F1(I)=F1I
1026                F2(I)=F2I
1027                F3(I)=F3I
1028                Z1(I)=T11*F1I+T12*F2I+T13*F3I
1029                Z2(I)=T21*F1I+T22*F2I+T23*F3I
1030                Z3(I)=T31*F1I+    F2I
1031             END DO
1032             IF (FACCON*DYNO.GT.FNEWT) GOTO 40
1033 C --- ERROR ESTIMATION
1034       CALL ESTRAD (N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
1035      &          H,DD1,DD2,DD3,FCN,NFCN,Y0,Y,IJOB,X,M1,M2,NM1,
1036      &          E1,LDE1,Z1,Z2,Z3,CONT,F1,F2,IP1,IPHES,SCAL,ERR,
1037      &          FIRST,REJECT,FAC1,RPAR,IPAR)
1038 C --- COMPUTATION OF HNEW
1039 C --- WE REQUIRE .2<=HNEW/H<=8.
1040       FAC=MIN(SAFE,CFAC/(NEWT+2*NIT))
1041       QUOT=MAX(FACR,MIN(FACL,ERR**.25D0/FAC))
1042       HNEW=H/QUOT
1043 C *** *** *** *** *** *** ***
1044 C  IS THE ERROR SMALL ENOUGH ?
1045 C *** *** *** *** *** *** ***
1046       IF (ERR.LT.1.D0) THEN
1047 C --- STEP IS ACCEPTED
1048          FIRST=.FALSE.
1049          NACCPT=NACCPT+1
1050          IF (PRED) THEN
1051 C       --- PREDICTIVE CONTROLLER OF GUSTAFSSON
1052             IF (NACCPT.GT.1) THEN
1053                FACGUS=(HACC/H)*(ERR**2/ERRACC)**0.25D0/SAFE
1054                FACGUS=MAX(FACR,MIN(FACL,FACGUS))
1055                QUOT=MAX(QUOT,FACGUS)
1056                HNEW=H/QUOT
1057             END IF
1058             HACC=H
1059             ERRACC=MAX(1.0D-2,ERR)
1060          END IF
1061          XOLD=X
1062          HOLD=H
1063          X=XPH
1064          DO I=1,N
1065             Y(I)=Y(I)+Z3(I)
1066             Z2I=Z2(I)
1067             Z1I=Z1(I)
1068             CONT(I+N)=(Z2I-Z3(I))/C2M1
1069             AK=(Z1I-Z2I)/C1MC2
1070             ACONT3=Z1I/C1
1071             ACONT3=(AK-ACONT3)/C2
1072             CONT(I+N2)=(AK-CONT(I+N))/C1M1
1073             CONT(I+N3)=CONT(I+N2)-ACONT3
1074          END DO
1075          IF (ITOL.EQ.0) THEN
1076              DO I=1,N
1077                 SCAL(I)=AbsTol(1)+RelTol(1)*ABS(Y(I))
1078              END DO
1079          ELSE
1080              DO I=1,N
1081                 SCAL(I)=AbsTol(I)+RelTol(I)*ABS(Y(I))
1082              END DO
1083          END IF
1084          IF (IOUT.NE.0) THEN
1085              NRSOL=NACCPT+1
1086              XSOL=X
1087              XOSOL=XOLD
1088              DO I=1,N
1089                 CONT(I)=Y(I)
1090              END DO
1091              NSOLU=N
1092              HSOL=HOLD
1093              CALL SOLOUT(NRSOL,XOSOL,XSOL,Y,CONT,LRC,NSOLU,
1094      &                   RPAR,IPAR,IRTRN)
1095              IF (IRTRN.LT.0) GOTO 179
1096          END IF
1097          CALJAC=.FALSE.
1098          IF (LAST) THEN
1099             H=HOPT
1100             IDID=1
1101             RETURN
1102          END IF
1103          CALL FCN(N,X,Y,Y0,RPAR,IPAR)
1104          NFCN=NFCN+1
1105          HNEW=POSNEG*MIN(ABS(HNEW),HMAXN)
1106          HOPT=HNEW
1107          HOPT=MIN(H,HNEW)
1108          IF (REJECT) HNEW=POSNEG*MIN(ABS(HNEW),ABS(H))
1109          REJECT=.FALSE.
1110          IF ((X+HNEW/QUOT1-XEND)*POSNEG.GE.0.D0) THEN
1111             H=XEND-X
1112             LAST=.TRUE.
1113          ELSE
1114             QT=HNEW/H
1115             HHFAC=H
1116             IF (THETA.LE.THET.AND.QT.GE.QUOT1.AND.QT.LE.QUOT2) GOTO 30
1117             H=HNEW
1118          END IF
1119          HHFAC=H
1120          IF (THETA.LE.THET) GOTO 20
1121          GOTO 10
1122       ELSE
1123 C --- STEP IS REJECTED
1124          REJECT=.TRUE.
1125          LAST=.FALSE.
1126          IF (FIRST) THEN
1127              H=H*0.1D0
1128              HHFAC=0.1D0
1129          ELSE
1130              HHFAC=HNEW/H
1131              H=HNEW
1132          END IF
1133          IF (NACCPT.GE.1) NREJCT=NREJCT+1
1134          IF (CALJAC) GOTO 20
1135          GOTO 10
1136       END IF
1137 C --- UNEXPECTED STEP-REJECTION
1138   78  CONTINUE
1139       IF (IER.NE.0) THEN
1140           NSING=NSING+1
1141           IF (NSING.GE.5) GOTO 176
1142       END IF
1143       H=H*0.5D0
1144       HHFAC=0.5D0
1145       REJECT=.TRUE.
1146       LAST=.FALSE.
1147       IF (CALJAC) GOTO 20
1148       GOTO 10
1149 C --- FAIL EXIT
1150  176  CONTINUE
1151       WRITE(6,979)X
1152       WRITE(6,*) ' MATRIX IS REPEATEDLY SINGULAR, IER=',IER
1153       IDID=-4
1154       RETURN
1155  177  CONTINUE
1156       WRITE(6,979)X
1157       WRITE(6,*) ' STEP SIZE T0O SMALL, H=',H
1158       IDID=-3
1159       RETURN
1160  178  CONTINUE
1161       WRITE(6,979)X
1162       WRITE(6,*) ' MORE THAN NMAX =',NMAX,'STEPS ARE NEEDED'
1163       IDID=-2
1164       RETURN
1165 C --- EXIT CAUSED BY SOLOUT
1166  179  CONTINUE
1167       WRITE(6,979)X
1168  979  FORMAT(' EXIT OF RADAU5 AT X=',E18.4)
1169       IDID=2
1170       RETURN
1171       END
1172 C
1173 C     END OF SUBROUTINE RADCOR
1174 C
1175 C ***********************************************************
1176 C
1177       KPP_REAL FUNCTION CONTR5(I,X,CONT,LRC)
1178 C ----------------------------------------------------------
1179 C     THIS FUNCTION CAN BE USED FOR CONINUOUS OUTPUT. IT PROVIDES AN
1180 C     APPROXIMATION TO THE I-TH COMPONENT OF THE SOLUTION AT X.
1181 C     IT GIVES THE VALUE OF THE COLLOCATION POLYNOMIAL, DEFINED FOR
1182 C     THE LAST SUCCESSFULLY COMPUTED STEP (BY RADAU5).
1183 C ----------------------------------------------------------
1184       IMPLICIT KPP_REAL (A-H,O-Z)
1185       DIMENSION CONT(LRC)
1186       COMMON /CONRA5/NN,NN2,NN3,NN4,XSOL,HSOL,C2M1,C1M1
1187       S=(X-XSOL)/HSOL
1188       CONTR5=CONT(I)+S*(CONT(I+NN)+(S-C2M1)*(CONT(I+NN2)
1189      &     +(S-C1M1)*CONT(I+NN3)))
1190       RETURN
1191       END
1192 C
1193 C     END OF FUNCTION CONTR5
1194 C
1195 C ***********************************************************
1196 C ******************************************
1197 C     VERSION OF SEPTEMBER 18, 1995
1198 C ******************************************
1199 C
1200       SUBROUTINE DECOMR(N,FJAC,LDJAC,FMAS,LDMAS,MLMAS,MUMAS,
1201      &            M1,M2,NM1,FAC1,E1,LDE1,IP1,IER,IJOB,CALHES,IPHES)
1202       IMPLICIT KPP_REAL (A-H,O-Z)
1203       DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),E1(LDE1,NM1),
1204      &          IP1(NM1),IPHES(N)
1205       LOGICAL CALHES
1206       COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
1207 C
1208       GOTO (1,2,3,4,5,6,7,55,55,55,11,12,13,14,15), IJOB
1209 C
1210 C -----------------------------------------------------------
1211 C
1212    1  CONTINUE
1213 C ---  B=IDENTITY, JACOBIAN A FULL MATRIX
1214       DO J=1,N
1215          DO  I=1,N
1216             E1(I,J)=-FJAC(I,J)
1217          END DO
1218          E1(J,J)=E1(J,J)+FAC1
1219       END DO
1220       CALL DEC (N,LDE1,E1,IP1,IER)
1221       RETURN
1222 C
1223 C -----------------------------------------------------------
1224 C
1225   11  CONTINUE
1226 C ---  B=IDENTITY, JACOBIAN A FULL MATRIX, SECOND ORDER
1227       DO J=1,NM1
1228          JM1=J+M1
1229          DO I=1,NM1
1230             E1(I,J)=-FJAC(I,JM1)
1231          END DO
1232          E1(J,J)=E1(J,J)+FAC1
1233       END DO
1234  45   MM=M1/M2
1235       DO J=1,M2
1236          DO I=1,NM1
1237             SUM=0.D0
1238             DO K=0,MM-1
1239                SUM=(SUM+FJAC(I,J+K*M2))/FAC1
1240             END DO
1241             E1(I,J)=E1(I,J)-SUM
1242          END DO
1243       END DO
1244       CALL DEC (NM1,LDE1,E1,IP1,IER)
1245       RETURN
1246 C
1247 C -----------------------------------------------------------
1248 C
1249    2  CONTINUE
1250 C ---  B=IDENTITY, JACOBIAN A BANDED MATRIX
1251       DO J=1,N
1252          DO I=1,MBJAC
1253             E1(I+MLE,J)=-FJAC(I,J)
1254          END DO
1255          E1(MDIAG,J)=E1(MDIAG,J)+FAC1
1256       END DO
1257       CALL DECB (N,LDE1,E1,MLE,MUE,IP1,IER)
1258       RETURN
1259 C
1260 C -----------------------------------------------------------
1261 C
1262   12  CONTINUE
1263 C ---  B=IDENTITY, JACOBIAN A BANDED MATRIX, SECOND ORDER
1264       DO J=1,NM1
1265          JM1=J+M1
1266          DO I=1,MBJAC
1267             E1(I+MLE,J)=-FJAC(I,JM1)
1268          END DO
1269          E1(MDIAG,J)=E1(MDIAG,J)+FAC1
1270       END DO
1271   46  MM=M1/M2
1272       DO J=1,M2
1273          DO I=1,MBJAC
1274             SUM=0.D0
1275             DO K=0,MM-1
1276                SUM=(SUM+FJAC(I,J+K*M2))/FAC1
1277             END DO
1278             E1(I+MLE,J)=E1(I+MLE,J)-SUM
1279          END DO
1280       END DO
1281       CALL DECB (NM1,LDE1,E1,MLE,MUE,IP1,IER)
1282       RETURN
1283 C
1284 C -----------------------------------------------------------
1285 C
1286    3  CONTINUE
1287 C ---  B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX
1288       DO J=1,N
1289          DO I=1,N
1290             E1(I,J)=-FJAC(I,J)
1291          END DO
1292          DO I=MAX(1,J-MUMAS),MIN(N,J+MLMAS)
1293             E1(I,J)=E1(I,J)+FAC1*FMAS(I-J+MBDIAG,J)
1294          END DO
1295       END DO
1296       CALL DEC (N,LDE1,E1,IP1,IER)
1297       RETURN
1298 C
1299 C -----------------------------------------------------------
1300 C
1301   13  CONTINUE
1302 C ---  B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER
1303       DO J=1,NM1
1304          JM1=J+M1
1305          DO I=1,NM1
1306             E1(I,J)=-FJAC(I,JM1)
1307          END DO
1308          DO I=MAX(1,J-MUMAS),MIN(NM1,J+MLMAS)
1309             E1(I,J)=E1(I,J)+FAC1*FMAS(I-J+MBDIAG,J)
1310          END DO
1311       END DO
1312       GOTO 45
1313 C
1314 C -----------------------------------------------------------
1315 C
1316    4  CONTINUE
1317 C ---  B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX
1318       DO J=1,N
1319          DO I=1,MBJAC
1320             E1(I+MLE,J)=-FJAC(I,J)
1321          END DO
1322          DO I=1,MBB
1323             IB=I+MDIFF
1324             E1(IB,J)=E1(IB,J)+FAC1*FMAS(I,J)
1325          END DO
1326       END DO
1327       CALL DECB (N,LDE1,E1,MLE,MUE,IP1,IER)
1328       RETURN
1329 C
1330 C -----------------------------------------------------------
1331 C
1332   14  CONTINUE
1333 C ---  B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX, SECOND ORDER
1334       DO J=1,NM1
1335          JM1=J+M1
1336          DO I=1,MBJAC
1337             E1(I+MLE,J)=-FJAC(I,JM1)
1338          END DO
1339          DO I=1,MBB
1340             IB=I+MDIFF
1341             E1(IB,J)=E1(IB,J)+FAC1*FMAS(I,J)
1342          END DO
1343       END DO
1344       GOTO 46
1345 C
1346 C -----------------------------------------------------------
1347 C
1348    5  CONTINUE
1349 C ---  B IS A FULL MATRIX, JACOBIAN A FULL MATRIX
1350       DO J=1,N
1351          DO I=1,N
1352             E1(I,J)=FMAS(I,J)*FAC1-FJAC(I,J)
1353          END DO
1354       END DO
1355       CALL DEC (N,LDE1,E1,IP1,IER)
1356       RETURN
1357 C
1358 C -----------------------------------------------------------
1359 C
1360   15  CONTINUE
1361 C ---  B IS A FULL MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER
1362       DO J=1,NM1
1363          JM1=J+M1
1364          DO I=1,NM1
1365             E1(I,J)=FMAS(I,J)*FAC1-FJAC(I,JM1)
1366          END DO
1367       END DO
1368       GOTO 45
1369 C
1370 C -----------------------------------------------------------
1371 C
1372    6  CONTINUE
1373 C ---  B IS A FULL MATRIX, JACOBIAN A BANDED MATRIX
1374 C ---  THIS OPTION IS NOT PROVIDED
1375       RETURN
1376 C
1377 C -----------------------------------------------------------
1378 C
1379    7  CONTINUE
1380 C ---  B=IDENTITY, JACOBIAN A FULL MATRIX, HESSENBERG-OPTION
1381       IF (CALHES) CALL Elmhes (LDJAC,N,1,N,FJAC,IPHES)
1382       CALHES=.FALSE.
1383       DO J=1,N-1
1384          J1=J+1
1385          E1(J1,J)=-FJAC(J1,J)
1386       END DO
1387       DO J=1,N
1388          DO I=1,J
1389             E1(I,J)=-FJAC(I,J)
1390          END DO
1391          E1(J,J)=E1(J,J)+FAC1
1392       END DO
1393       CALL DECH(N,LDE1,E1,1,IP1,IER)
1394       RETURN
1395 C
1396 C -----------------------------------------------------------
1397 C
1398   55  CONTINUE
1399       RETURN
1400       END
1401 C
1402 C     END OF SUBROUTINE DECOMR
1403 C
1404 C ***********************************************************
1405 C
1406       SUBROUTINE DECOMC(N,FJAC,LDJAC,FMAS,LDMAS,MLMAS,MUMAS,
1407      &            M1,M2,NM1,ALPHN,BETAN,E2R,E2I,LDE1,IP2,IER,IJOB)
1408       IMPLICIT KPP_REAL (A-H,O-Z)
1409       DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),
1410      &          E2R(LDE1,NM1),E2I(LDE1,NM1),IP2(NM1)
1411       COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
1412 C
1413       GOTO (1,2,3,4,5,6,7,55,55,55,11,12,13,14,15), IJOB
1414 C
1415 C -----------------------------------------------------------
1416 C
1417    1  CONTINUE
1418 C ---  B=IDENTITY, JACOBIAN A FULL MATRIX
1419       DO J=1,N
1420          DO I=1,N
1421             E2R(I,J)=-FJAC(I,J)
1422             E2I(I,J)=0.D0
1423          END DO
1424          E2R(J,J)=E2R(J,J)+ALPHN
1425          E2I(J,J)=BETAN
1426       END DO
1427       CALL DECC (N,LDE1,E2R,E2I,IP2,IER)
1428       RETURN
1429 C
1430 C -----------------------------------------------------------
1431 C
1432   11  CONTINUE
1433 C ---  B=IDENTITY, JACOBIAN A FULL MATRIX, SECOND ORDER
1434       DO J=1,NM1
1435          JM1=J+M1
1436          DO I=1,NM1
1437             E2R(I,J)=-FJAC(I,JM1)
1438             E2I(I,J)=0.D0
1439          END DO
1440          E2R(J,J)=E2R(J,J)+ALPHN
1441          E2I(J,J)=BETAN
1442       END DO
1443   45  MM=M1/M2
1444       ABNO=ALPHN**2+BETAN**2
1445       ALP=ALPHN/ABNO
1446       BET=BETAN/ABNO
1447       DO J=1,M2
1448          DO I=1,NM1
1449             SUMR=0.D0
1450             SUMI=0.D0
1451             DO K=0,MM-1
1452                SUMS=SUMR+FJAC(I,J+K*M2)
1453                SUMR=SUMS*ALP+SUMI*BET
1454                SUMI=SUMI*ALP-SUMS*BET
1455             END DO
1456             E2R(I,J)=E2R(I,J)-SUMR
1457             E2I(I,J)=E2I(I,J)-SUMI
1458          END DO
1459       END DO
1460       CALL DECC (NM1,LDE1,E2R,E2I,IP2,IER)
1461       RETURN
1462 C
1463 C -----------------------------------------------------------
1464 C
1465    2  CONTINUE
1466 C ---  B=IDENTITY, JACOBIAN A BANDED MATRIX
1467       DO J=1,N
1468          DO I=1,MBJAC
1469             IMLE=I+MLE
1470             E2R(IMLE,J)=-FJAC(I,J)
1471             E2I(IMLE,J)=0.D0
1472          END DO
1473          E2R(MDIAG,J)=E2R(MDIAG,J)+ALPHN
1474          E2I(MDIAG,J)=BETAN
1475       END DO
1476       CALL DECBC (N,LDE1,E2R,E2I,MLE,MUE,IP2,IER)
1477       RETURN
1478 C
1479 C -----------------------------------------------------------
1480 C
1481   12  CONTINUE
1482 C ---  B=IDENTITY, JACOBIAN A BANDED MATRIX, SECOND ORDER
1483       DO J=1,NM1
1484          JM1=J+M1
1485          DO I=1,MBJAC
1486             E2R(I+MLE,J)=-FJAC(I,JM1)
1487             E2I(I+MLE,J)=0.D0
1488          END DO
1489          E2R(MDIAG,J)=E2R(MDIAG,J)+ALPHN
1490          E2I(MDIAG,J)=E2I(MDIAG,J)+BETAN
1491       END DO
1492   46  MM=M1/M2
1493       ABNO=ALPHN**2+BETAN**2
1494       ALP=ALPHN/ABNO
1495       BET=BETAN/ABNO
1496       DO J=1,M2
1497          DO I=1,MBJAC
1498             SUMR=0.D0
1499             SUMI=0.D0
1500             DO K=0,MM-1
1501                SUMS=SUMR+FJAC(I,J+K*M2)
1502                SUMR=SUMS*ALP+SUMI*BET
1503                SUMI=SUMI*ALP-SUMS*BET
1504             END DO
1505             IMLE=I+MLE
1506             E2R(IMLE,J)=E2R(IMLE,J)-SUMR
1507             E2I(IMLE,J)=E2I(IMLE,J)-SUMI
1508          END DO
1509       END DO
1510       CALL DECBC (NM1,LDE1,E2R,E2I,MLE,MUE,IP2,IER)
1511       RETURN
1512 C
1513 C -----------------------------------------------------------
1514 C
1515    3  CONTINUE
1516 C ---  B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX
1517       DO  J=1,N
1518          DO  I=1,N
1519             E2R(I,J)=-FJAC(I,J)
1520             E2I(I,J)=0.D0
1521          END DO
1522       END DO
1523       DO J=1,N
1524          DO I=MAX(1,J-MUMAS),MIN(N,J+MLMAS)
1525             BB=FMAS(I-J+MBDIAG,J)
1526             E2R(I,J)=E2R(I,J)+ALPHN*BB
1527             E2I(I,J)=BETAN*BB
1528          END DO
1529       END DO
1530       CALL DECC(N,LDE1,E2R,E2I,IP2,IER)
1531       RETURN
1532 C
1533 C -----------------------------------------------------------
1534 C
1535   13  CONTINUE
1536 C ---  B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER
1537       DO J=1,NM1
1538          JM1=J+M1
1539          DO I=1,NM1
1540             E2R(I,J)=-FJAC(I,JM1)
1541             E2I(I,J)=0.D0
1542          END DO
1543          DO I=MAX(1,J-MUMAS),MIN(NM1,J+MLMAS)
1544             FFMA=FMAS(I-J+MBDIAG,J)
1545             E2R(I,J)=E2R(I,J)+ALPHN*FFMA
1546             E2I(I,J)=E2I(I,J)+BETAN*FFMA
1547          END DO
1548       END DO
1549       GOTO 45
1550 C
1551 C -----------------------------------------------------------
1552 C
1553    4  CONTINUE
1554 C ---  B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX
1555       DO J=1,N
1556          DO I=1,MBJAC
1557             IMLE=I+MLE
1558             E2R(IMLE,J)=-FJAC(I,J)
1559             E2I(IMLE,J)=0.D0
1560          END DO
1561          DO I=MAX(1,MUMAS+2-J),MIN(MBB,MUMAS+1-J+N)
1562             IB=I+MDIFF
1563             BB=FMAS(I,J)
1564             E2R(IB,J)=E2R(IB,J)+ALPHN*BB
1565             E2I(IB,J)=BETAN*BB
1566          END DO
1567       END DO
1568       CALL DECBC (N,LDE1,E2R,E2I,MLE,MUE,IP2,IER)
1569       RETURN
1570 C
1571 C -----------------------------------------------------------
1572 C
1573   14  CONTINUE
1574 C ---  B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX, SECOND ORDER
1575       DO J=1,NM1
1576          JM1=J+M1
1577          DO I=1,MBJAC
1578             E2R(I+MLE,J)=-FJAC(I,JM1)
1579             E2I(I+MLE,J)=0.D0
1580          END DO
1581          DO I=1,MBB
1582             IB=I+MDIFF
1583             FFMA=FMAS(I,J)
1584             E2R(IB,J)=E2R(IB,J)+ALPHN*FFMA
1585             E2I(IB,J)=E2I(IB,J)+BETAN*FFMA
1586          END DO
1587       END DO
1588       GOTO 46
1589 C
1590 C -----------------------------------------------------------
1591 C
1592    5  CONTINUE
1593 C ---  B IS A FULL MATRIX, JACOBIAN A FULL MATRIX
1594       DO J=1,N
1595          DO I=1,N
1596             BB=FMAS(I,J)
1597             E2R(I,J)=BB*ALPHN-FJAC(I,J)
1598             E2I(I,J)=BB*BETAN
1599          END DO
1600       END DO
1601       CALL DECC(N,LDE1,E2R,E2I,IP2,IER)
1602       RETURN
1603 C
1604 C -----------------------------------------------------------
1605 C
1606   15  CONTINUE
1607 C ---  B IS A FULL MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER
1608       DO J=1,NM1
1609          JM1=J+M1
1610          DO I=1,NM1
1611             E2R(I,J)=ALPHN*FMAS(I,J)-FJAC(I,JM1)
1612             E2I(I,J)=BETAN*FMAS(I,J)
1613          END DO
1614       END DO
1615       GOTO 45
1616 C
1617 C -----------------------------------------------------------
1618 C
1619    6  CONTINUE
1620 C ---  B IS A FULL MATRIX, JACOBIAN A BANDED MATRIX
1621 C ---  THIS OPTION IS NOT PROVIDED
1622       RETURN
1623 C
1624 C -----------------------------------------------------------
1625 C
1626    7  CONTINUE
1627 C ---  B=IDENTITY, JACOBIAN A FULL MATRIX, HESSENBERG-OPTION
1628       DO J=1,N-1
1629          J1=J+1
1630          E2R(J1,J)=-FJAC(J1,J)
1631          E2I(J1,J)=0.D0
1632       END DO
1633       DO J=1,N
1634          DO I=1,J
1635             E2I(I,J)=0.D0
1636             E2R(I,J)=-FJAC(I,J)
1637          END DO
1638          E2R(J,J)=E2R(J,J)+ALPHN
1639          E2I(J,J)=BETAN
1640       END DO
1641       CALL DECHC(N,LDE1,E2R,E2I,1,IP2,IER)
1642       RETURN
1643 C
1644 C -----------------------------------------------------------
1645 C
1646   55  CONTINUE
1647       RETURN
1648       END
1649 C
1650 C     END OF SUBROUTINE DECOMC
1651 C
1652 C ***********************************************************
1653 C
1654       SUBROUTINE SLVRAR(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
1655      &          M1,M2,NM1,FAC1,E1,LDE1,Z1,F1,IP1,IPHES,IER,IJOB)
1656       IMPLICIT KPP_REAL (A-H,O-Z)
1657       DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),E1(LDE1,NM1),
1658      &          IP1(NM1),IPHES(N),Z1(N),F1(N)
1659       COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
1660 C
1661       GOTO (1,2,3,4,5,6,7,55,55,55,11,12,13,13,15), IJOB
1662 C
1663 C -----------------------------------------------------------
1664 C
1665    1  CONTINUE
1666 C ---  B=IDENTITY, JACOBIAN A FULL MATRIX
1667       DO I=1,N
1668          Z1(I)=Z1(I)-F1(I)*FAC1
1669       END DO
1670       CALL SOL (N,LDE1,E1,Z1,IP1)
1671       RETURN
1672 C
1673 C -----------------------------------------------------------
1674 C
1675   11  CONTINUE
1676 C ---  B=IDENTITY, JACOBIAN A FULL MATRIX, SECOND ORDER
1677       DO I=1,N
1678          Z1(I)=Z1(I)-F1(I)*FAC1
1679       END DO
1680  48   CONTINUE
1681       MM=M1/M2
1682       DO J=1,M2
1683          SUM1=0.D0
1684          DO K=MM-1,0,-1
1685             JKM=J+K*M2
1686             SUM1=(Z1(JKM)+SUM1)/FAC1
1687             DO I=1,NM1
1688                IM1=I+M1
1689                Z1(IM1)=Z1(IM1)+FJAC(I,JKM)*SUM1
1690             END DO
1691          END DO
1692       END DO
1693       CALL SOL (NM1,LDE1,E1,Z1(M1+1),IP1)
1694  49   CONTINUE
1695       DO I=M1,1,-1
1696          Z1(I)=(Z1(I)+Z1(M2+I))/FAC1
1697       END DO
1698       RETURN
1699 C
1700 C -----------------------------------------------------------
1701 C
1702    2  CONTINUE
1703 C ---  B=IDENTITY, JACOBIAN A BANDED MATRIX
1704       DO I=1,N
1705          Z1(I)=Z1(I)-F1(I)*FAC1
1706       END DO
1707       CALL SOLB (N,LDE1,E1,MLE,MUE,Z1,IP1)
1708       RETURN
1709 C
1710 C -----------------------------------------------------------
1711 C
1712   12  CONTINUE
1713 C ---  B=IDENTITY, JACOBIAN A BANDED MATRIX, SECOND ORDER
1714       DO I=1,N
1715          Z1(I)=Z1(I)-F1(I)*FAC1
1716       END DO
1717   45  CONTINUE
1718       MM=M1/M2
1719       DO J=1,M2
1720          SUM1=0.D0
1721          DO K=MM-1,0,-1
1722             JKM=J+K*M2
1723             SUM1=(Z1(JKM)+SUM1)/FAC1
1724             DO I=MAX(1,J-MUJAC),MIN(NM1,J+MLJAC)
1725                IM1=I+M1
1726                Z1(IM1)=Z1(IM1)+FJAC(I+MUJAC+1-J,JKM)*SUM1
1727             END DO
1728          END DO
1729       END DO
1730       CALL SOLB (NM1,LDE1,E1,MLE,MUE,Z1(M1+1),IP1)
1731       GOTO 49
1732 C
1733 C -----------------------------------------------------------
1734 C
1735    3  CONTINUE
1736 C ---  B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX
1737       DO I=1,N
1738          S1=0.0D0
1739          DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS)
1740             S1=S1-FMAS(I-J+MBDIAG,J)*F1(J)
1741          END DO
1742          Z1(I)=Z1(I)+S1*FAC1
1743       END DO
1744       CALL SOL (N,LDE1,E1,Z1,IP1)
1745       RETURN
1746 C
1747 C -----------------------------------------------------------
1748 C
1749   13  CONTINUE
1750 C ---  B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER
1751       DO I=1,M1
1752          Z1(I)=Z1(I)-F1(I)*FAC1
1753       END DO
1754       DO I=1,NM1
1755          IM1=I+M1
1756          S1=0.0D0
1757          DO J=MAX(1,I-MLMAS),MIN(NM1,I+MUMAS)
1758             S1=S1-FMAS(I-J+MBDIAG,J)*F1(J+M1)
1759          END DO
1760          Z1(IM1)=Z1(IM1)+S1*FAC1
1761       END DO
1762       IF (IJOB.EQ.14) GOTO 45
1763       GOTO 48
1764 C
1765 C -----------------------------------------------------------
1766 C
1767    4  CONTINUE
1768 C ---  B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX
1769       DO I=1,N
1770          S1=0.0D0
1771          DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS)
1772             S1=S1-FMAS(I-J+MBDIAG,J)*F1(J)
1773          END DO
1774          Z1(I)=Z1(I)+S1*FAC1
1775       END DO
1776       CALL SOLB (N,LDE1,E1,MLE,MUE,Z1,IP1)
1777       RETURN
1778 C
1779 C -----------------------------------------------------------
1780 C
1781    5  CONTINUE
1782 C ---  B IS A FULL MATRIX, JACOBIAN A FULL MATRIX
1783       DO I=1,N
1784          S1=0.0D0
1785          DO J=1,N
1786             S1=S1-FMAS(I,J)*F1(J)
1787          END DO
1788          Z1(I)=Z1(I)+S1*FAC1
1789       END DO
1790       CALL SOL (N,LDE1,E1,Z1,IP1)
1791       RETURN
1792 C
1793 C -----------------------------------------------------------
1794 C
1795   15  CONTINUE
1796 C ---  B IS A FULL MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER
1797       DO I=1,M1
1798          Z1(I)=Z1(I)-F1(I)*FAC1
1799       END DO
1800       DO I=1,NM1
1801          IM1=I+M1
1802          S1=0.0D0
1803          DO J=1,NM1
1804             S1=S1-FMAS(I,J)*F1(J+M1)
1805          END DO
1806          Z1(IM1)=Z1(IM1)+S1*FAC1
1807       END DO
1808       GOTO 48
1809 C
1810 C -----------------------------------------------------------
1811 C
1812    6  CONTINUE
1813 C ---  B IS A FULL MATRIX, JACOBIAN A BANDED MATRIX
1814 C ---  THIS OPTION IS NOT PROVIDED
1815       RETURN
1816 C
1817 C -----------------------------------------------------------
1818 C
1819    7  CONTINUE
1820 C ---  B=IDENTITY, JACOBIAN A FULL MATRIX, HESSENBERG-OPTION
1821       DO I=1,N
1822          Z1(I)=Z1(I)-F1(I)*FAC1
1823       END DO
1824       DO MM=N-2,1,-1
1825           MP=N-MM
1826           MP1=MP-1
1827           I=IPHES(MP)
1828           IF (I.EQ.MP) GOTO 746
1829           ZSAFE=Z1(MP)
1830           Z1(MP)=Z1(I)
1831           Z1(I)=ZSAFE
1832  746      CONTINUE
1833           DO I=MP+1,N
1834              Z1(I)=Z1(I)-FJAC(I,MP1)*Z1(MP)
1835           END DO
1836        END DO
1837        CALL SOLH(N,LDE1,E1,1,Z1,IP1)
1838        DO MM=1,N-2
1839           MP=N-MM
1840           MP1=MP-1
1841           DO I=MP+1,N
1842              Z1(I)=Z1(I)+FJAC(I,MP1)*Z1(MP)
1843           END DO
1844           I=IPHES(MP)
1845           IF (I.EQ.MP) GOTO 750
1846           ZSAFE=Z1(MP)
1847           Z1(MP)=Z1(I)
1848           Z1(I)=ZSAFE
1849  750      CONTINUE
1850       END DO
1851       RETURN
1852 C
1853 C -----------------------------------------------------------
1854 C
1855   55  CONTINUE
1856       RETURN
1857       END
1858 C
1859 C     END OF SUBROUTINE SLVRAR
1860 C
1861 C ***********************************************************
1862 C
1863       SUBROUTINE SLVRAI(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
1864      &          M1,M2,NM1,ALPHN,BETAN,E2R,E2I,LDE1,Z2,Z3,
1865      &          F2,F3,CONT,IP2,IPHES,IER,IJOB)
1866       IMPLICIT KPP_REAL (A-H,O-Z)
1867       DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),
1868      &          IP2(NM1),IPHES(N),Z2(N),Z3(N),F2(N),F3(N)
1869       DIMENSION E2R(LDE1,NM1),E2I(LDE1,NM1)
1870       COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
1871 C
1872       GOTO (1,2,3,4,5,6,7,55,55,55,11,12,13,13,15), IJOB
1873 C
1874 C -----------------------------------------------------------
1875 C
1876    1  CONTINUE
1877 C ---  B=IDENTITY, JACOBIAN A FULL MATRIX
1878       DO I=1,N
1879          S2=-F2(I)
1880          S3=-F3(I)
1881          Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
1882          Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
1883       END DO
1884       CALL SOLC (N,LDE1,E2R,E2I,Z2,Z3,IP2)
1885       RETURN
1886 C
1887 C -----------------------------------------------------------
1888 C
1889   11  CONTINUE
1890 C ---  B=IDENTITY, JACOBIAN A FULL MATRIX, SECOND ORDER
1891       DO I=1,N
1892          S2=-F2(I)
1893          S3=-F3(I)
1894          Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
1895          Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
1896       END DO
1897  48   ABNO=ALPHN**2+BETAN**2
1898       MM=M1/M2
1899       DO J=1,M2
1900          SUM2=0.D0
1901          SUM3=0.D0
1902          DO K=MM-1,0,-1
1903             JKM=J+K*M2
1904             SUMH=(Z2(JKM)+SUM2)/ABNO
1905             SUM3=(Z3(JKM)+SUM3)/ABNO
1906             SUM2=SUMH*ALPHN+SUM3*BETAN
1907             SUM3=SUM3*ALPHN-SUMH*BETAN
1908             DO I=1,NM1
1909                IM1=I+M1
1910                Z2(IM1)=Z2(IM1)+FJAC(I,JKM)*SUM2
1911                Z3(IM1)=Z3(IM1)+FJAC(I,JKM)*SUM3
1912             END DO
1913          END DO
1914       END DO
1915       CALL SOLC (NM1,LDE1,E2R,E2I,Z2(M1+1),Z3(M1+1),IP2)
1916  49   CONTINUE
1917       DO I=M1,1,-1
1918          MPI=M2+I
1919          Z2I=Z2(I)+Z2(MPI)
1920          Z3I=Z3(I)+Z3(MPI)
1921          Z3(I)=(Z3I*ALPHN-Z2I*BETAN)/ABNO
1922          Z2(I)=(Z2I*ALPHN+Z3I*BETAN)/ABNO
1923       END DO
1924       RETURN
1925 C
1926 C -----------------------------------------------------------
1927 C
1928    2  CONTINUE
1929 C ---  B=IDENTITY, JACOBIAN A BANDED MATRIX
1930       DO I=1,N
1931          S2=-F2(I)
1932          S3=-F3(I)
1933          Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
1934          Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
1935       END DO
1936       CALL SOLBC (N,LDE1,E2R,E2I,MLE,MUE,Z2,Z3,IP2)
1937       RETURN
1938 C
1939 C -----------------------------------------------------------
1940 C
1941   12  CONTINUE
1942 C ---  B=IDENTITY, JACOBIAN A BANDED MATRIX, SECOND ORDER
1943       DO I=1,N
1944          S2=-F2(I)
1945          S3=-F3(I)
1946          Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
1947          Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
1948       END DO
1949   45  ABNO=ALPHN**2+BETAN**2
1950       MM=M1/M2
1951       DO J=1,M2
1952          SUM2=0.D0
1953          SUM3=0.D0
1954          DO K=MM-1,0,-1
1955             JKM=J+K*M2
1956             SUMH=(Z2(JKM)+SUM2)/ABNO
1957             SUM3=(Z3(JKM)+SUM3)/ABNO
1958             SUM2=SUMH*ALPHN+SUM3*BETAN
1959             SUM3=SUM3*ALPHN-SUMH*BETAN
1960             DO I=MAX(1,J-MUJAC),MIN(NM1,J+MLJAC)
1961                IM1=I+M1
1962                IIMU=I+MUJAC+1-J
1963                Z2(IM1)=Z2(IM1)+FJAC(IIMU,JKM)*SUM2
1964                Z3(IM1)=Z3(IM1)+FJAC(IIMU,JKM)*SUM3
1965             END DO
1966          END DO
1967       END DO
1968       CALL SOLBC (NM1,LDE1,E2R,E2I,MLE,MUE,Z2(M1+1),Z3(M1+1),IP2)
1969       GOTO 49
1970 C
1971 C -----------------------------------------------------------
1972 C
1973    3  CONTINUE
1974 C ---  B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX
1975       DO I=1,N
1976          S2=0.0D0
1977          S3=0.0D0
1978          DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS)
1979             BB=FMAS(I-J+MBDIAG,J)
1980             S2=S2-BB*F2(J)
1981             S3=S3-BB*F3(J)
1982          END DO
1983          Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
1984          Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
1985       END DO
1986       CALL SOLC(N,LDE1,E2R,E2I,Z2,Z3,IP2)
1987       RETURN
1988 C
1989 C -----------------------------------------------------------
1990 C
1991   13  CONTINUE
1992 C ---  B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER
1993       DO I=1,M1
1994          S2=-F2(I)
1995          S3=-F3(I)
1996          Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
1997          Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
1998       END DO
1999       DO I=1,NM1
2000          IM1=I+M1
2001          S2=0.0D0
2002          S3=0.0D0
2003          DO J=MAX(1,I-MLMAS),MIN(NM1,I+MUMAS)
2004             JM1=J+M1
2005             BB=FMAS(I-J+MBDIAG,J)
2006             S2=S2-BB*F2(JM1)
2007             S3=S3-BB*F3(JM1)
2008          END DO
2009          Z2(IM1)=Z2(IM1)+S2*ALPHN-S3*BETAN
2010          Z3(IM1)=Z3(IM1)+S3*ALPHN+S2*BETAN
2011       END DO
2012       IF (IJOB.EQ.14) GOTO 45
2013       GOTO 48
2014 C
2015 C -----------------------------------------------------------
2016 C
2017    4  CONTINUE
2018 C ---  B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX
2019       DO I=1,N
2020          S2=0.0D0
2021          S3=0.0D0
2022          DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS)
2023             BB=FMAS(I-J+MBDIAG,J)
2024             S2=S2-BB*F2(J)
2025             S3=S3-BB*F3(J)
2026          END DO
2027          Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
2028          Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
2029       END DO
2030       CALL SOLBC(N,LDE1,E2R,E2I,MLE,MUE,Z2,Z3,IP2)
2031       RETURN
2032 C
2033 C -----------------------------------------------------------
2034 C
2035    5  CONTINUE
2036 C ---  B IS A FULL MATRIX, JACOBIAN A FULL MATRIX
2037       DO I=1,N
2038          S2=0.0D0
2039          S3=0.0D0
2040          DO J=1,N
2041             BB=FMAS(I,J)
2042             S2=S2-BB*F2(J)
2043             S3=S3-BB*F3(J)
2044          END DO
2045          Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
2046          Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
2047       END DO
2048       CALL SOLC(N,LDE1,E2R,E2I,Z2,Z3,IP2)
2049       RETURN
2050 C
2051 C -----------------------------------------------------------
2052 C
2053   15  CONTINUE
2054 C ---  B IS A FULL MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER
2055       DO I=1,M1
2056          S2=-F2(I)
2057          S3=-F3(I)
2058          Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
2059          Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
2060       END DO
2061       DO I=1,NM1
2062          IM1=I+M1
2063          S2=0.0D0
2064          S3=0.0D0
2065          DO J=1,NM1
2066             JM1=J+M1
2067             BB=FMAS(I,J)
2068             S2=S2-BB*F2(JM1)
2069             S3=S3-BB*F3(JM1)
2070          END DO
2071          Z2(IM1)=Z2(IM1)+S2*ALPHN-S3*BETAN
2072          Z3(IM1)=Z3(IM1)+S3*ALPHN+S2*BETAN
2073       END DO
2074       GOTO 48
2075 C
2076 C -----------------------------------------------------------
2077 C
2078    6  CONTINUE
2079 C ---  B IS A FULL MATRIX, JACOBIAN A BANDED MATRIX
2080 C ---  THIS OPTION IS NOT PROVIDED
2081       RETURN
2082 C
2083 C -----------------------------------------------------------
2084 C
2085    7  CONTINUE
2086 C ---  B=IDENTITY, JACOBIAN A FULL MATRIX, HESSENBERG-OPTION
2087       DO I=1,N
2088          S2=-F2(I)
2089          S3=-F3(I)
2090          Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
2091          Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
2092       END DO
2093       DO MM=N-2,1,-1
2094           MP=N-MM
2095           MP1=MP-1
2096           I=IPHES(MP)
2097           IF (I.EQ.MP) GOTO 746
2098           ZSAFE=Z2(MP)
2099           Z2(MP)=Z2(I)
2100           Z2(I)=ZSAFE
2101           ZSAFE=Z3(MP)
2102           Z3(MP)=Z3(I)
2103           Z3(I)=ZSAFE
2104  746      CONTINUE
2105           DO I=MP+1,N
2106              E1IMP=FJAC(I,MP1)
2107              Z2(I)=Z2(I)-E1IMP*Z2(MP)
2108              Z3(I)=Z3(I)-E1IMP*Z3(MP)
2109           END DO
2110        END DO
2111        CALL SOLHC(N,LDE1,E2R,E2I,1,Z2,Z3,IP2)
2112        DO MM=1,N-2
2113           MP=N-MM
2114           MP1=MP-1
2115           DO I=MP+1,N
2116              E1IMP=FJAC(I,MP1)
2117              Z2(I)=Z2(I)+E1IMP*Z2(MP)
2118              Z3(I)=Z3(I)+E1IMP*Z3(MP)
2119           END DO
2120           I=IPHES(MP)
2121           IF (I.EQ.MP) GOTO 750
2122           ZSAFE=Z2(MP)
2123           Z2(MP)=Z2(I)
2124           Z2(I)=ZSAFE
2125           ZSAFE=Z3(MP)
2126           Z3(MP)=Z3(I)
2127           Z3(I)=ZSAFE
2128  750      CONTINUE
2129       END DO
2130       RETURN
2131 C
2132 C -----------------------------------------------------------
2133 C
2134   55  CONTINUE
2135       RETURN
2136       END
2137 C
2138 C     END OF SUBROUTINE SLVRAI
2139 C
2140 C ***********************************************************
2141 C
2142       SUBROUTINE SLVRAD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
2143      &          M1,M2,NM1,FAC1,ALPHN,BETAN,E1,E2R,E2I,LDE1,Z1,Z2,Z3,
2144      &          F1,F2,F3,CONT,IP1,IP2,IPHES,IER,IJOB)
2145       IMPLICIT KPP_REAL (A-H,O-Z)
2146       DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),E1(LDE1,NM1),
2147      &          E2R(LDE1,NM1),E2I(LDE1,NM1),IP1(NM1),IP2(NM1),
2148      &          IPHES(N),Z1(N),Z2(N),Z3(N),F1(N),F2(N),F3(N)
2149       COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
2150 C
2151       GOTO (1,2,3,4,5,6,7,55,55,55,11,12,13,13,15), IJOB
2152 C
2153 C -----------------------------------------------------------
2154 C
2155    1  CONTINUE
2156 C ---  B=IDENTITY, JACOBIAN A FULL MATRIX
2157       DO I=1,N
2158          S2=-F2(I)
2159          S3=-F3(I)
2160          Z1(I)=Z1(I)-F1(I)*FAC1
2161          Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
2162          Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
2163       END DO
2164       CALL SOL (N,LDE1,E1,Z1,IP1)
2165       CALL SOLC (N,LDE1,E2R,E2I,Z2,Z3,IP2)
2166       RETURN
2167 C
2168 C -----------------------------------------------------------
2169 C
2170   11  CONTINUE
2171 C ---  B=IDENTITY, JACOBIAN A FULL MATRIX, SECOND ORDER
2172       DO I=1,N
2173          S2=-F2(I)
2174          S3=-F3(I)
2175          Z1(I)=Z1(I)-F1(I)*FAC1
2176          Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
2177          Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
2178       END DO
2179  48   ABNO=ALPHN**2+BETAN**2
2180       MM=M1/M2
2181       DO J=1,M2
2182          SUM1=0.D0
2183          SUM2=0.D0
2184          SUM3=0.D0
2185          DO K=MM-1,0,-1
2186             JKM=J+K*M2
2187             SUM1=(Z1(JKM)+SUM1)/FAC1
2188             SUMH=(Z2(JKM)+SUM2)/ABNO
2189             SUM3=(Z3(JKM)+SUM3)/ABNO
2190             SUM2=SUMH*ALPHN+SUM3*BETAN
2191             SUM3=SUM3*ALPHN-SUMH*BETAN
2192             DO I=1,NM1
2193                IM1=I+M1
2194                Z1(IM1)=Z1(IM1)+FJAC(I,JKM)*SUM1
2195                Z2(IM1)=Z2(IM1)+FJAC(I,JKM)*SUM2
2196                Z3(IM1)=Z3(IM1)+FJAC(I,JKM)*SUM3
2197             END DO
2198          END DO
2199       END DO
2200       CALL SOL (NM1,LDE1,E1,Z1(M1+1),IP1)
2201       CALL SOLC (NM1,LDE1,E2R,E2I,Z2(M1+1),Z3(M1+1),IP2)
2202  49   CONTINUE
2203       DO I=M1,1,-1
2204          MPI=M2+I
2205          Z1(I)=(Z1(I)+Z1(MPI))/FAC1
2206          Z2I=Z2(I)+Z2(MPI)
2207          Z3I=Z3(I)+Z3(MPI)
2208          Z3(I)=(Z3I*ALPHN-Z2I*BETAN)/ABNO
2209          Z2(I)=(Z2I*ALPHN+Z3I*BETAN)/ABNO
2210       END DO
2211       RETURN
2212 C
2213 C -----------------------------------------------------------
2214 C
2215    2  CONTINUE
2216 C ---  B=IDENTITY, JACOBIAN A BANDED MATRIX
2217       DO I=1,N
2218          S2=-F2(I)
2219          S3=-F3(I)
2220          Z1(I)=Z1(I)-F1(I)*FAC1
2221          Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
2222          Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
2223       END DO
2224       CALL SOLB (N,LDE1,E1,MLE,MUE,Z1,IP1)
2225       CALL SOLBC (N,LDE1,E2R,E2I,MLE,MUE,Z2,Z3,IP2)
2226       RETURN
2227 C
2228 C -----------------------------------------------------------
2229 C
2230   12  CONTINUE
2231 C ---  B=IDENTITY, JACOBIAN A BANDED MATRIX, SECOND ORDER
2232       DO I=1,N
2233          S2=-F2(I)
2234          S3=-F3(I)
2235          Z1(I)=Z1(I)-F1(I)*FAC1
2236          Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
2237          Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
2238       END DO
2239   45  ABNO=ALPHN**2+BETAN**2
2240       MM=M1/M2
2241       DO J=1,M2
2242          SUM1=0.D0
2243          SUM2=0.D0
2244          SUM3=0.D0
2245          DO K=MM-1,0,-1
2246             JKM=J+K*M2
2247             SUM1=(Z1(JKM)+SUM1)/FAC1
2248             SUMH=(Z2(JKM)+SUM2)/ABNO
2249             SUM3=(Z3(JKM)+SUM3)/ABNO
2250             SUM2=SUMH*ALPHN+SUM3*BETAN
2251             SUM3=SUM3*ALPHN-SUMH*BETAN
2252             DO I=MAX(1,J-MUJAC),MIN(NM1,J+MLJAC)
2253                IM1=I+M1
2254                FFJA=FJAC(I+MUJAC+1-J,JKM)
2255                Z1(IM1)=Z1(IM1)+FFJA*SUM1
2256                Z2(IM1)=Z2(IM1)+FFJA*SUM2
2257                Z3(IM1)=Z3(IM1)+FFJA*SUM3
2258             END DO
2259          END DO
2260       END DO
2261       CALL SOLB (NM1,LDE1,E1,MLE,MUE,Z1(M1+1),IP1)
2262       CALL SOLBC (NM1,LDE1,E2R,E2I,MLE,MUE,Z2(M1+1),Z3(M1+1),IP2)
2263       GOTO 49
2264 C
2265 C -----------------------------------------------------------
2266 C
2267    3  CONTINUE
2268 C ---  B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX
2269       DO I=1,N
2270          S1=0.0D0
2271          S2=0.0D0
2272          S3=0.0D0
2273          DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS)
2274             BB=FMAS(I-J+MBDIAG,J)
2275             S1=S1-BB*F1(J)
2276             S2=S2-BB*F2(J)
2277             S3=S3-BB*F3(J)
2278          END DO
2279          Z1(I)=Z1(I)+S1*FAC1
2280          Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
2281          Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
2282       END DO
2283       CALL SOL (N,LDE1,E1,Z1,IP1)
2284       CALL SOLC(N,LDE1,E2R,E2I,Z2,Z3,IP2)
2285       RETURN
2286 C
2287 C -----------------------------------------------------------
2288 C
2289   13  CONTINUE
2290 C ---  B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER
2291       DO I=1,M1
2292          S2=-F2(I)
2293          S3=-F3(I)
2294          Z1(I)=Z1(I)-F1(I)*FAC1
2295          Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
2296          Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
2297       END DO
2298       DO I=1,NM1
2299          IM1=I+M1
2300          S1=0.0D0
2301          S2=0.0D0
2302          S3=0.0D0
2303          J1B=MAX(1,I-MLMAS)
2304          J2B=MIN(NM1,I+MUMAS)
2305          DO J=J1B,J2B
2306             JM1=J+M1
2307             BB=FMAS(I-J+MBDIAG,J)
2308             S1=S1-BB*F1(JM1)
2309             S2=S2-BB*F2(JM1)
2310             S3=S3-BB*F3(JM1)
2311          END DO
2312          Z1(IM1)=Z1(IM1)+S1*FAC1
2313          Z2(IM1)=Z2(IM1)+S2*ALPHN-S3*BETAN
2314          Z3(IM1)=Z3(IM1)+S3*ALPHN+S2*BETAN
2315       END DO
2316       IF (IJOB.EQ.14) GOTO 45
2317       GOTO 48
2318 C
2319 C -----------------------------------------------------------
2320 C
2321    4  CONTINUE
2322 C ---  B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX
2323       DO I=1,N
2324          S1=0.0D0
2325          S2=0.0D0
2326          S3=0.0D0
2327          DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS)
2328             BB=FMAS(I-J+MBDIAG,J)
2329             S1=S1-BB*F1(J)
2330             S2=S2-BB*F2(J)
2331             S3=S3-BB*F3(J)
2332          END DO
2333          Z1(I)=Z1(I)+S1*FAC1
2334          Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
2335          Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
2336       END DO
2337       CALL SOLB (N,LDE1,E1,MLE,MUE,Z1,IP1)
2338       CALL SOLBC(N,LDE1,E2R,E2I,MLE,MUE,Z2,Z3,IP2)
2339       RETURN
2340 C
2341 C -----------------------------------------------------------
2342 C
2343    5  CONTINUE
2344 C ---  B IS A FULL MATRIX, JACOBIAN A FULL MATRIX
2345       DO I=1,N
2346          S1=0.0D0
2347          S2=0.0D0
2348          S3=0.0D0
2349          DO J=1,N
2350             BB=FMAS(I,J)
2351             S1=S1-BB*F1(J)
2352             S2=S2-BB*F2(J)
2353             S3=S3-BB*F3(J)
2354          END DO
2355          Z1(I)=Z1(I)+S1*FAC1
2356          Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
2357          Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
2358       END DO
2359       CALL SOL (N,LDE1,E1,Z1,IP1)
2360       CALL SOLC(N,LDE1,E2R,E2I,Z2,Z3,IP2)
2361       RETURN
2362 C
2363 C -----------------------------------------------------------
2364 C
2365   15  CONTINUE
2366 C ---  B IS A FULL MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER
2367       DO I=1,M1
2368          S2=-F2(I)
2369          S3=-F3(I)
2370          Z1(I)=Z1(I)-F1(I)*FAC1
2371          Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
2372          Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
2373       END DO
2374       DO I=1,NM1
2375          IM1=I+M1
2376          S1=0.0D0
2377          S2=0.0D0
2378          S3=0.0D0
2379          DO J=1,NM1
2380             JM1=J+M1
2381             BB=FMAS(I,J)
2382             S1=S1-BB*F1(JM1)
2383             S2=S2-BB*F2(JM1)
2384             S3=S3-BB*F3(JM1)
2385          END DO
2386          Z1(IM1)=Z1(IM1)+S1*FAC1
2387          Z2(IM1)=Z2(IM1)+S2*ALPHN-S3*BETAN
2388          Z3(IM1)=Z3(IM1)+S3*ALPHN+S2*BETAN
2389       END DO
2390       GOTO 48
2391 C
2392 C -----------------------------------------------------------
2393 C
2394    6  CONTINUE
2395 C ---  B IS A FULL MATRIX, JACOBIAN A BANDED MATRIX
2396 C ---  THIS OPTION IS NOT PROVIDED
2397       RETURN
2398 C
2399 C -----------------------------------------------------------
2400 C
2401    7  CONTINUE
2402 C ---  B=IDENTITY, JACOBIAN A FULL MATRIX, HESSENBERG-OPTION
2403       DO I=1,N
2404          S2=-F2(I)
2405          S3=-F3(I)
2406          Z1(I)=Z1(I)-F1(I)*FAC1
2407          Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
2408          Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
2409       END DO
2410       DO MM=N-2,1,-1
2411           MP=N-MM
2412           MP1=MP-1
2413           I=IPHES(MP)
2414           IF (I.EQ.MP) GOTO 746
2415           ZSAFE=Z1(MP)
2416           Z1(MP)=Z1(I)
2417           Z1(I)=ZSAFE
2418           ZSAFE=Z2(MP)
2419           Z2(MP)=Z2(I)
2420           Z2(I)=ZSAFE
2421           ZSAFE=Z3(MP)
2422           Z3(MP)=Z3(I)
2423           Z3(I)=ZSAFE
2424  746      CONTINUE
2425           DO I=MP+1,N
2426              E1IMP=FJAC(I,MP1)
2427              Z1(I)=Z1(I)-E1IMP*Z1(MP)
2428              Z2(I)=Z2(I)-E1IMP*Z2(MP)
2429              Z3(I)=Z3(I)-E1IMP*Z3(MP)
2430           END DO
2431        END DO
2432        CALL SOLH(N,LDE1,E1,1,Z1,IP1)
2433        CALL SOLHC(N,LDE1,E2R,E2I,1,Z2,Z3,IP2)
2434        DO MM=1,N-2
2435           MP=N-MM
2436           MP1=MP-1
2437           DO I=MP+1,N
2438              E1IMP=FJAC(I,MP1)
2439              Z1(I)=Z1(I)+E1IMP*Z1(MP)
2440              Z2(I)=Z2(I)+E1IMP*Z2(MP)
2441              Z3(I)=Z3(I)+E1IMP*Z3(MP)
2442           END DO
2443           I=IPHES(MP)
2444           IF (I.EQ.MP) GOTO 750
2445           ZSAFE=Z1(MP)
2446           Z1(MP)=Z1(I)
2447           Z1(I)=ZSAFE
2448           ZSAFE=Z2(MP)
2449           Z2(MP)=Z2(I)
2450           Z2(I)=ZSAFE
2451           ZSAFE=Z3(MP)
2452           Z3(MP)=Z3(I)
2453           Z3(I)=ZSAFE
2454  750      CONTINUE
2455       END DO
2456       RETURN
2457 C
2458 C -----------------------------------------------------------
2459 C
2460   55  CONTINUE
2461       RETURN
2462       END
2463 C
2464 C     END OF SUBROUTINE SLVRAD
2465 C
2466 C ***********************************************************
2467 C
2468       SUBROUTINE ESTRAD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
2469      &          H,DD1,DD2,DD3,FCN,NFCN,Y0,Y,IJOB,X,M1,M2,NM1,
2470      &          E1,LDE1,Z1,Z2,Z3,CONT,F1,F2,IP1,IPHES,SCAL,ERR,
2471      &          FIRST,REJECT,FAC1,RPAR,IPAR)
2472       IMPLICIT KPP_REAL (A-H,O-Z)
2473       DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),E1(LDE1,NM1),IP1(NM1),
2474      &     SCAL(N),IPHES(N),Z1(N),Z2(N),Z3(N),F1(N),F2(N),Y0(N),Y(N)
2475       DIMENSION CONT(N),RPAR(1),IPAR(1)
2476       LOGICAL FIRST,REJECT
2477       COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
2478       EXTERNAL FCN
2479       HEE1=DD1/H
2480       HEE2=DD2/H
2481       HEE3=DD3/H
2482       GOTO (1,2,3,4,5,6,7,55,55,55,11,12,13,14,15), IJOB
2483 C
2484    1  CONTINUE
2485 C ------  B=IDENTITY, JACOBIAN A FULL MATRIX
2486       DO  I=1,N
2487          F2(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I)
2488          CONT(I)=F2(I)+Y0(I)
2489       END DO
2490       CALL SOL (N,LDE1,E1,CONT,IP1)
2491       GOTO 77
2492 C
2493   11  CONTINUE
2494 C ------  B=IDENTITY, JACOBIAN A FULL MATRIX, SECOND ORDER
2495       DO I=1,N
2496          F2(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I)
2497          CONT(I)=F2(I)+Y0(I)
2498       END DO
2499   48  MM=M1/M2
2500       DO J=1,M2
2501          SUM1=0.D0
2502          DO K=MM-1,0,-1
2503             SUM1=(CONT(J+K*M2)+SUM1)/FAC1
2504             DO I=1,NM1
2505                IM1=I+M1
2506                CONT(IM1)=CONT(IM1)+FJAC(I,J+K*M2)*SUM1
2507             END DO
2508          END DO
2509       END DO
2510       CALL SOL (NM1,LDE1,E1,CONT(M1+1),IP1)
2511       DO I=M1,1,-1
2512          CONT(I)=(CONT(I)+CONT(M2+I))/FAC1
2513       END DO
2514       GOTO 77
2515 C
2516    2  CONTINUE
2517 C ------  B=IDENTITY, JACOBIAN A BANDED MATRIX
2518       DO I=1,N
2519          F2(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I)
2520          CONT(I)=F2(I)+Y0(I)
2521       END DO
2522       CALL SOLB (N,LDE1,E1,MLE,MUE,CONT,IP1)
2523       GOTO 77
2524 C
2525   12  CONTINUE
2526 C ------  B=IDENTITY, JACOBIAN A BANDED MATRIX, SECOND ORDER
2527       DO I=1,N
2528          F2(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I)
2529          CONT(I)=F2(I)+Y0(I)
2530       END DO
2531   45  MM=M1/M2
2532       DO J=1,M2
2533          SUM1=0.D0
2534          DO K=MM-1,0,-1
2535             SUM1=(CONT(J+K*M2)+SUM1)/FAC1
2536             DO I=MAX(1,J-MUJAC),MIN(NM1,J+MLJAC)
2537                IM1=I+M1
2538                CONT(IM1)=CONT(IM1)+FJAC(I+MUJAC+1-J,J+K*M2)*SUM1
2539             END DO
2540          END DO
2541       END DO
2542       CALL SOLB (NM1,LDE1,E1,MLE,MUE,CONT(M1+1),IP1)
2543       DO I=M1,1,-1
2544          CONT(I)=(CONT(I)+CONT(M2+I))/FAC1
2545       END DO
2546       GOTO 77
2547 C
2548    3  CONTINUE
2549 C ------  B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX
2550       DO I=1,N
2551          F1(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I)
2552       END DO
2553       DO I=1,N
2554          SUM=0.D0
2555          DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS)
2556             SUM=SUM+FMAS(I-J+MBDIAG,J)*F1(J)
2557          END DO
2558          F2(I)=SUM
2559          CONT(I)=SUM+Y0(I)
2560       END DO
2561       CALL SOL (N,LDE1,E1,CONT,IP1)
2562       GOTO 77
2563 C
2564   13  CONTINUE
2565 C ------  B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER
2566       DO I=1,M1
2567          F2(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I)
2568          CONT(I)=F2(I)+Y0(I)
2569       END DO
2570       DO I=M1+1,N
2571          F1(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I)
2572       END DO
2573       DO I=1,NM1
2574          SUM=0.D0
2575          DO J=MAX(1,I-MLMAS),MIN(NM1,I+MUMAS)
2576             SUM=SUM+FMAS(I-J+MBDIAG,J)*F1(J+M1)
2577          END DO
2578          IM1=I+M1
2579          F2(IM1)=SUM
2580          CONT(IM1)=SUM+Y0(IM1)
2581       END DO
2582       GOTO 48
2583 C
2584    4  CONTINUE
2585 C ------  B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX
2586       DO I=1,N
2587          F1(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I)
2588       END DO
2589       DO I=1,N
2590          SUM=0.D0
2591          DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS)
2592             SUM=SUM+FMAS(I-J+MBDIAG,J)*F1(J)
2593          END DO
2594          F2(I)=SUM
2595          CONT(I)=SUM+Y0(I)
2596       END DO
2597       CALL SOLB (N,LDE1,E1,MLE,MUE,CONT,IP1)
2598       GOTO 77
2599 C
2600   14  CONTINUE
2601 C ------  B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX, SECOND ORDER
2602       DO I=1,M1
2603          F2(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I)
2604          CONT(I)=F2(I)+Y0(I)
2605       END DO
2606       DO I=M1+1,N
2607          F1(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I)
2608       END DO
2609       DO I=1,NM1
2610          SUM=0.D0
2611          DO J=MAX(1,I-MLMAS),MIN(NM1,I+MUMAS)
2612             SUM=SUM+FMAS(I-J+MBDIAG,J)*F1(J+M1)
2613          END DO
2614          IM1=I+M1
2615          F2(IM1)=SUM
2616          CONT(IM1)=SUM+Y0(IM1)
2617       END DO
2618       GOTO 45
2619 C
2620    5  CONTINUE
2621 C ------  B IS A FULL MATRIX, JACOBIAN A FULL MATRIX
2622       DO I=1,N
2623          F1(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I)
2624       END DO
2625       DO I=1,N
2626          SUM=0.D0
2627          DO J=1,N
2628             SUM=SUM+FMAS(I,J)*F1(J)
2629          END DO
2630          F2(I)=SUM
2631          CONT(I)=SUM+Y0(I)
2632       END DO
2633       CALL SOL (N,LDE1,E1,CONT,IP1)
2634       GOTO 77
2635 C
2636   15  CONTINUE
2637 C ------  B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER
2638       DO I=1,M1
2639          F2(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I)
2640          CONT(I)=F2(I)+Y0(I)
2641       END DO
2642       DO I=M1+1,N
2643          F1(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I)
2644       END DO
2645       DO I=1,NM1
2646          SUM=0.D0
2647          DO J=1,NM1
2648             SUM=SUM+FMAS(I,J)*F1(J+M1)
2649          END DO
2650          IM1=I+M1
2651          F2(IM1)=SUM
2652          CONT(IM1)=SUM+Y0(IM1)
2653       END DO
2654       GOTO 48
2655 C
2656    6  CONTINUE
2657 C ------  B IS A FULL MATRIX, JACOBIAN A BANDED MATRIX
2658 C ------  THIS OPTION IS NOT PROVIDED
2659       RETURN
2660 C
2661    7  CONTINUE
2662 C ------  B=IDENTITY, JACOBIAN A FULL MATRIX, HESSENBERG-OPTION
2663       DO I=1,N
2664          F2(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I)
2665          CONT(I)=F2(I)+Y0(I)
2666       END DO
2667       DO MM=N-2,1,-1
2668          MP=N-MM
2669          I=IPHES(MP)
2670          IF (I.EQ.MP) GOTO 310
2671          ZSAFE=CONT(MP)
2672          CONT(MP)=CONT(I)
2673          CONT(I)=ZSAFE
2674  310     CONTINUE
2675          DO I=MP+1,N
2676             CONT(I)=CONT(I)-FJAC(I,MP-1)*CONT(MP)
2677          END DO
2678       END DO
2679       CALL SOLH(N,LDE1,E1,1,CONT,IP1)
2680       DO MM=1,N-2
2681          MP=N-MM
2682          DO I=MP+1,N
2683             CONT(I)=CONT(I)+FJAC(I,MP-1)*CONT(MP)
2684          END DO
2685          I=IPHES(MP)
2686          IF (I.EQ.MP) GOTO 440
2687          ZSAFE=CONT(MP)
2688          CONT(MP)=CONT(I)
2689          CONT(I)=ZSAFE
2690  440     CONTINUE
2691       END DO
2692 C
2693 C --------------------------------------
2694 C
2695   77  CONTINUE
2696       ERR=0.D0
2697       DO  I=1,N
2698          ERR=ERR+(CONT(I)/SCAL(I))**2
2699       END DO
2700       ERR=MAX(SQRT(ERR/N),1.D-10)
2701 C
2702       IF (ERR.LT.1.D0) RETURN
2703       IF (FIRST.OR.REJECT) THEN
2704           DO I=1,N
2705              CONT(I)=Y(I)+CONT(I)
2706           END DO
2707           CALL FCN(N,X,CONT,F1,RPAR,IPAR)
2708           NFCN=NFCN+1
2709           DO I=1,N
2710              CONT(I)=F1(I)+F2(I)
2711           END DO
2712           GOTO (31,32,31,32,31,32,33,55,55,55,41,42,41,42,41), IJOB
2713 C ------ FULL MATRIX OPTION
2714   31      CONTINUE
2715           CALL SOL(N,LDE1,E1,CONT,IP1)
2716           GOTO 88
2717 C ------ FULL MATRIX OPTION, SECOND ORDER
2718  41      CONTINUE
2719          DO J=1,M2
2720             SUM1=0.D0
2721             DO K=MM-1,0,-1
2722                SUM1=(CONT(J+K*M2)+SUM1)/FAC1
2723                DO I=1,NM1
2724                   IM1=I+M1
2725                   CONT(IM1)=CONT(IM1)+FJAC(I,J+K*M2)*SUM1
2726                END DO
2727             END DO
2728          END DO
2729          CALL SOL(NM1,LDE1,E1,CONT(M1+1),IP1)
2730          DO I=M1,1,-1
2731             CONT(I)=(CONT(I)+CONT(M2+I))/FAC1
2732          END DO
2733          GOTO 88
2734 C ------ BANDED MATRIX OPTION
2735  32      CONTINUE
2736          CALL SOLB (N,LDE1,E1,MLE,MUE,CONT,IP1)
2737          GOTO 88
2738 C ------ BANDED MATRIX OPTION, SECOND ORDER
2739  42      CONTINUE
2740          DO J=1,M2
2741             SUM1=0.D0
2742             DO K=MM-1,0,-1
2743                SUM1=(CONT(J+K*M2)+SUM1)/FAC1
2744                DO I=MAX(1,J-MUJAC),MIN(NM1,J+MLJAC)
2745                   IM1=I+M1
2746                   CONT(IM1)=CONT(IM1)+FJAC(I+MUJAC+1-J,J+K*M2)*SUM1
2747                END DO
2748             END DO
2749          END DO
2750          CALL SOLB (NM1,LDE1,E1,MLE,MUE,CONT(M1+1),IP1)
2751          DO I=M1,1,-1
2752             CONT(I)=(CONT(I)+CONT(M2+I))/FAC1
2753          END DO
2754           GOTO 88
2755 C ------ HESSENBERG MATRIX OPTION
2756   33      CONTINUE
2757           DO MM=N-2,1,-1
2758              MP=N-MM
2759              I=IPHES(MP)
2760              IF (I.EQ.MP) GOTO 510
2761              ZSAFE=CONT(MP)
2762              CONT(MP)=CONT(I)
2763              CONT(I)=ZSAFE
2764  510         CONTINUE
2765              DO I=MP+1,N
2766                 CONT(I)=CONT(I)-FJAC(I,MP-1)*CONT(MP)
2767              END DO
2768           END DO
2769           CALL SOLH(N,LDE1,E1,1,CONT,IP1)
2770           DO MM=1,N-2
2771              MP=N-MM
2772              DO I=MP+1,N
2773                 CONT(I)=CONT(I)+FJAC(I,MP-1)*CONT(MP)
2774              END DO
2775              I=IPHES(MP)
2776              IF (I.EQ.MP) GOTO 640
2777              ZSAFE=CONT(MP)
2778              CONT(MP)=CONT(I)
2779              CONT(I)=ZSAFE
2780  640         CONTINUE
2781           END DO
2782 C -----------------------------------
2783    88     CONTINUE
2784           ERR=0.D0
2785           DO I=1,N
2786              ERR=ERR+(CONT(I)/SCAL(I))**2
2787           END DO
2788           ERR=MAX(SQRT(ERR/N),1.D-10)
2789        END IF
2790        RETURN
2791 C -----------------------------------------------------------
2792   55   CONTINUE
2793        RETURN
2794        END
2795 C
2796 C     END OF SUBROUTINE ESTRAD
2797 C
2798 C ***********************************************************
2799 C
2800       SUBROUTINE ESTRAV(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
2801      &          H,DD,FCN,NFCN,Y0,Y,IJOB,X,M1,M2,NM1,NS,NNS,
2802      &          E1,LDE1,ZZ,CONT,FF,IP1,IPHES,SCAL,ERR,
2803      &          FIRST,REJECT,FAC1,RPAR,IPAR)
2804       IMPLICIT KPP_REAL (A-H,O-Z)
2805       DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),E1(LDE1,NM1),IP1(NM1),
2806      &     SCAL(N),IPHES(N),ZZ(NNS),FF(NNS),Y0(N),Y(N)
2807       DIMENSION DD(NS),CONT(N),RPAR(1),IPAR(1)
2808       LOGICAL FIRST,REJECT
2809       COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
2810       EXTERNAL FCN
2811
2812       GOTO (1,2,3,4,5,6,7,55,55,55,11,12,13,14,15), IJOB
2813 C
2814    1  CONTINUE
2815 C ------  B=IDENTITY, JACOBIAN A FULL MATRIX
2816       DO  I=1,N
2817          SUM=0.D0
2818          DO K=1,NS
2819             SUM=SUM+DD(K)*ZZ(I+(K-1)*N)
2820          END DO
2821          FF(I+N)=SUM/H
2822          CONT(I)=FF(I+N)+Y0(I)
2823       END DO
2824       CALL SOL (N,LDE1,E1,CONT,IP1)
2825       GOTO 77
2826 C
2827   11  CONTINUE
2828 C ------  B=IDENTITY, JACOBIAN A FULL MATRIX, SECOND ORDER
2829       DO  I=1,N
2830          SUM=0.D0
2831          DO K=1,NS
2832             SUM=SUM+DD(K)*ZZ(I+(K-1)*N)
2833          END DO
2834          FF(I+N)=SUM/H
2835          CONT(I)=FF(I+N)+Y0(I)
2836       END DO
2837   48  MM=M1/M2
2838       DO J=1,M2
2839          SUM1=0.D0
2840          DO K=MM-1,0,-1
2841             SUM1=(CONT(J+K*M2)+SUM1)/FAC1
2842             DO I=1,NM1
2843                IM1=I+M1
2844                CONT(IM1)=CONT(IM1)+FJAC(I,J+K*M2)*SUM1
2845             END DO
2846          END DO
2847       END DO
2848       CALL SOL (NM1,LDE1,E1,CONT(M1+1),IP1)
2849       DO I=M1,1,-1
2850          CONT(I)=(CONT(I)+CONT(M2+I))/FAC1
2851       END DO
2852       GOTO 77
2853 C
2854    2  CONTINUE
2855 C ------  B=IDENTITY, JACOBIAN A BANDED MATRIX
2856       DO  I=1,N
2857          SUM=0.D0
2858          DO K=1,NS
2859             SUM=SUM+DD(K)*ZZ(I+(K-1)*N)
2860          END DO
2861          FF(I+N)=SUM/H
2862          CONT(I)=FF(I+N)+Y0(I)
2863       END DO
2864       CALL SOLB (N,LDE1,E1,MLE,MUE,CONT,IP1)
2865       GOTO 77
2866 C
2867   12  CONTINUE
2868 C ------  B=IDENTITY, JACOBIAN A BANDED MATRIX, SECOND ORDER
2869       DO  I=1,N
2870          SUM=0.D0
2871          DO K=1,NS
2872             SUM=SUM+DD(K)*ZZ(I+(K-1)*N)
2873          END DO
2874          FF(I+N)=SUM/H
2875          CONT(I)=FF(I+N)+Y0(I)
2876       END DO
2877   45  MM=M1/M2
2878       DO J=1,M2
2879          SUM1=0.D0
2880          DO K=MM-1,0,-1
2881             SUM1=(CONT(J+K*M2)+SUM1)/FAC1
2882             DO I=MAX(1,J-MUJAC),MIN(NM1,J+MLJAC)
2883                IM1=I+M1
2884                CONT(IM1)=CONT(IM1)+FJAC(I+MUJAC+1-J,J+K*M2)*SUM1
2885             END DO
2886          END DO
2887       END DO
2888       CALL SOLB (NM1,LDE1,E1,MLE,MUE,CONT(M1+1),IP1)
2889       DO I=M1,1,-1
2890          CONT(I)=(CONT(I)+CONT(M2+I))/FAC1
2891       END DO
2892       GOTO 77
2893 C
2894    3  CONTINUE
2895 C ------  B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX
2896       DO  I=1,N
2897          SUM=0.D0
2898          DO K=1,NS
2899             SUM=SUM+DD(K)*ZZ(I+(K-1)*N)
2900          END DO
2901          FF(I)=SUM/H
2902       END DO
2903       DO I=1,N
2904          SUM=0.D0
2905          DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS)
2906             SUM=SUM+FMAS(I-J+MBDIAG,J)*FF(J)
2907          END DO
2908          FF(I+N)=SUM
2909          CONT(I)=SUM+Y0(I)
2910       END DO
2911       CALL SOL (N,LDE1,E1,CONT,IP1)
2912       GOTO 77
2913 C
2914   13  CONTINUE
2915 C ------  B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER
2916       DO  I=1,M1
2917          SUM=0.D0
2918          DO K=1,NS
2919             SUM=SUM+DD(K)*ZZ(I+(K-1)*N)
2920          END DO
2921          FF(I+N)=SUM/H
2922          CONT(I)=FF(I+N)+Y0(I)
2923       END DO
2924       DO I=M1+1,N
2925          SUM=0.D0
2926          DO K=1,NS
2927             SUM=SUM+DD(K)*ZZ(I+(K-1)*N)
2928          END DO
2929          FF(I)=SUM/H
2930       END DO
2931       DO I=1,NM1
2932          SUM=0.D0
2933          DO J=MAX(1,I-MLMAS),MIN(NM1,I+MUMAS)
2934             SUM=SUM+FMAS(I-J+MBDIAG,J)*FF(J+M1)
2935          END DO
2936          IM1=I+M1
2937          FF(IM1+N)=SUM
2938          CONT(IM1)=SUM+Y0(IM1)
2939       END DO
2940       GOTO 48
2941 C
2942    4  CONTINUE
2943 C ------  B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX
2944       DO  I=1,N
2945          SUM=0.D0
2946          DO K=1,NS
2947             SUM=SUM+DD(K)*ZZ(I+(K-1)*N)
2948          END DO
2949          FF(I)=SUM/H
2950       END DO
2951       DO I=1,N
2952          SUM=0.D0
2953          DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS)
2954             SUM=SUM+FMAS(I-J+MBDIAG,J)*FF(J)
2955          END DO
2956          FF(I+N)=SUM
2957          CONT(I)=SUM+Y0(I)
2958       END DO
2959       CALL SOLB (N,LDE1,E1,MLE,MUE,CONT,IP1)
2960       GOTO 77
2961 C
2962   14  CONTINUE
2963 C ------  B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX, SECOND ORDER
2964       DO  I=1,M1
2965          SUM=0.D0
2966          DO K=1,NS
2967             SUM=SUM+DD(K)*ZZ(I+(K-1)*N)
2968          END DO
2969          FF(I+N)=SUM/H
2970          CONT(I)=FF(I+N)+Y0(I)
2971       END DO
2972       DO I=M1+1,N
2973          SUM=0.D0
2974          DO K=1,NS
2975             SUM=SUM+DD(K)*ZZ(I+(K-1)*N)
2976          END DO
2977          FF(I)=SUM/H
2978       END DO
2979       DO I=1,NM1
2980          SUM=0.D0
2981          DO J=MAX(1,I-MLMAS),MIN(NM1,I+MUMAS)
2982             SUM=SUM+FMAS(I-J+MBDIAG,J)*FF(J+M1)
2983          END DO
2984          IM1=I+M1
2985          FF(IM1+N)=SUM
2986          CONT(IM1)=SUM+Y0(IM1)
2987       END DO
2988       GOTO 45
2989 C
2990    5  CONTINUE
2991 C ------  B IS A FULL MATRIX, JACOBIAN A FULL MATRIX
2992       DO I=1,N
2993          SUM=0.D0
2994          DO K=1,NS
2995             SUM=SUM+DD(K)*ZZ(I+(K-1)*N)
2996          END DO
2997          FF(I)=SUM/H
2998       END DO
2999       DO I=1,N
3000          SUM=0.D0
3001          DO J=1,N
3002             SUM=SUM+FMAS(I,J)*FF(J)
3003          END DO
3004          FF(I+N)=SUM
3005          CONT(I)=SUM+Y0(I)
3006       END DO
3007       CALL SOL (N,LDE1,E1,CONT,IP1)
3008       GOTO 77
3009 C
3010   15  CONTINUE
3011 C ------  B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER
3012       DO  I=1,M1
3013          SUM=0.D0
3014          DO K=1,NS
3015             SUM=SUM+DD(K)*ZZ(I+(K-1)*N)
3016          END DO
3017          FF(I+N)=SUM/H
3018          CONT(I)=FF(I+N)+Y0(I)
3019       END DO
3020       DO I=M1+1,N
3021          SUM=0.D0
3022          DO K=1,NS
3023             SUM=SUM+DD(K)*ZZ(I+(K-1)*N)
3024          END DO
3025          FF(I)=SUM/H
3026       END DO
3027       DO I=1,NM1
3028          SUM=0.D0
3029          DO J=1,NM1
3030             SUM=SUM+FMAS(I,J)*FF(J+M1)
3031          END DO
3032          IM1=I+M1
3033          FF(IM1+N)=SUM
3034          CONT(IM1)=SUM+Y0(IM1)
3035       END DO
3036       GOTO 48
3037 C
3038    6  CONTINUE
3039 C ------  B IS A FULL MATRIX, JACOBIAN A BANDED MATRIX
3040 C ------  THIS OPTION IS NOT PROVIDED
3041       RETURN
3042 C
3043    7  CONTINUE
3044 C ------  B=IDENTITY, JACOBIAN A FULL MATRIX, HESSENBERG-OPTION
3045       DO  I=1,N
3046          SUM=0.D0
3047          DO K=1,NS
3048             SUM=SUM+DD(K)*ZZ(I+(K-1)*N)
3049          END DO
3050          FF(I+N)=SUM/H
3051          CONT(I)=FF(I+N)+Y0(I)
3052       END DO
3053       DO MM=N-2,1,-1
3054          MP=N-MM
3055          I=IPHES(MP)
3056          IF (I.EQ.MP) GOTO 310
3057          ZSAFE=CONT(MP)
3058          CONT(MP)=CONT(I)
3059          CONT(I)=ZSAFE
3060  310     CONTINUE
3061          DO I=MP+1,N
3062             CONT(I)=CONT(I)-FJAC(I,MP-1)*CONT(MP)
3063          END DO
3064       END DO
3065       CALL SOLH(N,LDE1,E1,1,CONT,IP1)
3066       DO MM=1,N-2
3067          MP=N-MM
3068          DO I=MP+1,N
3069             CONT(I)=CONT(I)+FJAC(I,MP-1)*CONT(MP)
3070          END DO
3071          I=IPHES(MP)
3072          IF (I.EQ.MP) GOTO 440
3073          ZSAFE=CONT(MP)
3074          CONT(MP)=CONT(I)
3075          CONT(I)=ZSAFE
3076  440     CONTINUE
3077       END DO
3078 C
3079 C --------------------------------------
3080 C
3081   77  CONTINUE
3082       ERR=0.D0
3083       DO  I=1,N
3084          ERR=ERR+(CONT(I)/SCAL(I))**2
3085       END DO
3086       ERR=MAX(SQRT(ERR/N),1.D-10)
3087 C
3088       IF (ERR.LT.1.D0) RETURN
3089       IF (FIRST.OR.REJECT) THEN
3090           DO I=1,N
3091              CONT(I)=Y(I)+CONT(I)
3092           END DO
3093           CALL FCN(N,X,CONT,FF,RPAR,IPAR)
3094           NFCN=NFCN+1
3095           DO I=1,N
3096              CONT(I)=FF(I)+FF(I+N)
3097           END DO
3098           GOTO (31,32,31,32,31,32,33,55,55,55,41,42,41,42,41), IJOB
3099 C ------ FULL MATRIX OPTION
3100  31      CONTINUE
3101          CALL SOL (N,LDE1,E1,CONT,IP1)
3102           GOTO 88
3103 C ------ FULL MATRIX OPTION, SECOND ORDER
3104  41      CONTINUE
3105          DO J=1,M2
3106             SUM1=0.D0
3107             DO K=MM-1,0,-1
3108                SUM1=(CONT(J+K*M2)+SUM1)/FAC1
3109                DO I=1,NM1
3110                   IM1=I+M1
3111                   CONT(IM1)=CONT(IM1)+FJAC(I,J+K*M2)*SUM1
3112                END DO
3113             END DO
3114          END DO
3115          CALL SOL (NM1,LDE1,E1,CONT(M1+1),IP1)
3116          DO I=M1,1,-1
3117             CONT(I)=(CONT(I)+CONT(M2+I))/FAC1
3118          END DO
3119           GOTO 88
3120 C ------ BANDED MATRIX OPTION
3121  32      CONTINUE
3122          CALL SOLB (N,LDE1,E1,MLE,MUE,CONT,IP1)
3123           GOTO 88
3124 C ------ BANDED MATRIX OPTION, SECOND ORDER
3125  42      CONTINUE
3126          DO J=1,M2
3127             SUM1=0.D0
3128             DO K=MM-1,0,-1
3129                SUM1=(CONT(J+K*M2)+SUM1)/FAC1
3130                DO I=MAX(1,J-MUJAC),MIN(NM1,J+MLJAC)
3131                   IM1=I+M1
3132                   CONT(IM1)=CONT(IM1)+FJAC(I+MUJAC+1-J,J+K*M2)*SUM1
3133                END DO
3134             END DO
3135          END DO
3136          CALL SOLB (NM1,LDE1,E1,MLE,MUE,CONT(M1+1),IP1)
3137          DO I=M1,1,-1
3138             CONT(I)=(CONT(I)+CONT(M2+I))/FAC1
3139          END DO
3140           GOTO 88
3141 C ------ HESSENBERG MATRIX OPTION
3142   33      CONTINUE
3143           DO MM=N-2,1,-1
3144              MP=N-MM
3145              I=IPHES(MP)
3146              IF (I.EQ.MP) GOTO 510
3147              ZSAFE=CONT(MP)
3148              CONT(MP)=CONT(I)
3149              CONT(I)=ZSAFE
3150  510         CONTINUE
3151              DO I=MP+1,N
3152                 CONT(I)=CONT(I)-FJAC(I,MP-1)*CONT(MP)
3153              END DO
3154           END DO
3155           CALL SOLH(N,LDE1,E1,1,CONT,IP1)
3156           DO MM=1,N-2
3157              MP=N-MM
3158              DO I=MP+1,N
3159                 CONT(I)=CONT(I)+FJAC(I,MP-1)*CONT(MP)
3160              END DO
3161              I=IPHES(MP)
3162              IF (I.EQ.MP) GOTO 640
3163              ZSAFE=CONT(MP)
3164              CONT(MP)=CONT(I)
3165              CONT(I)=ZSAFE
3166  640         CONTINUE
3167           END DO
3168 C -----------------------------------
3169   88      CONTINUE
3170           ERR=0.D0
3171           DO I=1,N
3172              ERR=ERR+(CONT(I)/SCAL(I))**2
3173           END DO
3174           ERR=MAX(SQRT(ERR/N),1.D-10)
3175        END IF
3176        RETURN
3177 C
3178 C -----------------------------------------------------------
3179 C
3180   55  CONTINUE
3181       RETURN
3182        END
3183 C
3184 C     END OF SUBROUTINE ESTRAV
3185 C
3186 C ***********************************************************
3187 C
3188       SUBROUTINE SLVROD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
3189      &          M1,M2,NM1,FAC1,E,LDE,IP,DY,AK,FX,YNEW,HD,IJOB,STAGE1)
3190       IMPLICIT KPP_REAL (A-H,O-Z)
3191       DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),E(LDE,NM1),
3192      &          IP(NM1),DY(N),AK(N),FX(N),YNEW(N)
3193       LOGICAL STAGE1
3194       COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
3195 C
3196       IF (HD.EQ.0.D0) THEN
3197          DO  I=1,N
3198            AK(I)=DY(I)
3199          END DO
3200       ELSE
3201          DO I=1,N
3202             AK(I)=DY(I)+HD*FX(I)
3203          END DO
3204       END IF
3205 C
3206       GOTO (1,2,3,4,5,6,55,55,55,55,11,12,13,13,15), IJOB
3207 C
3208 C -----------------------------------------------------------
3209 C
3210    1  CONTINUE
3211 C ---  B=IDENTITY, JACOBIAN A FULL MATRIX
3212       IF (STAGE1) THEN
3213          DO I=1,N
3214             AK(I)=AK(I)+YNEW(I)
3215          END DO
3216       END IF
3217       CALL SOL (N,LDE,E,AK,IP)
3218       RETURN
3219 C
3220 C -----------------------------------------------------------
3221 C
3222   11  CONTINUE
3223 C ---  B=IDENTITY, JACOBIAN A FULL MATRIX, SECOND ORDER
3224       IF (STAGE1) THEN
3225          DO I=1,N
3226             AK(I)=AK(I)+YNEW(I)
3227          END DO
3228       END IF
3229  48   MM=M1/M2
3230       DO J=1,M2
3231          SUM=0.D0
3232          DO K=MM-1,0,-1
3233             JKM=J+K*M2
3234             SUM=(AK(JKM)+SUM)/FAC1
3235             DO I=1,NM1
3236                IM1=I+M1
3237                AK(IM1)=AK(IM1)+FJAC(I,JKM)*SUM
3238             END DO
3239          END DO
3240       END DO
3241       CALL SOL (NM1,LDE,E,AK(M1+1),IP)
3242       DO I=M1,1,-1
3243          AK(I)=(AK(I)+AK(M2+I))/FAC1
3244       END DO
3245       RETURN
3246 C
3247 C -----------------------------------------------------------
3248 C
3249    2  CONTINUE
3250 C ---  B=IDENTITY, JACOBIAN A BANDED MATRIX
3251       IF (STAGE1) THEN
3252          DO I=1,N
3253             AK(I)=AK(I)+YNEW(I)
3254          END DO
3255       END IF
3256       CALL SOLB (N,LDE,E,MLE,MUE,AK,IP)
3257       RETURN
3258 C
3259 C -----------------------------------------------------------
3260 C
3261   12  CONTINUE
3262 C ---  B=IDENTITY, JACOBIAN A BANDED MATRIX, SECOND ORDER
3263       IF (STAGE1) THEN
3264          DO I=1,N
3265             AK(I)=AK(I)+YNEW(I)
3266          END DO
3267       END IF
3268   45  MM=M1/M2
3269       DO J=1,M2
3270          SUM=0.D0
3271          DO K=MM-1,0,-1
3272             JKM=J+K*M2
3273             SUM=(AK(JKM)+SUM)/FAC1
3274             DO I=MAX(1,J-MUJAC),MIN(NM1,J+MLJAC)
3275                IM1=I+M1
3276                AK(IM1)=AK(IM1)+FJAC(I+MUJAC+1-J,JKM)*SUM
3277             END DO
3278          END DO
3279       END DO
3280       CALL SOLB (NM1,LDE,E,MLE,MUE,AK(M1+1),IP)
3281       DO I=M1,1,-1
3282          AK(I)=(AK(I)+AK(M2+I))/FAC1
3283       END DO
3284       RETURN
3285 C
3286 C -----------------------------------------------------------
3287 C
3288    3  CONTINUE
3289 C ---  B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX
3290       IF (STAGE1) THEN
3291       DO  I=1,N
3292          SUM=0.D0
3293          DO  J=MAX(1,I-MLMAS),MIN(N,I+MUMAS)
3294             SUM=SUM+FMAS(I-J+MBDIAG,J)*YNEW(J)
3295          END DO
3296          AK(I)=AK(I)+SUM
3297       END DO
3298       END IF
3299       CALL SOL (N,LDE,E,AK,IP)
3300       RETURN
3301 C
3302 C -----------------------------------------------------------
3303 C
3304   13  CONTINUE
3305 C ---  B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER
3306       IF (STAGE1) THEN
3307          DO I=1,M1
3308             AK(I)=AK(I)+YNEW(I)
3309          END DO
3310          DO I=1,NM1
3311             SUM=0.D0
3312             DO J=MAX(1,I-MLMAS),MIN(NM1,I+MUMAS)
3313                 SUM=SUM+FMAS(I-J+MBDIAG,J)*YNEW(J+M1)
3314             END DO
3315             IM1=I+M1
3316             AK(IM1)=AK(IM1)+SUM
3317          END DO
3318       END IF
3319       IF (IJOB.EQ.14) GOTO 45
3320       GOTO 48
3321 C
3322 C -----------------------------------------------------------
3323 C
3324    4  CONTINUE
3325 C ---  B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX
3326       IF (STAGE1) THEN
3327       DO I=1,N
3328          SUM=0.D0
3329          DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS)
3330             SUM=SUM+FMAS(I-J+MBDIAG,J)*YNEW(J)
3331          END DO
3332          AK(I)=AK(I)+SUM
3333       END DO
3334       END IF
3335       CALL SOLB (N,LDE,E,MLE,MUE,AK,IP)
3336       RETURN
3337 C
3338 C -----------------------------------------------------------
3339 C
3340    5  CONTINUE
3341 C ---  B IS A FULL MATRIX, JACOBIAN A FULL MATRIX
3342       IF (STAGE1) THEN
3343       DO I=1,N
3344          SUM=0.D0
3345          DO J=1,N
3346             SUM=SUM+FMAS(I,J)*YNEW(J)
3347          END DO
3348          AK(I)=AK(I)+SUM
3349       END DO
3350       END IF
3351       CALL SOL (N,LDE,E,AK,IP)
3352       RETURN
3353 C
3354 C -----------------------------------------------------------
3355 C
3356   15  CONTINUE
3357 C ---  B IS A FULL MATRIX, JACOBIAN A FULL MATRIX, SECOND ORDER
3358       IF (STAGE1) THEN
3359          DO I=1,M1
3360             AK(I)=AK(I)+YNEW(I)
3361          END DO
3362          DO I=1,NM1
3363             SUM=0.D0
3364             DO J=1,NM1
3365                SUM=SUM+FMAS(I,J)*YNEW(J+M1)
3366             END DO
3367             IM1=I+M1
3368             AK(IM1)=AK(IM1)+SUM
3369          END DO
3370       END IF
3371       GOTO 48
3372 C
3373 C -----------------------------------------------------------
3374 C
3375    6  CONTINUE
3376 C ---  B IS A FULL MATRIX, JACOBIAN A BANDED MATRIX
3377 C ---  THIS OPTION IS NOT PROVIDED
3378       IF (STAGE1) THEN
3379       DO 624 I=1,N
3380          SUM=0.D0
3381          DO 623 J=1,N
3382   623       SUM=SUM+FMAS(I,J)*YNEW(J)
3383   624    AK(I)=AK(I)+SUM
3384       CALL SOLB (N,LDE,E,MLE,MUE,AK,IP)
3385       END IF
3386       RETURN
3387 C
3388 C -----------------------------------------------------------
3389 C
3390   55  CONTINUE
3391       RETURN
3392       END
3393 C
3394 C     END OF SUBROUTINE SLVROD
3395 C
3396 C
3397 C ***********************************************************
3398 C
3399       SUBROUTINE SLVSEU(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS,
3400      &          M1,M2,NM1,FAC1,E,LDE,IP,IPHES,DEL,IJOB)
3401       IMPLICIT KPP_REAL (A-H,O-Z)
3402       DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),E(LDE,NM1),DEL(N)
3403       DIMENSION IP(NM1),IPHES(N)
3404       COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
3405 C
3406       GOTO (1,2,1,2,1,55,7,55,55,55,11,12,11,12,11), IJOB
3407 C
3408 C -----------------------------------------------------------
3409 C
3410    1  CONTINUE
3411 C ---  B=IDENTITY, JACOBIAN A FULL MATRIX
3412       CALL SOL (N,LDE,E,DEL,IP)
3413       RETURN
3414 C
3415 C -----------------------------------------------------------
3416 C
3417   11  CONTINUE
3418 C ---  B=IDENTITY, JACOBIAN A FULL MATRIX, SECOND ORDER
3419       MM=M1/M2
3420       DO J=1,M2
3421          SUM=0.D0
3422          DO K=MM-1,0,-1
3423             JKM=J+K*M2
3424             SUM=(DEL(JKM)+SUM)/FAC1
3425             DO I=1,NM1
3426                IM1=I+M1
3427                DEL(IM1)=DEL(IM1)+FJAC(I,JKM)*SUM
3428             END DO
3429          END DO
3430       END DO
3431       CALL SOL (NM1,LDE,E,DEL(M1+1),IP)
3432       DO I=M1,1,-1
3433          DEL(I)=(DEL(I)+DEL(M2+I))/FAC1
3434       END DO
3435       RETURN
3436 C
3437 C -----------------------------------------------------------
3438 C
3439    2  CONTINUE
3440 C ---  B=IDENTITY, JACOBIAN A BANDED MATRIX
3441       CALL SOLB (N,LDE,E,MLE,MUE,DEL,IP)
3442       RETURN
3443 C
3444 C -----------------------------------------------------------
3445 C
3446   12  CONTINUE
3447 C ---  B=IDENTITY, JACOBIAN A BANDED MATRIX, SECOND ORDER
3448       MM=M1/M2
3449       DO J=1,M2
3450          SUM=0.D0
3451          DO K=MM-1,0,-1
3452             JKM=J+K*M2
3453             SUM=(DEL(JKM)+SUM)/FAC1
3454             DO I=MAX(1,J-MUJAC),MIN(NM1,J+MLJAC)
3455                IM1=I+M1
3456                DEL(IM1)=DEL(IM1)+FJAC(I+MUJAC+1-J,JKM)*SUM
3457             END DO
3458          END DO
3459       END DO
3460       CALL SOLB (NM1,LDE,E,MLE,MUE,DEL(M1+1),IP)
3461       DO I=M1,1,-1
3462          DEL(I)=(DEL(I)+DEL(M2+I))/FAC1
3463       END DO
3464       RETURN
3465 C
3466 C -----------------------------------------------------------
3467 C
3468    7  CONTINUE
3469 C ---  HESSENBERG OPTION
3470       DO MMM=N-2,1,-1
3471          MP=N-MMM
3472          MP1=MP-1
3473          I=IPHES(MP)
3474          IF (I.EQ.MP) GOTO 110
3475          ZSAFE=DEL(MP)
3476          DEL(MP)=DEL(I)
3477          DEL(I)=ZSAFE
3478  110     CONTINUE
3479          DO I=MP+1,N
3480             DEL(I)=DEL(I)-FJAC(I,MP1)*DEL(MP)
3481          END DO
3482       END DO
3483       CALL SOLH(N,LDE,E,1,DEL,IP)
3484       DO MMM=1,N-2
3485          MP=N-MMM
3486          MP1=MP-1
3487          DO I=MP+1,N
3488             DEL(I)=DEL(I)+FJAC(I,MP1)*DEL(MP)
3489          END DO
3490          I=IPHES(MP)
3491          IF (I.EQ.MP) GOTO 240
3492          ZSAFE=DEL(MP)
3493          DEL(MP)=DEL(I)
3494          DEL(I)=ZSAFE
3495  240     CONTINUE
3496       END DO
3497       RETURN
3498 C
3499 C -----------------------------------------------------------
3500 C
3501   55  CONTINUE
3502       RETURN
3503       END
3504 C
3505 C     END OF SUBROUTINE SLVSEU
3506 C
3507       SUBROUTINE DEC (N, NDIM, A, IP, IER)
3508 C VERSION REAL KPP_REAL
3509       INTEGER N,NDIM,IP,IER,NM1,K,KP1,M,I,J
3510       KPP_REAL A,T
3511       DIMENSION A(NDIM,N), IP(N)
3512 C-----------------------------------------------------------------------
3513 C  MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION.
3514 C  INPUT..
3515 C     N = ORDER OF MATRIX.
3516 C     NDIM = DECLARED DIMENSION OF ARRAY  A .
3517 C     A = MATRIX TO BE TRIANGULARIZED.
3518 C  OUTPUT..
3519 C     A(I,J), I.LE.J = UPPER TRIANGULAR FACTOR, U .
3520 C     A(I,J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I - L.
3521 C     IP(K), K.LT.N = INDEX OF K-TH PIVOT ROW.
3522 C     IP(N) = (-1)**(NUMBER OF INTERCHANGES) OR O .
3523 C     IER = 0 IF MATRIX A IS NONSINGULAR, OR K IF FOUND TO BE
3524 C           SINGULAR AT STAGE K.
3525 C  USE  SOL  TO OBTAIN SOLUTION OF LINEAR SYSTEM.
3526 C  DETERM(A) = IP(N)*A(1,1)*A(2,2)*...*A(N,N).
3527 C  IF IP(N)=O, A IS SINGULAR, SOL WILL DIVIDE BY ZERO.
3528 C
3529 C  REFERENCE..
3530 C     C. B. MOLER, ALGORITHM 423, LINEAR EQUATION KppSolveR,
3531 C     C.A.C.M. 15 (1972), P. 274.
3532 C-----------------------------------------------------------------------
3533       IER = 0
3534       IP(N) = 1
3535       IF (N .EQ. 1) GO TO 70
3536       NM1 = N - 1
3537       DO 60 K = 1,NM1
3538         KP1 = K + 1
3539         M = K
3540         DO 10 I = KP1,N
3541           IF (DABS(A(I,K)) .GT. DABS(A(M,K))) M = I
3542  10     CONTINUE
3543         IP(K) = M
3544         T = A(M,K)
3545         IF (M .EQ. K) GO TO 20
3546         IP(N) = -IP(N)
3547         A(M,K) = A(K,K)
3548         A(K,K) = T
3549  20     CONTINUE
3550         IF (T .EQ. 0.D0) GO TO 80
3551         T = 1.D0/T
3552         DO 30 I = KP1,N
3553  30       A(I,K) = -A(I,K)*T
3554         DO 50 J = KP1,N
3555           T = A(M,J)
3556           A(M,J) = A(K,J)
3557           A(K,J) = T
3558           IF (T .EQ. 0.D0) GO TO 45
3559           DO 40 I = KP1,N
3560  40         A(I,J) = A(I,J) + A(I,K)*T
3561  45       CONTINUE
3562  50       CONTINUE
3563  60     CONTINUE
3564  70   K = N
3565       IF (A(N,N) .EQ. 0.D0) GO TO 80
3566       RETURN
3567  80   IER = K
3568       IP(N) = 0
3569       RETURN
3570 C----------------------- END OF SUBROUTINE DEC -------------------------
3571       END
3572 C
3573 C
3574       SUBROUTINE SOL (N, NDIM, A, B, IP)
3575 C VERSION REAL KPP_REAL
3576       INTEGER N,NDIM,IP,NM1,K,KP1,M,I,KB,KM1
3577       KPP_REAL A,B,T
3578       DIMENSION A(NDIM,N), B(N), IP(N)
3579 C-----------------------------------------------------------------------
3580 C  SOLUTION OF LINEAR SYSTEM, A*X = B .
3581 C  INPUT..
3582 C    N = ORDER OF MATRIX.
3583 C    NDIM = DECLARED DIMENSION OF ARRAY  A .
3584 C    A = TRIANGULARIZED MATRIX OBTAINED FROM DEC.
3585 C    B = RIGHT HAND SIDE VECTOR.
3586 C    IP = PIVOT VECTOR OBTAINED FROM DEC.
3587 C  DO NOT USE IF DEC HAS SET IER .NE. 0.
3588 C  OUTPUT..
3589 C    B = SOLUTION VECTOR, X .
3590 C-----------------------------------------------------------------------
3591       IF (N .EQ. 1) GO TO 50
3592       NM1 = N - 1
3593       DO 20 K = 1,NM1
3594         KP1 = K + 1
3595         M = IP(K)
3596         T = B(M)
3597         B(M) = B(K)
3598         B(K) = T
3599         DO 10 I = KP1,N
3600  10       B(I) = B(I) + A(I,K)*T
3601  20     CONTINUE
3602       DO 40 KB = 1,NM1
3603         KM1 = N - KB
3604         K = KM1 + 1
3605         B(K) = B(K)/A(K,K)
3606         T = -B(K)
3607         DO 30 I = 1,KM1
3608  30       B(I) = B(I) + A(I,K)*T
3609  40     CONTINUE
3610  50   B(1) = B(1)/A(1,1)
3611       RETURN
3612 C----------------------- END OF SUBROUTINE SOL -------------------------
3613       END
3614 c
3615 c
3616       SUBROUTINE DECH (N, NDIM, A, LB, IP, IER)
3617 C VERSION REAL KPP_REAL
3618       INTEGER N,NDIM,IP,IER,NM1,K,KP1,M,I,J,LB,NA
3619       KPP_REAL A,T
3620       DIMENSION A(NDIM,N), IP(N)
3621 C-----------------------------------------------------------------------
3622 C  MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION OF A HESSENBERG
3623 C  MATRIX WITH LOWER BANDWIDTH LB
3624 C  INPUT..
3625 C     N = ORDER OF MATRIX A.
3626 C     NDIM = DECLARED DIMENSION OF ARRAY  A .
3627 C     A = MATRIX TO BE TRIANGULARIZED.
3628 C     LB = LOWER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED, LB.GE.1).
3629 C  OUTPUT..
3630 C     A(I,J), I.LE.J = UPPER TRIANGULAR FACTOR, U .
3631 C     A(I,J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I - L.
3632 C     IP(K), K.LT.N = INDEX OF K-TH PIVOT ROW.
3633 C     IP(N) = (-1)**(NUMBER OF INTERCHANGES) OR O .
3634 C     IER = 0 IF MATRIX A IS NONSINGULAR, OR K IF FOUND TO BE
3635 C           SINGULAR AT STAGE K.
3636 C  USE  SOLH  TO OBTAIN SOLUTION OF LINEAR SYSTEM.
3637 C  DETERM(A) = IP(N)*A(1,1)*A(2,2)*...*A(N,N).
3638 C  IF IP(N)=O, A IS SINGULAR, SOL WILL DIVIDE BY ZERO.
3639 C
3640 C  REFERENCE..
3641 C     THIS IS A SLIGHT MODIFICATION OF
3642 C     C. B. MOLER, ALGORITHM 423, LINEAR EQUATION KppSolveR,
3643 C     C.A.C.M. 15 (1972), P. 274.
3644 C-----------------------------------------------------------------------
3645       IER = 0
3646       IP(N) = 1
3647       IF (N .EQ. 1) GO TO 70
3648       NM1 = N - 1
3649       DO 60 K = 1,NM1
3650         KP1 = K + 1
3651         M = K
3652         NA = MIN0(N,LB+K)
3653         DO 10 I = KP1,NA
3654           IF (DABS(A(I,K)) .GT. DABS(A(M,K))) M = I
3655  10     CONTINUE
3656         IP(K) = M
3657         T = A(M,K)
3658         IF (M .EQ. K) GO TO 20
3659         IP(N) = -IP(N)
3660         A(M,K) = A(K,K)
3661         A(K,K) = T
3662  20     CONTINUE
3663         IF (T .EQ. 0.D0) GO TO 80
3664         T = 1.D0/T
3665         DO 30 I = KP1,NA
3666  30       A(I,K) = -A(I,K)*T
3667         DO 50 J = KP1,N
3668           T = A(M,J)
3669           A(M,J) = A(K,J)
3670           A(K,J) = T
3671           IF (T .EQ. 0.D0) GO TO 45
3672           DO 40 I = KP1,NA
3673  40         A(I,J) = A(I,J) + A(I,K)*T
3674  45       CONTINUE
3675  50       CONTINUE
3676  60     CONTINUE
3677  70   K = N
3678       IF (A(N,N) .EQ. 0.D0) GO TO 80
3679       RETURN
3680  80   IER = K
3681       IP(N) = 0
3682       RETURN
3683 C----------------------- END OF SUBROUTINE DECH ------------------------
3684       END
3685 C
3686 C
3687       SUBROUTINE SOLH (N, NDIM, A, LB, B, IP)
3688 C VERSION REAL KPP_REAL
3689       INTEGER N,NDIM,IP,NM1,K,KP1,M,I,KB,KM1,LB,NA
3690       KPP_REAL A,B,T
3691       DIMENSION A(NDIM,N), B(N), IP(N)
3692 C-----------------------------------------------------------------------
3693 C  SOLUTION OF LINEAR SYSTEM, A*X = B .
3694 C  INPUT..
3695 C    N = ORDER OF MATRIX A.
3696 C    NDIM = DECLARED DIMENSION OF ARRAY  A .
3697 C    A = TRIANGULARIZED MATRIX OBTAINED FROM DECH.
3698 C    LB = LOWER BANDWIDTH OF A.
3699 C    B = RIGHT HAND SIDE VECTOR.
3700 C    IP = PIVOT VECTOR OBTAINED FROM DEC.
3701 C  DO NOT USE IF DECH HAS SET IER .NE. 0.
3702 C  OUTPUT..
3703 C    B = SOLUTION VECTOR, X .
3704 C-----------------------------------------------------------------------
3705       IF (N .EQ. 1) GO TO 50
3706       NM1 = N - 1
3707       DO 20 K = 1,NM1
3708         KP1 = K + 1
3709         M = IP(K)
3710         T = B(M)
3711         B(M) = B(K)
3712         B(K) = T
3713         NA = MIN0(N,LB+K)
3714         DO 10 I = KP1,NA
3715  10       B(I) = B(I) + A(I,K)*T
3716  20     CONTINUE
3717       DO 40 KB = 1,NM1
3718         KM1 = N - KB
3719         K = KM1 + 1
3720         B(K) = B(K)/A(K,K)
3721         T = -B(K)
3722         DO 30 I = 1,KM1
3723  30       B(I) = B(I) + A(I,K)*T
3724  40     CONTINUE
3725  50   B(1) = B(1)/A(1,1)
3726       RETURN
3727 C----------------------- END OF SUBROUTINE SOLH ------------------------
3728       END
3729 C
3730       SUBROUTINE DECC (N, NDIM, AR, AI, IP, IER)
3731 C VERSION COMPLEX KPP_REAL
3732       IMPLICIT KPP_REAL (A-H,O-Z)
3733       INTEGER N,NDIM,IP,IER,NM1,K,KP1,M,I,J
3734       DIMENSION AR(NDIM,N), AI(NDIM,N), IP(N)
3735 C-----------------------------------------------------------------------
3736 C  MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION
3737 C  ------ MODIFICATION FOR COMPLEX MATRICES --------
3738 C  INPUT..
3739 C     N = ORDER OF MATRIX.
3740 C     NDIM = DECLARED DIMENSION OF ARRAYS  AR AND AI .
3741 C     (AR, AI) = MATRIX TO BE TRIANGULARIZED.
3742 C  OUTPUT..
3743 C     AR(I,J), I.LE.J = UPPER TRIANGULAR FACTOR, U ; REAL PART.
3744 C     AI(I,J), I.LE.J = UPPER TRIANGULAR FACTOR, U ; IMAGINARY PART.
3745 C     AR(I,J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I - L.
3746 C                                                    REAL PART.
3747 C     AI(I,J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I - L.
3748 C                                                    IMAGINARY PART.
3749 C     IP(K), K.LT.N = INDEX OF K-TH PIVOT ROW.
3750 C     IP(N) = (-1)**(NUMBER OF INTERCHANGES) OR O .
3751 C     IER = 0 IF MATRIX A IS NONSINGULAR, OR K IF FOUND TO BE
3752 C           SINGULAR AT STAGE K.
3753 C  USE  SOL  TO OBTAIN SOLUTION OF LINEAR SYSTEM.
3754 C  IF IP(N)=O, A IS SINGULAR, SOL WILL DIVIDE BY ZERO.
3755 C
3756 C  REFERENCE..
3757 C     C. B. MOLER, ALGORITHM 423, LINEAR EQUATION KppSolveR,
3758 C     C.A.C.M. 15 (1972), P. 274.
3759 C-----------------------------------------------------------------------
3760       IER = 0
3761       IP(N) = 1
3762       IF (N .EQ. 1) GO TO 70
3763       NM1 = N - 1
3764       DO 60 K = 1,NM1
3765         KP1 = K + 1
3766         M = K
3767         DO 10 I = KP1,N
3768           IF (DABS(AR(I,K))+DABS(AI(I,K)) .GT.
3769      &          DABS(AR(M,K))+DABS(AI(M,K))) M = I
3770  10     CONTINUE
3771         IP(K) = M
3772         TR = AR(M,K)
3773         TI = AI(M,K)
3774         IF (M .EQ. K) GO TO 20
3775         IP(N) = -IP(N)
3776         AR(M,K) = AR(K,K)
3777         AI(M,K) = AI(K,K)
3778         AR(K,K) = TR
3779         AI(K,K) = TI
3780  20     CONTINUE
3781         IF (DABS(TR)+DABS(TI) .EQ. 0.D0) GO TO 80
3782         DEN=TR*TR+TI*TI
3783         TR=TR/DEN
3784         TI=-TI/DEN
3785         DO 30 I = KP1,N
3786           PRODR=AR(I,K)*TR-AI(I,K)*TI
3787           PRODI=AI(I,K)*TR+AR(I,K)*TI
3788           AR(I,K)=-PRODR
3789           AI(I,K)=-PRODI
3790  30       CONTINUE
3791         DO 50 J = KP1,N
3792           TR = AR(M,J)
3793           TI = AI(M,J)
3794           AR(M,J) = AR(K,J)
3795           AI(M,J) = AI(K,J)
3796           AR(K,J) = TR
3797           AI(K,J) = TI
3798           IF (DABS(TR)+DABS(TI) .EQ. 0.D0) GO TO 48
3799           IF (TI .EQ. 0.D0) THEN
3800             DO 40 I = KP1,N
3801             PRODR=AR(I,K)*TR
3802             PRODI=AI(I,K)*TR
3803             AR(I,J) = AR(I,J) + PRODR
3804             AI(I,J) = AI(I,J) + PRODI
3805  40         CONTINUE
3806             GO TO 48
3807           END IF
3808           IF (TR .EQ. 0.D0) THEN
3809             DO 45 I = KP1,N
3810             PRODR=-AI(I,K)*TI
3811             PRODI=AR(I,K)*TI
3812             AR(I,J) = AR(I,J) + PRODR
3813             AI(I,J) = AI(I,J) + PRODI
3814  45         CONTINUE
3815             GO TO 48
3816           END IF
3817           DO 47 I = KP1,N
3818             PRODR=AR(I,K)*TR-AI(I,K)*TI
3819             PRODI=AI(I,K)*TR+AR(I,K)*TI
3820             AR(I,J) = AR(I,J) + PRODR
3821             AI(I,J) = AI(I,J) + PRODI
3822  47         CONTINUE
3823  48       CONTINUE
3824  50       CONTINUE
3825  60     CONTINUE
3826  70   K = N
3827       IF (DABS(AR(N,N))+DABS(AI(N,N)) .EQ. 0.D0) GO TO 80
3828       RETURN
3829  80   IER = K
3830       IP(N) = 0
3831       RETURN
3832 C----------------------- END OF SUBROUTINE DECC ------------------------
3833       END
3834 C
3835 C
3836       SUBROUTINE SOLC (N, NDIM, AR, AI, BR, BI, IP)
3837 C VERSION COMPLEX KPP_REAL
3838       IMPLICIT KPP_REAL (A-H,O-Z)
3839       INTEGER N,NDIM,IP,NM1,K,KP1,M,I,KB,KM1
3840       DIMENSION AR(NDIM,N), AI(NDIM,N), BR(N), BI(N), IP(N)
3841 C-----------------------------------------------------------------------
3842 C  SOLUTION OF LINEAR SYSTEM, A*X = B .
3843 C  INPUT..
3844 C    N = ORDER OF MATRIX.
3845 C    NDIM = DECLARED DIMENSION OF ARRAYS  AR AND AI.
3846 C    (AR,AI) = TRIANGULARIZED MATRIX OBTAINED FROM DEC.
3847 C    (BR,BI) = RIGHT HAND SIDE VECTOR.
3848 C    IP = PIVOT VECTOR OBTAINED FROM DEC.
3849 C  DO NOT USE IF DEC HAS SET IER .NE. 0.
3850 C  OUTPUT..
3851 C    (BR,BI) = SOLUTION VECTOR, X .
3852 C-----------------------------------------------------------------------
3853       IF (N .EQ. 1) GO TO 50
3854       NM1 = N - 1
3855       DO 20 K = 1,NM1
3856         KP1 = K + 1
3857         M = IP(K)
3858         TR = BR(M)
3859         TI = BI(M)
3860         BR(M) = BR(K)
3861         BI(M) = BI(K)
3862         BR(K) = TR
3863         BI(K) = TI
3864         DO 10 I = KP1,N
3865           PRODR=AR(I,K)*TR-AI(I,K)*TI
3866           PRODI=AI(I,K)*TR+AR(I,K)*TI
3867           BR(I) = BR(I) + PRODR
3868           BI(I) = BI(I) + PRODI
3869  10       CONTINUE
3870  20     CONTINUE
3871       DO 40 KB = 1,NM1
3872         KM1 = N - KB
3873         K = KM1 + 1
3874         DEN=AR(K,K)*AR(K,K)+AI(K,K)*AI(K,K)
3875         PRODR=BR(K)*AR(K,K)+BI(K)*AI(K,K)
3876         PRODI=BI(K)*AR(K,K)-BR(K)*AI(K,K)
3877         BR(K)=PRODR/DEN
3878         BI(K)=PRODI/DEN
3879         TR = -BR(K)
3880         TI = -BI(K)
3881         DO 30 I = 1,KM1
3882           PRODR=AR(I,K)*TR-AI(I,K)*TI
3883           PRODI=AI(I,K)*TR+AR(I,K)*TI
3884           BR(I) = BR(I) + PRODR
3885           BI(I) = BI(I) + PRODI
3886  30       CONTINUE
3887  40     CONTINUE
3888  50     CONTINUE
3889         DEN=AR(1,1)*AR(1,1)+AI(1,1)*AI(1,1)
3890         PRODR=BR(1)*AR(1,1)+BI(1)*AI(1,1)
3891         PRODI=BI(1)*AR(1,1)-BR(1)*AI(1,1)
3892         BR(1)=PRODR/DEN
3893         BI(1)=PRODI/DEN
3894       RETURN
3895 C----------------------- END OF SUBROUTINE SOLC ------------------------
3896       END
3897 C
3898 C
3899       SUBROUTINE DECHC (N, NDIM, AR, AI, LB, IP, IER)
3900 C VERSION COMPLEX KPP_REAL
3901       IMPLICIT KPP_REAL (A-H,O-Z)
3902       INTEGER N,NDIM,IP,IER,NM1,K,KP1,M,I,J
3903       DIMENSION AR(NDIM,N), AI(NDIM,N), IP(N)
3904 C-----------------------------------------------------------------------
3905 C  MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION
3906 C  ------ MODIFICATION FOR COMPLEX MATRICES --------
3907 C  INPUT..
3908 C     N = ORDER OF MATRIX.
3909 C     NDIM = DECLARED DIMENSION OF ARRAYS  AR AND AI .
3910 C     (AR, AI) = MATRIX TO BE TRIANGULARIZED.
3911 C  OUTPUT..
3912 C     AR(I,J), I.LE.J = UPPER TRIANGULAR FACTOR, U ; REAL PART.
3913 C     AI(I,J), I.LE.J = UPPER TRIANGULAR FACTOR, U ; IMAGINARY PART.
3914 C     AR(I,J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I - L.
3915 C                                                    REAL PART.
3916 C     AI(I,J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I - L.
3917 C                                                    IMAGINARY PART.
3918 C     LB = LOWER BANDWIDTH OF A (DIAGONAL NOT COUNTED), LB.GE.1.
3919 C     IP(K), K.LT.N = INDEX OF K-TH PIVOT ROW.
3920 C     IP(N) = (-1)**(NUMBER OF INTERCHANGES) OR O .
3921 C     IER = 0 IF MATRIX A IS NONSINGULAR, OR K IF FOUND TO BE
3922 C           SINGULAR AT STAGE K.
3923 C  USE  SOL  TO OBTAIN SOLUTION OF LINEAR SYSTEM.
3924 C  IF IP(N)=O, A IS SINGULAR, SOL WILL DIVIDE BY ZERO.
3925 C
3926 C  REFERENCE..
3927 C     C. B. MOLER, ALGORITHM 423, LINEAR EQUATION KppSolveR,
3928 C     C.A.C.M. 15 (1972), P. 274.
3929 C-----------------------------------------------------------------------
3930       IER = 0
3931       IP(N) = 1
3932       IF (LB .EQ. 0) GO TO 70
3933       IF (N .EQ. 1) GO TO 70
3934       NM1 = N - 1
3935       DO 60 K = 1,NM1
3936         KP1 = K + 1
3937         M = K
3938         NA = MIN0(N,LB+K)
3939         DO 10 I = KP1,NA
3940           IF (DABS(AR(I,K))+DABS(AI(I,K)) .GT.
3941      &          DABS(AR(M,K))+DABS(AI(M,K))) M = I
3942  10     CONTINUE
3943         IP(K) = M
3944         TR = AR(M,K)
3945         TI = AI(M,K)
3946         IF (M .EQ. K) GO TO 20
3947         IP(N) = -IP(N)
3948         AR(M,K) = AR(K,K)
3949         AI(M,K) = AI(K,K)
3950         AR(K,K) = TR
3951         AI(K,K) = TI
3952  20     CONTINUE
3953         IF (DABS(TR)+DABS(TI) .EQ. 0.D0) GO TO 80
3954         DEN=TR*TR+TI*TI
3955         TR=TR/DEN
3956         TI=-TI/DEN
3957         DO 30 I = KP1,NA
3958           PRODR=AR(I,K)*TR-AI(I,K)*TI
3959           PRODI=AI(I,K)*TR+AR(I,K)*TI
3960           AR(I,K)=-PRODR
3961           AI(I,K)=-PRODI
3962  30       CONTINUE
3963         DO 50 J = KP1,N
3964           TR = AR(M,J)
3965           TI = AI(M,J)
3966           AR(M,J) = AR(K,J)
3967           AI(M,J) = AI(K,J)
3968           AR(K,J) = TR
3969           AI(K,J) = TI
3970           IF (DABS(TR)+DABS(TI) .EQ. 0.D0) GO TO 48
3971           IF (TI .EQ. 0.D0) THEN
3972             DO 40 I = KP1,NA
3973             PRODR=AR(I,K)*TR
3974             PRODI=AI(I,K)*TR
3975             AR(I,J) = AR(I,J) + PRODR
3976             AI(I,J) = AI(I,J) + PRODI
3977  40         CONTINUE
3978             GO TO 48
3979           END IF
3980           IF (TR .EQ. 0.D0) THEN
3981             DO 45 I = KP1,NA
3982             PRODR=-AI(I,K)*TI
3983             PRODI=AR(I,K)*TI
3984             AR(I,J) = AR(I,J) + PRODR
3985             AI(I,J) = AI(I,J) + PRODI
3986  45         CONTINUE
3987             GO TO 48
3988           END IF
3989           DO 47 I = KP1,NA
3990             PRODR=AR(I,K)*TR-AI(I,K)*TI
3991             PRODI=AI(I,K)*TR+AR(I,K)*TI
3992             AR(I,J) = AR(I,J) + PRODR
3993             AI(I,J) = AI(I,J) + PRODI
3994  47         CONTINUE
3995  48       CONTINUE
3996  50       CONTINUE
3997  60     CONTINUE
3998  70   K = N
3999       IF (DABS(AR(N,N))+DABS(AI(N,N)) .EQ. 0.D0) GO TO 80
4000       RETURN
4001  80   IER = K
4002       IP(N) = 0
4003       RETURN
4004 C----------------------- END OF SUBROUTINE DECHC -----------------------
4005       END
4006 C
4007 C
4008       SUBROUTINE SOLHC (N, NDIM, AR, AI, LB, BR, BI, IP)
4009 C VERSION COMPLEX KPP_REAL
4010       IMPLICIT KPP_REAL (A-H,O-Z)
4011       INTEGER N,NDIM,IP,NM1,K,KP1,M,I,KB,KM1
4012       DIMENSION AR(NDIM,N), AI(NDIM,N), BR(N), BI(N), IP(N)
4013 C-----------------------------------------------------------------------
4014 C  SOLUTION OF LINEAR SYSTEM, A*X = B .
4015 C  INPUT..
4016 C    N = ORDER OF MATRIX.
4017 C    NDIM = DECLARED DIMENSION OF ARRAYS  AR AND AI.
4018 C    (AR,AI) = TRIANGULARIZED MATRIX OBTAINED FROM DEC.
4019 C    (BR,BI) = RIGHT HAND SIDE VECTOR.
4020 C    LB = LOWER BANDWIDTH OF A.
4021 C    IP = PIVOT VECTOR OBTAINED FROM DEC.
4022 C  DO NOT USE IF DEC HAS SET IER .NE. 0.
4023 C  OUTPUT..
4024 C    (BR,BI) = SOLUTION VECTOR, X .
4025 C-----------------------------------------------------------------------
4026       IF (N .EQ. 1) GO TO 50
4027       NM1 = N - 1
4028       IF (LB .EQ. 0) GO TO 25
4029       DO 20 K = 1,NM1
4030         KP1 = K + 1
4031         M = IP(K)
4032         TR = BR(M)
4033         TI = BI(M)
4034         BR(M) = BR(K)
4035         BI(M) = BI(K)
4036         BR(K) = TR
4037         BI(K) = TI
4038         DO 10 I = KP1,MIN0(N,LB+K)
4039           PRODR=AR(I,K)*TR-AI(I,K)*TI
4040           PRODI=AI(I,K)*TR+AR(I,K)*TI
4041           BR(I) = BR(I) + PRODR
4042           BI(I) = BI(I) + PRODI
4043  10       CONTINUE
4044  20     CONTINUE
4045  25     CONTINUE
4046       DO 40 KB = 1,NM1
4047         KM1 = N - KB
4048         K = KM1 + 1
4049         DEN=AR(K,K)*AR(K,K)+AI(K,K)*AI(K,K)
4050         PRODR=BR(K)*AR(K,K)+BI(K)*AI(K,K)
4051         PRODI=BI(K)*AR(K,K)-BR(K)*AI(K,K)
4052         BR(K)=PRODR/DEN
4053         BI(K)=PRODI/DEN
4054         TR = -BR(K)
4055         TI = -BI(K)
4056         DO 30 I = 1,KM1
4057           PRODR=AR(I,K)*TR-AI(I,K)*TI
4058           PRODI=AI(I,K)*TR+AR(I,K)*TI
4059           BR(I) = BR(I) + PRODR
4060           BI(I) = BI(I) + PRODI
4061  30       CONTINUE
4062  40     CONTINUE
4063  50     CONTINUE
4064         DEN=AR(1,1)*AR(1,1)+AI(1,1)*AI(1,1)
4065         PRODR=BR(1)*AR(1,1)+BI(1)*AI(1,1)
4066         PRODI=BI(1)*AR(1,1)-BR(1)*AI(1,1)
4067         BR(1)=PRODR/DEN
4068         BI(1)=PRODI/DEN
4069       RETURN
4070 C----------------------- END OF SUBROUTINE SOLHC -----------------------
4071       END
4072 C
4073       SUBROUTINE DECB (N, NDIM, A, ML, MU, IP, IER)
4074       KPP_REAL A,T
4075       DIMENSION A(NDIM,N), IP(N)
4076 C-----------------------------------------------------------------------
4077 C  MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION OF A BANDED
4078 C  MATRIX WITH LOWER BANDWIDTH ML AND UPPER BANDWIDTH MU
4079 C  INPUT..
4080 C     N       ORDER OF THE ORIGINAL MATRIX A.
4081 C     NDIM    DECLARED DIMENSION OF ARRAY  A.
4082 C     A       CONTAINS THE MATRIX IN BAND STORAGE.   THE COLUMNS
4083 C                OF THE MATRIX ARE STORED IN THE COLUMNS OF  A  AND
4084 C                THE DIAGONALS OF THE MATRIX ARE STORED IN ROWS
4085 C                ML+1 THROUGH 2*ML+MU+1 OF  A.
4086 C     ML      LOWER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED).
4087 C     MU      UPPER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED).
4088 C  OUTPUT..
4089 C     A       AN UPPER TRIANGULAR MATRIX IN BAND STORAGE AND
4090 C                THE MULTIPLIERS WHICH WERE USED TO OBTAIN IT.
4091 C     IP      INDEX VECTOR OF PIVOT INDICES.
4092 C     IP(N)   (-1)**(NUMBER OF INTERCHANGES) OR O .
4093 C     IER     = 0 IF MATRIX A IS NONSINGULAR, OR  = K IF FOUND TO BE
4094 C                SINGULAR AT STAGE K.
4095 C  USE  SOLB  TO OBTAIN SOLUTION OF LINEAR SYSTEM.
4096 C  DETERM(A) = IP(N)*A(MD,1)*A(MD,2)*...*A(MD,N)  WITH MD=ML+MU+1.
4097 C  IF IP(N)=O, A IS SINGULAR, SOLB WILL DIVIDE BY ZERO.
4098 C
4099 C  REFERENCE..
4100 C     THIS IS A MODIFICATION OF
4101 C     C. B. MOLER, ALGORITHM 423, LINEAR EQUATION KppSolveR,
4102 C     C.A.C.M. 15 (1972), P. 274.
4103 C-----------------------------------------------------------------------
4104       IER = 0
4105       IP(N) = 1
4106       MD = ML + MU + 1
4107       MD1 = MD + 1
4108       JU = 0
4109       IF (ML .EQ. 0) GO TO 70
4110       IF (N .EQ. 1) GO TO 70
4111       IF (N .LT. MU+2) GO TO 7
4112       DO 5 J = MU+2,N
4113       DO 5 I = 1,ML
4114   5   A(I,J) = 0.D0
4115   7   NM1 = N - 1
4116       DO 60 K = 1,NM1
4117         KP1 = K + 1
4118         M = MD
4119         MDL = MIN(ML,N-K) + MD
4120         DO 10 I = MD1,MDL
4121           IF (DABS(A(I,K)) .GT. DABS(A(M,K))) M = I
4122  10     CONTINUE
4123         IP(K) = M + K - MD
4124         T = A(M,K)
4125         IF (M .EQ. MD) GO TO 20
4126         IP(N) = -IP(N)
4127         A(M,K) = A(MD,K)
4128         A(MD,K) = T
4129  20     CONTINUE
4130         IF (T .EQ. 0.D0) GO TO 80
4131         T = 1.D0/T
4132         DO 30 I = MD1,MDL
4133  30       A(I,K) = -A(I,K)*T
4134         JU = MIN0(MAX0(JU,MU+IP(K)),N)
4135         MM = MD
4136         IF (JU .LT. KP1) GO TO 55
4137         DO 50 J = KP1,JU
4138           M = M - 1
4139           MM = MM - 1
4140           T = A(M,J)
4141           IF (M .EQ. MM) GO TO 35
4142           A(M,J) = A(MM,J)
4143           A(MM,J) = T
4144  35       CONTINUE
4145           IF (T .EQ. 0.D0) GO TO 45
4146           JK = J - K
4147           DO 40 I = MD1,MDL
4148             IJK = I - JK
4149  40         A(IJK,J) = A(IJK,J) + A(I,K)*T
4150  45       CONTINUE
4151  50       CONTINUE
4152  55     CONTINUE
4153  60     CONTINUE
4154  70   K = N
4155       IF (A(MD,N) .EQ. 0.D0) GO TO 80
4156       RETURN
4157  80   IER = K
4158       IP(N) = 0
4159       RETURN
4160 C----------------------- END OF SUBROUTINE DECB ------------------------
4161       END
4162 C
4163 C
4164       SUBROUTINE SOLB (N, NDIM, A, ML, MU, B, IP)
4165       KPP_REAL A,B,T
4166       DIMENSION A(NDIM,N), B(N), IP(N)
4167 C-----------------------------------------------------------------------
4168 C  SOLUTION OF LINEAR SYSTEM, A*X = B .
4169 C  INPUT..
4170 C    N      ORDER OF MATRIX A.
4171 C    NDIM   DECLARED DIMENSION OF ARRAY  A .
4172 C    A      TRIANGULARIZED MATRIX OBTAINED FROM DECB.
4173 C    ML     LOWER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED).
4174 C    MU     UPPER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED).
4175 C    B      RIGHT HAND SIDE VECTOR.
4176 C    IP     PIVOT VECTOR OBTAINED FROM DECB.
4177 C  DO NOT USE IF DECB HAS SET IER .NE. 0.
4178 C  OUTPUT..
4179 C    B      SOLUTION VECTOR, X .
4180 C-----------------------------------------------------------------------
4181       MD = ML + MU + 1
4182       MD1 = MD + 1
4183       MDM = MD - 1
4184       NM1 = N - 1
4185       IF (ML .EQ. 0) GO TO 25
4186       IF (N .EQ. 1) GO TO 50
4187       DO 20 K = 1,NM1
4188         M = IP(K)
4189         T = B(M)
4190         B(M) = B(K)
4191         B(K) = T
4192         MDL = MIN(ML,N-K) + MD
4193         DO 10 I = MD1,MDL
4194           IMD = I + K - MD
4195  10       B(IMD) = B(IMD) + A(I,K)*T
4196  20     CONTINUE
4197  25   CONTINUE
4198       DO 40 KB = 1,NM1
4199         K = N + 1 - KB
4200         B(K) = B(K)/A(MD,K)
4201         T = -B(K)
4202         KMD = MD - K
4203         LM = MAX0(1,KMD+1)
4204         DO 30 I = LM,MDM
4205           IMD = I - KMD
4206  30       B(IMD) = B(IMD) + A(I,K)*T
4207  40     CONTINUE
4208  50   B(1) = B(1)/A(MD,1)
4209       RETURN
4210 C----------------------- END OF SUBROUTINE SOLB ------------------------
4211       END
4212 C
4213       SUBROUTINE DECBC (N, NDIM, AR, AI, ML, MU, IP, IER)
4214       IMPLICIT KPP_REAL (A-H,O-Z)
4215       DIMENSION AR(NDIM,N), AI(NDIM,N), IP(N)
4216 C-----------------------------------------------------------------------
4217 C  MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION OF A BANDED COMPLEX
4218 C  MATRIX WITH LOWER BANDWIDTH ML AND UPPER BANDWIDTH MU
4219 C  INPUT..
4220 C     N       ORDER OF THE ORIGINAL MATRIX A.
4221 C     NDIM    DECLARED DIMENSION OF ARRAY  A.
4222 C     AR, AI     CONTAINS THE MATRIX IN BAND STORAGE.   THE COLUMNS
4223 C                OF THE MATRIX ARE STORED IN THE COLUMNS OF  AR (REAL
4224 C                PART) AND AI (IMAGINARY PART)  AND
4225 C                THE DIAGONALS OF THE MATRIX ARE STORED IN ROWS
4226 C                ML+1 THROUGH 2*ML+MU+1 OF  AR AND AI.
4227 C     ML      LOWER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED).
4228 C     MU      UPPER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED).
4229 C  OUTPUT..
4230 C     AR, AI  AN UPPER TRIANGULAR MATRIX IN BAND STORAGE AND
4231 C                THE MULTIPLIERS WHICH WERE USED TO OBTAIN IT.
4232 C     IP      INDEX VECTOR OF PIVOT INDICES.
4233 C     IP(N)   (-1)**(NUMBER OF INTERCHANGES) OR O .
4234 C     IER     = 0 IF MATRIX A IS NONSINGULAR, OR  = K IF FOUND TO BE
4235 C                SINGULAR AT STAGE K.
4236 C  USE  SOLBC  TO OBTAIN SOLUTION OF LINEAR SYSTEM.
4237 C  DETERM(A) = IP(N)*A(MD,1)*A(MD,2)*...*A(MD,N)  WITH MD=ML+MU+1.
4238 C  IF IP(N)=O, A IS SINGULAR, SOLBC WILL DIVIDE BY ZERO.
4239 C
4240 C  REFERENCE..
4241 C     THIS IS A MODIFICATION OF
4242 C     C. B. MOLER, ALGORITHM 423, LINEAR EQUATION KppSolveR,
4243 C     C.A.C.M. 15 (1972), P. 274.
4244 C-----------------------------------------------------------------------
4245       IER = 0
4246       IP(N) = 1
4247       MD = ML + MU + 1
4248       MD1 = MD + 1
4249       JU = 0
4250       IF (ML .EQ. 0) GO TO 70
4251       IF (N .EQ. 1) GO TO 70
4252       IF (N .LT. MU+2) GO TO 7
4253       DO 5 J = MU+2,N
4254       DO 5 I = 1,ML
4255       AR(I,J) = 0.D0
4256       AI(I,J) = 0.D0
4257   5   CONTINUE
4258   7   NM1 = N - 1
4259       DO 60 K = 1,NM1
4260         KP1 = K + 1
4261         M = MD
4262         MDL = MIN(ML,N-K) + MD
4263         DO 10 I = MD1,MDL
4264           IF (DABS(AR(I,K))+DABS(AI(I,K)) .GT.
4265      &          DABS(AR(M,K))+DABS(AI(M,K))) M = I
4266  10     CONTINUE
4267         IP(K) = M + K - MD
4268         TR = AR(M,K)
4269         TI = AI(M,K)
4270         IF (M .EQ. MD) GO TO 20
4271         IP(N) = -IP(N)
4272         AR(M,K) = AR(MD,K)
4273         AI(M,K) = AI(MD,K)
4274         AR(MD,K) = TR
4275         AI(MD,K) = TI
4276  20     IF (DABS(TR)+DABS(TI) .EQ. 0.D0) GO TO 80
4277         DEN=TR*TR+TI*TI
4278         TR=TR/DEN
4279         TI=-TI/DEN
4280         DO 30 I = MD1,MDL
4281           PRODR=AR(I,K)*TR-AI(I,K)*TI
4282           PRODI=AI(I,K)*TR+AR(I,K)*TI
4283           AR(I,K)=-PRODR
4284           AI(I,K)=-PRODI
4285  30       CONTINUE
4286         JU = MIN0(MAX0(JU,MU+IP(K)),N)
4287         MM = MD
4288         IF (JU .LT. KP1) GO TO 55
4289         DO 50 J = KP1,JU
4290           M = M - 1
4291           MM = MM - 1
4292           TR = AR(M,J)
4293           TI = AI(M,J)
4294           IF (M .EQ. MM) GO TO 35
4295           AR(M,J) = AR(MM,J)
4296           AI(M,J) = AI(MM,J)
4297           AR(MM,J) = TR
4298           AI(MM,J) = TI
4299  35       CONTINUE
4300           IF (DABS(TR)+DABS(TI) .EQ. 0.D0) GO TO 48
4301           JK = J - K
4302           IF (TI .EQ. 0.D0) THEN
4303             DO 40 I = MD1,MDL
4304             IJK = I - JK
4305             PRODR=AR(I,K)*TR
4306             PRODI=AI(I,K)*TR
4307             AR(IJK,J) = AR(IJK,J) + PRODR
4308             AI(IJK,J) = AI(IJK,J) + PRODI
4309  40         CONTINUE
4310             GO TO 48
4311           END IF
4312           IF (TR .EQ. 0.D0) THEN
4313             DO 45 I = MD1,MDL
4314             IJK = I - JK
4315             PRODR=-AI(I,K)*TI
4316             PRODI=AR(I,K)*TI
4317             AR(IJK,J) = AR(IJK,J) + PRODR
4318             AI(IJK,J) = AI(IJK,J) + PRODI
4319  45         CONTINUE
4320             GO TO 48
4321           END IF
4322           DO 47 I = MD1,MDL
4323             IJK = I - JK
4324             PRODR=AR(I,K)*TR-AI(I,K)*TI
4325             PRODI=AI(I,K)*TR+AR(I,K)*TI
4326             AR(IJK,J) = AR(IJK,J) + PRODR
4327             AI(IJK,J) = AI(IJK,J) + PRODI
4328  47         CONTINUE
4329  48       CONTINUE
4330  50       CONTINUE
4331  55     CONTINUE
4332  60     CONTINUE
4333  70   K = N
4334       IF (DABS(AR(MD,N))+DABS(AI(MD,N)) .EQ. 0.D0) GO TO 80
4335       RETURN
4336  80   IER = K
4337       IP(N) = 0
4338       RETURN
4339 C----------------------- END OF SUBROUTINE DECBC ------------------------
4340       END
4341 C
4342 C
4343       SUBROUTINE SOLBC (N, NDIM, AR, AI, ML, MU, BR, BI, IP)
4344       IMPLICIT KPP_REAL (A-H,O-Z)
4345       DIMENSION AR(NDIM,N), AI(NDIM,N), BR(N), BI(N), IP(N)
4346 C-----------------------------------------------------------------------
4347 C  SOLUTION OF LINEAR SYSTEM, A*X = B ,
4348 C                  VERSION BANDED AND COMPLEX-KPP_REAL.
4349 C  INPUT..
4350 C    N      ORDER OF MATRIX A.
4351 C    NDIM   DECLARED DIMENSION OF ARRAY  A .
4352 C    AR, AI TRIANGULARIZED MATRIX OBTAINED FROM DECB (REAL AND IMAG. PART).
4353 C    ML     LOWER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED).
4354 C    MU     UPPER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED).
4355 C    BR, BI RIGHT HAND SIDE VECTOR (REAL AND IMAG. PART).
4356 C    IP     PIVOT VECTOR OBTAINED FROM DECBC.
4357 C  DO NOT USE IF DECB HAS SET IER .NE. 0.
4358 C  OUTPUT..
4359 C    BR, BI SOLUTION VECTOR, X (REAL AND IMAG. PART).
4360 C-----------------------------------------------------------------------
4361       MD = ML + MU + 1
4362       MD1 = MD + 1
4363       MDM = MD - 1
4364       NM1 = N - 1
4365       IF (ML .EQ. 0) GO TO 25
4366       IF (N .EQ. 1) GO TO 50
4367       DO 20 K = 1,NM1
4368         M = IP(K)
4369         TR = BR(M)
4370         TI = BI(M)
4371         BR(M) = BR(K)
4372         BI(M) = BI(K)
4373         BR(K) = TR
4374         BI(K) = TI
4375         MDL = MIN(ML,N-K) + MD
4376         DO 10 I = MD1,MDL
4377           IMD = I + K - MD
4378           PRODR=AR(I,K)*TR-AI(I,K)*TI
4379           PRODI=AI(I,K)*TR+AR(I,K)*TI
4380           BR(IMD) = BR(IMD) + PRODR
4381           BI(IMD) = BI(IMD) + PRODI
4382  10     CONTINUE
4383  20     CONTINUE
4384  25     CONTINUE
4385       DO 40 KB = 1,NM1
4386         K = N + 1 - KB
4387         DEN=AR(MD,K)*AR(MD,K)+AI(MD,K)*AI(MD,K)
4388         PRODR=BR(K)*AR(MD,K)+BI(K)*AI(MD,K)
4389         PRODI=BI(K)*AR(MD,K)-BR(K)*AI(MD,K)
4390         BR(K)=PRODR/DEN
4391         BI(K)=PRODI/DEN
4392         TR = -BR(K)
4393         TI = -BI(K)
4394         KMD = MD - K
4395         LM = MAX0(1,KMD+1)
4396         DO 30 I = LM,MDM
4397           IMD = I - KMD
4398           PRODR=AR(I,K)*TR-AI(I,K)*TI
4399           PRODI=AI(I,K)*TR+AR(I,K)*TI
4400           BR(IMD) = BR(IMD) + PRODR
4401           BI(IMD) = BI(IMD) + PRODI
4402  30       CONTINUE
4403  40     CONTINUE
4404         DEN=AR(MD,1)*AR(MD,1)+AI(MD,1)*AI(MD,1)
4405         PRODR=BR(1)*AR(MD,1)+BI(1)*AI(MD,1)
4406         PRODI=BI(1)*AR(MD,1)-BR(1)*AI(MD,1)
4407         BR(1)=PRODR/DEN
4408         BI(1)=PRODI/DEN
4409  50   CONTINUE
4410       RETURN
4411 C----------------------- END OF SUBROUTINE SOLBC ------------------------
4412       END
4413 c
4414 C
4415       subroutine Elmhes(nm,n,low,igh,a,int)
4416 C
4417       integer i,j,m,n,la,nm,igh,kp1,low,mm1,mp1
4418       real*8 a(nm,n)
4419       real*8 x,y
4420       real*8 dabs
4421       integer int(igh)
4422 C
4423 C     this subroutine is a translation of the algol procedure elmhes,
4424 C     num. math. 12, 349-368(1968) by martin and wilkinson.
4425 C     handbook for auto. comp., vol.ii-linear algebra, 339-358(1971).
4426 C
4427 C     given a real general matrix, this subroutine
4428 C     reduces a submatrix situated in rows and columns
4429 C     low through igh to upper hessenberg form by
4430 C     stabilized elementary similarity transformations.
4431 C
4432 C     on input:
4433 C
4434 C      nm must be set to the row dimension of two-dimensional
4435 C        array parameters as declared in the calling program
4436 C        dimension statement;
4437 C
4438 C      n is the order of the matrix;
4439 C
4440 C      low and igh are integers determined by the balancing
4441 C        subroutine  balanc.      if  balanc  has not been used,
4442 C        set low=1, igh=n;
4443 C
4444 C      a contains the input matrix.
4445 C
4446 C     on output:
4447 C
4448 C      a contains the hessenberg matrix.  the multipliers
4449 C        which were used in the reduction are stored in the
4450 C        remaining triangle under the hessenberg matrix;
4451 C
4452 C      int contains information on the rows and columns
4453 C        interchanged in the reduction.
4454 C        only elements low through igh are used.
4455 C
4456 C     questions and comments should be directed to b. s. garbow,
4457 C     applied mathematics division, argonne national laboratory
4458 C
4459 C     ------------------------------------------------------------------
4460 C
4461       la = igh - 1
4462       kp1 = low + 1
4463       if (la .lt. kp1) go to 200
4464 C
4465       do 180 m = kp1, la
4466        mm1 = m - 1
4467        x = 0.0d0
4468        i = m
4469 C
4470        do 100 j = m, igh
4471           if (dabs(a(j,mm1)) .le. dabs(x)) go to 100
4472           x = a(j,mm1)
4473           i = j
4474   100   continue
4475 C
4476        int(m) = i
4477        if (i .eq. m) go to 130
4478 C    :::::::::: interchange rows and columns of a ::::::::::
4479        do 110 j = mm1, n
4480           y = a(i,j)
4481           a(i,j) = a(m,j)
4482           a(m,j) = y
4483   110   continue
4484 C
4485        do 120 j = 1, igh
4486           y = a(j,i)
4487           a(j,i) = a(j,m)
4488           a(j,m) = y
4489   120   continue
4490 C    :::::::::: end interchange ::::::::::
4491   130   if (x .eq. 0.0d0) go to 180
4492        mp1 = m + 1
4493 C
4494        do 160 i = mp1, igh
4495           y = a(i,mm1)
4496           if (y .eq. 0.0d0) go to 160
4497           y = y / x
4498           a(i,mm1) = y
4499 C
4500           do 140 j = m, n
4501   140      a(i,j) = a(i,j) - y * a(m,j)
4502 C
4503           do 150 j = 1, igh
4504   150      a(j,m) = a(j,m) + y * a(j,i)
4505 C
4506   160   continue
4507 C
4508   180 continue
4509 C
4510   200 return
4511 C    :::::::::: last card of elmhes ::::::::::
4512       end
4513
4514         SUBROUTINE SOLOUT (NR,XOLD,X,Y,CONT,LRC,N,RPAR,IPAR,IRTRN)
4515 C --- PRINTS SOLUTION AT EQUIDISTANT OUTPUT-POINTS BY USING "CONTR5"
4516         IMPLICIT KPP_REAL (A-H,O-Z)
4517         DIMENSION Y(N),CONT(LRC)
4518         COMMON /INTERN/XOUT
4519         IF (NR.EQ.1) THEN
4520            WRITE (6,99) X,Y(1),Y(2),NR-1
4521            XOUT=0.2D0
4522         ELSE
4523  10        CONTINUE
4524            IF (X.GE.XOUT) THEN
4525 C --- CONTINUOUS OUTPUT FOR RADAU5
4526               WRITE (6,99) XOUT,CONTR5(1,XOUT,CONT,LRC),
4527      &                     CONTR5(2,XOUT,CONT,LRC),NR-1
4528               XOUT=XOUT+0.2D0
4529               GOTO 10
4530            END IF
4531         END IF
4532  99     FORMAT(1X,'X =',F5.2,'    Y =',2E18.10,'    NSTEP =',I4)
4533         RETURN
4534         END
4535
4536         SUBROUTINE FUNC_CHEM (N, T, V, FCT, RPAR, IPAR)
4537         IMPLICIT NONE
4538         INCLUDE 'KPP_ROOT_Parameters.h'
4539         INCLUDE 'KPP_ROOT_Global.h'
4540         KPP_REAL V(NVAR), FCT(NVAR)
4541         KPP_REAL T, TOLD, RPAR(1)
4542         INTEGER N, IPAR(1)
4543         TOLD = TIME
4544         TIME = T
4545         CALL Update_SUN()
4546         CALL Update_RCONST()
4547         TIME = TOLD
4548         CALL Fun(V,  FIX, RCONST, FCT)
4549         RETURN
4550         END
4551
4552         SUBROUTINE JAC_CHEM (N, T, V, JV, NROWPD, RPAR, IPAR)
4553         IMPLICIT NONE
4554         INCLUDE 'KPP_ROOT_Parameters.h'
4555         INCLUDE 'KPP_ROOT_Global.h'
4556         KPP_REAL V(NVAR), JV(NVAR,NVAR)
4557         KPP_REAL T, TOLD, RPAR(1)
4558         INTEGER N, IPAR(1), NROWPD, i, j
4559         TOLD = TIME
4560         TIME = T
4561         CALL Update_SUN()
4562         CALL Update_RCONST()
4563         TIME = TOLD
4564         DO i=1,NVAR
4565            DO j=1,NVAR
4566               JV(i,j) = 0.D0
4567            END DO
4568         END DO
4569         CALL Jac(V,  FIX, RCONST, JV)
4570         RETURN
4571         END
4572