subroutine rk40(y,x,h,yout) c c======================================================================= c === c Given Y and its derivatives DYDX known at X, use the fourth-order === c Runge-Kutta method to advance the solution over an interval H and === c return the incremented variables as YOUT. Use DERIVS0 to compute === c the derivatives DYDX at X. === c === c======================================================================= c c----------------------------------------------------------------------- c Define global data. c----------------------------------------------------------------------- c #include #ifdef cyclic # include #endif #include c c----------------------------------------------------------------------- c Define local data. c----------------------------------------------------------------------- c integer i,index,n,nmax #ifdef cyclic * ,ncyc #endif parameter (nmax=3) FLOAT * h,hh,x,xh FLOAT * dydx(nmax),dyt(nmax),dym(nmax),y(8),yout(8),yt(nmax) c c----------------------------------------------------------------------- c Time origin is controlled from calling subroutine and DERIVS0 c----------------------------------------------------------------------- c index=0 n=3 call derivs0(x,y,dydx,index) hh=h*p5 xh=x+hh do 11 i=1,n yt(i)=y(i)+hh*dydx(i) 11 continue index=1 call derivs0(xh,yt,dyt,index) do 12 i=1,n yt(i)=y(i)+hh*dyt(i) 12 continue call derivs0(xh,yt,dym,index) do 13 i=1,n yt(i)=y(i)+h*dym(i) dym(i)=dyt(i)+dym(i) 13 continue index=2 call derivs0(x+h,yt,dyt,index) do 14 i=1,n yout(i+4)=(dydx(i)+dyt(i)+c2*dym(i))/c6 yout(i)=y(i)+h*yout(i+4) 14 continue #ifdef cyclic ncyc=int(yout(1)-1)/imtm2 if(yout(1).lt.c1) ncyc=ncyc-1 yout(1)=yout(1)-FLoaT(ncyc*imtm2) #endif return end