Module ks_diag_tests
  Use test_params
  Use mpi
  Implicit None
Contains

  Subroutine test_ks_array_diag()

    Use numbers_module , Only : wp
    Use ks_array_module, Only : ks_array, ks_array_init, ks_array_comm_to_base, ks_array_finalise, &
         K_POINT_REAL, K_POINT_COMPLEX, ks_array_replicated_1D
!!$    Use mpi            , Only : mpi_bcast, mpi_comm_world, mpi_double_complex, mpi_double_precision
    Use mpi            , Only : mpi_comm_world, mpi_double_complex, mpi_double_precision

    Integer, Parameter :: ns = 2
    Integer, Parameter :: nk = 3

    Type( ks_array ) :: A, A_split
    Type( ks_array ) :: Q
    Type( ks_array ) :: QT
    Type( ks_array ) :: B

    Type( ks_array_replicated_1D ), Dimension( 1:nk * ns ) :: E

    Type( ks_array ) :: base_k

    Complex( wp ), Dimension( :, :, :, : ), Allocatable :: A_c
    Complex( wp ), Dimension( :, :    ), Allocatable :: tmp_c

    Real( wp ), Dimension( :, :, :, : ), Allocatable :: A_r
    Real( wp ), Dimension( :, :    ), Allocatable :: tmp_r

    Real( wp ), Dimension( : ), Allocatable :: evals_ks

    Real( wp ) :: rand
    Real( wp ) :: max_diff

    Integer, Dimension( 1:3, 1:nk ) :: k_points
    Integer, Dimension(      1:nk ) :: k_types

    Integer :: k, s
    Integer :: i
    Integer :: n
    Integer :: error

    n = m

    Allocate( A_r( 1:n, 1:n, 1:nk, 1:ns ) )
    Allocate( A_c( 1:n, 1:n, 1:nk, 1:ns ) )
    Allocate( tmp_r( 1:n, 1:n ) )
    Allocate( tmp_c( 1:n, 1:n ) )

    A_r = Huge( A_r )
    A_c = Huge( Real( A_c, Kind( A_c ) ) )

    If( me == 0 ) Then

       Call Random_Number( A_r )
       Do s = 1, ns
          Do k = 1, nk
             A_r( :, :, k, s ) = A_r( :, :, k, s ) + Transpose( A_r( :, :, k, s  ) )
          End Do
       End Do

       Do s = 1, ns
          Do k = 1, nk
             Call Random_Number( tmp_r )
             A_c( :, :, k, s ) = tmp_r
             Call Random_Number( tmp_r )
             A_c( :, :, k, s ) = A_c( :, :, k, s ) + Cmplx( 0.0_wp, tmp_r, Kind = wp )
             A_c( :, :, k, s ) = A_c( :, :, k, s ) + Conjg( Transpose( A_c( :, :, k, s  ) ) )
          End Do
       End Do

       Do k = 1, nk
          k_points( :, k ) = [ k - 1, 0, 0 ]
          Call Random_Number( rand )
          k_types( k ) = Merge( K_POINT_REAL, K_POINT_COMPLEX, rand > 0.5_wp )
       End Do

    End If

    Call mpi_bcast( k_points, Size( k_points ), mpi_integer, 0, mpi_comm_world, error )
    Call mpi_bcast( k_types , Size( k_types  ), mpi_integer, 0, mpi_comm_world, error )

    Call mpi_bcast( A_r, Size( A_r ), mpi_double_precision, 0, mpi_comm_world, error )

    Call mpi_bcast( A_c, Size( A_c ), mpi_double_complex  , 0, mpi_comm_world, error )

    Call ks_array_init( n_block )
    Call ks_array_comm_to_base( MPI_COMM_WORLD, ns, k_types, k_points, base_k )

    Call A%create( n, n, base_k )
    If( verbose ) Then
       Call A%print_info( 'A', 200 )
    End If
    Do s = 1, ns
       Do k = 1, nk
          If( k_types( k ) == K_POINT_REAL ) Then
             Call A%set_by_global( k_points( :, k ), s, 1, n, 1, n, A_r( :, :, k, s ) )
          Else
             Call A%set_by_global( k_points( :, k ), s, 1, n, 1, n, A_c( :, :, k, s ) )
          End If
       End Do
    End Do

    Call A%split_ks( 2.0_wp, A_split )
    If( verbose ) Then
       Call A_split%print_info( 'A_split', 100 )
    End If

    Call A_split%diag( Q, E )
    QT = .Dagger. Q
    B = QT * A_split * Q
    If( verbose ) Then
       Call B%print_info( 'B', 100 )
    End If

    max_diff = -1.0_wp
    Do s = 1, ns
       Do k = 1, nk
          tmp_r = 0.0_wp
          evals_ks = E( k + ( s - 1 ) * nk )
          If( k_types( k ) == K_POINT_REAL ) Then
             Call B%get_by_global( k_points( :, k ), s, 1, n, 1, n, tmp_r )
             Do i = 1, n
                tmp_r( i, i ) = Abs( tmp_r( i, i ) - evals_ks( i ) )
             End Do
          Else
             Call B%get_by_global( k_points( :, k ), s, 1, n, 1, n, tmp_c )
             Do i = 1, n
                tmp_r( i, i ) = Abs( tmp_c( i, i ) - evals_ks( i ) )
             End Do
          End If
          max_diff = Max( max_diff, Maxval( tmp_r ) )
       End Do
    End Do
    If( me == 0 ) Then
       Write( *, error_format ) 'Error in ks_split diag ', max_diff, &
            Merge( passed, FAILED, max_diff < tol )
    End If

    Call ks_array_finalise

  End Subroutine test_ks_array_diag

  Subroutine test_ks_array_choleski()

    Use numbers_module , Only : wp
    Use ks_array_module, Only : ks_array, ks_array_init, ks_array_comm_to_base, ks_array_finalise, &
         K_POINT_REAL, K_POINT_COMPLEX
!!$    Use mpi            , Only : mpi_bcast, mpi_comm_world, mpi_double_complex, mpi_double_precision
    Use mpi            , Only : mpi_comm_world, mpi_double_complex, mpi_double_precision

    Integer, Parameter :: ns = 2
    Integer, Parameter :: nk = 3

    Type( ks_array ) :: A, A_split
    Type( ks_array ) :: L
    Type( ks_array ) :: LT
    Type( ks_array ) :: B

    Type( ks_array ) :: base_k

    Complex( wp ), Dimension( :, :, :, : ), Allocatable :: A_c
    Complex( wp ), Dimension( :, :    ), Allocatable :: tmp_c

    Real( wp ), Dimension( :, :, :, : ), Allocatable :: A_r
    Real( wp ), Dimension( :, :    ), Allocatable :: tmp_r

    Real( wp ) :: rand
    Real( wp ) :: max_diff

    Integer, Dimension( 1:3, 1:nk ) :: k_points
    Integer, Dimension(      1:nk ) :: k_types

    Integer :: k, s
    Integer :: n
    Integer :: error

    n = m

    Allocate( A_r( 1:n, 1:n, 1:nk, 1:ns ) )
    Allocate( A_c( 1:n, 1:n, 1:nk, 1:ns ) )
    Allocate( tmp_r( 1:n, 1:n ) )
    Allocate( tmp_c( 1:n, 1:n ) )

    A_r = Huge( A_r )
    A_c = Huge( Real( A_c, Kind( A_c ) ) )

    If( me == 0 ) Then

       Call Random_Number( A_r )
       Do s = 1, ns
          Do k = 1, nk
             ! Make sure matrix is positive definite
             A_r( :, :, k, s ) = Matmul( A_r( :, :, k, s ), Transpose( A_r( :, :, k, s  ) ) )
          End Do
       End Do

       Do s = 1, ns
          Do k = 1, nk
             Call Random_Number( tmp_r )
             A_c( :, :, k, s ) = tmp_r
             Call Random_Number( tmp_r )
             A_c( :, :, k, s ) = A_c( :, :, k, s ) + Cmplx( 0.0_wp, tmp_r, Kind = wp )
             ! Make sure matrix is positive definite
             A_c( :, :, k, s ) = Matmul( A_c( :, :, k, s ), Conjg( Transpose( A_c( :, :, k, s  ) ) ) )
          End Do
       End Do

       Do k = 1, nk
          k_points( :, k ) = [ k - 1, 0, 0 ]
          Call Random_Number( rand )
          k_types( k ) = Merge( K_POINT_REAL, K_POINT_COMPLEX, rand > 0.5_wp )
       End Do

    End If

    Call mpi_bcast( k_points, Size( k_points ), mpi_integer, 0, mpi_comm_world, error )
    Call mpi_bcast( k_types , Size( k_types  ), mpi_integer, 0, mpi_comm_world, error )

    Call mpi_bcast( A_r, Size( A_r ), mpi_double_precision, 0, mpi_comm_world, error )

    Call mpi_bcast( A_c, Size( A_c ), mpi_double_complex  , 0, mpi_comm_world, error )

    Call ks_array_init( n_block )
    Call ks_array_comm_to_base( MPI_COMM_WORLD, ns, k_types, k_points, base_k )

    Call A%create( n, n, base_k )
    If( verbose ) Then
       Call A%print_info( 'A', 200 )
    End If
    Do s = 1, ns
       Do k = 1, nk
          If( k_types( k ) == K_POINT_REAL ) Then
             Call A%set_by_global( k_points( :, k ), s, 1, n, 1, n, A_r( :, :, k, s ) )
          Else
             Call A%set_by_global( k_points( :, k ), s, 1, n, 1, n, A_c( :, :, k, s ) )
          End If
       End Do
    End Do

    Call A%split_ks( 2.0_wp, A_split )
    If( verbose ) Then
       Call A_split%print_info( 'A_split', 100 )
    End If

    L = .Choleski. A_split
    LT = .Dagger. L
    B = L * LT - A_split
    If( verbose ) Then
       Call B%print_info( 'B', 100 )
    End If

    max_diff = -1.0_wp
    Do s = 1, ns
       Do k = 1, nk
          tmp_r = 0.0_wp
          If( k_types( k ) == K_POINT_REAL ) Then
             Call B%get_by_global( k_points( :, k ), s, 1, n, 1, n, tmp_r )
             tmp_r = Abs( tmp_r )
          Else
             Call B%get_by_global( k_points( :, k ), s, 1, n, 1, n, tmp_c )
             tmp_r = Abs( tmp_c )
          End If
          max_diff = Max( max_diff, Maxval( tmp_r ) )
       End Do
    End Do
    If( me == 0 ) Then
       Write( *, error_format ) 'Error in ks_split choleski ', max_diff, &
            Merge( passed, FAILED, max_diff < tol )
    End If

    Call ks_array_finalise

  End Subroutine test_ks_array_choleski

  Subroutine test_ks_array_tr_inv()

    Use numbers_module , Only : wp
    Use ks_array_module, Only : ks_array, ks_array_init, ks_array_comm_to_base, ks_array_finalise, &
         K_POINT_REAL, K_POINT_COMPLEX
!!$    Use mpi            , Only : mpi_bcast, mpi_comm_world, mpi_double_complex, mpi_double_precision
    Use mpi            , Only : mpi_comm_world, mpi_double_complex, mpi_double_precision

    Integer, Parameter :: ns = 2
    Integer, Parameter :: nk = 3

    Type( ks_array ) :: A, A_split
    Type( ks_array ) :: L
    Type( ks_array ) :: L_inv
    Type( ks_array ) :: B

    Type( ks_array ) :: base_k

    Complex( wp ), Dimension( :, :, :, : ), Allocatable :: A_c
    Complex( wp ), Dimension( :, :    ), Allocatable :: tmp_c

    Real( wp ), Dimension( :, :, :, : ), Allocatable :: A_r
    Real( wp ), Dimension( :, :    ), Allocatable :: tmp_r

    Real( wp ) :: rand
    Real( wp ) :: max_diff

    Integer, Dimension( 1:3, 1:nk ) :: k_points
    Integer, Dimension(      1:nk ) :: k_types

    Integer :: k, s
    Integer :: i
    Integer :: n
    Integer :: error

    n = m

    Allocate( A_r( 1:n, 1:n, 1:nk, 1:ns ) )
    Allocate( A_c( 1:n, 1:n, 1:nk, 1:ns ) )
    Allocate( tmp_r( 1:n, 1:n ) )
    Allocate( tmp_c( 1:n, 1:n ) )

    A_r = Huge( A_r )
    A_c = Huge( Real( A_c, Kind( A_c ) ) )

    If( me == 0 ) Then

       Call Random_Number( A_r )
       Do s = 1, ns
          Do k = 1, nk
             ! Make sure matrix is positive definite
             A_r( :, :, k, s ) = Matmul( A_r( :, :, k, s ), Transpose( A_r( :, :, k, s  ) ) )
          End Do
       End Do

       Do s = 1, ns
          Do k = 1, nk
             Call Random_Number( tmp_r )
             A_c( :, :, k, s ) = tmp_r
             Call Random_Number( tmp_r )
             A_c( :, :, k, s ) = A_c( :, :, k, s ) + Cmplx( 0.0_wp, tmp_r, Kind = wp )
             ! Make sure matrix is positive definite
             A_c( :, :, k, s ) = Matmul( A_c( :, :, k, s ), Conjg( Transpose( A_c( :, :, k, s  ) ) ) )
          End Do
       End Do

       Do k = 1, nk
          k_points( :, k ) = [ k - 1, 0, 0 ]
          Call Random_Number( rand )
          k_types( k ) = Merge( K_POINT_REAL, K_POINT_COMPLEX, rand > 0.5_wp )
       End Do

    End If

    Call mpi_bcast( k_points, Size( k_points ), mpi_integer, 0, mpi_comm_world, error )
    Call mpi_bcast( k_types , Size( k_types  ), mpi_integer, 0, mpi_comm_world, error )

    Call mpi_bcast( A_r, Size( A_r ), mpi_double_precision, 0, mpi_comm_world, error )

    Call mpi_bcast( A_c, Size( A_c ), mpi_double_complex  , 0, mpi_comm_world, error )

    Call ks_array_init( n_block )
    Call ks_array_comm_to_base( MPI_COMM_WORLD, ns, k_types, k_points, base_k )

    Call A%create( n, n, base_k )
    If( verbose ) Then
       Call A%print_info( 'A', 200 )
    End If
    Do s = 1, ns
       Do k = 1, nk
          If( k_types( k ) == K_POINT_REAL ) Then
             Call A%set_by_global( k_points( :, k ), s, 1, n, 1, n, A_r( :, :, k, s ) )
          Else
             Call A%set_by_global( k_points( :, k ), s, 1, n, 1, n, A_c( :, :, k, s ) )
          End If
       End Do
    End Do

    Call A%split_ks( 2.0_wp, A_split )
    If( verbose ) Then
       Call A_split%print_info( 'A_split', 100 )
    End If

    ! Generate a traingular array by Choleski Decomposition
    L = .Choleski. A_split
    L_inv = .TrInv. L
    B = L * L_inv
    If( verbose ) Then
       Call B%print_info( 'B', 100 )
    End If

    max_diff = -1.0_wp
    Do s = 1, ns
       Do k = 1, nk
          tmp_r = 0.0_wp
          If( k_types( k ) == K_POINT_REAL ) Then
             ! B should be the unit matrix
             Call B%get_by_global( k_points( :, k ), s, 1, n, 1, n, tmp_r )
             tmp_r = Abs( tmp_r )
          Else
             Call B%get_by_global( k_points( :, k ), s, 1, n, 1, n, tmp_c )
             tmp_r = Abs( tmp_c )
          End If
          Do i = 1, n
             tmp_r( i, i ) = tmp_r( i, i ) - 1.0_wp
          End Do
          max_diff = Max( max_diff, Maxval( tmp_r ) )
       End Do
    End Do
    If( me == 0 ) Then
       Write( *, error_format ) 'Error in ks_split tr_inv ', max_diff, &
            Merge( passed, FAILED, max_diff < tol )
    End If

    Call ks_array_finalise

  End Subroutine test_ks_array_tr_inv

  Subroutine test_ks_array_tr_inv_with_iterator()

    Use numbers_module , Only : wp
    Use ks_array_module, Only : ks_array, ks_point_info, ks_array_init, ks_array_comm_to_base, ks_array_finalise, &
         K_POINT_REAL, K_POINT_COMPLEX, K_POINT_NOT_EXIST
!!$    Use mpi            , Only : mpi_bcast, mpi_comm_world, mpi_double_complex, mpi_double_precision
    Use mpi            , Only : mpi_comm_world, mpi_double_complex, mpi_double_precision

    Integer, Parameter :: ns = 2
    Integer, Parameter :: nk = 3

    Type( ks_array ) :: A, A_split
    Type( ks_array ) :: L
    Type( ks_array ) :: L_inv
    Type( ks_array ) :: B

    Type( ks_array ) :: base_k

    Type( ks_point_info ), Allocatable :: info

    Complex( wp ), Dimension( :, :, :, : ), Allocatable :: A_c
    Complex( wp ), Dimension( :, :    ), Allocatable :: tmp_c

    Real( wp ), Dimension( :, :, :, : ), Allocatable :: A_r
    Real( wp ), Dimension( :, :    ), Allocatable :: tmp_r

    Real( wp ) :: rand
    Real( wp ) :: max_diff

    Integer, Dimension( 1:3, 1:nk ) :: k_points
    Integer, Dimension(      1:nk ) :: k_types

    Integer :: k, s
    Integer :: i
    Integer :: n
    Integer :: error

    n = m

    Allocate( A_r( 1:n, 1:n, 1:nk, 1:ns ) )
    Allocate( A_c( 1:n, 1:n, 1:nk, 1:ns ) )
    Allocate( tmp_r( 1:n, 1:n ) )
    Allocate( tmp_c( 1:n, 1:n ) )

    A_r = Huge( A_r )
    A_c = Huge( Real( A_c, Kind( A_c ) ) )

    If( me == 0 ) Then

       Call Random_Number( A_r )
       Do s = 1, ns
          Do k = 1, nk
             ! Make sure matrix is positive definite
             A_r( :, :, k, s ) = Matmul( A_r( :, :, k, s ), Transpose( A_r( :, :, k, s  ) ) )
          End Do
       End Do

       Do s = 1, ns
          Do k = 1, nk
             Call Random_Number( tmp_r )
             A_c( :, :, k, s ) = tmp_r
             Call Random_Number( tmp_r )
             A_c( :, :, k, s ) = A_c( :, :, k, s ) + Cmplx( 0.0_wp, tmp_r, Kind = wp )
             ! Make sure matrix is positive definite
             A_c( :, :, k, s ) = Matmul( A_c( :, :, k, s ), Conjg( Transpose( A_c( :, :, k, s  ) ) ) )
          End Do
       End Do

       Do k = 1, nk
          k_points( :, k ) = [ k - 1, 0, 0 ]
          Call Random_Number( rand )
          k_types( k ) = Merge( K_POINT_REAL, K_POINT_COMPLEX, rand > 0.5_wp )
       End Do

    End If

    Call mpi_bcast( k_points, Size( k_points ), mpi_integer, 0, mpi_comm_world, error )
    Call mpi_bcast( k_types , Size( k_types  ), mpi_integer, 0, mpi_comm_world, error )

    Call mpi_bcast( A_r, Size( A_r ), mpi_double_precision, 0, mpi_comm_world, error )

    Call mpi_bcast( A_c, Size( A_c ), mpi_double_complex  , 0, mpi_comm_world, error )

    Call ks_array_init( n_block )
    Call ks_array_comm_to_base( MPI_COMM_WORLD, ns, k_types, k_points, base_k )

    Call A%create( n, n, base_k )
    If( verbose ) Then
       Call A%print_info( 'A', 200 )
    End If
!!$    Do s = 1, ns
!!$       Do k = 1, nk
!!$          If( k_types( k ) == K_POINT_REAL ) Then
!!$             Call A%set_by_global( k_points( :, k ), s, 1, n, 1, n, A_r( :, :, k, s ) )
!!$          Else
!!$             Call A%set_by_global( k_points( :, k ), s, 1, n, 1, n, A_c( :, :, k, s ) )
!!$          End If
!!$       End Do
!!$    End Do
    Call A%split_ks( 2.0_wp, A_split )
    If( verbose ) Then
       Call A_split%print_info( 'A_split', 100 )
    End If

    Call A_split%iterator_init()
    Do
       info = A_split%iterator_next()
!!$       Write( *, * ) info%k_type
       If( info%k_type == K_POINT_NOT_EXIST ) Exit
       Outer: Do s = 1, ns
          Do k = 1, nk
             If( All( info%k_indices == k_points( :, k ) .And. info%spin == s ) ) Then
                If( info%k_type == K_POINT_REAL ) Then
                   Call A_split%set_by_global( k_points( :, k ), s, 1, n, 1, n, A_r( :, :, k, s ) )
                Else
                   Call A_split%set_by_global( k_points( :, k ), s, 1, n, 1, n, A_c( :, :, k, s ) )
                End If
                Exit Outer
             End If
          End Do
       End Do Outer
    End Do
    Call A_split%iterator_reset()

    ! Generate a traingular array by Choleski Decomposition
    L = .Choleski. A_split
    L_inv = .TrInv. L
    B = L * L_inv
    If( verbose ) Then
       Call B%print_info( 'B', 100 )
    End If

    max_diff = -1.0_wp
    Do s = 1, ns
       Do k = 1, nk
          tmp_r = 0.0_wp
          If( k_types( k ) == K_POINT_REAL ) Then
             ! B should be the unit matrix
             Call B%get_by_global( k_points( :, k ), s, 1, n, 1, n, tmp_r )
             tmp_r = Abs( tmp_r )
          Else
             Call B%get_by_global( k_points( :, k ), s, 1, n, 1, n, tmp_c )
             tmp_r = Abs( tmp_c )
          End If
          Do i = 1, n
             tmp_r( i, i ) = tmp_r( i, i ) - 1.0_wp
          End Do
          max_diff = Max( max_diff, Maxval( tmp_r ) )
       End Do
    End Do
    If( me == 0 ) Then
       Write( *, error_format ) 'Error in ks_split tr_inv iterator ', max_diff, &
            Merge( passed, FAILED, max_diff < tol )
    End If

    Call ks_array_finalise

  End Subroutine test_ks_array_tr_inv_with_iterator

  Subroutine test_ks_array_tr_inv_with_iterator_and_extract()

    Use numbers_module , Only : wp
    Use ks_array_module, Only : ks_array, ks_point_info, ks_array_init, ks_array_comm_to_base, ks_array_finalise, &
         K_POINT_REAL, K_POINT_COMPLEX, K_POINT_NOT_EXIST
    Use mpi            , Only : mpi_comm_world, mpi_double_complex, mpi_double_precision

    Integer, Parameter :: ns = 2
    Integer, Parameter :: nk = 3

    Type( ks_array ) :: A, A_split
    Type( ks_array ) :: A_split_element
    Type( ks_array ) :: L
    Type( ks_array ) :: L_inv
    Type( ks_array ) :: B

    Type( ks_array ) :: base_k

    Type( ks_point_info ), Allocatable :: info

    Complex( wp ), Dimension( :, :, :, : ), Allocatable :: A_c
    Complex( wp ), Dimension( :, :    ), Allocatable :: tmp_c

    Real( wp ), Dimension( :, :, :, : ), Allocatable :: A_r
    Real( wp ), Dimension( :, :    ), Allocatable :: tmp_r

    Real( wp ) :: rand
    Real( wp ) :: max_diff

    Integer, Dimension( 1:3, 1:nk ) :: k_points
    Integer, Dimension(      1:nk ) :: k_types

    Integer :: k, s
    Integer :: i
    Integer :: n
    Integer :: error

    n = m

    Allocate( A_r( 1:n, 1:n, 1:nk, 1:ns ) )
    Allocate( A_c( 1:n, 1:n, 1:nk, 1:ns ) )
    Allocate( tmp_r( 1:n, 1:n ) )
    Allocate( tmp_c( 1:n, 1:n ) )

    A_r = Huge( A_r )
    A_c = Huge( Real( A_c, Kind( A_c ) ) )

    If( me == 0 ) Then

       Call Random_Number( A_r )
       Do s = 1, ns
          Do k = 1, nk
             ! Make sure matrix is positive definite
             A_r( :, :, k, s ) = Matmul( A_r( :, :, k, s ), Transpose( A_r( :, :, k, s  ) ) )
          End Do
       End Do

       Do s = 1, ns
          Do k = 1, nk
             Call Random_Number( tmp_r )
             A_c( :, :, k, s ) = tmp_r
             Call Random_Number( tmp_r )
             A_c( :, :, k, s ) = A_c( :, :, k, s ) + Cmplx( 0.0_wp, tmp_r, Kind = wp )
             ! Make sure matrix is positive definite
             A_c( :, :, k, s ) = Matmul( A_c( :, :, k, s ), Conjg( Transpose( A_c( :, :, k, s  ) ) ) )
          End Do
       End Do

       Do k = 1, nk
          k_points( :, k ) = [ k - 1, 0, 0 ]
          Call Random_Number( rand )
          k_types( k ) = Merge( K_POINT_REAL, K_POINT_COMPLEX, rand > 0.5_wp )
       End Do

    End If

    Call mpi_bcast( k_points, Size( k_points ), mpi_integer, 0, mpi_comm_world, error )
    Call mpi_bcast( k_types , Size( k_types  ), mpi_integer, 0, mpi_comm_world, error )

    Call mpi_bcast( A_r, Size( A_r ), mpi_double_precision, 0, mpi_comm_world, error )

    Call mpi_bcast( A_c, Size( A_c ), mpi_double_complex  , 0, mpi_comm_world, error )

    Call ks_array_init( n_block )
    Call ks_array_comm_to_base( MPI_COMM_WORLD, ns, k_types, k_points, base_k )

    Call A%create( n, n, base_k )
    If( verbose ) Then
       Call A%print_info( 'A', 200 )
    End If

    Call A%split_ks( 2.0_wp, A_split )
    If( verbose ) Then
       Call A_split%print_info( 'A_split', 100 )
    End If

    Call A_split%iterator_init()
    Do
       info = A_split%iterator_next()
       If( info%k_type == K_POINT_NOT_EXIST ) Exit
       Outer: Do s = 1, ns
          Do k = 1, nk
             If( All( info%k_indices == k_points( :, k ) .And. info%spin == s ) ) Then
                If( info%k_type == K_POINT_REAL ) Then
                   Call A_split%set_by_global( k_points( :, k ), s, 1, n, 1, n, A_r( :, :, k, s ) )
                Else
                   Call A_split%set_by_global( k_points( :, k ), s, 1, n, 1, n, A_c( :, :, k, s ) )
                End If
                Exit Outer
             End If
          End Do
       End Do Outer
    End Do
    Call A_split%iterator_reset()

    ! Generate a traingular array by Choleski Decomposition
    max_diff = -1.0_wp
    Call A_split%iterator_init()
    Do
       info = A_split%iterator_next()
       If( info%k_type == K_POINT_NOT_EXIST ) Exit

       A_split_element = A_split%extract_current_element()

       L = .Choleski. A_split_element
       L_inv = .TrInv. L
       B = L * L_inv
       If( verbose ) Then
          Call B%print_info( 'B', 100 )
       End If

       Do s = 1, ns
          Do k = 1, nk
             If( All( info%k_indices == k_points( :, k ) .And. info%spin == s ) ) Then
                tmp_r = 0.0_wp
                If( k_types( k ) == K_POINT_REAL ) Then
                   ! B should be the unit matrix
                   Call B%get_by_global( k_points( :, k ), s, 1, n, 1, n, tmp_r )
                   tmp_r = Abs( tmp_r )
                Else
                   Call B%get_by_global( k_points( :, k ), s, 1, n, 1, n, tmp_c )
                   tmp_r = Abs( tmp_c )
                End If
                Do i = 1, n
                   tmp_r( i, i ) = tmp_r( i, i ) - 1.0_wp
                End Do
                max_diff = Max( max_diff, Maxval( tmp_r ) )
             End If
          End Do
       End Do

    End Do

    If( me == 0 ) Then
       Write( *, error_format ) 'Error in ks_split tr_inv it + split ', max_diff, &
            Merge( passed, FAILED, max_diff < tol )
    End If

    Call ks_array_finalise

  End Subroutine test_ks_array_tr_inv_with_iterator_and_extract

  Subroutine test_ks_array_add_with_iterator_and_extract()

    Use numbers_module , Only : wp
    Use ks_array_module, Only : ks_array, ks_point_info, ks_array_init, ks_array_comm_to_base, ks_array_finalise, &
         K_POINT_REAL, K_POINT_COMPLEX, K_POINT_NOT_EXIST
    Use mpi            , Only : mpi_comm_world, mpi_double_complex, mpi_double_precision

    Integer, Parameter :: ns = 2
    Integer, Parameter :: nk = 3

    Type( ks_array ) :: A, C_split, A_split
    Type( ks_array ) :: A_split_element, C_split_element
    Type( ks_array ) :: B

    Type( ks_array ) :: base_k

    Type( ks_point_info ), Allocatable :: info

    Complex( wp ), Dimension( :, :, :, : ), Allocatable :: A_c
    Complex( wp ), Dimension( :, :    ), Allocatable :: tmp_c

    Real( wp ), Dimension( :, :, :, : ), Allocatable :: A_r
    Real( wp ), Dimension( :, :    ), Allocatable :: tmp_r

    Real( wp ) :: rand
    Real( wp ) :: max_diff, max_diff_ks

    Integer, Dimension( 1:3, 1:nk ) :: k_points
    Integer, Dimension(      1:nk ) :: k_types

    Integer :: k, s
    Integer :: n
    Integer :: error

    n = m

    Allocate( A_r( 1:n, 1:n, 1:nk, 1:ns ) )
    Allocate( A_c( 1:n, 1:n, 1:nk, 1:ns ) )
    Allocate( tmp_r( 1:n, 1:n ) )
    Allocate( tmp_c( 1:n, 1:n ) )

    A_r = Huge( A_r )
    A_c = Huge( Real( A_c, Kind( A_c ) ) )

    If( me == 0 ) Then

       Call Random_Number( A_r )

       Do s = 1, ns
          Do k = 1, nk
             Call Random_Number( tmp_r )
             A_c( :, :, k, s ) = tmp_r
             Call Random_Number( tmp_r )
             A_c( :, :, k, s ) = A_c( :, :, k, s ) + Cmplx( 0.0_wp, tmp_r, Kind = wp )
          End Do
       End Do

       Do k = 1, nk
          k_points( :, k ) = [ k - 1, 0, 0 ]
          Call Random_Number( rand )
          k_types( k ) = Merge( K_POINT_REAL, K_POINT_COMPLEX, rand > 0.5_wp )
       End Do

    End If

    Call mpi_bcast( k_points, Size( k_points ), mpi_integer, 0, mpi_comm_world, error )
    Call mpi_bcast( k_types , Size( k_types  ), mpi_integer, 0, mpi_comm_world, error )

    Call mpi_bcast( A_r, Size( A_r ), mpi_double_precision, 0, mpi_comm_world, error )

    Call mpi_bcast( A_c, Size( A_c ), mpi_double_complex  , 0, mpi_comm_world, error )

    Call ks_array_init( n_block )
    Call ks_array_comm_to_base( MPI_COMM_WORLD, ns, k_types, k_points, base_k )

    Call A%create( n, n, base_k )
    If( verbose ) Then
       Call A%print_info( 'A', 200 )
    End If

    Call A%split_ks( 2.0_wp, A_split )
    If( verbose ) Then
       Call A_split%print_info( 'A_split', 100 )
    End If

    Call A_split%iterator_init()
    Do
       info = A_split%iterator_next()
       If( info%k_type == K_POINT_NOT_EXIST ) Exit
       Outer: Do s = 1, ns
          Do k = 1, nk
             If( All( info%k_indices == k_points( :, k ) .And. info%spin == s ) ) Then
                If( info%k_type == K_POINT_REAL ) Then
                   Call A_split%set_by_global( k_points( :, k ), s, 1, n, 1, n, A_r( :, :, k, s ) )
                Else
                   Call A_split%set_by_global( k_points( :, k ), s, 1, n, 1, n, A_c( :, :, k, s ) )
                End If
                Exit Outer
             End If
          End Do
       End Do Outer
    End Do
    C_split = 3.0_wp * A_split

    Call A_split%iterator_reset()

    ! Generate a traingular array by Choleski Decomposition
    max_diff = -1.0_wp
    Call A_split%iterator_init()
    Call C_split%iterator_init()
    Do
       info = A_split%iterator_next()
       If( info%k_type == K_POINT_NOT_EXIST ) Exit
       info = C_split%iterator_next()

       A_split_element = A_split%extract_current_element()
       C_split_element = C_split%extract_current_element()

       B = A_split_element + C_split_element
       If( verbose ) Then
          Call B%print_info( 'B', 100 )
       End If

       Do s = 1, ns
          Do k = 1, nk
             If( All( info%k_indices == k_points( :, k ) .And. info%spin == s ) ) Then
                tmp_r = 0.0_wp
                If( k_types( k ) == K_POINT_REAL ) Then
                   ! B should be the unit matrix
                   Call B%get_by_global( k_points( :, k ), s, 1, n, 1, n, tmp_r )
                   max_diff_ks = Maxval( Abs( tmp_r - 4.0_wp * A_r( :, :, k, s ) ) )
                Else
                   Call B%get_by_global( k_points( :, k ), s, 1, n, 1, n, tmp_c )
                   max_diff_ks = Maxval( Abs( tmp_c - 4.0_wp * A_c( :, :, k, s ) ) )
                End If
                max_diff = Max( max_diff, max_diff_ks )
             End If
          End Do
       End Do

    End Do

    If( me == 0 ) Then
       Write( *, error_format ) 'Error in ks_split add it + split ', max_diff, &
            Merge( passed, FAILED, max_diff < tol )
    End If

    Call ks_array_finalise

  End Subroutine test_ks_array_add_with_iterator_and_extract

End Module ks_diag_tests


! Note: these implementations are not mine, but derive from ACM QUICKSORT - ALGORITHM #402 (W. H. Verity's implementation)

SUBROUTINE qsort_int (a,ord,n)
  IMPLICIT NONE

  INTEGER, intent(in)  :: a(n)
  INTEGER, intent(in)  :: n
  integer, INTENT(out) :: ord(n)

  !==============SORTS THE ARRAY A(I),I=1,2,...,N BY PUTTING THE
  !   ASCENDING ORDER VECTOR IN ORD.  THAT IS ASCENDING ORDERED A
  !   IS A(ORD(I)),I=1,2,...,N; DESCENDING ORDER A IS A(ORD(N-I+1)),
  !   I=1,2,...,N .  THIS SORT RUNS IN TIME PROPORTIONAL TO N LOG N .

  !     ACM QUICKSORT - ALGORITHM #402 - IMPLEMENTED IN FORTRAN 66 BY
  !                                 WILLIAM H. VERITY, WHV@PSUVM.PSU.EDU
  !                                 CENTER FOR ACADEMIC COMPUTING
  !                                 THE PENNSYLVANIA STATE UNIVERSITY
  !                                 UNIVERSITY PARK, PA.  16802

  INTEGER :: poplst(2,20)
  integer :: x,xx,z,zz,y
  INTEGER :: yp, ndeep, u1, l1, i, u, l, p, q, ix, iz, ip, iq

  !     TO SORT DIFFERENT INPUT TYPES, CHANGE THE FOLLOWING
  !     SPECIFICATION STATEMENTS; FOR EXAMPLE, FOR FORTRAN CHARACTER
  !     USE THE FOLLOWING:  CHARACTER *(*) A(N)

  ndeep=0
  u1=n
  l1=1
  DO   i=1,n
     ord(i)=i
  END DO
2 IF (u1 <= l1) then
     RETURN
  end if

3 l=l1
  u=u1

  ! PART

4 p=l
  q=u
  !     FOR CHARACTER SORTS, THE FOLLOWING 3 STATEMENTS WOULD BECOME
  !     X = ORD(P)
  !     Z = ORD(Q)
  !     IF (A(X) .LE. A(Z)) GO TO 2

  !     WHERE "CLE" IS A LOGICAL FUNCTION WHICH RETURNS "TRUE" IF THE
  !     FIRST ARGUMENT IS LESS THAN OR EQUAL TO THE SECOND, BASED ON "LEN"
  !     CHARACTERS.

  x=a(ord(p))
  z=a(ord(q))
  IF (x <= z) GO TO 5
  y=x
  x=z
  z=y
  yp=ord(p)
  ord(p)=ord(q)
  ord(q)=yp
5 IF (u-l <= 1) GO TO 15
  xx=x
  ix=p
  zz=z
  iz=q

  ! LEFT

6 p=p+1
  IF (p >= q) GO TO 7
  x=a(ord(p))
  IF (x >= xx) GO TO 8
  GO TO 6
7 p=q-1
  GO TO 13

  ! RIGHT

8 q=q-1
  IF (q <= p) GO TO 9
  z=a(ord(q))
  IF (z <= zz) GO TO 10
  GO TO 8
9 q=p
  p=p-1
  z=x
  x=a(ord(p))

  ! DIST

10 IF (x <= z) GO TO 11
  y=x
  x=z
  z=y
  ip=ord(p)
  ord(p)=ord(q)
  ord(q)=ip
11 IF (x <= xx) GO TO 12
  xx=x
  ix=p
12 IF (z >= zz) GO TO 6
  zz=z
  iz=q
  GO TO 6

  ! OUT

13 CONTINUE
  IF (.NOT.(p /= ix.AND.x /= xx)) GO TO 14
  ip=ord(p)
  ord(p)=ord(ix)
  ord(ix)=ip
14 CONTINUE
  IF (.NOT.(q /= iz.AND.z /= zz)) GO TO 15
  iq=ord(q)
  ord(q)=ord(iz)
  ord(iz)=iq
15 CONTINUE
  IF (u-q <= p-l) GO TO 16
  l1=l
  u1=p-1
  l=q+1
  GO TO 17
16 u1=u
  l1=q+1
  u=p-1
17 CONTINUE
  IF (u1 <= l1) GO TO 18

  ! START RECURSIVE CALL

  ndeep=ndeep+1
  poplst(1,ndeep)=u
  poplst(2,ndeep)=l
  GO TO 3
18 IF (u > l) GO TO 4

  ! POP BACK UP IN THE RECURSION LIST

  IF (ndeep == 0) GO TO 2
  u=poplst(1,ndeep)
  l=poplst(2,ndeep)
  ndeep=ndeep-1
  GO TO 18

END SUBROUTINE qsort_int

SUBROUTINE qsort_real (a, ord,n)
  IMPLICIT NONE

  real,    intent(in)  :: a(n)
  INTEGER, intent(in)  :: n
  integer, INTENT(out) :: ord(n)

  !==============SORTS THE ARRAY A(I),I=1,2,...,N BY PUTTING THE
  !   ASCENDING ORDER VECTOR IN ORD.  THAT IS ASCENDING ORDERED A
  !   IS A(ORD(I)),I=1,2,...,N; DESCENDING ORDER A IS A(ORD(N-I+1)),
  !   I=1,2,...,N .  THIS SORT RUNS IN TIME PROPORTIONAL TO N LOG N .

  !     ACM QUICKSORT - ALGORITHM #402 - IMPLEMENTED IN FORTRAN 66 BY
  !                                 WILLIAM H. VERITY, WHV@PSUVM.PSU.EDU
  !                                 CENTER FOR ACADEMIC COMPUTING
  !                                 THE PENNSYLVANIA STATE UNIVERSITY
  !                                 UNIVERSITY PARK, PA.  16802

  INTEGER :: poplst(2,20)
  REAL :: x,xx,z,zz,y
  INTEGER :: yp, ndeep, u1, l1, i, u, l, p, q, ix, iz, ip, iq

  !     TO SORT DIFFERENT INPUT TYPES, CHANGE THE FOLLOWING
  !     SPECIFICATION STATEMENTS; FOR EXAMPLE, FOR FORTRAN CHARACTER
  !     USE THE FOLLOWING:  CHARACTER *(*) A(N)

  ndeep=0
  u1=n
  l1=1
  DO   i=1,n
     ord(i)=i
  END DO
2 IF (u1 <= l1) then
     RETURN
  end if

3 l=l1
  u=u1

  ! PART

4 p=l
  q=u
  !     FOR CHARACTER SORTS, THE FOLLOWING 3 STATEMENTS WOULD BECOME
  !     X = ORD(P)
  !     Z = ORD(Q)
  !     IF (A(X) .LE. A(Z)) GO TO 2

  !     WHERE "CLE" IS A LOGICAL FUNCTION WHICH RETURNS "TRUE" IF THE
  !     FIRST ARGUMENT IS LESS THAN OR EQUAL TO THE SECOND, BASED ON "LEN"
  !     CHARACTERS.

  x=a(ord(p))
  z=a(ord(q))
  IF (x <= z) GO TO 5
  y=x
  x=z
  z=y
  yp=ord(p)
  ord(p)=ord(q)
  ord(q)=yp
5 IF (u-l <= 1) GO TO 15
  xx=x
  ix=p
  zz=z
  iz=q

  ! LEFT

6 p=p+1
  IF (p >= q) GO TO 7
  x=a(ord(p))
  IF (x >= xx) GO TO 8
  GO TO 6
7 p=q-1
  GO TO 13

  ! RIGHT

8 q=q-1
  IF (q <= p) GO TO 9
  z=a(ord(q))
  IF (z <= zz) GO TO 10
  GO TO 8
9 q=p
  p=p-1
  z=x
  x=a(ord(p))

  ! DIST

10 IF (x <= z) GO TO 11
  y=x
  x=z
  z=y
  ip=ord(p)
  ord(p)=ord(q)
  ord(q)=ip
11 IF (x <= xx) GO TO 12
  xx=x
  ix=p
12 IF (z >= zz) GO TO 6
  zz=z
  iz=q
  GO TO 6

  ! OUT

13 CONTINUE
  IF (.NOT.(p /= ix.AND.x /= xx)) GO TO 14
  ip=ord(p)
  ord(p)=ord(ix)
  ord(ix)=ip
14 CONTINUE
  IF (.NOT.(q /= iz.AND.z /= zz)) GO TO 15
  iq=ord(q)
  ord(q)=ord(iz)
  ord(iz)=iq
15 CONTINUE
  IF (u-q <= p-l) GO TO 16
  l1=l
  u1=p-1
  l=q+1
  GO TO 17
16 u1=u
  l1=q+1
  u=p-1
17 CONTINUE
  IF (u1 <= l1) GO TO 18

  ! START RECURSIVE CALL

  ndeep=ndeep+1
  poplst(1,ndeep)=u
  poplst(2,ndeep)=l
  GO TO 3
18 IF (u > l) GO TO 4

  ! POP BACK UP IN THE RECURSION LIST

  IF (ndeep == 0) GO TO 2
  u=poplst(1,ndeep)
  l=poplst(2,ndeep)
  ndeep=ndeep-1
  GO TO 18

END SUBROUTINE qsort_real


  program analytical_solution

! This program implements the analytical solution for the displacement vector in a 3D viscoelastic medium
! with a vertical force source located in (0,0,0).
! Implemented by Dimitri Komatitsch from the 2D plane-strain viscoelastic medium analytical solution
! of Appendix B of Carcione et al., Wave propagation simulation in a linear viscoelastic medium, GJI, vol. 95, p. 597-611 (1988)
! (note that that Appendix contains two typos, fixed in this code).
! The amplitude of the force is called F and is defined below.

!! DK DK Dimitri Komatitsch, CNRS Marseille, France, May 2018.

  implicit none

!! DK DK May 2018: the missing 1/L factor in older Carcione papers
!! DK DK May 2018: has been added to this code by Quentin Brissaud and by Etienne Bachmann
!! DK DK for the viscoacoustic code in directory EXAMPLES/attenuation/viscoacoustic,
!! DK DK it would be very easy to copy the changes from there to this viscoelastic version;
!! DK DK but then all the values of the tau_epsilon in the code below would need to change.

!! DK DK Dimitri Komatitsch, CNRS Marseille, France, April 2017: added the elastic reference calculation.

! compute the non-viscoacoustic case as a reference if needed, i.e. turn attenuation off
  logical, parameter :: TURN_ATTENUATION_OFF = .true.

! to see how small the contribution of the near-field term is,
! here the user can ask not to include it, to then compare with the full result obtained with this flag set to false
  logical, parameter :: DO_NOT_COMPUTE_THE_NEAR_FIELD = .false.

  integer, parameter :: iratio = 32

  integer, parameter :: nfreq = 524288
  integer, parameter :: nt = iratio * nfreq

  double precision, parameter :: freqmax = 80.d0
!! DK DK to print the velocity if we want to display the curve of how velocity varies with frequency
!! DK DK for instance to compute the unrelaxed velocity in the Zener model
! double precision, parameter :: freqmax = 20000.d0

  double precision, parameter :: freqseuil = 0.00005d0

  double precision, parameter :: pi = 3.141592653589793d0

! for the solution in time domain
  integer it,i
  real wsave(4*nt+15)
  complex c(nt)

!! DK DK for my slow inverse Discrete Fourier Transform using a double loop
  complex :: input(nt), i_imaginary_constant
  integer :: j,m

! density of the medium
  double precision, parameter :: rho = 2000.d0

! amplitude of the force source
  double precision, parameter :: F = 1.d0

! spatial dimension of the problem (this analytical code is valid only in three dimensions, thus do not change this)
  integer, parameter :: NDIM = 3

! definition position recepteur Carcione
  double precision, dimension(NDIM) :: x

! Definition source Dimitri
  double precision, parameter :: f0 = 18.d0
  double precision, parameter :: t0 = 1.2d0 / f0

! Definition source Carcione
! double precision f0,t0,eta,epsil
! parameter(f0 = 50.d0)
! parameter(t0 = 0.075d0)
! parameter(epsil = 1.d0)
! parameter(eta = 0.5d0)

! number of Zener standard linear solids in parallel
  integer, parameter :: Lnu = 3

! DK DK I implemented a very simple and slow inverse Discrete Fourier Transform
! DK DK at some point, for verification, using a double loop. I keep it just in case.
! DK DK For large number of points it is extremely slow because of the double loop.
! DK DK Thus there is no reason to turn this flag on.
  logical, parameter :: USE_SLOW_FOURIER_TRANSFORM = .false.

!! DK DK March 2018: this missing 1/L factor has been added to this code by Quentin Brissaud
!! DK DK for the viscoacoustic code in directory EXAMPLES/attenuation/viscoacoustic,
!! DK DK it would be very easy to copy the changes from there to this viscoelastic version;
!! DK DK but then all the values of the tau_epsilon below would need to change.

 double precision, dimension(Lnu) :: tau_sigma_kappa,tau_sigma_mu,tau_epsilon_kappa,tau_epsilon_mu

  integer :: ifreq
  double precision :: deltafreq,freq,omega,omega0,deltat,time,a
  double complex :: comparg

! Fourier transform of the Ricker wavelet source
  double complex fomega(0:nfreq)

! real and imaginary parts
  double precision ra(0:nfreq),rb(0:nfreq)

! spectral amplitude
  double precision ampli(0:nfreq)

! analytical solution for the three components
  double complex phi1(-nfreq:nfreq)
  double complex phi2(-nfreq:nfreq)
  double complex phi3(-nfreq:nfreq)

! external function
  double complex, external :: ui

! modules elastiques
  double complex :: Kappa_omega, Mu_omega, Vp_omega, Vs_omega, temp, Mu_relaxed, Kappa_relaxed

! ********** end of variable declarations ************

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!  Below are all the variables that need to be obtained from Specfem3D
!  These default values match the following simulation:
!
!  Vp(f0_ref) = 3297.849 m/s
!  Vs(f0_ref) = 2222.536 m/s
!  f0_ref = 18 Hz
!  QKappa = 9000
!  QMu = 10
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! We use Kappa_ref and Mu_ref to compute the case without attenuation
! They are chosen such that Kappa_ref = Kappa(2*pi*f0_ref) when there is
! attenuation
  double precision, parameter :: Vp_ref  = 3297.849
  double precision, parameter :: Vs_ref  = 2222.536

  double precision, parameter :: Mu_ref    =  rho * Vs_ref * Vs_ref
  double precision, parameter :: Kappa_ref =  rho * ( Vp_ref * Vp_ref - 4.d0/3 * Vs_ref * Vs_ref)

! add this for the case of an elastic medium, with viscoelasticity turned off
  double precision, parameter :: Mu_unrelaxed    =  Mu_ref
  double precision, parameter :: Kappa_unrelaxed =  Kappa_ref

! tau constants mimicking constant QKappa and KMu
  tau_epsilon_mu    = (/ 0.281966668348107  ,  3.607809663879573E-002 , 5.638875613224542E-003/)
  tau_sigma_mu      = (/ 0.186873539567019  ,  2.491998701168405E-002 , 3.323133676931235E-003/)
  tau_epsilon_kappa = (/ 0.186973292921151  ,  2.492998179955646E-002 , 3.324907855424433E-003/)
  tau_sigma_kappa   = (/ 0.186873539567019  ,  2.491998701168405E-002 , 3.323133676931235E-003/)

! Eq (32) of Jeroen's note, eq (2.199) of Carcione's book from 2014, third edition
  Kappa_relaxed = (Kappa_unrelaxed /(sum(tau_epsilon_kappa(:)/tau_sigma_kappa(:))/Lnu))
  Mu_relaxed    = (Mu_unrelaxed    /(sum(tau_epsilon_mu(:)/tau_sigma_mu(:))/Lnu))

! position of the receiver
  x(1) = +500.
  x(2) = +500.
  x(3) = +500.

  print *,'Force source located at the origin (0,0,0)'
  print *,'Receiver located in (x,y,z) = ',x(1),x(2),x(3)

  if (TURN_ATTENUATION_OFF) then
    print *,'BEWARE: computing the elastic reference solution (i.e., without attenuation) instead of the viscoelastic solution'
  else
    print *,'Computing the viscoelastic solution'
  endif

  if (DO_NOT_COMPUTE_THE_NEAR_FIELD) then
    print *,'BEWARE: computing the far-field solution only, rather than the full Green function'
  else
    print *,'Computing the full solution, including the near-field term of the Green function'
  endif

! step in frequency
  deltafreq = freqmax / dble(nfreq)

! define parameters for the Ricker source
  omega0 = 2.d0 * pi * f0
  a = pi**2 * f0**2

  deltat = 1.d0 / (freqmax*dble(iratio))

! define the spectrum of the source
  do ifreq=0,nfreq
      freq = deltafreq * dble(ifreq)
      omega = 2.d0 * pi * freq

! typo in equation (B7) of Carcione et al., Wave propagation simulation in a linear viscoelastic medium,
! Geophysical Journal, vol. 95, p. 597-611 (1988), the exponential should be of -i omega t0,
! fixed here by adding the minus sign
      comparg = dcmplx(0.d0,-omega*t0)

! definir le spectre du Ricker de Carcione avec cos()
! equation (B7) of Carcione et al., Wave propagation simulation in a linear viscoelastic medium,
! Geophysical Journal, vol. 95, p. 597-611 (1988)
!     fomega(ifreq) = pi * dsqrt(pi/eta) * (1.d0/omega0) * cdexp(comparg) * ( dexp(- (pi*pi/eta) * (epsil/2 - omega/omega0)**2) &
!         + dexp(- (pi*pi/eta) * (epsil/2 + omega/omega0)**2) )

! definir le spectre d'un Ricker classique (centre en t0)
      fomega(ifreq) = dsqrt(pi) * cdexp(comparg) * omega**2 * dexp(-omega**2/(4.d0*a)) / (2.d0 * dsqrt(a**3))

      ra(ifreq) = dreal(fomega(ifreq))
      rb(ifreq) = dimag(fomega(ifreq))
! prendre le module de l'amplitude spectrale
      ampli(ifreq) = dsqrt(ra(ifreq)**2 + rb(ifreq)**2)
  enddo

! sauvegarde du spectre d'amplitude de la source en Hz au format Gnuplot
  open(unit=10,file='spectrum_of_the_source_used.gnu',status='unknown')
  do ifreq = 0,nfreq
    freq = deltafreq * dble(ifreq)
    write(10,*) sngl(freq),sngl(ampli(ifreq))
  enddo
  close(10)

! ************** calcul solution analytique ****************

! d'apres Carcione GJI vol 95 p 611 (1988)
  do ifreq=0,nfreq
      freq = deltafreq * dble(ifreq)
      omega = 2.d0 * pi * freq

! critere ad-hoc pour eviter singularite en zero
  if (freq < freqseuil) omega = 2.d0 * pi * freqseuil

! use standard infinite frequency (unrelaxed) reference,
! in which waves slow down when attenuation is turned on.
  temp = dcmplx(0.d0,0.d0)
  do i=1,Lnu
    temp = temp + (1./Lnu)*dcmplx(1.d0,omega*tau_epsilon_kappa(i)) / dcmplx(1.d0,omega*tau_sigma_kappa(i))
  enddo

! Eq(31) of Jeroen's note
  Kappa_omega = Kappa_relaxed * temp

  temp = dcmplx(0.d0,0.d0)
  do i=1,Lnu
    temp = temp + (1./Lnu)*dcmplx(1.d0,omega*tau_epsilon_mu(i)) / dcmplx(1.d0,omega*tau_sigma_mu(i))
  enddo

! Eq (31) of Jeroen's note
  Mu_omega = Mu_relaxed * temp

  if (TURN_ATTENUATION_OFF) then
! wave speeds and moduli are the same for all omega when there is no attenuation
    Kappa_omega = Kappa_ref
    Mu_omega    = Mu_ref
  endif

  Vs_omega = cdsqrt( Mu_omega / rho)
  Vp_omega = cdsqrt( Kappa_omega / rho + (4.d0/3)*Vs_omega*Vs_omega )

!! DK DK print the velocity if we want to display the curve of how velocity varies with frequency
!! DK DK for instance to compute the unrelaxed velocity in the Zener model
! print *,freq,dsqrt(real(Vp_omega)**2 + imag(Vp_omega)**2)
!! DK DK print the velocity if we want to display the curve of how velocity varies with frequency
!! DK DK for instance to compute the unrelaxed velocity in the Zener model
! print *,freq,dsqrt(real(Vs_omega)**2 + imag(Vs_omega)**2)

! calcul de la solution analytique en frequence
  phi1(ifreq) = ui(1,omega,Vp_omega,Vs_omega,x,rho,NDIM,F,DO_NOT_COMPUTE_THE_NEAR_FIELD) * fomega(ifreq)
  phi2(ifreq) = ui(2,omega,Vp_omega,Vs_omega,x,rho,NDIM,F,DO_NOT_COMPUTE_THE_NEAR_FIELD) * fomega(ifreq)
  phi3(ifreq) = ui(3,omega,Vp_omega,Vs_omega,x,rho,NDIM,F,DO_NOT_COMPUTE_THE_NEAR_FIELD) * fomega(ifreq)

  enddo

! take the conjugate value for negative frequencies
  do ifreq=-nfreq,-1
      phi1(ifreq) = dconjg(phi1(-ifreq))
      phi2(ifreq) = dconjg(phi2(-ifreq))
      phi3(ifreq) = dconjg(phi3(-ifreq))
  enddo

! ***************************************************************************
! Calculation of the time domain solution (using routine "cfftb" from Netlib)
! ***************************************************************************

! **********
! Compute Ux
! **********

! initialize FFT arrays
  call cffti(nt,wsave)

! clear array of Fourier coefficients
  do it = 1,nt
    c(it) = cmplx(0.,0.)
  enddo

! use the Fourier values for Ux
  c(1) = cmplx(phi1(0))
  do ifreq=1,nfreq-2
      c(ifreq+1) = cmplx(phi1(ifreq))
      c(nt+1-ifreq) = conjg(cmplx(phi1(ifreq)))
  enddo

! perform the inverse FFT for Ux
  if (.not. USE_SLOW_FOURIER_TRANSFORM) then
    call cfftb(nt,c,wsave)
  else
! DK DK I implemented a very simple and slow inverse Discrete Fourier Transform here
! DK DK at some point, for verification, using a double loop. I keep it just in case.
! DK DK For large number of points it is extremely slow because of the double loop.
    input(:) = c(:)
!   imaginary constant "i"
    i_imaginary_constant = (0.,1.)
    do it = 1,nt
      if (mod(it,1000) == 0) print *,'FFT inverse it = ',it,' out of ',nt
      j = it
      c(j) = cmplx(0.,0.)
      do m = 1,nt
        c(j) = c(j) + input(m) * exp(2.d0 * PI * i_imaginary_constant * dble((m-1) * (j-1)) / nt)
      enddo
    enddo
  endif

! in the inverse Discrete Fourier transform one needs to divide by N, the number of samples (number of time steps here)
  c(:) = c(:) / nt

! value of a time step
  deltat = 1.d0 / (freqmax*dble(iratio))

! to get the amplitude right, we need to divide by the time step
  c(:) = c(:) / deltat

! save time result inverse FFT for Ux

  if (TURN_ATTENUATION_OFF) then
    if (DO_NOT_COMPUTE_THE_NEAR_FIELD) then
      open(unit=11,file='Ux_time_analytical_solution_elastic_without_near_field.dat',status='unknown')
    else
      open(unit=11,file='Ux_time_analytical_solution_elastic.dat',status='unknown')
    endif
  else
    if (DO_NOT_COMPUTE_THE_NEAR_FIELD) then
      open(unit=11,file='Ux_time_analytical_solution_viscoelastic_without_near_field.dat',status='unknown')
    else
      open(unit=11,file='Ux_time_analytical_solution_viscoelastic.dat',status='unknown')
    endif
  endif
  do it=1,nt
! DK DK Dec 2011: subtract t0 to be consistent with the SPECFEM2D code
        time = dble(it-1)*deltat - t0
! the seismograms are very long due to the very large number of FFT points used,
! thus keeping the useful part of the signal only (the first six seconds of the seismogram)
        if (time >= 0.d0 .and. time <= 6.d0) write(11,*) sngl(time),real(c(it))
  enddo
  close(11)

! **********
! Compute Uy
! **********

! clear array of Fourier coefficients
  do it = 1,nt
    c(it) = cmplx(0.,0.)
  enddo

! use the Fourier values for Uy
  c(1) = cmplx(phi2(0))
  do ifreq=1,nfreq-2
      c(ifreq+1) = cmplx(phi2(ifreq))
      c(nt+1-ifreq) = conjg(cmplx(phi2(ifreq)))
  enddo

! perform the inverse FFT for Uy
  if (.not. USE_SLOW_FOURIER_TRANSFORM) then
    call cfftb(nt,c,wsave)
  else
! DK DK I implemented a very simple and slow inverse Discrete Fourier Transform here
! DK DK at some point, for verification, using a double loop. I keep it just in case.
! DK DK For large number of points it is extremely slow because of the double loop.
    input(:) = c(:)
!   imaginary constant "i"
    i_imaginary_constant = (0.,1.)
    do it = 1,nt
      if (mod(it,1000) == 0) print *,'FFT inverse it = ',it,' out of ',nt
      j = it
      c(j) = cmplx(0.,0.)
      do m = 1,nt
        c(j) = c(j) + input(m) * exp(2.d0 * PI * i_imaginary_constant * dble((m-1) * (j-1)) / nt)
      enddo
    enddo
  endif

! in the inverse Discrete Fourier transform one needs to divide by N, the number of samples (number of time steps here)
  c(:) = c(:) / nt

! value of a time step
  deltat = 1.d0 / (freqmax*dble(iratio))

! to get the amplitude right, we need to divide by the time step
  c(:) = c(:) / deltat

! save time result inverse FFT for Uy
  if (TURN_ATTENUATION_OFF) then
    if (DO_NOT_COMPUTE_THE_NEAR_FIELD) then
      open(unit=11,file='Uy_time_analytical_solution_elastic_without_near_field.dat',status='unknown')
    else
      open(unit=11,file='Uy_time_analytical_solution_elastic.dat',status='unknown')
    endif
  else
    if (DO_NOT_COMPUTE_THE_NEAR_FIELD) then
      open(unit=11,file='Uy_time_analytical_solution_viscoelastic_without_near_field.dat',status='unknown')
    else
      open(unit=11,file='Uy_time_analytical_solution_viscoelastic.dat',status='unknown')
    endif
  endif
  do it=1,nt
! DK DK Dec 2011: subtract t0 to be consistent with the SPECFEM2D code
        time = dble(it-1)*deltat - t0
! the seismograms are very long due to the very large number of FFT points used,
! thus keeping the useful part of the signal only (the first six seconds of the seismogram)
        if (time >= 0.d0 .and. time <= 6.d0) write(11,*) sngl(time),real(c(it))
  enddo
  close(11)

! **********
! Compute Uz
! **********

! clear array of Fourier coefficients
  do it = 1,nt
    c(it) = cmplx(0.,0.)
  enddo

! use the Fourier values for Uz
  c(1) = cmplx(phi3(0))
  do ifreq=1,nfreq-2
      c(ifreq+1) = cmplx(phi3(ifreq))
      c(nt+1-ifreq) = conjg(cmplx(phi3(ifreq)))
  enddo

! perform the inverse FFT for Uz
  if (.not. USE_SLOW_FOURIER_TRANSFORM) then
    call cfftb(nt,c,wsave)
  else
! DK DK I implemented a very simple and slow inverse Discrete Fourier Transform here
! DK DK at some point, for verification, using a double loop. I keep it just in case.
! DK DK For large number of points it is extremely slow because of the double loop.
    input(:) = c(:)
!   imaginary constant "i"
    i_imaginary_constant = (0.,1.)
    do it = 1,nt
      if (mod(it,1000) == 0) print *,'FFT inverse it = ',it,' out of ',nt
      j = it
      c(j) = cmplx(0.,0.)
      do m = 1,nt
        c(j) = c(j) + input(m) * exp(2.d0 * PI * i_imaginary_constant * dble((m-1) * (j-1)) / nt)
      enddo
    enddo
  endif

! in the inverse Discrete Fourier transform one needs to divide by N, the number of samples (number of time steps here)
  c(:) = c(:) / nt

! value of a time step
  deltat = 1.d0 / (freqmax*dble(iratio))

! to get the amplitude right, we need to divide by the time step
  c(:) = c(:) / deltat

! save time result inverse FFT for Uz
  if (TURN_ATTENUATION_OFF) then
    if (DO_NOT_COMPUTE_THE_NEAR_FIELD) then
      open(unit=11,file='Uz_time_analytical_solution_elastic_without_near_field.dat',status='unknown')
    else
      open(unit=11,file='Uz_time_analytical_solution_elastic.dat',status='unknown')
    endif
  else
    if (DO_NOT_COMPUTE_THE_NEAR_FIELD) then
      open(unit=11,file='Uz_time_analytical_solution_viscoelastic_without_near_field.dat',status='unknown')
    else
      open(unit=11,file='Uz_time_analytical_solution_viscoelastic.dat',status='unknown')
    endif
  endif
  do it=1,nt
! DK DK Dec 2011: subtract t0 to be consistent with the SPECFEM2D code
        time = dble(it-1)*deltat - t0
! the seismograms are very long due to the very large number of FFT points used,
! thus keeping the useful part of the signal only (the first six seconds of the seismogram)
        if (time >= 0.d0 .and. time <= 6.d0) write(11,*) sngl(time),real(c(it))
  enddo
  close(11)

  end

! -----------

  double complex function ui(i,omega,v1,v2,x,rho,NDIM,F,DO_NOT_COMPUTE_THE_NEAR_FIELD)

  implicit none

  double precision, parameter :: pi = 3.141592653589793d0

  integer i,j
  double precision omega,F
  double complex :: v1,v2,Up_far_field,Us_far_field,near_field,i_imaginary_constant
  double complex :: fourier_transform_of_Heaviside_times_t_for_Vp,fourier_transform_of_Heaviside_times_t_for_Vs
  logical :: DO_NOT_COMPUTE_THE_NEAR_FIELD

  integer :: NDIM,kronecker_delta_symbol
  double precision, dimension(NDIM) :: x,gamma_vector
  double precision :: r,rho

! source-receiver distance
  r = dsqrt(x(1)**2 + x(2)**2 + x(3)**2)

! define the gamma vector of Aki and Richards (1980), which is the unit vector from the source to the receiver
! see Aki and Richards (1980) below Box 4.3 (continued) and above equation (4.23)
  gamma_vector(:) = x(:) / r

! the force is vertical in this analytical code, thus the j direction of the force is Z
  j = 3

! imaginary constant "i"
  i_imaginary_constant = (0.d0,1.d0)

! Kronecker delta symbol
  if (i == j) then
    kronecker_delta_symbol = 1
  else
    kronecker_delta_symbol = 0
  endif

! far-field P wave term
! see e.g. Aki and Richards (1980), equation (4.24) in the time domain; here we use the same equation but in the frequency domain
  Up_far_field = gamma_vector(i) * gamma_vector(j) * exp(-i_imaginary_constant*omega*r/v1) / (4.d0 * pi * rho * v1**2 * r)

! far-field S wave term
! see e.g. Aki and Richards (1980), equation (4.25) in the time domain; here we use the same equation but in the frequency domain
! note that Aki and Richards (1980) has a typo in equation (4.25), a term 1/r is missing; fixed here.
! The corrected equation is for instance in Jose Pujol, Elastic wave propagation and generation in seismology, equation (9.6.1)
  Us_far_field = (kronecker_delta_symbol - gamma_vector(i) * gamma_vector(j)) * exp(-i_imaginary_constant*omega*r/v2) &
               / (4.d0 * pi * rho * v2**2 * r)

! near field term (see e.g. Jose Pujol, Elastic wave propagation and generation in seismology, equation (9.6.1) in the time domain;
! here we use the same equation but in the frequency domain)
  fourier_transform_of_Heaviside_times_t_for_Vp = i_imaginary_constant * (-omega*r + i_imaginary_constant*v1) * &
                   exp(-i_imaginary_constant*omega*r/v1) / (v1 * omega**2)

  fourier_transform_of_Heaviside_times_t_for_Vs = i_imaginary_constant * (-omega*r + i_imaginary_constant*v2) * &
                   exp(-i_imaginary_constant*omega*r/v2) / (v2 * omega**2)

  near_field = (3 * gamma_vector(i) * gamma_vector(j) - kronecker_delta_symbol) * &
     (fourier_transform_of_Heaviside_times_t_for_Vp - fourier_transform_of_Heaviside_times_t_for_Vs) / (4.d0 * pi * rho * r**3)

  if (DO_NOT_COMPUTE_THE_NEAR_FIELD) near_field = dcmplx(0.d0,0.d0)

! the result is the sum of the three terms; we also multiply by the amplitude of the force
  ui = F * (Up_far_field + Us_far_field + near_field)

  end

! ***************** routine de FFT pour signal en temps ****************

! FFT routine taken from Netlib

  subroutine CFFTB (N,C,WSAVE)
  DIMENSION       C(1)       ,WSAVE(1)
  if (N == 1) return
  IW1 = N+N+1
  IW2 = IW1+N+N
  CALL CFFTB1 (N,C,WSAVE,WSAVE(IW1),WSAVE(IW2))
  return
  END
  subroutine CFFTB1 (N,C,CH,WA,IFAC)
  DIMENSION       CH(1)      ,C(1)       ,WA(1)      ,IFAC(1)
  NF = IFAC(2)
  NA = 0
  L1 = 1
  IW = 1
  DO 116 K1=1,NF
   IP = IFAC(K1+2)
   L2 = IP*L1
   IDO = N/L2
   IDOT = IDO+IDO
   IDL1 = IDOT*L1
   if (IP /= 4) goto 103
   IX2 = IW+IDOT
   IX3 = IX2+IDOT
   if (NA /= 0) goto 101
   CALL PASSB4 (IDOT,L1,C,CH,WA(IW),WA(IX2),WA(IX3))
   goto 102
  101    CALL PASSB4 (IDOT,L1,CH,C,WA(IW),WA(IX2),WA(IX3))
  102    NA = 1-NA
   goto 115
  103    if (IP /= 2) goto 106
   if (NA /= 0) goto 104
   CALL PASSB2 (IDOT,L1,C,CH,WA(IW))
   goto 105
  104    CALL PASSB2 (IDOT,L1,CH,C,WA(IW))
  105    NA = 1-NA
   goto 115
  106    if (IP /= 3) goto 109
   IX2 = IW+IDOT
   if (NA /= 0) goto 107
   CALL PASSB3 (IDOT,L1,C,CH,WA(IW),WA(IX2))
   goto 108
  107    CALL PASSB3 (IDOT,L1,CH,C,WA(IW),WA(IX2))
  108    NA = 1-NA
   goto 115
  109    if (IP /= 5) goto 112
   IX2 = IW+IDOT
   IX3 = IX2+IDOT
   IX4 = IX3+IDOT
   if (NA /= 0) goto 110
   CALL PASSB5 (IDOT,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4))
   goto 111
  110    CALL PASSB5 (IDOT,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4))
  111    NA = 1-NA
   goto 115
  112    if (NA /= 0) goto 113
   CALL PASSB (NAC,IDOT,IP,L1,IDL1,C,C,C,CH,CH,WA(IW))
   goto 114
  113    CALL PASSB (NAC,IDOT,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW))
  114    if (NAC /= 0) NA = 1-NA
  115    L1 = L2
   IW = IW+(IP-1)*IDOT
  116 continue
  if (NA == 0) return
  N2 = N+N
  DO 117 I=1,N2
   C(I) = CH(I)
  117 continue
  return
  END
  subroutine PASSB (NAC,IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA)
  DIMENSION       CH(IDO,L1,IP)          ,CC(IDO,IP,L1), &
                  C1(IDO,L1,IP)          ,WA(1)      ,C2(IDL1,IP), &
                  CH2(IDL1,IP)
  IDOT = IDO/2
  NT = IP*IDL1
  IPP2 = IP+2
  IPPH = (IP+1)/2
  IDP = IP*IDO
!
  if (IDO < L1) goto 106
  DO 103 J=2,IPPH
   JC = IPP2-J
   DO 102 K=1,L1
      DO 101 I=1,IDO
         CH(I,K,J) = CC(I,J,K)+CC(I,JC,K)
         CH(I,K,JC) = CC(I,J,K)-CC(I,JC,K)
  101       continue
  102    continue
  103 continue
  DO 105 K=1,L1
   DO 104 I=1,IDO
      CH(I,K,1) = CC(I,1,K)
  104    continue
  105 continue
  goto 112
  106 DO 109 J=2,IPPH
   JC = IPP2-J
   DO 108 I=1,IDO
      DO 107 K=1,L1
         CH(I,K,J) = CC(I,J,K)+CC(I,JC,K)
         CH(I,K,JC) = CC(I,J,K)-CC(I,JC,K)
  107       continue
  108    continue
  109 continue
  DO 111 I=1,IDO
   DO 110 K=1,L1
      CH(I,K,1) = CC(I,1,K)
  110    continue
  111 continue
  112 IDL = 2-IDO
  INC = 0
  DO 116 L=2,IPPH
   LC = IPP2-L
   IDL = IDL+IDO
   DO 113 IK=1,IDL1
      C2(IK,L) = CH2(IK,1)+WA(IDL-1)*CH2(IK,2)
      C2(IK,LC) = WA(IDL)*CH2(IK,IP)
  113    continue
   IDLJ = IDL
   INC = INC+IDO
   DO 115 J=3,IPPH
      JC = IPP2-J
      IDLJ = IDLJ+INC
      if (IDLJ > IDP) IDLJ = IDLJ-IDP
      WAR = WA(IDLJ-1)
      WAI = WA(IDLJ)
      DO 114 IK=1,IDL1
         C2(IK,L) = C2(IK,L)+WAR*CH2(IK,J)
         C2(IK,LC) = C2(IK,LC)+WAI*CH2(IK,JC)
  114       continue
  115    continue
  116 continue
  DO 118 J=2,IPPH
   DO 117 IK=1,IDL1
      CH2(IK,1) = CH2(IK,1)+CH2(IK,J)
  117    continue
  118 continue
  DO 120 J=2,IPPH
   JC = IPP2-J
   DO 119 IK=2,IDL1,2
      CH2(IK-1,J) = C2(IK-1,J)-C2(IK,JC)
      CH2(IK-1,JC) = C2(IK-1,J)+C2(IK,JC)
      CH2(IK,J) = C2(IK,J)+C2(IK-1,JC)
      CH2(IK,JC) = C2(IK,J)-C2(IK-1,JC)
  119    continue
  120 continue
  NAC = 1
  if (IDO == 2) return
  NAC = 0
  DO 121 IK=1,IDL1
   C2(IK,1) = CH2(IK,1)
  121 continue
  DO 123 J=2,IP
   DO 122 K=1,L1
      C1(1,K,J) = CH(1,K,J)
      C1(2,K,J) = CH(2,K,J)
  122    continue
  123 continue
  if (IDOT > L1) goto 127
  IDIJ = 0
  DO 126 J=2,IP
   IDIJ = IDIJ+2
   DO 125 I=4,IDO,2
      IDIJ = IDIJ+2
      DO 124 K=1,L1
         C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J)
         C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J)
  124       continue
  125    continue
  126 continue
  return
  127 IDJ = 2-IDO
  DO 130 J=2,IP
   IDJ = IDJ+IDO
   DO 129 K=1,L1
      IDIJ = IDJ
      DO 128 I=4,IDO,2
         IDIJ = IDIJ+2
         C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J)
         C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J)
  128       continue
  129    continue
  130 continue
  return
  END
  subroutine PASSB2 (IDO,L1,CC,CH,WA1)
  DIMENSION       CC(IDO,2,L1)           ,CH(IDO,L1,2), &
                  WA1(1)
  if (IDO > 2) goto 102
  DO 101 K=1,L1
   CH(1,K,1) = CC(1,1,K)+CC(1,2,K)
   CH(1,K,2) = CC(1,1,K)-CC(1,2,K)
   CH(2,K,1) = CC(2,1,K)+CC(2,2,K)
   CH(2,K,2) = CC(2,1,K)-CC(2,2,K)
  101 continue
  return
  102 DO 104 K=1,L1
   DO 103 I=2,IDO,2
      CH(I-1,K,1) = CC(I-1,1,K)+CC(I-1,2,K)
      TR2 = CC(I-1,1,K)-CC(I-1,2,K)
      CH(I,K,1) = CC(I,1,K)+CC(I,2,K)
      TI2 = CC(I,1,K)-CC(I,2,K)
      CH(I,K,2) = WA1(I-1)*TI2+WA1(I)*TR2
      CH(I-1,K,2) = WA1(I-1)*TR2-WA1(I)*TI2
  103    continue
  104 continue
  return
  END
  subroutine PASSB3 (IDO,L1,CC,CH,WA1,WA2)
  DIMENSION       CC(IDO,3,L1)           ,CH(IDO,L1,3), &
                  WA1(1)     ,WA2(1)
  DATA TAUR,TAUI /-.5,.866025403784439/
  if (IDO /= 2) goto 102
  DO 101 K=1,L1
   TR2 = CC(1,2,K)+CC(1,3,K)
   CR2 = CC(1,1,K)+TAUR*TR2
   CH(1,K,1) = CC(1,1,K)+TR2
   TI2 = CC(2,2,K)+CC(2,3,K)
   CI2 = CC(2,1,K)+TAUR*TI2
   CH(2,K,1) = CC(2,1,K)+TI2
   CR3 = TAUI*(CC(1,2,K)-CC(1,3,K))
   CI3 = TAUI*(CC(2,2,K)-CC(2,3,K))
   CH(1,K,2) = CR2-CI3
   CH(1,K,3) = CR2+CI3
   CH(2,K,2) = CI2+CR3
   CH(2,K,3) = CI2-CR3
  101 continue
  return
  102 DO 104 K=1,L1
   DO 103 I=2,IDO,2
      TR2 = CC(I-1,2,K)+CC(I-1,3,K)
      CR2 = CC(I-1,1,K)+TAUR*TR2
      CH(I-1,K,1) = CC(I-1,1,K)+TR2
      TI2 = CC(I,2,K)+CC(I,3,K)
      CI2 = CC(I,1,K)+TAUR*TI2
      CH(I,K,1) = CC(I,1,K)+TI2
      CR3 = TAUI*(CC(I-1,2,K)-CC(I-1,3,K))
      CI3 = TAUI*(CC(I,2,K)-CC(I,3,K))
      DR2 = CR2-CI3
      DR3 = CR2+CI3
      DI2 = CI2+CR3
      DI3 = CI2-CR3
      CH(I,K,2) = WA1(I-1)*DI2+WA1(I)*DR2
      CH(I-1,K,2) = WA1(I-1)*DR2-WA1(I)*DI2
      CH(I,K,3) = WA2(I-1)*DI3+WA2(I)*DR3
      CH(I-1,K,3) = WA2(I-1)*DR3-WA2(I)*DI3
  103    continue
  104 continue
  return
  END
  subroutine PASSB4 (IDO,L1,CC,CH,WA1,WA2,WA3)
  DIMENSION       CC(IDO,4,L1)           ,CH(IDO,L1,4), &
                  WA1(1)     ,WA2(1)     ,WA3(1)
  if (IDO /= 2) goto 102
  DO 101 K=1,L1
   TI1 = CC(2,1,K)-CC(2,3,K)
   TI2 = CC(2,1,K)+CC(2,3,K)
   TR4 = CC(2,4,K)-CC(2,2,K)
   TI3 = CC(2,2,K)+CC(2,4,K)
   TR1 = CC(1,1,K)-CC(1,3,K)
   TR2 = CC(1,1,K)+CC(1,3,K)
   TI4 = CC(1,2,K)-CC(1,4,K)
   TR3 = CC(1,2,K)+CC(1,4,K)
   CH(1,K,1) = TR2+TR3
   CH(1,K,3) = TR2-TR3
   CH(2,K,1) = TI2+TI3
   CH(2,K,3) = TI2-TI3
   CH(1,K,2) = TR1+TR4
   CH(1,K,4) = TR1-TR4
   CH(2,K,2) = TI1+TI4
   CH(2,K,4) = TI1-TI4
  101 continue
  return
  102 DO 104 K=1,L1
   DO 103 I=2,IDO,2
      TI1 = CC(I,1,K)-CC(I,3,K)
      TI2 = CC(I,1,K)+CC(I,3,K)
      TI3 = CC(I,2,K)+CC(I,4,K)
      TR4 = CC(I,4,K)-CC(I,2,K)
      TR1 = CC(I-1,1,K)-CC(I-1,3,K)
      TR2 = CC(I-1,1,K)+CC(I-1,3,K)
      TI4 = CC(I-1,2,K)-CC(I-1,4,K)
      TR3 = CC(I-1,2,K)+CC(I-1,4,K)
      CH(I-1,K,1) = TR2+TR3
      CR3 = TR2-TR3
      CH(I,K,1) = TI2+TI3
      CI3 = TI2-TI3
      CR2 = TR1+TR4
      CR4 = TR1-TR4
      CI2 = TI1+TI4
      CI4 = TI1-TI4
      CH(I-1,K,2) = WA1(I-1)*CR2-WA1(I)*CI2
      CH(I,K,2) = WA1(I-1)*CI2+WA1(I)*CR2
      CH(I-1,K,3) = WA2(I-1)*CR3-WA2(I)*CI3
      CH(I,K,3) = WA2(I-1)*CI3+WA2(I)*CR3
      CH(I-1,K,4) = WA3(I-1)*CR4-WA3(I)*CI4
      CH(I,K,4) = WA3(I-1)*CI4+WA3(I)*CR4
  103    continue
  104 continue
  return
  END
  subroutine PASSB5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4)
  DIMENSION       CC(IDO,5,L1)           ,CH(IDO,L1,5), &
                  WA1(1)     ,WA2(1)     ,WA3(1)     ,WA4(1)
  DATA TR11,TI11,TR12,TI12 /.309016994374947,.951056516295154, &
  -.809016994374947,.587785252292473/
  if (IDO /= 2) goto 102
  DO 101 K=1,L1
   TI5 = CC(2,2,K)-CC(2,5,K)
   TI2 = CC(2,2,K)+CC(2,5,K)
   TI4 = CC(2,3,K)-CC(2,4,K)
   TI3 = CC(2,3,K)+CC(2,4,K)
   TR5 = CC(1,2,K)-CC(1,5,K)
   TR2 = CC(1,2,K)+CC(1,5,K)
   TR4 = CC(1,3,K)-CC(1,4,K)
   TR3 = CC(1,3,K)+CC(1,4,K)
   CH(1,K,1) = CC(1,1,K)+TR2+TR3
   CH(2,K,1) = CC(2,1,K)+TI2+TI3
   CR2 = CC(1,1,K)+TR11*TR2+TR12*TR3
   CI2 = CC(2,1,K)+TR11*TI2+TR12*TI3
   CR3 = CC(1,1,K)+TR12*TR2+TR11*TR3
   CI3 = CC(2,1,K)+TR12*TI2+TR11*TI3
   CR5 = TI11*TR5+TI12*TR4
   CI5 = TI11*TI5+TI12*TI4
   CR4 = TI12*TR5-TI11*TR4
   CI4 = TI12*TI5-TI11*TI4
   CH(1,K,2) = CR2-CI5
   CH(1,K,5) = CR2+CI5
   CH(2,K,2) = CI2+CR5
   CH(2,K,3) = CI3+CR4
   CH(1,K,3) = CR3-CI4
   CH(1,K,4) = CR3+CI4
   CH(2,K,4) = CI3-CR4
   CH(2,K,5) = CI2-CR5
  101 continue
  return
  102 DO 104 K=1,L1
   DO 103 I=2,IDO,2
      TI5 = CC(I,2,K)-CC(I,5,K)
      TI2 = CC(I,2,K)+CC(I,5,K)
      TI4 = CC(I,3,K)-CC(I,4,K)
      TI3 = CC(I,3,K)+CC(I,4,K)
      TR5 = CC(I-1,2,K)-CC(I-1,5,K)
      TR2 = CC(I-1,2,K)+CC(I-1,5,K)
      TR4 = CC(I-1,3,K)-CC(I-1,4,K)
      TR3 = CC(I-1,3,K)+CC(I-1,4,K)
      CH(I-1,K,1) = CC(I-1,1,K)+TR2+TR3
      CH(I,K,1) = CC(I,1,K)+TI2+TI3
      CR2 = CC(I-1,1,K)+TR11*TR2+TR12*TR3
      CI2 = CC(I,1,K)+TR11*TI2+TR12*TI3
      CR3 = CC(I-1,1,K)+TR12*TR2+TR11*TR3
      CI3 = CC(I,1,K)+TR12*TI2+TR11*TI3
      CR5 = TI11*TR5+TI12*TR4
      CI5 = TI11*TI5+TI12*TI4
      CR4 = TI12*TR5-TI11*TR4
      CI4 = TI12*TI5-TI11*TI4
      DR3 = CR3-CI4
      DR4 = CR3+CI4
      DI3 = CI3+CR4
      DI4 = CI3-CR4
      DR5 = CR2+CI5
      DR2 = CR2-CI5
      DI5 = CI2-CR5
      DI2 = CI2+CR5
      CH(I-1,K,2) = WA1(I-1)*DR2-WA1(I)*DI2
      CH(I,K,2) = WA1(I-1)*DI2+WA1(I)*DR2
      CH(I-1,K,3) = WA2(I-1)*DR3-WA2(I)*DI3
      CH(I,K,3) = WA2(I-1)*DI3+WA2(I)*DR3
      CH(I-1,K,4) = WA3(I-1)*DR4-WA3(I)*DI4
      CH(I,K,4) = WA3(I-1)*DI4+WA3(I)*DR4
      CH(I-1,K,5) = WA4(I-1)*DR5-WA4(I)*DI5
      CH(I,K,5) = WA4(I-1)*DI5+WA4(I)*DR5
  103    continue
  104 continue
  return
  END

  subroutine CFFTI (N,WSAVE)
  DIMENSION       WSAVE(1)
  if (N == 1) return
  IW1 = N+N+1
  IW2 = IW1+N+N
  CALL CFFTI1 (N,WSAVE(IW1),WSAVE(IW2))
  return
  END
  subroutine CFFTI1 (N,WA,IFAC)
  DIMENSION       WA(1)      ,IFAC(1)    ,NTRYH(4)
  DATA NTRYH(1),NTRYH(2),NTRYH(3),NTRYH(4)/3,4,2,5/
  NL = N
  NF = 0
  J = 0
  101 J = J+1
  if (J-4) 102,102,103
  102 NTRY = NTRYH(J)
  goto 104
  103 NTRY = NTRY+2
  104 NQ = NL/NTRY
  NR = NL-NTRY*NQ
  if (NR) 101,105,101
  105 NF = NF+1
  IFAC(NF+2) = NTRY
  NL = NQ
  if (NTRY /= 2) goto 107
  if (NF == 1) goto 107
  DO 106 I=2,NF
   IB = NF-I+2
   IFAC(IB+2) = IFAC(IB+1)
  106 continue
  IFAC(3) = 2
  107 if (NL /= 1) goto 104
  IFAC(1) = N
  IFAC(2) = NF
  TPI = 6.28318530717959
  ARGH = TPI/FLOAT(N)
  I = 2
  L1 = 1
  DO 110 K1=1,NF
   IP = IFAC(K1+2)
   LD = 0
   L2 = L1*IP
   IDO = N/L2
   IDOT = IDO+IDO+2
   IPM = IP-1
   DO 109 J=1,IPM
      I1 = I
      WA(I-1) = 1.
      WA(I) = 0.
      LD = LD+L1
      FI = 0.
      ARGLD = FLOAT(LD)*ARGH
      DO 108 II=4,IDOT,2
         I = I+2
         FI = FI+1.
         ARG = FI*ARGLD
         WA(I-1) = COS(ARG)
         WA(I) = SIN(ARG)
  108       continue
      if (IP <= 5) goto 109
      WA(I1-1) = WA(I-1)
      WA(I1) = WA(I)
  109    continue
   L1 = L2
  110 continue
  return
  END

  subroutine CFFTF (N,C,WSAVE)
  DIMENSION       C(1)       ,WSAVE(1)
  if (N == 1) return
  IW1 = N+N+1
  IW2 = IW1+N+N
  CALL CFFTF1 (N,C,WSAVE,WSAVE(IW1),WSAVE(IW2))
  return
  END
  subroutine CFFTF1 (N,C,CH,WA,IFAC)
  DIMENSION       CH(1)      ,C(1)       ,WA(1)      ,IFAC(1)
  NF = IFAC(2)
  NA = 0
  L1 = 1
  IW = 1
  DO 116 K1=1,NF
   IP = IFAC(K1+2)
   L2 = IP*L1
   IDO = N/L2
   IDOT = IDO+IDO
   IDL1 = IDOT*L1
   if (IP /= 4) goto 103
   IX2 = IW+IDOT
   IX3 = IX2+IDOT
   if (NA /= 0) goto 101
   CALL PASSF4 (IDOT,L1,C,CH,WA(IW),WA(IX2),WA(IX3))
   goto 102
  101    CALL PASSF4 (IDOT,L1,CH,C,WA(IW),WA(IX2),WA(IX3))
  102    NA = 1-NA
   goto 115
  103    if (IP /= 2) goto 106
   if (NA /= 0) goto 104
   CALL PASSF2 (IDOT,L1,C,CH,WA(IW))
   goto 105
  104    CALL PASSF2 (IDOT,L1,CH,C,WA(IW))
  105    NA = 1-NA
   goto 115
  106    if (IP /= 3) goto 109
   IX2 = IW+IDOT
   if (NA /= 0) goto 107
   CALL PASSF3 (IDOT,L1,C,CH,WA(IW),WA(IX2))
   goto 108
  107    CALL PASSF3 (IDOT,L1,CH,C,WA(IW),WA(IX2))
  108    NA = 1-NA
   goto 115
  109    if (IP /= 5) goto 112
   IX2 = IW+IDOT
   IX3 = IX2+IDOT
   IX4 = IX3+IDOT
   if (NA /= 0) goto 110
   CALL PASSF5 (IDOT,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4))
   goto 111
  110    CALL PASSF5 (IDOT,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4))
  111    NA = 1-NA
   goto 115
  112    if (NA /= 0) goto 113
   CALL PASSF (NAC,IDOT,IP,L1,IDL1,C,C,C,CH,CH,WA(IW))
   goto 114
  113    CALL PASSF (NAC,IDOT,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW))
  114    if (NAC /= 0) NA = 1-NA
  115    L1 = L2
   IW = IW+(IP-1)*IDOT
  116 continue
  if (NA == 0) return
  N2 = N+N
  DO 117 I=1,N2
   C(I) = CH(I)
  117 continue
  return
  END
  subroutine PASSF (NAC,IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA)
  DIMENSION       CH(IDO,L1,IP)          ,CC(IDO,IP,L1), &
                  C1(IDO,L1,IP)          ,WA(1)      ,C2(IDL1,IP), &
                  CH2(IDL1,IP)
  IDOT = IDO/2
  NT = IP*IDL1
  IPP2 = IP+2
  IPPH = (IP+1)/2
  IDP = IP*IDO
!
  if (IDO < L1) goto 106
  DO 103 J=2,IPPH
   JC = IPP2-J
   DO 102 K=1,L1
      DO 101 I=1,IDO
         CH(I,K,J) = CC(I,J,K)+CC(I,JC,K)
         CH(I,K,JC) = CC(I,J,K)-CC(I,JC,K)
  101       continue
  102    continue
  103 continue
  DO 105 K=1,L1
   DO 104 I=1,IDO
      CH(I,K,1) = CC(I,1,K)
  104    continue
  105 continue
  goto 112
  106 DO 109 J=2,IPPH
   JC = IPP2-J
   DO 108 I=1,IDO
      DO 107 K=1,L1
         CH(I,K,J) = CC(I,J,K)+CC(I,JC,K)
         CH(I,K,JC) = CC(I,J,K)-CC(I,JC,K)
  107       continue
  108    continue
  109 continue
  DO 111 I=1,IDO
   DO 110 K=1,L1
      CH(I,K,1) = CC(I,1,K)
  110    continue
  111 continue
  112 IDL = 2-IDO
  INC = 0
  DO 116 L=2,IPPH
   LC = IPP2-L
   IDL = IDL+IDO
   DO 113 IK=1,IDL1
      C2(IK,L) = CH2(IK,1)+WA(IDL-1)*CH2(IK,2)
      C2(IK,LC) = -WA(IDL)*CH2(IK,IP)
  113    continue
   IDLJ = IDL
   INC = INC+IDO
   DO 115 J=3,IPPH
      JC = IPP2-J
      IDLJ = IDLJ+INC
      if (IDLJ > IDP) IDLJ = IDLJ-IDP
      WAR = WA(IDLJ-1)
      WAI = WA(IDLJ)
      DO 114 IK=1,IDL1
         C2(IK,L) = C2(IK,L)+WAR*CH2(IK,J)
         C2(IK,LC) = C2(IK,LC)-WAI*CH2(IK,JC)
  114       continue
  115    continue
  116 continue
  DO 118 J=2,IPPH
   DO 117 IK=1,IDL1
      CH2(IK,1) = CH2(IK,1)+CH2(IK,J)
  117    continue
  118 continue
  DO 120 J=2,IPPH
   JC = IPP2-J
   DO 119 IK=2,IDL1,2
      CH2(IK-1,J) = C2(IK-1,J)-C2(IK,JC)
      CH2(IK-1,JC) = C2(IK-1,J)+C2(IK,JC)
      CH2(IK,J) = C2(IK,J)+C2(IK-1,JC)
      CH2(IK,JC) = C2(IK,J)-C2(IK-1,JC)
  119    continue
  120 continue
  NAC = 1
  if (IDO == 2) return
  NAC = 0
  DO 121 IK=1,IDL1
   C2(IK,1) = CH2(IK,1)
  121 continue
  DO 123 J=2,IP
   DO 122 K=1,L1
      C1(1,K,J) = CH(1,K,J)
      C1(2,K,J) = CH(2,K,J)
  122    continue
  123 continue
  if (IDOT > L1) goto 127
  IDIJ = 0
  DO 126 J=2,IP
   IDIJ = IDIJ+2
   DO 125 I=4,IDO,2
      IDIJ = IDIJ+2
      DO 124 K=1,L1
         C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)+WA(IDIJ)*CH(I,K,J)
         C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)-WA(IDIJ)*CH(I-1,K,J)
  124       continue
  125    continue
  126 continue
  return
  127 IDJ = 2-IDO
  DO 130 J=2,IP
   IDJ = IDJ+IDO
   DO 129 K=1,L1
      IDIJ = IDJ
      DO 128 I=4,IDO,2
         IDIJ = IDIJ+2
         C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)+WA(IDIJ)*CH(I,K,J)
         C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)-WA(IDIJ)*CH(I-1,K,J)
  128       continue
  129    continue
  130 continue
  return
  END
  subroutine PASSF2 (IDO,L1,CC,CH,WA1)
  DIMENSION       CC(IDO,2,L1)           ,CH(IDO,L1,2), &
                  WA1(1)
  if (IDO > 2) goto 102
  DO 101 K=1,L1
   CH(1,K,1) = CC(1,1,K)+CC(1,2,K)
   CH(1,K,2) = CC(1,1,K)-CC(1,2,K)
   CH(2,K,1) = CC(2,1,K)+CC(2,2,K)
   CH(2,K,2) = CC(2,1,K)-CC(2,2,K)
  101 continue
  return
  102 DO 104 K=1,L1
   DO 103 I=2,IDO,2
      CH(I-1,K,1) = CC(I-1,1,K)+CC(I-1,2,K)
      TR2 = CC(I-1,1,K)-CC(I-1,2,K)
      CH(I,K,1) = CC(I,1,K)+CC(I,2,K)
      TI2 = CC(I,1,K)-CC(I,2,K)
      CH(I,K,2) = WA1(I-1)*TI2-WA1(I)*TR2
      CH(I-1,K,2) = WA1(I-1)*TR2+WA1(I)*TI2
  103    continue
  104 continue
  return
  END
  subroutine PASSF3 (IDO,L1,CC,CH,WA1,WA2)
  DIMENSION       CC(IDO,3,L1)           ,CH(IDO,L1,3), &
                  WA1(1)     ,WA2(1)
  DATA TAUR,TAUI /-.5,-.866025403784439/
  if (IDO /= 2) goto 102
  DO 101 K=1,L1
   TR2 = CC(1,2,K)+CC(1,3,K)
   CR2 = CC(1,1,K)+TAUR*TR2
   CH(1,K,1) = CC(1,1,K)+TR2
   TI2 = CC(2,2,K)+CC(2,3,K)
   CI2 = CC(2,1,K)+TAUR*TI2
   CH(2,K,1) = CC(2,1,K)+TI2
   CR3 = TAUI*(CC(1,2,K)-CC(1,3,K))
   CI3 = TAUI*(CC(2,2,K)-CC(2,3,K))
   CH(1,K,2) = CR2-CI3
   CH(1,K,3) = CR2+CI3
   CH(2,K,2) = CI2+CR3
   CH(2,K,3) = CI2-CR3
  101 continue
  return
  102 DO 104 K=1,L1
   DO 103 I=2,IDO,2
      TR2 = CC(I-1,2,K)+CC(I-1,3,K)
      CR2 = CC(I-1,1,K)+TAUR*TR2
      CH(I-1,K,1) = CC(I-1,1,K)+TR2
      TI2 = CC(I,2,K)+CC(I,3,K)
      CI2 = CC(I,1,K)+TAUR*TI2
      CH(I,K,1) = CC(I,1,K)+TI2
      CR3 = TAUI*(CC(I-1,2,K)-CC(I-1,3,K))
      CI3 = TAUI*(CC(I,2,K)-CC(I,3,K))
      DR2 = CR2-CI3
      DR3 = CR2+CI3
      DI2 = CI2+CR3
      DI3 = CI2-CR3
      CH(I,K,2) = WA1(I-1)*DI2-WA1(I)*DR2
      CH(I-1,K,2) = WA1(I-1)*DR2+WA1(I)*DI2
      CH(I,K,3) = WA2(I-1)*DI3-WA2(I)*DR3
      CH(I-1,K,3) = WA2(I-1)*DR3+WA2(I)*DI3
  103    continue
  104 continue
  return
  END
  subroutine PASSF4 (IDO,L1,CC,CH,WA1,WA2,WA3)
  DIMENSION       CC(IDO,4,L1)           ,CH(IDO,L1,4), &
                  WA1(1)     ,WA2(1)     ,WA3(1)
  if (IDO /= 2) goto 102
  DO 101 K=1,L1
   TI1 = CC(2,1,K)-CC(2,3,K)
   TI2 = CC(2,1,K)+CC(2,3,K)
   TR4 = CC(2,2,K)-CC(2,4,K)
   TI3 = CC(2,2,K)+CC(2,4,K)
   TR1 = CC(1,1,K)-CC(1,3,K)
   TR2 = CC(1,1,K)+CC(1,3,K)
   TI4 = CC(1,4,K)-CC(1,2,K)
   TR3 = CC(1,2,K)+CC(1,4,K)
   CH(1,K,1) = TR2+TR3
   CH(1,K,3) = TR2-TR3
   CH(2,K,1) = TI2+TI3
   CH(2,K,3) = TI2-TI3
   CH(1,K,2) = TR1+TR4
   CH(1,K,4) = TR1-TR4
   CH(2,K,2) = TI1+TI4
   CH(2,K,4) = TI1-TI4
  101 continue
  return
  102 DO 104 K=1,L1
   DO 103 I=2,IDO,2
      TI1 = CC(I,1,K)-CC(I,3,K)
      TI2 = CC(I,1,K)+CC(I,3,K)
      TI3 = CC(I,2,K)+CC(I,4,K)
      TR4 = CC(I,2,K)-CC(I,4,K)
      TR1 = CC(I-1,1,K)-CC(I-1,3,K)
      TR2 = CC(I-1,1,K)+CC(I-1,3,K)
      TI4 = CC(I-1,4,K)-CC(I-1,2,K)
      TR3 = CC(I-1,2,K)+CC(I-1,4,K)
      CH(I-1,K,1) = TR2+TR3
      CR3 = TR2-TR3
      CH(I,K,1) = TI2+TI3
      CI3 = TI2-TI3
      CR2 = TR1+TR4
      CR4 = TR1-TR4
      CI2 = TI1+TI4
      CI4 = TI1-TI4
      CH(I-1,K,2) = WA1(I-1)*CR2+WA1(I)*CI2
      CH(I,K,2) = WA1(I-1)*CI2-WA1(I)*CR2
      CH(I-1,K,3) = WA2(I-1)*CR3+WA2(I)*CI3
      CH(I,K,3) = WA2(I-1)*CI3-WA2(I)*CR3
      CH(I-1,K,4) = WA3(I-1)*CR4+WA3(I)*CI4
      CH(I,K,4) = WA3(I-1)*CI4-WA3(I)*CR4
  103    continue
  104 continue
  return
  END
  subroutine PASSF5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4)
  DIMENSION       CC(IDO,5,L1)           ,CH(IDO,L1,5), &
                  WA1(1)     ,WA2(1)     ,WA3(1)     ,WA4(1)
  DATA TR11,TI11,TR12,TI12 /.309016994374947,-.951056516295154, &
  -.809016994374947,-.587785252292473/
  if (IDO /= 2) goto 102
  DO 101 K=1,L1
   TI5 = CC(2,2,K)-CC(2,5,K)
   TI2 = CC(2,2,K)+CC(2,5,K)
   TI4 = CC(2,3,K)-CC(2,4,K)
   TI3 = CC(2,3,K)+CC(2,4,K)
   TR5 = CC(1,2,K)-CC(1,5,K)
   TR2 = CC(1,2,K)+CC(1,5,K)
   TR4 = CC(1,3,K)-CC(1,4,K)
   TR3 = CC(1,3,K)+CC(1,4,K)
   CH(1,K,1) = CC(1,1,K)+TR2+TR3
   CH(2,K,1) = CC(2,1,K)+TI2+TI3
   CR2 = CC(1,1,K)+TR11*TR2+TR12*TR3
   CI2 = CC(2,1,K)+TR11*TI2+TR12*TI3
   CR3 = CC(1,1,K)+TR12*TR2+TR11*TR3
   CI3 = CC(2,1,K)+TR12*TI2+TR11*TI3
   CR5 = TI11*TR5+TI12*TR4
   CI5 = TI11*TI5+TI12*TI4
   CR4 = TI12*TR5-TI11*TR4
   CI4 = TI12*TI5-TI11*TI4
   CH(1,K,2) = CR2-CI5
   CH(1,K,5) = CR2+CI5
   CH(2,K,2) = CI2+CR5
   CH(2,K,3) = CI3+CR4
   CH(1,K,3) = CR3-CI4
   CH(1,K,4) = CR3+CI4
   CH(2,K,4) = CI3-CR4
   CH(2,K,5) = CI2-CR5
  101 continue
  return
  102 DO 104 K=1,L1
   DO 103 I=2,IDO,2
      TI5 = CC(I,2,K)-CC(I,5,K)
      TI2 = CC(I,2,K)+CC(I,5,K)
      TI4 = CC(I,3,K)-CC(I,4,K)
      TI3 = CC(I,3,K)+CC(I,4,K)
      TR5 = CC(I-1,2,K)-CC(I-1,5,K)
      TR2 = CC(I-1,2,K)+CC(I-1,5,K)
      TR4 = CC(I-1,3,K)-CC(I-1,4,K)
      TR3 = CC(I-1,3,K)+CC(I-1,4,K)
      CH(I-1,K,1) = CC(I-1,1,K)+TR2+TR3
      CH(I,K,1) = CC(I,1,K)+TI2+TI3
      CR2 = CC(I-1,1,K)+TR11*TR2+TR12*TR3
      CI2 = CC(I,1,K)+TR11*TI2+TR12*TI3
      CR3 = CC(I-1,1,K)+TR12*TR2+TR11*TR3
      CI3 = CC(I,1,K)+TR12*TI2+TR11*TI3
      CR5 = TI11*TR5+TI12*TR4
      CI5 = TI11*TI5+TI12*TI4
      CR4 = TI12*TR5-TI11*TR4
      CI4 = TI12*TI5-TI11*TI4
      DR3 = CR3-CI4
      DR4 = CR3+CI4
      DI3 = CI3+CR4
      DI4 = CI3-CR4
      DR5 = CR2+CI5
      DR2 = CR2-CI5
      DI5 = CI2-CR5
      DI2 = CI2+CR5
      CH(I-1,K,2) = WA1(I-1)*DR2+WA1(I)*DI2
      CH(I,K,2) = WA1(I-1)*DI2-WA1(I)*DR2
      CH(I-1,K,3) = WA2(I-1)*DR3+WA2(I)*DI3
      CH(I,K,3) = WA2(I-1)*DI3-WA2(I)*DR3
      CH(I-1,K,4) = WA3(I-1)*DR4+WA3(I)*DI4
      CH(I,K,4) = WA3(I-1)*DI4-WA3(I)*DR4
      CH(I-1,K,5) = WA4(I-1)*DR5+WA4(I)*DI5
      CH(I,K,5) = WA4(I-1)*DI5-WA4(I)*DR5
  103    continue
  104 continue
  return
  END


text:
!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
!> \brief Dynamic Load Balancing

subroutine hecmw_transfer_data_f2c (hecMESHnew, adapRES)

  use  hecmw_util
  use  hecmw_io
  use  hecmw_dist_copy_f2c_f
  use  hecmw_result
  type (hecmwST_local_mesh) :: hecMESHnew
  type (hecmwST_result_data):: adapRES

  call hecmw_dlb_f2c_init
  call hecmw_dist_copy_f2c(hecMESHnew, ierr)
  call hecmw_dlb_f2c_finalize
  call test_mesh
  call hecmw_dist_result_copy_f2c(adapRES)
end subroutine hecmw_transfer_data_f2c

max_stars_repo_name: masae-hayashi/FrontISTR

max_issues_repo_name: masae-hayashi/FrontISTR

max_forks_repo_name: masae-hayashi/FrontISTR


text:
The McKeon Condominiums, also known as the La Buena Vida condominiums, are in a circular layout next to the Wildhorse development in East Davis.  There are four units per building, each unit gets one half of a twocar garage, and all four share a coinoperated washerdryer set.  Each unit is separately owned, and the condominium association maintains the building exterior and common areas. At the center of the circular layout is a shared pool and common room that can be reserved by residents.  These condos are generally among the least expensive housing in Davis. The majority of the condos are rented out, but there arent many students in the complex. The most common occupants seem to be families with young children and retirees, and the area is generally quiet and well maintained.

The La Buena Vida complex has special mention in Davis Municipal Code/39 due to the way it meters water usage.

Each building has a groundfloor flat (#1), overthegarage flat (#4), and two mirrorimage twostory units on the side.  Most reports indicate the #1 unit is most desirable, and the #4 is least desirable (according to selling price).  Note that they are not part of the new Wildhorse construction... they were constructed in the early 1970s, and until a few years ago were surrounded by open fields on three sides.  The new construction means the greenbelt, Nugget Fields, and Sandy Motley Park join the Oak Grove Park and nearby Oak Tree Plaza making this a very pleasant living space.

On the down side, the list of Codes, Covenants and Restrictions for this community is rather thick, so dont plan on doing things like changing your yard, or hanging clotheslines, and only the #4 unit has any significant porch.  All common areas (i.e. everything not inside a unit) are managed by the La Buena Vida Homeowners Association, which currently employs http://www.assocmc.com Association Management Concepts, Inc. for onsite management and fee collection. A 5person board of directors is elected annually from condo owners.

One rule that does seem to stay bent is the requirement to lock the laundry room in each unit... most are left unlocked to allow you to access more than one washer/dryer at a time to shorten laundry day, in spite of regular warnings from the management that this invites vandalism.
  If this was ever true, it no longer is. The current board makes every attempt to keep the locks functional and the doors secure.  Users/MattCzarnowski mc
    Does this mean someone goes around locking all the doors? Maybe my unit (which is the same as the original authors) is too far of a walk for them.  Users/AndrewChen

Unitrans L Bus Line stops near the center of the development in Fremont Court. Parking in the driveways is limited, though there is usually parking available on the street.

Computer and Internet Highspeed Internet is available but due to the distance from the central office, DSL speeds are limited to 384Kbps up/down. Comcast cable TV and Internet are available. Satellite TV is an option but the homeowners association must approve the installation of a dish. Finding a clear view of the south sky may be a problem as well.

Note: The units were built with wiki:wikipedia:Aluminium_wire aluminum wiring. When replacing electical sockets/switches or doing any other electrical work, you must use the proper CO/AL or CO/ALRrated parts unless your unit has already been retrofitted with copper connections.

20050823 16:31:03 nbsp While some things require board approval, its really not difficult to acquire. And although the official DSL speed is 384Kbps, I routinely get higher speeds. Users/MattCzarnowski

20060613 15:04:10 nbsp #4 is the least desirable unit of all 4 units in the complex and it is well known to the neighborhood and sold or rented with significantly lower price. Users/MoMostafavi

20060731 15:14:00 nbsp Although the selling prices of the #4 units seem to support the idea that they are less desireable, I find the layout of the #4 the best. Thats why I live in one. :) Users/MattCzarnowski

20070501 17:53:16 nbsp well i just moved in here and there are quite a bit of younger people who live here (probably because of the price). I also find that the #4 units have the best floorplan but the worst insulation. i dont mind because the 2 & 3 have only 1 bathroom and its upstairs.  Users/GregFurnish

20071201 12:52:50 nbsp I live in a #4 unit, and it is quite nice. It helped that the owner was able to install doublepaned windows before we moved in. I really like living in this complex; the proximity to Nugget, greenbelt, and empty space make for a convenient and comfortable place to live. However, BJ (the current home owners assoc. minion) is a major pain. In my opinion HOA should be banned on principle because they invariably lead to petty disagreements over little things.  Users/DylanBeaudette

20080413 09:53:55 nbsp Here is a fun experiment: leave your trash cans out for a couple extra days, and they will magically disappear. (Technically speaking, they dont really disappear rather the neighborhood helper BJ whisks them away to the dumpster dungeon over by the pool. You wont get a note or warning, and will have to navigate through her complicated system of message machines and pager if you want them back. ) Users/DylanBeaudette

20090408 22:02:15 nbsp I live here with my boyfriend and we have cable internet and its good. We live in a #1 unit and it is very spacious. 2 Bedrooms/ 1 Bath for $995 a month. We have to pay a set fee of $60 for water, sewage and garbage every month and every couple months they are suppose to send us another bill with the additional charges that the $60 didnt cover. We just got our bill and it was about $100 over 7 months, which isnt too bad.  Users/mowtov

20100305 05:11:25 nbsp Ya, nobody likes rules but BJ has greatly keep this complex from turning into a getto!  Users/Williams

20100912 18:44:57 nbsp The walls are horribly thin and my unit has singlepane windows. These factors combined means I can make out most conversations going on in other units and outside. If you value privacy or quiet, this is NOT the place for you! Users/JenniferF

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
Describe Users/PamelaTrokanski here.

Welcome to the wiki Pamela, its good to see you providing more info! Daubert
One of my buddies really enjoyed his class with you, I had fun as well


text:
cdis
cdis    Open Source License/Disclaimer, Forecast Systems Laboratory
cdis    NOAA/OAR/FSL, 325 Broadway Boulder, CO 80305
cdis
cdis    This software is distributed under the Open Source Definition,
cdis    which may be found at http://www.opensource.org/osd.html.
cdis
cdis    In particular, redistribution and use in source and binary forms,
cdis    with or without modification, are permitted provided that the
cdis    following conditions are met:
cdis
cdis    - Redistributions of source code must retain this notice, this
cdis    list of conditions and the following disclaimer.
cdis
cdis    - Redistributions in binary form must provide access to this
cdis    notice, this list of conditions and the following disclaimer, and
cdis    the underlying source code.
cdis
cdis    - All modifications to this software must be clearly documented,
cdis    and are solely the responsibility of the agent making the
cdis    modifications.
cdis
cdis    - If significant modifications or enhancements are made to this
cdis    software, the FSL Software Policy Manager
cdis    (softwaremgr@fsl.noaa.gov) should be notified.
cdis
cdis    THIS SOFTWARE AND ITS DOCUMENTATION ARE IN THE PUBLIC DOMAIN
cdis    AND ARE FURNISHED "AS IS."  THE AUTHORS, THE UNITED STATES
cdis    GOVERNMENT, ITS INSTRUMENTALITIES, OFFICERS, EMPLOYEES, AND
cdis    AGENTS MAKE NO WARRANTY, EXPRESS OR IMPLIED, AS TO THE USEFULNESS
cdis    OF THE SOFTWARE AND DOCUMENTATION FOR ANY PURPOSE.  THEY ASSUME
cdis    NO RESPONSIBILITY (1) FOR THE USE OF THE SOFTWARE AND
cdis    DOCUMENTATION; OR (2) TO PROVIDE TECHNICAL SUPPORT TO USERS.
cdis
cdis cdis
cdis
cdis
cdis
cdis
cdis
cdis
cdis
cdis
cdis
	Subroutine Precip_Water_Profile(Q,P,ZA,N,Pw)
C
C************************************************************************
C*									*
C*	Module Name:	Precip_Water_Profile				*
C*									*
C*	Language:	Fortran-77	   Library:			*
C*	Version.Rev:	1.0 30 April 96	Programmer: Kleespies		*
C*									*
C*	Calling Seq:	Call Precip_Water_Profile(Q,P,Za,N,Pw)		*
C*									*
C*	Description:	Computes precipitable water profilefrom a mixing*
C*			ratio profile and a pressure profile, along	*
C*			a slant path.					*
C*									*
C*			This code largely lifted from RADCOM.		*
C*									*
C*	Input Args:	R*4 	Q(N)	Mixing ratio profile		*
C*			R*4	P(N)	Pressure profile		*
C*			R*4	ZA	Zenith Angle			*
C*			R*4	N	Number of points in profile	*
C*									*
C*	Output Args:	R*4	Precipitable Water profile		*
C*									*
C*	Common Blks:	None						*
C*									*
C*	Include:	NOne						*
C*									*
C*	Externals:	None						*
C*									*
C*	Data Files:	None						*
C*									*
C*	Restrictions:	Q,P,Pw dimensioned in calling routine.		*
C*									*
C*	Error Codes:	None						*
C*									*
C*	Error Messages:	None						*
C*									*
C************************************************************************
C
	Implicit NONE

*	Input
	Integer N
	Real Q(N),P(N),ZA  ! modified by Dan Birkenheuer 8/10/98

*	Output
	Real Pw(N)
	Real Path, CGrav
	Data CGrav / 5.098581e-4 / ! .5/980.665
	Real Delp

	Integer i	! local variables

	Path = CGrav / Cos(ZA*acos(-1.)/180.)

	Pw(1) = 0.0

	Do i = 2 , N
	 Delp = Abs(P(i) - P(i-1))
	 Pw(i) = Pw(i-1) + Path*(Q(i) + Q(i-1))*Delp
	EndDo

	Return
	End

max_stars_repo_name: maxinye/laps-mirror

max_stars_repo_licenses: Intel Unlicense OLDAP-2.2.1 NetCDF

max_issues_repo_name: longwosion/laps-mirror

max_issues_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense

max_forks_repo_name: longwosion/laps-mirror

max_forks_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense


text:
      SUBROUTINE EEVMM (VWWNLM, PJVPLM, VWMPMT)

C note: no 3D version (eevmm3) because standard does not specify
C a unique transformation.

      REAL   VWMPMT(3,3), VWWNLM(4), PJVPLM(4)

      CALL EBLTM (VWWNLM(1),VWWNLM(3),
     1             PJVPLM(1)-VWWNLM(1), PJVPLM(3)-VWWNLM(3), 0.0,
     2            (PJVPLM(2)-PJVPLM(1))/(VWWNLM(2)-VWWNLM(1)),
     3            (PJVPLM(4)-PJVPLM(3))/(VWWNLM(4)-VWWNLM(3)),  VWMPMT)

      END


text:
!!  Copyright (C)  Stichting Deltares, 2012-2016.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License version 3,
!!  as published by the Free Software Foundation.
!!
!!  This program is distributed in the hope that it will be useful,
!!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!!  GNU General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!  contact: delft3d.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands
!!
!!  All indications and logos of, and references to registered trademarks
!!  of Stichting Deltares remain the property of Stichting Deltares. All
!!  rights reserved.

      SUBROUTINE RD_FILID ( DEFFDS,         FFORM , VFFORM, CONTEN,
     +                      VERSIO, SERIAL, RUNDAT, SOURCE, REMARK,
     +                      LUNREP, IERROR)
!
!     Deltares
!
!     CREATED            :  june 1999 by Jan van Beek
!
!     FUNCTION           :  Read File Identification group from NEFIS file
!
!     FILES              :  NEFIS file assumed opened
!
!     SUBROUTINES CALLED :
!
!     ARGUMENTS
!
!     NAME    TYPE     LENGTH     FUNCT.  DESCRIPTION
!     ----    -----    ------     ------- -----------
!     DEFFDS  INT      2993       I/O     Definition file descriptor
!     DATFDS  INT       999       I/O     Data file descriptor
!     FFORM   CH*40       1       O       File Format
!     VFFORM  REAL        1       O       Version File Format
!     CONTEN  CH*40       1       O       File contents
!     VERSIO  REAL        1       O       Version number process library
!     SERIAL  INT         1       O       File serial number
!     RUNDAT  CH*20       1       O       Creation date
!     SOURCE  CH*40       1       O       Source data
!     REMARK  CH*40       4       O       Remarks
!     LUNREP  INT         1       I       Unit number report file
!     IERROR  INT         1       O       Error
!
!     IMPLICIT NONE for extra compiler checks
!     SAVE to keep the group definition intact
!
      IMPLICIT NONE
      SAVE
!
!     declaration of arguments
!
      INTEGER       SERIAL      , LUNREP     ,
     +              IERROR
      INTEGER       DEFFDS
      REAL          VFFORM      , VERSIO
      CHARACTER*20  RUNDAT
      CHARACTER*40  FFORM       , CONTEN     ,
     +              SOURCE
      CHARACTER*40  REMARK(4)
!
!     Local variables
!
!     GRPNAM  CHAR*16     1       LOCAL   group name (table)
!     NELEMS  INTEGER     1       LOCAL   number of elements in group (=cell)
!     ELMNMS  CHAR*16  NELEMS     LOCAL   name of elements on file
!     ELMTPS  CHAR*16  NELEMS     LOCAL   type of elements
!     ELMDMS  INTEGER  6,NELEMS   LOCAL   dimension of elements
!     NBYTSG  INTEGER  NELEMS     LOCAL   length of elements (bytes)
!
      INTEGER       NELEMS
      PARAMETER   ( NELEMS = 8 )
!
      INTEGER       I               , IELM          ,
     +              BUFLEN
      INTEGER       ELMDMS(2,NELEMS), NBYTSG(NELEMS),
     +              UINDEX(3)
      CHARACTER*16  GRPNAM
      CHARACTER*16  ELMNMS(NELEMS)  , ELMTPS(NELEMS)
      CHARACTER*64  ELMDES(NELEMS)
!
!     External NEFIS Functions
!
      INTEGER   GETELS
     +         ,GETELT
      EXTERNAL  GETELS
     +         ,GETELT
!
!     element names
!
      DATA  GRPNAM  /'FILE_ID'/
      DATA
     + (ELMNMS(I),ELMTPS(I),NBYTSG(I),ELMDMS(1,I),ELMDMS(2,I),ELMDES(I),
     +  I = 1 , NELEMS)
     +/'FFORM ', 'CHARACTER', 40, 1, 1,'File Format                   ',
     + 'VFFORM', 'REAL'     ,  4, 1, 1,'Version File Format           ',
     + 'CONTEN', 'CHARACTER', 40, 1, 1,'File contents                 ',
     + 'VERSIO', 'REAL'     ,  4, 1, 1,'Version number process library',
     + 'SERIAL', 'INTEGER'  ,  4, 1, 1,'File serial number            ',
     + 'RUNDAT', 'CHARACTER', 20, 1, 1,'Creation date                 ',
     + 'SOURCE', 'CHARACTER', 40, 1, 1,'Source data                   ',
     + 'REMARK', 'CHARACTER', 40, 1, 4,'Remarks                       '/
!
!     Read group
!
!     WRITE(LUNREP,*) ' reading GROUP:',GRPNAM
      UINDEX(1) = 1
      UINDEX(2) = 1
      UINDEX(3) = 1
!     WRITE(LUNREP,*) ' reading ELEMENT:',ELMNMS(1)
      BUFLEN = NBYTSG(1)*ELMDMS(2,1)
      IERROR = GETELS (DEFFDS ,
     +                 GRPNAM , ELMNMS(1),
     +                 UINDEX , 1        ,
     +                 BUFLEN , FFORM    )
      IF ( IERROR .NE. 0 ) THEN
         WRITE(LUNREP,*) 'ERROR reading element',ELMNMS(1)
         WRITE(LUNREP,*) 'ERROR number:',IERROR
         GOTO 900
      ENDIF
!     WRITE(LUNREP,*) ' reading ELEMENT:',ELMNMS(2)
      BUFLEN = NBYTSG(2)*ELMDMS(2,2)
      IERROR = GETELT (DEFFDS ,
     +                 GRPNAM , ELMNMS(2),
     +                 UINDEX , 1        ,
     +                 BUFLEN , VFFORM   )
      IF ( IERROR .NE. 0 ) THEN
         WRITE(LUNREP,*) 'ERROR reading element',ELMNMS(2)
         WRITE(LUNREP,*) 'ERROR number:',IERROR
         GOTO 900
      ENDIF
!     WRITE(LUNREP,*) ' reading ELEMENT:',ELMNMS(3)
      BUFLEN = NBYTSG(3)*ELMDMS(2,3)
      IERROR = GETELS (DEFFDS ,
     +                 GRPNAM , ELMNMS(3),
     +                 UINDEX , 1        ,
     +                 BUFLEN , CONTEN   )
      IF ( IERROR .NE. 0 ) THEN
         WRITE(LUNREP,*) 'ERROR reading element',ELMNMS(3)
         WRITE(LUNREP,*) 'ERROR number:',IERROR
         GOTO 900
      ENDIF
!     WRITE(LUNREP,*) ' reading ELEMENT:',ELMNMS(4)
      BUFLEN = NBYTSG(4)*ELMDMS(2,4)
      IERROR = GETELT (DEFFDS ,
     +                 GRPNAM , ELMNMS(4),
     +                 UINDEX , 1        ,
     +                 BUFLEN , VERSIO   )
      IF ( IERROR .NE. 0 ) THEN
         WRITE(LUNREP,*) 'ERROR reading element',ELMNMS(4)
         WRITE(LUNREP,*) 'ERROR number:',IERROR
         GOTO 900
      ENDIF
!     WRITE(LUNREP,*) ' reading ELEMENT:',ELMNMS(5)
      BUFLEN = NBYTSG(5)*ELMDMS(2,5)
      IERROR = GETELT (DEFFDS ,
     +                 GRPNAM , ELMNMS(5),
     +                 UINDEX , 1        ,
     +                 BUFLEN , SERIAL   )
      IF ( IERROR .NE. 0 ) THEN
         WRITE(LUNREP,*) 'ERROR reading element',ELMNMS(5)
         WRITE(LUNREP,*) 'ERROR number:',IERROR
         GOTO 900
      ENDIF
!     WRITE(LUNREP,*) ' reading ELEMENT:',ELMNMS(6)
      BUFLEN = NBYTSG(6)*ELMDMS(2,6)
      IERROR = GETELS (DEFFDS ,
     +                 GRPNAM , ELMNMS(6),
     +                 UINDEX , 1        ,
     +                 BUFLEN , RUNDAT   )
      IF ( IERROR .NE. 0 ) THEN
         WRITE(LUNREP,*) 'ERROR reading element',ELMNMS(6)
         WRITE(LUNREP,*) 'ERROR number:',IERROR
         GOTO 900
      ENDIF
!     WRITE(LUNREP,*) ' reading ELEMENT:',ELMNMS(7)
      BUFLEN = NBYTSG(7)*ELMDMS(2,7)
      IERROR = GETELS (DEFFDS ,
     +                 GRPNAM , ELMNMS(7),
     +                 UINDEX , 1        ,
     +                 BUFLEN , SOURCE   )
      IF ( IERROR .NE. 0 ) THEN
         WRITE(LUNREP,*) 'ERROR reading element',ELMNMS(7)
         WRITE(LUNREP,*) 'ERROR number:',IERROR
         GOTO 900
      ENDIF
!     WRITE(LUNREP,*) ' reading ELEMENT:',ELMNMS(8)
      BUFLEN = NBYTSG(8)*ELMDMS(2,8)
      IERROR = GETELS (DEFFDS ,
     +                 GRPNAM , ELMNMS(8),
     +                 UINDEX , 1        ,
     +                 BUFLEN , REMARK   )
      IF ( IERROR .NE. 0 ) THEN
         WRITE(LUNREP,*) 'ERROR reading element',ELMNMS(8)
         WRITE(LUNREP,*) 'ERROR number:',IERROR
         GOTO 900
      ENDIF
!
  900 CONTINUE
      RETURN
!
      END

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
!***********************************************************************
!                                                                      *
      REAL(KIND(0.0D0)) FUNCTION CLRX (KAPPAA, K, KAPPAB)
!                                                                      *
!   The value of CLRX is the 3-j symbol:                               *
!                                                                      *
!                    ( JA        K        JB  )                        *
!                    ( 1/2       0       -1/2 )                        *
!                                                                      *
!   The  K'S are kappa angular quantum numbers. The formula is taken   *
!   from D M Brink and G R Satchler, <Angular Momentum>, second edi-   *
!   tion (Oxford: Clarendon press, 1968), p 138.   The logarithms of   *
!   the first  MFACT  factorials must be available in  COMMON/FACTS/   *
!   for this program to function correctly. Note that  N!  is stored   *
!   in FACT(N+1)                                                       *
!                                                                      *
!   No subroutines called.                                             *
!                                                                      *
!   Written by Farid A Parpia, at Oxford   Last updated: 06 Oct 1992   *
!                                                                      *
!***********************************************************************
!...Translated by Pacific-Sierra Research 77to90  4.3E  10:46:52   2/14/04
!...Modified by Charlotte Froese Fischer
!                     Gediminas Gaigalas  10/05/17
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  DOUBLE
      USE FACTS_C
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER, INTENT(IN) :: KAPPAA, K, KAPPAB
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: KA, KB, KAPKB, KABKP, KAMKB, KBMKA
      REAL(DOUBLE) :: EXPTRM
!-----------------------------------------------
!
!
!   Determine the absolute values of the kappas
!
      KA = ABS(KAPPAA)
      KB = ABS(KAPPAB)
!
!   Perform the triangularity check
!
      IF (ABS(KA - KB)<=K .AND. KA+KB-1>=K) THEN
!
!   Triangularity satisfied; compute the 3j coefficient
!
!   Begin with the logarithm of the square of the leading term
!
         EXPTRM = -LOG(DBLE(KA*KB))
!
!   Compute the logarithm of the square root of the leading term
!   and the factorial part that doesn't depend on the parity of
!   KA+KB+K (the delta factor)
!
         KAPKB = KA + KB
         KABKP = KAPKB + K
         KAMKB = KA - KB
         KBMKA = KB - KA
         EXPTRM = 0.5D00*(EXPTRM + GAM(KAPKB-K)+GAM(KAMKB+K+1)+GAM(KBMKA+K+1)-&
            GAM(KABKP+1))
!
!   The remainder depends on the parity of KA+KB+K
!
         IF (MOD(KABKP,2) == 0) THEN
!
!   Computation for even parity case
!
!   Include the phase factor: a minus sign if necessary
!
            IF (MOD(3*KABKP/2,2) == 0) THEN
               CLRX = 1.0D00
            ELSE
               CLRX = -1.0D00
            ENDIF
!
!   Include the contribution from the factorials
!
            EXPTRM = EXPTRM + GAM((KABKP+2)/2) - GAM((KAPKB-K)/2) - GAM((KAMKB+&
               K+2)/2) - GAM((KBMKA+K+2)/2)
!
         ELSE
!
!   Computation for odd parity case
!
!   Include the phase factor: a minus sign if necessary
!
            IF (MOD((3*KABKP - 1)/2,2) == 0) THEN
               CLRX = 1.0D00
            ELSE
               CLRX = -1.0D00
            ENDIF
!
!   Include the contribution from the factorials
!
            EXPTRM = EXPTRM + GAM((KABKP+1)/2) - GAM((KAPKB-K+1)/2) - GAM((&
               KAMKB+K+1)/2) - GAM((KBMKA+K+1)/2)
!
         ENDIF
!
!   Final assembly
!
         CLRX = CLRX*EXP(EXPTRM)
!
      ELSE
!
!   Triangularity violated; set the coefficient to zero
!
         CLRX = 0.0D00
!
      ENDIF
!
      RETURN
!
      END FUNCTION CLRX

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: sylas/grasp-continuum

max_forks_repo_name: sylas/grasp-continuum