!
!     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!     *                                                               *
!     *                  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_analysis_routines) vector_analysis_regular_grid_saved

contains

    ! Purpose:
    !
    !     subroutine vhaes(nlat, nlon, ityp, nt, v, w, idvw, jdvw, br, bi, cr, ci, &
    !                      mdab, ndab, wvhaes, lvhaes, work, lwork, ierror)
    !
    !     subroutine vhaes performs the vector spherical harmonic analysis
    !     on the vector field (v, w) and stores the result in the arrays
    !     br, bi, cr, and ci. v(i, j) and w(i, j) are the colatitudinal
    !     (measured from the north pole) and east longitudinal components
    !     respectively, located at colatitude theta(i) = (i-1)*pi/(nlat-1)
    !     and longitude phi(j) = (j-1)*2*pi/nlon. the spectral
    !     representation of (v, w) is given at output parameters v, w in
    !     subroutine vhses.
    !
    !     input parameters
    !
    !     nlat   the number of colatitudes on the full sphere including the
    !            poles. for example, nlat = 37 for a five degree grid.
    !            nlat determines the grid increment in colatitude as
    !            pi/(nlat-1).  if nlat is odd the equator is located at
    !            grid point i=(nlat + 1)/2. if nlat is even the equator is
    !            located half way between points i=nlat/2 and i=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 analysis
    !                 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 analysis
    !                 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 analysis
    !                 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 analysis is performed on the northern
    !                 hemisphere only.  i.e., if nlat is odd the analysis
    !                 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 analysis 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 analysis is performed on the northern
    !                 hemisphere only.  i.e., if nlat is odd the analysis
    !                 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 analysis 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 analysis is performed on the northern
    !                 hemisphere only.  i.e., if nlat is odd the analysis
    !                 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 analysis 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 analysis is performed on the northern
    !                 hemisphere only.  i.e., if nlat is odd the analysis
    !                 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 analysis 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 analysis is performed on the northern
    !                 hemisphere only.  i.e., if nlat is odd the analysis
    !                 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 analysis 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 analysis is performed on the northern
    !                 hemisphere only.  i.e., if nlat is odd the analysis
    !                 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 analysis 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 analyses.  in the program that calls vhaes,
    !            the arrays v, w, br, bi, cr, and ci can be three dimensional
    !            in which case multiple analyses will be performed.
    !            the third index is the analysis index which assumes the
    !            values k=1, ..., nt.  for a single analysis set nt=1. the
    !            discription of the remaining parameters is simplified
    !            by assuming that nt=1 or that all the arrays are two
    !            dimensional.
    !
    !     v, w    two or three dimensional arrays (see input parameter nt)
    !            that contain the vector function to be analyzed.
    !            v is the colatitudnal component and w is the east
    !            longitudinal component. v(i, j), w(i, j) contain the
    !            components at colatitude theta(i) = (i-1)*pi/(nlat-1)
    !            and longitude phi(j) = (j-1)*2*pi/nlon. the index ranges
    !            are defined above at the input parameter ityp.
    !
    !     idvw   the first dimension of the arrays v, w as it appears in
    !            the program that calls vhaes. 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 vhaes. jdvw must be at least nlon.
    !
    !     mdab   the first dimension of the arrays br, bi, cr, and ci as it
    !            appears in the program that calls vhaes. 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 vhaes. ndab must be at
    !            least nlat.
    !
    !     lvhaes an array which must be initialized by subroutine vhaesi.
    !            once initialized, wvhaes can be used repeatedly by vhaes
    !            as long as nlon and nlat remain unchanged.  wvhaes must
    !            not be altered between calls of vhaes.
    !
    !     lvhaes the dimension of the array wvhaes as it appears in the
    !            program that calls vhaes. 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 lvhaes must be at least
    !
    !            l1*l2(2*nlat-l1+1)+nlon+15
    !
    !     output parameters
    !
    !     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
    !            in the discription of subroutine vhses. br(mp1, np1),
    !            bi(mp1, np1), cr(mp1, np1), and ci(mp1, np1) are computed
    !            for mp1=1, ..., mmax and np1=mp1, ..., nlat except for np1=nlat
    !            and odd mp1. mmax=min(nlat, nlon/2) if nlon is even or
    !            mmax=min(nlat, (nlon + 1)/2) if nlon is odd.
    !
    !     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 lvhaes
    !
    module subroutine vhaes(nlat, nlon, ityp, nt, v, w, idvw, jdvw, br, bi, cr, ci, &
        mdab, ndab, wvhaes, 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(in)  :: v(idvw, jdvw, nt)
        real(wp),    intent(in)  :: w(idvw, jdvw, nt)
        integer(ip), intent(in)  :: idvw
        integer(ip), intent(in)  :: jdvw
        real(wp),    intent(out) :: br(mdab, ndab, nt)
        real(wp),    intent(out) :: bi(mdab, ndab, nt)
        real(wp),    intent(out) :: cr(mdab, ndab, nt)
        real(wp),    intent(out) :: ci(mdab, ndab, nt)
        integer(ip), intent(in)  :: mdab
        integer(ip), intent(in)  :: ndab
        real(wp),    intent(in)  :: wvhaes(:)
        integer(ip), intent(out) :: ierror

        ! Local variables
        integer(ip) :: idv, imid, idz, ist, lnl, lzimn, mmax
        integer(ip) :: required_wavetable_size
        integer(ip) :: lwork, workspace_indices(6)
        type(VectorAnalysisUtility) :: util

        imid = (nlat + 1)/2
        mmax = min(nlat, (nlon + 1)/2)
        idz = (mmax*(2*nlat-mmax+1))/2
        lzimn = idz*imid
        required_wavetable_size = (2 * lzimn) + nlon + 15

        ! Check calling arguments
        call util%check_vector_analysis_inputs(nlat, nlon, ityp, idvw, jdvw, &
            mdab, ndab, nt, required_wavetable_size, wvhaes, ierror)

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

        select case (ityp)
            case(0:2)
                idv = nlat
                ist = imid
            case default
                idv = imid
                ist = 0
        end select

        lnl = nt*idv*nlon

        ! Set required workspace size
        lwork = 2*lnl+idv*nlon

        block
            real(wp) :: work(lwork)

            !  Set workspace indices
            workspace_indices = get_vhaes_workspace_indices(ist, lnl, lzimn)

            associate (&
                iw1 => workspace_indices(1), &
                iw2 => workspace_indices(2), &
                iw3 => workspace_indices(3), &
                iw4 => workspace_indices(4), &
                jw1 => workspace_indices(5), &
                jw2 => workspace_indices(6) &
                )
                call vhaes_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:), idz, wvhaes, wvhaes(jw1:), wvhaes(jw2:))
            end associate
        end block

    end subroutine vhaes

    ! Purpose:
    !
    !     subroutine vhaesi(nlat, nlon, wvhaes, ierror)
    !
    !     subroutine vhaesi initializes the array wvhaes which can then be
    !     used repeatedly by subroutine vhaes until nlat or nlon is changed.
    !
    !     input parameters
    !
    !     nlat   the number of colatitudes on the full sphere including the
    !            poles. for example, nlat = 37 for a five degree grid.
    !            nlat determines the grid increment in colatitude as
    !            pi/(nlat-1).  if nlat is odd the equator is located at
    !            grid point i=(nlat + 1)/2. if nlat is even the equator is
    !            located half way between points i=nlat/2 and i=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.
    !
    !     lvhaes the dimension of the array wvhaes as it appears in the
    !            program that calls vhaes. 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 lvhaes must be at least
    !
    !               l1*l2*(2*nlat-l1+1)+nlon+15
    !
    !     output parameters
    !
    !     wvhaes an array which is initialized for use by subroutine vhaes.
    !            once initialized, wvhaes can be used repeatedly by vhaes
    !            as long as nlat or nlon remain unchanged.  wvhaes must not
    !            be altered between calls of vhaes.
    !
    !
    !     ierror = 0  no errors
    !            = 1  error in the specification of nlat
    !            = 2  error in the specification of nlon
    !            = 3  error in the specification of lvhaes
    !
    module subroutine vhaesi(nlat, nlon, wvhaes, ierror)

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

        ! Local variables
        integer(ip) :: imid, labc, lzimn, mmax
        integer(ip) :: required_wavetable_size
        integer(ip) :: workspace_indices(4)
        integer(ip) :: lwork, ldwork
        type(SpherepackUtility) :: util

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

        ! Check calling arguments
        if (nlat < 3) then
            ierror = 1
        else if (nlon < 1) then
            ierror = 2
        else if (size(wvhaes) < required_wavetable_size) then
            ierror = 3
        else
            ierror = 0
        end if

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

        !  Set workspace indices
        workspace_indices = get_vhaesi_workspace_indices(lzimn, nlat, imid, mmax)

        ! Set required workspace sizes
        lwork = 5*nlat*imid+labc
        ldwork = 2 * (nlat + 1)

        block
            real(wp) :: work(lwork), dwork(ldwork)

            associate (&
                jw1 => workspace_indices(1), &
                jw2 => workspace_indices(2), &
                iw1 => workspace_indices(3), &
                idz => workspace_indices(4) &
                )
                call vhaesi_lower_utility_routine( &
                    nlat, nlon, imid, wvhaes, wvhaes(jw1:), idz, work, work(iw1:), dwork)

                call util%hfft%initialize(nlon, wvhaes(jw2:))
            end associate
        end block

    end subroutine vhaesi

    pure function get_vhaes_workspace_indices(ist, lnl, lzimn) &
        result (return_value)

        ! Dummy arguments
        integer(ip), intent(in)  :: ist
        integer(ip), intent(in)  :: lnl
        integer(ip), intent(in)  :: lzimn
        integer(ip)              :: return_value(6)

        associate (i => return_value)
            i(1) = ist+1
            i(2) = lnl+1
            i(3) = i(2)+ist
            i(4) = i(2)+lnl
            i(5) = lzimn+1
            i(6) = i(5)+lzimn
        end associate

    end function get_vhaes_workspace_indices

    subroutine vhaes_lower_utility_routine(nlat, nlon, ityp, nt, imid, idvw, jdvw, v, w, mdab, &
        ndab, br, bi, cr, ci, idv, ve, vo, we, wo, work, idz, zv, zw, wrfft)

        ! Dummy arguments
        integer(ip), intent(in)  :: nlat
        integer(ip), intent(in)  :: nlon
        integer(ip), intent(in)  :: ityp
        integer(ip), intent(in)  :: nt
        integer(ip), intent(in)  :: imid
        integer(ip), intent(in)  :: idvw
        integer(ip), intent(in)  :: jdvw
        real(wp),    intent(in)  :: v(idvw, jdvw, nt)
        real(wp),    intent(in)  :: w(idvw, jdvw, nt)
        integer(ip), intent(in)  :: mdab
        integer(ip), intent(in)  :: ndab
        real(wp),    intent(out) :: br(mdab, ndab, nt)
        real(wp),    intent(out) :: bi(mdab, ndab, nt)
        real(wp),    intent(out) :: cr(mdab, ndab, nt)
        real(wp),    intent(out) :: ci(mdab, ndab, nt)
        integer(ip), intent(in)  :: idv
        real(wp),    intent(out) :: ve(idv, nlon, nt)
        real(wp),    intent(out) :: vo(idv, nlon, nt)
        real(wp),    intent(out) :: we(idv, nlon, nt)
        real(wp),    intent(out) :: wo(idv, nlon, nt)
        real(wp),    intent(out) :: work(*)
        integer(ip), intent(in)  :: idz
        real(wp),    intent(in)  :: zv(idz, *)
        real(wp),    intent(in)  :: zw(idz, *)
        real(wp),    intent(in)  :: wrfft(:)

        ! Local variables
        integer(ip) :: i, imm1, k, m, mb
        integer(ip) :: mmax, mp1, mp2, odd_stride, even_stride, np1
        type(VectorAnalysisUtility) :: util

        call util%analysis_setup(idvw, jdvw, mdab, ndab, &
            bi, br, ci, cr, even_stride, idv, imid, imm1, ityp, &
            mmax, nlat, nlon, nt, odd_stride,  v, ve, vo, w, we, wo, wrfft)

        !  Compute coefficients br, bi, cr, ci
        select case (ityp)
            case (0)

                !  case ityp=0,  no symmetries

                ! case m=0
                do k=1, nt
                    do i=1, imid
                        do np1=2, even_stride, 2
                            br(1, np1, k) = br(1, np1, k)+zv(np1, i)*ve(i, 1, k)
                            cr(1, np1, k) = cr(1, np1, k)-zv(np1, i)*we(i, 1, k)
                        end do
                    end do
                end do

                do k=1, nt
                    do i=1, imm1
                        do np1=3, odd_stride, 2
                            br(1, np1, k) = br(1, np1, k)+zv(np1, i)*vo(i, 1, k)
                            cr(1, np1, k) = cr(1, np1, k)-zv(np1, i)*wo(i, 1, k)
                        end do
                    end do
                end do

                ! case m = 1 through nlat-1
                if (mmax < 2) return

                main_loop_case_0: do mp1=2, mmax
                    m = mp1-1
                    mb = m*(nlat-1)-(m*(m-1))/2
                    mp2 = mp1+1

                    if (mp1 <= odd_stride) then
                        do k=1, nt
                            do i=1, imm1
                                do np1=mp1, odd_stride, 2
                                    br(mp1, np1, k) = br(mp1, np1, k)+zv(np1+mb, i)*vo(i, 2*mp1-2, k) &
                                        +zw(np1+mb, i)*we(i, 2*mp1-1, k)
                                    bi(mp1, np1, k) = bi(mp1, np1, k)+zv(np1+mb, i)*vo(i, 2*mp1-1, k) &
                                        -zw(np1+mb, i)*we(i, 2*mp1-2, k)
                                    cr(mp1, np1, k) = cr(mp1, np1, k)-zv(np1+mb, i)*wo(i, 2*mp1-2, k) &
                                        +zw(np1+mb, i)*ve(i, 2*mp1-1, k)
                                    ci(mp1, np1, k) = ci(mp1, np1, k)-zv(np1+mb, i)*wo(i, 2*mp1-1, k) &
                                        -zw(np1+mb, i)*ve(i, 2*mp1-2, k)
                                end do
                            end do
                        end do
                        if (mod(nlat, 2) /= 0) then
                            do k=1, nt
                                do np1=mp1, odd_stride, 2
                                    br(mp1, np1, k) = br(mp1, np1, k)+zw(np1+mb, imid)*we(imid, 2*mp1-1, k)
                                    bi(mp1, np1, k) = bi(mp1, np1, k)-zw(np1+mb, imid)*we(imid, 2*mp1-2, k)
                                    cr(mp1, np1, k) = cr(mp1, np1, k)+zw(np1+mb, imid)*ve(imid, 2*mp1-1, k)
                                    ci(mp1, np1, k) = ci(mp1, np1, k)-zw(np1+mb, imid)*ve(imid, 2*mp1-2, k)
                                end do
                            end do
                        end if
                    end if

                    if (mp2 > even_stride) cycle main_loop_case_0

                    do k=1, nt
                        do i=1, imm1
                            do np1=mp2, even_stride, 2
                                br(mp1, np1, k) = br(mp1, np1, k)+zv(np1+mb, i)*ve(i, 2*mp1-2, k) &
                                    +zw(np1+mb, i)*wo(i, 2*mp1-1, k)
                                bi(mp1, np1, k) = bi(mp1, np1, k)+zv(np1+mb, i)*ve(i, 2*mp1-1, k) &
                                    -zw(np1+mb, i)*wo(i, 2*mp1-2, k)
                                cr(mp1, np1, k) = cr(mp1, np1, k)-zv(np1+mb, i)*we(i, 2*mp1-2, k) &
                                    +zw(np1+mb, i)*vo(i, 2*mp1-1, k)
                                ci(mp1, np1, k) = ci(mp1, np1, k)-zv(np1+mb, i)*we(i, 2*mp1-1, k) &
                                    -zw(np1+mb, i)*vo(i, 2*mp1-2, k)
                            end do
                        end do
                    end do

                    if (mod(nlat, 2) == 0) cycle main_loop_case_0

                    do k=1, nt
                        do np1=mp2, even_stride, 2
                            br(mp1, np1, k) = br(mp1, np1, k)+zv(np1+mb, imid)*ve(imid, 2*mp1-2, k)
                            bi(mp1, np1, k) = bi(mp1, np1, k)+zv(np1+mb, imid)*ve(imid, 2*mp1-1, k)
                            cr(mp1, np1, k) = cr(mp1, np1, k)-zv(np1+mb, imid)*we(imid, 2*mp1-2, k)
                            ci(mp1, np1, k) = ci(mp1, np1, k)-zv(np1+mb, imid)*we(imid, 2*mp1-1, k)
                        end do
                    end do
                end do main_loop_case_0
            case (1)
                !
                !  case ityp=1 ,  no symmetries but cr and ci equal zero
                !
                !    case m=0
                !
                do k=1, nt
                    do i=1, imid
                        do np1=2, even_stride, 2
                            br(1, np1, k) = br(1, np1, k)+zv(np1, i)*ve(i, 1, k)
                        end do
                    end do
                end do

                do k=1, nt
                    do i=1, imm1
                        do np1=3, odd_stride, 2
                            br(1, np1, k) = br(1, np1, k)+zv(np1, i)*vo(i, 1, k)
                        end do
                    end do
                end do
                !
                !  case m = 1 through nlat-1
                !
                if (mmax < 2) return

                do mp1=2, mmax
                    m = mp1-1
                    mb = m*(nlat-1)-(m*(m-1))/2
                    mp2 = mp1+1

                    if (mp1 <= odd_stride) then
                        do k=1, nt
                            do i=1, imm1
                                do np1=mp1, odd_stride, 2
                                    br(mp1, np1, k) = br(mp1, np1, k)+zv(np1+mb, i)*vo(i, 2*mp1-2, k) &
                                        +zw(np1+mb, i)*we(i, 2*mp1-1, k)
                                    bi(mp1, np1, k) = bi(mp1, np1, k)+zv(np1+mb, i)*vo(i, 2*mp1-1, k) &
                                        -zw(np1+mb, i)*we(i, 2*mp1-2, k)
                                end do
                            end do
                        end do

                        if (mod(nlat, 2) /= 0) then
                            do k=1, nt
                                do np1=mp1, odd_stride, 2
                                    br(mp1, np1, k) = br(mp1, np1, k)+zw(np1+mb, imid)*we(imid, 2*mp1-1, k)
                                    bi(mp1, np1, k) = bi(mp1, np1, k)-zw(np1+mb, imid)*we(imid, 2*mp1-2, k)
                                end do
                            end do
                        end if
                    end if

                    if (mp2 > even_stride) return

                    do k=1, nt
                        do i=1, imm1
                            do np1=mp2, even_stride, 2
                                br(mp1, np1, k) = br(mp1, np1, k)+zv(np1+mb, i)*ve(i, 2*mp1-2, k) &
                                    +zw(np1+mb, i)*wo(i, 2*mp1-1, k)
                                bi(mp1, np1, k) = bi(mp1, np1, k)+zv(np1+mb, i)*ve(i, 2*mp1-1, k) &
                                    -zw(np1+mb, i)*wo(i, 2*mp1-2, k)
                            end do
                        end do
                    end do

                    if (mod(nlat, 2) == 0) return

                    do k=1, nt
                        do np1=mp2, even_stride, 2
                            br(mp1, np1, k) = br(mp1, np1, k)+zv(np1+mb, imid)*ve(imid, 2*mp1-2, k)
                            bi(mp1, np1, k) = bi(mp1, np1, k)+zv(np1+mb, imid)*ve(imid, 2*mp1-1, k)
                        end do
                    end do
                end do
            case (2)
                !
                !  case ityp=2 ,  no symmetries but br and bi equal zero
                !
                !    case m=0
                !
                do k=1, nt
                    do i=1, imid
                        do np1=2, even_stride, 2
                            cr(1, np1, k) = cr(1, np1, k)-zv(np1, i)*we(i, 1, k)
                        end do
                    end do
                end do

                do k=1, nt
                    do i=1, imm1
                        do np1=3, odd_stride, 2
                            cr(1, np1, k) = cr(1, np1, k)-zv(np1, i)*wo(i, 1, k)
                        end do
                    end do
                end do
                !
                !  case m = 1 through nlat-1
                !
                if (mmax < 2) return

                do mp1=2, mmax
                    m = mp1-1
                    mb = m*(nlat-1)-(m*(m-1))/2
                    mp2 = mp1+1

                    if (mp1 <= odd_stride) then
                        do k=1, nt
                            do i=1, imm1
                                do np1=mp1, odd_stride, 2
                                    cr(mp1, np1, k) = cr(mp1, np1, k)-zv(np1+mb, i)*wo(i, 2*mp1-2, k) &
                                        +zw(np1+mb, i)*ve(i, 2*mp1-1, k)
                                    ci(mp1, np1, k) = ci(mp1, np1, k)-zv(np1+mb, i)*wo(i, 2*mp1-1, k) &
                                        -zw(np1+mb, i)*ve(i, 2*mp1-2, k)
                                end do
                            end do
                        end do

                        if (mod(nlat, 2) /= 0) then
                            do k=1, nt
                                do np1=mp1, odd_stride, 2
                                    cr(mp1, np1, k) = cr(mp1, np1, k)+zw(np1+mb, imid)*ve(imid, 2*mp1-1, k)
                                    ci(mp1, np1, k) = ci(mp1, np1, k)-zw(np1+mb, imid)*ve(imid, 2*mp1-2, k)
                                end do
                            end do
                        end if
                    end if

                    if (mp2 > even_stride) return

                    do k=1, nt
                        do i=1, imm1
                            do np1=mp2, even_stride, 2
                                cr(mp1, np1, k) = cr(mp1, np1, k)-zv(np1+mb, i)*we(i, 2*mp1-2, k) &
                                    +zw(np1+mb, i)*vo(i, 2*mp1-1, k)
                                ci(mp1, np1, k) = ci(mp1, np1, k)-zv(np1+mb, i)*we(i, 2*mp1-1, k) &
                                    -zw(np1+mb, i)*vo(i, 2*mp1-2, k)
                            end do
                        end do
                    end do

                    if (mod(nlat, 2) == 0) return

                    do k=1, nt
                        do np1=mp2, even_stride, 2
                            cr(mp1, np1, k) = cr(mp1, np1, k)-zv(np1+mb, imid)*we(imid, 2*mp1-2, k)
                            ci(mp1, np1, k) = ci(mp1, np1, k)-zv(np1+mb, imid)*we(imid, 2*mp1-1, k)
                        end do
                    end do
                end do
            case (3)
                !
                !  case ityp=3 ,  v even , w odd
                !
                !    case m=0
                !
                do k=1, nt
                    do i=1, imid
                        do np1=2, even_stride, 2
                            br(1, np1, k) = br(1, np1, k)+zv(np1, i)*ve(i, 1, k)
                        end do
                    end do
                end do

                do k=1, nt
                    do i=1, imm1
                        do np1=3, odd_stride, 2
                            cr(1, np1, k) = cr(1, np1, k)-zv(np1, i)*wo(i, 1, k)
                        end do
                    end do
                end do
                !
                !  case m = 1 through nlat-1
                !
                if (mmax < 2) return

                do mp1=2, mmax
                    m = mp1-1
                    mb = m*(nlat-1)-(m*(m-1))/2
                    mp2 = mp1+1

                    if (mp1 <= odd_stride) then
                        do k=1, nt
                            do i=1, imm1
                                do np1=mp1, odd_stride, 2
                                    cr(mp1, np1, k) = cr(mp1, np1, k)-zv(np1+mb, i)*wo(i, 2*mp1-2, k) &
                                        +zw(np1+mb, i)*ve(i, 2*mp1-1, k)
                                    ci(mp1, np1, k) = ci(mp1, np1, k)-zv(np1+mb, i)*wo(i, 2*mp1-1, k) &
                                        -zw(np1+mb, i)*ve(i, 2*mp1-2, k)
                                end do
                            end do
                        end do

                        if (mod(nlat, 2) /= 0) then
                            do k=1, nt
                                do np1=mp1, odd_stride, 2
                                    cr(mp1, np1, k) = cr(mp1, np1, k)+zw(np1+mb, imid)*ve(imid, 2*mp1-1, k)
                                    ci(mp1, np1, k) = ci(mp1, np1, k)-zw(np1+mb, imid)*ve(imid, 2*mp1-2, k)
                                end do
                            end do
                        end if
                    end if

                    if (mp2 > even_stride) return

                    do k=1, nt
                        do i=1, imm1
                            do np1=mp2, even_stride, 2
                                br(mp1, np1, k) = br(mp1, np1, k)+zv(np1+mb, i)*ve(i, 2*mp1-2, k) &
                                    +zw(np1+mb, i)*wo(i, 2*mp1-1, k)
                                bi(mp1, np1, k) = bi(mp1, np1, k)+zv(np1+mb, i)*ve(i, 2*mp1-1, k) &
                                    -zw(np1+mb, i)*wo(i, 2*mp1-2, k)
                            end do
                        end do
                    end do

                    if (mod(nlat, 2) == 0) return

                    do  k=1, nt
                        do  np1=mp2, even_stride, 2
                            br(mp1, np1, k) = br(mp1, np1, k)+zv(np1+mb, imid)*ve(imid, 2*mp1-2, k)
                            bi(mp1, np1, k) = bi(mp1, np1, k)+zv(np1+mb, imid)*ve(imid, 2*mp1-1, k)
                        end do
                    end do
                end do
            case (4)
                !
                !  case ityp=4 ,  v even, w odd, and cr and ci equal 0.
                !
                !    case m=0
                !
                do k=1, nt
                    do i=1, imid
                        do np1=2, even_stride, 2
                            br(1, np1, k) = br(1, np1, k)+zv(np1, i)*ve(i, 1, k)
                        end do
                    end do
                end do
                !
                !  case m = 1 through nlat-1
                !
                if (mmax < 2) return

                do mp1=2, mmax
                    m = mp1-1
                    mb = m*(nlat-1)-(m*(m-1))/2
                    mp2 = mp1+1

                    if (mp2 > even_stride) exit

                    do k=1, nt
                        do i=1, imm1
                            do np1=mp2, even_stride, 2
                                br(mp1, np1, k) = br(mp1, np1, k)+zv(np1+mb, i)*ve(i, 2*mp1-2, k) &
                                    +zw(np1+mb, i)*wo(i, 2*mp1-1, k)
                                bi(mp1, np1, k) = bi(mp1, np1, k)+zv(np1+mb, i)*ve(i, 2*mp1-1, k) &
                                    -zw(np1+mb, i)*wo(i, 2*mp1-2, k)
                            end do
                        end do
                    end do

                    if (mod(nlat, 2) == 0) return

                    do k=1, nt
                        do np1=mp2, even_stride, 2
                            br(mp1, np1, k) = br(mp1, np1, k)+zv(np1+mb, imid)*ve(imid, 2*mp1-2, k)
                            bi(mp1, np1, k) = bi(mp1, np1, k)+zv(np1+mb, imid)*ve(imid, 2*mp1-1, k)
                        end do
                    end do
                end do
            case (5)
                !
                !  case ityp=5   v even, w odd, and br and bi equal zero
                !
                !    case m=0
                !
                do k=1, nt
                    do i=1, imm1
                        do np1=3, odd_stride, 2
                            cr(1, np1, k) = cr(1, np1, k)-zv(np1, i)*wo(i, 1, k)
                        end do
                    end do
                end do
                !
                !  case m = 1 through nlat-1
                !
                if (mmax < 2) return

                do mp1=2, mmax
                    m = mp1-1
                    mb = m*(nlat-1)-(m*(m-1))/2
                    mp2 = mp1+1

                    if (mp1 > odd_stride) exit

                    do k=1, nt
                        do i=1, imm1
                            do np1=mp1, odd_stride, 2
                                cr(mp1, np1, k) = cr(mp1, np1, k)-zv(np1+mb, i)*wo(i, 2*mp1-2, k) &
                                    +zw(np1+mb, i)*ve(i, 2*mp1-1, k)
                                ci(mp1, np1, k) = ci(mp1, np1, k)-zv(np1+mb, i)*wo(i, 2*mp1-1, k) &
                                    -zw(np1+mb, i)*ve(i, 2*mp1-2, k)
                            end do
                        end do
                    end do

                    if (mod(nlat, 2) == 0) exit

                    do k=1, nt
                        do np1=mp1, odd_stride, 2
                            cr(mp1, np1, k) = cr(mp1, np1, k)+zw(np1+mb, imid)*ve(imid, 2*mp1-1, k)
                            ci(mp1, np1, k) = ci(mp1, np1, k)-zw(np1+mb, imid)*ve(imid, 2*mp1-2, k)
                        end do
                    end do
                end do
            case (6)
                !
                !  case ityp=6 ,  v odd , w even
                !
                !    case m=0
                !
                do k=1, nt
                    do i=1, imid
                        do np1=2, even_stride, 2
                            cr(1, np1, k) = cr(1, np1, k)-zv(np1, i)*we(i, 1, k)
                        end do
                    end do
                end do

                do k=1, nt
                    do i=1, imm1
                        do np1=3, odd_stride, 2
                            br(1, np1, k) = br(1, np1, k)+zv(np1, i)*vo(i, 1, k)
                        end do
                    end do
                end do
                !
                !  case m = 1 through nlat-1
                !
                if (mmax < 2) return

                do mp1=2, mmax
                    m = mp1-1
                    mb = m*(nlat-1)-(m*(m-1))/2
                    mp2 = mp1+1

                    if (mp1 <= odd_stride) then

                        do  k=1, nt
                            do  i=1, imm1
                                do  np1=mp1, odd_stride, 2
                                    br(mp1, np1, k) = br(mp1, np1, k)+zv(np1+mb, i)*vo(i, 2*mp1-2, k) &
                                        +zw(np1+mb, i)*we(i, 2*mp1-1, k)
                                    bi(mp1, np1, k) = bi(mp1, np1, k)+zv(np1+mb, i)*vo(i, 2*mp1-1, k) &
                                        -zw(np1+mb, i)*we(i, 2*mp1-2, k)
                                end do
                            end do
                        end do

                        if (mod(nlat, 2) /= 0) then
                            do  k=1, nt
                                do  np1=mp1, odd_stride, 2
                                    br(mp1, np1, k) = br(mp1, np1, k)+zw(np1+mb, imid)*we(imid, 2*mp1-1, k)
                                    bi(mp1, np1, k) = bi(mp1, np1, k)-zw(np1+mb, imid)*we(imid, 2*mp1-2, k)
                                end do
                            end do
                        end if
                    end if

                    if (mp2 > even_stride) return

                    do k=1, nt
                        do  i=1, imm1
                            do  np1=mp2, even_stride, 2
                                cr(mp1, np1, k) = cr(mp1, np1, k)-zv(np1+mb, i)*we(i, 2*mp1-2, k) &
                                    +zw(np1+mb, i)*vo(i, 2*mp1-1, k)
                                ci(mp1, np1, k) = ci(mp1, np1, k)-zv(np1+mb, i)*we(i, 2*mp1-1, k) &
                                    -zw(np1+mb, i)*vo(i, 2*mp1-2, k)
                            end do
                        end do
                    end do

                    if (mod(nlat, 2) == 0) return

                    do k=1, nt
                        do np1=mp2, even_stride, 2
                            cr(mp1, np1, k) = cr(mp1, np1, k)-zv(np1+mb, imid)*we(imid, 2*mp1-2, k)
                            ci(mp1, np1, k) = ci(mp1, np1, k)-zv(np1+mb, imid)*we(imid, 2*mp1-1, k)
                        end do
                    end do
                end do
            case (7)
                !
                !  case ityp=7   v odd, w even, and cr and ci equal zero
                !
                !    case m=0
                !
                do k=1, nt
                    do i=1, imm1
                        do np1=3, odd_stride, 2
                            br(1, np1, k) = br(1, np1, k)+zv(np1, i)*vo(i, 1, k)
                        end do
                    end do
                end do
                !
                !  case m = 1 through nlat-1
                !
                if (mmax < 2) return

                do mp1=2, mmax
                    m = mp1-1
                    mb = m*(nlat-1)-(m*(m-1))/2
                    mp2 = mp1+1

                    if (mp1 > odd_stride) return

                    do  k=1, nt
                        do  i=1, imm1
                            do  np1=mp1, odd_stride, 2
                                br(mp1, np1, k) = br(mp1, np1, k)+zv(np1+mb, i)*vo(i, 2*mp1-2, k) &
                                    +zw(np1+mb, i)*we(i, 2*mp1-1, k)
                                bi(mp1, np1, k) = bi(mp1, np1, k)+zv(np1+mb, i)*vo(i, 2*mp1-1, k) &
                                    -zw(np1+mb, i)*we(i, 2*mp1-2, k)
                            end do
                        end do
                    end do

                    if (mod(nlat, 2) == 0) return

                    do k=1, nt
                        do np1=mp1, odd_stride, 2
                            br(mp1, np1, k) = br(mp1, np1, k)+zw(np1+mb, imid)*we(imid, 2*mp1-1, k)
                            bi(mp1, np1, k) = bi(mp1, np1, k)-zw(np1+mb, imid)*we(imid, 2*mp1-2, k)
                        end do
                    end do
                end do
            case (8)
                !
                !  case ityp=8   v odd, w even, and both br and bi equal zero
                !
                !    case m=0
                !
                do k=1, nt
                    do i=1, imid
                        do np1=2, even_stride, 2
                            cr(1, np1, k) = cr(1, np1, k)-zv(np1, i)*we(i, 1, k)
                        end do
                    end do
                end do
                !
                ! case m = 1 through nlat-1
                !
                if (mmax < 2) return

                do mp1=2, mmax
                    m = mp1-1
                    mb = m*(nlat-1)-(m*(m-1))/2
                    mp2 = mp1+1

                    if (mp2 > even_stride) return

                    do k=1, nt
                        do i=1, imm1
                            do np1=mp2, even_stride, 2
                                cr(mp1, np1, k) = cr(mp1, np1, k)-zv(np1+mb, i)*we(i, 2*mp1-2, k) &
                                    +zw(np1+mb, i)*vo(i, 2*mp1-1, k)
                                ci(mp1, np1, k) = ci(mp1, np1, k)-zv(np1+mb, i)*we(i, 2*mp1-1, k) &
                                    -zw(np1+mb, i)*vo(i, 2*mp1-2, k)
                            end do
                        end do
                    end do

                    if (mod(nlat, 2) == 0) return

                    do k=1, nt
                        do np1=mp2, even_stride, 2
                            cr(mp1, np1, k) = cr(mp1, np1, k)-zv(np1+mb, imid)*we(imid, 2*mp1-2, k)
                            ci(mp1, np1, k) = ci(mp1, np1, k)-zv(np1+mb, imid)*we(imid, 2*mp1-1, k)
                        end do
                    end do
                end do
        end select

    end subroutine vhaes_lower_utility_routine

    pure function get_vhaesi_workspace_indices(lzimn, nlat, imid, mmax) &
        result (return_value)

        ! Dummy arguments
        integer(ip), intent(in)  :: lzimn
        integer(ip), intent(in)  :: nlat
        integer(ip), intent(in)  :: imid
        integer(ip), intent(in)  :: mmax
        integer(ip)              :: return_value(4)

        associate (i => return_value)
            i(1) = lzimn+1
            i(2) = 2*lzimn+1
            i(3) = 3*nlat*imid+1
            i(4) = (mmax*(2*nlat-mmax+1))/2
        end associate

    end function get_vhaesi_workspace_indices

    subroutine vhaesi_lower_utility_routine(nlat, nlon, imid, zv, zw, idz, zin, wzvin, dwork)

        ! Dummy arguments
        integer(ip), intent(in)  :: nlat
        integer(ip), intent(in)  :: nlon
        integer(ip), intent(in)  :: imid
        real(wp),    intent(out) :: zv(idz, *)
        real(wp),    intent(out) :: zw(idz, *)
        integer(ip), intent(in)  :: idz
        real(wp),    intent(out) :: zin(imid, nlat, 3)
        real(wp),    intent(out) :: wzvin(*)
        real(wp),    intent(out) :: dwork(*)

        ! Local variables
        integer(ip)         :: i3, m, mn, mp1, np1, mmax
        type(SpherepackUtility) :: util

        mmax = min(nlat, (nlon + 1)/2)

        call util%zvinit(nlat, nlon, wzvin, dwork)

        do mp1=1, mmax
            m = mp1-1
            call util%zvin(0, nlat, nlon, m, zin, i3, wzvin)
            do np1=mp1, nlat
                mn = m*(nlat-1)-(m*(m-1))/2+np1
                zv(mn, 1:imid) = zin(1:imid, np1, i3)
            end do
        end do

        call util%zwinit(nlat, nlon, wzvin, dwork)

        do mp1=1, mmax
            m = mp1-1
            call util%zwin(0, nlat, nlon, m, zin, i3, wzvin)
            do np1=mp1, nlat
                mn = m*(nlat-1)-(m*(m-1))/2+np1
                zw(mn, 1:imid) = zin(1:imid, np1, i3)
            end do
        end do

    end subroutine vhaesi_lower_utility_routine

end submodule vector_analysis_regular_grid_saved


text:
      MODULE inid_I
      INTERFACE
!...Generated by Pacific-Sierra Research 77to90  4.4G  12:41:19  03/10/06
      SUBROUTINE inid
      END SUBROUTINE
      END INTERFACE
      END MODULE

max_stars_repo_name: openmopac/MOPAC-archive

max_issues_repo_name: openmopac/MOPAC-archive

max_forks_repo_name: openmopac/MOPAC-archive


text:
The Electric Trike Guy has been seen in all parts of Davis, but is more often spotted in Mace Ranch.  He seems to wear the same clothes and cap regardless of the weather.

20090915 11:42:39 nbsp Electric trike guy lives near Chestnut park, and can sometimes be rolling slow when his cart is dying, or pushing it when it finally has. Users/StevenDaubert
    Kind of disturbed why this guy is being spotlightedhe may simply have a physical disability that requires him to use this cart. Users/PeterBoulay
  not really, keep your political correctness on the East coast where it belongs. Dude is seen throughout town constantly and Im sure people wondered. I just shared my dealings with him via this tool of collaborative authorship Daubert
   Really Daubert, you have dealt with him? Could you tell us his name? Thanks.  Users/NickSchmalenberger

20111209 09:43:18 nbsp Didnt his trike used to have an awning? Striped red and white? It was quite delightful! Users/BonnieStewart
 Especially when the weather outside was frightful! Users/TomGarberson
I have a couple photos of him in traffic, I have been meaning to interview him Daubert

20111209 12:04 The comments herein have made my morning Brenton

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      DOUBLE PRECISION   FUNCTION PDLAMCH10( ICTXT, CMACH )
      include "mpif.h"
*
*  -- ScaLAPACK auxiliary routine (version 1.0) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     February 28, 1995
*
*     The name has been changed in order to avoid symbol collision with
*     the newer version of PDLAMCH in Scalapack which can not be used
*     with MPI context.
*
*     .. Scalar Arguments ..
      CHARACTER          CMACH
      INTEGER            ICTXT
*     ..
*
*  Purpose
*  =======
*
*  PDLAMCH determines double precision machine parameters.
*
*  Arguments
*  =========
*
*  ICTXT   (global input) INTEGER
*          The MPI context handle in which the computation takes
*          place.
*
*  CMACH   (global input) CHARACTER*1
*          Specifies the value to be returned by PDLAMCH:
*          = 'E' or 'e',   PDLAMCH := eps
*          = 'S' or 's ,   PDLAMCH := sfmin
*          = 'B' or 'b',   PDLAMCH := base
*          = 'P' or 'p',   PDLAMCH := eps*base
*          = 'N' or 'n',   PDLAMCH := t
*          = 'R' or 'r',   PDLAMCH := rnd
*          = 'M' or 'm',   PDLAMCH := emin
*          = 'U' or 'u',   PDLAMCH := rmin
*          = 'L' or 'l',   PDLAMCH := emax
*          = 'O' or 'o',   PDLAMCH := rmax
*
*          where
*
*          eps   = relative machine precision
*          sfmin = safe minimum, such that 1/sfmin does not overflow
*          base  = base of the machine
*          prec  = eps*base
*          t     = number of (base) digits in the mantissa
*          rnd   = 1.0 when rounding occurs in addition, 0.0 otherwise
*          emin  = minimum exponent before (gradual) underflow
*          rmin  = underflow threshold - base**(emin-1)
*          emax  = largest exponent before overflow
*          rmax  = overflow threshold  - (base**emax)*(1-eps)
*
*  =====================================================================
*
*     .. Local Scalars ..
      INTEGER            IDUMM
      DOUBLE PRECISION   TEMP, TEMP1
*     ..
*     .. External Subroutines ..
*      EXTERNAL           DGAMN2D, DGAMX2D
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      DOUBLE PRECISION   DLAMCH
      EXTERNAL           DLAMCH, LSAME
*     ..
*     .. Executable Statements ..
*
      TEMP1 = DLAMCH( CMACH )
*
      IF( LSAME( CMACH, 'E' ).OR.LSAME( CMACH, 'S' ).OR.
     $    LSAME( CMACH, 'M' ).OR.LSAME( CMACH, 'U' ) ) THEN
          CALL MPI_ALLREDUCE( TEMP1, TEMP, 1, MPI_DOUBLE_PRECISION,
     $                        MPI_MAX, ICTXT, IDUMM )
*         CALL DGAMX2D( ICTXT, 'All', ' ', 1, 1, TEMP, 1, IDUMM,
*     $                 IDUMM, 1, -1, IDUMM )
      ELSE IF( LSAME( CMACH, 'L' ).OR.LSAME( CMACH, 'O' ) ) THEN
          CALL MPI_ALLREDUCE( TEMP1, TEMP, 1, MPI_DOUBLE_PRECISION,
     $                        MPI_MIN, ICTXT, IDUMM )
*         CALL DGAMN2D( ICTXT, 'All', ' ', 1, 1, TEMP, 1, IDUMM,
*     $                 IDUMM, 1, -1, IDUMM )
      ELSE
          TEMP = TEMP1
      END IF
*
      PDLAMCH10 = TEMP
*
*     End of PDLAMCH10
*
      END

max_stars_repo_name: XpressAI/frovedis

max_issues_repo_name: XpressAI/frovedis

max_forks_repo_name: XpressAI/frovedis


!-------------------------------------------------------------------------------
!> Module vertical metrics
!!
!! @par Description
!!          In this module, the vertical metrics is calculated for the icoshaedral model
!!
!! @author NICAM developers, Team SCALE
!<
!-------------------------------------------------------------------------------
module mod_vmtr
  !-----------------------------------------------------------------------------
  !
  !++ Used modules
  !
  use scale_precision
  use scale_stdio
  use scale_prof

  use mod_adm, only: &
     ADM_LOG_FID
  use mod_adm, only: &
     ADM_TI,      &
     ADM_TJ,      &
     ADM_AI,      &
     ADM_AIJ,     &
     ADM_AJ,      &
     ADM_lall,    &
     ADM_lall_pl, &
     ADM_gall,    &
     ADM_gall_pl, &
     ADM_kall
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: VMTR_setup

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !

  ! index for VMTR_C2Wfact, W2Cfact
  integer, public, parameter :: I_a = 1
  integer, public, parameter :: I_b = 2

  integer, public, parameter :: I_c = 1
  integer, public, parameter :: I_d = 2

  integer, public, parameter :: I_a_GZXH = 1
  integer, public, parameter :: I_b_GZXH = 2
  integer, public, parameter :: I_a_GZYH = 3
  integer, public, parameter :: I_b_GZYH = 4
  integer, public, parameter :: I_a_GZZH = 5
  integer, public, parameter :: I_b_GZZH = 6

#ifdef _FIXEDINDEX_
  real(RP), public              :: VMTR_GAM2        (ADM_gall   ,ADM_kall,ADM_lall   )
  real(RP), public              :: VMTR_GAM2_pl     (ADM_gall_pl,ADM_kall,ADM_lall_pl)
  real(RP), public              :: VMTR_GAM2H       (ADM_gall   ,ADM_kall,ADM_lall   )
  real(RP), public              :: VMTR_GAM2H_pl    (ADM_gall_pl,ADM_kall,ADM_lall_pl)
  real(RP), public              :: VMTR_GSGAM2      (ADM_gall   ,ADM_kall,ADM_lall   )
  real(RP), public              :: VMTR_GSGAM2_pl   (ADM_gall_pl,ADM_kall,ADM_lall_pl)
  real(RP), public              :: VMTR_GSGAM2H     (ADM_gall   ,ADM_kall,ADM_lall   )
  real(RP), public              :: VMTR_GSGAM2H_pl  (ADM_gall_pl,ADM_kall,ADM_lall_pl)

  real(RP), public              :: VMTR_RGSQRTH     (ADM_gall   ,ADM_kall,ADM_lall   )
  real(RP), public              :: VMTR_RGSQRTH_pl  (ADM_gall_pl,ADM_kall,ADM_lall_pl)
  real(RP), public              :: VMTR_RGAM        (ADM_gall   ,ADM_kall,ADM_lall   )
  real(RP), public              :: VMTR_RGAM_pl     (ADM_gall_pl,ADM_kall,ADM_lall_pl)
  real(RP), public              :: VMTR_RGAMH       (ADM_gall   ,ADM_kall,ADM_lall   )
  real(RP), public              :: VMTR_RGAMH_pl    (ADM_gall_pl,ADM_kall,ADM_lall_pl)
  real(RP), public              :: VMTR_RGSGAM2     (ADM_gall   ,ADM_kall,ADM_lall   )
  real(RP), public              :: VMTR_RGSGAM2_pl  (ADM_gall_pl,ADM_kall,ADM_lall_pl)
  real(RP), public              :: VMTR_RGSGAM2H    (ADM_gall   ,ADM_kall,ADM_lall   )
  real(RP), public              :: VMTR_RGSGAM2H_pl (ADM_gall_pl,ADM_kall,ADM_lall_pl)

  real(RP), public              :: VMTR_W2Cfact     (ADM_gall,   ADM_kall,2,ADM_lall   )
  real(RP), public              :: VMTR_W2Cfact_pl  (ADM_gall_pl,ADM_kall,2,ADM_lall_pl)
  real(RP), public              :: VMTR_C2Wfact     (ADM_gall,   ADM_kall,2,ADM_lall   )
  real(RP), public              :: VMTR_C2Wfact_pl  (ADM_gall_pl,ADM_kall,2,ADM_lall_pl)
  real(RP), public              :: VMTR_C2WfactGz   (ADM_gall,   ADM_kall,6,ADM_lall   )
  real(RP), public              :: VMTR_C2WfactGz_pl(ADM_gall_pl,ADM_kall,6,ADM_lall_pl)

  real(RP), public              :: VMTR_VOLUME      (ADM_gall   ,ADM_kall,ADM_lall   )
  real(RP), public              :: VMTR_VOLUME_pl   (ADM_gall_pl,ADM_kall,ADM_lall_pl)

  real(RP), public              :: VMTR_PHI         (ADM_gall   ,ADM_kall,ADM_lall   )
  real(RP), public              :: VMTR_PHI_pl      (ADM_gall_pl,ADM_kall,ADM_lall_pl)
#else
  real(RP), public, allocatable :: VMTR_GAM2        (:,:,:)   ! Gamma^2 at the full level
  real(RP), public, allocatable :: VMTR_GAM2_pl     (:,:,:)
  real(RP), public, allocatable :: VMTR_GAM2H       (:,:,:)   ! Gamma^2 at the half level
  real(RP), public, allocatable :: VMTR_GAM2H_pl    (:,:,:)
  real(RP), public, allocatable :: VMTR_GSGAM2      (:,:,:)   ! G^1/2 X Gamma^2 at the full level
  real(RP), public, allocatable :: VMTR_GSGAM2_pl   (:,:,:)
  real(RP), public, allocatable :: VMTR_GSGAM2H     (:,:,:)   ! G^1/2 X Gamma^2 at the half level
  real(RP), public, allocatable :: VMTR_GSGAM2H_pl  (:,:,:)

  real(RP), public, allocatable :: VMTR_RGSQRTH     (:,:,:)   ! 1 / G^1/2 at the half level
  real(RP), public, allocatable :: VMTR_RGSQRTH_pl  (:,:,:)
  real(RP), public, allocatable :: VMTR_RGAM        (:,:,:)   ! 1 / Gamma at the integer level
  real(RP), public, allocatable :: VMTR_RGAM_pl     (:,:,:)
  real(RP), public, allocatable :: VMTR_RGAMH       (:,:,:)   ! 1 / Gamma at the half level
  real(RP), public, allocatable :: VMTR_RGAMH_pl    (:,:,:)
  real(RP), public, allocatable :: VMTR_RGSGAM2     (:,:,:)   ! 1 / (G^1/2 X Gamma^2) at the full level
  real(RP), public, allocatable :: VMTR_RGSGAM2_pl  (:,:,:)
  real(RP), public, allocatable :: VMTR_RGSGAM2H    (:,:,:)   ! 1 / (G^1/2 X Gamma^2) at the half level
  real(RP), public, allocatable :: VMTR_RGSGAM2H_pl (:,:,:)

  real(RP), public, allocatable :: VMTR_W2Cfact     (:,:,:,:) ! factor for half to full level
  real(RP), public, allocatable :: VMTR_W2Cfact_pl  (:,:,:,:)
  real(RP), public, allocatable :: VMTR_C2Wfact     (:,:,:,:) ! factor for full to half level
  real(RP), public, allocatable :: VMTR_C2Wfact_pl  (:,:,:,:)
  real(RP), public, allocatable :: VMTR_C2WfactGz   (:,:,:,:) ! factor for full to half level with Gz
  real(RP), public, allocatable :: VMTR_C2WfactGz_pl(:,:,:,:)

  real(RP), public, allocatable :: VMTR_VOLUME      (:,:,:)   ! volume at the full level
  real(RP), public, allocatable :: VMTR_VOLUME_pl   (:,:,:)

  real(RP), public, allocatable :: VMTR_PHI         (:,:,:)   ! geopotential at the full level
  real(RP), public, allocatable :: VMTR_PHI_pl      (:,:,:)
#endif

  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  logical, private :: deep = .false.

  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  !> Setup
  subroutine VMTR_setup
    use mod_adm, only: &
       ADM_CTL_FID,   &
       ADM_proc_stop, &
       ADM_have_pl,   &
       ADM_gmin,      &
       ADM_gmax,      &
       ADM_KNONE,     &
       ADM_kmin,      &
       ADM_kmax
    use scale_const, only: &
       CONST_GRAV
    use mod_comm, only: &
       COMM_data_transfer
    use mod_grd, only: &
       GRD_Z,        &
       GRD_ZH,       &
       GRD_dgz,      &
       GRD_dgzh,     &
       GRD_vz,       &
       GRD_vz_pl,    &
       GRD_afac,     &
       GRD_bfac,     &
       GRD_cfac,     &
       GRD_dfac,     &
       GRD_rscale,   &
       GRD_grid_type
    use mod_gmtr, only: &
       GMTR_area,   &
       GMTR_area_pl
    use mod_oprt, only: &
       OPRT_gradient,         &
       OPRT_horizontalize_vec
    implicit none

    integer, parameter :: var_max = 6

    integer, parameter :: JXH     = 1
    integer, parameter :: JYH     = 2
    integer, parameter :: JZH     = 3
    integer, parameter :: JX      = 4
    integer, parameter :: JY      = 5
    integer, parameter :: JZ      = 6

    real(RP) :: var   (ADM_gall,   ADM_kall,ADM_lall,   var_max)
    real(RP) :: var_pl(ADM_gall_pl,ADM_kall,ADM_lall_pl,var_max)

    !--- G^1/2
    real(RP) :: GSQRT    (ADM_gall   ,ADM_kall,ADM_lall   )
    real(RP) :: GSQRT_pl (ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP) :: GSQRTH   (ADM_gall   ,ADM_kall,ADM_lall   )
    real(RP) :: GSQRTH_pl(ADM_gall_pl,ADM_kall,ADM_lall_pl)
    !--- Gamma factor
    real(RP) :: GAM      (ADM_gall   ,ADM_kall,ADM_lall   )
    real(RP) :: GAM_pl   (ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP) :: GAMH     (ADM_gall   ,ADM_kall,ADM_lall   )
    real(RP) :: GAMH_pl  (ADM_gall_pl,ADM_kall,ADM_lall_pl)
    !--- vector G^z at the full level
    real(RP) :: GZX      (ADM_gall   ,ADM_kall,ADM_lall   )
    real(RP) :: GZX_pl   (ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP) :: GZY      (ADM_gall   ,ADM_kall,ADM_lall   )
    real(RP) :: GZY_pl   (ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP) :: GZZ      (ADM_gall   ,ADM_kall,ADM_lall   )
    real(RP) :: GZZ_pl   (ADM_gall_pl,ADM_kall,ADM_lall_pl)
    !--- vector G^z at the half level
    real(RP) :: GZXH     (ADM_gall   ,ADM_kall,ADM_lall   )
    real(RP) :: GZXH_pl  (ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP) :: GZYH     (ADM_gall   ,ADM_kall,ADM_lall   )
    real(RP) :: GZYH_pl  (ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP) :: GZZH     (ADM_gall   ,ADM_kall,ADM_lall   )
    real(RP) :: GZZH_pl  (ADM_gall_pl,ADM_kall,ADM_lall_pl)

    namelist / VMTRPARAM / &
       deep

    integer :: ierr
    integer :: g, k, l
    !---------------------------------------------------------------------------

    !--- read parameters
    write(ADM_LOG_FID,*)
    write(ADM_LOG_FID,*) '+++ Module[vmtr]/Category[common share]'
    rewind(ADM_CTL_FID)
    read(ADM_CTL_FID,nml=VMTRPARAM,iostat=ierr)
    if ( ierr < 0 ) then
       write(ADM_LOG_FID,*) '*** VMTRPARAM is not specified. use default.'
    elseif( ierr > 0 ) then
       write(*,          *) 'xxx Not appropriate names in namelist VMTRPARAM. STOP.'
       write(ADM_LOG_FID,*) 'xxx Not appropriate names in namelist VMTRPARAM. STOP.'
       call ADM_proc_stop
    endif
    write(ADM_LOG_FID,nml=VMTRPARAM)

#ifndef _FIXEDINDEX_
    allocate( VMTR_GAM2        (ADM_gall,   ADM_kall,ADM_lall   ) )
    allocate( VMTR_GAM2_pl     (ADM_gall_pl,ADM_kall,ADM_lall_pl) )
    allocate( VMTR_GAM2H       (ADM_gall,   ADM_kall,ADM_lall   ) )
    allocate( VMTR_GAM2H_pl    (ADM_gall_pl,ADM_kall,ADM_lall_pl) )
    allocate( VMTR_GSGAM2      (ADM_gall,   ADM_kall,ADM_lall   ) )
    allocate( VMTR_GSGAM2_pl   (ADM_gall_pl,ADM_kall,ADM_lall_pl) )
    allocate( VMTR_GSGAM2H     (ADM_gall,   ADM_kall,ADM_lall   ) )
    allocate( VMTR_GSGAM2H_pl  (ADM_gall_pl,ADM_kall,ADM_lall_pl) )

    allocate( VMTR_RGSQRTH     (ADM_gall,   ADM_kall,ADM_lall   ) )
    allocate( VMTR_RGSQRTH_pl  (ADM_gall_pl,ADM_kall,ADM_lall_pl) )
    allocate( VMTR_RGAM        (ADM_gall,   ADM_kall,ADM_lall   ) )
    allocate( VMTR_RGAM_pl     (ADM_gall_pl,ADM_kall,ADM_lall_pl) )
    allocate( VMTR_RGAMH       (ADM_gall,   ADM_kall,ADM_lall   ) )
    allocate( VMTR_RGAMH_pl    (ADM_gall_pl,ADM_kall,ADM_lall_pl) )
    allocate( VMTR_RGSGAM2     (ADM_gall,   ADM_kall,ADM_lall   ) )
    allocate( VMTR_RGSGAM2_pl  (ADM_gall_pl,ADM_kall,ADM_lall_pl) )
    allocate( VMTR_RGSGAM2H    (ADM_gall,   ADM_kall,ADM_lall   ) )
    allocate( VMTR_RGSGAM2H_pl (ADM_gall_pl,ADM_kall,ADM_lall_pl) )

    allocate( VMTR_W2Cfact     (ADM_gall,   ADM_kall,2,ADM_lall   ) )
    allocate( VMTR_W2Cfact_pl  (ADM_gall_pl,ADM_kall,2,ADM_lall_pl) )
    allocate( VMTR_C2Wfact     (ADM_gall,   ADM_kall,2,ADM_lall   ) )
    allocate( VMTR_C2Wfact_pl  (ADM_gall_pl,ADM_kall,2,ADM_lall_pl) )
    allocate( VMTR_C2WfactGz   (ADM_gall,   ADM_kall,6,ADM_lall   ) )
    allocate( VMTR_C2WfactGz_pl(ADM_gall_pl,ADM_kall,6,ADM_lall_pl) )

    allocate( VMTR_VOLUME      (ADM_gall,   ADM_kall,ADM_lall   ) )
    allocate( VMTR_VOLUME_pl   (ADM_gall_pl,ADM_kall,ADM_lall_pl) )

    allocate( VMTR_PHI         (ADM_gall,   ADM_kall,ADM_lall   ) )
    allocate( VMTR_PHI_pl      (ADM_gall_pl,ADM_kall,ADM_lall_pl) )
#endif

    !--- if 1 layer model( shallow water model ),
    if ( ADM_kall == ADM_KNONE ) then

       do k = 1, ADM_kall
          VMTR_VOLUME   (:,k,:) = GMTR_area   (:,:)
          VMTR_VOLUME_pl(:,k,:) = GMTR_area_pl(:,:)
       enddo

       return
    endif

    var   (:,:,:,:) = 0.0_RP
    var_pl(:,:,:,:) = 0.0_RP

    !--- calculation of Jxh, Jyh, and Jzh
    call OPRT_gradient( GRD_vz(:,:,:,GRD_ZH),  GRD_vz_pl(:,:,:,GRD_ZH), & !--- [IN]
                        var   (:,:,:,JXH:JZH), var_pl   (:,:,:,JXH:JZH) ) !--- [OUT]

    call OPRT_horizontalize_vec( var(:,:,:,JXH), var_pl(:,:,:,JXH), & !--- [INOUT]
                                 var(:,:,:,JYH), var_pl(:,:,:,JYH), & !--- [INOUT]
                                 var(:,:,:,JZH), var_pl(:,:,:,JZH)  ) !--- [INOUT]

    !--- calculation of Jx, Jy, and Jz
    call OPRT_gradient( GRD_vz(:,:,:,GRD_Z), GRD_vz_pl(:,:,:,GRD_Z), & !--- [IN]
                        var   (:,:,:,JX:JZ), var_pl   (:,:,:,JX:JZ)  ) !--- [OUT]

    call OPRT_horizontalize_vec( var(:,:,:,JX), var_pl(:,:,:,JX), & !--- [INOUT]
                                 var(:,:,:,JY), var_pl(:,:,:,JY), & !--- [INOUT]
                                 var(:,:,:,JZ), var_pl(:,:,:,JZ)  ) !--- [INOUT]

    !--- fill HALO
    call COMM_data_transfer( var, var_pl )

    var(suf(ADM_gmax+1,ADM_gmin-1),:,:,:) = var(suf(ADM_gmax+1,ADM_gmin),:,:,:)
    var(suf(ADM_gmin-1,ADM_gmax+1),:,:,:) = var(suf(ADM_gmin,ADM_gmax+1),:,:,:)

    !--- G^1/2 = dz/dgz
    do l = 1, ADM_lall
       !--- calculation of G^1/2 at full level
       do k = ADM_kmin, ADM_kmax
       do g = 1, ADM_gall
          GSQRT(g,k,l) = ( GRD_vz(g,k+1,l,GRD_ZH) - GRD_vz(g,k,l,GRD_ZH) ) / GRD_dgz(k)
       enddo
       enddo
       do g = 1, ADM_gall
          GSQRT(g,ADM_kmin-1,l) = GSQRT(g,ADM_kmin,l)
          GSQRT(g,ADM_kmax+1,l) = GSQRT(g,ADM_kmax,l)
       enddo

       !--- calculation of G^1/2 at half level
       do k = ADM_kmin, ADM_kmax+1
       do g = 1, ADM_gall
          GSQRTH(g,k,l) = ( GRD_vz(g,k,l,GRD_Z) - GRD_vz(g,k-1,l,GRD_Z) ) / GRD_dgzh(k)
       enddo
       enddo
       do g = 1, ADM_gall
          GSQRTH(g,ADM_kmin-1,l) = GSQRTH(g,ADM_kmin,l)
       enddo
    enddo

    !--- Gamma = (a+z) / a
    if ( deep ) then
       do l = 1, ADM_lall
       do k = 1, ADM_kall
       do g = 1, ADM_gall
          GAM (g,k,l) = 1.0_RP + GRD_vz(g,k,l,GRD_Z)  / GRD_rscale
          GAMH(g,k,l) = 1.0_RP + GRD_vz(g,k,l,GRD_ZH) / GRD_rscale
       enddo
       enddo
       enddo
    else
       do l = 1, ADM_lall
       do k = 1, ADM_kall
       do g = 1, ADM_gall
          GAM (g,k,l) = 1.0_RP
          GAMH(g,k,l) = 1.0_RP
       enddo
       enddo
       enddo
    endif

    do l = 1, ADM_lall
    do k = 1, ADM_kall
    do g = 1, ADM_gall
       VMTR_GAM2    (g,k,l) = GAM (g,k,l) * GAM (g,k,l)
       VMTR_GAM2H   (g,k,l) = GAMH(g,k,l) * GAMH(g,k,l)
       VMTR_GSGAM2  (g,k,l) = GAM (g,k,l) * GAM (g,k,l) * GSQRT (g,k,l)
       VMTR_GSGAM2H (g,k,l) = GAMH(g,k,l) * GAMH(g,k,l) * GSQRTH(g,k,l)

       VMTR_RGSQRTH (g,k,l) = 1.0_RP / GSQRTH(g,k,l)
       VMTR_RGAM    (g,k,l) = 1.0_RP / GAM (g,k,l)
       VMTR_RGAMH   (g,k,l) = 1.0_RP / GAMH(g,k,l)
       VMTR_RGSGAM2 (g,k,l) = 1.0_RP / VMTR_GSGAM2 (g,k,l)
       VMTR_RGSGAM2H(g,k,l) = 1.0_RP / VMTR_GSGAM2H(g,k,l)
    enddo
    enddo
    enddo

    ! full level <-> half level interpolation factor
    do l = 1, ADM_lall
       do k = ADM_kmin, ADM_kmax+1
       do g = 1, ADM_gall
          VMTR_C2Wfact(g,k,I_a,l) = 0.5_RP * GRD_afac(k) * VMTR_RGSGAM2(g,k  ,l) * VMTR_GSGAM2H(g,k,l)
          VMTR_C2Wfact(g,k,I_b,l) = 0.5_RP * GRD_bfac(k) * VMTR_RGSGAM2(g,k-1,l) * VMTR_GSGAM2H(g,k,l)
       enddo
       enddo
       do g = 1, ADM_gall
          VMTR_C2Wfact(g,ADM_kmin-1,I_a,l) = 0.0_RP
          VMTR_C2Wfact(g,ADM_kmin-1,I_b,l) = 0.0_RP
       enddo

       do k = ADM_kmin-1, ADM_kmax
       do g = 1, ADM_gall
          VMTR_W2Cfact(g,k,I_c,l) = 0.5_RP * GRD_cfac(k) * VMTR_GSGAM2(g,k,l) * VMTR_RGSGAM2H(g,k+1,l)
          VMTR_W2Cfact(g,k,I_d,l) = 0.5_RP * GRD_dfac(k) * VMTR_GSGAM2(g,k,l) * VMTR_RGSGAM2H(g,k  ,l)
       enddo
       enddo
       do g = 1, ADM_gall
          VMTR_W2Cfact(g,ADM_kmax+1,I_c,l) = 0.0_RP
          VMTR_W2Cfact(g,ADM_kmax+1,I_d,l) = 0.0_RP
       enddo
    enddo

    ! full level <-> half level interpolation factor with Gz

    !--- Gz(X) = - JX / G^1/2
    !--- Gz(Y) = - JY / G^1/2
    !--- Gz(Z) = - JZ / G^1/2
    do l = 1, ADM_lall
       do k = 1, ADM_kall
       do g = 1, ADM_gall
          GZXH(g,k,l) = -var(g,k,l,JXH) / GSQRTH(g,k,l)
          GZYH(g,k,l) = -var(g,k,l,JYH) / GSQRTH(g,k,l)
          GZZH(g,k,l) = -var(g,k,l,JZH) / GSQRTH(g,k,l)
          GZX (g,k,l) = -var(g,k,l,JX)  / GSQRT (g,k,l)
          GZY (g,k,l) = -var(g,k,l,JY)  / GSQRT (g,k,l)
          GZZ (g,k,l) = -var(g,k,l,JZ)  / GSQRT (g,k,l)
       enddo
       enddo

       do k = ADM_kmin, ADM_kmax+1
       do g = 1, ADM_gall
          VMTR_C2WfactGz(g,k,I_a_GZXH,l) = 0.5_RP * GRD_afac(k) * VMTR_RGSGAM2(g,k  ,l) * VMTR_GSGAM2H(g,k,l) &
                                         * GZXH(g,k,l)
          VMTR_C2WfactGz(g,k,I_b_GZXH,l) = 0.5_RP * GRD_bfac(k) * VMTR_RGSGAM2(g,k-1,l) * VMTR_GSGAM2H(g,k,l) &
                                         * GZXH(g,k,l)
          VMTR_C2WfactGz(g,k,I_a_GZYH,l) = 0.5_RP * GRD_afac(k) * VMTR_RGSGAM2(g,k  ,l) * VMTR_GSGAM2H(g,k,l) &
                                         * GZYH(g,k,l)
          VMTR_C2WfactGz(g,k,I_b_GZYH,l) = 0.5_RP * GRD_bfac(k) * VMTR_RGSGAM2(g,k-1,l) * VMTR_GSGAM2H(g,k,l) &
                                         * GZYH(g,k,l)
          VMTR_C2WfactGz(g,k,I_a_GZZH,l) = 0.5_RP * GRD_afac(k) * VMTR_RGSGAM2(g,k  ,l) * VMTR_GSGAM2H(g,k,l) &
                                         * GZZH(g,k,l)
          VMTR_C2WfactGz(g,k,I_b_GZZH,l) = 0.5_RP * GRD_bfac(k) * VMTR_RGSGAM2(g,k-1,l) * VMTR_GSGAM2H(g,k,l) &
                                         * GZZH(g,k,l)
       enddo
       enddo

       do g = 1, ADM_gall
          VMTR_C2WfactGz(g,ADM_kmin-1,I_a_GZXH,l) = 0.0_RP
          VMTR_C2WfactGz(g,ADM_kmin-1,I_b_GZXH,l) = 0.0_RP
          VMTR_C2WfactGz(g,ADM_kmin-1,I_a_GZYH,l) = 0.0_RP
          VMTR_C2WfactGz(g,ADM_kmin-1,I_b_GZYH,l) = 0.0_RP
          VMTR_C2WfactGz(g,ADM_kmin-1,I_a_GZZH,l) = 0.0_RP
          VMTR_C2WfactGz(g,ADM_kmin-1,I_b_GZZH,l) = 0.0_RP
       enddo
    enddo

    !--- calculation of volume, geopotential
    do l = 1, ADM_lall
    do k = 1, ADM_kall
    do g = 1, ADM_gall
       VMTR_VOLUME(g,k,l) = GMTR_area(g,l) * VMTR_GSGAM2(g,k,l) * GRD_dgz(k)

       VMTR_PHI   (g,k,l) = GRD_vz(g,k,l,GRD_Z) * CONST_GRAV
    enddo
    enddo
    enddo

    if ( ADM_have_pl ) then

       !---   G^1/2 = dz/dgz
       do l = 1, ADM_lall_pl
          !--- calculation of G^1/2 at full level
          do k = ADM_kmin, ADM_kmax
          do g = 1, ADM_gall_pl
             GSQRT_pl(g,k,l) = ( GRD_vz_pl(g,k+1,l,GRD_ZH) - GRD_vz_pl(g,k,l,GRD_ZH) ) / GRD_dgz(k)
          enddo
          enddo
          do g = 1, ADM_gall_pl
             GSQRT_pl(g,ADM_kmin-1,l) = GSQRT_pl(g,ADM_kmin,l)
             GSQRT_pl(g,ADM_kmax+1,l) = GSQRT_pl(g,ADM_kmax,l)
          enddo
          !--- calculation of G^1/2 at half level
          do k = ADM_kmin, ADM_kmax+1
          do g = 1, ADM_gall_pl
             GSQRTH_pl(g,k,l) = ( GRD_vz_pl(g,k,l,GRD_Z) - GRD_vz_pl(g,k-1,l,GRD_Z) ) / GRD_dgzh(k)
          enddo
          enddo
          do g = 1, ADM_gall_pl
             GSQRTH_pl(g,ADM_kmin-1,l) = GSQRTH_pl(g,ADM_kmin,l)
          enddo
       enddo

       !--- Gamma = (a+z) / a
       if ( deep ) then
          do l = 1, ADM_lall_pl
          do k = 1, ADM_kall
          do g = 1, ADM_gall_pl
             GAM_pl (g,k,l) = 1.0_RP + GRD_vz_pl(g,k,l,GRD_Z)  / GRD_rscale
             GAMH_pl(g,k,l) = 1.0_RP + GRD_vz_pl(g,k,l,GRD_ZH) / GRD_rscale
          enddo
          enddo
          enddo
       else
          do l = 1, ADM_lall_pl
          do k = 1, ADM_kall
          do g = 1, ADM_gall_pl
             GAM_pl (g,k,l) = 1.0_RP
             GAMH_pl(g,k,l) = 1.0_RP
          enddo
          enddo
          enddo
       endif

       do l = 1, ADM_lall_pl
       do k = 1, ADM_kall
       do g = 1, ADM_gall_pl
          VMTR_GAM2_pl    (g,k,l) = GAM_pl (g,k,l) * GAM_pl (g,k,l)
          VMTR_GAM2H_pl   (g,k,l) = GAMH_pl(g,k,l) * GAMH_pl(g,k,l)
          VMTR_GSGAM2_pl  (g,k,l) = GAM_pl (g,k,l) * GAM_pl (g,k,l) * GSQRT_pl (g,k,l)
          VMTR_GSGAM2H_pl (g,k,l) = GAMH_pl(g,k,l) * GAMH_pl(g,k,l) * GSQRTH_pl(g,k,l)

          VMTR_RGSQRTH_pl (g,k,l) = 1.0_RP / GSQRTH_pl(g,k,l)
          VMTR_RGAM_pl    (g,k,l) = 1.0_RP / GAM_pl (g,k,l)
          VMTR_RGAMH_pl   (g,k,l) = 1.0_RP / GAMH_pl(g,k,l)
          VMTR_RGSGAM2_pl (g,k,l) = 1.0_RP / VMTR_GSGAM2_pl (g,k,l)
          VMTR_RGSGAM2H_pl(g,k,l) = 1.0_RP / VMTR_GSGAM2H_pl(g,k,l)
       enddo
       enddo
       enddo

       ! full level <-> half level interpolation factor
       do l = 1, ADM_lall_pl
          do k = ADM_kmin, ADM_kmax+1
          do g = 1, ADM_gall_pl
             VMTR_C2Wfact_pl(g,k,I_a,l) = 0.5_RP * GRD_afac(k) * VMTR_RGSGAM2_pl(g,k  ,l) * VMTR_GSGAM2H_pl(g,k,l)
             VMTR_C2Wfact_pl(g,k,I_b,l) = 0.5_RP * GRD_bfac(k) * VMTR_RGSGAM2_pl(g,k-1,l) * VMTR_GSGAM2H_pl(g,k,l)
          enddo
          enddo
          do g = 1, ADM_gall_pl
             VMTR_C2Wfact_pl(g,ADM_kmin-1,I_a,l) = 0.0_RP
             VMTR_C2Wfact_pl(g,ADM_kmin-1,I_b,l) = 0.0_RP
          enddo

          do k = ADM_kmin-1, ADM_kmax
          do g = 1, ADM_gall_pl
             VMTR_W2Cfact_pl(g,k,I_c,l) = 0.5_RP * GRD_cfac(k) * VMTR_GSGAM2_pl(g,k,l) * VMTR_RGSGAM2H_pl(g,k+1,l)
             VMTR_W2Cfact_pl(g,k,I_d,l) = 0.5_RP * GRD_dfac(k) * VMTR_GSGAM2_pl(g,k,l) * VMTR_RGSGAM2H_pl(g,k  ,l)
          enddo
          enddo
          do g = 1, ADM_gall_pl
             VMTR_W2Cfact_pl(g,ADM_kmax+1,I_c,l) = 0.0_RP
             VMTR_W2Cfact_pl(g,ADM_kmax+1,I_d,l) = 0.0_RP
          enddo
       enddo

       ! full level <-> half level interpolation factor with Gz

       !--- Gz(X) = - JX / G^1/2
       !--- Gz(Y) = - JY / G^1/2
       !--- Gz(Z) = - JZ / G^1/2
       do l = 1, ADM_lall_pl
          do k = 1, ADM_kall
          do g = 1, ADM_gall_pl
             GZXH_pl(g,k,l) = -var_pl(g,k,l,JXH) / GSQRTH_pl(g,k,l)
             GZYH_pl(g,k,l) = -var_pl(g,k,l,JYH) / GSQRTH_pl(g,k,l)
             GZZH_pl(g,k,l) = -var_pl(g,k,l,JZH) / GSQRTH_pl(g,k,l)
             GZX_pl (g,k,l) = -var_pl(g,k,l,JX)  / GSQRT_pl (g,k,l)
             GZY_pl (g,k,l) = -var_pl(g,k,l,JY)  / GSQRT_pl (g,k,l)
             GZZ_pl (g,k,l) = -var_pl(g,k,l,JZ)  / GSQRT_pl (g,k,l)
          enddo
          enddo

          do k = ADM_kmin, ADM_kmax+1
          do g = 1, ADM_gall_pl
             VMTR_C2WfactGz_pl(g,k,I_a_GZXH,l) = 0.5_RP * GRD_afac(k) * VMTR_RGSGAM2_pl(g,k  ,l) * VMTR_GSGAM2H_pl(g,k,l) &
                                               * GZXH_pl(g,k,l)
             VMTR_C2WfactGz_pl(g,k,I_b_GZXH,l) = 0.5_RP * GRD_bfac(k) * VMTR_RGSGAM2_pl(g,k-1,l) * VMTR_GSGAM2H_pl(g,k,l) &
                                               * GZXH_pl(g,k,l)
             VMTR_C2WfactGz_pl(g,k,I_a_GZYH,l) = 0.5_RP * GRD_afac(k) * VMTR_RGSGAM2_pl(g,k  ,l) * VMTR_GSGAM2H_pl(g,k,l) &
                                               * GZYH_pl(g,k,l)
             VMTR_C2WfactGz_pl(g,k,I_b_GZYH,l) = 0.5_RP * GRD_bfac(k) * VMTR_RGSGAM2_pl(g,k-1,l) * VMTR_GSGAM2H_pl(g,k,l) &
                                               * GZYH_pl(g,k,l)
             VMTR_C2WfactGz_pl(g,k,I_a_GZZH,l) = 0.5_RP * GRD_afac(k) * VMTR_RGSGAM2_pl(g,k  ,l) * VMTR_GSGAM2H_pl(g,k,l) &
                                               * GZZH_pl(g,k,l)
             VMTR_C2WfactGz_pl(g,k,I_b_GZZH,l) = 0.5_RP * GRD_bfac(k) * VMTR_RGSGAM2_pl(g,k-1,l) * VMTR_GSGAM2H_pl(g,k,l) &
                                               * GZZH_pl(g,k,l)
          enddo
          enddo

          do g = 1, ADM_gall_pl
             VMTR_C2WfactGz_pl(g,ADM_kmin-1,I_a_GZXH,l) = 0.0_RP
             VMTR_C2WfactGz_pl(g,ADM_kmin-1,I_b_GZXH,l) = 0.0_RP
             VMTR_C2WfactGz_pl(g,ADM_kmin-1,I_a_GZYH,l) = 0.0_RP
             VMTR_C2WfactGz_pl(g,ADM_kmin-1,I_b_GZYH,l) = 0.0_RP
             VMTR_C2WfactGz_pl(g,ADM_kmin-1,I_a_GZZH,l) = 0.0_RP
             VMTR_C2WfactGz_pl(g,ADM_kmin-1,I_b_GZZH,l) = 0.0_RP
          enddo
       enddo

       !--- calculation of volume, geopotential
       do l = 1, ADM_lall_pl
       do k = 1, ADM_kall
       do g = 1, ADM_gall_pl
          VMTR_VOLUME_pl(g,k,l) = GMTR_area_pl(g,l) * VMTR_GSGAM2_pl(g,k,l) * GRD_dgz(k)

          VMTR_PHI_pl   (g,k,l) = GRD_vz_pl(g,k,l,GRD_Z) * CONST_GRAV
       enddo
       enddo
       enddo
    else
       VMTR_GAM2_pl     (:,:,:)   = 0.0_RP
       VMTR_GAM2H_pl    (:,:,:)   = 0.0_RP
       VMTR_GSGAM2_pl   (:,:,:)   = 0.0_RP
       VMTR_GSGAM2H_pl  (:,:,:)   = 0.0_RP
       VMTR_RGSQRTH_pl  (:,:,:)   = 0.0_RP
       VMTR_RGAM_pl     (:,:,:)   = 0.0_RP
       VMTR_RGAMH_pl    (:,:,:)   = 0.0_RP
       VMTR_RGSGAM2_pl  (:,:,:)   = 0.0_RP
       VMTR_RGSGAM2H_pl (:,:,:)   = 0.0_RP
       VMTR_W2Cfact_pl  (:,:,:,:) = 0.0_RP
       VMTR_C2Wfact_pl  (:,:,:,:) = 0.0_RP
       VMTR_C2WfactGz_pl(:,:,:,:) = 0.0_RP
       VMTR_VOLUME_pl   (:,:,:)   = 0.0_RP
       VMTR_PHI_pl      (:,:,:)   = 0.0_RP
    endif

    return
  end subroutine VMTR_setup

  !-----------------------------------------------------------------------------
  integer function suf(i,j)
    use mod_adm, only: &
       ADM_gall_1d
    implicit none

    integer :: i, j
    !---------------------------------------------------------------------------

    suf = ADM_gall_1d * (j-1) + i

  end function suf

end module mod_vmtr


text:
! { dg-do run }
! { dg-options "-finit-character=32" }

program init_flag_6
  call char_test
end program init_flag_6

! Test some initializations for both implicitly and
! explicitly declared local variables.
subroutine char_test
  character*1 c1
  character*8 c2, c3(5)
  character c4(10)
  if (c1 /= ' ') call abort
  if (c2 /= '        ') call abort
  if (c3(1) /= '        ') call abort
  if (c3(5) /= '        ') call abort
  if (c4(5) /= ' ') call abort
end subroutine char_test

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
      program main
      implicit none
! number of components
      integer ncom
      parameter (ncom=4)
! number of particles
      integer npart
      parameter (npart=2)
! Detail of particles
      integer partlst(npart)
      data partlst/-11,11/
! number of beam type
      integer nbeam
      parameter (nbeam=2)
! detail of beam id
      integer beamlst(nbeam)
      data beamlst/-11,11/
! nodes of the y grid
      integer nny
      parameter (nny=100)
! nodes of the z grid
      integer nnz
      parameter (nnz=20)
! the reference scale
      real*8 Qref
! the scale (temporary variable)
      real*8 Q
! the electron mass
      real*8 me
      parameter (me=0.511d-3)
! the electric coupling
      real*8 aem
! this is the Gmu scheme
      parameter (aem=1.d0/132.2332d0)
! grid nodes and grids
      real*8 yv(nny),zv(nnz)
      real*8,dimension(:,:,:,:,:), allocatable :: grids
c
      integer iny,inz,icom,ipart,ibeam
      integer i1,i2,i,j
c     integer ifakegrid
      real*8 ylow,yupp,zlow,zupp
c     real*8 alst,q0st,alxm,q0xm
      real*8 jkb,x
      real*8 eepdf_tilde_calc,eepdf_tilde_factor
!     temporary variables
      real*8 tv1,tv2,tvy
c
      ylow=1.d-6
      yupp=1d0-1.d-8
      zlow=log(1.d0/me)
      zupp=log(10000.d0/me)
      Qref=1d0
      allocate(grids(nny,nnz,ncom,npart,nbeam))
c     electron part
      do iny=1,nny
c     Using linear grid as the first step
c     More points at the end
      tvy=dfloat(iny-1)/dfloat(nny-1)
c     yv(iny)=(1.d0-(1.d0-tvy)**3.d0)*(yupp-ylow)+ylow
      yv(iny)=tvy*tvy*(6.d0+tvy*(-8.d0+3.d0*tvy))*(yupp-ylow)+ylow
c     yv(iny)=tvy*(3.d0-tvy*2.d0)*(yupp-ylow)+ylow
      do inz=1,nnz
      zv(inz)=dfloat(inz-1)/dfloat(nnz-1)*(zupp-zlow)+zlow
      do icom=1,ncom
      do ipart=1,npart
      do ibeam=1,nbeam
      Q=exp(zv(inz))*me
      tv1=eepdf_tilde_calc(yv(iny),Q*Q,icom
     #     ,partlst(ipart),beamlst(ibeam))
      tv2=eepdf_tilde_factor(yv(iny),Q*Q,icom
     #     ,partlst(ipart),beamlst(ibeam))
      grids(iny,inz,icom,ipart,ibeam)=tv1/tv2
!     print *,yv(iny),tv1,tv2,tv1/tv2
      enddo
      enddo
      enddo
      enddo
      enddo
c
      open(unit=10,file='eepdf.f',status='unknown')
c
      write(10,'(a)')
     #'      function eepdf_tilde(y,Q2,icom,ipart,ibeam)',
     #'      implicit none',
     #'      real*8 eepdf_tilde',
     #'      real*8 Q2,Qref,me',
     #'      integer icom,ipart,ibeam',
     #'      real*8 tmp,cstmin,cxmmin,cxmmax',
     #'      integer i,id0,listmin,lixmmin,lixmmax',
     #'      logical firsttime,check,T,F,grid(21)',
     #'      parameter (T=.true.)',
     #'      parameter (F=.false.)',
     #'      real*8 eepdf_tilde_factor'
      write(10,'(a)')
     #'      real*8 y,z',
     #'      real*8 ylow,yupp,zlow,zupp',
     #'      real*8 jkb'
      write(10,'(a,d15.8,a,d15.8,a)')
     #'      parameter (ylow=',ylow,',yupp=',yupp,')'
      write(10,'(a,d15.8,a,d15.8,a)')
     #'      parameter (zlow=',zlow,',zupp=',zupp,')'
      write(10,'(a,d15.8,a,d15.8,a)')
     #'      parameter (Qref=',Qref,',me=',me,')'
      do ipart=1,npart
      do ibeam=1,nbeam
      do icom=1,ncom
!     Use ipart instead of the real PID, as we have negative PID
        write(10,'(a,i0,a,i0,a,i0)')
     #'      real*8 eepdf_',icom,'_',ipart,'_',ibeam
      enddo
      enddo
      enddo
      write(10,'(a)')
     #'      z=0.5d0*log(Q2/me/me)'
!     No check now
      do icom=1,ncom
      if(icom.eq.1)then
        write(10,'(a)')
     #'      if(icom.eq.1)then'
      else
        write(10,'(a,i0,a)')
     #'      else if(icom.eq.',icom,')then'
      endif
        do ipart=1,npart
        if(ipart.eq.1)then
          write(10,'(a,i0,a)')
     #'        if(ipart.eq.',partlst(ipart),')then'
        else
          write(10,'(a,i0,a)')
     #'        else if(ipart.eq.',partlst(ipart),')then'
        endif
        do ibeam=1,nbeam
        if(ibeam.eq.1)then
          write(10,'(a,i0,a)')
     #'          if(ibeam.eq.',beamlst(ibeam),')then'
        else
          write(10,'(a,i0,a)')
     #'          else if(ibeam.eq.',beamlst(ibeam),')then'
        endif
        write(10,'(a,i0,a,i0,a,i0,a)')
     #'          tmp=eepdf_',icom,'_',ipart,'_',ibeam,'(y,z)'
        enddo
        write(10,'(a)')
     #'          else',
     #'            tmp=0d0',
     #'          endif'

        enddo
        write(10,'(a)')
     #'        else',
     #'          tmp=0d0',
     #'        endif'
        enddo
        write(10,'(a)')
     #'      else',
     #'        tmp=0d0',
     #'      endif',
     #'      eepdf_tilde=tmp*eepdf_tilde_factor(y,Q2,icom,ipart,ibeam)',
     #'      end'
        write(10,'(a)')
     #'c',
     #'c',
     #'cccc',
     #'c',
     #'c'

         do icom=1,ncom
         do ipart=1,npart
         do ibeam=1,nbeam

         write(10,'(a,i0,a,i0,a,i0,a)')
     #'      function eepdf_',icom,'_',ipart,'_',ibeam,'(y,z)'
         write(10,'(a)')
     #'      implicit none'
         write(10,'(a,i1,a,i1,a,i0,a)')
     #'      real*8 eepdf_',icom,'_',ipart,'_',ibeam,',y,z'
         write(10,'(a)')
     #'      integer narg,nny,nnz',
     #'      parameter (narg=2)'
        write(10,'(a,i0,a)')
     #'      parameter (nny=',nny,')'
        write(10,'(a,i0,a)')
     #'      parameter (nnz=',nnz,')'
        write(10,'(a)')
     #'      integer iny,inz,nent(narg)',
     #'      real*8 tmp,dfint,ymap,zmap',
     #'      real*8 arg(narg),ent(nny+nnz)',
     #'      real*8 yv(nny),zv(nnz),gridv(nny,nnz)',
     #'      logical firsttime',
     #'      external dfint,ymap,zmap'
c
        write(10,'(a)')
     #'      data yv/'
        i1=nny/4
        i2=mod(nny,4)
        if(i2.eq.0)i1=i1-1
        do i=1,i1
        write(10,'(a,4(d15.8,1h,))')
     #'     #',(yv(j),j=1+(i-1)*4,i*4)
        enddo
        write(10,'(a,4(d15.8,a))')
     #'     #',(yv(j),',',j=1+i1*4,nny-1),yv(nny),'/'
c
        write(10,'(a)')
     #'      data zv/'
        i1=nnz/4
        i2=mod(nnz,4)
        if(i2.eq.0)i1=i1-1
        do i=1,i1
        write(10,'(a,4(d15.8,1h,))')
     #'     #',(zv(j),j=1+(i-1)*4,i*4)
        enddo
        write(10,'(a,4(d15.8,a))')
     #'     #',(zv(j),',',j=1+i1*4,nnz-1),zv(nnz),'/'
c
        do inz=1,nnz
        write(10,'(a,i3,a,i3,a)')
     #'      data (gridv(iny,',inz,'),iny=1,',nny,')/'
        i1=nny/4
        i2=mod(nny,4)
        if(i2.eq.0)i1=i1-1
        do i=1,i1
        write(10,'(a,4(d15.8,1h,))')
     #'     #',(grids(j,inz,icom,ipart,ibeam),j=1+(i-1)*4,i*4)
        enddo
        write(10,'(a,4(d15.8,a))')
     #'     #',(grids(j,inz,icom,ipart,ibeam),',',j=1+i1*4,nny-1),
     #    grids(nny,inz,icom,ipart,ibeam),'/'
        enddo
c
        write(10,'(a)')
     #'      data firsttime/.true./',
     #'      save',
     #'c',
     #'      if(firsttime)then',
     #'        firsttime=.false.',
     #'        nent(1)=nny',
     #'        nent(2)=nnz',
     #'        do iny=1,nny',
     #'          ent(iny)=ymap(yv(iny))',
     #'        enddo',
     #'        do inz=1,nnz',
     #'          ent(nny+inz)=zmap(zv(inz))',
     #'        enddo',
     #'      endif',
     #'      arg(1)=ymap(y)',
     #'      arg(2)=zmap(z)',
     #'      tmp=dfint(narg,arg,nent,ent,gridv)'
        write(10,'(a,i0,a,i0,a,i0,a)')
     #'      eepdf_',icom,'_',ipart,'_',ibeam,'=tmp'
        write(10,'(a)')
     #'      return',
     #'      end'
c
        write(10,'(a)')
     #'c',
     #'c',
     #'cccc',
     #'c',
     #'c'
c end of loops over itype and ipart
        enddo
        enddo
        enddo
c
        write(10,'(a)')
     #'      function ymap(st)',
     #'c Use this function to interpolate by means of',
     #'c   stnode_i=ymap(stnode_stored_i).',
     #'c Example (to be used below): tmp=log10(st)',
     #'      implicit none',
     #'      real*8 ymap,st,tmp',
     #'c',
     #'      tmp=st',
     #'      ymap=tmp',
     #'      return',
     #'      end',
     #'  ',
     #'  ',
     #'      function zmap(xm)',
     #'c Use this function to interpolate by means of',
     #'c   xmnode_i=zmap(xmnode_stored_i).',
     #'c Example (to be used below): tmp=log10(xm)',
     #'      implicit none',
     #'      real*8 zmap,xm,tmp',
     #'c',
     #'      tmp=xm',
     #'      zmap=tmp',
     #'      return',
     #'      end'
c
        end program


!! unit tests and registration tests of HDF5 OO interface
use, intrinsic:: ieee_arithmetic, only: ieee_value, ieee_quiet_nan, ieee_is_nan
use, intrinsic:: iso_fortran_env, only: int64, int32, real32, real64, stderr=>error_unit
use, intrinsic:: iso_c_binding, only: c_null_char
use hdf5_interface, only: hdf5_file, toLower, strip_trailing_null

implicit none

type(hdf5_file) :: h5f
integer :: i1(4)
real(real32)    :: nan, r1(4), r2(4,4)

integer :: i

nan = ieee_value(1.0, ieee_quiet_nan)

do concurrent (i = 1:size(i1))
  i1(i) = i
enddo

r1 = i1

call test_lowercase()
print *,'PASSED: HDF5 character'
call test_strip_null()
print *,'PASSED: null strip'
call testNewHDF5()
print *,'PASSED: HDF5 scalar real / integer'
call testGroup()
print *,'PASSED: HDF5 group'
call testAddHDF5()
print *,'PASSED: HDF5 array'
call test_hdf5_deflate()
print *,'PASSED: HDF5 compression'
call test_write_attributes()
print *,'PASSED: HDF5 attributes'

call testrwHDF5(ng=69, nn=100, pn=5)
call test_string_rw()

contains

subroutine test_lowercase()

character(*), parameter :: hello = 'HeLl0 Th3rE !>? '
  !! Fortran 2003 allocatable string

if (.not.(toLower(hello)=='hell0 th3re !>? ')) error stop 'error: lowercase conversion'

if (.not.(trim(toLower(hello))=='hell0 th3re !>?')) error stop 'Allocatable lowercase conversion error'

if(.not.all(toLower(['Hi','hI'])==['hi','hi'])) error stop 'error on array conversion'

end subroutine test_lowercase

subroutine test_strip_null()

character(*), parameter :: hello = 'HeLl0 Th3rE !>? '

if (.not.strip_trailing_null(hello // c_null_char) == hello) error stop 'problem stripping trailing null'

end subroutine test_strip_null

subroutine testNewHDF5()

real(real32), allocatable :: rr1(:)
real(real32) :: rt
integer(int32) :: it
integer(int32), allocatable :: i1t(:)

call h5f%initialize('test.h5',status='new',action='w')

!! scalar tests
call h5f%add('/scalar_int', 42_int32)
call h5f%add('/scalar_real', 42._real32)

call h5f%get('/scalar_int', it)
call h5f%get('/scalar_real', rt)
if (.not.(rt==it .and. it==42)) then
  write(stderr,*) it,'/=',rt
  error stop 'scalar real / int: not equal 42'
endif

call h5f%add('/real1',r1)
call h5f%get('/real1',rr1)
if (.not.all(r1 == rr1)) error stop 'real: read does not match write'

call h5f%add('/ai1', i1)
call h5f%get('/ai1',i1t)

call h5f%finalize()

if (.not.all(i1==i1t)) error stop 'integer 1-D: read does not match write'

if (.not. h5f%filename == 'test.h5') error stop h5f%filename//' mismatch filename'

end subroutine testNewHDF5

subroutine testGroup

call h5f%initialize('test_groups.h5',status='new',action='rw')
call h5f%add('/test/')

call h5f%open('/test')
call h5f%add('group3/scalar', 1_int32)
call h5f%add('group3/scalar_real', 1._real32)
call h5f%close()

call h5f%finalize()

end subroutine testGroup

subroutine testAddHDF5()
!! tests that compression doesn't fail for very small datasets, where it really shouldn't be used (makes file bigger)

integer :: i2(4,4)
integer, allocatable :: i2t(:,:)
real(real32), allocatable :: rr2(:,:)
real(real32)  ::  nant

i2(1,:) = i1
do concurrent (i = 1:size(i2,2))
  i2(i,:) = i2(1,:) * i
enddo

r2 = i2

call h5f%initialize('test.h5',status='old',action='rw',comp_lvl=1)

call h5f%add('/test/group2/ai2', i2)
call h5f%get('/test/group2/ai2',i2t)
if (.not.all(i2==i2t)) error stop 'read does not match write'

call h5f%add('/test/real2',r2)
call h5f%get('/test/real2',rr2)
if (.not.all(r2 == rr2)) error stop 'real: read does not match write'

call h5f%add('/nan',nan)
call h5f%get('/nan',nant)
if (.not.ieee_is_nan(nant)) error stop 'failed storing or reading NaN'

call h5f%finalize()

end subroutine testAddHDF5

subroutine test_hdf5_deflate()

integer(int64), parameter :: N=1000
integer(int64) :: crat
integer ::  fsize, ibig2(N,N) = 0, ibig3(N,N,4) = 0
real(real32) :: big2(N,N) = 0., big3(N,N,4) = 0.

call h5f%initialize('test_deflate.h5',status='new',action='rw',comp_lvl=1)

call h5f%add('/big2', big2, chunk_size=[100,100])

call h5f%finalize()

inquire(file='test_deflate.h5', size=fsize)
crat = (N*N*storage_size(big2)/8) / fsize

print '(A,F6.2,A,I6)','filesize (Mbytes): ',fsize/1e6, '   2D compression ratio:',crat

if (h5f%comp_lvl > 0 .and. crat < 10) print *,'warning: 2D low compression'

!======================================
call h5f%initialize('test_deflate.h5',status='new',action='rw',comp_lvl=1)

call h5f%add('/big3', big3, chunk_size=[100,100,1])

call h5f%finalize()

inquire(file='test_deflate.h5', size=fsize)
crat = (N*N*storage_size(big3)/8) / fsize

print '(A,F6.2,A,I6)','filesize (Mbytes): ',fsize/1e6, '   3D compression ratio:',crat

if (h5f%comp_lvl > 0 .and. crat < 10) print *,'warning: 3D low compression'
!======================================
call h5f%initialize('test_deflate.h5',status='new',action='rw',comp_lvl=1)

call h5f%add('/ibig3', ibig3, chunk_size=[1000,100,1])

call h5f%finalize()

inquire(file='test_deflate.h5', size=fsize)
crat = (N*N*storage_size(ibig3)/8) / fsize

print '(A,F6.2,A,I6)','filesize (Mbytes): ',fsize/1e6, '   3D compression ratio:',crat

if (h5f%comp_lvl > 0 .and. crat < 10) print *,'warning: 3D low compression'
!======================================
call h5f%initialize('test_deflate.h5',status='new',action='rw',comp_lvl=1)

call h5f%add('/ibig2', ibig2, chunk_size=[100,100])

call h5f%finalize()

inquire(file='test_deflate.h5', size=fsize)
crat = (N*N*storage_size(ibig2)/8) / fsize

print '(A,F6.2,A,I6)','filesize (Mbytes): ',fsize/1e6, '   3D compression ratio:',crat

if (h5f%comp_lvl > 0 .and. crat < 10) print *,'warning: 3D low compression'

end subroutine test_hdf5_deflate

subroutine test_write_attributes()

call h5f%initialize('test_deflate.h5')

call h5f%writeattr('/little/','note','this is just a little number')

call h5f%finalize()

end subroutine test_write_attributes

subroutine test_string_rw()

character(2) :: value

call h5f%initialize('test_string.h5',status='new')

call h5f%add('/little','42')

call h5f%get('/little',value)

if (.not.value=='42') error stop 'string dataset read/write verification failure. Value: '// value

call h5f%finalize()

end subroutine test_string_rw

subroutine testrwHDF5(ng, nn, pn)
!! more group

integer, intent(in) :: ng, nn, pn

real(real32), allocatable :: flux(:,:),fo(:)
character(2) :: pnc,ic
integer :: i

allocate(flux(nn,ng),fo(nn))
flux = 1.0
write(pnc,'(I2)') pn

call h5f%initialize('p'//trim(adjustl(pnc))//'.h5',status='new',action='w')

do i = 1,ng
write(ic,'(I2)') i
call h5f%add('/group'//trim(adjustl(ic))//'/flux_node',flux(:ng,i))
enddo

call h5f%get('/group1/flux_node',fo)
if (.not.all(fo(:ng)==flux(:ng,1))) error stop 'test_read_write: read does not match write'

call h5f%finalize()

end subroutine testrwHDF5

end program


text:
SUBROUTINE h3kalk(k,omega,n,nl,rho,g,A33,B33,C33,L,Br,D,M,H3)
IMPLICIT NONE
!Erklærer variable og vektorer
INTEGER::n,nl,i,j
DOUBLE PRECISION,DIMENSION(1:n)::omega,k
REAL,DIMENSION(1:n,1:nl)::A33,B33
DOUBLE PRECISION,DIMENSION(1:n,1:nl)::eta3nevner,F3A1,F3A2,F3A
DOUBLE PRECISION,DIMENSION(1:n,1:nl)::eta3A,H3
REAL,DIMENSION(1:nl)::L,Br,D,M,C33
REAL rho,g,a,b

!Beregner integralene til kreftene i hiv, skriver kontroll til skjerm
DO i=1,n
   DO j=1,nl
      F3A1(i,j) = (2.0/k(i)) * SIN((k(i) * L(j)) / 2.0 )

      F3A2(i,j) = rho*g*Br(j)* EXP(- k(i)*D(j)) - (omega(i)**2) * (A33(i,j)/L(i,j)) * &
           EXP(- k(i)*D(j)/2.0)

      eta3nevner(i,j)= SQRT(  (C33(j) - (M(j)+A33(i,jy)) * omega(i)**2 )**2 + &
           (B33(i,j) * omega(i))**2  )

      F3A(i,j) = F3A1(i,j) * F3A2(i,j)
   END DO
END DO

!Beregner eta 3
DO i=1,n
   DO j=1,nl
      eta3A(i,j)=   F3A(i,j) / eta3nevner(i,j)
   END DO
END DO

!Beregner transfer-funksjoner
DO i=1,n
   DO j=1,nl
      H3(i,j)=ABS(eta3A(i,j))
   END DO
END DO

!Skriver resultatet til fil for sjekk/ plotting
OPEN(41,FILE="eta3plot.txt",ACTION="WRITE",STATUS="UNKNOWN")
DO i=1,n
   WRITE(41,*) omega(i),eta3A(i,1),eta3A(i,2)
END DO
CLOSE(41)

RETURN
END SUBROUTINE

max_stars_repo_name: torlarse/dmet

max_issues_repo_name: torlarse/dmet

max_forks_repo_name: torlarse/dmet