| 1 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! |
|---|
| 2 | ! SEULEX - Stiff extrapolation method based on linearly implicit Euler ! |
|---|
| 3 | ! By default the code employs the KPP sparse linear algebra routines ! |
|---|
| 4 | ! Compile with -DFULL_ALGEBRA to use full linear algebra (LAPACK) ! |
|---|
| 5 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~! |
|---|
| 6 | |
|---|
| 7 | MODULE KPP_ROOT_Integrator |
|---|
| 8 | |
|---|
| 9 | USE KPP_ROOT_Precision, ONLY: dp |
|---|
| 10 | USE KPP_ROOT_Parameters, ONLY: NVAR, LU_NONZERO |
|---|
| 11 | USE KPP_ROOT_Jacobian, ONLY: LU_DIAG |
|---|
| 12 | USE KPP_ROOT_LinearAlgebra |
|---|
| 13 | |
|---|
| 14 | IMPLICIT NONE |
|---|
| 15 | PUBLIC |
|---|
| 16 | SAVE |
|---|
| 17 | |
|---|
| 18 | ! variables from the former COMMON block /CONRA5/ are now here: |
|---|
| 19 | INTEGER :: NN, NN2, NN3, NN4 |
|---|
| 20 | KPP_REAL :: TSOL, HSOL |
|---|
| 21 | |
|---|
| 22 | ! Statistics |
|---|
| 23 | INTEGER :: Nfun,Njac,Nstp,Nacc,Nrej,Ndec,Nsol,Nsng |
|---|
| 24 | |
|---|
| 25 | ! Method parameters |
|---|
| 26 | |
|---|
| 27 | ! mz_rs_20050717: TODO: use strings of IERR_NAMES for error messages |
|---|
| 28 | ! description of the error numbers IERR |
|---|
| 29 | CHARACTER(LEN=50), PARAMETER, DIMENSION(-4:1) :: IERR_NAMES = (/ & |
|---|
| 30 | 'Matrix is repeatedly singular ', & ! -4 |
|---|
| 31 | 'Step size too small ', & ! -3 |
|---|
| 32 | 'More than Max_no_steps steps are needed ', & ! -2 |
|---|
| 33 | 'Insufficient storage for work or iwork ', & ! -1 |
|---|
| 34 | ' ', & ! 0 (not used) |
|---|
| 35 | 'Success ' /) ! 1 |
|---|
| 36 | |
|---|
| 37 | CONTAINS |
|---|
| 38 | |
|---|
| 39 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 40 | SUBROUTINE INTEGRATE( TIN, TOUT, & |
|---|
| 41 | ICNTRL_U, RCNTRL_U, ISTATUS_U, RSTATUS_U, IERR_U ) |
|---|
| 42 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 43 | |
|---|
| 44 | USE KPP_ROOT_Parameters, ONLY: nvar |
|---|
| 45 | USE KPP_ROOT_Global, ONLY: atol,rtol,var |
|---|
| 46 | |
|---|
| 47 | IMPLICIT NONE |
|---|
| 48 | |
|---|
| 49 | KPP_REAL :: TIN ! TIN - Start Time |
|---|
| 50 | KPP_REAL :: TOUT ! TOUT - End Time |
|---|
| 51 | INTEGER, INTENT(IN), OPTIONAL :: ICNTRL_U(20) |
|---|
| 52 | KPP_REAL, INTENT(IN), OPTIONAL :: RCNTRL_U(20) |
|---|
| 53 | INTEGER, INTENT(OUT), OPTIONAL :: ISTATUS_U(20) |
|---|
| 54 | KPP_REAL, INTENT(OUT), OPTIONAL :: RSTATUS_U(20) |
|---|
| 55 | INTEGER, INTENT(OUT), OPTIONAL :: IERR_U |
|---|
| 56 | |
|---|
| 57 | INTEGER :: Ncolumns, Ncolumns2, NRDENS |
|---|
| 58 | PARAMETER (Ncolumns=12,Ncolumns2=2+Ncolumns*(Ncolumns+3)/2,NRDENS=NVAR) |
|---|
| 59 | |
|---|
| 60 | KPP_REAL :: RCNTRL(20), RSTATUS(20) |
|---|
| 61 | INTEGER :: ICNTRL(20), ISTATUS(20) |
|---|
| 62 | INTEGER :: IERR |
|---|
| 63 | INTEGER, SAVE :: Ntotal = 0 |
|---|
| 64 | KPP_REAL, SAVE :: H |
|---|
| 65 | |
|---|
| 66 | H = 0.0_dp |
|---|
| 67 | |
|---|
| 68 | ICNTRL(1:20) = 0 |
|---|
| 69 | RCNTRL(1:20) = 0._dp |
|---|
| 70 | ICNTRL(10)=0 !~~~> OUTPUT ROUTINE IS NOT USED DURING INTEGRATION |
|---|
| 71 | |
|---|
| 72 | ! if optional parameters are given, and if they are >0, |
|---|
| 73 | ! they overwrite the default settings |
|---|
| 74 | IF (PRESENT(ICNTRL_U)) ICNTRL(:) = ICNTRL_U(:) |
|---|
| 75 | IF (PRESENT(RCNTRL_U)) RCNTRL(:) = RCNTRL_U(:) |
|---|
| 76 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~----- |
|---|
| 77 | |
|---|
| 78 | |
|---|
| 79 | CALL ATMSEULEX(NVAR,TIN,TOUT,VAR,H,RTOL,ATOL, & |
|---|
| 80 | RCNTRL,ICNTRL,RSTATUS,ISTATUS,IERR ) |
|---|
| 81 | Ntotal = Ntotal + Nstp |
|---|
| 82 | !!$ PRINT*,'NSTEPS=',Nstp,' (',Ntotal,') T=',TIN |
|---|
| 83 | |
|---|
| 84 | |
|---|
| 85 | Nfun = Nfun + ISTATUS(1) |
|---|
| 86 | Njac = Njac + ISTATUS(2) |
|---|
| 87 | Nstp = Nstp + ISTATUS(3) |
|---|
| 88 | Nacc = Nacc + ISTATUS(4) |
|---|
| 89 | Nrej = Nrej + ISTATUS(5) |
|---|
| 90 | Ndec = Ndec + ISTATUS(6) |
|---|
| 91 | Nsol = Nsol + ISTATUS(7) |
|---|
| 92 | |
|---|
| 93 | ! if optional parameters are given for output |
|---|
| 94 | ! use them to store information in them |
|---|
| 95 | IF (PRESENT(ISTATUS_U)) THEN |
|---|
| 96 | ISTATUS_U(:) = 0 |
|---|
| 97 | ISTATUS_U(1) = Nfun ! function calls |
|---|
| 98 | ISTATUS_U(2) = Njac ! jacobian calls |
|---|
| 99 | ISTATUS_U(3) = Nstp ! steps |
|---|
| 100 | ISTATUS_U(4) = Nacc ! accepted steps |
|---|
| 101 | ISTATUS_U(5) = Nrej ! rejected steps (except at the beginning) |
|---|
| 102 | ISTATUS_U(6) = Ndec ! LU decompositions |
|---|
| 103 | ISTATUS_U(7) = Nsol ! forward/backward substitutions |
|---|
| 104 | ENDIF |
|---|
| 105 | IF (PRESENT(RSTATUS_U)) THEN |
|---|
| 106 | RSTATUS_U(:) = 0. |
|---|
| 107 | RSTATUS_U(1) = TOUT ! final time |
|---|
| 108 | ENDIF |
|---|
| 109 | IF (PRESENT(IERR_U)) IERR_U = IERR |
|---|
| 110 | |
|---|
| 111 | ! mz_rs_20050716: IERR is returned to the user who then decides what to do |
|---|
| 112 | ! about it, i.e. either stop the run or ignore it. |
|---|
| 113 | !!$ IF (IERR < 0) THEN |
|---|
| 114 | !!$ PRINT *,'SEULEX: Unsuccessful exit at T=', TIN,' (IERR=',IERR,')' |
|---|
| 115 | !!$ STOP |
|---|
| 116 | !!$ ENDIF |
|---|
| 117 | |
|---|
| 118 | END SUBROUTINE INTEGRATE |
|---|
| 119 | |
|---|
| 120 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 121 | SUBROUTINE ATMSEULEX( N,Tinitial,Tfinal,Y,H,RelTol,AbsTol, & |
|---|
| 122 | RCNTRL,ICNTRL,RSTATUS,ISTATUS,IERR ) |
|---|
| 123 | |
|---|
| 124 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 125 | ! NUMERICAL SOLUTION OF A STIFF (OR DIFFERENTIAL ALGEBRAIC) |
|---|
| 126 | ! SYSTEM OF FIRST 0RDER ORDINARY DIFFERENTIAL EQUATIONS MY'=F(T,Y). |
|---|
| 127 | ! THIS IS AN EXTRAPOLATION-ALGORITHM, BASED ON THE |
|---|
| 128 | ! LINEARLY IMPLICIT EULER METHOD (WITH STEP SIZE CONTROL |
|---|
| 129 | ! AND ORDER SELECTION). |
|---|
| 130 | ! |
|---|
| 131 | ! AUTHORS: E. HAIRER AND G. WANNER |
|---|
| 132 | ! UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES |
|---|
| 133 | ! CH-1211 GENEVE 24, SWITZERLAND |
|---|
| 134 | ! E-MAIL: HAIRER@DIVSUN.UNIGE.CH, WANNER@DIVSUN.UNIGE.CH |
|---|
| 135 | ! INCLUSION OF DENSE OUTPUT BY E. HAIRER AND A. OSTERMANN |
|---|
| 136 | ! |
|---|
| 137 | ! THIS CODE IS PART OF THE BOOK: |
|---|
| 138 | ! E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL |
|---|
| 139 | ! EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS. |
|---|
| 140 | ! SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS 14, |
|---|
| 141 | ! SPRINGER-VERLAG (1991) |
|---|
| 142 | ! |
|---|
| 143 | ! VERSION OF SEPTEMBER 30, 1995 |
|---|
| 144 | ! |
|---|
| 145 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 146 | ! |
|---|
| 147 | ! INPUT PARAMETERS |
|---|
| 148 | ! ---------------- |
|---|
| 149 | ! N DIMENSION OF THE SYSTEM |
|---|
| 150 | ! |
|---|
| 151 | ! T INITIAL T-VALUE |
|---|
| 152 | ! |
|---|
| 153 | ! Y(N) INITIAL VALUES FOR Y |
|---|
| 154 | ! |
|---|
| 155 | ! Tend FINAL T-VALUE (Tend-T MAY BE POSITIVE OR NEGATIVE) |
|---|
| 156 | ! |
|---|
| 157 | ! H INITIAL STEP SIZE GUESS; |
|---|
| 158 | ! FOR STIFF EQUATIONS WITH INITIAL TRANSIENT, |
|---|
| 159 | ! H=1.D0/(NORM OF F'), USUALLY 1.D-2 OR 1.D-3, IS GOOD. |
|---|
| 160 | ! THIS CHOICE IS NOT VERY IMPORTANT, THE CODE QUICKLY |
|---|
| 161 | ! ADAPTS ITS STEP SIZE (IF H=0.D0, THE CODE PUTS H=1.D-6 |
|---|
| 162 | ! |
|---|
| 163 | ! RelTol,AbsTol RELATIVE AND ABSOLUTE ERROR TOLERANCES. THEY |
|---|
| 164 | ! CAN BE BOTH SCALARS OR ELSE BOTH VECTORS OF LENGTH N. |
|---|
| 165 | ! |
|---|
| 166 | ! JAC NAME (EXTERNAL) OF THE SUBROUTINE WHICH COMPUTES |
|---|
| 167 | ! THE PARTIAL DERIVATIVES OF F(T,Y) WITH RESPECT TO Y |
|---|
| 168 | ! |
|---|
| 169 | ! SOLOUT NAME (EXTERNAL) OF SUBROUTINE PROVIDING THE |
|---|
| 170 | ! NUMERICAL SOLUTION DURING INTEGRATION. |
|---|
| 171 | ! IF IOUT>=1, IT IS CALLED AFTER EVERY SUCCESSFUL STEP. |
|---|
| 172 | ! SUPPLY A DUMMY SUBROUTINE IF IOUT=0. |
|---|
| 173 | ! IT MUST HAVE THE FORM |
|---|
| 174 | ! SUBROUTINE SOLOUT (NR,TOLD,T,Y,RC,LRC,IC,LIC,N, |
|---|
| 175 | ! RPAR,IPAR,IRTRN) |
|---|
| 176 | ! KPP_REAL T,Y(N),RC(LRC),IC(LIC) |
|---|
| 177 | ! .... |
|---|
| 178 | ! SOLOUT FURNISHES THE SOLUTION "Y" AT THE NR-TH |
|---|
| 179 | ! GRID-POINT "T" (THEREBY THE INITIAL VALUE IS |
|---|
| 180 | ! THE FIRST GRID-POINT). |
|---|
| 181 | ! "TOLD" IS THE PRECEEDING GRID-POINT. |
|---|
| 182 | ! "IRTRN" SERVES TO INTERRUPT THE INTEGRATION. IF IRTRN |
|---|
| 183 | ! IS SET <0, SEULEX RETURNS TO THE CALLING PROGRAM. |
|---|
| 184 | ! DO NOT CHANGE THE ENTRIES OF RC(LRC),IC(LIC)! |
|---|
| 185 | ! |
|---|
| 186 | ! ----- CONTINUOUS OUTPUT (IF IOUT=2): ----- |
|---|
| 187 | ! DURING CALLS TO "SOLOUT", A CONTINUOUS SOLUTION |
|---|
| 188 | ! FOR THE INTERVAL [TOLD,T] IS AVAILABLE THROUGH |
|---|
| 189 | ! THE KPP_REAL FUNCTION |
|---|
| 190 | ! >>> CONTEX(I,S,RC,LRC,IC,LIC) <<< |
|---|
| 191 | ! WHICH PROVIDES AN APPROXIMATION TO THE I-TH |
|---|
| 192 | ! COMPONENT OF THE SOLUTION AT THE POINT S. THE VALUE |
|---|
| 193 | ! S SHOULD LIE IN THE INTERVAL [TOLD,T]. |
|---|
| 194 | ! |
|---|
| 195 | ! IOUT GIVES INFORMATION ON THE SUBROUTINE SOLOUT: |
|---|
| 196 | ! IOUT=0: SUBROUTINE IS NEVER CALLED |
|---|
| 197 | ! IOUT=1: SUBROUTINE IS USED FOR OUTPUT |
|---|
| 198 | ! IOUT=2: DENSE OUTPUT IS PERFORMED IN SOLOUT |
|---|
| 199 | ! |
|---|
| 200 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 201 | ! |
|---|
| 202 | ! SOPHISTICATED SETTING OF PARAMETERS |
|---|
| 203 | ! ----------------------------------- |
|---|
| 204 | ! SEVERAL PARAMETERS OF THE CODE ARE TUNED TO MAKE IT CNTRL |
|---|
| 205 | ! WELL. THEY MAY BE DEFINED BY SETTING CNTRL(1),..,CNTRL(13) |
|---|
| 206 | ! AS WELL AS ICNTRL(1),..,ICNTRL(4) DIFFERENT FROM ZERO. |
|---|
| 207 | ! FOR ZERO INPUT, THE CODE CHOOSES DEFAULT VALUES: |
|---|
| 208 | ! |
|---|
| 209 | ! |
|---|
| 210 | !~~~> INPUT PARAMETERS: |
|---|
| 211 | ! |
|---|
| 212 | ! Note: For input parameters equal to zero the default values of the |
|---|
| 213 | ! corresponding variables are used. |
|---|
| 214 | !~~~> |
|---|
| 215 | ! ICNTRL(1) = 1: F = F(y) Independent of T (autonomous) |
|---|
| 216 | ! = 0: F = F(t,y) Depends on T (non-autonomous) |
|---|
| 217 | ! |
|---|
| 218 | ! ICNTRL(2) = 0: AbsTol, RelTol are NVAR-dimensional vectors |
|---|
| 219 | ! = 1: AbsTol, RelTol are scalars |
|---|
| 220 | ! |
|---|
| 221 | ! ICNTRL(3) -> not used |
|---|
| 222 | ! |
|---|
| 223 | ! ICNTRL(4) -> maximum number of integration steps |
|---|
| 224 | ! For ICNTRL(4)=0 the default value of 100000 is used |
|---|
| 225 | ! |
|---|
| 226 | ! ICNTRL(11) THE MAXIMUM NUMBER OF COLUMNS IN THE EXTRAPOLATION |
|---|
| 227 | ! TABLE. THE DEFAULT VALUE (FOR ICNTRL(3)=0) IS 12. |
|---|
| 228 | ! IF ICNTRL(3).NE.0 THEN ICNTRL(3) SHOULD BE >= 3. |
|---|
| 229 | ! |
|---|
| 230 | ! ICNTRL(12) SWITCH FOR THE STEP SIZE SEQUENCE |
|---|
| 231 | ! IF ICNTRL(4) == 1 THEN 1,2,3,4,6,8,12,16,24,32,48,... |
|---|
| 232 | ! IF ICNTRL(4) == 2 THEN 2,3,4,6,8,12,16,24,32,48,64,... |
|---|
| 233 | ! IF ICNTRL(4) == 3 THEN 1,2,3,4,5,6,7,8,9,10,... |
|---|
| 234 | ! IF ICNTRL(4) == 4 THEN 2,3,4,5,6,7,8,9,10,11,... |
|---|
| 235 | ! THE DEFAULT VALUE (FOR ICNTRL(4)=0) IS ICNTRL(4)=2. |
|---|
| 236 | ! |
|---|
| 237 | ! ICNTRL(13) PARAMETER "LAMBDA" OF DENSE OUTPUT; POSSIBLE VALUES |
|---|
| 238 | ! ARE 0 AND 1; DEFAULT ICNTRL(5)=0. |
|---|
| 239 | ! |
|---|
| 240 | ! ICNTRL(14) = NRDENS = NUMBER OF COMPONENTS, FOR WHICH DENSE OUTPUT |
|---|
| 241 | ! IS REQUIRED |
|---|
| 242 | ! |
|---|
| 243 | ! ICNTRL(21),...,ICNTRL(NRDENS+20) INDICATE THE COMPONENTS, FOR WHICH |
|---|
| 244 | ! DENSE OUTPUT IS REQUIRED |
|---|
| 245 | ! |
|---|
| 246 | !~~~> Real parameters |
|---|
| 247 | ! |
|---|
| 248 | ! RCNTRL(1) -> Hmin, lower bound for the integration step size |
|---|
| 249 | ! It is strongly recommended to keep Hmin = ZERO |
|---|
| 250 | ! RCNTRL(2) -> Hmax, upper bound for the integration step size |
|---|
| 251 | ! RCNTRL(3) -> Hstart, starting value for the integration step size |
|---|
| 252 | ! RCNTRL(4) -> FacMin, lower bound on step decrease factor (default=0.2) |
|---|
| 253 | ! RCNTRL(5) -> FacMax, upper bound on step increase factor (default=6) |
|---|
| 254 | ! RCNTRL(6) -> FacRej, step decrease factor after multiple rejections |
|---|
| 255 | ! (default=0.1) |
|---|
| 256 | ! RCNTRL(7) -> FacSafe, by which the new step is slightly smaller |
|---|
| 257 | ! than the predicted value (default=0.9) |
|---|
| 258 | ! RCNTRL(8) -> ThetaMin. If Newton convergence rate smaller |
|---|
| 259 | ! than ThetaMin the Jacobian is not recomputed; |
|---|
| 260 | ! (default=0.001). Increase cntrl(3), to 0.01 say, when |
|---|
| 261 | ! Jacobian evaluations are costly. for small systems it |
|---|
| 262 | ! should be smaller. |
|---|
| 263 | ! RCNTRL(9) -> not used |
|---|
| 264 | ! RCNTRL(10,11) -> FAC1,FAC2 (parameters for step size selection) |
|---|
| 265 | ! RCNTRL(12,13) -> FAC3,FAC4 (parameters for order selection) |
|---|
| 266 | ! RCNTRL(14,15) -> FacSafe1, FacSafe2 |
|---|
| 267 | ! Safety factors for step size prediction |
|---|
| 268 | ! HNEW=H*FacSafe2*(FacSafe1*TOL/ERR)**(1/(J-1)) |
|---|
| 269 | ! RCNTRL(16:19) -> WorkFcn, WorkJac, WorkDec, WorkSol |
|---|
| 270 | ! estimated computational work |
|---|
| 271 | ! |
|---|
| 272 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 273 | ! |
|---|
| 274 | ! OUTPUT PARAMETERS |
|---|
| 275 | ! ----------------- |
|---|
| 276 | ! T T-VALUE WHERE THE SOLUTION IS COMPUTED |
|---|
| 277 | ! (AFTER SUCCESSFUL RETURN T=Tend) |
|---|
| 278 | ! |
|---|
| 279 | ! Y(N) SOLUTION AT T |
|---|
| 280 | ! |
|---|
| 281 | ! H PREDICTED STEP SIZE OF THE LAST ACCEPTED STEP |
|---|
| 282 | ! |
|---|
| 283 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 284 | ! DECLARATIONS |
|---|
| 285 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 286 | IMPLICIT NONE |
|---|
| 287 | INTEGER :: N, IERR, ITOL, Max_no_steps, Ncolumns, Nsequence, Lambda, & |
|---|
| 288 | NRDENS, i, Ncolumns2, NRD, IOUT |
|---|
| 289 | KPP_REAL :: Y(NVAR),AbsTol(*),RelTol(*) |
|---|
| 290 | KPP_REAL :: Tinitial, Tfinal, Roundoff, Hmin, Hmax, & |
|---|
| 291 | FacMin, FacMax, FAC1, FAC2, FAC3, FAC4, FacSafe1, & |
|---|
| 292 | FacSafe2, H, Hstart,WorkFcn,WorkJac, WorkDec, WorkSol,& |
|---|
| 293 | WorkRow, FacRej, FacSafe, ThetaMin, T |
|---|
| 294 | LOGICAL :: AUTNMS |
|---|
| 295 | KPP_REAL :: RCNTRL(20), RSTATUS(20) |
|---|
| 296 | INTEGER :: ICNTRL(20), ISTATUS(20) |
|---|
| 297 | KPP_REAL, PARAMETER :: ZERO = 0.0d0 |
|---|
| 298 | |
|---|
| 299 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 300 | ! SETTING THE PARAMETERS |
|---|
| 301 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 302 | Nfun=0 |
|---|
| 303 | Njac=0 |
|---|
| 304 | Nstp=0 |
|---|
| 305 | Nacc=0 |
|---|
| 306 | Nrej=0 |
|---|
| 307 | Ndec=0 |
|---|
| 308 | Nsol=0 |
|---|
| 309 | |
|---|
| 310 | IERR = 0 |
|---|
| 311 | |
|---|
| 312 | IF (ICNTRL(1) == 0) THEN |
|---|
| 313 | AUTNMS = .FALSE. |
|---|
| 314 | ELSE |
|---|
| 315 | AUTNMS = .TRUE. |
|---|
| 316 | END IF |
|---|
| 317 | |
|---|
| 318 | !~~~> For Scalar tolerances (ICNTRL(1)/=0) the code uses AbsTol(1) and RelTol(1) |
|---|
| 319 | !~~~> For Vector tolerances (ICNTRL(1)==0) the code uses AbsTol(1:NVAR) and RelTol(1:NVAR) |
|---|
| 320 | IF (ICNTRL(2) == 0) THEN |
|---|
| 321 | ITOL = 1 |
|---|
| 322 | ELSE |
|---|
| 323 | ITOL = 0 |
|---|
| 324 | END IF |
|---|
| 325 | |
|---|
| 326 | !~~~> Max_no_steps: the maximum number of time steps admitted |
|---|
| 327 | IF (ICNTRL(4) == 0) THEN |
|---|
| 328 | Max_no_steps = 100000 |
|---|
| 329 | ELSEIF (ICNTRL(4) > 0) THEN |
|---|
| 330 | Max_no_steps=ICNTRL(4) |
|---|
| 331 | ELSE |
|---|
| 332 | PRINT * ,'User-selected ICNTRL(4)=',ICNTRL(4) |
|---|
| 333 | CALL SEULEX_ErrorMsg(-1,Tinitial,ZERO,IERR); |
|---|
| 334 | END IF |
|---|
| 335 | |
|---|
| 336 | !~~~> IOUT = use (or not) the output routine |
|---|
| 337 | IOUT = ICNTRL(10) |
|---|
| 338 | IF ( IOUT<0 .OR. IOUT>2 ) THEN |
|---|
| 339 | PRINT * ,'User-selected ICNTRL(10)=',ICNTRL(10) |
|---|
| 340 | IOUT = 0 |
|---|
| 341 | END IF |
|---|
| 342 | |
|---|
| 343 | !~~~> Ncolumns: maximum number of columns in the extrapolation |
|---|
| 344 | IF (ICNTRL(11)==0) THEN |
|---|
| 345 | Ncolumns=12 |
|---|
| 346 | ELSEIF (ICNTRL(11) > 2) THEN |
|---|
| 347 | Ncolumns=ICNTRL(11) |
|---|
| 348 | ELSE |
|---|
| 349 | PRINT * ,'User-selected ICNTRL(11)=',ICNTRL(11) |
|---|
| 350 | CALL SEULEX_ErrorMsg(-2,Tinitial,ZERO,IERR); |
|---|
| 351 | END IF |
|---|
| 352 | |
|---|
| 353 | !~~~> Nsequence: choice of step size sequence |
|---|
| 354 | IF (ICNTRL(12)==0) THEN |
|---|
| 355 | Nsequence = 2 |
|---|
| 356 | ELSEIF ( (ICNTRL(12)>0).AND.(ICNTRL(12)<5) ) THEN |
|---|
| 357 | Nsequence = ICNTRL(4) |
|---|
| 358 | ELSE |
|---|
| 359 | PRINT * ,'User-selected ICNTRL(12)=',ICNTRL(12) |
|---|
| 360 | CALL SEULEX_ErrorMsg(-3,Tinitial,ZERO,IERR) |
|---|
| 361 | END IF |
|---|
| 362 | |
|---|
| 363 | !~~~> LAMBDA: parameter for dense output |
|---|
| 364 | LAMBDA = ICNTRL(13) |
|---|
| 365 | IF ( LAMBDA < 0 .OR. LAMBDA >= 2 ) THEN |
|---|
| 366 | PRINT * ,'User-selected ICNTRL(13)=',ICNTRL(13) |
|---|
| 367 | CALL SEULEX_ErrorMsg(-4,Tinitial,ZERO,IERR) |
|---|
| 368 | END IF |
|---|
| 369 | |
|---|
| 370 | !~~~>- NRDENS: number of dense output components |
|---|
| 371 | NRDENS=ICNTRL(14) |
|---|
| 372 | IF ( (NRDENS < 0) .OR. (NRDENS > N) ) THEN |
|---|
| 373 | PRINT * ,'User-selected ICNTRL(14)=',ICNTRL(14) |
|---|
| 374 | CALL SEULEX_ErrorMsg(-5,Tinitial,ZERO,IERR) |
|---|
| 375 | END IF |
|---|
| 376 | |
|---|
| 377 | |
|---|
| 378 | !~~~> Unit roundoff (1+Roundoff>1) |
|---|
| 379 | Roundoff = WLAMCH('E') |
|---|
| 380 | |
|---|
| 381 | !~~~> Lower bound on the step size: (positive value) |
|---|
| 382 | IF (RCNTRL(1) == ZERO) THEN |
|---|
| 383 | Hmin = ZERO |
|---|
| 384 | ELSEIF (RCNTRL(1) > ZERO) THEN |
|---|
| 385 | Hmin = RCNTRL(1) |
|---|
| 386 | ELSE |
|---|
| 387 | PRINT * , 'User-selected RCNTRL(1)=', RCNTRL(1) |
|---|
| 388 | CALL SEULEX_ErrorMsg(-7,Tinitial,ZERO,IERR) |
|---|
| 389 | RETURN |
|---|
| 390 | END IF |
|---|
| 391 | |
|---|
| 392 | !~~~> Upper bound on the step size: (positive value) |
|---|
| 393 | IF (RCNTRL(2) == ZERO) THEN |
|---|
| 394 | Hmax = ABS(Tfinal-Tinitial) |
|---|
| 395 | ELSEIF (RCNTRL(2) > ZERO) THEN |
|---|
| 396 | Hmax = MIN(ABS(RCNTRL(2)),ABS(Tfinal-Tinitial)) |
|---|
| 397 | ELSE |
|---|
| 398 | PRINT * , 'User-selected RCNTRL(2)=', RCNTRL(2) |
|---|
| 399 | CALL SEULEX_ErrorMsg(-7,Tinitial,ZERO,IERR) |
|---|
| 400 | RETURN |
|---|
| 401 | END IF |
|---|
| 402 | |
|---|
| 403 | !~~~> Starting step size: (positive value) |
|---|
| 404 | IF (RCNTRL(3) == ZERO) THEN |
|---|
| 405 | Hstart = MAX(Hmin,Roundoff) |
|---|
| 406 | ELSEIF (RCNTRL(3) > ZERO) THEN |
|---|
| 407 | Hstart = MIN(ABS(RCNTRL(3)),ABS(Tfinal-Tinitial)) |
|---|
| 408 | ELSE |
|---|
| 409 | PRINT * , 'User-selected Hstart: RCNTRL(3)=', RCNTRL(3) |
|---|
| 410 | CALL SEULEX_ErrorMsg(-7,Tinitial,ZERO,IERR) |
|---|
| 411 | RETURN |
|---|
| 412 | END IF |
|---|
| 413 | |
|---|
| 414 | |
|---|
| 415 | !~~~> Step size can be changed s.t. FacMin < Hnew/Hexit < FacMax |
|---|
| 416 | IF (RCNTRL(4) == ZERO) THEN |
|---|
| 417 | FacMin = 0.2_dp |
|---|
| 418 | ELSEIF (RCNTRL(4) > ZERO) THEN |
|---|
| 419 | FacMin = RCNTRL(4) |
|---|
| 420 | ELSE |
|---|
| 421 | PRINT * , 'User-selected FacMin: RCNTRL(4)=', RCNTRL(4) |
|---|
| 422 | CALL SEULEX_ErrorMsg(-8,Tinitial,ZERO,IERR) |
|---|
| 423 | RETURN |
|---|
| 424 | END IF |
|---|
| 425 | IF (RCNTRL(5) == ZERO) THEN |
|---|
| 426 | FacMax = 10.0_dp |
|---|
| 427 | ELSEIF (RCNTRL(5) > ZERO) THEN |
|---|
| 428 | FacMax = RCNTRL(5) |
|---|
| 429 | ELSE |
|---|
| 430 | PRINT * , 'User-selected FacMax: RCNTRL(5)=', RCNTRL(5) |
|---|
| 431 | CALL SEULEX_ErrorMsg(-8,Tinitial,ZERO,IERR) |
|---|
| 432 | RETURN |
|---|
| 433 | END IF |
|---|
| 434 | !~~~> FacRej: Factor to decrease step after 2 succesive rejections |
|---|
| 435 | IF (RCNTRL(6) == ZERO) THEN |
|---|
| 436 | FacRej = 0.1_dp |
|---|
| 437 | ELSEIF (RCNTRL(6) > ZERO) THEN |
|---|
| 438 | FacRej = RCNTRL(6) |
|---|
| 439 | ELSE |
|---|
| 440 | PRINT * , 'User-selected FacRej: RCNTRL(6)=', RCNTRL(6) |
|---|
| 441 | CALL SEULEX_ErrorMsg(-8,Tinitial,ZERO,IERR) |
|---|
| 442 | RETURN |
|---|
| 443 | END IF |
|---|
| 444 | !~~~> FacSafe: Safety Factor in the computation of new step size |
|---|
| 445 | IF (RCNTRL(7) == ZERO) THEN |
|---|
| 446 | FacSafe = 0.9_dp |
|---|
| 447 | ELSEIF (RCNTRL(7) > ZERO) THEN |
|---|
| 448 | FacSafe = RCNTRL(7) |
|---|
| 449 | ELSE |
|---|
| 450 | PRINT * , 'User-selected FacSafe: RCNTRL(7)=', RCNTRL(7) |
|---|
| 451 | CALL SEULEX_ErrorMsg(-8,Tinitial,ZERO,IERR) |
|---|
| 452 | RETURN |
|---|
| 453 | END IF |
|---|
| 454 | |
|---|
| 455 | !~~~> ThetaMin: DECIDES WHETHER THE JACOBIAN SHOULD BE RECOMPUTED; |
|---|
| 456 | ! INCREASE WORK(3), TO 0.01 SAY, WHEN JACOBIAN EVALUATIONS |
|---|
| 457 | ! ARE COSTLY. FOR SMALL SYSTEMS WORK(3) SHOULD BE SMALLER. |
|---|
| 458 | IF(RCNTRL(8) == 0.D0)THEN |
|---|
| 459 | ThetaMin = 1.0d-3 |
|---|
| 460 | ELSE |
|---|
| 461 | ThetaMin = RCNTRL(8) |
|---|
| 462 | END IF |
|---|
| 463 | |
|---|
| 464 | !~~~> FAC1,FAC2: PARAMETERS FOR STEP SIZE SELECTION |
|---|
| 465 | ! THE NEW STEP SIZE FOR THE J-TH DIAGONAL ENTRY IS |
|---|
| 466 | ! CHOSEN SUBJECT TO THE RESTRICTION |
|---|
| 467 | ! FACMIN/WORK(5) <= HNEW(J)/HOLD <= 1/FACMIN |
|---|
| 468 | ! WHERE FACMIN=WORK(4)**(1/(J-1)) |
|---|
| 469 | IF(RCNTRL(10) == 0.D0)THEN |
|---|
| 470 | FAC1=0.1D0 |
|---|
| 471 | ELSE |
|---|
| 472 | FAC1=RCNTRL(10) |
|---|
| 473 | END IF |
|---|
| 474 | IF(RCNTRL(11) == 0.D0)THEN |
|---|
| 475 | FAC2=4.0D0 |
|---|
| 476 | ELSE |
|---|
| 477 | FAC2=RCNTRL(11) |
|---|
| 478 | END IF |
|---|
| 479 | !~~~> FAC3, FAC4: PARAMETERS FOR THE ORDER SELECTION |
|---|
| 480 | ! ORDER IS DECREASED IF W(K-1) <= W(K)*WORK(6) |
|---|
| 481 | ! ORDER IS INCREASED IF W(K) <= W(K-1)*WORK(7) |
|---|
| 482 | IF(RCNTRL(12) == 0.D0)THEN |
|---|
| 483 | FAC3=0.7D0 |
|---|
| 484 | ELSE |
|---|
| 485 | FAC3=RCNTRL(12) |
|---|
| 486 | END IF |
|---|
| 487 | IF(RCNTRL(13) == 0.D0)THEN |
|---|
| 488 | FAC4=0.9D0 |
|---|
| 489 | ELSE |
|---|
| 490 | FAC4=RCNTRL(13) |
|---|
| 491 | END IF |
|---|
| 492 | !~~~>- FacSafe1, FacSafe2: safety factors for step size prediction |
|---|
| 493 | ! HNEW=H*WORK(9)*(WORK(8)*TOL/ERR)**(1/(J-1)) |
|---|
| 494 | IF(RCNTRL(14) == 0.D0)THEN |
|---|
| 495 | FacSafe1=0.6D0 |
|---|
| 496 | ELSE |
|---|
| 497 | FacSafe1=RCNTRL(14) |
|---|
| 498 | END IF |
|---|
| 499 | IF(RCNTRL(15) == 0.D0)THEN |
|---|
| 500 | FacSafe2=0.93D0 |
|---|
| 501 | ELSE |
|---|
| 502 | FacSafe2=RCNTRL(15) |
|---|
| 503 | END IF |
|---|
| 504 | |
|---|
| 505 | !~~~> WorkFcn: estimated computational work for a calls to FCN |
|---|
| 506 | IF(RCNTRL(16) == 0.D0)THEN |
|---|
| 507 | WorkFcn=1.D0 |
|---|
| 508 | ELSE |
|---|
| 509 | WorkFcn=RCNTRL(16) |
|---|
| 510 | END IF |
|---|
| 511 | !~~~> WorkJac: estimated computational work for calls to JAC |
|---|
| 512 | IF(RCNTRL(17) == 0.D0)THEN |
|---|
| 513 | WorkJac=5.D0 |
|---|
| 514 | ELSE |
|---|
| 515 | WorkJac=RCNTRL(17) |
|---|
| 516 | END IF |
|---|
| 517 | !~~~> WorkDec: estimated computational work for calls to DEC |
|---|
| 518 | IF(RCNTRL(18) == 0.D0)THEN |
|---|
| 519 | WorkDec=1.D0 |
|---|
| 520 | ELSE |
|---|
| 521 | WorkDec=RCNTRL(18) |
|---|
| 522 | END IF |
|---|
| 523 | !~~~> WorkSol: estimated computational work for calls to SOL |
|---|
| 524 | IF(RCNTRL(19) == 0.D0)THEN |
|---|
| 525 | WorkSol=1.D0 |
|---|
| 526 | ELSE |
|---|
| 527 | WorkSol=RCNTRL(19) |
|---|
| 528 | END IF |
|---|
| 529 | WorkRow=WorkFcn+WorkSol |
|---|
| 530 | |
|---|
| 531 | !~~~> Check if tolerances are reasonable |
|---|
| 532 | IF (ITOL == 0) THEN |
|---|
| 533 | IF (AbsTol(1) <= 0.D0.OR.RelTol(1) <= 10.D0*Roundoff) THEN |
|---|
| 534 | PRINT * , ' Scalar AbsTol = ',AbsTol(1) |
|---|
| 535 | PRINT * , ' Scalar RelTol = ',RelTol(1) |
|---|
| 536 | CALL SEULEX_ErrorMsg(-9,Tinitial,ZERO,IERR) |
|---|
| 537 | END IF |
|---|
| 538 | ELSE |
|---|
| 539 | DO i=1,N |
|---|
| 540 | IF (AbsTol(i) <= 0.D0.OR.RelTol(i) <= 10.D0*Roundoff) THEN |
|---|
| 541 | PRINT * , ' AbsTol(',i,') = ',AbsTol(i) |
|---|
| 542 | PRINT * , ' RelTol(',i,') = ',RelTol(i) |
|---|
| 543 | CALL SEULEX_ErrorMsg(-9,Tinitial,ZERO,IERR) |
|---|
| 544 | END IF |
|---|
| 545 | END DO |
|---|
| 546 | END IF |
|---|
| 547 | |
|---|
| 548 | IF (IERR < 0) RETURN |
|---|
| 549 | |
|---|
| 550 | !~~~>---- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN WORK ----- |
|---|
| 551 | Ncolumns2=(Ncolumns*(Ncolumns+1))/2 |
|---|
| 552 | NRD=MAX(1,NRDENS) |
|---|
| 553 | |
|---|
| 554 | T = Tinitial |
|---|
| 555 | !~~~> CALL TO CORE INTEGRATOR |
|---|
| 556 | CALL SEULEX_Integrator(N,T,Tfinal,Y,Hmax,H,Ncolumns,RelTol,AbsTol,ITOL, & |
|---|
| 557 | IOUT,IERR,Max_no_steps,Roundoff,Nsequence,AUTNMS, & |
|---|
| 558 | FAC1,FAC2,FAC3,FAC4,ThetaMin,FacSafe1,FacSafe2,WorkJac, & |
|---|
| 559 | WorkDec,WorkRow,Ncolumns2,NRD,LAMBDA,Nstp) |
|---|
| 560 | |
|---|
| 561 | ISTATUS(1)=Nfun |
|---|
| 562 | ISTATUS(2)=Njac |
|---|
| 563 | ISTATUS(3)=Nstp |
|---|
| 564 | ISTATUS(4)=Nacc |
|---|
| 565 | ISTATUS(5)=Nrej |
|---|
| 566 | ISTATUS(6)=Ndec |
|---|
| 567 | ISTATUS(7)=Nsol |
|---|
| 568 | |
|---|
| 569 | END SUBROUTINE ATMSEULEX |
|---|
| 570 | |
|---|
| 571 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 572 | SUBROUTINE SEULEX_ErrorMsg(Code,T,H,IERR) |
|---|
| 573 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 574 | ! Handles all error messages |
|---|
| 575 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 576 | |
|---|
| 577 | KPP_REAL, INTENT(IN) :: T, H |
|---|
| 578 | INTEGER, INTENT(IN) :: Code |
|---|
| 579 | INTEGER, INTENT(OUT) :: IERR |
|---|
| 580 | |
|---|
| 581 | IERR = Code |
|---|
| 582 | PRINT * , & |
|---|
| 583 | 'Forced exit from SEULEX due to the following error:' |
|---|
| 584 | |
|---|
| 585 | SELECT CASE (Code) |
|---|
| 586 | CASE (-1) |
|---|
| 587 | PRINT * , '--> Improper value for maximal no of steps' |
|---|
| 588 | CASE (-2) |
|---|
| 589 | PRINT * , '--> Improper value for maximum no of columns in extrapolation' |
|---|
| 590 | CASE (-3) |
|---|
| 591 | PRINT * , '--> Improper value for step size sequence' |
|---|
| 592 | CASE (-4) |
|---|
| 593 | PRINT * , '--> Improper value for Lambda (must be 0/1)' |
|---|
| 594 | CASE (-5) |
|---|
| 595 | PRINT * , '--> Improper number of dense output components' |
|---|
| 596 | CASE (-6) |
|---|
| 597 | PRINT * , '--> Improper parameters for second order equations' |
|---|
| 598 | CASE (-7) |
|---|
| 599 | PRINT * , '--> Hmin/Hmax/Hstart must be positive' |
|---|
| 600 | CASE (-8) |
|---|
| 601 | PRINT * , '--> FacMin/FacMax/FacRej must be positive' |
|---|
| 602 | CASE (-9) |
|---|
| 603 | PRINT * , '--> Improper tolerance values' |
|---|
| 604 | CASE (-10) |
|---|
| 605 | PRINT * , '--> No of steps exceeds maximum bound' |
|---|
| 606 | CASE (-11) |
|---|
| 607 | PRINT * , '--> Step size too small: T + 10*H = T', & |
|---|
| 608 | ' or H < Roundoff' |
|---|
| 609 | CASE (-12) |
|---|
| 610 | PRINT * , '--> Matrix is repeatedly singular' |
|---|
| 611 | CASE DEFAULT |
|---|
| 612 | PRINT *, 'Unknown Error code: ', Code |
|---|
| 613 | END SELECT |
|---|
| 614 | |
|---|
| 615 | PRINT *, "T=", T, "and H=", H |
|---|
| 616 | |
|---|
| 617 | END SUBROUTINE SEULEX_ErrorMsg |
|---|
| 618 | |
|---|
| 619 | |
|---|
| 620 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 621 | SUBROUTINE SEULEX_Integrator(N,T,Tend,Y,Hmax,H,Ncolumns,RelTol,AbsTol,ITOL,& |
|---|
| 622 | IOUT,IERR,Max_no_steps,Roundoff,Nsequence,AUTNMS, & |
|---|
| 623 | FAC1,FAC2,FAC3,FAC4,ThetaMin,FacSafe1,FacSafe2,WorkJac, & |
|---|
| 624 | WorkDec,WorkRow,Ncolumns2,NRD,LAMBDA,Nstp) |
|---|
| 625 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 626 | ! CORE INTEGRATOR FOR SEULEX |
|---|
| 627 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 628 | ! DECLARATIONS |
|---|
| 629 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 630 | USE KPP_ROOT_Parameters |
|---|
| 631 | USE KPP_ROOT_Jacobian |
|---|
| 632 | IMPLICIT KPP_REAL (A-H,O-Z) |
|---|
| 633 | IMPLICIT INTEGER (I-N) |
|---|
| 634 | |
|---|
| 635 | INTEGER :: N, Ncolumns, Ncolumns2, K, KC, KRIGHT, KLR, KK, KRN,& |
|---|
| 636 | KOPT, NRD |
|---|
| 637 | KPP_REAL :: Y(NVAR),DY(NVAR),FX(NVAR),YHH(NVAR) |
|---|
| 638 | KPP_REAL :: DYH(NVAR), DEL(NVAR), WH(NVAR) |
|---|
| 639 | KPP_REAL :: SCAL(NVAR), HH(Ncolumns), W(Ncolumns), A(Ncolumns) |
|---|
| 640 | #ifdef FULL_ALGEBRA |
|---|
| 641 | KPP_REAL :: FJAC(NVAR,NVAR) |
|---|
| 642 | #else |
|---|
| 643 | KPP_REAL :: FJAC(LU_NONZERO) |
|---|
| 644 | #endif |
|---|
| 645 | KPP_REAL Table(Ncolumns,N) |
|---|
| 646 | INTEGER IP(N),NJ(Ncolumns),IPHES(N),ICOMP(NRD) |
|---|
| 647 | KPP_REAL RelTol(*),AbsTol(*) |
|---|
| 648 | KPP_REAL FSAFE(Ncolumns2,NRD),FACUL(Ncolumns),E(N,N),DENS((Ncolumns+2)*NRD) |
|---|
| 649 | LOGICAL REJECT,LAST,ATOV,CALJAC,CALHES,AUTNMS |
|---|
| 650 | |
|---|
| 651 | KPP_REAL TOLDD,HHH,NNRD |
|---|
| 652 | COMMON /COSEU/TOLDD,HHH,NNRD,KRIGHT |
|---|
| 653 | |
|---|
| 654 | !~~~> COMPUTE COEFFICIENTS FOR DENSE OUTPUT |
|---|
| 655 | IF (IOUT == 2) THEN |
|---|
| 656 | NNRD=NRD |
|---|
| 657 | !~~~> COMPUTE THE FACTORIALS -------- |
|---|
| 658 | FACUL(1)=1.D0 |
|---|
| 659 | DO i=1,Ncolumns-1 |
|---|
| 660 | FACUL(i+1)=i*FACUL(i) |
|---|
| 661 | END DO |
|---|
| 662 | END IF |
|---|
| 663 | |
|---|
| 664 | !~~~> DEFINE THE STEP SIZE SEQUENCE |
|---|
| 665 | IF (Nsequence == 1) THEN |
|---|
| 666 | NJ(1)=1 |
|---|
| 667 | NJ(2)=2 |
|---|
| 668 | NJ(3)=3 |
|---|
| 669 | DO I=4,Ncolumns |
|---|
| 670 | NJ(i)=2*NJ(I-2) |
|---|
| 671 | END DO |
|---|
| 672 | END IF |
|---|
| 673 | IF (Nsequence == 2) THEN |
|---|
| 674 | NJ(1)=2 |
|---|
| 675 | NJ(2)=3 |
|---|
| 676 | DO I=3,Ncolumns |
|---|
| 677 | NJ(i)=2*NJ(I-2) |
|---|
| 678 | END DO |
|---|
| 679 | END IF |
|---|
| 680 | DO i=1,Ncolumns |
|---|
| 681 | IF (Nsequence == 3) NJ(i)=I |
|---|
| 682 | IF (Nsequence == 4) NJ(i)=I+1 |
|---|
| 683 | END DO |
|---|
| 684 | A(1)=WorkJac+NJ(1)*WorkRow+WorkDec |
|---|
| 685 | DO I=2,Ncolumns |
|---|
| 686 | A(i)=A(i-1)+(NJ(i)-1)*WorkRow+WorkDec |
|---|
| 687 | END DO |
|---|
| 688 | K=MAX0(3,MIN0(Ncolumns-2,INT(-DLOG10(RelTol(1)+AbsTol(1))*.6D0+1.5D0))) |
|---|
| 689 | |
|---|
| 690 | ! T = Tinitial |
|---|
| 691 | HmaxN = MIN(ABS(Hmax),ABS(Tend-T)) |
|---|
| 692 | IF (ABS(H) <= 10.D0*Roundoff) H=1.0D-6 |
|---|
| 693 | H=MIN(ABS(H),HmaxN) |
|---|
| 694 | Theta=2*ABS(ThetaMin) |
|---|
| 695 | ERR=0.D0 |
|---|
| 696 | W(1)=1.D30 |
|---|
| 697 | DO i=1,N |
|---|
| 698 | IF (ITOL == 0) THEN |
|---|
| 699 | SCAL(i)=AbsTol(1)+RelTol(1)*DABS(Y(i)) |
|---|
| 700 | ELSE |
|---|
| 701 | SCAL(i)=AbsTol(i)+RelTol(i)*DABS(Y(i)) |
|---|
| 702 | END IF |
|---|
| 703 | END DO |
|---|
| 704 | CALJAC=.FALSE. |
|---|
| 705 | REJECT=.FALSE. |
|---|
| 706 | LAST=.FALSE. |
|---|
| 707 | 10 CONTINUE |
|---|
| 708 | IF (REJECT) Theta=2*ABS(ThetaMin) |
|---|
| 709 | ATOV=.FALSE. |
|---|
| 710 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 711 | !~~~> IS Tend REACHED IN THE NEXT STEP? |
|---|
| 712 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 713 | H1=Tend-T |
|---|
| 714 | IF (H1 <= Roundoff) GO TO 110 |
|---|
| 715 | HOPT=H |
|---|
| 716 | H=MIN(H,H1,HmaxN) |
|---|
| 717 | IF (H >= H1-Roundoff) LAST=.TRUE. |
|---|
| 718 | IF (AUTNMS) THEN |
|---|
| 719 | CALL FUN_CHEM(T,Y,DY) |
|---|
| 720 | END IF |
|---|
| 721 | IF (Theta > ThetaMin.AND..NOT.CALJAC) THEN |
|---|
| 722 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 723 | ! COMPUTATION OF THE JACOBIAN |
|---|
| 724 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 725 | CALL JAC_CHEM(T,Y,FJAC) |
|---|
| 726 | CALJAC=.TRUE. |
|---|
| 727 | CALHES=.FALSE. |
|---|
| 728 | END IF |
|---|
| 729 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 730 | !~~~> THE FIRST AND LAST STEP |
|---|
| 731 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 732 | IF (Nstp == 0.OR.LAST) THEN |
|---|
| 733 | IPT=0 |
|---|
| 734 | Nstp=Nstp+1 |
|---|
| 735 | DO J=1,K |
|---|
| 736 | KC=J |
|---|
| 737 | CALL SEUL(J,N,T,Y,DY,FX,FJAC,LFJAC,E,LE,IP,H,Ncolumns, & |
|---|
| 738 | HmaxN,Table,SCAL,NJ,HH,W,A,YHH,DYH,DEL,WH,ERR,FacSafe1,FAC, & |
|---|
| 739 | FAC1,FAC2,FacSafe2,Theta,Nfun,Ndec,Nsol, & |
|---|
| 740 | ERROLD,IPHES,ICOMP,AUTNMS,REJECT, & |
|---|
| 741 | ATOV,FSAFE,Ncolumns2,NRD,IOUT,IPT,CALHES) |
|---|
| 742 | IF (ATOV) GOTO 10 |
|---|
| 743 | IF (J > 1 .AND. ERR <= 1.d0) GOTO 60 |
|---|
| 744 | END DO |
|---|
| 745 | GO TO 55 |
|---|
| 746 | END IF |
|---|
| 747 | !~~~> BASIC INTEGRATION STEP |
|---|
| 748 | 30 CONTINUE |
|---|
| 749 | IPT=0 |
|---|
| 750 | Nstp=Nstp+1 |
|---|
| 751 | IF (Nstp >= Max_no_steps) GOTO 120 |
|---|
| 752 | KC=K-1 |
|---|
| 753 | DO J=1,KC |
|---|
| 754 | CALL SEUL(J,N,T,Y,DY,FX,FJAC,LFJAC,E,LE,IP,H,Ncolumns,& |
|---|
| 755 | HmaxN,Table,SCAL,NJ,HH,W,A,YHH,DYH,DEL,WH,ERR,FacSafe1,& |
|---|
| 756 | FAC,FAC1,FAC2,FacSafe2,Theta,Nfun,Ndec,Nsol,& |
|---|
| 757 | ERROLD,IPHES,ICOMP,AUTNMS,REJECT,& |
|---|
| 758 | ATOV,FSAFE,Ncolumns2,NRD,IOUT,IPT,CALHES) |
|---|
| 759 | IF (ATOV) GOTO 10 |
|---|
| 760 | END DO |
|---|
| 761 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 762 | !~~~> CONVERGENCE MONITOR |
|---|
| 763 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 764 | IF (K == 2.OR.REJECT) GO TO 50 |
|---|
| 765 | IF (ERR <= 1.D0) GO TO 60 |
|---|
| 766 | IF (ERR > DBLE(NJ(K+1)*NJ(K))*4.D0) GO TO 100 |
|---|
| 767 | 50 CALL SEUL(K,N,T,Y,DY,FX,FJAC,LFJAC,E,LE,IP,H,Ncolumns,& |
|---|
| 768 | HmaxN,Table,SCAL,NJ,HH,W,A,YHH,DYH,DEL,WH,ERR,FacSafe1,& |
|---|
| 769 | FAC,FAC1,FAC2,FacSafe2,Theta,Nfun,Ndec,Nsol,& |
|---|
| 770 | ERROLD,IPHES,ICOMP,AUTNMS,REJECT,& |
|---|
| 771 | ATOV,FSAFE,Ncolumns2,NRD,IOUT,IPT,CALHES) |
|---|
| 772 | IF (ATOV) GOTO 10 |
|---|
| 773 | KC=K |
|---|
| 774 | IF (ERR <= 1.D0) GO TO 60 |
|---|
| 775 | !~~~> HOPE FOR CONVERGENCE IN LINE K+1 |
|---|
| 776 | 55 IF (ERR > DBLE(NJ(K+1))*2.D0) GO TO 100 |
|---|
| 777 | KC=K+1 |
|---|
| 778 | CALL SEUL(KC,N,T,Y,DY,FX,FJAC,LFJAC,E,LE,IP,H,Ncolumns,& |
|---|
| 779 | HmaxN,Table,SCAL,NJ,HH,W,A,YHH,DYH,DEL,WH,ERR,FacSafe1,& |
|---|
| 780 | FAC,FAC1,FAC2,FacSafe2,Theta,Nfun,Ndec,Nsol,& |
|---|
| 781 | ERROLD,IPHES,ICOMP,AUTNMS,REJECT,& |
|---|
| 782 | ATOV,FSAFE,Ncolumns2,NRD,IOUT,IPT,CALHES) |
|---|
| 783 | IF (ATOV) GOTO 10 |
|---|
| 784 | IF (ERR > 1.D0) GO TO 100 |
|---|
| 785 | !Adi IF ((ERR > 1.D0).and.(H.gt.Hmin)) GO TO 100 |
|---|
| 786 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 787 | !~~~> STEP IS ACCEPTED |
|---|
| 788 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 789 | 60 TOLD=T |
|---|
| 790 | T=T+H |
|---|
| 791 | IF (IOUT == 2) THEN |
|---|
| 792 | KRIGHT=KC |
|---|
| 793 | DO i=1,NRD |
|---|
| 794 | DENS(i)=Y(ICOMP(i)) |
|---|
| 795 | END DO |
|---|
| 796 | END IF |
|---|
| 797 | DO i=1,N |
|---|
| 798 | T1I=Table(1,I) |
|---|
| 799 | IF (ITOL == 0) THEN |
|---|
| 800 | SCAL(i)=AbsTol(1)+RelTol(1)*DABS(T1I) |
|---|
| 801 | ELSE |
|---|
| 802 | SCAL(i)=AbsTol(i)+RelTol(i)*DABS(T1I) |
|---|
| 803 | END IF |
|---|
| 804 | Y(i)=T1I |
|---|
| 805 | END DO |
|---|
| 806 | Nacc=Nacc+1 |
|---|
| 807 | CALJAC=.FALSE. |
|---|
| 808 | IF (IOUT == 2) THEN |
|---|
| 809 | TOLDD=TOLD |
|---|
| 810 | HHH=H |
|---|
| 811 | DO i=1,NRD |
|---|
| 812 | DENS(NRD+I)=Y(ICOMP(i)) |
|---|
| 813 | END DO |
|---|
| 814 | DO KLR=1,KRIGHT-1 |
|---|
| 815 | !~~~> COMPUTE DIFFERENCES |
|---|
| 816 | IF (KLR >= 2) THEN |
|---|
| 817 | DO KK=KLR,KC |
|---|
| 818 | LBEG=((KK+1)*KK)/2 |
|---|
| 819 | LEND=LBEG-KK+2 |
|---|
| 820 | DO L=LBEG,LEND,-1 |
|---|
| 821 | DO i=1,NRD |
|---|
| 822 | FSAFE(L,I)=FSAFE(L,I)-FSAFE(L-1,I) |
|---|
| 823 | END DO |
|---|
| 824 | END DO |
|---|
| 825 | END DO |
|---|
| 826 | END IF |
|---|
| 827 | !~~~> COMPUTE DERIVATIVES AT RIGHT END ---- |
|---|
| 828 | DO KK=KLR+LAMBDA,KC |
|---|
| 829 | FACNJ=NJ(KK) |
|---|
| 830 | FACNJ=FACNJ**KLR/FACUL(KLR+1) |
|---|
| 831 | IPT=((KK+1)*KK)/2 |
|---|
| 832 | DO I=1,NRD |
|---|
| 833 | KRN=(KK-LAMBDA+1)*NRD |
|---|
| 834 | DENS(KRN+I)=FSAFE(IPT,I)*FACNJ |
|---|
| 835 | END DO |
|---|
| 836 | END DO |
|---|
| 837 | DO J=KLR+LAMBDA+1,KC |
|---|
| 838 | DBLENJ=NJ(J) |
|---|
| 839 | DO L=J,KLR+LAMBDA+1,-1 |
|---|
| 840 | FACTOR=DBLENJ/NJ(L-1)-1.D0 |
|---|
| 841 | DO i=1,NRD |
|---|
| 842 | KRN=(L-LAMBDA+1)*NRD+I |
|---|
| 843 | DENS(KRN-NRD)=DENS(KRN)+(DENS(KRN)-DENS(KRN-NRD))/FACTOR |
|---|
| 844 | END DO |
|---|
| 845 | END DO |
|---|
| 846 | END DO |
|---|
| 847 | END DO |
|---|
| 848 | !~~~> COMPUTE THE COEFFICIENTS OF THE INTERPOLATION POLYNOMIAL |
|---|
| 849 | DO IN=1,NRD |
|---|
| 850 | DO J=1,KRIGHT |
|---|
| 851 | II=NRD*J+IN |
|---|
| 852 | DENS(II)=DENS(II)-DENS(II-NRD) |
|---|
| 853 | END DO |
|---|
| 854 | END DO |
|---|
| 855 | END IF |
|---|
| 856 | !~~~> COMPUTE OPTIMAL ORDER |
|---|
| 857 | IF (KC == 2) THEN |
|---|
| 858 | KOPT=3 |
|---|
| 859 | IF (REJECT) KOPT=2 |
|---|
| 860 | GO TO 80 |
|---|
| 861 | END IF |
|---|
| 862 | IF (KC <= K) THEN |
|---|
| 863 | KOPT=KC |
|---|
| 864 | IF (W(KC-1) < W(KC)*FAC3) KOPT=KC-1 |
|---|
| 865 | IF (W(KC) < W(KC-1)*FAC4) KOPT=MIN0(KC+1,Ncolumns-1) |
|---|
| 866 | ELSE |
|---|
| 867 | KOPT=KC-1 |
|---|
| 868 | IF (KC > 3.AND.W(KC-2) < W(KC-1)*FAC3) KOPT=KC-2 |
|---|
| 869 | IF (W(KC) < W(KOPT)*FAC4) KOPT=MIN0(KC,Ncolumns-1) |
|---|
| 870 | END IF |
|---|
| 871 | !~~~> AFTER A REJECTED STEP |
|---|
| 872 | 80 IF (REJECT) THEN |
|---|
| 873 | K=MIN0(KOPT,KC) |
|---|
| 874 | H=MIN(H,HH(K)) |
|---|
| 875 | REJECT=.FALSE. |
|---|
| 876 | GO TO 10 |
|---|
| 877 | END IF |
|---|
| 878 | !~~~> COMPUTE STEP SIZE FOR NEXT STEP |
|---|
| 879 | IF (KOPT <= KC) THEN |
|---|
| 880 | H=HH(KOPT) |
|---|
| 881 | ELSE |
|---|
| 882 | IF (KC < K.AND.W(KC) < W(KC-1)*FAC4) THEN |
|---|
| 883 | H=HH(KC)*A(KOPT+1)/A(KC) |
|---|
| 884 | ELSE |
|---|
| 885 | H=HH(KC)*A(KOPT)/A(KC) |
|---|
| 886 | END IF |
|---|
| 887 | END IF |
|---|
| 888 | K=KOPT |
|---|
| 889 | !Adi H = MAX(H, Hmin) |
|---|
| 890 | GO TO 10 |
|---|
| 891 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 892 | !~~~> STEP IS REJECTED |
|---|
| 893 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 894 | 100 K=MIN0(K,KC) |
|---|
| 895 | IF (K > 2.AND.W(K-1) < W(K)*FAC3) K=K-1 |
|---|
| 896 | Nrej=Nrej+1 |
|---|
| 897 | H=HH(K) |
|---|
| 898 | LAST=.FALSE. |
|---|
| 899 | REJECT=.TRUE. |
|---|
| 900 | IF (CALJAC) GOTO 30 |
|---|
| 901 | GO TO 10 |
|---|
| 902 | !~~~> SOLUTION EXIT |
|---|
| 903 | 110 CONTINUE |
|---|
| 904 | H=HOPT |
|---|
| 905 | IERR=1 |
|---|
| 906 | RETURN |
|---|
| 907 | !~~~> FAIL EXIT |
|---|
| 908 | 120 WRITE (6,979) T,H |
|---|
| 909 | 979 FORMAT(' EXIT OF SEULEX AT T=',D14.7,' H=',D14.7) |
|---|
| 910 | IERR=-1 |
|---|
| 911 | RETURN |
|---|
| 912 | |
|---|
| 913 | |
|---|
| 914 | END SUBROUTINE SEULEX_Integrator |
|---|
| 915 | |
|---|
| 916 | |
|---|
| 917 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 918 | SUBROUTINE SEUL(JJ,N,T,Y,DY,FX,FJAC,LFJAC,E,LE,IP,& |
|---|
| 919 | H,Ncolumns,HmaxN,Table,SCAL,NJ,HH,W,A,YH,DYH,DEL,WH,ERR,FacSafe1, & |
|---|
| 920 | FAC,FAC1,FAC2,FacSafe2,Theta,Nfun,Ndec,Nsol,& |
|---|
| 921 | ERROLD,IPHES,ICOMP, & |
|---|
| 922 | AUTNMS,REJECT,ATOV,FSAFE,Ncolumns2,NRD,IOUT, & |
|---|
| 923 | IPT,CALHES) |
|---|
| 924 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 925 | !~~~> THIS SUBROUTINE COMPUTES THE J-TH LINE OF THE |
|---|
| 926 | !~~~> EXTRAPOLATION TABLE AND PROVIDES AN ESTIMATE |
|---|
| 927 | !~~~> OF THE OPTIMAL STEP SIZE |
|---|
| 928 | USE KPP_ROOT_Parameters |
|---|
| 929 | USE KPP_ROOT_Jacobian |
|---|
| 930 | IMPLICIT KPP_REAL (A-H,O-Z) |
|---|
| 931 | IMPLICIT INTEGER (I-N) |
|---|
| 932 | INTEGER :: Ncolumns, Ncolumns2, N, NRD |
|---|
| 933 | KPP_REAL :: Y(NVAR),YH(NVAR),DY(NVAR),FX(NVAR),DYH(NVAR) |
|---|
| 934 | KPP_REAL :: DEL(NVAR),WH(NVAR),SCAL(NVAR),HH(Ncolumns),W(Ncolumns),A(Ncolumns) |
|---|
| 935 | #ifdef FULL_ALGEBRA |
|---|
| 936 | KPP_REAL :: FJAC(NVAR,NVAR), E(NVAR,NVAR) |
|---|
| 937 | #else |
|---|
| 938 | KPP_REAL :: FJAC(LU_NONZERO), E(LU_NONZERO) |
|---|
| 939 | #endif |
|---|
| 940 | KPP_REAL :: Table(Ncolumns,NVAR) |
|---|
| 941 | KPP_REAL :: FSAFE(Ncolumns2,NRD) |
|---|
| 942 | INTEGER :: IP(N),NJ(Ncolumns),IPHES(N),ICOMP(NRD) |
|---|
| 943 | LOGICAL ATOV,REJECT,AUTNMS,CALHES |
|---|
| 944 | |
|---|
| 945 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 946 | ! COMPUTE THE MATRIX E AND ITS DECOMPOSITION |
|---|
| 947 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 948 | HJ=H/NJ(JJ) |
|---|
| 949 | HJI=1.D0/HJ |
|---|
| 950 | #ifdef FULL_ALGEBRA |
|---|
| 951 | DO j=1,N |
|---|
| 952 | DO i=1,N |
|---|
| 953 | E(i,j)=-FJAC(i,j) |
|---|
| 954 | END DO |
|---|
| 955 | E(j,j)=E(j,j)+HJI |
|---|
| 956 | END DO |
|---|
| 957 | CALL DGETRF(N,N,E,N,IP,ISING) |
|---|
| 958 | #else |
|---|
| 959 | DO i=1,LU_NONZERO |
|---|
| 960 | E(i)=-FJAC(i) |
|---|
| 961 | END DO |
|---|
| 962 | DO j=1,N |
|---|
| 963 | E(LU_DIAG(j))=E(LU_DIAG(j))+HJI |
|---|
| 964 | END DO |
|---|
| 965 | CALL KppDecomp (E,ISING) |
|---|
| 966 | #endif |
|---|
| 967 | Ndec=Ndec+1 |
|---|
| 968 | IF (ISING.NE.0) GOTO 79 |
|---|
| 969 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 970 | !~~~> STARTING PROCEDURE |
|---|
| 971 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 972 | IF (.NOT.AUTNMS) THEN |
|---|
| 973 | CALL FUN_CHEM(T+HJ,Y,DY) |
|---|
| 974 | END IF |
|---|
| 975 | DO i=1,N |
|---|
| 976 | YH(i)=Y(i) |
|---|
| 977 | DEL(i)=DY(i) |
|---|
| 978 | END DO |
|---|
| 979 | #ifdef FULL_ALGEBRA |
|---|
| 980 | CALL DGETRS ('N',N,1,E,N,IP,DEL,N,ISING) |
|---|
| 981 | #else |
|---|
| 982 | CALL KppSolve (E,DEL) |
|---|
| 983 | #endif |
|---|
| 984 | Nsol=Nsol+1 |
|---|
| 985 | M=NJ(JJ) |
|---|
| 986 | IF (IOUT == 2.AND.M == JJ) THEN |
|---|
| 987 | IPT=IPT+1 |
|---|
| 988 | DO i=1,NRD |
|---|
| 989 | FSAFE(IPT,I)=DEL(ICOMP(i)) |
|---|
| 990 | END DO |
|---|
| 991 | END IF |
|---|
| 992 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 993 | !~~~> SEMI-IMPLICIT EULER METHOD |
|---|
| 994 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 995 | IF (M > 1) THEN |
|---|
| 996 | DO MM=1,M-1 |
|---|
| 997 | DO i=1,N |
|---|
| 998 | YH(i)=YH(i)+DEL(i) |
|---|
| 999 | END DO |
|---|
| 1000 | IF (AUTNMS) THEN |
|---|
| 1001 | CALL FUN_CHEM(T+HJ*MM,YH,DYH) |
|---|
| 1002 | ELSE |
|---|
| 1003 | CALL FUN_CHEM(T+HJ*(MM+1),YH,DYH) |
|---|
| 1004 | END IF |
|---|
| 1005 | |
|---|
| 1006 | IF (MM == 1.AND.JJ <= 2) THEN |
|---|
| 1007 | !~~~> STABILITY CHECK |
|---|
| 1008 | DEL1=0.D0 |
|---|
| 1009 | DO i=1,N |
|---|
| 1010 | DEL1=DEL1+(DEL(i)/SCAL(i))**2 |
|---|
| 1011 | END DO |
|---|
| 1012 | DEL1=SQRT(DEL1) |
|---|
| 1013 | IF (.NOT.AUTNMS) THEN |
|---|
| 1014 | CALL FUN_CHEM(T+HJ,YH,WH) |
|---|
| 1015 | |
|---|
| 1016 | DO i=1,N |
|---|
| 1017 | DEL(i)=WH(i)-DEL(i)*HJI |
|---|
| 1018 | END DO |
|---|
| 1019 | ELSE |
|---|
| 1020 | DO i=1,N |
|---|
| 1021 | DEL(i)=DYH(i)-DEL(i)*HJI |
|---|
| 1022 | END DO |
|---|
| 1023 | END IF |
|---|
| 1024 | #ifdef FULL_ALGEBRA |
|---|
| 1025 | CALL DGETRS ('N',N,1,E,N,IP,DEL,N,ISING) |
|---|
| 1026 | #else |
|---|
| 1027 | CALL KppSolve (E,DEL) |
|---|
| 1028 | #endif |
|---|
| 1029 | Nsol=Nsol+1 |
|---|
| 1030 | DEL2=0.D0 |
|---|
| 1031 | DO i=1,N |
|---|
| 1032 | DEL2=DEL2+(DEL(i)/SCAL(i))**2 |
|---|
| 1033 | END DO |
|---|
| 1034 | DEL2=SQRT(DEL2) |
|---|
| 1035 | Theta=DEL2/MAX(1.D0,DEL1) |
|---|
| 1036 | IF (Theta > 1.D0) GOTO 79 |
|---|
| 1037 | END IF |
|---|
| 1038 | #ifdef FULL_ALGEBRA |
|---|
| 1039 | CALL DGETRS ('N',N,1,E,N,IP,DYH,N,ISING) |
|---|
| 1040 | #else |
|---|
| 1041 | CALL KppSolve (E,DYH) |
|---|
| 1042 | #endif |
|---|
| 1043 | Nsol=Nsol+1 |
|---|
| 1044 | DO i=1,N |
|---|
| 1045 | DEL(i)=DYH(i) |
|---|
| 1046 | END DO |
|---|
| 1047 | IF (IOUT == 2.AND.MM >= M-JJ) THEN |
|---|
| 1048 | IPT=IPT+1 |
|---|
| 1049 | DO i=1,NRD |
|---|
| 1050 | FSAFE(IPT,i)=DEL(ICOMP(i)) |
|---|
| 1051 | END DO |
|---|
| 1052 | END IF |
|---|
| 1053 | END DO |
|---|
| 1054 | END IF |
|---|
| 1055 | DO i=1,N |
|---|
| 1056 | Table(JJ,I)=YH(i)+DEL(i) |
|---|
| 1057 | END DO |
|---|
| 1058 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 1059 | !~~~> POLYNOMIAL EXTRAPOLATION |
|---|
| 1060 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 1061 | IF (JJ == 1) RETURN |
|---|
| 1062 | DO L=JJ,2,-1 |
|---|
| 1063 | FAC=(DBLE(NJ(JJ))/DBLE(NJ(L-1)))-1.D0 |
|---|
| 1064 | DO i=1,N |
|---|
| 1065 | Table(L-1,I)=Table(L,I)+(Table(L,I)-Table(L-1,I))/FAC |
|---|
| 1066 | END DO |
|---|
| 1067 | END DO |
|---|
| 1068 | ERR=0.D0 |
|---|
| 1069 | DO i=1,N |
|---|
| 1070 | ERR=ERR+MIN(ABS((Table(1,I)-Table(2,I)))/SCAL(i),1.D15)**2 |
|---|
| 1071 | END DO |
|---|
| 1072 | IF (ERR >= 1.D30) GOTO 79 |
|---|
| 1073 | ERR=SQRT(ERR/DBLE(N)) |
|---|
| 1074 | IF (JJ > 2.AND.ERR >= ERROLD) GOTO 79 |
|---|
| 1075 | ERROLD=MAX(4*ERR,1.D0) |
|---|
| 1076 | !~~~> COMPUTE OPTIMAL STEP SIZES |
|---|
| 1077 | EXPO=1.D0/JJ |
|---|
| 1078 | FACMIN=FAC1**EXPO |
|---|
| 1079 | FAC=MIN(FAC2/FACMIN,MAX(FACMIN,(ERR/FacSafe1)**EXPO/FacSafe2)) |
|---|
| 1080 | FAC=1.D0/FAC |
|---|
| 1081 | HH(JJ)=MIN(H*FAC,HmaxN) |
|---|
| 1082 | W(JJ)=A(JJ)/HH(JJ) |
|---|
| 1083 | RETURN |
|---|
| 1084 | 79 ATOV=.TRUE. |
|---|
| 1085 | H=H*0.5D0 |
|---|
| 1086 | REJECT=.TRUE. |
|---|
| 1087 | RETURN |
|---|
| 1088 | END SUBROUTINE SEUL |
|---|
| 1089 | |
|---|
| 1090 | |
|---|
| 1091 | |
|---|
| 1092 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 1093 | SUBROUTINE FUN_CHEM( T, V, FCT ) |
|---|
| 1094 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 1095 | |
|---|
| 1096 | USE KPP_ROOT_Parameters |
|---|
| 1097 | USE KPP_ROOT_Global |
|---|
| 1098 | USE KPP_ROOT_Function, ONLY: Fun |
|---|
| 1099 | USE KPP_ROOT_Rates, ONLY: Update_SUN, Update_RCONST, Update_PHOTO |
|---|
| 1100 | |
|---|
| 1101 | IMPLICIT NONE |
|---|
| 1102 | |
|---|
| 1103 | KPP_REAL :: V(NVAR), FCT(NVAR) |
|---|
| 1104 | KPP_REAL :: T, TOLD |
|---|
| 1105 | |
|---|
| 1106 | !TOLD = TIME |
|---|
| 1107 | !TIME = T |
|---|
| 1108 | !CALL Update_SUN() |
|---|
| 1109 | !CALL Update_RCONST() |
|---|
| 1110 | !CALL Update_PHOTO() |
|---|
| 1111 | !TIME = TOLD |
|---|
| 1112 | CALL Fun(V, FIX, RCONST, FCT) |
|---|
| 1113 | Nfun=Nfun+1 |
|---|
| 1114 | |
|---|
| 1115 | END SUBROUTINE FUN_CHEM |
|---|
| 1116 | |
|---|
| 1117 | |
|---|
| 1118 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 1119 | SUBROUTINE JAC_CHEM ( T, V, Jcb ) |
|---|
| 1120 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 1121 | |
|---|
| 1122 | USE KPP_ROOT_Parameters |
|---|
| 1123 | USE KPP_ROOT_Global |
|---|
| 1124 | USE KPP_ROOT_Jacobian, ONLY: Jac_SP |
|---|
| 1125 | USE KPP_ROOT_Rates, ONLY: Update_SUN, Update_RCONST, Update_PHOTO |
|---|
| 1126 | |
|---|
| 1127 | IMPLICIT NONE |
|---|
| 1128 | |
|---|
| 1129 | KPP_REAL :: V(NVAR), T, TOLD |
|---|
| 1130 | #ifdef FULL_ALGEBRA |
|---|
| 1131 | KPP_REAL :: JV(LU_NONZERO), Jcb(NVAR,NVAR) |
|---|
| 1132 | INTEGER :: i,j |
|---|
| 1133 | #else |
|---|
| 1134 | KPP_REAL :: Jcb(LU_NONZERO) |
|---|
| 1135 | #endif |
|---|
| 1136 | |
|---|
| 1137 | !TOLD = TIME |
|---|
| 1138 | !TIME = T |
|---|
| 1139 | !CALL Update_SUN() |
|---|
| 1140 | !CALL Update_RCONST() |
|---|
| 1141 | !CALL Update_PHOTO() |
|---|
| 1142 | !TIME = TOLD |
|---|
| 1143 | |
|---|
| 1144 | #ifdef FULL_ALGEBRA |
|---|
| 1145 | CALL Jac_SP(V, FIX, RCONST, JV) |
|---|
| 1146 | DO j=1,NVAR |
|---|
| 1147 | DO i=1,NVAR |
|---|
| 1148 | Jcb(i,j) = 0.0D0 |
|---|
| 1149 | END DO |
|---|
| 1150 | END DO |
|---|
| 1151 | DO i=1,LU_NONZERO |
|---|
| 1152 | Jcb(LU_IROW(i),LU_ICOL(i)) = JV(i) |
|---|
| 1153 | END DO |
|---|
| 1154 | #else |
|---|
| 1155 | CALL Jac_SP(V, FIX, RCONST, Jcb) |
|---|
| 1156 | #endif |
|---|
| 1157 | Njac=Njac+1 |
|---|
| 1158 | |
|---|
| 1159 | END SUBROUTINE JAC_CHEM |
|---|
| 1160 | |
|---|
| 1161 | !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
|---|
| 1162 | |
|---|
| 1163 | END MODULE KPP_ROOT_Integrator |
|---|
