text:
!   nls.f90
!   Computational core of the project. Define routines that solve Non-Linear Schrodinger equation.
!   (c) Daniel Bershatsky, 2015-2016
!   See LISENCE for details.
module nls
    implicit none

    private

    integer, parameter :: sp = selected_real_kind(6, 37)
    integer, parameter :: dp = selected_real_kind(15, 307)

    public :: version
    public :: make_banded_matrix
    public :: clear_first_row_of_derivative
    public :: divide_derivative_on_radius
    public :: make_laplacian_o3, make_laplacian_o5, make_laplacian_o7, make_laplacian
    public :: make_laplacian_2d_o3, make_laplacian_2d_o5, make_laplacian_2d_o7, make_laplacian_2d
    public :: rbbmv_o3, rbbmv_o5, rbbmv_o7, rbbmv, rgbmv
    public :: revervoir, revervoir_2d
    public :: hamiltonian, hamiltonian_2d
    public :: runge_kutta, runge_kutta_2d
    public :: solve_nls_1d, solve_nls_2d, solve_nls
    public :: chemical_potential_1d, chemical_potential_2d

contains

    !   \brief Return current version of library in format `major.minor.patch`.
    subroutine version(major, minor, patch)
        implicit none

        integer, intent(out) :: major, minor, patch

        major = 0
        minor = 2
        patch = 0
    end subroutine version

    !   \brief Make banded matrix from a row. Banded matrix representation corresponds to BLAS-2 documentation. Memory
    !   usage O(nm), time complexity O(nm).
    !
    !   \param[in] n
    !   \verbatim
    !       Size of square banded matrix `mat`.
    !   \endverbatim
    !
    !   \param[in] m
    !   \verbatim
    !       Length of `row`.
    !   \endverbatim
    !
    !   \param[in] row
    !   \verbatim
    !       Row that composes banded matrix.
    !   \endverbatim
    !
    !   \param[out] mat
    !   \verbatim
    !       Banded matrix as a Fortran array of shape (m, n).
    !   \endverbatim
    pure subroutine make_banded_matrix(n, m, row, mat)
        implicit none

        integer, parameter :: sp = selected_real_kind(6, 37)
        integer, intent(in) :: n, m
        real(sp), intent(in), dimension(m) :: row
        real(sp), intent(out), dimension(m, n) :: mat

        integer :: j, k, l, r

        k = (m - 1) / 2
        l = k + 1
        r = n - k

        ! Fill left triangle part of band
        do j = 1, k
            mat(k + 2 - j + 0:, j) = row(k + 2 - j:)
            mat(:k + 2 - j - 1, j) = 0 ! can be removed
        end do

        ! Fill body
        do j = l, r
            mat(:, j) = row(:)
        end do

        ! Fill right triangle part of band
        do j = 1, k
            mat(:m - j + 0, n - k + j) = row(:m - j)
            mat(m - j + 1:, n - k + j) = 0 ! can be removed
        end do
    end subroutine make_banded_matrix

    pure subroutine clear_first_row_of_derivative(n, m, L1)
        implicit none

        integer, parameter :: sp = selected_real_kind(6, 37)
        integer, intent(in) :: n
        integer, intent(in) :: m
        real(sp), intent(out), dimension(m, n) :: L1

        integer :: i

        ! Here (m + 1) / 2 is middle row of banded materix L1
        do i = 1, (m + 1) / 2
            L1(i, (m + 1) / 2 - i + 1) = 0.0
        end do
    end subroutine clear_first_row_of_derivative

    pure subroutine divide_derivative_on_radius(n, m, h, L1)
        implicit none

        integer, parameter :: sp = selected_real_kind(6, 37)
        integer, intent(in) :: n
        integer, intent(in) :: m
        real(sp), intent(in) :: h
        real(sp), intent(out), dimension(m, n) :: L1

        integer :: i, j, row

        ! Row k = i + j - k0: k = 0,1,2,...,n - 1, where k0: m = 2 k0 - 3. Cause is sum (i + j) is invarian of row.
        do i = 1, m
            do j = 1, n
                row = i + j - (m + 3) / 2

                if (row > 0) then
                    L1(i, j) = L1(i, j) / (row * h)
                end if
            end do
        end do
    end subroutine divide_derivative_on_radius

    pure subroutine make_laplacian_o3(n, h, op)
        implicit none

        integer, parameter :: sp = selected_real_kind(6, 37)
        integer, intent(in) :: n
        integer, parameter :: m = 3
        real(sp), intent(in) :: h
        real(sp), intent(out), dimension(m, n) :: op

        real(sp), dimension(m) :: row1d, row2d
        real(sp), dimension(m, n) :: L1, L2
        real(sp) :: dx1, dx2

        dx1 = 2 * h ** 1
        dx2 = 1 * h ** 2

        row1d = (/ 1, 0, -1 /) / dx1
        row2d = (/ 1, -2, 1 /) / dx2

        call make_banded_matrix(n, m, row1d, L1)
        call make_banded_matrix(n, m, row2d, L2)

        L2(2, 1) = 2 * L2(2, 1)  ! node #0, row #0
        L2(1, 2) = 4 * L2(1, 2)  ! node #0, row #0

        call clear_first_row_of_derivative(n, m, L1)
        call divide_derivative_on_radius(n, m, h, L1)

        op = L1 + L2
    end subroutine make_laplacian_o3

    subroutine make_laplacian_o5(n, h, op)
        implicit none

        integer, parameter :: sp = selected_real_kind(6, 37)
        integer, intent(in) :: n
        integer, parameter :: m = 5
        real(sp), intent(in) :: h
        real(sp), intent(out), dimension(m, n) :: op

        real(sp), dimension(m) :: row1d, row2d
        real(sp), dimension(m, n) :: L1, L2
        real(sp) :: dx1, dx2

        dx1 = 12 * h ** 1
        dx2 = 12 * h ** 2

        row1d = (/ -1, 8, 0, -8, 1/) / dx1
        row2d = (/ -1, 16, -30, 16, -1 /) / dx2

        call make_banded_matrix(n, m, row1d, L1)
        call make_banded_matrix(n, m, row2d, L2)

        L1(3, 2) = L1(3, 2) + 1.0 / dx1
        L2(3, 2) = L2(3, 2) - 1.0 / dx2
        L2(3, 1) = 2 * L2(3, 1)
        L2(2, 2) = 4 * L2(2, 2)
        L2(1, 3) = 4 * L2(1, 3)

        call clear_first_row_of_derivative(n, m, L1)
        call divide_derivative_on_radius(n, m, h, L1)

        op = L1 + L2
    end subroutine make_laplacian_o5

    pure subroutine make_laplacian_o7(n, h, op)
        implicit none

        integer, parameter :: sp = selected_real_kind(6, 37)
        integer, intent(in) :: n
        integer, parameter :: m = 7
        real(sp), intent(in) :: h
        real(sp), intent(out), dimension(m, n) :: op

        real(sp), dimension(m) :: row1d, row2d
        real(sp), dimension(m, n) :: L1, L2
        real(sp) :: dx1, dx2

        dx1 = 60 * h ** 1
        dx2 = 180 * h ** 2

        row1d = (/ 1, -9, 45, 0, -45, 9, -1/) / dx1
        row2d = (/ 2, -27, 270, -490, 270, -27, 2 /) / dx2

        call make_banded_matrix(n, m, row1d, L1)
        call make_banded_matrix(n, m, row2d, L2)

        ! The second row
        L1(5, 1) = L1(5, 1)
        L1(4, 2) = L1(4, 2) + 9.0 / dx1
        L1(3, 3) = L1(3, 3) - 1.0 / dx1
        L1(2, 4) = L1(2, 4)
        L1(1, 5) = L1(1, 5)

        ! The third row
        L1(6, 1) = L1(6, 1)
        L1(5, 2) = L1(5, 2) - 1.0 / dx1
        L1(4, 3) = L1(4, 3)
        L1(3, 4) = L1(3, 4)
        L1(2, 5) = L1(2, 5)
        L1(1, 6) = L1(1, 6)

        ! The first row
        L2(4, 1) = 2 * L2(4, 1)
        L2(3, 2) = 4 * L2(3, 2)
        L2(2, 3) = 4 * L2(2, 3)
        L2(1, 4) = 4 * L2(1, 4)

        ! The second row
        L2(5, 1) = L2(5, 1)
        L2(4, 2) = L2(4, 2) - 27.0 / dx2
        L2(3, 3) = L2(3, 3) + 2.0 / dx2
        L2(2, 4) = L2(2, 4)
        L2(1, 5) = L2(1, 5)

        ! The third row
        L2(6, 1) = L2(6, 1)
        L2(5, 2) = L2(5, 2) + 2.0 / dx2
        L2(4, 3) = L2(4, 3)
        L2(3, 4) = L2(3, 4)
        L2(2, 5) = L2(2, 5)
        L2(1, 6) = L2(1, 6)

        call clear_first_row_of_derivative(n, m, L1)
        call divide_derivative_on_radius(n, m, h, L1)

        op = L1 + L2
    end subroutine make_laplacian_o7

    !   \brief Make laplacian matrix that approximates laplacian operator in axial symmetric case on a given grid and
    !   given approximation. TODO: implement approximation of order 3 and 7.
    !
    !   \param[in] n
    !   \verbatim
    !       Size of square banded matrix `op`.
    !   \endverbatim
    !
    !   \param[in] m
    !   \verbatim
    !       Order of approximation.
    !   \endverbatim
    !
    !   \param[out] op
    !   \verbatim
    !       Banded matrix of size `n` that approximates laplacian operator in order `m`.
    !   \endverbatim
    subroutine make_laplacian(n, m, h, op)
        implicit none

        integer, parameter :: sp = selected_real_kind(6, 37)
        integer, intent(in) :: n
        integer, intent(in) :: m ! order
        real(sp), intent(in) :: h
        real(sp), intent(out), dimension(m, n) :: op

        if (m == 3) then
            call make_laplacian_o3(n, h, op)
        else if (m == 5) then
            call make_laplacian_o5(n, h, op)
        else if (m == 7) then
            call make_laplacian_o7(n, h, op)
        else
        end if
    end subroutine make_laplacian

    subroutine make_laplacian_2d_o3(n, h, blocks, orders)
        implicit none

        integer, parameter :: m = 3
        integer, intent(in) :: n
        integer, intent(out), dimension(m) :: orders
        real(sp), intent(in) :: h
        real(sp), intent(out), dimension(n, 2 * m - 1) :: blocks

        real(sp) :: dx2
        real(sp), dimension(1) :: left, right
        real(sp), dimension(3) :: middle

        dx2 = 1.0 * h ** 2

        left = (/ 1 /) / dx2
        middle = (/ 1, -4, 1 /) / dx2
        right = (/ 1 /) / dx2

        call make_banded_matrix(n, 1, left, blocks(:, 1:1))
        call make_banded_matrix(n, 3, middle, blocks(:, 2:4))
        call make_banded_matrix(n, 1, right, blocks(:, 5:5))

        orders = (/ 0, 1, 0 /)
    end subroutine make_laplacian_2d_o3

    subroutine make_laplacian_2d_o5(n, h, blocks, orders)
        implicit none

        integer, parameter :: m = 5
        integer, intent(in) :: n
        integer, intent(out), dimension(m) :: orders
        real(sp), intent(in) :: h
        real(sp), intent(out), dimension(n, 2 * m - 1) :: blocks

        real(sp) :: dx2
        real(sp), dimension(1) :: left1, left2, right1, right2
        real(sp), dimension(m) :: middle

        dx2 = 12 * h ** 2

        left2 = (/ -1 /) / dx2
        left1 = (/ 16 /) / dx2
        middle = (/ -1, 16, -60, 16, -1 /) / dx2
        right1 = (/ 16 /) / dx2
        right2 = (/ -1 /) / dx2

        call make_banded_matrix(n, 1, left2, blocks(:, 1:1))
        call make_banded_matrix(n, 1, left1, blocks(:, 2:2))
        call make_banded_matrix(n, m, middle, blocks(:, 3:7))
        call make_banded_matrix(n, 1, right1, blocks(:, 8:8))
        call make_banded_matrix(n, 1, right2, blocks(:, 9:9))

        orders = (/ 0, 0, 2, 0, 0 /)
    end subroutine make_laplacian_2d_o5

    subroutine make_laplacian_2d_o7(n, h, blocks, orders)
        implicit none

        integer, parameter :: m = 7
        integer, intent(in) :: n
        integer, intent(out), dimension(m) :: orders
        real(sp), intent(in) :: h
        real(sp), intent(out), dimension(m, 2 * m - 1) :: blocks

        real(sp) :: dx2
        real(sp), dimension(1) :: left3, left2, left1, right1, right2, right3
        real(sp), dimension(m) :: middle

        dx2 = 180.0 * h ** 2

        left3 = (/ 2 /) / dx2
        left2 = (/ -27 /) / dx2
        left1 = (/ 270 /) / dx2
        middle = (/ 2, -27, 270, -980, 270, -27, 2 /) / dx2
        right1 = left1
        right2 = left2
        right3 = left3

        call make_banded_matrix(n, 1, left3, blocks(:, 1:1))
        call make_banded_matrix(n, 1, left2, blocks(:, 2:2))
        call make_banded_matrix(n, 1, left1, blocks(:, 3:3))
        call make_banded_matrix(n, m, middle, blocks(:, 4:10))
        call make_banded_matrix(n, 1, right1, blocks(:, 11:11))
        call make_banded_matrix(n, 1, right2, blocks(:, 12:12))
        call make_banded_matrix(n, 1, right3, blocks(:, 13:13))

        orders = (/ 0, 0, 0, 3, 0, 0, 0 /)
    end subroutine make_laplacian_2d_o7

    ! TODO: implement high order approximations.
    subroutine make_laplacian_2d(n, m, h, blocks, orders)
        implicit none

        integer, intent(in) :: n, m
        integer, intent(out), dimension(m) :: orders
        real(sp), intent(in) :: h
        real(sp), intent(out), dimension(n, 2 * m - 1) :: blocks

        if (m == 3) then
            call make_laplacian_2d_o3(n, h, blocks, orders)
        else if (m == 5) then
            call make_laplacian_2d_o5(n, h, blocks, orders)
        else if (m == 7) then
            call make_laplacian_2d_o7(n, h, blocks, orders)
        end if
    end subroutine make_laplacian_2d

    !   \brief Perform one of the matrix-vector operations   y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y in
    !   case of block band matrix.
    !   TODO: support high-order approximations.
    subroutine rbbmv_o3(x, y, sign, blocks, ms, n)
        implicit none

        integer, parameter :: m = 3
        integer, intent(in) :: n
        integer, intent(in), dimension(m) :: ms
        real(sp), intent(in) :: sign
        real(sp), intent(in), dimension(n * n) :: x
        real(sp), intent(inout), dimension(n * n) :: y
        real(sp), intent(in), dimension(n, 2 * m - 1) :: blocks

        integer :: i

        do i = 2, n
            call rgbmv(x((i - 2) * n + 1:(i - 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(:, 1:1), ms(1), n)
        end do

        do i = 1, n
            call rgbmv(x((i - 1) * n + 1:(i + 0) * n), y((i - 1) * n + 1:i * n), sign, blocks(:, 2:4), ms(2), n)
        end do

        do i = 1, n - 1
            call rgbmv(x((i + 0) * n + 1:(i + 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(:, 5:5), ms(3), n)
        end do
    end subroutine rbbmv_o3

    subroutine rbbmv_o5(x, y, sign, blocks, ms, n)
        implicit none

        integer, parameter :: m = 5
        integer, intent(in) :: n
        integer, intent(in), dimension(m) :: ms
        real(sp), intent(in) :: sign
        real(sp), intent(in), dimension(n * n) :: x
        real(sp), intent(inout), dimension(n * n) :: y
        real(sp), intent(in), dimension(n, 2 * m - 1) :: blocks

        integer :: i

        do i = 1, n - 2
            call rgbmv(x((i + 1) * n + 1:(i + 2) * n), y((i - 1) * n + 1:i * n), sign, blocks(:, 9:9), ms(5), n)
        end do

        do i = 1, n - 1
            call rgbmv(x((i - 0) * n + 1:(i + 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(:, 8:8), ms(4), n)
        end do

        do i = 1, n
            call rgbmv(x((i - 1) * n + 1:(i + 0) * n), y((i - 1) * n + 1:i * n), sign, blocks(:, 3:7), ms(3), n)
        end do

        do i = 2, n
            call rgbmv(x((i - 2) * n + 1:(i - 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(:, 2:2), ms(2), n)
        end do

        do i = 3, n
            call rgbmv(x((i - 3) * n + 1:(i - 2) * n), y((i - 1) * n + 1:i * n), sign, blocks(:, 1:1), ms(1), n)
        end do
    end subroutine rbbmv_o5

    subroutine rbbmv_o7(x, y, sign, blocks, ms, n)
        implicit none

        integer, parameter :: m = 7
        integer, intent(in) :: n
        integer, intent(in), dimension(m) :: ms
        real(sp), intent(in) :: sign
        real(sp), intent(in), dimension(n * n) :: x
        real(sp), intent(inout), dimension(n * n) :: y
        real(sp), intent(in), dimension(n, 2 * m - 1) :: blocks

        integer :: i

        do i = 4, n
            call rgbmv(x((i - 4) * n + 1:(i - 3) * n), y((i - 1) * n + 1:i * n), sign, blocks(:, 1:1), ms(1), n)
        end do

        do i = 3, n
            call rgbmv(x((i - 3) * n + 1:(i - 2) * n), y((i - 1) * n + 1:i * n), sign, blocks(:, 2:2), ms(2), n)
        end do

        do i = 2, n
            call rgbmv(x((i - 2) * n + 1:(i - 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(:, 3:3), ms(3), n)
        end do

        do i = 1, n
            call rgbmv(x((i - 1) * n + 1:(i + 0) * n), y((i - 1) * n + 1:i * n), sign, blocks(:, 4:10), ms(4), n)
        end do

        do i = 1, n - 1
            call rgbmv(x((i - 0) * n + 1:(i + 1) * n), y((i - 1) * n + 1:i * n), sign, blocks(:, 11:11), ms(5), n)
        end do

        do i = 1, n - 2
            call rgbmv(x((i + 1) * n + 1:(i + 2) * n), y((i - 1) * n + 1:i * n), sign, blocks(:, 12:12), ms(6), n)
        end do

        do i = 1, n - 3
            call rgbmv(x((i + 2) * n + 1:(i + 3) * n), y((i - 1) * n + 1:i * n), sign, blocks(:, 13:13), ms(7), n)
        end do
    end subroutine rbbmv_o7

    subroutine rbbmv(x, y, sign, blocks, ms, m, n)
        implicit none

        integer, intent(in) :: m, n
        integer, intent(in), dimension(m) :: ms
        real(sp), intent(in) :: sign
        real(sp), intent(in), dimension(n * n) :: x
        real(sp), intent(inout), dimension(n * n) :: y
        real(sp), intent(in), dimension(n, 2 * m - 1) :: blocks

        if (m == 3) then
            call rbbmv_o3(x, y, sign, blocks, ms, n)
        else if (m == 5) then
            call rbbmv_o5(x, y, sign, blocks, ms, n)
        else if (m == 7) then
            call rbbmv_o7(x, y, sign, blocks, ms, n)
        end if
    end subroutine rbbmv

    !   \brief Wrapper of BLAS-2 function `sgbmv` which is matvec implementation for banded matrix.
    subroutine rgbmv(x, u, sign, op, klu, n) ! reduced gbmv()
        implicit none

        integer, parameter :: sp = selected_real_kind(6, 37)
        integer, intent(in) :: n, klu
        real(sp), intent(in), dimension(n) :: x
        real(sp), intent(inout), dimension(n) :: u
        real(sp), intent(in), dimension(2 * klu + 1, n) :: op
        real(sp), intent(in) :: sign

        call sgbmv('N', n, n, klu, klu, sign, op, 2 * klu + 1, x, 1, 1.0, u, 1)
    end subroutine rgbmv

    !   \brief Calculate density of revervoir particles with given pumping profile and condensate density.
    !
    !   \param[in] pumping
    !   \verbatim
    !       Pumping profile on spacial grid.
    !   \endverbatim
    !
    !   \param[in] coeffs
    !   \verbatim
    !       Contains coefficients of NLS equation with reservoir. Each coefficient coresponds to each term of NLS
    !       equation or exciton-polariton revervoir. See `solve_nls`.
    !   \endverbatim
    !
    !   \param[in] u_sqr
    !   \verbatim
    !       Squared absolute value of psi-function.
    !   \endverbatim
    !
    !   \param[out] r
    !   \verbatim
    !       Density of revervoir particles.
    !   \endverbatim
    !
    !   \param[in] n
    !   \verbatim
    !       Number of nodes of spacial grid.
    !   \endverbatim
    subroutine revervoir(pumping, coeffs, u_sqr, r, n)
        implicit none

        integer, parameter :: sp = selected_real_kind(6, 37)
        integer, intent(in) :: n
        real(sp), intent(in), dimension(n) :: pumping
        real(sp), intent(in), dimension(23) :: coeffs
        real(sp), intent(in), dimension(n) :: u_sqr
        real(sp), intent(out), dimension(n) :: r ! actually n(r)

        r = coeffs(12) * pumping  / (coeffs(13) + coeffs(14) * u_sqr)
    end subroutine

    !   \brief Calculate action of hamiltonian operator on function defined on a grid with a given pumping and
    !   condensate density.
    !
    !   \param[in] pumping
    !   \verbatim
    !       Pumping profile on spacial grid.
    !   \endverbatim
    !
    !   \param[in] coeffs
    !   \verbatim
    !       Contains coefficients of NLS equation with reservoir. Each coefficient coresponds to each term of NLS
    !       equation or exciton-polariton revervoir. See `solve_nls`.
    !   \endverbatim
    !
    !   \param[in] u
    !   \verbatim
    !       Psi-function on which hamiltonian operator acts.
    !   \endverbatim
    !
    !   \param[out] v
    !   \verbatim
    !       Result of application of hamiltonian operator on the right hand side.
    !   \endverbatim
    !
    !   \param[out] op
    !   \verbatim
    !       Banded laplacian representation of order `klu` and size `n`. See `make_laplacian`.
    !   \endverbatim
    !
    !   \param[in] klu
    !   \verbatim
    !       Order of laplacian approximation.
    !   \endverbatim
    !
    !   \param[in] n
    !   \verbatim
    !       Number of nodes of spacial grid.
    !   \endverbatim
    subroutine hamiltonian(pumping, coeffs, u, v, op, klu, n)
        implicit none

        integer, parameter :: sp = selected_real_kind(6, 37)
        integer, intent(in) :: n, klu
        complex(sp), intent(in), dimension(n) :: u
        complex(sp), intent(out), dimension(n) :: v
        real(sp), intent(in), dimension(n) :: pumping
        real(sp), intent(in), dimension(23) :: coeffs
        real(sp), intent(in), dimension(2 * klu + 1, n) :: op

        complex(sp), parameter :: i = (0.0, 1.0)
        real(sp), parameter :: sign = 1.0
        real(sp), dimension(n) :: r, u_sqr
        real(sp), dimension(n) :: v_real, v_imag, u_real, u_imag

        u_real = real(u)
        u_imag = aimag(u)
        u_sqr = real(conjg(u) * u)

        call revervoir(pumping, coeffs, u_sqr, r, n)

        v_real = (coeffs(3) * r - coeffs(4)) * u_real + (coeffs(5) * u_sqr + coeffs(6) * r) * u_imag
        v_imag = (coeffs(3) * r - coeffs(4)) * u_imag - (coeffs(5) * u_sqr + coeffs(6) * r) * u_real

        call rgbmv(u_imag, v_real, -sign, op, klu, n)
        call rgbmv(u_real, v_imag, +sign, op, klu, n)

        v = cmplx(v_real, v_imag, sp)
    end subroutine

    !   \brief Yet another Runge-Kutta of order 4 implementation with matrix-valeud right hand side. The criterion of
    !   break is steady state solution which is not passed explicitly.
    !
    !   \param[in] dt
    !   \verbatim
    !       Time step.
    !   \endverbatim
    !
    !   \param[in] t0
    !   \verbatim
    !       Initial time.
    !   \endverbatim
    !
    !   \param[in] u0
    !   \verbatim
    !       Initial solution.
    !   \endverbatim
    !
    !   \param[in] op
    !   \verbatim
    !       Banded laplacian representation of order `order` and size `n`. See `make_laplacian`.
    !   \endverbatim
    !
    !   \param[in] n
    !   \verbatim
    !       Number of grid nodes.
    !   \endverbatim
    !
    !   \param[in] order
    !   \verbatim
    !       Order of laplacian approximation.
    !   \endverbatim
    !
    !   \param[in] iters
    !   \verbatim
    !       Number of iterations.
    !   \endverbatim
    !
    !   \param[out] u
    !   \verbatim
    !       Psi-function that solves NLS equation with reservoir.
    !   \endverbatim
    !
    !   \param[in] pumping
    !   \verbatim
    !       Pumping profile on spacial grid.
    !   \endverbatim
    !
    !   \param[in] coeffs
    !   \verbatim
    !       Contains coefficients of NLS equation with reservoir. Each coefficient coresponds to each term of NLS
    !       equation or exciton-polariton revervoir. See `solve_nls`.
    !   \endverbatim
    subroutine runge_kutta(dt, t0, u0, op, n, order, iters, u, pumping, coeffs)
        implicit none

        integer, parameter :: sp = selected_real_kind(6, 37)
        integer, intent(in) :: n, order, iters
        real(sp), intent(in) :: dt, t0
        real(sp), intent(in), dimension(order, n) :: op
        complex(sp), intent(in), dimension(n) :: u0
        complex(sp), intent(out), dimension(n) :: u
        real(sp), intent(in), dimension(n) :: pumping
        real(sp), intent(in), dimension(23) :: coeffs

        integer :: i, klu
        complex(sp) :: t
        complex(sp), dimension(n) :: k1, k2, k3, k4

        klu = (order - 1) / 2
        u = u0
        t = t0

        do i = 1, iters
            call hamiltonian(pumping, coeffs, u + 0. * dt / 2, k1, op, klu, n)
            call hamiltonian(pumping, coeffs, u + k1 * dt / 2, k2, op, klu, n)
            call hamiltonian(pumping, coeffs, u + k2 * dt / 2, k3, op, klu, n)
            call hamiltonian(pumping, coeffs, u + k3 * dt / 1, k4, op, klu, n)

            u = u + (k1 + 2 * k2 + 2 * k3 + k4) * dt / 6
            t = t + dt
        end do
    end subroutine runge_kutta

    !   \brief Solve Non-Linear Schrodinger equation in axial symmentric geometry. It is based on Runge-Kutta method of
    !   order 4 with laplacian approximation order that passed as argument. Initial solution is uniform everywhere on
    !   the grid.
    !
    !   \param[in] dt
    !   \verbatim
    !       Time step.
    !   \endverbatim
    !
    !   \param[in] dx
    !   \verbatim
    !       Spartial step.
    !   \endverbatim
    !
    !   \param[in] n
    !   \verbatim
    !       Number of grid nodes.
    !   \endverbatim
    !
    !   \param[in] order
    !   \verbatim
    !       Order of laplacian approximation.
    !   \endverbatim
    !
    !   \param[in] iters
    !   \verbatim
    !       Number of iterations.
    !   \endverbatim
    !
    !   \param[in] pumping
    !   \verbatim
    !       Pumping profile on spacial grid.
    !   \endverbatim
    !
    !   \param[in] coeffs
    !   \verbatim
    !       Contains coefficients of NLS equation with reservoir. Each coefficient coresponds to each term of NLS
    !       equation or exciton-polariton revervoir.
    !       Coefficient 1 corresponds to `i \partial_t` term.
    !       Coefficient 2 corresponds to `-\nabla^2` term.
    !       Coefficient 3 corresponds to `i n \psi` term.
    !       Coefficient 4 corresponds to `-i \psi` term.
    !       Coefficient 5 corresponds to `|\psi|^2 \psi` term.
    !       Coefficient 6 corresponds to `n \psi` term.
    !       Coefficient 11 corresponds to `\partial_t` term of reservoire equation.
    !       Coefficient 12 corresponds to `P term of reservoire equation.
    !       Coefficient 13 corresponds to `-n` term of reservoire equation.
    !       Coefficient 14 corresponds to `-n |\psi|^2` term of reservoire equation.
    !       Coefficient 15 corresponds to `\nabla^2` term of reservoire equation.
    !       The coresspondence is induced sequensially with equations in Wouters&Carusotto, 2007.
    !   \endverbatim
    !
    !   \param[in] u_0
    !   \verbatim
    !       Initial psi-function that solves NLS equation with reservoir.
    !   \endverbatim
    !
    !   \param[out] u
    !   \verbatim
    !       Psi-function that solves NLS equation with reservoir.
    !   \endverbatim
    subroutine solve_nls(dt, dx, n, order, iters, pumping, coeffs, u0, u)
        implicit none

        integer, parameter :: sp = selected_real_kind(6, 37)
        real(sp), intent(in) :: dt, dx
        integer, intent(in) :: n, order, iters
        real(sp), intent(in), dimension(n) :: pumping
        real(sp), intent(in), dimension(23) :: coeffs
        complex(sp), intent(in), dimension(n) :: u0
        complex(sp), intent(out), dimension(n) :: u

        real(sp), parameter :: t0 = 0.0
        real(sp), dimension(order, n) :: op

        call make_laplacian(n, order, dx, op)
        call runge_kutta(dt, t0, u0, op, n, order, iters, u, pumping, coeffs)
    end subroutine solve_nls

    subroutine solve_nls_1d(dt, dx, n, order, iters, pumping, coeffs, u0, u)
        implicit none

        integer, parameter :: sp = selected_real_kind(6, 37)
        real(sp), intent(in) :: dt, dx
        integer, intent(in) :: n, order, iters
        real(sp), intent(in), dimension(n) :: pumping
        complex(sp), intent(in), dimension(n) :: u0
        complex(sp), intent(out), dimension(n) :: u
        real(sp), intent(in), dimension(23) :: coeffs

        call solve_nls(dt, dx, n, order, iters, pumping, coeffs, u0, u)
    end subroutine solve_nls_1d

    subroutine revervoir_2d(pumping, coeffs, u_sqr, r, n)
        implicit none

        integer, intent(in) :: n
        real(sp), intent(in), dimension(n, n) :: pumping
        real(sp), intent(in), dimension(23) :: coeffs
        real(sp), intent(in), dimension(n, n) :: u_sqr
        real(sp), intent(out), dimension(n, n) :: r ! actually n(r)

        r = coeffs(12) * pumping  / (coeffs(13) + coeffs(14) * u_sqr)
    end subroutine revervoir_2d

    subroutine hamiltonian_2d(pumping, coeffs, u, v, blocks, orders, order, n)
        implicit none

        integer, intent(in) :: n, order
        integer, intent(in), dimension(order) :: orders
        complex(sp), intent(in), dimension(n, n) :: u
        complex(sp), intent(out), dimension(n, n) :: v
        real(sp), intent(in), dimension(n, n) :: pumping
        real(sp), intent(in), dimension(23) :: coeffs
        real(sp), intent(in), dimension(n, 2 * order - 1) :: blocks

        complex(sp), parameter :: i = (0.0, 1.0)
        real(sp), parameter :: sign = 1.0
        real(sp), dimension(n, n) :: r, u_sqr
        real(sp), dimension(n, n) :: v_real, v_imag, u_real, u_imag

        u_real = real(u)
        u_imag = aimag(u)
        u_sqr = real(conjg(u) * u)

        call revervoir_2d(pumping, coeffs, u_sqr, r, n)

        v_real = (coeffs(3) * r - coeffs(4)) * u_real + (coeffs(5) * u_sqr + coeffs(6) * r) * u_imag
        v_imag = (coeffs(3) * r - coeffs(4)) * u_imag - (coeffs(5) * u_sqr + coeffs(6) * r) * u_real

        call rbbmv(u_imag, v_real, -sign, blocks, orders, order, n)
        call rbbmv(u_real, v_imag, +sign, blocks, orders, order, n)

        v = cmplx(v_real, v_imag, sp)
    end subroutine hamiltonian_2d

    ! Code is almost the same as `runge_kutta`.
    subroutine runge_kutta_2d(dt, t0, u0, n, blocks, orders, order, iters, u, pumping, coeffs)
        implicit none

        integer, intent(in) :: n, order, iters
        integer, intent(in), dimension(order) :: orders
        real(sp), intent(in) :: dt, t0
        real(sp), intent(in), dimension(23) :: coeffs
        real(sp), intent(in), dimension(n, 2 * order - 1) :: blocks
        real(sp), intent(in), dimension(n, n) :: pumping
        complex(sp), intent(in), dimension(n, n) :: u0
        complex(sp), intent(out), dimension(n, n) :: u

        integer :: i
        real(sp) :: t
        complex(sp), dimension(n, n) :: k1, k2, k3, k4

        u = u0
        t = t0

        do i = 1, iters
            call hamiltonian_2d(pumping, coeffs, u + 0. * dt / 2, k1, blocks, orders, order, n)
            call hamiltonian_2d(pumping, coeffs, u + k1 * dt / 2, k2, blocks, orders, order, n)
            call hamiltonian_2d(pumping, coeffs, u + k2 * dt / 2, k3, blocks, orders, order, n)
            call hamiltonian_2d(pumping, coeffs, u + k3 * dt / 1, k4, blocks, orders, order, n)

            u = u + (k1 + 2 * k2 + 2 * k3 + k4) * dt / 6
            t = t + dt
        end do
    end subroutine runge_kutta_2d

    subroutine solve_nls_2d(dt, dx, n, order, iters, pumping, coeffs, u0, u)
        implicit none

        integer, intent(in) :: n, order, iters
        real(sp), intent(in) :: dt, dx
        real(sp), intent(in), dimension(n, n) :: pumping
        real(sp), intent(in), dimension(23) :: coeffs
        complex(sp), intent(in), dimension(n, n) :: u0
        complex(sp), intent(out), dimension(n, n) :: u

        integer, dimension(order) :: orders
        real(sp), parameter :: t0 = 0.0
        real(sp), dimension(n, 2 * order - 1) :: blocks

        call make_laplacian_2d(n, order, dx, blocks, orders)
        call runge_kutta_2d(dt, t0, u0, n, blocks, orders, order, iters, u, pumping, coeffs)
    end subroutine solve_nls_2d

    subroutine chemical_potential_1d(dx, n, pumping, coeffs, u0, mu)
        implicit none

        integer, intent(in) :: n
        real(sp), intent(in) :: dx
        real(sp), intent(in), dimension(n) :: pumping
        real(sp), intent(in), dimension(23) :: coeffs
        complex(sp), intent(out) :: mu  ! chemical potential
        complex(sp), intent(in), dimension(n) :: u0

        integer, parameter :: order = 5, klu = (order - 1) / 2
        integer :: i
        real(sp), dimension(order, n) :: op
        real(sp), dimension(n) :: r
        complex(sp), dimension(n) :: u
        complex(sp) :: M, E

        call make_laplacian(n, order, dx, op)
        call hamiltonian(pumping, coeffs, u0, u, op, klu, n)

        do i = 1, n
            r(i) = (i - 1.0) * dx
        end do

        M = dot_product(u0, u0 * r)
        E = (0.0, 1.0) * dot_product(u0, u * r)
        mu = E / M
    end subroutine chemical_potential_1d

    subroutine chemical_potential_2d(dx, n, pumping, coeffs, u0, mu)
        integer, intent(in) :: n
        real(sp), intent(in) :: dx
        real(sp), intent(in), dimension(n, n) :: pumping
        real(sp), intent(in), dimension(23) :: coeffs
        real(sp), intent(out) :: mu  ! chemical potential
        complex(sp), intent(in), dimension(n, n) :: u0

        integer, parameter :: order = 5, klu = (order - 1) / 2
        integer, dimension(order) :: orders
        real(sp), dimension(2 * order - 1, n) :: blocks
        complex(sp), dimension(n, n) :: u
        complex(sp) :: M, E

        call make_laplacian_2d(n, order, dx, blocks, orders)
        call hamiltonian_2d(pumping, coeffs, u0, u, blocks, orders, order, n)

        M = dot_product(reshape(u0, (/ n * n /)), reshape(u0, (/ n * n /)))
        E = (0.0, 1.0) * dot_product(reshape(u0, (/ n * n /)), reshape(u, (/ n * n /)))

        mu = E / M
    end subroutine chemical_potential_2d

end module nls

max_stars_repo_name: daskol/nls

max_issues_repo_name: daskol/nls

max_forks_repo_name: daskol/nls


text:
      double complex function pvfndd(n,x,iep)
C----Implementation of DD Eq. 4.11
      implicit none
      include 'TRonshellcutoff.f'
      integer j,n,infty
      double complex xm1,x,cln,cone
      double precision iep
      parameter(cone=(1d0,0d0),infty=16) ! number of terms in sum

      xm1=x-cone
      if (abs(x) .lt. 10d0) then
        if (abs(x-cone) .lt. onshellcutoff) then
          pvfndd=0d0
        else
          pvfndd=(cone-x**(n+1))*(cln(xm1,iep)-cln(x,iep))
        endif
        do j=0,n
          pvfndd=pvfndd-x**(n-j)/dfloat(j+1)
        enddo
      else
        pvfndd=cln(cone-cone/x,iep)
        do j=n+1,n+infty
          pvfndd=pvfndd+x**(n-j)/dfloat(j+1)
        enddo
      endif

      return
      end

ext: f

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
subroutine normal(pf)
!from physical to normalized ([0-1]) units

use share, only: dp, ndim, n_p, llimits, steps, ntie, indtie

implicit none

!params
real(dp), intent(inout) :: pf(ndim) ! full vector of parameters

!locals
real(dp)                :: ulimit
integer                 :: i

do i=1,ndim
  ulimit=llimits(i)+steps(i)*(n_p(i)-1)
  pf(i)=0.5_dp+(2._dp*pf(i)-llimits(i)-ulimit)/(2._dp*(ulimit-llimits(i)))
enddo

if (indtie(1) > 0) then
  do i=1,ntie
    if (pf(indtie(i)) < 0.0_dp) then
      pf(indtie(i))=0.0_dp
    else
        if (pf(indtie(i)) >= 1.0_dp) pf(indtie(i))=0.999999_dp
    endif
  enddo
endif

end

max_stars_repo_name: dnidever/apogee

max_issues_repo_name: dnidever/apogee

max_forks_repo_name: dnidever/apogee


text:
*DECK HW3CRT
      SUBROUTINE HW3CRT (XS, XF, L, LBDCND, BDXS, BDXF, YS, YF, M,
     +   MBDCND, BDYS, BDYF, ZS, ZF, N, NBDCND, BDZS, BDZF, ELMBDA,
     +   LDIMF, MDIMF, F, PERTRB, IERROR, W)
C***BEGIN PROLOGUE  HW3CRT
C***PURPOSE  Solve the standard seven-point finite difference
C            approximation to the Helmholtz equation in Cartesian
C            coordinates.
C***LIBRARY   SLATEC (FISHPACK)
C***CATEGORY  I2B1A1A
C***TYPE      SINGLE PRECISION (HW3CRT-S)
C***KEYWORDS  CARTESIAN, ELLIPTIC, FISHPACK, HELMHOLTZ, PDE
C***AUTHOR  Adams, J., (NCAR)
C           Swarztrauber, P. N., (NCAR)
C           Sweet, R., (NCAR)
C***DESCRIPTION
C
C     Subroutine HW3CRT solves the standard seven-point finite
C     difference approximation to the Helmholtz equation in Cartesian
C     coordinates:
C
C         (d/dX)(dU/dX) + (d/dY)(dU/dY) + (d/dZ)(dU/dZ)
C
C                    + LAMBDA*U = F(X,Y,Z) .
C
C    * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C
C    * * * * * * * *    Parameter Description     * * * * * * * * * *
C
C
C            * * * * * *   On Input    * * * * * *
C
C     XS,XF
C        The range of X, i.e. XS .LE. X .LE. XF .
C        XS must be less than XF.
C
C     L
C        The number of panels into which the interval (XS,XF) is
C        subdivided.  Hence, there will be L+1 grid points in the
C        X-direction given by X(I) = XS+(I-1)DX for I=1,2,...,L+1,
C        where DX = (XF-XS)/L is the panel width.  L must be at
C        least 5 .
C
C     LBDCND
C        Indicates the type of boundary conditions at X = XS and X = XF.
C
C        = 0  If the solution is periodic in X, i.e.
C             U(L+I,J,K) = U(I,J,K).
C        = 1  If the solution is specified at X = XS and X = XF.
C        = 2  If the solution is specified at X = XS and the derivative
C             of the solution with respect to X is specified at X = XF.
C        = 3  If the derivative of the solution with respect to X is
C             specified at X = XS and X = XF.
C        = 4  If the derivative of the solution with respect to X is
C             specified at X = XS and the solution is specified at X=XF.
C
C     BDXS
C        A two-dimensional array that specifies the values of the
C        derivative of the solution with respect to X at X = XS.
C        when LBDCND = 3 or 4,
C
C             BDXS(J,K) = (d/dX)U(XS,Y(J),Z(K)), J=1,2,...,M+1,
C                                                K=1,2,...,N+1.
C
C        When LBDCND has any other value, BDXS is a dummy variable.
C        BDXS must be dimensioned at least (M+1)*(N+1).
C
C     BDXF
C        A two-dimensional array that specifies the values of the
C        derivative of the solution with respect to X at X = XF.
C        When LBDCND = 2 or 3,
C
C             BDXF(J,K) = (d/dX)U(XF,Y(J),Z(K)), J=1,2,...,M+1,
C                                                K=1,2,...,N+1.
C
C        When LBDCND has any other value, BDXF is a dummy variable.
C        BDXF must be dimensioned at least (M+1)*(N+1).
C
C     YS,YF
C        The range of Y, i.e. YS .LE. Y .LE. YF.
C        YS must be less than YF.
C
C     M
C        The number of panels into which the interval (YS,YF) is
C        subdivided.  Hence, there will be M+1 grid points in the
C        Y-direction given by Y(J) = YS+(J-1)DY for J=1,2,...,M+1,
C        where DY = (YF-YS)/M is the panel width.  M must be at
C        least 5 .
C
C     MBDCND
C        Indicates the type of boundary conditions at Y = YS and Y = YF.
C
C        = 0  If the solution is periodic in Y, i.e.
C             U(I,M+J,K) = U(I,J,K).
C        = 1  If the solution is specified at Y = YS and Y = YF.
C        = 2  If the solution is specified at Y = YS and the derivative
C             of the solution with respect to Y is specified at Y = YF.
C        = 3  If the derivative of the solution with respect to Y is
C             specified at Y = YS and Y = YF.
C        = 4  If the derivative of the solution with respect to Y is
C             specified at Y = YS and the solution is specified at Y=YF.
C
C     BDYS
C        A two-dimensional array that specifies the values of the
C        derivative of the solution with respect to Y at Y = YS.
C        When MBDCND = 3 or 4,
C
C             BDYS(I,K) = (d/dY)U(X(I),YS,Z(K)), I=1,2,...,L+1,
C                                                K=1,2,...,N+1.
C
C        When MBDCND has any other value, BDYS is a dummy variable.
C        BDYS must be dimensioned at least (L+1)*(N+1).
C
C     BDYF
C        A two-dimensional array that specifies the values of the
C        derivative of the solution with respect to Y at Y = YF.
C        When MBDCND = 2 or 3,
C
C             BDYF(I,K) = (d/dY)U(X(I),YF,Z(K)), I=1,2,...,L+1,
C                                                K=1,2,...,N+1.
C
C        When MBDCND has any other value, BDYF is a dummy variable.
C        BDYF must be dimensioned at least (L+1)*(N+1).
C
C     ZS,ZF
C        The range of Z, i.e. ZS .LE. Z .LE. ZF.
C        ZS must be less than ZF.
C
C     N
C        The number of panels into which the interval (ZS,ZF) is
C        subdivided.  Hence, there will be N+1 grid points in the
C        Z-direction given by Z(K) = ZS+(K-1)DZ for K=1,2,...,N+1,
C        where DZ = (ZF-ZS)/N is the panel width.  N must be at least 5.
C
C     NBDCND
C        Indicates the type of boundary conditions at Z = ZS and Z = ZF.
C
C        = 0  If the solution is periodic in Z, i.e.
C             U(I,J,N+K) = U(I,J,K).
C        = 1  If the solution is specified at Z = ZS and Z = ZF.
C        = 2  If the solution is specified at Z = ZS and the derivative
C             of the solution with respect to Z is specified at Z = ZF.
C        = 3  If the derivative of the solution with respect to Z is
C             specified at Z = ZS and Z = ZF.
C        = 4  If the derivative of the solution with respect to Z is
C             specified at Z = ZS and the solution is specified at Z=ZF.
C
C     BDZS
C        A two-dimensional array that specifies the values of the
C        derivative of the solution with respect to Z at Z = ZS.
C        When NBDCND = 3 or 4,
C
C             BDZS(I,J) = (d/dZ)U(X(I),Y(J),ZS), I=1,2,...,L+1,
C                                                J=1,2,...,M+1.
C
C        When NBDCND has any other value, BDZS is a dummy variable.
C        BDZS must be dimensioned at least (L+1)*(M+1).
C
C     BDZF
C        A two-dimensional array that specifies the values of the
C        derivative of the solution with respect to Z at Z = ZF.
C        When NBDCND = 2 or 3,
C
C             BDZF(I,J) = (d/dZ)U(X(I),Y(J),ZF), I=1,2,...,L+1,
C                                                J=1,2,...,M+1.
C
C        When NBDCND has any other value, BDZF is a dummy variable.
C        BDZF must be dimensioned at least (L+1)*(M+1).
C
C     ELMBDA
C        The constant LAMBDA in the Helmholtz equation. If
C        LAMBDA .GT. 0, a solution may not exist.  However, HW3CRT will
C        attempt to find a solution.
C
C     F
C        A three-dimensional array that specifies the values of the
C        right side of the Helmholtz equation and boundary values (if
C        any).  For I=2,3,...,L, J=2,3,...,M, and K=2,3,...,N
C
C                   F(I,J,K) = F(X(I),Y(J),Z(K)).
C
C        On the boundaries F is defined by
C
C        LBDCND      F(1,J,K)         F(L+1,J,K)
C        ------   ---------------   ---------------
C
C          0      F(XS,Y(J),Z(K))   F(XS,Y(J),Z(K))
C          1      U(XS,Y(J),Z(K))   U(XF,Y(J),Z(K))
C          2      U(XS,Y(J),Z(K))   F(XF,Y(J),Z(K))   J=1,2,...,M+1
C          3      F(XS,Y(J),Z(K))   F(XF,Y(J),Z(K))   K=1,2,...,N+1
C          4      F(XS,Y(J),Z(K))   U(XF,Y(J),Z(K))
C
C        MBDCND      F(I,1,K)         F(I,M+1,K)
C        ------   ---------------   ---------------
C
C          0      F(X(I),YS,Z(K))   F(X(I),YS,Z(K))
C          1      U(X(I),YS,Z(K))   U(X(I),YF,Z(K))
C          2      U(X(I),YS,Z(K))   F(X(I),YF,Z(K))   I=1,2,...,L+1
C          3      F(X(I),YS,Z(K))   F(X(I),YF,Z(K))   K=1,2,...,N+1
C          4      F(X(I),YS,Z(K))   U(X(I),YF,Z(K))
C
C        NBDCND      F(I,J,1)         F(I,J,N+1)
C        ------   ---------------   ---------------
C
C          0      F(X(I),Y(J),ZS)   F(X(I),Y(J),ZS)
C          1      U(X(I),Y(J),ZS)   U(X(I),Y(J),ZF)
C          2      U(X(I),Y(J),ZS)   F(X(I),Y(J),ZF)   I=1,2,...,L+1
C          3      F(X(I),Y(J),ZS)   F(X(I),Y(J),ZF)   J=1,2,...,M+1
C          4      F(X(I),Y(J),ZS)   U(X(I),Y(J),ZF)
C
C        F must be dimensioned at least (L+1)*(M+1)*(N+1).
C
C        NOTE:
C
C        If the table calls for both the solution U and the right side F
C        on a boundary, then the solution must be specified.
C
C     LDIMF
C        The row (or first) dimension of the arrays F,BDYS,BDYF,BDZS,
C        and BDZF as it appears in the program calling HW3CRT. this
C        parameter is used to specify the variable dimension of these
C        arrays.  LDIMF must be at least L+1.
C
C     MDIMF
C        The column (or second) dimension of the array F and the row (or
C        first) dimension of the arrays BDXS and BDXF as it appears in
C        the program calling HW3CRT.  This parameter is used to specify
C        the variable dimension of these arrays.
C        MDIMF must be at least M+1.
C
C     W
C        A one-dimensional array that must be provided by the user for
C        work space.  The length of W must be at least 30 + L + M + 5*N
C        + MAX(L,M,N) + 7*(INT((L+1)/2) + INT((M+1)/2))
C
C
C            * * * * * *   On Output   * * * * * *
C
C     F
C        Contains the solution U(I,J,K) of the finite difference
C        approximation for the grid point (X(I),Y(J),Z(K)) for
C        I=1,2,...,L+1, J=1,2,...,M+1, and K=1,2,...,N+1.
C
C     PERTRB
C        If a combination of periodic or derivative boundary conditions
C        is specified for a Poisson equation (LAMBDA = 0), a solution
C        may not exist.  PERTRB is a constant, calculated and subtracted
C        from F, which ensures that a solution exists.  PWSCRT then
C        computes this solution, which is a least squares solution to
C        the original approximation.  This solution is not unique and is
C        unnormalized.  The value of PERTRB should be small compared to
C        the right side F.  Otherwise, a solution is obtained to an
C        essentially different problem.  This comparison should always
C        be made to insure that a meaningful solution has been obtained.
C
C     IERROR
C        An error flag that indicates invalid input parameters.  Except
C        for numbers 0 and 12, a solution is not attempted.
C
C        =  0  No error
C        =  1  XS .GE. XF
C        =  2  L .LT. 5
C        =  3  LBDCND .LT. 0 .OR. LBDCND .GT. 4
C        =  4  YS .GE. YF
C        =  5  M .LT. 5
C        =  6  MBDCND .LT. 0 .OR. MBDCND .GT. 4
C        =  7  ZS .GE. ZF
C        =  8  N .LT. 5
C        =  9  NBDCND .LT. 0 .OR. NBDCND .GT. 4
C        = 10  LDIMF .LT. L+1
C        = 11  MDIMF .LT. M+1
C        = 12  LAMBDA .GT. 0
C
C        Since this is the only means of indicating a possibly incorrect
C        call to HW3CRT, the user should test IERROR after the call.
C
C *Long Description:
C
C    * * * * * * *   Program Specifications    * * * * * * * * * * * *
C
C     Dimension of   BDXS(MDIMF,N+1),BDXF(MDIMF,N+1),BDYS(LDIMF,N+1),
C     Arguments      BDYF(LDIMF,N+1),BDZS(LDIMF,M+1),BDZF(LDIMF,M+1),
C                    F(LDIMF,MDIMF,N+1),W(see argument list)
C
C     Latest         December 1, 1978
C     Revision
C
C     Subprograms    HW3CRT,POIS3D,POS3D1,TRIDQ,RFFTI,RFFTF,RFFTF1,
C     Required       RFFTB,RFFTB1,COSTI,COST,SINTI,SINT,COSQI,COSQF,
C                    COSQF1,COSQB,COSQB1,SINQI,SINQF,SINQB,CFFTI,
C                    CFFTI1,CFFTB,CFFTB1,PASSB2,PASSB3,PASSB4,PASSB,
C                    CFFTF,CFFTF1,PASSF1,PASSF2,PASSF3,PASSF4,PASSF,
C                    PIMACH
C
C     Special        NONE
C     Conditions
C
C     Common         NONE
C     Blocks
C
C     I/O            NONE
C
C     Precision      Single
C
C     Specialist     Roland Sweet
C
C     Language       FORTRAN
C
C     History        Written by Roland Sweet at NCAR in July 1977
C
C     Algorithm      This subroutine defines the finite difference
C                    equations, incorporates boundary data, and
C                    adjusts the right side of singular systems and
C                    then calls POIS3D to solve the system.
C
C     Space          7862(decimal) = 17300(octal) locations on the
C     Required       NCAR Control Data 7600
C
C     Timing and        The execution time T on the NCAR Control Data
C     Accuracy       7600 for subroutine HW3CRT is roughly proportional
C                    to L*M*N*(log2(L)+log2(M)+5), but also depends on
C                    input parameters LBDCND and MBDCND.  Some typical
C                    values are listed in the table below.
C                       The solution process employed results in a loss
C                    of no more than three significant digits for L,M
C                    and N as large as 32.  More detailed information
C                    about accuracy can be found in the documentation
C                    for subroutine POIS3D which is the routine that
C                    actually solves the finite difference equations.
C
C
C                       L(=M=N)     LBDCND(=MBDCND=NBDCND)      T(MSECS)
C                       -------     ----------------------      --------
C
C                         16                  0                    300
C                         16                  1                    302
C                         16                  3                    348
C                         32                  0                   1925
C                         32                  1                   1929
C                         32                  3                   2109
C
C     Portability    American National Standards Institute FORTRAN.
C                    The machine dependent constant PI is defined in
C                    function PIMACH.
C
C     Required       COS,SIN,ATAN
C     Resident
C     Routines
C
C     Reference      NONE
C
C    * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  POIS3D
C***REVISION HISTORY  (YYMMDD)
C   801001  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C***END PROLOGUE  HW3CRT
C
C
      DIMENSION       BDXS(MDIMF,*)          ,BDXF(MDIMF,*)          ,
     1                BDYS(LDIMF,*)          ,BDYF(LDIMF,*)          ,
     2                BDZS(LDIMF,*)          ,BDZF(LDIMF,*)          ,
     3                F(LDIMF,MDIMF,*)       ,W(*)
C***FIRST EXECUTABLE STATEMENT  HW3CRT
      IERROR = 0
      IF (XF .LE. XS) IERROR = 1
      IF (L .LT. 5) IERROR = 2
      IF (LBDCND.LT.0 .OR. LBDCND.GT.4) IERROR = 3
      IF (YF .LE. YS) IERROR = 4
      IF (M .LT. 5) IERROR = 5
      IF (MBDCND.LT.0 .OR. MBDCND.GT.4) IERROR = 6
      IF (ZF .LE. ZS) IERROR = 7
      IF (N .LT. 5) IERROR = 8
      IF (NBDCND.LT.0 .OR. NBDCND.GT.4) IERROR = 9
      IF (LDIMF .LT. L+1) IERROR = 10
      IF (MDIMF .LT. M+1) IERROR = 11
      IF (IERROR .NE. 0) GO TO 188
      DY = (YF-YS)/M
      TWBYDY = 2./DY
      C2 = 1./(DY**2)
      MSTART = 1
      MSTOP = M
      MP1 = M+1
      MP = MBDCND+1
      GO TO (104,101,101,102,102),MP
  101 MSTART = 2
  102 GO TO (104,104,103,103,104),MP
  103 MSTOP = MP1
  104 MUNK = MSTOP-MSTART+1
      DZ = (ZF-ZS)/N
      TWBYDZ = 2./DZ
      NP = NBDCND+1
      C3 = 1./(DZ**2)
      NP1 = N+1
      NSTART = 1
      NSTOP = N
      GO TO (108,105,105,106,106),NP
  105 NSTART = 2
  106 GO TO (108,108,107,107,108),NP
  107 NSTOP = NP1
  108 NUNK = NSTOP-NSTART+1
      LP1 = L+1
      DX = (XF-XS)/L
      C1 = 1./(DX**2)
      TWBYDX = 2./DX
      LP = LBDCND+1
      LSTART = 1
      LSTOP = L
C
C     ENTER BOUNDARY DATA FOR X-BOUNDARIES.
C
      GO TO (122,109,109,112,112),LP
  109 LSTART = 2
      DO 111 J=MSTART,MSTOP
         DO 110 K=NSTART,NSTOP
            F(2,J,K) = F(2,J,K)-C1*F(1,J,K)
  110    CONTINUE
  111 CONTINUE
      GO TO 115
  112 DO 114 J=MSTART,MSTOP
         DO 113 K=NSTART,NSTOP
            F(1,J,K) = F(1,J,K)+TWBYDX*BDXS(J,K)
  113    CONTINUE
  114 CONTINUE
  115 GO TO (122,116,119,119,116),LP
  116 DO 118 J=MSTART,MSTOP
         DO 117 K=NSTART,NSTOP
            F(L,J,K) = F(L,J,K)-C1*F(LP1,J,K)
  117    CONTINUE
  118 CONTINUE
      GO TO 122
  119 LSTOP = LP1
      DO 121 J=MSTART,MSTOP
         DO 120 K=NSTART,NSTOP
            F(LP1,J,K) = F(LP1,J,K)-TWBYDX*BDXF(J,K)
  120    CONTINUE
  121 CONTINUE
  122 LUNK = LSTOP-LSTART+1
C
C     ENTER BOUNDARY DATA FOR Y-BOUNDARIES.
C
      GO TO (136,123,123,126,126),MP
  123 DO 125 I=LSTART,LSTOP
         DO 124 K=NSTART,NSTOP
            F(I,2,K) = F(I,2,K)-C2*F(I,1,K)
  124    CONTINUE
  125 CONTINUE
      GO TO 129
  126 DO 128 I=LSTART,LSTOP
         DO 127 K=NSTART,NSTOP
            F(I,1,K) = F(I,1,K)+TWBYDY*BDYS(I,K)
  127    CONTINUE
  128 CONTINUE
  129 GO TO (136,130,133,133,130),MP
  130 DO 132 I=LSTART,LSTOP
         DO 131 K=NSTART,NSTOP
            F(I,M,K) = F(I,M,K)-C2*F(I,MP1,K)
  131    CONTINUE
  132 CONTINUE
      GO TO 136
  133 DO 135 I=LSTART,LSTOP
         DO 134 K=NSTART,NSTOP
            F(I,MP1,K) = F(I,MP1,K)-TWBYDY*BDYF(I,K)
  134    CONTINUE
  135 CONTINUE
  136 CONTINUE
C
C     ENTER BOUNDARY DATA FOR Z-BOUNDARIES.
C
      GO TO (150,137,137,140,140),NP
  137 DO 139 I=LSTART,LSTOP
         DO 138 J=MSTART,MSTOP
            F(I,J,2) = F(I,J,2)-C3*F(I,J,1)
  138    CONTINUE
  139 CONTINUE
      GO TO 143
  140 DO 142 I=LSTART,LSTOP
         DO 141 J=MSTART,MSTOP
            F(I,J,1) = F(I,J,1)+TWBYDZ*BDZS(I,J)
  141    CONTINUE
  142 CONTINUE
  143 GO TO (150,144,147,147,144),NP
  144 DO 146 I=LSTART,LSTOP
         DO 145 J=MSTART,MSTOP
            F(I,J,N) = F(I,J,N)-C3*F(I,J,NP1)
  145    CONTINUE
  146 CONTINUE
      GO TO 150
  147 DO 149 I=LSTART,LSTOP
         DO 148 J=MSTART,MSTOP
            F(I,J,NP1) = F(I,J,NP1)-TWBYDZ*BDZF(I,J)
  148    CONTINUE
  149 CONTINUE
C
C     DEFINE A,B,C COEFFICIENTS IN W-ARRAY.
C
  150 CONTINUE
      IWB = NUNK+1
      IWC = IWB+NUNK
      IWW = IWC+NUNK
      DO 151 K=1,NUNK
         I = IWC+K-1
         W(K) = C3
         W(I) = C3
         I = IWB+K-1
         W(I) = -2.*C3+ELMBDA
  151 CONTINUE
      GO TO (155,155,153,152,152),NP
  152 W(IWC) = 2.*C3
  153 GO TO (155,155,154,154,155),NP
  154 W(IWB-1) = 2.*C3
  155 CONTINUE
      PERTRB = 0.
C
C     FOR SINGULAR PROBLEMS ADJUST DATA TO INSURE A SOLUTION WILL EXIST.
C
      GO TO (156,172,172,156,172),LP
  156 GO TO (157,172,172,157,172),MP
  157 GO TO (158,172,172,158,172),NP
  158 IF (ELMBDA) 172,160,159
  159 IERROR = 12
      GO TO 172
  160 CONTINUE
      MSTPM1 = MSTOP-1
      LSTPM1 = LSTOP-1
      NSTPM1 = NSTOP-1
      XLP = (2+LP)/3
      YLP = (2+MP)/3
      ZLP = (2+NP)/3
      S1 = 0.
      DO 164 K=2,NSTPM1
         DO 162 J=2,MSTPM1
            DO 161 I=2,LSTPM1
               S1 = S1+F(I,J,K)
  161       CONTINUE
            S1 = S1+(F(1,J,K)+F(LSTOP,J,K))/XLP
  162    CONTINUE
         S2 = 0.
         DO 163 I=2,LSTPM1
            S2 = S2+F(I,1,K)+F(I,MSTOP,K)
  163    CONTINUE
         S2 = (S2+(F(1,1,K)+F(1,MSTOP,K)+F(LSTOP,1,K)+F(LSTOP,MSTOP,K))/
     1                                                          XLP)/YLP
         S1 = S1+S2
  164 CONTINUE
      S = (F(1,1,1)+F(LSTOP,1,1)+F(1,1,NSTOP)+F(LSTOP,1,NSTOP)+
     1    F(1,MSTOP,1)+F(LSTOP,MSTOP,1)+F(1,MSTOP,NSTOP)+
     2                                   F(LSTOP,MSTOP,NSTOP))/(XLP*YLP)
      DO 166 J=2,MSTPM1
         DO 165 I=2,LSTPM1
            S = S+F(I,J,1)+F(I,J,NSTOP)
  165    CONTINUE
  166 CONTINUE
      S2 = 0.
      DO 167 I=2,LSTPM1
         S2 = S2+F(I,1,1)+F(I,1,NSTOP)+F(I,MSTOP,1)+F(I,MSTOP,NSTOP)
  167 CONTINUE
      S = S2/YLP+S
      S2 = 0.
      DO 168 J=2,MSTPM1
         S2 = S2+F(1,J,1)+F(1,J,NSTOP)+F(LSTOP,J,1)+F(LSTOP,J,NSTOP)
  168 CONTINUE
      S = S2/XLP+S
      PERTRB = (S/ZLP+S1)/((LUNK+1.-XLP)*(MUNK+1.-YLP)*
     1                                              (NUNK+1.-ZLP))
      DO 171 I=1,LUNK
         DO 170 J=1,MUNK
            DO 169 K=1,NUNK
               F(I,J,K) = F(I,J,K)-PERTRB
  169       CONTINUE
  170    CONTINUE
  171 CONTINUE
  172 CONTINUE
      NPEROD = 0
      IF (NBDCND .EQ. 0) GO TO 173
      NPEROD = 1
      W(1) = 0.
      W(IWW-1) = 0.
  173 CONTINUE
      CALL POIS3D (LBDCND,LUNK,C1,MBDCND,MUNK,C2,NPEROD,NUNK,W,W(IWB),
     1             W(IWC),LDIMF,MDIMF,F(LSTART,MSTART,NSTART),IR,W(IWW))
C
C     FILL IN SIDES FOR PERIODIC BOUNDARY CONDITIONS.
C
      IF (LP .NE. 1) GO TO 180
      IF (MP .NE. 1) GO TO 175
      DO 174 K=NSTART,NSTOP
         F(1,MP1,K) = F(1,1,K)
  174 CONTINUE
      MSTOP = MP1
  175 IF (NP .NE. 1) GO TO 177
      DO 176 J=MSTART,MSTOP
         F(1,J,NP1) = F(1,J,1)
  176 CONTINUE
      NSTOP = NP1
  177 DO 179 J=MSTART,MSTOP
         DO 178 K=NSTART,NSTOP
            F(LP1,J,K) = F(1,J,K)
  178    CONTINUE
  179 CONTINUE
  180 CONTINUE
      IF (MP .NE. 1) GO TO 185
      IF (NP .NE. 1) GO TO 182
      DO 181 I=LSTART,LSTOP
         F(I,1,NP1) = F(I,1,1)
  181 CONTINUE
      NSTOP = NP1
  182 DO 184 I=LSTART,LSTOP
         DO 183 K=NSTART,NSTOP
            F(I,MP1,K) = F(I,1,K)
  183    CONTINUE
  184 CONTINUE
  185 CONTINUE
      IF (NP .NE. 1) GO TO 188
      DO 187 I=LSTART,LSTOP
         DO 186 J=MSTART,MSTOP
            F(I,J,NP1) = F(I,J,1)
  186    CONTINUE
  187 CONTINUE
  188 CONTINUE
      RETURN
      END


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

      subroutine dlwq60 ( deriv  , conc   , notot  , noseg  , itfact ,
     *                    amass2 , isys   , nsys   , dmps   , intopt ,
     *                    isdmp  )

!     Deltares Software Centre

!>\File
!>           scales DERIV after the user quality processes, for steady state computation.

!     CREATED: april 3, 1988 by L.Postma

!     LOGICAL UNITNUMBERS : none

!     SUBROUTINES CALLED  : none

      use timers

      implicit none

!     Parameters          :

!     type     kind  function         name                      description
      integer   (4), intent(in   ) :: notot                   !< total number of substances
      integer   (4), intent(in   ) :: noseg                   !< number of computational volumes
      real      (4), intent(inout) :: deriv (notot ,noseg)    !< derivatives to be scaled
      real      (4), intent(inout) :: conc  (notot ,noseg)    !< concentrations per substance per volume
      integer   (4), intent(in   ) :: itfact                  !< scale factor between clocks
      real      (4), intent(inout) :: amass2(notot , 5   )    !< mass balance array
      integer   (4), intent(in   ) :: isys                    !< 'this' substance
      integer   (4), intent(in   ) :: nsys                    !< number of substances
      real      (4), intent(inout) :: dmps  (notot,*)         !< dumped fluxes is intopt > 7
      integer   (4), intent(in   ) :: intopt                  !< Integration suboptions
      integer   (4), intent(in   ) :: isdmp (noseg)           !< Pointer dumped segments

!     Local declarations

      integer(4) iseg, i, ip   ! Loop and help variables

      integer(4) ithandl /0/
      if ( timon ) call timstrt ( "dlwq60", ithandl )

!         loop accross deriv and conc

      do iseg = 1    , noseg
         conc  (isys,iseg) = conc  (isys,iseg)/itfact
         do i    = isys , isys +nsys-1
            deriv (i,iseg) = deriv (i,iseg)/itfact
            amass2(i,   2) = amass2(i,   2) + deriv(i,iseg)
         enddo
      enddo

      if ( mod(intopt,16) .ge. 8  ) then
         do iseg = 1 , noseg
            ip = isdmp(iseg)
            if ( ip .gt. 0 ) then
               do i    = isys , isys +nsys-1
                  dmps(i,ip)= dmps(i,ip) + deriv(i,iseg)
               enddo
            endif
         enddo
      endif

      if ( timon ) call timstop ( ithandl )

      return
      end

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


!
! Distributed under the OSI-approved Apache License, Version 2.0.  See
!  accompanying file Copyright.txt for details.
!
!  adios2_variable_min_mod.f90 : ADIOS2 Fortran bindings for overloaded
!                                adios2_variable_min subroutine
!   Created on: Nov 15, 2018
!       Author: William F Godoy godoywf@ornl.gov
!

module adios2_variable_min_mod
    use adios2_parameters_mod
    use adios2_variable_mod
    implicit none

    interface adios2_variable_min

        module procedure adios2_variable_min_real
        module procedure adios2_variable_min_dp
        module procedure adios2_variable_min_complex
        module procedure adios2_variable_min_complex_dp
        module procedure adios2_variable_min_integer1
        module procedure adios2_variable_min_integer2
        module procedure adios2_variable_min_integer4
        module procedure adios2_variable_min_integer8

    end interface
    external adios2_variable_min_f2c

contains

    subroutine adios2_variable_min_real(minimum, variable, ierr)
        real, intent(out) :: minimum
        type(adios2_variable), intent(in) :: variable
        integer, intent(out) :: ierr

        call adios2_variable_check_type(variable, adios2_type_real, &
                                        'variable_min', ierr)
        if (ierr == 0) then
            call adios2_variable_min_f2c(minimum, variable%f2c, ierr)
        end if

    end subroutine

    subroutine adios2_variable_min_dp(minimum, variable, ierr)
        real(kind=8), intent(out) :: minimum
        type(adios2_variable), intent(in) :: variable
        integer, intent(out) :: ierr

        call adios2_variable_check_type(variable, adios2_type_dp, &
                                        'variable_min', ierr)
        if (ierr == 0) then
            call adios2_variable_min_f2c(minimum, variable%f2c, ierr)
        end if

    end subroutine

    subroutine adios2_variable_min_complex(minimum, variable, ierr)
        complex, intent(out) :: minimum
        type(adios2_variable), intent(in) :: variable
        integer, intent(out) :: ierr

        call adios2_variable_check_type(variable, adios2_type_complex, &
                                        'variable_min', ierr)
        if (ierr == 0) then
            call adios2_variable_min_f2c(minimum, variable%f2c, ierr)
        end if

    end subroutine

    subroutine adios2_variable_min_complex_dp(minimum, variable, ierr)
        complex(kind=8), intent(out) :: minimum
        type(adios2_variable), intent(in) :: variable
        integer, intent(out) :: ierr

        call adios2_variable_check_type(variable, adios2_type_complex_dp, &
                                        'variable_min', ierr)
        if (ierr == 0) then
            call adios2_variable_min_f2c(minimum, variable%f2c, ierr)
        end if

    end subroutine

    subroutine adios2_variable_min_integer1(minimum, variable, ierr)
        integer(kind=1), intent(out) :: minimum
        type(adios2_variable), intent(in) :: variable
        integer, intent(out) :: ierr

        call adios2_variable_check_type(variable, adios2_type_integer1, &
                                        'variable_min', ierr)
        if (ierr == 0) then
            call adios2_variable_min_f2c(minimum, variable%f2c, ierr)
        end if

    end subroutine

    subroutine adios2_variable_min_integer2(minimum, variable, ierr)
        integer(kind=2), intent(out) :: minimum
        type(adios2_variable), intent(in) :: variable
        integer, intent(out) :: ierr

        call adios2_variable_check_type(variable, adios2_type_integer2, &
                                        'variable_min', ierr)
        if (ierr == 0) then
            call adios2_variable_min_f2c(minimum, variable%f2c, ierr)
        end if

    end subroutine

    subroutine adios2_variable_min_integer4(minimum, variable, ierr)
        integer(kind=4), intent(out) :: minimum
        type(adios2_variable), intent(in) :: variable
        integer, intent(out) :: ierr

        call adios2_variable_check_type(variable, adios2_type_integer4, &
                                        'variable_min', ierr)
        if (ierr == 0) then
            call adios2_variable_min_f2c(minimum, variable%f2c, ierr)
        end if

    end subroutine

    subroutine adios2_variable_min_integer8(minimum, variable, ierr)
        integer(kind=8), intent(out) :: minimum
        type(adios2_variable), intent(in) :: variable
        integer, intent(out) :: ierr

        call adios2_variable_check_type(variable, adios2_type_integer8, &
                                        'variable_min', ierr)
        if (ierr == 0) then
            call adios2_variable_min_f2c(minimum, variable%f2c, ierr)
        end if

    end subroutine

end module


text:
! RUN: bbc -emit-fir %s -o - | FileCheck %s

! EXPONENT
! CHECK-LABEL: exponent_test
subroutine exponent_test

    integer :: i1, i2, i3, i4
  ! CHECK: %[[i0:.*]] = fir.alloca i32 {bindc_name = "i1", uniq_name = "_QFexponent_testEi1"}
  ! CHECK: %[[i1:.*]] = fir.alloca i32 {bindc_name = "i2", uniq_name = "_QFexponent_testEi2"}
  ! CHECK: %[[i2:.*]] = fir.alloca i32 {bindc_name = "i3", uniq_name = "_QFexponent_testEi3"}
  ! CHECK: %[[i3:.*]] = fir.alloca i32 {bindc_name = "i4", uniq_name = "_QFexponent_testEi4"}

    real(kind = 4) :: x1
    real(kind = 8) :: x2
    real(kind = 10) :: x3
    real(kind = 16) :: x4
  ! CHECK: %[[x0:.*]] = fir.alloca f32 {bindc_name = "x1", uniq_name = "_QFexponent_testEx1"}
  ! CHECK: %[[x1:.*]] = fir.alloca f64 {bindc_name = "x2", uniq_name = "_QFexponent_testEx2"}
  ! CHECK: %[[x2:.*]] = fir.alloca f80 {bindc_name = "x3", uniq_name = "_QFexponent_testEx3"}
  ! CHECK: %[[x3:.*]] = fir.alloca f128 {bindc_name = "x4", uniq_name = "_QFexponent_testEx4"}

    i1 = exponent(x1)
  ! CHECK: %[[temp0:.*]] = fir.load %[[x0:.*]] : !fir.ref<f32>
  ! CHECK: %[[result0:.*]] = fir.call @_FortranAExponent4_4(%[[temp0:.*]]) : (f32) -> i32
  ! CHECK: fir.store %[[result0:.*]] to %[[i0:.*]] : !fir.ref<i32>

    i2 = exponent(x2)
  ! CHECK: %[[temp1:.*]] = fir.load %[[x1:.*]] : !fir.ref<f64>
  ! CHECK: %[[result1:.*]] = fir.call @_FortranAExponent8_4(%[[temp1:.*]]) : (f64) -> i32
  ! CHECK: fir.store %[[result1:.*]] to %[[i1:.*]] : !fir.ref<i32>

    i3 = exponent(x3)
  ! CHECK: %[[temp2:.*]] = fir.load %[[x2:.*]] : !fir.ref<f80>
  ! CHECK: %[[result2:.*]] = fir.call @_FortranAExponent10_4(%[[temp2:.*]]) : (f80) -> i32
  ! CHECK: fir.store %[[result2:.*]] to %[[i2:.*]] : !fir.ref<i32>

    i4 = exponent(x4)
  ! CHECK: %[[temp3:.*]] = fir.load %[[x3:.*]] : !fir.ref<f128>
  ! CHECK: %[[result3:.*]] = fir.call @_FortranAExponent16_4(%[[temp3:.*]]) : (f128) -> i32
  ! CHECK: fir.store %[[result3:.*]] to %[[i3:.*]] : !fir.ref<i32>
  end subroutine exponent_test

max_stars_repo_name: akashhansda/llvm-project

max_issues_repo_name: akashhansda/llvm-project

max_forks_repo_name: akashhansda/llvm-project


text:
( peephole optimizer )

: next-instruction ( codeaddr -- nextcodeaddr )
    dup @ swap cell+ swap
    find-bytecode ?hasarg if
	cell+
    then
;

: ?notendofword ( addr -- addr noteow? )
    dup @ ' eow <> ;

: ?isbranch ( addr -- addr branch? )
    dup @
    case
	' branch of 1 endof
	' 0branch of 1 endof
	' 1branch of 1 endof
	' <branch of 1 endof
	' >branch of 1 endof
	' <=branch of 1 endof
	' >=branch of 1 endof
	' 0<branch of 1 endof
	' 0>branch of 1 endof
	' 0<=branch of 1 endof
	' 0>=branch of 1 endof
	' <>branch of 1 endof
	' =branch of 1 endof
	0 swap
    endcase
;

: fixoffset-afterborder ( border offsetaddr offset -- border offsetaddr fixedoffset )
    dup 0< if           \ backward jump? will it go over the border?
	2dup +          ( border offsetaddr offset targetaddr )
	3 pick < if       \ jump over the border?
	    cell+
	then
    \ no, forward jump. no need to do anything
    then
;

: fixoffset-beforeborder ( border offsetaddr offset -- border offsetaddr fixedoffset )
    dup 0> if
	2dup +
	3 pick > if
	    cell-
	then
	\ forward jump and not over the border
    then
;

: fixoffset ( border offsetaddr offset -- border offsetaddr fixedoffset )
    2 pick 2 pick > if
	fixoffset-beforeborder
    else
	fixoffset-afterborder
    then
;

: fixbranches ( border currinst -- )
    begin
	?notendofword
    while
	    ?isbranch if
		cell+              ( border offsetaddr )
		dup @              ( border offsetaddr offset )
		fixoffset          ( border offsetaddr fixedoffset )
		over !
		cell+
	    else
		next-instruction
	    then
    repeat
    2drop
;

: save-code ( addr code -- addr+cell )
    over ! cell+
;

variable pattern-list

: patterns immediate
    here @   \ mark pattern and replacement list start. consists of cell pairs (pattern,replacement)
;

: end-patterns immediate
    , 0 , 0 ,          \ end of list marker pair
    pattern-list !   \ save pattern and replacement list into pattern list
;

: save-bytecode
    dup           ( wordname wordname )
    find          ( wordname word/0 )
    ?dup if
	nip >cfa @
    else
	number drop
    then
    save-code
;

hex
1badedda constant endpatternmarker
2dabafab constant wildcardmarker
decimal

: ?endofpattern endpatternmarker = ;
: save-wildcard wildcardmarker save-code ;
: terminate-list endpatternmarker save-code consthere ! ;
: wildcard? s" ?" str= ;

\ compile pattern into const space
: p{ immediate ( -- pattern )
    consthere @            ( patt_start )
    dup                    ( patt_start addr )
    begin
	word                 ( patt_start addr word )
	dup wildcard?        ( patt_start addr word iswildcard )
	if
	    drop save-wildcard
	else
	    dup s" }p" str= if
		drop terminate-list exit
	    else
		save-bytecode
	    then
	then
    again
;

\ replacement pattern
: r{ immediate ( -- replacement )
    consthere @
    dup                  ( replacement_start replacement )
    begin
	word             ( replacement_start replacement word )
	dup wildcard?
	if
	    drop save-wildcard
	else
	    dup s" }r" str=   ( replacement_start replacement word isendmarker )
	    if
		drop terminate-list exit
	    else
		save-bytecode
	    then
	then
    again
;

: -> immediate
    , ;

: advance ( addr1 addr2 -- addr1+cell addr2+cell )
    swap cell+ swap cell+ ;

: match ( codeaddr pattern -- true/false )
    begin
	over @ over @          ( codeaddr pattern bytecode patterncode )
	dup ?endofpattern if              \ if we reached pattern's end, that's a match!
	    2drop 2drop 1 exit
	else                   \ otherwise need to compare the codes
	    dup wildcardmarker = if
		\ patterncode wildcardmarker matches always
		2drop
	    else
		<> if             \ no match, exit with false
		    2drop 0 exit
		then
	    then
	    advance
	then
    again
;

: search-match ( codeaddr -- matchcodeaddr replacementaddr/0 )
    pattern-list @         ( codeaddr patternlist )
    begin
	2dup @             ( codeaddr patternlist codeaddr pattern )
	match if           ( codeaddr patternlist matched? )
	    cell+ @ exit
	else
	    cell+ cell+
	    dup @ 0= if
		drop 0 exit
	    then
	then
    again
;

: replace-pattern ( codeaddr replacementaddr -- )
    begin
	2dup @       ( codeaddr replacementaddr codeaddr replacement )
	dup ?endofpattern not while
	    dup wildcardmarker = if
		2drop   \ leave question marks intact
	    else
		swap !
	    then
	    advance
    repeat
    2drop 2drop
;

defer remove-noops

: scan-and-replace ( codeaddr -- )
    dup         ( beginaddr codeaddr )
    begin
	dup @         \ beginaddr codeaddr instruction
	' eow <>      \ beginaddr codeaddr iseow
    while             \ beginaddr codeaddr
	    search-match   \ beginaddr codeaddr match/0
	    ?dup if
		2dup replace-pattern drop          \ beginaddr codeaddr
		over remove-noops    ( remove noops now to make other patterns fire right away )
	    else
		next-instruction
	    then
    repeat
    2drop
;

: copy-instr ( dest src -- dest+cell src+cell )
    2dup @ swap !
    advance
;

: remove-noop ( startaddrofmove -- )
    dup dup cell+            ( startaddrofmove currdest nextinstraddr )
    begin
	dup @ find-bytecode ?hasarg if
	    copy-instr
	else
	    dup @ ' eow = if
		copy-instr 2drop drop exit
	    then
	then
	copy-instr
    again
;

: remove-noops ( wordaddr -- )
    dup                ( wordaddr currinst )
    begin
	dup @          ( wordaddr currinst opcode )
	' eow <>       ( wordaddr currinst isnoteow )
    while
	    dup @ ' noop = if    ( wordaddr currinst isnoop )
		dup 2 pick fixbranches
		dup remove-noop
	    else
		next-instruction
	    then
    repeat
    2drop
;

variable const-expr-pattern
variable const-expr-pattern2
variable const-expr-pattern3

p{ lit ? lit ? ? }p   const-expr-pattern !
p{ lit+ ? lit+ ? }p   const-expr-pattern2 !
p{ lit- ? lit- ? }p   const-expr-pattern3 !

: replace-const-expr ( matchedaddr op -- )
    swap
    dup cell+ @        ( op matchedaddr a )
    over 3 cells + @   ( op matchedaddr a b )
    3 pick exec-builtin ( op matchedaddr a_op_b )
    over cell+ !       ( op matchedaddr )
    2 cells +          ( op firstnoopaddr )
    dup ' noop swap !  ( op firstnoopaddr )
    cell+              ( op 2ndnoopaddr )
    dup ' noop swap !  ( op 2ndnoopaddr )
    cell+              ( op 3rdnoopaddr )
    ' noop swap !      ( op )
    drop
;

: replace-const-expr2 ( matchedaddr -- )
    dup cell+ @       ( matchedaddr a )
    over 3 cells + @  ( matchedaddr a b )
    +
    over cell+ !      ( matchedaddr )
    2 cells +
    dup ' noop swap !
    cell+
    ' noop swap !
;

: constant-fold-replace ( codeaddr -- nomatch )
    dup 4 cells + @ case                 ( codeaddr op )
	' + of ' + replace-const-expr 0 endof
	' - of ' - replace-const-expr 0 endof
	' / of ' / replace-const-expr 0 endof
	' * of ' * replace-const-expr 0 endof
	2drop 1 exit
    endcase
;

: constant-fold-search ( codeaddr -- )
    dup         ( wordaddr codeaddr )
    begin
	dup @                ( wordaddr codeaddr instr )
	' eow <>             ( wordaddr codeaddr noteow )
    while                          ( wordaddr codeaddr )
	    dup const-expr-pattern @ match if        ( wordaddr codeaddr isamatch? )
		dup constant-fold-replace if         ( wordaddr codeaddr notreplaced? )
		    next-instruction     \ not replaced, next instruction
		else
		    over remove-noops    \ replaced, remove noops now
		then
	    else
		dup const-expr-pattern2 @ match if
		    dup replace-const-expr2
		    over remove-noops
		else
		    dup const-expr-pattern3 @ match if
			dup replace-const-expr2
			over remove-noops
		    else
			next-instruction
		    then
		then
	    then
    repeat
    2drop
;

: optimize ( -- )
    latest @ >cfa
    dup constant-fold-search
    dup scan-and-replace
    constant-fold-search
;

: opt-word immediate ( -- )
    word find >cfa
    dup constant-fold-search
    dup scan-and-replace
    constant-fold-search
;

patterns
  p{ lit ? + }p               -> r{ lit+ ? noop }r ,
  p{ lit ? - }p               -> r{ lit- ? noop }r ,
  p{ lit 1 + }p               -> r{ 1+ noop noop }r  ,
  p{ lit 1 - }p               -> r{ 1- noop noop }r  ,
  p{ lit- 1 }p                -> r{ 1- noop }r ,
  p{ lit+ 1 }p                -> r{ 1+ noop }r ,
  p{ lit 1 * }p               -> r{ noop noop noop }r ,
  p{ lit 1 / }p               -> r{ noop noop noop }r ,
  p{ lit 0 + }p               -> r{ noop noop noop }r ,
  p{ lit 0 1- }p              -> r{ lit -1 noop }r ,
  p{ lit 0 1+ }p              -> r{ lit 1 noop }r ,
  p{ lit 0 - }p               -> r{ noop noop noop }r ,
  p{ lit -1 + }p              -> r{ 1- noop noop }r ,
  p{ lit -1 - }p              -> r{ 1+ noop noop }r ,
  p{ flit ? f+ }p             -> r{ flit+ ? noop }r ,
  p{ flit ? f- }p             -> r{ flit- ? noop }r ,
\  p{ swap dup }p              -> r{ swapdup noop }r ,
  p{ swap over }p             -> r{ tuck noop }r   ,
  p{ swap swap }p             -> r{ noop noop }r   ,
  p{ dup swap }p              -> r{ dup noop }r    ,
  p{ drop drop }p             -> r{ 2drop noop }r  ,
  p{ lit ? @ }p               -> r{ var@ ? noop }r  ,
  p{ lit ? ! }p               -> r{ var! ? noop }r  ,
  p{ lit ? f@ }p               -> r{ fvar@ ? noop }r  ,
  p{ lit ? f! }p               -> r{ fvar! ? noop }r  ,
  p{ lit ? s@ }p               -> r{ svar@ ? noop }r  ,
  p{ lit ? s! }p               -> r{ svar! ? noop }r  ,
  p{ lit ? + @ }p             -> r{ field@ ? noop noop }r ,
  p{ lit ? + ! }p             -> r{ field! ? noop noop }r ,
  p{ lit+ ? @ }p              -> r{ field@ ? noop }r ,
  p{ lit+ ? ! }p              -> r{ field! ? noop }r ,

  p{ lit ? + s@ }p             -> r{ sfield@ ? noop noop }r ,
  p{ lit ? + s! }p             -> r{ sfield! ? noop noop }r ,
  p{ lit+ ? s@ }p              -> r{ sfield@ ? noop }r ,
  p{ lit+ ? s! }p              -> r{ sfield! ? noop }r ,

  p{ lit ? + f@ }p             -> r{ ffield@ ? noop noop }r ,
  p{ lit ? + f! }p             -> r{ ffield! ? noop noop }r ,
  p{ lit+ ? f@ }p              -> r{ ffield@ ? noop }r ,
  p{ lit+ ? f! }p              -> r{ ffield! ? noop }r ,

  p{ swap drop swap drop }p   -> r{ 2nip noop noop noop }r ,
  p{ swap drop }p             -> r{ nip noop }r ,
  p{ nip nip nip }p           -> r{ 2nip nip noop }r ,
  p{ nip nip }p               -> r{ 2nip noop }r ,
  p{ over over }p             -> r{ 2dup noop }r ,
  p{ dup @ }p                 -> r{ dup@ noop }r ,
  p{ 0= 0branch ? }p          -> r{ noop 1branch ? }r ,
  p{ 0<> 0branch ? }p         -> r{ noop 0branch ? }r ,
  p{ r> r> }p                 -> r{ 2r> noop }r ,
  p{ >r >r }p                 -> r{ 2>r noop }r ,
  p{ >r r> }p                 -> r{ noop noop }r ,
  p{ r> >r }p                 -> r{ noop noop }r ,
  p{ rdrop rdrop }p           -> r{ 2rdrop noop }r ,
  p{ rot -rot }p              -> r{ noop noop }r ,
  p{ -rot rot }p              -> r{ noop noop }r ,
  p{ dup drop }p              -> r{ noop noop }r ,
  p{ over drop }p             -> r{ noop noop }r ,
  p{ + + }p                   -> r{ bi+ noop }r ,

  p{ 0> 0branch ? }p          -> r{ noop 0<=branch ? }r ,
  p{ 0< 0branch ? }p          -> r{ noop 0>=branch ? }r ,
  p{ 0>= 0branch ? }p          -> r{ noop 0<branch ? }r ,
  p{ 0<= 0branch ? }p          -> r{ noop 0>branch ? }r ,
  p{ 0= 0branch ? }p          -> r{ noop 1branch ? }r ,
  p{ 0<> 0branch ? }p          -> r{ noop 0branch ? }r ,
  p{ = 0branch ? }p           -> r{ noop <>branch ? }r ,
  p{ <> 0branch ? }p          -> r{ noop =branch ? }r ,
  p{ < 0branch ? }p          -> r{ noop >=branch ? }r ,
  p{ > 0branch ? }p          -> r{ noop <=branch ? }r ,
  p{ <= 0branch ? }p          -> r{ noop >branch ? }r ,
  p{ >= 0branch ? }p          -> r{ noop <branch ? }r ,
  p{ = 1branch ? }p           -> r{ noop =branch ? }r ,
  p{ <> 1branch ? }p          -> r{ noop <>branch ? }r ,
  p{ < 1branch ? }p          -> r{ noop <branch ? }r ,
  p{ > 1branch ? }p          -> r{ noop >branch ? }r ,
  p{ <= 1branch ? }p          -> r{ noop <=branch ? }r ,
  p{ >= 1branch ? }p          -> r{ noop >=branch ? }r ,
  p{ 0> 1branch ? }p          -> r{ noop 0>branch ? }r ,
  p{ 0< 1branch ? }p          -> r{ noop 0<branch ? }r ,
  p{ 0>= 1branch ? }p          -> r{ noop 0>=branch ? }r ,
  p{ 0<= 1branch ? }p          -> r{ noop 0<=branch ? }r ,
  p{ 0= 1branch ? }p          -> r{ noop 0branch ? }r ,
  p{ 0<> 1branch ? }p          -> r{ noop 1branch ? }r
end-patterns

s" ;" create immediate
' lit ,
' exit ,
' , ,
' lit ,
' eow ,
' , ,
' latest ,
' @ ,
' hidden ,
' call , ' optimize ,
' [ ,
' exit ,
' eow ,

hide next-instruction
hide ?notendofword
hide ?isbranch
hide fixoffset-afterborder
hide fixoffset-beforeborder
hide fixoffset
hide fixbranches
hide save-code
hide patterns
hide end-patterns
hide save-bytecode
hide save-wildcard
hide terminate-list
hide wildcard?
hide p{
hide r{
hide ->
hide advance
hide match
hide search-match
hide replace-pattern
hide scan-and-replace
hide copy-instr
hide remove-noop
hide remove-noops

max_stars_repo_name: slarba/simpleforth

max_issues_repo_name: slarba/simpleforth

max_forks_repo_name: slarba/simpleforth


text:
C
C      $Id: vc07f.f,v 1.7 2010-03-15 22:49:25 haley Exp $
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C                Copyright (C)  1996                                   C
C        University Corporation for Atmospheric Research               C
C                All Rights Reserved                                   C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C  File:       vc07f.f
C
C  Author:     David Brown (converted by Mary Haley)
C              National Center for Atmospheric Research
C              PO 3000, Boulder, Colorado
C
C  Date:       Wed Jul  3 9:25:32 MST 1996
C
C  Description: This example emulates the LLU example "fcover", overlaying
C               contours and vectors on a map plot.
C
C
      external NhlFAppClass
      external NhlFNcgmWorkstationClass
      external NhlFPSWorkstationClass
      external NhlFPDFWorkstationClass
      external NhlFCairoPSPDFWorkstationClass
      external NhlFCairoImageWorkstationClass
      external NhlFCairoWindowWorkstationClass
      external NhlFVectorFieldClass
      external NhlFVectorPlotClass
      external NhlFScalarFieldClass
      external NhlFContourPlotClass
      external NhlFMapPlotClass

      parameter(MSIZE=73,NSIZE=73,NROWS=11,NCOLORS=24)

      character*7  wks_type
      integer i, j
      integer appid, wid, cnid, vcid, mpid
      integer vfield, sfield
      integer grlist, rlist, len_dims(2)
      real U(NSIZE,MSIZE),V(NSIZE,MSIZE),P(NSIZE,MSIZE)
      integer ithin(NROWS)
      data ithin /90,15,5,5,4,4,3,3,2,2,2/
      character*256 filename
      integer flen
      real cmap(3,NCOLORS)
      data cmap/1.0,1.0,1.0,
     +          0.0,0.0,0.0,
     +          0.9,0.9,0.9,
     +          0.6,0.6,0.6,
     +          0.3,0.3,0.3,
     +          0.8,0.9,1.0,
     +          0.5,0.0,0.5,
     +          0.0,0.5,0.7,
     +          0.0,0.0,0.0,
     +          0.00000,1.00000,0.00000,
     +          0.14286,1.00000,0.00000,
     +          0.28571,1.00000,0.00000,
     +          0.42857,1.00000,0.00000,
     +          0.57143,1.00000,0.00000,
     +          0.71429,1.00000,0.00000,
     +          0.85714,1.00000,0.00000,
     +          1.00000,1.00000,0.00000,
     +          1.00000,0.85714,0.00000,
     +          1.00000,0.71429,0.00000,
     +          1.00000,0.57143,0.00000,
     +          1.00000,0.42857,0.00000,
     +          1.00000,0.28571,0.00000,
     +          1.00000,0.14286,0.00000,
     +          1.00000,0.00000,0.00000/

      len_dims(1) = 3
      len_dims(2) = NCOLORS

C
C Initialize the high level utility library
C
      call NhlFInitialize
C
C Create an application object.
C
      call NhlFRLCreate(rlist,'setrl')
      call NhlFRLCreate(grlist,'getrl')

      call NhlFRLClear(rlist)
      call NhlFRLSetString(rlist,'appUsrDir','./',ierr)
      call NhlFRLSetString(rlist,'appDefaultParent','True',ierr)
      call NhlFCreate(appid,'vc07',NhlFappClass,0,rlist,ierr)

      wks_type = "x11"

      if (wks_type.eq."ncgm".or.wks_type.eq."NCGM") then
C
C Create an NCGM workstation.
C
         call NhlFRLClear(rlist)
         call NhlFRLSetString(rlist,'wkMetaName','./vc07f.ncgm',ierr)
         call NhlFRLSetMDFloatArray(rlist,'wkColorMap',cmap,2,len_dims,
     +        ierr)
         call NhlFCreate(wid,'vc07Work',
     +        NhlFNcgmWorkstationClass,0,rlist,ierr)
      else if (wks_type.eq."x11".or.wks_type.eq."X11") then
C
C Create an X11 workstation.
C
         call NhlFRLClear(rlist)
         call NhlFRLSetString(rlist,'wkPause','True',ierr)
         call NhlFRLSetMDFloatArray(rlist,'wkColorMap',cmap,2,len_dims,
     +        ierr)
         call NhlFCreate(wid,'vc07Work',
     +        NhlFCairoWindowWorkstationClass,
     +        0,rlist,ierr)
      else if (wks_type.eq."oldps".or.wks_type.eq."OLDPS") then
C
C Create an older-style PostScript workstation.
C
         call NhlFRLClear(rlist)
         call NhlFRLSetString(rlist,'wkPSFileName','./vc07f.ps',ierr)
         call NhlFRLSetMDFloatArray(rlist,'wkColorMap',cmap,2,len_dims,
     +        ierr)
         call NhlFCreate(wid,'vc07Work',
     +        NhlFPSWorkstationClass,0,rlist,ierr)
      else if (wks_type.eq."oldpdf".or.wks_type.eq."OLDPDF") then
C
C Create an older-style PDF workstation.
C
         call NhlFRLClear(rlist)
         call NhlFRLSetString(rlist,'wkPDFFileName','./vc07f.pdf',ierr)
         call NhlFRLSetMDFloatArray(rlist,'wkColorMap',cmap,2,len_dims,
     +        ierr)
         call NhlFCreate(wid,'vc07Work',
     +        NhlFPDFWorkstationClass,0,rlist,ierr)
      else if (wks_type.eq."pdf".or.wks_type.eq."PDF".or.
     +         wks_type.eq."ps".or.wks_type.eq."PS") then
C
C Create a cairo PS/PDF workstation.
C
         call NhlFRLClear(rlist)
         call NhlFRLSetString(rlist,'wkFileName','./vc07f',ierr)
         call NhlFRLSetMDFloatArray(rlist,'wkColorMap',cmap,2,len_dims,
     +        ierr)
         call NhlFRLSetString(rlist,'wkFormat',wks_type,ierr)
         call NhlFCreate(wid,'vc07Work',
     +        NhlFCairoPSPDFWorkstationClass,0,rlist,ierr)
      else if (wks_type.eq."png".or.wks_type.eq."PNG") then
C
C Create a cairo PNG workstation.
C
         call NhlFRLClear(rlist)
         call NhlFRLSetString(rlist,'wkFileName','./vc07f',ierr)
         call NhlFRLSetMDFloatArray(rlist,'wkColorMap',cmap,2,len_dims,
     +        ierr)
         call NhlFRLSetString(rlist,'wkFormat',wks_type,ierr)
         call NhlFCreate(wid,'vc07Work',
     +        NhlFCairoImageWorkstationClass,0,rlist,ierr)
      endif
C
C Read the data file
C
      flen = 10
      call gngpat(filename,'examples',ierr)
      do 10 i=1,256
         if( filename(i:i).eq.char(0) ) then
            filename(i:i+flen) = '/fcover.dat'
            goto 15
         endif
 10   continue
 15   open(unit=10,file=filename,status='old')
      read (10,*) U
      read (10,*) V
      read (10,*) P
C
C Massage the data to eliminate surplus of vectors near the pole
C
      do 50 j=NSIZE,NSIZE-NROWS+1,-1
         do 49 i=1,MSIZE
            if (mod(i,ithin(NSIZE-j+1)).ne.0) then
               u(i,j) = -9999.0
            endif
 49      continue
 50   continue
C
C Create a MapPlot object.
C
      call NhlFCreate(mpid,'mapplot',NhlFmapPlotClass,wid,0,ierr)
C
C Create a ScalarField.
C
      len_dims(1) = NSIZE
      len_dims(2) = MSIZE
      call NhlFRLClear(rlist)
      call NhlFRLSetMDFloatArray(rlist,'sfDataArray',P,2,len_dims,ierr)
      call NhlFRLSetFloat(rlist,'sfXCStartV',-180.,ierr)
      call NhlFRLSetFloat(rlist,'sfYCStartV',-90.,ierr)
      call NhlFRLSetFloat(rlist,'sfXCEndV',180.,ierr)
      call NhlFRLSetFloat(rlist,'sfYCEndV',90.,ierr)
      call NhlFRLSetFloat(rlist,'sfMissingValueV', -9999.0,ierr)
      call NhlFCreate(sfield,'ScalarField',NhlFscalarFieldClass,appid,
     +     rlist,ierr)
C
C Create a VectorField.
C
      call NhlFRLClear(rlist)
      call NhlFRLSetMDFloatArray(rlist,'vfUDataArray',U,2,len_dims,ierr)
      call NhlFRLSetMDFloatArray(rlist,'vfVDataArray',V,2,len_dims,ierr)
      call NhlFRLSetFloat(rlist,'vfXCStartV',-180.,ierr)
      call NhlFRLSetFloat(rlist,'vfYCStartV',-90.,ierr)
      call NhlFRLSetFloat(rlist,'vfXCEndV', 180.,ierr)
      call NhlFRLSetFloat(rlist,'vfYCEndV',90.,ierr)
      call NhlFRLSetFloat(rlist,'vfMissingUValueV',-9999.0,ierr)
      call NhlFRLSetFloat(rlist,'vfMissingVValueV',-9999.0,ierr)
      call NhlFCreate(vfield,'VectorField',NhlFvectorFieldClass,appid,
     +     rlist,ierr)
C
C Create a VectorPlot object.
C
      call NhlFRLClear(rlist)
      call NhlFRLSetString(rlist,'vcUseScalarArray','true',ierr)
      call NhlFRLSetInteger(rlist,'vcVectorFieldData',vfield,ierr)
      call NhlFRLSetInteger(rlist,'vcScalarFieldData',sfield,ierr)
      call NhlFCreate(vcid,'vectorplot',NhlFvectorPlotClass,wid,rlist,
     +     ierr)

      call NhlFRLClear(grlist)
      call NhlFRLGetFloat(grlist,'vcMinMagnitudeF',vmin,ierr)
      call NhlFRLGetFloat(grlist,'vcMaxMagnitudeF',vmax,ierr)
      call NhlFGetValues(vcid,grlist,ierr)

      call NhlFRLClear(rlist)
      call NhlFRLSetFloat(rlist,'vcMinMagnitudeF',vmin+0.1*(vmax-vmin),
     +     ierr)
      call NhlFSetValues(vcid,rlist,ierr)
C
C Create a ContourPlot object.
C
      call NhlFRLClear(rlist)
      call NhlFRLSetInteger(rlist,'cnScalarFieldData',sfield,ierr)
      call NhlFCreate(cnid,'contourplot',NhlFcontourPlotClass,wid,
     +     rlist,ierr)
C
C Overlay the vectors and the contours on the map and draw
C everything.
C
      call NhlFDraw(mpid,ierr)
      call NhlFFrame(wid,ierr)

      call NhlFDraw(cnid,ierr)
      call NhlFFrame(wid,ierr)

      call NhlFDraw(vcid,ierr)
      call NhlFFrame(wid,ierr)

      call NhlFAddOverlay(mpid,vcid,-1,ierr)
      call NhlFAddOverlay(mpid,cnid,-1,ierr)
      call NhlFDraw(mpid,ierr)
      call NhlFFrame(wid,ierr)

      call NhlFDestroy(wid,ierr)
      call NhlFClose
      stop
      end

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


! { dg-do compile }
!
! Contributed by Hossein Talebi  <talebi.hossein@gmail.com>
!
!
Module part_base2_class

    implicit none

    type :: ty_moc1
        integer l
    end type ty_moc1
    integer,parameter ::  MAX_NUM_ELEMENT_TYPE=32

    type :: ty_element_index2

        class(ty_moc1),allocatable :: element
        class(ty_moc1),allocatable :: element_th(:)

    endtype ty_element_index2

    type :: ty_part_base2
        type(ty_element_index2)::element_index(MAX_NUM_ELEMENT_TYPE)
    end type ty_part_base2

    class(ty_part_base2),allocatable ::  part_tmp_obj

End Module part_base2_class

    use part_base2_class
    allocate (part_tmp_obj)
    allocate (part_tmp_obj%element_index(1)%element, source = ty_moc1(1))
    allocate (part_tmp_obj%element_index(1)%element_th(1), source = ty_moc1(99))
    allocate (part_tmp_obj%element_index(32)%element_th(1), source = ty_moc1(999))

    do i = 1, MAX_NUM_ELEMENT_TYPE
      if (allocated (part_tmp_obj%element_index(i)%element_th)) then
        print *, i, part_tmp_obj%element_index(i)%element_th(1)%l
      end if
    end do
    deallocate (part_tmp_obj)

end