!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!> This module allows one to convert a sparse matrix to a dense matrix. It also
!! supports dense the dense versions of core matrix routines. This module is
!! used in situations where matrices become too dense for good sparse matrix
!! performance.
MODULE DMatrixModule
  USE DataTypesModule, ONLY : NTREAL, NTCOMPLEX
  USE SMatrixModule, ONLY : Matrix_lsr, Matrix_lsc, &
       & ConstructMatrixFromTripletList
  USE TripletListModule, ONLY : TripletList_r, TripletList_c, &
       & AppendToTripletList, ConstructTripletList, DestructTripletList
  USE TripletModule, ONLY : Triplet_r, Triplet_c
  IMPLICIT NONE
  PRIVATE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> A datatype for storing a dense matrix.
  TYPE, PUBLIC :: Matrix_ldr
     REAL(NTREAL), DIMENSION(:,:), ALLOCATABLE :: DATA !< values of the matrix.
     INTEGER :: rows !< Matrix dimension: rows.
     INTEGER :: columns !< Matrix dimension: columns.
  END TYPE Matrix_ldr
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> A datatype for storing a dense matrix.
  TYPE, PUBLIC :: Matrix_ldc
     !> values of the matrix.
     COMPLEX(NTCOMPLEX), DIMENSION(:,:), ALLOCATABLE :: DATA
     INTEGER :: rows !< Matrix dimension: rows.
     INTEGER :: columns !< Matrix dimension: columns.
  END TYPE Matrix_ldc
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  PUBLIC :: ConstructEmptyMatrix
  PUBLIC :: ConstructMatrixDFromS
  PUBLIC :: ConstructMatrixSFromD
  PUBLIC :: CopyMatrix
  PUBLIC :: DestructMatrix
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  PUBLIC :: SplitMatrix
  PUBLIC :: ComposeMatrix
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  PUBLIC :: MatrixNorm
  PUBLIC :: IncrementMatrix
  PUBLIC :: MultiplyMatrix
  PUBLIC :: TransposeMatrix
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  INTERFACE Matrix_ldr
     MODULE PROCEDURE ConstructEmptyMatrix_ldr
  END INTERFACE Matrix_ldr
  INTERFACE Matrix_ldc
     MODULE PROCEDURE ConstructEmptyMatrix_ldc
  END INTERFACE Matrix_ldc
  INTERFACE ConstructEmptyMatrix
     MODULE PROCEDURE ConstructEmptyMatrixSup_ldr
     MODULE PROCEDURE ConstructEmptyMatrixSup_ldc
  END INTERFACE ConstructEmptyMatrix
  INTERFACE ConstructMatrixDFromS
     MODULE PROCEDURE ConstructMatrixDFromS_ldr
     MODULE PROCEDURE ConstructMatrixDFromS_ldc
  END INTERFACE ConstructMatrixDFromS
  INTERFACE ConstructMatrixSFromD
     MODULE PROCEDURE ConstructMatrixSFromD_ldr
     MODULE PROCEDURE ConstructMatrixSFromD_ldc
  END INTERFACE ConstructMatrixSFromD
  INTERFACE CopyMatrix
     MODULE PROCEDURE CopyMatrix_ldr
     MODULE PROCEDURE CopyMatrix_ldc
  END INTERFACE CopyMatrix
  INTERFACE DestructMatrix
     MODULE PROCEDURE DestructMatrix_ldr
     MODULE PROCEDURE DestructMatrix_ldc
  END INTERFACE DestructMatrix
  INTERFACE SplitMatrix
     MODULE PROCEDURE SplitMatrix_ldr
     MODULE PROCEDURE SplitMatrix_ldc
  END INTERFACE SplitMatrix
  INTERFACE ComposeMatrix
     MODULE PROCEDURE ComposeMatrix_ldr
     MODULE PROCEDURE ComposeMatrix_ldc
  END INTERFACE ComposeMatrix
  INTERFACE MatrixNorm
     MODULE PROCEDURE MatrixNorm_ldr
     MODULE PROCEDURE MatrixNorm_ldc
  END INTERFACE MatrixNorm
  INTERFACE IncrementMatrix
     MODULE PROCEDURE IncrementMatrix_ldr
     MODULE PROCEDURE IncrementMatrix_ldc
  END INTERFACE IncrementMatrix
  INTERFACE MultiplyMatrix
     MODULE PROCEDURE MultiplyMatrix_ldr
     MODULE PROCEDURE MultiplyMatrix_ldc
  END INTERFACE MultiplyMatrix
  INTERFACE TransposeMatrix
     MODULE PROCEDURE TransposeMatrix_ldr
     MODULE PROCEDURE TransposeMatrix_ldc
  END INTERFACE TransposeMatrix
CONTAINS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> A subroutine wrapper for the empty constructor.
  PURE SUBROUTINE ConstructEmptyMatrixSup_ldr(this, rows, columns)
    !> The matrix to construct
    TYPE(Matrix_ldr), INTENT(INOUT) :: this
    !> Rows of the matrix
    INTEGER, INTENT(IN) :: rows
    !> Columns of the matrix
    INTEGER, INTENT(IN) :: columns

    CALL DestructMatrix(this)
    this = ConstructEmptyMatrix_ldr(rows, columns)
  END SUBROUTINE ConstructEmptyMatrixSup_ldr
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> Construct an empty dense matrix with a set number of rows and columns
  PURE FUNCTION ConstructEmptyMatrix_ldr(rows, columns) RESULT(this)
    !> The matrix to construct.
    TYPE(Matrix_ldr) :: this
    !> Rows of the matrix
    INTEGER, INTENT(IN) :: rows
    !> Columns of the matrix.
    INTEGER, INTENT(IN) :: columns

  this%rows = rows
  this%columns = columns

  ALLOCATE(this%DATA(rows,columns))

  END FUNCTION ConstructEmptyMatrix_ldr
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> A function that converts a sparse matrix to a dense matrix.
  PURE SUBROUTINE ConstructMatrixDFromS_ldr(sparse_matrix, dense_matrix)
    !> The sparse matrix to convert.
    TYPE(Matrix_lsr), INTENT(IN) :: sparse_matrix
    !> Output. Must be preallocated.
    TYPE(Matrix_ldr), INTENT(INOUT) :: dense_matrix
    !! Helper Variables
    TYPE(Triplet_r) :: temporary

  !! Helper Variables
  INTEGER :: inner_counter, outer_counter
  INTEGER :: elements_per_inner
  INTEGER :: total_counter

  CALL ConstructEmptyMatrix(dense_matrix, sparse_matrix%rows, &
       & sparse_matrix%columns)

  !! Loop over elements.
  dense_matrix%DATA = 0
  total_counter = 1
  DO outer_counter = 1, sparse_matrix%columns
     elements_per_inner = sparse_matrix%outer_index(outer_counter+1) - &
          & sparse_matrix%outer_index(outer_counter)
     temporary%index_column = outer_counter
     DO inner_counter = 1, elements_per_inner
        temporary%index_row = sparse_matrix%inner_index(total_counter)
        temporary%point_value = sparse_matrix%values(total_counter)
        dense_matrix%DATA(temporary%index_row, temporary%index_column) = &
             & temporary%point_value
        total_counter = total_counter + 1
     END DO
  END DO

  END SUBROUTINE ConstructMatrixDFromS_ldr
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> A function that converts a dense matrix to a sparse matrix.
  PURE SUBROUTINE ConstructMatrixSFromD_ldr(dense_matrix, sparse_matrix, &
       & threshold_in)
    !> Matrix to convert.
    TYPE(Matrix_ldr), INTENT(IN) :: dense_matrix
    !> Output matrix.
    TYPE(Matrix_lsr), INTENT(INOUT) :: sparse_matrix
    !> Value for pruning values to zero.
    REAL(NTREAL), INTENT(IN), OPTIONAL :: threshold_in
    !! Local Variables
    TYPE(Triplet_r) :: temporary
    TYPE(TripletList_r) :: temporary_list

  !! Local Variables
  INTEGER :: inner_counter, outer_counter
  INTEGER :: columns, rows

  columns = dense_matrix%columns
  rows = dense_matrix%rows

  IF (PRESENT(threshold_in)) THEN
     CALL ConstructTripletList(temporary_list)
     DO outer_counter = 1, columns
        temporary%index_column = outer_counter
        DO inner_counter = 1, rows
           temporary%point_value = &
                & dense_matrix%DATA(inner_counter,outer_counter)
           IF (ABS(temporary%point_value) .GT. threshold_in) THEN
              temporary%index_row = inner_counter
              CALL AppendToTripletList(temporary_list,temporary)
           END IF
        END DO
     END DO
  ELSE
     CALL ConstructTripletList(temporary_list, rows*columns)
     DO outer_counter = 1, columns
        temporary%index_column = outer_counter
        DO inner_counter = 1, rows
           temporary%point_value = &
                & dense_matrix%DATA(inner_counter,outer_counter)
           temporary%index_row = inner_counter
           temporary_list%DATA(inner_counter+rows*(outer_counter-1)) = &
                & temporary
        END DO
     END DO
  END IF

  CALL ConstructMatrixFromTripletList(sparse_matrix, temporary_list, &
       & rows, columns)
  CALL DestructTripletList(temporary_list)

  END SUBROUTINE ConstructMatrixSFromD_ldr
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> Copy the matrix A into the B.
  PURE SUBROUTINE CopyMatrix_ldr(matA, matB)
    !> The matrix to copy.
    TYPE(Matrix_ldr), INTENT(IN) :: matA
    !> matB = matA
    TYPE(Matrix_ldr), INTENT(INOUT) :: matB

  CALL ConstructEmptyMatrix(matB, matA%rows, matA%columns)
  matB%DATA = matA%DATA

  END SUBROUTINE CopyMatrix_ldr
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> Deallocate the memory associated with this matrix.
  PURE SUBROUTINE DestructMatrix_ldr(this)
    !> The matrix to delete.
    TYPE(Matrix_ldr), INTENT(INOUT) :: this

  IF (ALLOCATED(this%DATA)) THEN
     DEALLOCATE(this%DATA)
  END IF

  END SUBROUTINE DestructMatrix_ldr
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> AXPY for dense matrices. B = B + alpha*A
  PURE SUBROUTINE IncrementMatrix_ldr(MatA,MatB,alpha_in)
    !> MatA is added
    TYPE(Matrix_ldr), INTENT(IN) :: MatA
    !> MatB is incremented.
    TYPE(Matrix_ldr), INTENT(INOUT) :: MatB
    !> A scaling parameter.
    REAL(NTREAL), OPTIONAL, INTENT(IN) :: alpha_in
    !! Temporary
    REAL(NTREAL) :: alpha

  !! Process Optional Parameters
  IF (.NOT. PRESENT(alpha_in)) THEN
     alpha = 1.0_NTREAL
  ELSE
     alpha = alpha_in
  END IF

  MatB%DATA = MatB%DATA + alpha*MatA%DATA

  END SUBROUTINE IncrementMatrix_ldr
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> Compute the norm of a dense matrix.
  !> Computes the Frobenius norm.
  FUNCTION MatrixNorm_ldr(this) RESULT(norm)
    !! Parameters
    !> The matrix to compute the norm of.
    TYPE(Matrix_ldr), INTENT(IN) :: this
    !> The norm of the matrix.
    REAL(NTREAL) :: norm
    INTEGER :: II, JJ

    norm = 0
    DO II =1, this%rows
       DO JJ = 1,  this%columns
          norm = norm + this%DATA(II,JJ)**2
       END DO
    END DO
  END FUNCTION MatrixNorm_ldr
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> Transpose a dense matrix.
  PURE SUBROUTINE TransposeMatrix_ldr(matA, matAT)
    !> matA the matrix to transpose.
    TYPE(Matrix_ldr), INTENT(IN) :: matA
    !> matAT = matA^T.
    TYPE(Matrix_ldr), INTENT(INOUT) :: matAT

  CALL ConstructEmptyMatrix(matAT, matA%columns, matA%rows)
  matAT%DATA = TRANSPOSE(matA%DATA)

  END SUBROUTINE TransposeMatrix_ldr
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> Create a big matrix from an array of matrices by putting them one next
  !> to another.
  PURE SUBROUTINE ComposeMatrix_ldr(mat_array, block_rows, block_columns, &
       & out_matrix)
    !> The number of rows of the array of blocks.
    INTEGER, INTENT(IN) :: block_rows
    !> The number of columns of the array of blocks.
    INTEGER, INTENT(IN) :: block_columns
    !> 2d array of matrices to compose.
    TYPE(Matrix_ldr), DIMENSION(block_rows, block_columns), INTENT(IN) :: &
         & mat_array
    !> The composed matrix.
    TYPE(Matrix_ldr), INTENT(INOUT) :: out_matrix

  !! Local Data
  INTEGER, DIMENSION(block_rows+1) :: row_offsets
  INTEGER, DIMENSION(block_columns+1) :: column_offsets
  INTEGER :: out_rows, out_columns
  INTEGER :: II, JJ

  !! Determine the size of the big matrix
  out_columns = 0
  column_offsets(1) = 1
  out_columns = 0
  DO JJ = 1, block_columns
     column_offsets(JJ+1) = column_offsets(JJ) + mat_array(1,JJ)%columns
     out_columns = out_columns + mat_array(1,JJ)%columns
  END DO
  row_offsets(1) = 1
  out_rows = 0
  DO II = 1, block_rows
     row_offsets(II+1) = row_offsets(II) + mat_array(II,1)%rows
     out_rows = out_rows + mat_array(II,1)%rows
  END DO

  !! Allocate Memory
  CALL ConstructEmptyMatrix(out_matrix, out_columns, out_rows)

  !! Copy
  DO JJ = 1, block_columns
     DO II = 1, block_rows
        out_matrix%DATA(row_offsets(II):row_offsets(II+1)-1, &
             & column_offsets(JJ):column_offsets(JJ+1)-1) = &
             & mat_array(II,JJ)%DATA
     END DO
  END DO

  END SUBROUTINE ComposeMatrix_ldr
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> Split a sparse matrix into an array of sparse matrices.
  PURE SUBROUTINE SplitMatrix_ldr(this, block_rows, block_columns, &
       & split_array, block_size_row_in, block_size_column_in)
    !> The matrix to split.
    TYPE(Matrix_ldr), INTENT(IN) :: this
    !> Number of rows to split the matrix into.
    INTEGER, INTENT(IN) :: block_rows
    !> Number of columns to split the matrix into.
    INTEGER, INTENT(IN) :: block_columns
    !> A block_columns x block_rows array for the output to go into.
    TYPE(Matrix_ldr), DIMENSION(:,:), INTENT(INOUT) :: split_array
    !> Specifies the size of the  rows.
    INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: block_size_row_in
    !> Specifies the size of the columns.
    INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: block_size_column_in

  !! Local Data
  INTEGER, DIMENSION(block_rows) :: block_size_row
  INTEGER, DIMENSION(block_columns) :: block_size_column
  INTEGER, DIMENSION(block_rows+1) :: row_offsets
  INTEGER, DIMENSION(block_columns+1) :: column_offsets
  !! Temporary Variables
  INTEGER :: divisor_row, divisor_column
  INTEGER :: II, JJ

  !! Calculate the split sizes
  IF (PRESENT(block_size_row_in)) THEN
     block_size_row = block_size_row_in
  ELSE
     divisor_row = this%rows/block_rows
     block_size_row = divisor_row
     block_size_row(block_rows) = this%rows - divisor_row*(block_rows-1)
  END IF
  IF (PRESENT(block_size_column_in)) THEN
     block_size_column = block_size_column_in
  ELSE
     divisor_column = this%columns/block_columns
     block_size_column = divisor_column
     block_size_column(block_columns) = this%columns - &
          & divisor_column*(block_columns-1)
  END IF

  !! Copy the block offsets
  row_offsets(1) = 1
  DO II = 1, block_rows
     row_offsets(II+1) = row_offsets(II) + block_size_row(II)
  END DO
  column_offsets(1) = 1
  DO JJ = 1, block_columns
     column_offsets(JJ+1) = column_offsets(JJ) + block_size_column(JJ)
  END DO

  !! Copy
  DO JJ = 1, block_columns
     DO II = 1, block_rows
        CALL ConstructEmptyMatrix(split_array(II,JJ), block_size_column(JJ), &
             & block_size_row(II))
        split_array(II,JJ)%DATA = &
             & this%DATA(row_offsets(II):row_offsets(II+1)-1, &
             & column_offsets(JJ):column_offsets(JJ+1)-1)
     END DO
  END DO

  END SUBROUTINE SplitMatrix_ldr
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> A wrapper for multiplying two dense matrices.
  SUBROUTINE MultiplyMatrix_ldr(MatA, MatB, MatC, IsATransposed_in, &
       & IsBTransposed_in)
    !> The first matrix.
    TYPE(Matrix_ldr), INTENT(IN) :: MatA
    !> The second matrix.
    TYPE(Matrix_ldr), INTENT(IN) :: MatB
    !> MatC = MatA*MatB.
    TYPE(Matrix_ldr), INTENT(INOUT) :: MatC
    !> True if A is already transposed.
    LOGICAL, OPTIONAL, INTENT(IN) :: IsATransposed_in
    !> True if B is already transposed.
    LOGICAL, OPTIONAL, INTENT(IN) :: IsBTransposed_in
    !! Local variables
    CHARACTER :: TRANSA
    CHARACTER :: TRANSB
    INTEGER :: M
    INTEGER :: N
    INTEGER :: K
    DOUBLE PRECISION, PARAMETER :: ALPHA = 1.0
    INTEGER :: LDA
    INTEGER :: LDB
    DOUBLE PRECISION, PARAMETER :: BETA = 0.0
    INTEGER :: LDC

    !! Optional Parameters
    TRANSA = 'N'
    IF (PRESENT(IsATransposed_in)) THEN
       IF (IsATransposed_in) THEN
          TRANSA = 'T'
       END IF
    END IF
    TRANSB = 'N'
    IF (PRESENT(IsBTransposed_in)) THEN
       IF (IsBTransposed_in) THEN
          TRANSB = 'T'
       END IF
    END IF

    !! Setup Lapack
    IF (TRANSA .EQ. 'T') THEN
       M = MatA%columns
    ELSE
       M = MatA%rows
    END IF

    IF (TRANSB .EQ. 'T') THEN
       N = MatB%rows
    ELSE
       N = MatB%columns
    END IF

    IF (TRANSA .EQ. 'T') THEN
       K = MatA%rows
    ELSE
       K = MatA%columns
    END IF

    IF (TRANSA .EQ. 'T') THEN
       LDA = K
    ELSE
       LDA = M
    END IF

    IF (TRANSB .EQ. 'T') THEN
       LDB = N
    ELSE
       LDB = K
    END IF

    LDC = M

    !! Multiply
    CALL ConstructEmptyMatrix(MatC, M, N)
    CALL DGEMM(TRANSA, TRANSB, M, N, K, ALPHA, MatA%DATA, LDA, MatB%DATA, &
         & LDB, BETA, MatC%DATA, LDC)

  END SUBROUTINE MultiplyMatrix_ldr
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> A subroutine style wrapper for the constructor.
  PURE SUBROUTINE ConstructEmptyMatrixSup_ldc(this, rows, columns)
    !> The matrix to construct.
    TYPE(Matrix_ldc), INTENT(INOUT) :: this
    !> The number of rows of the matrix.
    INTEGER, INTENT(IN) :: rows
    !> The number of columns o the matrix.
    INTEGER, INTENT(IN) :: columns

    CALL DestructMatrix(this)
    this = ConstructEmptyMatrix_ldc(rows, columns)
  END SUBROUTINE ConstructEmptyMatrixSup_ldc
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> Construct an empty dense matrix with a set number of rows and columns
  PURE FUNCTION ConstructEmptyMatrix_ldc(rows, columns) RESULT(this)
    !> The matrix to construct.
    TYPE(Matrix_ldc) :: this
    !> Rows of the matrix
    INTEGER, INTENT(IN) :: rows
    !> Columns of the matrix.
    INTEGER, INTENT(IN) :: columns

  this%rows = rows
  this%columns = columns

  ALLOCATE(this%DATA(rows,columns))

  END FUNCTION ConstructEmptyMatrix_ldc
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> A function that converts a sparse matrix to a dense matrix.
  PURE SUBROUTINE ConstructMatrixDFromS_ldc(sparse_matrix, dense_matrix)
    !> The sparse matrix to convert.
    TYPE(Matrix_lsc), INTENT(IN) :: sparse_matrix
    !> Dense matrix output. Must be preallocated.
    TYPE(Matrix_ldc), INTENT(INOUT) :: dense_matrix
    !! Helper Variables
    TYPE(Triplet_c) :: temporary

  !! Helper Variables
  INTEGER :: inner_counter, outer_counter
  INTEGER :: elements_per_inner
  INTEGER :: total_counter

  CALL ConstructEmptyMatrix(dense_matrix, sparse_matrix%rows, &
       & sparse_matrix%columns)

  !! Loop over elements.
  dense_matrix%DATA = 0
  total_counter = 1
  DO outer_counter = 1, sparse_matrix%columns
     elements_per_inner = sparse_matrix%outer_index(outer_counter+1) - &
          & sparse_matrix%outer_index(outer_counter)
     temporary%index_column = outer_counter
     DO inner_counter = 1, elements_per_inner
        temporary%index_row = sparse_matrix%inner_index(total_counter)
        temporary%point_value = sparse_matrix%values(total_counter)
        dense_matrix%DATA(temporary%index_row, temporary%index_column) = &
             & temporary%point_value
        total_counter = total_counter + 1
     END DO
  END DO

  END SUBROUTINE ConstructMatrixDFromS_ldc
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> A function that converts a dense matrix to a sparse matrix.
  PURE SUBROUTINE ConstructMatrixSFromD_ldc(dense_matrix, sparse_matrix, &
       & threshold_in)
    !> The matrix to convert.
    TYPE(Matrix_ldc), INTENT(IN) :: dense_matrix
    !> The sparse output matrix.
    TYPE(Matrix_lsc), INTENT(INOUT) :: sparse_matrix
    !> Value for pruning values to zero.
    REAL(NTREAL), INTENT(IN), OPTIONAL :: threshold_in
    !! Local Variables
    TYPE(Triplet_c) :: temporary
    TYPE(TripletList_c) :: temporary_list

  !! Local Variables
  INTEGER :: inner_counter, outer_counter
  INTEGER :: columns, rows

  columns = dense_matrix%columns
  rows = dense_matrix%rows

  IF (PRESENT(threshold_in)) THEN
     CALL ConstructTripletList(temporary_list)
     DO outer_counter = 1, columns
        temporary%index_column = outer_counter
        DO inner_counter = 1, rows
           temporary%point_value = &
                & dense_matrix%DATA(inner_counter,outer_counter)
           IF (ABS(temporary%point_value) .GT. threshold_in) THEN
              temporary%index_row = inner_counter
              CALL AppendToTripletList(temporary_list,temporary)
           END IF
        END DO
     END DO
  ELSE
     CALL ConstructTripletList(temporary_list, rows*columns)
     DO outer_counter = 1, columns
        temporary%index_column = outer_counter
        DO inner_counter = 1, rows
           temporary%point_value = &
                & dense_matrix%DATA(inner_counter,outer_counter)
           temporary%index_row = inner_counter
           temporary_list%DATA(inner_counter+rows*(outer_counter-1)) = &
                & temporary
        END DO
     END DO
  END IF

  CALL ConstructMatrixFromTripletList(sparse_matrix, temporary_list, &
       & rows, columns)
  CALL DestructTripletList(temporary_list)

  END SUBROUTINE ConstructMatrixSFromD_ldc
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> Copy the matrix A into the B.
  PURE SUBROUTINE CopyMatrix_ldc(matA, matB)
    !> The matrix to copy.
    TYPE(Matrix_ldc), INTENT(IN) :: matA
    !> matB = matA
    TYPE(Matrix_ldc), INTENT(INOUT) :: matB

  CALL ConstructEmptyMatrix(matB, matA%rows, matA%columns)
  matB%DATA = matA%DATA

  END SUBROUTINE CopyMatrix_ldc
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> Deallocate the memory associated with this matrix.
  PURE SUBROUTINE DestructMatrix_ldc(this)
    !> This the matrix to delete.
    TYPE(Matrix_ldc), INTENT(INOUT) :: this

  IF (ALLOCATED(this%DATA)) THEN
     DEALLOCATE(this%DATA)
  END IF

  END SUBROUTINE DestructMatrix_ldc
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> AXPY for dense matrices. B = B + alpha*A
  PURE SUBROUTINE IncrementMatrix_ldc(MatA,MatB,alpha_in)
    !> MatA is added
    TYPE(Matrix_ldc), INTENT(IN) :: MatA
    !> MatB is incremented.
    TYPE(Matrix_ldc), INTENT(INOUT) :: MatB
    !> A scaling parameter.
    REAL(NTREAL), OPTIONAL, INTENT(IN) :: alpha_in
    !! Temporary
    REAL(NTREAL) :: alpha

  !! Process Optional Parameters
  IF (.NOT. PRESENT(alpha_in)) THEN
     alpha = 1.0_NTREAL
  ELSE
     alpha = alpha_in
  END IF

  MatB%DATA = MatB%DATA + alpha*MatA%DATA

  END SUBROUTINE IncrementMatrix_ldc
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> Compute the norm of a dense matrix.
  !> Computes the Frobenius norm.
  FUNCTION MatrixNorm_ldc(this) RESULT(norm)
    !> The matrix to compute the norm of.
    TYPE(Matrix_ldc), INTENT(IN) :: this
    !> The norm of the matrix.
    REAL(NTREAL) :: norm
    !! Local Variables
    INTEGER :: II, JJ
    COMPLEX(NTCOMPLEX) :: val, conjval

    norm = 0
    DO II =1, this%rows
       DO JJ = 1,  this%columns
          val = this%DATA(II,JJ)
          conjval = CONJG(val)
          norm = norm + REAL(val*conjval,KIND=NTREAL)
       END DO
    END DO
  END FUNCTION MatrixNorm_ldc
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> Transpose a dense matrix.
  PURE SUBROUTINE TransposeMatrix_ldc(matA, matAT)
    !> The matrix to transpose.
    TYPE(Matrix_ldc), INTENT(IN) :: matA
    !> matAT = matA^T.
    TYPE(Matrix_ldc), INTENT(INOUT) :: matAT

  CALL ConstructEmptyMatrix(matAT, matA%columns, matA%rows)
  matAT%DATA = TRANSPOSE(matA%DATA)

  END SUBROUTINE TransposeMatrix_ldc
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> Create a big matrix from an array of matrices by putting them one next
  !> to another.
  PURE SUBROUTINE ComposeMatrix_ldc(mat_array, block_rows, block_columns, &
       & out_matrix)
    !> The number of rows of the array of blocks.
    INTEGER, INTENT(IN) :: block_rows
    !> The number of columns of the array of blocks.
    INTEGER, INTENT(IN) :: block_columns
    !> 2d array of matrices to compose.
    TYPE(Matrix_ldc), DIMENSION(block_rows, block_columns), INTENT(IN) :: &
         & mat_array
    !> The composed matrix.
    TYPE(Matrix_ldc), INTENT(INOUT) :: out_matrix

  !! Local Data
  INTEGER, DIMENSION(block_rows+1) :: row_offsets
  INTEGER, DIMENSION(block_columns+1) :: column_offsets
  INTEGER :: out_rows, out_columns
  INTEGER :: II, JJ

  !! Determine the size of the big matrix
  out_columns = 0
  column_offsets(1) = 1
  out_columns = 0
  DO JJ = 1, block_columns
     column_offsets(JJ+1) = column_offsets(JJ) + mat_array(1,JJ)%columns
     out_columns = out_columns + mat_array(1,JJ)%columns
  END DO
  row_offsets(1) = 1
  out_rows = 0
  DO II = 1, block_rows
     row_offsets(II+1) = row_offsets(II) + mat_array(II,1)%rows
     out_rows = out_rows + mat_array(II,1)%rows
  END DO

  !! Allocate Memory
  CALL ConstructEmptyMatrix(out_matrix, out_columns, out_rows)

  !! Copy
  DO JJ = 1, block_columns
     DO II = 1, block_rows
        out_matrix%DATA(row_offsets(II):row_offsets(II+1)-1, &
             & column_offsets(JJ):column_offsets(JJ+1)-1) = &
             & mat_array(II,JJ)%DATA
     END DO
  END DO

  END SUBROUTINE ComposeMatrix_ldc
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> Split a sparse matrix into an array of sparse matrices.
  PURE SUBROUTINE SplitMatrix_ldc(this, block_rows, block_columns, &
       & split_array, block_size_row_in, block_size_column_in)
    !> The matrix to split.
    TYPE(Matrix_ldc), INTENT(IN) :: this
    !> Number of rows to split the matrix into.
    INTEGER, INTENT(IN) :: block_rows
    !> Number of columns to split the matrix into.
    INTEGER, INTENT(IN) :: block_columns
    !> A COLUMNxROW array for the output to go into.
    TYPE(Matrix_ldc), DIMENSION(:,:), INTENT(INOUT) :: split_array
    !> Specifies the size of the  rows.
    INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: block_size_row_in
    !> Specifies the size of the columns.
    INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: block_size_column_in

  !! Local Data
  INTEGER, DIMENSION(block_rows) :: block_size_row
  INTEGER, DIMENSION(block_columns) :: block_size_column
  INTEGER, DIMENSION(block_rows+1) :: row_offsets
  INTEGER, DIMENSION(block_columns+1) :: column_offsets
  !! Temporary Variables
  INTEGER :: divisor_row, divisor_column
  INTEGER :: II, JJ

  !! Calculate the split sizes
  IF (PRESENT(block_size_row_in)) THEN
     block_size_row = block_size_row_in
  ELSE
     divisor_row = this%rows/block_rows
     block_size_row = divisor_row
     block_size_row(block_rows) = this%rows - divisor_row*(block_rows-1)
  END IF
  IF (PRESENT(block_size_column_in)) THEN
     block_size_column = block_size_column_in
  ELSE
     divisor_column = this%columns/block_columns
     block_size_column = divisor_column
     block_size_column(block_columns) = this%columns - &
          & divisor_column*(block_columns-1)
  END IF

  !! Copy the block offsets
  row_offsets(1) = 1
  DO II = 1, block_rows
     row_offsets(II+1) = row_offsets(II) + block_size_row(II)
  END DO
  column_offsets(1) = 1
  DO JJ = 1, block_columns
     column_offsets(JJ+1) = column_offsets(JJ) + block_size_column(JJ)
  END DO

  !! Copy
  DO JJ = 1, block_columns
     DO II = 1, block_rows
        CALL ConstructEmptyMatrix(split_array(II,JJ), block_size_column(JJ), &
             & block_size_row(II))
        split_array(II,JJ)%DATA = &
             & this%DATA(row_offsets(II):row_offsets(II+1)-1, &
             & column_offsets(JJ):column_offsets(JJ+1)-1)
     END DO
  END DO

  END SUBROUTINE SplitMatrix_ldc
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !> A wrapper for multiplying two dense matrices.
  SUBROUTINE MultiplyMatrix_ldc(MatA, MatB, MatC, IsATransposed_in, &
       & IsBTransposed_in)
    !> The first matrix.
    TYPE(Matrix_ldc), INTENT(IN) :: MatA
    !> The second matrix.
    TYPE(Matrix_ldc), INTENT(IN) :: MatB
    !> MatC = MatA*MatB.
    TYPE(Matrix_ldc), INTENT(INOUT) :: MatC
    !> True if A is already transposed.
    LOGICAL, OPTIONAL, INTENT(IN) :: IsATransposed_in
    !> True if B is already transposed.
    LOGICAL, OPTIONAL, INTENT(IN) :: IsBTransposed_in
    !! Local variables
    CHARACTER :: TRANSA
    CHARACTER :: TRANSB
    INTEGER :: M
    INTEGER :: N
    INTEGER :: K
    COMPLEX*16, PARAMETER :: ALPHA = 1.0
    INTEGER :: LDA
    INTEGER :: LDB
    COMPLEX*16, PARAMETER :: BETA = 0.0
    INTEGER :: LDC

    !! Optional Parameters
    TRANSA = 'N'
    IF (PRESENT(IsATransposed_in)) THEN
       IF (IsATransposed_in) THEN
          TRANSA = 'T'
       END IF
    END IF
    TRANSB = 'N'
    IF (PRESENT(IsBTransposed_in)) THEN
       IF (IsBTransposed_in) THEN
          TRANSB = 'T'
       END IF
    END IF

    !! Setup Lapack
    IF (TRANSA .EQ. 'T') THEN
       M = MatA%columns
    ELSE
       M = MatA%rows
    END IF

    IF (TRANSB .EQ. 'T') THEN
       N = MatB%rows
    ELSE
       N = MatB%columns
    END IF

    IF (TRANSA .EQ. 'T') THEN
       K = MatA%rows
    ELSE
       K = MatA%columns
    END IF

    IF (TRANSA .EQ. 'T') THEN
       LDA = K
    ELSE
       LDA = M
    END IF

    IF (TRANSB .EQ. 'T') THEN
       LDB = N
    ELSE
       LDB = K
    END IF

    LDC = M

    !! Multiply
    CALL ConstructEmptyMatrix(MatC, M, N)
    CALL ZGEMM(TRANSA, TRANSB, M, N, K, ALPHA, MatA%DATA, LDA, MatB%DATA, &
         & LDB, BETA, MatC%DATA, LDC)

  END SUBROUTINE MultiplyMatrix_ldc
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
END MODULE DMatrixModule


text:
! Copyright (C) 2015-2017, Simone Marsili
! All rights reserved.
! License: BSD 3 clause

module scrs
  use kinds
  implicit none
  private
  public :: compute_scores
  public :: print_mat

contains

  pure real(kflt) function frobenius(a)
    real(kflt), intent(in) :: a(:,:)

    frobenius = sqrt(sum(a**2))
  end function frobenius

  pure subroutine apc_correction(a)
    real(kflt), intent(inout) :: a(:,:)
    real(kflt), allocatable :: sums(:)
    real(kflt) :: totsum
    integer :: err
    integer :: i,j
    integer :: n

    n = size(a,1)
    allocate(sums(n),stat=err)
    do i = 1,n
       sums(i) = sum(a(:,i))
    end do
    totsum = sum(sums)
    do j = 1,n
       do i = j,n
          a(i,j) = a(i,j) - sums(i)*sums(j)/totsum
          a(j,i) = a(i,j)
       end do
    end do
    deallocate(sums)
  end subroutine apc_correction

  subroutine print_mat(mat,uscrs,frmt,nerr)
    ! https://people.sc.fsu.edu/~jburkardt/data/mm/mm.html
    real(kflt), intent(in) :: mat(:,:)
    integer, intent(in) :: uscrs
    character(1), intent(in) :: frmt
    integer, intent(out) :: nerr
    logical :: inds = .true.
    integer :: i,j,n1,n2

    nerr = 0
    ! check format
    if (.not. any(["r","c","a"] == frmt)) then
       write(0,*) "error ! unknown format for scores matrix file"
       nerr = nerr + 1
       return
    end if
    n1 = size(mat,1)
    n2 = size(mat,2)
    ! write header line
    select case(frmt)
    case("r")
       do j = 1,n1
          do i = j,n2
             write(uscrs,'(i5,1x,i5,1x,f8.4)') i, j, mat(i,j)
          end do
       end do
    case("a")
       write(uscrs,'(a)') "%%MatrixMarket matrix array real symmetric"
       write(uscrs,'(a)') "%Lower triangular, column-major order"
       write(uscrs,'(i5,1x,i5)') n1,n2
       do j = 1,n1
          do i = j,n2
             write(uscrs,'(f8.4)') mat(i,j)
          end do
       end do
    case("c")
       write(uscrs,'(a)') "%%MatrixMarket matrix coordinate real symmetric"
       write(uscrs,'(a)') "%Lower triangular, column-major order"
       write(uscrs,'(i5,1x,i5,1x,i5)') n1,n2,n1*n2
       do j = 1,n1
          do i = j,n2
             write(uscrs,'(i5,1x,i5,1x,f8.4)') i, j, mat(i,j)
          end do
       end do
    end select

  end subroutine print_mat

  subroutine compute_scores(nv,ns,couplings,ignore_pivot,scores)
    implicit none
    integer(kint), intent(in) :: nv,ns
    real(kflt), intent(in) :: couplings(ns,ns,nv,nv)
    integer, intent(in) :: ignore_pivot
    real(kflt), intent(out) :: scores(nv,nv)
    integer :: iv,jv,err,is,js
    real(kflt) :: sums(nv),rsum(ns),csum(ns),matij(ns,ns),matji(ns,ns)
    real(kflt) :: totsum
    integer :: inds(ns-1)

    scores = 0.0_kflt
    do iv = 1,nv-1
       do jv = iv+1,nv
          if (ignore_pivot >= 1 .and. ignore_pivot <= ns) then
             is = 0
             do js = 1,ns
                if (js /= ignore_pivot) then
                   is = is + 1
                   inds(is) = js
                end if
             end do
             matij = couplings(inds,inds,iv,jv)
             matji = couplings(inds,inds,jv,iv)
          else
             matij = couplings(:,:,iv,jv)
             matji = couplings(:,:,jv,iv)
          end if
          call double_center(matij)
          call double_center(matji)
          scores(iv,jv) = sum((matij + transpose(matji))**2)
          scores(jv,iv) = scores(iv,jv)
       end do
    end do
    scores = 0.5_kflt * sqrt(scores)
    do iv = 1,nv
       sums(iv) = sum(scores(:,iv))
    end do
    totsum = sum(sums)
    do jv = 1,nv
       do iv = jv,nv
          scores(iv,jv) = scores(iv,jv) - sums(iv)*sums(jv)/totsum
          scores(jv,iv) = scores(iv,jv)
       end do
    end do
  end subroutine compute_scores

  pure subroutine double_center(a)
    real(kflt), intent(inout) :: a(:,:)
    integer :: i,j,n,m
    real(kflt) :: totsum
    real(kflt), allocatable :: rsum(:),csum(:)
    integer :: err
    n = size(a,dim=1)
    m = size(a,dim=2)
    allocate(rsum(n),csum(m),stat=err)
    totsum = sum(a) / real(n*m)
    csum = sum(a,dim=2) / real(n)
    rsum = sum(a,dim=1) / real(m)
    do j = 1,m
       do i = 1,n
          a(i,j) =  a(i,j) - csum(i) - rsum(j) + totsum
       end do
    end do
    deallocate(rsum,csum)
  end subroutine double_center

end module scrs

max_stars_repo_name: simomarsili/fmpl

max_issues_repo_name: simomarsili/fmpl

max_forks_repo_name: simomarsili/fmpl


        !COMPILER-GENERATED INTERFACE MODULE: Wed Jul 21 10:02:25 2021
        ! This source file is for reference only and may not completely
        ! represent the generated interface used by the compiler.
        MODULE ERRORS__genmod
          INTERFACE
            SUBROUTINE ERRORS
            END SUBROUTINE ERRORS
          END INTERFACE
        END MODULE ERRORS__genmod


text:
subroutine DIFCO(nq, nz, mubar_z, T, den, edd, &
                 hscale, tauedd, DK, H_atm, bx1x2, scale_H)
  use photochem_data, only: grav_z, mass
  implicit none

  integer, intent(in) :: nq, nz
  real(8), intent(in) :: mubar_z(nz), T(nz), den(nz), edd(nz)

  real(8), intent(out) :: hscale(nz), tauedd(nz), DK(nz)
  real(8), intent(out) :: H_atm(nz)
  real(8), intent(out) :: bx1x2(nq,nz)
  real(8), intent(out) :: scale_H(nq,nz)

  ! local variables
  real*8 bkmg, eddav, denav
  real*8 h, tav
  integer i,j

  ! ***** DK(I) = K*N AT GRID STEP I+1/2 *****
  DO i = 1 , nz - 1
    eddav = SQRT(EDD(i)*EDD(i+1))
                                 !average eddy diffusion in grid center
    denav = SQRT(DEN(i)*DEN(i+1))
                               !average density at grid center
    DK(i) = eddav*denav
  ENDDO
  !
  !   COMPUTE DIFFUSION LIFETIME AT EACH HEIGHT (H*H/K)
  DO i = 1 , nz
    bkmg = 1.38E-16/(1.67E-24*mubar_z(i)*grav_z(i))
    h = bkmg*T(i)
    HSCALE(i) = h
    TAUEDD(i) = h*h/EDD(i)
  ENDDO

  DO i = 1 , nz - 1
    bkmg = 1.38E-16/(1.67E-24*mubar_z(i)*grav_z(i))
    tav = SQRT(T(i)*T(i+1)) !average temperature at grid center
    H_ATM(i) = bkmg*tav
    !  compute scale heights of all the species in the atmosphere at all
    !   heights
    DO j = 1 , nq
      SCALE_H(j,i) = bkmg*tav*mubar_z(i)/MASS(j)
    ENDDO
    do j = 1,nq
      bx1x2(j,i) = 1.52d18*(1.d0/mass(j)+1.d0/mubar_z(i))**0.5d0*(tav**0.5)
    enddo

  ENDDO

  ! top layer of atmosphere
  bkmg = 1.38E-16/(1.67E-24*mubar_z(nz)*grav_z(nz))
  H_ATM(nz) = bkmg*tav
  DO j = 1 , nq
     SCALE_H(j,nz) = bkmg*T(nz)*mubar_z(nz)/MASS(j)
  ENDDO
  do j = 1,nq
    bx1x2(j,NZ) = 1.52d18*(1.d0/mass(j)+1.d0/mubar_z(nz))**0.5d0*(T(nz)**0.5d0)
  enddo

end subroutine

max_stars_repo_name: Nicholaswogan/PhotochemPy

max_issues_repo_name: StAtmos/PhotochemPy

max_forks_repo_name: StAtmos/PhotochemPy


text:
program tsunami

  ! Tsunami simulator.
  !
  ! Solves the non-linear 2-d shallow water equation system:
  !
  !     du/dt + u du/dx + v du/dy + g dh/dx = 0
  !     dv/dt + u dv/dx + v dv/dy + g dh/dy = 0
  !     dh/dt + d(hu)/dx + d(hv)/dy = 0
  !
  ! This version is parallelized and uses derived types.

  use iso_fortran_env, only: int32, real32
  use mod_field, only: Field, diffx, diffy

  implicit none

  integer(int32) :: n

  integer(int32), parameter :: im = 201, jm = 201
  integer(int32), parameter :: num_time_steps = 1000

  real(real32), parameter :: dt = 0.02
  real(real32), parameter :: dx = 1, dy = 1
  real(real32), parameter :: g = 9.8

  integer(int32), parameter :: ic = im / 2 + 1
  integer(int32), parameter :: jc = jm / 2 + 1
  real(real32), parameter :: decay = 0.02

  type(Field) :: h, hm, u, v

  real(real32) :: hmin, hmax, hmean

  u = Field('u', [im, jm])
  v = Field('v', [im, jm])
  h = Field('h', [im, jm])
  hm = Field('h_mean', [im, jm])

  ! initialize a gaussian blob in the center
  call h % init_gaussian(decay, ic, jc)

  hm = 10.

  call h % write(0)

  time_loop: do n = 1, num_time_steps

    ! compute u at next time step
    u = u - (u * diffx(u) / dx + v * diffy(u) / dy &
      + g * diffx(h) / dx) * dt

    ! compute v at next time step
    v = v - (u * diffx(v) / dx + v * diffy(v) / dy &
      + g * diffy(h) / dy) * dt

    ! compute h at next time step
    h = h - (diffx(u * (hm + h)) / dx &
           + diffy(v * (hm + h)) / dy) * dt

    hmin = minval(h % data)
    call co_min(hmin, 1)

    hmax = maxval(h % data)
    call co_max(hmax, 1)

    hmean = sum(h % data(h % lb(1):h % ub(1),h % lb(2):h % ub(2))) &
          / size(h % data(h % lb(1):h % ub(1),h % lb(2):h % ub(2)))
    call co_sum(hmean, 1)
    hmean = hmean / num_images()

    if (this_image() == 1) &
      print '(a, i5, 3(f10.6))', &
        'step, min(h), max(h), mean(h):', &
         n, hmin, hmax, hmean

    call h % write(n)

  end do time_loop

end program tsunami

max_stars_repo_name: AAGAN/tsunami

max_issues_repo_name: AAGAN/tsunami

max_forks_repo_name: AAGAN/tsunami


text:
Campus Counsel provides inhouse legal advice for UC Davis on matters large and small.  Matters they have advised on include negotiations with Sodexho, and neverending Howard Zochlinski case.  Steve Drown is the head lawyer in the Campus Counsel office.

If wikizens would like to request information from the university that some lowerlevel official denied you, contact Lynnette Temple.  She handles information requests under the wiki:wikipedia:Freedom of Information Act (United States) Freedom of Information Act and the http://www.google.com/search?qCalifornia+Public+Records+Act California Public Records Act.

They are located on the second floor of Mrak Hall.  The office entrance is fairly unwelcoming  there are chairs in a tiny waiting room, and you cant even talk to a person inside unless they come out to speak to you.  Charming.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
	SUBROUTINE GDTXTM ( igdfln, gdatim, npts, timfnd, trange,
     +			    iret )
C************************************************************************
C* GDTXTM								*
C*									*
C* This subroutine sets time information for a contour series.		*
C*									*
C* GDTXTM ( IGDFLN, GDATIM, NPTS, TIMFND, TRANGE, IRET )		*
C*									*
C* Input parameters:							*
C*	IGDFLN		INTEGER		Grid file number		*
C*	GDATIM		CHAR*		User input date/time		*
C*									*
C* Output parameters:							*
C*	NPTS		INTEGER		Number of points		*
C*	TIMFND (NPTS)	CHAR*		Time values			*
C*	TRANGE		CHAR*		Actual range of grid times	*
C*	IRET		INTEGER		Return code			*
C*					  1 = user typed RETURN		*
C*					  0 = normal return		*
C*					 -4 = invalid grid point	*
C*					-13 = no grids found in GDATTIM	*
C**									*
C* Log:									*
C* G. Huffman/USRA	 5/89	Consolidated from GDTDTA, GDTTTL	*
C* M. desJardins/GSFC	12/89	Fixed // for Apollo			*
C* K. Brill/GSC          5/90   Made local character strings larger	*
C* S. Schotz/GSC	 6/90	Get respnd locally from IP_RESP		*
C* T.W.Barker/SS/SSD	 8/91	Created from GDTTIM (GDTSER)		*
C* K. Brill/NMC		11/92	Fix ER_WMSG call			*
C* G. Krueger/EAI	 5/93	Modified for use in GDPLOT		*
C* G. Krueger/EAI	 9/93	Modified for use in GDTHGT		*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	CHARACTER*(*)	gdatim, timfnd (*), trange
	LOGICAL		respnd
C*
	PARAMETER	( NTIME = 200 )
	CHARACTER*36	time, timlst (NTIME), tstrt, tstop, ttt, timout
C------------------------------------------------------------------------

	iret = 0
	npts = 0
C
C*	Get all times which might have data.
C
	CALL GD_GTIM ( igdfln, NTIME, timlst, ntimin, ier )
	IF ( ntimin .eq. 0 ) THEN
	    iret = -13
	    CALL ER_WMSG ( 'GDTHGT', iret, ' ', ier )
	    RETURN
	END IF
C
C*	If GDATIM is 'LIST' and the session is interactive, list the
C*	times and get a range.
C
	CALL ST_LCUC ( gdatim, time, ier )
	IF ( time .eq. 'LIST' )  THEN
C
	    CALL IP_RESP ( respnd, ier )
	    IF ( respnd )  THEN
		CALL TI_DSPL ( ntimin, timlst, time, ier )
		IF ( ier .eq. 1 )  THEN
		    iret = 1
		    RETURN
		ELSE
		    CALL ST_LCUC ( time, time, ier )
		END IF
C
	    ELSE
		iret   = -4
		CALL ER_WMSG ( 'TI', iret, ' ', ier )
		RETURN
	    END IF
	END IF
C
C*	Narrow the list to that specified in TIME.
C
	CALL TG_FIND ( time, ntimin, timlst, timout, npts, timfnd, ier )
	IF ( ier .ne. 0 )  THEN
	    CALL ER_WMSG ( 'TG', ier, ' ', ier1 )
	    iret = -13
	    RETURN
	END IF
C
C*	Encode the times into the "actual" grid time range (used in
C*	building the default title) . . . reduce the start/stop times
C*	to YYMMDD/HHthh if the other stuff is zero.  Drop A00 if both
C*	have it, and drop however much of the stop YYMMDD/HH matches
C*	the start values (starting at the left).
C
	tstrt = timfnd (1)
	tstop = timfnd (npts)
C
C*	The test on character 15 (last digit of 'hhh') ensures that
C*	the work hasn't already been done on 13:17 ('hhhmm').
C
	IF ( ( tstrt ( 15:15 ) .ne. ' ' ) .and.
     +       ( tstop ( 15:15 ) .ne. ' ' ) )  THEN
	    IF ( ( tstrt ( 16:17 ) .eq. '00' ) .and.
     +	         ( tstop ( 16:17 ) .eq. '00' ) )  THEN
	        tstrt ( 16:17 ) = '  '
	        tstop ( 16:17 ) = '  '
	    END IF
	    IF ( ( tstrt ( 13:13 )   .eq. '0' ) .and.
     +	         ( tstop ( 13:13 )   .eq. '0' ) )  THEN
		ttt   = tstrt
	        tstrt = ttt ( 1:12 ) // ttt ( 14: )
		ttt   = tstop
	        tstop = ttt ( 1:12 ) // ttt ( 14: )
	    END IF
	END IF
C
	IF ( ( tstrt ( 10:11 ) .eq. '00' ) .and.
     +	     ( tstop ( 10:11 ) .eq. '00' ) )  THEN
	    ttt   = tstrt
	    tstrt = ttt ( 1:9 ) // ttt ( 12: )
	    ttt   = tstop
	    tstop = ttt ( 1:9 ) // ttt ( 12: )
	END IF
	CALL ST_LSTR ( tstrt, lstrt, ier )
	IF ( ( tstrt ( lstrt-2:lstrt ) .eq. 'A00' ) .and.
     +	     ( tstop ( lstrt-2:lstrt ) .eq. 'A00' ) )  THEN
	    lstrt = lstrt - 3
	    tstop = tstop ( 1:lstrt )
	ELSE IF ( tstrt ( 1:9 ) .eq. tstop ( 1:9 ) )  THEN
	    tstop = tstop ( 10: )
	END IF
	IF ( tstrt ( 1:6 ) .eq. tstop ( 1:6 ) )  THEN
	    tstop = tstop ( 7: )
	ELSE IF ( tstrt ( 1:4 ) .eq. tstop ( 1:4 ) )  THEN
	    tstop = tstop ( 5: )
	ELSE IF ( tstrt ( 1:2 ) .eq. tstop ( 1:2 ) )  THEN
	    tstop = tstop ( 3: )
	END IF
	trange = tstrt ( 1:lstrt ) // '-' // tstop
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
      MODULE core_C
!...Created by Pacific-Sierra Research 77to90  4.3E  06:37:37  12/28/06
!...Modified by Charlotte Froese Fischer
!                     Gediminas Gaigalas  10/05/17
      INTEGER :: NCORE
      END MODULE core_C

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: sylas/grasp-continuum

max_forks_repo_name: sylas/grasp-continuum


	module reallocate
	  !     nur 1 (generischer) Name wird von aussen angesprochen
	  interface doreallocate
	    module procedure doreallocate_r8_d1
	    module procedure doreallocate_r8_d2
	    module procedure doreallocate_c16_d1
	    module procedure doreallocate_c16_d2
	    module procedure doreallocate_i4_d1
	    module procedure doreallocate_i4_d2
	    module procedure hugo     !   ;)
	  end interface
	contains

	  !     leider sind mehrere subroutines notwendig fuer verschiedene Typen
	  subroutine doreallocate_r8_d1(tf, newdimension)
	    real*8, pointer :: hilfsfeld(:), tf(:)
	    allocate(hilfsfeld(newdimension))
	    !     nur 1 mal kopieren reicht
	    !     auch fuer mehrdimensionale Felder schaut die Zuweisung gleich aus
            min1=min(newdimension,size(tf,1))
	    hilfsfeld(1:min1)=tf(1:min1)
	    deallocate(tf)
	    !     der Zeiger wird nur auf das neue Feld umgebogen, nicht neu alloziert
	    tf=>hilfsfeld
	  end subroutine

	  subroutine doreallocate_r8_d2(tf, newdimension1, newdimension2)
	    real*8, pointer :: hilfsfeld(:,:), tf(:,:)
	    allocate(hilfsfeld(newdimension1,newdimension2))
            min1=min(newdimension1,size(tf,1))
            min2=min(newdimension2,size(tf,2))
	    hilfsfeld(1:min1,1:min2)=tf(1:min1,1:min2)
	    deallocate(tf)
	    tf=>hilfsfeld
	  end subroutine

	  subroutine doreallocate_c16_d1(tf, newdimension)
	    complex*16, pointer :: hilfsfeld(:), tf(:)
	    allocate(hilfsfeld(newdimension))
            min1=min(newdimension,size(tf,1))
	    hilfsfeld(1:min1)=tf(1:min1)
	    deallocate(tf)
	    tf=>hilfsfeld
	  end subroutine

	  subroutine doreallocate_c16_d2(tf, newdimension1, newdimension2)
	    complex*16, pointer :: hilfsfeld(:,:), tf(:,:)
	    allocate(hilfsfeld(newdimension1,newdimension2))
            min1=min(newdimension1,size(tf,1))
            min2=min(newdimension2,size(tf,2))
	    hilfsfeld(1:min1,1:min2)=tf(1:min1,1:min2)
	    deallocate(tf)
	    tf=>hilfsfeld
	  end subroutine

	  subroutine doreallocate_i4_d1(tf, newdimension)
	    integer*4, pointer :: hilfsfeld(:), tf(:)
	    allocate(hilfsfeld(newdimension))
            min1=min(newdimension,size(tf,1))
	    hilfsfeld(1:min1)=tf(1:min1)
	    deallocate(tf)
	    tf=>hilfsfeld
	  end subroutine

	  subroutine doreallocate_i4_d2(tf, newdimension1, newdimension2)
	    integer*4, pointer :: hilfsfeld(:,:), tf(:,:)
	    allocate(hilfsfeld(newdimension1,newdimension2))
            min1=min(newdimension1,size(tf,1))
            min2=min(newdimension2,size(tf,2))
	    hilfsfeld(1:min1,1:min2)=tf(1:min1,1:min2)
	    deallocate(tf)
	    tf=>hilfsfeld
	  end subroutine

	!     Es gibt auch Methoden, um das Programm unleserlich zu machen :-)
	!     das sollten wir besser vermeiden ;-)
	  subroutine hugo
	    write(*,*) " Hier ist Hugo"
	  end subroutine
	end module


text:
program phaml_master
use phaml
implicit none
type(phaml_solution_type) :: soln
call phaml_create(soln,nproc=2)
call phaml_solve_pde(soln,max_vert=500,print_grid_when=FINAL, &
                     print_grid_who=MASTER)
call phaml_compress(soln)
call phaml_popen(soln,11,"savef.dat")
call phaml_store(soln,11)
call phaml_pclose(soln,11)
call phaml_destroy(soln,finalize_mpi=.false.)
call phaml_create(soln,nproc=2)
call phaml_popen(soln,12,"savef.dat")
call phaml_restore(soln,12)
call phaml_pclose(soln,12)
call phaml_solve_pde(soln,max_vert=500,print_grid_when=FINAL, &
                     print_grid_who=MASTER,inc_factor=1.1_my_real)
call phaml_compress(soln)
call phaml_popen(soln,13,"saveu.dat","UNFORMATTED")
call phaml_store(soln,13)
call phaml_pclose(soln,13)
call phaml_destroy(soln,finalize_mpi=.false.)
call phaml_create(soln,nproc=2)
call phaml_popen(soln,14,"saveu.dat","UNFORMATTED")
call phaml_restore(soln,14)
call phaml_pclose(soln,14)
call phaml_solve_pde(soln,max_vert=500,print_grid_when=FINAL, &
                     print_grid_who=MASTER,inc_factor=1.1_my_real)
call phaml_destroy(soln)

end program phaml_master

max_stars_repo_name: qsnake/phaml

max_issues_repo_name: qsnake/phaml

max_forks_repo_name: qsnake/phaml


text:
      program dp_example
      integer array_size
      parameter(array_size=10)
      dimension array(array_size)
      write(6, *) 'I can store up to', array_size, 'numbers'
      do i=1,array_size
        array(i) = i
      enddo
      end


! This module defines derived type ObsContainerType:
!
! ObsContainerType -- contains a pointer to an object of type
! ObserveType.  Its purpose is to allow ObserveType objects to be
! stored in an array.
!-----------------------------------------------------------------------
module ObsContainerModule

  use KindModule, only: DP, I4B
  use ObserveModule, only: ObserveType

  implicit none

  private
  public :: ObsContainerType

  type :: ObsContainerType
    ! -- Public members
    class(ObserveType), pointer, public :: obsrv => null()
  end type ObsContainerType

end module ObsContainerModule


text:
lock(a, stat = st)
sync all ()
syncimages (*)
syncmemory (errmsg = e)
exit
goto 123
goto (1,2,3) 4+5
if (a+c) 1,2,3
stop a+b
error stop
unlock(a)
endprogram

max_stars_repo_name: LChernigovskaya/fortran-plugin

max_issues_repo_name: LChernigovskaya/fortran-plugin

max_forks_repo_name: LChernigovskaya/fortran-plugin


text:
!
! Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
! See https://llvm.org/LICENSE.txt for license information.
! SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception

! 2d int array, minloc with dim, mask and back

program main
  integer :: res(4)
  integer, dimension(4,4) :: array
  integer :: expect(4)
  array = reshape((/4,2,9,-7,9,1,5,5,8,-1,-1,5,-7,5,9,-7/),shape(array))
  res = minloc(array, dim = 1, mask = array .lt. 7, back=.true.)
  expect = (/4,2,3,4/)
  call check(res, expect, 4)
  !print *, "base=true:", res
  res = minloc(array, dim = 1, mask = array .lt. 7, back=.false.)
  expect = (/4,2,2,1/)
  call check(res, expect, 4)
  !print *, "back=false:", res
end

max_stars_repo_name: abrahamtovarmob/flang

max_issues_repo_name: abrahamtovarmob/flang

max_forks_repo_name: abrahamtovarmob/flang


text:
!
!  Copyright 2017-2019 SALMON developers
!
!  Licensed under the Apache License, Version 2.0 (the "License");
!  you may not use this file except in compliance with the License.
!  You may obtain a copy of the License at
!
!      http://www.apache.org/licenses/LICENSE-2.0
!
!  Unless required by applicable law or agreed to in writing, software
!  distributed under the License is distributed on an "AS IS" BASIS,
!  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
!  See the License for the specific language governing permissions and
!  limitations under the License.
!
subroutine prep_ini
use salmon_parallel, only: nproc_id_global, nproc_group_global
use salmon_communication, only: comm_is_root, comm_bcast
use scf_data
implicit none
integer :: imol,jj
real(8),parameter :: rtmp=1.d-10

if(comm_is_root(nproc_id_global))then
  open(60,file=file_ini)
  read(60,*) file_OUT_ini
  read(60,*) H_ini
  read(60,*) (rLsize_ini(jj),jj=1,3)
  read(60,*) num_mol
  read(60,*) rlatcon
end if

call comm_bcast(file_OUT_ini, nproc_group_global)
call comm_bcast(H_ini,        nproc_group_global)
call comm_bcast(rLsize_ini,   nproc_group_global)
call comm_bcast(num_mol,      nproc_group_global)
call comm_bcast(rlatcon,      nproc_group_global)

allocate(coo_mol_ini(3,num_mol))

if(comm_is_root(nproc_id_global))then
  do imol=1,num_mol
    read(60,*) (coo_mol_ini(jj,imol),jj=1,3)
  end do
  close(60)
end if

call comm_bcast(coo_mol_ini,nproc_group_global)

H_ini=H_ini/a_B
rlatcon=rlatcon/a_B
rLsize_ini(:)=rLsize_ini(:)/a_B

coo_mol_ini(:,:)=rlatcon*coo_mol_ini(:,:)

lg_end_ini(:)=int((rLsize_ini(:)+rtmp)/2.d0/Hgs(:))
do jj=1,3
  select case(imesh_oddeven(jj))
    case(1)
      lg_sta_ini(jj)=-(int((rLsize_ini(jj)+rtmp)/2.d0/Hgs(jj)))
    case(2)
      lg_sta_ini(jj)=-(int((rLsize_ini(jj)+rtmp)/2.d0/Hgs(jj)))+1
  end select
end do
lg_num_ini(:)=lg_end_ini(:)-lg_sta_ini(:)+1

end subroutine prep_ini

max_stars_repo_name: shunsuke-sato/SALMON

max_issues_repo_name: shunsuke-sato/SALMON

max_forks_repo_name: shunsuke-sato/SALMON


text:
10 MODULE PSM
11 INTERFACE PSI
PROCEDURE :: a
12 PROCEDURE :: b, c
MODULE PROCEDURE :: d
13 MODULE PROCEDURE :: e, f
14 END INTERFACE

CONTAINS

SUBROUTINE a()
END SUBROUTINE

16 SUBROUTINE b(i)
END SUBROUTINE

SUBROUTINE c(i,j)
END SUBROUTINE

SUBROUTINE d(x)
END SUBROUTINE

SUBROUTINE e(x,y)
17 END SUBROUTINE

SUBROUTINE f(x,y,z)
END SUBROUTINE

18 END MODULE

max_stars_repo_name: ouankou/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
      SUBROUTINE ML5_0_MP_BORN_AMPS_AND_WFS(P)
C
C     Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
C     By the MadGraph5_aMC@NLO Development Team
C     Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
C
C     Computes all the AMP and WFS in quadruple precision for the
C     phase space point P(0:3,NEXTERNAL)
C
C     Process: g g > t t~ QCD<=2 QED=0 [ virt = QCD ]
C
      IMPLICIT NONE
C
C     CONSTANTS
C
      INTEGER NBORNAMPS
      PARAMETER (NBORNAMPS=3)
      INTEGER    NLOOPAMPS, NCTAMPS
      PARAMETER (NLOOPAMPS=129, NCTAMPS=85)
      INTEGER    NEXTERNAL
      PARAMETER (NEXTERNAL=4)
      INTEGER    NWAVEFUNCS
      PARAMETER (NWAVEFUNCS=10)
      INTEGER    NCOMB
      PARAMETER (NCOMB=16)
      REAL*16     ZERO
      PARAMETER (ZERO=0E0_16)
      COMPLEX*32 IMAG1
      PARAMETER (IMAG1=(0E0_16,1E0_16))

C
C     ARGUMENTS
C
      REAL*16 P(0:3,NEXTERNAL)
C
C     LOCAL VARIABLES
C
      INTEGER I,J,H
      INTEGER NHEL(NEXTERNAL), IC(NEXTERNAL)
      DATA IC/NEXTERNAL*1/
C
C     FUNCTIONS
C
      LOGICAL ML5_0_IS_HEL_SELECTED
C
C     GLOBAL VARIABLES
C
      INCLUDE 'mp_coupl_same_name.inc'

      INTEGER NTRY
      LOGICAL CHECKPHASE,HELDOUBLECHECKED
      REAL*8 REF
      COMMON/ML5_0_INIT/NTRY,CHECKPHASE,HELDOUBLECHECKED,REF

      LOGICAL GOODHEL(NCOMB)
      LOGICAL GOODAMP(NLOOPAMPS,NCOMB)
      COMMON/ML5_0_FILTERS/GOODAMP,GOODHEL

      INTEGER HELPICKED
      COMMON/ML5_0_HELCHOICE/HELPICKED

      COMPLEX*32 AMP(NBORNAMPS,NCOMB)
      COMMON/ML5_0_MP_AMPS/AMP
      COMPLEX*16 DPAMP(NBORNAMPS,NCOMB)
      COMMON/ML5_0_AMPS/DPAMP
      COMPLEX*32 W(20,NWAVEFUNCS,NCOMB)
      COMMON/ML5_0_MP_WFS/W

      COMPLEX*32 AMPL(3,NCTAMPS)
      COMMON/ML5_0_MP_AMPL/AMPL

      COMPLEX*16 DPW(20,NWAVEFUNCS,NCOMB)
      COMMON/ML5_0_WFCTS/DPW

      COMPLEX*16 DPAMPL(3,NLOOPAMPS)
      LOGICAL S(NLOOPAMPS)
      COMMON/ML5_0_AMPL/DPAMPL,S

      INTEGER HELC(NEXTERNAL,NCOMB)
      COMMON/ML5_0_HELCONFIGS/HELC

      LOGICAL MP_DONE_ONCE
      COMMON/ML5_0_MP_DONE_ONCE/MP_DONE_ONCE

C     This array specify potential special requirements on the
C      helicities to
C     consider. POLARIZATIONS(0,0) is -1 if there is not such
C      requirement.
      INTEGER POLARIZATIONS(0:NEXTERNAL,0:5)
      COMMON/ML5_0_BEAM_POL/POLARIZATIONS

C     ----------
C     BEGIN CODE
C     ---------

      MP_DONE_ONCE=.TRUE.

C     To be on the safe side, we always update the MP params here.
C     It can be redundant as this routine can be called a couple of
C      times for the same PS point during the stability checks.
C     But it is really not time consuming and I would rather be safe.
      CALL MP_UPDATE_AS_PARAM()

      DO H=1,NCOMB
        IF ((HELPICKED.EQ.H).OR.((HELPICKED.EQ.-1)
     $   .AND.((CHECKPHASE.OR..NOT.HELDOUBLECHECKED).OR.GOODHEL(H))))
     $    THEN
C         Handle the possible requirement of specific polarizations
          IF ((.NOT.CHECKPHASE)
     $     .AND.HELDOUBLECHECKED.AND.POLARIZATIONS(0,0)
     $     .EQ.0.AND.(.NOT.ML5_0_IS_HEL_SELECTED(H))) THEN
            CYCLE
          ENDIF
          DO I=1,NEXTERNAL
            NHEL(I)=HELC(I,H)
          ENDDO
          CALL MP_VXXXXX(P(0,1),ZERO,NHEL(1),-1*IC(1),W(1,1,H))
          CALL MP_VXXXXX(P(0,2),ZERO,NHEL(2),-1*IC(2),W(1,2,H))
          CALL MP_OXXXXX(P(0,3),MDL_MT,NHEL(3),+1*IC(3),W(1,3,H))
          CALL MP_IXXXXX(P(0,4),MDL_MT,NHEL(4),-1*IC(4),W(1,4,H))
          CALL MP_VVV1P0_1(W(1,1,H),W(1,2,H),GC_4,ZERO,ZERO,W(1,5,H))
C         Amplitude(s) for born diagram with ID 1
          CALL MP_FFV1_0(W(1,4,H),W(1,3,H),W(1,5,H),GC_5,AMP(1,H))
          CALL MP_FFV1_1(W(1,3,H),W(1,1,H),GC_5,MDL_MT,MDL_WT,W(1,6,H))
C         Amplitude(s) for born diagram with ID 2
          CALL MP_FFV1_0(W(1,4,H),W(1,6,H),W(1,2,H),GC_5,AMP(2,H))
          CALL MP_FFV1_2(W(1,4,H),W(1,1,H),GC_5,MDL_MT,MDL_WT,W(1,7,H))
C         Amplitude(s) for born diagram with ID 3
          CALL MP_FFV1_0(W(1,7,H),W(1,3,H),W(1,2,H),GC_5,AMP(3,H))
          CALL MP_FFV1P0_3(W(1,4,H),W(1,3,H),GC_5,ZERO,ZERO,W(1,8,H))
C         Counter-term amplitude(s) for loop diagram number 4
          CALL MP_R2_GG_1_R2_GG_2_0(W(1,5,H),W(1,8,H),R2_GGG_1
     $     ,R2_GGG_2,AMPL(1,1))
C         Counter-term amplitude(s) for loop diagram number 5
          CALL MP_FFV1_0(W(1,4,H),W(1,3,H),W(1,5,H),R2_GQQ,AMPL(1,2))
          CALL MP_FFV1_0(W(1,4,H),W(1,3,H),W(1,5,H),UV_GQQB_1EPS
     $     ,AMPL(2,3))
          CALL MP_FFV1_0(W(1,4,H),W(1,3,H),W(1,5,H),UV_GQQB_1EPS
     $     ,AMPL(2,4))
          CALL MP_FFV1_0(W(1,4,H),W(1,3,H),W(1,5,H),UV_GQQB_1EPS
     $     ,AMPL(2,5))
          CALL MP_FFV1_0(W(1,4,H),W(1,3,H),W(1,5,H),UV_GQQB_1EPS
     $     ,AMPL(2,6))
          CALL MP_FFV1_0(W(1,4,H),W(1,3,H),W(1,5,H),UV_GQQB_1EPS
     $     ,AMPL(2,7))
          CALL MP_FFV1_0(W(1,4,H),W(1,3,H),W(1,5,H),UV_GQQB_1EPS
     $     ,AMPL(2,8))
          CALL MP_FFV1_0(W(1,4,H),W(1,3,H),W(1,5,H),UV_GQQG_1EPS
     $     ,AMPL(2,9))
          CALL MP_FFV1_0(W(1,4,H),W(1,3,H),W(1,5,H),UV_GQQB,AMPL(1,10))
          CALL MP_FFV1_0(W(1,4,H),W(1,3,H),W(1,5,H),UV_GQQT,AMPL(1,11))
          CALL MP_FFV1_2(W(1,4,H),W(1,2,H),GC_5,MDL_MT,MDL_WT,W(1,9,H))
C         Counter-term amplitude(s) for loop diagram number 7
          CALL MP_R2_QQ_1_R2_QQ_2_0(W(1,9,H),W(1,6,H),R2_QQQ,R2_QQT
     $     ,AMPL(1,12))
          CALL MP_R2_QQ_2_0(W(1,9,H),W(1,6,H),UV_TMASS_1EPS,AMPL(2,13))
          CALL MP_R2_QQ_2_0(W(1,9,H),W(1,6,H),UV_TMASS,AMPL(1,14))
C         Counter-term amplitude(s) for loop diagram number 8
          CALL MP_FFV1_0(W(1,4,H),W(1,6,H),W(1,2,H),R2_GQQ,AMPL(1,15))
          CALL MP_FFV1_0(W(1,4,H),W(1,6,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,16))
          CALL MP_FFV1_0(W(1,4,H),W(1,6,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,17))
          CALL MP_FFV1_0(W(1,4,H),W(1,6,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,18))
          CALL MP_FFV1_0(W(1,4,H),W(1,6,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,19))
          CALL MP_FFV1_0(W(1,4,H),W(1,6,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,20))
          CALL MP_FFV1_0(W(1,4,H),W(1,6,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,21))
          CALL MP_FFV1_0(W(1,4,H),W(1,6,H),W(1,2,H),UV_GQQG_1EPS
     $     ,AMPL(2,22))
          CALL MP_FFV1_0(W(1,4,H),W(1,6,H),W(1,2,H),UV_GQQB,AMPL(1,23))
          CALL MP_FFV1_0(W(1,4,H),W(1,6,H),W(1,2,H),UV_GQQT,AMPL(1,24))
          CALL MP_FFV1_1(W(1,3,H),W(1,2,H),GC_5,MDL_MT,MDL_WT,W(1,10,H)
     $     )
C         Counter-term amplitude(s) for loop diagram number 10
          CALL MP_R2_QQ_1_R2_QQ_2_0(W(1,7,H),W(1,10,H),R2_QQQ,R2_QQT
     $     ,AMPL(1,25))
          CALL MP_R2_QQ_2_0(W(1,7,H),W(1,10,H),UV_TMASS_1EPS,AMPL(2,26)
     $     )
          CALL MP_R2_QQ_2_0(W(1,7,H),W(1,10,H),UV_TMASS,AMPL(1,27))
C         Counter-term amplitude(s) for loop diagram number 11
          CALL MP_FFV1_0(W(1,7,H),W(1,3,H),W(1,2,H),R2_GQQ,AMPL(1,28))
          CALL MP_FFV1_0(W(1,7,H),W(1,3,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,29))
          CALL MP_FFV1_0(W(1,7,H),W(1,3,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,30))
          CALL MP_FFV1_0(W(1,7,H),W(1,3,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,31))
          CALL MP_FFV1_0(W(1,7,H),W(1,3,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,32))
          CALL MP_FFV1_0(W(1,7,H),W(1,3,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,33))
          CALL MP_FFV1_0(W(1,7,H),W(1,3,H),W(1,2,H),UV_GQQB_1EPS
     $     ,AMPL(2,34))
          CALL MP_FFV1_0(W(1,7,H),W(1,3,H),W(1,2,H),UV_GQQG_1EPS
     $     ,AMPL(2,35))
          CALL MP_FFV1_0(W(1,7,H),W(1,3,H),W(1,2,H),UV_GQQB,AMPL(1,36))
          CALL MP_FFV1_0(W(1,7,H),W(1,3,H),W(1,2,H),UV_GQQT,AMPL(1,37))
C         Counter-term amplitude(s) for loop diagram number 13
          CALL MP_FFV1_0(W(1,4,H),W(1,10,H),W(1,1,H),R2_GQQ,AMPL(1,38))
          CALL MP_FFV1_0(W(1,4,H),W(1,10,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,39))
          CALL MP_FFV1_0(W(1,4,H),W(1,10,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,40))
          CALL MP_FFV1_0(W(1,4,H),W(1,10,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,41))
          CALL MP_FFV1_0(W(1,4,H),W(1,10,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,42))
          CALL MP_FFV1_0(W(1,4,H),W(1,10,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,43))
          CALL MP_FFV1_0(W(1,4,H),W(1,10,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,44))
          CALL MP_FFV1_0(W(1,4,H),W(1,10,H),W(1,1,H),UV_GQQG_1EPS
     $     ,AMPL(2,45))
          CALL MP_FFV1_0(W(1,4,H),W(1,10,H),W(1,1,H),UV_GQQB,AMPL(1,46)
     $     )
          CALL MP_FFV1_0(W(1,4,H),W(1,10,H),W(1,1,H),UV_GQQT,AMPL(1,47)
     $     )
C         Counter-term amplitude(s) for loop diagram number 14
          CALL MP_FFV1_0(W(1,9,H),W(1,3,H),W(1,1,H),R2_GQQ,AMPL(1,48))
          CALL MP_FFV1_0(W(1,9,H),W(1,3,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,49))
          CALL MP_FFV1_0(W(1,9,H),W(1,3,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,50))
          CALL MP_FFV1_0(W(1,9,H),W(1,3,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,51))
          CALL MP_FFV1_0(W(1,9,H),W(1,3,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,52))
          CALL MP_FFV1_0(W(1,9,H),W(1,3,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,53))
          CALL MP_FFV1_0(W(1,9,H),W(1,3,H),W(1,1,H),UV_GQQB_1EPS
     $     ,AMPL(2,54))
          CALL MP_FFV1_0(W(1,9,H),W(1,3,H),W(1,1,H),UV_GQQG_1EPS
     $     ,AMPL(2,55))
          CALL MP_FFV1_0(W(1,9,H),W(1,3,H),W(1,1,H),UV_GQQB,AMPL(1,56))
          CALL MP_FFV1_0(W(1,9,H),W(1,3,H),W(1,1,H),UV_GQQT,AMPL(1,57))
C         Counter-term amplitude(s) for loop diagram number 17
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,8,H),R2_3GG,AMPL(1,58))
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,8,H),UV_3GB_1EPS,AMPL(2
     $     ,59))
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,8,H),UV_3GB_1EPS,AMPL(2
     $     ,60))
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,8,H),UV_3GB_1EPS,AMPL(2
     $     ,61))
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,8,H),UV_3GB_1EPS,AMPL(2
     $     ,62))
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,8,H),UV_3GB_1EPS,AMPL(2
     $     ,63))
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,8,H),UV_3GB_1EPS,AMPL(2
     $     ,64))
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,8,H),UV_3GG_1EPS,AMPL(2
     $     ,65))
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,8,H),UV_3GB,AMPL(1,66))
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,8,H),UV_3GT,AMPL(1,67))
C         Counter-term amplitude(s) for loop diagram number 31
          CALL MP_R2_GG_1_0(W(1,5,H),W(1,8,H),R2_GGQ,AMPL(1,68))
          CALL MP_R2_GG_1_0(W(1,5,H),W(1,8,H),R2_GGQ,AMPL(1,69))
          CALL MP_R2_GG_1_0(W(1,5,H),W(1,8,H),R2_GGQ,AMPL(1,70))
          CALL MP_R2_GG_1_0(W(1,5,H),W(1,8,H),R2_GGQ,AMPL(1,71))
C         Counter-term amplitude(s) for loop diagram number 32
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,8,H),R2_3GQ,AMPL(1,72))
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,8,H),R2_3GQ,AMPL(1,73))
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,8,H),R2_3GQ,AMPL(1,74))
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,8,H),R2_3GQ,AMPL(1,75))
C         Counter-term amplitude(s) for loop diagram number 34
          CALL MP_R2_GG_1_R2_GG_3_0(W(1,5,H),W(1,8,H),R2_GGQ,R2_GGB
     $     ,AMPL(1,76))
C         Counter-term amplitude(s) for loop diagram number 35
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,8,H),R2_3GQ,AMPL(1,77))
C         Counter-term amplitude(s) for loop diagram number 37
          CALL MP_R2_GG_1_R2_GG_3_0(W(1,5,H),W(1,8,H),R2_GGQ,R2_GGT
     $     ,AMPL(1,78))
C         Counter-term amplitude(s) for loop diagram number 38
          CALL MP_VVV1_0(W(1,1,H),W(1,2,H),W(1,8,H),R2_3GQ,AMPL(1,79))
C         Amplitude(s) for UVCT diagram with ID 40
          CALL MP_FFV1_0(W(1,4,H),W(1,3,H),W(1,5,H),GC_5,AMPL(2,80))
          AMPL(2,80)=AMPL(2,80)*(4.0D0*UVWFCT_G_1_1EPS+2.0D0
     $     *UVWFCT_B_0_1EPS)
C         Amplitude(s) for UVCT diagram with ID 41
          CALL MP_FFV1_0(W(1,4,H),W(1,3,H),W(1,5,H),GC_5,AMPL(1,81))
          AMPL(1,81)=AMPL(1,81)*(2.0D0*UVWFCT_G_1+2.0D0*UVWFCT_G_2
     $     +2.0D0*UVWFCT_T_0)
C         Amplitude(s) for UVCT diagram with ID 42
          CALL MP_FFV1_0(W(1,4,H),W(1,6,H),W(1,2,H),GC_5,AMPL(2,82))
          AMPL(2,82)=AMPL(2,82)*(4.0D0*UVWFCT_G_1_1EPS+2.0D0
     $     *UVWFCT_B_0_1EPS)
C         Amplitude(s) for UVCT diagram with ID 43
          CALL MP_FFV1_0(W(1,4,H),W(1,6,H),W(1,2,H),GC_5,AMPL(1,83))
          AMPL(1,83)=AMPL(1,83)*(2.0D0*UVWFCT_G_1+2.0D0*UVWFCT_G_2
     $     +2.0D0*UVWFCT_T_0)
C         Amplitude(s) for UVCT diagram with ID 44
          CALL MP_FFV1_0(W(1,7,H),W(1,3,H),W(1,2,H),GC_5,AMPL(2,84))
          AMPL(2,84)=AMPL(2,84)*(4.0D0*UVWFCT_G_1_1EPS+2.0D0
     $     *UVWFCT_B_0_1EPS)
C         Amplitude(s) for UVCT diagram with ID 45
          CALL MP_FFV1_0(W(1,7,H),W(1,3,H),W(1,2,H),GC_5,AMPL(1,85))
          AMPL(1,85)=AMPL(1,85)*(2.0D0*UVWFCT_G_1+2.0D0*UVWFCT_G_2
     $     +2.0D0*UVWFCT_T_0)
C         Copy the qp wfs to the dp ones as they are used to setup the
C          CT calls.
          DO I=1,NWAVEFUNCS
            DO J=1,20
              DPW(J,I,H)=W(J,I,H)
            ENDDO
          ENDDO
C         Same for the counterterms amplitudes
          DO I=1,NCTAMPS
            DO J=1,3
              DPAMPL(J,I)=AMPL(J,I)
              S(I)=.TRUE.
            ENDDO
          ENDDO
          DO I=1,NBORNAMPS
            DPAMP(I,H)=AMP(I,H)
          ENDDO
        ENDIF
      ENDDO

      END


!
!  Copyright 2017 SALMON developers
!
!  Licensed under the Apache License, Version 2.0 (the "License");
!  you may not use this file except in compliance with the License.
!  You may obtain a copy of the License at
!
!      http://www.apache.org/licenses/LICENSE-2.0
!
!  Unless required by applicable law or agreed to in writing, software
!  distributed under the License is distributed on an "AS IS" BASIS,
!  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
!  See the License for the specific language governing permissions and
!  limitations under the License.
!
module read_pslfile_sub

use scf_data

integer,parameter :: Nlps=3, Nlmps=16

integer,allocatable :: Mlps0(:)

integer :: maxMr
real(8) :: rmin_step
real(8) :: rmaxRps

real(8),allocatable :: Zps(:)              ! Pseudo charge
real(8),allocatable :: Rps(:)              ! Core radius
real(8),allocatable :: Mass(:)             ! Atomic weight

real(8) :: rPC,rRC(0:Nlps)

integer,allocatable :: Mr(:)
real(8),allocatable :: step(:)

real(8), allocatable :: upp_f(:,:,:)
real(8), allocatable :: rhopp_f(:,:)
real(8), allocatable :: vpp_f(:,:,:)

real(8), allocatable :: rad_f(:,:)

contains
!==================================================================================================
subroutine read_pslfile
use salmon_parallel, only: nproc_id_global
use salmon_communication, only: comm_is_root
implicit none
integer :: ak
character(256) :: ps_file(MKI)
character(256) :: ps_file_tmp
integer :: ips_type,nlen_psf

allocate( Mlps0(MKI) )
allocate( Mr(MKI) )
allocate( step(MKI) )

allocate( Zps(MKI) )
allocate( Rps(MKI) )
allocate( Mass(MKI) )

maxMr=0
rmin_step=1.d8
rmaxRps=0

do ak=1,MKI

  ps_file_tmp = trim(pseudo_file(ak))
  nlen_psf = len_trim(ps_file_tmp)

  if(ps_file_tmp(nlen_psf+1-8:nlen_psf) == '_rps.dat')then
     ips_type = n_Yabana_Bertsch_psformat
     ps_format(ak) = 'KY'
  else if(ps_file_tmp(nlen_psf+1-6:nlen_psf) == '.pspnc')then
     ips_type = n_ABINIT_psformat
     ps_format(ak) = 'ABINIT'
  else if(ps_file_tmp(nlen_psf+1-4:nlen_psf) == '.cpi')then
     ips_type = n_FHI_psformat
     ps_format(ak) = 'FHI'
  else if(ps_file_tmp(nlen_psf+1-4:nlen_psf) == '.fhi')then
     ips_type = n_ABINITFHI_psformat
     ps_format(ak) = 'ABINITPHI'
  else
     stop 'Unprepared ps_format is required input_pseudopotential_YS'
  end if

  ipsfileform(ak) = ips_type

  if(ipsfileform(ak)==3.or.ipsfileform(ak)==4)then
    Mlps0(ak)=Mlps(ak)
  end if

  if(comm_is_root(nproc_id_global))then
    write(*,*) iZatom(ak)
  end if

  select case( iZatom(ak) )
    case(1) ; Atomname(ak)='H'  ; Mass(ak) = 1.008d0
    case(2) ; Atomname(ak)='He' ; Mass(ak) = 4.003d0
    case(3) ; Atomname(ak)='Li' ; Mass(ak) = 6.941d0
    case(4) ; Atomname(ak)='Be' ; Mass(ak) = 9.012d0
    case(5) ; Atomname(ak)='B'  ; Mass(ak) = 10.811d0
    case(6) ; Atomname(ak)='C'  ; Mass(ak) = 12.011d0
    case(7) ; Atomname(ak)='N'  ; Mass(ak) = 14.007d0
    case(8) ; Atomname(ak)='O'  ; Mass(ak) = 15.999d0
    case(9) ; Atomname(ak)='F'  ; Mass(ak) = 18.998d0
    case(10); Atomname(ak)='Ne' ; Mass(ak) = 20.180d0
    case(11); Atomname(ak)='Na' ; Mass(ak) = 22.990d0
    case(12); Atomname(ak)='Mg' ; Mass(ak) = 24.305d0
    case(13); Atomname(ak)='Al' ; Mass(ak) = 26.982d0
    case(14); Atomname(ak)='Si' ; Mass(ak) = 28.086d0
    case(15); Atomname(ak)='P'  ; Mass(ak) = 30.974d0
    case(16); Atomname(ak)='S'  ; Mass(ak) = 32.065d0
    case(17); Atomname(ak)='Cl' ; Mass(ak) = 35.453d0
    case(18); Atomname(ak)='Ar' ; Mass(ak) = 39.948d0
    case(19); Atomname(ak)='K'  ; Mass(ak) = 39.098d0
    case(20); Atomname(ak)='Ca' ; Mass(ak) = 40.078d0
    case(21); Atomname(ak)='Sc' ; Mass(ak) = 44.956d0
    case(22); Atomname(ak)='Ti' ; Mass(ak) = 47.867d0
    case(23); Atomname(ak)='V'  ; Mass(ak) = 50.942d0
    case(24); Atomname(ak)='Cr' ; Mass(ak) = 51.996d0
    case(25); Atomname(ak)='Mn' ; Mass(ak) = 54.938d0
    case(26); Atomname(ak)='Fe' ; Mass(ak) = 55.845d0
    case(27); Atomname(ak)='Co' ; Mass(ak) = 58.933d0
    case(28); Atomname(ak)='Ni' ; Mass(ak) = 58.693d0
    case(29); Atomname(ak)='Cu' ; Mass(ak) = 63.546d0
    case(30); Atomname(ak)='Zn' ; Mass(ak) = 65.38d0
    case(31); Atomname(ak)='Ga' ; Mass(ak) = 69.723d0
    case(32); Atomname(ak)='Ge' ; Mass(ak) = 72.64d0
    case(33); Atomname(ak)='As' ; Mass(ak) = 74.922d0
    case(34); Atomname(ak)='Se' ; Mass(ak) = 78.96d0
    case(35); Atomname(ak)='Br' ; Mass(ak) = 79.904d0
    case(36); Atomname(ak)='Kr' ; Mass(ak) = 83.798d0
    case(47); Atomname(ak)='Ag' ; Mass(ak) = 107.868d0
    case(79); Atomname(ak)='Au' ; Mass(ak) = 196.967d0
  end select

  select case( ipsfileform(ak) )
    case(1)
      select case( iZatom(ak) )
        case default ; ps_file(ak)=trim(ps_file_tmp) !trim(Atomname(ak))//'_rps.dat'
        case(10,18,36,47,79) ; ps_file(ak)= trim(ps_file_tmp)!trim(Atomname(ak))//'_crps.dat'
      end select
      call read_Mr_YB(ak,ps_file)
    case(2)
      ps_file(ak)=trim(ps_file_tmp) !trim(Atomname(ak))//'.pspnc'
      call read_Mr_ABINIT(ak,ps_file)
    case(3,4)
      ps_file(ak)=trim(ps_file_tmp) !trim(Atomname(ak))//'.cpi' or trim(Atomname(ak))//'.fhi'
      call read_Mr_fhi(ak,ps_file)
  end select
end do

maxMr=maxval(Mr)

allocate(upp_f(0:maxMr,0:Nlps,MKI))
allocate(rhopp_f(0:maxMr,MKI))
allocate(vpp_f(0:maxMr,0:Nlps,MKI))
allocate(rad_f(0:maxMr,MKI) )

do ak=1,MKI
  select case( ipsfileform(ak) )
    case(1)
      call read_psl_YB(ak,ps_file)
    case(2)
      call read_psl_ABINIT(ak,ps_file)
    case(3,4)
      call read_psl_fhi(ak,ps_file)
      call setRps_fhi(ak)
  end select
end do

maxlm=0
do ak=1,MKI
  if(Mlps(ak)>maxlm) maxlm=Mlps(ak)
end do
maxlm=(maxlm+1)**2

end subroutine read_pslfile

!======================================================================
subroutine read_Mr_YB(ak,ps_file)
implicit none
integer :: ak
character(256) :: ps_file(MKI)

open(4,file=ps_file(ak),status='old')
read(4,*) Mr(ak)
close(4)
return

end subroutine read_Mr_YB
!======================================================================
subroutine read_Mr_ABINIT(ak,ps_file)
implicit none
integer :: ak
character(256) :: ps_file(MKI)
real(8) :: zatom, zion, pspdat,pspcod,pspxc,lmaxabinit,lloc,mmax,r2well
character(1) :: dummy_text

open(4,file=ps_file(ak),status='old')
read(4,*) dummy_text
read(4,*) zatom, zion, pspdat
read(4,*) pspcod,pspxc,lmaxabinit,lloc,mmax,r2well
close(4)

Mr(ak)=mmax

end subroutine read_Mr_ABINIT

!======================================================================
subroutine read_Mr_fhi(ak,ps_file)
implicit none
integer :: ak
integer :: i
character(256) :: ps_file(MKI)
character(1) :: dummy_text

open(4,file=ps_file(ak),status='old')
select case( ipsfileform(ak) )
  case(3)
    do i=1,11
      read(4,*) dummy_text
    end do
  case(4)
    do i=1,18
      read(4,*) dummy_text
    end do
end select
read(4,*) Mr(ak)
close(4)

end subroutine read_Mr_fhi

!======================================================================
!======================================================================
subroutine read_psl_YB(ak,ps_file)
implicit none
integer :: ak
character(256) :: ps_file(MKI)
integer :: L
integer :: i
real(8) :: r

open(4,file=ps_file(ak),status='old')

read(4,*) Mr(ak),step(ak),Mlps0(ak),Zps(ak)

step(ak)=step(ak)/a_B

if( Mr(ak) > maxMr ) maxMr=Mr(ak)
if( step(ak) < rmin_step ) rmin_step=step(ak)

if ( Mlps0(ak) < Mlps(ak) ) then
  write(*,*) "too large Mlps value for ",ak," atom"
  stop
end if

rRC=0.d0

read(4,*) rPC, ( rRC(L) , L=0,Mlps0(ak) )
do L=0,Mlps0(ak)
  rRC(L)=rRC(L)/a_B
end do

Rps(ak)=maxval(rRC)
if(maxval( rRC ) > rmaxRps) rmaxRps=maxval( rRC )

do i=0,Mr(ak)
  read(4,*) r,rhopp_f(i,ak),(vpp_f(i,L,ak),L=0,Mlps0(ak))
end do
do i=0,Mr(ak)
  read(4,*) r,(upp_f(i,L,ak),L=0,Mlps0(ak))
end do

do i=0,Mr(ak)
  rhopp_f(i,ak)=rhopp_f(i,ak)*(a_B)**3d0
  do L=0,Mlps0(ak)
    vpp_f(i,L,ak)=vpp_f(i,L,ak)/2d0/Ry
  end do
end do

close(4)

return

end subroutine read_psl_YB

!==================================================================================================
!======================================================================
subroutine read_psl_ABINIT(ak,ps_file)
!See http://www.abinit.org/downloads/psp-links/psp-links/lda_tm
  implicit none
  integer,intent(in) :: ak
  character(256),intent(in) :: ps_file(MKI)
  integer :: i
  real(8) :: rZps
  integer :: ll
  real(8) :: zatom, zion, pspdat,pspcod,pspxc,lmaxabinit,lloc,mmax,r2well,dl
  real(8) :: e99_0,e99_9,nproj,rcpsp,rms,ekb1,ekb2,epsatm,rchrg,fchrg,qchrg
  character(1) :: dummy_text

  step(ak)=0.01d0/a_B

  open(4,file=ps_file(ak),status='old')
  read(4,*) dummy_text
  read(4,*) zatom, zion, pspdat
  rZps = zion
  Zps(ak)=int(rZps+1d-10)
  read(4,*) pspcod,pspxc,lmaxabinit,lloc,mmax,r2well
  Mlps0(ak)=lmaxabinit
  Mr(ak) = mmax - 1

  rRC=0.d0
  do ll=0,Mlps0(ak)
    read(4,*) dl,e99_0,e99_9,nproj,rcpsp
    read(4,*) rms,ekb1,ekb2,epsatm
    rRC(ll) = rcpsp
  end do

  Rps(ak)=maxval(rRC)
  if(maxval( rRC ) > rmaxRps) rmaxRps=maxval( rRC )

  read(4,*) rchrg,fchrg,qchrg
  do ll=0,Mlps0(ak)
    read(4,*) dummy_text
    do i=1,(Mr(ak)+1)/3
      read(4,*) vpp_f(3*(i-1),ll,ak),vpp_f(3*(i-1)+1,ll,ak),vpp_f(3*(i-1)+2,ll,ak)
    end do
  end do
  do ll=0,Mlps0(ak)
    read(4,*) dummy_text
    do i=1,(Mr(ak)+1)/3
      read(4,*) upp_f(3*(i-1),ll,ak),upp_f(3*(i-1)+1,ll,ak),upp_f(3*(i-1)+2,ll,ak)
    end do
  end do
  close(4)

  do i=1,Mr(ak)
    rhopp_f(i,ak)=0.d0
  end do

return

End Subroutine read_psl_ABINIT

!==================================================================================================
!======================================================================
subroutine read_psl_ABINIT_PBE(ak,ps_file)
!See http://www.abinit.org/downloads/psp-links/psp-links/gga-phi
implicit none
integer,intent(in) :: ak
character(256),intent(in) :: ps_file(MKI)
integer :: i,ibox
real(8) :: rZps
integer :: ll
real(8) :: zatom, zion, pspdat,pspcod,pspxc,lmaxabinit,lloc,mmax,r2well
real(8) :: rchrg,fchrg,qchrg
character(1) :: dummy_text
real(8) :: step_tmp(0:3)

step_tmp(0:3)=1.d8

open(4,file=ps_file(ak),status='old')
read(4,*) dummy_text
read(4,*) zatom, zion, pspdat
rZps = zion
Zps(ak)=int(rZps+1d-10)
read(4,*) pspcod,pspxc,lmaxabinit,lloc,mmax,r2well
Mlps0(ak)=lmaxabinit
Mr(ak) = mmax - 1
read(4,*) rchrg,fchrg,qchrg
do i=5,18
  read(4,*)
end do

do ll=0,Mlps0(ak)
  read(4,*) ibox,step_tmp(ll)
  do i=0,Mr(ak)
    read(4,*) ibox,rad_f(i,ak),upp_f(i,ll,ak),vpp_f(i,ll,ak)
  end do
  rad_f(0,ak)=0.d0
end do
close(4)

step(ak) = minval(step_tmp(0:Mlps0(ak)))

do i=1,Mr(ak)
  rhopp_f(i,ak)=0.d0
end do

return

End Subroutine read_psl_ABINIT_PBE

!==================================================================================================
!======================================================================
subroutine read_psl_fhi(ak,ps_file)
implicit none
integer,intent(in) :: ak
character(256),intent(in) :: ps_file(MKI)
integer :: i,ibox
integer :: ll
character(1) :: dummy_text
real(8) :: step_tmp(0:3)

step_tmp(0:3)=1.d8

open(4,file=ps_file(ak),status='old')

call set_Zps(ak)

select case( ipsfileform(ak) )
  case(3)
    do i=1,11
      read(4,*) dummy_text
    end do
  case(4)
    do i=1,18
      read(4,*) dummy_text
    end do
end select

do ll=0,Mlps0(ak)
  read(4,*) ibox,step_tmp(ll)
  Mr(ak)=ibox-1
  do i=0,Mr(ak)
    read(4,*) ibox,rad_f(i,ak),upp_f(i,ll,ak),vpp_f(i,ll,ak)
  end do
  rad_f(0,ak)=0.d0
end do
close(4)

step(ak) = minval(step_tmp(0:Mlps0(ak)))

do i=1,Mr(ak)
  rhopp_f(i,ak)=0.d0
end do

return

End Subroutine read_psl_fhi

!==================================================================================================
!======================================================================
Subroutine setRps_ABINIT_PBE(ak)
implicit none
integer,intent(in) :: ak
real(8) :: rcore(0:4,1:120)
integer :: ll

!rcore list in unit of bohr

rcore(0,1)=1.2763463d0
rcore(1,1)=1.2763463d0
rcore(2,1)=0.3502204d0
rcore(3,1)=1.2763463d0

rcore(0:3,2)=1.0916080d0

rcore(0,3)=2.1818685d0
rcore(1,3)=2.1818685d0
rcore(2,3)=2.4649735d0
rcore(3,3)=2.4649735d0

rcore(0,4)=2.0886093d0
rcore(1,4)=2.0886093d0
rcore(2,4)=2.3596136d0
rcore(3,4)=2.3596136d0

rcore(0:3,5)=1.6708874d0

rcore(0:3,6)=1.4981530d0

rcore(0:3,7)=1.4157818d0

rcore(0:3,8)=1.3995488d0

rcore(0:3,9)=1.3385227d0

rcore(0:3,10)=1.2961595d0

rcore(0,11)=2.7011879d0
rcore(1,11)=2.7011879d0
rcore(2,11)=2.9781169d0
rcore(3,11)=2.9781169d0

rcore(0,12)=1.7907575d0
rcore(1,12)=1.9743480d0
rcore(2,12)=2.1242906d0
rcore(3,12)=2.1242906d0

rcore(0,13)=2.0873202d0
rcore(1:3,13)=2.4760889d0

rcore(0,14)=1.7039185d0
rcore(1,14)=1.8786063d0
rcore(2,14)=2.0212776d0
rcore(3,14)=2.0212776d0

rcore(0,15)=1.6296049d0
rcore(1,15)=1.7966740d0
rcore(2,15)=1.9331230d0
rcore(3,15)=1.9331230d0

rcore(0,16)=1.5654902d0
rcore(1,16)=1.6843819d0
rcore(2,16)=1.8570667d0
rcore(3,16)=1.8570667d0

rcore(0,17)=1.5853006d0
rcore(1,17)=1.5853006d0
rcore(2,17)=1.7909988d0
rcore(3,17)=1.7909988d0

rcore(0,18)=1.4259180d0
rcore(1,18)=1.4259180d0
rcore(2,18)=1.7332789d0
rcore(3,18)=1.7332789d0

rcore(0,19)=2.8781425d0
rcore(1,19)=2.6749895d0
rcore(2,19)=3.1732130d0
rcore(3,19)=3.1732130d0

rcore(0,20)=2.8709747d0
rcore(1,20)=2.8709747d0
rcore(2,20)=2.5412400d0
rcore(3,20)=2.5412400d0

rcore(0,21)=2.7342616d0
rcore(1,21)=2.7342616d0
rcore(2,21)=2.4202286d0
rcore(3,21)=2.4202286d0

rcore(0,22)=2.6099770d0
rcore(1,22)=2.6099770d0
rcore(2,22)=2.3672806d0
rcore(3,22)=2.3672806d0

rcore(0,23)=2.5581633d0
rcore(1,23)=2.5581633d0
rcore(2,23)=2.3775959d0
rcore(3,23)=2.3775959d0

rcore(0,24)=2.5121270d0
rcore(1,24)=2.5121270d0
rcore(2,24)=2.3348091d0
rcore(3,24)=2.3348091d0

rcore(0,25)=2.4116419d0
rcore(1,25)=2.4712095d0
rcore(2,25)=2.2967797d0
rcore(3,25)=2.2967797d0

rcore(0,26)=2.3188865d0
rcore(1,26)=2.4348542d0
rcore(2,26)=2.2084421d0
rcore(3,26)=2.2084421d0

rcore(0,27)=2.2881569d0
rcore(1,27)=2.3446744d0
rcore(2,27)=2.1791761d0
rcore(3,27)=2.1791761d0

rcore(0,28)=2.6173915d0
rcore(1,28)=2.3167811d0
rcore(2,28)=2.1013484d0
rcore(3,28)=2.1013484d0

rcore(0,29)=2.6535190d0
rcore(1,29)=2.2921434d0
rcore(2,29)=2.0790016d0
rcore(3,29)=2.0790016d0

rcore(0,30)=2.3265479d0
rcore(1,30)=2.2704673d0
rcore(2,30)=2.0097016d0
rcore(3,30)=2.6933476d0

rcore(0,31)=2.0925765d0
rcore(1,31)=2.2514979d0
rcore(2,31)=2.0925765d0
rcore(3,31)=2.4823242d0

rcore(0,32)=1.9783190d0
rcore(1,32)=2.7167764d0
rcore(2,32)=2.4641489d0
rcore(3,32)=2.4641489d0

rcore(0,33)=1.9657537d0
rcore(1,33)=2.1672851d0
rcore(2,33)=2.1150435d0
rcore(3,33)=2.4484978d0

rcore(0,34)=1.8619473d0
rcore(1,34)=2.0528364d0
rcore(2,34)=1.8619473d0
rcore(3,34)=2.3764832d0

rcore(0,35)=1.7651496d0
rcore(1,35)=2.0939132d0
rcore(2,35)=1.6405570d0
rcore(3,35)=2.1986300d0

rcore(0,36)=1.7161177d0
rcore(1,36)=1.9866780d0
rcore(2,36)=1.4118000d0
rcore(3,36)=2.0860320d0

rcore(0,47)=2.4191908d0
rcore(1,47)=2.6029170d0
rcore(2,47)=2.4191908d0
rcore(3,47)=2.4191908d0

rcore(0,79)=2.5227012d0
rcore(1,79)=2.6488617d0
rcore(2,79)=2.5227012d0
rcore(3,79)=2.5227012d0

rRC=0.d0
do ll=0,Mlps0(ak)
  rRC(ll) = rcore(ll,iZatom(ak))
end do

Rps(ak)=maxval(rRC)
if(maxval( rRC ) > rmaxRps) rmaxRps=maxval( rRC )

End Subroutine setRps_ABINIT_PBE

!==================================================================================================
!======================================================================
Subroutine set_Zps(ak)
integer :: ak

select case( iZatom(ak) )
  case(1)   ; Zps(ak) = 1.d0
  case(2)   ; Zps(ak) = 2.d0
  case(3)   ; Zps(ak) = 1.d0
  case(4)   ; Zps(ak) = 2.d0
  case(5)   ; Zps(ak) = 3.d0
  case(6)   ; Zps(ak) = 4.d0
  case(7)   ; Zps(ak) = 5.d0
  case(8)   ; Zps(ak) = 6.d0
  case(9)   ; Zps(ak) = 7.d0
  case(10)  ; Zps(ak) = 8.d0
  case(11)  ; Zps(ak) = 1.d0
  case(12)  ; Zps(ak) = 2.d0
  case(13)  ; Zps(ak) = 3.d0
  case(14)  ; Zps(ak) = 4.d0
  case(15)  ; Zps(ak) = 5.d0
  case(16)  ; Zps(ak) = 6.d0
  case(17)  ; Zps(ak) = 7.d0
  case(18)  ; Zps(ak) = 8.d0
  case(19)  ; Zps(ak) = 1.d0
  case(20)  ; Zps(ak) = 2.d0
  case(21)  ; Zps(ak) = 3.d0
  case(22)  ; Zps(ak) = 4.d0
  case(23)  ; Zps(ak) = 5.d0
  case(24)  ; Zps(ak) = 6.d0
  case(25)  ; Zps(ak) = 7.d0
  case(26)  ; Zps(ak) = 8.d0
  case(27)  ; Zps(ak) = 9.d0
  case(28)  ; Zps(ak) = 10.d0
  case(29)  ; Zps(ak) = 11.d0
  case(30)  ; Zps(ak) = 12.d0
  case(31)  ; Zps(ak) = 3.d0
  case(32)  ; Zps(ak) = 4.d0
  case(33)  ; Zps(ak) = 5.d0
  case(34)  ; Zps(ak) = 6.d0
  case(35)  ; Zps(ak) = 7.d0
  case(36)  ; Zps(ak) = 8.d0
  case(37)  ; Zps(ak) = 1.d0
  case(38)  ; Zps(ak) = 2.d0
  case(39)  ; Zps(ak) = 3.d0
  case(40)  ; Zps(ak) = 4.d0
  case(41)  ; Zps(ak) = 5.d0
  case(42)  ; Zps(ak) = 6.d0
  case(43)  ; Zps(ak) = 7.d0
  case(44)  ; Zps(ak) = 8.d0
  case(45)  ; Zps(ak) = 9.d0
  case(46)  ; Zps(ak) = 10.d0
  case(47)  ; Zps(ak) = 11.d0
  case(48)  ; Zps(ak) = 12.d0
  case(49)  ; Zps(ak) = 3.d0
  case(50)  ; Zps(ak) = 4.d0
  case(51)  ; Zps(ak) = 5.d0
  case(52)  ; Zps(ak) = 6.d0
  case(53)  ; Zps(ak) = 7.d0
  case(54)  ; Zps(ak) = 8.d0
  case(55)  ; Zps(ak) = 1.d0
  case(56)  ; Zps(ak) = 2.d0
  case(57)  ; Zps(ak) = 3.d0
  case(58)  ; Zps(ak) = 4.d0
  case(59)  ; Zps(ak) = 5.d0
  case(60)  ; Zps(ak) = 6.d0
  case(61)  ; Zps(ak) = 7.d0
  case(62)  ; Zps(ak) = 8.d0
  case(63)  ; Zps(ak) = 9.d0
  case(64)  ; Zps(ak) = 10.d0
  case(65)  ; Zps(ak) = 11.d0
  case(66)  ; Zps(ak) = 12.d0
  case(67)  ; Zps(ak) = 13.d0
  case(68)  ; Zps(ak) = 14.d0
  case(69)  ; Zps(ak) = 15.d0
  case(70)  ; Zps(ak) = 16.d0
  case(71)  ; Zps(ak) = 17.d0
  case(72)  ; Zps(ak) = 4.d0
  case(73)  ; Zps(ak) = 5.d0
  case(74)  ; Zps(ak) = 6.d0
  case(75)  ; Zps(ak) = 7.d0
  case(76)  ; Zps(ak) = 8.d0
  case(77)  ; Zps(ak) = 9.d0
  case(78)  ; Zps(ak) = 10.d0
  case(79)  ; Zps(ak) = 11.d0
  case(80)  ; Zps(ak) = 12.d0
  case(81)  ; Zps(ak) = 3.d0
  case(82)  ; Zps(ak) = 4.d0
  case(83)  ; Zps(ak) = 5.d0
  case(84)  ; Zps(ak) = 6.d0
  case(85)  ; Zps(ak) = 7.d0
  case(86)  ; Zps(ak) = 8.d0
  case(87)  ; Zps(ak) = 1.d0
  case(88)  ; Zps(ak) = 2.d0
  case(89)  ; Zps(ak) = 3.d0
  case(90)  ; Zps(ak) = 4.d0
  case(91)  ; Zps(ak) = 5.d0
  case(92)  ; Zps(ak) = 6.d0
  case(93)  ; Zps(ak) = 7.d0
  case(94)  ; Zps(ak) = 8.d0
  case(95)  ; Zps(ak) = 9.d0
  case(96)  ; Zps(ak) = 10.d0
  case(97)  ; Zps(ak) = 11.d0
  case(98)  ; Zps(ak) = 12.d0
  case(99)  ; Zps(ak) = 13.d0
  case(100) ; Zps(ak) = 14.d0
  case(101) ; Zps(ak) = 15.d0
  case(102) ; Zps(ak) = 16.d0
  case(103) ; Zps(ak) = 17.d0
  case(104) ; Zps(ak) = 18.d0
  case(105) ; Zps(ak) = 19.d0
  case(106) ; Zps(ak) = 20.d0
  case(107) ; Zps(ak) = 21.d0
  case(108) ; Zps(ak) = 22.d0
  case(109) ; Zps(ak) = 23.d0
end select

End Subroutine set_Zps

!==================================================================================================
!======================================================================
Subroutine setRps_fhi(ak)
implicit none
integer,intent(in) :: ak
real(8) :: rcore(0:4,1:120)
integer :: ll

!rcore list in unit of bohr

rcore(0,1)=  1.2763463d0
rcore(1,1)=  1.2763463d0
rcore(2,1)=  1.2763463d0
rcore(3,1)=  1.2763463d0

rcore(0,2)=  1.0916080d0
rcore(1,2)=  1.0916080d0
rcore(2,2)=  1.0916080d0
rcore(3,2)=  1.0916080d0

rcore(0,3)=  2.1818685d0
rcore(1,3)=  2.1818685d0
rcore(2,3)=  2.4649735d0
rcore(3,3)=  2.4649735d0

rcore(0,4)=  2.0886093d0
rcore(1,4)=  2.0886093d0
rcore(2,4)=  2.3596136d0
rcore(3,4)=  2.3596136d0

rcore(0,5)=  1.6708874d0
rcore(1,5)=  1.6708874d0
rcore(2,5)=  1.6708874d0
rcore(3,5)=  1.6708874d0

rcore(0,6)=  1.4981530d0
rcore(1,6)=  1.4981530d0
rcore(2,6)=  1.4981530d0
rcore(3,6)=  1.4981530d0

rcore(0,7)=  1.4157818d0
rcore(1,7)=  1.4157818d0
rcore(2,7)=  1.4157818d0
rcore(3,7)=  1.4157818d0

rcore(0,8)=  1.3995488d0
rcore(1,8)=  1.3995488d0
rcore(2,8)=  1.3995488d0
rcore(3,8)=  1.3995488d0

rcore(0,9)=  1.3385227d0
rcore(1,9)=  1.3385227d0
rcore(2,9)=  1.3385227d0
rcore(3,9)=  1.3385227d0

rcore(0,10)=  1.2961595d0
rcore(1,10)=  1.2961595d0
rcore(2,10)=  1.2961595d0
rcore(3,10)=  1.2961595d0

rcore(0,11)=  2.7011879d0
rcore(1,11)=  2.7011879d0
rcore(2,11)=  2.9781169d0
rcore(3,11)=  2.9781169d0

rcore(0,12)=  2.0873202d0
rcore(1,12)=  2.4760889d0
rcore(2,12)=  2.4760889d0
rcore(3,12)=  2.4760889d0

rcore(0,13)=  1.7907575d0
rcore(1,13)=  1.9743480d0
rcore(2,13)=  2.1242906d0
rcore(3,13)=  2.1242906d0

rcore(0,14)=  1.7039185d0
rcore(1,14)=  1.8786063d0
rcore(2,14)=  2.0212776d0
rcore(3,14)=  2.0212776d0

rcore(0,15)=  1.6296049d0
rcore(1,15)=  1.7966740d0
rcore(2,15)=  1.9331230d0
rcore(3,15)=  1.9331230d0

rcore(0,16)=  1.5654902d0
rcore(1,16)=  1.6843819d0
rcore(2,16)=  1.8570667d0
rcore(3,16)=  1.8570667d0

rcore(0,17)=  1.5853006d0
rcore(1,17)=  1.5853006d0
rcore(2,17)=  1.7909988d0
rcore(3,17)=  1.7909988d0

rcore(0,18)=  1.4259180d0
rcore(1,18)=  1.4259180d0
rcore(2,18)=  1.7332789d0
rcore(3,18)=  1.7332789d0

rcore(0,19)=  3.2515914d0
rcore(1,19)=  3.1732130d0
rcore(2,19)=  3.1732130d0
rcore(3,19)=  3.1732130d0

rcore(0,20)=  2.7342353d0
rcore(1,20)=  2.9418878d0
rcore(2,20)=  2.5412400d0
rcore(3,20)=  2.5412400d0

rcore(0,21)=  2.7342616d0
rcore(1,21)=  2.7342616d0
rcore(2,21)=  2.4202286d0
rcore(3,21)=  2.4202286d0

rcore(0,22)=  2.6099770d0
rcore(1,22)=  2.6099770d0
rcore(2,22)=  2.3672806d0
rcore(3,22)=  2.3672806d0

rcore(0,23)=  2.5581633d0
rcore(1,23)=  2.5581633d0
rcore(2,23)=  2.3775959d0
rcore(3,23)=  2.3775959d0

rcore(0,24)=  2.5121270d0
rcore(1,24)=  2.5121270d0
rcore(2,24)=  2.3348091d0
rcore(3,24)=  2.3348091d0

rcore(0,25)=  2.4116419d0
rcore(1,25)=  2.4712095d0
rcore(2,25)=  2.2967797d0
rcore(3,25)=  2.2967797d0

rcore(0,26)=  2.3188865d0
rcore(1,26)=  2.4348542d0
rcore(2,26)=  2.2084421d0
rcore(3,26)=  2.2084421d0

rcore(0,27)=  2.2881569d0
rcore(1,27)=  2.3446744d0
rcore(2,27)=  2.1791761d0
rcore(3,27)=  2.1791761d0

rcore(0,28)=  2.2064370d0
rcore(1,28)=  2.3167811d0
rcore(2,28)=  2.1013484d0
rcore(3,28)=  2.1013484d0

rcore(0,29)=  2.0790016d0
rcore(1,29)=  2.2921434d0
rcore(2,29)=  2.0790016d0
rcore(3,29)=  2.0790016d0

rcore(0,30)=  2.0097016d0
rcore(1,30)=  2.2704673d0
rcore(2,30)=  2.0097016d0
rcore(3,30)=  2.6933476d0

rcore(0,31)=  2.0925765d0
rcore(1,31)=  2.2514979d0
rcore(2,31)=  2.0925765d0
rcore(3,31)=  2.4823242d0

rcore(0,32)=  1.9783190d0
rcore(1,32)=  2.1811386d0
rcore(2,32)=  2.4641489d0
rcore(3,32)=  2.4641489d0

rcore(0,33)=  1.9657537d0
rcore(1,33)=  2.1672851d0
rcore(2,33)=  2.4484978d0
rcore(3,33)=  2.4484978d0

rcore(0,34)=  1.8619473d0
rcore(1,34)=  2.0528364d0
rcore(2,34)=  2.3764832d0
rcore(3,34)=  2.3764832d0

rcore(0,35)=  1.7651496d0
rcore(1,35)=  2.0939132d0
rcore(2,35)=  2.1986300d0
rcore(3,35)=  2.1986300d0

rcore(0,36)=  1.7161177d0
rcore(1,36)=  1.9866780d0
rcore(2,36)=  2.0860320d0
rcore(3,36)=  2.0860320d0

rcore(0,38)=  3.4638956d0
rcore(1,38)=  3.9133480d0
rcore(2,38)=  3.4638956d0
rcore(3,38)=  3.4638956d0

rcore(0,39)=  3.1368487d0
rcore(1,39)=  3.1368487d0
rcore(2,39)=  2.6443347d0
rcore(3,39)=  2.6443347d0

rcore(0,40)=  2.9127601d0
rcore(1,40)=  3.0584275d0
rcore(2,40)=  2.5782263d0
rcore(3,40)=  2.5782263d0

rcore(0,41)=  2.7732187d0
rcore(1,41)=  2.9838317d0
rcore(2,41)=  2.4547114d0
rcore(3,41)=  2.4547114d0

rcore(0,42)=  2.6419339d0
rcore(1,42)=  2.7740572d0
rcore(2,42)=  2.3962659d0
rcore(3,42)=  2.3962659d0

rcore(0,43)=  2.5804936d0
rcore(1,43)=  2.6442318d0
rcore(2,43)=  2.6442318d0
rcore(3,43)=  2.6442318d0

rcore(0,44)=  2.4610579d0
rcore(1,44)=  2.5841356d0
rcore(2,44)=  2.4017350d0
rcore(3,44)=  2.4017350d0

rcore(0,45)=  2.4658050d0
rcore(1,45)=  2.5891201d0
rcore(2,45)=  2.3483631d0
rcore(3,45)=  2.3483631d0

rcore(0,46)=  2.4717819d0
rcore(1,46)=  2.5953959d0
rcore(2,46)=  2.3540553d0
rcore(3,46)=  2.3540553d0

rcore(0,47)=  2.4191908d0
rcore(1,47)=  2.6029170d0
rcore(2,47)=  2.4191908d0
rcore(3,47)=  2.4191908d0

rcore(0,48)=  2.3687909d0
rcore(1,48)=  2.4872544d0
rcore(2,48)=  2.3116922d0
rcore(3,48)=  2.3116922d0

rcore(0,49)=  2.4364941d0
rcore(1,49)=  2.4966755d0
rcore(2,49)=  2.3204483d0
rcore(3,49)=  2.3204483d0

rcore(0,50)=  2.3302081d0
rcore(1,50)=  2.3877642d0
rcore(2,50)=  2.6975849d0
rcore(3,50)=  2.6975849d0

rcore(0,51)=  2.2294503d0
rcore(1,51)=  2.3409453d0
rcore(2,51)=  2.6446911d0
rcore(3,51)=  2.6446911d0

rcore(0,52)=  2.2405847d0
rcore(1,52)=  2.2959271d0
rcore(2,52)=  2.6578993d0
rcore(3,52)=  2.6578993d0

rcore(0,53)=  2.1453201d0
rcore(1,53)=  2.1983095d0
rcore(2,53)=  2.6721617d0
rcore(3,53)=  2.6721617d0

rcore(0,54)=  2.0548375d0
rcore(1,54)=  2.0548375d0
rcore(2,54)=  2.4977638d0
rcore(3,54)=  2.4977638d0

rcore(0,55)=  3.9950035d0
rcore(1,55)=  4.4045760d0
rcore(2,55)=  3.9950035d0
rcore(3,55)=  3.9950035d0

rcore(0,56)=  4.2216481d0
rcore(1,56)=  4.8872263d0
rcore(2,56)=  3.9236641d0
rcore(3,56)=  3.9236641d0

rcore(0,58)=  3.4360930d0
rcore(1,58)=  3.9778218d0
rcore(2,58)=  2.9681408d0
rcore(3,58)=  3.1935572d0

rcore(0,59)=  2.4597067d0
rcore(1,59)=  2.4597067d0
rcore(2,59)=  2.4597067d0
rcore(3,59)=  2.4597067d0

rcore(0,61)=  2.4980376d0
rcore(1,61)=  2.4980376d0
rcore(2,61)=  2.4980376d0
rcore(3,61)=  2.4980376d0

rcore(0,62)=  2.4577467d0
rcore(1,62)=  2.4577467d0
rcore(2,62)=  2.4577467d0
rcore(3,62)=  2.4577467d0

rcore(0,63)=  2.4784776d0
rcore(1,63)=  2.4784776d0
rcore(2,63)=  2.4784776d0
rcore(3,63)=  2.4784776d0

rcore(0,64)=  2.4397514d0
rcore(1,64)=  2.4397514d0
rcore(2,64)=  2.4397514d0
rcore(3,64)=  2.4397514d0

rcore(0,65)=  2.4615515d0
rcore(1,65)=  2.4615515d0
rcore(2,65)=  2.4615515d0
rcore(3,65)=  2.4615515d0

rcore(0,66)=  2.4841343d0
rcore(1,66)=  2.4841343d0
rcore(2,66)=  2.4841343d0
rcore(3,66)=  2.4841343d0

rcore(0,68)=  2.4706250d0
rcore(1,68)=  2.4706250d0
rcore(2,68)=  2.4706250d0
rcore(3,68)=  2.4706250d0

rcore(0,72)=  2.7679638d0
rcore(1,72)=  2.8363325d0
rcore(2,72)=  2.4500600d0
rcore(3,72)=  2.4500600d0

rcore(0,73)=  2.6642398d0
rcore(1,73)=  2.7974787d0
rcore(2,73)=  2.4164976d0
rcore(3,73)=  2.4164976d0

rcore(0,74)=  2.6282365d0
rcore(1,74)=  2.6931540d0
rcore(2,74)=  2.4427231d0
rcore(3,74)=  2.4427231d0

rcore(0,75)=  2.5931934d0
rcore(1,75)=  2.6572453d0
rcore(2,75)=  2.4101535d0
rcore(3,75)=  2.4101535d0

rcore(0,76)=  2.5590724d0
rcore(1,76)=  2.6870519d0
rcore(2,76)=  2.4973870d0
rcore(3,76)=  2.4973870d0

rcore(0,77)=  2.5258377d0
rcore(1,77)=  2.6521551d0
rcore(2,77)=  2.5258377d0
rcore(3,77)=  2.5258377d0

rcore(0,78)=  2.4934552d0
rcore(1,78)=  2.6181531d0
rcore(2,78)=  2.4934552d0
rcore(3,78)=  2.4934552d0

rcore(0,79)=  2.5227012d0
rcore(1,79)=  2.6488617d0
rcore(2,79)=  2.5227012d0
rcore(3,79)=  2.5227012d0

rcore(0,80)=  2.4911674d0
rcore(1,80)=  2.4911674d0
rcore(2,80)=  2.3725176d0
rcore(3,80)=  2.3725176d0

rcore(0,81)=  2.4604123d0
rcore(1,81)=  2.4604123d0
rcore(2,81)=  2.4604123d0
rcore(3,81)=  2.4604123d0

rcore(0,82)=  2.3718232d0
rcore(1,82)=  2.3718232d0
rcore(2,82)=  2.6149854d0
rcore(3,82)=  2.6149854d0

rcore(0,83)=  2.2867640d0
rcore(1,83)=  2.2867640d0
rcore(2,83)=  2.4604330d0
rcore(3,83)=  2.4604330d0

rcore(0,84)=  2.4911914d0
rcore(1,84)=  2.4911914d0
rcore(2,84)=  2.4911914d0
rcore(3,84)=  2.4911914d0

rcore(0,85)=  2.4618832d0
rcore(1,85)=  2.4618832d0
rcore(2,85)=  2.4618832d0
rcore(3,85)=  2.4618832d0

rcore(0,86)=  2.4933581d0
rcore(1,86)=  2.4933581d0
rcore(2,86)=  2.4933581d0
rcore(3,86)=  2.4933581d0

rcore(0,96)=  2.4626281d0
rcore(1,96)=  2.4626281d0
rcore(2,96)=  2.4626281d0
rcore(3,96)=  2.4626281d0

rcore(0,97)=  2.4974400d0
rcore(1,97)=  2.4974400d0
rcore(2,97)=  2.4974400d0
rcore(3,97)=  2.4974400d0

rcore(0,99)=  2.4469867d0
rcore(1,99)=  2.4469867d0
rcore(2,99)=  2.4469867d0
rcore(3,99)=  2.4469867d0

rcore(0,100)=  2.4823529d0
rcore(1,100)=  2.4823529d0
rcore(2,100)=  2.4823529d0
rcore(3,100)=  2.4823529d0

rcore(0,103)=  2.4695797d0
rcore(1,103)=  2.4695797d0
rcore(2,103)=  2.4695797d0
rcore(3,103)=  2.4695797d0

rcore(0,104)=  2.4458337d0
rcore(1,104)=  2.4458337d0
rcore(2,104)=  2.4458337d0
rcore(3,104)=  2.4458337d0

rcore(0,105)=  2.4823768d0
rcore(1,105)=  2.4823768d0
rcore(2,105)=  2.4823768d0
rcore(3,105)=  2.4823768d0

rcore(0,106)=  2.4589582d0
rcore(1,106)=  2.4589582d0
rcore(2,106)=  2.4589582d0
rcore(3,106)=  2.4589582d0

rcore(0,107)=  2.4961459d0
rcore(1,107)=  2.4961459d0
rcore(2,107)=  2.4961459d0
rcore(3,107)=  2.4961459d0

rcore(0,108)=  2.4730334d0
rcore(1,108)=  2.4730334d0
rcore(2,108)=  2.4730334d0
rcore(3,108)=  2.4730334d0

rcore(0,109)=  2.4503450d0
rcore(1,109)=  2.4503450d0
rcore(2,109)=  2.4503450d0
rcore(3,109)=  2.4503450d0

select case( iZatom(ak) )
  case(37,57,60,67,69-71,87-95,98,101,102) ; write(*,*) "fhi file is not prepared yet." ; stop
end select

rRC=0.d0
do ll=0,Mlps(ak)
  rRC(ll) = rcore(ll,iZatom(ak))
end do

Rps(ak)=maxval(rRC)
if(maxval( rRC ) > rmaxRps) rmaxRps=maxval( rRC )

End Subroutine setRps_fhi

end module read_pslfile_sub