chem/KPP/kpp/kpp-2.1/int/oldies/twostepj.f

   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