text:
!
!*******************************************************************************
!
subroutine gauss_quad(nt,kind,alpha,beta,x,w)
!
!*******************************************************************************
!
!  Discussion:
!
!  FORTRAN 90 subroutine to calculate abscissas and weights for Gauss-quadruature
!  for a wide range of orthogonal polynomials ortogonalized over a given range
!  [a,b] and weight function w(x).
!
!  Inputs:
!         nt    = Number of Gauss-quadruature points, [integer(kind = 4)]
!         kind  = Type of Quadrature rule,  [integer(kind = 4)]
!
!               = 1,   Legendre,             [-1,1],     w(x) = 1.0
!               = 2,   Associated Laguerre,  [0,inf],    w(x) = exp(-x)*x^alpha
!               = 3,   Hermite,              [-inf,inf], w(x) = exp(-x^2)
!               = 4,   Jacobi,               [-1,1],     w(x) = (1-x)^alpha*(1+x)^beta
!               = 5,   Chebyshev (Type I),   [-1,1],     w(x) = (1-x^2)^(-0.5)
!               = 6,   Chebyshev (type II),  [-1,1],     w(x) = (1-x^2)^(0.5)
!               = 7,   Gegenbauer,           [-1,1],     w(x) = (1-x^2)^(alpha-0.5)
!
!         alpha = Value of alpha, if needed, [real(kind = 8)]
!         beta  = Value of beta, if needed, [real(kind = 8)]
!
!   Outputs:
!         x(1:nt) = abscissas for Gauss-quadrature rule, [real(kind=8)]
!         w(1:nt) = weights for Gauss-quadrature rule, [real(kind=8)]
!
!   Note that Chebyshev Type II is a special case of Gegenbauer with alpha = 1
!
!   The algorithm is based on the three-term recurrence relation for each
!   polynomial. The terms are found for the normalized polynonmials, i.e.,
!
!   F(n,alpha,beta,x) = p(n,alpha,beta,x)/sqrt(N(n,alpha,beta)),
!
!   where n is the order of the polynomial, and p(n,alpha,beta) denotes
!   the standard textbook definition with p(0,alpha,beta) = 1.0, and
!   N(n,alpha,beta) is the norm integral defined as
!
!   N(n,alpha,beta) = int[a,b] dx w(x)*[p(n,alpha,beta,x)]^2.
!
!   The three-term recurrence relations can be written as (with fixed
!   alpha and beta)
!
!   a(n-1)*F(n-1,x) + b(n)*F(n,x) + a(n)*F(n+1,x) = x*F(n,x)
!
!   Using the notation
!   <g,p> = int[a,b] dx w(x)*g(x)*F(n,x).
!
!   b(n) = <x*F(n,alpha,beta),F(n,alpha,beta)>
!
!   The coefficients a(n) = k(n)/k(n+1)
!
!   where k(n) is the leading coefficient for F(n), that is
!
!   F(n) = k(n)*x^n + k(n-1)*x^(n-1) + ...
!
!   The abscicca x(k) is the kth eigenvalue of the tridigonal matrix
!
!   b(0)   a(0)      0      0      0                  0       0
!   a(0)   b(1)   a(1)      0      0                  0       0
!      0   a(1)   b(2)   a(2)      0                  0       0
!      0      0   a(2)   b(3)   a(3)                  0       0
!
!                                                b(n-1)  a(n-1)
!      0      0      0      0      0             a(n-1)    b(n)
!
!   The weights are obtained from the 1st component of the kth eigenvector,
!   vec(1,k), via:
!
!   w(k) = N(0,alpha,beta)*vec(1,k)^2/(sum_k vec(1,k)^2)
!
!   The eigenvalues and eigenvectors are obtained from a modified version of
!   the EISPACK routine tql2, named tri_eig
!
!   Dependencies:
!
!     Modules:
!
!        None
!
!     Subroutines:
!
!        jacobi_norm
!        tri_eig
!
!     External functions:
!
!       real(kind=8) :: poly_a
!       real(kind=8) :: poly_b
!       real(kind=8) :: poly_norm
!
!     MPI routines:
!
!        None
!
!  Licensing:
!
!    SPDX-License-Identifier: MIT
!
!  Date:
!
!    11 May 2021
!
!  Author:
!
!      Erich Ormand, LLNL
!
!*******************************************************************************
!
   implicit none
   integer(kind=4), intent(in) :: nt, kind
   real(kind=8), intent(in) :: alpha, beta
   real(kind=8), intent(out) :: x(nt), w(nt)

   integer(kind=4) :: i, n
   real(kind=8) :: pi

   real(kind=8) :: xn, two_xn

   real(kind=8), allocatable :: d(:), e(:), z(:,:)
   real(kind=8) :: norm, int_norm
   integer(kind=4) :: ierror
!
!-----    External functions
!
   real(kind=8) :: poly_a
   real(kind=8) :: poly_b
   real(kind=8) :: poly_norm
!
!*******************************************************************************
!
   pi = 2.0d0*asin(1.0d0)
!
!----   Allocate and initalize arrays for the tri-diagonal matrix
!
   if(.not.allocated(d))allocate(d(nt))
   if(.not.allocated(e))allocate(e(nt))
   if(.not.allocated(z))allocate(z(nt,nt))
   d(1:nt) = 0.0d0
   e(1:nt) = 0.0d0
   z(1:nt,1:nt) = 0.0d0
!
!----   Set up tridiagonal matrix to be diagonalized for each integration type
!
   do i = 1, nt
      n = i - 1                           !   integer order of the polynomial, p(n,alpha,beta)
      xn = real(n,kind=8)
      two_xn = 2.0d0*xn
      d(i) = poly_b(n, kind, alpha, beta)
      if(i < nt)e(i+1)= poly_a(n, kind, alpha, beta)
      z(i,i) = 1.0d0
   end do
!
!   Find eigenvalues and eigenvectors of the tri-diagonal matrix
!
   call tri_eig(nt, nt, d, e, z, ierror)
!
!----    Find norm to set the weights
!
   norm = 0.0d0
   do i = 1, nt
      norm = norm + z(1,i)**2
   end do
!
!----    Set norm for the weights - integral with w(x)*[p(0,alpha,beta)]^2*dx
!
   int_norm = poly_norm(0,kind,alpha,beta)
!
!----   Put abscissas and weights into output arrays
!
   do i = 1, nt
      x(i) = d(i)
      w(i) = int_norm*z(1,i)**2/norm
   end do
!
!----   Delete temporary data
!
   if(allocated(d))deallocate(d)
   if(allocated(e))deallocate(e)
   if(allocated(z))deallocate(z)
!
!----   All done
!
   return
!
end subroutine gauss_quad
!
!*******************************************************************************
!
real(kind=8) function poly_b(n, kind, alpha, beta)
!
!*******************************************************************************
!
!  Discussion:
!
!  FORTRAN 90 function to calculate recurssion coefficient b(n), where
!   the three-term recurrence relations are written as (with fixed
!   alpha and beta) for the normalized orthogonal polynomials as:
!
!   a(n-1)*F(n-1,x) + b(n)*F(n,x) + a(n)*F(n+1,x) = x*F(n,x)
!
!  Inputs:
!         n     = order of the polynomial, [integer(kind = 4)]
!         kind  = Type of polynomial,  [integer(kind = 4)]
!
!               = 1,   Legendre,             [-1,1],     w(x) = 1.0
!               = 2,   Associated Laguerre,  [0,inf],    w(x) = exp(-x)*x^alpha
!               = 3,   Hermite,              [-inf,inf], w(x) = exp(-x^2)
!               = 4,   Jacobi,               [-1,1],     w(x) = (1-x)^alpha*(1+x)^beta
!               = 5,   Chebyshev (Type I),   [-1,1],     w(x) = (1-x^2)^(-0.5)
!               = 6,   Chebyshev (type II),  [-1,1],     w(x) = (1-x^2)^(0.5)
!               = 7,   Gegenbauer,           [-1,1],     w(x) = (1-x^2)^(alpha-0.5)
!
!         alpha = Value of alpha, if needed, [real(kind = 8)]
!         beta  = Value of beta, if needed, [real(kind = 8)]
!
!   Dependencies:
!
!     Modules:
!
!        None
!
!     Subroutines:
!
!        None
!
!     External functions:
!
!        None
!
!     MPI routines:
!
!        None
!
!  Licensing:
!
!    SPDX-License-Identifier: MIT
!
!  Date:
!
!    11 May 2021
!
!  Author:
!
!      Erich Ormand, LLNL
!
!*******************************************************************************
!
    implicit none
    integer(kind=4), intent(in) :: n, kind
    real(kind=8), intent(in) :: alpha, beta
!--------------------------------------------------------------------
    real(kind=8) :: xn, two_xn, xLn
!*************************************************************************
    poly_b = 0.0d0
    xn = real(n,kind=8)
    two_xn = 2.0d0*xn
    if(kind == 1)then                       !   Legendre
       poly_b = 0.0d0
    elseif(kind == 2)then                   !   Associated Laguerre
       poly_b = two_xn + alpha + 1.0d0
    elseif(kind == 3)then                   !   Hermite
       poly_b = 0.0d0
    elseif(kind == 4)then                   !   Jacobi
       poly_b = 0.0d0
       if(n == 0 .and. abs(beta**2 - alpha**2) < 1.0d-10)return
       xLn = two_xn + alpha + beta
       poly_b = (beta**2 - alpha**2)/(xLn*(xLn+2.0d0))
    elseif(kind == 5)then                   !   Chebyshev (Type I)
       poly_b = 0.0d0
    elseif(kind == 6)then                   !   Chebyshev (Type II)
       poly_b = 0.0d0
    elseif(kind == 7)then                   !   Gegenbauer
       if(alpha < 1.0d0)stop 'Error in Gengenbauer polynomilas, alpha < 1'
       poly_b = 0.0d0
    end if
    return
end function poly_b
!
!*******************************************************************************
!
real(kind=8) function poly_a(n, kind, alpha, beta)
!
!*******************************************************************************
!
!  Discussion:
!
!  FORTRAN 90 function to calculate recurssion coefficient a(n), where
!   the three-term recurrence relations are written as (with fixed
!   alpha and beta) for the normalized orthogonal polynomials as:
!
!   a(n-1)*F(n-1,x) + b(n)*F(n,x) + a(n)*F(n+1,x) = x*F(n,x)
!
!  Inputs:
!         n     = order of the polynomial, [integer(kind = 4)]
!         kind  = Type of polynomial,  [integer(kind = 4)]
!
!               = 1,   Legendre,             [-1,1],     w(x) = 1.0
!               = 2,   Associated Laguerre,  [0,inf],    w(x) = exp(-x)*x^alpha
!               = 3,   Hermite,              [-inf,inf], w(x) = exp(-x^2)
!               = 4,   Jacobi,               [-1,1],     w(x) = (1-x)^alpha*(1+x)^beta
!               = 5,   Chebyshev (Type I),   [-1,1],     w(x) = (1-x^2)^(-0.5)
!               = 6,   Chebyshev (type II),  [-1,1],     w(x) = (1-x^2)^(0.5)
!               = 7,   Gegenbauer,           [-1,1],     w(x) = (1-x^2)^(alpha-0.5)
!
!         alpha = Value of alpha, if needed, [real(kind = 8)]
!         beta  = Value of beta, if needed, [real(kind = 8)]
!
!   Dependencies:
!
!     Modules:
!
!        None
!
!     Subroutines:
!
!        jacobi_norm
!        tri_eig
!
!     External functions:
!
!        real(kind=8) :: poly_norm
!
!     MPI routines:
!
!        None
!
!  Licensing:
!
!    SPDX-License-Identifier: MIT
!
!  Date:
!
!    11 May 2021
!
!  Author:
!
!      Erich Ormand, LLNL
!
!*******************************************************************************
!
    implicit none
    integer(kind=4), intent(in) :: n, kind
    real(kind=8), intent(in) :: alpha, beta
!--------------------------------------------------------------------
    real(kind=8) :: xn, two_xn, XLn, temp
!---------   External functions
    real(kind=8) :: poly_norm
!*************************************************************************
!
    poly_a = 0.0d0
    xn = real(n,kind=8)
    two_xn = 2.0d0*xn
    if(kind == 1)then                       !   Legendre
       poly_a = (xn + 1.0d0)/sqrt((two_xn + 3.0d0)*(two_xn + 1.0d0))
    elseif(kind == 2)then                   !   Associated Laguerre
       poly_a = -sqrt((xn+1.0d0)*(xn+alpha+1.0d0))
    elseif(kind == 3)then                   !   Hermite
       poly_a = sqrt((xn + 1.0d0)/2.0d0)
    elseif(kind == 4)then                   !   Jacobi
       xLn = two_xn + alpha + beta
       temp = log_gamma(alpha + beta + two_xn + 1.0d0) + log_gamma(alpha + beta + xn + 2.0d0) -   &
              log_gamma(alpha + beta + xn + 1.0d0) - log_gamma(alpha + beta + two_xn + 3.0d0)
       temp = exp(temp)*2.0d0*(xn + 1.0d0)*sqrt(poly_norm(n+1,kind,alpha,beta)/poly_norm(n,kind,alpha,beta))
       poly_a = temp
    elseif(kind == 5)then                   !   Chebyshev (Type I)
       if(n == 0)then
          poly_a = 1.0d0/sqrt(2.0d0)
       else
          poly_a = 0.5d0
       end if
    elseif(kind == 6)then                   !   Chebyshev (Type II)
       poly_a = 0.5d0
    elseif(kind == 7)then                   !   Gegenbauer
       if(alpha < 1.0d0)stop 'Error in Gengenbauer polynomilas, alpha < 1'
       poly_a = 0.5d0*sqrt((xn+1.0d0)*(xn+2.0d0*alpha)/((xn+alpha)*(xn+alpha+1.0d0)))
    end if
    return
end function poly_a
!
!*******************************************************************************
!
real(kind=8) function poly_norm(n, kind, alpha, beta)
!
!*******************************************************************************
!
!  Discussion:
!
!  FORTRAN 90 function to calculate the normalization of orthogonal Polynomials
!
!   int[-1,1] w(x) [P(n,alpha,beta)]^2 dx
!  Inputs:
!         n     = order of the polynomial, [integer(kind = 4)]
!         kind  = Type of polynomial,  [integer(kind = 4)]
!
!               = 1,   Legendre,             [-1,1],     w(x) = 1.0
!               = 2,   Associated Laguerre,  [0,inf],    w(x) = exp(-x)*x^alpha
!               = 3,   Hermite,              [-inf,inf], w(x) = exp(-x^2)
!               = 4,   Jacobi,               [-1,1],     w(x) = (1-x)^alpha*(1+x)^beta
!               = 5,   Chebyshev (Type I),   [-1,1],     w(x) = (1-x^2)^(-0.5)
!               = 6,   Chebyshev (type II),  [-1,1],     w(x) = (1-x^2)^(0.5)
!               = 7,   Gegenbauer,           [-1,1],     w(x) = (1-x^2)^(alpha-0.5)
!
!         alpha = Value of alpha, if needed, [real(kind = 8)]
!         beta  = Value of beta, if needed, [real(kind = 8)]
!
!   Dependencies:
!
!     Modules:
!
!        None
!
!     Subroutines:
!
!        None
!
!     External functions:
!
!        None
!
!     MPI routines:
!
!        None
!
!  Licensing:
!
!    SPDX-License-Identifier: MIT
!
!  Date:
!
!    11 May 2021
!
!  Author:
!
!      Erich Ormand, LLNL
!
!*******************************************************************************
!
    implicit none
    integer(kind=4), intent(in) :: n, kind
    real(kind=8), intent(in) :: alpha, beta
!--------------------------------------------------------------------
    real(kind=8) :: factor, xn, two_xn, pi
!*************************************************************************
!
    pi = 2.0d0*asin(1.0d0)
    poly_norm = 0.0d0
    xn = real(n,kind=8)
    two_xn = 2.0d0*xn
    if(kind == 1)then                       !   Legendre
       poly_norm = 2.0d0/(two_xn + 1.0d0)
    elseif(kind == 2)then                   !   Associated Laguerre
       poly_norm = exp(log_gamma(xn + alpha + 1.0d0) - log_gamma(xn + 1.0d0))
    elseif(kind == 3)then                   !   Hermite
       poly_norm = 2.0d0**xn*sqrt(pi)*exp(log_gamma(xn + 1.0d0))
    elseif(kind == 4)then                   !   Jacobi
       factor = log_gamma(xn + alpha + 1.0d0) + log_gamma(xn + beta + 1.0d0) -  &
                log_gamma(xn + 1.0d0) - log_gamma(xn + alpha + beta + 1.0d0)
       factor = factor + (alpha + beta + 1.0d0)*log(2.0d0) -                    &
                log(2.0d0*xn + alpha + beta + 1.0d0)
       poly_norm = exp(factor)
    elseif(kind == 5)then                   !   Chebyshev (Type I)
       if(n == 0)then
          poly_norm = pi
       else
          poly_norm = 0.5d0*pi
       end if
    elseif(kind == 6)then                   !   Chebyshev (Type II)
       poly_norm = 0.5d0*pi
    elseif(kind == 7)then                   !   Gegenbauer
       factor = log_gamma(xn + 2.0d0*alpha) - log_gamma(xn + 1.0d0) - 2.0d0*log_gamma(alpha)
       poly_norm = pi*2.0d0**(1.0d0 - 2.0d0*alpha)*exp(factor)/(xn + alpha)
    end if
    return
end function poly_norm
!
!*******************************************************************************
!
real(kind=8) function poly(n, kind, alpha, beta, x)
!
!*******************************************************************************
!
!  Discussion:
!
!  FORTRAN 90 function to calculate values at a given value of x for the
!   orthogonal polynomial p(n,alpha,beta,x).
!   Uses the three-term recurrence relation (with fixed
!   alpha and beta) for the normalized orthogonal polynomials F(n,alpha,beta,x) as:
!
!   a(n-1)*F(n-1,x) + b(n)*F(n,x) + a(n)*F(n+1,x) = x*F(n,x)
!
!   or
!
!   F(n) = ((x - b(n-1)*F(n-1)) - a(n-2)*F(n-2))/a(n-1)
!
!   with
!
!   F(n,alpha,beta,x) = p(n,alpha,beta,x)/sqrt(N(n,alpha,beta)),
!
!  Inputs:
!         n     = order of the polynomial, [integer(kind = 4)]
!         kind  = Type of polynomial,  [integer(kind = 4)]
!
!               = 1,   Legendre,             [-1,1],     w(x) = 1.0
!               = 2,   Associated Laguerre,  [0,inf],    w(x) = exp(-x)*x^alpha
!               = 3,   Hermite,              [-inf,inf], w(x) = exp(-x^2)
!               = 4,   Jacobi,               [-1,1],     w(x) = (1-x)^alpha*(1+x)^beta
!               = 5,   Chebyshev (Type I),   [-1,1],     w(x) = (1-x^2)^(-0.5)
!               = 6,   Chebyshev (type II),  [-1,1],     w(x) = (1-x^2)^(0.5)
!               = 7,   Gegenbauer,           [-1,1],     w(x) = (1-x^2)^(alpha-0.5)
!
!         alpha = Value of alpha, if needed, [real(kind = 8)]
!         beta  = Value of beta, if needed, [real(kind = 8)]
!
!   Dependencies:
!
!     Modules:
!
!        None
!
!     Subroutines:
!
!        None
!
!     External functions:
!
!        real(kind=8) :: poly_a
!        real(kind=8) :: poly_b
!        real(kind=8) :: poly_norm
!
!     MPI routines:
!
!        None
!
!  Licensing:
!
!    SPDX-License-Identifier: MIT
!
!  Date:
!
!    11 May 2021
!
!  Author:
!
!      Erich Ormand, LLNL
!
!*******************************************************************************
!
   implicit none
   integer(kind=4), intent(in) :: n, kind
   real(kind=8), intent(in) :: alpha, beta
   real(kind=8), intent(in) :: x
!----------------------------------------------------------------
   integer(kind=4) :: nn
   real(kind=8) :: p0, p1, p2
!
!-----    External functions
!
   real(kind=8) :: poly_a
   real(kind=8) :: poly_b
   real(kind=8) :: poly_norm
!
!*******************************************************************************
!
   p0 = 1.0d0
   p1 = 0.0d0
!---   If n=1, trivial, and exit
   poly = p0
   if(n == 0)return
   if(kind == 1)then
      p1 = x
   elseif(kind == 2)then
      p1 = (-x + alpha + 1.0d0)
   elseif(kind == 3)then
      p1 = 2.0d0*x
   elseif(kind == 4)then
      p1 = (alpha + 1.0d0) + (alpha + beta + 2.0d0)*0.5d0*(x-1.0d0)
   elseif(kind == 5)then
      p1 = x
   elseif(kind == 6)then
      p1 = 2.0d0*x
   elseif(kind == 7)then
      p1 = 2.0d0*alpha*x
   end if
!---   If n=1, trivial, and exit
   poly = p1
   if(n == 1)return
!---   First and second normalized polynomials
   p0 = p0/sqrt(poly_norm(0,kind,alpha,beta))
   p1 = p1/sqrt(poly_norm(1,kind,alpha,beta))
   p2 = 0.0d0
!---   Build up to n via recurrence relation
   do nn = 2, n
      p2 = (x - poly_b(nn-1,kind,alpha,beta))*p1 -            &
           poly_a(nn-2,kind,alpha,beta)*p0
      p2 = p2/poly_a(nn-1,kind,alpha,beta)
      p0 = p1
      p1 = p2
   end do
!---   Undo normalization and return
   poly = p2*sqrt(poly_norm(n,kind,alpha,beta))
   return
end function poly

max_stars_repo_name: LLNL/Yet-Another-Hauser-Feshbach-Code

max_issues_repo_name: LLNL/Yet-Another-Hauser-Feshbach-Code

max_forks_repo_name: LLNL/Yet-Another-Hauser-Feshbach-Code


text:
The Blue Moon Literary & Art Review publishes fiction and poetry of all genres, from serious fiction to murder mystery. We are especially interested in short stories and excerpts from novels in progress, and we are dedicated to high quality, interesting and unusual works. We hope the Review will be entertaining and thought provoking; however, we have no political agenda, only a strong desire to promote and publish fine writing. In addition to works of fiction, we also invite for consideration reviews of art, books, and plays, as well as literary essays of any kind.

While we are a new publication, we hope to grow over the years as we promote and encourage new writers from all over the world. In these early years of production, we request that writers submit their work with a stamped, selfaddressed envelope (and a small Reader’s Fee of $10.00 if you would like a written critique, made payable to Blue Moon Literary & Art Review, in order to defray the expenses associated with reading and considering each submission).

The Fall 2009 Issue 4 is available in Davis at the Avid Reader, Newsbeat, John Natsoulas Center for the Arts, and Konditorei Austrian Pastry Cafe.  In Sacramento, copies may be found at Newsbeat located at 1050 20th Street. The next volume will be published in Spring 2010.

Our goal is to lead the literary life and to enable others to lead it with us. Good writing, like fine wine, is best shared and appreciated as a communal experience.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


submodule (complex_transform_routines) complex_initialization_routines

contains

    module subroutine cfft1i(n, wsave, lensav, ierror)
        !
        ! cfft1i: initialization for cfft1b and cfft1f.
        !
        !  Purpose:
        !
        !  cfft1i initializes array wsave for use in its companion routines
        !  cfft1b and cfft1f. Routine cfft1i must be called before the first
        !  call to cfft1b or cfft1f, and after whenever the value of integer
        !  n changes.
        !
        !  Parameters:
        !
        !  input,
        !  n, the length of the sequence to be
        !  transformed.  the transform is most efficient when n is a product
        !  of small primes.
        !
        !  input
        !  lensav, the dimension of the wsave array.
        !  lensav must be at least 2*n + int(log(real(n))) + 4.
        !
        !  output,
        !  wsave(lensav), containing the prime factors
        !  of n and  also containing certain trigonometric values which will be used
        !  in routines cfft1b or cfft1f.
        !
        !  output
        !  ierror, error_flag.
        !  0, successful exit;
        !  2, input parameter lensav not big enough.
        !
        !--------------------------------------------------------------
        ! Dummy arguments
        !--------------------------------------------------------------
        integer (ip), intent (in)  :: n
        real (wp),    intent (out) :: wsave(lensav)
        integer (ip), intent (in)  :: lensav
        integer (ip), intent (out) :: ierror
        !--------------------------------------------------------------
        ! Local variables
        !--------------------------------------------------------------
        integer (ip) :: iw1, iw2
        !--------------------------------------------------------------

        !
        !==> Check validity of input arguments
        !
        if ( size(wsave) < get_complex_1d_saved_workspace_length(n) ) then
            ierror = 2
            call fft_error_handler('cfftmi ', 3)
        else
            ierror = 0
        end if

        !
        !==> Perform transform
        !
        if (n == 1) return

        ! Set workspace indices
        iw1 = 2*n+1
        iw2 = iw1 + 1

        call initialize_factors_and_trig_tables(n, wsave, wsave(iw1), wsave(iw2))

    end subroutine cfft1i

    module subroutine cfft2i(l, m, wsave, lensav, ierror)
        !
        ! cfft2i: initialization for cfft2b and cfft2f.
        !
        !  purpose:
        !
        !  cfft2i initializes real array wsave for use in its companion
        !  routines cfft2f and cfft2b for computing two-dimensional fast
        !  fourier transforms of complex data.  prime factorizations of l and m,
        !  together with tabulations of the trigonometric functions, are
        !  computed and stored in array wsave.
        !
        !  on 10 may 2010, this code was modified by changing the value
        !  of an index into the wsave array.
        !
        !  input
        !  integer l, the number of elements to be transformed
        !  in the first dimension. the transform is most efficient when l is a
        !  product of small primes.
        !
        !  integer m, the number of elements to be transformed
        !  in the second dimension. the transform is most efficient when m is a
        !  product of small primes.
        !
        !  integer lensav, the dimension of the wsave array.
        !  lensav must be at least
        !
        !  2*(l+m) + int(log(real(l)))+ int(log(real(m))) + 8.
        !
        !  output
        !  real wsave(lensav), contains the prime factors of l
        !  and m, and also certain trigonometric values which will be used in
        !  routines cfft2b or cfft2f.
        !
        !  integer  ierror, error_flag.
        !  0, successful exit;
        !  2, input parameter lensav not big enough;
        !  20, input error returned by lower level routine.
        !
        !--------------------------------------------------------------
        ! Dummy arguments
        !--------------------------------------------------------------
        integer (ip), intent (in)  :: l
        integer (ip), intent (in)  :: m
        real (wp),    intent (out) :: wsave(lensav)
        integer (ip), intent (in)  :: lensav
        integer (ip), intent (out) :: ierror
        !--------------------------------------------------------------
        ! Local variables
        !--------------------------------------------------------------
        integer (ip) local_error_flag
        !--------------------------------------------------------------

        !
        !==> Check validity of input arguments
        !
        if (size(wsave) < get_complex_2d_saved_workspace_length(l, m)) then
            ierror = 2
            call fft_error_handler('cfft2i', 4)
            return
        else
            ierror = 0
        end if

        associate( lnsv => get_complex_1d_saved_workspace_length(l) )

            call cfftmi(l, wsave(1:), lnsv, local_error_flag)

        end associate

        ! Check error_flag
        if ( local_error_flag /= 0) then
            ierror = 20
            call fft_error_handler('cfft2i',-5)
            return
        end if

        associate( &
            iw => get_1d_saved_workspace_length(l) - 1, &
            lnsv => get_1d_saved_workspace_length(m) &
            )

            call cfftmi(m, wsave(iw:), lnsv, local_error_flag)

        end associate

        ! Check error_flag
        if (local_error_flag /= 0) then
            ierror = 20
            call fft_error_handler('cfft2i',-5)
        end if

    end subroutine cfft2i

    module subroutine cfftmi(n, wsave, lensav, ierror)
        !
        !  cfftmi: initialization for cfftmb and cfftmf.
        !
        !  purpose:
        !
        !  cfftmi initializes array wsave for use in its companion routines
        !  cfftmb and cfftmf.  cfftmi must be called before the first call
        !  to cfftmb or cfftmf, and after whenever the value of integer n changes.
        !
        !  parameters:
        !
        !  integer n, the length of each sequence to be
        !  transformed.  the transform is most efficient when n is a product of
        !  small primes.
        !
        !  integer lensav, the dimension of the wsave array.
        !  lensav must be at least 2*n + int(log(real(n))) + 4.
        !
        !  real wsave(lensav), containing the prime factors
        !  of n and also containing certain trigonometric values which will be used in
        !  routines cfftmb or cfftmf.
        !
        !  integer ierror, error_flag.
        !  0, successful exit;
        !  2, input parameter lensav not big enough.
        !
        !--------------------------------------------------------------
        ! Dummy arguments
        !--------------------------------------------------------------
        integer (ip), intent (in)  :: n
        real (wp),    intent (out) :: wsave(lensav)
        integer (ip), intent (in)  :: lensav
        integer (ip), intent (out) :: ierror
        !--------------------------------------------------------------
        ! Local variables
        !--------------------------------------------------------------
        integer (ip) :: iw1, iw2
        !--------------------------------------------------------------

        !
        !==> Check validity of input arguments
        !
        if (lensav < get_complex_nd_saved_workspace_length(n) ) then
            ierror = 2
            call fft_error_handler('cfftmi ', 3)
            return
        else
            ierror = 0
        end if

        !
        !==>  Perform transform
        !
        if (n == 1) return

        ! Set workspace index pointer
        iw1 = 2*n+1
        iw2 = iw1 + 1

        call initialize_factors_and_trig_tables(n, wsave, wsave(iw1), wsave(iw2:))

    end subroutine cfftmi

    pure subroutine initialize_factors_and_trig_tables(n, wa, fnf, fac)
        !
        ! Purpose:
        !
        ! Sets up factors and trig tables
        !
        !--------------------------------------------------
        ! Dummy arguments
        !--------------------------------------------------
        integer (ip), intent (in)  :: n
        real (wp),    intent (out) :: wa(*)
        real (wp),    intent (out) :: fnf
        real (wp),    intent (out) :: fac(*)
        !--------------------------------------------------
        ! Dummy arguments
        !--------------------------------------------------
        integer (ip) :: ido, iip, iw, k1, l1, l2, nf
        !--------------------------------------------------

        !
        !==> Get the factorization of n.
        !
        call compute_factorization(n, nf, fac)
        fnf = real(nf, kind=wp)
        iw = 1
        l1 = 1
        !
        !==> Set up the trigonometric tables.
        !
        do k1 = 1, nf
            iip = int(fac(k1), kind=ip)
            l2 = l1 * iip
            ido = n / l2
            call compute_trig_table(ido, iip, wa(iw))
            iw = iw + (iip - 1) * (2*ido)
            l1 = l2
        end do

    end subroutine initialize_factors_and_trig_tables

    pure subroutine compute_factorization(n, nf, fac)
        !
        ! Purpose:
        !
        ! Factors of an integer for wp-kind float precision computations.
        !
        !  Parameters:
        !
        !  n, the number for which factorization and other information is needed.
        !
        !  nf, the number of factors.
        !
        !  fac(*), a list of factors of n.
        !
        !--------------------------------------------------
        ! Dummy arguments
        !--------------------------------------------------
        integer (ip), intent (in)  :: n
        integer (ip), intent (out) :: nf
        real (wp),    intent (out) :: fac(*)
        !--------------------------------------------------
        ! Local variables
        !--------------------------------------------------
        integer (ip) :: j, nl, nq, nr, ntry
        !--------------------------------------------------

        ntry = 0
        nl = n
        nf = 0
        j = 0

        do while (1 < nl)
            j = j + 1
            select case (j)
                case (1)
                    ntry = 4
                case (2:3)
                    ntry = j
                case (4)
                    ntry = 5
                case default
                    ntry = ntry + 2
            end select

            inner_loop: do
                nq = nl / ntry
                nr = nl - ntry * nq

                if ( nr /= 0 ) then
                    exit inner_loop
                end if

                nf = nf + 1
                fac(nf) = real(ntry, kind=wp)
                nl = nq

            end do inner_loop
        end do

    end subroutine compute_factorization

    pure subroutine compute_trig_table(ido, iip, wa)
        !
        ! Purpose:
        !
        ! Computes trigonometric tables for wp-kind float precision arithmetic.
        !
        !--------------------------------------------------
        ! Dummy arguments
        !--------------------------------------------------
        integer (ip), intent (in)  :: ido
        integer (ip), intent (in)  :: iip
        real (wp),    intent (out) :: wa(ido,iip-1,2)
        !--------------------------------------------------
        ! Local variables
        !--------------------------------------------------
        integer (ip)         :: i, j !! Counters
        real (wp), parameter :: TWO_PI = TWO * acos(-ONE)
        real (wp)            :: argz, arg1, arg2, arg3, arg4
        !--------------------------------------------------

        argz = TWO_PI/iip
        arg1 = TWO_PI/( ido * iip)

        do j = 2, iip
            arg2 = real(j - 1, kind=wp) * arg1
            do i = 1, ido
                arg3 = real(i - 1, kind=wp) * arg2
                wa(i,j-1,1) = cos(arg3)
                wa(i,j-1,2) = sin(arg3)
            end do
            if (5 < iip) then
                arg4 = real(j - 1, kind=wp) * argz
                wa(1,j-1,1) = cos(arg4)
                wa(1,j-1,2) = sin(arg4)
            end if
        end do

    end subroutine compute_trig_table

end submodule complex_initialization_routines


text:
*	$Id: ffxf0.f,v 1.5 1996/02/07 09:28:49 gj Exp $
*###[ ffxf0:
	subroutine ffxf0(cf0,ce0i,cd0ij,xpi,ier)
***#[*comment:***********************************************************
*									*
*	calculate							*
*									*
*	      1  /   /						     \-1*
*	f0= -----\dq |(q^2-m_1^2)((q+p_1)^2-m_2^2)...((q-p_6)^2-m_6^2|	*
*	    ipi^2/   \						     /	*
*									*
*	following the six five-point-function method in ....		*
*	As an extra the ten fourpoint functions Dij are also returned	*
*	plus the six fivepoint functions Ei.				*
*									*
*	Input:	xpi = m_i^2	(real)	  i=1,6				*
*		xpi = p_i.p_i	(real)	  i=7,12 (note: B&D metric)	*
*		xpi = (p_i+p_{i+1})^2 (r) i=13,18			*
*		xpi = (p_i+p_{i+1}+p_{i+3})^2 (r)  i=19,21		*
*									*
*	Output:	cf0		(complex) F0				*
*		ce0i(6)		(complex) E0 with s_i missing		*
*		cd0ij(6,6)	(complex) D0 with s_i and s_j missing	*
*		ier		(integer) no of digits lost, >50 error	*
*									*
***#]*comment:***********************************************************
*  #[ declarations:
	implicit none
*
*	arguments
*
	DOUBLE PRECISION xpi(21)
	DOUBLE COMPLEX cf0,ce0i(6),cd0ij(6,6)
	integer ier
*
*	local variables
*
	integer i,j,ier0
	DOUBLE PRECISION dpipj(21,21)
*
*	common blocks:
*
	include 'ff.h'
*
*  #] declarations:
*  #[ get differences:
*
*	simulate the differences in the masses etc..
*
	if ( lwrite ) then
	    print *,'ffxf0: input xpi: '
	    print '(i3,e24.16)',(i,xpi(i),i=1,21)
	endif
*
*	no redundant input yet (may be necessary)
*
*
*	the differences
*
	ier0 = 0
	if ( lwarn ) then
	    do 20 i=1,21
		dpipj(i,i) = 0
		do 10 j=1,i-1
		    dpipj(j,i) = xpi(j) - xpi(i)
		    dpipj(i,j) = -dpipj(j,i)
		    if ( abs(dpipj(j,i)) .lt. xloss*abs(xpi(i))
     +			.and. xpi(i) .ne. xpi(j) ) then
			call ffwarn(193,ier0,dpipj(j,i),xpi(i))
			if ( lwrite ) print *,'between xpi(',i,
     +				') and xpi(',j,')'
		    endif
   10		continue
   20	    continue
	else
	    do 40 i=1,21
		do 30 j=1,21
		    dpipj(j,i) = xpi(j) - xpi(i)
   30		continue
   40	    continue
	endif
*  #] get differences:
*  #[ call ffxf0a:
	call ffxf0a(cf0,ce0i,cd0ij,xpi,dpipj,ier)
*  #] call ffxf0a:
*###] ffxf0:
	end
*###[ ffxf0a:
	subroutine ffxf0a(cf0,ce0i,cd0ij,xpi,dpipj,ier)
***#[*comment:***********************************************************
*									*
*	calculate							*
*									*
*	      1  /   /						     \-1*
*	f0= -----\dq |(q^2-m_1^2)((q+p_1)^2-m_2^2)...((q-p_5)^2-m_5^2|	*
*	    ipi^2/   \						     /	*
*									*
*	following the five four-point-function method in ....		*
*	As an extra the five fourpoint function Di are also reurned	*
*									*
*	Input:	xpi = m_i^2	   (real)  i=1,6			*
*		xpi = p_i.p_i	   (real)  i=7,12 (note: B&D metric)	*
*		xpi = (p_i+p_{i+1})^2 (r)  i=13,18			*
*		xpi = (p_i+p_{i+1}+p_{i+2})^2 (r)  i=19,21		*
*		dpipj(21,21)	   (real)  = pi(i) - pi(j)		*
*									*
*	Output:	cf0		  (complex)				*
*		ce0i(6)		  (complex) E0 with s_i missing		*
*		cd0ij(6,6)	  (complex) D0 with s_i,s_j missing	*
*		ier		  (integer) <50:lost # digits 100=error	*
*									*
***#]*comment:***********************************************************
*  #[ declarations:
	implicit none
*
*	arguments
*
	integer ier
	DOUBLE COMPLEX cf0,ce0i(6),cd0ij(6,6)
	DOUBLE PRECISION xpi(21),dpipj(21,21)
*
*	local variables
*
	integer i,j,k,l,m,ii(10),ier2,ier1,ier0,irota,itype,ndiv,idum,
     +		idone,ii4(6),is
	logical lwsav,ldel2s,lwhich
	DOUBLE COMPLEX c,cfac,cs,cd0i(5),csum,csi(7)
	DOUBLE PRECISION del6,xpi4(13),dpipj4(10,13),piDpj4(10,10),
     +		absc,xmax,piDpj(21,21),xqi4(13),dqiqj4(10,13),
     +		qiDqj4(10,10),del2s,xmx4(6,6),dl4rij(6,6),xpi5(20),
     +		dpipj5(15,20),piDpj5(15,15),dl4ri(5),dl5ri(6),xlosn,
     +		d5sp,dl4q(6),psum
	save ii4
*
*	common blocks:
*
	include 'ff.h'
*
*	statement function
*
	absc(c) = abs(DBLE(c)) + abs(DIMAG(c))
*
*	data
*
	data ii4 /5,6,7,8,9,10/
!$omp threadprivate(ii4)

*  #] declarations:
*  #[ initialisations:
*
	ndiv = 0
	idsub = 0
	cf0 = 0
	do 2 i=1,6
	    ce0i(i) = 0
	    do 1 j=1,6
		cd0ij(i,j) = 0
    1	    continue
    2	continue
*
*  #] initialisations:
*  #[ get dot products:
*
	idsub = idsub + 1
	call ffdot6(piDpj,xpi,dpipj,ier)
	if ( ldot ) then
	    do 6 i=1,21
		do 5 j=1,21
		    fpij6(j,i) = piDpj(j,i)
    5		continue
    6	    continue
	continue
	endif
	if ( ltest ) then
	    ii(1) = 7
	    ii(2) = 8
	    ii(3) = 9
	    ii(4) = 10
	    ii(5) = 11
	    call ffdl5p(xpi,piDpj,21,ii,ier)
	endif
	if ( lwrite ) print *,'After dotproducts ier = ',ier
*
*  #] get dot products:
*  #[ five and four point stuff:
*
	ier2 = ier
	do 100 i=1,6
*
*	    get the five-point momenta
*
	    ier1 = ier
	    call ffpi65(xpi5,dpipj5,piDpj5,xpi,dpipj,piDpj,i,ier1)
*
*	    get fourpoint functions
*
	    do 90 k=1,5
		j=k
		if ( lwrite ) print '(a,2i2,a)',
     +		    '####[ ffxf0a: fourpoint function nr ',i,j+1,': '
		if ( k.lt.i ) then
*		    we already have it
		else
		    j = j+1
		    ier0 = ier
*
*		    get four-point momenta
*
		    call ffpi54(xpi4,dpipj4,piDpj4,xpi5,dpipj5,piDpj5,k,
     +			ier0)
		    if ( ltest ) then
			idum = ier
			call ffpi64(xqi4,dqiqj4,qiDqj4,xpi,dpipj,piDpj,
     +				i,j,idum)
			xlosn = xloss*DBLE(10)**(-mod(ier0,50))
			do 12 l=1,13
			    if ( xlosn*abs(xpi4(l)-xqi4(l)).gt.precx*abs
     +			      (xpi4(l)) ) print*,'ffxf0a: error: xpi4(',
     +			      l,') != xqi4(',l,'): ',xpi4(l),xqi4(l)
			    do 11 m=1,10
				if ( xlosn*abs(dpipj4(m,l)-dqiqj4(m,l))
     +				    .gt.precx*abs(xpi4(l)) ) print *,
     +				    'ffxf0a: error: dpipj4(',m,l,') !=',
     +				    ' dqiqj4(',m,l,'): ',dpipj4(m,l),
     +				    dqiqj4(m,l),dpipj4(m,l)-dqiqj4(m,l)
   11			    continue
   12			continue
			do 14 l=1,10
			    do 13 m=1,10
				if ( piDpj4(m,l).ne.qiDqj4(m,l) ) print
     +				    *,'ffxf0a: error: piDpj4(',m,l,
     +				    ') != qiDqj4(',m,l,'): ',piDpj4(m,
     +				    l),qiDqj4(m,l)
   13			    continue
   14			continue
		    endif
		    ier1 = ier0
		    call ffxdir(cs,cfac,idone,xpi4,dpipj4,6,ndiv,ier1)
		    if ( idone .gt. 0 ) then
*			done
			xmax = abs(cs)*10d0**(-mod((ier1-ier0),50))
		    else
			ier1 = ier0
*
*			rotate to calculable posistion
*
			call ffrot4(irota,del2s,xqi4,dqiqj4,qiDqj4,xpi4,
     +				dpipj4,piDpj4,5,itype,ier0)
			if ( itype .lt. 0 ) then
			    print *,'ffxf0:  error:  Cannot handle '//
     +				'this 4point masscombination yet:'
			    print *,(xpi(j),j=1,20)
			    return
			endif
			if ( itype .eq. 1 ) then
			    ldel2s = .TRUE.
			    isgnal = +1
			    print *,'ffxf0a: Cannot handle del2s=0 yet'
			    stop
			else
			    ldel2s = .FALSE.
			endif
			if ( itype .eq. 2 ) then
			    print *,'ffxf0a: Cannot handle doubly IR ',
     +				'divergent function yet'
			    stop
			endif
*
*			get fourpoint function
*
			if ( lwrite ) then
			    print *,'xpi for ffxd0e: '
			    print '(i3,e24.16)',(m,xqi4(m),m=1,13)
			endif
			lwsav = lwrite
			lwrite = .FALSE.
			idsub = idsub + 1
			call ffxd0e(cs,cfac,xmax, .FALSE.,ndiv,xqi4,
     +				dqiqj4,qiDqj4,del2s,ldel2s,ier0)
			lwrite = lwsav
			ier1 = max(ier1,ier0)
		    endif
		    cd0ij(i,j) = cs*cfac
		    cd0ij(j,i) = cd0ij(i,j)
		    xmx4(i,j) = xmax*absc(cfac)
		    xmx4(j,i) = xmx4(i,j)
		    if ( ldot ) then
			call ffdl3p(fdl3ij(i,j),qiDqj4,10,ii4,ii4,
     +				ier0)
			fdl3ij(j,i) = fdl3ij(i,j)
			fd4sij(i,j) = fdel4s
			fd4sij(j,i) = fdel4s
*			let's check that these have been set by ffxd0e...
			if ( ltest ) then
			    ier0 = 0
			    call ffdel4(fdel4s,xpi4,piDpj4,10,ier0)
			    if ( xloss*10d0**(-ier0-1)*abs(fd4sij(i,j)-
     +			    	    fdel4s).gt.precx*abs(fdel4s) ) then
			    	print *,'ffxf0a: error: Del4s was not'//
     +			    		' correct',fd4sij(i,j),fdel4s,
     +			    		fd4sij(i,j)-fdel4s,ier0
			    endif
			endif
		    endif
		endif
*
*		get the coefficient determinant (not symmetric!)
*
		idsub = idsub + 1
		ier0 = ier
		call ffdl4r(dl4rij(i,j),xpi5,piDpj5,15,k,ier0)
		ier1 = max(ier1,ier0)
*
*		and fill the five-point linear arrays
*
		cd0i(k) = cd0ij(i,j)
		dl4ri(k) = dl4rij(i,j)
		if ( lwrite ) then
		    print '(a,2i2,a)',
     +			'####] ffxf0a: fourpoint function nr ',i,j,': '
		    print *,'dl4rij(',i,j,') = ',dl4rij(i,j)
		    print *,'cd0ij(',i,j,')  = ',cd0ij(i,j),xmx4(i,j),
     +			ier0
		endif
   90	    continue
*
*	    call ffxe00
*
	    if ( lwrite ) print '(a,i2,a)',
     +		'####[ ffxf0a: fivepoint function nr ',i,': '
	    call ffxe00(ce0i(i),cd0i,dl4ri,xpi5,piDpj5,ier1)
	    if ( lwrite ) print '(a,i2,a)',
     +		'####] ffxf0a: fivepoint function nr ',i,': '
	    if ( lwrite ) print *,'ce0i(',i,') = ',ce0i(i),ier1
	    if ( ldot ) fdl4i(i) = fdel4
	    ier2 = max(ier2,ier1)
  100	continue
	ier = ier2
	if ( lwrite ) print *,'after E0s ier = ',ier
*
*  #] five and four point stuff:
*  #[ six point stuff:
*
	ier1 = 0
	call ffdel6(del6,xpi,piDpj,21,ier1)
	csum = 0
	xmax = 0
	do i=1,6
	    if ( ce0i(i) .ne. 0 ) then
		ier0 = 0
		call ffdl5r(dl5ri(i),xpi,piDpj,21,i,ier0)
		csum = csum + DBLE(dl5ri(i))*ce0i(i)
		xmax = max(xmax,absc(csum))
		ier1 = max(ier1,ier0)
	    endif
	enddo
	ier = max(ier,ier1)
*
*	Check for cancellations in the final adding up
*
	if ( lwarn .and. 2*absc(csum) .lt. xloss*xmax )
     +		call ffwarn(191,ier,absc(csum),xmax)
*
*	Check for a sum close to the minimum of the range (underflow
*	problems)
*
	if ( lwarn .and. absc(csum).lt.xalogm/precc .and. csum.ne.0 )
     +		call ffwarn(192,ier,absc(csum),xalogm/precc)
*
*	If the imaginary part is very small it most likely is zero
*	(can be removed, just esthetically more pleasing)
*
	if ( abs(DIMAG(csum)) .lt. precc*abs(DBLE(csum)) )
     +		csum = DCMPLX(DBLE(csum))
*
*	Finally ...
*
	cf0 = csum*DBLE(-1/(2*del6))
*
*  #] six point stuff:
*  #[ print output:
	if ( lwrite ) then
	    print *,'ffxf0a: cf0 = ',cf0,ier
	endif
	if ( .FALSE. ) then
	    do i=1,6
		if ( xpi(i).eq.0 ) then
*		    assume it's a photon
		    psum = 0
		    do j=1,5
			k = i+j-1
			if ( k.gt.6 ) k = k-6
			psum = psum + xpi(k+6)
			do l=1,j-1
			    m = i+l-1
			    if ( m.gt.6 ) m = m-6
			    psum = psum + 2*piDpj(k+6,m+6)
			enddo
			k = i+j
			if ( k.gt.6 ) k = k-6
			if ( abs(piDpj(i,k)).gt.xloss*abs(xpi(k)) ) then
			    print *,'ratio coeffs ',k,' is ',dl5ri(k)/
     +				(-2*del6)
			    print *,'propagator   ',k,' is ',1/(psum-xpi
     +				(k))
			endif
		    enddo
		endif
	    enddo
	endif
*  #] print output:
*###] ffxf0a:
	end
*###[ ffxf0r:
	subroutine ffxf0r(cf0,ce0i,cd0ij,xpi,ier)
***#[*comment:***********************************************************
*									*
*	Tries all 12 easy permutations of the 6pointfunction		*
*									*
***#]*comment:***********************************************************
*  #[ declarations:
	implicit none
	integer ier,nrot
	parameter(nrot=12)
	DOUBLE PRECISION xpi(21),xqi(21)
	DOUBLE COMPLEX cf0,ce0i(6),cd0ij(6,6),cf0p,ce0ip(6),cd0ijp(6,6),
     +		ce0iq(6),cd0ijq(6,6)
	integer inew(21,nrot),irota,ier1,i,j,k,l,icon,ialsav,init
	parameter(icon=3)
	logical lcon
	save inew,init,lcon
	include 'ff.h'
	data inew
     +	/1,2,3,4,5,6, 7,8,9,10,11,12, 13,14,15,16,17,18, 19,20,21,
     +	 2,3,4,5,6,1, 8,9,10,11,12,7, 14,15,16,17,18,13, 20,21,19,
     +	 3,4,5,6,1,2, 9,10,11,12,7,8, 15,16,17,18,13,14, 21,19,20,
     +	 4,5,6,1,2,3, 10,11,12,7,8,9, 16,17,18,13,14,15, 19,20,21,
     +	 5,6,1,2,3,4, 11,12,7,8,9,10, 17,18,13,14,15,16, 20,21,19,
     +	 6,1,2,3,4,5, 12,7,8,9,10,11, 18,13,14,15,16,17, 21,19,20,
     +	 6,5,4,3,2,1, 11,10,9,8,7,12, 16,15,14,13,18,17, 21,20,19,
     +	 5,4,3,2,1,6, 10,9,8,7,12,11, 15,14,13,18,17,16, 20,19,21,
     +	 4,3,2,1,6,5, 9,8,7,12,11,10, 14,13,18,17,16,15, 19,21,20,
     +	 3,2,1,6,5,4, 8,7,12,11,10,9, 13,18,17,16,15,14, 21,20,19,
     +	 2,1,6,5,4,3, 7,12,11,10,9,8, 18,17,16,15,14,13, 20,19,21,
     +	 1,6,5,4,3,2, 12,11,10,9,8,7, 17,16,15,14,13,18, 19,21,20/
	data init /0/
!$omp threadprivate(inew,init,lcon)

*  #] declarations:
*  #[ open console for some activity on screen:
	if ( init .eq. 0 ) then
	    init = 1
	    if ( lwrite ) then
		open(icon,file='CON:',status='old',err=11)
		lcon = .TRUE.
		goto 13
	    endif
   11	    continue
	    lcon = .FALSE.
   13	    continue
	endif
*  #] open console for some activity on screen:
*  #[ calculations:
	cf0 = 0
	ier = 999
	ialsav = isgnal
	    do 20 irota=1,nrot
		do 10 i=1,21
		    xqi(inew(i,irota)) = xpi(i)
   10		continue
		print '(a,i2,a,i2)','---#[ rotation ',irota,': isgnal ',
     +			isgnal
		if (lcon) write(icon,'(a,i2,a,i2)')'rotation ',irota,',
     +			isgnal ',isgnal
		ier1 = 0
		ner = 0
		id = id + 1
		isgnal = ialsav
		call ffxf0(cf0p,ce0ip,cd0ijp,xqi,ier1)
		ier1 = ier1 + ner
		if ( ier.gt.5 ) call ffwarn(998,ier,x0,x0)
		print '(a,i2,a,i2)','---#] rotation ',irota,': isgnal ',
     +			isgnal
		print '(a,2g28.16,i3)','f0 = ',cf0p,ier1
		do 15 k=1,6
		    ce0iq(k) = ce0ip(inew(k,irota))
		    print '(a,2g28.16,i3)','e0 = ',ce0iq(k),k
   15		continue
		do 17 k=1,6
		    do 16 l=k+1,6
			cd0ijq(l,k)=cd0ijp(inew(l,irota),inew(k,irota))
			print '(a,2g28.16,2i3)','d0 = ',cd0ijq(l,k),l,k
   16		    continue
   17		continue
		if (lcon) write(icon,'(a,2g28.16,i3)')'f0 = ',cf0p,ier1
		if ( ier1 .lt. ier ) then
		    cf0 = cf0p
		    do 19 k=1,6
			ce0i(k) = ce0iq(k)
			do 18 l=k+1,6
			    cd0ij(l,k) =
     +				cd0ijp(inew(l,irota),inew(k,irota))
			    cd0ij(k,l) = cd0ij(l,k)
   18			continue
   19		    continue
		    ier = ier1
		endif
   20	    continue
*  #] calculations:
*###] ffxf0r:
	end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
!========================================================================
!                       LAMBDA OPERATOR MODULE
!
! For the Lambda Iteration scheme (and its sister, the Accelerated
! Lambda Iteration scheme) a Lambda Operator is required.
!
! A Lambda operator is simply an operator (subroutine) that computes
! the mean intensity at some wavelength using the integrals over a large
! number of rays. In contrast to the mean intensity subroutine in the
! montecarlo_module, for the Lambda operator we do not "follow a photon
! along its way" but we instead "integrate the formal transfer equation
! along a ray". If no an-isotropic dust scattering is included in the
! problem, and the scattering source function is known, this gives a
! more accurate answer. Also, depending on which method of Lambda
! operator you use, you get less noise. Finally, one can use the
! Accelerated Lambda Iteration scheme to speed up convergence of the
! Lamber Iteration.
!========================================================================
module lambdaoper_module
use mathroutines_module
use rtglobal_module
use montecarlo_module
use amr_module
use amrray_module
use constants_module

double precision :: lambdaoper_shift = 1.d90
integer :: lambdaoper_nfreq,lambdaoper_istar,lambdaoper_inu
double precision, allocatable :: lambdaoper_intensity(:),lambdaoper_bgint(:)
integer :: lambdaoper_order=1
logical :: lambdaoper_addtocell_meanint=.false.
logical :: lambdaoper_vertexbased=.false.
integer :: lambdaoper_random_ray_method=0
double precision :: lambdaoper_local_lamdiag

doubleprecision :: lambdaoper_snu_prev,lambdaoper_anu_prev,lambdaoper_jnu_prev
doubleprecision :: lambdaoper_snu_curr,lambdaoper_anu_curr,lambdaoper_jnu_curr

double precision, allocatable :: lambdaoper_meanint_cell(:)
double precision, allocatable :: lambdaoper_lambda_diag_cell(:)
double precision, allocatable :: lambdaoper_dstot_cell(:)
double precision, allocatable :: lambdaoper_meanint_vertex(:)
double precision, allocatable :: lambdaoper_lambda_diag_vertex(:)
double precision, allocatable :: lambdaoper_dstot_vertex(:)

contains

!------------------------------------------------------------------------
!                         INIT THIS MODULE
!------------------------------------------------------------------------
subroutine lambdaoper_init(nf,freq,vertexbased,bgint)
  implicit none
  integer :: nf
  double precision :: freq(1:nf)
  double precision, optional :: bgint(1:nf)
  logical :: vertexbased
  !
  ! First clean up
  !
  call lambdaoper_cleanup()
  !
  ! Then copy some stuff
  !
  lambdaoper_vertexbased = vertexbased
  lambdaoper_nfreq = nf
  !
  ! Allocate and init some arrays
  !
  allocate(lambdaoper_intensity(1:nf),lambdaoper_bgint(1:nf))
  lambdaoper_intensity(:) = 0.d0
  if(present(bgint)) then
     lambdaoper_bgint(:) = bgint(:)
  else
     lambdaoper_bgint(:) = 0.d0
  endif
  !
  ! Since we want to compute the mean intensity, here is the array for that
  !
  if(vertexbased) then
     allocate(lambdaoper_meanint_vertex(amr_nr_vertices_max))
     allocate(lambdaoper_lambda_diag_vertex(amr_nr_vertices_max))
     allocate(lambdaoper_dstot_vertex(amr_nr_vertices_max))
  else
     allocate(lambdaoper_meanint_cell(amr_nrleafs_max))
     allocate(lambdaoper_lambda_diag_cell(amr_nrleafs_max))
     allocate(lambdaoper_dstot_cell(amr_nrleafs_max))
  endif
end subroutine lambdaoper_init

!------------------------------------------------------------------------
!                        CLOSE THIS MODULE
!------------------------------------------------------------------------
subroutine lambdaoper_cleanup()
  implicit none
  if(allocated(lambdaoper_intensity)) deallocate(lambdaoper_intensity)
  if(allocated(lambdaoper_bgint)) deallocate(lambdaoper_bgint)
  if(allocated(lambdaoper_meanint_cell)) deallocate(lambdaoper_meanint_cell)
  if(allocated(lambdaoper_meanint_vertex)) deallocate(lambdaoper_meanint_vertex)
  if(allocated(lambdaoper_dstot_cell)) deallocate(lambdaoper_dstot_cell)
  if(allocated(lambdaoper_dstot_vertex)) deallocate(lambdaoper_dstot_vertex)
  if(allocated(lambdaoper_lambda_diag_cell)) deallocate(lambdaoper_lambda_diag_cell)
  if(allocated(lambdaoper_lambda_diag_vertex)) deallocate(lambdaoper_lambda_diag_vertex)
  lambdaoper_nfreq = 0
  lambdaoper_local_lamdiag = -1d99
end subroutine lambdaoper_cleanup

!------------------------------------------------------------------------
!              HELPER ROUTINE FOR RT INTEGRATION STEP
!------------------------------------------------------------------------
subroutine lambdaoper_integrate_element()
  implicit none
  integer :: inu
  double precision :: dtau1,theomax,e0,e1,ax,bx,xp,qdr
  double precision :: src(1),alp(1)
  double precision :: r,theta,phi
  !
  ! Check if we are/were in a star
  !
  if(lambdaoper_istar.le.0) then
     !
     ! Normal situation: we are not inside a star (istar=0), or we will
     ! hit the star but are not inside the star yet (istar<0).
     !
     ! Decide whether we do 1st order, 2nd order or 3rd order integration
     !
     select case(lambdaoper_order)
     case(-1)
        !
        ! Integrate only the optical depth (useful for escape probability)
        !
        ! Get the src and alp values
        !
        call sources_get_src_alp(lambdaoper_inu,lambdaoper_inu,1,src,alp,##########)
        !
        ! Now do the RT, including everything
        !
        lambdaoper_intensity(lambdaoper_inu) = lambdaoper_intensity(lambdaoper_inu) + alp(1) * ray_ds
        !
     case(1)
        !
        ! First order integration. Cell-based scheme.
        !
        ! Get the src and alp values
        !
        call sources_get_src_alp(lambdaoper_inu,lambdaoper_inu,1,src,alp,########)
        !
        ! Compute the exp(-tau)
        !
        xp = alp(1) * ray_ds
        if(xp.lt.0.d0) stop 7329
        if(xp.gt.1d-4) then
           xp  = exp(-xp)
           lambdaoper_local_lamdiag = 1.d0 - xp
        else
           lambdaoper_local_lamdiag = xp
           xp  = exp(-xp)
        endif
        !
        ! Now do the RT, including everything
        !
        lambdaoper_intensity(lambdaoper_inu) = xp * lambdaoper_intensity(lambdaoper_inu) + lambdaoper_local_lamdiag * src(1)
        !
     case(2)
        !
        ! Second order integration. Corner-based (=vertex-based) scheme.
        ! Here all sources have all been pre-calculated on
        ! the cell corners.
        !
        ! Backup snu and anu at previous crossing
        !
        lambdaoper_snu_prev = lambdaoper_snu_curr
        lambdaoper_anu_prev = lambdaoper_anu_curr
        lambdaoper_jnu_prev = lambdaoper_jnu_curr
        !
        ! Find snu and anu at new crossing
        !
        ! NOTE: If lambdaoper_interpol_jnu is set, then snu is actually
        !       j_nu (the emissivity). It will be translated below.
        !
        if(igrid_type.lt.100) then
           if(igrid_coord.lt.100) then
              !
              ! Cartesian
              !
              call lambdaoper_find_srcalp_interpol(       &
                   ray_cart_x,ray_cart_y,ray_cart_z,  &
                   lambdaoper_snu_curr,lambdaoper_anu_curr)
           elseif(igrid_coord.lt.200) then
              !
              ! Spherical
              !
              call amr_xyz_to_rthphi(ray_cart_x,ray_cart_y,ray_cart_z, &
                                           r,theta,phi)
              call lambdaoper_find_srcalp_interpol(    &
                   r,theta,phi,lambdaoper_snu_curr,lambdaoper_anu_curr)
           else
              stop 498
           endif
        else
           stop 499
        endif
        !
        ! If lambdaoper_interpol_jnu is set, then what we get back
        ! from the camera_find_srcalp_interpol() is in fact the
        ! emissivity j_nu instead of the source function S_nu.
        ! So translate either snu to jnu or vice versa.
        !
        if(lambdaoper_interpol_jnu) then
           jnu_curr = snu_curr
           snu_curr = jnu_curr / (anu_curr+1d-99)
        else
           jnu_curr = snu_curr * (anu_curr+1d-99)
        endif
        !
        ! Now do the integration step
        !
        if(ray_index.gt.0) then
           dtau1     = 0.5d0 * ( lambdaoper_anu_prev + lambdaoper_anu_curr ) * ray_ds
           theomax   = 0.5d0 * ( lambdaoper_jnu_prev + lambdaoper_jnu_curr ) * ray_ds
           if(dtau1.gt.1.d-6) then
              xp = exp(-dtau1)
              e0 = 1.d0 - xp
              e1 = dtau1 - e0
              bx = e1 / dtau1
              ax = e0 - bx
           else
              ax = 0.5d0 * dtau1
              bx = 0.5d0 * dtau1
              xp = 1.d0 - dtau1
           endif
           if(dtau1.gt.1e-9) then
              qdr = ax * snu_prev + bx * snu_curr
           else
              qdr = theomax
           endif
           qdr = min(qdr,theomax)
           lambdaoper_intensity(lambdaoper_inu) = xp * lambdaoper_intensity(lambdaoper_inu) + qdr
        endif
     case(3)
        !
        ! Third order integration
        !
        write(*,*) 'ERROR: Third-order integration not yet implemented'
        stop
     end select
  elseif(camera_incl_stars.ne.0) then
     !
     ! Special situation: we are inside a star.
     ! So set the intensity to that of the star.
     !
     lambdaoper_intensity(lambdaoper_inu) = find_starlight_interpol(sources_frequencies(lambdaoper_inu), &
                                                    lambdaoper_istar)
  endif
end subroutine lambdaoper_integrate_element

!------------------------------------------------------------------------
!            INTEGRATE FORMAL TRANSFER EQUATION ALONG A RAY
!------------------------------------------------------------------------
subroutine lambdaoper_integrate_ray()
  implicit none
  double precision :: s,ds
  doubleprecision :: snu_prev,anu_prev,jnu_prev,snu_curr,anu_curr,jnu_curr
  !
  ! Check
  !
  if(lambdaoper_shift.gt.0.d0) then
     write(stdo,*) 'ERROR in Lambda operator: Shift>0.'
     stop
  endif
  if(.not.allocated(lambdaoper_intensity)) then
     write(stdo,*) 'ERROR in Lambda operator: lambdaoper_intensity not allocated.'
     stop
  endif
  !
  ! Start the intensity
  !
  lambdaoper_intensity(inu0:inu1) = lambdaoper_bgint(inu0:inu1)
  !
  ! Set defaults
  !
  snu_prev = 0.d0
  anu_prev = 0.d0
  jnu_prev = 0.d0
  snu_curr = 0.d0
  anu_curr = 0.d0
  jnu_curr = 0.d0
  !
  ! Check the dust scattering mode, because this only works under
  ! certain conditions.
  !
  if(scattering_mode.ge.2) then
     write(stdo,*) 'ERROR in Lambda module: The Lambda Operator only'
     write(stdo,*) '      works with dust scattering if the scattering'
     write(stdo,*) '      is isotropic.'
     stop
  endif
  if(scattering_mode.ge.1) then
     if(.not.allocated(mcscat_scatsrc_iquv)) stop 3005
  endif
  !
  ! If mirror symmetry is switched on, then we must assure that
  ! z is above the midplane
  !
  if(amrray_mirror_equator) then
     if(ray_cart_z.eq.0.d0) then
        ray_cart_dirz = abs(ray_cart_dirz)
     elseif(ray_cart_z.lt.0.d0) then
        ray_cart_z    = -ray_cart_z
        ray_cart_dirz = -ray_cart_dirz
     endif
  endif
  !
  ! We start outside the domain, so all indices must be 0
  !
  ray_index     = 0
  ray_indexnext = 0
  !
  ! Now trace
  !
  arrived = .false.
  do while(.not.arrived)
     !
     ! Back up the current position
     !
     ray_prev_x     = ray_cart_x
     ray_prev_y     = ray_cart_y
     ray_prev_z     = ray_cart_z
     !
     ! Find new position
     !
     if(igrid_type.lt.100) then
        !
        ! We use a regular or AMR grid
        !
        if(igrid_coord.lt.100) then
           !
           ! We use cartesian coordinates
           !
           call amrray_find_next_location_cart(ray_dsend,            &
                ray_cart_x,ray_cart_y,ray_cart_z,                    &
                ray_cart_dirx,ray_cart_diry,ray_cart_dirz,           &
                ray_index,ray_indexnext,ray_ds,arrived,              &
                levelnext=levelnext)
           !
           ! If we went through one of the boundaries of the
           ! domain, we might need to implement the thermal
           ! boundary, if present
           !
           if(incl_thermbc.ne.0) then
              if((ray_index.eq.0).and.(ray_indexnext.ne.0)) then
                 !
                 ! First determine which boundary we just went
                 ! through
                 !
                 if(amrray_icross.ge.0) stop 3902  ! Self-consistency check
                 if(amrray_icross.eq.-1) then
                    bc_idir = 1
                    bc_ilr  = 1
                 elseif(amrray_icross.eq.-2) then
                    bc_idir = 1
                    bc_ilr  = 2
                 elseif(amrray_icross.eq.-3) then
                    bc_idir = 2
                    bc_ilr  = 1
                 elseif(amrray_icross.eq.-4) then
                    bc_idir = 2
                    bc_ilr  = 2
                 elseif(amrray_icross.eq.-5) then
                    bc_idir = 3
                    bc_ilr  = 1
                 elseif(amrray_icross.eq.-6) then
                    bc_idir = 3
                    bc_ilr  = 2
                 else
                    stop 4944
                 endif
                 !
                 ! Now set the intensity to the thermal
                 ! emission of the boundary, if the
                 ! thermal boundary is switched on
                 !
                 if(thermal_bc_active(bc_ilr,bc_idir)) then
                    do inu=1,lambdaoper_nfreq
                       lambdaoper_intensity(inu) = bplanck(thermal_bc_temp(bc_ilr,bc_idir), &
                                                lambdaoper_frequencies(inu))
                    enddo
                 endif
              endif
           endif
           !
        elseif(igrid_coord.lt.200) then
           !
           ! We use spherical coordinates
           !
           call amrray_find_next_location_spher(ray_dsend,           &
                ray_cart_x,ray_cart_y,ray_cart_z,                    &
                ray_cart_dirx,ray_cart_diry,ray_cart_dirz,           &
                ray_index,ray_indexnext,ray_ds,arrived)
        else
           write(stdo,*) 'ERROR: Cylindrical coordinates not yet implemented'
           stop
        endif
        !
        ! Set the lambdaoper_istar, which is used for treating finite-size
        ! stars (instead of point source stars)
        !
        lambdaoper_istar = amrray_ispherehit
        !
     else
        write(stdo,*) 'SORRY: Delaunay or Voronoi grids not yet implemented'
        stop
     endif
     !
     ! Path length
     !
     ds    = sqrt( (ray_cart_x-ray_prev_x)**2 + (ray_cart_y-ray_prev_y)**2 + (ray_cart_z-ray_prev_z)**2 )
     !
     ! Now do the integration along this path element
     !
     call lambdaoper_integrate_element()
     !
     ! Add information to the cell
     !
     ! NOTE: We add the final intensity to the cell, even though perhaps
     !       the average would be better. But the Lambda_diag approximate
     !       operator is calculated on the basis of the final intensity.
     !
     if(lambdaoper_addtocell_meanint) then
        lambdaoper_meanint_cell(ray_index) = lambdaoper_meanint_cell(ray_index) + &
             ds * lambdaoper_intensity(lambdaoper_inu)
        lambdaoper_lambda_diag_cell(ray_index) = lambdaoper_lambda_diag_cell(ray_index) + &
             ds * lambdaoper_local_lamdiag
        lambdaoper_dstot_cell(ray_index) = lambdaoper_dstot_cell(ray_index) + ds
     endif
     !
     ! Increase s
     !
     s = s + ds
     !
     ! Now make next cell the current cell
     !
     ray_index = ray_indexnext
     !
     ! Do the same for the pointers.
     ! This is grid type dependent.
     !
     if(igrid_type.lt.100) then
        !
        ! AMR type grid
        !
        if(ray_index.gt.0) then
           !
           ! We are in a cell
           !
           if(amr_tree_present) then
              amrray_cell => amr_index_to_leaf(ray_index)%link
           else
              call amr_regular_get_ixyz(ray_index,amrray_ix_curr,amrray_iy_curr,amrray_iz_curr)
           endif
        else
           nullify(amrray_cell)
           amrray_ix_curr = -1
           amrray_iy_curr = -1
           amrray_iz_curr = -1
        endif
     else
        write(stdo,*) 'ERROR: This grid type not yet implemented'
        stop
     endif
     !
  enddo
  !
end subroutine lambdaoper_integrate_ray

!------------------------------------------------------------------------
!                        FIND A RANDOM RAY
!------------------------------------------------------------------------
subroutine lambdaoper_draw_random_ray(doshift)
  implicit none
  integer :: icell
  logical :: doshift
  !
  ! We have different ways to choose random rays
  !
  select case(lambdaoper_random_ray_method)
  case(0)
     !
     ! Simplest method: Each cell has the same weight
     !
     ! Find cell
     !
     icell = floor(ran2(iseed)*nrcells)+1
     if(icell.gt.nrcells) stop 9231
     !
     ! Random location in cell
     !
     call montecarlo_random_pos_in_cell(icell)
     !
     ! Find the direction in which to emit
     !
     call montecarlo_randomdir(ray_cart_dirx,ray_cart_diry,ray_cart_dirz)
     !
!
! Note: When we use weighted cells (rays emanating from some cells are
!       more likely than from others), then we must include an extra
!       counter-weighting of the rays to prevent improper weighting in
!       the calculation of the mean intensity.
  case default
     write(stdo,*) 'ERROR: Do not know lambda random ray method ',lambdaoper_random_ray_method
     stop
  end select
  !
  ! If "doshift" is set, then we shift the starting point of the ray
  !
  if(doshift) then
     !
     ! Compute the base shift
     !
     if(grid_contsph_r.eq.0.d0) then
        write(stdo,*) 'ERROR: The grid_contsph_r is not set. Did you forget to'
        write(stdo,*) '       call the subroutine postprocess_grid_for_use() after'
        write(stdo,*) '       setting up the grid?'
        stop
     endif
     lambdaoper_shift = 2.01d0*grid_contsph_r
     !
     ! If star_maxexterndist>0
     !
     if(star_maxexterndist.gt.0.d0) then
        lambdaoper_shift = lambdaoper_shift + 2*star_maxexterndist
     endif
     !
     ! Now shift
     !
     ray_cart_x = ray_cart_x - lambdaoper_shift * ray_cart_dirx
     ray_cart_y = ray_cart_y - lambdaoper_shift * ray_cart_diry
     ray_cart_z = ray_cart_z - lambdaoper_shift * ray_cart_dirz
     !
     ! The end s-value is the shift
     !
     ray_dsend = lambdaoper_shift
  else
     ray_dsend = 1d99
  endif
end subroutine lambdaoper_draw_random_ray

!------------------------------------------------------------------------
!                      RANDOM RAY LAMBDA OPERATOR
!------------------------------------------------------------------------
subroutine lambdaoper_operator_cellbased_random_rays(nrays)
  implicit none
  integer :: iray,nrays
  !
  ! Check
  !
  if(.not.allocated(lambdaoper_meanint_cell)) then
     write(*,*) 'ERROR: Cannot use Random-Ray Lambda Operator without'
     write(*,*) '       first calling lambdaoper_init() with vertexbased=.false..'
     stop
  endif
  !
  ! Clear the mean intensity array
  !
  lambdaoper_meanint_cell(:) = 0.d0
  lambdaoper_dstot_cell(:) = 0.d0
  !
  ! Do a loop over all rays
  !
  do iray=1,nrays
     !
     ! Find a random ray: This sets the ray_cart_x,y,z and ray_cart_dirx,y,z
     ! according to some recipe
     !
     call lambdaoper_draw_random_ray(.true.)
     !
     ! Now integrate along the ray
     !
     call lambdaoper_integrate_ray()
  enddo
  !
  ! Compute the mean intensity
  !
  lambdaoper_meanint_cell(:) = lambdaoper_meanint_cell(:) / ( lambdaoper_dstot_cell(:) + 1d-99 )
  !
  ! Compute the diagonal elements of the Lambda Operator
  !
  lambdaoper_lambda_diag_cell(:) = lambdaoper_lambda_diag_cell(:) / ( lambdaoper_dstot_cell(:) + 1d-99 )
  !
end subroutine lambdaoper_operator_cellbased_random_rays

!========================================================================
!               LINE-TRANSFER LAMBDA OPERATOR STUFF
!========================================================================

end module lambdaoper_module

max_stars_repo_name: dlmatra/miao

max_issues_repo_name: dlmatra/miao

max_forks_repo_name: dlmatra/miao


text:
* ----------------- BEGIN gm_rv_drvr ---------------------------
      Program gm_rv_drvr

*  Obtains input parameters and calls random vibration simulation

* Notes regarding modification of driver:
*  1. All "include" statements assume all subprograms are in the same folder
*     as this program.
*  2. If values of stress other than those in the input parameter file are
*     to be used, set "stressc" equal to the new value, rather than "stress"
*     (and also set the other parameters if the stress changes with moment).

* Dates: 06/01/95 - Written by D.M. Boore; Modified from TD_DRVR
*        06/09/95 - Bells and whistles still being added
*        08/11/95 - Add frequency column to output
*        08/11/95 - Changed places where parameter-file and column-file names
*                   are requested
*        08/18/95 - Added call to Write_Params
*        10/17/95 - Added a flag to Get_Params and Write_Params to tell if using
*                   time domain or random vibration procedure
*        11/14/95 - Combined smsim.fi and rv.fi
*        11/16/95 - Minor changes to i/o; moved calculation of fup to get_params
*        12/06/95 - Added "patience" message on screen.
*        12/08/95 - Added switch for 91 standard periods.
*        12/14/95 - Print out frequency and duration of excitation
*        12/14/95 - Repeated computation of fup in this driver.
*        12/17/95 - Changed location of asking for col file and reordered
*                   some output
*        12/31/95 - Modified some of the *.sum output
*        01/03/96 - Added pk_rms_cl_eq68 to *.sum output
*        01/05/96 - Added column file stem name to column headings
*        02/06/96 - Minor formatting improvements
*        04/12/96 - changed names psv, psa to psv, psa
*        10/21/96 - added an override option for fup calculated in the driver
*        02/28/97 - added output of kappa (because I now allow it to be
*                   magnitude dependent)
*        04/20/97 - increased dimension of per from 91 to 400
*        01/14/99 - Added calculation of Arias intensity
*        01/22/99 - Modified output, including Sd, pga, pgv in column file
*        01/22/99 - changed names psv, psa back to psv, psa
*        02/09/99 - Added computation of displacement
*        02/10/99 - Include option to print pga, pgd, pgd and similarly
*                   pga/w, pga/w*w, etc in a column file to provide peak
*                   motion levels in plots of response spectra.  The file
*                   has so many columns that I decided to write the
*                   information to a separate file rather than include it
*                   in the response spectra file.
*        03/05/99 - Added include statements at end to bring in all needed
*                   subroutines, rather than including them in the
*                   cl_rv.bat batch file.  Also, I removed the include calls to
*                   rvtdsubs.for and recipes.for from smsim_rv.for.
*                   Write smsim version number
*        03/10/99 - Change from 'a\' to 'a' in format statements.  This is
*                   less pleasing aesthetically, but it allows
*                   for redirection of input/output, which did not
*                   work before (a test program showed that redirection
*                   works as long as the total of the sum of the characters
*                   in the query strings written to the screen do not exceed
*                   some undetermined number). Also, change "call banner(5)"
*                   to "call banner(6)" and remove the "cr to quit" in the
*                   "dist (cr to quit)" query
*        03/12/99 - Added check for nper = 1 in loop over computing log-spaced
*                   periods (this is "fool-proofing"; if only one period is
*                   desired the user should reply "y" to the individual
*                   periods query and the computations with nper-1 in the
*                   denominator should never be reached).
*        05/23/99 - Add an extra decimal to per output and add psa/pga (shape)
*                   to output
*        02/05/00 - Used trim_c for f_sum, f_col and corrected error in
*                   making column heads for f_sum stem less than 4 characters
*        06/13/02 - Renamed "rv_drvr" to "gmrvdrvr" and "smsim_rv" to "gm_rv"
*                   to indicate that various ground motion measures are being
*                   computed using random vibration theory.
*        02/09/03 - Only print stress variable if numsource = 1 or 2
*        07/16/05 - Suggest the name of output files
*        07/27/05 - Improve obtaining input (among other things, define "rply"
*                   with more characters, in case the user enters a space
*                   before the requested reply).
*        08/16/05 - Renamed from gmrvdrvr to gm_rv_drvr
*        05/11/07 - Removed "\smsim\" from include statements
*        08/19/08 - Change the construction of the suggested col file,
*                   because the stem name of the params file can contain periods
*                   before the last period.
*        02/26/09 - Updated calls to get date and time.
*        12/01/09 - Redo computation of fup, depending on the values of fmax and kappa
*        12/02/09 - Because the same code for computing fup appears in a number of drivers, I've
*                   consolidated it into a subroutine.
*        12/04/09 - Write calculated values of dursource and durpath to summary file.
*        12/20/09 - Use coefficients relating log10_H_eff to M to modify the distance.
*        12/22/09 - Add m, r, h_eff, rmod, numsource, stress to output
*        02/12/10 - As a result of feedback from John Douglas, I replaced "cr" and "psv"
*                   in the screen prompts with "Enter" and "response spectra", respectively.
!        08/25/10 - Shorten "patience" screen message.
!        09/30/10 - Use e format when writing psa/pga to avoid discretization effects in plots
!        01/15/11 - Add durex to output
!        01/17/11 - Add trms to output
!        01/21/11 - Add Der Kiureghain/Asfura estimates to summary file, as well
!                   as estimates using trms when calculation Nz and Ne.  This required
!                   changes to rv_subs and to smsim.fi
!        01/29/11 - Add amom0,1,2,4 to summary file
!        01/31/11 - Change nout to nu_sum, f_out to f_sum
!        02/06/11 - Explicitly set idva = 0 for an oscillator, even though this is also
!                   done in gm_rv in rv_subs (so this change is not needed--the output
!                   will be the same as before).  I am doing it here as well for clarity in
!                   reading the code, which has become quite convoluted.
!        04/05/11 - I replaced "h_eff" with "f_ff" (because the finite-fault factor is not
!                   really a pseudo-depth).
!                 - Call rmod_calc to compute rmod (rmod = rmod_calc(r, m, iflag, c1, c2, f_ff) )
!        05/23/11 - Print out avib (what Liu and Pezeshk call "k"; the value should be 1/3 for the
!                   BJ84 oscillator correction) to the summary file.
!        07/02/11 - Incorporate the Boore-Thompson modifications to the BJ84 oscillator correction
!                   for computing Trms.  The type of correction is specified in the params file, as follows:
!                   osc_crrctn(0=no correction;1=bj84;2=lp99; 3=bt12 WNA; 4=bt12 ENA; 5=bt12 avg. ena and wna)
!                 - Use generic names for pars files (ena250_bt12.pars and wna100_bt12.pars), so can
!                   change them without revising and recompiling the program.
!                 - Read in both ena and wna oscillator corrections.   The decision of whether to
!                   use one or the other or both is now made in rv_subs (as specified
!                   by osc_crrctn).
!        07/05/11 - Write smsim parameter file name in the first line of the output file.
!        08/02/11 - Changed "bt11" to "bt12", in anticipation of the paper being published in 2012.
!        08/08/11 - Change input order of BT12 coefficients
!        09/17/11 - The Trms4Osc coefficients are now read within get_params
!        12/02/11 - Add computation of am0 right after specification of amag (previously, it
!                   was contained in spect_scale, but it was not clear from where that routine was
!                   called).
!                 - I also corrected an error, in which I called rmod_calc, with amag as an
!                   argument, before specifying amag.  I corrected this by changing the order
!                   in which I ask for amag and r (amag first)
!        01/07/14 - Longer file names
!        10/13/14 - Add warning about using CL68 rather than DK80 R2PF
!                 - Allow for f_ff being defined by two lines
!        12/21/14 - Allow for f_ff being defined by a transition curve between two lines; this
!                   required a change to the calling arguments of rmod_calc
!        11/14/15 - Revised to use the simulations made using the DK80 rms-to-peak factors.
!                   The previous version, which used simulations made with the CL68 rms-to-peak
!                   factors, has been renamed gm_rv_drvr_cl

      character f_params*200, f_sum*200, buf*200, rply*10, f_pavd*200
      character f_suggest*120
      character colsd_head*12, colpsv_head*12, colpsa_head*12
      character col_shape*16
      character psvcalc*10, indvidper*10, f_col*200, message*80,
     :  standard_periods*10, write_pavd*10
      character date*8, time_start*10, time_stop*10
      logical tdflag, fparam_exist
      real t_pgd(2), t_pgv(2), t_pga(2)

      real*8 ane_d, xi_d, w

      REAL PER(400)
      DATA PER/0.040,0.042,0.044,0.046,0.048,
     * 0.050,0.055,0.060,0.065,0.070,
     * 0.075,0.080,0.085,0.090,0.095,0.10,0.11,0.12,0.13,0.14,
     * 0.15,0.16,0.17,0.18,0.19,0.20,0.22,0.24,0.26,0.28,
     * 0.30,0.32,0.34,0.36,0.38,0.40,0.42,0.44,0.46,0.48,
     * 0.50,0.55,0.60,0.65,0.70,0.75,0.80,0.85,0.90,0.95,
     * 1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,
     * 2.0,2.2,2.4,2.6,2.8,3.0,3.2,3.4,3.6,3.8,
     * 4.0,4.2,4.4,4.6,4.8,5.0,5.5,6.0,6.5,7.0,
     * 7.5,8.0,8.5,9.0,9.5,10.0,11.0,12.0,13.0,14.0,
     * 15.0,309*0.001/

      include 'bt12osc.fi'

      include 'smsim.fi'

      pi = 4.0 * atan(1.0)
      twopi = 2.0 * pi

      call banner(6)

      write(*, '(a)')
     :    ' Enter a message, if desired: '
      message = ' '
      read(*, '(a)') message

333   continue

      write(*,*)
      write(*,'(a)')
     :    ' ****** NOTE ******'
      write(*,'(a)')
     :    ' This program uses the Vanmarcke/Der Kiureghian'
      write(*,'(a)')
     :    ' rms-to-peak factor.'
      write(*,'(a)')
     :    ' ****** NOTE ******'
      write(*,*)

      f_params = ' '
      write(*, '(a)')
     :  ' Enter name of file with parameters (press "Enter" to quit): '
      read(*, '(a)') f_params
      if (f_params(1:4) == '    ') stop
      call trim_c(f_params,nc_f_params)
      inquire(file=f_params(1:nc_f_params), exist=fparam_exist)
      if (.not. fparam_exist) then
        write(*,'(a)') ' ******* FILE DOES NOT EXIST ******* '
        go to 333
      end if
      do i = nc_f_params, 1, -1
        if (f_params(i:i) == '.') then
           indexp = i
           go to 991
        end if
        indexp = 1
      end do
991   continue

      f_suggest = ' '
      f_suggest = f_params(1:indexp)//'gmrv.sum'
      call trim_c(f_suggest,nc_f_suggest)

      f_sum = ' '
      write(*, '(a)')
     :    ' Enter name of summary file ("Enter" = '//
     :     f_suggest(1:nc_f_suggest)//'):'
      read(*, '(a)') f_sum
      if (f_sum(1:4) == '    ') then
        f_sum = f_suggest(1:nc_f_suggest)
      end if
      call trim_c(f_sum,nc_f_sum)
      call get_lun(nu_sum)
      open(unit=nu_sum,file=f_sum(1:nc_f_sum),status='unknown')

      call banner(nu_sum)

      write(nu_sum, '(a)') message
      write(nu_sum, '(2a)') ' output file: ',
     :   f_sum(1:nc_f_sum)

      write(nu_sum, '(a)')
     :  ' *** Results computed using GM_RV_DRVR ***'

* Standard Fortran 90 intrinsic Subroutine DATE_AND_TIME
      call DATE_AND_TIME( date, time_start )
* Date is returned as 'CCYYMMDD' (character date*8)
* Time is returned as 'hhmmss.sss' (character time*10)
*     character datx*8, time_start*10
      write(nu_sum, *)
      write(nu_sum, '(a)')
     :   ' Date: '//date(1:4)//'/'//date(5:6)//'/'//date(7:8)
      write(nu_sum, '(a)')
     : ' Time Start: '//
     : time_start(1:2)//':'//time_start(3:4)//':'//time_start(5:10)
      write(nu_sum, *)

      write(nu_sum, '(2a)') ' file with parameters: ',
     :   f_params(1:20)

      write(*, '(a)') ' Compute response spectra? (y/n): '
      psvcalc = ' '
      read(*, '(a)') psvcalc
      call trim_c(psvcalc,nc_psvcalc)
      call upstr(psvcalc)
      if (psvcalc(1:1) == 'Y') then
        write(*, '(a)')
     :  ' Enter individual periods? (y/n): '
        indvidper = ' '
        read(*, '(a)') indvidper
        call trim_c(indvidper, nc_indvidper)
        call upstr(indvidper)
        if (indvidper(1:1) == 'Y') then
          write(*, '(a)')
     :       ' Enter fractional damping, nperiods: '
          read(*, *) damp, nper
          write(*, '(a)') '    Enter oscillator periods '//
     :        '(space or comma separated, single or multiple lines): '
          read(*,*) (per(i), i = 1, nper)
        else
          write(*, '(a)') ' Use standard set of 91 periods? (y/n): '
          standard_periods = ' '
          read(*, '(a)') standard_periods
          call trim_c(standard_periods, nc_standard_periods)
          call upstr(standard_periods)
          nper = 91
          if (standard_periods(1:1) == 'Y') then
            nper = 91
            write(*, '(a)') ' Enter fractional damping: '
            read(*, *) damp
          else
            write(*, '(a)')
     :       ' Enter fractional damping, perlow, perhigh, nper: '
            read(*,*) damp, perlow, perhigh, nper
            do i = 1, nper
              if (nper == 1) then
                per(i) = perlow
              else
                per(i) =
     :            perlow * (perhigh/perlow)**(float(i-1)/float(nper-1))
              end if
            end do
          end if
        end if
      else
        nper = 0
      end if

      tdflag = .false.  ! controls extent of input file read and written
      call get_params( f_params, tdflag )
      call write_params( nu_sum, tdflag)

! Compute fup.  Because the same code appears in a number of drivers, I've
! consolidated it into a subroutine.

      call get_fup(fm, akappa, amp_cutoff, nu_sum, fup)

100   continue

      write(*, '(a)') ' amag: '
      read(*, *) amag
      am0 = 10.**(1.5*amag + 16.05)

      write(*, '(a)') ' dist: '
      read(*, *) r

      rmod = rmod_calc(r, amag,
     :                  iflag_f_ff, nlines_f_ff,
     :                  c1_log10_f_ff, c2_log10_f_ff, mh_f_ff,
     :                  c3_log10_f_ff, c4_log10_f_ff,
     :                  m1t_f_ff, m2t_f_ff,
     :                  c0t_f_ff, c1t_f_ff, c2t_f_ff, c3t_f_ff,
     :                                                          f_ff)
      new_mr = .true.

      if (psvcalc(1:1) == 'Y') then
        f_suggest = ' '
        f_suggest = f_params(1:indexp-1)//'.mxxxxr00000_rs.rv.col'
        call trim_c(f_suggest,nc_f_suggest)
        write(f_suggest(indexp+2:indexp+5),'(f4.2)') amag
        if (rmod .lt. 10.0) then
          write(f_suggest(indexp+9:indexp+11),'(f3.1)') rmod
        else if (rmod .lt. 100.0) then
          write(f_suggest(indexp+8:indexp+11),'(f4.1)') rmod
        else if (rmod .lt. 1000.0) then
          write(f_suggest(indexp+7:indexp+11),'(f5.1)') rmod
        else if (rmod .lt. 10000.0) then
          write(f_suggest(indexp+8:indexp+11),'(i4)') int(rmod)
        else
          write(f_suggest(indexp+7:indexp+11),'(i5)') int(rmod)
        end if
        call trim_c(f_suggest,nc_f_suggest)

        f_col = ' '
        write(*, '(a)')
     :    ' Enter name of column file to contain'//
     :    ' response spectra ("Enter" = '//
     :     f_suggest(1:nc_f_suggest)//'):'
        read(*, '(a)') f_col
        call trim_c(f_col,nc_f_col)
        if (f_col(1:4) == '    ') then
          f_col = f_suggest(1:nc_f_suggest)
        end if
        call trim_c(f_col,nc_f_col)
        call get_lun(ncol)
        open(unit=ncol, file=f_col(1:nc_f_col), status='unknown')
        write(ncol,'(a)') ' Parameter file: '//f_params(1:nc_f_params)
        nc_f_col_stem = min0(8, index(f_col(1:nc_f_col),'.') - 1)
        colsd_head  =     ' sd:'//f_col(1:nc_f_col_stem)//'________'
        colpsv_head =     'psv:'//f_col(1:nc_f_col_stem)//'________'
        colpsa_head =     'psa:'//f_col(1:nc_f_col_stem)//'________'
        col_shape   = 'psa/pga:'//f_col(1:nc_f_col_stem)//'________'
        write(ncol, '(
     :    2x,a, 6x,a, 5x,a,
     :    6x,a, 7x,a, 4x,a, 4x,a,
     :    1x,a, 2x,a,
     :    1x,a, 1x,a, 1x,a, 1x,a,
     :    1x,a,
     :    2x,a, 4x,a, 6x,a, 7x,a
     :                            )')
     :   'damp', 'per', 'freq',
     :   'M', 'R', 'F_FF', 'Rmod',
     :   'numsource', 'stress',
     :    colsd_head, colpsv_head, colpsa_head, col_shape,
     :   'domfrqpsv',
     :   'dursource', 'durpath', 'durex', 'trms'

        write_pavd = ' '
        if (indvidper(1:1) /= 'Y') then
          write(*, '(a)')
     :      ' Write pga, pgv, pgd to a column file? (y/n): '
          read(*, '(a)') write_pavd
          call trim_c(write_pavd,nc_write_pavd)
          call upstr(write_pavd)
          if (write_pavd(1:1) == 'Y') then
            f_suggest= ' '
            f_suggest= f_params(1:indexp-1)//'.mxxxxr00000_pavd.rv.col'
            write(f_suggest(indexp+2:indexp+5),'(f4.2)') amag
            if (rmod .lt. 10.0) then
              write(f_suggest(indexp+9:indexp+11),'(f3.1)') rmod
            else if (rmod .lt. 100.0) then
              write(f_suggest(indexp+8:indexp+11),'(f4.1)') rmod
            else if (rmod .lt. 10000.0) then
              write(f_suggest(indexp+8:indexp+11),'(i4)') int(rmod)
            else
              write(f_suggest(indexp+7:indexp+11),'(i5)') int(rmod)
            end if
            call trim_c(f_suggest,nc_f_suggest)

            write(*, '(a)')
     :        ' Enter name of column file for to contain'//
     :        ' pga, pgv, pgd ("Enter" = '//
     :                       f_suggest(1:nc_f_suggest)//'):'
            f_pavd = ' '
            read( *, '(a)') f_pavd
            call trim_c(f_pavd, nc_f_pavd)
            if(f_pavd(1:4) == '    ') then
              f_pavd = f_suggest(1:nc_f_suggest)
            end if
            call trim_c(f_pavd,nc_f_pavd)

            call get_lun(npavd)
            open(unit=npavd, file=f_pavd(1:nc_f_pavd), status='unknown')
            tl = per(1)
            th = per(nper)
            tdiv = (th/tl)**(0.25)
            t_pga(1) = tl
            t_pga(2) = tl*tdiv
            t_pgv(1) = sqrt(th*tl/tdiv)
            t_pgv(2) = sqrt(th*tl*tdiv)
            t_pgd(1) = th/tdiv
            t_pgd(2) = th

            print *,' tl, th, tdiv, t_pgd(1,2)=',
     :                tl, th, tdiv, t_pgd(1),t_pgd(2)

            write(npavd,
     :            '(t4,a,    t17,a,    t26,a,    t36,a,
     :             t45,a,    t56,a,    t69,a,    t78,a,
     :             t86,a,    t96,a,   t108,a,   t121,a)')
     :           't_pgd',    'pgd',  'w*pgd', 'w2*pgd',
     :           't_pgv',  'pgv/w',    'pgv',  'w*pgv',
     :           't_pga', 'pga/w2',  'pga/w',    'pga'
          end if
        end if
      end if  ! end loop to compute response spectrum

      call DATE_AND_TIME( date, time_start )

      iaorins = 1

      idva = 0
      write(*,'(a)')
     : '  Computing peak displacement'
      call gm_rv(pgdsim)
      freq20pgd = freq20
      ane_pgd = ane
      anz_pgd = anz
      eps_pgd = eps_rv
      pk_rms_pgd = pk_rms_dk80_eq2
      amom0_pgd = amom0
      amom1_pgd = amom1
      amom2_pgd = amom2
      amom4_pgd = amom4

      idva = 1
      write(*,'(a)')
     : '  Computing peak velocity'
      call gm_rv(pgvsim)
      freq20pgv = freq20
      ane_pgv = ane
      anz_pgv = anz
      eps_pgv = eps_rv
      pk_rms_pgv = pk_rms_dk80_eq2
      amom0_pgv = amom0
      amom1_pgv = amom1
      amom2_pgv = amom2
      amom4_pgv = amom4

      idva = 2
      write(*,'(a)')
     : '  Computing peak acceleration'
      call gm_rv(pgasim)
      freq20pga = freq20
      ane_pga = ane
      anz_pga = anz
      eps_pga = eps_rv
      pk_rms_pga = pk_rms_dk80_eq2
      amom0_pga = amom0
      amom1_pga = amom1
      amom2_pga = amom2
      amom4_pga = amom4
      arias = arias_rv   ! must rename, because the variable arias_rv
                         ! is assigned a value for oscillator runs, but in
                         ! those cases it is not the true Arias intensity
                         ! (it is the zeroth moment of the oscillator response)

      write(nu_sum,*)
      write(nu_sum, '(a)') ' *********** NEW R AND M ********** '
      write(nu_sum, '(a, 1p5(1x,e10.3))')
     :      ' r, f_ff, rmod, amag, kappa = ',
     :        r, f_ff, rmod, amag, kappa_f(amag)
      if (psvcalc(1:1) == 'Y') then
        write(nu_sum, '(2x,2a)') 'File with response spectra: ',
     :                         f_col(1:30)
      else
        write(nu_sum, '(2x,a)') 'No response spectra computed'
      end if

      write(nu_sum, '(2x,a,1x,1pe10.3)') ' const= ', const

      if (numsource == 1 .or. numsource == 2) then
        write(nu_sum, '(2x,a, f6.3,
     :      1p, 1x,e9.2, 2(1x,e10.3), 3(1x,e9.2))')
     :    ' amag, stress, fa, fb, dursource, durpath, durex= ',
     :      amag, stress, fa, fb, dursource_calc, durpath_calc, durex
      else
        write(nu_sum, '(2x,a, f6.3,
     :      1p, 2(1x,e10.3), 3(1x,e9.2))')
     :    ' amag, fa, fb, dursource, durpath, durex= ',
     :      amag, fa, fb, dursource_calc, durpath_calc, durex
      end if
      write(nu_sum, '(2x,a,1p2(1x,e10.3))') ' am0, am0b_m0fa= ',
     :                                       am0, am0b_m0
      write(nu_sum,'( 30x,a,
     :                 2x,a, 4x,a,
     :                 6x,a, 6x,a, 1x,a,
     :                 8x,a, 8x,a, 8x,a, 8x,a
     :                                                     )')
     :             'pga(cm/s2)',
     :             'domfreq', 'eps',
     :             'nx', 'nz', 'pk_rms',
     :             'amom0', 'amom1', 'amom2', 'amom4'
      write(nu_sum,
     :      '( 30x,es10.2,
     :          2x,f7.2, 1x,f6.4,
     :          1x,f7.2, 1x,f7.2, 1x,f6.2,
     :        4(1x,es12.5) )')
     :      pgasim,
     :      freq20pga, eps_pga,
     :      ane_pga, anz_pga, pk_rms_pga,
     :      amom0_pga, amom1_pga, amom2_pga, amom4_pga

      write(nu_sum,'( 31x,a,
     :                 2x,a, 4x,a,
     :                 6x,a, 6x,a, 1x,a,
     :                 8x,a, 8x,a, 8x,a, 8x,a
     :                                                     )')
     :             'pgv(cm/s)',
     :             'domfreq', 'eps',
     :             'nx', 'nz', 'pk_rms',
     :             'amom0', 'amom1', 'amom2', 'amom4'
      write(nu_sum,
     :      '( 30x,es10.2,
     :          2x,f7.2, 1x,f6.4,
     :          1x,f7.2, 1x,f7.2, 1x,f6.2,
     :        4(1x,es12.5) )')
     :      pgvsim,
     :      freq20pgv, eps_pgv,
     :      ane_pgv, anz_pgv, pk_rms_pgv,
     :      amom0_pgv, amom1_pgv, amom2_pgv, amom4_pgv

      write(nu_sum,'( 33x,a,
     :                 2x,a, 4x,a,
     :                 6x,a, 6x,a, 1x,a,
     :                 8x,a, 8x,a, 8x,a, 8x,a
     :                                                     )')
     :             'pgd(cm)',
     :             'domfreq', 'eps',
     :             'nx', 'nz', 'pk_rms',
     :             'amom0', 'amom1', 'amom2', 'amom4'
      write(nu_sum,
     :      '( 30x,es10.2,
     :          2x,f7.2, 1x,f6.4,
     :          1x,f7.2, 1x,f7.2, 1x,f6.2,
     :        4(1x,es12.5) )')
     :      pgdsim,
     :      freq20pgd, eps_pgd,
     :      ane_pgd, anz_pgd, pk_rms_pgd,
     :      amom0_pgd, amom1_pgd, amom2_pgd, amom4_pgd

      write(nu_sum, '(t19,a)') 'Arias intensity (cm/s)'
      write(nu_sum, '(1p, t32, e9.2)') arias

      if (psvcalc(1:1) == 'Y') then
        write(nu_sum,'(2x,a,f6.3)')
     :        ' Fractional oscillator damping = ', damp
        write(nu_sum,'( 2x,a, 4x,a,
     :                5x,a, 2x,a, 1x,a,
     :               1x,a,
     :               4x,a, 5x,a, 5x,a, 1x,a, 1x,a,
     :               5x,a, 1x,a,
     :               1x,a,
     :               1x,a, 1x,a,
     :               1x,a, 1x,a, 1x,a,
     :               1x,a, 1x,a,
     :               1x,a, 1x,a,
     :               1x,a, 1x,a,
     :               2x,a, 1x,a,
     :               1x,a, 1x,a,
     :               4x,a, 1x,a,
     :               1x,a, 1x,a,
     :               2x,a, 1x,a,
     :               8x,a, 8x,a, 8x,a, 8x,a,
     :               7x,a,
     :               2x,a, 4x,a, 6x,a, 7x,a,
     :               4x,a,
     :               1x,a, 1x,a
     :                                              )')
     :             'per(s)', 'freq',
     :             'sd(cm)', 'psv(cm/s)', 'psa(cm/s2)',
     :             'domfreq',
     :             'eps', 'nx', 'nz', 'pk_rms', 'zupusd',
     :      'deltay',         'deltay_prime',
     :      'anz_prime',
     :      'ane_psv',        'ane_dur4n_trms',
     :      'anz_psv',        'anz_prime', 'anz_dur4n_trms',
     :      'pk_rms_cl_eq68', 'pk_rms_cl_eq68_dur4n_trms',
     :      'pk_rms_dk_2',    'pk_rms_dk_dur4n_trms',
     :      'pk_cl_eq68/w',   'pk_cl_eq68_dur4n_trms/w',
     :      'pk_dk_2/w',      'pk_dk_dur4n_trms/w',
     :      'pk_cl_eq68',     'pk_cl_eq68_dur4n_trms',
     :      'pk_dk_2',        'pk_dk_dur4n_trms',
     :      'pk_cl_eq68*w',   'pk_cl_eq68_dur4n_trms*w',
     :      'pk_dk_2*w',      'pk_dk_dur4n_trms*w',
     :      'amom0', 'amom1', 'amom2', 'amom4',
     :      'avib',
     :      'dursource', 'durpath', 'durex', 'trms',
     :      'per/durex',
     :      'sqrt(durex/trms)', 'sqrt(trms/durex)'

        iaorins = 2
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
! Assume that instrument response is relative to ground
! displacement.  In this case idva = 0, and to make sure
! that this is so, I include the following statement
! (this is also set in gm_rv in rv_subs; I am repeating here
! for clarity when reading this code later):

        idva = 0
!
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

        do i = 1, nper

          write(*,'(a,i3,a,i3,a)')
     :    '+ Computing response spectra number ',
     :     i, ' out of ', nper, ' periods.'

          perosc = per(i)
          call gm_rv(psvsim)
          freq20psv = freq20
          ane_psv = ane
          anz_psv = anz
          eps_psv = eps_rv
          pk_rms_psv = pk_rms_dk80_eq2

          w = twopi/perosc

          write(nu_sum,
     :      '( t2,f7.3, 1x,f7.3,
     :         1x,es10.3, 1x,es10.3, 1x,es10.3,
     :         1x,f7.2,
     :         1x,f6.4, 1x,f6.1, 1x,f6.1, 1x,f6.2, 2x,f5.2,
     :         1x,es10.3, 3x,es10.3,
     :         4x,f6.1,
     :         2x,f6.1, 9x,f6.1,
     :         2x,f6.1, 4x,f6.1, 9x,f6.1,
     :         9x,f6.2, 20x,f6.2,
     :         6x,f6.2, 15x,f6.2,
     :         3x,es10.3, 14x,es10.3,
     :         1x,es10.3, 9x,es10.3,
     :         1x,es10.3, 12x,es10.3,
     :         1x,es10.3, 7x,es10.3,
     :         3x,es10.3, 14x,es10.3,
     :         1x,es10.3, 9x,es10.3,
     :         1x,es12.5, 1x,es12.5, 1x,es12.5, 1x,es12.5,
     :         1x,es10.3,
     :         4(1x,es10.3),
     :         1x,es12.5,
     :         7x,es10.3, 7x,es10.3
     :                                                        )')
     :      perosc, 1.0/perosc,
     :      psvsim/w, psvsim, w*psvsim,
     :      freq20psv,
     :      eps_psv, ane_psv, anz_psv, pk_rms_psv, zup_used,
     :      deltay, deltay_prime,
     :      anz_prime,
     :      ane_psv, ane_dur4n_trms,
     :      anz_psv, anz_prime, anz_dur4n_trms,
     :      pk_rms_cl_eq68, pk_rms_cl_eq68_dur4n_trms,
     :      pk_rms_dk_2,   pk_rms_dk_2_dur4n_trms,
     :      pk_cl_eq68/w, pk_cl_eq68_dur4n_trms/w,
     :      pk_dk_2/w, pk_dk_2_dur4n_trms/w,
     :      pk_cl_eq68, pk_cl_eq68_dur4n_trms,
     :      pk_dk_2, pk_dk_2_dur4n_trms,
     :      pk_cl_eq68*w, pk_cl_eq68_dur4n_trms*w,
     :      pk_dk_2*w, pk_dk_2_dur4n_trms*w,
     :      amom0, amom1, amom2, amom4,
     :      avib,
     :      dursource_calc, durpath_calc, durex, trms,
     :      perosc/durex,
     :      sqrt(durex/trms), sqrt(trms/durex)

          if (numsource == 1 .or. numsource == 2 .or.
     :        numsource == 11) then
            stress_out = stress
          else
            stress_out = 0.0
          end if

          write(ncol,'(
     :      1x,f5.3, 1x,f8.4, 1x,f8.3,
     :      1x,f6.3, 1x,f7.2, 1x,f7.2, 1x,f7.2,
     :      8x,i2, 1x,f7.2,
     :      3x,es10.3, 3x,es10.3, 3x,es10.3,
     :      7x,es10.3,
     :      1x,f9.4,
     :      4(1x,es10.3)
     :                                 )')
     :      damp, perosc, 1.0/perosc,
     :      amag, r, f_ff, rmod,
     :      numsource, stress_out,
     :      psvsim/(twopi/perosc), psvsim, (twopi/perosc)*psvsim,
     :      (twopi/perosc)*psvsim/pgasim,
     :      freq20psv,
     :      dursource_calc, durpath_calc, durex, trms

        end do
        close(unit=ncol)
        if (write_pavd(1:1) == 'Y') then
          do i = 1, 2
            write(npavd, '(3(1x,f7.3,3(1x,es10.3)))')
     :       t_pgd(i),pgdsim,pgdsim*(twopi/t_pgd(i)),
     :                pgdsim*(twopi/t_pgd(i))**2,
     :       t_pgv(i),pgvsim/(twopi/t_pgv(i)),pgvsim,
     :                pgvsim*(twopi/t_pgv(i)),
     :       t_pga(i),pgasim/(twopi/t_pga(i))**2,
     :                pgasim/(twopi/t_pga(i)),pgasim
          end do
          close(unit=npavd)
        end if
      end if

      write(nu_sum, *)
      call DATE_AND_TIME( date, time_stop )
      call time_diff(time_start, time_stop, time_elapsed)
      write(nu_sum, '(a,1x,1pe10.3)')
     :     ' Elapsed time (sec): ', time_elapsed

      write(*, '(a)')
     : ' Compute results for another r and M? (y/n;"Enter"=quit): '
      read(*, '(a)') buf
      if (buf(1:4) == '    ') go to 999
      if (buf(1:1) == 'n' .or. buf(1:1) == 'N') go to 999

      goto 100

999   continue

      close(unit=nu_sum)
      stop
      end
* ----------------- END gm_rv_drvr ---------------------------

!      include 'rv_subs.for'
!      include 'rv_td_subs.for'
!      include 'recipes.for'
      include 'smsim_subs_for_rv_programs.for'
      include 'smsim_util_subs.for'

ext: for

max_stars_repo_name: cliftbar/FlaPyDisaster

max_issues_repo_name: cliftbar/FlaPyDisaster

max_forks_repo_name: cliftbar/FlaPyDisaster


!==============================================================================!
  subroutine Backup_Mod_Write_Variable(fh, disp, vc, var_name, var)
!------------------------------------------------------------------------------!
!   Writes a whole variable to backup file.                                    !
!------------------------------------------------------------------------------!
  implicit none
!---------------------------------[Arguments]----------------------------------!
  integer          :: fh, disp, vc
  character(len=*) :: var_name
  type(Var_Type)   :: var
!-----------------------------------[Locals]-----------------------------------!
  type(Comm_Type), pointer :: comm
  character(len=80)        :: vn
  integer                  :: vs  ! variable size
!==============================================================================!

  ! Take alias
  comm => var % pnt_grid % comm

  if(this_proc < 2) print *, '# Writing variable: ', trim(var_name)

  ! Increase variable count
  vc = vc + 1

  ! Vector without boundaries
  vn = var_name
  call Comm_Mod_Write_Text(fh, vn, disp)
  vs = (2*comm % nc_t + 2*comm % nb_t) * SIZE_REAL
  call Comm_Mod_Write_Int (fh, vs, disp)

  call Comm_Mod_Write_Cell_Real(comm, fh, var % n(1:comm % nc_s),  disp)
  call Comm_Mod_Write_Bnd_Real (comm, fh, var % n(-comm % nb_f:  &
                                                  -comm % nb_l),   disp)
  call Comm_Mod_Write_Bnd_Real (comm, fh, var % q(-comm % nb_f:  &
                                                  -comm % nb_l),   disp)
  call Comm_Mod_Write_Cell_Real(comm, fh, var % o(1:comm % nc_s),  disp)

  end subroutine