!
!     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!     *                                                               *
!     *                  copyright (c) 1998 by UCAR                   *
!     *                                                               *
!     *       University Corporation for Atmospheric Research         *
!     *                                                               *
!     *                      all rights reserved                      *
!     *                                                               *
!     *                         Spherepack                            *
!     *                                                               *
!     *       A Package of Fortran Subroutines and Programs           *
!     *                                                               *
!     *              for Modeling Geophysical Processes               *
!     *                                                               *
!     *                             by                                *
!     *                                                               *
!     *                  John Adams and Paul Swarztrauber             *
!     *                                                               *
!     *                             of                                *
!     *                                                               *
!     *         the National Center for Atmospheric Research          *
!     *                                                               *
!     *                Boulder, Colorado  (80307)  U.S.A.             *
!     *                                                               *
!     *                   which is sponsored by                       *
!     *                                                               *
!     *              the National Science Foundation                  *
!     *                                                               *
!     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!
module type_AssociatedLegendrePolynomialGenerator

    use spherepack_precision, only: &
        wp, & ! working precision
        ip, & ! integer precision
        PI, HALF_PI

    use type_FastFourierTransform, only: &
        FastFourierTransform

    ! Explicit typing only
    implicit none

    ! Everything is private unless stated otherwise
    private
    public :: AssociatedLegendrePolynomialGenerator
    public :: alfk, lfp, lfpt, lfim, lfin, get_legendre_function

    type, public :: AssociatedLegendrePolynomialGenerator
    contains
        ! Type-bound procedures
        procedure, nopass :: alfk
        procedure, nopass :: lfim
        procedure, nopass :: lfin
        procedure, nopass :: lfp
        procedure, nopass :: lfpt
        procedure, nopass :: get_legendre_function
    end type AssociatedLegendrePolynomialGenerator

contains

    subroutine get_legendre_function(lat, ntrunc, legfunc)

        ! Dummy arguments
        real(wp),              intent(in)  :: lat
        integer(ip),           intent(in)  :: ntrunc
        real(wp), allocatable, intent(out) :: legfunc(:)

        ! Local variables
        integer(ip)  :: alloc_stat, required_size, workspace_size

        !  Allocate memory
        required_size = (ntrunc + 1) * (ntrunc + 2)/2
        allocate (legfunc(required_size), stat=alloc_stat)

        ! Check allocation status
        if (alloc_stat /= 0) then
            error stop "Failed to allocate legfunc in get_legendre_function"
        end if

        ! Compute required workspace size
        workspace_size = (ntrunc/2)+1

        block
            real(wp), parameter :: DEGREE_TO_RADIAN = PI/180
            real(wp)            :: cp(workspace_size), theta
            integer(ip)         :: n, m, nm, nmstrt ! Counters

            theta = HALF_PI - (DEGREE_TO_RADIAN * lat)
            nmstrt = 0

            do m=1, ntrunc+1
                do n=m, ntrunc+1
                    nm = nmstrt + n - m + 1

                    !  Compute normalized associate Legendre function at theta
                    call alfk(n - 1, m - 1, cp)
                    call lfpt(n - 1, m - 1, theta, cp, legfunc(nm))

                end do
                nmstrt = nmstrt + ntrunc - m + 2
            end do
        end block

    end subroutine get_legendre_function

    subroutine alfk(n, m, cp)
        ! subroutine alfk (n, m, cp)
        !
        ! dimension of           real cp(n/2 + 1)
        ! arguments
        !
        ! purpose                routine alfk computes double precision fourier
        !                        coefficients in the trigonometric series
        !                        representation of the normalized associated
        !                        legendre function pbar(n, m, theta) for use by
        !                        routines lfp and lfpt in calculating double
        !                        precision pbar(n, m, theta).
        !
        !                        first define the normalized associated
        !                        legendre functions
        !
        !                        pbar(m, n, theta) = sqrt((2*n+1)*factorial(n-m)
        !                        /(2*factorial(n+m)))*sin(theta)**m/(2**n*
        !                        factorial(n)) times the (n+m)th derivative of
        !                        (x**2-1)**n with respect to x=cos(theta)
        !
        !                        where theta is colatitude.
        !
        !                        then subroutine alfk computes the coefficients
        !                        cp(k) in the following trigonometric
        !                        expansion of pbar(m, n, theta).
        !
        !                        1) for n even and m even, pbar(m, n, theta) =
        !                           .5*cp(1) plus the sum from k=1 to k=n/2
        !                           of cp(k+1)*cos(2*k*th)
        !
        !                        2) for n even and m odd, pbar(m, n, theta) =
        !                           the sum from k=1 to k=n/2 of
        !                           cp(k)*sin(2*k*th)
        !
        !                        3) for n odd and m even, pbar(m, n, theta) =
        !                           the sum from k=1 to k=(n+1)/2 of
        !                           cp(k)*cos((2*k-1)*th)
        !
        !                        4) for n odd and m odd,  pbar(m, n, theta) =
        !                           the sum from k=1 to k=(n+1)/2 of
        !                           cp(k)*sin((2*k-1)*th)
        !
        !
        ! usage                  call alfk(n, m, cp)
        !
        ! arguments
        !
        ! on input               n
        !                          nonnegative integer specifying the degree of
        !                          pbar(n, m, theta)
        !
        !                        m
        !                          is the order of pbar(n, m, theta). m can be
        !                          any integer however cp is computed such that
        !                          pbar(n, m, theta) = 0 if abs(m) is greater
        !                          than n and pbar(n, m, theta) = (-1)**m*
        !                          pbar(n, -m, theta) for negative m.
        !
        ! on output              cp
        !                          double precision array of length (n/2)+1
        !                          which contains the fourier coefficients in
        !                          the trigonometric series representation of
        !                          pbar(n, m, theta)
        !
        !
        ! special conditions     none
        !
        ! precision              64-bit double precision
        !
        ! algorithm              the highest order coefficient is determined in
        !                        closed form and the remainig coefficients are
        !                        determined as the solution of a backward
        !                        recurrence relation.
        !
        ! accuracy               comparison between routines alfk and double
        !                        precision dalfk on the cray1 indicates
        !                        greater accuracy for smaller values
        !                        of input parameter n.  agreement to 14
        !                        places was obtained for n=10 and to 13
        !                        places for n=100.
        !

        ! Dummy arguments
        integer(ip), intent(in)   :: n
        integer(ip), intent(in)   :: m
        real(wp),    intent(out)  :: cp(n/2+1)

        ! Local variables
        integer(ip)         :: i, l, ma, nex, nmms2
        real(wp), parameter :: sc10 = 1024.0_wp
        real(wp), parameter :: sc20 = sc10**2
        real(wp), parameter :: sc40 = sc20**2
        real(wp)            :: a1, b1, c1, t1, t2
        real(wp)            :: fk, cp2, pm1
        real(wp)            :: fden, fnmh, fnum, fnnp1, fnmsq

        cp(1) = 0.0_wp
        ma = abs(m)

        if (ma > n) return

        if (n < 1) then
            !
            !  n less than 1
            !
            cp(1) = sqrt(2.0_wp)
        else if (n == 1) then
            !
            !  n equals 1
            !
            if (ma == 0) then
                cp(1) = sqrt(1.5_wp)
            else
                cp(1) = sqrt(0.75_wp)
                select case (m)
                    case (-1)
                        cp(1) = -cp(1)
                end select
            end if
        else
            !
            !  n greater than 1
            !
            select case (mod(n+ma, 2))
                case (0)
                    nmms2 = (n-ma)/2
                    fnum = n+ma+1
                    fnmh = n-ma+1
                    pm1 = 1.0_wp
                case default
                    nmms2 = (n-ma-1)/2
                    fnum = n+ma+2
                    fnmh = n-ma+2
                    pm1 = -1.0_wp
            end select

            t1 = 1.0_wp/sc20
            nex = 20
            fden = 2.0_wp

            if (nmms2 >= 1) then
                do  i=1, nmms2
                    t1 = fnum*t1/fden
                    if (t1 > sc20) then
                        t1 = t1/sc40
                        nex = nex+40
                    end if
                    fnum = fnum+2.0_wp
                    fden = fden+2.0_wp
                end do
            end if

            t1 = t1/2.0_wp**(n-1-nex)

            if (mod(ma/2, 2) /= 0) t1 = -t1

            t2 = 1.0_wp

            if (ma /= 0) then
                do  i=1, ma
                    t2 = fnmh*t2/(fnmh+pm1)
                    fnmh = fnmh+2.0_wp
                end do
            end if

            cp2 = t1*sqrt((real(n, kind=wp)+0.5_wp)*t2)
            fnnp1 = real(n*(n+1), kind=wp)
            fnmsq = fnnp1-2.0_wp*(ma**2)
            l = (n+1)/2

            if (mod(n, 2) == 0 .and. mod(ma, 2) == 0) l = l+1

            cp(l) = cp2

            if ((m < 0) .and. (mod(ma, 2) /= 0)) cp(l) = -cp(l)

            if (l <= 1) return

            fk = n
            a1 = (fk-2.0_wp)*(fk-1.0_wp)-fnnp1
            b1 = 2.0_wp*(fk**2-fnmsq)
            cp(l-1) = b1*cp(l)/a1

            do
                l = l-1

                if (l <= 1) return

                fk = fk-2.0_wp
                a1 = (fk-2.0_wp)*(fk-1.0_wp)-fnnp1
                b1 = -2.0_wp*(fk*fk-fnmsq)
                c1 = (fk+1.0_wp)*(fk+2.0_wp)-fnnp1
                cp(l-1) = -(b1*cp(l)+c1*cp(l+1))/a1
                cycle
            end do
        end if

    end subroutine alfk

    subroutine lfim(init, theta, l, n, nm, pb, id, wlfim)
        !
        ! subroutine lfim (init, theta, l, n, nm, pb, id, wlfim)
        !
        ! dimension of           theta(l),  pb(id, nm+1),  wlfim(4*l*(nm+1))
        ! arguments
        !
        ! purpose                given n and l, routine lfim calculates
        !                        the normalized associated legendre functions
        !                        pbar(n, m, theta) for m=0, ..., n and theta(i)
        !                        for i=1, ..., l where
        !
        !                        pbar(m, n, theta) = sqrt((2*n+1)*factorial(n-m)
        !                        /(2*factorial(n+m)))*sin(theta)**m/(2**n*
        !                        factorial(n)) times the (n+m)th derivative of
        !                        (x**2-1)**n with respect to x=cos(theta)
        !
        ! usage                  call lfim (init, theta, l, n, nm, pb, id, wlfim)
        !
        ! arguments
        ! on input               init
        !                        = 0
        !                            initialization only - using parameters
        !                            l, nm and array theta, subroutine lfim
        !                            initializes array wlfim for subsequent
        !                            use in the computation of the associated
        !                            legendre functions pb. initialization
        !                            does not have to be repeated unless
        !                            l, nm, or array theta are changed.
        !                        = 1
        !                            subroutine lfim uses the array wlfim that
        !                            was computed with init = 0 to compute pb.
        !
        !                        theta
        !                          an array that contains the colatitudes
        !                          at which the associated legendre functions
        !                          will be computed. the colatitudes must be
        !                          specified in radians.
        !
        !                        l
        !                          the length of the theta array. lfim is
        !                          vectorized with vector length l.
        !
        !                        n
        !                          nonnegative integer, less than nm, specifying
        !                          degree of pbar(n, m, theta). subroutine lfim
        !                          must be called starting with n=0. n must be
        !                          incremented by one in subsequent calls and
        !                          must not exceed nm.
        !
        !                        nm
        !                          the maximum value of n and m
        !
        !                        id
        !                          the first dimension of the two dimensional
        !                          array pb as it appears in the program that
        !                          calls lfim. (see output parameter pb)
        !
        !                        wlfim
        !                          an array with length 4*l*(nm+1) which
        !                          must be initialized by calling lfim
        !                          with init=0 (see parameter init)  it
        !                          must not be altered between calls to
        !                          lfim.
        !
        !
        ! on output              pb
        !                          a two dimensional array with first
        !                          dimension id in the program that calls
        !                          lfim. the second dimension of pb must
        !                          be at least nm+1. starting with n=0
        !                          lfim is called repeatedly with n being
        !                          increased by one between calls. on each
        !                          call, subroutine lfim computes
        !                          = pbar(m, n, theta(i)) for m=0, ..., n and
        !                          i=1, ...l.
        !
        !                        wlfim
        !                          array containing values which must not
        !                          be altered unless l, nm or the array theta
        !                          are changed in which case lfim must be
        !                          called with init=0 to reinitialize the
        !                          wlfim array.
        !
        ! special conditions     n must be increased by one between calls
        !                        of lfim in which n is not zero.
        !
        ! precision              64-bit double precision
        !
        !
        ! algorithm              routine lfim calculates pbar(n, m, theta) using
        !                        a four term recurrence relation. (unpublished
        !                        notes by paul n. swarztrauber)
        !

        ! Dummy arguments

        integer(ip), intent(in)  :: init
        real(wp),    intent(in)  :: theta(*)
        integer(ip), intent(in)  :: l
        integer(ip), intent(in)  :: n
        integer(ip), intent(in)  :: nm
        real(wp),    intent(out) :: pb(1)
        integer(ip), intent(in)  :: id
        real(wp),    intent(out) :: wlfim(1)

        ! Local variables

        integer(ip) :: workspace_indices(3)

        !
        !   total length of wlfim is 4*l*(nm+1)
        !
        workspace_indices = get_workspace_indices(l, nm)

        associate (&
            iw1 => workspace_indices(1), &
            iw2 => workspace_indices(2), &
            iw3 => workspace_indices(3) &
            )

            call lfim_lower_utility_routine(init, theta, l, n, nm, id, pb, wlfim, wlfim(iw1), &
                wlfim(iw2), wlfim(iw3), wlfim(iw2))

        end associate

    contains

        pure function get_workspace_indices(l, nm) result (return_value)

            ! Dummy arguments

            integer(ip), intent(in) :: l
            integer(ip), intent(in) :: nm
            integer(ip)              :: return_value(3)

            ! Local variables

            integer(ip) :: lnx

            associate (i => return_value)

                lnx = l*(nm + 1)
                i(1) = lnx + 1
                i(2) = 2*lnx + 1
                i(3) = 3*lnx + 1

            end associate

        end function get_workspace_indices

        subroutine lfim_lower_utility_routine(init, theta, l, n, nm, id, p3, phz, ph1, p1, p2, cp)

            ! Dummy arguments

            integer(ip), intent(in)  :: init
            real(wp),    intent(in)  :: theta(*)
            integer(ip), intent(in)  :: l
            integer(ip), intent(in)  :: n
            integer(ip), intent(in)  :: nm
            integer(ip), intent(in)  :: id
            real(wp),    intent(out) :: p3(id, *)
            real(wp),    intent(out) :: phz(l, *)
            real(wp),    intent(out) :: ph1(l, *)
            real(wp),    intent(out) :: p1(l, *)
            real(wp),    intent(out) :: p2(l, *)
            real(wp),    intent(out) :: cp(*)

            ! Local variables

            integer(ip)         :: i, m, nm1, nh, mp1, np1, nmp1
            real(wp)            :: cc, dd, ee, cn, fm, fn, fnmm, fnpm
            real(wp)            :: tn, temp
            real(wp), parameter :: SQRT2 = sqrt(2.0_wp)
            real(wp), parameter :: SQRT5 = sqrt(5.0_wp)
            real(wp), parameter :: SQRT6 = sqrt(6.0_wp)
            real(wp), parameter :: ONE_OVER_SQRT2 = 1.0_wp/SQRT2
            real(wp), parameter :: ONE_OVER_SQRT6 = 1.0_wp/SQRT6
            real(wp), parameter :: SQRT5_OVER_SQRT6 = SQRT5/SQRT6

            nmp1 = nm+1

            select case (init)
                case (0)
                    phz(:, 1) = ONE_OVER_SQRT2
                    do np1=2, nmp1
                        nh = np1-1
                        call alfk(nh, 0, cp)
                        do i=1, l
                            call lfpt(nh, 0, theta(i), cp, phz(i, np1))
                        end do
                        call alfk(nh, 1, cp)
                        do  i=1, l
                            call lfpt(nh, 1, theta(i), cp, ph1(i, np1))
                        end do
                    end do
                case default
                    if (n <= 2) then
                        if (n < 1) then
                            p3(:, 1) = phz(:, 1)
                        else if (n == 1) then
                            p3(:, 1) = phz(:, 2)
                            p3(:, 2) = ph1(:, 2)
                        else
                            p3(:, 1) = phz(:, 3)
                            p3(:, 2) = ph1(:, 3)
                            p3(:, 3) = SQRT5_OVER_SQRT6 * phz(:, 1) - ONE_OVER_SQRT6 * p3(:, 1)
                            p1(:, 1) = phz(:, 2)
                            p1(:, 2) = ph1(:, 2)
                            p2(:, 1) = phz(:, 3)
                            p2(:, 2) = ph1(:, 3)
                            p2(:, 3) = p3(:, 3)
                        end if
                    else
                        nm1 = n-1
                        np1 = n+1
                        fn = real(n, kind=wp)
                        tn = 2.0_wp*fn
                        cn = (tn+1.0_wp)/(tn-3.0_wp)
                        p3(:, 1) = phz(:, np1)
                        p3(:, 2) = ph1(:, np1)

                        if (nm1 >= 3) then
                            do mp1=3, nm1
                                m = mp1-1
                                fm = real(m, kind=wp)
                                fnpm = fn+fm
                                fnmm = fn-fm
                                temp = fnpm*(fnpm-1.0_wp)
                                cc = sqrt(cn*(fnpm-3.0_wp)*(fnpm-2.0_wp)/temp)
                                dd = sqrt(cn*fnmm*(fnmm-1.0_wp)/temp)
                                ee = sqrt((fnmm+1.0_wp)*(fnmm+2.0_wp)/temp)
                                p3(:, mp1) = cc*p1(i, mp1-2)+dd*p1(:, mp1)-ee*p3(:, mp1-2)
                            end do
                        end if
                        fnpm = 2.0_wp*fn-1.0_wp
                        temp = fnpm*(fnpm-1.0_wp)
                        cc = sqrt(cn*(fnpm-3.0_wp)*(fnpm-2.0_wp)/temp)
                        ee = sqrt(6.0_wp/temp)
                        p3(:, n) = cc*p1(:, n-2)-ee*p3(:, n-2)
                        fnpm = 2.0_wp*fn
                        temp = fnpm*(fnpm-1.0_wp)
                        cc = sqrt(cn*(fnpm-3.0_wp)*(fnpm-2.0_wp)/temp)
                        ee = sqrt(2.0_wp/temp)
                        p3(:, n+1) = cc*p1(:, n-1)-ee*p3(:, n-1)
                        p1(:, 1:np1) = p2(:, 1:np1)
                        p2(:, :np1) = p3(:, 1:np1)
                    end if
            end select

        end subroutine lfim_lower_utility_routine

    end subroutine lfim

    subroutine lfin(init, theta, l, m, nm, pb, id, wlfin)
        !
        ! subroutine lfin (init, theta, l, m, nm, pb, id, wlfin)
        !
        ! dimension of           theta(l),  pb(id, nm+1),  wlfin(4*l*(nm+1))
        ! arguments
        !
        ! purpose                given m and l, routine lfin calculates
        !                        the normalized associated legendre functions
        !                        pbar(n, m, theta) for n=m, ..., nm and theta(i)
        !                        for i=1, ..., l where
        !
        !                        pbar(m, n, theta) = sqrt((2*n+1)*factorial(n-m)
        !                        /(2*factorial(n+m)))*sin(theta)**m/(2**n*
        !                        factorial(n)) times the (n+m)th derivative of
        !                        (x**2-1)**n with respect to x=cos(theta)
        !
        ! usage                  call lfin (init, theta, l, m, nm, pb, id, wlfin)
        !
        ! arguments
        ! on input               init
        !                        = 0
        !                            initialization only - using parameters
        !                            l, nm and the array theta, subroutine lfin
        !                            initializes the array wlfin for subsequent
        !                            use in the computation of the associated
        !                            legendre functions pb. initialization does
        !                            not have to be repeated unless l, nm or
        !                            the array theta are changed.
        !                        = 1
        !                            subroutine lfin uses the array wlfin that
        !                            was computed with init = 0 to compute pb
        !
        !                        theta
        !                          an array that contains the colatitudes
        !                          at which the associated legendre functions
        !                          will be computed. the colatitudes must be
        !                          specified in radians.
        !
        !                        l
        !                          the length of the theta array. lfin is
        !                          vectorized with vector length l.
        !
        !                        m
        !                          nonnegative integer, less than nm, specifying
        !                          degree of pbar(n, m, theta). subroutine lfin
        !                          must be called starting with n=0. n must be
        !                          incremented by one in subsequent calls and
        !                          must not exceed nm.
        !
        !                        nm
        !                          the maximum value of n and m
        !
        !                        id
        !                          the first dimension of the two dimensional
        !                          array pb as it appears in the program that
        !                          calls lfin. (see output parameter pb)
        !
        !                        wlfin
        !                          an array with length 4*l*(nm+1) which
        !                          must be initialized by calling lfin
        !                          with init=0 (see parameter init)  it
        !                          must not be altered between calls to
        !                          lfin.
        !
        !
        ! on output              pb
        !                          a two dimensional array with first
        !                          dimension id in the program that calls
        !                          lfin. the second dimension of pb must
        !                          be at least nm+1. starting with m=0
        !                          lfin is called repeatedly with m being
        !                          increased by one between calls. on each
        !                          call, subroutine lfin computes pb(i, n+1)
        !                          = pbar(m, n, theta(i)) for n=m, ..., nm and
        !                          i=1, ...l.
        !
        !                        wlfin
        !                          array containing values which must not
        !                          be altered unless l, nm or the array theta
        !                          are changed in which case lfin must be
        !                          called with init=0 to reinitialize the
        !                          wlfin array.
        !
        ! special conditions     m must be increased by one between calls
        !                        of lfin in which m is not zero.
        !
        ! precision              single
        !
        ! algorithm              routine lfin calculates pbar(n, m, theta) using
        !                        a four term recurrence relation. (unpublished
        !                        notes by paul n. swarztrauber)
        !

        ! Dummy arguments
        integer(ip), intent(in)  :: init
        real(wp),    intent(in)  :: theta(l)
        integer(ip), intent(in)  :: l
        integer(ip), intent(in)  :: m
        integer(ip), intent(in)  :: nm
        real(wp),    intent(out) :: pb(id, nm+1)
        integer(ip), intent(in)  :: id
        real(wp),    intent(out) :: wlfin(4*l*(nm+1))

        ! Local variables
        integer(ip) :: workspace_indices(3)

        !  total length of wlfin is 4*l*(nm+1)
        workspace_indices = get_workspace_indices(l, nm)

        associate (&
            iw1=> workspace_indices(1), &
            iw2 => workspace_indices(2), &
            iw3 => workspace_indices(3) &
            )
            call lfin_lower_utility_routine(init, theta, l, m, nm, id, pb, wlfin, wlfin(iw1), &
                wlfin(iw2), wlfin(iw3), wlfin(iw2))
        end associate

    contains

        pure function get_workspace_indices(l, nm) result (return_value)

            ! Dummy arguments
            integer(ip), intent(in) :: l
            integer(ip), intent(in) :: nm
            integer(ip)              :: return_value(3)

            ! Local variables
            integer(ip) :: lnx

            associate (i => return_value)
                lnx = l*(nm + 1)
                i(1) = lnx + 1
                i(2) = 2*lnx + 1
                i(3) = 3*lnx + 1
            end associate

        end function get_workspace_indices

        subroutine lfin_lower_utility_routine(init, theta, l, m, nm, id, p3, phz, ph1, p1, p2, cp)

            ! Dummy arguments
            integer(ip), intent(in)  :: init
            real(wp),    intent(in)  :: theta(l)
            integer(ip), intent(in)  :: l
            integer(ip), intent(in)  :: m
            integer(ip), intent(in)  :: nm
            integer(ip), intent(in)  :: id
            real(wp),    intent(out) :: p3(id, *)
            real(wp),    intent(out) :: phz(l, *)
            real(wp),    intent(out) :: ph1(l, *)
            real(wp),    intent(out) :: p1(l, *)
            real(wp),    intent(out) :: p2(l, *)
            real(wp),    intent(out) :: cp(*)

            ! Local variables
            integer(ip)         :: i, n, nh, mp1, np1, mp3, nmp1
            real(wp)            :: cc, dd, ee, cn, fm, fn, fnmm, fnpm
            real(wp)            :: tm, tn, temp
            real(wp), parameter :: SQRT2 = sqrt(2.0_wp)
            real(wp), parameter :: ONE_OVER_SQRT2 = 1.0_wp/sqrt2

            nmp1 = nm+1

            select case (init)
                case (0)

                    phz(:, 1) = ONE_OVER_SQRT2

                    do np1=2, nmp1
                        nh = np1-1
                        call alfk(nh, 0, cp)
                        do i=1, l
                            call lfpt(nh, 0, theta(i), cp, phz(i, np1))
                        end do
                        call alfk(nh, 1, cp)
                        do i=1, l
                            call lfpt(nh, 1, theta(i), cp, ph1(i, np1))
                        end do
                    end do

                case default
                    mp1 = m+1
                    fm = real(m)
                    tm = fm+fm

                    if (m < 1) then

                        p3(:, 1:nmp1) = phz(:, 1:nmp1)
                        p1(:, 1:nmp1) = phz(:, 1:nmp1)

                    else if (m == 1) then

                        p3(:, 2:nmp1) = ph1(:, 2:nmp1)
                        p2(:, 2:nmp1) = ph1(:, 2:nmp1)

                    else
                        temp = tm*(tm-1.0_wp)
                        cc = sqrt((tm+1.0_wp)*(tm-2.0_wp)/temp)
                        ee = sqrt(2.0_wp/temp)
                        p3(:, m+1) = cc*p1(:, m-1)-ee*p1(:, m+1)

                        if (m == nm) return

                        temp = tm*(tm+1.0_wp)
                        cc = sqrt((tm+3.0_wp)*(tm-2.0_wp)/temp)
                        ee = sqrt(6.0_wp/temp)
                        p3(:, m+2) = cc*p1(:, m)-ee*p1(:, m+2)
                        mp3 = m+3

                        if (nmp1 >= mp3) then
                            do np1=mp3, nmp1
                                n = np1-1
                                fn = real(n, kind=wp)
                                tn = fn+fn
                                cn = (tn+1.0_wp)/(tn-3.0_wp)
                                fnpm = fn+fm
                                fnmm = fn-fm
                                temp = fnpm*(fnpm-1.0_wp)
                                cc = sqrt(cn*(fnpm-3.0_wp)*(fnpm-2.0_wp)/temp)
                                dd = sqrt(cn*fnmm*(fnmm-1.0_wp)/temp)
                                ee = sqrt((fnmm+1.0_wp)*(fnmm+2.0_wp)/temp)
                                p3(:, np1) = cc*p1(:, np1-2)+dd*p3(:, np1-2)-ee*p1(:, np1)
                            end do
                        end if

                        p1(:, m:nmp1) = p2(:, m:nmp1)
                        p2(:, m:nmp1) = p3(:, m:nmp1)

                    end if
            end select

        end subroutine lfin_lower_utility_routine

    end subroutine lfin

    subroutine lfp (init, n, m, l, cp, pb, w)
        !
        !
        ! subroutine lfp (init, n, m, l, cp, pb, w)
        !
        ! dimension of           cp((n/2)+1), pb(l), w(5*l+41)
        ! arguments
        !
        ! purpose                routine lfp uses coefficients computed by
        !                        routine alfk to calculate the 64-bit double precision
        !                        normalized associated legendre function pbar(n,
        !                        m, theta) at colatitudes theta=(i-1)*pi/(l-1),
        !                        i=1, ..., l. subroutine lfp evaluates pbar
        !                        using one of the following trigonometric
        !                        expansions
        !
        !                        1) for n even and m even, pbar(m, n, theta) =
        !                           .5*cp(1) plus the sum from k=1 to k=n/2
        !                           of cp(k)*cos(2*k*th)
        !
        !                        2) for n even and m odd, pbar(m, n, theta) =
        !                           the sum from k=1 to k=n/2 of
        !                           cp(k)*sin(2*k*th)
        !
        !                        3) for n odd and m even, pbar(m, n, theta) =
        !                           the sum from k=1 to k=(n+1)/2 of
        !                           cp(k)*cos((2*k-1)*th)
        !
        !                        4) for n odd and m odd,  pbar(m, n, theta) =
        !                           the sum from k=1 to k=(n+1)/2 of
        !                           cp(k)*sin((2*k-1)*th)
        !
        !
        ! usage                  call lfp(init, n, m, l, cp, pb, w)
        !
        ! arguments
        !
        ! on input               init
        !                          = 0 initialization only
        !                          = 1 compute pbar(n, m, theta)
        !
        !                          lfp call with init = 0 initializes array w;
        !                          no values of pbar(n, m, theta) are computed.
        !                          init=0 should be used on the first call, or
        !                          if l or w values differ from those in the
        !                          previous call.
        !
        !                        n
        !                          nonnegative integer, less than l, specifying
        !                          the degree of pbar(n, m, theta)
        !
        !                        m
        !                          is the order of pbar(n, m, theta). m can be
        !                          any integer however pbar(n, m, theta) = 0
        !                          if abs(m) is greater than n and
        !                          pbar(n, m, theta) = (-1)**m*pbar(n, -m, theta)
        !                          for negative m.
        !
        !                        l
        !                          number of colatitudes theta=(i-1)*pi/(l-1)
        !                          for i=1, ..., l where l is greater than 1.
        !                          l must be an odd integer.
        !
        !                        cp
        !                          64-bit double precision array of length (n/2)+1
        !                          containing coefficients computed by routine
        !                          alfk
        !
        !                        w
        !                          a 64-bit double precision work array with at
        !                          least 5*l+41 locations
        !
        ! on output              pb
        !                          64-bit double precision array of length l containing
        !                          pbar(n, m, theta), theta=(i-1)*pi/(l-1) for i=1
        !                          , ..., l.
        !
        !                        w
        !                          a 64-bit double precision array containing values
        !                          which must not be destroyed if the next call
        !                          will have the same value of input parameter n
        !
        ! special conditions     calls to routine lfp must be preceded by an
        !                        appropriate call to routine alfk.
        !
        ! precision              64-bit double precision
        !
        ! algorithm              the trigonometric series formula used by
        !                        routine lfp to calculate pbar(n, m, theta) for
        !                        theta=(i-1)*pi/(l-1), i=1, ..., n, depends on
        !                        m and n as follows:
        !
        !                           1) for n even and m even, the formula is
        !                              .5*cp(1) plus the sum from k=1 to k=n/2
        !                              of cp(k)*cos(2*k*theta)
        !                           2) for n even and m odd. the formula is
        !                              the sum from k=1 to k=n/2 of
        !                              cp(k)*sin(2*k*theta)
        !                           3) for n odd and m even, the formula is
        !                              the sum from k=1 to k=(n+1)/2 of
        !                              cp(k)*cos((2*k-1)*theta)
        !                           4) for n odd and m odd, the formula is
        !                              the sum from k=1 to k=(n+1)/2 of
        !                              cp(k)*sin((2*k-1)*theta)
        !
        ! accuracy               comparison between routines lfp and double
        !                        precision dlfp on the cray1 indicates greater
        !                        accuracy for smaller values of input parameter
        !                        n.  agreement to 12 places was obtained for
        !                        n=10 and to 11 places for n=100.
        !
        ! timing                 time per call to routine lfp is dependent on
        !                        the input parameters l and n.
        !

        ! Dummy arguments

        integer(ip), intent(in)  :: init
        integer(ip), intent(in)  :: n
        integer(ip), intent(in)  :: m
        integer(ip), intent(in)  :: l
        real(wp)                  :: cp(n/2 + 1)
        real(wp)                  :: pb(l)
        real(wp)                  :: w(5*l+41)

        ! Dummy arguments

        integer(ip) :: ma, iw1, iw2

        pb = 0.0_wp

        ma = abs(m)

        if (ma > n) return

        !
        !  Set workspace indices
        !
        iw1 = 2*l+12
        iw2 = iw1+3*(l+1)/2+15

        call lfp_lower_utility_routine(init, n, ma, l, cp, pb, w, w(iw1), w(iw2))

    contains

        subroutine lfp_lower_utility_routine(init, n, m, l, cp, p, wsave1, wsave2, wsave3)

            ! Dummy arguments

            integer(ip), intent(in) :: init
            integer(ip), intent(in) :: n
            integer(ip), intent(in) :: m
            integer(ip), intent(in) :: l
            real(wp)                 :: cp(*)
            real(wp)                 :: p(*)
            real(wp)                 :: wsave1(*)
            real(wp)                 :: wsave2(*)
            real(wp)                 :: wsave3(*)

            ! Local variables
            integer(ip), save   :: lc, lq, ls
            integer(ip)         :: i
            integer(ip)         :: lm1, np1, ls2, kdp, lmi
            real(wp)            :: dt
            real(wp), parameter :: ONE_OVER_SQRT2 = 1.0_wp/sqrt(2.0_wp)
            type(FastFourierTransform)       :: fft

            select case (init)
                case (0)

                    lc=(l+1)/2
                    ls=lc-2
                    lq=lc-1
                    call fft%sinti(ls, wsave1)
                    call fft%costi(lc, wsave2)
                    call fft%cosqi(lq, wsave3)

                case default

                    if ((n <= 0) .and. (m <= 0)) then
                        p(1:l) = ONE_OVER_SQRT2
                    else
                        ls2 = (l+1)/2
                        lm1 = l-1
                        np1 = n+1

                        dt = PI/lm1

                        if (mod(n, 2) <= 0) then
                            if (mod(m, 2) <= 0) then
                                kdp = n/2+1
                                p(1:kdp)=0.5_wp*cp(1:kdp)
                                p(lc)= 2.0_wp*p(lc)

                                call fft%cost(lc, p, wsave2)

                                do i=1, lc
                                    lmi=l-i
                                    p(lmi+1)=p(i)
                                end do
                            else
                                kdp=n/2
                                p(2:kdp+1)=0.5_wp*cp(1:kdp)
                                p(ls+2)=0.0_wp

                                call fft%sint(ls, p(2), wsave1)

                                do i=1, ls
                                    lmi=l-i
                                    p(lmi)=-p(i+1)
                                end do

                                p(l)=0.0_wp
                            end if
                        else
                            kdp=(n+1)/2

                            if (mod(m, 2) <= 0) then

                                p(1:kdp)=0.25_wp*cp(1:kdp)

                                call fft%cosqb(lq, p, wsave3)

                                do i=1, lq
                                    lmi=l-i
                                    p(lmi+1)=-p(i)
                                end do

                            else
                                p(2:kdp+1)=0.25_wp*cp(1:kdp)

                                call fft%sinqb(lq, p(2), wsave3)

                                do i=1, lq
                                    lmi=l-i
                                    p(lmi)=p(i+1)
                                end do

                                p(l)=0.0_wp
                            end if
                        end if
                    end if
            end select

        end subroutine lfp_lower_utility_routine

    end subroutine lfp

    subroutine lfpt(n, m, theta, cp, pb)
        !
        ! subroutine lfpt (n, m, theta, cp, pb)
        !
        ! dimension of
        ! arguments
        !                        cp((n/2)+1)
        !
        ! purpose                routine lfpt uses coefficients computed by
        !                        routine alfk to compute the single precision
        !                        normalized associated legendre function pbar(n,
        !                        m, theta) at colatitude theta.
        !
        ! usage                  call lfpt(n, m, theta, cp, pb)
        !
        ! arguments
        !
        ! on input               n
        !                          nonnegative integer specifying the degree of
        !                          pbar(n, m, theta)
        !                        m
        !                          is the order of pbar(n, m, theta). m can be
        !                          any integer however pbar(n, m, theta) = 0
        !                          if abs(m) is greater than n and
        !                          pbar(n, m, theta) = (-1)**m*pbar(n, -m, theta)
        !                          for negative m.
        !
        !                        theta
        !                          single precision colatitude in radians
        !
        !                        cp
        !                          single precision array of length (n/2)+1
        !                          containing coefficients computed by routine
        !                          alfk
        !
        ! on output              pb
        !                          single precision variable containing
        !                          pbar(n, m, theta)
        !
        ! special conditions     calls to routine lfpt must be preceded by an
        !                        appropriate call to routine alfk.
        !
        ! precision              single
        !
        ! algorithm              the trigonometric series formula used by
        !                        routine lfpt to calculate pbar(n, m, th) at
        !                        colatitude th depends on m and n as follows:
        !
        !                           1) for n even and m even, the formula is
        !                              .5*cp(1) plus the sum from k=1 to k=n/2
        !                              of cp(k)*cos(2*k*th)
        !                           2) for n even and m odd. the formula is
        !                              the sum from k=1 to k=n/2 of
        !                              cp(k)*sin(2*k*th)
        !                           3) for n odd and m even, the formula is
        !                              the sum from k=1 to k=(n+1)/2 of
        !                              cp(k)*cos((2*k-1)*th)
        !                           4) for n odd and m odd, the formula is
        !                              the sum from k=1 to k=(n+1)/2 of
        !                              cp(k)*sin((2*k-1)*th)
        !
        ! accuracy               comparison between routines lfpt and double
        !                        precision dlfpt on the cray1 indicates greater
        !                        accuracy for greater values on input parameter
        !                        n.  agreement to 13 places was obtained for
        !                        n=10 and to 12 places for n=100.
        !
        ! timing                 time per call to routine lfpt is dependent on
        !                        the input parameter n.
        !

        ! Dummy arguments

        integer(ip), intent(in)  :: n
        integer(ip), intent(in)  :: m
        real(wp),    intent(in)  :: theta
        real(wp),    intent(in)  :: cp(n/2+1)
        real(wp),    intent(out) :: pb

        ! Local variables

        integer(ip) :: ma, np1, k, kdo, kp1
        real(wp)    :: cos2t, sin2t, cost, sint, temp, summation

        pb = 0.0_wp
        ma = abs(m)

        if (ma <= n) then
            if (n <= 0) then
                if (ma <= 0) pb= sqrt(0.5_wp)
            else
                np1 = n+1
                if (mod(n, 2) <= 0) then
                    if (mod(ma, 2) <= 0) then
                        kdo = n/2+1
                        cos2t = cos(2.0_wp*theta)
                        sin2t = sin(2.0_wp*theta)
                        cost = 1.0_wp
                        sint = 0.0_wp
                        summation = 0.5_wp * cp(1)

                        do kp1=2, kdo
                            temp = cos2t*cost-sin2t*sint
                            sint = sin2t*cost+cos2t*sint
                            cost = temp
                            summation = summation+cp(kp1)*cost
                        end do

                        pb = summation
                    else
                        kdo = n/2
                        cos2t = cos(2.0_wp*theta)
                        sin2t = sin(2.0_wp*theta)
                        cost = 1.0_wp
                        sint = 0.0_wp
                        summation = 0.0_wp

                        do k=1, kdo
                            temp = cos2t*cost-sin2t*sint
                            sint = sin2t*cost+cos2t*sint
                            cost = temp
                            summation = summation+cp(k)*sint
                        end do

                        pb = summation
                    end if
                else
                    kdo = (n+1)/2
                    if (mod(ma, 2) <= 0) then
                        cos2t = cos(2.0*theta)
                        sin2t = sin(2.0*theta)
                        cost = cos(theta)
                        sint = -sin(theta)
                        summation = 0.0_wp

                        do k=1, kdo
                            temp = cos2t*cost-sin2t*sint
                            sint = sin2t*cost+cos2t*sint
                            cost = temp
                            summation = summation+cp(k)*cost
                        end do

                        pb= summation
                    else
                        cos2t = cos(2.0*theta)
                        sin2t = sin(2.0*theta)
                        cost = cos(theta)
                        sint = -sin(theta)
                        summation = 0.0_wp

                        do k=1, kdo
                            temp = cos2t*cost-sin2t*sint
                            sint = sin2t*cost+cos2t*sint
                            cost = temp
                            summation = summation+cp(k)*sint
                        end do

                        pb = summation
                    end if
                end if
            end if
        end if

    end subroutine lfpt

end module type_AssociatedLegendrePolynomialGenerator


text:
Noise permits allow you to have a temporary exception to the noise ordinance. You can download a noise permit http://www.cityofdavis.org/police/pdfs/forms/Noise_Permit_Application_and_Instructions.pdf from the City website. Permits cost $25 even if you are denied the permit. You can also pick up a noise permit application from the Davis Police Department.  You must turn in your application to the Police Department. The permits allow you to have noise above the noise ordinance until midnight on weekends at 10pm during weekdays. You cannot receive a noise permit if you have been cited with a noise violation. You cannot also receive a noise permit if one has been given to you or a property within 200 yards in the past 3 months. Noise permits are not needed on campus as there is no noise ordinance, although excessive loud noise will probably bring the UC police.

The police do not need to verify the intensity of the sound to cite you with a violation of the noise ordinance, all they need is a complaint. The majority of the time, Community Service Officers are the ones giving out the violations while a the real police stand by as backup.

If you are served 3 or more noise violations it is considered a misdeamnor and you may be arrested.

Receiving a noise violation does not give the police permission to enter your residence unless they cannot immediately find a resident or if a resident or partygoer gives them permission. Howevere, the police may use other reasons for entry, such as seeing underage drinkers or the smell of illegal substances. Police may also enter a residence without permission to break up a fight or to stop a crime in progress.

You do not have to, but it is always smart to be courteous to the police as it may make the difference between a $160+ fine and a warning.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
Emmanuel DiazOrdaz is a current ASUCD ASUCD Senate Senator. He was elected in the Fall 2010 ASUCD Election as an independent.

Support Emmanuel on Facebook http://www.facebook.com/home.php?#!/event.php?eid156028331099935 here

UCD Involvement :
P.E.A.C.E. TrainerCross Cultural Center,
P.E.A.C.E. CoCoordinatorCross Cultural Center
Admin. Director for Yik’al Kuyum Student Recruitment & Retention Center
ASUCD Gender And Sexualities Commissioner, Spring 2009Present
REACE Participant Winter 2010, Participant
Queer Leadership Retreat Spring 2010, Participant
Chican@/Latin@ Leadership Retreat Fall 2010, Workshop Facilitator
Native Leadership Retreat Fall 2010, Workshop Facilitator
Outreach to multiple Yolo/Sacramento County High Schools
Aggies of Color

 Candidate Statement

Hello beautiful Aggies!! My name is Emmanuel DiazOrdaz and I am running for ASUCD Senate!!! (vote for me :)

In case youve never voted in an ASUCD election before (and thats most of the student body) heres a little bit of infothere are twelve seats in the ASUCD Senate, and six of them are up for election this quarter! Because ASUCD elections use choice voting, where you RANK the candidates instead of just picking one, just ...remember to vote Emmanuel DiazOrdaz #1 and Cameron Brown #2!!

Campaign Platform
AB540 AWARENESS WEEK
A weeklong event dedicated to the issues around AB540, and undocumented students ending with a scholarship giveaway open to all students.

TACO TRUCK ON CAMPUS
In collaboration with the Retention Coordinator for Yik’al Kuyum at the Student Recruitment & Retention Center, I want to bring a Taco Truck to campus at least twice a month during La Raza Tuesdays.

GENERATION SEX WEEK PROGRAMING
Gender Soliloqueers
The second annual Gender Soliloqueers invites students of all genders and sexualities to express themselves through performances, spoken word, and radial monologues.
Foreskin Awareness Day
A day dedicated to foreskin awareness specifically it’s effects on men physically, mentally, emotionally, and creatively.
Rape Culture In Greek Life
A program advocating for victims of violence and rape, specifically within Greek communities.

INCREASE LEADERSHIP RETREAT/CONFERENCE ALLOCATION
As of now, some money is being allocated to certain leadership retreats on campus. However, it is not nearly enough. There needs to be more institutionalized money for these very crucial community building, ally developing retreats and conferences in order for us to learn from each other and other UC’s to make our own campus a safer space for all communities.

Short and Simple, mostly.
If you have any questions, ask me :


text:
      subroutine view(x,y,z)
      common/gfmatr/t(32),nt,nt1
       common/gfview/mac(20),xp1,yp1,zp1
      nw=nt
      nt1=3
      xp1=x
      yp1=y
      zp1=z
      if(nt.ne.1)
     1call hcinv(x,y,z,xp1,yp1,zp1)
      call hcrot1(x,y,z)
       r=sqrt(x**2+y**2+z**2)
      call hcprsp(r)
      nt=nw
      return
       end

max_stars_repo_name: sergev/grafor

max_issues_repo_name: sergev/grafor

max_forks_repo_name: sergev/grafor


text:
C    Copyright (c) 2014, Sandia Corporation.
C    Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
C    the U.S. Government retains certain rights in this software.
C
C    Redistribution and use in source and binary forms, with or without
C    modification, are permitted provided that the following conditions are
C    met:
C
C        * Redistributions of source code must retain the above copyright
C          notice, this list of conditions and the following disclaimer.
C
C        * Redistributions in binary form must reproduce the above
C          copyright notice, this list of conditions and the following
C          disclaimer in the documentation and/or other materials provided
C          with the distribution.
C
C        * Neither the name of Sandia Corporation nor the names of its
C          contributors may be used to endorse or promote products derived
C          from this software without specific prior written permission.
C
C    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C    A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C    OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C    THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
C

C $Id: inhole.f,v 1.1 1990/11/30 11:09:43 gdsjaar Exp $
C $Log: inhole.f,v $
C Revision 1.1  1990/11/30 11:09:43  gdsjaar
C Initial revision
C
C
CC* FILE: [.MAIN]INHOLE.FOR
CC* MODIFIED BY: TED BLACKER
CC* MODIFICATION DATE: 7/6/90
CC* MODIFICATION: COMPLETED HEADER INFORMATION
C
      SUBROUTINE INHOLE (MR, N7, N29, JPNTR, IIN, IFOUND, IFHOLE, NHPR,
     &   IHLIST, MERGE, NOROOM)
C***********************************************************************
C
C  SUBROUTINE INHOLE  =  INPUTS A REGION'S HOLES INTO THE DATABASE
C
C***********************************************************************
C
      DIMENSION IFHOLE(MR), NHPR(MR), IHLIST(MR*2)
      DIMENSION IIN(IFOUND)
C
      LOGICAL NOROOM, MERGE
C
      NOROOM = .TRUE.
C
C  ADD THE REGION INTO THE DATABASE
C
      J = JPNTR
      IFHOLE(J) = N29 + 1
      DO 100 I = 1, IFOUND
         JJ = IIN(I)
         IF (JJ .EQ. 0) GO TO 110
         N29 = N29 + 1
         IF (N29 .GT. MR*2) RETURN
         IHLIST(N29) = JJ
  100 CONTINUE
C
  110 CONTINUE
      NHPR(J) = N29 - IFHOLE(J) + 1
      IF (NHPR(J) .LT. 1) THEN
         WRITE(*, 10000) J
         N29 = IFHOLE(J) - 1
      END IF
C
      NOROOM = .FALSE.
      RETURN
C
10000 FORMAT(' REGION:', I5, ' HAS LESS THAN ONE HOLE', /,
     &   ' THE HOLES FOR THIS REGION WILL NOT BE INPUT INTO DATABASE')
C
      END

max_stars_repo_name: milljm/seacas

max_issues_repo_name: milljm/seacas

max_forks_repo_name: milljm/seacas


text:
C
C $Id: aggtch.f,v 1.6 2008-07-27 00:14:34 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE AGGTCH (IDCS,CHST,LNCS)
C
      CHARACTER*(*) CHST
C
C This routine gets character strings previously stored by the routine
C AGSTCH (which see).  It has the following arguments:
C
C -- IDCS is the identifying integer returned by AGSTCH when the string
C    was stored.
C
C -- CHST is the character string returned.
C
C -- LNCS is the length of the character string returned in CHST.
C
C The following common blocks contain variables which are required for
C the character-storage-and-retrieval scheme of AUTOGRAPH.
C
      COMMON /AGCHR1/ LNIC,INCH(2,50),LNCA,INCA
      SAVE /AGCHR1/
C
      COMMON /AGCHR2/ CHRA(2000)
      CHARACTER*1 CHRA
      SAVE /AGCHR2/
C
C First, blank-fill the character variable to be returned.
C
      CHST=' '
C
C If the identifier is less than -LNIC, the (one-character) string is
C retrieved from it.
C
      IF (IDCS.LT.(-LNIC)) THEN
        CHST=CHAR(-IDCS-LNIC-1)
        LNCS=1
C
C If the identifier is between -LNIC and -1, its absolute value is the
C index, in INCH, of the descriptor of the character string stored in
C CHRA.
C
      ELSE IF (IDCS.LE.(-1)) THEN
        I=-IDCS
        J=INCH(1,I)-1
        IF (J.GE.0) THEN
          LNCS=MIN(LEN(CHST),INCH(2,I))
          DO 101 K=1,LNCS
            J=J+1
            CHST(K:K)=CHRA(J)
  101     CONTINUE
        ELSE
          LNCS=0
        END IF
C
C In all other cases, return a single blank.
C
      ELSE
        LNCS=1
C
      END IF
C
C Done.
C
      RETURN
C
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


! RUN: %S/test_errors.sh %s %t %flang -fopenmp
! REQUIRES: shell
! OpenMP Version 4.5
! 2.15.3.3 private Clause
! Variables that appear in namelist statements may not appear in a private clause.

module test
  integer :: a, b, c
  namelist /nlist1/ a, b
end module

program omp_private
  use test

  integer :: p(10) ,q(10)
  namelist /nlist2/ c, d

  a = 5
  b = 10
  c = 100

  !ERROR: Variable 'a' in NAMELIST cannot be in a PRIVATE clause
  !ERROR: Variable 'c' in NAMELIST cannot be in a PRIVATE clause
  !$omp parallel private(a, c)
  d = a + b
  !$omp end parallel

  call sb()

  contains
    subroutine sb()
      namelist /nlist3/ p, q

      !ERROR: Variable 'p' in NAMELIST cannot be in a PRIVATE clause
      !ERROR: Variable 'd' in NAMELIST cannot be in a PRIVATE clause
      !$omp parallel private(p, d)
      p = c * b
      q = p * d
      !$omp end parallel

      write(*, nlist1)
      write(*, nlist2)
      write(*, nlist3)

    end subroutine

end program omp_private


text:
In order to make driving in Davis easier for you it will be helpful to know about the local traffic patterns. Or avoid the whole thing by Telecommute telecommuting. The City of Davis also has a program called Street Smarts that provides information on using transportation in town.

 Bicycle
 on the UCD campus

Busy times include 9AM and lunchtime.  Watch where youre going, and look out for people who arent paying attention!  Pedestrians should be courteous and walk on sidewalks rather than in the bike lanes.  If youre forced to walk in a bike lane, please try to keep to the rightmost side.  Such common courtesy is especially important when there are a lot of folks on the move.

 in town

Swarms of kids on bikes converge on the various Schools around town, especially the http://daviswiki.org/Schools#head0842148f0ca04c936405774f0e6ee4a76998e459 Junior High Schools.  The busiest times are 7AM8AM and ~3PM.

 Auto
 in town

The numerous seemingly untimed stoplights on 5th Street, which turns into Russell Boulevard on the West side of town, contributes to what is probably the worst street traffic in Davis during commute hours.  If youre trying to get across town during 8AM9PM and/or 5PM6PM, you might be better off using Covell Boulevard or Interstate 80 instead.

 Interstate 80

The most predictable traffic jams in the Davis area on Interstate 80 I80:

  7AM9AM Eastbound (folks going to work in Sacramento)
  7AM9AM Westbound near Vacaville (folks going to work in the San Francisco SF Bay Area)
  Fridays 3PM7PM Eastbound (folks going to Tahoe, etc. for the weekend)

 Highway 113

This short stretch of Highway 113 highway rarely sees any traffic jams, except when theres a wreck.


text:
C
C
*+ image_print

       subroutine image_print(minirt,data,status)
C      ------------------------------------------
C
C Print values from an image
C
C Given:
C   redtape items defining image
       integer    minirt(*)
C   image data array
       real*4     data(*)
C Returned:
C   status variable
       integer    status
C
C To produce a formatted printed output of a region of the specified image.
C The user will be requested to define the format of the output printing.
C
C Output will be sent to the user defined output device.
C All user IO is handled using the IOLIB routine library.
C
C STATUS should be zero on entry to this routine and is not changed
C explicitly.
C
C The routine is NON-STANDARD a call is made to PLT_GETUV routine.
C
C-
C Version 1.0 10/10/89  [P. Alexander, MRAO]
C

C loop counters
       integer       i, ii, n, iu, iv, iuv(2), iu1, iu2
C local UV-range of map
       integer       uv_range(4)
C local image statistics
       real*4        image_stats_r(4), image_max
       integer       image_stats_i(4)
C parameters controlling type of output
       logical       real_output
       real*4        scale_factor, value(50)
       integer       field_size
       character     output_format*80, title_format*80, title*80
       character*1   dash(132)
       integer       len_format, ivalue(50)
       integer       npanel, ipanel, width
C output unit for terminal IO
       integer       iout
C output unit for printed output
       integer       iprint

C Functions
       include '/mrao/include/iolib_functions.inc'
       include '/mrao/include/chrlib_functions.inc'
       logical   cmd_dblev

C check status on entry
       if (status.ne.0) return

C initialise character variables
       iprint = 0
       title_format = ' '
       output_format= ' '
       do i=1,132
         dash(i) = '-'
       end do

C find the range of image to display and max/min
       do i = 1,4
         uv_range(i) = minirt(i)
       end do
       call plot_getuv('UV-range : ','*',uv_range,status)
       call scnmap2(minirt,data,uv_range,image_stats_r,image_stats_i,
     *              status)
       image_max = max(abs(image_stats_r(1)),abs(image_stats_r(2)))
       if (status.ne.0) goto 999

C determine type of output
       call io_enqout(iout)
       write(iout,1) image_max
  1    format(1x,'Absolute maximum in range = ',1PE11.3)
       real_output = io_yesno(
     *      'Print values in floating point format ? ','no',status)
       if (status.ne.0) goto 999
       if (real_output) then
         call io_geti('Number of decimal places : ','2',
     *                field_size,status)
       else
         call io_geti('Size of each output field : ','8',
     *                field_size,status)
         scale_factor = 10.0**(nint(field_size - 2 - log10(image_max)))
         call io_getr('Scale-factor : ','*',scale_factor,status)
       end if

C open the output unit for image display
2      call io_geti('Output width (40<=width<=132) : ','80',
     *              width,status)
       if (status.ne.0) goto 999
       if (width.lt.40 .or. width.gt.132) then
         call cmd_wrerr('PRINT-IMAGE','Illegal device width')
         goto 2
       end if
       call io_getstr('Title : ',' ',title,status)
       call io_opeout(iprint,status)
       write(iprint,'(1X/1X,''PRINT-IMAGE Title: '',a)')
     *       title(1:chr_lenb(title))
       if (status.ne.0) goto 999

C loop through the output data printing the values in the requested
C format
       if (real_output) then

C real output
         n = (width - 6)/(field_size + 8)
         if (n.gt.50) n = 50
         npanel = (uv_range(2)-uv_range(1)-1)/n + 1
         if (iprint.ne.iout) then
           write(iprint,'(a,i2)')'Image displayed as ',npanel,' panels'
         end if
         print *,'.. Image displayed as ',npanel,' panels'
         if (field_size.ge.2) then
           write(output_format,3) 8+field_size,field_size
   3       format('(1X,I5,'':'',50(1PE',i2,'.',i1,'))')
           write(title_format,31) 8+field_size
   31      format('(1X/1X/1X,''Panel number: '',i3/',
     *            '1X,'' V/U :'',50I',i2,')')
         else
           write(output_format,4) 8+field_size,field_size
   4       format('(1X,I5,'':'',50(1PE',i1,'.',i1,'))')
           write(title_format,41) 8+field_size
   41      format('(1X/1X/1X,''Panel number: '',i3/',
     *            '1X,'' V/U :'',50I',i1,')')
         end if
         if (cmd_dblev(5)) then
           print *,output_format(1:chr_lenb(output_format))
           print *,title_format(1:chr_lenb(title_format))
         end if
         do ipanel = 1,npanel
           iu1 = uv_range(1) + (ipanel-1)*n
           iu2 = min(uv_range(2),iu1+n-1)
           ii = 0
           if (cmd_dblev(2)) then
             print *,'.. panel: ',ipanel
             print *,'.. range: ',iu1,iu2,uv_range(3), uv_range(4)
           end if
           do iu=iu1,iu2
             ii = ii + 1
             ivalue(ii) = iu
           end do
           len_format = chr_lenb(title_format)
           write(iprint,fmt=title_format(1:len_format))
     *           ipanel,(ivalue(i), i=1,ii)
           write(iprint,'(1X,132A)') (dash(i), i=1,width)
           do iv=uv_range(3),uv_range(4),-1
             ii = 0
             do iu=iu1,iu2
               ii = ii + 1
               iuv(1) = iu
               iuv(2) = iv
               call iuvval2(minirt,data,iuv,value(ii),status)
             end do
             len_format = chr_lenb(output_format)
             write(iprint,fmt=output_format(1:len_format))
     *             iv, (value(i), i=1,ii)
             if (status.ne.0) goto 999
           end do
           write(iprint,'(1X,132A)') (dash(i), i=1,width)
         end do
       else

C integer output
         n = (width - 6)/(field_size + 1)
         npanel = (uv_range(2)-uv_range(1)-1)/n + 1
         if (iprint.ne.iout) then
           write(iprint,'(a,i2)')'Image displayed as ',npanel,' panels'
         end if
         print *,'.. Image displayed as ',npanel,' panels'
         if ((field_size+1).lt.10) then
           write(output_format,5) field_size+1
   5       format('(1X,I5,'':'',50I',i1,')')
           write(title_format,51) field_size+1
   51      format(
     *      '(1X/1X/1X,''Panel number :'',i3,''  Scale: '',1PE12.3/',
     *      '1X,'' V/U :'',50I',i1,')')
         else
           write(output_format,6) field_size+1
   6       format('(1X,I5,'':'',50I',i2,')')
           write(title_format,61) field_size+1
   61      format(
     *      '(1X/1X/1X,''Panel number :'',i3,''  Scale: '',1PE12.3/',
     *      '1X,'' V/U :'',50I',i2,')')
         end if
         if (cmd_dblev(5)) then
           print *,output_format(1:chr_lenb(output_format))
           print *,title_format(1:chr_lenb(title_format))
         end if
         do ipanel = 1,npanel
           len_format = chr_lenb(title_format)
           iu1 = uv_range(1) + (ipanel-1)*n
           iu2 = min(uv_range(2),iu1+n-1)
           ii = 0
           do iu=iu1,iu2
             ii = ii + 1
             ivalue(ii) = iu
           end do
           if (cmd_dblev(2)) then
             print *,'.. panel: ',ipanel, n
             print *,'.. range: ',iu1,iu2,uv_range(3), uv_range(4)
           end if
           len_format = chr_lenb(title_format)
           write(iprint,fmt=title_format(1:len_format))
     *           ipanel,scale_factor,(ivalue(i), i=1,ii)
           write(iprint,'(1X,132A)') (dash(i), i=1,width)
           do iv=uv_range(3),uv_range(4),-1
             ii = 0
             do iu=iu1,iu2
               ii = ii + 1
               iuv(1) = iu
               iuv(2) = iv
               call iuvval2(minirt,data,iuv,value(ii),status)
             end do
             len_format = chr_lenb(output_format)
             write(iprint,fmt=output_format(1:len_format))
     *             iv, (nint(value(i)*scale_factor), i=1,ii)
             if (status.ne.0) goto 999
           end do
           write(iprint,'(1X,132A)') (dash(i), i=1,width)
         end do

       end if

999    if ((iprint.ne.iout).and.iprint.ne.0) then
         close (iprint)
         call io_setout(iout)
       end if
       call cmd_err(status,'IMAGE_PRINT','Failed')
       end

max_stars_repo_name: CavendishAstrophysics/anmap

max_issues_repo_name: CavendishAstrophysics/anmap

max_forks_repo_name: CavendishAstrophysics/anmap


!
!  Copyright 2019-2020 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 parallelization
  implicit none

  !!! multi process
  integer, public :: nproc_group_global
  integer, public :: nproc_id_global
  integer, public :: nproc_size_global

  !!! elapsed-time gap between compute node.
  real(8), private :: gap_time_system

  ! call once
  public :: setup_parallel
  public :: end_parallel

  ! util
  public :: adjust_elapse_time
  public :: get_thread_id
  public :: get_nthreads
  public :: is_distributed_parallel

  private :: find_elapse_time_gap

contains
  subroutine setup_parallel
    use communication
    implicit none
    call comm_init
    call comm_get_globalinfo(nproc_group_global, nproc_id_global, nproc_size_global)
    call find_elapse_time_gap
  end subroutine

  subroutine end_parallel
    use communication
    implicit none
    call comm_finalize
  end subroutine

  subroutine find_elapse_time_gap
    use communication
    use misc_routines, only: get_wtime
    implicit none
    real(8) :: now, me
    call comm_sync_all
    now = get_wtime()
    me  = now
    call comm_sync_all
    call comm_get_min(now, nproc_group_global)
    gap_time_system = me - now
  end subroutine

  function adjust_elapse_time(time)
    implicit none
    real(8), intent(in) :: time
    real(8)             :: adjust_elapse_time
    adjust_elapse_time = time - gap_time_system
  end function

  function get_thread_id() result(nid)
#ifdef _OPENMP
    use omp_lib
#endif
    implicit none
    integer :: nid
#ifdef _OPENMP
    nid = omp_get_thread_num()
#else
    nid = 0
#endif
  end function

  function get_nthreads() result(nsize)
#ifdef _OPENMP
    use omp_lib
#endif
    implicit none
    integer :: nsize
#ifdef _OPENMP
    nsize = omp_get_max_threads()
#else
    nsize = 1
#endif
  end function

  function is_distributed_parallel() result(ret)
    implicit none
    logical :: ret
    ret = (nproc_size_global > 1)
  end function
end module


text:
	program matriz_un_medio

	implicit none

	real,allocatable,dimension(:,:)::mat
	integer:: fila,colum,fila2,colum2,i

	Write(*,*) "Ingrese las dimenciones de la matriz"
	Write(*,*) "Filas"
	read(*,*) fila
	Write(*,*) "Columnas"
	read(*,*) colum

	allocate(mat(fila,colum))

	write(*,*) "Ingrese los elementos de la matriz"

	 fila2=1
	 colum2=1

	 do while(fila2<=fila)

		do while(fila2<=fila.and.colum2<=colum)

	Write(*,*)"Ingrese el elemento de la matriz",fila2,",",colum2

		read(*,*)mat(fila2,colum2)
		colum2=colum2+1
		enddo

		fila2=fila2+1
		colum2=1

	  enddo

	Write(*,*)"La matriz ingresada fue :"

	do i=1,fila,1

	write(*,*) mat(i,:)
	enddo

!-----------------------------------------
	  deallocate(mat)
	  pause

	  end program

max_stars_repo_name: patricio-c/fortran

max_issues_repo_name: patricio-c/fortran

max_forks_repo_name: patricio-c/fortran


text:
! This Source Code Form is subject to the terms of the Mozilla Public
! License, v. 2.0. If a copy of the MPL was not distributed with this
! file, You can obtain one at http://mozilla.org/MPL/2.0/.
!> \brief The Molecular Orbital module. Handles exclusively the MOs
!> \author David Krisiloff
!> \date 9/1/2011
!>
!> Designed in the hopes of making the code just slightly more modular.
!> This should reduced the number of copies of the MOs (we need just 1 !!)
!>

! A couple key points:
! (1) This is in no way compatabile with anything PAO ... mostly because of the headache
!     that rewriting so much of the code ( which I'm not familiar with ) would be. I'm
!     terribly sorry if this causes someone else a headache
!
! (2) Read but DO NOT WRITE TO molecular_orbitals. Unless you really know what you are
!     doing. If you randomly alter the molecular_orbitals, weird things are going to
!     happen.

module molecular_orbital_mod
    use fortran_timing
    use global_var_mod
    use molecule_var_mod
    use utilities_mod
    use get_integrals
    use math_utils

    implicit none

    real(kind=real8), allocatable, dimension(:,:) :: molecular_orbitals

    contains

    !> \brief Does all the molecular orbital setup, reading them off this and possibly
    !> calculating the LOVOs
    subroutine get_molecular_orbitals()
        implicit none

        ! local variables
        integer :: nOrb                         ! number of molecular orbitals
        integer :: nBas                         ! number of atom centered basis functions
        integer :: ierr
        integer :: iCoeff                                   ! a pointer for MOLCAS

        ! Initialize just for safety
        ierr = 0
        iCoeff = 0

        write(ioOutput,130)
    130 format(/,1x,"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!",/&
                ,1x,"!//                              ",/&
                ,1x,"!// ENTERED MOLECULAR ORBITALS   ",/&
                ,1X,"!//                              ",/&
                ,1x,"!//                              ",/&
                ,1x,"!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!",/)

        !> \todo We should have as a user input the name of the orbital file
        nOrb = num_orbitals
        nBas = number_bas

        ! Allocate space for the molecular_orbitals. Notice the typical FORTRAN column ordering.
        ! Each column is a separate orbital
        allocate( molecular_orbitals(nBas, nOrb), stat=ierr )
        call allocatecheck(ierr,'molecular_orbitals')
        molecular_orbitals = 0
        ! Read the MOs from disk (doesn't matter if it's local or nonlocal anymore)
        call get_coefficients(nBas, nOrb, molecular_orbitals)

        ! This next part shouldn't be necessary. However because ever once in a while
        ! something bad happens to (i.e. I screw up) the orbitals file ... lets just double check that
        ! everything is normalized.
        call check_norms(nOrb, nBas)

    end subroutine get_molecular_orbitals

    !> \brief Double checks that all molecular orbitals are normalized to 1
    subroutine check_norms(nOrb, nBas)
        implicit none

        ! input variables
        integer, intent(in) :: nOrb                         ! number of molecular orbitals
        integer, intent(in) :: nBas                         ! number of atom centered basis functions

        ! local variables
        real(kind=real8), allocatable, dimension(:,:) :: AOoverlap
        real(kind=real8) :: norm
        integer :: iOrb, jOrb, ierr
        logical :: error

        ! First allocate and obtain the AOoverlap matrix
        allocate( AOoverlap(nBas, nBas), stat=ierr)
        call allocatecheck(ierr, "AOoverlap")
        call get_overlap(nBas, AOoverlap)

        ! Go through each orbital and check the norm
        error = .false.
        write(*,*) "Checking orbital norms ..."
        do iOrb = 1 , nOrb
            norm = calc_norm(molecular_orbitals(:,iOrb), molecular_orbitals(:,iOrb), AOoverlap, nBas )
            if ( abs(norm -1.0) > 5.0D-6 ) then
               100 format("Orbital ", i5, " has an unusual norm " , E20.10)
               write(*,100) iOrb, norm
               error = .true.
            end if
            do jOrb = iOrb + 1, nOrb
               norm = calc_norm(molecular_orbitals(:,iOrb), molecular_orbitals(:,jOrb), AOoverlap, nBas )
               if ( abs(norm) > 5.0D-6 ) then
101               format("Orbitals ", i5, " and ", i5, " are not orthogonal, with overlap " , E20.10)
                  write(*,101) iOrb, jOrb, norm
                  error = .true.
               end if
            end do
        end do
        write(*,*)

        if (error) then
           write(*,*) "******************************************"
           write(*,*) "                                          "
           write(*,*) "SERIOUS WARNING FROM MOLECULAR_ORBITAL_MOD"
           write(*,*) "NORMALIZATION ISSUES !!!! "
           write(*,*) "                                          "
           write(*,*) "******************************************"
           stop
        else
           write(*,*) "Molecular orbitals are normalized and orthogonal!"
        end if

    end subroutine check_norms

    !> \brief Calculates the norm of a molecular orbital, (AO basis is not orthogonal)
    real function calc_norm (vec1, vec2, overlap, nBas )
        implicit none

        ! input variables
        real(kind=real8), dimension(:), intent(in)   :: vec1, vec2
        real(kind=real8), dimension(:,:), intent(in) :: overlap
        integer, intent(in) :: nBas

        ! local variables
        real(kind=real8) :: tmp
        integer :: i,j

        tmp = 0.0
        do i = 1, nBas
            do j = 1, nBas
                tmp = tmp + vec1(i) * vec2(j) * overlap(i,j)
            end do
        end do
        calc_norm = tmp

    end function calc_norm

!#ifndef TIGER_GAMESS
!    !> \brief Writes out the orbital to disk in MOLCAS' format (useful for visualization)
!    !> \param nBas Number of AO orbitals
!    !> \param nOrb Number of MO orbitals to write to disk
!    subroutine write_orbs_MOLCAS(nBas,nOrb)
!        implicit none
!
!        ! inputs
!        integer, intent(in) :: nBas, nOrb
!
!        ! local variables
!        character(len=8) :: fName, Title
!        integer :: unit_num, nSym, ierr
!        integer, dimension(:), allocatable :: indt
!        real(kind=real8), dimension(:), allocatable :: eOrb, Occ
!
!        ! allocate some memory
!        allocate(eOrb(nOrb), Occ(nOrb), indt(nOrb), stat=ierr)
!        call allocatecheck(ierr, "in write_orbs_MOLCAS")
!        eOrb = 0
!        Occ  = 0
!        indt = 0
!
!        !call MOLCAS routine to write the orbitals
!        fName = 'LOVOs'
!        Title = 'debug'
!        unit_num = 768
!        call WrVec( fName, unit_num, 'C', nSym, nBas, nOrb, molecular_orbitals, Occ, eOrb, indt, Title)
!
!        deallocate(Occ, eOrb, indt)
!
!    end subroutine write_orbs_MOLCAs
!#endif

end module molecular_orbital_mod

max_stars_repo_name: EACcodes/TigerCI

max_issues_repo_name: EACcodes/TigerCI

max_forks_repo_name: EACcodes/TigerCI


text:
!
!     CalculiX - A 3-dimensional finite element program
!     Copyright (C) 1998-2018 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     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, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
!
!     subroutine to find the right node for different element types
!     based on face number jface and node number ii
!
      integer function getnodel(ii,jface,nope)
!
!     author: Saskia Sitzmann
!
      implicit none
!
      integer ii,jface,nope,
     &        ifaceq(8,6),ifacet(6,4),ifacew1(4,5),ifacew2(8,5)
!
      include "gauss.f"
!
      data ifaceq /4,3,2,1,11,10,9,12,
     &            5,6,7,8,13,14,15,16,
     &            1,2,6,5,9,18,13,17,
     &            2,3,7,6,10,19,14,18,
     &            3,4,8,7,11,20,15,19,
     &            4,1,5,8,12,17,16,20/
      data ifacet /1,3,2,7,6,5,
     &             1,2,4,5,9,8,
     &             2,3,4,6,10,9,
     &             1,4,3,8,10,7/
!
      data ifacew1 /1,3,2,0,
     &             4,5,6,0,
     &             1,2,5,4,
     &             2,3,6,5,
     &             3,1,4,6/
!
!     nodes per face for quadratic wedge elements
!
      data ifacew2 /1,3,2,9,8,7,0,0,
     &             4,5,6,10,11,12,0,0,
     &             1,2,5,4,7,14,10,13,
     &             2,3,6,5,8,15,11,14,
     &             3,1,4,6,9,13,12,15/
!
      getnodel=0
!
      if((nope.eq.20).or.(nope.eq.8)) then
               getnodel=ifaceq(ii,jface)
      elseif((nope.eq.10).or.(nope.eq.4)) then
         getnodel=ifacet(ii,jface)
      elseif(nope.eq.6) then
         getnodel=ifacew1(ii,jface)
      else
         getnodel=ifacew2(ii,jface)
      endif
!
      end
!

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop