!
!     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!     *                                                               *
!     *                  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                  *
!     *                                                               *
!     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!
submodule(vector_synthesis_routines) vector_synthesis_gaussian_grid

contains
    !     subroutine vhsgc(nlat, nlon, ityp, nt, v, w, idvw, jdvw, br, bi, cr, ci, &
    !                      mdab, ndab, wvhsgc, ierror)
    !
    !     subroutine vhsgc performs the vector spherical harmonic synthesis
    !     of the arrays br, bi, cr, and ci and stores the result in the
    !     arrays v and w. v(i, j) and w(i, j) are the colatitudinal
    !     (measured from the north pole) and east longitudinal components
    !     respectively, located at the gaussian colatitude point theta(i)
    !     and longitude phi(j) = (j-1)*2*pi/nlon. the spectral
    !     representation of (v, w) is given below at output parameters v, w.
    !
    !     input parameters
    !
    !     nlat   the number of points in the gaussian colatitude grid on the
    !            full sphere. these lie in the interval (0, pi) and are computed
    !            in radians in theta(1) <...< theta(nlat) by subroutine compute_gaussian_latitudes_and_weights.
    !            if nlat is odd the equator will be included as the grid point
    !            theta((nlat + 1)/2).  if nlat is even the equator will be
    !            excluded as a grid point and will lie half way between
    !            theta(nlat/2) and theta(nlat/2+1). nlat must be at least 3.
    !            note: on the half sphere, the number of grid points in the
    !            colatitudinal direction is nlat/2 if nlat is even or
    !            (nlat + 1)/2 if nlat is odd.
    !
    !     nlon   the number of distinct londitude points.  nlon determines
    !            the grid increment in longitude as 2*pi/nlon. for example
    !            nlon = 72 for a five degree grid. nlon must be greater
    !            than zero. the axisymmetric case corresponds to nlon=1.
    !            the efficiency of the computation is improved when nlon
    !            is a product of small prime numbers.
    !
    !     ityp   = 0  no symmetries exist about the equator. the synthesis
    !                 is performed on the entire sphere.  i.e. on the
    !                 arrays v(i, j), w(i, j) for i=1, ..., nlat and
    !                 j=1, ..., nlon.
    !
    !            = 1  no symmetries exist about the equator. the synthesis
    !                 is performed on the entire sphere.  i.e. on the
    !                 arrays v(i, j), w(i, j) for i=1, ..., nlat and
    !                 j=1, ..., nlon. the curl of (v, w) is zero. that is,
    !                 (d/dtheta (sin(theta) w) - dv/dphi)/sin(theta) = 0.
    !                 the coefficients cr and ci are zero.
    !
    !            = 2  no symmetries exist about the equator. the synthesis
    !                 is performed on the entire sphere.  i.e. on the
    !                 arrays v(i, j), w(i, j) for i=1, ..., nlat and
    !                 j=1, ..., nlon. the divergence of (v, w) is zero. i.e.,
    !                 (d/dtheta (sin(theta) v) + dw/dphi)/sin(theta) = 0.
    !                 the coefficients br and bi are zero.
    !
    !            = 3  v is symmetric and w is antisymmetric about the
    !                 equator. the synthesis is performed on the northern
    !                 hemisphere only.  i.e., if nlat is odd the synthesis
    !                 is performed on the arrays v(i, j), w(i, j) for
    !                 i=1, ..., (nlat + 1)/2 and j=1, ..., nlon. if nlat is
    !                 even the synthesis is performed on the the arrays
    !                 v(i, j), w(i, j) for i=1, ..., nlat/2 and j=1, ..., nlon.
    !
    !            = 4  v is symmetric and w is antisymmetric about the
    !                 equator. the synthesis is performed on the northern
    !                 hemisphere only.  i.e., if nlat is odd the synthesis
    !                 is performed on the arrays v(i, j), w(i, j) for
    !                 i=1, ..., (nlat + 1)/2 and j=1, ..., nlon. if nlat is
    !                 even the synthesis is performed on the the arrays
    !                 v(i, j), w(i, j) for i=1, ..., nlat/2 and j=1, ..., nlon.
    !                 the curl of (v, w) is zero. that is,
    !                 (d/dtheta (sin(theta) w) - dv/dphi)/sin(theta) = 0.
    !                 the coefficients cr and ci are zero.
    !
    !            = 5  v is symmetric and w is antisymmetric about the
    !                 equator. the synthesis is performed on the northern
    !                 hemisphere only.  i.e., if nlat is odd the synthesis
    !                 is performed on the arrays v(i, j), w(i, j) for
    !                 i=1, ..., (nlat + 1)/2 and j=1, ..., nlon. if nlat is
    !                 even the synthesis is performed on the the arrays
    !                 v(i, j), w(i, j) for i=1, ..., nlat/2 and j=1, ..., nlon.
    !                 the divergence of (v, w) is zero. i.e.,
    !                 (d/dtheta (sin(theta) v) + dw/dphi)/sin(theta) = 0.
    !                 the coefficients br and bi are zero.
    !
    !            = 6  v is antisymmetric and w is symmetric about the
    !                 equator. the synthesis is performed on the northern
    !                 hemisphere only.  i.e., if nlat is odd the synthesis
    !                 is performed on the arrays v(i, j), w(i, j) for
    !                 i=1, ..., (nlat + 1)/2 and j=1, ..., nlon. if nlat is
    !                 even the synthesis is performed on the the arrays
    !                 v(i, j), w(i, j) for i=1, ..., nlat/2 and j=1, ..., nlon.
    !
    !            = 7  v is antisymmetric and w is symmetric about the
    !                 equator. the synthesis is performed on the northern
    !                 hemisphere only.  i.e., if nlat is odd the synthesis
    !                 is performed on the arrays v(i, j), w(i, j) for
    !                 i=1, ..., (nlat + 1)/2 and j=1, ..., nlon. if nlat is
    !                 even the synthesis is performed on the the arrays
    !                 v(i, j), w(i, j) for i=1, ..., nlat/2 and j=1, ..., nlon.
    !                 the curl of (v, w) is zero. that is,
    !                 (d/dtheta (sin(theta) w) - dv/dphi)/sin(theta) = 0.
    !                 the coefficients cr and ci are zero.
    !
    !            = 8  v is antisymmetric and w is symmetric about the
    !                 equator. the synthesis is performed on the northern
    !                 hemisphere only.  i.e., if nlat is odd the synthesis
    !                 is performed on the arrays v(i, j), w(i, j) for
    !                 i=1, ..., (nlat + 1)/2 and j=1, ..., nlon. if nlat is
    !                 even the synthesis is performed on the the arrays
    !                 v(i, j), w(i, j) for i=1, ..., nlat/2 and j=1, ..., nlon.
    !                 the divergence of (v, w) is zero. i.e.,
    !                 (d/dtheta (sin(theta) v) + dw/dphi)/sin(theta) = 0.
    !                 the coefficients br and bi are zero.
    !
    !
    !     nt     the number of syntheses.  in the program that calls vhsgc,
    !            the arrays v, w, br, bi, cr, and ci can be three dimensional
    !            in which case multiple syntheses will be performed.
    !            the third index is the synthesis index which assumes the
    !            values k=1, ..., nt.  for a single synthesis set nt=1. the
    !            discription of the remaining parameters is simplified
    !            by assuming that nt=1 or that all the arrays are two
    !            dimensional.
    !
    !     idvw   the first dimension of the arrays v, w as it appears in
    !            the program that calls vhsgc. if ityp <= 2 then idvw
    !            must be at least nlat.  if ityp > 2 and nlat is
    !            even then idvw must be at least nlat/2. if ityp > 2
    !            and nlat is odd then idvw must be at least (nlat + 1)/2.
    !
    !     jdvw   the second dimension of the arrays v, w as it appears in
    !            the program that calls vhsgc. jdvw must be at least nlon.
    !
    !     br, bi  two or three dimensional arrays (see input parameter nt)
    !     cr, ci  that contain the vector spherical harmonic coefficients
    !            in the spectral representation of v(i, j) and w(i, j) given
    !            below at the discription of output parameters v and w.
    !
    !     mdab   the first dimension of the arrays br, bi, cr, and ci as it
    !            appears in the program that calls vhsgc. mdab must be at
    !            least min(nlat, nlon/2) if nlon is even or at least
    !            min(nlat, (nlon + 1)/2) if nlon is odd.
    !
    !     ndab   the second dimension of the arrays br, bi, cr, and ci as it
    !            appears in the program that calls vhsgc. ndab must be at
    !            least nlat.
    !
    !     wvhsgc an array which must be initialized by subroutine vhsgci.
    !            once initialized, wvhsgc can be used repeatedly by vhsgc
    !            as long as nlon and nlat remain unchanged.  wvhsgc must
    !            not be altered between calls of vhsgc.
    !
    !     lvhsgc the dimension of the array wvhsgc as it appears in the
    !            program that calls vhsgc. define
    !
    !               l1 = min(nlat, nlon/2) if nlon is even or
    !               l1 = min(nlat, (nlon + 1)/2) if nlon is odd
    !
    !            and
    !
    !               l2 = nlat/2        if nlat is even or
    !               l2 = (nlat + 1)/2    if nlat is odd
    !
    !            then lvhsgc must be at least
    !
    !               4*nlat*l2+3*max(l1-2, 0)*(2*nlat-l1-1)+nlon+15
    !
    !     output parameters
    !
    !     v, w    two or three dimensional arrays (see input parameter nt)
    !            in which the synthesis is stored. v is the colatitudinal
    !            component and w is the east longitudinal component.
    !            v(i, j), w(i, j) contain the components at the gaussian
    !            colatitude theta(i) and longitude phi(j) = (j-1)*2*pi/nlon.
    !            the index ranges are defined above at the input parameter
    !            ityp. v and w are computed from the formulas given below.
    !
    !     define
    !
    !     1.  theta is colatitude and phi is east longitude
    !
    !     2.  the normalized associated legendre funnctions
    !
    !         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)
    !
    !     3.  vbar(m, n, theta) = the derivative of pbar(m, n, theta) with
    !                           respect to theta divided by the square
    !                           root of n(n+1).
    !
    !         vbar(m, n, theta) is more easily computed in the form
    !
    !         vbar(m, n, theta) = (sqrt((n+m)*(n-m+1))*pbar(m-1, n, theta)
    !         -sqrt((n-m)*(n+m+1))*pbar(m+1, n, theta))/(2*sqrt(n*(n+1)))
    !
    !     4.  wbar(m, n, theta) = m/(sin(theta))*pbar(m, n, theta) divided
    !                           by the square root of n(n+1).
    !
    !         wbar(m, n, theta) is more easily computed in the form
    !
    !         wbar(m, n, theta) = sqrt((2n+1)/(2n-1))*(sqrt((n+m)*(n+m-1))
    !         *pbar(m-1, n-1, theta)+sqrt((n-m)*(n-m-1))*pbar(m+1, n-1, theta))
    !         /(2*sqrt(n*(n+1)))
    !
    !
    !    the colatitudnal dependence of the normalized surface vector
    !                spherical harmonics are defined by
    !
    !     5.    bbar(m, n, theta) = (vbar(m, n, theta), i*wbar(m, n, theta))
    !
    !     6.    cbar(m, n, theta) = (i*wbar(m, n, theta), -vbar(m, n, theta))
    !
    !
    !    the coordinate to index mappings
    !
    !     7.   phi(j) = (j-1)*2*pi/nlon, theta(i) is the i(th) guassian
    !          point (see nlat as an input parameter).
    !
    !     the maximum (plus one) longitudinal wave number
    !
    !     8.     mmax = min(nlat, nlon/2) if nlon is even or
    !            mmax = min(nlat, (nlon + 1)/2) if nlon is odd.
    !
    !    if we further define the output vector as
    !
    !     9.    h(i, j) = (v(i, j), w(i, j))
    !
    !    and the complex coefficients
    !
    !     10.   b(m, n) = cmplx(br(m+1, n+1), bi(m+1, n+1))
    !
    !     11.   c(m, n) = cmplx(cr(m+1, n+1), ci(m+1, n+1))
    !
    !
    !    then for i=1, ..., nlat and  j=1, ..., nlon
    !
    !        the expansion for real h(i, j) takes the form
    !
    !     h(i, j) = the sum from n=1 to n=nlat-1 of the real part of
    !
    !         .5*(b(0, n)*bbar(0, n, theta(i))+c(0, n)*cbar(0, n, theta(i)))
    !
    !     plus the sum from m=1 to m=mmax-1 of the sum from n=m to
    !     n=nlat-1 of the real part of
    !
    !              b(m, n)*bbar(m, n, theta(i))*exp(i*m*phi(j))
    !             +c(m, n)*cbar(m, n, theta(i))*exp(i*m*phi(j))
    !
    !   *************************************************************
    !
    !   in terms of real variables this expansion takes the form
    !
    !             for i=1, ..., nlat and  j=1, ..., nlon
    !
    !     v(i, j) = the sum from n=1 to n=nlat-1 of
    !
    !               .5*br(1, n+1)*vbar(0, n, theta(i))
    !
    !     plus the sum from m=1 to m=mmax-1 of the sum from n=m to
    !     n=nlat-1 of the real part of
    !
    !       (br(m+1, n+1)*vbar(m, n, theta(i))-ci(m+1, n+1)*wbar(m, n, theta(i)))
    !                                          *cos(m*phi(j))
    !      -(bi(m+1, n+1)*vbar(m, n, theta(i))+cr(m+1, n+1)*wbar(m, n, theta(i)))
    !                                          *sin(m*phi(j))
    !
    !    and for i=1, ..., nlat and  j=1, ..., nlon
    !
    !     w(i, j) = the sum from n=1 to n=nlat-1 of
    !
    !              -.5*cr(1, n+1)*vbar(0, n, theta(i))
    !
    !     plus the sum from m=1 to m=mmax-1 of the sum from n=m to
    !     n=nlat-1 of the real part of
    !
    !      -(cr(m+1, n+1)*vbar(m, n, theta(i))+bi(m+1, n+1)*wbar(m, n, theta(i)))
    !                                          *cos(m*phi(j))
    !      +(ci(m+1, n+1)*vbar(m, n, theta(i))-br(m+1, n+1)*wbar(m, n, theta(i)))
    !                                          *sin(m*phi(j))
    !
    !
    !      br(m+1, nlat), bi(m+1, nlat), cr(m+1, nlat), and ci(m+1, nlat) are
    !      assumed zero for m even.
    !
    !
    !     ierror = 0  no errors
    !            = 1  error in the specification of nlat
    !            = 2  error in the specification of nlon
    !            = 3  error in the specification of ityp
    !            = 4  error in the specification of nt
    !            = 5  error in the specification of idvw
    !            = 6  error in the specification of jdvw
    !            = 7  error in the specification of mdab
    !            = 8  error in the specification of ndab
    !            = 9  error in the specification of lvhsgc
    !
    module subroutine vhsgc(nlat, nlon, ityp, nt, v, w, idvw, jdvw, br, bi, cr, ci, &
        mdab, ndab, wvhsgc, ierror)

        ! Dummy arguments
        integer(ip), intent(in)  :: nlat
        integer(ip), intent(in)  :: nlon
        integer(ip), intent(in)  :: ityp
        integer(ip), intent(in)  :: nt
        real(wp),    intent(out) :: v(idvw, jdvw, nt)
        real(wp),    intent(out) :: w(idvw, jdvw, nt)
        integer(ip), intent(in)  :: idvw
        integer(ip), intent(in)  :: jdvw
        real(wp),    intent(in)  :: br(mdab, ndab, nt)
        real(wp),    intent(in)  :: bi(mdab, ndab, nt)
        real(wp),    intent(in)  :: cr(mdab, ndab, nt)
        real(wp),    intent(in)  :: ci(mdab, ndab, nt)
        integer(ip), intent(in)  :: mdab
        integer(ip), intent(in)  :: ndab
        real(wp),    intent(in)  :: wvhsgc(:)
        integer(ip), intent(out) :: ierror

        ! Local variables
        integer(ip) :: idv, imid, ist, n1, n2
        integer(ip) :: labc, lnl, required_wavetable_size
        integer(ip) :: mmax, lzz1, lwork
        type(SpherepackUtility) :: util

        ! Check input arguments
        required_wavetable_size = util%get_lvhsgc(nlat, nlon)

        call util%check_vector_transform_inputs(ityp, idvw, jdvw, &
            mdab, ndab, nlat, nlon, nt, required_wavetable_size, &
            wvhsgc, ierror)

        ! Check error flag
        if (ierror /= 0) return

        imid = (nlat + 1)/2
        mmax = min(nlat, (nlon + 1)/2)
        lzz1 = 2*nlat*imid
        labc = 3*(max(mmax-2, 0)*(2*nlat-mmax-1))/2
        n1 = min(nlat, (nlon + 1)/2)
        n2 = (nlat + 1)/2

        ! Set required workspace size
        select case (ityp)
            case (0:2)
                idv = nlat
                ist = imid
                lwork = nlat*(2*nt*nlon+max(6*imid, nlon))
            case default
                idv = imid
                ist = 0
                lwork = imid*(2*nt*nlon+max(6*nlat, nlon))
        end select

        lnl = nt*idv*nlon

        block
            real(wp)    :: work(lwork)
            integer(ip) :: iw1, iw2, iw3, iw4, iw5
            integer(ip) :: jw1, jw2, lwzvin

            iw1 = ist+1
            iw2 = lnl+1
            iw3 = iw2+ist
            iw4 = iw2+lnl
            iw5 = iw4+3*imid*nlat
            lzz1 = 2*nlat*imid
            labc = 3*(max(mmax-2, 0)*(2*nlat-mmax-1))/2
            lwzvin = lzz1+labc
            jw1 = lwzvin+1
            jw2 = jw1+lwzvin

            call vhsgc_lower_utility_routine(nlat, nlon, ityp, nt, imid, &
                idvw, jdvw, v, w, mdab, ndab, br, bi, cr, ci, idv, work, &
                work(iw1:), work(iw2:), work(iw3:), work(iw4:), work(iw5:), &
                wvhsgc, wvhsgc(jw1:), wvhsgc(jw2:))
        end block

    end subroutine vhsgc

    !     subroutine vhsgci(nlat, nlon, wvhsgc, ierror)
    !
    !     subroutine vhsgci initializes the array wvhsgc which can then be
    !     used repeatedly by subroutine vhsgc until nlat or nlon is changed.
    !
    !     input parameters
    !
    !     nlat   the number of points in the gaussian colatitude grid on the
    !            full sphere. these lie in the interval (0, pi) and are computed
    !            in radians in theta(1) <...< theta(nlat) by subroutine compute_gaussian_latitudes_and_weights.
    !            if nlat is odd the equator will be included as the grid point
    !            theta((nlat + 1)/2).  if nlat is even the equator will be
    !            excluded as a grid point and will lie half way between
    !            theta(nlat/2) and theta(nlat/2+1). nlat must be at least 3.
    !            note: on the half sphere, the number of grid points in the
    !            colatitudinal direction is nlat/2 if nlat is even or
    !            (nlat + 1)/2 if nlat is odd.
    !
    !     nlon   the number of distinct londitude points.  nlon determines
    !            the grid increment in longitude as 2*pi/nlon. for example
    !            nlon = 72 for a five degree grid. nlon must be greater
    !            than zero. the axisymmetric case corresponds to nlon=1.
    !            the efficiency of the computation is improved when nlon
    !            is a product of small prime numbers.
    !
    !     lvhsgc the dimension of the array wvhsgc as it appears in the
    !            program that calls vhsgc. define
    !
    !               l1 = min(nlat, nlon/2) if nlon is even or
    !               l1 = min(nlat, (nlon + 1)/2) if nlon is odd
    !
    !            and
    !
    !               l2 = nlat/2        if nlat is even or
    !               l2 = (nlat + 1)/2    if nlat is odd
    !
    !            then lvhsgc must be at least
    !
    !               4*nlat*l2+3*max(l1-2, 0)*(2*nlat-l1-1)+nlon+15
    !
    !
    !     output parameters
    !
    !     wvhsgc an array which is initialized for use by subroutine vhsgc.
    !            once initialized, wvhsgc can be used repeatedly by vhsgc
    !            as long as nlat and nlon remain unchanged.  wvhsgc must not
    !            be altered between calls of vhsgc.
    !
    !
    !     ierror = 0  no errors
    !            = 1  error in the specification of nlat
    !            = 2  error in the specification of nlon
    !            = 3  error in the specification of lvhsgc
    !
    module subroutine vhsgci(nlat, nlon, wvhsgc, ierror)

        ! Dummy arguments
        integer(ip), intent(in)  :: nlat
        integer(ip), intent(in)  :: nlon
        real(wp),    intent(out) :: wvhsgc(:)
        integer(ip), intent(out) :: ierror

        ! Local variables
        type(SpherepackUtility) :: util
        integer(ip) :: imid
        integer(ip) :: iw1
        integer(ip) :: iw2
        integer(ip) :: iwrk
        integer(ip) :: jw1
        integer(ip) :: jw2
        integer(ip) :: jw3
        integer(ip) :: labc
        integer(ip) :: lwvbin
        integer(ip) :: lzz1
        integer(ip) :: mmax, ldwork

        associate (lvhsgc => size(wvhsgc))

            imid = (nlat + 1)/2
            lzz1 = 2*nlat*imid
            mmax = min(nlat, (nlon + 1)/2)
            labc = 3*(max(mmax-2, 0)*(2*nlat-mmax-1))/2

            ! Check calling arguments
            if (nlat < 3) then
                ierror = 1
            else if (nlon < 1) then
                ierror = 2
            else if (lvhsgc < 2*(lzz1+labc)+nlon+15) then
                ierror = 3
            else
                ierror = 0
            end if

            ! Check error flag
            if (ierror /= 0) return

            ! Set workspace index pointers
            jw1 = 1
            jw2 = jw1+nlat
            jw3 = jw2+nlat
            iwrk = (nlat + 1)/2 + 1
            lwvbin = lzz1+labc
            iw1 = lwvbin+1
            iw2 = iw1+lwvbin

            ! Set required workspace size
            ldwork = 2 * nlat * (nlat + 1) + 1

            block
                real(wp) :: dwork(ldwork)

                call compute_gaussian_latitudes_and_weights( &
                    nlat, dwork(jw1:), dwork(jw2:), ierror)

                call util%initialize_polar_components_gaussian_grid(nlat, nlon, dwork, &
                    wvhsgc, dwork(iwrk:))

                call util%initialize_azimuthal_components_gaussian_grid(nlat, nlon, &
                    dwork, wvhsgc(iw1:), dwork(iwrk:))

                call util%hfft%initialize(nlon, wvhsgc(iw2:))
            end block
        end associate

    end subroutine vhsgci

    subroutine vhsgc_lower_utility_routine(nlat, nlon, ityp, nt, imid, idvw, jdvw, v, w, mdab, &
        ndab, br, bi, cr, ci, idv, ve, vo, we, wo, vb, wb, wvbin, wwbin, wrfft)

        real(wp) :: bi
        real(wp) :: br
        real(wp) :: ci
        real(wp) :: cr
        integer(ip) :: i
        integer(ip) :: idv
        integer(ip) :: idvw
        integer(ip) :: imid
        integer(ip) :: imm1
        integer(ip) :: ityp
        integer(ip) :: iv
        integer(ip) :: iw

        integer(ip) :: jdvw
        integer(ip) :: k
        integer(ip) :: m
        integer(ip) :: mdab

        integer(ip) :: mmax
        integer(ip) :: mp1
        integer(ip) :: mp2
        integer(ip) :: ndab
        integer(ip) :: odd_stride
        integer(ip) :: even_stride
        integer(ip) :: nlat
        integer(ip) :: nlon
        integer(ip) :: np1
        integer(ip) :: nt
        real(wp) :: v
        real(wp) :: vb
        real(wp) :: ve
        real(wp) :: vo
        real(wp) :: w
        real(wp) :: wb
        real(wp) :: we
        real(wp) :: wo
        real(wp) :: wrfft
        real(wp) :: wvbin
        real(wp) :: wwbin
        dimension v(idvw, jdvw, nt), w(idvw, jdvw, nt), br(mdab, ndab, nt), &
            bi(mdab, ndab, nt), cr(mdab, ndab, nt), ci(mdab, ndab, nt), &
            ve(idv, nlon, nt), vo(idv, nlon, nt), we(idv, nlon, nt), &
            wo(idv, nlon, nt), wvbin(*), wwbin(*), wrfft(*), &
            vb(imid, nlat, 3), wb(imid, nlat, 3)

        type(VectorSynthesisUtility) :: util

        call util%synthesis_setup(even_stride, imid, imm1, mmax, nlat, &
            odd_stride, ve, vo, we, wo)

        vector_symmetry_cases: select case (ityp)
            case (0)
                ! case ityp=0   no symmetries
                call util%compute_polar_component(0, nlat, nlon, 0, vb, iv, wvbin)

                ! case m = 0
                do k=1, nt
                    do np1=2, even_stride, 2
                        do i=1, imid
                            ve(i, 1, k)=ve(i, 1, k)+br(1, np1, k)*vb(i, np1, iv)
                            we(i, 1, k)=we(i, 1, k)-cr(1, np1, k)*vb(i, np1, iv)
                        end do
                    end do
                end do
                do k=1, nt
                    do np1=3, odd_stride, 2
                        do i=1, imm1
                            vo(i, 1, k)=vo(i, 1, k)+br(1, np1, k)*vb(i, np1, iv)
                            wo(i, 1, k)=wo(i, 1, k)-cr(1, np1, k)*vb(i, np1, iv)
                        end do
                    end do
                end do

                ! case m = 1 through nlat-1
                if (mmax < 2) exit vector_symmetry_cases
                do mp1=2, mmax
                    m = mp1-1
                    mp2 = mp1+1
                    call util%compute_polar_component(0, nlat, nlon, m, vb, iv, wvbin)
                    call util%compute_azimuthal_component(0, nlat, nlon, m, wb, iw, wwbin)
                    if (mp1 <= odd_stride) then
                        do k=1, nt
                            do np1=mp1, odd_stride, 2
                                do i=1, imm1
                                    vo(i, 2*mp1-2, k) = vo(i, 2*mp1-2, k)+br(mp1, np1, k)*vb(i, np1, iv)
                                    ve(i, 2*mp1-2, k) = ve(i, 2*mp1-2, k)-ci(mp1, np1, k)*wb(i, np1, iw)
                                    vo(i, 2*mp1-1, k) = vo(i, 2*mp1-1, k)+bi(mp1, np1, k)*vb(i, np1, iv)
                                    ve(i, 2*mp1-1, k) = ve(i, 2*mp1-1, k)+cr(mp1, np1, k)*wb(i, np1, iw)
                                    wo(i, 2*mp1-2, k) = wo(i, 2*mp1-2, k)-cr(mp1, np1, k)*vb(i, np1, iv)
                                    we(i, 2*mp1-2, k) = we(i, 2*mp1-2, k)-bi(mp1, np1, k)*wb(i, np1, iw)
                                    wo(i, 2*mp1-1, k) = wo(i, 2*mp1-1, k)-ci(mp1, np1, k)*vb(i, np1, iv)
                                    we(i, 2*mp1-1, k) = we(i, 2*mp1-1, k)+br(mp1, np1, k)*wb(i, np1, iw)
                                end do
                                if (odd(nlat)) then
                                    ve(imid, 2*mp1-2, k) = ve(imid, 2*mp1-2, k) &
                                        -ci(mp1, np1, k)*wb(imid, np1, iw)
                                    ve(imid, 2*mp1-1, k) = ve(imid, 2*mp1-1, k) &
                                        +cr(mp1, np1, k)*wb(imid, np1, iw)
                                    we(imid, 2*mp1-2, k) = we(imid, 2*mp1-2, k) &
                                        -bi(mp1, np1, k)*wb(imid, np1, iw)
                                    we(imid, 2*mp1-1, k) = we(imid, 2*mp1-1, k) &
                                        +br(mp1, np1, k)*wb(imid, np1, iw)
                                end if
                            end do
                        end do
                    end if
                    if (mp2 <= even_stride) then
                        do k=1, nt
                            do np1=mp2, even_stride, 2
                                do i=1, imm1
                                    ve(i, 2*mp1-2, k) = ve(i, 2*mp1-2, k)+br(mp1, np1, k)*vb(i, np1, iv)
                                    vo(i, 2*mp1-2, k) = vo(i, 2*mp1-2, k)-ci(mp1, np1, k)*wb(i, np1, iw)
                                    ve(i, 2*mp1-1, k) = ve(i, 2*mp1-1, k)+bi(mp1, np1, k)*vb(i, np1, iv)
                                    vo(i, 2*mp1-1, k) = vo(i, 2*mp1-1, k)+cr(mp1, np1, k)*wb(i, np1, iw)
                                    we(i, 2*mp1-2, k) = we(i, 2*mp1-2, k)-cr(mp1, np1, k)*vb(i, np1, iv)
                                    wo(i, 2*mp1-2, k) = wo(i, 2*mp1-2, k)-bi(mp1, np1, k)*wb(i, np1, iw)
                                    we(i, 2*mp1-1, k) = we(i, 2*mp1-1, k)-ci(mp1, np1, k)*vb(i, np1, iv)
                                    wo(i, 2*mp1-1, k) = wo(i, 2*mp1-1, k)+br(mp1, np1, k)*wb(i, np1, iw)
                                end do
                                if (odd(nlat)) then
                                    ve(imid, 2*mp1-2, k) = ve(imid, 2*mp1-2, k) &
                                        +br(mp1, np1, k)*vb(imid, np1, iv)
                                    ve(imid, 2*mp1-1, k) = ve(imid, 2*mp1-1, k) &
                                        +bi(mp1, np1, k)*vb(imid, np1, iv)
                                    we(imid, 2*mp1-2, k) = we(imid, 2*mp1-2, k) &
                                        -cr(mp1, np1, k)*vb(imid, np1, iv)
                                    we(imid, 2*mp1-1, k) = we(imid, 2*mp1-1, k) &
                                        -ci(mp1, np1, k)*vb(imid, np1, iv)
                                end if
                            end do
                        end do
                    end if
                end do
            case(1)
                ! case ityp=1   no symmetries,  cr and ci equal zero
                call util%compute_polar_component(0, nlat, nlon, 0, vb, iv, wvbin)

                ! case m = 0
                do k=1, nt
                    do np1=2, even_stride, 2
                        do i=1, imid
                            ve(i, 1, k)=ve(i, 1, k)+br(1, np1, k)*vb(i, np1, iv)
                        end do
                    end do
                end do
                do k=1, nt
                    do np1=3, odd_stride, 2
                        do i=1, imm1
                            vo(i, 1, k)=vo(i, 1, k)+br(1, np1, k)*vb(i, np1, iv)
                        end do
                    end do
                end do

                ! case m = 1 through nlat-1
                if (mmax < 2) exit vector_symmetry_cases
                do mp1=2, mmax
                    m = mp1-1
                    mp2 = mp1+1
                    call util%compute_polar_component(0, nlat, nlon, m, vb, iv, wvbin)
                    call util%compute_azimuthal_component(0, nlat, nlon, m, wb, iw, wwbin)
                    if (mp1 <= odd_stride) then
                        do k=1, nt
                            do np1=mp1, odd_stride, 2
                                do i=1, imm1
                                    vo(i, 2*mp1-2, k) = vo(i, 2*mp1-2, k)+br(mp1, np1, k)*vb(i, np1, iv)
                                    vo(i, 2*mp1-1, k) = vo(i, 2*mp1-1, k)+bi(mp1, np1, k)*vb(i, np1, iv)
                                    we(i, 2*mp1-2, k) = we(i, 2*mp1-2, k)-bi(mp1, np1, k)*wb(i, np1, iw)
                                    we(i, 2*mp1-1, k) = we(i, 2*mp1-1, k)+br(mp1, np1, k)*wb(i, np1, iw)
                                end do
                                if (odd(nlat)) then
                                    we(imid, 2*mp1-2, k) = we(imid, 2*mp1-2, k) &
                                        -bi(mp1, np1, k)*wb(imid, np1, iw)
                                    we(imid, 2*mp1-1, k) = we(imid, 2*mp1-1, k) &
                                        +br(mp1, np1, k)*wb(imid, np1, iw)
                                end if
                            end do
                        end do
                    end if
                    if (mp2 <= even_stride) then
                        do k=1, nt
                            do np1=mp2, even_stride, 2
                                do i=1, imm1
                                    ve(i, 2*mp1-2, k) = ve(i, 2*mp1-2, k)+br(mp1, np1, k)*vb(i, np1, iv)
                                    ve(i, 2*mp1-1, k) = ve(i, 2*mp1-1, k)+bi(mp1, np1, k)*vb(i, np1, iv)
                                    wo(i, 2*mp1-2, k) = wo(i, 2*mp1-2, k)-bi(mp1, np1, k)*wb(i, np1, iw)
                                    wo(i, 2*mp1-1, k) = wo(i, 2*mp1-1, k)+br(mp1, np1, k)*wb(i, np1, iw)
                                end do
                                if (odd(nlat)) then
                                    ve(imid, 2*mp1-2, k) = ve(imid, 2*mp1-2, k) &
                                        +br(mp1, np1, k)*vb(imid, np1, iv)
                                    ve(imid, 2*mp1-1, k) = ve(imid, 2*mp1-1, k) &
                                        +bi(mp1, np1, k)*vb(imid, np1, iv)
                                end if
                            end do
                        end do
                    end if
                end do
            case(2)
                ! case ityp=2   no symmetries,  br and bi are equal to zero
                call util%compute_polar_component(0, nlat, nlon, 0, vb, iv, wvbin)

                ! case m = 0
                do k=1, nt
                    do np1=2, even_stride, 2
                        do i=1, imid
                            we(i, 1, k)=we(i, 1, k)-cr(1, np1, k)*vb(i, np1, iv)
                        end do
                    end do
                end do
                do k=1, nt
                    do np1=3, odd_stride, 2
                        do i=1, imm1
                            wo(i, 1, k)=wo(i, 1, k)-cr(1, np1, k)*vb(i, np1, iv)
                        end do
                    end do
                end do

                ! case m = 1 through nlat-1
                if (mmax < 2) exit vector_symmetry_cases
                do mp1=2, mmax
                    m = mp1-1
                    mp2 = mp1+1
                    call util%compute_polar_component(0, nlat, nlon, m, vb, iv, wvbin)
                    call util%compute_azimuthal_component(0, nlat, nlon, m, wb, iw, wwbin)
                    if (mp1 <= odd_stride) then
                        do k=1, nt
                            do np1=mp1, odd_stride, 2
                                do i=1, imm1
                                    ve(i, 2*mp1-2, k) = ve(i, 2*mp1-2, k)-ci(mp1, np1, k)*wb(i, np1, iw)
                                    ve(i, 2*mp1-1, k) = ve(i, 2*mp1-1, k)+cr(mp1, np1, k)*wb(i, np1, iw)
                                    wo(i, 2*mp1-2, k) = wo(i, 2*mp1-2, k)-cr(mp1, np1, k)*vb(i, np1, iv)
                                    wo(i, 2*mp1-1, k) = wo(i, 2*mp1-1, k)-ci(mp1, np1, k)*vb(i, np1, iv)
                                end do
                                if (odd(nlat)) then
                                    ve(imid, 2*mp1-2, k) = ve(imid, 2*mp1-2, k) &
                                        -ci(mp1, np1, k)*wb(imid, np1, iw)
                                    ve(imid, 2*mp1-1, k) = ve(imid, 2*mp1-1, k) &
                                        +cr(mp1, np1, k)*wb(imid, np1, iw)
                                end if
                            end do
                        end do
                    end if
                    if (mp2 <= even_stride) then
                        do k=1, nt
                            do np1=mp2, even_stride, 2
                                do i=1, imm1
                                    vo(i, 2*mp1-2, k) = vo(i, 2*mp1-2, k)-ci(mp1, np1, k)*wb(i, np1, iw)
                                    vo(i, 2*mp1-1, k) = vo(i, 2*mp1-1, k)+cr(mp1, np1, k)*wb(i, np1, iw)
                                    we(i, 2*mp1-2, k) = we(i, 2*mp1-2, k)-cr(mp1, np1, k)*vb(i, np1, iv)
                                    we(i, 2*mp1-1, k) = we(i, 2*mp1-1, k)-ci(mp1, np1, k)*vb(i, np1, iv)
                                end do
                                if (odd(nlat)) then
                                    we(imid, 2*mp1-2, k) = we(imid, 2*mp1-2, k) &
                                        -cr(mp1, np1, k)*vb(imid, np1, iv)
                                    we(imid, 2*mp1-1, k) = we(imid, 2*mp1-1, k) &
                                        -ci(mp1, np1, k)*vb(imid, np1, iv)
                                end if
                            end do
                        end do
                    end if
                end do
            case(3)
                ! case ityp=3   v even,  w odd
                call util%compute_polar_component(0, nlat, nlon, 0, vb, iv, wvbin)

                ! case m = 0
                do k=1, nt
                    do np1=2, even_stride, 2
                        do i=1, imid
                            ve(i, 1, k)=ve(i, 1, k)+br(1, np1, k)*vb(i, np1, iv)
                        end do
                    end do
                end do

                do k=1, nt
                    do np1=3, odd_stride, 2
                        do i=1, imm1
                            wo(i, 1, k)=wo(i, 1, k)-cr(1, np1, k)*vb(i, np1, iv)
                        end do
                    end do
                end do

                ! case m = 1 through nlat-1
                if (mmax < 2) exit vector_symmetry_cases
                do mp1=2, mmax
                    m = mp1-1
                    mp2 = mp1+1
                    call util%compute_polar_component(0, nlat, nlon, m, vb, iv, wvbin)
                    call util%compute_azimuthal_component(0, nlat, nlon, m, wb, iw, wwbin)
                    if (mp1 <= odd_stride) then
                        do k=1, nt
                            do np1=mp1, odd_stride, 2
                                do i=1, imm1
                                    ve(i, 2*mp1-2, k) = ve(i, 2*mp1-2, k)-ci(mp1, np1, k)*wb(i, np1, iw)
                                    ve(i, 2*mp1-1, k) = ve(i, 2*mp1-1, k)+cr(mp1, np1, k)*wb(i, np1, iw)
                                    wo(i, 2*mp1-2, k) = wo(i, 2*mp1-2, k)-cr(mp1, np1, k)*vb(i, np1, iv)
                                    wo(i, 2*mp1-1, k) = wo(i, 2*mp1-1, k)-ci(mp1, np1, k)*vb(i, np1, iv)
                                end do
                                if (odd(nlat)) then
                                    ve(imid, 2*mp1-2, k) = ve(imid, 2*mp1-2, k) &
                                        -ci(mp1, np1, k)*wb(imid, np1, iw)
                                    ve(imid, 2*mp1-1, k) = ve(imid, 2*mp1-1, k) &
                                        +cr(mp1, np1, k)*wb(imid, np1, iw)
                                end if
                            end do
                        end do
                    end if
                    if (mp2 <= even_stride) then
                        do k=1, nt
                            do np1=mp2, even_stride, 2
                                do i=1, imm1
                                    ve(i, 2*mp1-2, k) = ve(i, 2*mp1-2, k)+br(mp1, np1, k)*vb(i, np1, iv)
                                    ve(i, 2*mp1-1, k) = ve(i, 2*mp1-1, k)+bi(mp1, np1, k)*vb(i, np1, iv)
                                    wo(i, 2*mp1-2, k) = wo(i, 2*mp1-2, k)-bi(mp1, np1, k)*wb(i, np1, iw)
                                    wo(i, 2*mp1-1, k) = wo(i, 2*mp1-1, k)+br(mp1, np1, k)*wb(i, np1, iw)
                                end do
                                if (odd(nlat)) then
                                    ve(imid, 2*mp1-2, k) = ve(imid, 2*mp1-2, k) &
                                        +br(mp1, np1, k)*vb(imid, np1, iv)
                                    ve(imid, 2*mp1-1, k) = ve(imid, 2*mp1-1, k) &
                                        +bi(mp1, np1, k)*vb(imid, np1, iv)
                                end if
                            end do
                        end do
                    end if
                end do
            case(4)
                ! case ityp=4   v even,  w odd, and both cr and ci equal zero
                call util%compute_polar_component(1, nlat, nlon, 0, vb, iv, wvbin)

                ! case m = 0
                do k=1, nt
                    do np1=2, even_stride, 2
                        do i=1, imid
                            ve(i, 1, k)=ve(i, 1, k)+br(1, np1, k)*vb(i, np1, iv)
                        end do
                    end do
                end do

                ! case m = 1 through nlat-1
                if (mmax < 2) exit vector_symmetry_cases
                do mp1=2, mmax
                    m = mp1-1
                    mp2 = mp1+1
                    call util%compute_polar_component(1, nlat, nlon, m, vb, iv, wvbin)
                    call util%compute_azimuthal_component(1, nlat, nlon, m, wb, iw, wwbin)
                    if (mp2 <= even_stride) then
                        do k=1, nt
                            do np1=mp2, even_stride, 2
                                do i=1, imm1
                                    ve(i, 2*mp1-2, k) = ve(i, 2*mp1-2, k)+br(mp1, np1, k)*vb(i, np1, iv)
                                    ve(i, 2*mp1-1, k) = ve(i, 2*mp1-1, k)+bi(mp1, np1, k)*vb(i, np1, iv)
                                    wo(i, 2*mp1-2, k) = wo(i, 2*mp1-2, k)-bi(mp1, np1, k)*wb(i, np1, iw)
                                    wo(i, 2*mp1-1, k) = wo(i, 2*mp1-1, k)+br(mp1, np1, k)*wb(i, np1, iw)
                                end do
                                if (odd(nlat)) then
                                    ve(imid, 2*mp1-2, k) = ve(imid, 2*mp1-2, k) &
                                        +br(mp1, np1, k)*vb(imid, np1, iv)
                                    ve(imid, 2*mp1-1, k) = ve(imid, 2*mp1-1, k) &
                                        +bi(mp1, np1, k)*vb(imid, np1, iv)
                                end if
                            end do
                        end do
                    end if
                end do
            case(5)
                ! case ityp=5   v even,  w odd, br and bi equal zero
                call util%compute_polar_component(2, nlat, nlon, 0, vb, iv, wvbin)

                ! case m = 0
                do k=1, nt
                    do np1=3, odd_stride, 2
                        do i=1, imm1
                            wo(i, 1, k)=wo(i, 1, k)-cr(1, np1, k)*vb(i, np1, iv)
                        end do
                    end do
                end do

                ! case m = 1 through nlat-1
                if (mmax < 2) exit vector_symmetry_cases
                do mp1=2, mmax
                    m = mp1-1
                    mp2 = mp1+1
                    call util%compute_polar_component(2, nlat, nlon, m, vb, iv, wvbin)
                    call util%compute_azimuthal_component(2, nlat, nlon, m, wb, iw, wwbin)
                    if (mp1 <= odd_stride) then
                        do k=1, nt
                            do np1=mp1, odd_stride, 2
                                do i=1, imm1
                                    ve(i, 2*mp1-2, k) = ve(i, 2*mp1-2, k)-ci(mp1, np1, k)*wb(i, np1, iw)
                                    ve(i, 2*mp1-1, k) = ve(i, 2*mp1-1, k)+cr(mp1, np1, k)*wb(i, np1, iw)
                                    wo(i, 2*mp1-2, k) = wo(i, 2*mp1-2, k)-cr(mp1, np1, k)*vb(i, np1, iv)
                                    wo(i, 2*mp1-1, k) = wo(i, 2*mp1-1, k)-ci(mp1, np1, k)*vb(i, np1, iv)
                                end do
                                if (odd(nlat)) then
                                    ve(imid, 2*mp1-2, k) = ve(imid, 2*mp1-2, k) &
                                        -ci(mp1, np1, k)*wb(imid, np1, iw)
                                    ve(imid, 2*mp1-1, k) = ve(imid, 2*mp1-1, k) &
                                        +cr(mp1, np1, k)*wb(imid, np1, iw)
                                end if
                            end do
                        end do
                    end if
                end do
            case(6)
                ! case ityp=6   v odd  ,  w even
                call util%compute_polar_component(0, nlat, nlon, 0, vb, iv, wvbin)

                ! case m = 0
                do k=1, nt
                    do np1=2, even_stride, 2
                        do i=1, imid
                            we(i, 1, k)=we(i, 1, k)-cr(1, np1, k)*vb(i, np1, iv)
                        end do
                    end do
                end do
                do k=1, nt
                    do np1=3, odd_stride, 2
                        do i=1, imm1
                            vo(i, 1, k)=vo(i, 1, k)+br(1, np1, k)*vb(i, np1, iv)
                        end do
                    end do
                end do

                ! case m = 1 through nlat-1
                if (mmax < 2) exit vector_symmetry_cases
                do mp1=2, mmax
                    m = mp1-1
                    mp2 = mp1+1
                    call util%compute_polar_component(0, nlat, nlon, m, vb, iv, wvbin)
                    call util%compute_azimuthal_component(0, nlat, nlon, m, wb, iw, wwbin)
                    if (mp1 <= odd_stride) then
                        do k=1, nt
                            do np1=mp1, odd_stride, 2
                                do i=1, imm1
                                    vo(i, 2*mp1-2, k) = vo(i, 2*mp1-2, k)+br(mp1, np1, k)*vb(i, np1, iv)
                                    vo(i, 2*mp1-1, k) = vo(i, 2*mp1-1, k)+bi(mp1, np1, k)*vb(i, np1, iv)
                                    we(i, 2*mp1-2, k) = we(i, 2*mp1-2, k)-bi(mp1, np1, k)*wb(i, np1, iw)
                                    we(i, 2*mp1-1, k) = we(i, 2*mp1-1, k)+br(mp1, np1, k)*wb(i, np1, iw)
                                end do
                                if (odd(nlat)) then
                                    we(imid, 2*mp1-2, k) = we(imid, 2*mp1-2, k) &
                                        -bi(mp1, np1, k)*wb(imid, np1, iw)
                                    we(imid, 2*mp1-1, k) = we(imid, 2*mp1-1, k) &
                                        +br(mp1, np1, k)*wb(imid, np1, iw)
                                end if
                            end do
                        end do
                    end if
                    if (mp2 <= even_stride) then
                        do k=1, nt
                            do np1=mp2, even_stride, 2
                                do i=1, imm1
                                    vo(i, 2*mp1-2, k) = vo(i, 2*mp1-2, k)-ci(mp1, np1, k)*wb(i, np1, iw)
                                    vo(i, 2*mp1-1, k) = vo(i, 2*mp1-1, k)+cr(mp1, np1, k)*wb(i, np1, iw)
                                    we(i, 2*mp1-2, k) = we(i, 2*mp1-2, k)-cr(mp1, np1, k)*vb(i, np1, iv)
                                    we(i, 2*mp1-1, k) = we(i, 2*mp1-1, k)-ci(mp1, np1, k)*vb(i, np1, iv)
                                end do
                                if (odd(nlat)) then
                                    we(imid, 2*mp1-2, k) = we(imid, 2*mp1-2, k) &
                                        -cr(mp1, np1, k)*vb(imid, np1, iv)
                                    we(imid, 2*mp1-1, k) = we(imid, 2*mp1-1, k) &
                                        -ci(mp1, np1, k)*vb(imid, np1, iv)
                                end if
                            end do
                        end do
                    end if
                end do
            case(7)
                ! case ityp=7   v odd, w even   cr and ci equal zero
                call util%compute_polar_component(2, nlat, nlon, 0, vb, iv, wvbin)

                ! case m = 0
                do k=1, nt
                    do np1=3, odd_stride, 2
                        do i=1, imm1
                            vo(i, 1, k)=vo(i, 1, k)+br(1, np1, k)*vb(i, np1, iv)
                        end do
                    end do
                end do

                ! case m = 1 through nlat-1
                if (mmax < 2) exit vector_symmetry_cases
                do mp1=2, mmax
                    m = mp1-1
                    mp2 = mp1+1
                    call util%compute_polar_component(2, nlat, nlon, m, vb, iv, wvbin)
                    call util%compute_azimuthal_component(2, nlat, nlon, m, wb, iw, wwbin)
                    if (mp1 <= odd_stride) then
                        do k=1, nt
                            do np1=mp1, odd_stride, 2
                                do i=1, imm1
                                    vo(i, 2*mp1-2, k) = vo(i, 2*mp1-2, k)+br(mp1, np1, k)*vb(i, np1, iv)
                                    vo(i, 2*mp1-1, k) = vo(i, 2*mp1-1, k)+bi(mp1, np1, k)*vb(i, np1, iv)
                                    we(i, 2*mp1-2, k) = we(i, 2*mp1-2, k)-bi(mp1, np1, k)*wb(i, np1, iw)
                                    we(i, 2*mp1-1, k) = we(i, 2*mp1-1, k)+br(mp1, np1, k)*wb(i, np1, iw)
                                end do
                                if (odd(nlat)) then
                                    we(imid, 2*mp1-2, k) = we(imid, 2*mp1-2, k) &
                                        -bi(mp1, np1, k)*wb(imid, np1, iw)
                                    we(imid, 2*mp1-1, k) = we(imid, 2*mp1-1, k) &
                                        +br(mp1, np1, k)*wb(imid, np1, iw)
                                end if
                            end do
                        end do
                    end if
                end do
            case(8)
                ! case ityp=8   v odd,  w even   br and bi equal zero
                call util%compute_polar_component(1, nlat, nlon, 0, vb, iv, wvbin)

                ! case m = 0
                do k=1, nt
                    do np1=2, even_stride, 2
                        do i=1, imid
                            we(i, 1, k)=we(i, 1, k)-cr(1, np1, k)*vb(i, np1, iv)
                        end do
                    end do
                end do

                ! case m = 1 through nlat-1
                if (mmax < 2) exit vector_symmetry_cases
                do mp1=2, mmax
                    m = mp1-1
                    mp2 = mp1+1
                    call util%compute_polar_component(1, nlat, nlon, m, vb, iv, wvbin)
                    call util%compute_azimuthal_component(1, nlat, nlon, m, wb, iw, wwbin)
                    if (mp2 <= even_stride) then
                        do k=1, nt
                            do np1=mp2, even_stride, 2
                                do i=1, imm1
                                    vo(i, 2*mp1-2, k) = vo(i, 2*mp1-2, k)-ci(mp1, np1, k)*wb(i, np1, iw)
                                    vo(i, 2*mp1-1, k) = vo(i, 2*mp1-1, k)+cr(mp1, np1, k)*wb(i, np1, iw)
                                    we(i, 2*mp1-2, k) = we(i, 2*mp1-2, k)-cr(mp1, np1, k)*vb(i, np1, iv)
                                    we(i, 2*mp1-1, k) = we(i, 2*mp1-1, k)-ci(mp1, np1, k)*vb(i, np1, iv)
                                end do
                                if (odd(nlat)) then
                                    we(imid, 2*mp1-2, k) = we(imid, 2*mp1-2, k) &
                                        -cr(mp1, np1, k)*vb(imid, np1, iv)
                                    we(imid, 2*mp1-1, k) = we(imid, 2*mp1-1, k) &
                                        -ci(mp1, np1, k)*vb(imid, np1, iv)
                                end if
                            end do
                        end do
                    end if
                end do
        end select vector_symmetry_cases

        call util%assemble_transform(idvw, jdvw, idv, imid, &
            imm1, ityp, nlat, nlon, nt, v, ve, vo, w, we, wo, wrfft)

    end subroutine vhsgc_lower_utility_routine

end submodule vector_synthesis_gaussian_grid


text:
The Causeway Classic, named for the Yolo Causeway that connects Davis to Sacramento, is an annual UC Davis Football football game where UC Davis plays against their rival, Sacramento State University Sacramento State. The game is analagous to the Big Game where Cal plays Stanford.

2010 Causeway Classic
 In 2010, the Sac State Hornets were heavily favored to win and the game was played in a downpour from the 2nd quarter on. Sac quickly went up 130, including a score due to a blocked punt... But a fumble by Sac State let to a 137 deficit for the Ags. Then a Touchdown strike made it 1413. The conditions at Aggie Stadium were very bad; there were several missed field goals. Finally, Sac State got a FG to end the first half at 1614 Sac State. Then, at halftime, lightning and hail moved in. The lightning came very close to the stadium (possibly even as close as a 1/2 or 1/4 of a mile). The game was delayed at the end of halftime for 47 minutes due to the lightning danger. Everyone was cleared out of the seats due to the threat of lightning plus wet aluminum seats. Most of the crowd left to go home and/or watch the rest on tv at this time due to the lightning and heavy rain but those who stayed were rewarded in the end. The 3rd Quarter was sloppy and nothing happened. In the 4th, The Ags punted to Sac after blowing a scoring chance, but a freak play saved the Ags. The the ball bounced off the arms of the returner from Sac State and it was caught PERFECTLY by the aggies who ran into the end zone (but were given possession at the 7 by rule). The ags kicked a fg to make it 1716 UCD with about 6 mins left. A few drives later, with less than three minutes left, the Ags were pinned down inside their own 10 yrd line and were forced to punt from their own end zone. The Aggies needed a long, good punt in order to prevent the Hornets from chewing up the clock and scoring the final field goal or touchdown without time for the Aggies to respond. They needed every yard. Sacramento State brought in their quarterback from last year as a sub (He led the Hornets to a comeback win at the Hornets Nest in 2009 against UCD in a game that is controversial b/c of a Chris Carter catch that was not called, which led to a collapse in that game. It is recalled with angst among many aggie faithful who witnessed the game last year). This year, he threw many long, successful passes and it looked bad for Davis. Sac State was within the 30 yard line of Davis with about a minute to go and a timeout or two remaining. However, the Aggies forced a Sac State sack at the 37 that put them out of field goal range. State called time out and was forced to pass on 2nd and 17. On 2nd and 3rd down, the Ags nearly picked off the pass both times. Finally on 4th and 17, UC Davis prevailed on an INTERCEPTION that sealed the deal for UC Davis first win since 2007 and first win over Sac State at Aggie Stadium. The win crushed Sac States hopes for the playoffs.

There is also a track dual meet between UCD and CSUS that is called the Causeway Classic but the two events do not take place at the same time.

To add to the confusion there is also a pizza served at Woodstocks Pizza Woodstocks called the Causeway Classic.

Some more history of the Causeway Classic can be found on the http://www.csus.edu/pubaf/journal/fall2003/26causeway.htm Sac State webpage

The State Hornet, Sac States student newspaper posted an timeline with clickable game capsules for every regular season meeting of the two teams up to 2010 as a preview for the 2011 Causeway Classic at Hornet Stadium. It can be accessed here: http://www.statehornet.com/app/Causeway%20Classic%202011%20Timeline/index.html

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
PROGRAM HBOND

IMPLICIT NONE

INTEGER, PARAMETER :: ik8 = SELECTED_INT_KIND (8)
INTEGER, PARAMETER :: ik4 = SELECTED_INT_KIND (4)
INTEGER, PARAMETER :: rk8 = SELECTED_REAL_KIND (8)
INTEGER, PARAMETER :: prec = SELECTED_REAL_KIND (8) !for the sorting routine

LOGICAL check
INTEGER (KIND = ik4) :: dt, id, n_ion, n_atom, n_header, mol, typ, nn
INTEGER (KIND = ik8) :: i, j, k, m, l, time_start, time_stop, mino
INTEGER (KIND = ik8) :: n_step, counter, ind, a_count, a_end
INTEGER (KIND = ik8) :: b_start, a_start, b_end
INTEGER (KIND = ik4), ALLOCATABLE, DIMENSION (:, :) :: Neigh_O !O neigh list
INTEGER (KIND = ik4), ALLOCATABLE, DIMENSION (:, :) :: NH4 !ammonium atom list
REAL (KIND = rk8) :: db, da, max_bond, min_bond, max_angle, min_angle, dist
REAL (KIND = rk8) :: dx, dy, dz, vx1, vx2, vy1, vy2, vz1, vz2, v1mag, v2mag
REAL (KIND = rk8) :: dot, angle, bx, by, bz, x, y, z, x1, y1, z1, x2, y2, z2
REAL (KIND = rk8) :: mindist, x3, y3, z3
REAL (KIND = rk8), DIMENSION (4) :: Bavg !Average hydrogen bond values
REAL (KIND = rk8), DIMENSION (4) :: Aavg !Average hydrogen bond angle values
REAL (KIND = rk8), DIMENSION (4) :: Bavgsort !sorted avg bond values
REAL (KIND = rk8), DIMENSION (4) :: Aavgsort !sorted avg hbond values
REAL (KIND = rk8), DIMENSION (4) :: Btemp1 !temp bond values
REAL (KIND = rk8), DIMENSION (4) :: Btemp2 !sorted temp bond values
REAL (KIND = rk8), DIMENSION (4) :: Atemp1 !temp angle values
REAL (KIND = rk8), DIMENSION (4) :: Atemp2 !sorted temp angle values
REAL (KIND = rk8), ALLOCATABLE, DIMENSION (:) :: Olist !sorted oxygen list
REAL (KIND = rk8), ALLOCATABLE, DIMENSION (:, :) :: Bonds !bond histograms
REAL (KIND = rk8), ALLOCATABLE, DIMENSION (:, :) :: Angles !angle histograms
REAL (KIND = rk8), ALLOCATABLE, DIMENSION (:, :, :) :: O_Dist !N-O distances
REAL (KIND = rk8), ALLOCATABLE, DIMENSION (:, :, :) :: Coords !atom info
REAL (KIND = rk8), ALLOCATABLE, DIMENSION (:, :, :) :: Box !box size info
CHARACTER filename1*40, filename2*40

nn = 10 !number of nearest oxygens for neighborlist
n_header = 5 !header lines to skip in trajectory
n_atom = 15120 !number of atoms
n_ion = 1512 !number of each type of ion
time_start = 100000 !starting timestep
time_stop = 200000 !ending timestep
dt = 10000 !frequency of trajectory dump
n_step = (time_stop - time_start) / dt + 1 !snapshots in trajectory
db = 0.00001 !resolution for bond histogram
da = 0.001 !resolution for angle histogram
max_bond = 5.0 !maximum hydrogen bond length for histogram - Changed from 3.5 on
!10/15
min_bond = 0.10 !minimum hydrogen bond length for histogram - Changed from 0.3
!on 10/15
max_angle = 180.0 !maximum hydrogen bond angle for histogram
min_angle = 1.0 !minimum hydrogen bond angle for histogram - Changed from 30 on
!10/15
b_start = INT (min_bond / db) - 1 !lower index for bond histogram
b_end = INT (max_bond / db) + 1 !upper index for bond histogram
a_start = INT (min_angle / da) - 1 !lower index for angle histogram
a_end = INT (max_angle / da) + 1 !upper index for angle histogram

!WRITE (*, *) 'Allocating'

ALLOCATE (Olist (n_ion * 4)) !total number of O atoms per snapshot
ALLOCATE (Neigh_O (n_ion, 10)) !number of O neigh per NH4
ALLOCATE (NH4 (n_ion, 5)) !number of atoms per NH4
ALLOCATE (Coords (5, n_step, n_atom), Box (2, 3, n_step)) !atom + box info
ALLOCATE (Bonds (4, b_start:b_end), Angles (4, a_start:b_end)) !histogramming
ALLOCATE (O_Dist (n_ion, 2, n_ion * 4)) !for each NH4, all O-N dist + O atom id

!Initialize some arrays
NH4 (:, :) = 0
Neigh_O (:, :) = 0
Coords (:, :, :) = 0.0
Box (:, :, :) = 0.0
O_Dist (:, :, :) = 0.0

!Open trajectory file
OPEN (UNIT = 10, FILE = 'out.lammpstrj', STATUS = 'OLD')

!WRITE (*, *) 'Reading'

!Loop over snap shots of trajectory file
DO i = 1, n_step

  !Skip over header lines in trajectory file
  DO j = 1, n_header

    READ (10, *)

  ENDDO

  !Read box dimensions into an array
  DO j = 1, 3

    READ (10, *) Box (1, j, i), Box (2, j, i)

  ENDDO

  !Skip line describing elements listed per line in trajectory
  READ (10, *)

  !Loop over total number of atoms, saving per atom information
  DO j = 1, n_atom

    READ (10, *) id, mol, typ, x, y, z
    Coords (1, i, id) = x
    Coords (2, i, id) = y
    Coords (3, i, id) = z
    Coords (4, i, id) = mol
    Coords (5, i, id) = typ

  ENDDO

ENDDO

!WRITE (*, *) 'Done Reading'

!Create separate list for NH4 atom ids, sorted by mol id for easy reference
!Loop over total number of atoms in first snapshot
!NH4 array = (1:n_ion, 5) where (:, 1) are N and (:, 2:5) are H
DO j = 1, n_atom

  !Find Nitrogen atoms, type = 3
  IF (Coords (5, 1, j) .EQ. 3) THEN

      NH4 (Coords (4, 1, j) - n_ion, 1) = j

    !Find Hydrogen atoms, type = 2
    ELSE IF (Coords (5, 1, j) .EQ. 2) THEN

      IF (NH4 (Coords (4, 1, j) - n_ion, 2) .EQ. 0) THEN

          NH4 (Coords (4, 1, j) - n_ion, 2) = j

        ELSE IF (NH4 (Coords (4, 1, j) - n_ion, 3) .EQ. 0) THEN

          NH4 (Coords (4, 1, j) - n_ion, 3) = j

        ELSE IF (NH4 (Coords (4, 1, j) - n_ion, 4) .EQ. 0) THEN

          NH4 (Coords (4, 1, j) - n_ion, 4) = j

        ELSE

          NH4 (Coords (4, 1, j) - n_ion, 5) = j

      END IF

    ELSE

      CYCLE

  END IF

END DO

!WRITE (*, *) 'Assigned NH4 array'

!Calculate box lengths for first snapshot and initialize counter
bx = Box (2, 1, 1) - Box (1, 1, 1)
by = Box (2, 2, 1) - Box (1, 2, 1)
bz = Box (2, 3, 1) - Box (1, 3, 1)
counter = 0 !keeps track of number of oxygen atoms found per NH4

!Calculate and save all N-O distances
!Loop over all atoms
DO i = 1, n_atom

  !Search for Oxygen, type = 4, load coordinates, increment oxygen counter
  IF (Coords (5, 1, i) .EQ. 4) THEN

      counter = counter + 1
      x1 = Coords (1, 1, i)
      y1 = Coords (2, 1, i)
      z1 = Coords (3, 1, i)

      !Loop over all Nitrogen, calculate O-N distance, save dist and O atom id
      DO j = 1, n_ion

        x2 = Coords (1, 1, NH4 (j, 1))
        y2 = Coords (2, 1, NH4 (j, 1))
        z2 = Coords (3, 1, NH4 (j, 1))
        dx = x2 - x1
        dy = y2 - y1
        dz = z2 - z1
        dx = dx - bx * ANINT (dx / bx)
        dy = dy - by * ANINT (dy / by)
        dz = dz - bz * ANINT (dz / bz)
        O_Dist (j, 1, counter) = i !saving oxygen atom id
        O_Dist (j, 2, counter) = dx ** 2 + dy ** 2 + dz ** 2 !saving O distance

      END DO

    ELSE

      CYCLE

  END IF

END DO

!WRITE (*, *) 'Calculated Oxygen distances to Nitrogens'

!Sort O-N distances and find associated O atom ids
!Only save desired number of nearest neighboring O atoms "nn"
!Loop over all NH4 ions
DO i = 1, n_ion

  !Create temporary list to send to sorting subroutine
  Olist (:) = O_Dist (i, 2, :)
  CALL dpquicksort (Olist)

  !Only save O atom ids of the "nn" quantity of shortest N-O distances
  DO j = 1, nn

    !Initialize check and k
    check = .FALSE. !if TRUE, break out of do while loop
    k = 1 !atom id index, incremented in the do while loop

    !Search for oxygen atom corresponding to the jth shortest O-N distance
    DO WHILE (check .EQ. .FALSE.)

      !If jth shortest O-N distance found, save O atom id in neighborlist
      IF (Olist (j) .EQ. O_Dist (i, 2, k)) THEN

          check = .TRUE. !breaks out of do while loop
          Neigh_O (i, j) = O_dist (i, 1, k) !i=ion, j=O neigh, (i, 1, k)=O id

        ELSE

          !increment to next O-N distance comparison
          k = k + 1

      END IF

    END DO

  END DO

END DO

!WRITE (*, *) 'Created Neighborlists'

!No longer need temp Olist or full list of O-N distances of first snapshot
DEALLOCATE (Olist, O_Dist)
CLOSE (UNIT = 10, STATUS = 'KEEP') !Close trajectory file
OPEN (UNIT = 40, FILE = 'HBondstats.txt', STATUS = 'UNKNOWN') !Bond averages
OPEN (UNIT = 50, FILE = 'HAnglestats.txt', STATUS = 'UNKNOWN') !Angle averages

!WRITE (*, *) 'Starting Bond and Angle Calculations'

!Initialize arrays
Bonds (:, :) = 0.0
Angles (:, :) = 0.0
Bavg (:) = 0.0
Aavg (:) = 0.0
Bavgsort (:) = 0.0
Aavgsort (:) = 0.0

!Loop over all snapshots
DO i = 1, n_step

!  WRITE (*, '(A, I0)') 'snapshot=', i

  !Calculate box dimensions for snapshot
  bx = Box (2, 1, i) - Box (1, 1, i)
  by = Box (2, 2, i) - Box (1, 2, i)
  bz = Box (2, 3, i) - Box (1, 3, i)

  !Loop over all NH4 ions
  DO j = 1, n_ion

    !Nitrogen coordinates of the jth NH4 ion
    x1 = Coords (1, i, NH4 (j, 1))
    y1 = Coords (2, i, NH4 (j, 1))
    z1 = Coords (3, i, NH4 (j, 1))

    !Loop over the hydrogen atoms of the jth NH4 ion
    DO k = 2, 5

      !kth Hydrogen coordinates of the jth NH4 ion
      x2 = Coords (1, i, NH4 (j, k))
      y2 = Coords (2, i, NH4 (j, k))
      z2 = Coords (3, i, NH4 (j, k))
      mindist = 10.0 !initialize min. hydrogen bond dist (arb. large value)

      !Loop over the "nn" nearest oxygen neighbors of the jth NH4 ion
      !Attempt to find minimum O-H distance possible from provided neighlist
      DO l = 1, nn

        !Oxygen coordinates of nnth Oxygen in neighborlist for the jth NH4 ion
        !Calculate O-H distance
        x3 = Coords (1, i, Neigh_O (j, l))
        y3 = Coords (2, i, Neigh_O (j, l))
        z3 = Coords (3, i, Neigh_O (j, l))
        dx = x3 - x2
        dy = y3 - y2
        dz = z3 - z2
        dx = dx - bx * ANINT (dx / bx)
        dy = dy - by * ANINT (dy / by)
        dz = dz - bz * ANINT (dz / bz)
        dist = SQRT (dx ** 2 + dy ** 2 + dz ** 2)

        !Keep the information of the smallest O-H distance
        IF (dist .LT. mindist) THEN

            mindist = dist !update minimum O-H distance for kth H of jth NH4
            mino = l !save associated oxygen in neighborlist

          ELSE

            CYCLE

        END IF

      END DO

      !Add hydrogen bond information to histogram and average arrays
      ind = INT (mindist / db) + 1 !calculate index for bond histogram
      Bonds (k - 1, ind) = Bonds (k - 1, ind) + 1 !increment bond histogram
      Btemp1 (k - 1) = mindist !save hbond dist to temp array to be sorted

      !Oxygen coordinates of hydrogen bond + calculate angle
      x3 = Coords (1, i, Neigh_O (j, mino))
      y3 = Coords (2, i, Neigh_O (j, mino))
      z3 = Coords (3, i, Neigh_O (j, mino))
      vx1 = x2 - x1
      vy1 = y2 - y1
      vz1 = z2 - z1
      vx2 = x3 - x2
      vy2 = y3 - y2
      vz2 = z3 - z2
      vx1 = vx1 - bx * ANINT (vx1 / bx)
      vy1 = vy1 - by * ANINT (vy1 / by)
      vz1 = vz1 - bz * ANINT (vz1 / bz)
      vx2 = vx2 - bx * ANINT (vx2 / bx)
      vy2 = vy2 - by * ANINT (vy2 / by)
      vz2 = vz2 - bz * ANINT (vz2 / bz)
      v1mag = SQRT ((vx1 ** 2) + (vy1 ** 2) + (vz1 ** 2))
      v2mag = SQRT ((vx2 ** 2) + (vy2 ** 2) + (vz2 ** 2))
      vx1 = vx1 / v1mag
      vy1 = vy1 / v1mag
      vz1 = vz1 / v1mag
      vx2 = vx2 / v2mag
      vy2 = vy2 / v2mag
      vz2 = vz2 / v2mag
      dot = vx1 * vx2 + vy1 * vy2 + vz1 * vz2
      angle = 180.0 - (DACOS (dot)) * (180.0 / 3.141592654)

      !Add hydrogen bond angle to histogram and average arrays
      ind = INT (angle / da) + 1 !calculate index for angle histogram
      Angles (k - 1, ind) = Angles (k - 1, ind) + 1 !increment angle histogram
      Atemp1 (k - 1) = angle !save hbond angle to temp array to be sorted

    END DO

      !Create copy of hydrogen bond list for jth NH4 ion, sort, add to avg tot.
      Btemp2 (:) = Btemp1 (:) !copy hbond list
      CALL dpquicksort (Btemp2) !sort hbond list
      Bavgsort (:) = Bavgsort (:) + Btemp2 (:) !add hbond values for running total
      Bavg (:) = Bavg (:) + Btemp1 (:) !running total of unsorted hbonds

      !sort hydrogen bond angles to keep initial association with hbonds
      DO l = 1, 4

        !Compares original Hbond list to sorted Hbond list
        !Perform same sorting of the indicies for hydrogen bond angles
        IF (Btemp1 (l) .EQ. Btemp2 (1)) THEN

            Atemp2 (1) = Atemp1 (l)

          ELSE IF (Btemp1 (l) .EQ. Btemp2 (2)) THEN

            Atemp2 (2) = Atemp1 (l)

          ELSE IF (Btemp1 (l) .EQ. Btemp2 (3)) THEN

            Atemp2 (3) = Atemp1 (l)

          ELSE

            Atemp2 (4) = Atemp1 (l)

        END IF

      END DO

      !add hbond angle values to running total to be averaged at the end
      Aavg (:) = Aavg (:) + Atemp1 (:) !unsorted
      Aavgsort (:) = Aavgsort (:) + Atemp2 (:) !sorted

  END DO !end loop over jth NH4 ion

END DO !end loop over ith snapshort

!open file for Hbond histogram
!WRITE (filename1, '(A)') "HBondhist.txt"
!OPEN (UNIT = 20, FILE = filename1, STATUS = 'UNKNOWN')

!time average hbond histogram
!Bonds (:, :) = Bonds (:, :) / REAL (n_step)

!write hbond histogram to file
!DO i = b_start, b_end

!  WRITE (20, *) (i - 1) * db, Bonds (:, i)

!ENDDO

!CLOSE (UNIT = 20, STATUS = 'keep') !close Hbond histogram file

!open file for Hbond angle histogram
!WRITE (filename2, '(A)') "HAnglehist.txt"
!OPEN (UNIT = 30, FILE = filename2, STATUS = 'UNKNOWN')

!time average hbond angle histogram
!Angles (:, :) = Angles (:, :) / REAL (n_step)

!write hbond angle histogram to file
!DO i = a_start, a_end

!  WRITE (30, *) (i - 1) * da, Angles (:, i)

!ENDDO

!CLOSE (UNIT = 30, STATUS = 'KEEP') !close Hbond angle histogram file

!time and ion average the different hydrogen bonds and angles
Bavg (:) = Bavg (:) / REAL (n_step * n_ion) !unsorted bond averaging
Bavgsort (:) = Bavgsort (:) / REAL (n_step * n_ion) !sorted bond averaging
Aavg (:) = Aavg (:) / REAL (n_step * n_ion) !unsorted angle averaging
Aavgsort (:) = Aavgsort (:) / REAL (n_step * n_ion) !sorted angle averaging
WRITE (40, *) 'Unsorted', Bavg (:) !write unsorted average bonds
WRITE (40, *) 'Sorted', Bavgsort (:) !write sorted average bonds
WRITE (50, *) 'Unsorted', Aavg (:) !write unsorted average angles
WRITE (50, *) 'Sorted', Aavgsort (:) !write sorted average angles

CLOSE (UNIT = 40, STATUS = 'KEEP') !close average hbond file
CLOSE (UNIT = 50, STATUS = 'KEEP') !close average hbond angle file

!subroutine included for sorting algorithm
CONTAINS

  ! ripped from https://www.mjr19.org.uk/IT/sorts/sorts.f90
  ! dual pivot quicksort
  recursive subroutine dpquicksort(array)
    real(prec), intent(inout)::array(:)
    real(prec) :: temp,p1,p2
    integer :: i,j,last,l,k,g

    last=size(array)

    if (last.lt.40) then ! use insertion sort on small arrays
       do i=2,last
          temp=array(i)
          do j=i-1,1,-1
             if (array(j).le.temp) exit
             array(j+1)=array(j)
          enddo
          array(j+1)=temp
       enddo
       return
    endif
    p1=array(last/3)
    p2=array(2*last/3)
    if (p2.lt.p1) then
       temp=p1
       p1=p2
       p2=temp
    endif
    array(last/3)=array(1)
    array(1)=p1
    array(2*last/3)=array(last)
    array(last)=p2

    g=last
    l=2
    do while (array(l).lt.p1)
       l=l+1
    enddo
    k=l

    do while(k.lt.g)
       temp=array(k)
       if (temp.lt.p1) then
          array(k)=array(l)
          array(l)=temp
          l=l+1
       else if (temp.gt.p2) then
          do while(array(g-1).gt.p2)
             g=g-1
          enddo
          if (k.ge.g) exit
          g=g-1
          if (array(g).lt.p1) then
             array(k)=array(l)
             array(l)=array(g)
             array(g)=temp
             l=l+1
          else
             array(k)=array(g)
             array(g)=temp
          endif
       endif
       k=k+1
    enddo
    if (l.gt.2) then
       array(1)=array(l-1)
       array(l-1)=p1
       call dpquicksort(array(1:l-2))
    endif
    call dpquicksort(array(l:g-1))
    if (g.lt.last) then
       array(last)=array(g)
       array(g)=p2
       call dpquicksort(array(g+1:last))
    endif

  end subroutine dpquicksort

END PROGRAM HBOND

max_stars_repo_name: dowlinglab/ap-fffit

max_issues_repo_name: dowlinglab/ap-fffit

max_forks_repo_name: dowlinglab/ap-fffit


text:
      program mdslib_test

      implicit none

      integer descr, MdsValue, MdsOpen, MdsPut

      integer IDTYPE_LONG, IDTYPE_FLOAT, IDTYPE_CSTRING

      parameter (IDTYPE_LONG=8, IDTYPE_FLOAT=10, IDTYPE_CSTRING=14)

      integer dsc
      integer sts
      integer i,j
      integer iresult
      integer nx,ny
      integer size
      parameter (nx=2, ny=13)
      real    result
      real    resultarr(20)
      real    array2d(nx,ny)
      character cresult*12

      do i=1,nx
         do j=1,ny
            array2d(i,j) = i*10 + j
         enddo
      enddo

      dsc = descr(IDTYPE_FLOAT,result,0,0)
      sts = MdsValue("1.23"//CHAR(0),dsc,0,size)
      if (abs(result-1.23) .gt. 1e-5) then
        write (6,*) "MdsValue('1.23') : ",result,sts
        stop 1
      end if

      dsc = descr(IDTYPE_FLOAT,resultarr,5,0)
      sts = MdsValue("2.2 : 8. : 2."//CHAR(0),dsc,0,size)
      if (abs(resultarr(1)-2.2) .gt. 1e-5) then
        write (6,*) "MdsValue('2. : 40. : 2.') : ",resultarr,sts,size
        stop 1
      end if

      dsc = descr(IDTYPE_LONG,iresult,0)
      sts = MdsValue("_=SetEnv('test_path=.')"//CHAR(0),dsc,0,size)
      if (and(and(sts, iresult), 1) .ne. 1) then
        write (6,*) "setenv(): ",iresult,sts,size
        stop 1
      end if
      sts = MdsValue("TreeOpenNew('test',1)"//CHAR(0),dsc,0,size)
      if (and(and(sts, iresult), 1) .ne. 1) then
        write (6,*) "TreeOpenNew(): ",iresult,sts,size
        stop 1
      end if
      sts = MdsValue("TreeAddNode('A',_,'ANY')"//CHAR(0),dsc,0,size)
      if (and(and(sts, iresult), 1) .ne. 1) then
        write (6,*) "TreeAddNode: ",iresult,sts,size
        stop 1
      end if
      sts = MdsValue("TreeWrite()"//CHAR(0),dsc,0,size)
      if (and(and(sts, iresult), 1) .ne. 1) then
        write (6,*) "TreeWrite(): ",iresult,sts,size
        stop 1
      end if
      sts = MdsValue("TreeClose()"//CHAR(0),dsc,0,size)
      if (and(and(sts, iresult), 1) .ne. 1) then
        write (6,*) "TreeClose(): ",iresult,sts,size
        stop 1
      end if

      dsc = descr(IDTYPE_CSTRING,cresult,0,12)
      sts = MdsValue("$EXPT"//CHAR(0),dsc,0,size)
      if (and(sts, 1) .eq. 1) then
        write (6,*) "MdsValue($EXPT)/*before*/: ",cresult,sts,size
        stop 1
      end if

      sts = MdsOpen("test"//CHAR(0),1)
      if (and(sts, 1) .ne. 1) then
        write (6,*) "MdsOpen('test',1) :", sts
        stop 1
      end if

      dsc = descr(IDTYPE_CSTRING,cresult,0,12)
      sts = MdsValue("$EXPT"//CHAR(0),dsc,0,size)
      if (and(sts, 1) .ne. 1) then
        write (6,*) "MdsValue($EXPT)/*after*/: ",cresult,sts,size
        stop 1
      end if

      dsc = descr(IDTYPE_LONG,42042,0)
      sts = MdsPut("A"//CHAR(0),"$"//CHAR(0),dsc,0)
      if (and(sts, 1) .ne. 1) then
        write (6,*) "MdsPut('A',$,42042): ",sts
        stop 1
      end if

      dsc = descr(IDTYPE_LONG,iresult,0)
      sts = MdsValue("A"//CHAR(0),dsc,0,size)
      if (and(sts, 1) .ne. 1) then
        write (6,*) "MdsValue('A')/*i*/: ",iresult,sts
        stop 1
      end if

      dsc = descr(IDTYPE_FLOAT,resultarr,20,0)
      do i=1,20
         resultarr(i) = 0.
      enddo
      sts = MdsPut("A"//CHAR(0),"$"//CHAR(0),dsc,0)
      if (and(sts, 1) .ne. 1) then
        write (6,*) "MdsPut('A','$',array) : ",sts
        stop 1
      end if

      sts = MdsValue("A"//CHAR(0),dsc,0,size)
      if (and(sts, 1) .ne. 1) then
        write (6,*) "MdsValue('A')/*/array*/: ",resultarr,sts,size
        stop 1
      end if

      dsc = descr(IDTYPE_FLOAT,array2d,nx,ny,0)
      sts = MdsPut("A"//CHAR(0),"$"//CHAR(0),dsc,0)
      if (and(sts, 1) .ne. 1) then
        write (6,*) "MdsPut('A',$,array2d) : ",sts
        stop 1
      end if

      sts = MdsValue("A"//CHAR(0),dsc,0,size)
      if (and(sts, 1) .ne. 1) then
        write (6,*) "MdsValue('A')/*/array2d*/: ",resultarr,sts,size
        stop 1
      end if

      end

max_stars_repo_name: fluffynukeit/mdsplus

max_issues_repo_name: fluffynukeit/mdsplus

max_forks_repo_name: fluffynukeit/mdsplus


text:
! Fortran 95 version

function f95add(a, b) result(ret)
  implicit none
  integer :: a
  integer :: b
  integer :: ret
  ret = a + b
end function f95add

max_stars_repo_name: plicease/FFI-Platypus-Lang-Fortran

max_issues_repo_name: PerlFFI/FFI-Platypus-Lang-Fortran


!--------+--------+--------+--------+--------+--------+--------+------!
! main program of a set of tools based on Wannier90 TB
! constructed by Q.S.Wu on 4/9/2010
! change      by Q.S.Wu on 4/22/2010
! changed     by Q.S.wu on July/15/2010
! version     HmnR.data  contains soc
! mpi-version is not test yet, take you own risk.
! mpi-version is  tested , please report bugs to QSWU
! Jan 25 2015 by Q.S.Wu at ETH Zurich
! wuquansheng@gmail.com
! Copyright (c) 2010 QuanSheng Wu. All rights reserved.
!--------+--------+--------+--------+--------+--------+--------+------!

  program main

     use wmpi
     use para
     implicit none

     !> file existence
     logical  :: exists
     integer :: ierr
     character(8) :: cht

     !> time measure
     real(dp) :: time_start, time_end, time_init

     ierr = 0
     cpuid= 0
     num_cpu= 1
     !> initial the environment of mpi
#if defined (MPI)
     call mpi_init(ierr)
     call mpi_comm_rank(mpi_cmw,cpuid,ierr)
     call mpi_comm_size(mpi_cmw,num_cpu,ierr)
#endif

     if (cpuid==0) open(unit=stdout, file='WT.out')

     !> if mpi initial wrong, alarm
     if (cpuid==0.and.ierr.ne.0)then
        write(stdout,*)'mpi initialize wrong'
        stop
     endif

     call now(time_init)
     call header

     !> print information for mpi
     if (cpuid==0) then
        write(stdout, '(1x, a, i5, a)')'You are using ', num_cpu, ' CPU cores'
        write(stdout, *)' '
     endif

     !> readin the control parameters for this program
     call readinput

     !> open file Hmn_R.data to get Num_wann and Nrpts
     if (cpuid.eq.0)then
        write(stdout,*)''
        inquire (file =Hrfile, EXIST = exists)
        if (exists)then
           if (index(Hrfile, 'HWR')==0) then
              write(stdout,'(2x,a,a,a)')'File ',trim(Hrfile), &
                 ' exist, We are using HmnR from wannier90'
              open(unit=1001,file=Hrfile,status='old')
              read(1001,*)
              read(1001,*) Num_wann
              read(1001,*) Nrpts
              write(stdout,*)'>> Num_wann', Num_wann
              write(stdout,*)'>> NRPTS', NRPTS
              close(1001)
           else
              write(stdout,'(2x, 3a)')'File ',trim(Hrfile), &
                 ' exist, We are using HmnR from HWR'
              open(unit=1001,file=Hrfile,status='old')
              read(1001,*)
              read(1001,'(a26, i3)')cht, Num_wann
              read(1001,'(a32,i10)')cht,Nrpts
              write(stdout,*)'>> Num_wann', Num_wann
              write(stdout,*)'>> NRPTS', NRPTS
              close(1001)
           endif ! hwr or not
        else
           write(stdout,'(2x,a25)')'>>> Error : no HmnR input'
           stop
        endif ! exists or not
        if ((soc==0 .and. sum(nprojs)/=Num_wann) .or. &
           (soc>0 .and. sum(nprojs)/=Num_wann/2))then
           print *, sum(nprojs), num_wann, num_wann/2
           stop 'projectors are wrong'
        endif
     endif ! cpuid

     !> broadcast and Nrpts to every cpu
#if defined (MPI)
     call MPI_bcast(Num_wann,1,mpi_in,0,mpi_cmw,ierr)
     call MPI_bcast(Nrpts,1,mpi_in,0,mpi_cmw,ierr)
#endif
     !> dimension for surface green's function
     Ndim= Num_wann* Np

     !> allocate necessary arrays for tight binding hamiltonians
     allocate(irvec(3,nrpts))
     allocate(ndegen(nrpts))
     allocate(HmnR(num_wann,num_wann,nrpts))

     if (cpuid==0)then
        write(stdout,*) ' >> Begin to read Hmn_R.data'
     endif
     call readHmnR()
     if (cpuid==0)then
        write(stdout,*) ' << Read Hmn_R.data successfully'
     endif

     !> broadcast data to every cpu
#if defined (MPI)
     call MPI_bcast(irvec,size(irvec),mpi_in,0,mpi_cmw,ierr)
     call MPI_bcast(HmnR,size(HmnR),mpi_dc,0,mpi_cmw,ierr)
     call MPI_bcast(ndegen,size(ndegen),mpi_in,0,mpi_cmw,ierr)
#endif

     !> import symmetry
    !call symmetry

     !> bulk band
     if (BulkBand_calc) then
        if(cpuid.eq.0)write(stdout, *)' '
        if(cpuid.eq.0)write(stdout, *)'>> Start of calculating bulk band'
        call now(time_start)
        call ek_bulk
       !call ek_bulk_mirror_x
       !call ek_bulk_mirror_z
        call now(time_end)
        call print_time_cost(time_start, time_end, 'BulkBand')
        if(cpuid.eq.0)write(stdout, *)'<< End of calculating bulk band'
     endif

     if (BulkFS_calc) then
        if(cpuid.eq.0)write(stdout, *)' '
        if(cpuid.eq.0)write(stdout, *)'>> Start of calculating the bulk FS'
        call now(time_start)
       !call fermisurface3D
        call fermisurface
        call now(time_end)
        call print_time_cost(time_start, time_end, 'BulkFS')
        if(cpuid.eq.0)write(stdout, *)'<< End of calculating the bulk FS'
     endif

     if (JDos_calc.and.Dos_calc) then
        if(cpuid.eq.0)write(stdout, *)' '
        if(cpuid.eq.0)write(stdout, *)'>> Start of calculating DOS and Jdos for bulk system'
        call now(time_start)
        call dos_joint_dos
        call now(time_end)
        call print_time_cost(time_start, time_end, 'Dos_calc and Jdos_calc')
        if(cpuid.eq.0)write(stdout, *)'<< End of calculating the DOS and Jdos for bulk system'
     else
        if (Dos_calc) then
           if(cpuid.eq.0)write(stdout, *)' '
           if(cpuid.eq.0)write(stdout, *)'>> Start of calculating DOS for bulk system'
           call now(time_start)
           call dos_sub
           call now(time_end)
           call print_time_cost(time_start, time_end, 'Dos_calc')
           if(cpuid.eq.0)write(stdout, *)'<< End of calculating the DOS for bulk system'
        endif

        if (JDos_calc) then
           if(cpuid.eq.0)write(stdout, *)' '
           if(cpuid.eq.0)write(stdout, *)'>> Start of calculating JDOS for bulk system'
           call now(time_start)
           call Joint_dos
           call now(time_end)
           call print_time_cost(time_start, time_end, 'JDos_calc')
           if(cpuid.eq.0)write(stdout, *)'<< End of calculating the JDOS for bulk system'
        endif
     endif

     !> effective mass
     if (EffectiveMass_calc) then
        if(cpuid.eq.0)write(stdout, *)' '
        if(cpuid.eq.0)write(stdout, *)'>> Start of calculating the effective mass'
        call now(time_start)
        call effective_mass_calc
        call now(time_end)
        call print_time_cost(time_start, time_end, 'EffectiveMass_calc')
        if(cpuid.eq.0)write(stdout, *)'<< End of calculating the effective mass'
     endif

     if (BulkGap_plane_calc) then
        if(cpuid.eq.0)write(stdout, *)' '
        if(cpuid.eq.0)write(stdout, *)'>> Start of calculating the bulk gap in plane'
        call now(time_start)
       !call psik_bulk
       !call ek_bulk_polar
       !call ek_bulk_fortomas
       !call ek_bulk_spin
       !call dos_calc
       !call ek_bulk2D
       !call ek_bulk2D_spin
        call gapshape
        call now(time_end)
        call print_time_cost(time_start, time_end, 'BulkGap_plane')
        if(cpuid.eq.0)write(stdout, *)'<< End of calculating the bulk gap in plane'
     endif

     if (BulkGap_Cube_calc) then
        if(cpuid.eq.0)write(stdout, *)' '
        if(cpuid.eq.0)write(stdout, *)'>> start of calculating the bulk gap in Cube'
        call now(time_start)
        call gapshape3D
       !call landau_level_k
       !call landau_level_B
        call now(time_end)
        call print_time_cost(time_start, time_end, 'BulkGap_Cube')
        if(cpuid.eq.0)write(stdout, *)'<< End of calculating the bulk gap in Cube'
     endif

     !> slab band
     if (SlabBand_calc)then
        if(cpuid.eq.0)write(stdout, *)' '
        if(cpuid.eq.0)write(stdout, *)'>> Start of calculating the slab band structure'
        call now(time_start)
        call ek_slab
       !call psik
        call now(time_end)
        call print_time_cost(time_start, time_end, 'SlabBand')
        if(cpuid.eq.0)write(stdout, *)'<< End of calculating the slab band structure'
     endif

     if (BerryCurvature_calc)then
        if(cpuid.eq.0)write(stdout, *)' '
        if(cpuid.eq.0)write(stdout, *)'>> Start of calculating the Berry curvature'
        call now(time_start)
        call berry_curvarture
        call now(time_end)
        call print_time_cost(time_start, time_end, 'BerryCurvature')
        if(cpuid.eq.0)write(stdout, *)'End of calculating the Berry curvature'
     endif

     if (WireBand_calc) then
        if(cpuid.eq.0)write(stdout, *)' '
        if(cpuid.eq.0)write(stdout, *)'>> Start of calculating the wire band'
        call now(time_start)
        call ek_ribbon
        call now(time_end)
        call print_time_cost(time_start, time_end, 'WireBand')
        if(cpuid.eq.0)write(stdout, *)'End of calculating the wire band'
     endif

     !> wannier center calculate
     if (wanniercenter_calc)then
        if(cpuid.eq.0)write(stdout, *)' '
        if(cpuid.eq.0)write(stdout, *)'>> Start of calculating the Wilson loop'
        call now(time_start)
        call wannier_center3D_plane
       !call wannier_center2D
       !call wannier_center2D_alt
       !call wannier_center3D
       !call wannier_center3D_plane_mirror_plus
       !call wannier_center3D_plane_mirror_minus
        call now(time_end)
        call print_time_cost(time_start, time_end, 'WannierCenter')
        if(cpuid.eq.0)write(stdout, *)'<< End of calculating the Wilson loop'
     endif

     if (BerryPhase_calc) then
        if(cpuid.eq.0)write(stdout, *)' '
        if(cpuid.eq.0)write(stdout, *)'>> Start of calculating the Berry phase'
        call now(time_start)
        call berryphase
        call now(time_end)
        call print_time_cost(time_start, time_end, 'BerryPhase')
        if(cpuid.eq.0)write(stdout, *)'End of calculating the Berry phase'
     endif

     !> surface state
     if (SlabSS_calc) then
        if(cpuid.eq.0)write(stdout, *)' '
        if(cpuid.eq.0)write(stdout, *)'>> Start of calculating the surface state'
        call now(time_start)
        call surfstat
        call now(time_end)
        call print_time_cost(time_start, time_end, 'BerryCurvature')
        if(cpuid.eq.0)write(stdout, *)'End of calculating the surface state'
     endif

     !> fermi arc
     if (SlabArc_calc) then
        if(cpuid.eq.0)write(stdout, *)' '
        if(cpuid.eq.0)write(stdout, *)'>> Start of calculating the surface arc'
        call now(time_start)
        call fermiarc
        call now(time_end)
        call print_time_cost(time_start, time_end, 'SlabArc')
        if(cpuid.eq.0)write(stdout, *)'End of calculating the surface arc'
     endif

     !> fermi arc QPI
     if (SlabQPI_calc) then
        if(cpuid.eq.0)write(stdout, *)' '
        if(cpuid.eq.0)write(stdout, *)'>> Start of calculating the surface QPI'
        call now(time_start)
        call fermiarc_jdos
        call now(time_end)
        call print_time_cost(time_start, time_end, 'SlabQPI')
        if(cpuid.eq.0)write(stdout, *)'End of calculating the surface QPI'
     endif

     !> calculate spin-texture
     if (SlabSpintexture_calc)then
        if(cpuid.eq.0)write(stdout, *)' '
        if(cpuid.eq.0)write(stdout, *)'>> Start of calculating the spin texture for surface'
        call now(time_start)
        call spintext
        call now(time_end)
        call print_time_cost(time_start, time_end, 'SlabSpintexture')
        if(cpuid.eq.0)write(stdout, *)'End of calculating the spin texture for surface'
     endif

     call now(time_end)

     if(cpuid.eq.0)write(stdout, *)' '
     call print_time_cost(time_init, time_end, 'whole program')
     if (cpuid.eq.0)write(stdout,*)'Congratulations! you finished the calculation.'
     if (cpuid.eq.0)write(stdout,*)'See you next time :)'

#if defined (MPI)
     call mpi_finalize(ierr)
#endif

  end  !<< end of main program


text:
      SUBROUTINE QLOADL (IOPT)
C
C     THIS ROUTINE CALCULATES THERMAL LOADS FROM QBDY1, QBDY2, OR
C     QVECT DATA. THE INPUT DATA, READ FROM FILE SLT, IS -
C
C     ENTRY       QBDY1         QBDY2          QVECT
C     -----       -----         -----          -----
C       1          TYPE         EL.ID.          SIL1
C       2         EL.ID.         TYPE           SIL2
C       3          SIL1          SIL1           SIL3
C       4          SIL2          SIL2           SIL4
C       5          SIL3          SIL3          EL.ID.
C       6          SIL4          SIL4           TYPE
C      7-10      C1,C2,C3,C4    -SAME          -SAME
C     11-13       E VECTOR       NONE           NONE
C     14-16      V1 VECTOR         *             *
C     17-19      V2 VECTOR         *             *
C
      LOGICAL         NOGO,TRANST
      INTEGER         SLT,OLD,BG,SUBR(2),IGRIDS(6),TYPE,SILS(4),IE(3),
     1                MINUS(2)
      REAL            E(3),COEF(4),V1(3),V2(3),CARD(19)
      CHARACTER       UFM*23
      COMMON /XMSSG / UFM
      COMMON /LOADX / LC,SLT,BG,OLD,NXN(12),IFM,NYN(2),ILID
      COMMON /ZZZZZZ/ CORE(1)
      COMMON /SYSTEM/ SYSBUF,IOUT
      COMMON /QVECT / ITRAN,IQVECT
      EQUIVALENCE     (TYPE,CARD(1)),(ID,CARD(2)),(SILS(1),CARD(3)),
     1                (COEF(1),CARD(7)),(E(1),IE(1),CARD(11)),
     2                (V1(1),CARD(14)),(V2(1),CARD(17))
      DATA    IGRIDS/ 1,2,2,3,4,2   /
      DATA    SUBR  / 4HQLOA,4HDL   /
      DATA    ITRAN1, IOLD,MINUS    / 4HTRAN,0,-1,-1 /
C
      TRANST = .FALSE.
      IF (ITRAN .EQ. ITRAN1) TRANST = .TRUE.
      NWORDS = 10
      IF (IOPT .EQ. 3) NWORDS = 19
C
      CALL READ (*100,*110,SLT,CARD(1),NWORDS,0,FLAG)
C
C     REARRANGE CARD ARRAY FOR UNIFORMITY.
C
      GO TO (20,10,40), IOPT
   10 DOT = CARD(1)
      CARD(1) = CARD(2)
      CARD(2) = DOT
   20 N = IGRIDS(TYPE)
C
C     QBDY1 OR QBDY2
C
      DO 30 I = 1,N
      ISIL = SILS(I)
      CORE(ISIL) = CORE(ISIL) + COEF(I)
   30 CONTINUE
      RETURN
C
C     QVECT LOADS
C
   40 DOT     = CARD(5)
      DOT2    = CARD(6)
      CARD(6) = CARD(4)
      CARD(5) = CARD(3)
      CARD(4) = CARD(2)
      CARD(3) = CARD(1)
      CARD(2) = DOT
      CARD(1) = DOT2
      N       = IGRIDS(TYPE)
      DOT     = 0.0
      INT     = 0
      IF (TYPE .EQ. 6) GO TO 70
      DO 50 I = 1,3
      IF (NUMTYP(IE(I)) .EQ. 1) GO TO 51
      DOT = DOT + E(I)*V1(I)
      GO TO 50
   51 INT = INT + 1
   50 CONTINUE
      IF (INT .GT.   0) GO TO 90
      IF (DOT .GE. 0.0) RETURN
      DO 60 I = 1,N
      ISIL = SILS(I)
   60 CORE(ISIL) = CORE(ISIL) - DOT*COEF(I)
      RETURN
C
C     QVECT ON ELCYL ELEMENT
C
   70 DOT2 = 0.0
      DO 80 I = 1,3
      IF (NUMTYP(IE(I)) .EQ. 1) GO TO 81
      DOT  = DOT  + E(I)*V1(I)
      DOT2 = DOT2 + E(I)*V2(I)
      GO TO 80
   81 INT = INT + 1
   80 CONTINUE
      IF (INT .GT. 0) GO TO 90
      COEF(1) = COEF(1)*SQRT(DOT**2 + DOT2**2)
      COEF(2) = COEF(1)
      ISIL = SILS(1)
      CORE(ISIL) = CORE(ISIL) + COEF(1)
      ISIL = SILS(2)
      CORE(ISIL) = CORE(ISIL) + COEF(2)
      RETURN
C
C     GOES HERE IF INTEGERS ARE FOUND IN E VECTOR
C
   90 IF (.NOT. TRANST) GO TO 120
C
C     BUILD QVECT RECORDS FOR TRANSIENT
C
      IF (ILID .EQ. IOLD) GO TO 91
      IF (IOLD .EQ.    0) GO TO 92
C
C     TERMINATE OLD RECORD
C
      CALL WRITE (IQVECT,MINUS,2,0)
   92 IOLD = ILID
      CALL WRITE (IQVECT,ILID,1,0)
C
C     DUMP DATA ON IQVECT
C
   91 CALL WRITE (IQVECT,N,1,0)
      DO 93 I = 1,N
      CALL WRITE (IQVECT,SILS(I),1,0)
      CALL WRITE (IQVECT,COEF(I),1,0)
   93 CONTINUE
      CALL WRITE (IQVECT,IE,3,0)
      CALL WRITE (IQVECT,V1,6,0)
      RETURN
C
  100 CALL MESAGE (-2,SLT,SUBR)
  110 CALL MESAGE (-3,SLT,SUBR)
  120 NOGO = .TRUE.
      WRITE  (IOUT,130) UFM,ID
  130 FORMAT (A23,' 3080, ERROR IN QVECT DATA, INTEGER VALUES SPECIFIED'
     1,      ' FOR THERMAL FLUX VECTOR COMPONENTS', /30X,
     2       'IN A NON-TRANSIENT ANALYSIS.', /30X,'ELEMENT ID = ',I9)
      CALL MESAGE (-61,0,SUBR)
      RETURN
      END


text:
      subroutine compute_nxyz(ijob)
C***********************************************************************
!
      use comparttr
      use comparttr_sg
      use comunits, only : error_unit_number
      use comgrid,    only : n_grid_points
      use comgrid_sg, only : n_grid_points_sg
!
      implicit none
C
C     Variables passed into subroutine
C
      integer ijob
C
C     ijob = 0 : Assume nxyz was filled by reading from input.
C                Only compute the number of grid points and
C                verify valid (.ge. 1) values for nxyz.
C
C     ijob = 1 : divide cells into N equal increments so that
C                each increment is less than or equal to the
C                specified length scale.
C
C     ijob = 2 : same as one with the added constraint that N
C                must be a power of ibase = 2. If the value of
C                ibase is modified other powers can be easily
C                implemented.
C
C
C     Subroutine variables
C
      integer i,j,k, lu, icounter
      integer n_power, ibase
      integer nearest_power_int
      integer n_warning_value, if_error
!
      integer local_debug
      data    local_debug / 0 /
!
      n_warning_value = 32
      ibase = 2
C
C***********************************************************************

      if(ijob .gt. 0)then
C
C     Compute number of subdivisions (nxyz), based on
C     cell size (dxyz) and desired length scale (scale_xyz)
C
      icounter = 0
      do i = 1, n_touched_cells
         j = touched_cells(i)
         do k = 1,3
            nxyz(i,k) =
     1                   int(dxyz(j,k)/scale_xyz(i,k) + 1.0d0 -
     2                 1.d-8*dxyz(j,k)/scale_xyz(i,k))
         enddo
      enddo
      endif
C
C     Promote integers to nearest larger power of 2
C
      if(ijob .eq. 2)then
      do i = 1, n_touched_cells
         do k = 1,3
           n_power = nearest_power_int(ibase,nxyz(i,k))
           j = ibase**n_power
           if(j .lt. nxyz(i,k)) then
              nxyz(i,k) = j*ibase
           else
              nxyz(i,k) = j
           endif
         enddo
      enddo
      endif
C
C     Print error if nx,ny,nz is less than 1
C
      if_error = 0
      do i = 1, n_touched_cells
         if((nxyz(i,1).lt. 1).or.
     1      (nxyz(i,2).lt. 1).or.
     2      (nxyz(i,3).lt. 1))then
            if_error = if_error + 1
            write(error_unit_number,*)
     1       'ERROR Subgrid Cell Refinement ERROR'
            write(error_unit_number,100)
     1       dxyz(i,1),scale_xyz(i,1), nxyz(i,1)
            write(error_unit_number,110)
     1       dxyz(i,2),scale_xyz(i,2), nxyz(i,2)
            write(error_unit_number,120)
     1       dxyz(i,3),scale_xyz(i,3), nxyz(i,3)
          endif
          if(if_error .ne. 0)then
            write(error_unit_number,*)
     1       'ERROR-STOP in compute_nxyz-ERROR'
            stop
          endif
      enddo
C
C     Compute the number of grid points in the subgrid.
C
      n_grid_points_sg = 0
      do i = 1, n_touched_cells
         n_grid_points_sg = n_grid_points_sg +
     1                      nxyz(i,1)*nxyz(i,2)*nxyz(i,3)
      enddo
      n_grid_points_sg =
     1  n_grid_points_sg + n_grid_points - n_touched_cells
C
C     Print warnings if nx,ny,nz is larger than n_warning_value
C
      do i = 1, n_touched_cells
         if((nxyz(i,1).gt. n_warning_value).or.
     1      (nxyz(i,2).gt. n_warning_value).or.
     2      (nxyz(i,3).gt. n_warning_value))then
            write(error_unit_number,*)
     1       'Subgrid Cell Refinement WARNING'
            write(error_unit_number,*)
     1       'Input will result in a large number of '
     2           ,'subgrid refined cells.'
            write(error_unit_number,100)
     1       dxyz(i,1),scale_xyz(i,1), nxyz(i,1)
            write(error_unit_number,110)
     1       dxyz(i,2),scale_xyz(i,2), nxyz(i,2)
            write(error_unit_number,120)
     1       dxyz(i,3),scale_xyz(i,3), nxyz(i,3)
  100       format
     1  ('cell dx =',e14.7,' scale factor =',e14.7,' sg nx=',i6)
  110       format
     1  ('cell dy =',e14.7,' scale factor =',e14.7,' sg ny=',i6)
  120       format
     1  ('cell dz =',e14.7,' scale factor =',e14.7,' sg nz=',i6)
          endif
         enddo
!
!     Output when local_debug is not 0
!
      if(local_debug .ne. 0)then

      write(error_unit_number, *) 'compute_nxyz: local_debug .ne. 0'
      write(error_unit_number, *)
     1       'index, cell_id, nxyz(x), dxyz(x), scale_xyz(x)...y...z'
      do i = 1, n_touched_cells
         j = touched_cells(i)
         write(error_unit_number,900)i,j,
     1      nxyz(i,1),dxyz(j,1),scale_xyz(i,1),
     2      nxyz(i,2),dxyz(j,2),scale_xyz(i,2),
     3      nxyz(i,3),dxyz(j,3),scale_xyz(i,3)
 900        format(i8,1x,i8,1x,3(i8,1x,e14.6,1x,e14.6))
      enddo
      write(error_unit_number, *) 'compute_nxyz: local_debug .ne. 0'

      endif

      return
      end
      function nearest_power_int(ibase,int_value)
      integer ibase, int_value
      real*8 base, real_value

      base = float(ibase)
      real_value = float(int_value)
      nearest_power_int  = nint(log10(real_value)/log10(base))

      return
      end

max_stars_repo_name: lanl/PLUMECALC

max_issues_repo_name: lanl/PLUMECALC

max_forks_repo_name: lanl/PLUMECALC


        !COMPILER-GENERATED INTERFACE MODULE: Fri May 22 16:58:43 2020
        ! This source file is for reference only and may not completely
        ! represent the generated interface used by the compiler.
        MODULE COMPUTE_VECTOR_WHOLE_MEDIUM__genmod
          INTERFACE
            SUBROUTINE COMPUTE_VECTOR_WHOLE_MEDIUM(FIELD_ACOUSTIC,      &
     &FIELD_GRAVITOACOUSTIC,FIELD_GRAVITO,FIELD_ELASTIC,                &
     &FIELDS_POROELASTIC)
              USE SPECFEM_PAR
              REAL(KIND=4) :: FIELD_ACOUSTIC(NGLOB_ACOUSTIC)
              REAL(KIND=4) :: FIELD_GRAVITOACOUSTIC(                    &
     &NGLOB_GRAVITOACOUSTIC)
              REAL(KIND=4) :: FIELD_GRAVITO(NGLOB_GRAVITOACOUSTIC)
              REAL(KIND=4) :: FIELD_ELASTIC(3,NGLOB_ELASTIC)
              REAL(KIND=4) :: FIELDS_POROELASTIC(2,NGLOB_POROELASTIC)
            END SUBROUTINE COMPUTE_VECTOR_WHOLE_MEDIUM
          END INTERFACE
        END MODULE COMPUTE_VECTOR_WHOLE_MEDIUM__genmod