text:
!
!     (c) 2019 Guide Star Engineering, LLC
!     This Software was developed for the US Nuclear Regulatory Commission (US NRC)
!     under contract "Multi-Dimensional Physics Implementation into Fuel Analysis under
!     Steady-state and Transients (FAST)", contract # NRC-HQ-60-17-C-0007
!
!
!    NEMO - Numerical Engine (for) Multiphysics Operators
! Copyright (c) 2007, Stefano Toninel
!                     Gian Marco Bianchi  University of Bologna
!              David P. Schmidt    University of Massachusetts - Amherst
!              Salvatore Filippone University of Rome Tor Vergata
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without modification,
! are permitted provided that the following conditions are met:
!
!     1. Redistributions of source code must retain the above copyright notice,
!        this list of conditions and the following disclaimer.
!     2. Redistributions in binary form must reproduce the above copyright notice,
!        this list of conditions and the following disclaimer in the documentation
!        and/or other materials provided with the distribution.
!     3. Neither the name of the NEMO project nor the names of its contributors
!        may be used to endorse or promote products derived from this software
!        without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
! ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
! WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
! DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
! ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!
!---------------------------------------------------------------------------------
SUBMODULE (tools_mesh) cmp_mesh_implementation
    USE class_face
    IMPLICIT NONE

    CONTAINS

        MODULE PROCEDURE cmp_moving_surf
        USE class_psblas!, ONLY : mypnum_
        USE class_connectivity
        USE class_surface!, ONLY : alloc_surface
        USE class_vertex
        IMPLICIT NONE
        !! $Id: cmp_moving_surf.f90 8157 2014-10-09 13:02:44Z sfilippo $
        !!
        !! Description: Allocates and calculates shapes of moving surfaces
        !!
        INTEGER :: ib, info

        IF(mypnum_() == 0) THEN
            WRITE(*,*) 'Autodetecting virtual surfaces:'

            ! Not necessary since INTENT(OUT)
    !!$  if (associated (surf) ) then
    !!$     call abort_psblas
    !!$  endif

            ! Allocates pointer of SURFACE objects
            ALLOCATE(surf(nbc),stat=info)
            IF (info /= 0) THEN
                WRITE(*,100)
                CALL abort_psblas
            ENDIF

            ! Identifies surfaces
            DO ib = 1, nbc
                CALL alloc_surface(v2b,ib,verts,surf(ib))
            END DO

            WRITE(*,*)

        ENDIF

100     FORMAT(' ERROR! Failure to allocate memory in CMP_MOVING_SURF')

        END PROCEDURE cmp_moving_surf

        MODULE PROCEDURE cmp_mesh_c2c
        USE class_psblas
        USE class_connectivity
        IMPLICIT NONE
        !! $Id: cmp_mesh_c2c.f90 8157 2014-10-09 13:02:44Z sfilippo $
        !!
        !! Description:
        !!
        INTEGER :: ic, IF, im, is, info
        INTEGER :: n, ncells, nfl_faces
        INTEGER, ALLOCATABLE :: itab(:,:), kconn(:)

        IF(mypnum_() == 0) THEN
            WRITE(*,*) 'Computing cells adjacency graph C2C'

            n         = f2c%max_conn()            ! Maximum connectivity degree
            ncells    = f2c%nel_()                ! Number of cells
            nfl_faces = COUNT(faces%flag_() == 0) ! Number of fluid faces

            ALLOCATE(itab(n,ncells),kconn(ncells),stat=info)
            IF(info /= 0) THEN
                WRITE(*,100)
                CALL abort_psblas
            END IF

            itab = 0
            kconn = 0

            ! WARNING! CMP_MESH_C2C is called before the global-to-local reallocation.
            ! Thus faces are still numbered in such a way that fluid ones come first

            DO IF = 1, nfl_faces
            im = faces(IF)%master_()
            is = faces(IF)%slave_()
                kconn(im) = kconn(im) + 1
                kconn(is) = kconn(is) + 1
                itab(kconn(im),im) = is
                itab(kconn(is),is) = im
            END DO

            ! Counts total amount of connectivities
            n = 0
            DO ic = 1, ncells
                n = n + kconn(ic)
            END DO

            ! Creates C2C connectivity object
            CALL alloc_conn(c2c,nel=ncells,nconn=n)

            ! Sets C2C connectivity
            DO ic = 1, ncells
                n = kconn(ic)
                CALL c2c%set_ith_conn(ic,itab(1:n,ic))
            END DO

            DEALLOCATE(itab,kconn,stat=info)
            IF(info /= 0) THEN
                WRITE(*,200)
                CALL abort_psblas
            END IF

        ENDIF

100     FORMAT(' ERROR! Memory allocation failure in CMP_MESH_C2C')
200     FORMAT(' ERROR! Memory deallocation failure in CMP_MESH_C2C')

        END PROCEDURE cmp_mesh_c2c

        MODULE PROCEDURE cmp_mesh_f2b
        USE class_psblas
        USE class_connectivity
        IMPLICIT NONE
        !! $Id: cmp_mesh_f2b.f90 3175 2008-06-13 12:59:07Z sfilippo $
        !!
        !! Description:
        !!
        INTEGER :: ib, IF, info, maxconn, n, nfaces,j
        INTEGER :: flag(SIZE(faces)), iface(SIZE(faces)), kf(-1:nbc)
        INTEGER, ALLOCATABLE :: iconn(:)

        IF(mypnum_() == 0) &
            & WRITE(*,*) 'Computing face to boundary connectivity F2B'

        nfaces = SIZE(faces)
        iface = (/(IF, IF = 1, nfaces)/)

        flag = faces%flag_()

        ! Computes maximum connectivity degree
        IF (.FALSE.) THEN
            maxconn = 0
            DO ib = -1, nbc
                kf(ib) = COUNT(flag == ib)
                maxconn = MAX(maxconn,kf(ib))
            END DO
        ELSE
            CALL psb_msort(flag(1:nfaces),ix=iface(1:nfaces),flag=psb_sort_keep_idx_)
            IF (nfaces > 0) THEN
                IF ((flag(1) < -1).OR.(flag(nfaces)>nbc)) THEN
                    WRITE(0,*) 'Error in cmp_mesh_f2b: flags out of bounds ',flag(1),flag(nfaces)
                    CALL abort_psblas
                END IF
            END IF
            maxconn = 0
            kf(:)   = 0
            DO IF = 1, nfaces
                ib  = flag(IF)
                kf(ib) = kf(ib) + 1
                maxconn = MAX(maxconn,kf(ib))
            END DO
        END IF

        CALL alloc_conn(f2b,nel=nbc+2,nconn=nfaces,lb=-1)

        ALLOCATE(iconn(maxconn),stat=info)
        IF(info /= 0) THEN
            WRITE(*,100)
            CALL abort_psblas
        END IF

        j = 0
        DO ib = -1, nbc
            n = kf(ib)
            iconn(1:n) = iface(j+1:j+n)
    !!$     iconn(1:n) = pack(iface,flag == ib)
            CALL f2b%set_ith_conn(ib,iconn(1:n))
            j = j + n
        END DO

        DEALLOCATE(iconn)

100     FORMAT(' ERROR! Memory allocation failure in CMP_MESH_F2B')

        END PROCEDURE cmp_mesh_f2b

        MODULE PROCEDURE cmp_mesh_f2f
        USE class_psblas
        USE class_connectivity
        IMPLICIT NONE
        !! $Id: cmp_mesh_f2f.f90 8157 2014-10-09 13:02:44Z sfilippo $
        !!
        !! Description:
        !!    Provides the Face To Face adjacency graph.
        !!    The criterium for associating the graph nodes is arbitrarly defined by the
        !!    user. In this case, using a co-located cell-centered strategy for variable
        !!    displacement, the F2F connectivity is obtained by linking each face to all
        !!    the ones belonging to its master and slave cells.
        !!
        INTEGER :: i, ic, IF, if1, if2, info, j, n, ncells, ncf
        INTEGER, POINTER :: if2c(:) => NULL()
        INTEGER, ALLOCATABLE :: kconn(:), itab(:,:)

        IF(mypnum_() == 0) THEN
            WRITE(*,*) 'Computing faces adjacency graph F2F'

            n = (f2c%max_conn() - 1) * 2 ! F2F maximum connectivity degree
            ncells = f2c%nel_()          ! Number of cells

            ALLOCATE(itab(n,nfaces),kconn(nfaces),stat=info)
            IF(info /= 0) THEN
                WRITE(*,100)
                CALL abort_psblas
            END IF

            kconn = 0
            DO ic = 1, ncells
                CALL f2c%get_ith_conn(if2c,ic)
                ncf = SIZE(if2c)
                DO i = 1, ncf - 1
                    if1 = if2c(i)
                    DO j = i + 1, ncf
                        if2 = if2c(j)

                        ! Add direct arc
                        kconn(if1) = kconn(if1) + 1
                        itab(kconn(if1),if1) = if2

                        ! Add indirect arc
                        kconn(if2) = kconn(if2) + 1
                        itab(kconn(if2),if2) = if1
                    END DO
                END DO
            END DO

            ! Computes total number of connectivities
            n = 0
            DO IF = 1, nfaces
                n = n + kconn(IF)
            END DO
            CALL alloc_conn(f2f,nel=nfaces,nconn=n)

            ! Sets F2F connectivity
            DO IF = 1, nfaces
                n = kconn(IF)
                CALL f2f%set_ith_conn(IF,itab(1:n,IF))
            END DO

            ! Frees memory
            NULLIFY(if2c)
            DEALLOCATE(itab,kconn)

        ENDIF

100     FORMAT(' ERROR! Memory allocation failure in CMP_MESH_F2F')

        END PROCEDURE cmp_mesh_f2f

        MODULE PROCEDURE cmp_mesh_part
        USE class_psblas
        USE class_connectivity
        USE tools_part
        USE part_graph, ONLY : bld_part_graph, get_part_graph
        USE part_block, ONLY : bld_part_block
        USE part_random, ONLY : bld_part_random
        IMPLICIT NONE
        !! $Id: cmp_mesh_part.f90 8157 2014-10-09 13:02:44Z sfilippo $
        !!
        !! Description:
        !!
        INTEGER :: info, nprocs, ncells
        INTEGER, ALLOCATABLE :: xadj_glob(:)
        INTEGER, ALLOCATABLE :: adjncy_glob(:)

        IF (mypnum_() == 0) THEN
            CALL sw_par%tic()

            nprocs = nprocs_()

            IF(nprocs > 1 ) WRITE(*,*)

            ncells = c2c%nel_()

            ! Checks status of PART_CELLS vector
            IF(ALLOCATED(part_cells)) DEALLOCATE(part_cells)
            ALLOCATE(part_cells(ncells),stat=info)
            IF(info /= 0) THEN
                WRITE(*,100)
                CALL abort_psblas
            END IF

            ! Retrieves connectivity data in CSR-like format
            IF(ipart > 0) THEN
                ! ParMetis requires global connectivity in CSR format, therefore
                ! the call to GET_CONN_CSR is done here, outside BLD_PART_GRAPH.
                CALL c2c%get_conn_csr(xadj_glob,adjncy_glob)
            END IF

            SELECT CASE(ipart)
            CASE(-1)
                ! Random partitioning
                CALL bld_part_random(ncells,nprocs,part_cells)
            CASE(0)
                ! HPF Block partitioning
                CALL bld_part_block(ncells,nprocs,part_cells)
            CASE(1)
                ! METIS Partitioning
                CALL bld_part_graph(xadj_glob,adjncy_glob,ipart=1)
                CALL get_part_graph(part_cells)
            CASE DEFAULT
                WRITE(*,200)
                CALL abort_psblas
            END SELECT
            WRITE(*,*)

            IF(ipart > 0) THEN
                DEALLOCATE(xadj_glob,adjncy_glob,stat=info)
                IF(info /= 0) THEN
                    WRITE(*,300)
                    CALL abort_psblas
                END IF
            END IF

            CALL sw_par%toc()
        ENDIF

100     FORMAT(' ERROR! Memory allocation failure in CMP_MESH_PART')
200     FORMAT(' ERROR! Unknown partitioning method in CMP_MESH_PART')
300     FORMAT(' ERROR! Memory deallocation failure in CMP_MESH_PART')

        END PROCEDURE cmp_mesh_part

        MODULE PROCEDURE cmp_mesh_renum
        USE class_psblas
        USE class_cell
        USE class_connectivity
        USE renum
        IMPLICIT NONE
        !! $Id: cmp_mesh_renum.f90 8157 2014-10-09 13:02:44Z sfilippo $
        !!
        !! Description:
        !!
        INTEGER :: icontxt, mypnum

        icontxt = icontxt_()
        mypnum = mypnum_()

        CALL sw_ord%tic()

        SELECT CASE(irenum)
        CASE(0)
            IF(mypnum == 0) THEN
                WRITE(*,*) 'Matrix renumbering: none'
                WRITE(*,*)
            END IF
        CASE(1)
            ! Builds permutation array on P0, according to IRENUM value.
            IF(mypnum == 0) THEN
                CALL start_renum(irenum,c2c)
                !!     call print_renum(0)
                ! Applies permutation to cell-related structures

                CALL apply_renum(c2c,to_a_and_b_)
                CALL apply_renum(v2c,to_b_)
                CALL apply_renum(f2c,to_b_)
                CALL apply_renum(c2g,to_a_)
                CALL apply_renum(cells)
                CALL apply_renum(faces)

                ! Deallocates private storage of renum module
                CALL stop_renum
            ENDIF
        CASE DEFAULT
            WRITE(*,100)
            CALL abort_psblas
        END SELECT

        CALL sw_ord%toc()

100     FORMAT(' ERROR! Unknown renumbering method in CMP_MESH_RENUM')

        END PROCEDURE cmp_mesh_renum

        MODULE PROCEDURE cmp_mesh_v2b
        USE class_psblas
        USE class_connectivity
        IMPLICIT NONE
        !! $Id: cmp_mesh_v2b.f90 8157 2014-10-09 13:02:44Z sfilippo $
        !!
        !! Description:  Sets up vertex to boundary connectivity, so that given a boundary, the
        !! GET_ITH_CONN method can be used to obtain the list of vertices that belong to the boundary.
        !! Boundary ID's correspond to those used to flag faces:
        !!
        !!---------------------------|
        !! Vertex Flag | Vertex Type |
        !!---------------------------|
        !!      0      | interior    |
        !!    < 0      | proc. bndry |
        !!    > 0      | phys. bndry |
        !!---------------------------|
        !!
        LOGICAL :: inspected(0:nbc)
        INTEGER :: flag, IF, info, iv, k, n, nverts
        INTEGER :: kconn(0:nbc)
        INTEGER, ALLOCATABLE :: itab(:,:)
        INTEGER, POINTER :: if2v(:) => NULL()
        TYPE(connectivity) :: f2v

        IF(mypnum_() == 0) THEN
            WRITE(*,*) 'Computing vertex to boundary connectivity V2B'

            CALL v2f%get_dual_conn(f2v)

            nverts = f2v%nel_()

            ALLOCATE(itab(nverts,0:nbc),stat=info)
            IF(info /= 0) THEN
                WRITE(*,100)
                CALL abort_psblas
            END IF

            itab = 0
            kconn = 0

            DO iv = 1, nverts
                inspected(:) = .FALSE.

                CALL f2v%get_ith_conn(if2v,iv)
                n = SIZE(if2v)

                DO k = 1, n
                    IF = if2v(k)
                    flag = faces(IF)%flag_()
                    IF(flag > 0 .AND. .NOT.inspected(flag)) THEN
                        inspected(flag) = .TRUE.
                        kconn(flag) = kconn(flag) + 1
                        itab(kconn(flag),flag) = iv
                    END IF
                END DO

                IF(ALL(inspected .EQV. .FALSE.)) THEN
                    kconn(0) = kconn(0) + 1
                    itab(kconn(0),0) = iv
                END IF
            END DO

            ! Computes number of V2B connectivity
            n = 0
            DO k = 0, nbc
                n = n + kconn(k)
            END DO

            CALL alloc_conn(v2b,nel=nbc+1,nconn=n,lb=0)

            DO k = 0, nbc
                n = kconn(k)
                CALL v2b%set_ith_conn(k,itab(1:n,k))
            END DO

            DEALLOCATE(itab)
            NULLIFY(if2v)
            CALL free_conn(f2v)

        ENDIF

100     FORMAT(' ERROR! Memory allocation failure in CMP_MESH_V2B')

        END PROCEDURE cmp_mesh_v2b

        MODULE PROCEDURE  cmp_mesh_v2e
        USE class_psblas
        USE class_connectivity
        IMPLICIT NONE
        !! $Id: cmp_mesh_v2e.f90 8157 2014-10-09 13:02:44Z sfilippo $
        !!
        !! Description:
        !!    Provides the Vertex To Edge connectivity.
        !!
        LOGICAL, ALLOCATABLE :: inspected(:), twice(:)
        INTEGER :: ie, ie1, ie2
        INTEGER :: iv, iv1, iv2, iv3, iv4
        INTEGER :: i, IF, info, j, n, nedges, nfaces, nverts
        INTEGER, POINTER :: iv2f(:) => NULL()
        INTEGER, POINTER :: iv2e(:) => NULL(), ie2v(:) => NULL()
        TYPE(connectivity) :: work, e2v
        TYPE(stopwatch) :: sw_v2e

        INTEGER :: d1, d2, s1, s2

        ! V2E makes sense only for 3D meshes!
        IF(ncd == 2) RETURN

        IF(mypnum_() == 0) &
            WRITE(*,*) 'Computing vertex to edge connectivity V2E'

        ! Constructs stopwatch
        sw_v2e = stopwatch_(icontxt_())

        ! Start timing
        CALL sw_v2e%tic()

        IF(mypnum_() == 0) THEN

            nfaces = v2f%nel_()

            nedges = 0
            DO IF = 1, nfaces
                CALL v2f%get_ith_conn(iv2f,IF)
                nedges = nedges + SIZE(iv2f)
            END DO

            CALL alloc_conn(work,nel=nedges,nconn=2 * nedges)

            ie = 1
            DO IF = 1, nfaces

                CALL v2f%get_ith_conn(iv2f,IF)

                n = SIZE(iv2f)
                ! # of face vertices is equal to # of face edges N

                DO i = 1, n - 1
                    CALL work%set_ith_conn(ie,(/ iv2f(i), iv2f(i+1) /))
                    ie = ie + 1
                END DO
                CALL work%set_ith_conn(ie,(/ iv2f(n), iv2f(1) /))
                ie = ie + 1
            END DO

            CALL work%get_dual_conn(e2v)

            ALLOCATE(inspected(nedges),twice(nedges),stat=info)
            IF(info /= 0) THEN
                WRITE(*,100)
                CALL abort_psblas
            END IF

            nverts = e2v%nel_()

            inspected = .FALSE.
            twice     = .FALSE.
            DO iv = 1,  nverts
                CALL e2v%get_ith_conn(ie2v,iv)
                n = SIZE(ie2v)

                edge1: DO i = 1, n - 1
                    ie1 = ie2v(i)
                    IF(twice(ie1) .OR. inspected(ie1)) CYCLE edge1
                    inspected(ie1) = .TRUE.

                    CALL work%get_ith_conn(iv2e,ie1)
                    iv1 = iv2e(1)
                    iv2 = iv2e(2)
                    s1 = iv1 + iv2
                    d1 = ABS(iv1 - iv2)

                    edge2: DO j = i + 1, n
                        ie2 = ie2v(j)
                        IF(twice(ie2) .OR. inspected(ie2)) CYCLE edge2

                        CALL work%get_ith_conn(iv2e,ie2)
                        iv3 = iv2e(1)
                        iv4 = iv2e(2)
                        s2 = iv3 + iv4
                        d2 = ABS(iv3 - iv4)

                        IF(s1 == s2 .AND. d1 == d2) twice(ie2) = .TRUE.
    !!$           if((iv4 == iv1 .and. iv3 == iv2) .or.&
    !!$               iv4 == iv2 .and. iv3 == iv1 ) then
    !!$              twice(ie2) = .true.
    !!$           end if
                    END DO edge2

                END DO edge1
            END DO

            nedges = COUNT(.NOT.twice)
            CALL alloc_conn(v2e,nel=nedges,nconn=nedges*2)

            ie = 0
            n = work%nel_()
            DO i = 1, n
                IF(twice(i)) CYCLE
                ie = ie + 1
                CALL work%get_ith_conn(iv2e,i)
                CALL v2e%set_ith_conn(ie,iv2e)
            END DO

            ! Frees storage
            NULLIFY(iv2e,ie2v,iv2f)
            DEALLOCATE(inspected,twice)
            CALL free_conn(e2v)
            CALL free_conn(work)

        ENDIF

        ! Stop timing
        CALL sw_v2e%toc()

        ! Synchronization
        CALL sw_v2e%synchro()

        IF(mypnum_() == 0) WRITE(*,'(a,es13.6)') &
            &'  - Elapsed time for V2E building:', sw_v2e%partial_()

100     FORMAT(' ERROR! Memory allocation failure in CMP_MESH_V2E')

        END PROCEDURE cmp_mesh_v2e

        MODULE PROCEDURE cmp_mesh_v2v
        USE class_psblas
        USE class_connectivity
        IMPLICIT NONE
        !! $Id: cmp_mesh_v2v.f90 3093 2008-04-22 14:51:09Z sfilippo $
        !!
        !! Description:
        !!    Provides the Vertex To Vertex adjacency graph.
        !!    The criterium for associating the graph nodes is arbitrarly defined by the
        !!    user. In this case, using a co-located cell-centered strategy for variable
        !!    displacement, the V2V connectivity is obtained by linking each vertex to all
        !!    the ones belonging to the cells which share it.
        !!
        INTEGER, PARAMETER :: ncv_max = 50 ! Max # of cells sharing a vertex
        INTEGER :: i, ic, inb, info, iv, iv1, iv2, j, n, ncv, ncells
        INTEGER, POINTER :: iv2c(:) => NULL()
        INTEGER, ALLOCATABLE :: kconn(:), itab(:,:)
        LOGICAL :: inspected(nverts)
        TYPE(stopwatch) :: sw_v2v

        IF(mypnum_() == 0) WRITE(*,*) 'Computing vertices adjacency graph V2V'

        ! Constructs stopwatch
        sw_v2v = stopwatch_(icontxt_())

        ! Start timing
        CALL sw_v2v%tic()

        n = (v2c%max_conn() - 1) * ncv_max ! Maximum conn. degree (overestimated)
        ncells = v2c%nel_()

        ALLOCATE(itab(n,nverts),kconn(nverts),stat=info)
        IF(info /= 0) THEN
            WRITE(*,100)
            CALL abort_psblas
        END IF

        kconn = 0
        DO ic = 1, ncells
            CALL v2c%get_ith_conn(iv2c,ic)
            ncv = SIZE(iv2c)
            DO i = 1, ncv - 1
                iv1 = iv2c(i)
                DO j = i + 1, ncv
                    iv2 = iv2c(j)

                    ! Add direct arc
                    kconn(iv1) = kconn(iv1) + 1
                    itab(kconn(iv1),iv1) = iv2

                    ! Add indirect arc
                    kconn(iv2) = kconn(iv2) + 1
                    itab(kconn(iv2),iv2) = iv1
                END DO
            END DO
        END DO

        ! Eliminates duplicates
        inspected = .FALSE.
        DO iv = 1, nverts
            j = 0

            n = kconn(iv)
            DO i = 1, n
                inb = itab(i,iv)
                IF(.NOT.inspected(inb)) THEN
                    j = j + 1
                    itab(j,iv) = inb
                    inspected(inb) = .TRUE.
                END IF
            END DO
            kconn(iv) = j

            ! Reset inspected
            DO i = 1, j
                inb = itab(i,iv)
                inspected(inb) = .FALSE.
            END DO
        END DO

        ! Computes total number of connectivities
        n = 0
        DO iv = 1, nverts
            n = n + kconn(iv)
        END DO
        CALL alloc_conn(v2v,nel=nverts,nconn=n)

        ! Sets V2V connectivity
        DO iv = 1, nverts
            n = kconn(iv)
            CALL v2v%set_ith_conn(iv,itab(1:n,iv))
        END DO

        ! Frees memory
        NULLIFY(iv2c)
        DEALLOCATE(kconn,itab)

        ! Stop timing
        CALL sw_v2v%toc()

        ! Synchronization
        CALL sw_v2v%synchro()

        IF(mypnum_() == 0) WRITE(*,'(a,es13.6)') &
            & '  - Elapsed time for V2V building:', sw_v2v%partial_()

100     FORMAT(' ERROR! Memory allocation failure in CMP_MESH_V2V')

        END PROCEDURE cmp_mesh_v2v

        MODULE PROCEDURE cmp_mesh_v2ve
        USE class_psblas
        USE class_connectivity
        USE tools_mesh, ONLY : cmp_mesh_v2e
        IMPLICIT NONE
        !! $Id: cmp_mesh_v2ve.f90 8157 2014-10-09 13:02:44Z sfilippo $
        !!
        !! Description:
        !!    Builds a vertex to vertex connectivity based on edge connections.
        !!
        INTEGER :: i, ie, info, iv, iv_nb, nconn, nmax, nverts
        INTEGER, POINTER :: ie2v(:) => NULL(), iv2e(:) => NULL()
        INTEGER, ALLOCATABLE :: iv2v(:,:), kconn(:)
        TYPE(connectivity) :: v2e, e2v
        TYPE(stopwatch) :: sw_v2v

        IF(mypnum_() == 0) WRITE(*,*) 'Computing vertices adjacency graph V2VE'

        ! Constructs stopwatch
        sw_v2v = stopwatch_(icontxt_())

        ! Start timing
        CALL sw_v2v%tic()

        SELECT CASE(ncd)
        CASE(2)
            IF(mypnum_() == 0)    CALL v2f%get_dual_conn(e2v)

        CASE(3)
            CALL cmp_mesh_v2e(ncd,v2f,v2e)
            IF(mypnum_() == 0)    CALL v2e%get_dual_conn(e2v)
        END SELECT

        IF(mypnum_() == 0) THEN

            nverts = e2v%nel_()
            nmax   = e2v%max_conn()

            ALLOCATE(iv2v(nverts,nmax),kconn(nverts),stat=info)
            IF(info /= 0) THEN
                WRITE(*,100)
                CALL abort_psblas
            END IF

            SELECT CASE(ncd)
            CASE(2) ! 2D

                DO iv = 1, nverts

                    CALL e2v%get_ith_conn(ie2v,iv)
                    nconn = SIZE(ie2v)

                    DO i = 1, nconn
                        ie = ie2v(i)

                        ! In 2D edges = faces => use V2F
                        CALL v2f%get_ith_conn(iv2e,ie)
                        iv_nb = iv2e(1) + iv2e(2) - iv

                        iv2v(iv,i) = iv_nb
                    END DO

                    kconn(iv) = nconn
                END DO

            CASE(3) ! 3D

                DO iv = 1, nverts

                    CALL e2v%get_ith_conn(ie2v,iv)
                    nconn = SIZE(ie2v)

                    DO i = 1, nconn
                        ie = ie2v(i)

                        ! In 3D edges /= faces => use V2E
                        CALL v2e%get_ith_conn(iv2e,ie)
                        iv_nb = iv2e(1) + iv2e(2) - iv

                        iv2v(iv,i) = iv_nb
                    END DO

                    kconn(iv) = nconn
                END DO

            END SELECT

            ! Computes total number of connectivities
            nconn = 0
            DO iv = 1, nverts
                nconn = nconn + kconn(iv)
            END DO

            ! Allocate and sets V2V
            CALL alloc_conn(v2v,nel=nverts,nconn=nconn)
            DO iv = 1, nverts
                nconn = kconn(iv)
                CALL v2v%set_ith_conn(iv,iv2v(iv,1:nconn))
            END DO

            ! Free memory storage
            NULLIFY(iv2e,ie2v)
            DEALLOCATE(iv2v,kconn)
            CALL free_conn(e2v)
            IF(ncd == 3) CALL free_conn(v2e)

        ENDIF

        ! Stop timing
        CALL sw_v2v%toc()

        ! Synchronization
        CALL sw_v2v%synchro()

        IF(mypnum_() == 0) WRITE(*,'(a,es13.6)') &
            &'  - Elapsed time for V2VE building:', sw_v2v%partial_()

100     FORMAT(' ERROR! Memory allocation failure in CMP_MESH_V2VE')

        END PROCEDURE cmp_mesh_v2ve

END SUBMODULE cmp_mesh_implementation

max_stars_repo_name: sourceryinstitute/MORFEUS-Source

max_issues_repo_name: sourceryinstitute/MORFEUS-Source

max_forks_repo_name: sourceryinstitute/MORFEUS-Source


text:
!> When a fine grid is updated, this subroutine is called to save fluxes
!! going into an adjacent coarse cell in **svdflx**, which is for conservative
!! fix later.
!!
c -------------------------------------------------------
c
      subroutine fluxad(xfluxm,xfluxp,yfluxm,yfluxp,
     1                  svdflx,mptr,mitot,mjtot,
     2                   nvar,lenbc,lratiox,lratioy,ng,dtf,dx,dy)
c

      use amr_module
      implicit double precision (a-h,o-z)

c :::::::::::::::::::: FLUXAD ::::::::::::::::::::::::::::::::::
c  save fine grid fluxes  at the border of the grid, for fixing
c  up the adjacent coarse cells. at each edge of the grid, only
c  save the plus or minus fluxes, as necessary. For ex., on
c  left edge of fine grid, it is the minus xfluxes that modify the
c  coarse cell.
c :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

      dimension xfluxm(nvar,mitot,mjtot), yfluxm(nvar,mitot,mjtot)
      dimension xfluxp(nvar,mitot,mjtot), yfluxp(nvar,mitot,mjtot)
      dimension svdflx(nvar,lenbc)

      nx  = mitot-2*ng
      ny  = mjtot-2*ng
      nyc = ny/lratioy
      nxc = nx/lratiox

c ::::: left side saved first
      lind = 0

      do 100 j=1,nyc
         lind = lind + 1
         jfine = (j-1)*lratioy + ng
         do 110 ivar = 1, nvar
            do 120 l=1,lratioy
               svdflx(ivar,lind) = svdflx(ivar,lind) +
     1                             xfluxm(ivar,ng+1,jfine+l)*dtf*dy
c              write(dbugunit,900)lind,xfluxm(ivar,1,jfine+l),
c     .                           xfluxp(ivar,1,jfine+l)
 900           format(' lind ', i4,' m & p ',2e15.7,' svd ',e15.7)
120         continue
110      continue
100   continue

c ::::: top side
c      write(dbugunit,*)" saving top side "
      do 200 i=1,nxc
         lind = lind + 1
         ifine = (i-1)*lratiox + ng
         do 210 ivar = 1, nvar
            do 220 l=1,lratiox
               svdflx(ivar,lind) = svdflx(ivar,lind) +
     1                     yfluxp(ivar,ifine+l,mjtot-ng+1)*dtf*dx
c              write(dbugunit,900)lind,yfluxm(ivar,ifine+l,mjtot-ng+1,
c     .                           ),yfluxp(ivar,ifine+l,mjtot-ng+1),
c     .                           svdflx(ivar,lind)
220         continue
210      continue
200   continue

c ::::: right side
      do 300 j=1,nyc
         lind = lind + 1
         jfine = (j-1)*lratioy + ng
         do 310 ivar = 1, nvar
            do 320 l=1,lratioy
               svdflx(ivar,lind) = svdflx(ivar,lind) +
     1                     xfluxp(ivar,mitot-ng+1,jfine+l)*dtf*dy
c              write(dbugunit,900)lind,xfluxm(ivar,mitot-ng+1,jfine+l,
c                                 ),xfluxp(ivar,mitot-ng+1,jfine+l)
320         continue
310      continue
300   continue

c ::::: bottom side
c      write(dbugunit,*)" saving bottom side "
      do 400 i=1,nxc
         lind = lind + 1
         ifine = (i-1)*lratiox + ng
         do 410 ivar = 1, nvar
            do 420 l=1,lratiox
               svdflx(ivar,lind) = svdflx(ivar,lind) +
     1                             yfluxm(ivar,ifine+l,ng+1)*dtf*dx
c              write(dbugunit,900)lind,yfluxm(ivar,ifine+l,ng+1),
c     .                      yfluxp(ivar,ifine+l,ng+1),svdflx(ivar,lind)
420         continue
410      continue
400   continue

      return
      end

ext: f

max_stars_repo_name: navravi/amrclaw


program main
    implicit none
    real :: l
    real, external :: period
    open(unit=1, file='out.txt')
    write (*, "(A)") 'Please input the length of pendulum:'
    read  (*, *) l
    write (1, "(A, 2(F0.6, A))") 'The period of the pendulum with length ', l, 'm is ', period(l), 's.'
    close(1)
    stop
end program main

function period(l)
    implicit none
    real :: l
    real :: period
    real, parameter :: pi = 3.1415926535
    real, parameter :: g = 9.80665
    period = 2 * pi * (l / g) ** 0.5
    return
end function


text:
      SUBROUTINE SHADEK(DIZ,PERSUN,PERSND,INDS,K,ISATD)
C
C           'SHADEK' COMPUTES SHADOWING OF ELEMENTS AND SWITCHES INDICES
C
      IMPLICIT REAL*8(A-H,O-Z)
C
      COMMON/IPOOL1/ IGRAV,IDAMP,IK,K1,ITIM,IAB,IAPS,IBB,IBPS,NK(10),
     .               LK(10),LLK(10)
C
      COMMON/RPOOL1/ RHOK(10),TIME,SA(3,3),FM1(3,3),ZLK(10),OMEG(3),
     .               ZLKP(10),ZLKDP(10),CMAT(3,3),GBAR(3,3),YBCM(3),
     .               ZBZK(3,10),FCM(3,3),DTO,PHID,PHI
C
      COMMON/RPOOL6/ FM(3,3),CIY(3,3),CIZ(3,3),SAT(3,3),SZ1,SZ2,SZ3
C
      COMMON/SATLSH/ TAUK(10),OCULTK(10),RADSH
C
      DIMENSION DIZ(3),OMK(3),ISATD(10)
C
       INDS=0
      TEST=DIZ(1)
      IF(TEST.LT.0.0D0) GO TO 6
      DTRAN=DSQRT(DIZ(2)**2+DIZ(3)**2)
      IF(DTRAN.LE.0.001D0) GO TO 5
      SHADL=RADSH*DIZ(1)/DTRAN
      PERSUN=1.0D0-SHADL/ZLK(K)
      IF(PERSUN.GT.0.0D0) GO TO 10
    5 CONTINUE
      PERSUN=0.0D0
      GO TO 10
    6 CONTINUE
      PERSUN=1.0D0
   10 CONTINUE
      IF(PERSUN.EQ.1.0D0.OR.PERSUN.EQ.0.0D0) GO TO 30
C
      DO 11 I=1,3
      OMK(I)=FCM(1,I)*OMEG(1)+FCM(2,I)*OMEG(2)+FCM(3,I)*OMEG(3)
      IF(IDAMP.EQ.0.OR.K.GT.K1) GO TO 11
      OMK(I)=OMK(I)+FM(2,I)*PHID
   11 CONTINUE
      DIKDT=OMK(3)*DIZ(2)-OMK(2)*DIZ(3)
C
      INDS=ISATD (K)
      SHADLD=RADSH*DIKDT/DTRAN**3
      PERSND=-SHADLD/ZLK(K)
      MIND=0
      TEST=PERSND*TAUK(K)
      TEST=DABS(TEST)
      IF(TEST.GT.OCULTK(K)) MIND=1
      IF(PERSND.GT.0.0D0) GO TO 20
C     ENTERING SHADOW
      IF(ISATD(K).NE.1.AND.MIND.EQ.1) INDS=1
      GO TO 40
   20 CONTINUE
      IF(ISATD(K).NE.2.AND.MIND.EQ.1) INDS=2
      GO TO 40
   30 CONTINUE
      PERSND=0.0D0
   40 CONTINUE
      RETURN
      END

ext: for


text:
module gemread
  implicit none

contains

  subroutine ggi(gemFilename, mxgrd, n, nx, ny, proj, ang, lllat, lllon, urlat, urlon,&
       gdattm, level, ivcord, vcord, parm)
    character(*),  intent(in)  :: gemFilename
    integer,       intent(out) :: mxgrd, n, nx, ny
    character(20), intent(out) :: proj
    real,          intent(out) :: ang(3)
    real,          intent(out) :: lllat, lllon, urlat, urlon
    integer(1), intent(inout) :: gdattm(:,:,:)  ! First dim 20 chars
    integer,       intent(inout) :: level(:,:)
    integer,       intent(inout) :: ivcord(:)
    integer(1),  intent(inout) :: vcord(:,:)  ! First dim 4 chars
    integer(1), intent(inout) :: parm(:,:)  ! First dim 12 chars

    integer, parameter :: LLNNAV = 256, LLNANL = 128

    real :: bknav(LLNNAV), bkanl(LLNANL)
    integer :: status, iacss, navsz, ianlsz, ihdrsz, projtype, i
    character(20) :: firstm, lasttm

    character(20), allocatable :: gdattmloc(:,:)
    character(4), allocatable :: vcordloc(:)
    character(12), allocatable :: parmloc(:)

    integer :: c1, c2, j, k
    c1 = size(gdattm,2)
    c2 = size(gdattm,3)

    allocate(gdattmloc(c1,c2), vcordloc(c1), parmloc(c1))

    ! Initialize to bad values in case of error
    mxgrd = -1
    n = -1
    nx = -1
    ny = -1
    proj = "ERROR"
    ang = -9999.
    lllat = -9999.
    lllon = -9999.
    urlat = -9999.
    urlon = -9999.

    call in_bdta(status)
    if (status == 0) then
       call gd_init(status)
    else
       print *, 'Inexplicable error with in_bdta!'
       return
    end if

    if (status /= 0) then
       print *, 'Error initializing GEMPAK interface!'
       return
    end if

    call gd_ofil(gemFilename, .false., .false., iacss, navsz, bknav, ianlsz,  &
         bkanl, ihdrsz, mxgrd, status)
    if (status /= 0) then
       print *, 'Error opening GEMPAK file!'
       return
    end if

    call gd_ngrd(iacss, n, firstm, lasttm, status)
    if (status /= 0) then
       print *, 'Error getting number of grids!'
       return
    end if

    do i = 1, n
       call gd_gidn(iacss, i, gdattmloc(i,:), level(i,:), ivcord(i), parmloc(i), &
            status)
       if (status /= 0) then
          print *, 'Error reading grid info!'
          exit
       end if
       call lv_ccrd(ivcord(i), vcordloc(i), status)
       if (status /= 0) print *, 'Oopsie!'
    end do

    call gd_clos(iacss, status)
    if (status /= 0) print *, 'Error closing GEMPAK file!'

    call st_itoc(bknav(2), 1, proj, status)
    if (status /= 0) print *, 'Could not determine projection!'

    ! Grab the data we need
    nx = bknav(5)
    ny = bknav(6)
    projtype = bknav(1)
    select case (projtype)
    case (1)
       print *, 'Case 1'
       ang = 0.
    case (2)
       ang = [bknav(11), bknav(12), bknav(13)]
    case (3)
       print *, 'Case 3'
       ang = 0.
    case default
       print *, 'Unimplemented projection'
       return
    end select

    lllat = bknav(7)
    lllon = bknav(8)
    urlat = bknav(9)
    urlon = bknav(10)

    !print *, gdattmloc(:100,:), vcordloc(:100), parmloc(:100)
    do j = 1, c2
       do i = 1, c1
          do k = 1, 20
             gdattm(k,i,j) = iachar(gdattmloc(i,j)(k:k))
          end do
       end do
    end do
    do i = 1, c1
       do j = 1, 4
          vcord(j,i) = iachar(vcordloc(i)(j:j))
       end do
       do j = 1, 12
          parm(j,i) = iachar(parmloc(i)(j:j))
       end do
    end do
  end subroutine ggi

  subroutine get_grid_info(gemFilename, mxgrd, n, nx, ny, proj, ang, lllat,  &
       lllon, urlat, urlon, gdattm, level, ivcord, vcord, parm)
    character(*),  intent(in)    :: gemFilename
    integer,       intent(out)   :: mxgrd, n, nx, ny
    character(20), intent(out)   :: proj
    real,          intent(out)   :: ang(3)
    real,          intent(out)   :: lllat, lllon, urlat, urlon
    character(*), intent(inout) :: gdattm(:,:)
    integer,       intent(inout) :: level(:,:)
    integer,       intent(inout) :: ivcord(:)
    character(*),  intent(inout) :: vcord(:)
    character(*), intent(inout) :: parm(:)

    integer, parameter :: LLNNAV = 256, LLNANL = 128

    real :: bknav(LLNNAV), bkanl(LLNANL)
    integer :: status, iacss, navsz, ianlsz, ihdrsz, projtype, i
    character(20) :: firstm, lasttm

    ! Initialize to bad values in case of error
    mxgrd = -1
    n = -1
    nx = -1
    ny = -1
    proj = "ERROR"
    ang = -9999.
    lllat = -9999.
    lllon = -9999.
    urlat = -9999.
    urlon = -9999.

    call in_bdta(status)
    if (status == 0) then
       call gd_init(status)
    else
       print *, 'Inexplicable error with in_bdta!'
       return
    end if

    if (status /= 0) then
       print *, 'Error initializing GEMPAK interface!'
       return
    end if

    call gd_ofil(gemFilename, .false., .false., iacss, navsz, bknav, ianlsz,  &
         bkanl, ihdrsz, mxgrd, status)
    if (status /= 0) then
       print *, 'Error opening GEMPAK file!'
       return
    end if

    call gd_ngrd(iacss, n, firstm, lasttm, status)
    if (status /= 0) then
       print *, 'Error getting number of grids!'
       return
    end if

    do i = 1, n
       call gd_gidn(iacss, i, gdattm(i,:), level(i,:), ivcord(i), parm(i), &
            status)
       if (status /= 0) then
          print *, 'Error reading grid info!'
          exit
       end if
       call lv_ccrd(ivcord(i), vcord(i), status)
       if (status /= 0) print *, 'Oopsie!'
    end do

    call gd_clos(iacss, status)
    if (status /= 0) print *, 'Error closing GEMPAK file!'

    call st_itoc(bknav(2), 1, proj, status)
    if (status /= 0) print *, 'Could not determine projection!'

    ! Grab the data we need
    nx = bknav(5)
    ny = bknav(6)
    projtype = bknav(1)
    select case (projtype)
    case (1)
       print *, 'Case 1'
       ang = 0.
    case (2)
       ang = [bknav(11), bknav(12), bknav(13)]
    case (3)
       print *, 'Case 3'
       ang = 0.
    case default
       print *, 'Unimplemented projection'
       return
    end select

    lllat = bknav(7)
    lllon = bknav(8)
    urlat = bknav(9)
    urlon = bknav(10)
  end subroutine get_grid_info

  subroutine get_max_grids(gemFilename, mxgrd)
    character(*), intent(in) :: gemFilename
    integer, intent(out) :: mxgrd

    integer, parameter :: LLNNAV = 256, LLNANL = 128

    real :: bknav(LLNNAV), bkanl(LLNANL)
    integer :: status, iacss, navsz, ianlsz, ihdrsz

    call in_bdta(status)
    if (status == 0) call gd_init(status)
    if (status /= 0) then
       print *, 'Error initializing GEMPAK interface!'
       mxgrd = -1
       return
    end if

    call gd_ofil(gemFilename, .false., .false., iacss, navsz, bknav, ianlsz,  &
         bkanl, ihdrsz, mxgrd, status)
    if (status /= 0) then
       print *, 'Error opening GEMPAK file!'
       mxgrd = -1
       return
    end if

    call gd_clos(iacss, status)
    if (status /=0) print *, 'Error closing GEMPAK file!'
  end subroutine get_max_grids

  subroutine get_num_grids(gemFilename, n)
    character(*), intent(in) :: gemFilename
    integer, intent(out) :: n

    integer, parameter :: LLNNAV = 256, LLNANL = 128

    real :: bknav(LLNNAV), bkanl(LLNANL)
    integer :: status, iacss, navsz, ianlsz, ihdrsz, mxgrd
    character(20) :: firstm, lasttm

    call in_bdta(status)
    if (status == 0) call gd_init(status)
    if (status /= 0) then
       print *, 'Error initializing GEMPAK interface!'
       n = -1
       return
    end if

    call gd_ofil(gemFilename, .false., .false., iacss, navsz, bknav, ianlsz,  &
         bkanl, ihdrsz, mxgrd, status)
    if (status /= 0) then
       print *, 'Error opening GEMPAK file!'
       n = -1
       return
    end if

    call gd_ngrd(iacss, n, firstm, lasttm, status)
    if (status /= 0) then
       print *, 'Error getting number of grids!'
       n = -1
       return
    end if

    call gd_clos(iacss, status)
    if (status /=0) print *, 'Error closing GEMPAK file!'
  end subroutine get_num_grids

  subroutine get_grid_dims(gemFilename, kx, ky)
    character(*), intent(in) :: gemFilename
    integer, intent(out) :: kx, ky

    integer, parameter :: LLNNAV = 256, LLNANL = 128

    real :: bknav(LLNNAV), bkanl(LLNANL)
    integer :: status, iacss, navsz, ianlsz, ihdrsz, mxgrd
    character(4) :: proj

    kx = -1 ; ky = -1

    call in_bdta(status)
    if (status == 0) call gd_init(status)
    if (status /= 0) then
       print *, 'Error initializing GEMPAK interface!'
       return
    end if

    call gd_ofil(gemFilename, .false., .false., iacss, navsz, bknav, ianlsz,  &
         bkanl, ihdrsz, mxgrd, status)
    if (status /= 0) then
       print *, 'Error opening GEMPAK file!'
       return
    end if

    call gr_rnav(bknav, proj, kx, ky, status)
    if (status /= 0) print *, 'Error reading GEMPAK navigation block!'

    call gd_clos(iacss, status)
    if (status /=0) print *, 'Error closing GEMPAK file!'
  end subroutine get_grid_dims

  subroutine read_grid(gemFilename, gdattm1, gdattm2, level1, level2, vcord, &
       parmin, grid)
    character(*), intent(in) :: gemFilename, gdattm1, gdattm2, vcord, parmin
    integer, intent(in) :: level1, level2
    real, dimension(:,:), intent(inout) :: grid

    integer, parameter :: LLNNAV = 256, LLNANL = 128

    integer, dimension(:), allocatable :: ighdr

    real :: bknav(LLNNAV), bkanl(LLNANL)
    integer, dimension(2) :: level
    integer :: ivcord, status, iacss, navsz, ianlsz, ihdrsz, mxgrd, igx, igy
    character(20) :: gdattm(2)
    character(12) :: parm
    character(4) :: vcoord, junk

    grid = -999.

    gdattm = (/ gdattm1, gdattm2 /)

    level = (/ level1, level2 /)

    call st_lcuc(vcord, vcoord, status)
    call lv_cord(vcoord, junk, ivcord, status)
    if (status /= 0) then
       print *, 'Error with vertical coordinate!'
       return
    end if

    parm = parmin

    call in_bdta(status)
    if (status == 0) call gd_init(status)
    if (status /= 0) then
       print *, 'Error initializing GEMPAK interface!'
       return
    end if

    call gd_ofil(gemFilename, .false., .false., iacss, navsz, bknav, ianlsz,  &
         bkanl, ihdrsz, mxgrd, status)
    if (status /= 0) then
       print *, 'Error opening GEMPAK file!'
       return
    end if

    allocate(ighdr(ihdrsz))

    call gd_rdat(iacss, gdattm, level, ivcord, parm, grid, igx, igy, ighdr,  &
         status)
    if (status /=0) print *, 'Error reading GEMPAK file!'

    call gd_clos(iacss, status)
    if (status /=0) print *, 'Error closing GEMPAK file!'
  end subroutine read_grid

  subroutine create_gemfile(newfname, copyfname, mxgrd, status)
    character(*), intent(in) :: newfname, copyfname
    integer, intent(in) :: mxgrd
    integer, intent(out) :: status

    integer, parameter :: LLNNAV = 256, LLNANL = 128

    real :: bknav(LLNNAV), bkanl(LLNANL)
    integer :: iacss, navsz, ianlsz, ihdrsz, origmxgrd
    logical :: exists

    status = -1

    inquire(file=newfname, exist=exists)
    if (exists) then
       print *, 'File to create already exists!'
       return
    end if
    inquire(file=copyfname, exist=exists)
    if (.not. exists) then
       print *, 'File to copy from does not exist!'
       return
    end if

    call in_bdta(status)
    if (status == 0) call gd_init(status)
    if (status /= 0) then
       print *, 'Error initializing GEMPAK interface!'
       return
    end if

    call gd_ofil(copyfname, .false., .false., iacss, navsz, bknav, ianlsz,  &
         bkanl, ihdrsz, origmxgrd, status)
    if (status /= 0) then
       print *, 'Error opening GEMPAK file!'
       return
    end if

    call gd_clos(iacss, status)
    if (status /= 0) print *, 'Error closing GEMPAK file!'

    call gd_cref(newfname, LLNNAV, bknav, LLNANL, bkanl, 2, mxgrd, iacss,  &
         status)
    if (status /= 0) then
       print *, 'Error creating GEMPAK file!'
       return
    end if

    call gd_clos(iacss, status)
    if (status /= 0) print *, 'Error closing GEMPAK file!'
  end subroutine create_gemfile

  subroutine write_grid(gemFilename, gdattm1, gdattm2, level1, level2, vcord, &
       parmin, grid, status)
    character(*), intent(in) :: gemFilename, gdattm1, gdattm2, vcord, parmin
    integer, intent(in) :: level1, level2
    real, dimension(:,:), intent(in) :: grid
    integer, intent(out) :: status

    integer, parameter :: LLNNAV = 256, LLNANL = 128

    integer, dimension(2) :: ighdr = 0

    real :: bknav(LLNNAV), bkanl(LLNANL)
    integer, dimension(2) :: level
    integer :: ivcord, iacss, navsz, ianlsz, ihdrsz, mxgrd, igx, igy
    character(20) :: gdattm(2)
    character(12) :: parm
    character(4) :: vcoord, junk

    igx = size(grid, 1)
    igy = size(grid, 2)

    gdattm = (/ gdattm1, gdattm2 /)

    level = (/ level1, level2 /)

    call st_lcuc(vcord, vcoord, status)
    call lv_cord(vcoord, junk, ivcord, status)
    if (status /= 0) then
       print *, 'Error with vertical coordinate!'
       return
    end if

    parm = parmin

    call in_bdta(status)
    if (status == 0) call gd_init(status)
    if (status /= 0) then
       print *, 'Error initializing GEMPAK interface!'
       return
    end if

    call gd_ofil(gemFilename, .true., .false., iacss, navsz, bknav, ianlsz,  &
         bkanl, ihdrsz, mxgrd, status)
    if (status /= 0) then
       print *, 'Error opening GEMPAK file!'
       return
    end if

    call gd_wdat(iacss, grid, igx, igy, ighdr, gdattm, level, ivcord, parm,  &
         .false., status)
    if (status /= 0) then
       print *, 'Error writing GEMPAK file!'
       print *, status
    end if

    call gd_clos(iacss, status)
    if (status /= 0) print *, 'Error closing GEMPAK file!'
  end subroutine write_grid
end module gemread

program test
  use gemread
  implicit none

  character(*), parameter :: fname = 'nam211.gem', fname2 = 'newfile.gem'
  integer,      parameter :: lun = 19

  real,          dimension(:,:), allocatable :: grid
  integer,       dimension(:,:), allocatable :: level
  character(20), dimension(:,:), allocatable :: gdattm
  integer,       dimension(:),   allocatable :: ivcord
  character(4),  dimension(:),   allocatable :: vcord
  character(12), dimension(:),   allocatable :: parm
  real,          dimension(3) :: ang
  real          :: lllat, lllon, urlat, urlon
  integer       :: i, n, nx, ny, stat
  character(20) :: proj

  call get_num_grids(fname, n)
  allocate(gdattm(n,2), level(n,2), ivcord(n), vcord(n), parm(n))
  call get_grid_info(fname, i, n, nx, ny, proj, ang, lllat, lllon, urlat, &
       urlon, gdattm, level, ivcord, vcord, parm)
!  call get_max_grids(fname, i)
  print *, i
!  call get_num_grids(fname, n)
  print *, n
!  call get_grid_dims(fname, nx, ny)
  print *, nx, ny
  print *, proj, ang, lllat, lllon, urlat, urlon
  print *, gdattm(:100,:)
  print *, level(:100,:)
  print *, ivcord(:100)
  print *, vcord(:100)
  print *, parm(:100)
  allocate(grid(nx,ny))
  call read_grid(fname, '150610/1200F000', '', 850, -1, 'PRES', 'HGHT', grid)
  print *, grid(10,20)
  call create_gemfile(fname2, fname, 1000, i)
  print *, i
  call write_grid(fname2, '121025/0000F000', '', 850, -1, 'PRES', 'HGHT',  &
       grid, i)
  print *, i
  open(unit=19, iostat=stat, file=fname2, status='old')
  if (stat == 0) close(19, status='delete')
end program test

max_stars_repo_name: sgdecker/pygempak

max_issues_repo_name: sgdecker/pygempak

max_forks_repo_name: sgdecker/pygempak


text:
! Program to test initialization expressions involving subobjects
program parameter_2
  implicit none
  type :: SS
     integer :: I
     integer :: J
  end type SS
  type :: TT
    integer :: N
    type (SS), dimension(2) :: o
  end type

  type (SS), parameter :: s =  SS (1, 2)
  type (TT), parameter :: t = TT(42, (/ SS(3, 4), SS(8, 9) /))

  integer, parameter :: a(2) = (/5, 10/)
  integer, parameter :: d1 = s%i
  integer, parameter :: d2 = a(2)
  integer, parameter :: d4 = t%o(2)%j

  integer q1, q2, q3, q4
  common /c1/q1(d1), q2(d2), q3(a(1)), q4(d4) ! legal
end

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
A cheesecake is a dessert consisting of soft cheese normally on top of crushed cookies, a pastry, or a sponge. Cheesecake styles can vary quite a bit depending on what region the style comes from. If you ever get in the mood for cheesecake in Davis, you should check out the following places:

  Our House
  Burgers and Brew  They have a cheesecake with Nutella spread on it.
  Ciocolat
  Jack In The Box
  Safeway  They normally have a New York style cheesecake but can be quite dry. Good place to go if you are in a hurry or in the mood for something cheap.

Outside of Davis
  The Cheesecake Factory
  Kitchen 428/ Mojos Lounge and Bar  Offers a rotating cheesecake selection made by pastry chef Danny Gutierez.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
What better way to express your love for Davis Wikiand spread Wikivangelists The Word than putting a wiki sticker on your laptop, bike, car, binder or child? The original sticker reserve has since run dry, so rather than make a second run of the same sticker, Users/ArlenAbraham I thought Id give the community a chance to design the next wiki sticker.

 What

Were having a sticker design contest. The Wikizen who designs the best sticker will have their design printed en masse and receive a pile of free stickers. Please limit your stickers to smallish squares or rectangles. Two colors is preferable (one vinyl color, one ink color), but we might make exceptions if the design is super awesome. Designs should be in a vector format or bitmapped at 300 dpi. Leave margins of at least 1/8. Dont try to do any halftoning or any other pixel tricks as the sticker process doesnt take kindly to that kind of thing.

 Who

You design it, users/philipneustrom Philip picks a winner (hes paying for them).

 When

Please upload your designs to this page by midnight on May 21.

 How

Here is some source material:

If you dont have Photoshop, Flash or Illustrator, you can use opensource alternatives like http://www.gimp.org/ The GIMP or http://www.inkscape.org/ Inkscape.

 Late Entries

Image(URL.png, thumbnail, 400, Per Dauberts suggestion) {{{<}}} You cant tell from this, but everything that is white would be transparent, and the sticker itself would be those resuable static cling to glass type stickers, Afixable in a store window or a car etc ...  Users/StevenDaubert
Sadly, looking into it, static cling stickers are a bit pricey, but then again, I only checked google. ~DavePoole

Something in the same vein as the original. The tiffs are the correct size (4.25 x 1.38) at 300 dpi and the designs have sufficient margins.

Inverted

20070507 19:48:23 nbsp I couldnt open that .svg extension (maybe its just me, though) Users/MaryLieth
 Its in a format thats very convenient to a graphic designer. Its all good. The same thing is in the PDF. Users/WilliamLewis wl

20070508 01:02:40 nbsp Basically the same as before, but transparent except for the logo and the URL which would be black, it would also be the kind that sticks to glass I belive by static cling or something. it should say daviswiki.org/.............
where people can put something with sharpie or whatnot, ideal for business etc. Im will whip up this idea sooner or later, please not it wouldnt be periods but a underline, but the wiki hates me so it underlines the rest of my commentUsers/StevenDaubert
  Post Script: If I win I would use my pile of stickers to go around to business and have them pimp for themselves  and the wiki as noted in my wiki proposal. Meybe I should link my wiki proposal but I will do so later
   oooo, i like that idea. fillintheblank is great. Users/michelleaccurso michelle

20070509 12:14:30 nbsp I just want a sticker of the old design for my bike. :\ Users/WilliamLewis
    20070510 08:33:17 nbsp go to the sticker request page and make a seemingly futile request for one!  Users/StevenDaubert

20070510 14:15:15 nbsp How about.... (First line) What dont you know about Davis?
(Second line) DavisWiki logo Users/stevie13

20070510 14:42:54 nbsp What about some kind of Read about it on the Davis Wiki campaign?  In other words, when people go to a new business or organization because of something they read on the DW, they should tell that business or organization.  I did this recently with a local hotel, and the response was, more or less, hunh?  That shouldnt happen.  A read about it on the Davis Wiki campaign would encourage people to talk up the wiki via word of mouth;  businesses and organizations would definitely take notice, since they are themselves trying to figure how best to spread the word.. Users/CovertProfessor

20070523 12:56:06 nbsp Nobody submitted a design!  Even if somebody would recreate the old sticker with the proper dimensions itd help get some out faster! Users/PhilipNeustrom
   I was going to have David Poole whip up my concept sometime in the near future, but JA pretty much hit the nail on the head...

20070523 19:50:05 nbsp Can we combine the second line on Dauberts with the first entry, it is more of a request for feedback that way maybe?  Users/DavidPoole

20070525 10:06:47 nbsp any feedback on the designs I submitted? I hope theyre used! ) Users/EliseKane

20070525 13:36:23 nbsp I like the cow. Users/JabberWokky

20070526 00:56:07 nbsp Needs some color, bolder text, sharper picture lines... Users/Jedron
  Excellent  could you upload your suggested changes?

20070526 11:19:57 nbsp Suggestion:  A tandem bike (a la Elises excellent bike idea) with the motto Davis Wiki: Ride Together or Davis Wiki: The More the Merrier.  The bike idea is really fantastic, Elise.  I just thought a tandem bike might capture the multiple people working together idea.  Single or tandem, though, the bike logo is great. Users/JabberWokky
  I have this vision of an edit war being depicted as a bunch of bikes riding in a constant circle around a UCD bike traffic circle. Users/WesHardaker

20070608 22:55:06 nbsp So has a design been chosen? Whens it bein printed? Users/EliseKane

20070617 22:06:07 nbsp Ill say that http://daviswiki.org/Wiki_Sticker/StickiWiki_Design_Contest?actionFiles&doview&targetShare.png this is the best here, stickerwise.  It gives a hint of what Davis Wiki is.  Ill start looking into printing it.   I also like the bike one, but it needs a lot more polish before it could be printed.  Users/PhilipNeustrom

20070617 22:45:22 nbsp Philip, what kind of polish are we talking? Users/EliseKane

20070618 13:49:19 nbsp More color would make it stand out and be noticed much more. Users/Jedron

20070618 14:29:35 nbsp I still think mine is the best. FU B.R.! Users/JamesSchwab

20070619  What do people think of this color combination? Users/EliseKane

  I prefer the the black and white that you did originally  I think it will be more visible and more striking.  I like your design quite a bit, actually.  The only change I would make would be to smooth out the bike graphic a little bit.  Maybe add a dash symbol, i.e: FIND IT FAST – DAVISWIKI.ORG (with the dash on the first line).  Users/CovertProfessor

  The graphic needs smoothing out and the bike and text dont really go together crispness wise.
    Yes... a distressed font would work well.  Or something like a schoolbook simple rounded serif font.  I flashed to a version that reads Davis and Wiki above, and underneath Goes together like a Bicycle and your Ass.  But I dont think that works well.  http://cheshirehall.net/link/ws/wink.png  Still, something drawing a comparison to how the wiki fits Davis like Bicycles fit Davis would work well. Users/JabberWokky

20070625 19:12:54 nbsp Can somebody upload a simple replication of the original sticker design (see Wiki Sticker) that is completely ready for print (see http://www.123stickers.com/index.asp?PageActionVIEWCATS&Category202 for some details, perhaps)?  If someone does this Ill send it off to be printed.. Users/PhilipNeustrom

20070625 19:14:26 nbsp How do people feel about a square sticker with just the roundabout?  Perhaps some text someplace around / near it..? Users/PhilipNeustrom

20070626 13:24:43 nbsp I think its a waste that a few people, myself included, spent time designing stickers and now a decision is being made to print the original boring sticker. If decisionmakers are going to go with their own ideas, then why create a design contest in the first place? Users/EliseKane

   Because none of the images uploaded on this page are print quality.  Upload a print quality sticker to the specs of the sticker site (thats linked here) and I will immediately go and order stickers.  Nobody has uploaded a sticker that is printable, even as of now.   The requirements are roughly:  a single color (black and white is absolutely the cheapest  going up to two colors costs about $100 more), crisp, no borders (looks odd on a sticker) and looks good, and of the proper dimensions.  I could spend time to fix up some of these designs to make them printable, but thats why we are having a contest  so I dont have to.  :/  Short of that, I am thinking we should just print up the old design as none of the more exciting ones are printable right now and demand for stickers is through the roof.  Though, reading through this page, the requirements were not made clear.  Ill hold off reordering the old sticker for a few days to see if anyone uploads something print quality before then.  Users/PhilipNeustrom

20070628 21:34:54 nbsp I like Dauberts fill in the blanks suggestion the best. And I feel that it could work on a sticker.  It probably isnt the only one we should make now, though. Users/BrentLaabs

20070628 22:19:47 nbsp Davis Wiki  Why everyone knows more about Davis than you do.
DavisWiki.Org Users/EdWins

20070628 22:29:02 nbsp Remember that were making stickers, not billboards. Even a simple sentence would be too much for these stickers as wed have to shrink the text down to the point where youd have to be rather close to the sticker in order to read it. Also, some of the designs are incredibly intricate. At a distance, its all mush and hard to make sense of. Also, the sticker process that would be used doesnt take kindly to pixel manipulation to make textures.

Remember: The simpler, the better. We want these things to do their work, even at a distance. Users/WilliamLewis

20070629 08:52:23 nbsp I vote for the orginalishinvereted.
 Users/JamesSchwab

20070629 14:06:47 nbsp I was originally thought of it being like http://Daviswiki.org/ with the {{{_______}}} being below it, and when combined makes the URL for the specific page that sticker has been customized. This way it can be slapped randomly with nothing, or put on a pay phone with the payphone project DW url, or on a store with that stores wiki page etc. Users/StevenDaubert

20070629 16:56:41 nbsp Hiya. Can anyone show me the link for the proper dimensions of the sticker? I would be grateful. Users/EliseKane


text:
      SUBROUTINE GRAMAR
C     READ AND INTERPRETE THE META-GRAMMAR
C     2022-02-07: no XML
C     GF 01.08.1980: ']' TERMINATES 'RULE' TOO
C     GF 27.07.1980: WITH 'EOSTMT,AXIOM'
C     GF 19.07.1980: WITHOUT 'SETS', 'GRASYB', 'FMEM'-CHECK FOR 'SEMWRI'
C           PARAMETER 'GRANOU' SETS NO 'UNCH'
C     GF 12.07.1980 : INT PARASK, TEST,MEMHIB=... DELETED
C
      INCLUDE 'PARS.f'
      INCLUDE 'MEMS.f'
      INCLUDE 'PROS.f'
      INCLUDE 'SEMS.f'
      INCLUDE 'STRS.f'
      INCLUDE 'SYMS.f'
      INTEGER*2 PARASK
      INTEGER*2 I
     = ,ADDSUB   ! = 1(0) : ADD (SUBTRACT) A PRODUCTION
     = ,ENTITY   ! RESULT OF 'SCAN'
     = ,LEFT     ! THE LEFT SIDE OF A PRODUCTION
     = ,LNG      ! THE LENGTH OF A PRODUCTION
     = ,MON      ! -> MEMBER NUMBER ONE OF A PRODUCTION
     = ,PROD     ! -> A PRODUCTION
     = ,SEMA     ! -> THE SEMANTICS FOR A PRODUCTION
     = ,SYMBOL   ! RESULT OF 'SCAN'
     = ,SINGLE   ! = 1 (0) IF EVERY PRODUCTION IS (NOT) PROCESSED
      INTEGER*2 TEST
C
      SINGLE = PARASK('SINGLE',1,6,0)
      TEST   = PARASK('GRAMAR',1,6,0)

CAXI  CALL LRINIT

      CALL SCAINI
      CALL SCANR(I,EOFILE)
      CALL SCANR(I,IDENT)
      CALL SCANR(I,NUMBER)
      CALL SCANR(I,STRIN)
      CALL SPEINS('=>',1,2,ARROW)
      CALL SPEINS('|' ,1,1,BAR)
      CALL SPEINS('=' ,1,1,EQUALS)
      CALL SPEINS('-' ,1,1,MINUS)
      CALL SPEINS('.' ,1,1,PERIOD)
      CALL SPEINS(']' ,1,1,BUS)
      CALL SPEINS('[' ,1,1,SUB)
      CALL SPEINS('#' ,1,1,NUMSIG)
      CALL SPEINS('(' ,1,1,OPENP)
      CALL SPEINS(')' ,1,1,CLOSP)
      CALL SPEINS('+' ,1,1,PLUS)
      CALL SPEINS('*' ,1,1,TIMES)
      CALL SPEINS('/' ,1,1,DIVIDE)
      CALL SPEINS(';' ,1,1,I)
      CALL SCANR (I,EOSTMT)
      AXIOM = FSYM
C
CAXI  CALL LRAXIO ()
      SYMBOL = I ! was I 2005-04-07
      SYMBOL = 0 ! was I 2005-04-07
2     IF(SYMBOL .EQ. SUB .OR. SYMBOL .EQ. EOFILE) GOTO 1
C       PROCESS THE REST OF THE INTERFACE TO THE SCANNER
        CALL SCANR(SYMBOL,ENTITY)
        IF (SYMBOL .NE. STRIN) GOTO 3
          CALL SPEINS(STRNG,SYMPOS(ENTITY),SYMEND(ENTITY),I)
3       CONTINUE ! STRIN
      GOTO 2
1     CONTINUE ! SKIP INTERFACE
      GOTO 100
101   CONTINUE
        CALL ASSERT(50,ENTITY,-ENTITY)
100   CONTINUE
        CALL SCANR (SYMBOL,ENTITY)
C       CALL XML2 ('GRAMAR', 1, 6, SYMBOL,ENTITY)
130   CONTINUE
5       IF (SYMBOL .GT. BUS) GOTO 107 ! LIDEN
          GOTO (101,101,101,101,101
     =     ,106,107,108,109,110,111,112,113,114,115),SYMBOL
C           EOF IDE NUM STR  =>  |   =   -   .   ]
109   CONTINUE
        CALL SPEINS(STRNG,SYMPOS(ENTITY),SYMEND(ENTITY),I)
        ENTITY = I
108   CONTINUE
107   CONTINUE
        ! STORE THE ENTITY IN 'MEM'
        IF (FMEM .GE. MEMHIB) GOTO 7
          MEM(FMEM) = ENTITY
          FMEM = FMEM + 1
        GOTO 8 ! NO 'MEM'-OVERFLOW
7       CONTINUE
          CALL ASSERT(57,FMEM,ENTITY)
8       CONTINUE
        GOTO 100
C
112   CONTINUE
        ADDSUB = 1
131   CONTINUE
        FMEM = FMEM - 1 ! DELETE LEFT SIDE
        LEFT = MEM(FMEM)
        MON = FMEM
        SEMA = 0
        GOTO 100
113   CONTINUE
        ADDSUB = 0
        GOTO 131
110   CONTINUE
        CALL SCANR(SYMBOL,ENTITY)
        IF (SYMBOL .NE. NUMSIG) GOTO 9
          CALL SCANR(SYMBOL,ENTITY)
          IF (SYMBOL .NE. NUMBER) GOTO 11
            CALL SYMNUM(ENTITY,SEMA)
11        CONTINUE ! NUMBER
        GOTO 10 ! = NUMSIG
9       CONTINUE
          CALL ASSERT(55,SYMBOL,NUMSIG)
10      CONTINUE
12      IF(SYMBOL .EQ. PERIOD .OR. SYMBOL .EQ. BAR
     =    .OR. SYMBOL .EQ. BUS
     =  ) GOTO 13
          CALL SCANR(SYMBOL,ENTITY)
        GOTO 12
13      CONTINUE
        GOTO 130
111   CONTINUE
114   CONTINUE
115   CONTINUE ! ']'
        LNG = FMEM - MON
        CALL PROCHA(ADDSUB,LEFT,MON,LNG,PROD)
        IF (SINGLE .NE. 1) GOTO 14
C         CALL REORG
14      CONTINUE
        IF (SEMA .LE. 0) GOTO 15
          PROSEM(PROD) = FSEM
          CALL SEMWRI(0,SEPR,SEMA)
          CALL SEMWRI(0,EOS,0)
          SEMA = 0
        GOTO 18 ! SEMA .LE. 0
15      CONTINUE
          PROSEM(PROD) = PNOSEM + PARASK('GRANOU',1,6,0)
18      CONTINUE
        MON = FMEM
        IF (SYMBOL .NE. BUS)
     =    GOTO 100 ! FOR ']' FALL THROUGH
106   CONTINUE
        RETURN
      END

max_stars_repo_name: gfis/jextra

max_issues_repo_name: gfis/jextra

max_forks_repo_name: gfis/jextra


text:
! --------------------------------------------------------------------
!
! Copyright (C) 2015 Rocco Meli
!
! This program is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
! ---------------------------------------------------------------------

PROGRAM input_test

    USE INPUT
    USE OUTPUT

    IMPLICIT NONE

    ! System an basis set sizes
    INTEGER :: Ne   ! Number of electrons
    INTEGER :: Nn   ! Number of nuclei
    INTEGER :: K    ! Basis set size
    INTEGER :: c    ! Contractions (TODO: split valence basis set)

    ! System and basis set informations
    REAL*8, allocatable, dimension(:,:) :: R            ! Atomic potisions
    INTEGER, allocatable, dimension(:) :: Z             ! Atomic charges
    REAL*8, allocatable, dimension(:,:) :: basis_R      ! Basis functions' centers
    INTEGER, allocatable, dimension(:,:) :: basis_L     ! Basis functions' angular momenta
    REAL*8, allocatable, dimension(:,:) :: basis_A      ! Contraction exponential coefficients
    REAL*8, allocatable, dimension(:,:) :: basis_D      ! Conttaction linear coefficients

    CALL load("tests/N2_f.in",Ne,Nn,K,c,R,Z,basis_R,basis_L,basis_A,basis_D)

    CALL print_real_matrix(Nn,3,R)
    CALL print_integer_matrix(Nn,1,Z)
    CALL print_real_matrix(K,3,basis_R)
    CALL print_integer_matrix(K,3,basis_L)
    CALL print_real_matrix(K,c,basis_A)
    CALL print_real_matrix(K,c,basis_D)

    DEALLOCATE(R,Z,basis_R,basis_L,basis_A,basis_D)

END PROGRAM input_test

max_stars_repo_name: RMeli/Hartree-Fock

max_issues_repo_name: RMeli/Hartree-Fock

max_forks_repo_name: RMeli/Hartree-Fock


text:
Yolo Atari Club was a retrogaming and Atari enthusiast society that formed in the 1990s and apparently disbanded at some point within the following decade. The club held sporadic gaming and programming events. A paid membership included access to a software lending library, a newsletter, and other Atari resources.

Although seemingly unconnected, the ongoing series of Atari Party Atari Parties fills some of the void left behind by the old clubs demise.


text:
${FPIO}/rtl/fpio_fifo_if.sv
${FPIO}/rtl/fpio_fifo.sv
${FPIO}/rtl/fpio_if.sv
${FPIO}/rtl/fpio_tx.sv
${FPIO}/rtl/fpio.sv

max_stars_repo_name: mballance/fpio

max_issues_repo_name: mballance/fpio

max_forks_repo_name: mballance/fpio


text:
!-----------------------------------------------------------------------
!     Function:      dflux_simpson
!     Authors:       S. Lazerson (lazerson@pppl.gov)
!     Date:          03/02/2012
!     Description:   This subroutine calculates the flux due to a line
!                    segment using simpson formula.
!-----------------------------------------------------------------------
      REAL FUNCTION dflux_simpson(x0,y0,z0,x1,y1,z1,cg)
!-----------------------------------------------------------------------
!     Libraries
!-----------------------------------------------------------------------
      USE stel_kinds, ONLY: rprec
      USE virtual_casing_mod
      USE biotsavart
      USE bsc_T
!-----------------------------------------------------------------------
!     Local Variables
!          ier            Error Flag
!          iunit          File ID Number
!-----------------------------------------------------------------------
      IMPLICIT NONE
      INTEGER, intent(in) :: cg
      REAL(rprec)   , intent(in) :: x0,y0,z0,x1,y1,z1
      REAL(rprec) :: delf,delt,dx,dy,dz,xp,yp,zp,ax,ay,az
      REAL(rprec) :: xvec(3),avec(3)
      INTEGER     :: i,ier
      integer, parameter :: nop=3
      real(rprec), dimension(nop), parameter :: ci=(/1./6.,2./3.,1./6./)
!-----------------------------------------------------------------------
!     Begin Function
!-----------------------------------------------------------------------

      delf = 0
      delt = REAL(1)/REAL(nop-1)
      dx = x1-x0
      dy = y1-y0
      dz = z1-z0
      IF (cg > 0) THEN
         DO i = 1, nop
            xvec(1) = x0+(i-1)*dx*delt
            xvec(2) = y0+(i-1)*dy*delt
            xvec(3) = z0+(i-1)*dz*delt
            CALL bsc_a(coil_group(cg),xvec,avec)
            delf = delf+ ci(i)*(avec(1)*dx+avec(2)*dy+avec(3)*dz)
         END DO
      ELSE
         DO i = 1, nop
            ier = 1; ax = 0; ay = 0; az = 0
            xp = x0+(i-1)*dx*delt
            yp = y0+(i-1)*dy*delt
            zp = z0+(i-1)*dz*delt
            CALL vecpot_vc(xp,yp,zp,ax,ay,az,ier)
            delf = delf+ ci(i)*(ax*dx+ay*dy+az*dz)
         END DO
      END IF
      dflux_simpson = delf
      RETURN
!-----------------------------------------------------------------------
!     End Function
!-----------------------------------------------------------------------
      END FUNCTION dflux_simpson

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT

max_forks_repo_name: joseluisvelasco/STELLOPT


text:
      SUBROUTINE speval(n,xk,c,t,f,ifail)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE stel_constants
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER :: n, ifail
      REAL(rprec), DIMENSION(0:*) :: xk, c
      REAL(rprec), DIMENSION(*) :: f
      REAL(rprec) :: t
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: i, j, l, ip1
      REAL(rprec) :: a10, a11, a12, a20, a21
      REAL(rprec) :: b02, b12, b03, b13, b23, b04, b14, b24, b34
      REAL(rprec) :: d10, d11, d12, d20, d21, d30
      REAL(rprec) :: tm0, tm1, tm2, tp1, tp2, tp3

!     evaluate the cubic spline with knots xk and b-spline
!     coefficients c at the point t. return function value
!     in f(1) and derivatives in f(2)-f(4).

      ifail=1
      DO j=1,4
         f(j)=0
      END DO
      IF((t.lt.xk(4)) .or. (t.gt.xk(n+1))) RETURN
      i=4
      ip1=n+1
   20 l=(i+ip1)/2
      IF(ip1-i.le.1) GO TO 40
      IF(t.lt.xk(l)) GO TO 30
      i=l
      GO TO 20
   30 ip1=l
      GO TO 20
   40 tm2=t-xk(i-2)
      tm1=t-xk(i-1)
      tm0=t-xk(i)
      tp1=xk(i+1)-t
      tp2=xk(i+2)-t
      tp3=xk(i+3)-t
      d10=tp1+tm0
      d11=tp1+tm1
      d20=tp2+tm0
      d12=tp1+tm2
      d21=tp2+tm1
      d30=tp3+tm0
      b12=tp1/d10
      b02=tm0/d10
      b23=tp1*b12/d11
      b13=tm1*b12/d11+tp2*b02/d20
      b03=tm0*b02/d20
      b34=tp1*b23/d12
      b24=tm2*b23/d12+tp2*b13/d21
      b14=tm1*b13/d21+tp3*b03/d30
      b04=tm0*b03/d30
      f(1)=c(i-3)*b34+c(i-2)*b24+c(i-1)*b14+c(i)*b04
      a12=(c(i-2)-c(i-3))/d12
      a11=(c(i-1)-c(i-2))/d21
      a10=(c(i)-c(i-1))/d30
      f(2)=3*(a12*b23+a11*b13+a10*b03)
      a21=(a11-a12)/d11
      a20=(a10-a11)/d20
      f(3)=6*(a21*b12+a20*b02)
      f(4)=6*(a20-a21)/d10
      ifail=0

      END SUBROUTINE speval


text: #redirect American Hotel


! Make sure that frontend driver options that require arguments are
! correctly rejected when the argument value is missing.

! REQUIRES: new-flang-driver

!-----------
! RUN lines
!-----------
! RUN: not %flang_fc1 -E %s -o 2>&1 | FileCheck %s
! RUN: not %flang_fc1 -E %s -U 2>&1 | FileCheck %s
! RUN: not %flang_fc1 -E %s -D 2>&1 | FileCheck %s
! RUN: not %flang_fc1 -E %s -I 2>&1 | FileCheck %s
! RUN: not %flang_fc1 -E %s -J 2>&1 | FileCheck %s
! RUN: not %flang_fc1 -E %s -module-dir 2>&1 | FileCheck %s
! RUN: not %flang_fc1 -E %s -module-suffix 2>&1 | FileCheck %s
! RUN: not %flang_fc1 -E %s -fintrinsic-modules-path 2>&1 | FileCheck %s

!-----------------------
! EXPECTED OUTPUT
!-----------------------
! CHECK: error: argument to '-{{.*}}' is missing (expected 1 value)


text:
      SUBROUTINE DEOFTS7(X,NROW,NCOL,NOBS,MSTA,XMSG,NEVAL,EVEC,JOPT,
     +                   IFLAG,XX,WRK,EVECTS,EVTSAV,IER)
      DOUBLE PRECISION XBAR
      DOUBLE PRECISION XVAR
      DOUBLE PRECISION XSTD
      DOUBLE PRECISION CON

c f77

c takes output of prneof.f and creates a time series of the
c .   eof amplitudes

c nomenclature :
c .   x         - matrix containing the data.  it contains n
c .               observations for each of m stations or grid pts.
c .   nrow,ncol - exact row (observation) and column (station)
c .               dimensions of x in the calling routine.
c .   nobs      - actual number of observations (nobs <= nrow)
c .   msta      - actual number of stations     (msta <= ncol)
c .   xmsg      - missing code (if no obs are missing set to some
c .               number which will not be encountered)
c .   neval     - no. of eigenvalues and eigen vectors computed by
c .               prneof.
c .   evec      - array created by prneof
c .               this must be dimensioned at least (ncol,neval) in the
c .               calling routine.
c .   jopt      - =0 means covariance  was used to generate evec
c .               =1 means correlation was used to generate evec
c .   xx        - work vector of the same dimensions as x
c .   wrk       - work vector of length nobs
c .   evects    - time series of eof amplitudes for each eigenvalue
c .   ier       - error code

C INPUT

C dimensions of x
      INTEGER NROW,NCOL,NOBS,MSTA

C monthly data from station/grid pt
C work array
C missing code (if any)
      DOUBLE PRECISION X(NROW,NCOL),EVEC(NCOL,NEVAL),XX(NROW,NCOL),
     +                 WRK(NOBS),XMSG
      INTEGER NEVAL,IFLAG

C OUTPUT

C monthly anomalies from long term monthly mean
      DOUBLE PRECISION EVECTS(NROW,NEVAL)
      DOUBLE PRECISION EVTSAV(NEVAL)
C error code
      INTEGER IER

      IER = 0
      IF (NROW.LE.0 .OR. NCOL.LE.0) IER = IER + 1
      IF (NOBS.LE.0 .OR. MSTA.LE.0) IER = IER + 10
      IF (IER.NE.0) THEN
          WRITE (*,FMT='(/,'' sub eofts7: ier='',5i5)') IER,NROW,NCOL,
     +      NOBS,MSTA
c          stop
          RETURN
      END IF

c set to msg as default

      DO M = 1,NCOL
          DO N = 1,NROW
              XX(N,M) = X(N,M)
          END DO
      END DO

c the following was added in Sept 2003 so that the
c 'normalization' that follows the "if" would
c not be done.

c c c if (iflag.eq.1 .and. jopt.eq.0) go to 10
      IF (JOPT.EQ.1) THEN

c calculate mean/stddev at each station/grid-pt

          DO M = 1,MSTA
              DO N = 1,NOBS
                  WRK(N) = X(N,M)
              END DO

              CALL DSTAT2(WRK,NOBS,XMSG,XBAR,XVAR,XSTD,KNTX,IER)
              CON = 1.0D0
              IF (JOPT.EQ.1 .AND. XSTD.GT.0.D0) THEN
                  CON = 1.D0/XSTD
              END IF

c         write (*,"(' eofts7: m,nobs,kntx,xbar,xstd=',3i5,2f8.2)")
c     *                        m,nobs,kntx,xbar,xstd

              DO N = 1,NOBS
                  IF (X(N,M).NE.XMSG .AND. XBAR.NE.XMSG) THEN
                      XX(N,M) = (X(N,M)-XBAR)*CON
                  ELSE
                      XX(N,M) = XMSG
                  END IF
              END DO

C end msta
          END DO

c c c write (*,"(//,'ANOMALIES',/)")
c c c do n=1,nobs
c c c     write (*,"(i5 , 10(1x,f9.3) )") n, (xx(n,m),m=1,msta)
c c c enddo

      END IF

c calculate the amplitude time series

   10 DO K = 1,NEVAL
          DO N = 1,NOBS
              KNTX = 0
              EVECTS(N,K) = 0.0D0
              DO M = 1,MSTA
                  IF (XX(N,M).NE.XMSG .AND. EVEC(M,K).NE.XMSG) THEN
                      KNTX = KNTX + 1
                      EVECTS(N,K) = EVECTS(N,K) + EVEC(M,K)*XX(N,M)
                  END IF
              END DO
              IF (KNTX.EQ.0) EVECTS(N,K) = XMSG
          END DO
          CALL DSTAT2(EVECTS(1,K),NOBS,XMSG,EVTSAV(K),XVAR,XSTD,KNTX,
     +                IER)
          IF (IFLAG.EQ.0 .AND. EVTSAV(K).NE.XMSG) THEN
              DO N = 1,NOBS
                  IF (EVECTS(N,K).NE.XMSG) THEN
                      EVECTS(N,K) = EVECTS(N,K) - EVTSAV(K)
                  END IF
              END DO
          END IF
      END DO

      RETURN
      END

      SUBROUTINE DEOFTSCA(XX,NROW,NCOL,NOBS,MSTA,XMSG
     +                   ,NEVAL,EVEC,JOPT,IFLAG,EVECTS,PCRIT,IER,X,WRK)
      IMPLICIT NONE

c takes output of prneof_ca.f and creates a time series of the
c .   eof amplitudes. It is essential that xx and pcrit not change
c .   from the call to prneof_ca.f

c use f90 to dynamically allocate/deallocate memory (work space)
c     also some array syntax

c nomenclature :
c .   xx        - matrix containing the data.  it contains n observations
c .               for each of m stations or grid pts.
c .   nrow,ncol - exact row (observation) and column (station)
c .               dimensions of x in the calling routine.
c .   nobs      - actual number of observations (nobs <= nrow)
c .   msta      - actual number of stations     (msta <= ncol)
c .   xmsg      - missing code (if no obs are missing set to some
c .               number which will not be encountered)
c .   neval     - no. of eigenvalues and eigen vectors computed by
c .               prneof.
c .   evec      - array created by prneof_ca
c .               this must be dimensioned at least (ncol,neval) in the
c .               calling routine.
c .   evects    - time series of eof amplitudes for each eigenvalue
c .   pcrit     - minimum % of non-missing values required for the
c .               station/grid-pt be used in the covariance/correlation
c .               matrix, must be between 0 and 100 inclusive
c .   ier       - error code
c INPUT

c dimensions of xx
      INTEGER NROW,NCOL,NOBS,MSTA
      INTEGER NEVAL, JOPT, IFLAG

c monthly data from station/grid pt
c missing code (if any)
      DOUBLE PRECISION XX(NROW,NCOL),EVEC(NCOL,NEVAL),XMSG
c OUTPUT

c monthly anomalies from long term monthly mean
      DOUBLE PRECISION EVECTS(NROW,NEVAL)
c error code
      INTEGER IER
c temporaray arrays
      DOUBLE PRECISION X(NROW,NCOL),WRK(NOBS)

      DOUBLE PRECISION PCRIT
      DOUBLE PRECISION PCRITX
      DOUBLE PRECISION XBAR
      DOUBLE PRECISION XVAR
      DOUBLE PRECISION XSTD
      DOUBLE PRECISION CON
      INTEGER NSTA, NC, NR, KNTX, M, N, K

      IER = 0

      PCRITX = PCRIT*0.01D0

c counts the total number of
      NSTA = 0
      DO NC = 1,MSTA
c counter for this location
          KNTX = 0
          DO NR = 1,NOBS
              IF (XX(NR,NC).NE.XMSG) THEN
                  KNTX = KNTX + 1
              END IF
          END DO
          IF (DBLE(KNTX)/DBLE(NOBS).GE.PCRITX) THEN
              NSTA = NSTA + 1
              DO N = 1,NOBS
                  X(N,NSTA) = XX(N,NC)
              END DO
          END IF
      END DO

c this if added Sept 2003 to allow computation without
c removing the mean.

      IF (JOPT.EQ.0 .AND. IFLAG.EQ.1) GO TO 10

c calculate each station/grid-pt  long-term mean and standard deviation
c use the "nsta" limit; this pertains to array x

      DO M = 1,NSTA
          DO N = 1,NOBS
              WRK(N) = X(N,M)
          END DO

          CALL DSTAT2(WRK,NOBS,XMSG,XBAR,XVAR,XSTD,KNTX,IER)
          CON = 1.0D0
          IF (JOPT.EQ.1 .AND. XSTD.GT.0.D0) THEN
              CON = 1.D0/XSTD
          END IF

          DO N = 1,NOBS
              IF (X(N,M).NE.XMSG .AND. XBAR.NE.XMSG) THEN
                  X(N,M) = (X(N,M)-XBAR)*CON
              ELSE
                  X(N,M) = XMSG
              END IF
          END DO
c end nsta
c
      END DO
   10 CONTINUE

c calculate the amplitude time series

      DO K = 1,NEVAL
          DO N = 1,NOBS
              KNTX = 0
              EVECTS(N,K) = 0.0D0
              DO M = 1,NSTA
                  IF (X(N,M).NE.XMSG) THEN
                      KNTX = KNTX + 1
                      EVECTS(N,K) = EVECTS(N,K) + EVEC(M,K)*X(N,M)
                  END IF
              END DO
              if (KNTX.EQ.0) EVECTS(N,K) = XMSG
          END DO
      END DO

      RETURN
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
      call suba(5hjunky)
      end
      subroutine suba(i)
      write(*,'(a)') i
      end

max_stars_repo_name: kurtsansom/ftnchek


text:
!!# MODULE: <USR_SimpleSet>
MODULE USR_SimpleSet

!!## PURPOSE
!! The scalar set <SimpleSet> is a simple pointer array of values
!! which are kept sorted so that duplicate values are
!! never added.  A user can <ADD> to the SET,
!! <REMOVE> from the SET and <FINALIZE> the SET.

!!## DETAILS
!! @ Items are added to the set in log_2(n) time,
!! @ removed from the set in log_2(n) time, and
!! @ values are found in log_2(n) time
!! @ the set is finalized in constant time.

!!## MODULES
!! @ function to return the index of a string in an array of strings
!! @ function to return the trimmed length (really trimmed SIZE) of an array of strings
!! @ reallocation routine
!! @ string clearing routine
USE FUN_INDEXa            !!((08-B-FUN_INDEXa.f90))
USE FUN_LEN_TRIMa         !!((03-A-FUN_LEN_TRIMa.f90))
USE SUB_Reallocate        !!((04-B-SUB_Reallocate.f90))
USE SUB_CLEAR             !!((04-A-SUB_CLEAR.f90))
USE FUN_SIZEa             !!((06-B-FUN_SIZEa.f90))
USE FUN_Error             !!((04-A-FUN_Error.f90))
USE FUN_EQUILOC           !!((03-A-FUN_EQUILOC.f90))
USE FUN_IntervalBisection !!((03-A-FUN_IntervalBisection.f90))
USE FUN_Default           !!((04-A-FUN_Default.f90))
USE VAR_Units             !!((03-A-VAR_Units.f90))
USE PRN_Table             !!((11-B-PRN_Table.f90))
USE FUN_STR               !!((05-B-FUN_STR.f90))
USE PAR_IntrinsicLengths  !!((02-A-PAR_IntrinsicLengths.f90))
USE KND_IntrinsicTypes    !!((01-A-KND_IntrinsicTypes.f90))

!!### DEFAULT IMPLICIT
IMPLICIT NONE

!!### DEFAULT ACCESS
PRIVATE

!!### INTERFACES
INTERFACE ADD_TO_SET
 MODULE PROCEDURE ADD_TO_SET_S
 MODULE PROCEDURE ADD_TO_SET_I
 MODULE PROCEDURE ADD_TO_SET_Rdp
END INTERFACE

INTERFACE ADDMANY_TO_SET
 MODULE PROCEDURE ADDMANY_TO_SET_I
END INTERFACE

INTERFACE REMOVE_FROM_SET
 MODULE PROCEDURE REMOVE_FROM_SET_S
 MODULE PROCEDURE REMOVE_FROM_SET_I
 MODULE PROCEDURE REMOVE_FROM_SET_rdp
END INTERFACE

INTERFACE FINALIZE_SET
 MODULE PROCEDURE FINALIZE_SET_S
 MODULE PROCEDURE FINALIZE_SET_I
 MODULE PROCEDURE FINALIZE_SET_Rdp
END INTERFACE

INTERFACE PRINT_SET
 MODULE PROCEDURE PRINT_SET_S
 MODULE PROCEDURE PRINT_SET_I
 MODULE PROCEDURE PRINT_SET_Rdp
END INTERFACE

INTERFACE FIND_IN_SET
 MODULE PROCEDURE FIND_IN_SET_S
 MODULE PROCEDURE FIND_IN_SET_I
 MODULE PROCEDURE FIND_IN_SET_Rdp
END INTERFACE

!!### PUBLIC ACCESS SET
PUBLIC :: PRINT_SET
PUBLIC :: ADD_TO_SET
PUBLIC :: ADDMANY_TO_SET
PUBLIC :: REMOVE_FROM_SET
PUBLIC :: FINALIZE_SET
PUBLIC :: FIND_IN_SET

!!### PARAMETERS
INTEGER,PARAMETER :: DEFAULT_dN=5

!!## CONTAINED ROUTINES
CONTAINS

!!### SUBROUTINE: ADD_TO_SET_S
SUBROUTINE ADD_TO_SET_S( SET , Entry )
!!#### PURPOSE
!! Add an Entry to the SET.

!!#### REQUIRED OUTPUT
CHARACTER(*),POINTER :: SET(:)

!!#### REQUIRED INPUT
CHARACTER(*),INTENT(IN) :: Entry

!!#### LOCAL VARIABLES
INTEGER :: l(1:2),N

!!--begin--
IF( .NOT.ASSOCIATED(SET) )THEN
 ALLOCATE( SET(1:DEFAULT_dN) )
 CALL CLEAR(SET)
 N = 0
 SET(1) = Entry

ELSE
 !get number of elements in set
 N = LEN_TRIMa(SET)

 !easy setup if N=0
 IF( N==0 )THEN
  SET(1) = Entry
  RETURN
 END IF

 !get interval
 l = IntervalBisection(Entry,SET(1:N))

 !first check if we have a duplicate value
 IF( l(1)/=0 )THEN
  IF( SET(l(1))==Entry )RETURN
 END IF

 !increment
 N = N + 1

 !reallocate set if full
 IF( N>=SIZE(SET) )THEN
  CALL Reallocate( SET , DEFAULT_dN , fill=REPEAT(" ",LEN(SET)))
 END IF

 !bump
 SET(l(2)+1:N) = SET(l(2):N-1)

 !set
 SET(l(2)) = Entry

END IF

!!--end--
END SUBROUTINE

!!### FUNCTION: FIND_IN_SET_S
FUNCTION FIND_IN_SET_S(SET,Entry) RESULT(I)
!!#### PURPOSE
!! Find an entry in the set.

!!#### REQUIRED INPUT
CHARACTER(*),POINTER    :: SET(:)
CHARACTER(*),INTENT(IN) :: Entry

!!#### REQUIRED OUTPUT
INTEGER :: I

!!#### LOCAL VARIABLES
INTEGER :: l(1:2),N

!!--begin--
I = 0

IF( ASSOCIATED(SET) )THEN

 N = LEN_TRIMa(SET)

 IF( N==0 )RETURN

 l = IntervalBisection(Entry,SET(1:N))

 IF( l(1)<1 )RETURN

 IF( SET(l(1))==Entry )THEN
  I = l(1)
 END IF

END IF

!!--end--
END FUNCTION

!!### SUBROUTINE: REMOVE_FROM_SET_S
SUBROUTINE REMOVE_FROM_SET_S(SET,Entry)
!!#### PURPOSE
!! Remove an entry from the SET.

!!#### REQUIRED OUTPUT
CHARACTER(*),POINTER    :: SET(:)

!!#### REQUIRED INPUT
CHARACTER(*),INTENT(IN) :: Entry

!!#### LOCAL VARIABLES
INTEGER :: N,I

!!--begin--
IF( ASSOCIATED(SET) )THEN
 I = FIND_IN_SET(Set,Entry)

 IF( I==0 )RETURN

 N = SIZE(SET)

 SET(I:N-1) = SET(I+1:N)

 CALL Clear(SET(N))

END IF

!!--end--
END SUBROUTINE

!!### SUBROUTINE: PRINT_SET_S
SUBROUTINE PRINT_SET_S( SET , Unit )
!!#### PURPOSE
!! Print out the contents of a set.

!!#### REQUIRED INPUT/OUTPUT
CHARACTER(*),POINTER :: SET(:)

!!#### OPTIONAL INPUT
INTEGER,OPTIONAL,INTENT(IN) :: Unit

!!#### LOCAL VARIABLES
INTEGER :: i,Unit_
CHARACTER(LEN(SET)+LEN_I),POINTER :: D(:,:)

!!--begin--
Unit_ = DEFAULT(DEFAULT_OUTPUT_UNIT,Unit)

IF( .NOT.ASSOCIATED(SET) )THEN
 WRITE(Unit_,"(a)")" N U L L   S E T"
 RETURN
ELSE
 ALLOCATE( D(1+SIZE(SET),2) )
 CALL CLEAR(D)
 D(1,:) = (/"index","value"/)
 DO i=1,SIZE(SET)
  D(i+1,1) = TRIM(STR(i))
  D(i+1,2) = TRIM(SET(i))
 END DO
 CALL PRINT_Table(D,Unit=Unit_)
END IF

!!--end--
END SUBROUTINE

!!### SUBROUTINE: FINALIZE_SET_S
SUBROUTINE FINALIZE_SET_S( SET )
!!#### PURPOSE
!! Finalize the SET by allocation to its LEN_TRIMa length.

!!#### REQUIRED INPUT/OUTPUT
CHARACTER(*),POINTER :: SET(:)

!!#### LOCAL VARIABLES
INTEGER :: dn,N

!!--begin--
IF( ASSOCIATED(SET) )THEN
 N = LEN_TRIMa(SET)
 IF( N==0 )THEN
  DEALLOCATE(SET)
  NULLIFY(SET)
 ELSE
  dn=LEN_TRIMa(SET)-SIZE(SET)
  CALL Reallocate( SET , dn )
 END IF
ENDIF

!!--end--
END SUBROUTINE

!!### SUBROUTINE: ADD_TO_SET_I
SUBROUTINE ADD_TO_SET_I( SET , Entry )
!!#### PURPOSE
!! Add an Entry to the SET.

!!#### REQUIRED OUTPUT
INTEGER,POINTER :: SET(:)

!!#### REQUIRED INPUT
INTEGER,INTENT(IN) :: Entry

!!#### LOCAL VARIABLES
INTEGER :: l(1:2),N

!!--begin--

IF( .NOT.ASSOCIATED(SET) )THEN
 ALLOCATE( SET(1:DEFAULT_dN) )
 SET = ERROR(1)
 N = 0
 SET(1) = Entry

ELSE
 !get number of elements in set
 N = SIZEa(SET)

 !easy setup if N=0
 IF( N==0 )THEN
  SET(1) = Entry
  RETURN
 END IF

 !get interval
 l = IntervalBisection(Entry,SET(1:N))

 !first check if we have a duplicate value
 IF( l(1)/=0 )THEN
  IF( SET(l(1))==Entry )RETURN
 END IF

 !increment
 N = N + 1

 !reallocate set if full
 IF( N>=SIZE(SET) )THEN
  CALL Reallocate( SET , DEFAULT_dN , fill=ERROR(1) )
 END IF

 !bump
 SET(l(2)+1:N) = SET(l(2):N-1)

 !set
 SET(l(2)) = Entry

END IF

!!--end--
END SUBROUTINE

!!### SUBROUTINE: ADDMANY_TO_SET_I
SUBROUTINE ADDMANY_TO_SET_I( SET , Entries )
!!#### PURPOSE
!! Add many Entries to the SET.

!!#### REQUIRED OUTPUT
INTEGER,POINTER :: SET(:)

!!#### REQUIRED INPUT
INTEGER,INTENT(IN) :: Entries(:)

!!#### LOCAL VARIABLES
INTEGER :: I

!!--begin--

DO I=1,SIZE(Entries)
 CALL ADD_TO_SET( SET , Entries(I) )
END DO

!!--end--
END SUBROUTINE

!!### FUNCTION: FIND_IN_SET_I
FUNCTION FIND_IN_SET_I(SET,Entry) RESULT(I)
!!#### PURPOSE
!! Find an entry in the set.

!!#### REQUIRED INPUT
INTEGER,POINTER    :: SET(:)
INTEGER,INTENT(IN) :: Entry

!!#### REQUIRED OUTPUT
INTEGER :: I

!!#### LOCAL VARIABLES
INTEGER :: l(1:2),N

!!--begin--
I = 0

IF( ASSOCIATED(SET) )THEN

 N = SIZEa(SET)

 IF( N==0 )RETURN

 l = IntervalBisection(Entry,SET(1:N))

 IF( l(1)<1 )RETURN

 IF( SET(l(1))==Entry )THEN
  I = l(1)
 END IF

END IF

!!--end--
END FUNCTION

!!### SUBROUTINE: REMOVE_FROM_SET_I
SUBROUTINE REMOVE_FROM_SET_I(SET,Entry)
!!#### PURPOSE
!! Remove an entry from the SET.

!!#### REQUIRED OUTPUT
INTEGER,POINTER    :: SET(:)

!!#### REQUIRED INPUT
INTEGER,INTENT(IN) :: Entry

!!#### LOCAL VARIABLES
INTEGER :: N,I

!!--begin--
IF( ASSOCIATED(SET) )THEN
 I = FIND_IN_SET(Set,Entry)

 IF( I==0 )RETURN

 N = SIZE(SET)

 SET(I:N-1) = SET(I+1:N)

 SET(N) = ERROR(1)

END IF

!!--end--
END SUBROUTINE

!!### SUBROUTINE: PRINT_SET_I
SUBROUTINE PRINT_SET_I( SET , Unit )
!!#### PURPOSE
!! Print out the contents of a set.

!!#### REQUIRED INPUT/OUTPUT
INTEGER,POINTER :: SET(:)

!!#### OPTIONAL INPUT
INTEGER,OPTIONAL,INTENT(IN) :: Unit

!!#### LOCAL VARIABLES
INTEGER :: i,Unit_
CHARACTER(LEN_I),POINTER :: D(:,:)

!!--begin--
Unit_ = DEFAULT(DEFAULT_OUTPUT_UNIT,Unit)

IF( .NOT.ASSOCIATED(SET) )THEN
 WRITE(Unit_,"(a)")" N U L L   S E T"
 RETURN
ELSE
 ALLOCATE( D(1+SIZE(SET),2) )
 CALL CLEAR(D)
 D(1,:) = (/"index","value"/)
 DO i=1,SIZE(SET)
  D(i+1,1) = TRIM(STR(i))
  D(i+1,2) = TRIM(STR(SET(i)))
 END DO
 CALL PRINT_Table(D,Unit=Unit_)
END IF

!!--end--
END SUBROUTINE

!!### SUBROUTINE: FINALIZE_SET_I
SUBROUTINE FINALIZE_SET_I( SET )
!!#### PURPOSE
!! Finalize the SET by allocation to its SIZEa length.

!!#### REQUIRED INPUT/OUTPUT
INTEGER,POINTER :: SET(:)

!!#### LOCAL VARIABLES
INTEGER :: dn,N

!!--begin--
IF( ASSOCIATED(SET) )THEN
 N = SIZEa(SET)
 IF( N==0 )THEN
  DEALLOCATE(SET)
  NULLIFY(SET)
 ELSE
  dn=SIZEa(SET)-SIZE(SET)
  CALL Reallocate( SET , dn )
 END IF
ENDIF

!!--end--
END SUBROUTINE

!!### SUBROUTINE: ADD_TO_SET_Rdp
SUBROUTINE ADD_TO_SET_Rdp( SET , Entry )
!!#### PURPOSE
!! Add an Entry to the SET.

!!#### REQUIRED OUTPUT
REAL(KIND_Rdp),POINTER :: SET(:)

!!#### REQUIRED INPUT
REAL(KIND_Rdp),INTENT(IN) :: Entry

!!#### LOCAL VARIABLES
INTEGER :: l(1:2),N

!!--begin--
IF( .NOT.ASSOCIATED(SET) )THEN
 ALLOCATE( SET(1:DEFAULT_dN) )
 SET = ERROR(1._KIND_Rdp)
 N = 0
 SET(1) = Entry

ELSE
 !get number of elements in set
 N = SIZEa(SET)

 !easy setup if N=0
 IF( N==0 )THEN
  SET(1) = Entry
  RETURN
 END IF

 !get interval
 l = IntervalBisection(Entry,SET(1:N))

 !first check if we have a duplicate value
 IF( l(1)/=0 )THEN
  IF( SET(l(1))==Entry )RETURN
 END IF

 !increment
 N = N + 1

 !reallocate set if full
 IF( N>=SIZE(SET) )THEN
  CALL Reallocate( SET , DEFAULT_dN , fill=ERROR(1._KIND_Rdp) )
 END IF

 !bump
 SET(l(2)+1:N) = SET(l(2):N-1)

 !set
 SET(l(2)) = Entry

END IF

!!--end--
END SUBROUTINE

!!### FUNCTION: FIND_IN_SET_Rdp
FUNCTION FIND_IN_SET_Rdp(SET,Entry) RESULT(I)
!!#### PURPOSE
!! Find an entry in the set.

!!#### REQUIRED INPUT
REAL(KIND_Rdp),POINTER    :: SET(:)
REAL(KIND_Rdp),INTENT(IN) :: Entry

!!#### REQUIRED OUTPUT
INTEGER :: I

!!#### LOCAL VARIABLES
INTEGER :: l(1:2),N

!!--begin--
I = 0

IF( ASSOCIATED(SET) )THEN

 N = SIZEa(SET)

 IF( N==0 )RETURN

 l = IntervalBisection(Entry,SET(1:N))

 IF( l(1)<1 )RETURN

 IF( SET(l(1))==Entry )THEN
  I = l(1)
 END IF

END IF

!!--end--
END FUNCTION

!!### SUBROUTINE: REMOVE_FROM_SET_Rdp
SUBROUTINE REMOVE_FROM_SET_Rdp(SET,Entry)
!!#### PURPOSE
!! Remove an entry from the SET.

!!#### REQUIRED OUTPUT
REAL(KIND_Rdp),POINTER    :: SET(:)

!!#### REQUIRED INPUT
REAL(KIND_Rdp),INTENT(IN) :: Entry

!!#### LOCAL VARIABLES
INTEGER :: N,I

!!--begin--
IF( ASSOCIATED(SET) )THEN
 I = FIND_IN_SET(Set,Entry)

 IF( I==0 )RETURN

 N = SIZE(SET)

 SET(I:N-1) = SET(I+1:N)

 SET(N) = ERROR(1._KIND_Rdp)

END IF

!!--end--
END SUBROUTINE

!!### SUBROUTINE: PRINT_SET_Rdp
SUBROUTINE PRINT_SET_Rdp( SET , Unit )
!!#### PURPOSE
!! Print out the contents of a set.

!!#### REQUIRED INPUT/OUTPUT
REAL(KIND_Rdp),POINTER :: SET(:)

!!#### OPTIONAL INPUT
INTEGER,OPTIONAL,INTENT(IN) :: Unit

!!#### LOCAL VARIABLES
INTEGER :: i,Unit_
CHARACTER(LEN_Rdp),POINTER :: D(:,:)

!!--begin--
Unit_ = DEFAULT(DEFAULT_OUTPUT_UNIT,Unit)

IF( .NOT.ASSOCIATED(SET) )THEN
 WRITE(Unit_,"(a)")" N U L L   S E T"
 RETURN
ELSE
 ALLOCATE( D(1+SIZE(SET),2) )
 CALL CLEAR(D)
 D(1,:) = (/"index","value"/)
 DO i=1,SIZE(SET)
  D(i+1,1) = TRIM(STR(i))
  D(i+1,2) = TRIM(STR(SET(i)))
 END DO
 CALL PRINT_Table(D,Unit=Unit_)
END IF

!!--end--
END SUBROUTINE

!!### SUBROUTINE: FINALIZE_SET_Rdp
SUBROUTINE FINALIZE_SET_Rdp( SET )
!!#### PURPOSE
!! Finalize the SET by allocation to its SIZEa length.

!!#### REQUIRED INPUT/OUTPUT
REAL(KIND_Rdp),POINTER :: SET(:)

!!#### LOCAL VARIABLES
INTEGER :: dn,N

!!--begin--
IF( ASSOCIATED(SET) )THEN
 N = SIZEa(SET)
 IF( N==0 )THEN
  DEALLOCATE(SET)
  NULLIFY(SET)
 ELSE
  dn=SIZEa(SET)-SIZE(SET)
  CALL Reallocate( SET , dn )
 END IF
ENDIF

!!--end--
END SUBROUTINE

END MODULE

max_stars_repo_name: wawiesel/tapack

max_issues_repo_name: wawiesel/tapack

max_forks_repo_name: wawiesel/tapack


text:
*> \brief \b SGER
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*  Definition:
*  ===========
*
*       SUBROUTINE SGER(M,N,ALPHA,X,INCX,Y,INCY,A,LDA)
*
*       .. Scalar Arguments ..
*       REAL ALPHA
*       INTEGER INCX,INCY,LDA,M,N
*       ..
*       .. Array Arguments ..
*       REAL A(LDA,*),X(*),Y(*)
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> SGER   performs the rank 1 operation
*>
*>    A := alpha*x*y**T + A,
*>
*> where alpha is a scalar, x is an m element vector, y is an n element
*> vector and A is an m by n matrix.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] M
*> \verbatim
*>          M is INTEGER
*>           On entry, M specifies the number of rows of the matrix A.
*>           M must be at least zero.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>           On entry, N specifies the number of columns of the matrix A.
*>           N must be at least zero.
*> \endverbatim
*>
*> \param[in] ALPHA
*> \verbatim
*>          ALPHA is REAL
*>           On entry, ALPHA specifies the scalar alpha.
*> \endverbatim
*>
*> \param[in] X
*> \verbatim
*>          X is REAL array of dimension at least
*>           ( 1 + ( m - 1 )*abs( INCX ) ).
*>           Before entry, the incremented array X must contain the m
*>           element vector x.
*> \endverbatim
*>
*> \param[in] INCX
*> \verbatim
*>          INCX is INTEGER
*>           On entry, INCX specifies the increment for the elements of
*>           X. INCX must not be zero.
*> \endverbatim
*>
*> \param[in] Y
*> \verbatim
*>          Y is REAL array of dimension at least
*>           ( 1 + ( n - 1 )*abs( INCY ) ).
*>           Before entry, the incremented array Y must contain the n
*>           element vector y.
*> \endverbatim
*>
*> \param[in] INCY
*> \verbatim
*>          INCY is INTEGER
*>           On entry, INCY specifies the increment for the elements of
*>           Y. INCY must not be zero.
*> \endverbatim
*>
*> \param[in,out] A
*> \verbatim
*>          A is REAL array of DIMENSION ( LDA, n ).
*>           Before entry, the leading m by n part of the array A must
*>           contain the matrix of coefficients. On exit, A is
*>           overwritten by the updated matrix.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>           On entry, LDA specifies the first dimension of A as declared
*>           in the calling (sub) program. LDA must be at least
*>           max( 1, m ).
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date December 2016
*
*> \ingroup single_blas_level2
*
*> \par Further Details:
*  =====================
*>
*> \verbatim
*>
*>  Level 2 Blas routine.
*>
*>  -- Written on 22-October-1986.
*>     Jack Dongarra, Argonne National Lab.
*>     Jeremy Du Croz, Nag Central Office.
*>     Sven Hammarling, Nag Central Office.
*>     Richard Hanson, Sandia National Labs.
*> \endverbatim
*>
*  =====================================================================
      SUBROUTINE SGER(M,N,ALPHA,X,INCX,Y,INCY,A,LDA)
*
*  -- Reference BLAS level2 routine (version 3.7.0) --
*  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     December 2016
*
*     .. Scalar Arguments ..
      REAL ALPHA
      INTEGER INCX,INCY,LDA,M,N
*     ..
*     .. Array Arguments ..
      REAL A(LDA,*),X(*),Y(*)
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL ZERO
      PARAMETER (ZERO=0.0E+0)
*     ..
*     .. Local Scalars ..
      REAL TEMP
      INTEGER I,INFO,IX,J,JY,KX
*     ..
*     .. External Subroutines ..
      EXTERNAL XERBLA
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC MAX
*     ..
*
*     Test the input parameters.
*
      INFO = 0
      IF (M.LT.0) THEN
          INFO = 1
      ELSE IF (N.LT.0) THEN
          INFO = 2
      ELSE IF (INCX.EQ.0) THEN
          INFO = 5
      ELSE IF (INCY.EQ.0) THEN
          INFO = 7
      ELSE IF (LDA.LT.MAX(1,M)) THEN
          INFO = 9
      END IF
      IF (INFO.NE.0) THEN
          CALL XERBLA('SGER  ',INFO)
          RETURN
      END IF
*
*     Quick return if possible.
*
      IF ((M.EQ.0) .OR. (N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN
*
*     Start the operations. In this version the elements of A are
*     accessed sequentially with one pass through A.
*
      IF (INCY.GT.0) THEN
          JY = 1
      ELSE
          JY = 1 - (N-1)*INCY
      END IF
      IF (INCX.EQ.1) THEN
          DO 20 J = 1,N
              IF (Y(JY).NE.ZERO) THEN
                  TEMP = ALPHA*Y(JY)
                  DO 10 I = 1,M
                      A(I,J) = A(I,J) + X(I)*TEMP
   10             CONTINUE
              END IF
              JY = JY + INCY
   20     CONTINUE
      ELSE
          IF (INCX.GT.0) THEN
              KX = 1
          ELSE
              KX = 1 - (M-1)*INCX
          END IF
          DO 40 J = 1,N
              IF (Y(JY).NE.ZERO) THEN
                  TEMP = ALPHA*Y(JY)
                  IX = KX
                  DO 30 I = 1,M
                      A(I,J) = A(I,J) + X(IX)*TEMP
                      IX = IX + INCX
   30             CONTINUE
              END IF
              JY = JY + INCY
   40     CONTINUE
      END IF
*
      RETURN
*
*     End of SGER  .
*
      END

max_stars_repo_name: lingochamp/OpenBLAS

max_issues_repo_name: lingochamp/OpenBLAS

max_forks_repo_name: lingochamp/OpenBLAS