text:
*DECK XERMSG
      SUBROUTINE XERMSG (LIBRAR, SUBROU, MESSG, NERR, LEVEL)
C***BEGIN PROLOGUE  XERMSG
C***PURPOSE  Process error messages for SLATEC and other libraries.
C***LIBRARY   SLATEC (XERROR)
C***CATEGORY  R3C
C***TYPE      ALL (XERMSG-A)
C***KEYWORDS  ERROR MESSAGE, XERROR
C***AUTHOR  Fong, Kirby, (NMFECC at LLNL)
C***DESCRIPTION
C
C   XERMSG processes a diagnostic message in a manner determined by the
C   value of LEVEL and the current value of the library error control
C   flag, KONTRL.  See subroutine XSETF for details.
C
C    LIBRAR   A character constant (or character variable) with the name
C             of the library.  This will be 'SLATEC' for the SLATEC
C             Common Math Library.  The error handling package is
C             general enough to be used by many libraries
C             simultaneously, so it is desirable for the routine that
C             detects and reports an error to identify the library name
C             as well as the routine name.
C
C    SUBROU   A character constant (or character variable) with the name
C             of the routine that detected the error.  Usually it is the
C             name of the routine that is calling XERMSG.  There are
C             some instances where a user callable library routine calls
C             lower level subsidiary routines where the error is
C             detected.  In such cases it may be more informative to
C             supply the name of the routine the user called rather than
C             the name of the subsidiary routine that detected the
C             error.
C
C    MESSG    A character constant (or character variable) with the text
C             of the error or warning message.  In the example below,
C             the message is a character constant that contains a
C             generic message.
C
C                   CALL XERMSG ('SLATEC', 'MMPY',
C                  *'THE ORDER OF THE MATRIX EXCEEDS THE ROW DIMENSION',
C                  *3, 1)
C
C             It is possible (and is sometimes desirable) to generate a
C             specific message--e.g., one that contains actual numeric
C             values.  Specific numeric values can be converted into
C             character strings using formatted WRITE statements into
C             character variables.  This is called standard Fortran
C             internal file I/O and is exemplified in the first three
C             lines of the following example.  You can also catenate
C             substrings of characters to construct the error message.
C             Here is an example showing the use of both writing to
C             an internal file and catenating character strings.
C
C                   CHARACTER*5 CHARN, CHARL
C                   WRITE (CHARN,10) N
C                   WRITE (CHARL,10) LDA
C                10 FORMAT(I5)
C                   CALL XERMSG ('SLATEC', 'MMPY', 'THE ORDER'//CHARN//
C                  *   ' OF THE MATRIX EXCEEDS ITS ROW DIMENSION OF'//
C                  *   CHARL, 3, 1)
C
C             There are two subtleties worth mentioning.  One is that
C             the // for character catenation is used to construct the
C             error message so that no single character constant is
C             continued to the next line.  This avoids confusion as to
C             whether there are trailing blanks at the end of the line.
C             The second is that by catenating the parts of the message
C             as an actual argument rather than encoding the entire
C             message into one large character variable, we avoid
C             having to know how long the message will be in order to
C             declare an adequate length for that large character
C             variable.  XERMSG calls XERPRN to print the message using
C             multiple lines if necessary.  If the message is very long,
C             XERPRN will break it into pieces of 72 characters (as
C             requested by XERMSG) for printing on multiple lines.
C             Also, XERMSG asks XERPRN to prefix each line with ' *  '
C             so that the total line length could be 76 characters.
C             Note also that XERPRN scans the error message backwards
C             to ignore trailing blanks.  Another feature is that
C             the substring '$$' is treated as a new line sentinel
C             by XERPRN.  If you want to construct a multiline
C             message without having to count out multiples of 72
C             characters, just use '$$' as a separator.  '$$'
C             obviously must occur within 72 characters of the
C             start of each line to have its intended effect since
C             XERPRN is asked to wrap around at 72 characters in
C             addition to looking for '$$'.
C
C    NERR     An integer value that is chosen by the library routine's
C             author.  It must be in the range -99 to 999 (three
C             printable digits).  Each distinct error should have its
C             own error number.  These error numbers should be described
C             in the machine readable documentation for the routine.
C             The error numbers need be unique only within each routine,
C             so it is reasonable for each routine to start enumerating
C             errors from 1 and proceeding to the next integer.
C
C    LEVEL    An integer value in the range 0 to 2 that indicates the
C             level (severity) of the error.  Their meanings are
C
C            -1  A warning message.  This is used if it is not clear
C                that there really is an error, but the user's attention
C                may be needed.  An attempt is made to only print this
C                message once.
C
C             0  A warning message.  This is used if it is not clear
C                that there really is an error, but the user's attention
C                may be needed.
C
C             1  A recoverable error.  This is used even if the error is
C                so serious that the routine cannot return any useful
C                answer.  If the user has told the error package to
C                return after recoverable errors, then XERMSG will
C                return to the Library routine which can then return to
C                the user's routine.  The user may also permit the error
C                package to terminate the program upon encountering a
C                recoverable error.
C
C             2  A fatal error.  XERMSG will not return to its caller
C                after it receives a fatal error.  This level should
C                hardly ever be used; it is much better to allow the
C                user a chance to recover.  An example of one of the few
C                cases in which it is permissible to declare a level 2
C                error is a reverse communication Library routine that
C                is likely to be called repeatedly until it integrates
C                across some interval.  If there is a serious error in
C                the input such that another step cannot be taken and
C                the Library routine is called again without the input
C                error having been corrected by the caller, the Library
C                routine will probably be called forever with improper
C                input.  In this case, it is reasonable to declare the
C                error to be fatal.
C
C    Each of the arguments to XERMSG is input; none will be modified by
C    XERMSG.  A routine may make multiple calls to XERMSG with warning
C    level messages; however, after a call to XERMSG with a recoverable
C    error, the routine should return to the user.  Do not try to call
C    XERMSG with a second recoverable error after the first recoverable
C    error because the error package saves the error number.  The user
C    can retrieve this error number by calling another entry point in
C    the error handling package and then clear the error number when
C    recovering from the error.  Calling XERMSG in succession causes the
C    old error number to be overwritten by the latest error number.
C    This is considered harmless for error numbers associated with
C    warning messages but must not be done for error numbers of serious
C    errors.  After a call to XERMSG with a recoverable error, the user
C    must be given a chance to call NUMXER or XERCLR to retrieve or
C    clear the error number.
C***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
C                 Error-handling Package, SAND82-0800, Sandia
C                 Laboratories, 1982.
C***ROUTINES CALLED  FDUMP, J4SAVE, XERCNT, XERHLT, XERPRN, XERSVE
C***REVISION HISTORY  (YYMMDD)
C   880101  DATE WRITTEN
C   880621  REVISED AS DIRECTED AT SLATEC CML MEETING OF FEBRUARY 1988.
C           THERE ARE TWO BASIC CHANGES.
C           1.  A NEW ROUTINE, XERPRN, IS USED INSTEAD OF XERPRT TO
C               PRINT MESSAGES.  THIS ROUTINE WILL BREAK LONG MESSAGES
C               INTO PIECES FOR PRINTING ON MULTIPLE LINES.  '$$' IS
C               ACCEPTED AS A NEW LINE SENTINEL.  A PREFIX CAN BE
C               ADDED TO EACH LINE TO BE PRINTED.  XERMSG USES EITHER
C               ' ***' OR ' *  ' AND LONG MESSAGES ARE BROKEN EVERY
C               72 CHARACTERS (AT MOST) SO THAT THE MAXIMUM LINE
C               LENGTH OUTPUT CAN NOW BE AS GREAT AS 76.
C           2.  THE TEXT OF ALL MESSAGES IS NOW IN UPPER CASE SINCE THE
C               FORTRAN STANDARD DOCUMENT DOES NOT ADMIT THE EXISTENCE
C               OF LOWER CASE.
C   880708  REVISED AFTER THE SLATEC CML MEETING OF JUNE 29 AND 30.
C           THE PRINCIPAL CHANGES ARE
C           1.  CLARIFY COMMENTS IN THE PROLOGUES
C           2.  RENAME XRPRNT TO XERPRN
C           3.  REWORK HANDLING OF '$$' IN XERPRN TO HANDLE BLANK LINES
C               SIMILAR TO THE WAY FORMAT STATEMENTS HANDLE THE /
C               CHARACTER FOR NEW RECORDS.
C   890706  REVISED WITH THE HELP OF FRED FRITSCH AND REG CLEMENS TO
C           CLEAN UP THE CODING.
C   890721  REVISED TO USE NEW FEATURE IN XERPRN TO COUNT CHARACTERS IN
C           PREFIX.
C   891013  REVISED TO CORRECT COMMENTS.
C   891214  Prologue converted to Version 4.0 format.  (WRB)
C   900510  Changed test on NERR to be -9999999 < NERR < 99999999, but
C           NERR .ne. 0, and on LEVEL to be -2 < LEVEL < 3.  Added
C           LEVEL=-1 logic, changed calls to XERSAV to XERSVE, and
C           XERCTL to XERCNT.  (RWC)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  XERMSG
      CHARACTER*(*) LIBRAR, SUBROU, MESSG
      CHARACTER*8 XLIBR, XSUBR
      CHARACTER*72  TEMP
      CHARACTER*20  LFIRST
C***FIRST EXECUTABLE STATEMENT  XERMSG
      LKNTRL = J4SAVE (2, 0, .FALSE.)
      MAXMES = J4SAVE (4, 0, .FALSE.)
C
C       LKNTRL IS A LOCAL COPY OF THE CONTROL FLAG KONTRL.
C       MAXMES IS THE MAXIMUM NUMBER OF TIMES ANY PARTICULAR MESSAGE
C          SHOULD BE PRINTED.
C
C       WE PRINT A FATAL ERROR MESSAGE AND TERMINATE FOR AN ERROR IN
C          CALLING XERMSG.  THE ERROR NUMBER SHOULD BE POSITIVE,
C          AND THE LEVEL SHOULD BE BETWEEN 0 AND 2.
C
      IF (NERR.LT.-9999999 .OR. NERR.GT.99999999 .OR. NERR.EQ.0 .OR.
     *   LEVEL.LT.-1 .OR. LEVEL.GT.2) THEN
         CALL XERPRN (' ***', -1, 'FATAL ERROR IN...$$ ' //
     *      'XERMSG -- INVALID ERROR NUMBER OR LEVEL$$ '//
     *      'JOB ABORT DUE TO FATAL ERROR.', 72)
         CALL XERSVE (' ', ' ', ' ', 0, 0, 0, KDUMMY)
         CALL XERHLT (' ***XERMSG -- INVALID INPUT')
         RETURN
      ENDIF
C
C       RECORD THE MESSAGE.
C
      I = J4SAVE (1, NERR, .TRUE.)
      CALL XERSVE (LIBRAR, SUBROU, MESSG, 1, NERR, LEVEL, KOUNT)
C
C       HANDLE PRINT-ONCE WARNING MESSAGES.
C
      IF (LEVEL.EQ.-1 .AND. KOUNT.GT.1) RETURN
C
C       ALLOW TEMPORARY USER OVERRIDE OF THE CONTROL FLAG.
C
      XLIBR  = LIBRAR
      XSUBR  = SUBROU
      LFIRST = MESSG
      LERR   = NERR
      LLEVEL = LEVEL
      CALL XERCNT (XLIBR, XSUBR, LFIRST, LERR, LLEVEL, LKNTRL)
C
      LKNTRL = MAX(-2, MIN(2,LKNTRL))
      MKNTRL = ABS(LKNTRL)
C
C       SKIP PRINTING IF THE CONTROL FLAG VALUE AS RESET IN XERCNT IS
C       ZERO AND THE ERROR IS NOT FATAL.
C
      IF (LEVEL.LT.2 .AND. LKNTRL.EQ.0) GO TO 30
      IF (LEVEL.EQ.0 .AND. KOUNT.GT.MAXMES) GO TO 30
      IF (LEVEL.EQ.1 .AND. KOUNT.GT.MAXMES .AND. MKNTRL.EQ.1) GO TO 30
      IF (LEVEL.EQ.2 .AND. KOUNT.GT.MAX(1,MAXMES)) GO TO 30
C
C       ANNOUNCE THE NAMES OF THE LIBRARY AND SUBROUTINE BY BUILDING A
C       MESSAGE IN CHARACTER VARIABLE TEMP (NOT EXCEEDING 66 CHARACTERS)
C       AND SENDING IT OUT VIA XERPRN.  PRINT ONLY IF CONTROL FLAG
C       IS NOT ZERO.
C
      IF (LKNTRL .NE. 0) THEN
         TEMP(1:21) = 'MESSAGE FROM ROUTINE '
         I = MIN(LEN(SUBROU), 16)
         TEMP(22:21+I) = SUBROU(1:I)
         TEMP(22+I:33+I) = ' IN LIBRARY '
         LTEMP = 33 + I
         I = MIN(LEN(LIBRAR), 16)
         TEMP(LTEMP+1:LTEMP+I) = LIBRAR (1:I)
         TEMP(LTEMP+I+1:LTEMP+I+1) = '.'
         LTEMP = LTEMP + I + 1
         CALL XERPRN (' ***', -1, TEMP(1:LTEMP), 72)
      ENDIF
C
C       IF LKNTRL IS POSITIVE, PRINT AN INTRODUCTORY LINE BEFORE
C       PRINTING THE MESSAGE.  THE INTRODUCTORY LINE TELLS THE CHOICE
C       FROM EACH OF THE FOLLOWING THREE OPTIONS.
C       1.  LEVEL OF THE MESSAGE
C              'INFORMATIVE MESSAGE'
C              'POTENTIALLY RECOVERABLE ERROR'
C              'FATAL ERROR'
C       2.  WHETHER CONTROL FLAG WILL ALLOW PROGRAM TO CONTINUE
C              'PROG CONTINUES'
C              'PROG ABORTED'
C       3.  WHETHER OR NOT A TRACEBACK WAS REQUESTED.  (THE TRACEBACK
C           MAY NOT BE IMPLEMENTED AT SOME SITES, SO THIS ONLY TELLS
C           WHAT WAS REQUESTED, NOT WHAT WAS DELIVERED.)
C              'TRACEBACK REQUESTED'
C              'TRACEBACK NOT REQUESTED'
C       NOTICE THAT THE LINE INCLUDING FOUR PREFIX CHARACTERS WILL NOT
C       EXCEED 74 CHARACTERS.
C       WE SKIP THE NEXT BLOCK IF THE INTRODUCTORY LINE IS NOT NEEDED.
C
      IF (LKNTRL .GT. 0) THEN
C
C       THE FIRST PART OF THE MESSAGE TELLS ABOUT THE LEVEL.
C
         IF (LEVEL .LE. 0) THEN
            TEMP(1:20) = 'INFORMATIVE MESSAGE,'
            LTEMP = 20
         ELSEIF (LEVEL .EQ. 1) THEN
            TEMP(1:30) = 'POTENTIALLY RECOVERABLE ERROR,'
            LTEMP = 30
         ELSE
            TEMP(1:12) = 'FATAL ERROR,'
            LTEMP = 12
         ENDIF
C
C       THEN WHETHER THE PROGRAM WILL CONTINUE.
C
         IF ((MKNTRL.EQ.2 .AND. LEVEL.GE.1) .OR.
     *       (MKNTRL.EQ.1 .AND. LEVEL.EQ.2)) THEN
            TEMP(LTEMP+1:LTEMP+14) = ' PROG ABORTED,'
            LTEMP = LTEMP + 14
         ELSE
            TEMP(LTEMP+1:LTEMP+16) = ' PROG CONTINUES,'
            LTEMP = LTEMP + 16
         ENDIF
C
C       FINALLY TELL WHETHER THERE SHOULD BE A TRACEBACK.
C
         IF (LKNTRL .GT. 0) THEN
            TEMP(LTEMP+1:LTEMP+20) = ' TRACEBACK REQUESTED'
            LTEMP = LTEMP + 20
         ELSE
            TEMP(LTEMP+1:LTEMP+24) = ' TRACEBACK NOT REQUESTED'
            LTEMP = LTEMP + 24
         ENDIF
         CALL XERPRN (' ***', -1, TEMP(1:LTEMP), 72)
      ENDIF
C
C       NOW SEND OUT THE MESSAGE.
C
      CALL XERPRN (' *  ', -1, MESSG, 72)
C
C       IF LKNTRL IS POSITIVE, WRITE THE ERROR NUMBER AND REQUEST A
C          TRACEBACK.
C
      IF (LKNTRL .GT. 0) THEN
         WRITE (TEMP, '(''ERROR NUMBER = '', I8)') NERR
         DO 10 I=16,22
            IF (TEMP(I:I) .NE. ' ') GO TO 20
   10    CONTINUE
C
   20    CALL XERPRN (' *  ', -1, TEMP(1:15) // TEMP(I:23), 72)
         CALL FDUMP
      ENDIF
C
C       IF LKNTRL IS NOT ZERO, PRINT A BLANK LINE AND AN END OF MESSAGE.
C
      IF (LKNTRL .NE. 0) THEN
         CALL XERPRN (' *  ', -1, ' ', 72)
         CALL XERPRN (' ***', -1, 'END OF MESSAGE', 72)
         CALL XERPRN ('    ',  0, ' ', 72)
      ENDIF
C
C       IF THE ERROR IS NOT FATAL OR THE ERROR IS RECOVERABLE AND THE
C       CONTROL FLAG IS SET FOR RECOVERY, THEN RETURN.
C
   30 IF (LEVEL.LE.0 .OR. (LEVEL.EQ.1 .AND. MKNTRL.LE.1)) RETURN
C
C       THE PROGRAM WILL BE STOPPED DUE TO AN UNRECOVERED ERROR OR A
C       FATAL ERROR.  PRINT THE REASON FOR THE ABORT AND THE ERROR
C       SUMMARY IF THE CONTROL FLAG AND THE MAXIMUM ERROR COUNT PERMIT.
C
      IF (LKNTRL.GT.0 .AND. KOUNT.LT.MAX(1,MAXMES)) THEN
         IF (LEVEL .EQ. 1) THEN
            CALL XERPRN
     *         (' ***', -1, 'JOB ABORT DUE TO UNRECOVERED ERROR.', 72)
         ELSE
            CALL XERPRN(' ***', -1, 'JOB ABORT DUE TO FATAL ERROR.', 72)
         ENDIF
         CALL XERSVE (' ', ' ', ' ', -1, 0, 0, KDUMMY)
         CALL XERHLT (' ')
      ELSE
         CALL XERHLT (MESSG)
      ENDIF
      RETURN
      END


text:
c-----------------------------------------------------------------------
      subroutine conaa(iret)
c-----------------------------------------------------------------------
c  determines interaction configuration
c-----------------------------------------------------------------------

      include 'epos.inc'
      include 'epos.incems'
      include 'epos.incsem'
      include 'epos.incpar'

      common/geom/rmproj,rmtarg,bmax,bkmx
      common/nucl3/phi,bimp
      common /cncl/xproj(mamx),yproj(mamx),zproj(mamx)
     *            ,xtarg(mamx),ytarg(mamx),ztarg(mamx)
      common/cfacmss/facmss
      double precision yyrmax
      common/scrangle/ phik3(kollmx),thetak3(kollmx)

      call utpri('conaa ',ish,ishini,4)

      iret=0

c     initialisations
c     ---------------

      vel=tanh(ypjtl-yhaha)+tanh(yhaha)

c     determine phi, bimp, coll, iproj, itarg, x/y/zproj, x/y/ztarg
c     ---------------------------------------------------------------

           if(iokoll.eq.1)then

      koll=matarg
      do k=1,koll
        do n=1,4
          coord(n,k)=0.
        enddo
      enddo
      bimp=0
      phi=0
      xproj(1)=0
      yproj(1)=0
      zproj(1)=0
      lproj(1)=koll
      lproj3(1)=0
      do k=1,koll
        bij=bkmx*sqrt(rangen())
        bk(k)=bij
        iproj(k)=1
        itarg(k)=k
        phi=2.*pi*rangen()
        xtarg(k)=bij*cos(phi)
        ytarg(k)=bij*sin(phi)
        ztarg(k)=0
        ltarg(k)=1
        kproj(1,k)=k
        ktarg(k,1)=k
        ltarg3(k)=0
        ktarg3(k,1)=0
        kproj3(k,1)=0
        if(iscreen.ne.0.and.bij.le.bkmxndif)then
          if(zbrmax.le.0..or. bij.lt.zbcutx+zbrmax*rangen())then
            lproj3(1)=lproj3(1)+1
            ltarg3(k)=1
            kproj3(1,lproj3(1))=k
            ktarg3(k,1)=k
          endif
        endif
      enddo

           elseif(maproj.eq.1.and.matarg.eq.1)then

      b1=bminim
      b2=amin1(bkmx,bmaxim)
      if(b1.gt.b2)call utstop('conaa: bmin > bmax&')
      bimp=sqrt(b1*b1+(b2*b2-b1*b1)*rangen())
      koll=1
      do n=1,4
        coord(n,1)=0.
      enddo
      bk(1)=bimp
      iproj(1)=1
      itarg(1)=1
      phi=2.*pi*rangen()
      xproj(1)=bimp*cos(phi)
      yproj(1)=bimp*sin(phi)
      zproj(1)=0
      xtarg(1)=0
      ytarg(1)=0
      ztarg(1)=0
      lproj(1)=1
      ltarg(1)=1
      lproj3(1)=1
      ltarg3(1)=1
      kproj3(1,1)=1
      ktarg3(1,1)=1
      kproj(1,1)=1
      ktarg(1,1)=1

           else

      call conxyz('p',mamx,xproj,yproj,zproj,ypjtl-yhaha)
      call conxyz('t',mamx,xtarg,ytarg,ztarg,yhaha)

      bx=0
      by=0
      if(maproj.gt.0)then
      if(bimevt.lt.0)then
        b1=bminim
        b2=amin1(rmproj+rmtarg,bmaxim)
        if(b1.gt.b2)call utstop('conaa: bmin > bmax&')
        bimp=sqrt(b1**2+(b2**2-b1**2)*rangen())
        if(nbarray.gt.0)bimp=barray(mod(nrevt,nbarray)+1)
        if(jpsi.gt.0)then
          bimp=b1+(b2-b1)*(float(mod(nrevt,12))+rangen())/12.
          bimevt=bimp
        endif
        phi=phimin+rangen()*(phimax-phimin)
      else
        phi=phievt
        bimp=bimevt
      endif
      bx=cos(phi)*bimp
      by=sin(phi)*bimp
      endif
      if(jpsi.lt.0)then !modify b
        bx=xtarg(1)
        by=ytarg(1)
      endif
      if(maproj.eq.0)goto1000
      koll=0
      do i=1,maproj
      lproj(i)=0
      lproj3(i)=0
      enddo
      do j=1,matarg
      ltarg(j)=0
      ltarg3(j)=0
      enddo
      do 12 i=1,maproj
      do 11 j=1,matarg
      if(jpsi.lt.0.and.ztarg(j).le.ztarg(1))goto11
      bij=sqrt((xproj(i)+bx-xtarg(j))**2+(yproj(i)+by-ytarg(j))**2)
      if(ish.ge.7)write(ifch,*)'i_p:',i,' i_t:',j,' b_ij:',bij
      if(bij.gt.bkmx)goto 11

      koll=koll+1
      if(koll.gt.kollmx)call utstop('conaa: kollmx too small&')
      bk(koll)=bij
      bkx(koll)=xproj(i)+bx-xtarg(j)
      bky(koll)=yproj(i)+by-ytarg(j)
      iproj(koll)=i
      itarg(koll)=j
      lproj(i)=lproj(i)+1
      ltarg(j)=ltarg(j)+1
      kproj(i,lproj(i))=koll
      ktarg(j,ltarg(j))=koll
      phik3(koll)=0.
      thetak3(koll)=0.
      if(iscreen.ne.0.and.bij.le.bkmxndif)then
        if(zbrmax.gt.0..and.bij.gt.zbcutx+zbrmax*rangen())goto 11
        lproj3(i)=lproj3(i)+1
        ltarg3(j)=ltarg3(j)+1
        kproj3(i,lproj3(i))=koll
        ktarg3(j,ltarg3(j))=koll
c define angle for anti-shadowing
        if(abs(bky(koll)).gt.1.e-6)then
          if(abs(bkx(koll)).gt.1.e-6)then
            phik3(koll)=atan(bky(koll)/bkx(koll))
          else
            phik3(koll)=sign(0.5*pi,bky(koll))
          endif
        elseif(bkx(koll).lt.0.)then
          phik3(koll)=pi
        endif
        if(bk(koll).gt.0.)then
          thetak3(koll)=atan(bglaubx/bk(koll))
        else
          thetak3(koll)=0.5*pi
        endif
      endif

11    continue
12    continue

      do k=1,koll
        do n=1,4
          coord(n,k)=0.
        enddo
      enddo

          endif

      if(ish.ge.3)write(ifch,*)'koll=',koll
      if(koll.eq.0)goto 1001

c     determine coord
c     ---------------
      do kl=1,koll
      i=iproj(kl)
      j=itarg(kl)
      dist=ztarg(j)-zproj(i)
      coord(1,kl)=(xproj(i)+xtarg(j))*0.5
      coord(2,kl)=(yproj(i)+ytarg(j))*0.5
      coord(3,kl)=(zproj(i)+ztarg(j))*0.5
      coord(4,kl)=dist/vel
      enddo

      if(iscreen.ne.0)call CalcScrPair(bimp)

    !~~~~~redefine energy in case of imposed radial flow~~~~~~~~~~~~~~~~
      yrmaxi=max(0.,yradmx+yradmi*log(1.+engy/sqrt(float(koll))))
      if(yrmaxi.gt.1e-5)then
        yyrmax=dble(yrmaxi)
        fradflii=sngl(1d0/
     &  ((sinh(yyrmax)*yyrmax-cosh(yyrmax)+1d0)/(yyrmax**2/2d0)))
      else
        fradflii=1.
      endif
      if(ish.ge.3)write(ifch,*)'yrmaxi=',yrmaxi

c     exit
c     ----
1000  continue
      if(ish.ge.5)then
      write(ifch,*)'ia,x/y/zproj:'
      do mm=1,maproj
      write(ifch,*)mm,xproj(mm),yproj(mm),zproj(mm)
      enddo
      write(ifch,*)'ia,x/y/ztarg:'
      do mm=1,matarg
      write(ifch,*)mm,xtarg(mm),ytarg(mm),ztarg(mm)
      enddo
      write(ifch,*)'iret',iret
      endif
      call utprix('conaa ',ish,ishini,4)
      return

1001  continue !iret=1 causes redo of whole collision
      iret=1
      if(ish.ge.3)then
      write(ifch,*)
      write(ifch,*)'***** subroutine conaa:'
      write(ifch,*)'***** no nucleon pair found --> no interaction'
      write(ifch,*)
      endif
      goto 1000

      end

c-----------------------------------------------------------------------
      function frdmzcut(b)
c-----------------------------------------------------------------------
c  determines interaction configuration
c-----------------------------------------------------------------------

      include 'epos.inc'
      include 'epos.incpar'

      frdmzcut=zbcutx*exp(-(b*rangen()))
c      frdmzcut=zbcutx*rangen()
c      bdum=b
      return
      end

c-----------------------------------------------------------------------
      subroutine CalcScrPair(b)
c-----------------------------------------------------------------------
c  determines interaction configuration
c-----------------------------------------------------------------------

      include 'epos.inc'
      include 'epos.incems'
      include 'epos.incpar'

      common/scrangle/ phik3(kollmx),thetak3(kollmx)
      logical lascr(kollmx),cont

c     order pairs for splitting as a function of b
c     ---------------
      do i=1,maproj
        if(lproj3(i).gt.1)then
          do l=2,lproj3(i)
            m=l
 100        continue
              kp=kproj3(i,m-1)
              kl=kproj3(i,m)
              if(bk(kl).lt.bk(kp))then
                kproj3(i,m-1)=kl
                kproj3(i,m)=kp
                m=m-1
            if(m.gt.1)goto 100
              endif
          enddo
        endif
      enddo
      do j=1,matarg
        if(ltarg3(j).gt.1)then
          do l=2,ltarg3(j)
            m=l
 200        continue
              kp=ktarg3(j,m-1)
              kl=ktarg3(j,m)
              if(bk(kl).lt.bk(kp))then
                ktarg3(j,m-1)=kl
                ktarg3(j,m)=kp
                m=m-1
            if(m.gt.1)goto 200
              endif
          enddo
        endif
      enddo

      if(koll.gt.1)then
c Define anti-shadowing as a consequence of geometrical screening of
c nucleons at large b by the one at small b

c Projectile
c Check phi not to be in the range of a nucleon with smaller b
        do i=1,maproj
          if(lproj3(i).gt.1)then
            kl=kproj3(i,1)
            lascr(kl)=.true.
            do 300 l=lproj3(i),2,-1
              kl=kproj3(i,l)
              lascr(kl)=.true.
              if(bk(kl).ge.frdmzcut(b))then
                do m=1,l-1
                  km=kproj3(i,m)
                  if(kl.ne.km.and.bk(km).ge.frdmzcut(b)
     &          .and.phik3(kl).ge.phik3(km)-thetak3(km)
     &          .and.phik3(kl).le.phik3(km)+thetak3(km))then
                    lascr(kl)=.false.
                    goto 300
                  endif
                enddo
              endif
 300        continue
          endif
        enddo
c suppress screened pair from the list
        do i=1,maproj
          if(lproj3(i).gt.1)then
            do l=1,lproj3(i)
              kl=kproj3(i,l)
            enddo
          endif
        enddo
        do i=1,maproj
          if(lproj3(i).gt.1)then
            n=2
            cont=lproj3(i).gt.1
            do while(cont)
              l=lproj3(i)
              kl=kproj3(i,l)
c suppress end of the list in order to have the last pair active
              do while(.not.lascr(kl).and.l.gt.1)
                kproj3(i,l)=0
                lproj3(i)=lproj3(i)-1
                l=lproj3(i)
                kl=kproj3(i,l)
              enddo
              cont=lproj3(i).gt.n
              if(cont)then
c compress list
                kn=kproj3(i,n)
                if(.not.lascr(kn))then
                  m=lproj3(i)   !last pair always active
                  km=kproj3(i,m)
                  kproj3(i,n)=km
                  kproj3(i,m)=0
                  lproj3(i)=lproj3(i)-1
                endif
              endif
              n=min(lproj3(i)-1,n)+1
              cont=lproj3(i).ne.n
            enddo
          endif
        enddo

c Target
c Check phi not to be in the range of a nucleon with smaller b
        do j=1,matarg
          if(ltarg3(j).gt.1)then
            kl=ktarg3(j,1)
            lascr(kl)=.true.
            do 400 l=ltarg3(j),2,-1
              kl=ktarg3(j,l)
              lascr(kl)=.true.
              if(bk(kl).ge.frdmzcut(b))then
                do m=1,l-1
                  km=ktarg3(j,m)
                  if(km.ne.kl.and.bk(km).ge.frdmzcut(b)
     &          .and.phik3(kl).ge.phik3(km)-thetak3(km)
     &          .and.phik3(kl).le.phik3(km)+thetak3(km))then
                    lascr(kl)=.false.
                    goto 400
                  endif
                enddo
              endif
 400        continue
          endif
        enddo
c suppress screened pair from the list
        do j=1,matarg
          if(ltarg3(j).gt.1)then
            n=2
            cont=ltarg3(j).gt.1
            do while(cont)
              l=ltarg3(j)
              kl=ktarg3(j,l)
c suppress end of the list in order to have the last pair active
              do while(.not.lascr(kl).and.l.gt.1)
                ktarg3(j,l)=0
                ltarg3(j)=ltarg3(j)-1
                l=ltarg3(j)
                kl=ktarg3(j,l)
              enddo
              cont=ltarg3(j).gt.n
              if(cont)then
c compress list
                kn=ktarg3(j,n)
                if(.not.lascr(kn))then
                  m=ltarg3(j)   !last pair always active
                  km=ktarg3(j,m)
                  ktarg3(j,n)=km
                  ktarg3(j,m)=0
                  ltarg3(j)=ltarg3(j)-1
                endif
              endif
              n=min(ltarg3(j)-1,n)+1
              cont=ltarg3(j).ne.n
            enddo
          endif
        enddo
      endif

      end

c-----------------------------------------------------------------------
      subroutine xGeometry(iii)
c-----------------------------------------------------------------------
      include 'epos.inc'
      include 'epos.incems'
      include 'epos.incsem'
      common/xgeom/nnn,naa(kollmx),nbb(kollmx)
      character*5 fmt1,fmt2

      if(iii.eq.0)then

      do k=1,kollmx
        naa(k)=0
      enddo
      nnn=0

      elseif(iii.eq.1)then

      ngl=0
      do k=1,koll
        r=bk(k)
        if(r.le.sqrt(sigine/10./pi))ngl=ngl+1
      enddo
      if(ngl.ne.0)then
        nnn=nnn+1
        naa(ngl)=naa(ngl)+1
      endif

      elseif(iii.eq.2)then

      if(xpar1.eq.0..and.xpar2.eq.0.)then
       print*,'---------- xGeometry -----------------------'
       return
      endif
      x1=1-0.01*xpar2
      x2=1-0.01*xpar1
      kmx=0
      nbb(1)=naa(1)
      do k=2,kollmx
        if(naa(k).ne.0.)kmx=k
        nbb(k)=nbb(k-1)+naa(k)
      enddo
      k1=0
      k2=0
      do k=1,kmx
        x=nbb(k)/float(nnn)
        if(x.lt.x1)k1=k
        if(x.lt.x2)k2=k
      enddo
      k1=k1+1
      k2=k2+1
      x=0
      av=0
      su=0
      p1=0.
      p2=0.
      do k=1,kmx
        xb=x
        x=nbb(k)/float(nnn)
        y=naa(k)/float(nnn)
        dx=x-xb
        p=0.
        if(k.eq.k1)then
          p=(x-x1)/dx
          p1=p
        elseif(k.eq.k2)then
          p=(x2-xb)/dx
          p2=p
        elseif(k.gt.k1.and.k.lt.k2)then
          p=1
        endif
      av=av+y*p*k
      su=su+y*p
      enddo
      av=av/su
      n1=nint(100*p1)
      n2=nint(100*p2)
      if(n1.eq.0)then
        k1=k1+1
        n1=100
      endif
      if(k1.le.9)fmt1='i1,4x'
      if(k1.gt.9.and.k1.le.99)fmt1='i2,3x'
      if(k1.gt.99.and.k1.le.999)fmt1='i3,2x'
      if(k1.gt.999.and.k1.le.9999)fmt1='i4,1x'
      if(k2.le.9)fmt2='i1,4x'
      if(k2.gt.9.and.k2.le.99)fmt2='i2,3x'
      if(k2.gt.99.and.k2.le.999)fmt2='i3,2x'
      if(k2.gt.999.and.k2.le.9999)fmt2='i4,1x'
      write(6,'(i4,a,i5,a,'//fmt1//',i6,a,i5,a,'//fmt2//',5x,a,f8.2)')
     &       nint(xpar2),':MIN',n1,'%',k1
     &      ,nint(xpar1),':MAX',n2,'%',k2   ,'av:',av
      endif

      end

c-----------------------------------------------------------------------
      function conbmx()
c-----------------------------------------------------------------------
      double precision om1intbc,p,eps
      include 'epos.inc'
      include 'epos.incsem'

      conbmx=0.
      b1=0.
      b2=7.
      eps=5.0d-3
      p=1.d0-dexp(-om1intbc(b2))
      if(p.gt.2.d0*eps)return

      ntry=0

10    ntry=ntry+1
      b=b1+.5*(b2-b1)

      p=(1.d0-dexp(-om1intbc(b)))

      if(p.gt.eps)then
       if(p.gt.2.d0*eps)then
        b1=b
       else
        conbmx=b
       return
       endif
      else
       if(p.lt.eps/5.d0)then
        b2=b
       else
        conbmx=b
        return
       endif
      endif

      if(ntry.le.1000)goto 10
      write(ifmt,*)'Too much try in conbmx ... bmax=',b
      conbmx=b
      return

      end

c-----------------------------------------------------------------------
      function conbmxndif()
c-----------------------------------------------------------------------
      double precision om1intbc,p,eps
      include 'epos.inc'
      include 'epos.incsem'

      iomegasave=iomega
      iomega=2
      conbmxndif=0.
      b1=0.
      b2=7.
      conbmxndif=b2
      eps=1d-10
      p=1.d0-dexp(-om1intbc(b2))
      if(p.gt.2.d0*eps)goto 100

      ntry=0

10    ntry=ntry+1
      b=b1+.5*(b2-b1)

      p=(1.d0-dexp(-om1intbc(b)))

      if(p.gt.eps)then
       if(p.gt.2.d0*eps)then
        b1=b
       else
        conbmxndif=b
       goto 100
       endif
      else
       if(p.lt.eps/5.d0)then
        b2=b
       else
        conbmxndif=b
        goto 100
       endif
      endif

      if(ntry.le.1000)goto 10
      write(ifmt,*)'Too much try in conbmxndif ... bkmxndif=',b
      conbmxndif=b
 100  iomega=iomegasave
      return

      end

c-----------------------------------------------------------------------
      function conbmxdif()
c-----------------------------------------------------------------------
c find b to have (1-exp(-om))pmax=pdiff
c-----------------------------------------------------------------------
      double precision om1intbc,pmax,drootom,pdiff
      include 'epos.inc'
      include 'epos.incsem'

      conbmxdif=0.
      b1=0.
      bmax=7.
      iomegasave=iomega
      iomega=2

      eps=1.e-5
      pmax=1.d0-dexp(-om1intbc(b1))
      pdiff=facdif
      if(pmax.gt.eps)then
        conbmxdif=drootom(pdiff,pmax,bmax,eps)
      endif
      iomega=iomegasave

      return

      end

c-----------------------------------------------------------------------
      subroutine conre
c-----------------------------------------------------------------------
c  initializes remnants
c-----------------------------------------------------------------------
      include "epos.incems"
      include 'epos.inc'

      call utpri('conre ',ish,ishini,6)

c     proj
c     ----
      la=laproj
      ma=iabs(maproj)
      las=0
      mas=0
      do l=1,ma
      if(la.lt.0)then
        if(iabs(idproj).lt.20)then
          id=idproj
        else
          ia=iabs(idproj/10)
          is=idproj/iabs(idproj)
          if(ia.ne.111.and.ia.ne.222.and.ia.ne.333)id=idproj/10*10
          if(ia.eq.111.or. ia.eq.222.or. ia.eq.333)id=idproj/10*10+is
          if(ia.eq.213)id=1230*is
        endif
      else
        id=1220
        if(rangen().le.(la-las)*1./(ma-mas))id=1120
        if(id.eq.1120)las=las+1
        mas=mas+1
      endif
      ic1=idtrai(1,id,1)
      ic2=idtrai(2,id,1)
      icproj(1,l)=ic1
      icproj(2,l)=ic2
      enddo

c     targ
c     ----
      la=latarg
      ma=iabs(matarg)
      las=0
      mas=0
      do l=1,ma
      if(la.lt.0)then
        if(iabs(idtarg).lt.20)then
          id=idtarg
        else
          ia=iabs(idtarg/10)
          is=idtarg/iabs(idtarg)
          if(ia.ne.111.and.ia.ne.222.and.ia.ne.333)id=idtarg/10*10
          if(ia.eq.111.or. ia.eq.222.or. ia.eq.333)id=idtarg/10*10+is
          if(ia.eq.213)id=1230*is
        endif
      else
        id=1220
        if(rangen().le.(la-las)*1./(ma-mas))id=1120
        if(id.eq.1120)las=las+1
        mas=mas+1
      endif
      ic1=idtrai(1,id,1)
      ic2=idtrai(2,id,1)
      ictarg(1,l)=ic1
      ictarg(2,l)=ic2
      enddo

      call utprix('conre ',ish,ishini,6)
      return
      end

c-----------------------------------------------------------------------
      subroutine conrl
c-----------------------------------------------------------------------
c  initializes target remnant in case of appl lepton
c-----------------------------------------------------------------------
      include "epos.incems"
      common/nucl1/laproj,maproj,latarg,matarg,core,fctrmx
      common/hadr2/iomodl,idproj,idtarg,wexcit

c     targ
c     ----
      la=latarg
      ma=iabs(matarg)
      las=0
      mas=0
      do l=1,ma
      if(la.lt.0)then
      id=idtarg
      else
      id=1220
      if(rangen().le.(la-las)*1./(ma-mas))id=1120
      if(id.eq.1120)las=las+1
      mas=mas+1
      endif
      ic1=idtrai(1,id,1)
      ic2=idtrai(2,id,1)
      ictarg(1,l)=ic1
      ictarg(2,l)=ic2
      enddo

      return
      end

c-----------------------------------------------------------------------
      subroutine conwr
c-----------------------------------------------------------------------
c     writes /cptl/
c-----------------------------------------------------------------------
      include "epos.inc"
      include "epos.incems"
      double precision XA(64,3),XB(64,3),BQGS,BMAXQGS,BMAXNEX,BMINNEX
      COMMON /Q_QGSNEX1/ XA,XB,BQGS,BMAXQGS,BMAXNEX,BMINNEX
      common/nucl3/phi,bimp
      common /cncl/xproj(mamx),yproj(mamx),zproj(mamx)
     *,xtarg(mamx),ytarg(mamx),ztarg(mamx)
      parameter(iapmax=208)
      double precision bqgs2,bmaxqgs2,bmaxnex2,bminnex2,xan,xbn
      common /qgsIInex1/xan(iapmax,3),xbn(iapmax,3)
     *,bqgs2,bmaxqgs2,bmaxnex2,bminnex2
      common/photrans/phoele(4),ebeam
      integer ic(2)

      call utpri('conwr ',ish,ishini,6)

      bx=cos(phi)*bimp
      by=sin(phi)*bimp

c     write /cptl/
c     ------------
      nptl=0

           if(iokoll.eq.1)then   ! precisely matarg collisions

      nptl=nptl+1
      do 3 i=1,4
3     xorptl(i,nptl)=0
      tivptl(1,nptl)=-ainfin
      tivptl(2,nptl)=0
      istptl(nptl)=1
      iorptl(nptl)=-1
      jorptl(nptl)=0
      do 1 k=1,koll
      nptl=nptl+1
      do 4 i=1,4
4     xorptl(i,nptl)=0
      tivptl(1,nptl)=-ainfin
      tivptl(2,nptl)=0
      istptl(nptl)=1
      iorptl(nptl)=-1
      jorptl(nptl)=0
1     continue

           elseif(iappl.ne.7)then

c             print *,'proj'

      do 6 i=1,maproj
      nptl=nptl+1
      istptl(nptl)=0
      iorptl(nptl)=0
      jorptl(nptl)=0
      if(model.eq.2)then       !QGSJet
      xproj(i)=XA(i,1)
      yproj(i)=XA(i,2)
      zproj(i)=XA(i,3)
      istptl(nptl)=1
      iorptl(nptl)=-1
      elseif(model.eq.7.or.model.eq.11)then       !QGSJetII
      xproj(i)=xan(i,1)
      yproj(i)=xan(i,2)
      zproj(i)=xan(i,3)
      istptl(nptl)=1
      iorptl(nptl)=-1
      elseif(model.ge.3)then       !Gheisha, ...
      istptl(nptl)=1
      iorptl(nptl)=-1
      endif
      xorptl(1,nptl)=xproj(i)+bx/2
      xorptl(2,nptl)=yproj(i)+by/2
      xorptl(3,nptl)=zproj(i)
      xorptl(4,nptl)=0
      tivptl(1,nptl)=-ainfin
c for visualisation uncomment
c-c   tivptl(1,nptl)=-100
      tivptl(2,nptl)= ainfin
c             print *,i,xorptl(1,nptl),xorptl(2,nptl),xorptl(3,nptl)
6     continue
c             print *,'targ'
      do 7 i=1,matarg
      nptl=nptl+1
      istptl(nptl)=0
      iorptl(nptl)=0
      jorptl(nptl)=0
      if(model.eq.2)then       !QGSJet
      xtarg(i)=XB(i,1)
      ytarg(i)=XB(i,2)
      ztarg(i)=XB(i,3)
      istptl(nptl)=1
      iorptl(nptl)=-1
      elseif(model.eq.7.or.model.eq.11)then       !QGSJetII
      xtarg(i)=xbn(i,1)
      ytarg(i)=xbn(i,2)
      ztarg(i)=xbn(i,3)
      istptl(nptl)=1
      iorptl(nptl)=-1
      elseif(model.ge.3)then       !Gheisha, ...
      istptl(nptl)=1
      iorptl(nptl)=-1
      endif
      xorptl(1,nptl)=xtarg(i)-bx/2
      xorptl(2,nptl)=ytarg(i)-by/2
      xorptl(3,nptl)=ztarg(i)
      xorptl(4,nptl)=0
      tivptl(1,nptl)=-ainfin
c for visualisation uncomment
c-c   tivptl(1,nptl)=-100
      tivptl(2,nptl)= ainfin
c             print *,i,xorptl(1,nptl),xorptl(2,nptl),xorptl(3,nptl)
7     continue
      if(abs(idprojin).eq.12)then   !electron for fake DIS
c electron projectile
        nptl=nptl+1
        istptl(nptl)=41
        iorptl(nptl)=-1
        jorptl(nptl)=-1
        iorptl(1)=nptl         !pi0 (porjectile) coming from lepton
        xorptl(1,nptl)=bx/2
        xorptl(2,nptl)=by/2
        xorptl(3,nptl)=0.
        xorptl(4,nptl)=0.
        tivptl(1,nptl)=-ainfin
        tivptl(2,nptl)=0.
c target nucleons (in lab frame)
        do i=1,matarg
          nptl=nptl+1
          istptl(nptl)=41
          iorptl(nptl)=-1
          jorptl(nptl)=-1
          xorptl(1,nptl)=xtarg(i)-bx/2
          xorptl(2,nptl)=ytarg(i)-by/2
          xorptl(3,nptl)=ztarg(i)
          xorptl(4,nptl)=0
          tivptl(1,nptl)=-ainfin
c         for visualisation uncomment
c         -c   tivptl(1,nptl)=-100
          tivptl(2,nptl)= ainfin
        enddo
c electron remnant
        nptl=nptl+1
        istptl(nptl)=0
        iorptl(nptl)=maproj+matarg+1
        jorptl(nptl)=-1
        xorptl(1,nptl)=bx/2
        xorptl(2,nptl)=by/2
        xorptl(3,nptl)=0.
        xorptl(4,nptl)=0.
        tivptl(1,nptl)=0.
        tivptl(2,nptl)= ainfin
      endif

          endif

      nptl=0
      if(iappl.le.2)then
      do i=1,maproj
      nptl=nptl+1
      ic(1)=icproj(1,i)
      ic(2)=icproj(2,i)
      id=idtra(ic,0,0,0)
c      id=idtra(ic,0,0,3)      !tp071107 imix=3 ??????????
      call idmass(id,ams)
      idptl(nptl)=id
      pptl(1,nptl)=0.
      pptl(2,nptl)=0.
      pptl(3,nptl)=pnullx
      pptl(4,nptl)=sqrt(pnullx**2+ams**2)
      pptl(5,nptl)=ams
      ifrptl(1,nptl)=0
      ifrptl(2,nptl)=0
      ityptl(nptl)=0
      enddo
      endif
      if(iappl.ne.7)then
      do i=1,matarg
      nptl=nptl+1
      ic(1)=ictarg(1,i)
      ic(2)=ictarg(2,i)
      id=idtra(ic,0,0,0)
c      id=idtra(ic,0,0,3)      !tp071107 imix=3 ??????????
      call idmass(id,ams)
      idptl(nptl)=id
      pptl(1,nptl)=0.
      pptl(2,nptl)=0.
      pptl(3,nptl)=-pnullx
      pptl(4,nptl)=sqrt(pnullx**2+ams**2)
      pptl(5,nptl)=ams
      ifrptl(1,nptl)=0
      ifrptl(2,nptl)=0
      ityptl(nptl)=0
      enddo
      if(abs(idprojin).eq.12)then   !electron for fake DIS
c electron projectile
        nptl=nptl+1
        id=idprojin
        call idmass(id,ams)
        idptl(nptl)=id
        pptl(1,nptl)=0.
        pptl(2,nptl)=0.
        pptl(3,nptl)=sqrt(max(0.,(elepti+ams)*(elepti-ams)))
        pptl(4,nptl)=elepti
        pptl(5,nptl)=ams
        ifrptl(1,nptl)=1
        ifrptl(2,nptl)=1
        ityptl(nptl)=40
c target nucleons (in lab frame)
        do i=1,matarg
          nptl=nptl+1
          idptl(nptl)=idptl(maproj+i)
          pptl(1,nptl)=0.
          pptl(2,nptl)=0.
          pptl(3,nptl)=-pnll
          pptl(4,nptl)=ebeam
          pptl(5,nptl)=pptl(5,maproj+i)
          ifrptl(1,nptl)=maproj+i
          ifrptl(2,nptl)=maproj+i
          ityptl(nptl)=50
        enddo
c electron remnant
        nptl=nptl+1
        idptl(nptl)=id
        pptl(1,nptl)=phoele(1)
        pptl(2,nptl)=phoele(2)
        pptl(3,nptl)=phoele(3)
        pptl(4,nptl)=phoele(4)
        pptl(5,nptl)=ams
        ifrptl(1,nptl)=0
        ifrptl(2,nptl)=0
        ityptl(nptl)=40
      endif

      else

      nptl=nptl+1
      id=idproj
      call idmass(id,ams)
      idptl(nptl)=id
      pptl(1,nptl)=0.
      pptl(2,nptl)=0.
      pptl(3,nptl)=pnullx
      pptl(4,nptl)=sqrt(pnullx**2+ams**2)
      pptl(5,nptl)=ams
      ifrptl(1,nptl)=0
      ifrptl(2,nptl)=0
      ityptl(nptl)=0
      iorptl(nptl)=-1
      jorptl(nptl)=0
      istptl(nptl)=0
      do 5 i=1,4
 5      xorptl(i,nptl)=0
      tivptl(1,nptl)=0
      tivptl(2,nptl)=0
      endif

c     exit
c     ----

      call utprix('conwr ',ish,ishini,6)
      return
      end

c------------------------------------------------------------------------
      subroutine conxyz(ch,n,x,y,z,ynuc)
c-----------------------------------------------------------------------
      include 'epos.inc'

      real x(n),y(n),z(n)
      character ch*1

      massnr=0
      iii=0
      if(ch.eq.'p')then
      massnr=maproj
      iii=1
      elseif(ch.eq.'t')then
      massnr=matarg
      iii=2
      else
      call utstop('conxyz: nucleus neither proj nor targ&')
      endif

      if(massnr.eq.0)return
      if(massnr.gt.n)call utstop('conxyz: massnr.gt.n&')
      if(massnr.eq.1)then
      x(1)=0
      y(1)=0
      z(1)=0
      return
      endif

      rad=radnuc(massnr)

      if(massnr.ge.10)then !---wood-saxon density---

        rad=rad/difnuc(massnr)
        cr1=1.+3./rad+6./rad**2+6./rad**3
        cr2=3./rad
        cr3=3./rad+6./rad**2
        do i=1,massnr
   1      zuk=rangen()*cr1-1.
          if(zuk.le.0.)then
            tt=rad*(rangen()**.3333-1.)
          elseif(zuk.le.cr2 )then
            tt=-log(rangen())
          elseif(zuk.lt.cr3 )then
            tt=-log(rangen())-log(rangen())
          else
            tt=-log(rangen())-log(rangen())-log(rangen())
          endif
          if(rangen().gt.1./(1.+exp(-abs(tt))))goto 1
          rim=tt+rad
          zz=rim*(2.*rangen()-1.)
          rim=sqrt(rim*rim-zz*zz)
          z(i)=zz*difnuc(massnr)
          call pscs(c,s)
          x(i)=rim*c*difnuc(massnr)
          y(i)=rim*s*difnuc(massnr)
        enddo

      elseif(massnr.ge.3)then  ! ---gaussian density---

        rad=rad*sqrt(2.*massnr/(massnr-1.))   !van hove simulation
        do l=1,3
          summ=0.
          do i=1,massnr-1
            j=massnr-i
            aks=rad *(rangen()+rangen()+rangen()-1.5)
            k=j+1
            if(l.eq.1)x(k)=summ-aks*sqrt(float(j)/k)
            if(l.eq.2)y(k)=summ-aks*sqrt(float(j)/k)
            if(l.eq.3)z(k)=summ-aks*sqrt(float(j)/k)
            summ=summ+aks/sqrt(float(j*k))
          enddo
          if(l.eq.1)x(1)=summ
          if(l.eq.2)y(1)=summ
          if(l.eq.3)z(1)=summ
        enddo

      elseif(massnr.eq.2)then  ! ---deuteron---

        !.........r=t*difnuc(massnr), t~exp(-2*t)*(1-exp(-a*t))
        a=radnuc(massnr)
  2     t=-0.5*alog(rangen())  !~exp(-2*t)
        if(rangen().gt.(1-exp(-a*t))**2)goto2
        r=t*difnuc(massnr)
        zz=r*(2.*rangen()-1.)
        call pscs(c,s)
        rxy=sqrt(r*r-zz*zz)
        z(1)=0.5*zz
        x(1)=0.5*rxy*c
        y(1)=0.5*rxy*s
        z(2)=-z(1)
        x(2)=-x(1)
        y(2)=-y(1)

      else

        stop'conxyz: wrong massnr.     '

      endif

c...plot preparation

      rmax=(radnuc(massnr)+3)
      drnucl(iii)=rmax/mxnucl
      nrnucl(iii)=nrnucl(iii)+1
      do i=1,massnr
        r=sqrt(x(i)**2+y(i)**2+z(i)**2)
        k=1+int(r/drnucl(iii))
        if(k.le.mxnucl)rnucl(k,iii)=rnucl(k,iii)+1
      enddo

c...lorentz trafo

      do i=1,massnr
      z(i)=z(i)/cosh(ynuc)
      enddo

      return
      end

c-----------------------------------------------------------------------
      subroutine conini
c-----------------------------------------------------------------------
      include 'epos.inc'

      imax=max(maproj,matarg)
      if(idtargin.eq.0)imax=max(imax,40)
      do massnr=1,mamxx
        dif=0.54
        rad=0.
        if(massnr.gt.imax.or.massnr.eq.1)then
          dif=0
          rad=0
        elseif(massnr.eq.197)then
          dif=0.562
          rad=6.5
        elseif(massnr.ge.10)then
          rad=1.12*massnr**0.33333-0.86*massnr**(-0.33333)
        elseif(massnr.ge.3)then
          rad=.9*float(massnr)**.3333
        elseif(massnr.eq.2)then
          dif=4.316
          rad=4.68
        endif
        difnuc(massnr)=dif
        radnuc(massnr)=rad
      enddo

      end

c-----------------------------------------------------------------------
      subroutine xConNuclDens(iii)
c-----------------------------------------------------------------------
c plots distribution of nucleons in nuclei
c  iii = 1 (proj) or 2 (targ)
c-----------------------------------------------------------------------
      include 'epos.inc'
      if(model.ne.1)return
      massnr=1
      if(iii.eq.1)then
        massnr=maproj
      elseif(iii.eq.2)then
        massnr=matarg
      endif
      if(massnr.eq.1)return
      a=1./4.316
      b=4.68
      write(ifhi,'(a)') '!-----------------------------------------'
      write(ifhi,'(a)') '!          nuclear density          '
      write(ifhi,'(a)') '!-----------------------------------------'
      write(ifhi,'(a)')       'openhisto'
      if(massnr.ge.10)write(ifhi,'(a)')'htyp lin xmod lin ymod lin'
      if(massnr.lt.10)write(ifhi,'(a)')'htyp lin xmod lin ymod log'
      write(ifhi,'(a)')       'text 0 0 "title nuclear density"'
      write(ifhi,'(a)')       'text 0.99 -0.15 " r (fm) "'
      write(ifhi,'(a)')       'text 0 0 "yaxis rho(r)"'
      write(ifhi,'(a,2e11.3)')'xrange',0.,mxnucl*drnucl(iii)
      write(ifhi,'(3a)')'yrange',' 0 ',' auto'
      write(ifhi,'(a)')       'array 2'
      do j=1,mxnucl
      r=(j-0.5)*drnucl(iii)
      d=0.5*drnucl(iii)
      write(ifhi,'(2e12.4)')  r,rnucl(j,iii)/nrnucl(iii)/
     *                     (4./3.*pi*((r+d)**3-(r-d)**3))
      enddo
      write(ifhi,'(a)')       '  endarray'
      write(ifhi,'(a)')       'closehisto plot 0-'
      write(ifhi,'(a)')       'openhisto'
      write(ifhi,'(a)')       'htyp lbo '
      write(ifhi,'(a)')       'array 2'
      do j=1,mxnucl
      r=(j-0.5)*drnucl(iii)
      rr=2*r
      rho=1
      if(massnr.eq.2)then
        rho=1.00*((1-exp(-b*a*rr))*exp(-a*rr)/rr)**2
      elseif(massnr.eq.197)then
        rho=0.16/(1+exp((r-6.5)/0.562))
      elseif(massnr.ge.10)then
        rho=0.16/(1+exp((r-radnuc(massnr))/difnuc(massnr)))
      endif
      write(ifhi,'(2e12.4)')  r,rho
      enddo
      write(ifhi,'(a)')       '  endarray'
      write(ifhi,'(a)')       'closehisto plot 0'
      end

c-----------------------------------------------------------------------
      subroutine xConThick(iii)
c-----------------------------------------------------------------------
      ! plots sigma_pp *T_A (b)  (=average number of collisions)
      ! T_A = thickness function
      !  iii = 1 (proj) or 2 (targ)
      !----------------------------------------------------------------
      include 'epos.inc'
      parameter(iconimax=20,iconkmax=100)
      real thick(2,0:iconimax)
      imax=iconimax
      kmax=iconkmax
      if(model.ne.1)return
      ramx=mxnucl*drnucl(iii)
      do i=0,imax
        x=i/float(imax)*ramx
        sum=0
        rho0=conrho(iii,0.)
        h=ramx/kmax
        do k=1,kmax
          z=k*h
          r=sqrt(x*x+z*z)
          rho2=conrho(iii,r)
          z=(k-0.5)*h
          r=sqrt(x*x+z*z)
          rho1=conrho(iii,r)
          sum=sum+h/6.*(rho0+4*rho1+rho2)
          rho0=rho2
        enddo
        sum=sum*2 ! integral fro -infty to +infty
        thick(1,i)=x
        thick(2,i)=sum
      enddo
      write(ifhi,'(a)')       'openhisto'
      write(ifhi,'(a)')       'htyp lbo '
      write(ifhi,'(a)')       'txt "xaxis b (fm)" '
      write(ifhi,'(a)')       'txt "yaxis [s]?pp! T?A! (b) " '
      write(ifhi,'(a)')       'array 2'
      do i=0,imax
        write(ifhi,'(2e12.4)') thick(1,i),sigine/10.*thick(2,i)
      enddo
      write(ifhi,'(a)')       '  endarray'
      write(ifhi,'(a)')       'closehisto plot 0'

      end

c-----------------------------------------------------------------------
      function conrho(iii,r)
c-----------------------------------------------------------------------
      ! nuclear density
      !  iii = 1 (proj) or 2 (targ)
      !----------------------------------------------------------------
      include 'epos.inc'
      conrho=1.
      if(model.ne.1)return
      massnr=1
      if(iii.eq.1)then
        massnr=maproj
      elseif(iii.eq.2)then
        massnr=matarg
      endif
      if(massnr.eq.1)return
      a=1./4.316
      b=4.68
      rr=2*r
      rho=1
      if(massnr.eq.2.and.rr.gt.0.)then
        rho=1.00*((1-exp(-b*a*rr))*exp(-a*rr)/rr)**2
      elseif(massnr.eq.197)then
        rho=0.16/(1+exp((r-6.5)/0.562))
      elseif(massnr.ge.10)then
        rho=0.16/(1+exp((r-radnuc(massnr))/difnuc(massnr)))
      endif
      conrho=rho
      end

ext: f

max_stars_repo_name: kotania/impy

max_issues_repo_name: kotania/impy

max_forks_repo_name: kotania/impy


text:
+incdir+..
../riscv_vip_pkg.sv
../riscv_vip_class_pkg.sv
riscv_vip_uvc_pkg.sv

max_stars_repo_name: kencycheng/riscv-vip

max_issues_repo_name: kencycheng/riscv-vip

max_forks_repo_name: kencycheng/riscv-vip


text:
! (C) Copyright 2009-2016 ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation nor
! does it submit to any jurisdiction.

!> Linked list implementation

!> Linked list subroutines

!> Initialize the linked list
subroutine init_(self)
 class(registry_t) :: self

 !set count to zero and allocate the head of the list
 if(.not.self%l_init.or..not.associated(self%head)) then
  self%count = 0
  allocate(self%head)
  nullify(self%head%next)
  self%l_init=.true.
 endif
end subroutine

!> Add element to the linked list
subroutine add_(self,key)
 class(registry_t) :: self
 integer           :: key

 type(node_t), pointer :: next
 type(node_t), pointer :: current

 !increase global counter and assign key
 self%count = self%count+1
 key = self%count

 !allocate next element and assign key
 allocate(next)
 next%key = key

 !move the head to the front of the list
 next%next => self%head%next
 self%head%next => next
end subroutine

!> Fetch element of the linked list by key
subroutine get_(self,key,ptr)
 class(registry_t)           :: self
 integer                     :: key
 type (LISTED_TYPE), pointer :: ptr

 type(node_t), pointer :: next
 integer               :: lkey

 !note that the list starts from self%head%next
 next => self%head
 ptr => NULL()

 !sweep the linked list to find matching key
 do while(associated(next))
  next=>next%next
  if(key.eq.next%key) then
   ptr => next%element
   exit
  endif
 enddo
end subroutine

!> Remove element of the linked list
subroutine remove_(self,key)
 class(registry_t) :: self
 integer           :: key

 type(node_t), pointer :: prev
 type(node_t), pointer :: next
 integer               :: lkey

 next => self%head%next
 prev => NULL()

 !sweep the linked list to find matching key,
 do while(associated(next))
  if(key.eq.next%key) exit
  prev => next
  next => next%next
 enddo
 !reconnect the list
 if(associated(next%next)) then
  if(associated(prev)) then
   prev%next => next%next
  else
   self%head%next=>next%next
  endif
 endif
 !remove the node and set key to 0
 if(associated(next)) deallocate(next)
 key=0
 return
end subroutine

!> Finalize the linked list, deallocate all nodes
subroutine finalize_(self)
 class(registry_t) :: self
 type(node_t), pointer :: current
 type(node_t), pointer :: next

 !sweep the linked list and deallocate all nodes
 next => self%head
 do while(associated(next))
  current => next
  next => next%next
  deallocate(current)
 enddo
end subroutine


text:
! Bug report from Rice: 15-generic-typebound-procedure.f90
! A generic type-bound procedure name causes an assertion failure.

module m101
  type  :: t
  contains
    generic :: p => q  ! assertion failure: name stack not empty
  end type
end module

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
!***********************************************************************
!                                                                      *
      SUBROUTINE BESSEL(IA, IB, IK, IW, K)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
!                                                                      *
!   This routine evaluates the functions                               *
!                                                                      *
!                  (2K+1)!!                                            *
!      BESSJ  =    --------  J  (w   r) - 1  =  PHI  (w   r) - 1       *
!                         K   K   ab               K   ab              *
!                  (w   r)                                             *
!                    ab                                                *
!   and                                                                *
!                         K+1                                          *
!                  (w   r)                                             *
!                    ab                                                *
!      BESSN  =  - --------  N  (w   r) - 1  =  PSI  (w   r) - 1       *
!                  (2K-1)!!   K   ab               K   ab              *
!                                                                      *
!                                                                      *
!   where J and N ARE spherical Bessel functions, and                  *
!                                                                      *
!                w   = ABS ( E(IA) - E(IB) )/c                         *
!                 ab                                                   *
!                                                                      *
!   where  E(I)  is the eigenvalue for orbital  I . The writeup (B J   *
!   McKenzie, I P Grant, and P H Norrington, Computer Phys Commun 21   *
!   (1980) 233-246) is incorrect in its description of the output of   *
!   this routine.                                                      *
!                                                                      *
!   The  routine uses equations given in M Abramowitz and I A STegun   *
!   to evaluate the functions. Devices are used to reduce the number   *
!   of actual evaluations of these functions.                          *
!                                                                      *
!                                           Last update: 09 Oct 1992   *
!                                                                      *
!***********************************************************************
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:04:58   1/ 3/07
!...Modified by Charlotte Froese Fischer
!                     Gediminas Gaigalas  10/05/17
!-----------------------------------------------
!   M O D U L E S
!-----------------------------------------------
      USE vast_kind_param, ONLY: DOUBLE
      USE parameter_def,   ONLY: KEYORB
      USE bess_C
      USE debug_C
      USE def_C
      USE grid_C
      USE orb_C
      USE stor_C
      USE wfaC_c
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER, INTENT(IN) :: IA
      INTEGER, INTENT(IN) :: IB
      INTEGER, INTENT(IN) :: IK
      INTEGER, INTENT(IN) :: IW
      INTEGER, INTENT(IN) :: K
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
      INTEGER, PARAMETER :: KEY = KEYORB
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: ICODE, IWKP, IKKP, L, I, NN, J, JCHAN, IREM, ISWAP
      REAL(DOUBLE) :: EPSI, W, WA, XBESS1, XBESS2, S1, S2, DFNM, DFN, SSN, SCN&
         , SN, CN, OBWA, B, SKEEP
!-----------------------------------------------
!
      EPSI = SQRT(0.1D00*ACCY)
!
!   Form unique label symmetric in IA, IB
!
      ICODE = MAX(IA,IB) + KEY*(MIN(IA,IB) + KEY*K)
!
!   Function in position; return
!
      IF (ICODE == KEEP(IK,IW)) RETURN
!
!   Function not in position; is it available in BESS arrays?
!
      W = WFACT*ABS(E(IA)-E(IB))/C
      WIJ(IW) = W
!
      DO IWKP = 1, 2
         DO IKKP = 1, 2
            IF (KEEP(IKKP,IWKP) /= ICODE) CYCLE
!
!   Function found move into position
!
            KEEP(IK,IW) = ICODE
            IF (LDBPR(7)) WRITE (99, 302) NP(IA), NH(IA), NP(IB), NH(IB), K, &
               IKKP, IWKP, IK, IW
            BESSJ(IK,IW,:N) = BESSJ(IKKP,IWKP,:N)
            BESSN(IK,IW,:N) = BESSN(IKKP,IWKP,:N)
            RETURN
         END DO
      END DO
!
!   Function not found; evaluate it
!
      IF (LDBPR(7)) WRITE (99, 303) NP(IA), NH(IA), NP(IB), NH(IB), K, IK, IW
!
      KEEP(IK,IW) = ICODE
!
      IF (W < EPSI**2) THEN
!
!   Negligible w
!
         BESSJ(IK,IW,:N) = 0.0D00
         BESSN(IK,IW,:N) = 0.0D00
         RETURN
!
      ENDIF
!
      NN = K
!
      BESSJ(IK,IW,1) = 0.0D00
      BESSN(IK,IW,1) = 0.0D00
!
!   Use a four-term power series for low w*r
!
      L5: DO J = 2, N
         WA = -0.5D00*(R(J)*W)**2
         XBESS1 = 1.0D00
         XBESS2 = 1.0D00
         S1 = 0.0D00
         S2 = 0.0D00
         DO I = 1, 4
            XBESS1 = XBESS1*WA/DBLE(I*(2*(NN + I) + 1))
            XBESS2 = XBESS2*WA/DBLE(I*(2*(I - NN) - 1))
            S1 = S1 + XBESS1
            S2 = S2 + XBESS2
            IF (ABS(XBESS1)>=ABS(S1)*EPSI .OR. ABS(XBESS2)>=ABS(S2)*EPSI) &
               CYCLE
            BESSJ(IK,IW,J) = S1
            BESSN(IK,IW,J) = S2
            CYCLE  L5
         END DO
         JCHAN = J
         GO TO 6
      END DO L5
!
!   If here then calculated whole array using four-term power
!   series.  Hence return
!
      RETURN
!
!   Use sin/cos expansion when power series requires more than
!   four terms terms to converge
!
    6 CONTINUE
      IF (NN == 0) THEN
         DFNM = 1.0D00
         DFN = 1.0D00
      ELSE
         DFNM = 1.0D00
         DO I = 3, 2*NN - 1, 2
            DFNM = DFNM*DBLE(I)
         END DO
         DFN = DFNM*DBLE(2*NN + 1)
      ENDIF
      DFNM = 1.0D00/DFNM
!
      IREM = MOD(NN,4)
!
      SELECT CASE (IREM)
      CASE (1)
!
!   NN = 1, 5, 9, ...
!
         SSN = -1.0D00
         SCN = 1.0D00
         ISWAP = 1
!
      CASE (2)
!
!   N = 2, 6, 10, ....
!
         SSN = -1.0D00
         SCN = -1.0D00
         ISWAP = 0
!
      CASE (3)
!
!   NN = 3, 7, 11,...
!
         SSN = 1.0D00
         SCN = -1.0D00
         ISWAP = 1
!
      CASE DEFAULT
!
!   NN = 0, 4, 8,...
!
         SSN = 1.0D00
         SCN = 1.0D00
         ISWAP = 0
!
      END SELECT
!
      DO J = JCHAN, N
         WA = W*R(J)
         IF (ISWAP == 0) THEN
            SN = SSN*SIN(WA)
            CN = SCN*COS(WA)
         ELSE
            SN = SSN*COS(WA)
            CN = SCN*SIN(WA)
         ENDIF
         OBWA = 1.0D00/WA
         B = OBWA
         S1 = B*SN
         S2 = B*CN
         DO I = 1, NN
            SKEEP = SN
            SN = CN
            CN = -SKEEP
            B = B*OBWA*DBLE((NN + I)*(NN - I + 1))/DBLE(2*I)
            S1 = S1 + B*SN
            S2 = S2 + B*CN
         END DO
         S1 = S1*DFN/WA**NN - 1.0D00
         S2 = S2*WA**(NN + 1)*DFNM - 1.0D00
         BESSJ(IK,IW,J) = S1
         BESSN(IK,IW,J) = S2
      END DO
      RETURN
!
  303 FORMAT(93X,I2,A2,2X,I2,A2,2X,I2,2X,'New',6X,'(',I2,',',I2,')')
  302 FORMAT(93X,I2,A2,2X,I2,A2,2X,I2,2X,'(',I2,',',I2,')',2X,'(',I2,',',I2,')'&
         )
      RETURN
!
      END SUBROUTINE BESSEL

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: sylas/grasp-continuum

max_forks_repo_name: sylas/grasp-continuum


text:
subroutine lmn_bouwal(a,e,gpvel,visac,denac,delta,&
                      wmatr,tract)

!-----------------------------------------------------------------------
!****f* Nstinc/lmn_bouwal
! NAME
!    lmn_bouwal
! DESCRIPTION
!    This routine computes the surface traction for the NS equations at
!    a given integration point of a boundary IBOUN received by argument
!    due to the use of a turbulent wall law. The algorithm is:
!    - Compute the tangential velocity u at the integration point x.
!      In fact, u is not necessarily tangential at x. For example:
!      -> u    -> u
!      o---x---o---x---o
!                      |\ u
!                      |
!    - Given the distance to the wall y, compute (U*) (friction velocity)
!    - Compute y+=(y)(U*)/(nu)
!      if(y+>5) then
!        t=-(rho)*(U*)^2*(u/|u|)
!      else
!        u+=y+ => U*^2=u*nu/y so that
!        t=-mu*u/y
!      end if
! USES
!    vecnor
!    frivel
! USED BY
!    lmn_bouope
!***
!-----------------------------------------------------------------------
   use typre
   use Mod_Element
   use Mod_LowMach
   implicit none
   type(LowMachProblem) :: a
   class(FiniteElement)  :: e
   real(rp) :: gpvel(e%ndime),visac,denac,delta
   real(rp) :: wmatr(a%ndofn,e%mnode,a%ndofn,e%mnode),tract(e%ndime)

   integer(ip)                :: jlins
   integer(ip)                :: ldime,idime,ldimb

   real(rp)                   :: vikin,velfr,yplus    ! nu, U*, y+, y
   real(rp)                   :: tveno,tvelo(e%ndime),ovelo(e%ndime)             ! |u|, u
   real(rp)                   :: trano,xmuit,tcoef,dtrac
   integer(ip)                :: inodb,jnodb,jdime,kdime,inode,jnode
   real(rp)                   :: relax

   jlins = 1
   if(delta<zelmn) return

   !normal velocity, projection of a over b: a_b = (a . b) b
   ovelo = dot_product(e%baloc(:,e%ndime),gpvel)*e%baloc(:,e%ndime)

   tvelo = gpvel - ovelo
   call vecnor(tvelo,e%ndime,tveno,2)                       ! |u|
   vikin=visac/denac                                      ! nu
   if(tveno<=zelmn) then
      velfr=0.0_rp
      tract=0.0_rp
      return
   else
      call frivel(delta,tveno,vikin,velfr)                  ! U*
   end if

   ! Compute prescribed traction
   if(jlins==0) then                                       ! RHS

       yplus=delta*velfr/vikin
       if(yplus<5.0_rp) then
          tract = -denac*vikin*tvelo/delta                  ! t=-mu*u/y
       else
          tract = -denac*velfr*velfr*tvelo/tveno            ! t=-rho*U*^2*(u/|u|)
       end if

   else if(jlins==1) then                                  ! Assembly (the sign changes on the left)

      yplus=delta*velfr/vikin
      if(yplus<5.0_rp) then
         trano=denac*vikin/delta                              ! t=-mu*u/y
      else
         trano=denac*velfr*velfr/tveno                        ! t=-rho*U*^2*(u/|u|)
      end if
      do inodb=1,e%pnodb
         inode = e%lboel(inodb)
         do idime=1,e%ndime
            do jnodb=1,e%pnodb
               jnode = e%lboel(jnodb)
               xmuit=e%shapb(inodb,e%igaub)*e%shapb(jnodb,e%igaub)
               do jdime=1,e%ndime
                  do ldimb=1,e%ndime-1
                     wmatr(idime,inode,jdime,jnode) = &
                        wmatr(idime,inode,jdime,jnode) + &
                         xmuit*e%baloc(idime,ldimb)*e%baloc(jdime,ldimb)*trano
                  end do
               end do
            end do
         end do
      end do

   elseif(jlins==2) then                                   ! Newton on the traction

   else
      call runend('lmn_bouwal: jlins not defined')
   end if

end subroutine lmn_bouwal

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
!> Author:   Jabir Ali Ouassou
!> Category: System
!>
!> This file renames the ISO input/output units to the standard UNIX names,
!> defines the ANSI escape codes for colored output, and defines a number of
!> auxiliary subroutines for e.g. writing out error and warning messages.

module stdio_m
  use, intrinsic :: iso_fortran_env
  public

  ! Declare standard input/output units
  integer :: stdin  = input_unit
  integer :: stdout = output_unit
  integer :: stderr = error_unit

  ! Define escape codes for terminal colors
  character(*), parameter :: color_none   = ''
  character(*), parameter :: color_red    = ''
  character(*), parameter :: color_green  = ''
  character(*), parameter :: color_yellow = ''
  character(*), parameter :: color_blue   = ''
  character(*), parameter :: color_purple = ''
  character(*), parameter :: color_cyan   = ''
  character(*), parameter :: color_white  = ''

  ! Declare public interfaces
  interface dump
    !! Public interface for functions that dump results to files
    module procedure dump_arrays, dump_scalar
  end interface
contains
  impure subroutine message(msg)
    !! This subroutine provides a way to report a status message.
    character(*), intent(in) :: msg

    write(stderr,'(a)') color_green  // ' >> INFO: '    // color_none // msg
    flush(stderr)
  end subroutine

  impure subroutine warning(msg)
    !! This subroutine provides a way to report a warning message.
    character(*), intent(in) :: msg

    write(stderr,'(a)') color_yellow // ' >> WARNING: ' // color_none // msg
    flush(stderr)
  end subroutine

  impure subroutine error(msg)
    !! This subroutine provides a way to report an error message and halt the program.
    character(*), intent(in) :: msg

    write(stderr,'(a)') color_red    // ' >> ERROR: '   // color_none // msg
    flush(stderr)
    stop
  end subroutine

  impure subroutine status_head(title)
    !! This subroutine is used to write boxed status messages to standard out;
    !! in particular, this routine writes out a boxed title with a timestamp.
    character(len=*), intent(in) :: title
    character(len=33)            :: title_
    real                         :: time
    integer                      :: hh, mm, ss

    ! Calculate the current time
    call cpu_time(time)
    hh = int(time/3600.0)
    mm = int(mod(time,3600.0)/60.0)
    ss = int(mod(time,60.0))

    ! adjust the provided title
    title_ = ''
    title_((len(title_)-len(title)+1)/2:) = title

    ! Write out the boxed header
    write(stdout,*)
    write(stdout,'(a)') &
      '╒═══════════════════════════════════╕'
    write(stdout,'(a)') &
      '│ '         // title_ //          ' │'
    write(stdout,'(a)') &
      '├───────────────────────────────────┤'
    write(stdout,'(a,3x,a,7x,i3.2,a,i2.2,a,i2.2,3x,a)') &
      '│', 'Elapsed time:', hh, ':', mm, ':', ss, '│'
  end subroutine

  impure subroutine status_body(title, value)
    !! This subroutine is used to write boxed status messages to standard out;
    !! in particular, this routine writes out the name and value of a variable.
    character(len=*), intent(in) :: title
    character(len=20)            :: title_
    class(*),         intent(in) :: value

    ! Adjust the provided title
    title_ = trim(title) // ':'

    ! Print out the title and value
    select type(value)
      type is (integer)
        write(stdout,'(a,3x,a,i10  ,2x,a)') '│', title_, value, '│'
      type is (real)
        write(stdout,'(a,3x,a,f10.8,2x,a)') '│', title_, value, '│'
      type is (double precision)
        write(stdout,'(a,3x,a,f10.8,2x,a)') '│', title_, value, '│'
    end select
  end subroutine

  impure subroutine status_foot()
    !! This subroutine is used to write boxed status messages to standard out;
    !! in particular, this routine writes out the bottom edge of such a box.

    ! Write out the boxed footer
    write(stdout,'(a)') &
      '╘═══════════════════════════════════╛'

    ! Flush the information to standard out
    flush(unit=stdout)
  end subroutine

  impure subroutine status_box(title)
    !! This subroutine is used to write boxed status messages to standard out.
    character(len=*), intent(in) :: title
    character(len=33)            :: title_

    ! adjust the provided title
    title_ = ''
    title_((len(title_)-len(title)+1)/2:) = title

    ! Write out the boxed message
    write(stdout,*)
    write(stdout,'(a)') &
      '╒═══════════════════════════════════╕'
    write(stdout,'(a)') &
      '│ '         // title_ //          ' │'
    write(stdout,'(a)') &
      '╘═══════════════════════════════════╛'
  end subroutine

  impure function input(file) result(unit)
    !! This function is used to open an input file for reading.
    character(len=*), intent(in) :: file
    integer                      :: iostat
    integer                      :: unit

    ! Open the output file
    open(newunit = unit, file = file, iostat = iostat, action = 'read', status = 'old')
    if (iostat /= 0) then
      call error('Failed to open input file "' // file // '"!')
    end if
  end function

  impure function output(file) result(unit)
    !! This function is used to open an output file for writing.
    character(len=*), intent(in) :: file
    integer                      :: iostat
    integer                      :: unit

    ! Open the output file
    open(newunit = unit, file = file, iostat = iostat, action = 'write', status = 'replace')
    if (iostat /= 0) then
      call error('Failed to open output file "' // file // '"!')
    end if
  end function

  impure subroutine dump_arrays(filename, arrays, header)
    !! This subroutine is used to dump numerical arrays to an output file.
    use :: iso_fortran_env

    character(len=*),               intent(in) :: filename
    real(real64),     dimension(:), intent(in) :: arrays
    character(len=*), dimension(:), intent(in) :: header

    real(real64), dimension(size(arrays)/size(header),size(header)) :: matrix
    integer :: unit
    integer :: n

    ! Reshape the data
    matrix = reshape(arrays, shape(matrix))

    ! Open the output file
    unit = output(filename)

    ! Write the header line
    write(unit,'(*(a20,:,"	"))') '# ' // header(1), header(2:)

    ! Loop over the matrix rows
    do n=1,size(matrix,1)
      ! Write the matrix column to file
      write(unit,'(*(es20.12e3,:,"	"))') matrix(n,:)
    end do

    ! Close the output file
    close(unit = unit)
  end subroutine

  impure subroutine dump_scalar(filename, scalar)
    !! This subroutine is used to dump a numerical result to an output file.
    use :: iso_fortran_env

    character(len=*), intent(in) :: filename
    real(real64),     intent(in) :: scalar

    character(len=2048) :: str

    integer :: unit
    integer :: n

    ! Open the output file
    unit = output(filename)

    ! Write out command-line arguments
    do n=1,command_argument_count()
      call get_command_argument(n, str)
      write(unit,'(a,"	")',advance='no') trim(str)
    end do

    ! Write the scalar value to file
    write(unit,'(*(es20.12e3,:,"	"))') scalar

    ! Close the output file
    close(unit = unit)
  end subroutine
end module

max_stars_repo_name: jabirali/NEUS

max_issues_repo_name: jabirali/NEUS

max_forks_repo_name: jabirali/NEUS


text:
! This program provides an elementary check of some of the parts of the
!   Fortran 90 interface to netCDF 3.5. It is a Fortran 90 implementation
!   of the nctst.cpp program provided with the C++ interface to netcdf
!   (in the src/cxx directory of the netcdf distribution).
!
program netcdfTest
  use typeSizes
  use netcdf
  implicit none

  ! netcdf related variables
  integer :: ncFileID,                                   &
             latDimID, lonDimID, frTimeDimID, timeDimID, &
             pressVarID, latVarID, lonVarID, frTimeVarID, refTimeVarID, scalarVarID

  ! Local variables
  integer, parameter :: numLats = 4, numLons = 3, &
                        numFrTimes = 2, timeStringLen = 20
  character (len = *), parameter :: fileName = "example.nc"
  integer :: counter
  real, dimension(numLons, numLats, numFrTimes) :: pressure
  integer (kind = FourByteInt), dimension(numFrTimes) :: frTimeVals
  real (kind = FourByteReal) fillVal;
  real (kind = FourByteReal), dimension(2) :: validRange;

  ! --------------------
  ! Code begins
  ! --------------------
  if(.not. byteSizesOK()) then
    print *, "Compiler does not appear to support required kinds of variables."
    stop
  end if

  ! Create the file
  call check(nf90_create(path = trim(fileName), cmode = nf90_clobber, ncid = ncFileID))

  ! Define the dimensions
  call check(nf90_def_dim(ncid = ncFileID, name = "lat",     len = numLats,        dimid = latDimID))
  call check(nf90_def_dim(ncid = ncFileID, name = "lon",     len = numLons,        dimid = lonDimID))
  call check(nf90_def_dim(ncid = ncFileID, name = "frtime",  len = nf90_unlimited, dimid = frTimeDimID))
  call check(nf90_def_dim(ncid = ncFileID, name = "timelen", len = timeStringLen,  dimid = timeDimID))

  ! Create variables and attributes
  call check(nf90_def_var(ncid = ncFileID, name = "P", xtype = nf90_float,     &
                     dimids = (/ lonDimID, latDimID, frTimeDimID /), varID = pressVarID) )
  call check(nf90_put_att(ncFileID, pressVarID, "long_name",   "pressure at maximum wind"))
  call check(nf90_put_att(ncFileID, pressVarID, "units",       "hectopascals") )
  ! Use 4-byte reals explicitly, to match 4-byte attribute type in test file
  validRange(1) = 0.
  validRange(2) = 1500
  call check(nf90_put_att(ncFileID, pressVarID, "valid_range", validRange))
  ! Use a 4-byte float constant, to match variable type
  fillVal = -9999.0
  call check(nf90_put_att(ncFileID, pressVarID,  "_FillValue", fillVal ) )

  call check(nf90_def_var(ncFileID, "lat", nf90_float, dimids = latDimID, varID = latVarID) )
  call check(nf90_put_att(ncFileID, latVarID, "long_name", "latitude"))
  call check(nf90_put_att(ncFileID, latVarID, "units", "degrees_north"))

  call check(nf90_def_var(ncFileID, "lon", nf90_float, lonDimID, lonVarID) )
  call check(nf90_put_att(ncFileID, lonVarID, "long_name", "longitude"))
  call check(nf90_put_att(ncFileID, lonVarID, "units",     "degrees_east"))

  call check(nf90_def_var(ncFileID, "frtime", nf90_int, frTimeDimID, frTimeVarID) )
  call check(nf90_put_att(ncFileID, frTimeVarID, "long_name", "forecast time"))
  call check(nf90_put_att(ncFileID, frTimeVarID, "units",     "hours"))

  call check(nf90_def_var(ncFileID, "reftime", nf90_char, timeDimID, refTimeVarID) )
  call check(nf90_put_att(ncFileID, refTimeVarID, "long_name", "reference time"))
  call check(nf90_put_att(ncFileID, refTimeVarID, "units",     "text_time"))

  ! In the C++ interface the define a scalar variable - do we know how to do this?
  call check(nf90_def_var(ncFileID, "ScalarVariable", nf90_real, scalarVarID))

  ! Global attributes
  call check(nf90_put_att(ncFileID, nf90_global, "history", &
                     "created by Unidata LDM from NPS broadcast"))
  call check(nf90_put_att(ncFileID, nf90_global, "title", &
                     "NMC Global Product Set: Pressure at Maximum Wind"))

  ! Leave define mode
  call check(nf90_enddef(ncfileID))

  ! Write the dimension variables
  call check(nf90_put_var(ncFileID, latVarId,     (/ -90., -87.5, -85., -82.5 /)) )
  call check(nf90_put_var(ncFileID, lonVarId,     (/ -180, -175, -170 /)      ) )
  ! Don't use anonymous array here, in case platform has 8-byte integers
  frTimeVals(1) = 12
  frTimeVals(2) = 18
  call check(nf90_put_var(ncFileID, frTimeVarId,  frTimeVals                  ) )
  call check(nf90_put_var(ncFileID, reftimeVarID, "1992-3-21 12:00"           ) )

  ! Write the pressure variable. Write a slab at a time to check incrementing.
  pressure = 949. + real(reshape( (/ (counter, counter = 1, numLats * numLons * numFrTimes) /),  &
                                    (/ numLons, numLats, numFrTimes /) ) )
  call check(nf90_put_var(ncFileID, pressVarID, pressure(:, :, 1:1)) )
  call check(nf90_put_var(ncFileID, pressVarID, pressure(:, :, 2:2), start = (/ 1, 1, 2 /)) )

  call check(nf90_put_var(ncFileID, scalarVarID, 10))
  call check(nf90_close(ncFileID))

contains
  ! Internal subroutine - checks error status after each netcdf, prints out text message each time
  !   an error code is returned.
  subroutine check(status)
    integer, intent ( in) :: status

    if(status /= nf90_noerr) then
      print *, trim(nf90_strerror(status))
    end if
  end subroutine check
end program netcdfTest

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
      subroutine chksrr(delxi,dlxmin,drer0,dxval0,eps100,iodb,jreac,
     $ nodbmx,noutpt,nord,nrct,nrctmx,nrd1mx,nttyo,rrelr0,rrelrp,
     $ tolsrr,ureac,xi0,xi1,xval0)
c
c     This subroutine checks the signs of the relative rates. It finds
c     the point of reaction progress at which any relative rate of an
c     irreversible reaction becomes zero. The relative rates are
c     tracked using finite differences.
c
c     This subroutine is called by:
c
c       EQ6/path.f
c
c-----------------------------------------------------------------------
c
c     Principal input:
c
c
c     Principal output:
c
c
c-----------------------------------------------------------------------
c
      implicit none
c
c-----------------------------------------------------------------------
c
c     Calling sequence variable declarations.
c
      integer nodbmx,nrctmx,nrd1mx
c
      integer noutpt,nttyo
c
      integer iodb(nodbmx),jreac(nrctmx)
c
      integer nord,nrct
c
      character*24 ureac(nrctmx)
c
      real*8 drer0(nrd1mx,nrctmx),dxval0(nrd1mx),rrelr0(nrctmx),
     $ rrelrp(nrctmx)
c
      real*8 delxi,dlxmin,eps100,tolsrr,xi0,xi1,xval0
c
c-----------------------------------------------------------------------
c
c     Local variable declarations.
c
      integer icount,ier,ilsign,j2,krzero,n,nrc,nrzero
c
      integer ilnobl
c
      character*48 usearch
      character*24 unam24
c
      real*8 arrxp,atgsrr,dxp,rrx,rrx0,rrxp,tolsx,xtargv,xval
c
      real*8 fctrl
c
c-----------------------------------------------------------------------
c
      data usearch /'a relative rate changes sign                    '/
c
c-----------------------------------------------------------------------
c
      if (nord .le. 0) go to 999
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     The search target is not where a relative rate is zero. The
c     magnitude of the search target is atgsrr = 0.5*tolsrr. The search
c     tolerance has the same value. Thus, in the case of a relative rate
c     going from positive to negative values, the search target is
c     -atgsrr and the interval for convergence is (-tolsrr,0.). In the
c     case of a relative rate going from negative to positive values,
c     the search target is atgsrr, and the interval for convergence is
c     (0.,tolsrr).
c
      atgsrr = 0.5*tolsrr
      tolsx = atgsrr
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
      icount = -1
  100 xi1 = xi0 + delxi
      icount = icount + 1
      if (icount .gt. 50) then
        if (iodb(1).gt.0 .or. iodb(5).gt.0) then
          write (noutpt,1000)
 1000     format(/3x,'A scan for reactant relative rates changing sign',
     $    ' has failed.',/3x,'Dropping to the minimum step size.')
        endif
        delxi = dlxmin
        go to 999
      endif
c
c     Estimate the relative rates from Taylor's series expansions.
c
      do nrc = 1,nrct
        rrx = 0.
        if (jreac(nrc).eq.0 .or. jreac(nrc).eq.-1) then
          rrx = rrelr0(nrc)
          dxp = 1.
          do n = 1,nord
            dxp = dxp*delxi
            rrx = rrx + ( drer0(n,nrc)/fctrl(n) )*dxp
          enddo
        endif
        rrelrp(nrc) = rrx
      enddo
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c     Find any cross-overs. Note that there are two kinds, positive to
c     negative, and negative to positive.
c
      xval = 0.
      nrzero = 0
      krzero = 0
      unam24 = ' '
      do nrc = 1,nrct
        if (jreac(nrc).eq.0 .or. jreac(nrc).eq.-1) then
          rrx0 = rrelr0(nrc)
          rrxp = rrelrp(nrc)
          if (rrx0 .gt. 0.) then
            if (rrxp .lt. -tolsrr) then
c
c             A case has been found in which the cross over tolerance
c             is exceeded for the positive to negative case.
c
              krzero = krzero + 1
              arrxp = -rrxp
              if (arrxp .gt. xval) then
                xval = arrxp
                nrzero = nrc
                unam24 = ureac(nrc)
                ilsign = 1
                xtargv = -atgsrr
              endif
            endif
          elseif (rrx0 .lt. 0.) then
            if (rrxp .gt. tolsrr) then
c
c             A case has been found in which the cross over tolerance
c             is exceeded for the negative to positive case.
c
              krzero = krzero + 1
              arrxp = rrxp
              if (arrxp .gt. xval) then
                xval = arrxp
                nrzero = nrc
                unam24 = ureac(nrc)
                ilsign = -1
                xtargv = atgsrr
              endif
            endif
          endif
        endif
      enddo
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
      if (krzero .gt. 0) then
        if (delxi .gt. dlxmin) then
          if (abs(xval - xtargv) .gt. tolsx) then
c
            xval0 = rrelr0(nrzero)
            do n = 1,nord
              dxval0(n) = drer0(n,nrzero)
            enddo
c
            call search(delxi,dlxmin,dxval0,eps100,ier,ilsign,iodb,
     $      nodbmx,nord,noutpt,nrd1mx,nttyo,tolsx,unam24,usearch,
     $      xtargv,xval0)
c
            if (ier .le. 0) go to 100
            if (ier .ge. 2) then
c
c             Note: if ier = 1, the returned "safe" value of delxi
c             is used.
c
              delxi = dlxmin
            endif
            go to 999
          endif
        endif
      endif
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
      if (iodb(1).gt.0 .or. iodb(5).gt.0) then
        if (krzero .gt. 0) then
          if (abs(xval - xtargv) .le. tolsx) then
            j2 = ilnobl(unam24)
            write (noutpt,1010) unam24(1:j2),xi1,delxi
 1010       format(/" --- Taylor's series predict a zero relative",
     $      /7x,'rate for ',a,' at Xi= ',1pe11.4,',',
     $      /7x,'delxi= ',1pe11.4,' ---')
          endif
        endif
      endif
c
c* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
  999 continue
      end


text:
C$Procedure            NCPOS ( NOT character position )

      INTEGER FUNCTION NCPOS ( STR, CHARS, START )

C$ Abstract
C
C     Find the first occurrence in a string of a character NOT belonging
C     to a collection of characters, starting at a specified location,
C     searching forwards.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     SCANNING
C
C$ Keywords
C
C     CHARACTER
C     SEARCH
C     UTILITY
C
C$ Declarations

      CHARACTER*(*)         STR
      CHARACTER*(*)         CHARS
      INTEGER               START

C$ Brief_I/O
C
C     VARIABLE  I/O  DESCRIPTION
C     --------  ---  --------------------------------------------------
C     STR        I   Any character string.
C     CHARS      I   A collection of characters.
C     START      I   Position to begin looking for one not in CHARS
C
C     The function returns the index of the first character of STR
C     at or following index START that is not in the collection CHARS.
C
C$ Detailed_Input
C
C     STR        is any character string.
C
C     CHARS      is a character string containing a collection of
C                characters.  Spaces in CHARS are significant.
C
C     START      is the position in STR to begin looking for
C                characters not in CHARS.
C
C$ Detailed_Output
C
C     The function returns the index of the first character of STR (at
C     or following index START) that is not one of the characters in the
C     string CHARS.  If no such character is found, the function returns
C     zero.
C
C$ Parameters
C
C     None.
C
C$ Exceptions
C
C     Error Free.
C
C     1) If START is less than 1, the search begins at the first
C        character of the string.
C
C     2) If START is greater than the length of the string, NCPOS
C        returns zero.
C
C$ Files
C
C     None.
C
C$ Particulars
C
C     NCPOS is case sensitive.
C
C     An entire family of related SPICELIB routines (POS, CPOS, NCPOS,
C     POSR, CPOSR, NCPOSR) is described in the Required Reading.
C
C     Those familiar with the True BASIC language should note that
C     these functions are equivalent to the True BASIC intrinsic
C     functions with the same names.
C
C$ Examples
C
C     Let STRING = 'BOB, JOHN, TED, AND MARTIN    '
C                   123456789012345678901234567890
C
C     Let CHAR   = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
C
C     Normal (Sequential) Searching:
C     ------------------------------
C
C           NCPOS( STRING, CHAR,  1 ) = 4
C           NCPOS( STRING, CHAR,  5 ) = 5
C           NCPOS( STRING, CHAR,  6 ) = 10
C           NCPOS( STRING, CHAR, 11 ) = 11
C           NCPOS( STRING, CHAR, 12 ) = 15
C           NCPOS( STRING, CHAR, 16 ) = 16
C           NCPOS( STRING, CHAR, 17 ) = 20
C           NCPOS( STRING, CHAR, 21 ) = 27
C           NCPOS( STRING, CHAR, 28 ) = 28
C           NCPOS( STRING, CHAR, 29 ) = 29
C           NCPOS( STRING, CHAR, 30 ) = 30
C           NCPOS( STRING, CHAR, 31 ) =  0
C
C     START out of bounds:
C     --------------------
C
C           NCPOS( STRING, CHAR, -12 ) = 4
C           NCPOS( STRING, CHAR,   0 ) = 4
C           NCPOS( STRING, CHAR,  31 ) = 0
C           NCPOS( STRING, CHAR, 123 ) = 0
C
C     Order within CHARS:
C     -------------------
C
C           NCPOS( STRING, 'JOHN', 7 ) = 10
C           NCPOS( STRING, 'OHJN', 7 ) = 10
C
C$ Restrictions
C
C     None.
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     W.L. Taber      (JPL)
C     H.A. Neilan     (JPL)
C     B.V. Semenov    (JPL)
C
C$ Version
C
C-    SPICELIB Version 1.0.3, 31-JAN-2008 (BVS)
C
C        Removed non-standard end-of-declarations marker
C        'C%&END_DECLARATIONS' from comments.
C
C-    SPICELIB Version 1.0.2, 10-MAR-1992 (WLT)
C
C        Comment section for permuted index source lines was added
C        following the header.
C
C-    SPICELIB Version 1.0.1, 26-MAR-1991 (HAN)
C
C        The Required Reading file POSITION was renamed to SCANNING.
C        This header was updated to reflect the change.
C
C-    SPICELIB Version 1.0.0, 31-JAN-1990 (WLT)
C
C-&

C$ Index_Entries
C
C     forward search for position of unlisted character
C
C-&

C
C     Local variables
C
      LOGICAL               FOUND
      INTEGER               LENSTR
      INTEGER               B

      LENSTR = LEN ( STR )

      B      = MAX(1,START)
      FOUND  = .FALSE.
      NCPOS   = 0

      DO WHILE ( .NOT. FOUND )

         IF ( B .GT. LENSTR ) THEN
            RETURN
         ELSE IF ( INDEX( CHARS, STR(B:B) ) .EQ. 0 ) THEN
            NCPOS = B
            RETURN
         ELSE
            B = B + 1
         END IF

      END DO

      RETURN
      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
!!  Copyright (C)  Stichting Deltares, 2012-2016.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License version 3,
!!  as published by the Free Software Foundation.
!!
!!  This program is distributed in the hope that it will be useful,
!!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!!  GNU General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!  contact: delft3d.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands
!!
!!  All indications and logos of, and references to registered trademarks
!!  of Stichting Deltares remain the property of Stichting Deltares. All
!!  rights reserved.

!
!  *********************************************************************
!  *         SUBROUTINE TO PERFORM SPECIAL INTERPOLATION               *
!  *********************************************************************
!
      SUBROUTINE EBCALC(X,F,FPR,NUMGR)
      IMPLICIT REAL*8 (A-H,O-Z)
      INCLUDE 'blmdim.inc'
      INCLUDE 'arran.inc'
!
!  CHECK WHETHER X IS TOO LOW OR TOO HIGH
!
      IF (X .GT. ZVEC(1)) GO TO 20
      F=FUN(1,NUMGR)
      FPR=0.0
      RETURN
   20 IF (X .LT. ZVEC(NZ)) GO TO 40
      F=FUN(NZ,NUMGR)
      FPR=0.0
      RETURN
   40 DO 60 I=2,NZ
      IF (X .LE. ZVEC(I)) GO TO 80
   60 CONTINUE
   80 EX=DEXP(-X)
      EI=DEXP(-ZVEC(I))
      EI1=DEXP(-ZVEC(I-1))
      ALAM=(EX-EI)/(EI1-EI)
      FPR=ALAM*DER(I-1,NUMGR)+(1.0-ALAM)*DER(I,NUMGR)
      C0=((EX/EI1)-1.0+X-ZVEC(I-1))/((EX/EI1)-1.0)
      C1=((EI1/EX)-1.0+ZVEC(I-1)-X)/((EI1/EX)-1.0)
      F=FUN(I-1,NUMGR)-C0*FPR+C1*DER(I-1,NUMGR)
      RETURN
      END