[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 | EXTERNAL ITER |
---|
| 11 | KPP_REAL T |
---|
| 12 | KPP_REAL V(NVAR), VOLD(NVAR), VNEW(NVAR) |
---|
| 13 | KPP_REAL startdt, hmin, hmax, h |
---|
| 14 | |
---|
| 15 | INTEGER INFO(5) |
---|
| 16 | |
---|
| 17 | INFO(1) = Autonomous |
---|
| 18 | h = hmin |
---|
| 19 | |
---|
| 20 | c Number of Jacobi-Seidel iterations |
---|
| 21 | numit = 3 |
---|
| 22 | |
---|
| 23 | |
---|
| 24 | DO i=1,NVAR |
---|
| 25 | RTOL(i) = 1.e-2 |
---|
| 26 | ENDDO |
---|
| 27 | |
---|
| 28 | CALL twostepj(NVAR,TIN,TOUT,h,hmin,hmax, |
---|
| 29 | + VOLD,VAR,VNEW, |
---|
| 30 | + ATOL,RTOL,numit, |
---|
| 31 | + nfcn,naccpt,nrejec,nstart,startdt,ITER) |
---|
| 32 | |
---|
| 33 | |
---|
| 34 | RETURN |
---|
| 35 | END |
---|
| 36 | |
---|
| 37 | |
---|
| 38 | |
---|
| 39 | SUBROUTINE ITER(n,T,y,yp,yl) |
---|
| 40 | INCLUDE 'KPP_ROOT_params.h' |
---|
| 41 | INCLUDE 'KPP_ROOT_global.h' |
---|
| 42 | REAL*8 T, y(NVAR), yp(NVAR), yl(NVAR) |
---|
| 43 | TOLD = TIME |
---|
| 44 | TIME = T |
---|
| 45 | CALL Update_SUN() |
---|
| 46 | CALL Update_RCONST() |
---|
| 47 | CALL FunSPLIT_VAR(y,Rad,yp,yl) |
---|
| 48 | TIME = TOLD |
---|
| 49 | RETURN |
---|
| 50 | END |
---|
| 51 | |
---|
| 52 | |
---|
| 53 | subroutine twostepj(n,t,te,dt,dtmin,dtmax, |
---|
| 54 | + yold,y,ynew, |
---|
| 55 | + atol,rtol,numit, |
---|
| 56 | + nfcn,naccpt,nrejec,nstart,startdt,ITER) |
---|
| 57 | implicit real*8 (a-h,o-z) |
---|
| 58 | external ITER |
---|
| 59 | integer n,numit,nfcn,naccpt,nrejec,nstart,i,j |
---|
| 60 | real*8 t,te,dt,dtmin,dtmax,startdt,ytol, |
---|
| 61 | + ratio,dtold,a1,a2,c,cp1,dtg,errlte,dy |
---|
| 62 | real*8 yold(n),y(n),ynew(n),yp(n),yl(n), |
---|
| 63 | + work(n),sum(n),atol(n),rtol(n) |
---|
| 64 | logical accept,failer,restart |
---|
| 65 | |
---|
| 66 | c |
---|
| 67 | c Initialization of counters, etc. |
---|
| 68 | |
---|
| 69 | naccpt=0 |
---|
| 70 | nrejec=0 |
---|
| 71 | nfcn=0 |
---|
| 72 | nstart=0 |
---|
| 73 | failer=.false. |
---|
| 74 | restart=.false. |
---|
| 75 | accept=.true. |
---|
| 76 | |
---|
| 77 | c Initial stepsize computation. |
---|
| 78 | |
---|
| 79 | 10 if (dtmin.eq.dtmax) then |
---|
| 80 | nstart=1 |
---|
| 81 | dt=min(dtmin,(te-t)/2) |
---|
| 82 | goto 28 |
---|
| 83 | endif |
---|
| 84 | CALL ITER(n,t,y,yp,yl) |
---|
| 85 | nfcn=nfcn+1 |
---|
| 86 | dt=te-t |
---|
| 87 | do 20 i=1,n |
---|
| 88 | ytol=atol(i)+rtol(i)*abs(y(i)) |
---|
| 89 | dy=yp(i)-y(i)*yl(i) |
---|
| 90 | if (dy.ne.0.0) dt=min(dt,ytol/abs(dy)) |
---|
| 91 | 20 continue |
---|
| 92 | 25 nstart=nstart+1 |
---|
| 93 | if (restart) dt=dt/10.0 |
---|
| 94 | restart=.true. |
---|
| 95 | dt=max(dtmin,min(dt,dtmax)) |
---|
| 96 | CALL FIT(t,te,dt) |
---|
| 97 | dt=min(dt,(te-t)/2) |
---|
| 98 | startdt=dt |
---|
| 99 | |
---|
| 100 | c The starting step is carried out, using the implicit Euler method. |
---|
| 101 | |
---|
| 102 | 28 do 30 i=1,n |
---|
| 103 | ynew(i)=y(i) |
---|
| 104 | yold(i)=y(i) |
---|
| 105 | sum(i)=y(i) |
---|
| 106 | 30 continue |
---|
| 107 | do 40 i=1,numit |
---|
| 108 | CALL ITER(n,t+dt,ynew,yp,yl) |
---|
| 109 | do i2i=1,n |
---|
| 110 | ynew(i2i) = (sum(i2i) + dt*yp(i2i))/(1.+dt*yl(i2i)) |
---|
| 111 | end do |
---|
| 112 | nfcn=nfcn+1 |
---|
| 113 | 40 continue |
---|
| 114 | naccpt=naccpt+1 |
---|
| 115 | t=t+dt |
---|
| 116 | do 50 j=1,n |
---|
| 117 | y(j)=ynew(j) |
---|
| 118 | 50 continue |
---|
| 119 | |
---|
| 120 | c Subsequent steps are carried out with the two-step BDF method. |
---|
| 121 | |
---|
| 122 | dtold=dt |
---|
| 123 | ratio=1.0 |
---|
| 124 | 60 continue |
---|
| 125 | c=1.0/ratio |
---|
| 126 | cp1=c+1.0 |
---|
| 127 | a1=((c+1.0)**2)/(c*c+2.0*c) |
---|
| 128 | a2=-1.0/(c*c+2.0*c) |
---|
| 129 | dtg=dt*(1.0+c)/(2.0+c) |
---|
| 130 | do 70 j=1,n |
---|
| 131 | sum(j)=a1*y(j)+a2*yold(j) |
---|
| 132 | ynew(j)=max(0.0,y(j)+ratio*(y(j)-yold(j))) |
---|
| 133 | 70 continue |
---|
| 134 | do 80 i=1,numit |
---|
| 135 | CALL ITER(n,t+dt,ynew,yp,yl) |
---|
| 136 | do i2i=1,n |
---|
| 137 | ynew(i2i) = (sum(i2i) + dtg*yp(i2i))/(1.+dtg*yl(i2i)) |
---|
| 138 | end do |
---|
| 139 | nfcn=nfcn+1 |
---|
| 140 | 80 continue |
---|
| 141 | |
---|
| 142 | c If stepsizes should remain equal, stepsize control is omitted. |
---|
| 143 | |
---|
| 144 | if (dtmin.eq.dtmax) then |
---|
| 145 | t=t+dtold |
---|
| 146 | naccpt=naccpt+1 |
---|
| 147 | do 85 j=1,n |
---|
| 148 | yold(j)=y(j) |
---|
| 149 | y(j)=ynew(j) |
---|
| 150 | 85 continue |
---|
| 151 | if (dt.ne.dtold) then |
---|
| 152 | t=t-dtold+dt |
---|
| 153 | goto 120 |
---|
| 154 | endif |
---|
| 155 | dt=min(dtold,te-t) |
---|
| 156 | ratio=dt/dtold |
---|
| 157 | if (t.ge.te) goto 120 |
---|
| 158 | goto 60 |
---|
| 159 | endif |
---|
| 160 | |
---|
| 161 | c Otherwise stepsize control is carried out. |
---|
| 162 | |
---|
| 163 | errlte=0.0 |
---|
| 164 | do 90 i=1,n |
---|
| 165 | ytol=atol(i)+rtol(i)*abs(y(i)) |
---|
| 166 | errlte=max(errlte,abs(c*ynew(i)-cp1*y(i)+yold(i))/ytol) |
---|
| 167 | 90 continue |
---|
| 168 | errlte=2.0*errlte/(c+c*c) |
---|
| 169 | CALL NEWDT(t,te,dt,dtold,ratio,errlte,accept, |
---|
| 170 | + dtmin,dtmax) |
---|
| 171 | |
---|
| 172 | c Here the step has been accepted. |
---|
| 173 | |
---|
| 174 | if (accept) then |
---|
| 175 | 201 format(2(E24.16,1X)) |
---|
| 176 | failer=.false. |
---|
| 177 | restart=.false. |
---|
| 178 | t=t+dtold |
---|
| 179 | naccpt=naccpt+1 |
---|
| 180 | do 100 j=1,n |
---|
| 181 | yold(j)=y(j) |
---|
| 182 | y(j)=ynew(j) |
---|
| 183 | 100 continue |
---|
| 184 | if (t.ge.te) goto 120 |
---|
| 185 | goto 60 |
---|
| 186 | endif |
---|
| 187 | |
---|
| 188 | c A restart check is carried out. |
---|
| 189 | |
---|
| 190 | if (failer) then |
---|
| 191 | nrejec=nrejec+1 |
---|
| 192 | failer=.false. |
---|
| 193 | naccpt=naccpt-1 |
---|
| 194 | t=t-dtold |
---|
| 195 | do 110 j=1,n |
---|
| 196 | y(j)=yold(j) |
---|
| 197 | 110 continue |
---|
| 198 | goto 25 |
---|
| 199 | endif |
---|
| 200 | |
---|
| 201 | c Here the step has been rejected. |
---|
| 202 | |
---|
| 203 | nrejec=nrejec+1 |
---|
| 204 | failer=.true. |
---|
| 205 | goto 60 |
---|
| 206 | |
---|
| 207 | c End of TWOSTEP. |
---|
| 208 | 120 end |
---|
| 209 | c===================================================================== |
---|
| 210 | |
---|
| 211 | subroutine NEWDT(t,te,dt,dtold,ratio,errlte, |
---|
| 212 | + accept,dtmin,dtmax) |
---|
| 213 | real*8 t,te,dt,dtold,ratio,errlte,ts,dtmin,dtmax |
---|
| 214 | logical accept |
---|
| 215 | if (errlte.gt.1.0.and.dt.gt.dtmin) then |
---|
| 216 | accept=.false. |
---|
| 217 | ts=t |
---|
| 218 | else |
---|
| 219 | accept=.true. |
---|
| 220 | dtold=dt |
---|
| 221 | ts=t+dtold |
---|
| 222 | endif |
---|
| 223 | dt=max(0.5,min(2.0,0.8/sqrt(errlte)))*dt |
---|
| 224 | dt=max(dtmin,min(dt,dtmax)) |
---|
| 225 | CALL FIT(ts,te,dt) |
---|
| 226 | ratio=dt/dtold |
---|
| 227 | end |
---|
| 228 | |
---|
| 229 | subroutine FIT(t,te,dt) |
---|
| 230 | real*8 t,te,dt,rns |
---|
| 231 | integer ns |
---|
| 232 | rns=(te-t)/dt |
---|
| 233 | if (rns.gt.10.0) goto 10 |
---|
| 234 | ns=int(rns)+1 |
---|
| 235 | dt=(te-t)/ns |
---|
| 236 | dt=(dt+t)-t |
---|
| 237 | 10 return |
---|
| 238 | end |
---|
| 239 | |
---|
| 240 | |
---|
| 241 | |
---|
| 242 | C End of MAIN function |
---|
| 243 | C **************************************************************** |
---|
| 244 | |
---|