text:
C***********************************************************************
C                                                                      *
C           LEGEND EDITOR FOR PROGRAMABLE DISPLAY PUSHBUTTONS          *
C                                  BY                                  *
C                            WILLIAM S. LUCK                           *
C                                                                      *
C                             MAIN PROGRAM                             *
C                                                                      *
C***********************************************************************
C
      CHARACTER*9 RFILE
C  RFILE is the name of the file to which micro switch messages are sent
      BYTE RESP
C  RESP is the numbered response to the main menu.
      BYTE TOPLEG(44), BOTLEG(44)
C  TOPLEG and BOTLEG are dummy arrays to be passed when calling PATMAP
C  subroutine.
      CHARACTER*4 TERM
C  TERM is the name of the RS-232 terminal port to which the PDPs are
C  connected.
      COMMON TERM
C  TERM is stored in a common block shared with the TALK subroutine.
      READ (1,5) TERM
5       FORMAT (A4)
      WRITE (6,*) 'ENTER FILE NAME'
      READ (5,10) RFILE
10      FORMAT (A9)
      OPEN (UNIT=4, FILE=RFILE, DEFAULTFILE='.PDP',
     $STATUS='UNKNOWN', ORGANIZATION='INDEXED', ACCESS='KEYED',
     $RECORDTYPE='VARIABLE', FORM='FORMATTED', RECL=15000,
     $KEY=(2:21:CHARACTER))
C  The indexed disk data file, RFILE, is opened and given the logic #4.
      OPEN (UNIT=2, FILE=RFILE, DEFAULTFILE='.TXT',
     $STATUS='UNKNOWN', ORGANIZATION='INDEXED', ACCESS='KEYED',
     $RECORDTYPE='VARIABLE', FORM='FORMATTED', RECL=15000,
     $KEY=(2:21:CHARACTER))
C  The indexed disk text file, RFILE, is opened and given the logic #2.
15    WRITE (6,20) RFILE, RFILE, RFILE
20      FORMAT (/////' ENTER NUMBER CORRESPONDING TO DESIRED COMMAND'//
     $'  1.  GENERAL CONTROL COMMAND'/
     $'  2.  SPECIFIC CONTROL COMMAND'/
     $'  3.  CHARACTER COMMAND'/
     $'  4.  BIT MAP COMMAND'/
     $'  5.  PATTERN MAP COMMAND'/
     $'  6.  LUMINANCE COMMAND'/
     $'  7.  EXECUTE MESSAGE FROM ',A9/
     $'  8.  DELETE MESSAGE FROM ',A9/
     $'  9.  LIST CONTENTS OF ',A9/
     $' 10.  END EDITING SESSION')
      READ (5,*,ERR=15) RESP
      IF (RESP .EQ. 1) THEN
        CALL GENCOM
C  GENCOM is the subroutine for general control commands.
      ELSE IF (RESP .EQ. 2) THEN
        CALL SPECCOM
C  SPECCOM is the subroutine for specific control commands.
      ELSE IF (RESP .EQ. 3) THEN
        CALL CHAR
C  CHAR is the subroutine for character control commands.
      ELSE IF (RESP .EQ. 4) THEN
        CALL BITMAP
C  BITMAP is the subroutine for bit map commands.
      ELSE IF (RESP .EQ. 5) THEN
        CALL PATMAP(TOPLEG,BOTLEG)
C  PATMAP is the subroutine for pattern map commands.
      ELSE IF (RESP .EQ. 6) THEN
        CALL LUMIN
C  LUMIN is the subroutine for luminance commands.
      ELSE IF (RESP .EQ. 7) THEN
        CALL EXECMES
C  EXECMES is the subroutine for executing messages from file.
      ELSE IF (RESP .EQ. 8) THEN
        CALL DELMES
C  DELMES is the subroutine for deleting messages from file.
      ELSE IF (RESP .EQ. 9) THEN
        CALL FILELIST(RFILE)
C  FILELIST is the subroutine for listing the contents of a file.
      ELSE IF (RESP .EQ. 10) THEN
        STOP
C  Program execution is terminated.
      END IF
      GO TO 15
C  Display main menu again and wait for response.
      END
C
C
C***********************************************************************
C                                                                      *
C                           SUBROUTINE GENCOM                          *
C                                                                      *
C***********************************************************************
C
      SUBROUTINE GENCOM
C  Subroutine for general control commands.
      BYTE LEGEND(44), NUMBER
C  LEGEND is the array that holds each byte of every message sent to a
C  micro switch.  The maximum number of bytes in a message is 44.
C  NUMBER is the number of bytes in a message.
      BYTE SUBRESP
C  SUBRESP is the numbered response to a subroutine menu.
5     NUMBER=1
C  NUMBER is 1 for all general command messages.
      WRITE (6,10)
10    FORMAT (///' ENTER NUMBER CORRESPONDING TO DESIRED GENERAL CONTROL
     $ COMMAND'//
     $' 1.  SELF TEST COMMAND'/
     $' 2.  CLEAR DISPLAY COMMAND'/
     $' 3.  CLEAR ERROR COMMAND'/
     $' 4.  END OF TRANSMISSION COMMAND'/
     $' 5.  RETRANSMIT COMMAND'/
     $' 6.  RETURN TO MAIN MENU')
      READ (5,*,ERR=5) SUBRESP
      IF (SUBRESP .EQ. 1) THEN
        LEGEND(1) = '88'X
C  Assign hex value of self test command to first byte of message.
      ELSE IF (SUBRESP .EQ. 2) THEN
        LEGEND(1) = '8C'X
C  Assign hex value of clear display command to first byte of message.
      ELSE IF (SUBRESP .EQ. 3) THEN
        LEGEND(1) = 'D0'X
C  Assign hex value of clear error command to first byte of message.
      ELSE IF (SUBRESP .EQ. 4) THEN
        LEGEND(1) = 'C0'X
C  Assign hex value of end of transmission command to first byte of mes.
      ELSE IF (SUBRESP .EQ. 5) THEN
        LEGEND(1) = 'C4'X
C  Assign hex value of retransmit command to first byte of message.
      ELSE IF (SUBRESP .EQ. 6) THEN
        RETURN
C  Return to main menu.
      ELSE
        GO TO 5
C  Rewrite subroutine menu if response number is out of bounds.
      END IF
      CALL SUMCHECK(LEGEND,NUMBER)
C  Call SUMCHECK subroutine to compute sumcheck byte, to send message to
C  micro switch, and to send message to data file.
      GO TO 5
C  Rewrite subroutine menu after return from SUMCHECK.
      END
C
C
C***********************************************************************
C                                                                      *
C                          SUBROUTINE SPECCOM                          *
C                                                                      *
C***********************************************************************
C
      SUBROUTINE SPECCOM
C  Subroutine for specific control commands.
      BYTE LEGEND(44), NUMBER
C  LEGEND is the array that holds each byte of every message sent to a
C  micro switch.  The maximum number of bytes in a message is 44.
C  NUMBER is the number of bytes in a message.
      BYTE SUBRESP, SWITCH
C  SUBRESP is the numbered response to the subroutine menu.
C  SWITCH is the number of the particular switch to which the message is
C  sent.
5     NUMBER=1
C  NUMBER is 1 for all specific control commmands.
      WRITE (6,10)
10    FORMAT (///' ENTER NUMBER CORRESPONDING TO DESIRED SPECIFIC CONTRO
     $L COMMAND AND <RETURN>,'/' THEN ENTER NUMBER FROM 0 TO 3 OF SWITCH
     $ TO RECEIVE COMMAND'//
     $' 1.  BLINK SWITCH COMMAND'/
     $' 2.  CLEAR SWITCH COMMAND'/
     $' 3.  STOP BLINK SWITCH COMMAND'/
     $' 4.  RETURN TO MAIN MENU')
      READ (5,*,ERR=5) SUBRESP
      IF (SUBRESP .EQ. 4) RETURN
C  Return to main menu.
      READ (5,*,ERR=5) SWITCH
      IF (SWITCH.LT. 0 .OR. SWITCH.GT.3) GO TO 5
C  If switch number is out of bounds, rewrite subroutine menu.
      IF (SUBRESP .EQ. 1) THEN
        LEGEND(1) = '90'X + SWITCH
C  Add blink switch command to switch number.
      ELSE IF (SUBRESP .EQ. 2) THEN
	LEGEND(1) = '98'X + SWITCH
C  Add clear switch command to switch number.
      ELSE IF (SUBRESP .EQ. 3) THEN
        LEGEND(1) = 'D4'X + SWITCH
C  Add stop blink switch command to switch number.
      ELSE
        GO TO 5
C  Rewrite subroutine menu if response number is out of bounds.
      END IF
      CALL SUMCHECK(LEGEND,NUMBER)
C  Call SUMCHECK subroutine to compute sumcheck byte, to send message to
C  micro switches, and to send message to data file.
      GO TO 5
C  Rewrite subroutine menu after return from SUMCHECK.
      END
C
C
C***********************************************************************
C                                                                      *
C                            SUBROUTINE CHAR                           *
C                                                                      *
C***********************************************************************
C
      SUBROUTINE CHAR
C  Subroutine for character control commands.
      BYTE LEGEND(44), NUMBER
C  LEGEND is the array that holds each byte of every message sent to a
C  micro switch.  The maximum number of bytes in a message is 44.
C  NUMBER is the number of bytes in a message.
      BYTE SUBRESP, SWITCH, MAXN, N
C  SUBRESP is the numbered response to a subroutine menu.
C  SWITCH is the number of the particular switch to which the message is
C  sent.
C  MAXN is the maximum number of characters for a particular FONT size.
C  N is the number of characters in the character string to be sent.
      CHARACTER*1 ASCII(6)
C  ASCII is the array that holds the character string to be sent.
5     WRITE (6,10)
10    FORMAT (///' ENTER NUMBER CORRESPONDING TO DESIRED CHARACTER COMMA
     $ND'//
     $' 1.  5X7 FONT CHARACTER COMMAND'/
     $' 2.  10X14 FONT CHARACTER COMMAND'/
     $' 3.  RETURN TO MAIN MENU')
      READ (5,*,ERR=5) SUBRESP
      IF (SUBRESP .EQ. 1) THEN
        LEGEND(1) = '84'X
C  Assign hex value of 5x7 FONT character command to first byte of
C  message.
        MAXN=7
C  Maximum number of 5x7 characters is less than 7.
      ELSE IF (SUBRESP .EQ. 2) THEN
        LEGEND(1) = 'A8'X
C  Assign hex value of 10x14 Font character command to first byte of
C  message.
        MAXN=4
C  Maximum number of 10x14 characters is less than 4.
      ELSE IF (SUBRESP .EQ. 3) THEN
        RETURN
C  Return to main menu.
      ELSE
        GO TO 5
C  Rewrite subroutine menu if response number is out of bounds.
      END IF
20    WRITE (6,25) MAXN
25    FORMAT (' ENTER SWITCH # FROM 0-3, X COORD, Y COORD, AND # OF CHAR
     $S < ',I1)
      READ (5,*,ERR=20) SWITCH, LEGEND(2), LEGEND(3), N
      IF (SWITCH.LT.0 .OR. SWITCH.GT.3) GO TO 20
C  Prompt for new input if switch number is out of bounds.
      IF (LEGEND(2).LT.0 .OR. LEGEND(2).GT.34) GO TO 20
C  Prompt for new input if x coordinate is out of bounds.
      IF (LEGEND(3).LT.0 .OR. LEGEND(3).GT.15) GO TO 20
C  Prompt for new input if y coordinate is out of bounds.
      IF (N.LT.1 .OR. N.GE.MAXN) GO TO 20
C  Prompt for new input if number of characters is out of bounds.
      LEGEND(1) = LEGEND(1) + SWITCH
C  Add FONT command to switch number and assign to first byte of
C  message.
      WRITE (6,30) N
30      FORMAT (' ENTER ',I1,' CHARACTER STRING')
      READ (5,40) (ASCII(I), I=1,N)
C  Load each character of the string into the array ASCII.
40      FORMAT (<N>(A1))
      DO 50 I=1,N
        LEGEND(I+3) = ICHAR(ASCII(I))
C  The decimal equivalent of each character in the string is assigned to
C  a particular byte in the message.
50    CONTINUE
      NUMBER=N+3
C  The number of bytes in the message is equal to the number of
C  characters in the string plus three.
      CALL SUMCHECK(LEGEND,NUMBER)
C  Call SUMCHECK subroutine to compute sumcheck byte, to send message to
C  micro switches, and to send message to data file.
      GO TO 5
C  Rewrite subroutine menu after return from SUMCHECK.
      END
C
C
C***********************************************************************
C                                                                      *
C                           SUBROUTINE BITMAP                          *
C                                                                      *
C***********************************************************************
C
      SUBROUTINE BITMAP
C  Subroutine for bit map commands.
      BYTE LEGEND(44), NUMBER
C  LEGEND is the array that holds each byte of every message sent to a
C  micro switch.  The maximum number of bytes in a message is 44.
C  NUMBER is the number of bytes in a message.
      BYTE SUBRESP, SWITCH, PAIRS, X, Y
C  SUBRESP is the numbered response to a subroutine menu.
C  SWITCH is the number of the particular switch to which the message is
C  sent.
C  PAIRS is the number of lines to be drawn by the message.
C  X is the x coordinate for a particular line.
C  Y is the y coordinate for a particular line.
5     WRITE (6,10)
10    FORMAT (///' ENTER NUMBER CORRESPONDING TO DESIRED BIT MAP COMMAND
     $'//
     $' 1.  HORIZONTAL BIT MAP COMMAND'/
     $' 2.  VERTICAL BIT MAP COMMAND'/
     $' 3.  RETURN TO MAIN MENU')
      READ (5,*,ERR=5) SUBRESP
      IF (SUBRESP.EQ.1 .OR. SUBRESP.EQ.2) THEN
C  If either a horizontal or a vertical bit map is selected, the
C  following block is executed.
30      WRITE (6,*) 'ENTER SWITCH # FROM 0-3, LINE LENGTH, AND # OF X,Y
     $COORD PAIRS < 16'
        READ (5,*,ERR=30) SWITCH, LEGEND(2), PAIRS
C  The line length is assigned to LEGEND(2) which is the second byte of
C  the message.
        IF (SWITCH.LT.0 .OR. SWITCH.GT.3) GO TO 30
C  Prompt for new input if switch number is out of bounds.
        IF (LEGEND(2).LT.1 .OR. LEGEND(2).GT.63) GO TO 30
C  Prompt for new input if line length is out of bounds.
        IF (PAIRS.LT.1 .OR. PAIRS.GT.15) GO TO 30
C  Prompt for new input if number of lines is out of bounds.
        DO 40 I=1,PAIRS
C  The DO loop reads the x and y coordinates of each line and assigns
C  each coordinate to a byte in the message.
32        WRITE (6,35) I
35        FORMAT (' ENTER X AND Y COORDINATES OF START POINT FOR PAIR ',
     $I2)
          READ (5,*,ERR=32) X, Y
          IF (X.LT.0 .OR. X.GT.34 .OR. Y.LT.0 .OR. Y.GT.15) GO TO 32
C  Prompt for new input if x or y coordinate is out of bounds.
          LEGEND(2*I+1)=X
          LEGEND(2*I+2)=Y
C  Increment the position of the LEGEND array and assign the x and y
C  coordinates to the LEGEND array.
40      CONTINUE
      ELSE IF (SUBRESP .EQ. 3) THEN
        RETURN
C  Return to main menu.
      ELSE
        GO TO 5
C  Rewrite subroutine menu if response number is out of bounds.
      END IF
      IF (SUBRESP .EQ. 1) THEN
        LEGEND(1) = '9C'X + SWITCH
C  Add horizontal bit map command to switch number and assign to first
C  byte of message.
      ELSE IF (SUBRESP .EQ. 2) THEN
        LEGEND(1) = 'A0'X + SWITCH
C  Add vertical bit map command to switch number and assign to first
C  byte of message.
      END IF
        NUMBER=2*(PAIRS+1)
C  The number of bytes in the message is equal to twice the number of
C  lines plus two.
      CALL SUMCHECK(LEGEND,NUMBER)
C  Call SUMCHECK subroutine to compute sumcheck byte, to send message to
C  micro switches, and to send message to data file.
      GO TO 5
C  Rewrite subroutine menu after return from SUMCHECK.
      END
C
C
C***********************************************************************
C                                                                      *
C                           SUBROUTINE PATMAP                          *
C                                                                      *
C***********************************************************************
C
      SUBROUTINE PATMAP(TOPLEG,BOTLEG)
C  Subroutine for pattern map commands.
      INTEGER*4 SYS$ASSIGN
C  SYS$ASSIGN assigns a channel to an I/O process.
      INTEGER*4 SYS$QIOW
C  SYS$QIOW is queue I/O request and wait for event flag.
      BYTE TOPLEG(44), BOTLEG(44)
C  TOPLEG is the array that holds each byte of every message sent to the
C  top half of the display.
C  BOTLEG is the array that holds each byte of every message sent to the
C  bottom half of the display.
      INTEGER*2 ITOPLEG(44), IBOTLEG(44)
C  ITOPLEG and IBOTLEG are INTEGER arrays set equal to the BYTE arrays
C  TOPLEG and BOTLEG respectively for the purpose of performing BIT
C  operations in the subroutines MAPSET and MAPINP.
      BYTE NUMBER, SWITCH, FLAG
C  NUMBER is the number of bytes in a message.
C  SWITCH is the number of the particular switch to which the message is
C  sent.
C  FLAG is 1 when displaying previous legend and is 0 when clearing the
C  screen.
      INTEGER*2 X, Y
C  X is the x coordinate for a particular LED.
C  Y is the y coordinate for a particular LED.
      CHARACTER*3 CURSOR
C  CURSOR holds response from keyboard while editing legend.
      INCLUDE '($IODEF)'
C  Access $IODEF module in the FORSYSDEF library.
      FLAG=1
C  Initialize FLAG to 1.
1     IFLAG=LIB$ERASE_PAGE(%REF(1),%REF(1))
      IF(.NOT.IFLAG) WRITE (6,*) 'ERASE PAGE FAILURE'
C  Clear screen and check status.
      KFLAG=LIB$PUT_LINE('   0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1
     $ 2 3 4 5 6 7 8 9 0 1 2 3 4',%REF(1),)
      IF(.NOT.KFLAG) WRITE (6,*) 'PUT LINE FAILURE'
C  Display x coordinates and check status.
      WRITE (6,*)
      JFLAG=LIB$SET_CURSOR(%REF(1),%REF(1))
      IF(.NOT.JFLAG) WRITE (6,*) 'SET CURSOR FAILURE'
C  The previous 3 lines are needed for transition from PUT_LINE to WRITE
      DO 10 I=1,16
        WRITE (6,5) I-1
5         FORMAT (X,I2)
C  Display y coordinates.
10    CONTINUE
      WRITE (6,*)
      WRITE (6,*) 'USE ARROW KEYS TO MOVE CURSOR TO DESIRED COORDINATES'
      WRITE (6,*) 'STRIKE THE RETURN KEY TO TOGGLE A COORDINATE PAIR'
      WRITE (6,*) 'PF1 = SEND MESSAGE; PF2 = CLEAR SCREEN; PF3 = RETURN
     $ TO MAIN MENU'
      JERRFLAG=SYS$ASSIGN('TT:',JOHN,,)
C  Assign channel to I/O process and check status.
      IF(.NOT.JERRFLAG) WRITE (6,*) 'ASSIGN FAILURE'
C  JOHN is the channel number assigned by the system service.
      IF (FLAG .EQ. 1) THEN
C  If the flag is 1, then equate the BYTE arrays to their respective
C  INTEGER arrays.
        DO 15 I=1,44
          ITOPLEG(I)=TOPLEG(I)
          IBOTLEG(I)=BOTLEG(I)
15      CONTINUE
        CALL MAPSET(ITOPLEG,0)
C  Call MAPSET subroutine to display the top half of the legend as
C  specified by the array ITOPLEG.
        CALL MAPSET(IBOTLEG,8)
C  Call MAPSET subroutine to display the bottom half of the legend as
C  specified by the array IBOTLEG.
      END IF
      FLAG=1
C  Reset value of FLAG to 1.
      X=4
      Y=2
C  X and Y are initialized to 4 and 2 respectively in order to set the
C  cursor at its beginning position.
17    JFLAG=LIB$SET_CURSOR(%REF(Y),%REF(X))
      IF(.NOT.JFLAG) WRITE (6,*) 'SET CURSOR FAILURE'
C  Set cursor to present X and Y position and check status.
      KERRFLG=SYS$QIOW(%VAL(0),%VAL(JOHN),%VAL(IO$_READVBLK),,,,
     $%REF(CURSOR),%VAL(3),,,,)
C  Read response from keyboard and check status.
C  0 is the flag to be set upon completion of the service.
C  JOHN is the channel number for the service.
C  IO$_READVBLK means to read a virtual block.
C  CURSOR is the address that is to receive the data read.
C  3 is the number of bytes to be read.
      IF(.NOT.KERRFLG) WRITE (6,*) 'READ FAILURE'
      IF (CURSOR .EQ. CHAR(27)//CHAR(91)//CHAR(65)) THEN
C  Decrement Y if the UP key is pressed; but do not let Y become < 2.
        Y=Y-1
        IF (Y .LT. 2) Y=2
      ELSE IF (CURSOR .EQ. CHAR(27)//CHAR(91)//CHAR(66)) THEN
C  Increment Y if the DOWN key is pressed; but do not let Y become > 17.
        Y=Y+1
        IF (Y .GT. 17) Y=17
      ELSE IF (CURSOR .EQ. CHAR(27)//CHAR(91)//CHAR(68)) THEN
C  Decrement X if the LEFT key is pressed; but do not let X become < 4.
        X=X-2
        IF (X .LT. 4) X=4
      ELSE IF (CURSOR .EQ. CHAR(27)//CHAR(91)//CHAR(67)) THEN
C  Increment X if the RIGHT key is pressed; but do not let X become > 72
        X=X+2
        IF (X .GT. 72) X=72
      ELSE IF (CURSOR .EQ. CHAR(27)//CHAR(79)//CHAR(80)) THEN
C  If PF 1 key is pressed, this block sends a message to micro switch.
        WRITE (6,*)
        JFLAG=LIB$SET_CURSOR(%REF(21),%REF(1))
        IF(.NOT.JFLAG) WRITE (6,*) 'SET CURSOR FAILURE'
C  The previous 3 lines are needed for transition from PUT_LINE to WRITE
20      WRITE (6,*) 'ENTER NUMBER FROM 0 TO 3 OF SWITCH TO RECEIVE MESSA
     $GE'
      READ (5,*,ERR=20) SWITCH
      IF (SWITCH.LT.0 .OR. SWITCH.GT.3) GO TO 20
C  Prompt for new input if switch number is out of bounds.
        DO 30 I=1,44
          TOPLEG(I)=ITOPLEG(I)
          BOTLEG(I)=IBOTLEG(I)
C  Equate the BYTE arrays to their respective INTEGER arrays.
30      CONTINUE
        TOPLEG(1) = 'AC'X + SWITCH
C  Add pattern map command to switch number and assign to first byte of
C  message.
        TOPLEG(2) = '00'X
C  Assign the value of the starting y position for the top half of the
C  display to the second byte of message.
        NUMBER=42
C  NUMBER is 42 for all pattern map commands.
        CALL SUMCHECK(TOPLEG,NUMBER)
C  Call SUMCHECK subroutine to compute sumcheck byte, to send message to
C  micro switch, and to send message to data file.
        BOTLEG(1) = 'AC'X + SWITCH
C  Add pattern map command to switch number and assign to first byte of
C  message.
        BOTLEG(2) = '08'X
C  Assign the value of the starting y position for the bottom half of
C  the display to the second byte of message.
        NUMBER=42
C  NUMBER is 42 for all pattern map commands.
        CALL SUMCHECK(BOTLEG,NUMBER)
C  Call SUMCHECK subroutine to compute sumcheck byte, to send message to
C  micro switch, and to send message to data file.
        GO TO 1
C  Go to clear screen.
      ELSE IF (CURSOR .EQ. CHAR(27)//CHAR(79)//CHAR(81)) THEN
C  If PF 2 key is pressed, this block clears the screen by setting the
C  value of the ITOPLEG and IBOTLEG arrays to 0.
        DO 40 I=1,44
          ITOPLEG(I)=0
          IBOTLEG(I)=0
40      CONTINUE
        FLAG=0
C  FLAG is set to 0 so the bits of the ITOPLEG and IBOTLEG arrays will
C  not be polled.
        GO TO 1
C  Go to clear screen.
      ELSE IF (CURSOR .EQ. CHAR(27)//CHAR(79)//CHAR(82)) THEN
C  If PF 3 key is pressed, this block returns control to the calling
C  routine.
        DO 50 I=1,44
          TOPLEG(I)=ITOPLEG(I)
          BOTLEG(I)=IBOTLEG(I)
C  Equate the INTEGER arrays to their respective BYTE arrays.
50      CONTINUE
      IFLAG=LIB$ERASE_PAGE(%REF(1),%REF(1))
      IF(.NOT.IFLAG) WRITE (6,*) 'ERASE PAGE FAILURE'
C  Clear screen and check status.
        RETURN
C  Return to calling routine.
      ELSE IF (CURSOR(1:1) .EQ. CHAR(13)) THEN
C  If RETURN key is pressed, this block toggles the display of the
C  present x,y position with either an * or a space.
        JFLAG=LIB$SET_CURSOR(%REF(Y),%REF(X))
        IF(.NOT.JFLAG) WRITE (6,*) 'SET CURSOR FAILURE'
C  Set cursor and check status.
        IF (Y .LT. 10) THEN
C  Call MAPINP subroutine to decide whether to put an * or a space at
C  the x and y coordinates passed to it.
          CALL MAPINP(ITOPLEG,X/2-2,Y-2)
        ELSE
          CALL MAPINP(IBOTLEG,X/2-2,Y-10)
        END IF
      END IF
      GO TO 17
C  Go and wait for another key to be pressed.
      END
C
C
C***********************************************************************
C                                                                      *
C                           SUBROUTINE MAPSET                          *
C                                                                      *
C***********************************************************************
C
      SUBROUTINE MAPSET(LEGEND,Z)
C  Subroutine MAPSET puts an * on the screen whenever a bit of the
C  LEGEND array is set.
      INTEGER*2 LEGEND(44)
C  LEGEND is the array on which BIT operations are performed.
      INTEGER*2 I, REM
C  I is the array element of LEGEND to be operated on.
C  REM holds the remainder and specifies the bit position to be operated
C  on.
      INTEGER*2 X, Y, Z
C  X is the x coordinate for a particular LED.
C  Y is the y coordinate for a particular LED.
C  Z is 0 when LEGEND represents the top half of the screen and is 8
C  when LEGEND represents the bottom half of the screen.
      DO 20 Y=0,7
        DO 10 X=0,34
C  The purpose of the two DO loops is to poll every bit of the LEGEND
C  array
          I = 5*Y + 3 + X/7
C  Compute the value of I.  Every 35 position line is stored in 5 bytes
C  representing 7 positions each.  The data bytes begin in the third
C  byte of a message.
          REM=MOD(X+7,7)
C  Compute the remainder of (X+7)/7
          IF (BTEST(LEGEND(I),REM)) THEN
C  See if particular bit is set.  If the bit is set, set the cursor to
C  the specified x,y position and put an * there.
            JFLAG=LIB$SET_CURSOR(%REF(Y+Z+2),%REF(2*X+4))
            IF(.NOT.JFLAG) WRITE (6,*) 'SET CURSOR FAILURE'
C  Set cursor and check status.
            KFLAG=LIB$PUT_LINE('*',%REF(0),)
            IF(.NOT.KFLAG) WRITE (6,*) 'PUT LINE FAILURE'
C  Put an * at the cursor position and check status.
          END IF
10      CONTINUE
20    CONTINUE
      RETURN
C  Return to calling routine.
      END
C
C
C***********************************************************************
C                                                                      *
C                           SUBROUTINE MAPINP                          *
C                                                                      *
C***********************************************************************
C
      SUBROUTINE MAPINP(LEGEND,X,Y)
C  Subroutine MAPINP decides whether to put an * or a space at the x and
C  y coordinates passed to it.  It then increments or decrements the
C  array element of LEGEND accordingly.
      INTEGER*2 LEGEND(44)
C  LEGEND is the array on which BIT operations are performed.
      INTEGER*2 X, Y
C  X is the x coordinate for a particular LED.
C  Y is the y coordinate for a particular LED.
      INTEGER*2 I, REM
C  I is the array element of LEGEND to be operated on.
C  REM holds the remainder and specifies the bit position to be operated
C  on.
      I = 5*Y + 3 + X/7
C  Compute the value of I.  Every 35 position line is stored in 5 bytes
C  representing 7 positions each.  The data bytes begin in the third
C  byte of a message.
      REM=MOD(X+7,7)
C  Compute the remainder of (X+7)/7.
      IF(.NOT. BTEST(LEGEND(I),REM)) THEN
C  See if particular bit is set. If bit is not set, set the bit and put
C  an * in the x,y position.
        LEGEND(I)=IBSET(LEGEND(I),REM)
        KFLAG=LIB$PUT_LINE('*',%REF(0),)
      ELSE
C  If the bit is set, clear the bit and put a space in the x,y position.
        LEGEND(I)=IBCLR(LEGEND(I),REM)
        KFLAG=LIB$PUT_LINE(' ',%REF(0),)
      END IF
      IF(.NOT.KFLAG) WRITE (6,*) 'PUT LINE FAILURE'
      RETURN
C  Return to calling routine.
      END
C
C
C***********************************************************************
C                                                                      *
C                           SUBROUTINE LUMIN                           *
C                                                                      *
C***********************************************************************
C
      SUBROUTINE LUMIN
C  Subroutine for luminance commands.
      BYTE LEGEND(44), NUMBER
C  LEGEND is the array that holds each byte of every message sent to a
C  micro switch.  The maximum number of bytes in a message is 44.
C  NUMBER is the number of bytes in a message.
      BYTE SUBRESP
C  SUBRESP is the numbered response to a subroutine menu.
5     NUMBER=2
C  NUMBER is 2 for all luminance command messages.
      WRITE (6,10)
10    FORMAT (///' CHOOSE A RELATIVE LUMINANCE, RANGING FROM 0 TO 36,'/
     $' WITH 0 BEING MINIMUM LUMINANCE AND 36 BEING MAXIMUM LUMINANCE.'/
     $' ENTER ANY OTHER BYTE NUMBER TO RETURN TO MAIN MENU')
      READ (5,*,ERR=5) SUBRESP
      IF (SUBRESP.LT.0 .OR. SUBRESP.GT.36) RETURN
C  Return to main menu if response number is out of bounds.
      LEGEND(1)= 'BC'X
C  Assign the hexadecimal value of a luminance command to the first
C  byte of the message.
      LEGEND(2) = '36'X + 2*SUBRESP
C  The second byte of the message is assigned a relative value depending
C  on the desired brightness.
      CALL SUMCHECK(LEGEND,NUMBER)
C  Call SUMCHECK subroutine to compute sumcheck byte, to send message to
C  micro switch, and to send message to data file.
      GO TO 5
C  Rewrite subroutine menu after return from SUMCHECK.
      END
C
C
C***********************************************************************
C                                                                      *
C                          SUBROUTINE EXECMES                          *
C                                                                      *
C***********************************************************************
C
      SUBROUTINE EXECMES
C  Subroutine for executing messages from file.
      BYTE LEGEND(44), SECLEG(44), NUMBER
C  LEGEND is the array that holds each byte of every message sent to a
C  micro switch.  The maximum number of bytes in a message is 44.
C  SECLEG is the array that holds the second half of a pattern map
C  legend.
C  NUMBER is number of bytes in a message.
      CHARACTER*20 NAME
C  NAME is an up to 20 character string assigned by the user to be
C  executed from the data file.
      WRITE (6,10)
10      FORMAT (///' ENTER NAME OF MESSAGE TO BE EXECUTED')
      READ (5,20) NAME
20      FORMAT (A20)
      READ (4,30,KEY=NAME,ERR=40) NUMBER, (LEGEND(I), I=1,NUMBER)
C  Read the number of bytes and values of the bytes of the output
C  message in the data file with the specified name.
30      FORMAT (22X,I2,2X,<NUMBER>(X,I4))
      IF (LEGEND(1).GE.-84 .AND. LEGEND(1).LE.-81) THEN
C  If the first byte of a message has a value between -84 and -81, read
C  the other half of the pattern map legend.
        WRITE (6,*) 'ENTER NAME OF OTHER HALF OF PATTERN MAP LEGEND'
        READ (5,35) NAME
35        FORMAT (A20)
        READ (4,30,KEY=NAME,ERR=40) NUMBER, (SECLEG(I), I=1,NUMBER)
        IF (SECLEG(1).LT.-84 .OR. SECLEG(1).GT.-81) GO TO 50
C  Make sure second message read is a pattern map command.
        IF (LEGEND(2).EQ.0) THEN
C  Decide which array is the top half of the legend and which is the
C  bottom.
          CALL PATMAP(LEGEND,SECLEG)
        ELSE IF (LEGEND(2).EQ.8) THEN
          CALL PATMAP(SECLEG,LEGEND)
C  Call PATMAP subroutine to display the legend display on the terminal
C  screen.
        END IF
      ELSE
        CALL TALK(LEGEND,NUMBER)
C  Call TALK subroutine to send to and receive messages from the micro
C  switches using the QIOW system service.
      END IF
      RETURN
C  Return to main menu.
40    WRITE (6,*) 'NAME DOES NOT EXIST ON FILE'
C  Action to be taken on error from indexed READ statement.
      RETURN
C  Return to main menu.
50    WRITE (6,*) 'NAME DOES NOT REPRESENT A PATTERN MAP MESSAGE'
C  Action to be taken when message read is not second half of pattern
C  map legend.
      RETURN
C  Return to main menu.
      END
C
C
C***********************************************************************
C                                                                      *
C                           SUBROUTINE DELMES                          *
C                                                                      *
C***********************************************************************
C
      SUBROUTINE DELMES
C  Subroutine for deleting messages from file.
      CHARACTER*20 NAME
C  NAME is an up to 20 character string assigned by the user to be
C  deleted from the data file.
      WRITE (6,10)
10      FORMAT (///' ENTER NAME OF MESSAGE TO BE DELETED')
      READ (5,20) NAME
20      FORMAT (A20)
      READ (4,30,KEY=NAME,ERR=40)
C  Before a record can be deleted from an indexed data file, it must
C  first be read.
      DELETE (4)
C  Delete the specified record from the data file.
      READ (2,30,KEY=NAME)
      DELETE (2)
C  Delete the specified record from the text file.
30    FORMAT ()
C  Dummy FORMAT statement.
      RETURN
C  Return to main menu.
40    WRITE (6,*) 'NAME DOES NOT EXIST ON FILE.'
C  Action to be taken on error from indexed READ statement.
      RETURN
C  Return to main menu.
      END
C
C
C***********************************************************************
C                                                                      *
C                          SUBROUTINE FILELIST                         *
C                                                                      *
C***********************************************************************
C
      SUBROUTINE FILELIST(RFILE)
C  Subroutine for listing the contents of a file.
      CHARACTER*9 RFILE
C  RFILE is the name of the file to which micro switch messages are sent
      CHARACTER*20 NAME
C  NAME is an up to 20 character string assigned by the user to a
C  message.
      WRITE (6,10) RFILE
10      FORMAT (///X,A9/)
      READ (2,20,KEYGE=' ',ERR=30) NAME
C  Go to beginning of file.
      WRITE (6,*) NAME
15    READ (2,20,ERR=30) NAME
20      FORMAT (X,A20)
      WRITE (6,*) NAME
      GO TO 15
C  Go to read another record from file.
30    RETURN
C  Return to main menu after reading entire file.
      END
C
C
C***********************************************************************
C                                                                      *
C                          SUBROUTINE SUMCHECK                         *
C                                                                      *
C***********************************************************************
C
      SUBROUTINE SUMCHECK(LEGEND,NUMBER)
C  Subroutine SUMCHECK computes the sumcheck byte, sends message to
C  micro switch, and sends message to data file.
      BYTE NUMBER
C  NUMBER is number of bytes in a message.
      BYTE LEGEND(NUMBER+2)
C  Redimension the array LEGEND to be the size of the number of bytes
C  sent by the calling routine plus the end of message byte and sumcheck
C  byte.  LEGEND is the array that holds each byte of every message sent
C  to a micro switch.  The maximum number of bytes in a message is 44.
      INTEGER*4 TOTAL
C  TOTAL is the sum of each byte of a message.
      CHARACTER*20 NAME
C  NAME is an up to 20 character string assigned by the user to be the
C  title of the message sent to the data file.
      LEGEND(NUMBER+1) = '80'X
C  The second to the last byte of a message is assign the hexadecimal
C  value for an end of message byte.
      TOTAL=0
C  TOTAL is initialized at 0.
      DO 10 I=1,NUMBER+1
        TOTAL=TOTAL+LEGEND(I)
C  The sum of the values of each byte of a message is assigned to TOTAL.
10    CONTINUE
      NUMBER=NUMBER+2
C  The number of bytes in a message is increased by two in order to
C  include the end of message byte and the sumcheck byte.
11    IF (TOTAL .LT. -128) THEN
        TOTAL=TOTAL+256
      ELSE IF (TOTAL .GT. 127) THEN
        TOTAL=TOTAL-256
      END IF
      IF (TOTAL.LT.-128 .OR. TOTAL.GT.127) GO TO 11
C  TOTAL is adjusted so that its value can be represented by an eight
C  bit byte.  Therefore the carry bit is ignored.
      LEGEND(NUMBER)=TOTAL
C  The last byte of a message is the sumcheck byte and is assigned the
C  same value of TOTAL.
C
      CALL TALK(LEGEND,NUMBER)
C  Call TALK subroutine to send to and receive messages from the micro
C  switches using the QIOW system service.
C
15    WRITE (6,*) 'GIVE 20 CHARACTER OR LESS NAME TO MESSAGE'
      WRITE (6,*) 'ENTER <RETURN> TO NOT SEND MESSAGE TO DATA FILE'
      READ (5,20) NAME
20      FORMAT (A20)
      IF (NAME .EQ. ' ') RETURN
C  If user does not wish to send message to data file, control is
C  returned to the calling routine.
      WRITE (4,30,ERR=40) NAME, NUMBER, (LEGEND(I), I=1,NUMBER)
C  Write the title, the number of bytes in the message, and the decimal
C  equivalent of each byte of the message to the data file opened in
C  the main program.  If an error occurs, branch to specified line.
30      FORMAT (X,A20,X,I2,2X,<NUMBER>(X,I4))
      WRITE (2,35) NAME
C  Write the title of the message to the text file opened in the main
C  program.
35      FORMAT (X,A20)
      RETURN
C  Return to calling routine.
40    WRITE (6,*) CHAR(7), 'DUPLICATE NAME; TRY AGAIN.'
C  Action to be taken on error from indexed WRITE statement.
      GO TO 15
C  Go to ask for another name.
      END
C
C
C***********************************************************************
C                                                                      *
C                            SUBROUTINE TALK                           *
C                                                                      *
C***********************************************************************
C
      SUBROUTINE TALK(LEGEND,NUMBER)
C  Subroutine TALK sends to and receives messages from the micro
C  switches using the QIOW system service.
      INTEGER*4 SYS$ASSIGN
C  SYS$ASSIGN assigns a channel to an I/O process.
      INTEGER*4 SYS$QIOW
C  SYS$QIOW is queue I/O request and wait for event flag.
      BYTE NUMBER
C  NUMBER is number of bytes in a message.
      BYTE LEGEND(NUMBER)
C  Redimension the array LEGEND to be the size of the number of bytes
C  sent by the calling routine plus the end of message byte and sumcheck
C  byte.  LEGEND is the array that holds each byte of every message sent
C  to a micro switch.  The maximum number of bytes in a message is 44.
      BYTE OUTMES
C  OUTMES is the response byte sent by the LRCU.
      CHARACTER*4 TERM
C  TERM is the name of the RS-232 terminal port to which the PDPs are
C  connected.
      COMMON TERM
C  TERM is stored in a common block shared with the main program.
C
      INCLUDE '($IODEF)'
C  Access $IODEF module in the FORSYSDEF library.
C
      JERRFLAG=SYS$ASSIGN(%DESCR(TERM),JOHN,,)
C  Assign channel to I/O process and check status.
C  JOHN is the channel number assigned by the system service.
      IF(.NOT.JERRFLAG) WRITE (6,*) 'ASSIGN FAILURE'
C
      KERRFLAG=SYS$QIOW(%VAL(0),%VAL(JOHN),%VAL(IO$_WRITEVBLK),,,,
     $%REF(LEGEND),%VAL(NUMBER),,,,)
C  Write command message to LRCU and check status.
C  0 is the flag to be set upon completion of the service.
C  JOHN is the channel number for the service.
C  IO$_WRITEVBLK means to write a virtual block.
C  LEGEND represents the block to be written.
C  NUMBER is the number of bytes to be written.
      IF(.NOT.KERRFLAG) WRITE (6,*) 'WRITE FAILURE'
C
10    LERRFLAG=SYS$QIOW(%VAL(1),%VAL(JOHN),%VAL(IO$_READVBLK),,,,
     $%REF(OUTMES),%VAL(1),,,,)
C  Read response byte from LRCU and check status.
C  1 is the flag to be set upon completion of the service.
C  JOHN is the channel number for the service.
C  IO$_READVBLK means to read a virtual block.
C  OUTMES is the address that is to receive the data read.
C  1 is the number of bytes to be read.
      IF(.NOT.LERRFLAG) WRITE (6,*) 'READ FAILURE'
C
      IF (OUTMES .EQ. 16) THEN
        WRITE (6,*) 'MESSAGE ACKNOWLEDGED'
        RETURN
C  Return to calling routine.
      ELSE IF (OUTMES .EQ. 4) THEN
        WRITE (6,*) 'SELF TEST PASSED'
      ELSE IF (LEGEND(1) .EQ. '88'X) THEN
        CALL TESTFAIL(OUTMES)
C  Call TESTFAIL subroutine to process self test failure messages.
      ELSE
        WRITE (6,*) OUTMES
      END IF
      GO TO 10
C  Go to get more output.
      END
C
C
C***********************************************************************
C                                                                      *
C                          SUBROUTINE TESTFAIL                         *
C                                                                      *
C***********************************************************************
C
      SUBROUTINE TESTFAIL(OUTMES)
C  Subroutine TESTFAIL processes self test failure messages.
      BYTE OUTMES
C  OUTMES is the response byte sent by the LRCU.
      IF (OUTMES .EQ. 'A1'X) WRITE (6,*) 'ROM FAILURE'
      IF (OUTMES .EQ. 'A2'X) WRITE (6,*) 'RAM FAILURE'
      IF (OUTMES .EQ. '23'X) WRITE (6,*) 'REFRESH CONTROLLER CHIP FAILUR
     $E'
      IF (OUTMES .EQ. 'A4'X) WRITE (6,*) 'SWITCH 0 OUTPUT FAILURE'
      IF (OUTMES .EQ. '25'X) WRITE (6,*) 'SWITCH 1 OUTPUT FAILURE'
      IF (OUTMES .EQ. '26'X) WRITE (6,*) 'SWITCH 2 OUTPUT FAILURE'
      IF (OUTMES .EQ. 'A7'X) WRITE (6,*) 'SWITCH 3 OUTPUT FAILURE'
      IF (OUTMES .EQ. 'A8'X) WRITE (6,*) 'REFRESH SWITCH 0 MEMORY FAILUR
     $E'
      IF (OUTMES .EQ. '29'X) WRITE (6,*) 'REFRESH SWITCH 1 MEMORY FAILUR
     $E'
      IF (OUTMES .EQ. '2A'X) WRITE (6,*) 'REFRESH SWITCH 2 MEMORY FAILUR
     $E'
      IF (OUTMES .EQ. 'AB'X) WRITE (6,*) 'REFRESH SWITCH 3 MEMORY FAILUR
     $E'
      RETURN
C  Return to calling routine.
      END
C
C

ext: for


! { dg-do compile }
! { dg-options "-Winteger-division" }
program main
  print *,10**(-3) ! { dg-warning "Negative exponent of integer has zero result" }
end program main


! = = = = = = = = = = = = = = = =
!
!  Diffusion proccess for SSC
!
! = = = = = = = = = = = = = = = =

!    SUBROUTINE SSC_Diffusion_2D_Model(SSC1,SSC2)
!
!        IMPLICIT NONE
!
!        REAL(fp) , DIMENSION(:,:) , POINTER :: SSC1 , SSC2 , TEMP
!        INTEGER :: I , J , i1 , i2 , j1 , j2
!
!
!! Calculate local erosion
!
!        CALL LocalErosion
!
!!- - - - - - - - - - - - - -
!!    ADI method is used
!!- - - - - - - - - - - - - -
!
!! Y-Direcsion
!
!        LX   => DY
!        LY   => DX
!        LT   => DT
!        LXX = LX**2
!        LYY = LY**2
!
!        ILOW => ML
!        IUP  => MU
!        JLOW => NL
!        JUP  => NU
!
!        DO J = JLOW , JUP
!
!            j1 = J+1
!            j2 = J-1
!            IF(J==JLOW) j2=J
!            IF(J==JUP) j1=J
!
!! results of last step
!
!            VARS  => SSC1(:,j2)
!            VARP  => SSC1(:,J)
!            VARN  => SSC1(:,j1)
!
!! wet AND dry information
!
!            GIDS  => HNDD(:,j2)
!            GIDP  => HNDD(:,J)
!            GIDN  => HNDD(:,j1)
!
!! diffusion coefficient -Y
!
!            VSCP  => SSC_DIFFU_HWL_Y(:,J)
!
!! diffusion coefficient -X
!
!            VSCS  => SSC_DIFFU_HWL_X(:,j2)
!            VSCN  => SSC_DIFFU_HWL_X(:,J)
!
!! erosion term: local source term
!
!            CEP   => ERO(:,J)
!
!! local water Depth
!
!            HP    => HTDEP0(:,J)
!
!! result
!
!            FRUIT => SSC2(:,J)
!
!! ADI solver
!
!            CALL Tridiagonal_ADI('Y',J)
!
!        END DO
!
!! Exchange pointers
!
!        TEMP=>SSC2
!        SSC2=>SSC1
!        SSC1=>TEMP
!        NULLIFY(VARS,VARP,VARN,FRUIT,GIDS,GIDP,GIDN,VSCP,VSCS,VSCN,TEMP,CEP)
!        NULLIFY(LX,LY,LT,ILOW,IUP,JLOW,JUP)
!
! ! - - - - - - - - - - - - - - - - -
!
! ! X-Direcsion
!
!        LX   => DX
!        LY   => DY
!        LT   => DT
!        LXX = LX**2
!        LYY = LY**2
!
!        ILOW => NL
!        IUP  => NU
!        JLOW => ML
!        JUP  => MU
!
!        DO J = JLOW , JUP
!
!            j1 = J+1
!            j2 = J-1
!            IF(J==JLOW) j2=J
!            IF(J==JUP ) j1=J
!
!! results of last step
!
!            VARS  => SSC1(j2,:)
!            VARP  => SSC1(J,:)
!            VARN  => SSC1(j1,:)
!
!! wet AND dry information
!
!            GIDS  => HNDD(j2,:)
!            GIDP  => HNDD(J,:)
!            GIDN  => HNDD(j1,:)
!
!! diffusion coefficient -Y
!
!            VSCP  => SSC_DIFFU_HWL_X(J,:)
!
!! diffusion coefficient -X
!
!            VSCS  => SSC_DIFFU_HWL_Y(j2,:)
!            VSCN  => SSC_DIFFU_HWL_Y(J,:)
!
!! erosion term: local source term
!
!            CEP   => ERO(J,:)
!
!! local water Depth
!
!            HP    => HTDEP0(J,:)
!
!! result
!
!            FRUIT => SSC2(J,:)
!
!! ADI solver
!
!            CALL Tridiagonal_ADI('X',J)
!
!        END DO
!
!! Exchange pointers
!
!        TEMP=>SSC2
!        SSC2=>SSC1
!        SSC1=>TEMP
!        NULLIFY(VARS,VARP,VARN,FRUIT,GIDS,GIDP,GIDN,VSCP,VSCS,VSCN,TEMP,CEP)
!        NULLIFY(LX,LY,LT,ILOW,IUP,JLOW,JUP)
!
!    END SUBROUTINE SSC_Diffusion_2D_Model

! = = = = = = = = = = = = = = = =
!
!     SSC constant model
!
! = = = = = = = = = = = = = = = =

    !SUBROUTINE SSC_Constant_Model(SSC1)
    !    IMPLICIT NONE
    !
    !    REAL(fp) , DIMENSION(:,:) , POINTER :: SSC1
    !
    !    CALL LocalErosion
    !
    !END SUBROUTINE SSC_Constant_Model

    ! 2D routine

!    SUBROUTINE Bed_Evolution(SSC1)
!
!        IMPLICIT NONE
!
!        REAL(fp) , DIMENSION(:,:) , POINTER :: SSC1
!        INTEGER :: I , J
!
!! Erosion is positive and Deposition is negative
!
!        DO I = ML , MU
!            DO J = NL , NU
!                EVOSTEP(I,J) = ERO(I,J) - SSC1(I,J)*WS
!                Evolution_TDCYC(I,J) = Evolution_TDCYC(I,J) + EVOSTEP(I,J)
!            END DO
!        END DO
!
!    END SUBROUTINE Bed_Evolution

! = = = = = = = = = = = = = = = = = = = =
!
!   ADI ROUTINE: solve tridiagonal matrix
!
! = = = = = = = = = = = = = = = = = = = =

!    SUBROUTINE Tridiagonal_ADI(DIR,ROW)
!
!        IMPLICIT NONE
!
!        CHARACTER(LEN=1) :: DIR
!        INTEGER          :: ROW
!
!        REAL(fp) :: AMT(6,3001)
!        REAL(fp) :: DCE , DCW , DCN , DCS
!        INTEGER  :: I , IC , K
!        INTEGER  :: IDE, IDW, IDS, IDN
!        IC = 0
!
!        DO I = ILOW , IUP
!
!            IC = IC+1
!
!! Local point is dry
!
!            IF(GIDP(I)==1) THEN
!
!! set zero
!
!                amt(1,IC) = 0.D0
!                amt(2,IC) = 1.D0
!                amt(3,IC) = 0.D0
!                amt(4,IC) = 0.D0
!                CYCLE
!
!! Local point is wet
!
!            ELSEIF(GIDP(I)==0) THEN
!
!!- - - - - - - - - - - - - -
!! Set diffusion coefficients
!!- - - - - - - - - - - - - -
!
!!SOUTH
!
!                IF(ROW==JLOW) THEN
!                    IDS = 1
!                ELSE
!                    IDS = GIDS(I)
!                END IF
!
!                IF(IDS==0) THEN
!                    DCS = VSCS(I)  * LT / LYY * 0.5D0 / HP(I)
!                END IF
!
!!NORTH
!
!                IF(ROW==JUP) THEN
!                    IDN = 1
!                ELSE
!                    IDN = GIDN(I)
!                END IF
!
!                IF(IDN==0) THEN
!                    DCN = VSCN(I)  * LT / LYY * 0.5D0 / HP(I)
!                END IF
!
!!WEST
!
!                IF(I==ILOW) THEN
!                    IDW = 1
!                ELSE
!                    IDW = GIDP(I-1)
!                END IF
!
!                IF(IDW==0) THEN
!                    DCW = VSCP(I-1) * LT / LXX * 0.5D0 / HP(I)
!                END IF
!
!!EAST
!
!                IF(I==IUP) THEN
!                    IDE = 1
!                ELSE
!                    IDE = GIDP(I+1)
!                END IF
!
!                IF(IDE==0) THEN
!                    DCE = VSCP(I)  * LT / LXX * 0.5D0 / HP(I)
!                END IF
!
!!- - - - - - - - - - - - - - - -
!! Left hand side of the equation
!!- - - - - - - - - - - - - - - -
!
!                IF(IDE==1.AND.IDW==0) THEN
!
!                    amt(1,IC) = - 2.D0*DCW
!                    amt(2,IC) = + 2.D0*DCW + 1.D0
!                    amt(3,IC) = 0.D0
!
!                ELSEIF(IDW==1.AND.IDE==0) THEN
!
!                    amt(1,IC) =   0.D0
!                    amt(2,IC) = + 2.D0*DCE + 1.D0
!                    amt(3,IC) = - 2.D0*DCE
!
!                ELSEIF(IDW==1.AND.IDE==1) THEN
!
!                    amt(1,IC) =   0.D0
!                    amt(2,IC) =   1.D0
!                    amt(3,IC) =   0.D0
!
!                ELSEIF(IDW==0.AND.IDE==0)  THEN
!
!                    amt(1,IC) = - DCW
!                    amt(2,IC) = + DCW + DCE + 1.D0
!                    amt(3,IC) = - DCE
!
!                ELSE
!
!                    WRITE(*,*) ROW, GIDP(I-1) , GIDP(I+1)
!                    PAUSE 'tri_matrix_error_1'
!
!                END IF
!
!!- - - - - - - - - - - - - - - -
!! Left hand side of the equation
!!- - - - - - - - - - - - - - - -
!
!                IF(    IDN==1.AND.IDS==0) THEN
!
!                    amt(4,IC) = VARP(I) + 2.D0*DCS*(VARS(I)-VARP(I))
!
!                ELSEIF(IDS==1.AND.IDN==0) THEN
!
!                    amt(4,IC) = VARP(I) + 2.D0*DCN*(VARN(I)-VARP(I))
!
!                ELSEIF(IDS==1.AND.IDN==1) THEN
!
!                    amt(4,IC) = VARP(I)
!
!                ELSEIF(IDS==0.AND.IDN==0)  THEN
!
!                    amt(4,IC) = VARP(I) + DCN*(VARN(I)-VARP(I)) + DCS*(VARS(I)-VARP(I))
!
!                ELSE
!
!                    WRITE(*,*) ROW, IDS , IDN
!                    PAUSE 'tri_matrix_error_2'
!
!                END IF
!
!        ! sourse term
!
!                amt(4,IC) = amt(4,IC) + CEP(I)/HP(I)*0.5D0*LT
!                amt(2,IC) = amt(2,IC) + WS/HP(I)*0.5D0*LT
!
!            END IF
!
!        END DO
!
!    !    OPEN(5,file='amt.txt')
!    !    do I = 1, IC
!    !        WRITE(5,'(10e25.10)') (amt(K,I),K=1,4)
!    !    end do
!    !    CLOSE(5)
!
!        CALL solve3_LU(IC)
!
!        fruit(:) = amt(4,:)
!
!        !IF(pros=='sedi') THEN
!        !IF(dir=='y') THEN
!        !OPEN(5,file='amt2.txt')
!        !do I = 1, IC
!        !    WRITE(5,'(10e25.10)') (amt(K,I),K=1,4)
!        !end do
!        !CLOSE(5)
!        !pause 'xg'
!        !end IF
!
!    END SUBROUTINE Tridiagonal_ADI

!- - - - - - - - - - - - - -
!
!   Bed evolution per step
!
!- - - - - - - - - - - - - -

! 1D routine


module module_mint
  use, intrinsic :: iso_fortran_env
  implicit none
  integer(int64), parameter :: modulus = 998244353_int64
  type mint
    integer(int64), private :: val = 0_int64
  contains
    procedure, public :: get => get
    procedure, public :: inv => inv
  end type
  interface mint
    module procedure :: new_mint_int64, new_mint_int32, new_mint_int16, new_mint_int8
  end interface
  interface assignment(=)
    module procedure :: asgn_int64, asgn_int32, asgn_int16, asgn_int8
  end interface
  interface operator(+)
    module procedure :: add, add_int64, add_int64_rev, add_int32, add_int32_rev, add_int16, add_int16_rev, add_int8, add_int8_rev
  end interface
  interface operator(-)
    module procedure :: sub, sub_int64, sub_int64_rev, sub_int32, sub_int32_rev, sub_int16, sub_int16_rev, sub_int8, sub_int8_rev
  end interface
  interface operator(*)
    module procedure :: mul, mul_int64, mul_int64_rev, mul_int32, mul_int32_rev, mul_int16, mul_int16_rev, mul_int8, mul_int8_rev
  end interface
  interface operator(/)
    module procedure :: div, div_int64, div_int64_rev, div_int32, div_int32_rev, div_int16, div_int16_rev, div_int8, div_int8_rev
  end interface
  interface operator(**)
    module procedure :: pow_int64, pow_int32, pow_int16, pow_int8
  end interface
  interface dot_product
    module procedure :: dot_prod
  end interface dot_product
  interface matmul
    module procedure :: matmul1, matmul2, matmul3
  end interface matmul
  public :: mint, assignment(=), operator(+), operator(-), operator(*), operator(/), operator(**), dot_product, matmul
  private :: new_mint_int64, new_mint_int32, new_mint_int16, new_mint_int8
  private :: asgn_int64, asgn_int32, asgn_int16, asgn_int8
  private :: add, add_int64, add_int64_rev, add_int32, add_int32_rev, add_int16, add_int16_rev, add_int8, add_int8_rev
  private :: sub, sub_int64, sub_int64_rev, sub_int32, sub_int32_rev, sub_int16, sub_int16_rev, sub_int8, sub_int8_rev
  private :: mul, mul_int64, mul_int64_rev, mul_int32, mul_int32_rev, mul_int16, mul_int16_rev, mul_int8, mul_int8_rev
  private :: div, div_int64, div_int64_rev, div_int32, div_int32_rev, div_int16, div_int16_rev, div_int8, div_int8_rev
  private :: pow_int64, pow_int32, pow_int16, pow_int8
  private :: dot_prod, matmul1, matmul2, matmul3
  private :: get, inv
contains
  pure elemental integer(int64) function get(this) result(res)
    class(mint), intent(in) :: this
    res = this%val
  end
  pure elemental subroutine asgn_int64(this, other)
    type(mint), intent(inout) :: this
    integer(int64), intent(in) :: other
    this%val = merge(other, modulo(other, modulus), 0_int64 <= other .and. other < modulus)
  end
  pure elemental subroutine asgn_int32(this, other)
    type(mint), intent(inout) :: this
    integer(int32), intent(in) :: other
    call asgn_int64(this, int(other, int64))
  end
  pure elemental subroutine asgn_int16(this, other)
    type(mint), intent(inout) :: this
    integer(int16), intent(in) :: other
    call asgn_int64(this, int(other, int64))
  end
  pure elemental subroutine asgn_int8(this, other)
    type(mint), intent(inout) :: this
    integer(int8), intent(in) :: other
    call asgn_int64(this, int(other, int64))
  end
  pure elemental type(mint) function add(this, other) result(res)
    type(mint), intent(in) :: this
    type(mint), intent(in) :: other
    integer(int64) :: tmp
    tmp = this%val + other%val
    res%val = merge(tmp, tmp - modulus, 0_int64 <= tmp .and. tmp < modulus)
  end
  pure elemental type(mint) function sub(this, other) result(res)
    type(mint), intent(in) :: this
    type(mint), intent(in) :: other
    integer(int64) :: tmp
    tmp = this%val - other%val
    res%val = merge(tmp, tmp + modulus, 0_int64 <= tmp .and. tmp < modulus)
  end
  pure elemental type(mint) function mul(this, other) result(res)
    type(mint), intent(in) :: this
    type(mint), intent(in) :: other
    integer(int64) :: tmp
    tmp = this%val * other%val
    res%val = merge(tmp, modulo(tmp, modulus), 0_int64 <= tmp .and. tmp < modulus)
  end
  pure elemental type(mint) function div(this, other) result(res)
    type(mint), intent(in) :: this
    type(mint), intent(in) :: other
    res = mul(this, inv(other))
  end
  pure elemental type(mint) function inv(this) result(res)
    class(mint), intent(in) :: this
    integer(int64) :: a, b, c, n
    integer(int64) :: s, t
    a = this%val
    b = modulus
    c = 0_int64
    n = 1_int64
    do while (b /= 0_int64)
      s = mod(a, b)
      t = n - a / b * c
      a = b
      b = s
      n = c
      c = t
    end do
    res%val = modulo(n, modulus)
  end
  pure elemental type(mint) function pow_int64(this, other) result(res)
    type(mint), intent(in) :: this
    integer(int64), intent(in) :: other
    integer(int64) :: a, k, ak
    a = get(this)
    k = other
    ak = 1_int64
    do while (k > 0_int64)
      if (btest(k, 0)) then
        ak = ak * a
        if (ak >= modulus) ak = mod(ak, modulus)
      end if
      a = a * a
      if (a >= modulus) a = mod(a, modulus)
      k = rshift(k, 1)
    end do
    res%val = modulo(ak, modulus)
  end
  pure elemental type(mint) function pow_int32(this, other) result(res)
    type(mint), intent(in) :: this
    integer(int32), intent(in) :: other
    res = pow_int64(this, int(other, int64))
  end
  pure elemental type(mint) function pow_int16(this, other) result(res)
    type(mint), intent(in) :: this
    integer(int16), intent(in) :: other
    res = pow_int64(this, int(other, int64))
  end
  pure elemental type(mint) function pow_int8(this, other) result(res)
    type(mint), intent(in) :: this
    integer(int8), intent(in) :: other
    res = pow_int64(this, int(other, int64))
  end
  pure type(mint) function dot_prod(x, y) result(res)
    type(mint), intent(in) :: x(:), y(:)
    integer(int32) :: l, k
    l = size(x, 1)
    do k = 1, l
      res%val = mod(res%val + x(k)%val * y(k)%val, modulus)
    end do
  end
  pure function matmul1(x, y) result(res)
    type(mint), intent(in) :: x(:, :), y(:)
    type(mint) :: res(size(x, 1))
    integer(int32) :: n, l, i, k
    integer(int64) :: tmp
    n = size(x, 1)
    l = size(x, 2)
    do k = 1, l
      tmp = y(k)%val
      do i = 1, n
        res(i)%val = mod(res(i)%val + x(i, k)%val * tmp, modulus)
      end do
    end do
  end
  pure function matmul2(x, y) result(res)
    type(mint), intent(in) :: x(:), y(:, :)
    type(mint) :: res(size(y, 2))
    integer(int32) :: m, l, j, k
    m = size(y, 2)
    l = size(y, 1)
    do j = 1, m
      do k = 1, l
        res(j)%val = mod(res(j)%val + x(k)%val * y(k, j)%val, modulus)
      end do
    end do
  end
  pure function matmul3(x, y) result(res)
    type(mint), intent(in) :: x(:, :), y(:, :)
    type(mint) :: res(size(x, 1), size(y, 2))
    integer(int32) :: n, m, l, i, j, k
    integer(int64) :: tmp
    n = size(x, 1)
    m = size(y, 2)
    l = size(x, 2)
    do j = 1, m
      do k = 1, l
        tmp = y(k, j)%val
        do i = 1, n
          res(i, j)%val = mod(res(i, j)%val + x(i, k)%val * tmp, modulus)
        end do
      end do
    end do
  end
  pure function matpow(x, n_) result(res)
    type(mint), intent(in) :: x(:, :)
    integer(int64), intent(in) :: n_
    type(mint) :: res(size(x, 1), size(x, 2)), tmp(size(x, 1), size(x, 2))
    integer(int32) :: i, n
    do i = 1, size(x, 1)
      res(i, i)%val = 1_int64
    end do
    tmp = x
    n = n_
    do while (n > 0_int64)
      if (btest(n, 0)) res = matmul3(res, tmp)
      tmp = matmul3(tmp, tmp)
      n = rshift(n, 1)
    end do
  end
  pure elemental type(mint) function new_mint_int64(val) result(res)
    integer(int64), intent(in) :: val
    call asgn_int64(res, val)
  end
  pure elemental type(mint) function add_int64(this, other) result(res)
    type(mint), intent(in) :: this
    integer(int64), intent(in) :: other
    res = add(this, new_mint_int64(other))
  end
  pure elemental type(mint) function sub_int64(this, other) result(res)
    type(mint), intent(in) :: this
    integer(int64), intent(in) :: other
    res = sub(this, new_mint_int64(other))
  end
  pure elemental type(mint) function mul_int64(this, other) result(res)
    type(mint), intent(in) :: this
    integer(int64), intent(in) :: other
    res = mul(this, new_mint_int64(other))
  end
  pure elemental type(mint) function div_int64(this, other) result(res)
    type(mint), intent(in) :: this
    integer(int64), intent(in) :: other
    res = div(this, new_mint_int64(other))
  end
  pure elemental type(mint) function add_int64_rev(other, this) result(res)
    type(mint), intent(in) :: this
    integer(int64), intent(in) :: other
    res = add(new_mint_int64(other), this)
  end
  pure elemental type(mint) function sub_int64_rev(other, this) result(res)
    type(mint), intent(in) :: this
    integer(int64), intent(in) :: other
    res = sub(new_mint_int64(other), this)
  end
  pure elemental type(mint) function mul_int64_rev(other, this) result(res)
    type(mint), intent(in) :: this
    integer(int64), intent(in) :: other
    res = mul(new_mint_int64(other), this)
  end
  pure elemental type(mint) function div_int64_rev(other, this) result(res)
    type(mint), intent(in) :: this
    integer(int64), intent(in) :: other
    res = div(new_mint_int64(other), this)
  end
  pure elemental type(mint) function new_mint_int32(val) result(res)
    integer(int32), intent(in) :: val
    call asgn_int32(res, val)
  end
  pure elemental type(mint) function add_int32(this, other) result(res)
    type(mint), intent(in) :: this
    integer(int32), intent(in) :: other
    res = add(this, new_mint_int32(other))
  end
  pure elemental type(mint) function sub_int32(this, other) result(res)
    type(mint), intent(in) :: this
    integer(int32), intent(in) :: other
    res = sub(this, new_mint_int32(other))
  end
  pure elemental type(mint) function mul_int32(this, other) result(res)
    type(mint), intent(in) :: this
    integer(int32), intent(in) :: other
    res = mul(this, new_mint_int32(other))
  end
  pure elemental type(mint) function div_int32(this, other) result(res)
    type(mint), intent(in) :: this
    integer(int32), intent(in) :: other
    res = div(this, new_mint_int32(other))
  end
  pure elemental type(mint) function add_int32_rev(other, this) result(res)
    type(mint), intent(in) :: this
    integer(int32), intent(in) :: other
    res = add(new_mint_int32(other), this)
  end
  pure elemental type(mint) function sub_int32_rev(other, this) result(res)
    type(mint), intent(in) :: this
    integer(int32), intent(in) :: other
    res = sub(new_mint_int32(other), this)
  end
  pure elemental type(mint) function mul_int32_rev(other, this) result(res)
    type(mint), intent(in) :: this
    integer(int32), intent(in) :: other
    res = mul(new_mint_int32(other), this)
  end
  pure elemental type(mint) function div_int32_rev(other, this) result(res)
    type(mint), intent(in) :: this
    integer(int32), intent(in) :: other
    res = div(new_mint_int32(other), this)
  end
  pure elemental type(mint) function new_mint_int16(val) result(res)
    integer(int16), intent(in) :: val
    call asgn_int16(res, val)
  end
  pure elemental type(mint) function add_int16(this, other) result(res)
    type(mint), intent(in) :: this
    integer(int16), intent(in) :: other
    res = add(this, new_mint_int16(other))
  end
  pure elemental type(mint) function sub_int16(this, other) result(res)
    type(mint), intent(in) :: this
    integer(int16), intent(in) :: other
    res = sub(this, new_mint_int16(other))
  end
  pure elemental type(mint) function mul_int16(this, other) result(res)
    type(mint), intent(in) :: this
    integer(int16), intent(in) :: other
    res = mul(this, new_mint_int16(other))
  end
  pure elemental type(mint) function div_int16(this, other) result(res)
    type(mint), intent(in) :: this
    integer(int16), intent(in) :: other
    res = div(this, new_mint_int16(other))
  end
  pure elemental type(mint) function add_int16_rev(other, this) result(res)
    type(mint), intent(in) :: this
    integer(int16), intent(in) :: other
    res = add(new_mint_int16(other), this)
  end
  pure elemental type(mint) function sub_int16_rev(other, this) result(res)
    type(mint), intent(in) :: this
    integer(int16), intent(in) :: other
    res = sub(new_mint_int16(other), this)
  end
  pure elemental type(mint) function mul_int16_rev(other, this) result(res)
    type(mint), intent(in) :: this
    integer(int16), intent(in) :: other
    res = mul(new_mint_int16(other), this)
  end
  pure elemental type(mint) function div_int16_rev(other, this) result(res)
    type(mint), intent(in) :: this
    integer(int16), intent(in) :: other
    res = div(new_mint_int16(other), this)
  end
  pure elemental type(mint) function new_mint_int8(val) result(res)
    integer(int8), intent(in) :: val
    call asgn_int8(res, val)
  end
  pure elemental type(mint) function add_int8(this, other) result(res)
    type(mint), intent(in) :: this
    integer(int8), intent(in) :: other
    res = add(this, new_mint_int8(other))
  end
  pure elemental type(mint) function sub_int8(this, other) result(res)
    type(mint), intent(in) :: this
    integer(int8), intent(in) :: other
    res = sub(this, new_mint_int8(other))
  end
  pure elemental type(mint) function mul_int8(this, other) result(res)
    type(mint), intent(in) :: this
    integer(int8), intent(in) :: other
    res = mul(this, new_mint_int8(other))
  end
  pure elemental type(mint) function div_int8(this, other) result(res)
    type(mint), intent(in) :: this
    integer(int8), intent(in) :: other
    res = div(this, new_mint_int8(other))
  end
  pure elemental type(mint) function add_int8_rev(other, this) result(res)
    type(mint), intent(in) :: this
    integer(int8), intent(in) :: other
    res = add(new_mint_int8(other), this)
  end
  pure elemental type(mint) function sub_int8_rev(other, this) result(res)
    type(mint), intent(in) :: this
    integer(int8), intent(in) :: other
    res = sub(new_mint_int8(other), this)
  end
  pure elemental type(mint) function mul_int8_rev(other, this) result(res)
    type(mint), intent(in) :: this
    integer(int8), intent(in) :: other
    res = mul(new_mint_int8(other), this)
  end
  pure elemental type(mint) function div_int8_rev(other, this) result(res)
    type(mint), intent(in) :: this
    integer(int8), intent(in) :: other
    res = div(new_mint_int8(other), this)
  end
end module


PROGRAM exercise13_4
IMPLICIT NONE
REAL :: x = 12., y = -3., result
INTEGER :: i = 6, j = 4
WRITE (*,100) ' Before call: x, y, i, j = ', x, y, i, j
100 FORMAT (A,2F6.1,2I6)
result = exec(y,i)
WRITE (*,*) 'The result is ', result
WRITE (*,100) ' After call:  x, y, i, j = ', x, y, i, j

CONTAINS
   REAL FUNCTION exec(x,i)
   REAL, INTENT(IN) :: x
   INTEGER, INTENT(IN) :: i
   WRITE (*,100) ' In exec:     x, y, i, j = ', x, y, i, j
   100 FORMAT (A,2F6.1,2I6)
   exec = ( x + y ) / REAL ( i + j )
   j = i
   END FUNCTION exec
END PROGRAM exercise13_4


text:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!   EVB-QMDFF - RPMD molecular dynamics and rate constant calculations on
!               black-box generated potential energy surfaces
!
!   Copyright (c) 2021 by Julien Steffen (steffen@pctc.uni-kiel.de)
!                         Stefan Grimme (grimme@thch.uni-bonn.de) (QMDFF code)
!
!   Permission is hereby granted, free of charge, to any person obtaining a
!   copy of this software and associated documentation files (the "Software"),
!   to deal in the Software without restriction, including without limitation
!   the rights to use, copy, modify, merge, publish, distribute, sublicense,
!   and/or sell copies of the Software, and to permit persons to whom the
!   Software is furnished to do so, subject to the following conditions:
!
!   The above copyright notice and this permission notice shall be included in
!   all copies or substantial portions of the Software.
!
!   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
!   IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
!   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
!   THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
!   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
!   FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
!   DEALINGS IN THE SOFTWARE.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!
!     ##################################################################
!     ##                                                              ##
!     ##  module evb_mod  --  parameters for evb_qmdff-calculations   ##
!     ##                                                              ##
!     ##################################################################
!
!
!     "evb_mod" contains all important global informations that are
!     used in the evb-qmdff-parts: qmdffgen,evbopt,dynamic etc.
!

module evb_mod
implicit none
save

integer::nqmdff  !the number of used QMDFF´s
integer n_one,natoms
real(kind=8),allocatable :: xyz (:,:)
real(kind=8),allocatable :: q   (:)
real(kind=8),allocatable :: g_one (:,:)
real(kind=8),allocatable :: c6xy(:,:)
real(kind=8),allocatable :: cn  (:)
integer,allocatable :: at (:)
integer,allocatable :: at2(:)
integer,allocatable :: imass(:)
!integer,allocatable :: molnum(:)  ! for solvent box simulations
!real(kind=8)::dens  ! density (i.e. number of molecules in box)

integer n_two
real(kind=8) ,allocatable :: xyz_two (:,:)
real(kind=8) ,allocatable :: q_two(:)
real(kind=8) ,allocatable :: g_two(:,:)
real(kind=8) ,allocatable :: c6xy_two(:,:)
real(kind=8) ,allocatable :: cn_two(:)
integer,allocatable :: at_two(:)
integer,allocatable :: imass_two(:)

integer n_three
real(kind=8) ,allocatable :: xyz_three (:,:)
real(kind=8) ,allocatable :: q_three(:)
real(kind=8) ,allocatable :: g_three(:,:)
real(kind=8) ,allocatable :: c6xy_three(:,:)
real(kind=8) ,allocatable :: cn_three(:)
integer,allocatable :: at_three(:)
integer,allocatable :: imass_three(:)

real(kind=8)::offa,offb,offc,offd,offe,offf,offg,offh,offi,offj,&
        &offk,offl,offm,offn
! The shifted single-point energy of the QMDFFs
real(kind=8)::E_zero1,E_zero2,E_zero3

real(kind=8),dimension(94,94)::r0ab,zab,sr42,r094_mod
real(kind=8),dimension(12) :: start_params
real(kind=8),dimension(14) :: evb3_params
!integer::mode,off_choice
character(len=40)::mode1,off_basis,calc
character(len=40)::orca_inp
character(len=:),allocatable::orca_egrad
! the simple dE coupling term
logical::evb_de
! decides, if evb_qmdff or orca_calculation
logical::orca_calculation
! shall plumed be used or not!
logical::use_plumed
! if you want to control that the energy is conserved
logical::energycon
! separate printout of covalent and noncovalent energies for a single QMDFF
logical::energysplit
real(kind=8)::e_cov_split,e_noncov_split
! maximal iterations for levmarq
integer::optsteps
! start-value for levmarq parameter between newton and CG and the alternation
real(kind=8)::lm_par
real(kind=8)::lm_par_change
! threshold-parameter for convergence in levmarq
real(kind=8)::lm_threshold
! numercial differentation-step in levmarq for 3x3-EVB
real(kind=8)::diff_step
! for the evb_dq-gradient
real(kind=8),dimension(:),allocatable::ts_coordinates,ts_coordinates2
logical::use_dq,evb_dq
! for mopac reference-calculations: number of commandlines after coordsection
integer::mopac_maxline
! numerical gradient for 3x3-evb for comparison, sets 1-3 elements not zero
logical num_grad, full
real(kind=8)::num_grad_step   ! stepsize for numerical gradients
!  parameters for DG-EVB with more than one reference point
logical::dg_evb
! Plumed Inputfile
character(len=80)::plumed_input
integer::dg_evb_points ! number of reference points for DG-EVB calculation
integer::nat6 ! NUMBER OF INTERNAL COORDINATES
integer::num_coord ! number of internal coordinates defined by user
real(kind=8),dimension(:),allocatable::all_ens,allf1_ens,allf2_ens,f_vec,b_vec
real(kind=8),dimension(:,:),allocatable::all_xyz,d_mat,all_int,point_int
real(kind=8),dimension(:,:),allocatable::all_grad,allf1_grad,allf2_grad
real(kind=8),dimension(:,:,:),allocatable::all_hess,allf1_hess,allf2_hess
logical::dg_evb_opt_multi,diff_ana ! if the alpha values should be optimized separately or at once
integer::maxstart !maximal number of local starts for optimization
integer::dg_mode       ! for DG-EVB (1,2 or 3)
integer::ref_input_unit  ! unit for the ref.input file in DG-EVB calculations
real(kind=8),dimension(:),allocatable::alph_opt,b_opt
! start interval for random number
real(kind=8)::lower_bond,upper_bond
! if the ordering of internal coordinates should be read in manually
logical::read_coord
! if the internal coordinates shall be defined as distance matrix!
logical::dist_matrix
! global array of internal coordinate definitions (D:1; Z:4)
integer,dimension(:,:),allocatable::coord_nr,coord_def
integer,dimension(:,:),allocatable::bonds,angles,ooplanes,diheds
! global integers for numbers of single coordinates (bonds,angles,ooplane,diheds)
integer::nbonds,nangles,noops,ndiheds
! for normal mode analysis
logical::print_nm
! for geometry optimization
integer::geomax
! global parameter for additional info output during programs
logical::more_info
! if only a fragment of the structure shall be considered for calculation
logical::calc_frag
integer,dimension(:),allocatable::indi
integer,dimension(:),allocatable::indi3
integer::nats
! for IRC optimizations
character(len=50):: irc_method
real(kind=8)::irc_steplen
real(kind=8)::irc_eulerlen
integer::irc_maxstep
! test case: mueller brown surface
! name of the DG-EVB reference file
character(len=80)::dg_ref_file
! if the DG-EVB coefficient matrix D shall be calculated numerically
logical::dg_mat_num
! minimum value below that distant Gaussians are neglected
real(kind=8)::g_thres
! if for DG-EVB mode 3 two alphas per reference point shall be used
logical::double_alpha
integer::add_alph
! for the TREQ / RP-EVB coupling method
logical::treq
integer::rp_evb_points  ! number of reference points with energies, gradients and hessians
integer::rp_irc_points  ! number of structures with energies along the IRC
integer::rp_spl2_dim   ! total number of reference informations per gradient/hessian point
real(kind=8)::lambda  ! the inverse average disctance between two reference frames
real(kind=8)::pareta  ! shape parameter for the damping of V12
real(kind=8)::interp_tol  ! how small the interval for s-determin. shall become
real(kind=8)::rp_ana_step  ! for analytical gradient: step of internal coordinate elongation
real(kind=8),allocatable::rp_spl_ref(:,:) ! reference information for spline interpolation
real(kind=8),allocatable::rp_spl_d2(:,:) ! second derivatives for spline interpolation
real(kind=8),allocatable::rp_spl_s(:)  ! the parameter reference values for spline interp.
real(kind=8),allocatable::rp_spl_ref2(:,:) ! same as above, but for gradients/hessians
real(kind=8),allocatable::rp_spl_d22(:,:) ! same as above, but for gradients/hessians
real(kind=8),allocatable::rp_spl_s2(:) ! same as above, but for gradients/hessians
real(kind=8),allocatable::rp_point_s(:) ! s values of the points with gradients and hessians
real(kind=8),allocatable::rp_point_v12(:)  ! coupling strengths on these points
real(kind=8),dimension(:,:),allocatable::dv12 ! the coupling gradients at the TREQ points
real(kind=8),dimension(:,:,:),allocatable::d2v12  ! the coupling hessians at the TREQ points
integer::qmdff_order    ! if QMDFF1=educt, QMDFF2=product or the other way
real(kind=8)::s_transfer  ! the actual value of the reaction coordinate s
real(kind=8)::s_ts    ! position of the TS along the path (in terms of the s-parameter)
real(kind=8)::s_tot    ! the total length of the IRC for rescaling
real(kind=8)::par_epsi  ! slope of the correction potential parabola
real(kind=8)::pre_exp   ! exponential prefactor of the gaussian damping for RP-EVB
real(kind=8)::corr_max   ! maximal value of the parabola prefactor for the QMDFF corrections
real(kind=8)::deltp   !  half size of the whole interval for which QMDFF shall be corrected
real(kind=8)::ddeltp  !  shift of the correction interval borders for the QMDFF hessians
logical::s_bord_man  ! if the borders along the IRC shall be read in manually
real(kind=8)::rp_mid_tot,rp_mid_trans   ! NEW 04.09.2018: only energy region near the TS
real(kind=8)::trans_r_lo,trans_r_hi  ! transition borders right to the ts (s values)
real(kind=8)::trans_l_lo,trans_l_hi   ! transition borders left to the ts (s values)
real(kind=8),dimension(:,:),allocatable::xyz_r_hi,xyz_l_lo  ! xyz structures of RP-EVB borders
real(kind=8),dimension(:),allocatable::int_r_hi,int_l_lo  ! int proj structures of RP-EVB borders
logical::shift_man   !  turn off automatical optimization of QMDFF minimum shift energies
logical::int_grad_plot  ! option to activate printing of internal gradient components to file..
logical::int_coord_plot  ! option to activate printing of internal coordinate components
logical::num_wilson  ! option, to numerically calculate the wilson b matrix for coord. transf.
logical::use_internals ! option to use internal coordinates for geometry optimizations
logical::spline_all  ! option to activate full IRC scan for Xi value for RP-EVB points
real(kind=8)::path_dist_limit  ! if the z-value is larger than this limit, QMDFF1/2 will be assumed
logical::no_evb  ! if the RP-EVB part shall be deactivated and only direct interpolation/QMDFF be used
real(kind=8)::iq_l_lo,iq_l_hi,iq_r_lo,iq_r_hi  ! borders between direct interpolation and QMDFFs
! for Mechanochemistry
logical::add_force
integer::force1_at,force2_at
real(kind=8)::force1_k,force2_k
real(kind=8)::force1_v(3),force2_v(3)
logical::afm_run
integer::afm_fix_at,afm_move_at
real(kind=8)::afm_move_dist
real(kind=8)::afm_move_v(3)
real(kind=8)::afm_move_first(3)
real(kind=8)::afm_k
integer::afm_steps
! for the new RPMD program
integer::nbeads   ! the number of ring polymer beads to be sampled
integer::npaths  ! number of equivalent reaction paths
real(kind=8),dimension(:,:,:),allocatable::p_i,q_i  ! positions and momenta of the system
logical::use_rpmd  ! if the rpmd.x program is used, allow also 1 internal coordinate..
integer,allocatable::rpmd_atoms(:) ! if not all atoms shall be treated with RPMD, define the active ones
real(kind=8),dimension(:,:),allocatable::ts_ref,ed_ref
real(kind=8),dimension(:),allocatable::int_ed,int_ts,xi_scale ! scale int. coordinates to 0..1
real(kind=8),dimension(:),allocatable::form_ref,break_ref  ! reference values for the TS bonds
real(kind=8),dimension(:),allocatable::form_ed,break_ed ! reference values of educts bonds (unimol.)
real(kind=8)::k_force ! umbrella force constant
real(kind=8),allocatable::struc_equi(:,:,:)
real(kind=8)::umbr_lo,umbr_hi ! bonds for umbrella samplings
real(kind=8)::umbr_dist ! distance between single umbrella samplings
character(len=20)::umbr_type  ! type of coordinate system to be used for umbrella sampling
integer::form_num,break_num  ! number of forming and breaking bonds for the reac.type
integer,dimension(:,:),allocatable::bond_form,bond_break  ! forming and breaking bonds
real(kind=8),allocatable::r_refs(:,:)  ! the needed educt-educt distances for the pre-equilibration
integer::shift_atom,shift_coord  ! for the ATOM_SHIFT coordinate
real(kind=8)::shift_lo,shift_hi  ! for the ATOM_SHIFT coordinate
logical::fix_atoms  ! if some atoms shall be fixed in dynamics calculations
integer,allocatable::fix_list(:) ! file with list of fixed atoms
integer::fix_num  ! number of fixed atoms
integer::umbr_traj ! number of umbrella trajectories per Xi value
real(kind=8)::t_avg ! TEST for temperature equilibrium
integer::n_samplings,n_over  ! number of sampling windows
logical::writestat ! decides if MD statistics shall be written out
character(len=20)::pmf_method  ! if WHAM or umbrella integration shall be used!
real(kind=8),allocatable::average(:),variance(:)   ! for umbrella integration!
integer::sum_eds   ! total number of educt molcules (1 up to 4)
integer,allocatable::at_ed(:,:)  ! atom numbers of all educts (up to four for ADDITION4)
real(kind=8)::R_inf   ! reference distance of both/all educts
integer,allocatable::n_ed(:) ! number of atoms in the educts (up to four for ADDITION4)
real(kind=8),allocatable::mass_ed(:)  ! masses of the educts (up to four for ADDITION4)
real(kind=8)::beta  ! inverse temperature
real(kind=8)::andersen_time ! application interval for Andersen thermostat
character(len=40)::thermo ! which thermostat shall be used (character)?
!integer::thermostat   ! switch for thermostat: 0: Anderson, 1: GLE
integer::andersen_step  ! after how many steps the Andersen thermostat shall be applied
integer::recr_equi,child_tot,child_interv,child_point,child_evol ! recrossing calc. parameters
integer::child_times
logical::recross_calc  ! signalize of a recrossing calculation is done currently
real(kind=8),allocatable::gle_mat_A(:,:),gle_mat_C(:,:)  ! the matrices used for the GLE thermostat
integer::gle_ns  ! number of dimensions for GLE thermostat
real(kind=8),allocatable::gle_S(:,:),gle_T(:,:),gle_p(:,:,:,:)  ! GLE thermostat arrays
real(kind=8),allocatable::gle_np(:,:,:,:) !  GLE thermostat arrays
logical::gen_test  ! perform only initial start structure generations
integer::gen_pr_frac,gen_pr_act   ! only each N'th structure, energy etc shall be written out
real(kind=8)::gen_energies  ! total sum of energies during structure generation
! For box simulations without periodic boundary conditions
logical::box_nonp
real(kind=8)::box_center(3)
real(kind=8)::box_size
integer::box_fixed(1000)
integer::num_fixed
! For constant pressure npt simulations with a barostat
logical::npt
real(kind=8)::pressure
logical::verbose  ! activates verbose output of several programs

!   If the starting velocity shall be read in
logical::read_vel
real(kind=8),allocatable::vel_start(:,:)

!   For QMDFF parameter optimization
real(kind=8),allocatable::xyz_ref_lm(:,:,:)

!   the averaged kinetic energy
real(kind=8)::ekin_avg
real(kind=8)::ekin2_avg
logical::calc_ekin
integer::ekin_num
integer::ekin_atoms(100000)

!    change mass of one element/isotope
logical::change_mass
real(kind=8)::newmass
character(len=2)::elem_mass
! for orca calculations (ab initio MD)
logical::orca
character(len=100)::orca_com
! for the water SPC model
logical::water_spc
real(kind=8),allocatable::water_pars(:) ! all needed parameters
integer::nwater  ! number of water molecules
integer,allocatable::water_act(:)  ! to which water molecule the atoms belong
!  for boxes with hard walls
logical::box_walls
real(kind=8)::walldim_x,walldim_y,walldim_z
!   For the Nose-Hoover/Berendsen barostat
integer::barostat   ! switch: 1: Nose-Hoover, 2: Berendsen
real(kind=8)::vel_baro
real(kind=8)::temp_test  ! for andersen testing..
real(kind=8)::press_act   ! the actual pressure to be calculated
real(kind=8)::mass_tot   ! for density calculations (mass of the system)
! for acceleration of calculations: store previous infos
real(kind=8),allocatable::inter_old(:,:),change_tol
integer,allocatable::i_best_old(:)
integer::irc_local  ! readin: how many structures left/right shall be investigated
! for geometry optimization (minima, TS)
logical::newton_raphson  ! if the simpler Newton-Raphson method shall be used instead of P-RFO
integer::maxiter  ! maximal number of geoopt steps
real(kind=8)::stepmax ! maximal step size
real(kind=8)::ethr,gthr,dthr,gmaxthr,dmaxthr  ! convergence criteria (geoopt)
real(kind=8)::irc_ethr,irc_gthr  ! convergence criteria (IRC)
integer::rp_evb_mode  ! if only energies and gradients shall be used for RP-EVB coupling
! for error handling
logical::act_check    ! if error handling is activated at all..
logical::recross_check  ! if checking is activated for recrossing
logical::loose_check   ! samplings are not canceled after single error
integer::err_max,err_act_max
real(kind=8)::energy_tol  ! tolerance for energy in relation to TS energy per bead
!    TEST TEST TEST
real(kind=8),allocatable::error_vec_ens(:)  ! vector with energy parts
real(kind=8),allocatable::error_vec_v12(:)  ! vector with coupling parts
real(kind=8)::xi_test  ! for manual shutdown
real(kind=8)::no_evb_xi  ! maximum xi value for pure QMDFF for no-evb
! for time measuring
real(kind=8)::time1,time2
real(kind=8)::time(10),duration
end module evb_mod

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: chrinide/EVB-QMDFF

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: chrinide/EVB-QMDFF

max_forks_repo_licenses: MIT Unlicense


text:
      SUBROUTINE CALPSER (ISTL_)
C
C CHANGE RECORD
C ** SUBROUTINE CALPSER UPDATES TIME VARIABLE SURFACE ELEVATION
C ** BOUNDARY CONDITIONS
C
      USE GLOBAL
      PSERT(0)=0.
      DO NS=1,NPSER
        IF(ISDYNSTP.EQ.0)THEN
          TIME=DT*FLOAT(N)/TCPSER(NS)+TBEGIN*(TCON/TCPSER(NS))
        ELSE
          TIME=TIMESEC/TCPSER(NS)
        ENDIF
        M1=MPTLAST(NS)
  100   CONTINUE
        M2=M1+1
        IF(TIME.GT.TPSER(M2,NS))THEN
          M1=M2
          GOTO 100
        ELSE
          MPTLAST(NS)=M1
        ENDIF
        TDIFF=TPSER(M2,NS)-TPSER(M1,NS)
        WTM1=(TPSER(M2,NS)-TIME)/TDIFF
        WTM2=(TIME-TPSER(M1,NS))/TDIFF
        PSERT(NS)=WTM1*PSER(M1,NS)+WTM2*PSER(M2,NS)
      ENDDO
 6000 FORMAT('N, PSERT = ',I6,4X,F12.4)
      RETURN
      END

ext: f


text:
C
C $Id: mapstr.f,v 1.20 2008-09-18 00:42:16 kennison Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE MAPSTR (WHCH,RVAL)
        CHARACTER*(*) WHCH
        REAL          RVAL
        IF (ICFELL('MAPSTR - UNCLEARED PRIOR ERROR',1).NE.0) RETURN
        CALL MDSETR (WHCH,RVAL)
        IF (ICFELL('MAPSTR',2).NE.0) RETURN
        RETURN
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
C    @(#)gtge_num.f	20.12 5/3/00
C****************************************************************
C
C     File: gtge_num.f
C
C     Purpose: Routine to obtain DEI area, zone, and bus numbers
C
C     Input parameters:
C
C             nb       - bus external number
C             nptia    - GE area index
C             nptiz    - GE zone index
C             nptib    - GE bus index
C
C     Author: Walt Powell  Date: 21 May 1996
C     Called by: saveptid.f
C
C****************************************************************
      integer function gtge_num (nb, iptia, iptiz, iptib)
      integer nb, iptia, iptiz, iptib

      include 'ipfinc/parametr.inc'

      include 'ipfinc/blank.inc'
      include 'ipfinc/pti_data.inc'
      include 'ipfinc/filnam.inc'
      include 'ipfinc/lfiles.inc'
      include 'ipfinc/bus.inc'
      include 'ipfinc/prt.inc'
      include 'ipfinc/area.inc'
      include 'ipfinc/arcntl.inc'
      include 'ipfinc/pseudo_b.inc'
      include 'ipfinc/qsdup.inc'
      include 'ipfinc/qksrt.inc'

      integer MAXPTIRECORDS
      parameter (MAXPTIRECORDS = 16000)
      common /scratch/ count, array(4,MAXBUS), htable_2(MAXBUS),
     &                 nextptr_2(MAXBUS), count_newbus,
     &                 newbusno(MAXBUS), count_newzone,
     &                 newzoneno(MAXCZN), count_newown,
     &                 newownno(MAXOWN), tempc(MAXPTIRECORDS),
     &                 sort_tempc(MAXPTIRECORDS)
      integer array, count, htable_2, count_newbus, newbusno,
     &        count_newzone, newzoneno, count_newown, newownno,
     &        sort_tempc
      character tempc*80

      common /bpa_num / user_rule, num_area_rule, num_zone_rule,
     &                  num_owner_rule, num_default_rule,
     &                  area_rule(3,MAXCAR), zone_rule(3,MAXZON),
     &                  owner_rule(3,MAXOWN), default_rule(2,100),
     &                  owner_code_ge(MAXOWN)
      integer user_rule, num_area_rule, num_zone_rule, num_owner_rule,
     &        num_default_rule, area_rule, zone_rule, owner_rule,
     &        default_rule
      character owner_code_ge*4

      common /verify_pti/ num_verify, hash_verify(PTI_HASHSIZE),
     &                    next_verify(PTI_MAXBUS),
     &                    ptinum_verify(PTI_MAXBUS),
     &                    bpanum_verify(PTI_MAXBUS)
      integer num_verify, hash_verify, next_verify, ptinum_verify,
     &                    bpanum_verify

      logical finished, build_area, found
      integer fnd_ptib, fnd_ptiy, error, add_ptin, fnd_ptia, fnd_ptin,
     &        fnd_ptiz, fnd_ptiq, fnd_ptio, h, p, areanos(0:MAXCAR),
     &        rule, add_ptiz, add_ptio
      character tempc2*10
      external kmp_ptiz, swp_ptiz

      data build_area / .false. /

      save build_area, areanos

      error = 0           ! initialize error count
      iptia = 0
      iptiz = 0
      iptib = 0

      if (.not. build_area) then
        do i = 1, ntot
          bus_number(i) = 0
        enddo
        do i = 1, MAXCAR
          areanos(i) = 0
          array(1,i) = i
          array(2,i) = 0
        enddo

        key = 5
        call qiksrt (1, num_anam, kmp_ptiz, swp_ptiz)
        do i = 1, num_anam
          areanos(i) = pti_anum(array(1,i))
        enddo
        areanos(num_anam+1) = areanos(num_anam) + 10
        build_area = .true.
      endif
c
c     Get GE area number
c
      ia = jarzn(nb)
      tempc2 = arcnam(ia)
      iptia = fnd_ptia (tempc2)
      if (iptia .le. 0) then
        write (errbuf(1), 10041) bus(nb), base(nb),
     &    arcnam(jarzn(nb))
10041   format (' Bus (', a8, f7.1,
     &          ') has a non-existant GE area (', a, ')')
        call prterx ('E',1)
        error = 1
        go to 900
      endif
c
c     Get default GE zone number
c
      iptiz = fnd_ptiy (zone(nb))
      if (iptiz .eq. 0) then
        last = 0
        do i = 1, num_znam
          last = max0 (pti_znum(i), last)
        enddo
        iptiz = add_ptiz (zone(nb), last+1)
      endif
      iptiz_def = iptiz
c
c     Get default GE owner number
c
      iptio = fnd_ptiq (owner(nb))
      if (iptio .eq. 0) then
        last = 0
        do i = 1, num_onam
          last = max0 (pti_onum(i), last)
        enddo
        iptio = add_ptio (last+1, owner(nb))
      endif
      iptio_def = iptio
c
c     Get new GE bus number
c
      if (bus_number(nb) .eq. 0) then
        iptib = fnd_ptib (bus(nb), base(nb), arcnam(ia))
        if (iptib .gt. 0) then
          iptiz = fnd_ptiz (pti_zone(iptib))
          if (iptiz .eq. 0) then
            write (errbuf(1), 10000) pti_zone(iptib), bus(nb), base(nb)
10000       format (' Zone number (', i4,
     &        ') in translation file for bus (', a8, f7.1,
     &        ') is missing a corresponding name')
            call prterx ('W',1)
            iptiz = iptiz_def
          endif
          iptio = fnd_ptio (pti_owner(iptib))
          if (iptio .eq. 0) then
            write (errbuf(1), 10010) pti_owner(iptib), bus(nb), base(nb)
10010       format (' Owner number (', i4,
     &        ') in translation file for bus (', a8, f7.1,
     &        ') is missing a corresponding name')
            call prterx ('W',1)
            iptio = iptio_def
          endif
        endif
      else
        numptib = bus_number(nb)
        iptib = fnd_ptin (numptib)
        iptiz = fnd_ptiz (pti_zone(iptib))
        iptio = fnd_ptio (pti_owner(iptib))
      endif
      if (iptib .le. 0) then
        rule = 1               ! Step through all rules
        if (rule .eq. 1) then
          indx = 0
          do i = 1, num_area_rule
            if (area_rule(1,i) .eq. pti_anum(iptia)) indx = i
          enddo
          if (indx .gt. 0) then
            finished = .false.
            numptib = area_rule(2,indx)
            do while (.not. finished)
              iptib = fnd_ptin (numptib)
              if (iptib .le. 0) then
                area_rule(2,indx) = numptib + 1
                iptib = add_ptin (numptib)
                pti_name(iptib) = bus(nb)
                pti_base(iptib) = base(nb)
                pti_zone(iptib) = pti_znum(iptiz)
                pti_area(iptib) = pti_anum(iptia)
                pti_owner(iptib) = pti_onum(iptio)
                bus_number(nb) = numptib
                finished = .true.
              else
                numptib = numptib + 1
                finished = (numptib .gt. area_rule(3,indx))
              endif
            enddo
          endif
        endif
        if (iptib .eq. 0 .and. rule .eq. 1) rule = 2
        if (iptib .eq. 0 .and. rule .eq. 2) then
          indx = 0
          do i = 1, num_zone_rule
            if (zone_rule(1,i) .eq. pti_znum(iptiz)) indx = i
          enddo
          if (indx .gt. 0) then
            finished = .false.
            numptib = zone_rule(2,indx)
            do while (.not. finished)
              iptib = fnd_ptin (numptib)
              if (iptib .le. 0) then
                zone_rule(2,indx) = numptib + 1
                iptib = add_ptin (numptib)
                pti_name(iptib) = bus(nb)
                pti_base(iptib) = base(nb)
                pti_zone(iptib) = pti_znum(iptiz)
                pti_area(iptib) = pti_anum(iptia)
                pti_owner(iptib) = pti_onum(iptio)
                bus_number(nb) = numptib
                finished = .true.
              else
                numptib = numptib + 1
                finished = (numptib .gt. zone_rule(3,indx))
              endif
            enddo
          endif
        endif
        if (iptib .eq. 0 .and. rule .eq. 2) rule = 3
        if (iptib .eq. 0 .and. rule .eq. 3) then
          if (iptio .gt. 0) then
            indx = 0
            do i = 1, num_owner_rule
              if (owner_rule(1,i) .eq. pti_onum(iptio)) indx = i
            enddo
          else
            indx = 0
          endif
          if (indx .gt. 0) then
            finished = .false.
            numptib = owner_rule(2,indx)
            do while (.not. finished)
              iptib = fnd_ptin (numptib)
              if (iptib .le. 0) then
                owner_rule(2,indx) = numptib + 1
                iptib = add_ptin (numptib)
                pti_name(iptib) = bus(nb)
                pti_base(iptib) = base(nb)
                pti_zone(iptib) = pti_znum(iptiz)
                pti_area(iptib) = pti_anum(iptia)
                pti_owner(iptib) = pti_onum(iptio)
                bus_number(nb) = numptib
                finished = .true.
              else
                numptib = numptib + 1
                finished = (numptib .gt. owner_rule(3,indx))
              endif
            enddo
          endif
        endif
        if (iptib .eq. 0 .and. rule .eq. 3) rule = 4
        if (iptib .eq. 0 .and. rule .eq. 4) then
          indx = 0
          if (num_default_rule .gt. 0) indx = 1
          if (indx .gt. 0) then
            finished = .false.
            numptib = default_rule(1,indx)
            do while (.not. finished)
              iptib = fnd_ptin (numptib)
              if (iptib .le. 0) then
                default_rule(1,indx) = numptib + 1
                iptib = add_ptin (numptib)
                pti_name(iptib) = bus(nb)
                pti_base(iptib) = base(nb)
                pti_zone(iptib) = pti_znum(iptiz)
                pti_area(iptib) = pti_anum(iptia)
                pti_owner(iptib) = pti_onum(iptio)
                bus_number(nb) = numptib
                finished = .true.
              else
                numptib = numptib + 1
                finished = (numptib .gt. default_rule(2,indx))
              endif
            enddo
          endif
        endif
        if (iptib .le. 0) then
          write (errbuf(1), 10045) bus(nb), base(nb)
10045     format (' Bus (', a8, f7.1,
     &        ') could not be assigned a viable GE number')
          call prterx ('E',1)
          error = 1
          go to 900
        endif
c
c       Store bus number.  Future calls will use this shortcut.
c
        bus_number(nb) = pti_num(iptib)

        count_newbus = count_newbus + 1
        newbusno(count_newbus) = iptib

      endif
      if (iptib .gt. 0) then

        h = mod (pti_num(iptib), PTI_HASHSIZE) + 1
        p = hash_verify(h)
        found = .false.
        do while (p .gt. 0 .and. .not. found)
          if (ptinum_verify(p) .eq. pti_num(iptib)) then
            found = .true.
          else
            p = next_verify(p)
          endif
        enddo
        if (p .eq. 0) then
          num_verify = num_verify + 1
          p = num_verify
          next_verify(p) = hash_verify(h)
          hash_verify(h) = p
          ptinum_verify(p) = pti_num(iptib)
          bpanum_verify(p) = nb
        else if (bpanum_verify(p) .ne. nb) then
          nbx = bpanum_verify(p)
          write (errbuf(1), 10050) pti_num(iptib), bus(nb), base(nb),
     &      bus(nbx), base(nbx)
10050     format (' Duplicate number (', i5,
     &      ') shared by two buses (', a8, f7.1, ') and (', a8, f7.1,
     &      ')')
          call prterx ('W',1)
        endif
      endif

  900 gtge_num = error   ! return status (0/1) = (successful/failure)

      return
      end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
!============================================================
      SUBROUTINE TRMM_WINDOW_TRMM_WNFLT(COSVZA)
      INCLUDE 'fl_trmm_window.inc'
      REAL TRMM_WINDOW_TRMM_FULIOU_FILTER
!fu_sr      = FU-liou Radiance in each window band
!fu_sf      = FU-liou Flux in each window band

! sat_flt_sr = emmulation of TRMM FILTERED RADIANCE in each window band
! sat_unf_sr = emmulation of TRMM UNFILTERED RADIANCE in each window band

! sat_flt_r = emmulation of TRMM FILTERED RADIANCE
! sat_unf_r = emmulation of TRMM UNFILTERED RADIANCE
! sat_f     = emmulation of TRMM Window Flux

      SAT_FLT_R = 0.0
      SAT_UNF_R = 0.0
      SAT_F = 0.0
C *** Do construct BAND ***
      DO IB = 11, 13

         ADM(IB) = FU_SF(IB)/(3.14159*FU_SR(IB))
         FF = TRMM_WINDOW_TRMM_FULIOU_FILTER(IB,0,IWNCLD(IB),FU_SR(IB),
     1      COSVZA)
         UFF = TRMM_WINDOW_TRMM_FULIOU_FILTER(IB,1,IWNCLD(IB),FU_SR(IB),
     1      COSVZA)

         SAT_FLT_SR(IB) = FU_SR(IB)*FF
         SAT_UNF_SR(IB) = SAT_FLT_SR(IB)*UFF
         SAT_FLT_R = SAT_FLT_R + SAT_FLT_SR(IB)
         SAT_UNF_R = SAT_UNF_R + SAT_UNF_SR(IB)
! Flux method 1 ( SPECTRAL ADM  )
         SAT_F = SAT_F + SAT_UNF_SR(IB)*(3.14159*ADM(IB))

!  FLUX METHOD 2 ( DIFFUSE ANGLE Aprox)
!     ff  = TRMM_FULIOU_FILTER(ib,0,iwncld(ib),fu_sr(ib),0.6)
!     uff = TRMM_FULIOU_FILTER(ib,1,iwncld(ib),fu_sr(ib),0.6)
!     sat_flt_sf(ib) =  fu_sf(ib) * ff
!     sat_unf_sf(ib) = sat_flt_sf(ib) * uff
!     sat_f     = sat_f     + sat_unf_sf(ib)
! !

      END DO

!     print* , 'SCCWN',sum( sat_unf_sr ) / sum( sat_flt_sr )

      RETURN
      END
!============================================================
      REAL FUNCTION TRMM_WINDOW_TRMM_FULIOU_FILTER (IB, IDIR, ICLD,
     1   UNFILTERED_RADIANCE, COSVZA)
      PARAMETER (IBW = 1)
      INCLUDE 'fl_trmm_window.inc'
      REAL COEFS(4,3,0:1,0:1,2)
      INTEGER IBT(11:13)
      DATA IBT/ 3, 2, 1/
      DATA COEFS/ 5.137E-01, 5.002E-04, -4.887E-03, 3.724E-05, 7.156E-01
     1   , -2.305E-04, 9.184E-04, 1.272E-05, 6.067E-01, 6.400E-04,
     2   -9.066E-03, 3.633E-05, 1.382E+00, -8.756E-05, 7.637E-04,
     3   -3.077E-06, 1.398E+00, 4.567E-04, -1.826E-03, -2.524E-05,
     4   1.463E+00, -3.787E-04, 1.142E-03, -4.245E-05, 4.946E-01,
     5   3.123E-03, -2.532E-04, -6.721E-05, 7.187E-01, -1.150E-03,
     6   2.028E-04, 9.825E-05, 6.101E-01, -2.080E-03, -2.730E-03,
     7   2.223E-04, 1.384E+00, -3.153E-04, 1.790E-04, 6.385E-06,
     8   1.391E+00, 2.233E-03, -3.908E-04, -1.896E-04, 1.460E+00,
     9   8.689E-04, -2.329E-04, -1.286E-04, 5.137E-01, 5.002E-04,
     .   -4.887E-03, 3.724E-05, 7.156E-01, -2.305E-04, 9.184E-04,
     1   1.272E-05, 6.067E-01, 6.400E-04, -9.066E-03, 3.633E-05,
     2   1.346E+00, 3.661E-05, 5.024E-04, -1.379E-06, 1.398E+00,
     3   4.567E-04, -1.826E-03, -2.524E-05, 1.390E+00, 4.193E-04,
     4   -8.215E-03, 5.476E-05, 4.946E-01, 3.123E-03, -2.532E-04,
     5   -6.721E-05, 7.187E-01, -1.150E-03, 2.028E-04, 9.825E-05,
     6   6.101E-01, -2.080E-03, -2.730E-03, 2.223E-04, 1.344E+00,
     7   2.939E-04, 1.753E-04, -1.015E-05, 1.391E+00, 2.233E-03,
     8   -3.908E-04, -1.896E-04, 1.395E+00, -2.624E-03, -2.174E-03,
     9   2.697E-04/

      IF (IB.LT.11 .OR. IB.GT.13) THEN
         TRMM_WINDOW_TRMM_FULIOU_FILTER = 1.0000
         RETURN
      ENDIF

      IF((IDIR.EQ.0.OR.IDIR.EQ.1).AND.(ICLD.EQ.0.OR.ICLD.EQ.1))THEN
      ELSE
         STOP ' ERROR idir or icld NE to 0 or 1'
      ENDIF

      ICLD_USE = ICLD
      IF (ICLD_USE .EQ. 1) THEN
         IF (IB.EQ.11 .AND. UNFILTERED_RADIANCE.GT.8.0 .OR. IB.EQ.12
     1       .AND. UNFILTERED_RADIANCE.GT.8.0 .OR. IB.EQ.13 .AND.
     2      UNFILTERED_RADIANCE.GT.20.0) ICLD_USE = 0
      ENDIF

      TRMM_WINDOW_TRMM_FULIOU_FILTER = COEFS(1,IBT(IB),IDIR,ICLD_USE,1)
     1    + UNFILTERED_RADIANCE*COEFS(2,IBT(IB),IDIR,ICLD_USE,1) +
     2   COSVZA*COEFS(3,IBT(IB),IDIR,ICLD_USE,1) + UNFILTERED_RADIANCE**
     3   2*COEFS(4,IBT(IB),IDIR,ICLD_USE,1)

      RETURN

      END