[2696] | 1 | SUBROUTINE INTEGRATE( TIN, TOUT ) |
---|
| 2 | |
---|
| 3 | INCLUDE 'KPP_ROOT_params.h' |
---|
| 4 | INCLUDE 'KPP_ROOT_global.h' |
---|
| 5 | |
---|
| 6 | C TIN - Start Time |
---|
| 7 | KPP_REAL TIN |
---|
| 8 | C TOUT - End Time |
---|
| 9 | KPP_REAL TOUT |
---|
| 10 | |
---|
| 11 | INTEGER INFO(5) |
---|
| 12 | |
---|
| 13 | EXTERNAL FUNC_CHEM, JAC_CHEM |
---|
| 14 | |
---|
| 15 | INFO(1) = Autonomous |
---|
| 16 | |
---|
| 17 | CALL ROS1(NVAR,TIN,TOUT,STEPMIN,VAR, |
---|
| 18 | + Info,FUNC_CHEM,JAC_CHEM) |
---|
| 19 | |
---|
| 20 | |
---|
| 21 | RETURN |
---|
| 22 | END |
---|
| 23 | |
---|
| 24 | |
---|
| 25 | |
---|
| 26 | |
---|
| 27 | SUBROUTINE ROS1(N,T,Tnext,Hstart, |
---|
| 28 | + y,Info,FUNC_CHEM,JAC_CHEM) |
---|
| 29 | |
---|
| 30 | INCLUDE 'KPP_ROOT_params.h' |
---|
| 31 | INCLUDE 'KPP_ROOT_sparse.h' |
---|
| 32 | C |
---|
| 33 | C Linearly Implicit Euler |
---|
| 34 | C A method of theoretical interest but of no practical value |
---|
| 35 | C |
---|
| 36 | C INPUT ARGUMENTS: |
---|
| 37 | C y = Vector of (NVAR) concentrations, contains the |
---|
| 38 | C initial values on input |
---|
| 39 | C [T, Tnext] = the integration interval |
---|
| 40 | C Hmin, Hmax = lower and upper bounds for the selected step-size. |
---|
| 41 | C Note that for Step = Hmin the current computed |
---|
| 42 | C solution is unconditionally accepted by the error |
---|
| 43 | C control mechanism. |
---|
| 44 | C AbsTol, RelTol = (NVAR) dimensional vectors of |
---|
| 45 | C componentwise absolute and relative tolerances. |
---|
| 46 | C FUNC_CHEM = name of routine of derivatives. KPP syntax. |
---|
| 47 | C See the header below. |
---|
| 48 | C JAC_CHEM = name of routine that computes the Jacobian, in |
---|
| 49 | C sparse format. KPP syntax. See the header below. |
---|
| 50 | C Info(1) = 1 for Autonomous system |
---|
| 51 | C = 0 for nonAutonomous system |
---|
| 52 | C |
---|
| 53 | C OUTPUT ARGUMENTS: |
---|
| 54 | C y = the values of concentrations at Tend. |
---|
| 55 | C T = equals TENDon output. |
---|
| 56 | C Info(2) = # of FUNC_CHEM CALLs. |
---|
| 57 | C Info(3) = # of JAC_CHEM CALLs. |
---|
| 58 | C Info(4) = # of accepted steps. |
---|
| 59 | C Info(5) = # of rejected steps. |
---|
| 60 | C Hstart = The last accepted stepsize |
---|
| 61 | C |
---|
| 62 | C Adrian Sandu, December 2001 |
---|
| 63 | C |
---|
| 64 | KPP_REAL Fv(NVAR) |
---|
| 65 | KPP_REAL JAC(LU_NONZERO) |
---|
| 66 | KPP_REAL H, Hstart |
---|
| 67 | KPP_REAL y(NVAR) |
---|
| 68 | KPP_REAL T, Tnext, Tplus |
---|
| 69 | KPP_REAL elo,ghinv,uround |
---|
| 70 | |
---|
| 71 | INTEGER n,nfcn,njac,Naccept,Nreject,i,j |
---|
| 72 | INTEGER Info(5) |
---|
| 73 | LOGICAL IsReject, Autonomous |
---|
| 74 | EXTERNAL FUNC_CHEM, JAC_CHEM |
---|
| 75 | |
---|
| 76 | |
---|
| 77 | H = Hstart |
---|
| 78 | Tplus = T |
---|
| 79 | Nfcn = 0 |
---|
| 80 | Njac = 0 |
---|
| 81 | |
---|
| 82 | C === Starting the time loop === |
---|
| 83 | 10 CONTINUE |
---|
| 84 | |
---|
| 85 | Tplus = T + H |
---|
| 86 | IF ( Tplus .gt. Tnext ) THEN |
---|
| 87 | H = Tnext - T |
---|
| 88 | Tplus = Tnext |
---|
| 89 | END IF |
---|
| 90 | |
---|
| 91 | C Initial Function and Jacobian values |
---|
| 92 | CALL FUNC_CHEM(NVAR, T, y, Fv) |
---|
| 93 | Nfcn = Nfcn+1 |
---|
| 94 | CALL JAC_CHEM(NVAR, T, y, JAC) |
---|
| 95 | Njac = Njac+1 |
---|
| 96 | |
---|
| 97 | C Form the Prediction matrix and compute its LU factorization |
---|
| 98 | DO 40 j=1,NVAR |
---|
| 99 | JAC(LU_DIAG(j)) = JAC(LU_DIAG(j)) - 1.0d0/H |
---|
| 100 | 40 CONTINUE |
---|
| 101 | CALL KppDecomp (JAC, ier) |
---|
| 102 | C |
---|
| 103 | IF (ier.ne.0) THEN |
---|
| 104 | PRINT *,'ROS1: Singular factorization at T=',T,'; H=',H |
---|
| 105 | STOP |
---|
| 106 | END IF |
---|
| 107 | |
---|
| 108 | C ------------ STAGE 1------------------------- |
---|
| 109 | CALL KppSolve (JAC, Fv) |
---|
| 110 | |
---|
| 111 | C ---- The Solution --- |
---|
| 112 | DO 160 j = 1,NVAR |
---|
| 113 | y(j) = y(j) - Fv(j) |
---|
| 114 | 160 CONTINUE |
---|
| 115 | T = T + H |
---|
| 116 | |
---|
| 117 | C ======= End of the time loop =============================== |
---|
| 118 | IF ( T .lt. Tnext ) GO TO 10 |
---|
| 119 | |
---|
| 120 | C ======= Output Information ================================= |
---|
| 121 | Info(2) = Nfcn |
---|
| 122 | Info(3) = Njac |
---|
| 123 | Info(4) = Njac |
---|
| 124 | Info(5) = 0 |
---|
| 125 | Hstart = H |
---|
| 126 | |
---|
| 127 | RETURN |
---|
| 128 | END |
---|
| 129 | |
---|
| 130 | |
---|
| 131 | SUBROUTINE FUNC_CHEM(N, T, Y, P) |
---|
| 132 | INCLUDE 'KPP_ROOT_params.h' |
---|
| 133 | INCLUDE 'KPP_ROOT_global.h' |
---|
| 134 | INTEGER N |
---|
| 135 | KPP_REAL T, Told |
---|
| 136 | KPP_REAL Y(NVAR), P(NVAR) |
---|
| 137 | Told = TIME |
---|
| 138 | TIME = T |
---|
| 139 | CALL Update_SUN() |
---|
| 140 | CALL Update_RCONST() |
---|
| 141 | CALL Fun( Y, FIX, RCONST, P ) |
---|
| 142 | TIME = Told |
---|
| 143 | RETURN |
---|
| 144 | END |
---|
| 145 | |
---|
| 146 | |
---|
| 147 | SUBROUTINE JAC_CHEM(N, T, Y, J) |
---|
| 148 | INCLUDE 'KPP_ROOT_params.h' |
---|
| 149 | INCLUDE 'KPP_ROOT_global.h' |
---|
| 150 | INTEGER N |
---|
| 151 | KPP_REAL Told, T |
---|
| 152 | KPP_REAL Y(NVAR), J(LU_NONZERO) |
---|
| 153 | Told = TIME |
---|
| 154 | TIME = T |
---|
| 155 | CALL Update_SUN() |
---|
| 156 | CALL Update_RCONST() |
---|
| 157 | CALL Jac_SP( Y, FIX, RCONST, J ) |
---|
| 158 | TIME = Told |
---|
| 159 | RETURN |
---|
| 160 | END |
---|
| 161 | |
---|
| 162 | |
---|
| 163 | |
---|
| 164 | |
---|
| 165 | |
---|
| 166 | |
---|