text:
*DECK XERCNT
      SUBROUTINE XERCNT (LIBRAR, SUBROU, MESSG, NERR, LEVEL, KONTRL)
C***BEGIN PROLOGUE  XERCNT
C***SUBSIDIARY
C***PURPOSE  Allow user control over handling of errors.
C***LIBRARY   SLATEC (XERROR)
C***CATEGORY  R3C
C***TYPE      ALL (XERCNT-A)
C***KEYWORDS  ERROR, XERROR
C***AUTHOR  Jones, R. E., (SNLA)
C***DESCRIPTION
C
C     Abstract
C        Allows user control over handling of individual errors.
C        Just after each message is recorded, but before it is
C        processed any further (i.e., before it is printed or
C        a decision to abort is made), a call is made to XERCNT.
C        If the user has provided his own version of XERCNT, he
C        can then override the value of KONTROL used in processing
C        this message by redefining its value.
C        KONTRL may be set to any value from -2 to 2.
C        The meanings for KONTRL are the same as in XSETF, except
C        that the value of KONTRL changes only for this message.
C        If KONTRL is set to a value outside the range from -2 to 2,
C        it will be moved back into that range.
C
C     Description of Parameters
C
C      --Input--
C        LIBRAR - the library that the routine is in.
C        SUBROU - the subroutine that XERMSG is being called from
C        MESSG  - the first 20 characters of the error message.
C        NERR   - same as in the call to XERMSG.
C        LEVEL  - same as in the call to XERMSG.
C        KONTRL - the current value of the control flag as set
C                 by a call to XSETF.
C
C      --Output--
C        KONTRL - the new value of KONTRL.  If KONTRL is not
C                 defined, it will remain at its original value.
C                 This changed value of control affects only
C                 the current occurrence of the current message.
C
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  (NONE)
C***REVISION HISTORY  (YYMMDD)
C   790801  DATE WRITTEN
C   861211  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900206  Routine changed from user-callable to subsidiary.  (WRB)
C   900510  Changed calling sequence to include LIBRARY and SUBROUTINE
C           names, changed routine name from XERCTL to XERCNT.  (RWC)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  XERCNT
      CHARACTER*(*) LIBRAR, SUBROU, MESSG
C***FIRST EXECUTABLE STATEMENT  XERCNT
      RETURN
      END
*DECK XERHLT
      SUBROUTINE XERHLT (MESSG)
C***BEGIN PROLOGUE  XERHLT
C***SUBSIDIARY
C***PURPOSE  Abort program execution and print error message.
C***LIBRARY   SLATEC (XERROR)
C***CATEGORY  R3C
C***TYPE      ALL (XERHLT-A)
C***KEYWORDS  ABORT PROGRAM EXECUTION, ERROR, XERROR
C***AUTHOR  Jones, R. E., (SNLA)
C***DESCRIPTION
C
C     Abstract
C        ***Note*** machine dependent routine
C        XERHLT aborts the execution of the program.
C        The error message causing the abort is given in the calling
C        sequence, in case one needs it for printing on a dayfile,
C        for example.
C
C     Description of Parameters
C        MESSG is as in XERMSG.
C
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  (NONE)
C***REVISION HISTORY  (YYMMDD)
C   790801  DATE WRITTEN
C   861211  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900206  Routine changed from user-callable to subsidiary.  (WRB)
C   900510  Changed calling sequence to delete length of character
C           and changed routine name from XERABT to XERHLT.  (RWC)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  XERHLT
      CHARACTER*(*) MESSG
C***FIRST EXECUTABLE STATEMENT  XERHLT
      STOP
      END
*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
*DECK XERPRN
      SUBROUTINE XERPRN (PREFIX, NPREF, MESSG, NWRAP)
C***BEGIN PROLOGUE  XERPRN
C***SUBSIDIARY
C***PURPOSE  Print error messages processed by XERMSG.
C***LIBRARY   SLATEC (XERROR)
C***CATEGORY  R3C
C***TYPE      ALL (XERPRN-A)
C***KEYWORDS  ERROR MESSAGES, PRINTING, XERROR
C***AUTHOR  Fong, Kirby, (NMFECC at LLNL)
C***DESCRIPTION
C
C This routine sends one or more lines to each of the (up to five)
C logical units to which error messages are to be sent.  This routine
C is called several times by XERMSG, sometimes with a single line to
C print and sometimes with a (potentially very long) message that may
C wrap around into multiple lines.
C
C PREFIX  Input argument of type CHARACTER.  This argument contains
C         characters to be put at the beginning of each line before
C         the body of the message.  No more than 16 characters of
C         PREFIX will be used.
C
C NPREF   Input argument of type INTEGER.  This argument is the number
C         of characters to use from PREFIX.  If it is negative, the
C         intrinsic function LEN is used to determine its length.  If
C         it is zero, PREFIX is not used.  If it exceeds 16 or if
C         LEN(PREFIX) exceeds 16, only the first 16 characters will be
C         used.  If NPREF is positive and the length of PREFIX is less
C         than NPREF, a copy of PREFIX extended with blanks to length
C         NPREF will be used.
C
C MESSG   Input argument of type CHARACTER.  This is the text of a
C         message to be printed.  If it is a long message, it will be
C         broken into pieces for printing on multiple lines.  Each line
C         will start with the appropriate prefix and be followed by a
C         piece of the message.  NWRAP is the number of characters per
C         piece; that is, after each NWRAP characters, we break and
C         start a new line.  In addition the characters '$$' embedded
C         in MESSG are a sentinel for a new line.  The counting of
C         characters up to NWRAP starts over for each new line.  The
C         value of NWRAP typically used by XERMSG is 72 since many
C         older error messages in the SLATEC Library are laid out to
C         rely on wrap-around every 72 characters.
C
C NWRAP   Input argument of type INTEGER.  This gives the maximum size
C         piece into which to break MESSG for printing on multiple
C         lines.  An embedded '$$' ends a line, and the count restarts
C         at the following character.  If a line break does not occur
C         on a blank (it would split a word) that word is moved to the
C         next line.  Values of NWRAP less than 16 will be treated as
C         16.  Values of NWRAP greater than 132 will be treated as 132.
C         The actual line length will be NPREF + NWRAP after NPREF has
C         been adjusted to fall between 0 and 16 and NWRAP has been
C         adjusted to fall between 16 and 132.
C
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  I1MACH, XGETUA
C***REVISION HISTORY  (YYMMDD)
C   880621  DATE WRITTEN
C   880708  REVISED AFTER THE SLATEC CML SUBCOMMITTEE MEETING OF
C           JUNE 29 AND 30 TO CHANGE THE NAME TO XERPRN AND TO REWORK
C           THE HANDLING OF THE NEW LINE SENTINEL TO BEHAVE LIKE THE
C           SLASH CHARACTER IN FORMAT STATEMENTS.
C   890706  REVISED WITH THE HELP OF FRED FRITSCH AND REG CLEMENS TO
C           STREAMLINE THE CODING AND FIX A BUG THAT CAUSED EXTRA BLANK
C           LINES TO BE PRINTED.
C   890721  REVISED TO ADD A NEW FEATURE.  A NEGATIVE VALUE OF NPREF
C           CAUSES LEN(PREFIX) TO BE USED AS THE LENGTH.
C   891013  REVISED TO CORRECT ERROR IN CALCULATING PREFIX LENGTH.
C   891214  Prologue converted to Version 4.0 format.  (WRB)
C   900510  Added code to break messages between words.  (RWC)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  XERPRN
      CHARACTER*(*) PREFIX, MESSG
      INTEGER NPREF, NWRAP
      CHARACTER*148 CBUFF
      INTEGER IU(5), NUNIT
      CHARACTER*2 NEWLIN
      PARAMETER (NEWLIN = '$$')
C***FIRST EXECUTABLE STATEMENT  XERPRN
      CALL XGETUA(IU,NUNIT)
C
C       A ZERO VALUE FOR A LOGICAL UNIT NUMBER MEANS TO USE THE STANDARD
C       ERROR MESSAGE UNIT INSTEAD.  I1MACH(4) RETRIEVES THE STANDARD
C       ERROR MESSAGE UNIT.
C
      N = I1MACH(4)
      DO 10 I=1,NUNIT
         IF (IU(I) .EQ. 0) IU(I) = N
   10 CONTINUE
C
C       LPREF IS THE LENGTH OF THE PREFIX.  THE PREFIX IS PLACED AT THE
C       BEGINNING OF CBUFF, THE CHARACTER BUFFER, AND KEPT THERE DURING
C       THE REST OF THIS ROUTINE.
C
      IF ( NPREF .LT. 0 ) THEN
         LPREF = LEN(PREFIX)
      ELSE
         LPREF = NPREF
      ENDIF
      LPREF = MIN(16, LPREF)
      IF (LPREF .NE. 0) CBUFF(1:LPREF) = PREFIX
C
C       LWRAP IS THE MAXIMUM NUMBER OF CHARACTERS WE WANT TO TAKE AT ONE
C       TIME FROM MESSG TO PRINT ON ONE LINE.
C
      LWRAP = MAX(16, MIN(132, NWRAP))
C
C       SET LENMSG TO THE LENGTH OF MESSG, IGNORE ANY TRAILING BLANKS.
C
      LENMSG = LEN(MESSG)
      N = LENMSG
      DO 20 I=1,N
         IF (MESSG(LENMSG:LENMSG) .NE. ' ') GO TO 30
         LENMSG = LENMSG - 1
   20 CONTINUE
   30 CONTINUE
C
C       IF THE MESSAGE IS ALL BLANKS, THEN PRINT ONE BLANK LINE.
C
      IF (LENMSG .EQ. 0) THEN
         CBUFF(LPREF+1:LPREF+1) = ' '
         DO 40 I=1,NUNIT
            WRITE(IU(I), '(A)') CBUFF(1:LPREF+1)
   40    CONTINUE
         RETURN
      ENDIF
C
C       SET NEXTC TO THE POSITION IN MESSG WHERE THE NEXT SUBSTRING
C       STARTS.  FROM THIS POSITION WE SCAN FOR THE NEW LINE SENTINEL.
C       WHEN NEXTC EXCEEDS LENMSG, THERE IS NO MORE TO PRINT.
C       WE LOOP BACK TO LABEL 50 UNTIL ALL PIECES HAVE BEEN PRINTED.
C
C       WE LOOK FOR THE NEXT OCCURRENCE OF THE NEW LINE SENTINEL.  THE
C       INDEX INTRINSIC FUNCTION RETURNS ZERO IF THERE IS NO OCCURRENCE
C       OR IF THE LENGTH OF THE FIRST ARGUMENT IS LESS THAN THE LENGTH
C       OF THE SECOND ARGUMENT.
C
C       THERE ARE SEVERAL CASES WHICH SHOULD BE CHECKED FOR IN THE
C       FOLLOWING ORDER.  WE ARE ATTEMPTING TO SET LPIECE TO THE NUMBER
C       OF CHARACTERS THAT SHOULD BE TAKEN FROM MESSG STARTING AT
C       POSITION NEXTC.
C
C       LPIECE .EQ. 0   THE NEW LINE SENTINEL DOES NOT OCCUR IN THE
C                       REMAINDER OF THE CHARACTER STRING.  LPIECE
C                       SHOULD BE SET TO LWRAP OR LENMSG+1-NEXTC,
C                       WHICHEVER IS LESS.
C
C       LPIECE .EQ. 1   THE NEW LINE SENTINEL STARTS AT MESSG(NEXTC:
C                       NEXTC).  LPIECE IS EFFECTIVELY ZERO, AND WE
C                       PRINT NOTHING TO AVOID PRODUCING UNNECESSARY
C                       BLANK LINES.  THIS TAKES CARE OF THE SITUATION
C                       WHERE THE LIBRARY ROUTINE HAS A MESSAGE OF
C                       EXACTLY 72 CHARACTERS FOLLOWED BY A NEW LINE
C                       SENTINEL FOLLOWED BY MORE CHARACTERS.  NEXTC
C                       SHOULD BE INCREMENTED BY 2.
C
C       LPIECE .GT. LWRAP+1  REDUCE LPIECE TO LWRAP.
C
C       ELSE            THIS LAST CASE MEANS 2 .LE. LPIECE .LE. LWRAP+1
C                       RESET LPIECE = LPIECE-1.  NOTE THAT THIS
C                       PROPERLY HANDLES THE END CASE WHERE LPIECE .EQ.
C                       LWRAP+1.  THAT IS, THE SENTINEL FALLS EXACTLY
C                       AT THE END OF A LINE.
C
      NEXTC = 1
   50 LPIECE = INDEX(MESSG(NEXTC:LENMSG), NEWLIN)
      IF (LPIECE .EQ. 0) THEN
C
C       THERE WAS NO NEW LINE SENTINEL FOUND.
C
         IDELTA = 0
         LPIECE = MIN(LWRAP, LENMSG+1-NEXTC)
         IF (LPIECE .LT. LENMSG+1-NEXTC) THEN
            DO 52 I=LPIECE+1,2,-1
               IF (MESSG(NEXTC+I-1:NEXTC+I-1) .EQ. ' ') THEN
                  LPIECE = I-1
                  IDELTA = 1
                  GOTO 54
               ENDIF
   52       CONTINUE
         ENDIF
   54    CBUFF(LPREF+1:LPREF+LPIECE) = MESSG(NEXTC:NEXTC+LPIECE-1)
         NEXTC = NEXTC + LPIECE + IDELTA
      ELSEIF (LPIECE .EQ. 1) THEN
C
C       WE HAVE A NEW LINE SENTINEL AT MESSG(NEXTC:NEXTC+1).
C       DON'T PRINT A BLANK LINE.
C
         NEXTC = NEXTC + 2
         GO TO 50
      ELSEIF (LPIECE .GT. LWRAP+1) THEN
C
C       LPIECE SHOULD BE SET DOWN TO LWRAP.
C
         IDELTA = 0
         LPIECE = LWRAP
         DO 56 I=LPIECE+1,2,-1
            IF (MESSG(NEXTC+I-1:NEXTC+I-1) .EQ. ' ') THEN
               LPIECE = I-1
               IDELTA = 1
               GOTO 58
            ENDIF
   56    CONTINUE
   58    CBUFF(LPREF+1:LPREF+LPIECE) = MESSG(NEXTC:NEXTC+LPIECE-1)
         NEXTC = NEXTC + LPIECE + IDELTA
      ELSE
C
C       IF WE ARRIVE HERE, IT MEANS 2 .LE. LPIECE .LE. LWRAP+1.
C       WE SHOULD DECREMENT LPIECE BY ONE.
C
         LPIECE = LPIECE - 1
         CBUFF(LPREF+1:LPREF+LPIECE) = MESSG(NEXTC:NEXTC+LPIECE-1)
         NEXTC  = NEXTC + LPIECE + 2
      ENDIF
C
C       PRINT
C
      DO 60 I=1,NUNIT
         WRITE(IU(I), '(A)') CBUFF(1:LPREF+LPIECE)
   60 CONTINUE
C
      IF (NEXTC .LE. LENMSG) GO TO 50
      RETURN
      END
*DECK XERSVE
      SUBROUTINE XERSVE (LIBRAR, SUBROU, MESSG, KFLAG, NERR, LEVEL,
     +   ICOUNT)
C***BEGIN PROLOGUE  XERSVE
C***SUBSIDIARY
C***PURPOSE  Record that an error has occurred.
C***LIBRARY   SLATEC (XERROR)
C***CATEGORY  R3
C***TYPE      ALL (XERSVE-A)
C***KEYWORDS  ERROR, XERROR
C***AUTHOR  Jones, R. E., (SNLA)
C***DESCRIPTION
C
C *Usage:
C
C        INTEGER  KFLAG, NERR, LEVEL, ICOUNT
C        CHARACTER * (len) LIBRAR, SUBROU, MESSG
C
C        CALL XERSVE (LIBRAR, SUBROU, MESSG, KFLAG, NERR, LEVEL, ICOUNT)
C
C *Arguments:
C
C        LIBRAR :IN    is the library that the message is from.
C        SUBROU :IN    is the subroutine that the message is from.
C        MESSG  :IN    is the message to be saved.
C        KFLAG  :IN    indicates the action to be performed.
C                      when KFLAG > 0, the message in MESSG is saved.
C                      when KFLAG=0 the tables will be dumped and
C                      cleared.
C                      when KFLAG < 0, the tables will be dumped and
C                      not cleared.
C        NERR   :IN    is the error number.
C        LEVEL  :IN    is the error severity.
C        ICOUNT :OUT   the number of times this message has been seen,
C                      or zero if the table has overflowed and does not
C                      contain this message specifically.  When KFLAG=0,
C                      ICOUNT will not be altered.
C
C *Description:
C
C   Record that this error occurred and possibly dump and clear the
C   tables.
C
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  I1MACH, XGETUA
C***REVISION HISTORY  (YYMMDD)
C   800319  DATE WRITTEN
C   861211  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900413  Routine modified to remove reference to KFLAG.  (WRB)
C   900510  Changed to add LIBRARY NAME and SUBROUTINE to calling
C           sequence, use IF-THEN-ELSE, make number of saved entries
C           easily changeable, changed routine name from XERSAV to
C           XERSVE.  (RWC)
C   910626  Added LIBTAB and SUBTAB to SAVE statement.  (BKS)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  XERSVE
      PARAMETER (LENTAB=10)
      INTEGER LUN(5)
      CHARACTER*(*) LIBRAR, SUBROU, MESSG
      CHARACTER*8  LIBTAB(LENTAB), SUBTAB(LENTAB), LIB, SUB
      CHARACTER*20 MESTAB(LENTAB), MES
      DIMENSION NERTAB(LENTAB), LEVTAB(LENTAB), KOUNT(LENTAB)
      SAVE LIBTAB, SUBTAB, MESTAB, NERTAB, LEVTAB, KOUNT, KOUNTX, NMSG
      DATA KOUNTX/0/, NMSG/0/
C***FIRST EXECUTABLE STATEMENT  XERSVE
C
      IF (KFLAG.LE.0) THEN
C
C        Dump the table.
C
         IF (NMSG.EQ.0) RETURN
C
C        Print to each unit.
C
         CALL XGETUA (LUN, NUNIT)
         DO 20 KUNIT = 1,NUNIT
            IUNIT = LUN(KUNIT)
            IF (IUNIT.EQ.0) IUNIT = I1MACH(4)
C
C           Print the table header.
C
            WRITE (IUNIT,9000)
C
C           Print body of table.
C
            DO 10 I = 1,NMSG
               WRITE (IUNIT,9010) LIBTAB(I), SUBTAB(I), MESTAB(I),
     *            NERTAB(I),LEVTAB(I),KOUNT(I)
   10       CONTINUE
C
C           Print number of other errors.
C
            IF (KOUNTX.NE.0) WRITE (IUNIT,9020) KOUNTX
            WRITE (IUNIT,9030)
   20    CONTINUE
C
C        Clear the error tables.
C
         IF (KFLAG.EQ.0) THEN
            NMSG = 0
            KOUNTX = 0
         ENDIF
      ELSE
C
C        PROCESS A MESSAGE...
C        SEARCH FOR THIS MESSG, OR ELSE AN EMPTY SLOT FOR THIS MESSG,
C        OR ELSE DETERMINE THAT THE ERROR TABLE IS FULL.
C
         LIB = LIBRAR
         SUB = SUBROU
         MES = MESSG
         DO 30 I = 1,NMSG
            IF (LIB.EQ.LIBTAB(I) .AND. SUB.EQ.SUBTAB(I) .AND.
     *         MES.EQ.MESTAB(I) .AND. NERR.EQ.NERTAB(I) .AND.
     *         LEVEL.EQ.LEVTAB(I)) THEN
                  KOUNT(I) = KOUNT(I) + 1
                  ICOUNT = KOUNT(I)
                  RETURN
            ENDIF
   30    CONTINUE
C
         IF (NMSG.LT.LENTAB) THEN
C
C           Empty slot found for new message.
C
            NMSG = NMSG + 1
            LIBTAB(I) = LIB
            SUBTAB(I) = SUB
            MESTAB(I) = MES
            NERTAB(I) = NERR
            LEVTAB(I) = LEVEL
            KOUNT (I) = 1
            ICOUNT    = 1
         ELSE
C
C           Table is full.
C
            KOUNTX = KOUNTX+1
            ICOUNT = 0
         ENDIF
      ENDIF
      RETURN
C
C     Formats.
C
 9000 FORMAT ('0          ERROR MESSAGE SUMMARY' /
     +   ' LIBRARY    SUBROUTINE MESSAGE START             NERR',
     +   '     LEVEL     COUNT')
 9010 FORMAT (1X,A,3X,A,3X,A,3I10)
 9020 FORMAT ('0OTHER ERRORS NOT INDIVIDUALLY TABULATED = ', I10)
 9030 FORMAT (1X)
      END
*DECK XGETUA
      SUBROUTINE XGETUA (IUNITA, N)
C***BEGIN PROLOGUE  XGETUA
C***PURPOSE  Return unit number(s) to which error messages are being
C            sent.
C***LIBRARY   SLATEC (XERROR)
C***CATEGORY  R3C
C***TYPE      ALL (XGETUA-A)
C***KEYWORDS  ERROR, XERROR
C***AUTHOR  Jones, R. E., (SNLA)
C***DESCRIPTION
C
C     Abstract
C        XGETUA may be called to determine the unit number or numbers
C        to which error messages are being sent.
C        These unit numbers may have been set by a call to XSETUN,
C        or a call to XSETUA, or may be a default value.
C
C     Description of Parameters
C      --Output--
C        IUNIT - an array of one to five unit numbers, depending
C                on the value of N.  A value of zero refers to the
C                default unit, as defined by the I1MACH machine
C                constant routine.  Only IUNIT(1),...,IUNIT(N) are
C                defined by XGETUA.  The values of IUNIT(N+1),...,
C                IUNIT(5) are not defined (for N .LT. 5) or altered
C                in any way by XGETUA.
C        N     - the number of units to which copies of the
C                error messages are being sent.  N will be in the
C                range from 1 to 5.
C
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  J4SAVE
C***REVISION HISTORY  (YYMMDD)
C   790801  DATE WRITTEN
C   861211  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  XGETUA
      DIMENSION IUNITA(5)
C***FIRST EXECUTABLE STATEMENT  XGETUA
      N = J4SAVE(5,0,.FALSE.)
      DO 30 I=1,N
         INDEX = I+4
         IF (I.EQ.1) INDEX = 3
         IUNITA(I) = J4SAVE(INDEX,0,.FALSE.)
   30 CONTINUE
      RETURN
      END
*DECK FDUMP
      SUBROUTINE FDUMP
C***BEGIN PROLOGUE  FDUMP
C***PURPOSE  Symbolic dump (should be locally written).
C***LIBRARY   SLATEC (XERROR)
C***CATEGORY  R3
C***TYPE      ALL (FDUMP-A)
C***KEYWORDS  ERROR, XERMSG
C***AUTHOR  Jones, R. E., (SNLA)
C***DESCRIPTION
C
C        ***Note*** Machine Dependent Routine
C        FDUMP is intended to be replaced by a locally written
C        version which produces a symbolic dump.  Failing this,
C        it should be replaced by a version which prints the
C        subprogram nesting list.  Note that this dump must be
C        printed on each of up to five files, as indicated by the
C        XGETUA routine.  See XSETUA and XGETUA for details.
C
C     Written by Ron Jones, with SLATEC Common Math Library Subcommittee
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  (NONE)
C***REVISION HISTORY  (YYMMDD)
C   790801  DATE WRITTEN
C   861211  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C***END PROLOGUE  FDUMP
C***FIRST EXECUTABLE STATEMENT  FDUMP
      RETURN
      END
*DECK J4SAVE
      FUNCTION J4SAVE (IWHICH, IVALUE, ISET)
C***BEGIN PROLOGUE  J4SAVE
C***SUBSIDIARY
C***PURPOSE  Save or recall global variables needed by error
C            handling routines.
C***LIBRARY   SLATEC (XERROR)
C***TYPE      INTEGER (J4SAVE-I)
C***KEYWORDS  ERROR MESSAGES, ERROR NUMBER, RECALL, SAVE, XERROR
C***AUTHOR  Jones, R. E., (SNLA)
C***DESCRIPTION
C
C     Abstract
C        J4SAVE saves and recalls several global variables needed
C        by the library error handling routines.
C
C     Description of Parameters
C      --Input--
C        IWHICH - Index of item desired.
C                = 1 Refers to current error number.
C                = 2 Refers to current error control flag.
C                = 3 Refers to current unit number to which error
C                    messages are to be sent.  (0 means use standard.)
C                = 4 Refers to the maximum number of times any
C                     message is to be printed (as set by XERMAX).
C                = 5 Refers to the total number of units to which
C                     each error message is to be written.
C                = 6 Refers to the 2nd unit for error messages
C                = 7 Refers to the 3rd unit for error messages
C                = 8 Refers to the 4th unit for error messages
C                = 9 Refers to the 5th unit for error messages
C        IVALUE - The value to be set for the IWHICH-th parameter,
C                 if ISET is .TRUE. .
C        ISET   - If ISET=.TRUE., the IWHICH-th parameter will BE
C                 given the value, IVALUE.  If ISET=.FALSE., the
C                 IWHICH-th parameter will be unchanged, and IVALUE
C                 is a dummy parameter.
C      --Output--
C        The (old) value of the IWHICH-th parameter will be returned
C        in the function value, J4SAVE.
C
C***SEE ALSO  XERMSG
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  (NONE)
C***REVISION HISTORY  (YYMMDD)
C   790801  DATE WRITTEN
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900205  Minor modifications to prologue.  (WRB)
C   900402  Added TYPE section.  (WRB)
C   910411  Added KEYWORDS section.  (WRB)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  J4SAVE
      LOGICAL ISET
      INTEGER IPARAM(9)
      SAVE IPARAM
      DATA IPARAM(1),IPARAM(2),IPARAM(3),IPARAM(4)/0,2,0,10/
      DATA IPARAM(5)/1/
      DATA IPARAM(6),IPARAM(7),IPARAM(8),IPARAM(9)/0,0,0,0/
C***FIRST EXECUTABLE STATEMENT  J4SAVE
      J4SAVE = IPARAM(IWHICH)
      IF (ISET) IPARAM(IWHICH) = IVALUE
      RETURN
      END

max_issues_repo_name: bslakman/RMG-Java

max_forks_repo_name: bslakman/RMG-Java


text:
module rojff_string_sink_m
    implicit none
    private
    public :: string_sink_t

    type, abstract :: string_sink_t
    contains
        procedure(append_i), deferred :: append
    end type

    abstract interface
        subroutine append_i(self, part)
            import :: string_sink_t

            implicit none

            class(string_sink_t), intent(inout) :: self
            character(len=*), intent(in) :: part
        end subroutine
    end interface
end module

max_stars_repo_name: everythingfunctional/rojff

max_issues_repo_name: everythingfunctional/rojff

max_forks_repo_name: everythingfunctional/rojff


text:
!*******************************************************************
!                                                                  *
      SUBROUTINE RMEAJJ(LL,IT,LQ,J,ITS,LQS,J1S,COEF)
!                                                                  *
!   Written by  G. Gaigalas                                        *
!   Transform to fortran 90/95 by G. Gaigalas       December 2012  *
!   The last modification made by G. Gaigalas       October  2017  *
!                                                                  *
!*******************************************************************
!
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY: DOUBLE
      USE CONS_C,          ONLY: ZERO
      USE ribojj_C
!-----------------------------------------------
!   I n t e r f a c e   B l o c k s
!-----------------------------------------------
      USE rmeajj11_I
      USE rmeajj9_I
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER,      INTENT(IN)  :: LL, IT, LQ, J, ITS, LQS, J1S
      REAL(DOUBLE), INTENT(OUT) :: COEF
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: I1, I2, IE1,  JT, JTS
      INTEGER, DIMENSION(1,2) :: IDS3, IDV3
      INTEGER, DIMENSION(3,3) :: IDS5, IDV5
      INTEGER, DIMENSION(6,8) :: IDS7, IDV7
!-----------------------------------------------
      DATA IDS3/12,20/
      DATA IDV3/1,1/
      DATA IDS5/-24,2*0,-30,120,-90,-54,-48,330/
      DATA IDV5/4*1,2*7,1,2*7/
      DATA IDS7/2*0,-40,3*0,-54,-33,-40,-65,30,0,88,-1,0,-1755,-40,  &
      0,198,-36,-72,4500,-234,-360,-78,156,0,108,-30,96,66,49,0,27,  &
      130,136,0,195,-104,-245,-624,1224,3*0,2720,-2448,-7752/
      DATA IDV7/6*1,7,2*1,7,2*1,7,2*1,77,11,1,7,11,1,77,2*11,35,5,1, &
      91,1,13,2*5,2*1,2*7,1,11,1,3*11,3*1,143,77,91/
!
      COEF=ZERO
      IF(LL > 37) RETURN
      IF(LL <= 9) THEN
        IF(IT < 300) THEN
          IF(IMPTJJ(IT) /= IMPNJJ(ITS)) RETURN
          IF(IT < ITS) THEN
            JT=IT
            JTS=ITS
          ELSE
            JT=ITS
            JTS=IT
          ENDIF
        ENDIF
      ENDIF
      IF(LL == 1) THEN
        COEF=-2
      ELSEIF(LL == 3)THEN
        I1=JT-2
        I2=JTS-3
        COEF=-DSQRT(DBLE(IDS3(I1,I2))/DBLE(IDV3(I1,I2)))
      ELSEIF(LL == 5)THEN
        I1=JT-5
        I2=JTS-8
        IF(IDS5(I1,I2) >= 0) THEN
          COEF=DSQRT(DBLE(IDS5(I1,I2))/DBLE(IDV5(I1,I2)))
        ELSE
          COEF=-DSQRT(-DBLE(IDS5(I1,I2))/DBLE(IDV5(I1,I2)))
        ENDIF
      ELSEIF(LL == 7)THEN
        I1=JT-11
        I2=JTS-17
        IF(IDS7(I1,I2) >= 0) THEN
          COEF=DSQRT(DBLE(IDS7(I1,I2))/DBLE(IDV7(I1,I2)))
        ELSE
          COEF=-DSQRT(-DBLE(IDS7(I1,I2))/DBLE(IDV7(I1,I2)))
        ENDIF
      ELSEIF(LL == 9) THEN
        IF(IT > 300) THEN
          IF(ITS < 300) THEN
            WRITE(6,'(A)') ' ERROR IN  RMEAJJ '
            STOP
          ENDIF
          IF(LL == J1S) THEN
            CALL RMEAJJ11(IT,ITS,LL,COEF)
          ELSE
            CALL RMEAJJ11(ITS,IT,LL,COEF)
            IE1=LQ-LQS+J-J1S+LL-1
            IF((IE1/4)*4 /= IE1)COEF=-COEF
          ENDIF
        ELSE
          CALL RMEAJJ9(IT,LQ,J,ITS,LQS,J1S,COEF)
          IF(IT > ITS) THEN
            IF(MOD(LQ+J-LQS-J1S+LL-1,4) /= 0) COEF=-COEF
          ENDIF
          WRITE(0,'(A)') ' KLAIDA RMEAJJ SUB. '
          STOP
        END IF
      ELSE
        IF(LL == J1S) THEN
          CALL RMEAJJ11(IT,ITS,LL,COEF)
        ELSE
          CALL RMEAJJ11(ITS,IT,LL,COEF)
          IE1=LQ-LQS+J-J1S+LL-1
          IF((IE1/4)*4 /= IE1)COEF=-COEF
        END IF
      ENDIF
      IF(LL < 9) THEN
        IF(IT > ITS) THEN
          IF(MOD(LQ+J-LQS-J1S+LL-1,4) /= 0) COEF=-COEF
        ENDIF
      ENDIF
      RETURN
      END SUBROUTINE RMEAJJ

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: sylas/grasp-continuum

max_forks_repo_name: sylas/grasp-continuum


text:
c     DRDINT1F
c>> 2000-12-03 DRDINT1F Krogh  Declared WORK as WORK(1) in DINTF.
c>> 1994-11-02 DRDINT1F Krogh  Changes to use M77CON
c>> 1994-08-08 DRDINT1F Snyder Took '0' out of formats for C conversion
c>> 1993-05-05 DRDINT1F Krogh  Adjusted to simplify conversion to C.
c>> 1987-12-09 DRDINT1F Snyder Initial Code.
c--D replaces "?": DR?INT1f, ?INT1, ?INTF
c
c     DEMO DRIVER for 1 dimensional quadrature subprogram DINT1.
c
c     Compute the integral for X = 0.0 to X = PI of SIN(X), then
c     subtract 2.0 from the ANSWER.  The result should be zero.
c
c     The integrand is evaluated by the subprogram DINTF.
c
c
      double precision A, B, ANSWER, WORK(1)
      integer IOPT(10)
10    format (/
     1' DRDINT1F:'/
     2' Compute the integral for X = 0.0 to X = PI of SIN(X), then'/
     3' subtract 2.0 from the ANSWER.  The result should be zero.')
20    format (/' ANSWER          =',G15.8/
     1         ' ERROR ESTIMATE  =',G15.8/
     2         ' STATUS FLAG     =',I6/
     3         ' FUNCTION VALUES =',I6)
c
      print 10
      A = 0.0D0
      B = 4.0D0 * ATAN(1.0D0)
      IOPT(2) = 10
      IOPT(3) = 0
      IOPT(4) = 0
c
      call DINT1 (A, B, ANSWER, WORK, IOPT)
      ANSWER = ANSWER - 2.0d0
      print 20, ANSWER, WORK(1), IOPT(1), IOPT(3)
      stop
      end
c   End of DRDINT1F

      subroutine DINTF (ANSWER, WORK, IFLAG)
c
c     Subroutine to provide integrand for DINT1.
c
      double precision ANSWER, WORK(1)
      integer IFLAG
c
c     IFLAG is not used in this example.
c
      answer = sin(work(1))
      return
c
      end

ext: f


text:
! Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
! See https://llvm.org/LICENSE.txt for license information.
! SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
!
!  array constructors of character and assumed-length character
!
module stopwatch

 interface watch
  module procedure watch_a, watch_s
 end interface

 integer n
 integer result(100)
contains

 subroutine create()
  call watch(clock=(/"cpu ","user","sys ","wall"/))
 end subroutine

 subroutine watch_a(clock)
  character(len=*), intent(in), dimension(:) :: clock
  !print *,'in watch_a, character len is ', len(clock(1))
  n = n + 1
  result(n) = len(clock(1))
  do i = 1,ubound(clock,1)
   !print *, i, '->', clock(i)
   n = n + 1
   result(n) = iachar(clock(i)(1:1))
  enddo
 end subroutine

 subroutine watch_s(clock)
  character(len=*), optional, intent(in) :: clock
  if (present(clock)) then
   call watch_a((/clock/))
  else
   call watch_a((/"none"/))
  end if
 end subroutine
end module

 use stopwatch
 integer expect(15)
 data expect/4,99,117,115,119,4,110,10,117,4,100,98,97,102,-1/

 n = 0
 call create
 call watch_s()
 call watch_s('ugly clock')
 call watch_a( (/ 'dont','be','a','fool' /) )
 !print *,n
 !print *,result(1:n)
 n = n + 1
 result(n) = -1
 call check(result,expect,15)
end

max_stars_repo_name: abrahamtovarmob/flang

max_issues_repo_name: abrahamtovarmob/flang

max_forks_repo_name: abrahamtovarmob/flang


text:
program regex_test
use iso_fortran_env
use iso_c_binding
use regex_mod
implicit none

  character(len=:),allocatable :: string
  character(len=:),allocatable :: string2

  character(len=:),allocatable :: my_match
  character(len=:),allocatable :: my_regex

  integer(4) :: i

  type(re_posix) :: re
  type(re_match) :: match_scalar
  type(re_match),allocatable :: match_array(:)

  string = "asdf ghjkl asdf ghjkl"
  string2 = "assssdf ghjkkkkl asddddf ghjjjjkl"

  write(*,*) "scan example:"
  write(*,*) scan(string,"aj")
  write(*,*) scan_all(string,"aj")
  write(*,*) string
  write(*,*) ( merge(string(i:i)," ", any(scan_all(string,"aj")==i) ) ,i = 1, len(string))

  write(*,*) "---------------------------------------------"

  write(*,*) "verify example:"
  write(*,*) verify(string,"adhjk")
  write(*,*) verify_all(string,"adhjk")
  write(*,*) string
  write(*,*) ( merge(string(i:i)," ", any(verify_all(string,"adhjk")==i) ) ,i = 1, len(string))

  write(*,*) "---------------------------------------------"

  write(*,*) "index example:"
  write(*,*) index(string,"f g")
  write(*,*) index_all(string,"f g")
  write(*,*) string
  write(*,*) ( merge(string(i:i)," ", any(index_all(string,"f g")==i) ) ,i = 1, len(string))

  write(*,*) "---------------------------------------------"

  write(*,*) "replace example:"
  write(*,*) string
  write(*,*) replace(string,"f g","12 345")
  write(*,*) replace_all(string,"f g","123 45")

  write(*,*) "---------------------------------------------"

  write(*,*) "regex crete example"
  re = re_posix("[skdj]+","xn")
  write(*,*) re

  write(*,*) "---------------------------------------------"

  write(*,*) "regex match example:"
  write(*,*) string
  match_scalar = re%match(string)
  write(*,*) match_scalar

  match_scalar = match(re,string)
  write(*,*) match_scalar

  match_scalar = match(string, "[skdj]+","xn")
  write(*,*) match_scalar

  write(*,*) "---------------------------------------------"

  write(*,*) "regex match_N example:"
  write(*,*) string
  match_array = re%match_N(string, 5)
  write(*,*) match_array

  match_array = match_N(re,string, 5)
  write(*,*) match_array

  match_array = match_N(string, "[skdj]+",5,"xn")
  write(*,*) match_array

  write(*,*) "---------------------------------------------"

  write(*,*) "regex match_all example:"
  write(*,*) string
  match_array = re%match_all(string)
  write(*,*) match_array

  match_array = match_all(re,string)
  write(*,*) match_array

  match_array = match_all(string, "[skdj]+","xn")
  write(*,*) match_array

  write(*,*) "---------------------------------------------"

  write(*,*) "regex replace example:"
  write(*,*) string2
  write(*,*) re%replace(string2, "<......>")
  write(*,*) re_replace(re, string2, "<......>")
  write(*,*) re_replace(string2, "[skdj]+","<......>","xn")

  write(*,*) "---------------------------------------------"

  write(*,*) "regex replace_all example:"
  write(*,*) string2
  write(*,*) re%replace_all(string2, "<......>")
  write(*,*) re_replace_all(string2, "[skdj]+","<......>","xn")
  write(*,*) re_replace_all(re,string2, "<......>")

  write(*,*) "---------------------------------------------"

  write(*,*) "regex assignment (=) example:"
  re = re_posix("[skdj]+","xn")
  write(*,*) re
  write(*,*) string
  match_scalar = re%match(string)
  write(*,*) match_scalar

  my_regex = re
  write(*,*) my_regex

  my_match = match_scalar
  write(*,*) my_match

  write(*,*) "---------------------------------------------"

  write(*,*) "regex operator (==) example:"
  re = re_posix("[skdj]+","xn")
  match_scalar = re%match(string)

  write(*,*) match_scalar == match(string, "[skdj]+","xn")
  write(*,*) re == re_posix("[skdj]+","xn")
  write(*,*)
  write(*,*) match_scalar == match(string, "[0-9]+","xn")
  write(*,*) re == re_posix("[0-9]+","xn")
  write(*,*)
  write(*,*) match_scalar == match(string, "[skdj]+","xni")
  write(*,*) re == re_posix("[skdj]+","xni")

  write(*,*) "---------------------------------------------"

  write(*,*) "regex match group example:"

  match_scalar = match("asdf@jklp.test", "([a-zA-Z]+)@(([a-zA-Z]+).([a-zA-Z]+))","xn")
  write(*,*) match_scalar
  write(*,*) "group: 1",match_scalar%group(1)
  write(*,*) "group: 2",match_scalar%group(2)
  write(*,*) "group: 3",match_scalar%group(3)
  write(*,*) "group: 4",match_scalar%group(4)
  write(*,*) "group: 5",match_scalar%group(5)
  write(*,*) "group: 6",match_scalar%group(6)
  write(*,*) "group: 7",match_scalar%group(7)

end program

max_stars_repo_name: msvbd/Regex-for-Fortran

max_issues_repo_name: msvbd/Regex-for-Fortran

max_forks_repo_name: msvbd/Regex-for-Fortran


text:
      subroutine sphere_tag4refinement(mx,my,mbc,
     &      meqn, xlower,ylower,dx,dy,t, blockno,
     &      q, tag_threshold, init_flag,tag_patch)
      implicit none

      integer mx,my, mbc, meqn, tag_patch, init_flag
      integer blockno
      double precision xlower, ylower, dx, dy, t
      double precision tag_threshold
      double precision q(1-mbc:mx+mbc,1-mbc:my+mbc,meqn)

      integer example
      common /example_comm/ example

      integer refine_pattern
      common /refine_comm/ refine_pattern

      integer initchoice
      common /initchoice_comm/ initchoice

      integer*8 cont, get_context

      integer i,j, mq
      double precision qmin, qmax, xc, yc, lap
      logical refine

      cont = get_context()

c     # Refine based only on first variable in system.
      refine = .false.
      mq = 1
      do j = 1,my
         do i = 1,mx
c            xc = xlower + (i-0.5)*dx
c            yc = ylower + (j-0.5)*dy
c            lap = (q(i+1,j,1) + q(i-1,j,1) + q(i,j+1,1) + q(i,j-1,1)
c     &            - 4*q(i,j,1))/dx**2
c            if (abs(lap) .gt. tag_threshold) then
c                  refine = .true.
c            endif
            if (init choice .le. 3) then
                refine = q(i,j,mq) .gt.  tag_threshold
            endif
            if (refine) then
                tag_patch = 1
                return
            endif
         enddo
      enddo

      end

max_stars_repo_name: ECLAIRWaveS/ForestClaw

max_issues_repo_name: ECLAIRWaveS/ForestClaw

max_forks_repo_name: ECLAIRWaveS/ForestClaw


text:
      SUBROUTINE DT_MASHEL(Pa1,Pa2,Xm1,Xm2,P1,P2,Irej)

C***********************************************************************
C                                                                      *
C    rescaling of momenta of two partons to put both                   *
C                                       on mass shell                  *
C                                                                      *
C    input:       PA1,PA2   input momentum vectors                     *
C                 XM1,2     desired masses of particles afterwards     *
C                 P1,P2     changed momentum vectors                   *
C                                                                      *
C The original version is written by R. Engel.                         *
C This version dated 12.12.94 is modified by S. Roesler.               *
C***********************************************************************

      IMPLICIT NONE
      DOUBLE PRECISION anorf , bgx , bgy , bgz , cod , cof , del ,
     &                 DT_YLAMB , ee , ee1 , ee2 , gam , ONE , P1 , P2 ,
     &                 Pa1 , Pa2 , pcmp , ppt , ptot1
      DOUBLE PRECISION ptot2 , px , py , pz , sid , sif , ss , TINY10 ,
     &                 Xm1 , xm12 , Xm2 , xm22 , xms , xptot , xx , yy ,
     &                 ZERO , zz
      INTEGER idev , Irej , k , mode
      SAVE

      INCLUDE 'inc/dtflka'

      PARAMETER (TINY10=1.0D-10,ONE=1.0D0,ZERO=0.0D0)

      DIMENSION Pa1(4) , Pa2(4) , P1(4) , P2(4)

      Irej = 0

C Lorentz transformation into system CMS
      px = Pa1(1) + Pa2(1)
      py = Pa1(2) + Pa2(2)
      pz = Pa1(3) + Pa2(3)
      ee = Pa1(4) + Pa2(4)
      xptot = SQRT(px**2+py**2+pz**2)
      xms = (ee-xptot)*(ee+xptot)
C        WRITE(LOUT,'(3E12.4)')XMS,XM1,XM2
      IF ( xms.LT.(Xm1+Xm2)**2 ) THEN

         Irej = 1
         GOTO 99999
      END IF
      xms = SQRT(xms)
      bgx = px/xms
      bgy = py/xms
      bgz = pz/xms
      gam = ee/xms
      CALL DT_DALTRA(gam,-bgx,-bgy,-bgz,Pa1(1),Pa1(2),Pa1(3),Pa1(4),
     &               ptot1,P1(1),P1(2),P1(3),P1(4))
C rotation angles
      cod = P1(3)/ptot1
C     SID = SQRT((ONE-COD)*(ONE+COD))
      ppt = SQRT(P1(1)**2+P1(2)**2)
      sid = ppt/ptot1
      cof = ONE
      sif = ZERO
      IF ( ptot1*sid.GT.TINY10 ) THEN
         cof = P1(1)/(sid*ptot1)
         sif = P1(2)/(sid*ptot1)
         anorf = SQRT(cof*cof+sif*sif)
         cof = cof/anorf
         sif = sif/anorf
      END IF
C new CM momentum and energies (for masses XM1,XM2)
      xm12 = SIGN(Xm1**2,Xm1)
      xm22 = SIGN(Xm2**2,Xm2)
      ss = xms**2
      pcmp = DT_YLAMB(ss,xm12,xm22)/(2.D0*xms)
      ee1 = SQRT(xm12+pcmp**2)
      ee2 = xms - ee1
C back rotation
      mode = 1
      CALL DT_MYTRAN(mode,ZERO,ZERO,pcmp,cod,sid,cof,sif,xx,yy,zz)
      CALL DT_DALTRA(gam,bgx,bgy,bgz,xx,yy,zz,ee1,ptot1,P1(1),P1(2),
     &               P1(3),P1(4))
      CALL DT_DALTRA(gam,bgx,bgy,bgz,-xx,-yy,-zz,ee2,ptot2,P2(1),P2(2),
     &               P2(3),P2(4))
C check consistency
      del = xms*0.0001D0
      IF ( ABS(px-P1(1)-P2(1)).GT.del ) THEN
         idev = 1
      ELSE IF ( ABS(py-P1(2)-P2(2)).GT.del ) THEN
         idev = 2
      ELSE IF ( ABS(pz-P1(3)-P2(3)).GT.del ) THEN
         idev = 3
      ELSE IF ( ABS(ee-P1(4)-P2(4)).GT.del ) THEN
         idev = 4
      ELSE
         idev = 0
      END IF
      IF ( idev.NE.0 ) THEN

         IF ( LPRi.GT.4 ) WRITE (LOUt,'(/1X,A,I3)')
     &         'MASHEL: inconsistent transformation' , idev

         IF ( LPRi.GT.4 ) WRITE (LOUt,'(1X,A)')
     &         'MASHEL: input momenta/masses:'

         IF ( LPRi.GT.4 ) WRITE (LOUt,'(1X,5E12.5)') (Pa1(k),k=1,4) ,
     &        Xm1

         IF ( LPRi.GT.4 ) WRITE (LOUt,'(1X,5E12.5)') (Pa2(k),k=1,4) ,
     &        Xm2

         IF ( LPRi.GT.4 ) WRITE (LOUt,'(1X,A)')
     &         'MASHEL: output momenta:'

         IF ( LPRi.GT.4 ) WRITE (LOUt,'(5X,4E12.5)') (P1(k),k=1,4)

         IF ( LPRi.GT.4 ) WRITE (LOUt,'(5X,4E12.5)') (P2(k),k=1,4)
      END IF
99999 END SUBROUTINE

max_stars_repo_name: pzhristov/DPMJET

max_issues_repo_name: pzhristov/DPMJET

max_forks_repo_name: pzhristov/DPMJET


text:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!                                                                   !!
!!                   GNU General Public License                      !!
!!                                                                   !!
!! This file is part of the Flexible Modeling System (FMS).          !!
!!                                                                   !!
!! FMS is free software; you can redistribute it and/or modify       !!
!! it and are expected to follow the terms of the GNU General Public !!
!! License as published by the Free Software Foundation.             !!
!!                                                                   !!
!! FMS 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 FMS; if not, write to:                                 !!
!!          Free Software Foundation, Inc.                           !!
!!          59 Temple Place, Suite 330                               !!
!!          Boston, MA  02111-1307  USA                              !!
!! or see:                                                           !!
!!          http://www.gnu.org/licenses/gpl.txt                      !!
!!                                                                   !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

module bgrid_vert_mod

!-----------------------------------------------------------------------
!
!     allocates storage and initializes vertical grid constants
!
!     contains several interfaces for compute pressure and height
!
!-----------------------------------------------------------------------

use types_mod, only : r8
use constants_mod, only: GRAV, RDGAS, RVGAS
use       fms_mod, only: error_mesg, FATAL

implicit none
private

!-----------------------------------------------------------------------
!  public defined data type (vert_grid_type)
!  -----------------------------------------
!    nlev  = number of vertical levels (integer)
!    nplev = number of pure pressure levels at the top of the model,
!            will equal zero when not using hybrid coordinate (integer)
!    deta  = vertical eta thickness/depth of model layers
!    aeta  = eta values at model levels (full levels)
!    eta   = eta values at model layer interfaces (half levels)
!    dpeta = vertical pressure thickness/depth of model layers
!    apeta = pressure values at model levels (full levels)
!    peta  = pressure values at model layer interfaces (half levels)
!    dfl   = reference values of geopotental height at half levels
!
!    psmin  = minimum allowable surface pressure to avoid negative
!             mass in a model layer (important for hybrid coord)
!    hybrid = logical flag (true for hybrid coordinate)
!    pzero  = logical flag (true for pres = 0 at top of model)
!
!    pref,tref,gamma = reference values for computing dfl
!
!
!  public interfaces
!  -----------------
!    vert_grid_init      - initializes the vert_grid_type
!    compute_height_half - computes the geopotential height (in meters)
!                            at half model levels
!    compute_pres_depth  - computes the pressure weight (mass) of a
!                            model layer
!    compute_pres_full   - computes the pressure at full model levels
!    compute_pres_half   - computes the pressure at half model levels
!
!-----------------------------------------------------------------------
!------- interfaces -------

public  vert_grid_init, compute_geop_height, compute_height,   &
        compute_pres_depth, compute_pres_full, compute_pres_half, &
        compute_pres_weights, compute_pressures, compute_height_bottom

!------- public defined data type -------

public  vert_grid_type

type vert_grid_type
   integer                     :: nlev, nplev
   real(r8), pointer, dimension(:) :: deta,  aeta,  eta,  dfl,  &
                                  dpeta, apeta, peta, wta, wtb
   real(r8)                        :: pref, tref, gamma, psmin
   logical                     :: hybrid, pzero
end type vert_grid_type

!-----------------------------------------------------------------------

  real(r8), parameter :: d608 = (RVGAS-RDGAS)/RDGAS
  real(r8), parameter :: ginv = 1./GRAV

!------ parameters for eta coordinate reference surface heights --------

  real(r8), parameter :: pref = 101325., tref = 288., gamma = 0.0065

!------ parameters for performance timing of code sections -----

  integer, parameter :: timlev = 9
  logical :: do_clock_init = .true.
  integer, dimension(7) :: id
  character(len=16), dimension(7) :: &
            names = (/ 'comp_pres_depth ', &
                       'comp_pres_full  ', &
                       'comp_pressures  ', &
                       'comp_pres_half  ', &
                       'comp_pres_wghts ', &
                       'comp_geop_hght  ', &
                       'comp_height_btm ' /)

!-----------------------------------------------------------------------

contains

!#######################################################################

function vert_grid_init (eta, peta, verbose) result (Vgrid)

!-----------------------------------------------------------------------
!
!    all arguments are input only
!    ----------------------------
!
!    deta_in  =  vertical grid spacing, an array of order one
!

     real(r8), intent (in)           ::  eta(:)
     real(r8), intent (in), optional :: peta(:)
  integer, intent (in), optional :: verbose

     type(vert_grid_type) :: Vgrid

!-----------------------------------------------------------------------
integer :: k, kx, lverbose
real(r8)    :: rgog
real(r8), dimension(size(eta))   :: phalf, lphalf, pres
real(r8), dimension(size(eta)-1) ::        lpfull

!-----------------------------------------------------------------------

   lverbose = 0;  if (present(verbose)) lverbose = verbose

!--------------derived vertical constants-------------------------------

      kx = size(eta) - 1
      allocate (Vgrid% deta(kx), Vgrid% aeta(kx), Vgrid% eta(kx+1),  &
                Vgrid%dpeta(kx), Vgrid%apeta(kx), Vgrid%peta(kx+1),  &
                Vgrid%wta  (kx), Vgrid%wtb  (kx),                    &
                Vgrid% dfl(kx+1))

      Vgrid % nlev = kx

!--------- set-up eta values and hybrid pressure levels ----------
!--------- note: eta(1) and eta(kx+1) have set values -----
!--------- also note: peta(kx+1) = 0.0 -----

        Vgrid % eta(1)    = 0.0
        Vgrid % eta(kx+1) = 1.0
        Vgrid % eta(2:kx) = eta(2:kx)

        Vgrid % peta = 0.0
        if (present(peta)) Vgrid % peta(1:kx) = peta(1:kx)

      do k = 1, kx
        Vgrid %  deta(k) =  Vgrid %  eta(k+1) - Vgrid %  eta(k)
        Vgrid %  aeta(k) = (Vgrid %  eta(k+1) + Vgrid %  eta(k)) * 0.5
        Vgrid % dpeta(k) =  Vgrid % peta(k+1) - Vgrid % peta(k)
        Vgrid % apeta(k) =  0.0
        Vgrid % wta  (k) =  0.0
        Vgrid % wtb  (k) =  0.0
      enddo

!----------- is this a hybrid coordinate ??? -----

      Vgrid % hybrid = .false.

      do k = 1, kx+1
         if ( Vgrid % peta(k) > 0.0 ) then
              Vgrid % hybrid = .true.
              exit
         endif
      enddo

!----------- find lowest pure pressure level --------------

      Vgrid % nplev = 0

      do k = 1, kx
         if ( Vgrid % deta(k) > 0.0 ) exit
         Vgrid % nplev = k
      enddo

!     ---- need average pressure in these layers ----

      Vgrid % pzero = .true.

      if ( Vgrid % nplev >= 1 ) then
         phalf(:) = Vgrid%peta(:) + Vgrid%eta(:)*pref
         if ( phalf(1) <= epsilon(phalf) ) then
            lphalf(1) = 0.0
            lphalf(2:) = log(phalf(2:))
         else
            lphalf(:) = log(phalf(:))
            Vgrid % pzero = .false.
         endif

         do k = 1, kx
            lpfull(k) = (phalf(k+1)*lphalf(k+1) - phalf(k)*lphalf(k)) &
                        / (phalf(k+1)-phalf(k)) - 1.0
            Vgrid % apeta(k) = exp(lpfull(k))
            Vgrid % wtb  (k) = lphalf(k+1) - lpfull(k)
            Vgrid % wta  (k) = lpfull(k)   - lphalf(k)
         enddo
         if (Vgrid % pzero) Vgrid % wta(1) =Vgrid % wtb(1)

      endif

!----------- find the minimum allowable surface pressure ------

      Vgrid % psmin = 0.0

      do k = 1, kx
         if ( Vgrid % deta(k) > 0.0 )  Vgrid % psmin =  &
               max ( Vgrid % psmin, -Vgrid % dpeta(k)/Vgrid % deta(k) )
      enddo

!--- optional output of coordinate values ----

   if (lverbose > 1) then
       print *,  ' eta=',Vgrid %   eta
       print *,  'deta=',Vgrid %  deta
       print *,  'aeta=',Vgrid %  aeta
       print *, ' peta=',Vgrid %  peta
       print *, 'dpeta=',Vgrid % dpeta
       print *, 'apeta=',Vgrid % apeta
       print *, 'nplev=',Vgrid % nplev
       print *, 'minimum allowable surface pressure = ', Vgrid % psmin
   endif

!---------- set-up eta coordinate geopotential heights -------------

      rgog = RDGAS*gamma/GRAV

      do k=1,kx
        pres(k) = Vgrid%peta(k) + Vgrid%eta(k)*pref
        Vgrid % dfl(k) = GRAV*tref*(1.0-(pres(k)/pref)**rgog)/gamma
      enddo
        Vgrid % dfl(kx+1) = 0.0

        Vgrid % pref  = pref
        Vgrid % tref  = tref
        Vgrid % gamma = gamma

! initialize code sections for performance timing

 !if (do_clock_init) then
 !  do k = 1, size(id)
 !     id(k) = mpp_clock_init ('BGRID: vert ('//trim(names(k))//')', timlev, flags=MPP_CLOCK_SYNC)
 !  enddo
 !  do_clock_init = .false.
 !endif

!-----------------------------------------------------------------------

end function vert_grid_init

!#######################################################################

subroutine compute_pres_depth (Vgrid, pssl, pdepth)

  type(vert_grid_type), intent(in)  :: Vgrid
       real(r8)           , intent(in)  :: pssl(:,:)
       real(r8)           , intent(out) :: pdepth(:,:,:)

       integer :: k, kp, ke

!-----------------------------------------------------------------------
!   compute the pressure weight (depth) for model layers
!-----------------------------------------------------------------------
  !call mpp_clock_begin (id(1))

  kp = Vgrid % nplev
  ke = Vgrid % nlev

  if (size(pdepth,3) /= ke) call error_mesg (  &
                              'compute_pres_depth in bgrid_vert_mod', &
                              'incorrect dimension 3 for pdepth', FATAL)

! --- check for zero/negative depth layers ---

  if (Vgrid % hybrid) then
      if (minval(pssl) <= Vgrid % psmin) call error_mesg  &
                 ('compute_pres_depth in bgrid_vert_mod',  &
                  'pressure depth <= 0.0', FATAL)
  endif

! --- compute depth ---

  do k = 1, kp
    pdepth(:,:,k) = Vgrid % dpeta(k)
  enddo

  do k = kp+1, ke
    pdepth(:,:,k) = Vgrid % dpeta(k) + Vgrid % deta(k) * pssl(:,:)
  enddo

  !call mpp_clock_end (id(1))

end subroutine compute_pres_depth

!#######################################################################

subroutine compute_pres_full (Vgrid, pssl, pfull, phalf, dpde)

  type(vert_grid_type), intent(in)  :: Vgrid
  real(r8)                , intent(in)  :: pssl(:,:)
  real(r8)                , intent(out) :: pfull(:,:,:)
  real(r8), optional      , intent(in)  :: phalf(:,:,:), dpde(:,:,:)

   real(r8), dimension(size(pfull,1),size(pfull,2),size(pfull,3)+1) :: ph
   real(r8), dimension(size(pfull,1),size(pfull,2),size(pfull,3))   :: dp
   integer :: k, kp, ke
!-----------------------------------------------------------------------
!      compute the pressure at full model levels
!-----------------------------------------------------------------------
  !call mpp_clock_begin (id(2))

  kp = Vgrid % nplev
  ke = Vgrid % nlev

  if (size(pfull,3) /= ke) call error_mesg (  &
                               'compute_pres_full in bgrid_vert_mod', &
                               'incorrect dimension 3 for pfull', FATAL)

!--- set or compute optional arguments ---

  if (present(phalf)) then
      ph = phalf
  else
      call compute_pres_half (Vgrid, pssl, ph)
  endif

  if (present(dpde)) then
      dp = dpde
  else
      call compute_pres_depth (Vgrid, pssl, dp)
  endif

!--- compute p*logp at half levels ---

  if ( Vgrid % pzero ) then
       ph(:,:,1) = 0.0
       ph(:,:,2:ke+1) = ph(:,:,2:ke+1) * log(ph(:,:,2:ke+1))
  else
       ph(:,:,:) = ph(:,:,:) * log(ph(:,:,:))
  endif

!--- compute pressure at full levels ---

  do k = 1, kp
    pfull(:,:,k) = Vgrid % apeta(k)
  enddo

  do k = kp+1, ke
    pfull(:,:,k) = exp( (ph(:,:,k+1)-ph(:,:,k))/dp(:,:,k) - 1.0 )
  enddo

  !call mpp_clock_end (id(2))

end subroutine compute_pres_full

!#######################################################################

subroutine compute_pressures (Vgrid, pssl, phalf, pfull, dpde, wta, wtb)

  type(vert_grid_type), intent(in)  :: Vgrid
  real(r8)                , intent(in)  :: pssl(:,:)
  real(r8)                , intent(out) :: phalf(:,:,:), pfull(:,:,:)
  real(r8), optional      , intent(out) :: dpde(:,:,:)
  real(r8), optional      , intent(out) :: wta(:,:,:), wtb(:,:,:)

  real(r8), dimension(size(pfull,1),size(pfull,2),size(pfull,3)+1) :: ph,lph
  real(r8), dimension(size(pfull,1),size(pfull,2),size(pfull,3))   :: dp,lpf
  integer :: k, kp, ke
!-----------------------------------------------------------------------
!      compute the pressure at full model levels
!-----------------------------------------------------------------------

  kp = Vgrid % nplev
  ke = Vgrid % nlev

  if (size(pfull,3) /= ke) call error_mesg (  &
                               'compute_pressures in bgrid_vert_mod', &
                               'incorrect dimension 3 for pfull', FATAL)

!--- set or compute optional arguments ---

  call compute_pres_half (Vgrid, pssl, ph)
  phalf = ph

  call compute_pres_depth (Vgrid, pssl, dp)
  if (present(dpde))  dpde = dp

! do not include time for previous calls
  !call mpp_clock_begin (id(3))

!--- compute p*logp at half levels ---

  if ( Vgrid % pzero ) then
      lph(:,:,1) = 0.0
       ph(:,:,1) = 0.0
      lph(:,:,2:ke+1) = log(ph(:,:,2:ke+1))
       ph(:,:,2:ke+1) = ph(:,:,2:ke+1) * lph(:,:,2:ke+1)
  else
      lph(:,:,:) = log(ph(:,:,:))
       ph(:,:,:) = ph(:,:,:) * lph(:,:,:)
  endif

!--- compute pressure at full levels ---

  do k = 1, kp
    pfull(:,:,k) = Vgrid % apeta(k)
  enddo

! if (present(wta) .or. present(wtb)) then
!   do k = 1, kp
!     lpf(:,:,k) = log(pfull(:,:,k))
!   enddo
! endif

  do k = kp+1, ke
    lpf(:,:,k) = (ph(:,:,k+1)-ph(:,:,k))/dp(:,:,k) - 1.0
    pfull(:,:,k) = exp( lpf(:,:,k) )
  enddo

!--- compute weights at full levels ---

   if (present(wtb)) then
      do k = 1, kp
        wtb(:,:,k) = Vgrid % wtb(k)
      enddo
      do k = kp+1, size(wtb,3)
        wtb(:,:,k) = lph(:,:,k+1) - lpf(:,:,k)
      enddo
   endif

   if (present(wta)) then
      do k = 1, kp
        wta(:,:,k) = Vgrid % wta(k)
      enddo
      do k = kp+1, size(wta,3)
        wta(:,:,k) = lpf(:,:,k) - lph(:,:,k)
      enddo
      if (Vgrid % pzero .and. kp == 0) wta(:,:,1) = wtb(:,:,1)
   endif

  !call mpp_clock_end (id(3))

end subroutine compute_pressures

!#######################################################################

subroutine compute_pres_half (Vgrid, pssl, phalf)

  type(vert_grid_type), intent(in)  :: Vgrid
       real(r8)           , intent(in)  :: pssl(:,:)
       real(r8)           , intent(out) :: phalf(:,:,:)

       integer :: k, kp, ke

!-----------------------------------------------------------------------
!    compute the pressure at the interface between model layers
!                      (half model levels)
!-----------------------------------------------------------------------
  !call mpp_clock_begin (id(4))

  kp = Vgrid % nplev + 1
  ke = Vgrid % nlev + 1

  if (size(phalf,3) /= ke) call error_mesg (  &
                               'compute_pres_half in bgrid_vert_mod', &
                               'incorrect dimension 3 for phalf', FATAL)

  do k = 1, kp
    phalf(:,:,k) = Vgrid % peta(k)
  enddo

  do k = kp+1, ke
    phalf(:,:,k) = Vgrid % peta(k) + Vgrid % eta(k) * pssl(:,:)
  enddo

  !call mpp_clock_end (id(4))

end subroutine compute_pres_half

!#######################################################################

subroutine compute_pres_weights ( Vgrid, phalf, pfull, wta, wtb )

 type(vert_grid_type), intent(in)  :: Vgrid
 real(r8), intent(in)  :: phalf(:,:,:), pfull(:,:,:)
 real(r8), intent(out) :: wta(:,:,:), wtb(:,:,:)

 real(r8), dimension(size(phalf,1),size(phalf,2),size(phalf,3)) :: logph
 real(r8), dimension(size(pfull,1),size(pfull,2),size(pfull,3)) :: logpf
 integer :: k, kp, kx, ks

  !call mpp_clock_begin (id(5))

   kp = Vgrid % nplev
   kx = size(pfull,3)

   if (Vgrid%pzero) then
       ks = max(2,kp+1)
   else
       ks = kp+1
   endif

   logph(:,:,ks  :kx+1) = log(phalf(:,:,ks  :kx+1))
   logpf(:,:,kp+1:kx  ) = log(pfull(:,:,kp+1:kx  ))

   do k = 1, kp
     wtb(:,:,k) = Vgrid % wtb(k)
   enddo
   do k = kp+1, kx
     wtb(:,:,k) = logph(:,:,k+1) - logpf(:,:,k)
   enddo

   if (Vgrid%pzero .and. kp == 0) wta(:,:,1) = wtb(:,:,1)
   do k = 1, kp
     wta(:,:,k) = Vgrid % wta(k)
   enddo
   do k = ks, kx
     wta(:,:,k) = logpf(:,:,k) - logph(:,:,k)
   enddo

  !call mpp_clock_end (id(5))

end subroutine compute_pres_weights

!#######################################################################

subroutine compute_geop_height (Vgrid, fssl, vtemp, wta, wtb, &
                                zfull, zhalf, mask)

    type(vert_grid_type), intent(in)   :: Vgrid
    real(r8), intent(in), dimension(:,:)   :: fssl
    real(r8), intent(in), dimension(:,:,:) :: vtemp, wta, wtb
    real(r8), intent(out)                  :: zfull(:,:,:)
    real(r8), intent(out), optional        :: zhalf(:,:,:)
    real(r8), intent(in),  optional        :: mask(:,:,:)

       integer :: k, klev
       real(r8), dimension(size(vtemp,1),size(vtemp,2)) :: zb, zt, rt

!-----------------------------------------------------------------------
!    compute the height (in meters) at the interface between
!              model layers (half model levels)
!-----------------------------------------------------------------------
  !call mpp_clock_begin (id(6))

  klev = Vgrid % nlev

  if (size(zfull,3) /= klev) call error_mesg (  &
                               'compute_geop_height in bgrid_vert_mod', &
                               'incorrect dimension 3 for zfull', FATAL)

  if (present(zhalf)) then
     if (size(zhalf,3) /= klev+1) call error_mesg (  &
                               'compute_geop_height in bgrid_vert_mod', &
                               'incorrect dimension 3 for zhalf', FATAL)
  endif

      zb(:,:) = fssl(:,:)
      if (present(zhalf)) zhalf(:,:,klev+1) = zb(:,:)

!------- vertical integration loop (bottom to top) ----------

      do k = klev, 1, -1

         rt(:,:) = RDGAS * vtemp(:,:,k)
         zt(:,:) = zb(:,:) + rt(:,:) * (wta(:,:,k)+wtb(:,:,k))
         zfull(:,:,k) = zb(:,:) + rt(:,:) * wtb(:,:,k)
         if (present(mask)) then
            where (mask(:,:,k) < 0.5)
                zt(:,:) = Vgrid % dfl(k)
                zfull(:,:,k) = 0.5*(Vgrid % dfl(k)+Vgrid % dfl(k+1))
            endwhere
         endif
         zb(:,:) = zt(:,:)
         if (present(zhalf)) zhalf(:,:,k) = zb(:,:)

      enddo

  !call mpp_clock_end (id(6))

end subroutine compute_geop_height

!#######################################################################

subroutine compute_height (Vgrid, fssl, temp, sphum, pfull, phalf, &
                           zfull, zhalf, mask)

   type(vert_grid_type), intent(in)    :: Vgrid
   real(r8), intent(in),  dimension(:,:)   :: fssl
   real(r8), intent(in),  dimension(:,:,:) :: temp, sphum, pfull, phalf
   real(r8), intent(out), dimension(:,:,:) :: zfull, zhalf
   real(r8), intent(in),  optional         :: mask(:,:,:)

   real(r8), dimension(size(temp,1),size(temp,2),size(temp,3)) ::  &
                                                   wta, wtb, vtemp

!-----------------------------------------------------------------------
!    compute the height (in meters) at the interface between
!              model layers (half model levels)
!    assumes that specific humidity will be used to compute
!                the virtual temperature
!-----------------------------------------------------------------------

  call compute_pres_weights ( Vgrid, phalf, pfull, wta, wtb )

  vtemp = temp * (1.0+d608*sphum)

  if (present(mask)) then
     call compute_geop_height ( Vgrid, fssl, vtemp, wta, wtb, &
                                zfull, zhalf, mask )
  else
     call compute_geop_height ( Vgrid, fssl, vtemp, wta, wtb, &
                                zfull, zhalf )
  endif

  zfull = zfull * ginv
  zhalf = zhalf * ginv

end subroutine compute_height

!#######################################################################

subroutine compute_height_bottom ( Vgrid, pssl, tbot, qbot,  &
                                   zbot, pbot, kbot )

   type(vert_grid_type), intent(in)  :: Vgrid
   real(r8), intent(in),  dimension(:,:) :: pssl, tbot, qbot
   real(r8), intent(out), dimension(:,:) :: zbot, pbot
   integer, intent(in),  optional    :: kbot(:,:)

   real(r8), dimension(size(pssl,1),size(pssl,2)) :: rt, dp, phb, pht, &
                                                 lphb, lpht, lpf
   integer :: i, j, kb

!  ----- pressure at top and bottom interface of bottom level -----

   !call mpp_clock_begin (id(7))

   if (present(kbot)) then
      do j = 1, size(pssl,2)
      do i = 1, size(pssl,1)
         kb = kbot(i,j)
         pht(i,j) = Vgrid%peta(kb  ) + Vgrid%eta(kb  )*pssl(i,j)
         phb(i,j) = Vgrid%peta(kb+1) + Vgrid%eta(kb+1)*pssl(i,j)
      enddo
      enddo
   else
         kb = Vgrid%nlev
         pht(:,:) = Vgrid%peta(kb+1) + Vgrid%eta(kb+1)*pssl(:,:)
         phb(:,:) = Vgrid%peta(kb  ) + Vgrid%eta(kb  )*pssl(:,:)
   endif

   rt = ginv*RDGAS * (tbot * (1.+d608*qbot))
   dp = phb - pht
   lphb = log(phb)
   lpht = log(pht)
   lpf = (phb*lphb-pht*lpht)/dp -1

   zbot = rt * (lphb-lpf)
   pbot = exp(lpf)

   !call mpp_clock_end (id(7))

end subroutine compute_height_bottom

!#######################################################################

end module bgrid_vert_mod

max_stars_repo_name: hkershaw-brown/feature-preprocess

max_issues_repo_name: hkershaw-brown/feature-preprocess

max_forks_repo_name: hkershaw-brown/feature-preprocess


text:
*----------------------------------------------------------------------*
      subroutine optc_orthvec(nadd,prenorm,
     &     ff_ssbsp,iord_ssbsp,sred,
     &     ff_sbsp,iord_sbsp,ndim_sbsp,mxsbsp,zero_vec,
     &     use_s,ioff_snew,ffsnew,ffnew,nnew,nopt,
     &     nwfpar,nincore,xbuf1,xbuf2,xbuf3,lenbuf)
*----------------------------------------------------------------------*
*     orthogonalize and add nnew new vectors (on ff_new) to subspace
*     (on ff_sbsp); actual number of new vectors returned on nadd
*----------------------------------------------------------------------*
      implicit none

      include 'stdunit.h'
      include 'def_filinf.h'
      include 'def_file_array.h'

      integer, parameter ::
     &     ntest = 00

      integer, intent(out) ::
     &     nadd
      integer, intent(in) ::
     &     nnew, mxsbsp, lenbuf, nincore, nopt, nwfpar(nopt), ioff_snew
      real(8), intent(in) ::
     &     sred(mxsbsp,ndim_sbsp)
      integer, intent(inout) ::
     &     ndim_sbsp, iord_sbsp(mxsbsp), iord_ssbsp(mxsbsp)
      type(file_array), intent(in) ::
     &     ff_sbsp(nopt), ff_ssbsp(nopt)
      type(file_array), intent(in) ::
     &     ffnew(nopt), ffsnew(nopt)
      real(8), intent(inout) ::
     &     xbuf1(*), xbuf2(*), xbuf3(*)
      logical, intent(in) ::
     &     zero_vec(ndim_sbsp), prenorm, use_s(*)

      integer ::
     &     nold, inew, jnew, iold, ii, jj, idx, ivec, irec, iopt,
     &     jold
      real(8) ::
     &     xdum, xnorm, xnorm_
      logical ::
     &     previous_zero

      integer ::
     &     iordnew(nnew)
      real(8), allocatable ::
     &     smat(:,:), xmat(:,:), scr(:)

      integer, external ::
     &     ioptc_get_sbsp_rec
      real(8), external ::
     &     ddot, da_ddot

      if (ntest.ge.100) then
        call write_title(lulog,wst_dbg_subr,'optc_orthvec')
        write(lulog,*) 'ndim, nnew, nopt: ',ndim_sbsp, nnew, nopt
        write(lulog,*) 'use_s: ',use_s(1:nopt)
      end if

      ! look for zero-vectors (in initial iteration); they should
      ! reside at the end of the initial subspace
      nold = 0
      previous_zero = .false.
      do ivec = 1, ndim_sbsp
        if (zero_vec(ivec)) then
          previous_zero = .true.
        else
          if (previous_zero) then
            write(lulog,*) 'zero_vec: ',zero_vec(1:ndim_sbsp)
            call quit(1,'optc_orthvec',
     &           'zero vectors should reside at the end of the'//
     &           ' initial subspace')
          end if
          nold = nold+1
        end if
      end do

      ! calculate overlap matrix
      allocate(
     &     smat(nold+nnew,nold+nnew),
     &     xmat(nold+nnew,nold+nnew),
     &     scr(2*(nold+nnew)))

      if (ntest.ge.100)
     &     write(lulog,*) 'nold, nnew: ',nold, nnew

      smat(1:nold+nnew,1:nold+nnew) = 0d0
      do iold = 1, nold
        smat(iold,iold) = 1d0
      end do
      if (any(use_s(1:nopt))) smat(1:nold,1:nold) = sred(1:nold,1:nold)

      do iopt = 1, nopt

c      ! 0) <old|old> (lower triangle in case of asymmetric metric)
c      if (use_s(iopt)) then
c        if (nincore.lt.2) call quit(1,'optc_orthvec',
c     &                              'not yet : nincore<2')
c        do jold = 1, nold
c
c          jj = iord_sbsp(jold)
c            call vec_from_da(ff_sbsp(iopt)%fhand,jold,
c     &                                          xbuf1,nwfpar(iopt))
c
c          do iold = 1, nold
c
c              ii = iord_ssbsp(iold)
c              if (jj.le.ii) cycle ! only lower triangle
c              call vec_from_da(ff_ssbsp(iopt)%fhand,iold,
c     &                                             xbuf2,nwfpar(iopt))
c              smat(jj,ii) = smat(jj,ii)
c     &             + ddot(nwfpar(iopt),xbuf1,1,xbuf2,1)
c
c          end do
c
c        end do
c
c      end if

      ! a) <old|new>
      if (nold.gt.0.and.nold.ge.nnew) then
        do inew = 1, nnew

          ii = inew+nold
          if (nincore.ge.2.and.use_s(iopt)) then
c dbg
c           print *,'1 loading metric record',inew+ioff_snew
c dbgend
            call vec_from_da(ffsnew(iopt)%fhand,inew+ioff_snew,xbuf1,
     &                       nwfpar(iopt))
          else if (nincore.ge.2) then
            call vec_from_da(ffnew(iopt)%fhand,inew,xbuf1,nwfpar(iopt))
          end if

          do iold = 1, nold

            jj = iord_sbsp(iold)
            if (nincore.ge.2.and.use_s(iopt)) then
              call vec_from_da(ff_sbsp(iopt)%fhand,iold,
     &                                             xbuf2,nwfpar(iopt))
              smat(jj,ii) = smat(jj,ii)
     &             + ddot(nwfpar(iopt),xbuf1,1,xbuf2,1)
c              smat(ii,jj) = smat(jj,ii)
            else if (nincore.ge.2) then
              call vec_from_da(ff_sbsp(iopt)%fhand,iold,
     &                                             xbuf2,nwfpar(iopt))
              smat(jj,ii) = smat(jj,ii)
     &                    + ddot(nwfpar(iopt),xbuf1,1,xbuf2,1)
              smat(ii,jj) = smat(jj,ii)
            else if (use_s(iopt)) then
              smat(jj,ii) = smat(jj,ii)
     &            + da_ddot(ff_sbsp(iopt)%fhand,iold,
     &                      ffsnew(iopt)%fhand,inew+ioff_snew,
     &                      nwfpar(iopt),xbuf1,xbuf2,lenbuf)
c              smat(ii,jj) = smat(jj,ii)
            else
              smat(jj,ii) = smat(jj,ii) +
     &            da_ddot(ff_sbsp(iopt)%fhand,iold,
     &                    ffnew(iopt)%fhand,  inew,
     &                    nwfpar(iopt),xbuf1,xbuf2,lenbuf)
              smat(ii,jj) = smat(jj,ii)
            end if

          end do

        end do

      else if (nold.gt.0.and.nnew.gt.nold) then

        do iold = 1, nold

          ii = iord_sbsp(iold)
          if (nincore.ge.2) then
            call vec_from_da(ff_sbsp(iopt)%fhand,iold,
     &                                            xbuf1,nwfpar(iopt))
          end if

          do inew = 1, nnew

            jj = inew+nold
            if (nincore.ge.2.and.use_s(iopt)) then
c dbg
c               print *,'2 loading metric record',inew+ioff_snew
c dbgend
              call vec_from_da(ffsnew(iopt)%fhand,inew+ioff_snew,xbuf2,
     &                         nwfpar(iopt))
              smat(jj,ii) = smat(jj,ii)
     &                + ddot(nwfpar(iopt),xbuf1,1,xbuf2,1)
c              smat(ii,jj) = smat(jj,ii)
            else if (nincore.ge.2) then
              call vec_from_da(ffnew(iopt)%fhand,inew,
     &                                          xbuf2,nwfpar(iopt))
              smat(jj,ii) = smat(jj,ii) +
     &             ddot(nwfpar(iopt),xbuf1,1,xbuf2,1)
              smat(ii,jj) = smat(jj,ii)
            else if (use_s(iopt)) then
              smat(jj,ii) = smat(jj,ii)
     &           + da_ddot(ff_sbsp(iopt)%fhand,iold,
     &                     ffsnew(iopt)%fhand,inew+ioff_snew,
     &                     nwfpar(iopt),xbuf1,xbuf2,lenbuf)
c              smat(ii,jj) = smat(jj,ii)
            else
              smat(jj,ii) = smat(jj,ii) +
     &            da_ddot(ff_sbsp(iopt)%fhand,iold,
     &                    ffnew(iopt)%fhand,  inew,
     &                    nwfpar(iopt),xbuf1,xbuf2,lenbuf)
              smat(ii,jj) = smat(jj,ii)
            end if

          end do

        end do

      end if

      ! a) <new|old> (needed seperately for asymmetric metric)
      if (use_s(iopt).and.nold.gt.0.and.nold.ge.nnew) then
        if (nincore.lt.2) call quit(1,'optc_orthvec',
     &                              'not yet : nincore<2')

        do inew = 1, nnew
          jj = inew+nold
            call vec_from_da(ffnew(iopt)%fhand,inew,xbuf1,nwfpar(iopt))
          do iold = 1, nold

              ii = iord_ssbsp(iold)
              call vec_from_da(ff_ssbsp(iopt)%fhand,iold,
     &                                             xbuf2,nwfpar(iopt))
              smat(jj,ii) = smat(jj,ii)
     &             + ddot(nwfpar(iopt),xbuf1,1,xbuf2,1)
          end do
        end do

      else if (use_s(iopt).and.nold.gt.0.and.nnew.gt.nold) then
        if (nincore.lt.2) call quit(1,'optc_orthvec',
     &                              'not yet : nincore<2')

        do iold = 1, nold
            jj = iord_ssbsp(iold)
            call vec_from_da(ff_ssbsp(iopt)%fhand,iold,
     &                                            xbuf1,nwfpar(iopt))
          do inew = 1, nnew
            ii = inew+nold
              call vec_from_da(ffnew(iopt)%fhand,inew,xbuf2,
     &                         nwfpar(iopt))
              smat(jj,ii) = smat(jj,ii)
     &                + ddot(nwfpar(iopt),xbuf1,1,xbuf2,1)
          end do
        end do

      end if

      ! b) <new|new>
      do inew = 1, nnew
        if (nincore.ge.2) then
          call vec_from_da(ffnew(iopt)%fhand,inew,xbuf1,
     &                     nwfpar(iopt))
          if (use_s(iopt)) then
c dbg
c           print *,'3 loading metric record',inew+ioff_snew
c dbgend
            call vec_from_da(ffsnew(iopt)%fhand,inew+ioff_snew,xbuf2,
     &                       nwfpar(iopt))
            smat(nold+inew,nold+inew) = smat(nold+inew,nold+inew) +
     &           ddot(nwfpar(iopt),xbuf1,1,xbuf2,1)
          else
            smat(nold+inew,nold+inew) = smat(nold+inew,nold+inew) +
     &           ddot(nwfpar(iopt),xbuf1,1,xbuf1,1)
          end if
c dbg
c           print *,'|new|: ',sqrt(smat(nold+inew,nold+inew))
c dbg
        else if (use_s(iopt)) then
          smat(nold+inew,nold+inew) = smat(nold+inew,nold+inew) +
     &         da_ddot(ffnew(iopt)%fhand,inew,
     &                 ffsnew(iopt)%fhand,inew+ioff_snew,
     &                 nwfpar(iopt),xbuf1,xbuf2,lenbuf)
        else
          smat(nold+inew,nold+inew) = smat(nold+inew,nold+inew) +
     &         da_ddot(ffnew(iopt)%fhand,inew,
     &                 ffnew(iopt)%fhand,inew,
     &         nwfpar(iopt),xbuf1,xbuf2,lenbuf)
        end if

c        do jnew = inew+1, nnew
        do jnew = 1, nnew
          if (inew.eq.jnew) cycle
          if (nincore.ge.2.and.use_s(iopt)) then
c dbg
c             print *,'4 loading metric record',jnew+ioff_snew
c dbgend
            call vec_from_da(ffsnew(iopt)%fhand,jnew+ioff_snew,xbuf2,
     &                         nwfpar(iopt))
            smat(nold+inew,nold+jnew) = smat(nold+inew,nold+jnew) +
     &           ddot(nwfpar(iopt),xbuf1,1,xbuf2,1)
c            smat(nold+jnew,nold+inew) = smat(nold+inew,nold+jnew)
          else if (nincore.ge.2) then
            if (jnew.le.inew) cycle
            call vec_from_da(ffnew(iopt)%fhand,jnew,xbuf2,nwfpar(iopt))
            smat(nold+inew,nold+jnew) = smat(nold+inew,nold+jnew) +
     &           ddot(nwfpar(iopt),xbuf1,1,xbuf2,1)
            smat(nold+jnew,nold+inew) = smat(nold+inew,nold+jnew)
          else
            if (jnew.le.inew) cycle ! revise this for asymmetric metric
            if (use_s(iopt)) then
              smat(nold+inew,nold+jnew) = smat(nold+inew,nold+jnew) +
     &           da_ddot(ffnew(iopt)%fhand,inew,
     &                   ffsnew(iopt)%fhand,jnew+ioff_snew,
     &                   nwfpar(iopt),xbuf1,xbuf2,lenbuf)
            else
              smat(nold+inew,nold+jnew) = smat(nold+inew,nold+jnew) +
     &             da_ddot(ffnew(iopt)%fhand,inew,
     &                     ffnew(iopt)%fhand,jnew,
     &             nwfpar(iopt),xbuf1,xbuf2,lenbuf)
            end if
            smat(nold+jnew,nold+inew) = smat(nold+inew,nold+jnew)
          end if
        end do

      end do

      if (ntest.ge.100) then
        write(lulog,*) 'the overlap matrix after iopt = ',iopt
        call wrtmat2(smat,nold+nnew,nold+nnew,nold+nnew,nold+nnew)
      end if

      end do ! iopt

c dbg
c     allocate larger scr!!
c      print *,'eigenvalues of smat:'
c      xmat = smat
c      call rs(nold+nnew,nold+nnew,xmat,scr,0,
c     &        scr(2*(nold+nnew)+1),scr(3*(nold+nnew)+1),
c     &     scr(4*(nold+nnew)+1),idx)
c      print *,scr(1:nold+nnew)
c dbg

      ! call modified Gram-Schmidt:
      call mgs(xmat,smat,nold+nnew,scr)

      if (ntest.ge.100) then
        write(lulog,*) 'the reorth. matrix: '
        call wrtmat2(xmat,nold+nnew,nold+nnew,nold+nnew,nold+nnew)
      end if

      ! assemble new vectors

      ! dummy iord for ffnew:
      do inew = 1, nnew
        iordnew(inew) = inew
      end do

      do iopt = 1, nopt

      idx = nold
c dbg
c      print *,'iopt = ',iopt,'  idx = ',idx
c dbg
      do inew = 1, nnew
        ! check whether vector is nonzero
c dbg
c        print *,'test = ',ddot(nold+nnew,xmat(1,nold+inew),1,
c     &                     xmat(1,nold+inew),1),nold,nnew,inew
c dbg
        if (ddot(nold+nnew,xmat(1,nold+inew),1,
     &                     xmat(1,nold+inew),1).lt.1d-10)
     &       cycle
        idx = idx+1
c dbg
c        print *,'idx, ndim_sbsp: ',idx, ndim_sbsp
c dbg
        if (idx.le.ndim_sbsp) then
          ! overwrite previous (zero) record
          irec = ioptc_get_sbsp_rec(idx,iord_sbsp,ndim_sbsp,mxsbsp)
        else
          ! get new record
          irec = ioptc_get_sbsp_rec(0,iord_sbsp,ndim_sbsp,mxsbsp)
        end if
c dbg
c        print *,'idx = ',idx,'   irec = ',irec
c dbg

        if (ntest.ge.100)
     &       write(lulog,*) 'next new record: ',irec

        ! add contributions from new space
        call optc_expand_vec(xmat(nold+1,nold+inew),nnew,xdum,.false.,
     &                       ff_sbsp(iopt)%fhand,irec,0d0,
     &                                   ffnew(iopt)%fhand,iordnew,
     &                       nincore,nwfpar(iopt),lenbuf,xbuf1,xbuf2)
        ! add contributions from old space
        call optc_expand_vec(xmat(1,nold+inew),nold,xdum,.false.,
     &                      ff_sbsp(iopt)%fhand,irec,1d0,
     &                                  ff_sbsp(iopt)%fhand,iord_sbsp,
     &                       nincore,nwfpar(iopt),lenbuf,xbuf1,xbuf2)

      end do

      end do ! iopt

      nadd = idx-nold

      deallocate(smat,xmat,scr)

      return
      end

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public