c5ca3916411862ba9abe4fc23d815c0098f4a4f8
[scilab.git] / scilab / modules / differential_equations / sci_gateway / fortran / sci_f_ode.f
1 c     Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 c     Copyright (C) INRIA
3 c     ...
4 c     
5 c     This file must be used under the terms of the CeCILL.
6 c     This source file is licensed as described in the file COPYING,
7 c     which
8 c     you should have received as part of this distribution.  The terms
9 c     are also available at
10 c     http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
11 c     
12       subroutine sciode
13 c     
14 c     ode
15 c     
16       include 'stack.h'
17       integer iadr,sadr
18 c      
19       character buftmp*(bsiz)
20 c     
21 c     common de lsode,lsoda,lsodar
22       double precision xxxx,yyyy,rlsr
23       integer ilsr
24       common/ls0001/xxxx(219),iiii(39)
25       common/lsa001/yyyy(22),jjjj(9)
26       common/lsr001/ rlsr(5),ilsr(9)
27       common/eh0001/kkkk(2)
28       save /ls0001/,/lsa001/,/lsr001/,/eh0001/
29 c     
30 c     commons avec bydot,bjac,....
31 c     
32       character*(nlgh+1) namef,namej,names
33       common/cydot/namef
34       common/cjac/namej
35       common/csurf/names
36       integer       iero
37       common/ierajf/iero
38 c     
39       double precision atol,rtol,t0,tout,dir
40       double precision h0,hmax,hmin,tcrit,tmax
41       integer top1,top2,tope,hsize
42 c     meth is simulator number, and jactyp the jacobian type used
43       integer meth,jactyp
44       logical jaco,achaud,withw,single
45       external bydot,bjac,bsurf
46       integer raide(2),root(2),adams,discre,rkf(3),rk(2),fix(2)
47       integer params(nsiz)
48       data raide/28,29/,adams/10/,root/27,24/
49       data discre/13/,rkf/27,20,15/,rk/27,20/,fix/15,18/
50       data params/-235739080,-303896856,669247720,3*673720360/
51 c     
52       iadr(l)=l+l-1
53       sadr(l)=(l/2)+1
54
55       iflagcr=0
56 c     
57 c     get %ODEOPTIONS variable
58       ifin=fin
59       fin=-1
60       istate=1
61       call stackg(params)
62       if(fin.eq.0) then
63 c     call msgs(72,0)
64          iopt=0
65          itask=1
66          jactyp=2
67          ml=-1
68          mu=-1
69       else
70          iopt=1
71          il=iadr(lstk(fin))
72          l=sadr(il+4)
73 c     %ODEOPTIONS=[itask,tcrit,h0,hmax,hmin,jactyp,mxstep,..
74 c     mxordn,mxords,ixpr, ml,mu]
75          itask=int(stk(l))
76          if(itask.lt.1.or.itask.gt.5) then
77             buf=' invalid option (first entry)'
78             call error(9999)
79             return
80          endif
81          tcrit=stk(l+1)
82          h0=stk(l+2)
83          hmax=stk(l+3)
84          hmin=stk(l+4)
85          if(hmin.gt.hmax) then
86             buf=' what? hmin greater than hmax?'
87             call error(9999)
88             return
89          endif
90          jactyp=int(stk(l+5))
91          if(jactyp.lt.0.or.jactyp.gt.5) then
92             buf=' invalid option (entry no 4)'
93             call error(9999)
94             return
95          endif
96          mxstep=int(stk(l+6))
97          mxordn=int(stk(l+7))
98          mxords=int(stk(l+8))
99          ixpr=int(stk(l+9))
100          ml=int(stk(l+10))
101          mu=int(stk(l+11))
102       endif
103 c     .....
104 c     
105       fin=ifin
106       withw=.false.
107
108       tope=top
109       itype=0
110
111       if(rhs.lt.4) then
112          call error(39)
113          return
114       endif
115 c     
116 c     lw=premiere adresse libre dans la pile
117       lw=lstk(top+1)
118 c     
119 c     test demarrage a chaud
120       ifin=iadr(lstk(top))
121       achaud=istk(ifin).eq.1
122       if(achaud) then
123 c     ilc=adresse of lsod* integer work space
124          top=top-2
125          il=iadr(lstk(top+2))
126          if(istk(il).ne.1) then
127             err=rhs
128             call error(53)
129             return
130          endif
131          liwp=istk(il+2)*istk(il+1)
132          lci=sadr(il+4)
133          ilc=iadr(lci)
134 c     lc=adresse of lsod* real work space
135          il=iadr(lstk(top+1))
136          if(istk(il).ne.1) then
137             err=rhs-1
138             call error(53)
139             return
140          endif
141          lc=sadr(il+4)
142          lrwp=istk(il+1)*istk(il+2)
143       endif
144 c     
145       top2=top-rhs+1
146       if(achaud) top2=top2+2
147       ile=iadr(lstk(top2))
148 c     
149       if(istk(ile).eq.10) then
150          top2=top2+1
151          if(abs(istk(ile+6)).eq.adams) then
152 c     lsode (adams)
153             meth=1
154          elseif(abs(istk(ile+6)).eq.raide(1) .and.
155      $           abs(istk(ile+7)).eq.raide(2)) then
156 c     lsode (gear)
157             meth=2
158          elseif(abs(istk(ile+6)).eq.root(1) .and.
159      $           abs(istk(ile+7)).eq.root(2)) then
160 c     lsodar
161             meth=3
162          elseif(abs(istk(ile+6)).eq.discre) then
163 c     ldisc
164             meth=4
165          elseif(abs(istk(ile+6)).eq.rkf(1) .and.
166      $           abs(istk(ile+7)).eq.rkf(2) .and.
167      $           abs(istk(ile+8)).eq.rkf(3)) then
168 c     rkf45
169             meth=6
170          elseif(abs(istk(ile+6)).eq.rk(1) .and.
171      $           abs(istk(ile+7)).eq.rk(2)) then
172 c     rk4
173             meth=5
174          elseif(abs(istk(ile+6)).eq.fix(1) .and.
175      $           abs(istk(ile+7)).eq.fix(2)) then
176 c     rksimp
177             meth=7
178          else
179             call error(42)
180             return
181          endif
182       else
183 c     lsoda
184          meth=0
185       endif
186 c     
187       if(meth.lt.3) then
188          if(lhs.ne.3.and.lhs.ne.1) then
189             call error(41)
190             return
191          endif
192       elseif(meth.eq.3) then
193          if(lhs.eq.3.or.lhs.gt.4) then
194             call error(41)
195             return
196          endif
197       elseif(meth.ge.4) then
198          if(lhs.ne.1) then
199             call error(41)
200             return
201          endif
202       endif
203 c     
204       top1=top
205 c     
206       if(meth.eq.3) then
207 c     on recupere le simulateur des equations des surfaces
208          ilsurf=iadr(lstk(top1))
209          if(istk(ilsurf).ne.10.and.istk(ilsurf).ne.15.and.
210      $        istk(ilsurf).ne.11.and.istk(ilsurf).ne.13) then
211             err=rhs-(tope-top1)
212             call error(80)
213             return
214          endif
215          if(istk(ilsurf).eq.10) then
216             names=' '
217             call cvstr(istk(ilsurf+5)-1,istk(ilsurf+6),names,1)
218             names(istk(ilsurf+5):istk(ilsurf+5))=char(0)
219             call setfsurf(names,irep)
220             if ( irep.eq.1) then
221                buf = names
222                call error(50)
223                return
224             endif
225          elseif(istk(ilsurf).eq.15) then
226             le1=sadr(ilsurf+istk(ilsurf+1)+3)
227             if(istk(iadr(le1)).ne.11.and.istk(iadr(le1)).ne.13) then
228                err=rhs-(tope-top1)
229                call error(80)
230                return
231             endif
232          endif
233          ksurf=top1
234          top1=top1-1
235 c     ... et le nombre d'equations
236          il=iadr(lstk(top1))
237          if(istk(il).ne.1) then
238             err=rhs-(tope-top1)
239             call error(53)
240             return
241          endif
242          if(istk(il+1)*istk(il+2).ne.1) then
243             err=rhs-(tope-top1)
244             call error(89)
245             return
246          endif
247          l=sadr(il+4)
248          nsurf=stk(l)
249          top1=top1-1
250       else
251          ksurf=0
252       endif
253
254       il=iadr(lstk(top1-1))
255       if(istk(il).eq.10.or.istk(il).eq.15.or.
256      $     istk(il).eq.11.or.istk(il).eq.13) then
257 c     JACOBIAN IS GIVEN (variable top1)
258          ilj=iadr(lstk(top1))
259          islj=istk(ilj)
260          if(islj.lt.10.or.islj.gt.15.or.islj.eq.12) then
261             err=rhs-(tope-top1)
262             call error(80)
263             return
264          endif
265          if(islj.eq.10) then
266             namej=' '
267             call cvstr(istk(ilj+5)-1,istk(ilj+6),namej,1)
268             namej(istk(ilj+5):istk(ilj+5))=char(0)
269             call setfjac(namej,irep)
270             if ( irep.eq.1) then
271                buf = namej
272                call error(50)
273                return
274             endif
275          elseif(islj.eq.15) then
276             le1=sadr(ilj+istk(ilj+1)+3)
277             if(istk(iadr(le1)).ne.11.and.istk(iadr(le1)).ne.13) then
278                err=rhs-(tope-top1)
279                call error(80)
280                return
281             endif
282          endif
283          if (meth.ge.4) then
284             call msgs(75,0)
285          endif
286          if(jactyp.eq.0.and.(meth.eq.2.or.meth.eq.1)) then
287             call msgs(75,0)
288          endif
289          jaco=.true.
290          if(iopt.eq.0) then
291 c     set jactyp (jacobian is supposed full)
292             jactyp=1
293          else
294 c     check jactyp
295             if(jactyp.eq.2.or.jactyp.eq.5) then
296                call msgs(75,0)
297             endif
298          endif
299          kjac=top1
300          top1=top1-1
301       else
302 c     JACOBIAN NOT GIVEN
303          if(iopt.eq.0) then
304 c     set jactyp (estimated jacobian is supposed full)
305             jactyp=2
306          else
307 c     check jactyp
308             if(jactyp.eq.1.or.jactyp.eq.4) then
309 c     %ODEOPTIONS requires the jacobian
310                call msgs(76,0)
311                jactyp=jactyp+1
312             endif
313          endif
314          jaco=.false.
315          kjac=0
316       endif
317
318       kytop=top1
319
320 c     
321 c     rhs
322       ilf=iadr(lstk(top1))
323       islf=istk(ilf)
324       if(islf.ne.10.and.islf.ne.15.and.islf.ne.11.and.islf.ne.13) then
325          err=rhs-(tope-top1)
326          call error(80)
327          return
328       endif
329       kydot=top1
330       if(islf.eq.10) then
331          namef=' '
332          call cvstr(istk(ilf+5)-1,istk(ilf+6),namef,1)
333          namef(istk(ilf+5):istk(ilf+5))=char(0)
334          call setfydot(namef,irep)
335          if ( irep.eq.1) then
336             buf = namef
337             call error(50)
338             return
339          endif
340 c     test list('fex', ...) or list(fun,...)
341       elseif(islf.eq.15) then
342          le1=sadr(ilf+istk(ilf+1)+3)
343          if(istk(iadr(le1)).eq.10) then
344             withw=.true.
345 c     .     next line just to tell to bydot that external is in fortran
346             istk(ilf)=10
347             if(istk(ilf+1).ne.2) then
348                buf='wrong list passed: needs two elts in list'
349                call error(9999)
350                return
351             endif
352             long1=istk(ilf+3)
353             lf=lstk(top1)
354             illl=iadr(lf+istk(ilf+3))
355             nblett=istk(illl-1)-1
356             namef=' '
357             call cvstr(istk(ilf+11)-1,istk(ilf+12),namef,1)
358             namef(istk(ilf+11):istk(ilf+11))=char(0)
359             call setfydot(namef,irep)
360             if ( irep.eq.1) then
361                buf = namef
362                call error(50)
363                return
364             endif
365             ll1=sadr(ilf+5)
366             ll2=ll1+long1-1
367             il2=iadr(ll2)
368             nbw=istk(il2+1)*istk(il2+2)
369             if(istk(il2+3).ne.0) then
370                buf='working array must be real'
371                call error(9999)
372                return
373             endif
374             lww=sadr(il2+4)
375 c     .     lww = adr w , nbw = size (w)
376          elseif(istk(iadr(le1)).ne.11.and.istk(iadr(le1)).ne.13) then
377             err=rhs-(tope-top1)
378             call error(80)
379             return
380          endif
381       endif
382 c     
383 c     jaco,type and meth initialized ...
384 c     top2 point on y0
385 c     
386 c     y0
387       kynew=top2
388       il=iadr(lstk(top2))
389       it=istk(il+3)
390 c     
391       if(istk(il).eq.1) then
392          hsize=4
393          ny=istk(il+1)*istk(il+2)*(istk(il+3)+1)
394       elseif(istk(il).eq.2) then
395          mn=istk(il+1)*istk(il+2)
396          hsize=9+mn
397          ny=(istk(il+8+mn)-1)*(istk(il+3)+1)
398       else
399          err=rhs-(tope-top2)
400          call error(44)
401          return
402       endif
403       if(it.eq.1) nys2=ny/2
404       ly=sadr(il+hsize)
405
406 c     list('fex',w)
407       if(withw) then
408          if(top+1.ge.bot) then
409             call error(18)
410             return
411          endif
412          top=top+1
413 c     .  kynew=top
414          ily=iadr(lstk(top))
415          err=sadr(ily+4)+ny+nbw-lstk(bot)
416          if(err.gt.0) then
417             call error(17)
418             return
419          endif
420          istk(ily+1)=1
421          istk(ily+2)=ny+nbw
422          istk(ily+3)=1
423          istk(ily+4)=0
424          ly1=sadr(ily+4)
425          call unsfdcopy(ny,stk(ly),1,stk(ly1),1)
426          call unsfdcopy(nbw,stk(lww),1,stk(ly1+ny),1)
427          lstk(top+1)=ly1+ny+nbw
428          ly=ly1
429          lw=lstk(top+1)
430       endif
431       lw1=lw
432       lw=sadr(iadr(lw1)+13)
433 c     
434 c     t0
435       top2=top2+1
436       kttop=top2
437       il=iadr(lstk(top2))
438       if(istk(il).ne.1) then
439          err=rhs-(tope-top2)
440          call error(53)
441          return
442       endif
443       l=sadr(il+4)
444       t0=stk(l)
445 c     t1
446       top2=top2+1
447       il=iadr(lstk(top2))
448       if(istk(il).ne.1) then
449          err=rhs-(tope-top2)
450          call error(53)
451          return
452       endif
453
454 c     number of output points
455       nn=istk(il+1)*istk(il+2)
456 c     pointer on  output time vector
457       lt1=sadr(il+4)
458 c     
459 c     optionnal parameters rtol et atol
460       top2=top2+1
461 c     default values
462       if(meth.eq.6.or.meth.eq.7) then
463          rtol=1.d-3
464          atol=1.d-4
465       else
466          rtol=1.0d-7
467          atol=1.0d-9
468       endif
469       nr=1
470       na=1
471       jobtol=kytop-top2+1
472 c     jobtol=(nothing ,rtol only,rtol and atol)
473 c     
474       if(jobtol.eq.1) then
475 c     default tolerances
476          lr=lw
477          la=lr+1
478          stk(la)=atol
479          stk(lr)=rtol
480       else
481 c     rtol given
482          lr=lw
483 c     rtol
484          il=iadr(lstk(top2))
485          if(istk(il).ne.1) then
486             err=rhs-(tope-top2)
487             call error(53)
488             return
489          endif
490          nr=istk(il+1)*istk(il+2)
491          if(nr.ne.1.and.nr.ne.ny) then
492             err=rhs-(tope-top2)
493             call error(89)
494             return
495          endif
496          lrt=sadr(il+4)
497          call unsfdcopy(nr,stk(lrt),1,stk(lr),1)
498          la=lr+nr
499 c     .  atol
500          if(jobtol.eq.2) then
501 c     .     default
502             stk(la)=atol
503          else
504 c     .     atol given
505             top2=top2+1
506             il=iadr(lstk(top2))
507             if(istk(il).ne.1) then
508                err=rhs-(tope-top2)
509                call error(53)
510                return
511             endif
512             na=istk(il+1)*istk(il+2)
513             if(na.ne.1.and.na.ne.ny) then
514                err=rhs-(tope-top2)
515                call error(89)
516                return
517             endif
518             lat=sadr(il+4)
519             call unsfdcopy(na,stk(lat),1,stk(la),1)
520          endif
521       endif
522       lw=la+na
523
524 c     set input top value
525       if(achaud) top=top+2
526 c     
527       if(nr.eq.1.and.na.eq.1) itol=1
528       if(nr.eq.1.and.na.gt.1) itol=2
529       if(nr.gt.1.and.na.eq.1) itol=3
530       if(nr.gt.1.and.na.gt.1) itol=4
531 c     
532 c     compute integrator workspace  sizes
533       if(meth.eq.0) then
534 c     lsoda
535          lrw=22+ny*max(16,ny+9)
536          liw=20+ny
537          nsizd=241
538          nsizi=50
539          if(jactyp.eq.4.or.jactyp.eq.5) then
540             lrn=20+16*ny
541             lrs=22+10*ny+(2*ml+mu)*ny
542             lrw=max(lrn,lrs)
543          endif
544       elseif(meth.eq.1) then
545 c     lsode - adams
546          if(jactyp.eq.1.or.jactyp.eq.2) then
547             lrw=22+16*ny+ny*ny
548          elseif(jactyp.eq.4.or.jactyp.eq.5) then
549             lrw=22+16*ny+(2*ml+mu+1)*ny
550          else
551             lrw=20+16*ny
552          endif
553          liw=20+ny
554          nsizd=219
555          nsizi=41
556       elseif(meth.eq.2) then
557 c     lsode gear
558          if(jactyp.eq.1.or.jactyp.eq.2) then
559             lrw=22+9*ny+ny*ny
560          elseif(jactyp.eq.4.or.jactyp.eq.5) then
561             lrw=22+9*ny+(2*ml+mu+1)*ny
562          else
563             lrw=20+9*ny
564          endif
565          liw=20+ny
566          nsizd=219
567          nsizi=41
568       elseif(meth.eq.3) then
569 c     lsodar
570          ilroot=iadr(lw)
571          lw=sadr(ilroot+nsurf)
572          lrw= 22 + ny * max(16, ny + 9) + 3*nsurf
573          liw=20+ny
574          nsizd=246
575          nsizi=59
576       elseif(meth.eq.4) then
577 c     lsdisc
578          lrw=ny
579          liw=1
580       elseif(meth.eq.5) then
581 c     lsrgk
582          lrw=3*ny
583          liw=1
584       elseif(meth.eq.6) then
585 c     rkf45
586          lrw=3+8*ny
587          liw=5
588       elseif(meth.eq.7) then
589 c     rksimp
590          lrw=3+8*ny
591          liw=1
592       endif
593 c     
594 c     hot start
595 c     
596       if(achaud) then
597          if(meth.le.3) then
598 c     check for input hot start table consistency
599             if(lrwp.ne.lrw+nsizd) then
600                buf=' Real hot start table has incorrect size'
601                call error(9999)
602                return
603             endif
604             if(liwp.ne.liw+nsizi) then
605                buf=' Input hot start table has incorrect size'
606                call error(9999)
607                return
608             endif
609          endif
610          istate=2
611 c     commons retrieval from hot start tables
612          if(meth.eq.0) then
613 c     lsoda
614             lsavs=lc+lrwp-nsizd
615             lsavi=lci+liwp-nsizi
616             call rscma1(stk(lsavs),stk(lsavi))
617          elseif(meth.eq.1.or.meth.eq.2) then
618 c     lsode
619             lsavs=lc+lrwp-nsizd
620             lsavi=lci+liwp-nsizi
621             call rscom1(stk(lsavs),stk(lsavi))
622          elseif(meth.eq.3) then
623 c     lsodar
624             lsavs=lc+lrwp-nsizd
625             lsavi=lci+liwp-nsizi
626             call rscar1(stk(lsavs),stk(lsavi))
627          endif
628 c     integer workspace retrieval
629          do 40 k=1,liw
630             istk(ilc+k-1)=int(stk(lci+k-1))
631  40      continue
632       endif
633 c     
634 c     
635 c     compute pointer on ode real and integer work spaces
636       lc0=lw
637       li=lc0+lrw
638 c     
639       ili=iadr(li)
640       lw=sadr(ili+liw)
641 c     
642 c     get memory to store results
643       lyp=lw
644       if(itask.eq.2.or.itask.eq.3.or.itask.eq.5) then
645 c     unknown number of output points.  space for  points
646 c     will be allocated later
647          single=.true.
648          lw=lyp
649          if(nn.ne.1) then
650             call msgs(77,0)
651             stk(lt1)=stk(lt1+nn-1)
652             nn=1
653          endif
654          if(it.ne.0) then
655             buf='itask=2,3 or 5: y0 must be a real vector'
656             call error(9999)
657             return
658          endif
659       else
660 c     number of output points is equal to number of t points all
661 c     space allocated here
662          single=.false.
663          lw=lw+nn*ny
664       endif
665
666 c     top points on external workspace
667       top=top+1
668       lstk(top+1)=lw
669       err=lstk(top+1)-lstk(bot)
670       if(err.gt.0) then
671          call error(17)
672          return
673       endif
674 c     
675       call xsetf(1)
676       call xsetun(wte)
677 c     
678       if(.not.achaud) then
679          lc=lc0
680          ilc=ili
681       endif
682 c     
683 c     data structure passed to externals, it contains pointer
684 c     to externals parameters
685 c     
686       ilw1=iadr(lw1)
687       istk(ilw1)=3
688       istk(ilw1+1)=ilw1+4
689       istk(ilw1+2)=ilw1+7
690       istk(ilw1+3)=ilw1+10
691       istk(ilw1+4)=kydot
692       istk(ilw1+5)=kttop
693       istk(ilw1+6)=kynew
694       istk(ilw1+7)=kjac
695       istk(ilw1+8)=kttop
696       istk(ilw1+9)=kynew
697       istk(ilw1+10)=ksurf
698       istk(ilw1+11)=kttop
699       istk(ilw1+12)=kynew
700 c     
701       if(iopt.eq.1) then
702 c     copy integration options in lsod* workspace
703          if(itask.ge.4) then
704             stk(lc)=tcrit
705          endif
706          stk(lc+4)=h0
707          stk(lc+5)=hmax
708          stk(lc+6)=hmin
709          if(meth.eq.0.or.meth.eq.3) then
710 c     lsoda/lsodar
711             if(jactyp.eq.4.or.jactyp.eq.5) then
712                istk(ilc)=ml
713                istk(ilc+1)=mu
714             endif
715             istk(ilc+4)=ixpr
716             istk(ilc+5)=mxstep
717             istk(ilc+6)=0
718             istk(ilc+7)=mxordn
719             istk(ilc+8)=mxords
720          elseif(meth.lt.3) then
721 c     lsode
722             if(jactyp.eq.4.or.jactyp.eq.5) then
723                istk(ilc)=ml
724                istk(ilc+1)=mu
725             endif
726             if(meth.lt.2) then
727                istk(ilc+4)=mxordn
728             else
729                istk(ilc+4)=mxords
730             endif
731             istk(ilc+5)=mxstep
732             istk(ilc+6)=0
733          endif
734       endif
735       tmax=stk(lt1+nn-1)
736       niter=nn
737       if(ixpr.eq.1.and.iopt.eq.1) then
738          call writebufodea(buf,itask,meth,jactyp,ml,mu,iopt)
739          call basout(io,wte,buf(1:80))
740          call writebufodeb(buf,tcrit,stk(lc+4),stk(lc+5),stk(lc+6))
741          call basout(io,wte,buf(1:80))
742       endif
743       if(single) then
744 c     loop til t=tout
745 c     --------------
746          dir=sign(1.0D0,tmax-t0)
747          tout=tmax
748          k=0
749          nn=0
750  50      k=k+1
751          nn=nn+1
752          if(meth.eq.0) then
753             call lsoda(bydot,ny,stk(ly),t0,tout,itol,stk(lr),
754      1           stk(la),itask,istate,iopt,stk(lc),lrw,istk(ilc),
755      2           liw,bjac,jactyp)
756          elseif(meth.eq.1.or.meth.eq.2) then
757             call lsode(bydot,ny,stk(ly),t0,tout,itol,stk(lr),
758      1           stk(la),itask,istate,iopt,stk(lc),lrw,istk(ilc),
759      2           liw,bjac,10*meth+jactyp)
760          elseif(meth.eq.3) then
761             call lsodar(bydot,ny,stk(ly),t0,tout,itol,stk(lr),
762      1           stk(la),itask,istate,iopt,stk(lc),lrw,istk(ilc),
763      2           liw,bjac,jactyp,bsurf,nsurf,istk(ilroot))
764          elseif(meth.eq.4) then
765             call lsdisc(bydot,ny,stk(ly),t0,tout, stk(lc),lrw,
766      1           istate)
767          elseif(meth.eq.5) then
768             call lsrgk(bydot,ny,stk(ly),t0,tout,itol,stk(lr),
769      1           stk(la),itask,istate,iopt,stk(lc),lrw,istk(ilc),
770      2           liw,bjac,meth)
771             if(iero.eq.-1) then
772                write(buftmp,'(e10.3)') tout
773                buf = buftmp
774                call msgs(70,0)
775             endif
776          elseif(meth.eq.6) then
777             call rkf45(bydot,ny,stk(ly),t0,tout,itol,rtol,
778      1           atol,itask,istate,iopt,stk(lc),lrw,istk(ilc),
779      2           liw,bjac,meth)
780          elseif(meth.eq.7) then
781             call rksimp(bydot,ny,stk(ly),t0,tout,itol,rtol,
782      1           atol,itask,istate,iopt,stk(lc),lrw,istk(ilc),
783      2           liw,bjac,meth)
784          endif
785          if(err.gt.0.or.err1.gt.0) return
786          if(istate.lt.0) then
787             if(meth.le.3) then
788                if(istate.eq.-3) then
789                   buf='illegal input'
790                   call error(9999)
791                   return
792                endif
793             elseif(meth.eq.5) then
794                call msgs(71,0)
795             elseif(meth.eq.4) then
796                if(istate.eq.-3) then
797                   buf ='ode discrete : a requested k is smaller '
798      $                 // ' than initial one'
799                   call error(999)
800                   return
801                else
802                   return
803                endif
804             endif
805             call msgs(4,ierr)
806             nn=k-1
807             goto 500
808          endif
809          if((meth.eq.6.or.meth.eq.7).and.istate.ne.2) then
810             nn=k-1
811             call msgs(74,0)
812             goto 500
813          endif
814 c     store intermediate result
815          lys=lyp+(k-1)*(ny+1)
816          lstk(top+1)=lys+(ny+1)
817          err=lstk(top+1)-lstk(bot)
818          if(err.gt.0) then
819             call error(17)
820             return
821          endif
822          stk(lys)=t0
823          call unsfdcopy(ny,stk(ly),1,stk(lys+1),1)
824          if((t0-tout)*dir.ge.0.0d0)  then
825 c     tout reached
826             nn=k
827             goto 500
828          endif
829          if(meth.eq.3.and.istate.eq.3) then
830 c     lsodar: a root found
831             nn=k
832             goto 500
833          endif
834          if(itask.eq.4.and.iflagcr.eq.1) then
835 c     tcrit reached
836             nn=k
837
838             call msgs(73,0)
839             goto 500
840          endif
841          goto 50
842 c     
843       else
844 c     
845 c     loop on t points
846 c--------------------
847          do 60 k=1,niter
848             tout=stk(lt1+k-1)
849             if(itask.ge.4.and.tout.gt.tcrit) then
850                tout=tcrit
851                iflagcr=1
852             endif
853             if(meth.eq.0) then
854                call lsoda(bydot,ny,stk(ly),t0,tout,itol,stk(lr),
855      1              stk(la),itask,istate,iopt,stk(lc),lrw,istk(ilc),
856      2              liw,bjac,jactyp)
857             elseif(meth.eq.1.or.meth.eq.2) then
858                call lsode(bydot,ny,stk(ly),t0,tout,itol,stk(lr),
859      1              stk(la),itask,istate,iopt,stk(lc),lrw,istk(ilc),
860      2              liw,bjac,10*meth+jactyp)
861             elseif(meth.eq.3) then
862                call lsodar(bydot,ny,stk(ly),t0,tout,itol,stk(lr),
863      1              stk(la),itask,istate,iopt,stk(lc),lrw,istk(ilc),
864      2              liw,bjac,jactyp,bsurf,nsurf,istk(ilroot))
865             elseif(meth.eq.4) then
866                call lsdisc(bydot,ny,stk(ly),t0,tout, stk(lc),lrw,
867      1              istate)
868             elseif(meth.eq.5) then
869                call lsrgk(bydot,ny,stk(ly),t0,tout,itol,stk(lr),
870      1              stk(la),itask,istate,iopt,stk(lc),lrw,istk(ilc),
871      2              liw,bjac,meth)
872                if(iero.eq.-1) then
873                   write(buftmp,'(e10.3)') tout
874                   buf = buftmp
875                   call msgs(70,0)
876                endif
877             elseif(meth.eq.6) then
878                call rkf45(bydot,ny,stk(ly),t0,tout,itol,rtol,
879      1              atol,itask,istate,iopt,stk(lc),lrw,istk(ilc),
880      2              liw,bjac,meth)
881                istk(ilc)=0
882                istk(ilc+1)=0
883             elseif(meth.eq.7) then
884                call rksimp(bydot,ny,stk(ly),t0,tout,itol,rtol,
885      1              atol,itask,istate,iopt,stk(lc),lrw,istk(ilc),
886      2              liw,bjac,meth)
887             endif
888             if(err.gt.0.or.err1.gt.0) return
889
890             if(istate.lt.0) then
891                if(meth.le.3) then
892                   if(istate.eq.-3) then
893                      buf='illegal input'
894                      call error(9999)
895                      return
896                   endif
897                endif
898                if(meth.eq.5) call msgs(71,0)
899                call msgs(4,ierr)
900                nn=k-1
901                goto 500
902             endif
903             if((meth.eq.6.or.meth.eq.7).and.istate.ne.2) then
904                nn=k-1
905                call msgs(74,0)
906                goto 500
907             endif
908 c     store intermediate result
909             if(it.eq.0) then
910                lys=lyp+(k-1)*ny
911                call unsfdcopy(ny,stk(ly),1,stk(lys),1)
912             else
913                lys=lyp+(k-1)*nys2
914                call unsfdcopy(nys2,stk(ly),1,stk(lys),1)
915                call unsfdcopy(nys2,stk(ly+nys2),1,stk(lys+nn*nys2),1)
916             endif
917             if(meth.eq.3.and.istate.eq.3) then
918 c     lsodar: a root found
919                nn=k
920                goto 500
921             endif
922             if(itask.eq.4.and.iflagcr.eq.1) then
923 c     tcrit reached
924                nn=k
925                call msgs(73,0)
926                goto 500
927             endif
928  60      continue
929       endif
930 c     
931 c     form results for output
932  500  continue
933 c      if(lhs.ge.3) then
934 c     preserve lsod* working spaces
935 c         lw=lyp+nn*(ny+1)
936 c         ilw=iadr(lw+lrw)
937 c         err=sadr(ilw+liw)-lstk(bot)
938 c         if(err.gt.0) then
939 c            call error(17)
940 c            return
941 c         endif
942 c        call unsfdcopy(lrw,stk(lc),1,stk(lw),1)
943 c        call icopy(liw,istk(ilc),1,istk(ilw),1)
944 c     endif
945 c     form state output
946       ils=iadr(lstk(kynew))
947       top=tope-rhs+1
948       call icopy(hsize,istk(ils),1,istk(ile),1)
949       ly=sadr(ile+hsize)
950       nel=istk(ile+1)*istk(ile+2)
951       if(nn.eq.0) then
952          istk(ile)=1
953          istk(ile+1)=0
954          istk(ile+2)=0
955          istk(ile+3)=0
956          lstk(top+1)=sadr(ile+4)
957       elseif(single) then
958          istk(ile+1)=istk(ile+1)+1
959          istk(ile+2)=nn*istk(ile+2)
960          inc=1
961          if(ly.gt.lyp) inc=-1
962          call unsfdcopy((ny+1)*nn,stk(lyp),inc,stk(ly),inc)
963          lstk(top+1)=ly+(ny+1)*nn
964       else
965          istk(ile+2)=nn*istk(ile+2)
966          if(istk(ile).eq.2) ly=sadr(ile+9+nel*nn)
967          inc=1
968          if(ly.gt.lyp) inc=-1
969
970          if(meth.eq.3.and.iadr(ly+ny*nn).gt.ilroot) then
971 c     preserve jroot table
972             ilr1=iadr(lyp+ny*nn)
973             err=sadr(ilr1)-lstk(bot)
974             if(err.gt.0) then
975                call error(17)
976                return
977             endif
978             call icopy(nsurf,istk(ilroot),1,istk(ilr1),1)
979             ilroot=ilr1
980          endif
981
982          call unsfdcopy(ny*nn,stk(lyp),inc,stk(ly),inc)
983          lstk(top+1)=ly+ny*nn
984          if(istk(ile).eq.2) then
985 c     on defini la table des pointeurs
986             il=ile+7
987             do 502 i=2,nn
988                do 501 j=1,nel
989                   istk(il+nel+j+1)=istk(ile+8+j)-1+istk(il+nel+1)
990  501           continue
991                il=il+nel
992  502        continue
993          endif
994       endif
995
996       if(meth.eq.3) then
997          if(lhs.lt.2) return
998 c     lsodar: form roots output
999          top=top+1
1000          il=iadr(lstk(top))
1001          istk(il)=1
1002          istk(il+3)=0
1003          l=sadr(il+4)
1004          if(istate.eq.3) then
1005             istk(il+1)=1
1006             istk(il+2)=1
1007             stk(l)=t0
1008             do 503 i=0,nsurf-1
1009                if(istk(ilroot+i).ne.0) then
1010                   l=l+1
1011                   istk(il+2)=istk(il+2)+1
1012                   stk(l)=i+1
1013                endif
1014  503        continue
1015          else
1016             istk(il+1)=0
1017             istk(il+2)=0
1018          endif
1019          lstk(top+1)=l+1
1020       endif
1021 c     form w and iw output
1022       if(lhs.lt.3) return
1023 c     w
1024       top=top+1
1025       il=iadr(lstk(top))
1026       istk(il)=1
1027       istk(il+1)=1
1028       istk(il+2)=lrw+nsizd
1029       istk(il+3)=0
1030       l=sadr(il+4)
1031       lstk(top+1)=l+lrw+nsizd
1032       call unsfdcopy(lrw,stk(lc),1,stk(l),1)
1033       lsvs=l+lrw
1034 c     iw
1035       top=top+1
1036       il=iadr(lstk(top))
1037       istk(il)=1
1038       istk(il+1)=1
1039       istk(il+2)=liw+nsizi
1040       istk(il+3)=0
1041       l=sadr(il+4)
1042       lstk(top+1)=l+liw+nsizi
1043       do 506 k=1,liw
1044          stk(l+k-1)=dble(istk(ilc+k-1))
1045  506  continue
1046       lsvi=l+liw
1047       if(meth.eq.0) then
1048          call svcma1(stk(lsvs),stk(lsvi))
1049       elseif(meth.lt.3) then
1050          call svcom1(stk(lsvs),stk(lsvi))
1051       else
1052          call svcar1(stk(lsvs),stk(lsvi))
1053       endif
1054       return
1055       end
1056