*-----------------------------------------------------------------------
*     Runge-Kutta-Gill algorithm routine.
*                                                 Oct. 5, 1990  S.Sakai
*-----------------------------------------------------------------------
*     Copyright (C) 2000-2016 GFD Dennou Club. All rights reserved.
*-----------------------------------------------------------------------
      SUBROUTINE ODRKG(N, FCN, T, DT, X, DX, XOUT, WORK)
      DIMENSION  X(N), DX(N), XOUT(N), WORK(N)

      DT2 = DT/2.
      C1 =  (1. - 1./SQRT(2.))*DT
      C2 =  -2. + 3./SQRT(2.)
      C3 =   2. -    SQRT(2.)
      C4 =  (1. + 1./SQRT(2.))*DT
      C5 =  (2. + 3./SQRT(2.))/3.
      C6 = -(2. +    SQRT(2.))/3.

*--------------------------- 1st step ----------------------------------

      DO 10 I=1,N
        XOUT(I) = X(I) + DT2*DX(I)
   10 CONTINUE

*--------------------------- 2nd step ----------------------------------

      TT = T+DT2
      CALL FCN(N, TT, XOUT, WORK)
      DO 20 I=1,N
        XOUT(I) = XOUT(I)  + C1*(WORK(I)-DX(I))
        DX(I)   = C2*DX(I) + C3*WORK(I)
   20 CONTINUE

*--------------------------- 3rd step ----------------------------------

      CALL FCN(N, TT, XOUT, WORK)
      DO 30 I=1,N
        XOUT(I) = XOUT(I)  + C4*(WORK(I)-DX(I))
        DX(I)   = C5*DX(I) + C6*WORK(I)
   30 CONTINUE

*--------------------------- 4th step ----------------------------------

      TT = T+DT
      CALL FCN(N, TT, XOUT, WORK)
      DO 40 I=1,N
        XOUT(I) = XOUT(I) + DT*(WORK(I)/6.+DX(I))
   40 CONTINUE

      RETURN
      END
