utils.f90

Spectral Element Method for wave propagation and rupture dynamics.

! DSORT sorts array DX and optionally makes the same interchanges in
! array DY.  The array DX may be sorted in increasing order or
! decreasing order.  A slightly modified quicksort algorithm is used.

! Description of Parameters
! DX - array of values to be sorted   (usually abscissas)
! DY - array to be (optionally) carried along
! N  - number of values in array DX to be sorted
! KFLAG - control parameter
! =  2  means sort DX in increasing order and carry DY along.
! =  1  means sort DX in increasing order (ignoring DY)
! = -1  means sort DX in decreasing order (ignoring DY)
! = -2  means sort DX in decreasing order and carry DY along.

!***REFERENCES  R. C. Singleton, Algorithm 347, An efficient algorithm
! for sorting with minimal storage, Communications of
! the ACM, 12, 3 (1969), pp. 185-187.
!***ROUTINES CALLED  XERMSG
!***REVISION HISTORY  (YYMMDD)
! 761101  DATE WRITTEN
! 761118  Modified to use the Singleton quicksort algorithm.  (JAW)
! 890531  Changed all specific intrinsics to generic.  (WRB)
! 890831  Modified array declarations.  (WRB)
! 891009  Removed unreferenced statement labels.  (WRB)
! 891024  Changed category.  (WRB)
! 891024  REVISION DATE from Version 3.2
! 891214  Prologue converted to Version 4.0 format.  (BAB)
! 900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
! 901012  Declared all variables; changed X,Y to DX,DY; changed
! code to parallel SSORT. (M. McClain)
! 920501  Reformatted the REFERENCES section.  (DWL, WRB)
! 920519  Clarified error messages.  (DWL)
! 920801  Declarations section rebuilt and code restructured to use
! IF-THEN-ELSE-ENDIF.  (RWC, WRB)
!***END PROLOGUE  DSORT
! .. Array Arguments ..
    double precision :: DX(:), DY(:)
! .. Local Scalars ..
    double precision :: R, T, TT, TTY, TY
    INTEGER :: I, IJ, J, K, L, M, NN
! .. Local Arrays ..
    INTEGER :: IL(21), IU(21)
! .. External Subroutines ..
!    EXTERNAL XERMSG
! .. Intrinsic Functions ..
!    INTRINSIC ABS, INT
!***FIRST EXECUTABLE STATEMENT  DSORT
    NN = size(DX)

! Alter array DX to get decreasing order if needed
!
!    IF (KFLAG <= -1) THEN
!        DO 10 I=1,NN
!            DX(I) = -DX(I)
!        10 end do
!    ENDIF

! Sort DX and carry DY along

    M = 1
    I = 1
    J = NN
    R = 0.375D0

    110 IF (I == J) GO TO 150
    IF (R <= 0.5898437D0) THEN
        R = R+3.90625D-2
    ELSE
        R = R-0.21875D0
    ENDIF

    120 K = I

! Select a central element of the array and save it in location T

    IJ = I + INT((J-I)*R)
    T = DX(IJ)
    TY = DY(IJ)

! If first element of array is greater than T, interchange with T

    IF (DX(I) > T) THEN
        DX(IJ) = DX(I)
        DX(I) = T
        T = DX(IJ)
        DY(IJ) = DY(I)
        DY(I) = TY
        TY = DY(IJ)
    ENDIF
    L = J

! If last element of array is less than T, interchange with T

    IF (DX(J) < T) THEN
        DX(IJ) = DX(J)
        DX(J) = T
        T = DX(IJ)
        DY(IJ) = DY(J)
        DY(J) = TY
        TY = DY(IJ)
    
    ! If first element of array is greater than T, interchange with T
    
        IF (DX(I) > T) THEN
            DX(IJ) = DX(I)
            DX(I) = T
            T = DX(IJ)
            DY(IJ) = DY(I)
            DY(I) = TY
            TY = DY(IJ)
        ENDIF
    ENDIF

! Find an element in the second half of the array which is smaller
! than T

    130 L = L-1
    IF (DX(L) > T) GO TO 130

! Find an element in the first half of the array which is greater
! than T

    140 K = K+1
    IF (DX(K) < T) GO TO 140

! Interchange these elements

    IF (K <= L) THEN
        TT = DX(L)
        DX(L) = DX(K)
        DX(K) = TT
        TTY = DY(L)
        DY(L) = DY(K)
        DY(K) = TTY
        GO TO 130
    ENDIF

! Save upper and lower subscripts of the array yet to be sorted

    IF (L-I > J-K) THEN
        IL(M) = I
        IU(M) = L
        I = K
        M = M+1
    ELSE
        IL(M) = K
        IU(M) = J
        J = L
        M = M+1
    ENDIF
    GO TO 160

! Begin again on another portion of the unsorted array

    150 M = M-1
    IF (M == 0) GO TO 190
    I = IL(M)
    J = IU(M)

    160 IF (J-I >= 1) GO TO 120
    IF (I == 1) GO TO 110
    I = I-1

    170 I = I+1
    IF (I == J) GO TO 150
    T = DX(I+1)
    TY = DY(I+1)
    IF (DX(I) <= T) GO TO 170
    K = I

    180 DX(K+1) = DX(K)
    DY(K+1) = DY(K)
    K = K-1
    IF (T < DX(K)) GO TO 180
    DX(K+1) = T
    DY(K+1) = TY
    GO TO 170

! Clean up

!    190 IF (KFLAG <= -1) THEN
!        DO 200 I=1,NN
!            DX(I) = -DX(I)
!        200 end do
!    ENDIF
190    RETURN
    end SUBROUTINE DSORT

!-----------------------------------------------------------------------
!  HUNT is adapted from Numerical Recipes for double precision, F90
!
!  jlo = index in xx right before x
!      = 0 or length(xx) if x is out of xx bounds 
!
subroutine hunt(xx,x,jlo)
  integer, intent(inout) :: jlo
  double precision, intent(in) :: x,xx(:)
  integer :: inc,jhi,jm,n
  logical :: ascnd
  n = size(xx)
  ascnd=xx(n)>xx(1)
  if(jlo<=0.or.jlo>n)then
    jlo=0
    jhi=n+1
  else
   inc=1
   if(x>=xx(jlo).eqv.ascnd)then
     do
        jhi=jlo+inc
        if(jhi>n)then
          jhi=n+1
          exit
        else if(x>=xx(jhi).eqv.ascnd)then
          jlo=jhi
          inc=inc+inc
        else
          exit
        endif
     enddo
   else
    jhi=jlo
    do
        jlo=jhi-inc
        if(jlo<1)then
          jlo=0
          exit
        else if(x<xx(jlo).eqv.ascnd)then
          jhi=jlo
          inc=inc+inc
        else
          exit
        endif
    enddo
   endif
  endif
  do while(jhi-jlo/=1)
    jm=(jhi+jlo)/2
    if(x>xx(jm).eqv.ascnd)then
        jlo=jm
    else
        jhi=jm
    endif
  enddo
 
end subroutine hunt

!***********************************************************************
!  SPLINE INTERPOLATION
!  adapted from Numerical Recipes for double precision, Fortran 90

  SUBROUTINE spline(x,y,n,yp1,ypn,y2)

  INTEGER, intent(in) :: n
  double precision, intent(in) :: yp1,ypn,x(n),y(n)
  double precision, intent(out) :: y2(n)

  INTEGER i,k
  double precision :: p,qn,sig,un
  double precision :: u(n)

  if (yp1 > .99d30) then
    y2(1)=0d0
    u(1)=0d0
  else
    y2(1)=-0.5d0
    u(1)=(3d0/(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
  endif
  do i=2,n-1
    sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
    p=sig*y2(i-1)+2d0
    y2(i)=(sig-1d0)/p
    u(i)=(6d0*((y(i+1)-y(i)) &
         /(x(i+1)-x(i))-(y(i)-y(i-1)) &
         /(x(i)-x(i-1))) &
         /(x(i+1)-x(i-1)) &
         -sig*u(i-1))/p
  enddo
  if (ypn > .99e30) then
    qn=0d0
    un=0d0
  else
    qn=0.5d0
    un=(3d0/(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
  endif
    y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1d0)
  do k=n-1,1,-1
    y2(k)=y2(k)*y2(k+1)+u(k)
  enddo
  
  END SUBROUTINE spline

!---------------------------------------------------------------
!  adapted from Numerical Recipes for double precision, Fortran 90
!  and to speed up table lookup

  SUBROUTINE splint(xa,ya,y2a,n,x,y)

  INTEGER, intent(in) :: n
  double precision, intent(in) :: x,xa(n),y2a(n),ya(n)
  double precision, intent(out) :: y

  integer :: k,khi
  integer, save :: klo=1
  double precision :: a,b,h

! original : find bracketing indices bissection
!  klo=1
!  khi=n
!  do while (khi-klo > 1)
!    k=(khi+klo)/2
!    if(xa(k).gt.x)then
!      khi=k
!    else
!      klo=k
!    endif
!  enddo

! modified: 
  call hunt(xa,x,klo)
  khi=klo+1

  h=xa(khi)-xa(klo)
  a=(xa(khi)-x)/h
  b=(x-xa(klo))/h
  y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**2)/6d0
  
  END SUBROUTINE splint

end module utils