!
!                Parallel Sparse BLAS  version 3.5
!      (C) Copyright 2006-2018
!        Salvatore Filippone
!        Alfredo Buttari
!
!    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. The name of the PSBLAS group or the names of its contributors may
!         not be used to endorse or promote products derived from this
!         software without specific 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 PSBLAS GROUP OR ITS 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.
!
!
! Author: Gibbs-Poole-Stockmeyer (revised by Stefano Toninel)
!
! Routines for Gibbs-Poole-Stockmeyer matrix bandwidth and profile reduction.
! Originally released in ACM-TOMS no. 508 written in Fortran77.
! Now revised and ported to Fortran90.
! Further revised and ported into the PSBLAS environment.
!
module psb_gps_mod
  !
  use psb_base_mod, only : psb_ipk_
  public psb_gps_reduce
  !
  ! COMMON /GRA/ N, IDPTH, IDEG
  !
  private
  ! common /CC/ XCC,SIZEG,STPT
  integer(psb_ipk_), save            :: xcc,n,idpth,ideg
  integer(psb_ipk_),allocatable,SAVE :: SIZEG(:),STPT(:)
  !
  ! COMMON /LVLW/ NHIGH,NLOW,NACUM
  integer(psb_ipk_),allocatable,target,save :: NHIGH(:),NLOW(:),NACUM(:),AUX(:)
  integer(psb_ipk_),PARAMETER :: INIT=500
  !
CONTAINS
  !
!!$  SUBROUTINE psb_gps_reduce(NDSTK, NR, IOLD, RENUM, NDEG, LVL, LVLS1, LVLS2,&
!!$       & CCSTOR, IBW2, IPF2,NE,IDPTHE,IDEGE)
  SUBROUTINE psb_gps_reduce(NDSTK, NR, IDEGE, IOLD, RENUM, NDEG,ibw2,ipf2,IDPTHE)
    !  SUBROUTINE REDUCE DETERMINES A ROW AND COLUMN PERMUTATION WHICH,
    !  WHEN APPLIED TO A GIVEN SPARSE MATRIX, PRODUCES A PERMUTED
    !  MATRIX WITH A SMALLER BANDWIDTH AND PROFILE.
    !  THE INPUT ARRAY IS A CONNECTION TABLE WHICH REPRESENTS THE
    !  INDICES OF THE NONZERO ELEMENTS OF THE MATRIX, A.  THE ALGO-
    !  RITHM IS DESCRIBED IN TERMS OF THE ADJACENCY GRAPH WHICH
    !  HAS THE CHARACTERISTIC THAT THERE IS AN EDGE (CONNECTION)
    !  BETWEEN NODES I AND J IF A(I,J) /= 0 AND I /= J.
    !  DIMENSIONING INFORMATION--THE FOLLOWING INTEGER ARRAYS MUST BE
    !  DIMENSIONED IN THE CALLING ROUTINE.
    !    NDSTK(NR,D1)        D1 IS >= MAXIMUM DEGREE OF ALL NODES.
    !    IOLD(D2)            D2 AND NR ARE >= THE TOTAL NUMBER OF
    !    RENUM(D2+1)         NODES IN THE GRAPH.
    !    NDEG(D2)            STORAGE REQUIREMENTS CAN BE SIGNIFICANTLY
    !    LVL(D2)             DECREASED FOR IBM 360 AND 370 COMPUTERS
    !    LVLS1(D2)           BY REPLACING INTEGER NDSTK BY
    !    LVLS2(D2)           INTEGER*2 NDSTK IN SUBROUTINES REDUCE,
    !    CCSTOR(D2)          DGREE, FNDIAM, TREE AND NUMBER.
    !  COMMON INFORMATION--THE FOLLOWING COMMON BLOCK MUST BE IN THE
    !  CALLING ROUTINE.
    !    COMMON/GRA/N,IDPTH,IDEG
    !  EXPLANATION OF INPUT VARIABLES--
    !    NDSTK-     CONNECTION TABLE REPRESENTING GRAPH.
    !               NDSTK(I,J)=NODE NUMBER OF JTH CONNECTION TO NODE
    !               NUMBER I.  A CONNECTION OF A NODE TO ITSELF IS NOT
    !               LISTED.  EXTRA POSITIONS MUST HAVE ZERO FILL.
    !    NR-        ROW DIMENSION ASSIGNED NDSTK IN CALLING PROGRAM.
    !    IOLD(I)-   NUMBERING OF ITH NODE UPON INPUT.
    !               IF NO NUMBERING EXISTS THEN IOLD(I)=I.
    !    N-         NUMBER OF NODES IN GRAPH (EQUAL TO ORDER OF MATRIX).
    !    IDEG-      MAXIMUM DEGREE OF ANY NODE IN THE GRAPH.
    !  EXPLANATION OF OUTPUT VARIABLES--
    !    RENUM(I)-  THE NEW NUMBER FOR THE ITH NODE.
    !    NDEG(I)-   THE DEGREE OF THE ITH NODE.
    !    IBW2-      THE BANDWIDTH AFTER RENUMBERING.
    !    IPF2-      THE PROFILE AFTER RENUMBERING.
    !    IDPTH-     NUMBER OF LEVELS IN REDUCE LEVEL STRUCTURE.
    !  THE FOLLOWING ONLY HAVE MEANING IF THE GRAPH WAS CONNECTED--
    !    LVL(I)-    INDEX INTO LVLS1 TO THE FIRST NODE IN LEVEL I.
    !               LVL(I+1)-LVL(I)= NUMBER OF NODES IN ITH LEVEL
    !    LVLS1-     NODE NUMBERS LISTED BY LEVEL.
    !    LVLS2(I)-  THE LEVEL ASSIGNED TO NODE I BY REDUCE.
    !  WORKING STORAGE VARIABLE--
    !    CCSTOR
    !  LOCAL STORAGE--
    !    COMMON/CC/-SUBROUTINES REDUCE, SORT2 AND PIKLVL ASSUME THAT
    !               THE GRAPH HAS AT MOST 50 CONNECTED COMPONENTS.
    !               SUBROUTINE FNDIAM ASSUMES THAT THERE ARE AT MOST
    !               100 NODES IN THE LAST LEVEL.
    !    COMMON/LVLW/-SUBROUTINES SETUP AND PIKLVL ASSUME THAT THERE
    !               ARE AT MOST 100 LEVELS.
    ! USE INTEGER*2 NDSTK  WITH AN IBM 360 OR 370.
    use psb_base_mod
    implicit none
    INTEGER(psb_ipk_) :: NR, IDEGE, IBW2, IPF2, IDPTHE
    ! COMMON /GRA/ N, IDPTH, IDEG
    ! IT IS ASSUMED THAT THE GRAPH HAS AT MOST 50 CONNECTED COMPONENTS.
    ! COMMON /CC/ XCC, SIZEG(50), STPT(50)
    ! COMMON /LVLW/ NHIGH(100), NLOW(100), NACUM(100)
    integer(psb_ipk_) :: stnode, rvnode, stnum, sbnum
    integer(psb_ipk_) :: ndstk(nr,idege), iold(nr), renum(nr+1), ndeg(nr)
    integer(psb_ipk_) :: lvl(nr), lvls1(nr), lvls2(nr), ccstor(nr)
    integer(psb_ipk_) :: nflg, info, i, ibw1, ipf1, idflt, isdir, lroot, lowdg
    integer(psb_ipk_) :: lvlbot, lvln, lvlwth, maxlw, num
    n     = nr
    ideg  = idege
    idpth = 0

    ALLOCATE(SIZEG(NR),STPT(NR), STAT=INFO)
    IF(INFO /= psb_success_) THEN
      write(psb_out_unit,*) 'ERROR! MEMORY ALLOCATION # 1 FAILED IN GPS'
      STOP
    END IF
    !
    ALLOCATE(NHIGH(INIT), NLOW(INIT), NACUM(INIT), AUX(INIT), STAT=INFO)
    IF(INFO /= psb_success_) THEN
      write(psb_out_unit,*) 'ERROR! MEMORY ALLOCATION # 2 FAILED IN GPS'
      STOP
    END IF
    !
    IBW2 = 0
    IPF2 = 0
    ! SET RENUM(I)=0 FOR ALL I TO INDICATE NODE I IS UNNUMBERED
    DO  I=1,N
      RENUM(I) = 0
    END DO
    ! COMPUTE DEGREE OF EACH NODE AND ORIGINAL BANDWIDTH AND PROFILE
    CALL DGREE(NDSTK, NR, NDEG, IOLD, IBW1, IPF1)
    ! SBNUM= LOW END OF AVAILABLE NUMBERS FOR RENUMBERING
    ! STNUM= HIGH END OF AVAILABLE NUMBERS FOR RENUMBERING
    SBNUM = 1
    STNUM = N
    ! NUMBER THE NODES OF DEGREE ZERO
    DO I=1,N
      IF (NDEG(I) > 0) CYCLE
      RENUM(I) = STNUM
      STNUM = STNUM - 1
    END DO
    ! FIND AN UNNUMBERED NODE OF MIN DEGREE TO START ON
    do
      LOWDG = IDEG + 1
      NFLG = 1
      ISDIR = 1
      DO  I=1,N
        IF (NDEG(I) >= LOWDG) CYCLE
        IF (RENUM(I) > 0) CYCLE
        LOWDG = NDEG(I)
        STNODE = I
      END DO
      ! FIND PSEUDO-DIAMETER AND ASSOCIATED LEVEL STRUCTURES.
      ! STNODE AND RVNODE ARE THE ENDS OF THE DIAM AND LVLS1 AND LVLS2
      ! ARE THE RESPECTIVE LEVEL STRUCTURES.
      CALL FNDIAM(STNODE, RVNODE, NDSTK, NR, NDEG, LVL, LVLS1,LVLS2, CCSTOR, IDFLT)
      IF (.not.(ndeg(stnode) <= ndeg(rvnode))) then
        ! NFLG INDICATES THE END TO BEGIN NUMBERING ON
        NFLG = -1
        STNODE = RVNODE
      endif
      CALL SETUP(LVL, LVLS1, LVLS2)
      ! FIND ALL THE CONNECTED COMPONENTS  (XCC COUNTS THEM)
      XCC = 0
      LROOT = 1
      LVLN = 1
      DO  I=1,N
        IF (LVL(I) /= 0) CYCLE
        XCC = XCC + 1
        STPT(XCC) = LROOT
        CALL TREE(I, NDSTK, NR, LVL, CCSTOR, NDEG, LVLWTH, LVLBOT,LVLN, MAXLW, N)
        SIZEG(XCC) = LVLBOT + LVLWTH - LROOT
        LROOT = LVLBOT + LVLWTH
        LVLN = LROOT
      END DO
      if (sort2() /= 0) then
        CALL PIKLVL(LVLS1, LVLS2, CCSTOR, IDFLT, ISDIR)
      endif
      ! ON RETURN FROM PIKLVL, ISDIR INDICATES THE DIRECTION THE LARGEST
      ! COMPONENT FELL.  ISDIR IS MODIFIED NOW TO INDICATE THE NUMBERING
      ! DIRECTION.  NUM IS SET TO THE PROPER VALUE FOR THIS DIRECTION.
      ISDIR = ISDIR*NFLG
      NUM = SBNUM
      IF (ISDIR < 0) NUM = STNUM
      CALL NUMBER(STNODE, NUM, NDSTK, LVLS2, NDEG, RENUM, LVLS1,LVL,&
           & NR, NFLG, IBW2, IPF2, CCSTOR, ISDIR)
      ! UPDATE STNUM OR SBNUM AFTER NUMBERING
      IF (ISDIR < 0) STNUM = NUM
      IF (ISDIR > 0) SBNUM = NUM
      IF (.not.(sbnum <= stnum)) exit
    end do
    IF (IBW2 > IBW1) then
      ! IF ORIGINAL NUMBERING IS BETTER THAN NEW ONE, SET UP TO RETURN IT
      DO I=1,N
        RENUM(I) = IOLD(I)
      END DO
      IBW2 = IBW1
      IPF2 = IPF1
      !
    endif
    DEALLOCATE(SIZEG,STPT,NHIGH,NLOW,AUX,NACUM)
    idpthe = idpth
    RETURN
  END SUBROUTINE PSB_GPS_REDUCE
  !
  SUBROUTINE DGREE(NDSTK, NR, NDEG, IOLD, IBW1, IPF1)
    !  DGREE COMPUTES THE DEGREE OF EACH NODE IN NDSTK AND STORES
    !  IT IN THE ARRAY NDEG.  THE BANDWIDTH AND PROFILE FOR THE ORIGINAL
    !  OR INPUT RENUMBERING OF THE GRAPH IS COMPUTED ALSO.
    ! USE INTEGER*2 NDSTK  WITH AN IBM 360 OR 370.
    implicit none
    INTEGER(psb_ipk_) :: NR, IBW1, IPF1, NDSTK(NR,IDEG), NDEG(N), IOLD(N)
    ! COMMON /GRA/ N, IDPTH, IDEG
    integer(psb_ipk_) :: i, itst, j, idif, irw

    IBW1 = 0
    IPF1 = 0
    DO  I=1,N
      NDEG(I) = 0
      IRW = 0
      DO J=1,IDEG
        ITST = NDSTK(I,J)
        IF(ITST <= 0) EXIT
        NDEG(I) = NDEG(I) + 1
        IDIF = IOLD(I) - IOLD(ITST)
        IF (IRW < IDIF) IRW = IDIF
      END DO
      IPF1 = IPF1 + IRW
      IF (IRW > IBW1) IBW1 = IRW
    END DO
    RETURN
  END SUBROUTINE DGREE
  !
  SUBROUTINE FNDIAM(SND1, SND2, NDSTK, NR, NDEG, LVL, LVLS1,LVLS2, IWK, IDFLT)
    !  FNDIAM IS THE CONTROL PROCEDURE FOR FINDING THE PSEUDO-DIAMETER OF
    !  NDSTK AS WELL AS THE LEVEL STRUCTURE FROM EACH END
    !  SND1-        ON INPUT THIS IS THE NODE NUMBER OF THE FIRST
    !               ATTEMPT AT FINDING A DIAMETER.  ON OUTPUT IT
    !               CONTAINS THE ACTUAL NUMBER USED.
    !  SND2-        ON OUTPUT CONTAINS OTHER END OF DIAMETER
    !  LVLS1-       ARRAY CONTAINING LEVEL STRUCTURE WITH SND1 AS ROOT
    !  LVLS2-       ARRAY CONTAINING LEVEL STRUCTURE WITH SND2 AS ROOT
    !  IDFLT-       FLAG USED IN PICKING FINAL LEVEL STRUCTURE, SET
    !               =1 IF WIDTH OF LVLS1 <= WIDTH OF LVLS2, OTHERWISE =2
    !  LVL,IWK-     WORKING STORAGE
    ! USE INTEGER*2 NDSTK  WITH AN IBM 360 OR 370.
    implicit none
    INTEGER(psb_ipk_) :: FLAG, SND, SND1, SND2, NR, idflt
    ! COMMON /GRA/ N, IDPTH, IDEG
    ! IT IS ASSUMED THAT THE LAST LEVEL HAS AT MOST 100 NODES.
    ! COMMON /CC/ NDLST(100)
    integer(psb_ipk_),POINTER :: NDLST(:)
    integer(psb_ipk_) :: NDSTK(NR,IDEG), NDEG(1), LVL(N), LVLS1(N), LVLS2(N),IWK(N)
    integer(psb_ipk_) :: i, j, mtw2, ndxn, ndxl, inow, lvlbot,lvln,lvlwth
    integer(psb_ipk_) :: k,mtw1, maxlw
    !
    NDLST => AUX
    !
    FLAG = 0
    MTW2 = N
    SND = SND1
    ! ZERO LVL TO INDICATE ALL NODES ARE AVAILABLE TO TREE
10  DO 20 I=1,N
       LVL(I) = 0
20  END DO
    LVLN = 1
    ! DROP A TREE FROM SND
    CALL TREE(SND, NDSTK, NR, LVL, IWK, NDEG, LVLWTH, LVLBOT,LVLN, MAXLW, MTW2)
    IF (FLAG >= 1) GO TO 50
    FLAG = 1
30  IDPTH = LVLN - 1
    MTW1 = MAXLW
    ! COPY LEVEL STRUCTURE INTO LVLS1
    DO 40 I=1,N
       LVLS1(I) = LVL(I)
40  END DO
    NDXN = 1
    NDXL = 0
    MTW2 = N
    ! SORT LAST LEVEL BY DEGREE  AND STORE IN NDLST
    CALL SORTDG(NDLST, IWK(LVLBOT), NDXL, LVLWTH, NDEG)
    SND = NDLST(1)
    GO TO 10
50  IF (IDPTH >= LVLN-1) GO TO 60
    ! START AGAIN WITH NEW STARTING NODE
    SND1 = SND
    GO TO 30
60  IF (MAXLW >= MTW2) GO TO 80
    MTW2 = MAXLW
    SND2 = SND
    ! STORE NARROWEST REVERSE LEVEL STRUCTURE IN LVLS2
    DO 70 I=1,N
       LVLS2(I) = LVL(I)
70  END DO
80  IF (NDXN == NDXL) GO TO 90
    ! TRY NEXT NODE IN NDLST
    NDXN = NDXN + 1
    SND = NDLST(NDXN)
    GO TO 10
90  IDFLT = 1
    IF (MTW2 <= MTW1) IDFLT = 2
    NULLIFY(NDLST)
    RETURN
  END SUBROUTINE FNDIAM
  !
  SUBROUTINE TREE(IROOT, NDSTK, NR, LVL, IWK, NDEG, LVLWTH, LVLBOT, LVLN, MAXLW, IBORT)
    !  TREE DROPS A TREE IN NDSTK FROM IROOT
    !  LVL-         ARRAY INDICATING AVAILABLE NODES IN NDSTK WITH ZERO
    !               ENTRIES. TREE ENTERS LEVEL NUMBERS ASSIGNED
    !               DURING EXECUTION OF THIS PROCEDURE
    !  IWK-         ON OUTPUT CONTAINS NODE NUMBERS USED IN TREE
    !               ARRANGED BY LEVELS (IWK(LVLN) CONTAINS IROOT
    !               AND IWK(LVLBOT+LVLWTH-1) CONTAINS LAST NODE ENTERED)
    !  LVLWTH-      ON OUTPUT CONTAINS WIDTH OF LAST LEVEL
    !  LVLBOT-      ON OUTPUT CONTAINS INDEX INTO IWK OF FIRST
    !               NODE IN LAST LEVEL
    !  MAXLW-       ON OUTPUT CONTAINS THE MAXIMUM LEVEL WIDTH
    !  LVLN-        ON INPUT THE FIRST AVAILABLE LOCATION IN IWK
    !               USUALLY ONE BUT IF IWK IS USED TO STORE PREVIOUS
    !               CONNECTED COMPONENTS, LVLN IS NEXT AVAILABLE LOCATION.
    !               ON OUTPUT THE TOTAL NUMBER OF LEVELS + 1
    !  IBORT-       INPUT PARAM WHICH TRIGGERS EARLY RETURN IF
    !               MAXLW BECOMES >= IBORT
    ! USE INTEGER*2 NDSTK  WITH AN IBM 360 OR 370.
    implicit none
    integer(psb_ipk_) :: IROOT, NR, NDSTK(NR,*), LVL(*), IWK(*), NDEG(*)
    integer(psb_ipk_) :: LVLWTH, LVLBOT, LVLN, MAXLW, IBORT
    integer(psb_ipk_) :: itest, iwknow, itop, lvltop,j , inow, ndrow
    MAXLW = 0
    ITOP = LVLN
    INOW = LVLN
    LVLBOT = LVLN
    LVLTOP = LVLN + 1
    LVLN = 1
    LVL(IROOT) = 1
    IWK(ITOP) = IROOT
10  LVLN = LVLN + 1
20  IWKNOW = IWK(INOW)
    NDROW = NDEG(IWKNOW)
    DO 30 J=1,NDROW
       ITEST = NDSTK(IWKNOW,J)
       IF (LVL(ITEST) /= 0) CYCLE
       LVL(ITEST) = LVLN
       ITOP = ITOP + 1
       IWK(ITOP) = ITEST
30  END DO
    INOW = INOW + 1
    IF (INOW < LVLTOP) GO TO 20
    LVLWTH = LVLTOP - LVLBOT
    IF (MAXLW < LVLWTH) MAXLW = LVLWTH
    IF (MAXLW >= IBORT) RETURN
    IF (ITOP < LVLTOP) RETURN
    LVLBOT = INOW
    LVLTOP = ITOP + 1
    GO TO 10
  END SUBROUTINE TREE
  !
  SUBROUTINE SORTDG(STK1, STK2, X1, X2, NDEG)
    ! SORTDG SORTS STK2 BY DEGREE OF THE NODE AND ADDS IT TO THE END
    ! OF STK1 IN ORDER OF LOWEST TO HIGHEST DEGREE.  X1 AND X2 ARE THE
    ! NUMBER OF NODES IN STK1 AND STK2 RESPECTIVELY.
    implicit none
    INTEGER(psb_ipk_) :: X1, X2, STK1, STK2, TEMP,NDEG
    ! COMMON /GRA/ N, IDPTH, IDEG
    DIMENSION NDEG(N), STK1(X1+X2), STK2(X2)
    integer(psb_ipk_) :: ind,itest,i,j,istk2,jstk2
    IND = X2
10  ITEST = 0
    IND = IND - 1
    IF (IND < 1) GO TO 30
    DO 20 I=1,IND
       J = I + 1
       ISTK2 = STK2(I)
       JSTK2 = STK2(J)
       IF (NDEG(ISTK2) <= NDEG(JSTK2)) CYCLE
       ITEST = 1
       TEMP = STK2(I)
       STK2(I) = STK2(J)
       STK2(J) = TEMP
20  END DO
    IF (ITEST == 1) GO TO 10
30  DO 40 I=1,X2
       X1 = X1 + 1
       STK1(X1) = STK2(I)
40  END DO
    RETURN
  END SUBROUTINE SORTDG
  !
  SUBROUTINE SETUP(LVL, LVLS1, LVLS2)
    ! SETUP COMPUTES THE REVERSE LEVELING INFO FROM LVLS2 AND STORES
    ! IT INTO LVLS2.  NACUM(I) IS INITIALIZED TO NODES/ITH LEVEL FOR NODES
    ! ON THE PSEUDO-DIAMETER OF THE GRAPH.  LVL IS INITIALIZED TO NON-
    ! ZERO FOR NODES ON THE PSEUDO-DIAM AND NODES IN A DIFFERENT
    ! COMPONENT OF THE GRAPH.
    !      COMMON /GRA/ N, IDPTH, IDEG
    ! IT IS ASSUMED THAT THERE ARE AT MOST 100 LEVELS.
    ! COMMON /LVLW/ NHIGH(100), NLOW(100), NACUM(100)
    use psb_base_mod
    implicit none
    integer(psb_ipk_) :: LVL(N), LVLS1(N), LVLS2(N)
    integer(psb_ipk_) :: SZ,i,itemp
    !-----------------------------------------------------
    SZ=SIZE(NACUM)
    IF(SZ < IDPTH) THEN
       write(psb_out_unit,*) 'GPS_SETUP: on fly reallocation of NACUM'
       CALL REALLOC(NACUM,SZ,IDPTH)
    END IF
    !-----------------------------------------------------
    DO 10 I=1,IDPTH
       NACUM(I) = 0
10  END DO
    DO 30 I=1,N
       LVL(I) = 1
       LVLS2(I) = IDPTH + 1 - LVLS2(I)
       ITEMP = LVLS2(I)
       IF (ITEMP > IDPTH) CYCLE
       IF (ITEMP /= LVLS1(I)) GO TO 20
       NACUM(ITEMP) = NACUM(ITEMP) + 1
       CYCLE
20     LVL(I) = 0
30  END DO
    RETURN
  END SUBROUTINE SETUP
  !
  FUNCTION SORT2() result(val)
    implicit none
    INTEGER(psb_ipk_) ::  val
    ! SORT2 SORTS SIZEG AND STPT INTO DESCENDING ORDER ACCORDING TO
    ! VALUES OF SIZEG. XCC=NUMBER OF ENTRIES IN EACH ARRAY
    INTEGER(psb_ipk_) :: TEMP,ind,itest,i,j
    ! IT IS ASSUMED THAT THE GRAPH HAS AT MOST 50 CONNECTED COMPONENTS.
    !COMMON /CC/ XCC, SIZEG(50), STPT(50)

    VAL = 0
    IF (XCC == 0) RETURN
    VAL = 1
    IND = XCC
10  ITEST = 0
    IND = IND - 1
    IF (IND < 1) RETURN
    DO 20 I=1,IND
       J = I + 1
       IF (SIZEG(I) >= SIZEG(J)) CYCLE
       ITEST = 1
       TEMP = SIZEG(I)
       SIZEG(I) = SIZEG(J)
       SIZEG(J) = TEMP
       TEMP = STPT(I)
       STPT(I) = STPT(J)
       STPT(J) = TEMP
20  END DO
    IF (ITEST == 1) GO TO 10
    RETURN
  END FUNCTION SORT2
  !
  SUBROUTINE PIKLVL(LVLS1, LVLS2, CCSTOR, IDFLT, ISDIR)
    use psb_base_mod
    implicit none
    ! PIKLVL CHOOSES THE LEVEL STRUCTURE  USED IN NUMBERING GRAPH
    ! LVLS1-       ON INPUT CONTAINS FORWARD LEVELING INFO
    ! LVLS2-       ON INPUT CONTAINS REVERSE LEVELING INFO
    !              ON OUTPUT THE FINAL LEVEL STRUCTURE CHOSEN
    ! CCSTOR-      ON INPUT CONTAINS CONNECTED COMPONENT INFO
    ! IDFLT-       ON INPUT =1 IF WDTH LVLS1 <= WDTH LVLS2, =2 OTHERWISE
    ! NHIGH        KEEPS TRACK OF LEVEL WIDTHS FOR HIGH NUMBERING
    ! NLOW-        KEEPS TRACK OF LEVEL WIDTHS FOR LOW NUMBERING
    ! NACUM-       KEEPS TRACK OF LEVEL WIDTHS FOR CHOSEN LEVEL STRUCTURE
    ! XCC-         NUMBER OF CONNECTED COMPONENTS
    ! SIZEG(I)-    SIZE OF ITH CONNECTED COMPONENT
    ! STPT(I)-     INDEX INTO CCSTORE OF 1ST NODE IN ITH CON COMPT
    ! ISDIR-       FLAG WHICH INDICATES WHICH WAY THE LARGEST CONNECTED
    !              COMPONENT FELL.  =+1 IF LOW AND -1 IF HIGH
    ! COMMON /GRA/ N, IDPTH, IDEG
    ! IT IS ASSUMED THAT THE GRAPH HAS AT MOST 50 COMPONENTS AND
    ! THAT THERE ARE AT MOST 100 LEVELS.
    ! COMMON /LVLW/ NHIGH(100), NLOW(100), NACUM(100)
    ! COMMON /CC/ XCC, SIZEG(50), STPT(50)
    INTEGER(psb_ipk_) :: LVLS1(N), LVLS2(N), CCSTOR(N)
    integer(psb_ipk_) :: SZ, ENDC,i,j,max1,max2,inode
    integer(psb_ipk_) :: lvlnh, it, k, lvlnl,idflt,isdir
    ! FOR EACH CONNECTED COMPONENT DO
    DO 80 I=1,XCC
       J = STPT(I)
       ENDC= SIZEG(I) + J - 1
       ! SET NHIGH AND NLOW EQUAL TO NACUM
       !-----------------------------------------------------
       SZ=SIZE(NHIGH)
       IF(SZ < IDPTH) THEN
          write(psb_out_unit,*) 'GPS_PIKLVL: on fly reallocation of NHIGH'
          CALL REALLOC(NHIGH,SZ,IDPTH)
       END IF
       SZ=SIZE(NLOW)
       IF(SZ < IDPTH) THEN
          write(psb_out_unit,*) 'GPS_PIKLVL: on fly reallocation of NLOW'
          CALL REALLOC(NLOW,SZ,IDPTH)
       END IF
       !-----------------------------------------------------
       DO 10 K=1,IDPTH
          NHIGH(K) = NACUM(K)
          NLOW(K) = NACUM(K)
10     END DO
       ! UPDATE NHIGH AND NLOW FOR EACH NODE IN CONNECTED COMPONENT
       DO 20 K=J,ENDC
          INODE = CCSTOR(K)
          LVLNH = LVLS1(INODE)
          NHIGH(LVLNH) = NHIGH(LVLNH) + 1
          LVLNL = LVLS2(INODE)
          NLOW(LVLNL) = NLOW(LVLNL) + 1
20     END DO
       MAX1 = 0
       MAX2 = 0
       ! SET MAX1=LARGEST NEW NUMBER IN NHIGH
       ! SET MAX2=LARGEST NEW NUMBER IN NLOW
       DO 30 K=1,IDPTH
          IF (2*NACUM(K) == NLOW(K)+NHIGH(K)) CYCLE
          IF (NHIGH(K) > MAX1) MAX1 = NHIGH(K)
          IF (NLOW(K) > MAX2) MAX2 = NLOW(K)
30     END DO
       ! SET IT= NUMBER OF LEVEL STRUCTURE TO BE USED
       IT = 1
       IF (MAX1 > MAX2) IT = 2
       IF (MAX1 == MAX2) IT = IDFLT
       IF (IT == 2) GO TO 60
       IF (I == 1) ISDIR = -1
       ! COPY LVLS1 INTO LVLS2 FOR EACH NODE IN CONNECTED COMPONENT
       DO 40 K=J,ENDC
          INODE = CCSTOR(K)
          LVLS2(INODE) = LVLS1(INODE)
40     END DO
       ! UPDATE NACUM TO BE THE SAME AS NHIGH
       DO 50 K=1,IDPTH
          NACUM(K) = NHIGH(K)
50     END DO
       CYCLE
       ! UPDATE NACUM TO BE THE SAME AS NLOW
60     DO 70 K=1,IDPTH
          NACUM(K) = NLOW(K)
70     END DO
80  END DO
    RETURN
  END SUBROUTINE PIKLVL
  !
  SUBROUTINE NUMBER(SND, NUM, NDSTK, LVLS2, NDEG, RENUM, LVLST,LSTPT,&
       & NR, NFLG, IBW2, IPF2, IPFA, ISDIR)
    use psb_base_mod
    implicit none
    !  NUMBER PRODUCES THE NUMBERING OF THE GRAPH FOR MIN BANDWIDTH
    !  SND-         ON INPUT THE NODE TO BEGIN NUMBERING ON
    !  NUM-         ON INPUT AND OUTPUT, THE NEXT AVAILABLE NUMBER
    !  LVLS2-       THE LEVEL STRUCTURE TO BE USED IN NUMBERING
    !  RENUM-       THE ARRAY USED TO STORE THE NEW NUMBERING
    !  LVLST-       ON OUTPUT CONTAINS LEVEL STRUCTURE
    !  LSTPT(I)-    ON OUTPUT, INDEX INTO LVLST TO FIRST NODE IN ITH LVL
    !               LSTPT(I+1) - LSTPT(I) = NUMBER OF NODES IN ITH LVL
    !  NFLG-        =+1 IF SND IS FORWARD END OF PSEUDO-DIAM
    !               =-1 IF SND IS REVERSE END OF PSEUDO-DIAM
    !  IBW2-        BANDWIDTH OF NEW NUMBERING COMPUTED BY NUMBER
    !  IPF2-        PROFILE OF NEW NUMBERING COMPUTED BY NUMBER
    !  IPFA-        WORKING STORAGE USED TO COMPUTE PROFILE AND BANDWIDTH
    !  ISDIR-       INDICATES STEP DIRECTION USED IN NUMBERING(+1 OR -1)
    ! USE INTEGER*2 NDSTK  WITH AN IBM 360 OR 370.
    INTEGER(psb_ipk_) :: SND, NUM, XA, XB, XC, XD, CX, ENDC, TEST, NR, ISDIR
    ! COMMON /GRA/ N, IDPTH, IDEG
    ! THE STORAGE IN COMMON BLOCKS CC AND LVLW IS NOW FREE AND CAN
    ! BE USED FOR STACKS.
    !COMMON /LVLW/ STKA(100), STKB(100), STKC(100)
    !COMMON /CC/ STKD(100)
    INTEGER(psb_ipk_) :: IPFA(N), NDSTK(NR,IDEG), LVLS2(N),&
         & NDEG(N), RENUM(N+1), LVLST(N),LSTPT(N),ipf2,ibw2,nflg, nbw
    integer(psb_ipk_),POINTER :: STKA(:),STKB(:),STKC(:),STKD(:)
    integer(psb_ipk_) :: SZ1,SZ2,i,j,nstpt,lvln, lst, lnd, inx, max, ipro,&
         & lvlnl, k, it
    !
    STKA => NHIGH
    STKB => NLOW
    STKC => NACUM
    STKD => AUX
    !
    ! SET UP LVLST AND LSTPT FROM LVLS2
    DO 10 I=1,N
       IPFA(I) = 0
10  END DO
    NSTPT = 1
    DO 30 I=1,IDPTH
       LSTPT(I) = NSTPT
       DO 20 J=1,N
          IF (LVLS2(J) /= I) CYCLE
          LVLST(NSTPT) = J
          NSTPT = NSTPT + 1
20     END DO
30  END DO
    LSTPT(IDPTH+1) = NSTPT
    ! STKA, STKB, STKC AND STKD ARE STACKS WITH POINTERS
    ! XA,XB,XC, AND XD.  CX IS A SPECIAL POINTER INTO STKC WHICH
    ! INDICATES THE PARTICULAR NODE BEING PROCESSED.
    ! LVLN KEEPS TRACK OF THE LEVEL WE ARE WORKING AT.
    ! INITIALLY STKC CONTAINS ONLY THE INITIAL NODE, SND.
    LVLN = 0
    IF (NFLG < 0) LVLN = IDPTH + 1
    XC = 1
    STKC(XC) = SND
40  CX = 1
    XD = 0
    LVLN = LVLN + NFLG
    LST = LSTPT(LVLN)
    LND = LSTPT(LVLN+1) - 1
    ! BEGIN PROCESSING NODE STKC(CX)
50  IPRO = STKC(CX)
    RENUM(IPRO) = NUM
    NUM = NUM + ISDIR
    ENDC = NDEG(IPRO)
    XA = 0
    XB = 0
    ! CHECK ALL ADJACENT NODES
    DO 80 I=1,ENDC
       TEST = NDSTK(IPRO,I)
       INX = RENUM(TEST)
       ! ONLY NODES NOT NUMBERED OR ALREADY ON A STACK ARE ADDED
       IF (INX == 0) GO TO 60
       IF (INX < 0) CYCLE
       ! DO PRELIMINARY BANDWIDTH AND PROFILE CALCULATIONS
       NBW = (RENUM(IPRO)-INX)*ISDIR
       IF (ISDIR > 0) INX = RENUM(IPRO)
       IF (IPFA(INX) < NBW) IPFA(INX) = NBW
       CYCLE
60     RENUM(TEST) = -1
       ! PUT NODES ON SAME LEVEL ON STKA, ALL OTHERS ON STKB
       IF (LVLS2(TEST) == LVLS2(IPRO)) GO TO 70
       XB = XB + 1
       STKB(XB) = TEST
       CYCLE
70     XA = XA + 1
       STKA(XA) = TEST
80  END DO
    ! SORT STKA AND STKB INTO INCREASING DEGREE AND ADD STKA TO STKC
    ! AND STKB TO STKD
    IF (XA == 0) GO TO 100
    IF (XA == 1) GO TO 90
    !-----------------------------------------------------------------
    SZ1=SIZE(STKC)
    SZ2=XC+XA
    IF(SZ1 < SZ2) THEN
       write(psb_out_unit,*) 'GPS_NUMBER - Check #1: on fly reallocation of STKC'
       CALL REALLOC(NACUM,SZ1,SZ2)
       STKC => NACUM
    END IF
    !-----------------------------------------------------------------
    CALL SORTDG(STKC, STKA, XC, XA, NDEG)
    GO TO 100
90  XC = XC + 1
    !-----------------------------------------------------------------
    SZ1=SIZE(STKC)
    SZ2=XC
    IF(SZ1 < SZ2) THEN
       write(psb_out_unit,*) 'GPS_NUMBER - Check #2: on fly reallocation of STKC'
       SZ2=SZ2+INIT
       CALL REALLOC(NACUM,SZ1,SZ2)
       STKC => NACUM
    END IF
    !-----------------------------------------------------------------
    STKC(XC) = STKA(XA)
100 IF (XB == 0) GO TO 120
    IF (XB == 1) GO TO 110
    !-----------------------------------------------------------------
    SZ1=SIZE(STKD)
    SZ2=XD+XB
    IF(SZ1 < SZ2) THEN
       write(psb_out_unit,*) 'GPS_NUMBER - Check #3: on fly reallocation of STKD'
       CALL REALLOC(AUX,SZ1,SZ2)
       STKD => AUX
    END IF
    !-----------------------------------------------------------------
    CALL SORTDG(STKD, STKB, XD, XB, NDEG)
    GO TO 120
110 XD = XD + 1
    !-----------------------------------------------------------------
    SZ1=SIZE(STKD)
    SZ2=XD
    IF(SZ1 < SZ2) THEN
       write(psb_out_unit,*) 'GPS_NUMBER - Check #4: on fly reallocation of STKD'
       SZ2=SZ2+INIT
       CALL REALLOC(AUX,SZ1,SZ2)
       STKD => AUX
    END IF
    !-----------------------------------------------------------------
    STKD(XD) = STKB(XB)
    ! BE SURE TO PROCESS ALL NODES IN STKC
120 CX = CX + 1
    IF (XC >= CX) GO TO 50
    ! WHEN STKC IS EXHAUSTED LOOK FOR MIN DEGREE NODE IN SAME LEVEL
    ! WHICH HAS NOT BEEN PROCESSED
    MAX = IDEG + 1
    SND = N + 1
    DO 130 I=LST,LND
       TEST = LVLST(I)
       IF (RENUM(TEST) /= 0) CYCLE
       IF (NDEG(TEST) >= MAX) CYCLE
       RENUM(SND) = 0
       RENUM(TEST) = -1
       MAX = NDEG(TEST)
       SND = TEST
130 END DO
    IF (SND == N+1) GO TO 140
    XC = XC + 1
    !-----------------------------------------------------------------
    SZ1=SIZE(STKC)
    SZ2=XC
    IF(SZ1 < SZ2) THEN
       write(psb_out_unit,*) 'GPS_NUMBER - Check #5: on fly reallocation of STKC'
       SZ2=SZ2+INIT
       CALL REALLOC(NACUM,SZ1,SZ2)
       STKC => NACUM
    END IF
    !-----------------------------------------------------------------
    STKC(XC) = SND
    GO TO 50
    ! IF STKD IS EMPTY WE ARE DONE, OTHERWISE COPY STKD ONTO STKC
    ! AND BEGIN PROCESSING NEW STKC
140 IF (XD == 0) GO TO 160
    !-----------------------------------------------------------------
    SZ1=SIZE(STKC)
    SZ2=XD
    IF(SZ1 < SZ2) THEN
       write(psb_out_unit,*) 'GPS_NUMBER - Check #6: on fly reallocation of STKC'
       SZ2=SZ2+INIT
       CALL REALLOC(NACUM,SZ1,SZ2)
       STKC => NACUM
    END IF
    !-----------------------------------------------------------------
    DO 150 I=1,XD
       STKC(I) = STKD(I)
150 END DO
    XC = XD
    GO TO 40
    ! DO FINAL BANDWIDTH AND PROFILE CALCULATIONS
160 DO 170 I=1,N
       IF (IPFA(I) > IBW2) IBW2 = IPFA(I)
       IPF2 = IPF2 + IPFA(I)
170 END DO
    !
    RETURN
  END SUBROUTINE NUMBER
  !
  ! ---------------------------------------------------------
  SUBROUTINE REALLOC(VET,SZ1,SZ2)
    use psb_base_mod
    ! PERFORM ON FLY REALLOCATION OF POINTER VET INCREASING
    ! ITS SIZE FROM SZ1 TO SZ2
    IMPLICIT NONE
    integer(psb_ipk_),allocatable  :: VET(:)
    integer(psb_ipk_) :: SZ1,SZ2,INFO

    call psb_realloc(sz2,vet,info)
    IF(INFO /= psb_success_) THEN
       write(psb_out_unit,*) 'Error! Memory allocation failure in REALLOC'
       STOP
    END IF
    RETURN
  END SUBROUTINE REALLOC
  !
end module psb_gps_mod

module psb_lgps_mod
  !
  use psb_base_mod, only : psb_lpk_
  public psb_lgps_reduce
  !
  ! COMMON /GRA/ N, IDPTH, IDEG
  !
  private
  ! common /CC/ XCC,SIZEG,STPT
  integer(psb_lpk_), save            :: xcc,n,idpth,ideg
  integer(psb_lpk_),allocatable,SAVE :: SIZEG(:),STPT(:)
  !
  ! COMMON /LVLW/ NHIGH,NLOW,NACUM
  integer(psb_lpk_),allocatable,target,save :: NHIGH(:),NLOW(:),NACUM(:),AUX(:)
  integer(psb_lpk_),PARAMETER :: INIT=500
  !
CONTAINS
  !
!!$  SUBROUTINE psb_gps_reduce(NDSTK, NR, IOLD, RENUM, NDEG, LVL, LVLS1, LVLS2,&
!!$       & CCSTOR, IBW2, IPF2,NE,IDPTHE,IDEGE)
  SUBROUTINE psb_lgps_reduce(NDSTK, NR, IDEGE, IOLD, RENUM, NDEG,ibw2,ipf2,IDPTHE)
    !  SUBROUTINE REDUCE DETERMINES A ROW AND COLUMN PERMUTATION WHICH,
    !  WHEN APPLIED TO A GIVEN SPARSE MATRIX, PRODUCES A PERMUTED
    !  MATRIX WITH A SMALLER BANDWIDTH AND PROFILE.
    !  THE INPUT ARRAY IS A CONNECTION TABLE WHICH REPRESENTS THE
    !  INDICES OF THE NONZERO ELEMENTS OF THE MATRIX, A.  THE ALGO-
    !  RITHM IS DESCRIBED IN TERMS OF THE ADJACENCY GRAPH WHICH
    !  HAS THE CHARACTERISTIC THAT THERE IS AN EDGE (CONNECTION)
    !  BETWEEN NODES I AND J IF A(I,J) /= 0 AND I /= J.
    !  DIMENSIONING INFORMATION--THE FOLLOWING INTEGER ARRAYS MUST BE
    !  DIMENSIONED IN THE CALLING ROUTINE.
    !    NDSTK(NR,D1)        D1 IS >= MAXIMUM DEGREE OF ALL NODES.
    !    IOLD(D2)            D2 AND NR ARE >= THE TOTAL NUMBER OF
    !    RENUM(D2+1)         NODES IN THE GRAPH.
    !    NDEG(D2)            STORAGE REQUIREMENTS CAN BE SIGNIFICANTLY
    !    LVL(D2)             DECREASED FOR IBM 360 AND 370 COMPUTERS
    !    LVLS1(D2)           BY REPLACING INTEGER NDSTK BY
    !    LVLS2(D2)           INTEGER*2 NDSTK IN SUBROUTINES REDUCE,
    !    CCSTOR(D2)          DGREE, FNDIAM, TREE AND NUMBER.
    !  COMMON INFORMATION--THE FOLLOWING COMMON BLOCK MUST BE IN THE
    !  CALLING ROUTINE.
    !    COMMON/GRA/N,IDPTH,IDEG
    !  EXPLANATION OF INPUT VARIABLES--
    !    NDSTK-     CONNECTION TABLE REPRESENTING GRAPH.
    !               NDSTK(I,J)=NODE NUMBER OF JTH CONNECTION TO NODE
    !               NUMBER I.  A CONNECTION OF A NODE TO ITSELF IS NOT
    !               LISTED.  EXTRA POSITIONS MUST HAVE ZERO FILL.
    !    NR-        ROW DIMENSION ASSIGNED NDSTK IN CALLING PROGRAM.
    !    IOLD(I)-   NUMBERING OF ITH NODE UPON INPUT.
    !               IF NO NUMBERING EXISTS THEN IOLD(I)=I.
    !    N-         NUMBER OF NODES IN GRAPH (EQUAL TO ORDER OF MATRIX).
    !    IDEG-      MAXIMUM DEGREE OF ANY NODE IN THE GRAPH.
    !  EXPLANATION OF OUTPUT VARIABLES--
    !    RENUM(I)-  THE NEW NUMBER FOR THE ITH NODE.
    !    NDEG(I)-   THE DEGREE OF THE ITH NODE.
    !    IBW2-      THE BANDWIDTH AFTER RENUMBERING.
    !    IPF2-      THE PROFILE AFTER RENUMBERING.
    !    IDPTH-     NUMBER OF LEVELS IN REDUCE LEVEL STRUCTURE.
    !  THE FOLLOWING ONLY HAVE MEANING IF THE GRAPH WAS CONNECTED--
    !    LVL(I)-    INDEX INTO LVLS1 TO THE FIRST NODE IN LEVEL I.
    !               LVL(I+1)-LVL(I)= NUMBER OF NODES IN ITH LEVEL
    !    LVLS1-     NODE NUMBERS LISTED BY LEVEL.
    !    LVLS2(I)-  THE LEVEL ASSIGNED TO NODE I BY REDUCE.
    !  WORKING STORAGE VARIABLE--
    !    CCSTOR
    !  LOCAL STORAGE--
    !    COMMON/CC/-SUBROUTINES REDUCE, SORT2 AND PIKLVL ASSUME THAT
    !               THE GRAPH HAS AT MOST 50 CONNECTED COMPONENTS.
    !               SUBROUTINE FNDIAM ASSUMES THAT THERE ARE AT MOST
    !               100 NODES IN THE LAST LEVEL.
    !    COMMON/LVLW/-SUBROUTINES SETUP AND PIKLVL ASSUME THAT THERE
    !               ARE AT MOST 100 LEVELS.
    ! USE INTEGER*2 NDSTK  WITH AN IBM 360 OR 370.
    use psb_base_mod
    implicit none
    INTEGER(psb_lpk_) :: NR, IDEGE, IBW2, IPF2, IDPTHE
    ! COMMON /GRA/ N, IDPTH, IDEG
    ! IT IS ASSUMED THAT THE GRAPH HAS AT MOST 50 CONNECTED COMPONENTS.
    ! COMMON /CC/ XCC, SIZEG(50), STPT(50)
    ! COMMON /LVLW/ NHIGH(100), NLOW(100), NACUM(100)
    integer(psb_lpk_) :: stnode, rvnode, stnum, sbnum
    integer(psb_lpk_) :: ndstk(nr,idege), iold(nr), renum(nr+1), ndeg(nr)
    integer(psb_lpk_) :: lvl(nr), lvls1(nr), lvls2(nr), ccstor(nr)
    integer(psb_lpk_) :: nflg, info, i, ibw1, ipf1, idflt, isdir, lroot, lowdg
    integer(psb_lpk_) :: lvlbot, lvln, lvlwth, maxlw, num
    n     = nr
    ideg  = idege
    idpth = 0

    ALLOCATE(SIZEG(NR),STPT(NR), STAT=INFO)
    IF(INFO /= psb_success_) THEN
      write(psb_out_unit,*) 'ERROR! MEMORY ALLOCATION # 1 FAILED IN GPS'
      STOP
    END IF
    !
    ALLOCATE(NHIGH(INIT), NLOW(INIT), NACUM(INIT), AUX(INIT), STAT=INFO)
    IF(INFO /= psb_success_) THEN
      write(psb_out_unit,*) 'ERROR! MEMORY ALLOCATION # 2 FAILED IN GPS'
      STOP
    END IF
    !
    IBW2 = 0
    IPF2 = 0
    ! SET RENUM(I)=0 FOR ALL I TO INDICATE NODE I IS UNNUMBERED
    DO  I=1,N
      RENUM(I) = 0
    END DO
    ! COMPUTE DEGREE OF EACH NODE AND ORIGINAL BANDWIDTH AND PROFILE
    CALL DGREE(NDSTK, NR, NDEG, IOLD, IBW1, IPF1)
    ! SBNUM= LOW END OF AVAILABLE NUMBERS FOR RENUMBERING
    ! STNUM= HIGH END OF AVAILABLE NUMBERS FOR RENUMBERING
    SBNUM = 1
    STNUM = N
    ! NUMBER THE NODES OF DEGREE ZERO
    DO I=1,N
      IF (NDEG(I) > 0) CYCLE
      RENUM(I) = STNUM
      STNUM = STNUM - 1
    END DO
    ! FIND AN UNNUMBERED NODE OF MIN DEGREE TO START ON
    do
      LOWDG = IDEG + 1
      NFLG = 1
      ISDIR = 1
      DO  I=1,N
        IF (NDEG(I) >= LOWDG) CYCLE
        IF (RENUM(I) > 0) CYCLE
        LOWDG = NDEG(I)
        STNODE = I
      END DO
      ! FIND PSEUDO-DIAMETER AND ASSOCIATED LEVEL STRUCTURES.
      ! STNODE AND RVNODE ARE THE ENDS OF THE DIAM AND LVLS1 AND LVLS2
      ! ARE THE RESPECTIVE LEVEL STRUCTURES.
      CALL FNDIAM(STNODE, RVNODE, NDSTK, NR, NDEG, LVL, LVLS1,LVLS2, CCSTOR, IDFLT)
      IF (.not.(ndeg(stnode) <= ndeg(rvnode))) then
        ! NFLG INDICATES THE END TO BEGIN NUMBERING ON
        NFLG = -1
        STNODE = RVNODE
      endif
      CALL SETUP(LVL, LVLS1, LVLS2)
      ! FIND ALL THE CONNECTED COMPONENTS  (XCC COUNTS THEM)
      XCC = 0
      LROOT = 1
      LVLN = 1
      DO  I=1,N
        IF (LVL(I) /= 0) CYCLE
        XCC = XCC + 1
        STPT(XCC) = LROOT
        CALL TREE(I, NDSTK, NR, LVL, CCSTOR, NDEG, LVLWTH, LVLBOT,LVLN, MAXLW, N)
        SIZEG(XCC) = LVLBOT + LVLWTH - LROOT
        LROOT = LVLBOT + LVLWTH
        LVLN = LROOT
      END DO
      if (sort2() /= 0) then
        CALL PIKLVL(LVLS1, LVLS2, CCSTOR, IDFLT, ISDIR)
      endif
      ! ON RETURN FROM PIKLVL, ISDIR INDICATES THE DIRECTION THE LARGEST
      ! COMPONENT FELL.  ISDIR IS MODIFIED NOW TO INDICATE THE NUMBERING
      ! DIRECTION.  NUM IS SET TO THE PROPER VALUE FOR THIS DIRECTION.
      ISDIR = ISDIR*NFLG
      NUM = SBNUM
      IF (ISDIR < 0) NUM = STNUM
      CALL NUMBER(STNODE, NUM, NDSTK, LVLS2, NDEG, RENUM, LVLS1,LVL,&
           & NR, NFLG, IBW2, IPF2, CCSTOR, ISDIR)
      ! UPDATE STNUM OR SBNUM AFTER NUMBERING
      IF (ISDIR < 0) STNUM = NUM
      IF (ISDIR > 0) SBNUM = NUM
      IF (.not.(sbnum <= stnum)) exit
    end do
    IF (IBW2 > IBW1) then
      ! IF ORIGINAL NUMBERING IS BETTER THAN NEW ONE, SET UP TO RETURN IT
      DO I=1,N
        RENUM(I) = IOLD(I)
      END DO
      IBW2 = IBW1
      IPF2 = IPF1
      !
    endif
    DEALLOCATE(SIZEG,STPT,NHIGH,NLOW,AUX,NACUM)
    idpthe = idpth
    RETURN
  end subroutine psb_lgps_reduce
  !
  SUBROUTINE DGREE(NDSTK, NR, NDEG, IOLD, IBW1, IPF1)
    !  DGREE COMPUTES THE DEGREE OF EACH NODE IN NDSTK AND STORES
    !  IT IN THE ARRAY NDEG.  THE BANDWIDTH AND PROFILE FOR THE ORIGINAL
    !  OR INPUT RENUMBERING OF THE GRAPH IS COMPUTED ALSO.
    ! USE INTEGER*2 NDSTK  WITH AN IBM 360 OR 370.
    implicit none
    INTEGER(psb_lpk_) :: NR, IBW1, IPF1, NDSTK(NR,IDEG), NDEG(N), IOLD(N)
    ! COMMON /GRA/ N, IDPTH, IDEG
    integer(psb_lpk_) :: i, itst, j, idif, irw

    IBW1 = 0
    IPF1 = 0
    DO  I=1,N
      NDEG(I) = 0
      IRW = 0
      DO J=1,IDEG
        ITST = NDSTK(I,J)
        IF(ITST <= 0) EXIT
        NDEG(I) = NDEG(I) + 1
        IDIF = IOLD(I) - IOLD(ITST)
        IF (IRW < IDIF) IRW = IDIF
      END DO
      IPF1 = IPF1 + IRW
      IF (IRW > IBW1) IBW1 = IRW
    END DO
    RETURN
  END SUBROUTINE DGREE
  !
  SUBROUTINE FNDIAM(SND1, SND2, NDSTK, NR, NDEG, LVL, LVLS1,LVLS2, IWK, IDFLT)
    !  FNDIAM IS THE CONTROL PROCEDURE FOR FINDING THE PSEUDO-DIAMETER OF
    !  NDSTK AS WELL AS THE LEVEL STRUCTURE FROM EACH END
    !  SND1-        ON INPUT THIS IS THE NODE NUMBER OF THE FIRST
    !               ATTEMPT AT FINDING A DIAMETER.  ON OUTPUT IT
    !               CONTAINS THE ACTUAL NUMBER USED.
    !  SND2-        ON OUTPUT CONTAINS OTHER END OF DIAMETER
    !  LVLS1-       ARRAY CONTAINING LEVEL STRUCTURE WITH SND1 AS ROOT
    !  LVLS2-       ARRAY CONTAINING LEVEL STRUCTURE WITH SND2 AS ROOT
    !  IDFLT-       FLAG USED IN PICKING FINAL LEVEL STRUCTURE, SET
    !               =1 IF WIDTH OF LVLS1 <= WIDTH OF LVLS2, OTHERWISE =2
    !  LVL,IWK-     WORKING STORAGE
    ! USE INTEGER*2 NDSTK  WITH AN IBM 360 OR 370.
    implicit none
    INTEGER(psb_lpk_) :: FLAG, SND, SND1, SND2, NR, idflt
    ! COMMON /GRA/ N, IDPTH, IDEG
    ! IT IS ASSUMED THAT THE LAST LEVEL HAS AT MOST 100 NODES.
    ! COMMON /CC/ NDLST(100)
    integer(psb_lpk_),POINTER :: NDLST(:)
    integer(psb_lpk_) :: NDSTK(NR,IDEG), NDEG(1), LVL(N), LVLS1(N), LVLS2(N),IWK(N)
    integer(psb_lpk_) :: i, j, mtw2, ndxn, ndxl, inow, lvlbot,lvln,lvlwth
    integer(psb_lpk_) :: k,mtw1, maxlw
    !
    NDLST => AUX
    !
    FLAG = 0
    MTW2 = N
    SND = SND1
    ! ZERO LVL TO INDICATE ALL NODES ARE AVAILABLE TO TREE
10  DO 20 I=1,N
       LVL(I) = 0
20  END DO
    LVLN = 1
    ! DROP A TREE FROM SND
    CALL TREE(SND, NDSTK, NR, LVL, IWK, NDEG, LVLWTH, LVLBOT,LVLN, MAXLW, MTW2)
    IF (FLAG >= 1) GO TO 50
    FLAG = 1
30  IDPTH = LVLN - 1
    MTW1 = MAXLW
    ! COPY LEVEL STRUCTURE INTO LVLS1
    DO 40 I=1,N
       LVLS1(I) = LVL(I)
40  END DO
    NDXN = 1
    NDXL = 0
    MTW2 = N
    ! SORT LAST LEVEL BY DEGREE  AND STORE IN NDLST
    CALL SORTDG(NDLST, IWK(LVLBOT), NDXL, LVLWTH, NDEG)
    SND = NDLST(1)
    GO TO 10
50  IF (IDPTH >= LVLN-1) GO TO 60
    ! START AGAIN WITH NEW STARTING NODE
    SND1 = SND
    GO TO 30
60  IF (MAXLW >= MTW2) GO TO 80
    MTW2 = MAXLW
    SND2 = SND
    ! STORE NARROWEST REVERSE LEVEL STRUCTURE IN LVLS2
    DO 70 I=1,N
       LVLS2(I) = LVL(I)
70  END DO
80  IF (NDXN == NDXL) GO TO 90
    ! TRY NEXT NODE IN NDLST
    NDXN = NDXN + 1
    SND = NDLST(NDXN)
    GO TO 10
90  IDFLT = 1
    IF (MTW2 <= MTW1) IDFLT = 2
    NULLIFY(NDLST)
    RETURN
  END SUBROUTINE FNDIAM
  !
  SUBROUTINE TREE(IROOT, NDSTK, NR, LVL, IWK, NDEG, LVLWTH, LVLBOT, LVLN, MAXLW, IBORT)
    !  TREE DROPS A TREE IN NDSTK FROM IROOT
    !  LVL-         ARRAY INDICATING AVAILABLE NODES IN NDSTK WITH ZERO
    !               ENTRIES. TREE ENTERS LEVEL NUMBERS ASSIGNED
    !               DURING EXECUTION OF THIS PROCEDURE
    !  IWK-         ON OUTPUT CONTAINS NODE NUMBERS USED IN TREE
    !               ARRANGED BY LEVELS (IWK(LVLN) CONTAINS IROOT
    !               AND IWK(LVLBOT+LVLWTH-1) CONTAINS LAST NODE ENTERED)
    !  LVLWTH-      ON OUTPUT CONTAINS WIDTH OF LAST LEVEL
    !  LVLBOT-      ON OUTPUT CONTAINS INDEX INTO IWK OF FIRST
    !               NODE IN LAST LEVEL
    !  MAXLW-       ON OUTPUT CONTAINS THE MAXIMUM LEVEL WIDTH
    !  LVLN-        ON INPUT THE FIRST AVAILABLE LOCATION IN IWK
    !               USUALLY ONE BUT IF IWK IS USED TO STORE PREVIOUS
    !               CONNECTED COMPONENTS, LVLN IS NEXT AVAILABLE LOCATION.
    !               ON OUTPUT THE TOTAL NUMBER OF LEVELS + 1
    !  IBORT-       INPUT PARAM WHICH TRIGGERS EARLY RETURN IF
    !               MAXLW BECOMES >= IBORT
    ! USE INTEGER*2 NDSTK  WITH AN IBM 360 OR 370.
    implicit none
    integer(psb_lpk_) :: IROOT, NR, NDSTK(NR,*), LVL(*), IWK(*), NDEG(*)
    integer(psb_lpk_) :: LVLWTH, LVLBOT, LVLN, MAXLW, IBORT
    integer(psb_lpk_) :: itest, iwknow, itop, lvltop,j , inow, ndrow
    MAXLW = 0
    ITOP = LVLN
    INOW = LVLN
    LVLBOT = LVLN
    LVLTOP = LVLN + 1
    LVLN = 1
    LVL(IROOT) = 1
    IWK(ITOP) = IROOT
10  LVLN = LVLN + 1
20  IWKNOW = IWK(INOW)
    NDROW = NDEG(IWKNOW)
    DO 30 J=1,NDROW
       ITEST = NDSTK(IWKNOW,J)
       IF (LVL(ITEST) /= 0) CYCLE
       LVL(ITEST) = LVLN
       ITOP = ITOP + 1
       IWK(ITOP) = ITEST
30  END DO
    INOW = INOW + 1
    IF (INOW < LVLTOP) GO TO 20
    LVLWTH = LVLTOP - LVLBOT
    IF (MAXLW < LVLWTH) MAXLW = LVLWTH
    IF (MAXLW >= IBORT) RETURN
    IF (ITOP < LVLTOP) RETURN
    LVLBOT = INOW
    LVLTOP = ITOP + 1
    GO TO 10
  END SUBROUTINE TREE
  !
  SUBROUTINE SORTDG(STK1, STK2, X1, X2, NDEG)
    ! SORTDG SORTS STK2 BY DEGREE OF THE NODE AND ADDS IT TO THE END
    ! OF STK1 IN ORDER OF LOWEST TO HIGHEST DEGREE.  X1 AND X2 ARE THE
    ! NUMBER OF NODES IN STK1 AND STK2 RESPECTIVELY.
    implicit none
    INTEGER(psb_lpk_) :: X1, X2, STK1, STK2, TEMP,NDEG
    ! COMMON /GRA/ N, IDPTH, IDEG
    DIMENSION NDEG(N), STK1(X1+X2), STK2(X2)
    integer(psb_lpk_) :: ind,itest,i,j,istk2,jstk2
    IND = X2
10  ITEST = 0
    IND = IND - 1
    IF (IND < 1) GO TO 30
    DO 20 I=1,IND
       J = I + 1
       ISTK2 = STK2(I)
       JSTK2 = STK2(J)
       IF (NDEG(ISTK2) <= NDEG(JSTK2)) CYCLE
       ITEST = 1
       TEMP = STK2(I)
       STK2(I) = STK2(J)
       STK2(J) = TEMP
20  END DO
    IF (ITEST == 1) GO TO 10
30  DO 40 I=1,X2
       X1 = X1 + 1
       STK1(X1) = STK2(I)
40  END DO
    RETURN
  END SUBROUTINE SORTDG
  !
  SUBROUTINE SETUP(LVL, LVLS1, LVLS2)
    ! SETUP COMPUTES THE REVERSE LEVELING INFO FROM LVLS2 AND STORES
    ! IT INTO LVLS2.  NACUM(I) IS INITIALIZED TO NODES/ITH LEVEL FOR NODES
    ! ON THE PSEUDO-DIAMETER OF THE GRAPH.  LVL IS INITIALIZED TO NON-
    ! ZERO FOR NODES ON THE PSEUDO-DIAM AND NODES IN A DIFFERENT
    ! COMPONENT OF THE GRAPH.
    !      COMMON /GRA/ N, IDPTH, IDEG
    ! IT IS ASSUMED THAT THERE ARE AT MOST 100 LEVELS.
    ! COMMON /LVLW/ NHIGH(100), NLOW(100), NACUM(100)
    use psb_base_mod
    implicit none
    integer(psb_lpk_) :: LVL(N), LVLS1(N), LVLS2(N)
    integer(psb_lpk_) :: SZ,i,itemp
    !-----------------------------------------------------
    SZ=SIZE(NACUM)
    IF(SZ < IDPTH) THEN
       write(psb_out_unit,*) 'GPS_SETUP: on fly reallocation of NACUM'
       CALL REALLOC(NACUM,SZ,IDPTH)
    END IF
    !-----------------------------------------------------
    DO 10 I=1,IDPTH
       NACUM(I) = 0
10  END DO
    DO 30 I=1,N
       LVL(I) = 1
       LVLS2(I) = IDPTH + 1 - LVLS2(I)
       ITEMP = LVLS2(I)
       IF (ITEMP > IDPTH) CYCLE
       IF (ITEMP /= LVLS1(I)) GO TO 20
       NACUM(ITEMP) = NACUM(ITEMP) + 1
       CYCLE
20     LVL(I) = 0
30  END DO
    RETURN
  END SUBROUTINE SETUP
  !
  FUNCTION SORT2() result(val)
    implicit none
    INTEGER(psb_lpk_) ::  val
    ! SORT2 SORTS SIZEG AND STPT INTO DESCENDING ORDER ACCORDING TO
    ! VALUES OF SIZEG. XCC=NUMBER OF ENTRIES IN EACH ARRAY
    INTEGER(psb_lpk_) :: TEMP,ind,itest,i,j
    ! IT IS ASSUMED THAT THE GRAPH HAS AT MOST 50 CONNECTED COMPONENTS.
    !COMMON /CC/ XCC, SIZEG(50), STPT(50)

    VAL = 0
    IF (XCC == 0) RETURN
    VAL = 1
    IND = XCC
10  ITEST = 0
    IND = IND - 1
    IF (IND < 1) RETURN
    DO 20 I=1,IND
       J = I + 1
       IF (SIZEG(I) >= SIZEG(J)) CYCLE
       ITEST = 1
       TEMP = SIZEG(I)
       SIZEG(I) = SIZEG(J)
       SIZEG(J) = TEMP
       TEMP = STPT(I)
       STPT(I) = STPT(J)
       STPT(J) = TEMP
20  END DO
    IF (ITEST == 1) GO TO 10
    RETURN
  END FUNCTION SORT2
  !
  SUBROUTINE PIKLVL(LVLS1, LVLS2, CCSTOR, IDFLT, ISDIR)
    use psb_base_mod
    implicit none
    ! PIKLVL CHOOSES THE LEVEL STRUCTURE  USED IN NUMBERING GRAPH
    ! LVLS1-       ON INPUT CONTAINS FORWARD LEVELING INFO
    ! LVLS2-       ON INPUT CONTAINS REVERSE LEVELING INFO
    !              ON OUTPUT THE FINAL LEVEL STRUCTURE CHOSEN
    ! CCSTOR-      ON INPUT CONTAINS CONNECTED COMPONENT INFO
    ! IDFLT-       ON INPUT =1 IF WDTH LVLS1 <= WDTH LVLS2, =2 OTHERWISE
    ! NHIGH        KEEPS TRACK OF LEVEL WIDTHS FOR HIGH NUMBERING
    ! NLOW-        KEEPS TRACK OF LEVEL WIDTHS FOR LOW NUMBERING
    ! NACUM-       KEEPS TRACK OF LEVEL WIDTHS FOR CHOSEN LEVEL STRUCTURE
    ! XCC-         NUMBER OF CONNECTED COMPONENTS
    ! SIZEG(I)-    SIZE OF ITH CONNECTED COMPONENT
    ! STPT(I)-     INDEX INTO CCSTORE OF 1ST NODE IN ITH CON COMPT
    ! ISDIR-       FLAG WHICH INDICATES WHICH WAY THE LARGEST CONNECTED
    !              COMPONENT FELL.  =+1 IF LOW AND -1 IF HIGH
    ! COMMON /GRA/ N, IDPTH, IDEG
    ! IT IS ASSUMED THAT THE GRAPH HAS AT MOST 50 COMPONENTS AND
    ! THAT THERE ARE AT MOST 100 LEVELS.
    ! COMMON /LVLW/ NHIGH(100), NLOW(100), NACUM(100)
    ! COMMON /CC/ XCC, SIZEG(50), STPT(50)
    INTEGER(psb_lpk_) :: LVLS1(N), LVLS2(N), CCSTOR(N)
    integer(psb_lpk_) :: SZ, ENDC,i,j,max1,max2,inode
    integer(psb_lpk_) :: lvlnh, it, k, lvlnl,idflt,isdir
    ! FOR EACH CONNECTED COMPONENT DO
    DO 80 I=1,XCC
       J = STPT(I)
       ENDC= SIZEG(I) + J - 1
       ! SET NHIGH AND NLOW EQUAL TO NACUM
       !-----------------------------------------------------
       SZ=SIZE(NHIGH)
       IF(SZ < IDPTH) THEN
          write(psb_out_unit,*) 'GPS_PIKLVL: on fly reallocation of NHIGH'
          CALL REALLOC(NHIGH,SZ,IDPTH)
       END IF
       SZ=SIZE(NLOW)
       IF(SZ < IDPTH) THEN
          write(psb_out_unit,*) 'GPS_PIKLVL: on fly reallocation of NLOW'
          CALL REALLOC(NLOW,SZ,IDPTH)
       END IF
       !-----------------------------------------------------
       DO 10 K=1,IDPTH
          NHIGH(K) = NACUM(K)
          NLOW(K) = NACUM(K)
10     END DO
       ! UPDATE NHIGH AND NLOW FOR EACH NODE IN CONNECTED COMPONENT
       DO 20 K=J,ENDC
          INODE = CCSTOR(K)
          LVLNH = LVLS1(INODE)
          NHIGH(LVLNH) = NHIGH(LVLNH) + 1
          LVLNL = LVLS2(INODE)
          NLOW(LVLNL) = NLOW(LVLNL) + 1
20     END DO
       MAX1 = 0
       MAX2 = 0
       ! SET MAX1=LARGEST NEW NUMBER IN NHIGH
       ! SET MAX2=LARGEST NEW NUMBER IN NLOW
       DO 30 K=1,IDPTH
          IF (2*NACUM(K) == NLOW(K)+NHIGH(K)) CYCLE
          IF (NHIGH(K) > MAX1) MAX1 = NHIGH(K)
          IF (NLOW(K) > MAX2) MAX2 = NLOW(K)
30     END DO
       ! SET IT= NUMBER OF LEVEL STRUCTURE TO BE USED
       IT = 1
       IF (MAX1 > MAX2) IT = 2
       IF (MAX1 == MAX2) IT = IDFLT
       IF (IT == 2) GO TO 60
       IF (I == 1) ISDIR = -1
       ! COPY LVLS1 INTO LVLS2 FOR EACH NODE IN CONNECTED COMPONENT
       DO 40 K=J,ENDC
          INODE = CCSTOR(K)
          LVLS2(INODE) = LVLS1(INODE)
40     END DO
       ! UPDATE NACUM TO BE THE SAME AS NHIGH
       DO 50 K=1,IDPTH
          NACUM(K) = NHIGH(K)
50     END DO
       CYCLE
       ! UPDATE NACUM TO BE THE SAME AS NLOW
60     DO 70 K=1,IDPTH
          NACUM(K) = NLOW(K)
70     END DO
80  END DO
    RETURN
  END SUBROUTINE PIKLVL
  !
  SUBROUTINE NUMBER(SND, NUM, NDSTK, LVLS2, NDEG, RENUM, LVLST,LSTPT,&
       & NR, NFLG, IBW2, IPF2, IPFA, ISDIR)
    use psb_base_mod
    implicit none
    !  NUMBER PRODUCES THE NUMBERING OF THE GRAPH FOR MIN BANDWIDTH
    !  SND-         ON INPUT THE NODE TO BEGIN NUMBERING ON
    !  NUM-         ON INPUT AND OUTPUT, THE NEXT AVAILABLE NUMBER
    !  LVLS2-       THE LEVEL STRUCTURE TO BE USED IN NUMBERING
    !  RENUM-       THE ARRAY USED TO STORE THE NEW NUMBERING
    !  LVLST-       ON OUTPUT CONTAINS LEVEL STRUCTURE
    !  LSTPT(I)-    ON OUTPUT, INDEX INTO LVLST TO FIRST NODE IN ITH LVL
    !               LSTPT(I+1) - LSTPT(I) = NUMBER OF NODES IN ITH LVL
    !  NFLG-        =+1 IF SND IS FORWARD END OF PSEUDO-DIAM
    !               =-1 IF SND IS REVERSE END OF PSEUDO-DIAM
    !  IBW2-        BANDWIDTH OF NEW NUMBERING COMPUTED BY NUMBER
    !  IPF2-        PROFILE OF NEW NUMBERING COMPUTED BY NUMBER
    !  IPFA-        WORKING STORAGE USED TO COMPUTE PROFILE AND BANDWIDTH
    !  ISDIR-       INDICATES STEP DIRECTION USED IN NUMBERING(+1 OR -1)
    ! USE INTEGER*2 NDSTK  WITH AN IBM 360 OR 370.
    INTEGER(psb_lpk_) :: SND, NUM, XA, XB, XC, XD, CX, ENDC, TEST, NR, ISDIR
    ! COMMON /GRA/ N, IDPTH, IDEG
    ! THE STORAGE IN COMMON BLOCKS CC AND LVLW IS NOW FREE AND CAN
    ! BE USED FOR STACKS.
    !COMMON /LVLW/ STKA(100), STKB(100), STKC(100)
    !COMMON /CC/ STKD(100)
    INTEGER(psb_lpk_) :: IPFA(N), NDSTK(NR,IDEG), LVLS2(N),&
         & NDEG(N), RENUM(N+1), LVLST(N),LSTPT(N),ipf2,ibw2,nflg, nbw
    integer(psb_lpk_),POINTER :: STKA(:),STKB(:),STKC(:),STKD(:)
    integer(psb_lpk_) :: SZ1,SZ2,i,j,nstpt,lvln, lst, lnd, inx, max, ipro,&
         & lvlnl, k, it
    !
    STKA => NHIGH
    STKB => NLOW
    STKC => NACUM
    STKD => AUX
    !
    ! SET UP LVLST AND LSTPT FROM LVLS2
    DO 10 I=1,N
       IPFA(I) = 0
10  END DO
    NSTPT = 1
    DO 30 I=1,IDPTH
       LSTPT(I) = NSTPT
       DO 20 J=1,N
          IF (LVLS2(J) /= I) CYCLE
          LVLST(NSTPT) = J
          NSTPT = NSTPT + 1
20     END DO
30  END DO
    LSTPT(IDPTH+1) = NSTPT
    ! STKA, STKB, STKC AND STKD ARE STACKS WITH POINTERS
    ! XA,XB,XC, AND XD.  CX IS A SPECIAL POINTER INTO STKC WHICH
    ! INDICATES THE PARTICULAR NODE BEING PROCESSED.
    ! LVLN KEEPS TRACK OF THE LEVEL WE ARE WORKING AT.
    ! INITIALLY STKC CONTAINS ONLY THE INITIAL NODE, SND.
    LVLN = 0
    IF (NFLG < 0) LVLN = IDPTH + 1
    XC = 1
    STKC(XC) = SND
40  CX = 1
    XD = 0
    LVLN = LVLN + NFLG
    LST = LSTPT(LVLN)
    LND = LSTPT(LVLN+1) - 1
    ! BEGIN PROCESSING NODE STKC(CX)
50  IPRO = STKC(CX)
    RENUM(IPRO) = NUM
    NUM = NUM + ISDIR
    ENDC = NDEG(IPRO)
    XA = 0
    XB = 0
    ! CHECK ALL ADJACENT NODES
    DO 80 I=1,ENDC
       TEST = NDSTK(IPRO,I)
       INX = RENUM(TEST)
       ! ONLY NODES NOT NUMBERED OR ALREADY ON A STACK ARE ADDED
       IF (INX == 0) GO TO 60
       IF (INX < 0) CYCLE
       ! DO PRELIMINARY BANDWIDTH AND PROFILE CALCULATIONS
       NBW = (RENUM(IPRO)-INX)*ISDIR
       IF (ISDIR > 0) INX = RENUM(IPRO)
       IF (IPFA(INX) < NBW) IPFA(INX) = NBW
       CYCLE
60     RENUM(TEST) = -1
       ! PUT NODES ON SAME LEVEL ON STKA, ALL OTHERS ON STKB
       IF (LVLS2(TEST) == LVLS2(IPRO)) GO TO 70
       XB = XB + 1
       STKB(XB) = TEST
       CYCLE
70     XA = XA + 1
       STKA(XA) = TEST
80  END DO
    ! SORT STKA AND STKB INTO INCREASING DEGREE AND ADD STKA TO STKC
    ! AND STKB TO STKD
    IF (XA == 0) GO TO 100
    IF (XA == 1) GO TO 90
    !-----------------------------------------------------------------
    SZ1=SIZE(STKC)
    SZ2=XC+XA
    IF(SZ1 < SZ2) THEN
       write(psb_out_unit,*) 'GPS_NUMBER - Check #1: on fly reallocation of STKC'
       CALL REALLOC(NACUM,SZ1,SZ2)
       STKC => NACUM
    END IF
    !-----------------------------------------------------------------
    CALL SORTDG(STKC, STKA, XC, XA, NDEG)
    GO TO 100
90  XC = XC + 1
    !-----------------------------------------------------------------
    SZ1=SIZE(STKC)
    SZ2=XC
    IF(SZ1 < SZ2) THEN
       write(psb_out_unit,*) 'GPS_NUMBER - Check #2: on fly reallocation of STKC'
       SZ2=SZ2+INIT
       CALL REALLOC(NACUM,SZ1,SZ2)
       STKC => NACUM
    END IF
    !-----------------------------------------------------------------
    STKC(XC) = STKA(XA)
100 IF (XB == 0) GO TO 120
    IF (XB == 1) GO TO 110
    !-----------------------------------------------------------------
    SZ1=SIZE(STKD)
    SZ2=XD+XB
    IF(SZ1 < SZ2) THEN
       write(psb_out_unit,*) 'GPS_NUMBER - Check #3: on fly reallocation of STKD'
       CALL REALLOC(AUX,SZ1,SZ2)
       STKD => AUX
    END IF
    !-----------------------------------------------------------------
    CALL SORTDG(STKD, STKB, XD, XB, NDEG)
    GO TO 120
110 XD = XD + 1
    !-----------------------------------------------------------------
    SZ1=SIZE(STKD)
    SZ2=XD
    IF(SZ1 < SZ2) THEN
       write(psb_out_unit,*) 'GPS_NUMBER - Check #4: on fly reallocation of STKD'
       SZ2=SZ2+INIT
       CALL REALLOC(AUX,SZ1,SZ2)
       STKD => AUX
    END IF
    !-----------------------------------------------------------------
    STKD(XD) = STKB(XB)
    ! BE SURE TO PROCESS ALL NODES IN STKC
120 CX = CX + 1
    IF (XC >= CX) GO TO 50
    ! WHEN STKC IS EXHAUSTED LOOK FOR MIN DEGREE NODE IN SAME LEVEL
    ! WHICH HAS NOT BEEN PROCESSED
    MAX = IDEG + 1
    SND = N + 1
    DO 130 I=LST,LND
       TEST = LVLST(I)
       IF (RENUM(TEST) /= 0) CYCLE
       IF (NDEG(TEST) >= MAX) CYCLE
       RENUM(SND) = 0
       RENUM(TEST) = -1
       MAX = NDEG(TEST)
       SND = TEST
130 END DO
    IF (SND == N+1) GO TO 140
    XC = XC + 1
    !-----------------------------------------------------------------
    SZ1=SIZE(STKC)
    SZ2=XC
    IF(SZ1 < SZ2) THEN
       write(psb_out_unit,*) 'GPS_NUMBER - Check #5: on fly reallocation of STKC'
       SZ2=SZ2+INIT
       CALL REALLOC(NACUM,SZ1,SZ2)
       STKC => NACUM
    END IF
    !-----------------------------------------------------------------
    STKC(XC) = SND
    GO TO 50
    ! IF STKD IS EMPTY WE ARE DONE, OTHERWISE COPY STKD ONTO STKC
    ! AND BEGIN PROCESSING NEW STKC
140 IF (XD == 0) GO TO 160
    !-----------------------------------------------------------------
    SZ1=SIZE(STKC)
    SZ2=XD
    IF(SZ1 < SZ2) THEN
       write(psb_out_unit,*) 'GPS_NUMBER - Check #6: on fly reallocation of STKC'
       SZ2=SZ2+INIT
       CALL REALLOC(NACUM,SZ1,SZ2)
       STKC => NACUM
    END IF
    !-----------------------------------------------------------------
    DO 150 I=1,XD
       STKC(I) = STKD(I)
150 END DO
    XC = XD
    GO TO 40
    ! DO FINAL BANDWIDTH AND PROFILE CALCULATIONS
160 DO 170 I=1,N
       IF (IPFA(I) > IBW2) IBW2 = IPFA(I)
       IPF2 = IPF2 + IPFA(I)
170 END DO
    !
    RETURN
  END SUBROUTINE NUMBER
  !
  ! ---------------------------------------------------------
  SUBROUTINE REALLOC(VET,SZ1,SZ2)
    use psb_base_mod
    ! PERFORM ON FLY REALLOCATION OF POINTER VET INCREASING
    ! ITS SIZE FROM SZ1 TO SZ2
    IMPLICIT NONE
    integer(psb_lpk_),allocatable  :: VET(:)
    integer(psb_lpk_) :: SZ1,SZ2
    integer(psb_ipk_) :: info
    call psb_realloc(sz2,vet,info)
    IF(INFO /= psb_success_) THEN
       write(psb_out_unit,*) 'Error! Memory allocation failure in REALLOC'
       STOP
    END IF
    RETURN
  END SUBROUTINE REALLOC
  !
end module psb_lgps_mod


text:
        subroutine bc3per (iBC,  res, iper, ilwork,nQs)
c
c----------------------------------------------------------------------
c
c This routine satisfies the BC of the periodic nodes after Ap product
c
c input:
c  iBC   (nshg)        : Boundary Condition Code
c  iper  (nshg)        : partners of periodic nodes
c  res   (nshg,nQs)   : residual before BC is applied
c
c output:
c  res   (nshg,nQs)   : residual after satisfaction of BC
c
c
c Kenneth Jansen,  Winter 1998.
c----------------------------------------------------------------------
c
        use periodicity  ! this gives you rcount(1:nshg) (real*8)
        include "common.h"
c
        dimension iBC(nshg),
     &            res(nshg,nQs),           ilwork(nlwork),
     &            iper(nshg)
c
c.... local periodic (and axisymmetric) boundary conditions (no communications)
c
           do j = 1,nshg
              if ((btest(iBC(j),10)) .or. (btest(iBC(j),12))) then
                 i = iper(j)
                 res(i,:) = res(i,:) + res(j,:)
                 res(j,:) = zero
              endif
           enddo

        if(numpe.gt.1) then
c
c.... nodes treated on another processor are eliminated
c
           numtask = ilwork(1)
           itkbeg = 1

           do itask = 1, numtask

              iacc   = ilwork (itkbeg + 2)
              numseg = ilwork (itkbeg + 4)

              if (iacc .eq. 0) then
                 do is = 1,numseg
                    isgbeg = ilwork (itkbeg + 3 + 2*is)
                    lenseg = ilwork (itkbeg + 4 + 2*is)
                    isgend = isgbeg + lenseg - 1
                    res(isgbeg:isgend,:) = zero
                 enddo
              endif

              itkbeg = itkbeg + 4 + 2*numseg

           enddo
        endif
c
c.... return
c
        return
        end

ext: f


program tsp_recursive
    implicit none
    integer :: n, min_distance, i
    integer, allocatable :: d(:, :)
    integer, allocatable :: min_route(:), route(:)

    print '(a)', "How many towns?"
    read *, n

    allocate(d(n, n))
    print '(a)', "Distances between the towns are"
    do i = 1, n
        read *, d(i, :)
    end do

    min_distance = huge(0)
    allocate(min_route(n))
    allocate(route(n))

    route(:) = [ (i, i = 1, n) ] ! [1, 2, 3, ..., n]
    call permute(route, 1)

    print '(a, i0, a)', "The shortest distance ", min_distance, " is for route:"
    print '(*(i0, :, 2x))', min_route(:)

    deallocate(d)
    deallocate(min_route)
    deallocate(route)
contains
    subroutine swap(a, b)
        integer, intent(inout) :: a, b
        integer :: temp
        temp = a; a = b; b = temp
    end subroutine swap

    recursive subroutine permute(a, m)
        integer, intent(inout) :: a(:)
        integer, intent(in)    :: m
        integer :: asize, j, l

        asize = size(a)
        if (m == asize) then
            !print '(*(i2, 2x))', a(:)
            l = d(a(asize), a(1))
            do j = 1, asize - 1
                l = l + d(a(j), a(j + 1))
            end do
            if (l < min_distance) then
                min_distance = l
                min_route(:) = a(:)
            end if
        else
            do j = m, asize
                call swap(a(m), a(j))
                call permute(a, m + 1)
                call swap(a(m), a(j))
            end do
        end if
    end subroutine permute
end program tsp_recursive


!Forward Flux Sampling with the Rosenbluth algorithm
program ffs
use mpi
implicit none

integer,parameter:: Nw=20000 !number of initial conditions at R_A
double precision,parameter:: vo=9.0d0
integer, parameter:: N=80
double precision, parameter:: tau=1d-5 !timestep
double precision:: x0(N,2), theta0(N),xcont(N,2),thetacont(N)

!MPI
integer :: ierr,tag,cores,id,my_id,status(MPI_STATUS_SIZE),iexit !mpi stuff
double precision :: start,finish !keeps time

!initialization
integer, parameter:: frames=10000
integer:: iframe
double precision:: init(frames,N,3)
double precision:: rtemp !random number

!counters
integer:: i,j,lamcount,kcount,tempcounter

!ffs
integer,parameter:: nlam=12 !number of interfaces
double precision,parameter,dimension(nlam):: lams=(/1.25d0,1.29d0,1.33d0,&
        & 1.38d0,1.43d0,1.50d0,1.57d0,1.65d0,1.72d0,1.77d0,1.81d0,1.85d0/)
integer,parameter:: kvals=100 !each interface starts kvals trajs
integer:: Sarray(nlam-1) !number of successful trajs
integer:: indicator(kvals) !within one branching step lists which trajectories make it
double precision:: ttotA !total time spent in A during initialization
integer:: trajinit,trajlens(2) !time spent per individual trajectory doing initiatlization and propagation
integer,parameter:: dtmax=2d4 !maximum length of each segment
double precision:: trajx(kvals,N,2),trajtheta(kvals,N)
integer:: trajlen(kvals) !length of the trajectory segments
integer:: iinit !random number for which endpoint to start the next trajectories from
double precision:: xlam1(N,2),thetalam1(N) !coordinates at lam1
integer,parameter:: maxlam1=1d7 !maximum time spent in A
integer:: ind1,ind2 !indicator for crossing lam1 or reaching lamend
integer:: breakind,sumind !temporary variables

!random number stuff
integer :: sizerand
integer, allocatable :: putrand(:)

!mpi stuff
call MPI_Init(ierr)
call MPI_Comm_rank(MPI_COMM_WORLD,my_id,ierr)
call MPI_Comm_size(MPI_COMM_WORLD,cores,ierr)

if (cores/=kvals) then
    print *, 'number of threads must be same as kvals'
    stop
end if

!seed random numbers
call random_seed(size=sizerand)
allocate(putrand(sizerand))
putrand(:)=my_id+10
call random_seed(put=putrand)
do i=1,50
    call random_number(rtemp)
end do

if (my_id==0) then
    !read initial conditions
    open(unit=140,file='v9_frames.txt')
    do i=1,frames
      do j=1,N
          read(140,*)init(i,j,:)
      end do
    end do

    !output files
    open(unit=10,file='N0_ttot.txt')
    open(unit=20,file='Stot.txt')
    open(unit=30,file='dttot.txt')
end if

ttotA=0.0d0
tempcounter=0
ind1=2
do j=1,Nw

start=MPI_Wtime()
if (my_id==0) then !initialize only in main node

!start the initial segment loop here
ind2=0
trajinit=0
do while(ind1==2 .or. ind2==0) !denotes it ended in B well or that old traj has to be continued
if (ind1==2) then!choose a new frame
    call random_number(rtemp)
    iframe=floor(rtemp*frames)+1 !an integer between [1,frames]
    x0=init(iframe,:,1:2)
    theta0=init(iframe,:,3)
end if !else x0,theta0 already has previous info
!call a trajectory
call inittraj(N,tau,x0,theta0,vo,my_id,j,lams(1),lams(nlam),maxlam1,&
&       xlam1,thetalam1,trajlen(1),ind1,ind2,tempcounter)
!initial flux calculation, outputs total time instants spent, the initial traj and its length
if (ind1==2) write (*,*) 'ended in B for my_id,Nwcount',my_id,j

ttotA=ttotA+(trajlen(1)-1)*tau
trajinit=trajinit+trajlen(1)-1
tempcounter=tempcounter+1

end do !end initialization loop

if (j<Nw) then
    x0=xlam1
    theta0=thetalam1
end if

!AD
write(*,*) 'Initialization done for Nwcount,tempcounter,iframe',j,tempcounter,iframe
!call flush()
trajlens(:)=0

end if !other cores don't need to do initialization

breakind=0
do lamcount=2,nlam
if (breakind==1) exit
!trajs start from interface at lams(lamcount-1)

if (my_id==0) then !choose the starting point at main node
!choose a starting point randomly
if (lamcount==2) then
    xcont=xlam1
    thetacont=thetalam1
else
    do while(.true.)
        call random_number(rtemp)
        iinit=floor(rtemp*kvals)+1 !an integer between [1,kvals]
        if (indicator(iinit)==1) then
            xcont=trajx(iinit,:,:)
            thetacont=trajtheta(iinit,:)
            trajlens(2)=trajlens(2)+trajlen(iinit)-1
            exit
        end if
    end do
end if

end if !main node condition ends
!broadcast the init cond
call MPI_Barrier(MPI_COMM_WORLD,ierr)
call MPI_Bcast(xcont,N*2,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr)
call MPI_Bcast(thetacont,N,MPI_DOUBLE,0,MPI_COMM_WORLD,ierr)
call MPI_Barrier(MPI_COMM_WORLD,ierr)

indicator=0
call proptraj(N,tau,xcont,thetacont,&
        &       vo,my_id,j,lams(lamcount),lams(1),dtmax,&
        &       trajx(my_id+1,:,:),trajtheta(my_id+1,:),&
        &       trajlen(my_id+1),indicator(my_id+1),&
        &       lamcount,my_id+1)
if (trajlen(my_id+1)==dtmax .and. indicator(my_id+1)==0) then
    write(*,*) 'Short segment at my_id,Nwcount,lamcount,kcount',j,lamcount,my_id+1
end if
!print *, 'propagation done for Nw,lamcount,kcount',j,lamcount,my_id+1
!receive in main node trajx,trajtheta,trajlen,indicator from other nodes
if (my_id>0) then
    id=0
    call MPI_Send(trajx(my_id+1,:,:),N*2,MPI_DOUBLE,id,100*my_id+4*(lamcount-2)+1,MPI_COMM_WORLD,ierr)
    call MPI_Send(trajtheta(my_id+1,:),N,MPI_DOUBLE,id,100*my_id+4*(lamcount-2)+2,MPI_COMM_WORLD,ierr)
    call MPI_Send(indicator(my_id+1),1,MPI_INTEGER,id,100*my_id+4*(lamcount-2)+3,MPI_COMM_WORLD,ierr)
    call MPI_Send(trajlen(my_id+1),1,MPI_INTEGER,id,100*my_id+4*(lamcount-2)+4,MPI_COMM_WORLD,ierr)
else
    do i=1,cores-1
        id=i
        call MPI_Recv(trajx(id+1,:,:),N*2,MPI_DOUBLE,id,100*id+4*(lamcount-2)+1,MPI_COMM_WORLD,status,ierr)
        call MPI_Recv(trajtheta(id+1,:),N,MPI_DOUBLE,id,100*id+4*(lamcount-2)+2,MPI_COMM_WORLD,status,ierr)
        call MPI_Recv(indicator(id+1),1,MPI_INTEGER,id,100*id+4*(lamcount-2)+3,MPI_COMM_WORLD,status,ierr)
        call MPI_Recv(trajlen(id+1),1,MPI_INTEGER,id,100*id+4*(lamcount-2)+4,MPI_COMM_WORLD,status,ierr)
        !print *, 'received for core', i
    end do
    trajlens(1)=trajlens(1)+sum(trajlen)-kvals
    !write the results of this step
    sumind=sum(indicator)
    write(20,'(I5.1)',advance='no') sumind
    !print *, j,lamcount,sumind
    !if no trajectories make it, break the lamcount loop and start fresh
    if (sumind==0) then
        do i=lamcount+1,nlam
            write(20,'(I5.1)',advance='no') 0
        end do
        breakind=1
    end if

end if !printed results of this lamcount loop
!broadcast breakind
call MPI_Barrier(MPI_COMM_WORLD,ierr)
call MPI_Bcast(breakind,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
call MPI_Barrier(MPI_COMM_WORLD,ierr)

!print *, 'kloop done for lamcount,my_id is', lamcount,my_id
!call flush()
end do !end lamcount loop

finish=MPI_Wtime()
if (my_id==0) then
    print *, 'Nwcount, time taken is',j,finish-start
    write(20,*) !add newline
    write(30,*) trajinit,trajlens(:)
    !call flush()
end if

end do !end Nw loop
if (my_id==0) write(10,*) Nw,ttotA
!call flush()

call MPI_Finalize(ierr)

end program ffs
!*************************************************************************!
!*************************************************************************!
!*************************************************************************!
!subroutine for running a single trajectory
subroutine inittraj(N,tau,x0,theta0,vo,my_id,Nwcount,lam1,lamend,maxlam1,&
        &       xlam1,thetalam1,trajlen,ind1,ind2,tempcounter)

implicit none

!arguments
integer, intent(in):: N
double precision, intent(in):: x0(N,2),theta0(N)
double precision, intent(in):: lam1,lamend,vo,tau
integer, intent(in):: my_id,Nwcount,maxlam1,tempcounter
double precision:: xlam1(N,2),thetalam1(N)
integer:: trajlen,ind1,ind2

!trajectory parameters
!integer, parameter:: N=80
double precision, parameter:: sigma=1.0d0
double precision, parameter:: sigmap6=sigma**6, sigmap12=sigmap6**2
double precision, parameter:: cutwca=2**(1./6) *sigma
double precision, parameter:: rthick=0.5d0*sigma !rskin-rcut
double precision, parameter:: rskin=cutwca+rthick
double precision, parameter:: Dd=0.25d0,Db=0.5d0,D_theta=1.0d0 !diffusion constant for dimer and bath
double precision, parameter:: fd=2.0d0,fb=1.0d0 !friction coefficient for dimer and bath
double precision:: etavard,etavarb,eta_theta_var!prefactor for gaussian noise
double precision,parameter:: pi=4.0d0*atan(1.0d0)
double precision, parameter:: rho=0.6 !A density
!double precision, parameter:: L=(NA*pi* (sigmaA)**2/(4.0d0*rho))**(1.0d0/2.0d0)
double precision:: L!=(N/rho)**(1.0d0/2.0d0)

!force calc
double precision,dimension(N,2):: delx !holder for jump in every timestep
double precision,dimension(N):: deltheta
double precision,dimension(N,2):: F0,F0wca,F0dw,Fa !force
double precision:: potdw,rint,rintold
double precision:: frc

!random number stuff
integer :: sizerand
integer, allocatable :: putrand(:)
!double precision:: rand1(N,2),rand2(N,2),r(N,2),theta(N,2),eta(N,2)
double precision::rand1(N,2),rand2(N,2),rand3(N),rand4(N),r(N,2),thetaBM1(N,2),thetaBM2(N),eta(N,2),eta_theta(N)

!neighbor list
double precision, dimension(N,2):: dxnei
double precision:: drnei,drneimax,drneimax2
integer:: nnflag
integer,dimension(N,N):: nlist

!counters and temporary variables
integer:: i,j,k,p
double precision:: xr(2),absr
integer:: indicator,counter,boxx

!determine length of box
L=(N/rho)**(1.0d0/2.0d0)
etavard=(2.0d0*Dd*tau)**0.5d0
etavarb=(2.0d0*Db*tau)**0.5d0
eta_theta_var=(2.0d0*D_theta*tau)**0.5d0

!seed random numbers
call random_seed(size=sizerand)
allocate(putrand(sizerand))
putrand(:)=1d6*my_id+1d3*Nwcount+tempcounter
call random_seed(put=putrand)
do i=1,50
    call random_number(rand1(1,1))
end do
!x0(1,:)=0.0d0!0.5*L !first particle in the middle of the box
!x0(2,1)=0.0d0+0.25*2.0d0
!x0(2,2)=0.0d0
!do i=3,N
!!    print *, 'starting particle', i
!    indicator=1
!    do while (indicator==1)
!        call random_number(rand1(i,:)) !three random coordinates for ith particle
!        rand1(i,:)=(rand1(i,:)-0.5d0)*L !rescales the coordinates to the whole box
!        indicator=0
!        do j=1,i-1  !Check overlap with previous particles
!            xr=rand1(i,:)-x0(j,:)
!            !!periodic boundaries
!            xr(1)=xr(1)-L*nint(xr(1)/L)
!            xr(2)=xr(2)-L*nint(xr(2)/L)
!            absr=dot_product(xr,xr)**0.5
!            if (absr<0.95*sigma) indicator=1  !overlap!
!        end do
!    end do
!    x0(i,:)=rand1(i,:)
!    if (my_id==0) print *, 'ending particle', i
!end do

!boxx=floor(L/sigma) !side length of box
!counter=0
!do i=1,boxx
!    do j=1,boxx
!        counter=counter+1
!        x0(counter,1)=-0.5*L+sigma*(i-1)
!        x0(counter,2)=-0.5*L+sigma*(j-1)
!        if (my_id==0) print *, counter, x0(counter,:)
!        if (counter==N) exit
!    end do
!    if (counter==N) exit
!end do

frc=24.0d0*(2.0d0*sigmap12/cutwca**13-sigmap6/cutwca**7)
dxnei=1000.0d0
xlam1=x0
thetalam1=theta0
ind1=0
counter=0

!begin actual transition path trajectory
do while (.true.)

    counter=counter+1

    !generate random noise using Box Mueller algo
    if (mod(counter,2)==1) then
        call random_number(rand1)
        call random_number(rand2)
        r=(-2.0d0*log(rand2))**0.5d0
        thetaBM1=2.0d0*pi*rand1
        rand1=r*cos(thetaBM1)
        rand2=r*sin(thetaBM1)
        eta=rand1

        call random_number(rand3)
        call random_number(rand4)
        r(:,1)=(-2.0d0*log(rand4))**0.5d0
        thetaBM2=2.0d0*pi*rand3
        rand3=r(:,1)*cos(thetaBM2)
        rand4=r(:,1)*sin(thetaBM2)
        eta_theta=rand3
    else
        eta=rand2
        eta_theta=rand4
    end if

    !Updating neighbor list if required by calculating two maximum displacements
    drneimax=0.0
    drneimax2=0.0
    do j=1,N
        drnei=sqrt(dot_product(dxnei(j,:),dxnei(j,:)))
        if (drnei > drneimax) then
            drneimax2=drneimax
            drneimax=drnei
        else
            if (drnei > drneimax2) then
                drneimax2=drnei
            endif
        endif
    enddo
    if (drneimax+drneimax2 > rthick) then
        nnflag=1 !update neighbor list in the next call for force
    !print *,'Updating neighbor list at step', i
    end if

    if (nnflag==1) dxnei=0

    if (counter>1) rintold=rint

    !force calc
    call force(N,sigma,sigmap6,sigmap12,cutwca,&
    &   xlam1,F0,F0wca,F0dw,rint,L,frc,&
    &   rskin,nlist,nnflag)

    if (counter>1) then
        if (rint>lam1 .and. rintold<lam1) then
            ind1=1 !crossed the first interface left to right
            ind2=1
            exit
        end if
    end if

    if (rint>lamend) then
        ind1=2 !reached lam2 in previous loop
        exit
    end if

    if (counter==maxlam1+1) then
        write(*,*) 'my_id,Nwcount',my_id,Nwcount,'initial segment unsuccessful'
        ind1=2 !need to use a different initial condition !!this biases the trajectory ensemble
        exit
    end if

    Fa(:,1)=vo*cos(thetalam1)
    Fa(:,2)=vo*sin(thetalam1)
    Fa(1,:)=0.0d0
    Fa(2,:)=0.0d0
    F0=F0+Fa

    !propagate position
    delx(3:N,:)=tau*F0(3:N,:)/fb+etavarb*eta(3:N,:)
    delx(1:2,:)=tau*F0(1:2,:)/fd+etavard*eta(1:2,:)
    xlam1=xlam1+delx
    dxnei=dxnei+delx

    deltheta=eta_theta_var*eta_theta
    thetalam1=thetalam1+deltheta

    do j=1,N
        xlam1(j,1)=xlam1(j,1)-L*nint(xlam1(j,1)/L)
        xlam1(j,2)=xlam1(j,2)-L*nint(xlam1(j,2)/L)
    end do

end do !trajectory has reached lam1
trajlen=counter

end subroutine inittraj
!*************************************************************************!
!*************************************************************************!
!*************************************************************************!
!subroutine for running a single trajectory
subroutine proptraj(N,tau,x0,theta0,vo,my_id,Nwcount,lam,lam1,dtmax,&
        &       trajx,trajtheta,trajlen,indicator,lamcount,kcount)

implicit none

!arguments
integer, intent(in):: N
double precision, intent(in):: x0(N,2),theta0(N)
double precision, intent(in):: lam,lam1,vo,tau
integer, intent(in):: my_id,Nwcount,dtmax,lamcount,kcount
double precision:: trajx(N,2),trajtheta(N)
integer:: trajlen,indicator

!trajectory parameters
!integer, parameter:: N=80
double precision, parameter:: sigma=1.0d0
double precision, parameter:: sigmap6=sigma**6, sigmap12=sigmap6**2
double precision, parameter:: cutwca=2**(1./6) *sigma
double precision, parameter:: rthick=0.5d0*sigma !rskin-rcut
double precision, parameter:: rskin=cutwca+rthick
double precision, parameter:: Dd=0.25d0,Db=0.5d0,D_theta=1.0d0 !diffusion constant for dimer and bath
double precision, parameter:: fd=2.0d0,fb=1.0d0 !friction coefficient for dimer and bath
double precision :: etavard,etavarb,eta_theta_var!prefactor for gaussian noise
double precision,parameter:: pi=4.0d0*atan(1.0d0)
double precision, parameter:: rho=0.6 !A density
!double precision, parameter:: L=(NA*pi* (sigmaA)**2/(4.0d0*rho))**(1.0d0/2.0d0)
double precision:: L!=(N/rho)**(1.0d0/2.0d0)

!force calc
double precision,dimension(N,2):: delx !holder for jump in every timestep
double precision,dimension(N):: deltheta
double precision,dimension(N,2):: F0,F0wca,F0dw,Fa !force
double precision:: potdw,rint
double precision:: frc

!random number stuff
integer :: sizerand
integer, allocatable :: putrand(:)
!double precision:: rand1(N,2),rand2(N,2),r(N,2),theta(N,2),eta(N,2)
double precision::rand1(N,2),rand2(N,2),rand3(N),rand4(N),r(N,2),thetaBM1(N,2),thetaBM2(N),eta(N,2),eta_theta(N)

!neighbor list
double precision, dimension(N,2):: dxnei
double precision:: drnei,drneimax,drneimax2
integer:: nnflag
integer,dimension(N,N):: nlist

!counters and temporary variables
integer:: i,j,k,p
double precision:: xr(2),absr
integer:: ind1,ind2,counter,boxx

!determine length of box
L=(N/rho)**(1.0d0/2.0d0)
etavard=(2.0d0*Dd*tau)**0.5d0
etavarb=(2.0d0*Db*tau)**0.5d0
eta_theta_var=(2.0d0*D_theta*tau)**0.5d0

!seed random numbers
call random_seed(size=sizerand)
allocate(putrand(sizerand))
putrand(:)=1d6*my_id+1d4*Nwcount+1d2*lamcount+kcount !AD makes sense when kvals=100
call random_seed(put=putrand)
do i=1,50
    call random_number(rand1(1,1))
end do

!x0(1,:)=0.0d0!0.5*L !first particle in the middle of the box
!x0(2,1)=0.0d0+0.25*2.0d0
!x0(2,2)=0.0d0
!do i=3,N
!!    print *, 'starting particle', i
!    indicator=1
!    do while (indicator==1)
!        call random_number(rand1(i,:)) !three random coordinates for ith particle
!        rand1(i,:)=(rand1(i,:)-0.5d0)*L !rescales the coordinates to the whole box
!        indicator=0
!        do j=1,i-1  !Check overlap with previous particles
!            xr=rand1(i,:)-x0(j,:)
!            !!periodic boundaries
!            xr(1)=xr(1)-L*nint(xr(1)/L)
!            xr(2)=xr(2)-L*nint(xr(2)/L)
!            absr=dot_product(xr,xr)**0.5
!            if (absr<0.95*sigma) indicator=1  !overlap!
!        end do
!    end do
!    x0(i,:)=rand1(i,:)
!    if (my_id==0) print *, 'ending particle', i
!end do

!boxx=floor(L/sigma) !side length of box
!counter=0
!do i=1,boxx
!    do j=1,boxx
!        counter=counter+1
!        x0(counter,1)=-0.5*L+sigma*(i-1)
!        x0(counter,2)=-0.5*L+sigma*(j-1)
!        if (my_id==0) print *, counter, x0(counter,:)
!        if (counter==N) exit
!    end do
!    if (counter==N) exit
!end do

frc=24.0d0*(2.0d0*sigmap12/cutwca**13-sigmap6/cutwca**7)
dxnei=1000.0d0
trajx=x0
trajtheta=theta0
ind1=0 !this notes whether lam1 has been reached
indicator=0 !this notes whether lam has been reached
counter=0

!begin actual transition path trajectory
do while (.true.) !hasn't reached lam1 or lam

    counter=counter+1

    !generate random noise using Box Mueller algo
    if (mod(counter,2)==1) then
        call random_number(rand1)
        call random_number(rand2)
        r=(-2.0d0*log(rand2))**0.5d0
        thetaBM1=2.0d0*pi*rand1
        rand1=r*cos(thetaBM1)
        rand2=r*sin(thetaBM1)
        eta=rand1

        call random_number(rand3)
        call random_number(rand4)
        r(:,1)=(-2.0d0*log(rand4))**0.5d0
        thetaBM2=2.0d0*pi*rand3
        rand3=r(:,1)*cos(thetaBM2)
        rand4=r(:,1)*sin(thetaBM2)
        eta_theta=rand3
    else
        eta=rand2
        eta_theta=rand4
    end if

    !Updating neighbor list if required by calculating two maximum displacements
    drneimax=0.0
    drneimax2=0.0
    do j=1,N
        drnei=sqrt(dot_product(dxnei(j,:),dxnei(j,:)))
        if (drnei > drneimax) then
            drneimax2=drneimax
            drneimax=drnei
        else
            if (drnei > drneimax2) then
                drneimax2=drnei
            endif
        endif
    enddo
    if (drneimax+drneimax2 > rthick) then
        nnflag=1 !update neighbor list in the next call for force
    !print *,'Updating neighbor list at step', i
    end if

    if (nnflag==1) dxnei=0

    !force calc
    call force(N,sigma,sigmap6,sigmap12,cutwca,&
    &   trajx,F0,F0wca,F0dw,rint,L,frc,&
    &   rskin,nlist,nnflag)

    if (rint<lam1) then
        ind1=1 !reached lam1 in previous loop
        exit
    else if (rint>lam) then
        indicator=1 !reached lam in previous loop
        exit
    end if

    if (counter==dtmax+1) exit

    Fa(:,1)=vo*cos(trajtheta)
    Fa(:,2)=vo*sin(trajtheta)
    Fa(1,:)=0.0d0
    Fa(2,:)=0.0d0
    F0=F0+Fa

    !propagate position
    delx(3:N,:)=tau*F0(3:N,:)/fb+etavarb*eta(3:N,:)
    delx(1:2,:)=tau*F0(1:2,:)/fd+etavard*eta(1:2,:)
    trajx=trajx+delx
    dxnei=dxnei+delx

    deltheta=eta_theta_var*eta_theta
    trajtheta=trajtheta+deltheta

    do j=1,N
        trajx(j,1)=trajx(j,1)-L*nint(trajx(j,1)/L)
        trajx(j,2)=trajx(j,2)-L*nint(trajx(j,2)/L)
    end do

end do !trajectory has reached lam or lam1
trajlen=counter

end subroutine proptraj
!*************************************************************************!
!*************************************************************************!
!*************************************************************************!
!*************************************************************************!
!subroutine for force calculation
subroutine force(N,sigma,sigmap6,sigmap12,cutwca,&
    &   x,F0,F0wca,F0dw,rint,L,frc,&
    &   rskin,nlist,nnflag)
implicit none

integer, intent(in):: N
double precision, intent(in):: x(N,2),L,frc,rskin
integer:: nnflag, nlist(N,N)
double precision, intent(in):: sigma,sigmap6,sigmap12,cutwca
double precision:: F0(N,2),F0wca(N,2),F0dw(N,2)
double precision:: potdw,rint
integer:: i,j,k

double precision:: xr(2),rcap(2),absr,rprod1,rprod2,Fint(2),fac1,fac2
double precision, parameter:: ewca=1.0d0,hdw=3.5d0,wdw=0.45d0

F0=0.0d0
F0wca=0.d0
F0dw=0.0d0

!generic forces common for all particles
if (nnflag==0) then !do calculation according to neighbor list

do i=1,N
!put 1body force here

!do j=i+1,N
do k=1,nlist(i,1)
  j=nlist(i,k+1)
  if (j<=i) cycle

!!define absr and xr right away and work with those
xr=x(i,:)-x(j,:)
!!periodic boundaries
xr(1)=xr(1)-nint(xr(1)/L)*L
xr(2)=xr(2)-nint(xr(2)/L)*L
absr=dot_product(xr,xr)**0.5
rcap=xr/absr

!!WCA evaluation
if (absr<cutwca) then
rprod1=absr*absr
rprod1=rprod1*rprod1
rprod1=rprod1*rprod1/absr !7th power
rprod2=rprod1*rprod1/absr !13th power
Fint=ewca*(24.0d0*(2.0d0*sigmap12/rprod2-sigmap6/rprod1)-frc)*rcap
F0wca(i,:)=F0wca(i,:)+Fint
F0wca(j,:)=F0wca(j,:)-Fint
end if
end do

end do

else !make new neighbor list
nlist(:,1)=0

do i=1,N
!put 1body force here

do j=i+1,N

!!define absr and xr right away and work with those
xr=x(i,:)-x(j,:)
!!periodic boundaries
xr(1)=xr(1)-nint(xr(1)/L)*L
xr(2)=xr(2)-nint(xr(2)/L)*L
absr=dot_product(xr,xr)**0.5
rcap=xr/absr

if (absr<=rskin) then !i and j are neighbors
nlist(i,nlist(i,1)+2)=j
nlist(j,nlist(j,1)+2)=i
nlist(i,1)=nlist(i,1)+1
nlist(j,1)=nlist(j,1)+1

!!WCA evaluation
if (absr<cutwca) then
rprod1=absr*absr
rprod1=rprod1*rprod1
rprod1=rprod1*rprod1/absr !7th power
rprod2=rprod1*rprod1/absr !13th power
Fint=ewca*(24.0d0*(2.0d0*sigmap12/rprod2-sigmap6/rprod1)-frc)*rcap
F0wca(i,:)=F0wca(i,:)+Fint
F0wca(j,:)=F0wca(j,:)-Fint
end if
end if !end calculation for i and j being neighbors

end do !loop over 2body forces

end do !loop over 1body forces
nnflag=0
end if

!double well forces for first two particles
xr=x(1,:)-x(2,:)
!!periodic boundaries
xr(1)=xr(1)-nint(xr(1)/L)*L
xr(2)=xr(2)-nint(xr(2)/L)*L
absr=dot_product(xr,xr)**0.5
rcap=xr/absr
fac1=absr-cutwca-wdw
fac2=1.0d0-(fac1)**2/wdw/wdw
rint=absr
potdw=hdw*fac2*fac2
Fint=4.0d0*hdw/wdw/wdw *fac1*fac2*rcap
F0dw(1,:)=Fint
F0dw(2,:)=-Fint

F0=F0wca+F0dw

end subroutine force
!*************************************************************************!
!*************************************************************************!
!*************************************************************************!
!*************************************************************************!
!subroutine for calculating total potential
subroutine potential(N,sigma,sigmap6,sigmap12,cutwca,x0,L,frc,pot0)
implicit none
integer, intent(in):: N
double precision, intent(in)::sigma,sigmap6,sigmap12,cutwca,x0(N,2),L,frc
double precision:: urc,pot0
integer:: i,j,k

double precision:: xr(2),absr,rprod1,rprod2,fac1,fac2
double precision, parameter:: ewca=1.0d0,hdw=3.5d0,wdw=0.45d0

urc=-1.0d0 !for wca potential
pot0=0.0d0

do i=1,N
do j=i+1,N

!!define absr and xr right away and work with those
xr=x0(i,:)-x0(j,:)
!!periodic boundaries
xr(1)=xr(1)-nint(xr(1)/L)*L
xr(2)=xr(2)-nint(xr(2)/L)*L
absr=dot_product(xr,xr)**0.5

!!WCA evaluation
if (absr<cutwca) then
rprod1=absr**6
rprod2=rprod1*rprod1
pot0=pot0+ewca*(4.0d0*(sigmap12/rprod2-sigmap6/rprod1)+(absr-cutwca)*frc-urc)
end if
end do
end do

!double well forces for first two particles
xr=x0(1,:)-x0(2,:)
!!periodic boundaries
xr(1)=xr(1)-nint(xr(1)/L)*L
xr(2)=xr(2)-nint(xr(2)/L)*L
absr=dot_product(xr,xr)**0.5
fac1=absr-cutwca-wdw
fac2=1.0d0-(fac1)**2/wdw/wdw
pot0=pot0+hdw*fac2*fac2

end subroutine potential
!*************************************************************************!
!*************************************************************************!
!*************************************************************************!
!*************************************************************************!


! { dg-do compile }
!
! PR fortran/38718
!
  implicit none
  real(kind=8), parameter :: r = kind(0) + 0.2
  complex(kind=8), parameter :: c = (r, -9.3)
  integer, parameter :: k = nint(dreal(c))
  integer, parameter :: l = nint(realpart(c))
  integer(kind=k) :: i
  integer(kind=l) :: j
  i = 42
  j = 42
  print *, k, i, j, r
  end


text:
subroutine Xsections_general(species,reactions_rx,kj,nz,nw,kw,wavl,T,jn,sq)
  implicit none

  ! module variables

  ! inputs
  character(len=8), intent(in) :: species
  character(len=*), intent(in) :: reactions_rx
  integer, intent(in) :: kj, nz, nw, kw
  double precision, dimension(nw+1), intent(in) :: wavl
  double precision, dimension(nz), intent(in) :: T

  ! in/out
  integer, intent(inout) :: jn
  double precision, dimension(kj,nz,kw), intent(inout) :: sq

  ! other vars
  character(len=50) :: fmt
  character(len=10) :: reac1, reac2, prod1, prod2, prod3
  character(len=25) :: reacs, prods
  character(len=5) :: label
  character(len=500) :: xsname, qyname
  character(len=1000) :: temperatureline, line
  character(len=8), dimension(100) :: temperatures
  character(len=8) :: temp_dumb
  integer :: io, stat, numtempcols, kdata, i, ierr, iw
  double precision, allocatable, dimension(:,:) :: yy, yy_grid
  double precision, allocatable, dimension(:) :: temperature_nums
  double precision, allocatable, dimension(:) :: x1, y1, x2, qy
  double precision, allocatable, dimension(:) :: dumby
  double precision :: deltax = 1.d-4, biggest=1.d+36, zero=0.0d0
  double precision, dimension(nw) :: yg1, yq1
  double precision, dimension(nw,nz) :: yg2
  double precision, dimension(nz) :: p_interp

  character(len=500) :: rootdir
  rootdir = '../PhotochemPy/'

  ! initialize
  ierr = 0

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!!! Read in cross sections !!!!
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  ! open xs file
  xsname = trim(rootdir)//'DATA/XSECTIONS/'//trim(species)//'/'// &
            trim(species)//'.XS.dat'
  open(11, file=trim(xsname), status='OLD')

  ! determine number of temperature columns
  read(11,*)
  read(11,*)
  read(11,'(A)') temperatureline
  do i=1,100
    read(temperatureline,*,iostat=io) temperatures(1:i)
    if (io == -1) exit
  enddo
  if (i == 101) then
    print*,'More temperature data than allowed for ',species
  endif
  ! now i-2 is number of temperature columns
  numtempcols = max(i-2,1)

  ! determine how many xsections points there are
  read(11,*)
  read(11,*)
  io = 0
  kdata = 0
  do while(io == 0)
    kdata = kdata + 1
    READ(11, *, IOSTAT=io)
  enddo

  ! allocate memory
  allocate(yy(kdata+4,numtempcols))
  allocate(yy_grid(nw,numtempcols+2))
  allocate(temperature_nums(numtempcols+2))
  allocate(y1(kdata+4))
  allocate(x1(kdata+4))
  allocate(dumby(numtempcols+1))

  ! now read in data
  rewind(11)
  do i=1,4
    read(11,*)
  enddo
  do i=1,kdata
    read(11,'(A)') line
    read(line,*) dumby(1:numtempcols+1)
    x1(i) = dumby(1)
    yy(i,:) = dumby(2:)
  enddo
  close(11)

  ! re-bin data to the grid
  y1 = 0.d0
  y1 = yy(:,1)
  call addpnt(x1,y1,kdata+4,kdata,x1(1)*(1.d0-deltax),zero)
  call addpnt(x1,y1,kdata+4,kdata,               zero,zero)
  call addpnt(x1,y1,kdata+4,kdata,x1(kdata)*(1.d0+deltax),zero)
  call addpnt(x1,y1,kdata+4,kdata,            biggest,zero)
  call inter2(nw+1,wavl,yg1,kdata,x1,y1,ierr)
  kdata = kdata - 4
  yy_grid(:,1) = yg1
  do i=1,numtempcols
    x1(:kdata) = x1(3:kdata+2)
    x1(kdata+1:kdata+4) = 0.d0
    y1 = 0.d0
    y1 = yy(:,i)
    call addpnt(x1,y1,kdata+4,kdata,x1(1)*(1.d0-deltax),zero)
    call addpnt(x1,y1,kdata+4,kdata,               zero,zero)
    call addpnt(x1,y1,kdata+4,kdata,x1(kdata)*(1.d0+deltax),zero)
    call addpnt(x1,y1,kdata+4,kdata,            biggest,zero)
    call inter2(nw+1,wavl,yg1,kdata,x1,y1,ierr)
    kdata = kdata-4
    yy_grid(:,i+1) = yg1
  enddo
  x1(:kdata) = x1(3:kdata+2)
  x1(kdata+1:kdata+4) = 0.d0
  y1 = 0.d0
  y1 = yy(:,numtempcols)
  call addpnt(x1,y1,kdata+4,kdata,x1(1)*(1.d0-deltax),zero)
  call addpnt(x1,y1,kdata+4,kdata,               zero,zero)
  call addpnt(x1,y1,kdata+4,kdata,x1(kdata)*(1.d0+deltax),zero)
  call addpnt(x1,y1,kdata+4,kdata,            biggest,zero)
  call inter2(nw+1,wavl,yg1,kdata,x1,y1,ierr)
  kdata = kdata-4
  yy_grid(:,numtempcols+2) = yg1
  if (ierr .NE. 0) then
    print*, "failed while binning cross sections for ", species
  endif

  ! find temperatures of each column
  if (numtempcols == 1) then
    temperature_nums(2) = 300.d0 ! if 1 temp column then assume its 300K
  else
    do i=2,numtempcols+1
      temp_dumb = temperatures(i)
      read(temp_dumb(:index(temp_dumb,'K')-1),*) temperature_nums(i)
    enddo
  endif
  temperature_nums(1) = 0.d0
  temperature_nums(numtempcols+2) = 3000.d0
  do i=1,nz
    if ((T(i) > 3000.d0) .or. (T(i) < 0.d0)) then
      print*, 'Atmospheric temperature can not be above 3000 K or below 0 K'
      stop
    endif
  enddo

  ! now interpolate to atmospheric temperature for each bin
  ! and at each temperature
  do i = 1, nw
    call interp_linear(1, numtempcols+2, temperature_nums, yy_grid(i,:), nz, &
                       T, p_interp)
    yg2(i,:) = p_interp
  enddo

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!!!!!! Quantum Yields !!!!!!!
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  ! should use this intead of reading reactions again
  ! j is the index of sq(j,i,iw)
  ! spec_ind = index(ispec,species)
  ! do j=1,kj
  !   if (photoreac(j) == spec_ind) then
  !     reac1 = chemj(1,photonums(j))
  !     reac2 = chemj(2,photonums(j))
  !     prod1 = chemj(3,photonums(j))
  !     prod2 = chemj(4,photonums(j))
  !     prod3 = chemj(5,photonums(j))
  !     print*,reac1,species
  !   endif
  ! enddo

  fmt = '(A10,A10,A10,A10,A8,A5)'
  open(12, file=trim(reactions_rx),status='OLD')
  io = 0
  do while(io == 0)
    read(12,trim(fmt),IOSTAT=io) reac1,reac2,prod1,prod2,prod3,label

    if ((label.eq.'PHOTO').and.(trim(species).eq.trim(reac1))) then

      reacs = trim(reac1) // " " // trim(reac2)
      prods = trim(prod1) // " " // trim(prod2) // " " &
             // trim(prod3)

      do i=1, len(trim(reacs))
        if(reacs(i:i) == ' ') reacs(i:i) = '_'
      end do

      do i=1, len(trim(prods))
        if(prods(i:i) == ' ') prods(i:i) = '_'
      enddo

      ! open QY file
      qyname = trim(rootdir)//'DATA/XSECTIONS/'//trim(species)//'/'// &
                trim(reacs)//"_"//trim(prods)//".QY.dat"

      ! open qy
      open(13, file=trim(qyname), status='OLD')
      ! count lines
      kdata = 0
      stat = 0
      do i=1,5
        read(13,*)
      enddo
      do while(stat == 0)
        read(13, *, IOSTAT=stat)
        kdata = kdata + 1
      end do

      allocate(qy(kdata+4))
      allocate(x2(kdata+4))

      ! rewind and read in data
      rewind(13)
      do i=1,4
        read(13,*)
      enddo
      do i=1,kdata
        read(13,*) x2(i), qy(i)
      enddo
      close(13)

      ! re-bin data to grid
      call addpnt(x2,qy,kdata+4,kdata,x2(1)*(1.d0-deltax),zero)
      call addpnt(x2,qy,kdata+4,kdata,               zero,zero)
      call addpnt(x2,qy,kdata+4,kdata,x2(kdata)*(1.d0+deltax),zero)
      call addpnt(x2,qy,kdata+4,kdata,            biggest,zero)
      call inter2(nw+1,wavl,yq1,kdata,x2,qy,ierr)
      if (ierr .ne. 0) THEN
        print*, "failed while reading in quantum yields for ", species
      endif

      do iw = 1, nw
        do i = 1, nz
          sq(jn,i,iw) = yg2(iw,i)*yq1(iw)
        enddo
      enddo

      jn = jn+1

      deallocate(qy)
      deallocate(x2)

    endif
  enddo
  close(12)

  ! deallocate
  ! all should auto-deallocate after going out of scope.
  deallocate(yy,yy_grid)
  deallocate(temperature_nums)
  deallocate(x1, y1)
  deallocate(dumby)

end subroutine

max_stars_repo_name: Nicholaswogan/PhotochemPy

max_issues_repo_name: StAtmos/PhotochemPy

max_forks_repo_name: StAtmos/PhotochemPy


text:
subroutine blbopo(a)
   use typre
   use Mod_Mesh
   implicit none
   class(FemMesh) :: a

   integer(ip) :: iaux,iaux2,iboun,ielem,inodb,inode,ipoin,ispos,ispos2,pblty

   call a%Timer%Blbopo%Tic

   pblty=1
   call a%Memor%alloc(a%npoin+1,a%pbopo,'pbopo','blbopo')
   !first loop to count
   do iboun = 1,a%nboun
      ispos2 = (a%nnodb(1)+1)*(iboun-1)
      if (a%nelty > 1) then
         ispos2 = a%pboel(iboun)-1
         pblty = a%ltypb(a%pboel(iboun+1)-a%pboel(iboun)-1)
      endif
      ielem = a%lboel(ispos2+a%nnodb(pblty)+1)
      ispos = a%nnode(1)*(ielem-1)
      if (a%nelty > 1) ispos = a%pnods(ielem)-1
      do inodb = 1,a%nnodb(pblty)
         inode = a%lboel(ispos2+inodb)
         ipoin = a%lnods(ispos+inode)
         a%pbopo(ipoin) = a%pbopo(ipoin)+1
      enddo
   enddo

   !second loop to fill
   iaux = a%pbopo(1)
   a%pbopo(1) = 1
   do ipoin = 2,a%npoin
      iaux2 = a%pbopo(ipoin)
      a%pbopo(ipoin) = a%pbopo(ipoin-1) + iaux
      iaux = iaux2
   enddo
   if (a%npoin > 0) then
      a%pbopo(a%npoin+1) = a%pbopo(a%npoin) +iaux
   endif

   call a%Memor%alloc(a%pbopo(a%npoin+1),a%lbopo,'lbopo','blbopo')

   do iboun = 1,a%nboun
      ispos2 = (a%nnodb(1)+1)*(iboun-1)
      if (a%nelty > 1) then
         ispos2 = a%pboel(iboun)-1
         pblty = a%ltypb(a%pboel(iboun+1)-a%pboel(iboun)-1)
      endif
      ielem = a%lboel(ispos2+a%nnodb(pblty)+1)
      ispos = a%nnode(1)*(ielem-1)
      if (a%nelty > 1) ispos = a%pnods(ielem)-1
      do inodb = 1,a%nnodb(pblty)
         inode = a%lboel(ispos2+inodb)
         ipoin = a%lnods(ispos+inode)
         a%lbopo(a%pbopo(ipoin)) = iboun
         a%pbopo(ipoin) = a%pbopo(ipoin)+1
      enddo
   enddo
   do ipoin = a%npoin,2,-1
      a%pbopo(ipoin) = a%pbopo(ipoin-1)
   enddo
   a%pbopo(1) = 1

   call a%Timer%Blbopo%Toc

end subroutine

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
SUBROUTINE CFFT( DATA, N, IFORW )

  ! complex FFT

  IMPLICIT NONE
  INTEGER I1, I2A, I2B, I3, I3Rev, IP1, IP2, N, ISign, IForw
  REAL    theta, sinth
  COMPLEX DATA( N ), TEMP, W, WSTP

  ISIGN = -IFORW
  I3REV = 1

  DO I3 = 1, N
     IF ( I3 < I3REV ) THEN   ! switch values
        TEMP          = DATA( I3 )
        DATA( I3 )    = DATA( I3REV )
        DATA( I3REV ) = TEMP
     ENDIF
     ! following loop is just to compute I3REV
     IP1 = N / 2
     DO WHILE (  I3REV > IP1 )
        IF ( IP1 <= 1 ) EXIT
        I3REV = I3REV - IP1
        IP1   = IP1 / 2
     END DO
     I3REV = I3REV + IP1
  END DO

  IP1 = 1

  DO WHILE  ( IP1 < N )
     IP2   = IP1 * 2
     THETA = 6.283185307 / FLOAT( ISIGN * IP2 )
     SINTH = SIN( THETA / 2.)
     WSTP  = CMPLX( -2. * SINTH * SINTH, SIN( THETA ) )
     W     = 1.

     DO I1 = 1, IP1
        DO I3 = I1, N, IP2
           I2A = I3
           I2B = I2A + IP1
           TEMP      = W * DATA( I2B )
           DATA( I2B ) = DATA( I2A ) - TEMP
           DATA( I2A ) = DATA( I2A ) + TEMP
        END DO

        W = W * WSTP + W
     END DO

     IP1 = IP2
  END DO

  RETURN
END SUBROUTINE CFFT

max_stars_repo_name: IvanaEscobar/sandbox

max_issues_repo_name: IvanaEscobar/sandbox

max_forks_repo_licenses: Naumen Condor-1.1 MS-PL


text:
*
*-----------------------------------------------------------------------
*	A C U R A T
*-----------------------------------------------------------------------
*
*	Coefficients of Slater integrals are square roots of
*   rational numbers.  To improve the accuaracy, certain commonly
*   occuring coefficients are improved to machine accuracy.
*
	DOUBLE PRECISION FUNCTION ACURAT(C)
        IMPLICIT DOUBLE PRECISION(A-H,O-Z)
        INTEGER NUM
        INTEGER DEN(11)
        DATA DEN/2,3,7,9,15,35,49,175,189,315,441/
        DATA D1/1.D0/
*
        C2 = C*C
        ACURAT = C
        DO 1 I = 1,11
           PROD = DEN(I)*C2
           NUM = NINT(PROD)
           EPS = DABS(NUM-PROD)/DEN(I)
           IF (EPS .LE. 1.D-8) THEN
              IF (EPS .NE. 0.) THEN
                 ACURAT = DSQRT((NUM*D1)/DEN(I))
                 IF (C .LT. 0.) ACURAT = -ACURAT
                 RETURN
              END IF
           END IF
  1     CONTINUE
        END


text:
      PROGRAM WMEX04
C
C  Examples of parameter control of fronts.
C
C Define the error file, the Fortran unit number, the workstation type,
C and the workstation ID to be used in calls to GKS routines.
C
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=1,  IWKID=1)   ! NCGM
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=8,  IWKID=1)   ! X Windows
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=11, IWKID=1)   ! PDF
C     PARAMETER (IERRF=6, LUNIT=2, IWTYPE=20, IWKID=1)   ! PostScript
C
      PARAMETER (IERRF=6, LUNIT=2, IWTYPE=1,  IWKID=1)
C
      PARAMETER (NS=2)
      DIMENSION XS(NS),YS(NS)
      DATA XS/ 0.10, 0.90/
C
C  Open GKS, open and activate a workstation.
C
      CALL GOPKS (IERRF, ISZDM)
      CALL GOPWK (IWKID, LUNIT, IWTYPE)
      CALL GACWK (IWKID)
C
C  Define a color table.
C
      CALL GSCR(IWKID, 0, 1.0, 1.0, 1.0)
      CALL GSCR(IWKID, 1, 0.0, 0.0, 0.0)
      CALL GSCR(IWKID, 2, 1.0, 0.0, 0.0)
      CALL GSCR(IWKID, 3, 0.0, 0.0, 1.0)
      CALL GSCR(IWKID, 4, 0.4, 0.0, 0.4)
C
C  Plot title.
C
      CALL PLCHHQ(0.50,0.94,
     + ':F26:Parameter control of front attributes', 0.03,0.,0.)
C
C  Various fronts with different attributes.
C
      CALL PCSETI('CC',4)
      FSIZE = .021
      YS(1) = .76
      YS(2) = .76
      CALL PCSETC('FC','%')
      CALL PLCHHQ(XS(1),YS(1)+.06,'%F22%Starting with FRO=''STA'', WFC=2
     + (red), CFC=3 (blue):',FSIZE,0.,-1.)
      CALL WMSETC('FRO','STA')
      CALL WMSETI('WFC',2)
      CALL WMSETI('CFC',3)
      CALL WMDRFT(NS,XS,YS)
C
      YS(1) = .60
      YS(2) = .60
      CALL PLCHHQ(XS(1),YS(1)+.06,'%F22%then setting BEG=0., END=.05, BE
     +T=.03 gives:',FSIZE,0.,-1.)
      CALL WMSETR('BEG',0.00)
      CALL WMSETR('END',0.05)
      CALL WMSETR('BET',0.03)
      CALL WMDRFT(NS,XS,YS)
C
      YS(1) = .44
      YS(2) = .44
      CALL PLCHHQ(XS(1),YS(1)+.06,'%F22%then setting NMS=5 and STY=-1,2,
     +1,-2,2 gives:',FSIZE,0.,-1.)
      CALL WMSETI('NMS',5)
      CALL WMSETI('PAI',1)
      CALL WMSETI('STY',-1)
      CALL WMSETI('PAI',2)
      CALL WMSETI('STY',2)
      CALL WMSETI('PAI',3)
      CALL WMSETI('STY',1)
      CALL WMSETI('PAI',4)
      CALL WMSETI('STY',-2)
      CALL WMSETI('PAI',5)
      CALL WMSETI('STY',2)
      CALL WMDRFT(NS,XS,YS)
C
      YS(1) = .27
      YS(2) = .27
      CALL PLCHHQ(XS(1),YS(1)+.07,'%F22%then setting SWI=.05 and LIN=12.
     + gives:',FSIZE,0.,-1.)
      CALL WMSETR('SWI',0.05)
      CALL WMSETR('LIN',12.)
      CALL WMDRFT(NS,XS,YS)
C
      YS(1) = .10
      YS(2) = .10
      CALL PLCHHQ(XS(1),YS(1)+.07,'%F22%then setting REV=1 gives:',
     +            .022,0.,-1.)
      CALL WMSETI('REV',1)
      CALL WMDRFT(NS,XS,YS)
      CALL FRAME
C
      CALL GDAWK(IWKID)
      CALL GCLWK(IWKID)
      CALL GCLKS
      STOP
C
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl