!************************************************************************
!*                                                                      *
!* Solve an ordinary system of first order differential equations using *
!* -------------------------------------------------------------------- *
!* automatic step size control                                          *
!* ----------------------------                                         *
!*                                                                      *
!* Programming language: ANSI C                                         *
!* Author:               Klaus Niederdrenk (FORTRAN)                    *
!* Adaptation:           Juergen Dietel, Computer Center, RWTH Aachen   *
!* Source:               existing C, Pascal, QuickBASIC and FORTRAN     *
!*                       codes                                          *
!* Date:                 6.2.1992, 10.2.1995                            *
!*                                                                      *
!*                       F90 Release By J-P Moreau, Paris.              *
!*                       (www.jpmoreau.fr)                              *
!* -------------------------------------------------------------------- *
!* REF.: "Numerical Algorithms with C, By Gisela Engeln-Muellges        *
!*        and Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11].           *
!************************************************************************

!debug only
Subroutine PrintVec(name, n, V)
character*(*) name
integer n,  i
real*8 V(0:n-1)
  print *, name
  write(*,10) (V(i), i=0,n-1)
  print *,' '
  return
10 format(10E12.6)
End


! Maximum norm of a difference vector ........
real*8 Function dist_max(vector1,vector2, n)
       real*8  vector1(0:n-1),  &
               vector2(0:n-1)
       integer n
!************************************************************************
!* Compute the maximum norm of the difference of two [0..n-1] vectors   *
!*                                                                      *
!* global name used:                                                    *
!* ================                                                     *
!*   None                                                               *
!************************************************************************
  real*8  abstand,   &     !reference value for computation of distance
          hilf             !distance of two vector elements
  integer i

  abstand = 0.d0
  do i=n-1, 0, -1
    hilf = DABS(vector1(i)-vector2(i))
    if (hilf > abstand)  abstand = hilf
  end do
  dist_max = abstand
  return
End


! Runge-Kutta embedding formulas of 2nd, 3rd degree ....................
Subroutine ruku23(x,y,bspn,n,h,y2,y3,yhilf,k1,k2,k3)
   real*8 x, y(0:n-1)
   integer bspn, n
   real*8 h, y2(0:n-1), y3(0:n-1)
! auxiliary vectors
   real*8 yhilf(0:n-1),k1(0:n-1),k2(0:n-1),k3(0:n-1)
!************************************************************************
!* Compute 2nd and 3rd order approximates y2, y3 at x + h starting with *
!* a solution y at x by using Runge-Kutta embedding formulas on the     *
!* first order system of n differential equations   y' = f(x,y) , as    *
!* supplied by  dgl().                                                  *
!*                                                                      *
!* Input parameters:                                                    *
!* =================                                                    *
!* x    x-value of left end point                                       *
!* y    y-values at x                                                   *
!* bspn # example                                                       *
!* n    number of differential equations                                *
!* h    step size                                                       *
!*                                                                      *
!* yhilf,k1,k2,k3: auxiliary vectors defined in awp.                    *
!*                                                                      *
!* Output parameters:                                                   *
!* ==================                                                   *
!* y2   2nd order approximation for y at x + h                          *
!* y3   3rd order approximation for y at x + h                          *
!*                                                                      *
!* External function:                                                   *
!* =================                                                    *
!* dgl  function that evaluates the right hand side of the system       *
!*      y' = f(x,y)                                                     *
!*      (See file t_dgls.f90).                                          * 
!************************************************************************
  integer i               !loop variable

  call dgl(bspn,n,x,y,k1)
  do i=0, n-1
    yhilf(i) = y(i) + h * k1(i)
  end do
  call dgl(bspn,n,x+h,yhilf,k2)
  do i=0, n-1
    yhilf(i) = y(i) + 0.25D0 * h * (k1(i) + k2(i))
  end do
  call dgl(bspn,n,x+0.5*h,yhilf,k3)
  do i=0, n-1
    y2(i) = y(i) + 0.5D0 * h * (k1(i) + k2(i))
    y3(i) = y(i) + h / 6.D0 * (k1(i) + k2(i) + 4.D0 * k3(i))
  end do
  return
End


! England's Einbettungs formulas of 4th and 5th degree ................
Subroutine engl45(x,y,bspn,n,h,y4,y5,yhilf,k1,k2,k3,k4,k5,k6)
  real*8    x                   !starting point of integration ........
  real*8    y(0:n-1)            !initial value at x ...................
  integer   bspn,       &       !# example
            n                   !number of differential equations .....
  real*8    h                   !step size ............................
  real*8   y4(0:n-1),   &       !4th order approximation for y at x + h
           y5(0:n-1)            !5th order approximation for y at x + h
  !auxiliary vectors
  real*8 yhilf(0:n-1),k1(0:n-1),k2(0:n-1),k3(0:n-1),k4(0:n-1),k5(0:n-1),k6(0:n-1)
!************************************************************************
!* Compute 4th and 5th order approximates y4, y5 at x + h starting with *
!* a solution y at x by using the England embedding formulas on the     *
!* first order system of n differential equations   y' = f(x,y) , as    *
!* supplied by  dgl().                                                  *
!*                                                                      *
!* Input parameters:                                                    *
!* =================                                                    *
!* x    initial x-value                                                 *
!* y    y-values at x, type pVEC                                        *
!* n    number of differential equations                                *
!* h    step size                                                       *
!*                                                                      *
!* yhilf, k1..K6: auxiliary vectors                                     *
!*                                                                      *
!* Output parameters:                                                   *
!* ==================                                                   *
!* y4   4th order approximation for y at x + h (pVEC)                   *
!* y5   5th order approximation for y at x + h (pVEC)                   *
!*                                                                      *
!* External function:                                                   *
!* =================                                                    *
!* dgl  function that evaluates the right hand side of the system       *
!*      y' = f(x,y)                                                     *
!*      (See file t_dgls.f90).                                          * 
!************************************************************************
  integer i               !loop variable
  call dgl(bspn,n,x,y,k1)
  do i=0, n-1
    yhilf(i) = y(i) + 0.5d0 * h * k1(i)
  end do
  call dgl(bspn,n,x+0.5d0*h,yhilf,k2)
  do i=0, n-1
    yhilf(i) = y(i) + (0.25D0 * h * (k1(i) + k2(i)))
  end do
  call dgl(bspn,n,x+0.5d0*h,yhilf,k3)
  do i=0, n-1
    yhilf(i) = y(i) + h * (-k2(i) + 2.D0 * k3(i))
  end do
  call dgl(bspn,n,x+h,yhilf,k4)
  do i=0, n-1
    yhilf(i) = y(i) + h / 27.D0 * (7.D0 * k1(i) + 10.D0 * k2(i) + k4(i))
  end do
  call dgl(bspn,n,x+(2.d0/3.d0)*h,yhilf,k5)
  do i=0, n-1
    yhilf(i) = y(i) + h / 625.d0 * (28.d0 * k1(i) - 125.d0 * k2(i) + &
               546.d0 * k3(i) + 54.d0 * k4(i) - 378.d0 * k5(i))
  end do
  call dgl(bspn,n,x+h/5.d0,yhilf,k6)
  do i=0, n-1
    y4(i) = y(i) + h / 6.d0 * (k1(i) + 4.d0 * k3(i) + k4(i))
    y5(i) = y(i) + h / 336.d0 * (14.d0 * k1(i) + 35.d0 * k4(i) +  &
            162.d0 * k5(i) + 125.d0 * k6(i))
  end do
  return
End


! embedding formulas of Prince-Dormand of 4./5. order .........................
Subroutine prdo45(x,y,bspn,n,h,y4,y5,steif1,steifanz,steif2, &
                  yhilf,k1,k2,k3,k4,k5,k6,k7, g6,g7)         
  real*8    x,           &          !starting point of integration .....
            y(0:n-1)                !initial value at x ................
  integer   bspn,        &          !# example .........................
            n                       !number of DEs .....................
  real*8    h,           &          !step size .........................
            y4(0:n-1),   &          !solution of 4th order at x+h ......
            y5(0:n-1)               !solution of 5th order at x+h ......
  !auxiliary flags
  integer steif1,steifanz,steif2
  !auxiliary vectors
  real*8 yhilf(0:n-1)
  real*8 k1(0:n-1),k2(0:n-1),k3(0:n-1),k4(0:n-1),k5(0:n-1),k6(0:n-1),k7(0:n-1)
  real*8 g6(0:n-1),g7(0:n-1)
!************************************************************************
!* Compute 4th and 5th order approximates y4, y5 at x + h starting with *
!* a solution y at x by using the Prince-Dormand embedding formulas on  *
!* the first order system of n differential equations y' = f(x,y) , as  *
!* supplied by  dgl().                                                  *
!* Simultaneously we perform two tests for stiffness whose results are  *
!* stored in steif1 and steif2.                                         *
!*                                                                      *
!* Input parameters:                                                    *
!* =================                                                    *
!* x    initial x-value                                                 *
!* y    y-values at x (pVEC)                                            *
!* n    number of differential equations                                *
!* h    step size                                                       *
!*                                                                      *
!* yhilf, k1..k7,g6,g7: auxiliary vectors.                              *
!*                                                                      *
!* Output parameters:                                                   *
!* ==================                                                   *
!* y4   4th order approximation for y at x + h                          *
!* y5   5th order approximation for y at x + h                          *
!*                                                                      *
!* External function:                                                   *
!* =================                                                    *
!* dgl  function that evaluates the right hand side of the system       *
!*      y' = f(x,y)                                                     *
!*      (See file t_dgls.f90).                                          * 
!***********************************************************************}
  integer i            !loop variable
  integer steifa       !Flag which is set if the second test for stiffness
                       !Shampine und Hiebert) is positive; otherwise the
                       !flag is erased.
  real*8 dist_max

  call dgl(bspn,n,x,y,k1)
  do i=0, n-1
    yhilf(i) = y(i) + 0.2D0 * h * k1(i)
  end do
  call dgl(bspn,n,x+0.2D0*h,yhilf,k2)
  do i=0, n-1
    yhilf(i) = y(i) + 0.075D0 * h * (k1(i) + 3.D0 * k2(i))
  end do
  call dgl(bspn,n,x+0.3D0*h,yhilf,k3)
  do i=0, n-1
    yhilf(i) = y(i) + h / 45.D0 * (44.D0 * k1(i) - 168.D0 * k2(i) + 160.D0 * k3(i))
  end do
  call dgl(bspn,n,x+0.8D0*h,yhilf,k4)
  do i=0, n-1
    yhilf(i) = y(i) + h / 6561.D0 * (19372.D0 * k1(i) - 76080.D0 * k2(i)  &
               + 64448.D0 * k3(i) - 1908.D0 * k4(i))
  end do
  call dgl(bspn,n,x+(8.D0/9.D0)*h,yhilf,k5)
  do i=0, n-1
    g6(i) = y(i) + h / 167904.D0 * (477901.D0 * k1(i) - 1806240.D0 * k2(i)  &
            + 1495424.D0 * k3(i) + 46746.D0 * k4(i) - 45927.D0 * k5(i))
  end do
  call dgl(bspn,n,x+h,g6,k6)
  do i=0, n-1
    g7(i) = y(i) + h / 142464.D0 * (12985.D0 * k1(i) + 64000.D0 * k3(i)  &
            + 92750.D0 * k4(i) - 45927.D0 * k5(i) + 18656.D0 * k6(i))
  end do
  call dgl(bspn,n,x+h,g7,k7)
  do i=0, n-1
    y5(i) = g7(i)
    y4(i) = y(i) + h / 21369600.D0 * (1921409.D0 * k1(i) + 9690880.D0 * k3(i)  &
            + 13122270.D0 * k4(i)- 5802111.D0 * k5(i) + 1902912.D0 * k6(i)     &
            + 534240.D0 * k7(i))
  end do

! Test for stiffness via dominant eigenvalue

  if (dist_max(k7, k6, n) > 3.3D0 * dist_max(g7, g6, n))  steif1 = 1

! one step in steffness test of Shampine & Hiebert

  do i=0, n-1
    g6(i) = h * (2.2D0 * k2(i) + 0.13D0 * k4(i) + 0.144D0 * k5(i))
    g7(i) = h * (2.134D0 * k1(i) + 0.24D0 * k3(i) + 0.1D0 * k6(i))
  end do

  if (dist_max(g6, g7, n) < dist_max(y4, y5, n))  then
    steifa = 1
  else  
    steifa = 0
  end if

  if (steifa > 0) then
    steifanz = steifanz + 1
    if (steifanz >= 3) steif2 = 1
  else
    steifanz = 0
  end if
  return

End


! Find the maximum norm of a REAL vector ............................
real*8 Function norm_max(vektor, n)    
       real*8  vektor(0:n-1)               ! vector .................
       integer n                           ! length of vector .......
! ************************************************************************
! *          Return the maximum norm of a [0..n-1] vector v              *
! *                                                                      *
! ************************************************************************
    real*8 norm,  &                                         ! local max.
    betrag                                    ! magnitude of a component
    integer i
  
    norm=0.D0
    do i=0, n-1
      betrag = dabs(vektor(i))
      if (betrag > norm)  norm = betrag
    end do   
    norm_max = norm
	return
End

! 1st order DESs with automatic step size control ...........................
Subroutine awp(x,xend,bspn,n,y,epsabs,epsrel,h,methode,fmax,aufrufe,fehler, &
               yhilf, k1, k2, k3, k4, k5, k6)
   real*8 x,             &        !initial/final x value ..............
          xend                    !desired end point ..................
  integer bspn,          &        !# example 
          n                       !number of DEs ......................
   real*8 y(0:n-1),      &        !initial/final y value ..............
          epsabs,        &        !absolute error bound ...............
          epsrel,        &        !relative error bound ...............
          h                       !initial/final step size ............
  integer methode,       &        !desired method (3, 6, 7) ...........
          fmax,          &        !maximal # of calls of  dgl() .......
          aufrufe,       &        !actual # of calls of  dgl() ........
          fehler                  !error code .........................
  !auxiliary vectors
  real*8 yhilf(0:9),k1(0:9),k2(0:9),k3(0:9),k4(0:9),k5(0:9),k6(0:9)
!************************************************************************
!* Compute the solution y of a system of first order ordinary           *
!* differential equations       y' = f(x,y)   at xend from the given    *
!* initial data (x0, y0).                                               *
!* We use automatic step size control internally so that the error of   *
!* y (absolutely or relatively) lies within the given error bounds      *
!* epsabs and epsrel.                                                   *
!*                                                                      *
!* Input parameters:                                                    *
!* =================                                                    *
!* x        initial value for x                                         *
!* y        initial values for y(0:n-1)                                 *
!* bspn     # example                                                   *
!* n        number of differential equations                            *
!* dgl      function that evaluates the right hand side of the system   *
!*          y' = f(x,y)  (see t_dgls) (Removed from list of parameters) *
!* xend     end of integration; xend > x0                               *
!* h        initial step size                                           *
!* epsabs   absolute error bound; >= 0; if = 0 we only check the        *
!*          relative error.                                             *
!* epsrel   relative error bound; >= 0; if = 0 we check only the        *
!*          absolute eror.                                              *
!* fmax     max number of evaluations of right hand side in dgl()       *
!* methode  chooses the method                                          *
!*          = 3: Runge-Kutta method of 2nd/3rd order                    *
!*          = 6: England formula of 4th/5th order                       *
!*          = 7: Formula of Prince-Dormand of 4th/5th order             *
!*                                                                      *
!* Output parameters:                                                   *
!* ==================                                                   *
!* x        final x-value for iteration. If fehler = 0  we usually have *
!*            x = xend.                                                 *
!* y        final y-values for the solution at x                        *
!* h        final step size used; leave for subsequent calls            *
!* aufrufe  actual number of calls of dgl()                             *
!*                                                                      *
!* Return value (fehler):                                               *
!* =====================                                                *
!* = 0: all ok                                                          *
!* = 1: both error bounds chosen too small for the given mach. constant *
!* = 2: xend <= x0                                                      *
!* = 3: h <= 0                                                          *
!* = 4: n <= 0                                                          *
!* = 5: more right hand side calls than allowed: aufrufe > fmax,        *
!*      x and h contain the current values when stop occured.           *
!* = 6: improper input for embedding formula                            *
!* = 7: lack of available memory (not used here)                        *
!* = 8: Computations completed, but the Prince Dormand formula stiff-   *
!*      ness test indicates possible stiffness.                         *
!* = 9: Computations completed, but both Prince Dormand formula stiff-  *
!*      ness tests indicate possible stiffness. Use method for stiff    *
!*      systems instead !                                               *
!* =10: aufrufe > fmax, see error code 5; AND the Prince Dormand formula*
!*      indicates stiffness; retry using a stiff DE solver !            *
!*                                                                      *
!************************************************************************
  real*8, parameter :: MACH_EPS=1.2d-16
  real*8, parameter :: MACH_2 = 100.d0 * MACH_EPS
  real*8   &    
  xend_h,  &       !|xend| - MACH_2, carrying same sign as xend
  ymax,    &       !Maximum norm of newest approximation of max
                   !order.
  hhilf,   &       !aux storage for the latest value of h
                   !produced by step size control. It is saved
                   !here in order to avoid to return a `h' that
                   !resulted from an arbitrary reduction at the
                   !end of the interval.
  diff,    &       !distance of the two approximations from the
                   !embedding formula.
  s                !indicates acceptance level for results from
                   !embeding formula.

  !approximate solution of low order
  real*8, pointer :: y_bad(:)
  !ditto of high order
  real*8, pointer ::y_good(:)
   
  real*8 mach_1    !machine constant dependent variable which
                   !avoids using too little steps near xend.
  integer i,  &    !Loop variable
  amEnde,     &    !flag that shows if the end of the interval
                   !can be reached with the actual step size.
  fertig           !flag indicating end of iterations.

  integer    &
  steif1,    &     !Flag, that is set in prdo45() if its
                   !stiffness test (dominant eigenvalue)
                   !indicates so. Otherwise no changes.
  steifanz,  &     !counter for number of successive successes
                   !of stiffness test of Shampine and Hiebert in
                   !prdo45().
  steif2,    &     !Flag, set in prdo45(), when the stiffness
                   !test of  Shampine and Hiebert wa successful
                   !three times in a row; otherwise no changes.
  ialloc           !allocation status (must be zero).
   
  !additional auxiliary vectors for prdo45 method (not used here)
  real*8 k7(0:n-1), g6(0:n-1), g7(0:n-1)

  real*8 dist_max, norm_max, Min, Max   !external functions
 
  !allocate auxiliary vectors
  allocate(y_bad(0:n-1), stat=ialloc)
  allocate(y_good(0:n-1), stat=ialloc)

  if (ialloc.ne.0) then
    print *,' awp: memory full !'!
    return
  end if

  !initialize some variables
  fehler   = 0                           
  mach_1   = MACH_EPS**0.75d0
  amEnde   = 0
  fertig   = 0
  steif1   = 0
  steif2   = 0
  steifanz = 0
  aufrufe  = 1
  ymax     = norm_max(y, n)

  if (xend >= 0.d0) then 
    xend_h = xend * (1.d0 - MACH_2)
  else 
    xend_h = xend * (1.d0 + MACH_2)
  end if

! ----------------------- check inputs ----------------------
  if (epsabs <= MACH_2 * ymax.and.epsrel <= MACH_2) then
    fehler=1
    return
  end if
  if (xend_h < x) then
    fehler= 2
    return
  end if
  if (h < MACH_2 * DABS(x)) then
    fehler=3
    return
  end if
  if (n <= 0) then
    fehler=4
    return
  end if
  if (methode.ne.3.and.methode.ne.6.and.methode.ne.7) then
    fehler=6
    return
  end if
 
! **********************************************************************
! *                                                                    *
! *                       I t e r a t i o n s                          *
! *                                                                    *
! **********************************************************************
  if (x + h > xend_h) then            
                                    !almost at end point ?
    hhilf  = h                      !A shortened step might be
    h      = xend - x               !enough.
    amEnde = 1
  end if

  do while(fertig.eq.0)             !solve DE system by integrating from
                                    !x0 to xend by suitable steps.
    Select Case(methode)
      Case(3) 
	    call ruku23(x, y, bspn,n, h, y_bad, y_good, yhilf,k1,k2,k3)
      Case(6)
	    call engl45(x, y, bspn, n, h, y_bad, y_good, yhilf,k1,k2,k3,k4,k5,k6)
      Case(7)
	    call prdo45(x, y, bspn, n, h, y_bad, y_good, steif1,steifanz,steif2, &
		            yhilf,k1,k2,k3,k4,k5,k6,k7, g6,g7)
    End Select

    aufrufe = aufrufe + methode

	diff = dist_max(y_bad, y_good, n)

    if (diff < MACH_2) then         !compute s
      s = 2.d0
    else
      ymax = norm_max(y_good, n)
      s    = DSQRT(h * (epsabs + epsrel * ymax) / diff)
      if (methode.ne.3)   s = DSQRT(s)
    end if

    if (s > 1.d0) then              !integration acceptable ?
      do i=0, n-1                   !accept highest order solution
        y(i) = y_good(i)            !move x
      end do
      
	  x = x + h

      if (amEnde.ne.0) then         !at end of interval ?
        fertig = 1                  !stop iteration
        if (methode.eq.7) then
          if (steif1>0.or.steif2>0)  fehler = 8
          if (steif1>0.and.steif2>0) fehler = 9
        end if
      else if (aufrufe > fmax) then !too many calls of   dgl() ?
        hhilf  = h                  !save actual step size
        fehler = 5                  !report error and stop
        fertig = 1
        if (methode.eq.7.and.(steif1>0.or.steif2>0))  fehler = 10
      else                          !Integration was successful
                                    !not at the interval end ?
        h = h * Min(2.d0, 0.98d0*s) !increase step size for next
                                    !step properly, at most by
                                    !factor two. Value `0.98*s' is
                                    !used in order to avoid that
                                    !the theoretical value s is
                                    !exceeded by accidental
                                    !rounding errors.
        if (x + h > xend_h) then    !nearly reached xend ?
          hhilf  = h                !=> One further step with
          h      = xend - x         !reduced step size might be
          amEnde = 1                !enough.
          if (h < mach_1 * DABS(xend)) then  !very close to xend ?
            fertig = 1                       !finish iteration.
          end if
        end if
      end if
    else                            !step unsuccessful ?
                                    !before repeating this step
      h = h * Max(0.5d0, 0.98d0*s)  !reduce step size properly, at
                                    !most by factor 1/2 (for factor
      amEnde = 0                    !0.98: see above).
    end if

  end do
 
  h = hhilf           !return the latest step size computed by step
                      !size control and  error code to the caller.

  return
End

! -------------------------- END  uawp.f90 ----------------------------