1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59
|
do i = 0, ndisc
if(have_to_init) then
h = h0
tin = t(i)
end if
have_to_init = .true.
tout = t(i+1)
if (type.eq.'IDE') then
call mebdfi(n,tin,h,y,dy,tout,t(i+1),mf,idid,lout,lwork,work,liwork, &
iwork,mbnd,maxder,itol,rtol,atol,rpar,ipar,pdervide,residide,ierr)
else if (type.eq.'ODE') then
call mebdfi(n,tin,h,y,dy,tout,t(i+1),mf,idid,lout,lwork,work,liwork, &
iwork,mbnd,maxder,itol,rtol,atol,rpar,ipar,pdervode,residode,ierr)
else if (type.eq.'DAE') then
call mebdfi(n,tin,h,y,dy,tout,t(i+1),mf,idid,lout,lwork,work,liwork, &
iwork,mbnd,maxder,itol,rtol,atol,rpar,ipar,pdervdae,residdae,ierr)
end if
if (ierr.eq.-1) then
h = h / 2d0
idid = -1
ierr = 0
have_to_init = .false.
cycle
elseif (idid.eq.1) then
if (printsolout) then
idid = 3
! The approximation at t0 + h0 is printed in the output file
write(90,formatout) tin, (y(indsol(it)),it=1,nindsol)
else
idid = 0
end if
have_to_init = .false.
cycle
elseif (idid.eq.3) then
! We print only the point in the considered interval [t(i) t(i+1)]
if ((printsolout).and.(tin.le.tout)) then
write(90,formatout) tin, (y(indsol(it)),it=1,nindsol)
end if
have_to_init = .false.
cycle
elseif (idid.ne.0) then
print *, 'MEBDFID: ERROR: ', 'MEBDFI returned IDID = ', idid
stop
endif
! The approximation in t(i+1) is printed in the output file
if (printsolout) then
write(90,formatout) tin, (y(indsol(it)),it=1,nindsol)
idid = 3
end if
end do |
Partager