subroutine rka( x, nX, t, tau, err, derivsRK, param ) integer*4 MAXnX, MAXnparam parameter( MAXnX = 50, MAXnparam = 1000 ) integer*4 nX real*8 x(MAXnX), t, tau, err, param(MAXnparam) external derivsRK ! Adaptive Runge-Kutta routine ! Inputs ! x Current value of the dependent variable ! nX Number of elements in dependent variable x ! t Independent variable (usually time) ! tau Step size (usually time step) ! err Desired fractional local truncation error ! derivsRK Right hand side of the ODE; derivsRK is the ! name of the function which returns dx/dt ! Calling format derivsRK(x,t,param). ! param Extra parameters passed to derivsRK ! Outputs ! x New value of the dependent variable ! t New value of the independent variable ! tau Suggested step size for next call to rka integer*4 i, iTry, maxTry real*8 tSave, safe1, safe2, xSmall(MAXnX), xBig(MAXnX) real*8 errorRatio, eps, scale, xDiff, ratio, tau_old, half_tau !* Set initial variables tSave = t ! Save initial value safe1 = 0.9 safe2 = 4.0 ! Safety factors !* Loop over maximum number of attempts to satisfy error bound maxTry = 100 do iTry=1,maxTry !* Take the two small time steps half_tau = 0.5 * tau do i=1,nX xSmall(i) = x(i) enddo call rk4(xSmall,nX,tSave,half_tau,derivsRK,param) t = tSave + half_tau call rk4(xSmall,nX,t,half_tau,derivsRK,param) !* Take the single big time step t = tSave + tau do i=1,nX xBig(i) = x(i) enddo call rk4(xBig,nX,tSave,tau,derivsRK,param) !* Compute the estimated truncation error errorRatio = 0.0 eps = 1.0e-16 do i=1,nX scale = err * (abs(xSmall(i)) + abs(xBig(i)))/2.0 xDiff = xSmall(i) - xBig(i) ratio = abs(xDiff)/(scale + eps) if( ratio .gt. errorRatio ) then errorRatio = ratio endif enddo !* Estimate new tau value (including safety factors) tau_old = tau tau = safe1*tau_old*errorRatio**(-0.20) if( tau .lt. tau_old/safe2 ) then tau = tau_old/safe2 else if( tau .gt. safe2*tau_old ) then tau = safe2*tau_old endif !* If error is acceptable, return computed values if (errorRatio .lt. 1) then do i=1,nX x(i) = xSmall(i) enddo return endif enddo !* Issue error message if error bound never satisfied write(*,*) 'ERROR: Adaptive Runge-Kutta routine failed' return end