text:
      SUBROUTINE READMO
      IMPLICIT DOUBLE PRECISION (A-H, O-Z)
      INCLUDE 'SIZES'
C
C MODULE TO READ IN GEOMETRY FILE, OUTPUT IT TO THE USER,
C AND CHECK THE DATA TO SEE IF IT IS REASONABLE.
C EXIT IF NECESSARY.
C
C
C
C  ON EXIT NATOMS    = NUMBER OF ATOMS PLUS DUMMY ATOMS (IF ANY).
C          KEYWRD    = KEYWORDS TO CONTROL CALCULATION
C          KOMENT    = COMMENT CARD
C          TITLE     = TITLE CARD
C          LABELS    = ARRAY OF ATOMIC LABELS INCLUDING DUMMY ATOMS.
C          GEO       = ARRAY OF INTERNAL COORDINATES.
C          LOPT      = FLAGS FOR OPTIMIZATION OF MOLECULE
C          NA        = ARRAY OF LABELS OF ATOMS, BOND LENGTHS.
C          NB        = ARRAY OF LABELS OF ATOMS, BOND ANGLES.
C          NC        = ARRAY OF LABELS OF ATOMS, DIHEDRAL ANGLES.
C          LATOM     = LABEL OF ATOM OF REACTION COORDINATE.
C          LPARAM    = RC: 1 FOR LENGTH, 2 FOR ANGLE, AND 3 FOR DIHEDRAL
C          REACT(200)= REACTION COORDINATE PARAMETERS
C          LOC(1,I)  = LABEL OF ATOM TO BE OPTIMIZED.
C          LOC(2,I)  = 1 FOR LENGTH, 2 FOR ANGLE, AND 3 FOR DIHEDRAL.
C          NVAR      = NUMBER OF PARAMETERS TO BE OPTIMIZED.
C          XPARAM    = STARTING VALUE OF PARAMETERS TO BE OPTIMIZED.
C
************************************************************************
C *** INPUT THE TRIAL GEOMETRY  \IE.  KGEOM=0\
C   LABEL(I) = THE ATOMIC NUMBER OF ATOM\I\.
C            = 99, THEN THE I-TH ATOM IS A DUMMY ATOM USED ONLY TO
C              SIMPLIFY THE DEFINITION OF THE MOLECULAR GEOMETRY.
C   GEO(1,I) = THE INTERNUCLEAR SEPARATION \IN ANGSTROMS\ BETWEEN ATOMS
C              NA(I) AND (I).
C   GEO(2,I) = THE ANGLE NB(I):NA(I):(I) INPUT IN DEGREES; STORED IN
C              RADIANS.
C   GEO(3,I) = THE ANGLE BETWEEN THE VECTORS NC(I):NB(I) AND NA(I):(I)
C              INPUT IN DEGREES - STORED IN RADIANS.
C  LOPT(J,I) = -1 IF GEO(J,I) IS THE REACTION COORDINATE.
C            = +1 IF GEO(J,I) IS A PARAMETER TO BE OPTIMIZED
C            =  0 OTHERWISE.
C *** NOTE:    MUCH OF THIS DATA IS NOT INCLUDED FOR THE FIRST 3 ATOMS.
C     ATOM1  INPUT LABELS(1) ONLY.
C     ATOM2  INPUT LABELS(2) AND GEO(1,2) SEPARATION BETWEEN ATOMS 1+2
C     ATOM3  INPUT LABELS(3), GEO(1,3)    SEPARATION BETWEEN ATOMS 2+3
C              AND GEO(2,3)              ANGLE ATOM1 : ATOM2 : ATOM3
C
************************************************************************
C
      DIMENSION LOPT(3,NUMATM)
      CHARACTER KEYWRD*241, KOMENT*81, TITLE*81, LINE*80, BANNER*80
      CHARACTER KEYS(80)*1, SPACE*1, SPACE2*2, CH*1, CH2*2
      CHARACTER ELEMNT*2, IDATE*24, GETNAM*80, NAME*4
      COMMON /KEYWRD/ KEYWRD
      COMMON /TITLES/ KOMENT,TITLE
      COMMON /GEOVAR/ NVAR, LOC(2,MAXPAR), IDUMY, XPARAM(MAXPAR)
      COMMON /PATH  / LATOM,LPARAM,REACT(200)
      COMMON /MESH  / LATOM1, LPARA1, LATOM2, LPARA2
      COMMON /ELEMTS/ ELEMNT(107)
      COMMON /MOLKST/ NUMAT,NAT(NUMATM),NFIRST(NUMATM),
     * NMIDLE(NUMATM),NLAST(NUMATM),NORBS,NELECS,NALPHA,NBETA,
     * NCLOSE,NOPEN,NDUMY,FRACT
      COMMON /OKMANY/ ISOK
      PARAMETER (MXDIM=MAXORB+NUMATM)
      COMMON /SYMRES/ TRANS,RTR,SIG,NAME,NAMO(MXDIM),INDX(MXDIM),
     * ISTA(2)
      COMMON /ISTOPE/ AMS(107)
      COMMON /GEOM  / GEO(3,NUMATM)
      COMMON /NUMCAL/ NUMCAL
      COMMON /GEOKST/ NATOMS,LABELS(NUMATM),
     1NA(NUMATM),NB(NUMATM),NC(NUMATM)
      COMMON /GEOSYM/ NDEP, LOCPAR(MAXPAR), IDEPFN(MAXPAR),
     1                      LOCDEP(MAXPAR)
      LOGICAL INT, AIGEO, ISOK
      SAVE SPACE, SPACE2, IREACT, INT
      DIMENSION COORD(3,NUMATM),VALUE(40)
      EQUIVALENCE (KEYS(1),KEYWRD)
      DATA SPACE, SPACE2/' ','  '/
      CONVTR=2.D0*ASIN(1.D0)/180.D0
      AIGEO=.FALSE.
   10 CONTINUE
C
      CALL GETTXT
      IF(INDEX(KEYWRD,'ECHO').NE.0)THEN
         REWIND 5
         ISOK=.FALSE.
         DO 50 I=1,1000
            READ(5,'(A)',END=60)KEYWRD
            DO 20 J=80,2,-1
   20       IF(KEYWRD(J:J).NE.' ')GOTO 30
            J=1
   30       DO 40 K=1,J
   40       IF(ICHAR(KEYWRD(K:K)).LT.32)KEYWRD(K:K)='*'
            WRITE(6,'(1X,A)')KEYWRD(1:J)
   50    CONTINUE
   60    CONTINUE
         REWIND 5
         CALL GETTXT
      ENDIF
      IF(INDEX(KEYWRD,'ECHO').NE.0)WRITE(6,'(''1'')')
      IF(KEYWRD(1:1) .NE. SPACE) THEN
         CH=KEYWRD(1:1)
         KEYWRD(1:1)=SPACE
         DO 70 I=2,239
            CH2=KEYWRD(I:I)
            KEYWRD(I:I)=CH
            CH=CH2
            IF(KEYWRD(I+1:I+2) .EQ. SPACE2) THEN
               KEYWRD(I+1:I+1)=CH
               GOTO 80
            ENDIF
   70    CONTINUE
         CH2=KEYWRD(240:240)
         KEYWRD(240:240)=CH
         KEYWRD(241:241)=CH2
   80    CONTINUE
      ENDIF
      IF(KOMENT(1:1) .NE. SPACE) THEN
         CH=KOMENT(1:1)
         KOMENT(1:1)=SPACE
         DO 90 I=2,79
            CH2=KOMENT(I:I)
            KOMENT(I:I)=CH
            CH=CH2
            IF(KOMENT(I+1:I+2) .EQ. SPACE2) THEN
               KOMENT(I+1:I+1)=CH
               GOTO 100
            ENDIF
   90    CONTINUE
         CH2=KOMENT(80:80)
         KOMENT(80:80)=CH
         KOMENT(81:81)=CH2
  100    CONTINUE
      ENDIF
      IF(TITLE(1:1) .NE. SPACE) THEN
         CH=TITLE(1:1)
         TITLE(1:1)=SPACE
         DO 110 I=2,79
            CH2=TITLE(I:I)
            TITLE(I:I)=CH
            CH=CH2
            IF(TITLE(I+1:I+2) .EQ. SPACE2) THEN
               TITLE(I+1:I+1)=CH
               GOTO 120
            ENDIF
  110    CONTINUE
         CH2=TITLE(80:80)
         TITLE(80:80)=CH
         TITLE(81:81)=CH2
  120    CONTINUE
      ENDIF
      DO 121 I=1,200
  121 REACT(I)=0.D0
      LATOM=0
      LPARAM=0
      IF(INDEX(KEYWRD,'OLDGEO').EQ.0) THEN
         NVAR=0
         NDEP=0
         IF(AIGEO.OR.INDEX(KEYWRD,'AIGIN').NE.0)THEN
            CALL GETGEG(5,LABELS,GEO,     NA,NB,NC,AMS,NATOMS,INT)
         IF(NVAR.EQ.0)THEN
         DO 122 J=1,3
         DO 122 I=1,NATOMS
  122    LOPT(J,I)=0
         ENDIF
         ELSE
            CALL GETGEO(5,LABELS,GEO,LOPT,NA,NB,NC,AMS,NATOMS,INT)
            IF(NATOMS.LT.0)THEN
               REWIND 5
               IF(NUMCAL.NE.1)THEN
                  WRITE(6,'(//,A)')'   GAUSSIAN INPUT REQUIRES STAND-ALO
     1NE JOB'
                  WRITE(6,'(/,A)')'   OR KEYWORD "AIGIN"'
                  STOP
               ENDIF
               AIGEO=.TRUE.
               GOTO 10
            ENDIF
         ENDIF
         IF(NATOMS.EQ.0)STOP
      ELSE
      DEGREE=90.D0/ASIN(1.D0)
      IF(NA(1).EQ.99)THEN
      DO 128 I=1,NATOMS
      DO 128 J=1,3
      LOPT(J,I)=1
  128 COORD(J,I)=GEO(J,I)
      LOPT(1,1)=0
      LOPT(2,1)=0
      LOPT(3,1)=0
      LOPT(2,2)=0
      LOPT(3,2)=0
      LOPT(3,3)=0
      CALL XYZINT(COORD,NATOMS,NA,NB,NC,DEGREE,GEO)
      ELSE
      DO 130 I=1,NATOMS
      DO 130 J=2,3
 130  GEO(J,I)=GEO(J,I)*DEGREE
      ENDIF
      ENDIF
      IF(INDEX(KEYWRD,'FORCE').NE.0 .AND. LABELS(NATOMS).EQ.107) THEN
      DO 131 I=1,NA(NATOMS)
      IF(LABELS(I).EQ.99)THEN
      WRITE(6,'(A)')' NO DUMMY ATOMS ALLOWED BEFORE TRANSLATION'
      WRITE(6,'(A)')' ATOM IN A FORCE CALCULATION'
      STOP
      ENDIF
  131 CONTINUE
      ENDIF
C
C
C OUTPUT FILE TO UNIT 6
C
C    WRITE HEADER
      IDATE=' '
      CALL fdate(IDATE)
      WRITE(6,'(1X,15(''*****''),''****'')')
C
C     CHANGE THE FOLLOWING LINE TO SUIT LOCAL ENVIRONMENT, IF DESIRED
C
      BANNER=' ** MOPAC (PUBLIC DOMAIN) FOR DEVELOPMENT USE '//
     1'ONLY.  NOT FOR PRODUCTION WORK  **'
      WRITE(6,'(A)')BANNER
C
C    THE BANNER DOES NOT APPEAR ANYWHERE ELSE.
C
      WRITE(6,'(1X,79(''*''))')
      LINE='   MNDO'
      IF(INDEX(KEYWRD,'MINDO') .NE. 0) LINE='MINDO/3'
      IF(INDEX(KEYWRD,'AM1')   .NE. 0) LINE='    AM1'
      IF(INDEX(KEYWRD,'PM3')   .NE. 0) LINE='    PM3'
      WRITE(6,'(/29X,A,'' CALCULATION RESULTS'',28X,///1X,15(''*****'')
     1,''****'' )')LINE(:7)
      WRITE(6,'('' *'',10X,''MOPAC:  VERSION '',F5.2,
     115X,''CALC''''D. '',A)') VERSON, IDATE
C
C CONVERT ANGLES TO RADIANS
      DO 140 J=2,3
C$DOIT VBEST
         DO 140 I=1,NATOMS
            GEO(J,I) = GEO(J,I) * CONVTR
  140 CONTINUE
C
C CHECK DATA
C
      NA(1)=0
      NB(1)=0
      NC(1)=0
      DO 150 I=1,NATOMS
         IF (LABELS(I) .LE. 0 ) THEN
            WRITE(6,'('' ATOMIC NUMBER OF '',I3,'' ?'')') LABELS(I)
            IF(I.EQ.1) THEN
               WRITE(6,'(A)')' THIS WAS THE FIRST ATOM'
            ELSE
               WRITE(6,'(A)')'    GEOMETRY UP TO, BUT NOT INCLUDING, THE
     1'//' FAULTY ATOM'
               NATOMS=I-1
               CALL GEOUT(6)
            ENDIF
            STOP
         ENDIF
         IF (  NA(I).GE.I.OR. NB(I).GE.I.OR. NC(I).GE.I
     1  .OR. (NA(I).EQ.NB(I))   .AND. I.GT.1
     2  .OR. (NA(I).EQ.NC(I).OR.NB(I).EQ.NC(I))  .AND. I.GT.2
     3  .OR.  NA(I)*NB(I)*NC(I).EQ.0  .AND. I.GT.3) THEN
            WRITE(6,'('' ATOM NUMBER '',I3,'' IS ILL-DEFINED'')') I
            IF(I.EQ.1)STOP
            WRITE(6,'(/,''  GEOMETRY READ IN'',/)')
            CALL GEOUT(6)
            STOP
         ENDIF
  150 CONTINUE
C
C WRITE KEYWORDS BACK TO USER AS FEEDBACK
      CALL WRTKEY(KEYWRD)
      WRITE(6,'(1X,14(''*****''),''*'',I3.3,''BY'',I3.3)')MAXHEV,MAXLIT
C
C FILL IN GEO MATRIX IF NEEDED
      IF(INDEX(KEYWRD,'OLDGEO').EQ.0.AND.INDEX(KEYWRD,'SYM') .NE. 0
     1.AND. NDEP.EQ.0) CALL GETSYM
      IF(NDEP.NE.0) CALL SYMTRY
C
C INITIALIZE FLAGS FOR OPTIMIZE AND PATH
      IFLAG = 0
      LATOM = 0
      NUMAT=0
      IF(NVAR.NE.0)THEN
         NUMAT=NATOMS
      ELSE
         DO 180 I=1,NATOMS
            IF(LABELS(I).NE.99.AND.LABELS(I).NE.107)NUMAT=NUMAT+1
            DO 180 J=1,3
               IF (LOPT(J,I) ) 160, 180, 170
C    FLAG FOR PATH
  160          CONVRT=1.D0
               IF ( IFLAG .NE. 0 ) THEN
                  IF(INDEX(KEYWRD,'STEP1').NE.0)THEN
                     LPARA1=LPARAM
                     LATOM1=LATOM
                     LPARA2=J
                     LATOM2=I
                     LATOM=0
                     IFLAG=0
                     GOTO 180
                  ELSE
                     WRITE(6,'('' ONLY ONE REACTION COORDINATE PERMITTED
     1'')')
                     STOP
                  ENDIF
               ENDIF
               LATOM  = I
               LPARAM = J
               IF(J.GT.1) CONVRT=0.01745329252D00
               REACT(1)  = GEO(J,I)
               IREACT=1
               IFLAG = 1
               GO TO 180
C    FLAG FOR OPTIMIZE
  170          NVAR = NVAR + 1
               LOC(1,NVAR) = I
               LOC(2,NVAR) = J
               XPARAM(NVAR)   = GEO(J,I)
  180    CONTINUE
      ENDIF
C READ IN PATH VALUES
      IF(IFLAG.EQ.0) GO TO 220
      IF(INDEX(KEYWRD,'NLLSQ').NE.0)THEN
         WRITE(6,'(A)')' NLLSQ USED WITH REACTION PATH; '//
     1'THIS OPTION IS NOT ALLOWED'
         STOP
      ENDIF
      IF(INDEX(KEYWRD,'SIGMA').NE.0)THEN
         WRITE(6,'(A)')' SIGMA USED WITH REACTION PATH; '//
     1'THIS OPTION IS NOT ALLOWED'
         STOP
      ENDIF
      IF(INDEX(KEYWRD,'STEP')+INDEX(KEYWRD,'POINTS').NE.0)THEN
         STEP=READA(KEYWRD,INDEX(KEYWRD,'STEP=')+5)
         NPTS=READA(KEYWRD,INDEX(KEYWRD,'POINT=')+6)
         IF(NPTS.GT.200)THEN
            WRITE(6,'(///,''    ONLY TWO HUNDRED POINTS ALLOWED IN REACT
     1'',''ION COORDINATE'')')
            STOP
         ENDIF
         IF(LPARAM.EQ.1.AND.STEP.LE.0)THEN
            WRITE(6,'(///,''    STEP FOR BOND LENGTH SHOULD BE SET POSIT
     1IVE '',''TO PREVENT TWO ATOMS COLLAPSE'')')
            STOP
         ENDIF
         GO TO 220
      ENDIF
  190 READ(5,'(A)',END=210) LINE
      CALL NUCHAR(LINE,VALUE,NREACT)
      IF(NREACT.EQ.0)GOTO 210
      DO 200 I=1,NREACT
         IJ=IREACT+I
         IF(IJ.GT.200)THEN
            WRITE(6,'(///,''    ONLY TWO HUNDRED POINTS ALLOWED IN REACT
     1ION'','' COORDINATE'')')
            STOP
         ENDIF
         REACT(IJ)=VALUE(I)*CONVRT
         IF(ABS(REACT(IJ)-REACT(IJ-1)).LT.1.D-5)THEN
            DUM1 = REACT(IJ)/CONVRT
            DUM2 = REACT(IJ-1)/CONVRT
            WRITE(6,'(///,'' TWO ADJACENT POINTS ARE IDENTICAL:  '',
     1 F7.3,2X,F7.3,/,'' THIS IS NOT ALLOWED IN A PATH CALCULATION'')')
     2 DUM1,DUM2
            STOP
         ENDIF
  200 CONTINUE
      IREACT=IREACT+NREACT
      GO TO 190
  210 CONTINUE
      DEGREE=1.D0
      IF(LPARAM.GT.1)DEGREE=90.D0/ASIN(1.D0)
      IF(IREACT.LE.1) THEN
         WRITE(6,'(//10X,'' NO POINTS SUPPLIED FOR REACTION PATH'')')
         WRITE(6,'(//10X,'' GEOMETRY AS READ IN IS AS FOLLOWS'')')
         CALL GEOUT(1)
         STOP
      ELSE
         WRITE(6,'(//10X,'' POINTS ON REACTION COORDINATE'')')
         WRITE(6,'(10X,8F8.2)')(REACT(I)*DEGREE,I=1,IREACT)
      ENDIF
      IEND=IREACT+1
      REACT(IEND)=-1.D12
C
C OUTPUT GEOMETRY AS FEEDBACK
C
  220 CALL WRTTXT(6)
      IF(INDEX(KEYWRD,'NOLOG').EQ.0)THEN
         OPEN(UNIT=11, FORM='FORMATTED', STATUS='UNKNOWN',
     +FILE=GETNAM('FOR011'))
         CALL WRTTXT(11)
      ENDIF
      CALL GEOUT(1)
      CALL GMETRY(GEO,COORD)
C
C  IF A POLYMER, EXPAND TO MERS
C
      IF(INDEX(KEYWRD,' MERS').NE.0)CALL MAKPOL(COORD)
      IF (INDEX(KEYWRD,'NOXYZ') .EQ. 0) THEN
         IF(INDEX(KEYWRD,'0SCF').NE.0)THEN
C
C  WRITE OUT CARTESIAN COORDINATES FOR USE AS A DATA-SET
C
            WRITE(6,'(A)')'   GEOMETRY IN CARTESIAN COORDINATE FORMAT'
            CALL WRTTXT(6)
            J=0
            DO 230 I=1,NATOMS
               IF(LABELS(I).NE.99)THEN
                  J=J+1
                  WRITE(6,'(2X,A,3(F19.13,I3))')
     1    ELEMNT(LABELS(I)),(COORD(K,J),1,K=1,3)
               ENDIF
  230       CONTINUE
         ELSE
            WRITE(6,'(//10X,''CARTESIAN COORDINATES '',/)')
            WRITE(6,'(4X,''NO.'',7X,''ATOM'',9X,''X'',
     1  9X,''Y'',9X,''Z'',/)')
            L=0
            DO 240 I=1,NATOMS
               IF(LABELS(I) .EQ. 99.OR.LABELS(I).EQ.107) GOTO 240
               L=L+1
               WRITE(6,'(I6,8X,A2,4X,3F10.4)')
     1  L,ELEMNT(LABELS(I)),(COORD(J,L),J=1,3)
  240       CONTINUE
         ENDIF
      ENDIF
      CALL SYMTRZ(COORD,C,NORBS,NORBS,.FALSE.,.FALSE.)
      WRITE(6,'(//''     MOLECULAR POINT GROUP   :   '',A4)') NAME
      IF(   INDEX(KEYWRD,' XYZ') .NE. 0 )THEN
         IF( NVAR .NE. 0 .AND.
     1 INT.AND.(NDEP .NE. 0 .OR.  NVAR.LT.3*NUMAT-6)) THEN
            IF(NDEP.NE.0)
     1WRITE(6,'(//10X,'' INTERNAL COORDINATES READ IN, AND SYMMETRY''
     2,/10X,'' SPECIFIED, BUT CALCULATION TO BE RUN IN CARTESIAN ''
     3,''COORDINATES'')')
            IF(NVAR.LT.3*NUMAT-6)
     1WRITE(6,'(//10X,'' INTERNAL COORDINATES READ IN, AND'',
     2'' CALCULATION '',/10X,''TO BE RUN IN CARTESIAN COORDINATES, '',
     3/10X,''BUT NOT ALL COORDINATES MARKED FOR OPTIMISATION'')')
            WRITE(6,'(//10X,'' THIS INVOLVES A LOGICALLLY ABSURD CHOICE'
     1',/10X,'' SO THE CALCULATION IS TERMINATED AT THIS POINT'')')
            STOP
         ENDIF
         SUMX=0.D0
         SUMY=0.D0
         SUMZ=0.D0
         DO 250 J=1,NUMAT
            SUMX=SUMX+COORD(1,J)
            SUMY=SUMY+COORD(2,J)
  250    SUMZ=SUMZ+COORD(3,J)
         SUMX=SUMX/NUMAT
         SUMY=SUMY/NUMAT
         SUMZ=SUMZ/NUMAT
         DO 260 J=1,NUMAT
            GEO(1,J)=COORD(1,J)-SUMX
            GEO(2,J)=COORD(2,J)-SUMY
  260    GEO(3,J)=COORD(3,J)-SUMZ
         NA(1)=99
         J=0
         NVAR=0
         DO 280 I=1,NATOMS
            IF(LABELS(I).NE.99)THEN
               J=J+1
               IF(J.EQ.1)THEN
                  K=0
               ELSEIF(J.LT.4)THEN
                  K=MIN(3,I-1)
               ELSE
                  K=3
               ENDIF
               DO 270 L=1,K
                  NVAR=NVAR+1
                  LOC(1,NVAR)=J
                  LOC(2,NVAR)=L
  270          XPARAM(NVAR)=GEO(L,J)
               LABELS(J)=LABELS(I)
            ENDIF
  280    CONTINUE
         NATOMS=NUMAT
      ELSE
         IF(NVAR.EQ.0) RETURN
         IF( .NOT. INT.AND.(NDEP .NE. 0 .OR.  NVAR.LT.3*NUMAT-6)) THEN
            IF(NDEP.NE.0)
     1WRITE(6,'(//10X,'' CARTESIAN COORDINATES READ IN, AND SYMMETRY''
     2,/10X,'' SPECIFIED, BUT CALCULATION TO BE RUN IN INTERNAL ''
     3,''COORDINATES'')')
            IF(NVAR.LT.3*NUMAT-6)
     1WRITE(6,'(//10X,'' CARTESIAN COORDINATES READ IN, AND'',
     2'' CALCULATION '',/10X,''TO BE RUN IN INTERNAL COORDINATES, '',
     3/10X,''BUT NOT ALL COORDINATES MARKED FOR OPTIMISATION'')')
            WRITE(6,'(//10X,''MOPAC, BY DEFAULT, USES INTERNAL COORDINAT
     1ES'',/10X,''TO SPECIFY CARTESIAN COORDINATES USE KEY-WORD :XYZ:'')
     2')
            WRITE(6,'(10X,''YOUR CURRENT CHOICE OF KEY-WORDS INVOLVES''
     1,'' A LOGICALLLY'',/10X,''ABSURD CHOICE SO THE CALCULATION IS'',
     2'' TERMINATED AT THIS POINT'')')
            STOP
         ENDIF
      ENDIF
      RETURN
      END

max_stars_repo_name: openmopac/MOPAC-archive

max_issues_repo_name: openmopac/MOPAC-archive

max_forks_repo_name: openmopac/MOPAC-archive


text:
      subroutine f
        print *, "Hello World"
      end

ext: f

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
Nails by Tam operates out of a converted house in the old east Davis area downtown.

$12 for a manicure, $24 for a pedicure. No deal for both, so its $36 for a manicure/pedicure.  $10 for a mens haircut.

  To learn about other businesses in Davis that provide professional aesthetic services, visit our Beauty Salons and Barber Shops page.

20070613 09:25:51 nbsp for a mans haircut, the best in town.  I have lived here for 6 years and this place is the only place I go to now that I have found it.  The guy takes his time and cuts your hair very well.  Plus it is only $9 for a haircut, cant beat that! Users/DanielB

20080109 08:11:47 nbsp Haircuts are $9, as advertised.  They dont wash, dry or style  just spray water and cut, which is why its still $9 for a womans haircut.  I thought he did a good job  I had him take off some length, thin it with a razor, and put a few layers in front.  He cut off more than I expected/asked for, so it might be a good idea to err on the conservative side when saying how many inches you want to trim. Users/ruthelizabeth

20080127 12:23:24 nbsp Got to agree with DanielB, The nice guy who cuts your hair takes his time and does excellent edge work. I have a very simple hair cut just 2 different razor lengths and yet he still takes the time to make sure its completely even. I always walk away with a better haircut than I ever got from supercuts where I usually shell out like 15 bucks for a 3 minute buzz. Thats why I usually leave a $6 tip. Users/ArmyDecoy

20100425 13:19:45 nbsp Went there today, the policy must of changed. They told me they did haircuts for men only. :( Users/Adalwolf

20110217 10:49:22 nbsp Hue always takes his time giving me a great hair cut! I highly recommend him. Users/BrianNeal

20121202 09:46:11 nbsp I went here based on a recommendation from a friend, and I was not disappointed.  Hue did a wonderful job cutting my hair and beard with high attention to detail, and it was very affordable as well.  I really recommend this place. Users/DavidLoyola

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
************************************************************************
* Read parameters for nonlinear solver of CC2D
*
* This routine reads the parameters for the nonlinear solver used
* in the nonlinear CC2D-iteration from a given file and stores
* them in the corresponding COMMON blockdiscretization.
*
* In:
*   MDATA  - Unit number to use for reading process
*   MSHOW  - Level of output in the initialization phase
*   CFNAME - Name of the file
************************************************************************

      SUBROUTINE RDNLCC (MDATA,MSHOW,CFNAME)

      IMPLICIT NONE

      INCLUDE 'cout.inc'
      INCLUDE 'cfileout.inc'
      INCLUDE 'cerr.inc'

      INCLUDE 'cnonlinsol_cc2d.inc'

C     parameters

      CHARACTER CFNAME*(*)
      INTEGER MSHOW,MDATA

C     local variables

      INTEGER I,IFMTS
      CHARACTER CSTR*(255)

C     Open DAT file

      IFMTS = 1
      CALL OF0 (MDATA,CFNAME,IFMTS)

C     Ignore the first 6 lines

      DO I=1,6
        READ (MDATA,*)
      END DO

C     Check version number

      READ (MDATA,*) I
      IF (I.NE.100) THEN
        WRITE (*,'(A)') 'Error in reading nonlinear-solver parameters'
        WRITE (*,'(A)') 'Version number incorrect'
        STOP
      END IF

      WRITE (CSTR,'(A)') ' Stationary nonlinear solver parameters'
      CALL CNOUTS (MSHOW,MTERM,MFILE,.TRUE.,CSTR)
      WRITE (CSTR,'(A)') '----------------------------------------'
      CALL CNOUTS (MSHOW,MTERM,MFILE,.TRUE.,CSTR)

C     Ignore separator

      READ (MDATA,*)

C     Minimum/Maximum number of iterations

      READ (MDATA,*) INLMIN
      READ (MDATA,*) INLMAX

      INLMIN = ABS(INLMIN)
      INLMAX = ABS(INLMAX)
      IF (INLMAX.LT.INLMIN) INLMAX=INLMIN

      WRITE (CSTR,1000)
     *      'Min. # of nonlin. iterations  : INLMIN = ', INLMIN
      CALL CNOUTS (MSHOW,MTERM,MFILE,.TRUE.,CSTR)

      WRITE (CSTR,1000)
     *      'Max. # of nonlin. iterations  : INLMAX = ', INLMAX
      CALL CNOUTS (MSHOW,MTERM,MFILE,.TRUE.,CSTR)

C     Error bounds, stopping criteria

      READ (MDATA,*) EPSD
      WRITE (CSTR,1001)
     *      'Limit for U-defects           : EPSD   = ', EPSD
      CALL CNOUTS (MSHOW,MTERM,MFILE,.TRUE.,CSTR)

      READ (MDATA,*) EPSDIV
      WRITE (CSTR,1001)
     *      'Limit for DIV-defects         : EPSDIV = ', EPSDIV
      CALL CNOUTS (MSHOW,MTERM,MFILE,.TRUE.,CSTR)

      READ (MDATA,*) EPSUR
      WRITE (CSTR,1001)
     *      'Limit for U-changes           : EPSUR  = ', EPSUR
      CALL CNOUTS (MSHOW,MTERM,MFILE,.TRUE.,CSTR)

      READ (MDATA,*) EPSPR
      WRITE (CSTR,1001)
     *      'Limit for P-changes           : EPSPR  = ', EPSPR
      CALL CNOUTS (MSHOW,MTERM,MFILE,.TRUE.,CSTR)

      READ (MDATA,*) DMPD
      WRITE (CSTR,1001)
     *      'Limit for defect improvement  : DMPD   = ', DMPD
      CALL CNOUTS (MSHOW,MTERM,MFILE,.TRUE.,CSTR)

C     Separator

      READ (MDATA,*)
      WRITE (MTERM,9001)

C     Parameters for optimal correction

      READ (MDATA,*) OMGMIN
      READ (MDATA,*) OMGMAX

      IF (OMGINI.LT.OMGMIN) OMGINI=OMGMIN
      IF (OMGINI.GT.OMGMAX) OMGINI=OMGMAX

      WRITE (CSTR,1001)
     *      'Lower limit for optimal OMEGA : OMGMIN = ', OMGMIN
      CALL CNOUTS (MSHOW,MTERM,MFILE,.TRUE.,CSTR)

      WRITE (CSTR,1001)
     *      'Upper limit for optimal OMEGA : OMGMAX = ', OMGMAX
      CALL CNOUTS (MSHOW,MTERM,MFILE,.TRUE.,CSTR)

      READ (MDATA,*) OMGINI
      WRITE (CSTR,1001)
     *      'Start value for optimal OMEGA : OMGINI = ', OMGINI
      CALL CNOUTS (MSHOW,MTERM,MFILE,.TRUE.,CSTR)

C     That's it.

      WRITE (MTERM,*)

      CLOSE (MDATA)

1000  FORMAT (A,I4)
1001  FORMAT (A,E16.8)
1002  FORMAT (A,2I4)
1003  FORMAT (A,A)

9000  FORMAT(79('-'))
9001  FORMAT(60('-'))

      END

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


text:
	SUBROUTINE CCBOXT (cval, clabel, ilabel, start, ill, jll,
     +			   idirc, scflag, zgrid, jhline, jvline, iret)
C************************************************************************
C* CCBOXT								*
C*									*
C* This subroutine draws a contour line through a box defined by four	*
C* grid points.								*
C*									*
C* CCBOXT  ( CVAL, CLABEL, ILABEL, START, ILL, JLL, IDIRC, SCFLAG, 	*
C*	     ZGRID, JHLINE, JVLINE, IRET )				*
C*									*
C* Input parameters:							*
C*	CVAL		REAL		Contour level			*
C*	CLABEL		CHAR*		Contour label			*
C*	ILABEL		INTEGER		Label type			*
C*	START		LOGICAL		Start flag			*
C*	SCFLAG		LOGICAL		Suppress small contour flag	*
C*                                                                      *
C* Work parameters:                                                     *
C*      ZGRID (KX*KY)   REAL            Array for grid points           *
C*      JHLINE (KX*KY)  INTEGER         Array for Horiz crossings       *
C*      JVLINE (KX*KY)  INTEGER         Array for Vert crossings        *
C*									*
C* Input and output parameters:						*
C*	ILL		INTEGER		Lower left x value for box	*
C*	JLL		INTEGER		Lower left y value for box	*
C*	IDIRC		INTEGER		Direction entering box		*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C*					 -3 = edge reached		*
C**									*
C* Log:									*
C* WOLFPLOT			Original code for SACC			*
C* M. desJardins/GSFC	 7/85	Adapted from AOIPS code for GEMPAK 3.1	*
C* M. desJardins/GSFC	 6/89	Rewritten for GEMPAK 5			*
C* J. Whistler/SSAI	 6/91	Set internal grids to size LLMXDG	*
C* M. desJardins/NMC	12/91	Renamed: CBOXIT-->CCBOXT		*
C* K. Brill/NMC		02/92	Shorten some too-long lines of code	*
C* G. Krueger/EAI	12/94	Conditional CCOEFF; CCGTBX call->CCCRSS	*
C* S. Jacobs/NCEP	 5/99	Modified check for multiple directions	*
C* C. Bailey/HPC	 6/06	Added contour label to calling sequence	*
C* C. Bailey/HPC	10/06	Added suppress small contour flag	*
C* S. Gilbert/NCEP      03/07   Added work arrays to call sequence      *
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	INCLUDE		'CONTUR.CMN'
C*
	CHARACTER*24	clabel(LLCLEV)
	LOGICAL		start, scflag
C*
	REAL		box (4), boxmod (4), zgrid (*)
        INTEGER         jhline (*), jvline (*)
	LOGICAL		tstr, tsts, tstl, edge
C*
	INTEGER		ippt1 (4), ippt2 (4), ippt3 (4), ippt4 (4)
	DATA		ippt1 / 4, 1, 2, 3 /, ippt2 / 3, 4, 1, 2 /,
     +			ippt3 / 2, 3, 4, 1 /, ippt4 / 1, 2, 3, 4 /
C
C*	This statement function gets the fractional distance along a
C*	side to the value 0.
C
	FIND0 ( x, y ) = x / ( x - y )
C------------------------------------------------------------------------
	iret = 0
	xll  = FLOAT ( ill )
	yll  = FLOAT ( jll )
C
C*	Set flag indicating we've been here.
C
	IF  ( idirc .eq. 1 )  THEN
	    jhline ( (jll-1)*isize+ill ) = jhline ( (jll-1)*isize+ill )
     +					   + 1
	  ELSE IF  ( idirc .eq. 2 )  THEN
	    jvline ( (jll-1)*isize+ill ) = jvline ( (jll-1)*isize+ill )
     +					   + 1
	  ELSE IF  ( idirc .eq. 3 )  THEN
	    jhline ( (jll)*isize+ill ) =
     +				   jhline ( (jll)*isize+ill ) + 1
	  ELSE IF  ( idirc .eq. 4 )  THEN
	    jvline ( (jll-1)*isize+ill+1 ) = jvline ( (jll-1)*isize+ill
     +							+1 ) + 1
	END IF
C
C*	First check for missing data.
C
	CALL CCCRSS ( 0, 0, ill, jll, cval, zgrid, iret )
	IF  ( iret .eq. -2 .or. iret .eq. -1 )  THEN
	    iret = -3
	    RETURN
	END IF
C
C*	Compute the coefficients to evaluate points inside box, if we
C*	are subboxing.
C
	IF ( numsub .gt. 2 ) CALL CCOEFF ( ill, jll, zgrid, ier )
C
C*	For a new contour line, find the point to start.
C
	IF  ( start )  THEN
	    CALL CCSTRT  ( cval, ill, jll, idirc, 1, illbox, jllbox,
     +			   zgrid, iret )
	    IF  ( iret .ne. 0 )  RETURN
C
C*	    Otherwise, check if point left by last box needs to be
C*	    modified for new coefficients.
C
	  ELSE
	    CALL CCCHCK  ( cval, ill, jll, idirc, illbox, jllbox,
     +                     zgrid, iret )
	    IF  ( iret .ne. 0 )  RETURN
	END IF
C
C*	Loop through all the grid boxes until an edge is reached.
C
	knt  =  1
	maxk =  4 * numsub
	edge = .false.
	DO WHILE  ( ( iret .eq. 0 ) .and. ( .not. edge ) )
C
C*	    Get corners of this subbox.
C
	    CALL CCGTSB  ( illbox, jllbox, ill, jll, cval, box, boxmod,
     +			   zgrid, iret )
C
C*	    Get points to look to right, straight and left.
C
	    ip1 = ippt1 ( idirc )
	    ip2 = ippt2 ( idirc )
	    ip3 = ippt3 ( idirc )
	    ip4 = ippt4 ( idirc )
C
C*	    Check for crossings to right, straight and left.
C
	    tstr = ( boxmod ( ip1 ) * boxmod ( ip2 ) .lt. 0. )
	    tsts = ( boxmod ( ip2 ) * boxmod ( ip3 ) .lt. 0. )
	    tstl = ( boxmod ( ip3 ) * boxmod ( ip4 ) .lt. 0. )
C
C*	    If there is more than one possible direction, the top
C*	    and bottom sides turn left, and the left and right
C*	    sides turn right.
C
	    IF  ( tstr .and. tstl )  THEN
		IF  ( ( idirc .eq. 1 ) .or. ( idirc .eq. 3 ) )  THEN
		    tstr = .false.
		  ELSE
		    tstl = .false.
		END IF
		tsts = .false.
	    END IF
C
C*	    Get new direction and points to interpolate.
C
	    IF  ( tstr )  THEN
		idirn = idirc + 1
		IF  ( idirn .eq. 5 )  idirn = 1
		dist = FIND0 ( box (ip1), box (ip2) )
		IF  ( ( idirc .eq. 3 ) .or. ( idirc .eq. 4 ) )
     +					dist = 1. - dist
	      ELSE IF  ( tsts )  THEN
		idirn = idirc
		dist  = FIND0 ( box (ip2), box (ip3) )
		IF  ( ( idirc .eq. 1 ) .or. ( idirc .eq. 4 ) )
     +					dist = 1. - dist
	      ELSE IF  ( tstl )  THEN
		idirn = idirc - 1
		IF  ( idirn .eq. 0 )  idirn = 4
		dist = FIND0 ( box (ip3), box (ip4) )
		IF  ( ( idirc .eq. 1 ) .or. ( idirc .eq. 2 ) )
     +					dist = 1. - dist
	      ELSE
		iret = -3
		RETURN
	    END IF
C
C*	    Compute value based on new direction.
C
	    IF  ( idirn .eq. 1 )  THEN
		x1 = xll + ( illbox - 1. + dist ) * fincxy
		y1 = yll + ( jllbox ) * fincxy
		jllbox = jllbox + 1
		IF  ( jllbox .eq. numsub )  THEN
		    edge = .true.
		    jll  = jll + 1
		    IF  ( jll .ge. jsize )  THEN
			iret = -3
			jhline ( (jll-1)*isize+ill ) =
     +					jhline ( (jll-1)*isize+ill ) + 1
		      ELSE IF  ( jhline ( (jll-1)*isize+ill ) .gt. 1 )
     +								    THEN
			iret = -5
		    END IF
		END IF
	      ELSE IF  ( idirn .eq. 2 )  THEN
		x1 = xll + ( illbox ) * fincxy
		y1 = yll + ( jllbox - 1. + dist ) * fincxy
		illbox = illbox + 1
		IF  ( illbox .eq. numsub )  THEN
		    edge = .true.
		    ill  = ill + 1
		    IF  ( ill .ge. isize )  THEN
			iret = -3
			jvline ( (jll-1)*isize+ill ) =
     +				      jvline ( (jll-1)*isize+ill ) + 1
		      ELSE IF  ( jvline ( (jll-1)*isize+ill ) .gt. 1 )
     +								    THEN
			iret = -5
		    END IF
		END IF
	      ELSE IF  ( idirn .eq. 3 )  THEN
		x1 = xll + ( illbox - 1. + dist ) * fincxy
		y1 = yll + ( jllbox - 1. ) * fincxy
		jllbox = jllbox - 1
		IF  ( jllbox .lt. 1 )  THEN
		    edge = .true.
		    jll  = jll - 1
		    IF  ( jll .lt. 1 )  THEN
			iret = -3
			jhline ( ill ) =
     +					  jhline ( ill ) + 1
		      ELSE IF  ( jhline ( (jll)*isize+ill ) .gt.
     +							      1 )  THEN
			iret = -5
		    END IF
		END IF
	      ELSE IF  ( idirn .eq. 4 )  THEN
		x1 = xll + ( illbox - 1. ) * fincxy
		y1 = yll + ( jllbox - 1. + dist ) * fincxy
		illbox = illbox - 1
		IF  ( illbox .lt. 1 )  THEN
		    edge = .true.
		    ill  = ill - 1
		    IF  ( ill .lt. 1 )  THEN
			iret = -3
			jvline ( (jll-1)*isize+1 ) =
     +					jvline ( (jll-1)*isize+1 )
     +						     + 1
		      ELSE IF  ( jvline ( (jll-1)*isize+ill+1 ) .gt. 1 )
     +								    THEN
			iret = -5
		    END IF
		END IF
	    END IF
C
C*	    Add point to buffer.
C
	    number = number + 1
	    IF  ( number .eq. 1024 )  THEN
		number = 1023
		xval1  = xval (1023)
		yval1  = yval (1023)
		CALL CCPLOT  ( cval, clabel, ilabel, scflag, iret )
		xval (1) = xval1
		yval (1) = yval1
		number   = 2
	    END IF
	    xval (number) = x1
	    yval (number) = y1
C
C*	    Set new direction.
C
	    idirc = idirn
C
C*	    Check for looping within box.
C
	    knt = knt + 1
	    IF  ( knt .gt. maxk )  THEN
		iret = -5
		edge = .true.
	    END IF
	END DO
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


!     =======================================
!     mpimod_dummy.f90
!     ----------------
!     This module replaces <mpimod.f90> for
!     single CPU runs
!
!     The module is shared by PUMA and PlaSim
!     =======================================

      subroutine mpbci(k) ! broadcast 1 integer
      implicit none
      integer k
      k=k
      return
      end

      subroutine mpbcin(k,n) ! broadcast n integer
      implicit none
      integer n
      integer :: k(n)
      n=n
      k(1)=k(1)
      return
      end

      subroutine mpbcr(p) ! broadcast 1 real
      implicit none
      real p
      p=p
      return
      end

      subroutine mpbcrn(p,n) ! broadcast n real
      implicit none
      integer n
      real :: p(n)
      p=p
      return
      end

      subroutine mpscin(k,n) ! scatter n integer
      implicit none
      integer n
      integer :: k(n)
      k=k
      return
      end

      subroutine mpscrn(p,n) ! scatter n real
      implicit none
      integer n
      real :: p(n)
      p=p
      return
      end

      subroutine mpscdn(p,n) ! scatter n double precision
      implicit none
      integer n
      real (kind=8) :: p(n)
      p=p
      return
      end

      subroutine mpscsp(pf,pp,klev) ! scatter spectral fields
      use pumamod
      implicit none
      integer klev
      real pf(NESP,klev)
      real pp(NSPP,klev)
      pp(1:NSPP,1:klev) = pf(1:NSPP,1:klev)
      return
      end

      subroutine mpscgp(pf,pp,klev) ! scatter gridpoint fields
      use pumamod
      implicit none
      integer klev
      real pf(NLON*NLAT,klev)
      real pp(NHOR,klev)
      pp(1:NHOR,1:klev) = pf(1:NHOR,1:klev)
      return
      end

      subroutine mpgasp(pf,pp,klev) ! gather spectral fields
      use pumamod
      implicit none
      integer klev
      real pf(NESP,klev)
      real pp(NSPP,klev)
      pf(1:NSPP,1:klev) = pp(1:NSPP,1:klev)
      return
      end

      subroutine mpgagp(pf,pp,klev) ! gather gridpoint fields
      use pumamod
      implicit none
      integer klev
      real pf(NHOR,klev)
      real pp(NHOR,klev)
      pf = pp
      return
      end

      subroutine mpgacs(pcs) ! gather cross sections
      use pumamod
      implicit none
      real pcs
      dimension pcs(NLAT,NLEV)
      pcs=pcs
      return
      end

      subroutine mpgallsp(pf,pp,klev) ! gather spectral to all
      use pumamod
      implicit none
      integer klev
      real pf(NESP,klev)
      real pp(NSPP,klev)
      pf(1:NSPP,1:klev) = pp(1:NSPP,1:klev)
      return
      end

      subroutine mpsum(psp,klev) ! sum spectral fields
      implicit none
      integer klev
      real :: psp(klev)
      psp=psp
      return
      end

      subroutine mpsumsc(psf,psp,klev) ! sum & scatter spectral
      use pumamod
      implicit none
      integer klev
      real psf(NESP,klev)
      real psp(NSPP,klev)
      psp(1:NSPP,1:klev) = psf(1:NSPP,1:klev)
      return
      end

      subroutine mpsumr(pr,kdim) ! sum kdim reals
      implicit none
      integer kdim
      real :: pr(kdim)
      pr=pr
      return
      end subroutine mpsumr

      subroutine mpsumbcr(pr,kdim) ! sum & broadcast kdim reals
      implicit none
      integer kdim
      real :: pr(kdim)
      pr=pr
      return
      end

!PBH dummy routine for scalars
      subroutine mpsumbcr1(pr) ! sum & broadcast kdim reals
      implicit none
      real pr
      pr=pr
      return
      end

!     ==================
!     SUBROUTINE MPABORT
!     ==================

      subroutine mpabort(ym)
      implicit none
      integer ilmess,ilabor,ilen,j,ioff

      character (len=* ) :: ym
      character (len=64) :: ystar = ' '
      character (len=64) :: yline = ' '
      character (len=64) :: yabor = 'Program aborted'
      character (len=64) :: ymess = ' '
      character (len=64) :: yhead = ' '

      ilmess = len_trim(ym)
      ilabor = len_trim(yabor)
      ilen   = 60

      do j = 1 , ilen+4
         ystar(j:j) = '*'
         yline(j:j) = '-'
      enddo

      ioff = 2
      if (ilmess < ilen-1) ioff = ioff + (ilen - ilmess) / 2
      ymess(1+ioff:ilmess+ioff) = trim(ym)

      ioff = 2
      if (ilabor < ilen-1) ioff = ioff + (ilen - ilabor) / 2
      yhead(1+ioff:ilabor+ioff) = trim(yabor)

      yline(1:1) = '*'
      ymess(1:1) = '*'
      yhead(1:1) = '*'
      yline(2:2) = ' '
      ymess(2:2) = ' '
      yhead(2:2) = ' '
      j = ilen + 4
      yline(j:j) = '*'
      ymess(j:j) = '*'
      yhead(j:j) = '*'
      j = ilen + 3
      yline(j:j) = ' '
      ymess(j:j) = ' '
      yhead(j:j) = ' '

      open (44,file='Abort_Message')
      write(44,'(A)') trim(ystar)
      write(44,'(A)') trim(yhead)
      write(44,'(A)') trim(yline)
      write(44,'(A)') trim(ymess)
      write(44,'(A)') trim(ystar)
      close(44)
      write( *,'(/,A)') trim(ystar)
      write( *,'(A)') trim(yhead)
      write( *,'(A)') trim(yline)
      write( *,'(A)') trim(ymess)
      write( *,'(A,/)') trim(ystar)

      stop
      end

      subroutine mpstart ! initialization
      use pumamod
      implicit none
      if (NPRO > 1) then
        print*,'error : scalar version compiled with NPRO > 1!'
        stop
      endif
      return
      end

      subroutine mpstop
      implicit none
      return
      end

      subroutine mpreadsp(ktape,p,kdim,klev)
      implicit none
      integer ktape,kdim,klev
      real p(kdim,klev)
      read (ktape) p
      return
      end

      subroutine mpreadgp(ktape,p,kdim,klev)
      implicit none
      integer ktape,kdim,klev
      real p(kdim,klev)
      read (ktape) p
      return
      end

      subroutine mpwritesp(ktape,p,kdim,klev)
      implicit none
      integer ktape,kdim,klev
      real p(kdim,klev)
      write (ktape) p
      return
      end

      subroutine mpwritegp(ktape,p,kdim,klev)
      implicit none
      integer ktape,kdim,klev
      real p(kdim,klev)
      write (ktape) p
      return
      end

      subroutine mpwritegph(ktape,p,kdim,klev,ihead)
      implicit none
      integer ktape,kdim,klev
      real p(kdim,klev)
!
      real(kind=4) zp(kdim,klev)
!
      integer ihead(8)
      write(ktape) ihead
      zp(:,:)=p(:,:)
      write(ktape) zp

      return
      end

      subroutine mpi_info(nprocess,pid)    ! get nproc and pid
      implicit none
      integer nprocess, pid
      nprocess = 1
      pid = 0
      return
      end subroutine mpi_info

      subroutine mpgetsp(yn,p,kdim,klev)
      implicit none
      character (len=*) :: yn
      integer kdim,klev
      real :: p(kdim,klev)
      call get_restart_array(yn,p,kdim,kdim,klev)
      return
      end subroutine mpgetsp

      subroutine mpgetgp(yn,p,kdim,klev)
      implicit none
      character (len=*) :: yn
      integer kdim,klev
      real :: p(kdim,klev)
      call get_restart_array(yn,p,kdim,kdim,klev)
      return
      end subroutine mpgetgp

      subroutine mpputsp(yn,p,kdim,klev)
      implicit none
      character (len=*) :: yn
      integer kdim,klev
      real :: p(kdim,klev)
      call put_restart_array(yn,p,kdim,kdim,klev)
      return
      end subroutine mpputsp

      subroutine mpputgp(yn,p,kdim,klev)
      implicit none
      character (len=*) :: yn
      integer kdim,klev
      real :: p(kdim,klev)
      call put_restart_array(yn,p,kdim,kdim,klev)
      return
      end subroutine mpputgp

      subroutine mpsurfgp(yn,p,kdim,klev)
      implicit none
      character (len=*) :: yn
      integer kdim,klev,iread
      real :: p(kdim,klev)
      call get_surf_array(yn,p,kdim,klev,iread)
      return
      end subroutine mpsurfgp

      subroutine mpmaxval(p,kdim,klev,pmax)
      implicit none
      integer kdim,klev
      real pmax
      real :: p(kdim,klev)
      pmax = maxval(p(:,:))
      return
      end subroutine mpmaxval

      subroutine mpsumval(p,kdim,klev,psum)
      implicit none
      integer klev,kdim
      real psum
      real :: p(kdim,klev)
      psum = sum(p(:,:))
      return
      end subroutine mpsumval


text:
C Copyright(C) 1999-2020 National Technology & Engineering Solutions
C of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C See packages/seacas/LICENSE for details

C=======================================================================
      SUBROUTINE MEMERR
C=======================================================================

C   --*** MEMERR *** (ETCLIB) Flag dynamic memory error
C   --   Written by Amy Gilkey - revised 02/23/88
C   --
C   --MEMERR prints an error message for a dynamic memory error.

      CALL MDSTAT (NERR, MEM)
      CALL MDERPT (2, NOVER)
      IF (NERR .LE. NOVER) THEN
         CALL PRTERR ('FATAL', 'Too much dynamic memory requested')
      ELSE
         CALL PRTERR ('PROGRAM', 'Dynamic allocation problem')
         CALL MDEROR (6)
      END IF

      RETURN
      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


! MIT License
!
! Copyright (c) 2020 SHEMAT-Suite
!
! 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.

!>    @brief initialisation routine for reinjection
!>    @param[in] ismpl local sample index
      SUBROUTINE user_init(ismpl)
        use arrays
        use mod_genrl
        use mod_linfos
        IMPLICIT NONE
        integer :: ismpl
        integer :: i

        IF (linfos(3)>=2) WRITE(*,*) ' ... user_init (REINJECT-3D)'
!     avoid side effects
        DO i = 1, nbc_data
!        init from the first set
          dbc_dataold(i) = dbc_data(i,1,1)
        END DO
        RETURN
      END

!>    @brief reinjection routine
!>    @param[in] ismpl local sample index
!>    @details
!>injection 1   13       8         5\n
!>injection 2   13       8        13\n
!>injection 3   13       8        21\n
!>production 1  13       18        5\n
!>production 2  13       18       13\n
!>production 3  13       18       21\n
      SUBROUTINE calc_user(ismpl)
        use arrays
        use mod_genrl
        use mod_time
        use mod_linfos
        use mod_wells3d
        IMPLICIT NONE
        integer :: ismpl
        integer :: i, j

        DOUBLE PRECISION tempc, pumpc, tempt, deltt
        LOGICAL lbc_found, lbc_found1, lbc_found2
        DOUBLE PRECISION deltat
        EXTERNAL deltat
        INTRINSIC dabs

        IF (linfos(3)>=2) WRITE(*,*) ' ... calc_user (REINJECT-3D)'

!VR --- special case, danger !!! ---

        deltt = deltat(simtime(ismpl),ismpl)
        IF ((simtime(ismpl)+deltt)/tunit>=stoptime) THEN
          tempc = 0.D0
          pumpc = 0.D0
          tempt = 0.D0
!       for all productions
          DO j = 1, num_pro
            lbc_found = .FALSE.
!          search in all boundary condition
            DO i = 1, nbc_data
              IF (ibc_data(i,cbc_i)==ipro(j) .AND. &
                  ibc_data(i,cbc_j)==jpro(j) .AND. &
                  ibc_data(i,cbc_k)==kpro(j) .AND. &
                  ibc_data(i,cbc_pv)==pv_head .AND. &
                  ibc_data(i,cbc_bt)==bt_neum) THEN
                IF (dbc_data(i,1,ismpl)>=0.D0) THEN
                  WRITE(*,'(2A,3I8,1A)') 'error: HEAD Neumann-boundary&
                    & point for production', &
                    ' has a value >= 0.0 at [', ipro(j), jpro(j), &
                    kpro(j), '] !'
                  STOP
                END IF
                lbc_found = .TRUE.
                tempc = tempc + conc(ipro(j),jpro(j),kpro(j),1,ismpl)* &
                  dabs(dbc_data(i,1,ismpl))
                tempt = tempt + temp(ipro(j),jpro(j),kpro(j),ismpl)* &
                  dabs(dbc_data(i,1,ismpl))
                pumpc = pumpc + dabs(dbc_data(i,1,ismpl))
                WRITE(*,'(1A,3(e12.4))') '!!!!! GPK2 center ', &
                  conc(ipro(j),jpro(j),kpro(j),1,ismpl), &
                  temp(ipro(j),jpro(j),kpro(j),ismpl), &
                  dbc_data(i,1,ismpl)
              END IF
            END DO
            IF ( .NOT. lbc_found) THEN
              WRITE(*,'(1A,3I8,1A)') 'error: no HEAD Neumann-boundary &
                &point found for production [', ipro(j), jpro(j), &
                kpro(j), '] !'
              STOP
            END IF
          END DO

          IF (pumpc/=0.D0) THEN
            tempc = tempc/pumpc
            tempt = tempt/pumpc
          END IF

          WRITE(*,'(4(1A,1e12.4))') '!!!!! reinjection rate = ', &
            pumpc, ' l/s,   concentration = ', tempc, &
            ' mmol/l, temperature = ', tempt, ' C, time =', &
            simtime(ismpl)

!       for all sources
          DO j = 1, num_in
            lbc_found = .FALSE.
            lbc_found1 = .FALSE.
            lbc_found2 = .FALSE.
!          search in all boundary condition
            DO i = 1, nbc_data
              IF (ibc_data(i,cbc_i)==iin(j) .AND. &
                  ibc_data(i,cbc_j)==jin(j) .AND. &
                  ibc_data(i,cbc_k)==kin(j) .AND. &
                  ibc_data(i,cbc_pv)==pv_conc .AND. &
                  ibc_data(i,cbc_bt)==bt_diri .AND. &
                  ibc_data(i,cbc_si)==1) THEN
                dbc_data(i,1,ismpl) = tempc
                lbc_found = .TRUE.
              END IF
              IF (ibc_data(i,cbc_i)==iin(j) .AND. &
                  ibc_data(i,cbc_j)==jin(j) .AND. &
                  ibc_data(i,cbc_k)==kin(j) .AND. &
                  ibc_data(i,cbc_pv)==pv_temp .AND. &
                  ibc_data(i,cbc_bt)==bt_diri) THEN
                dbc_data(i,1,ismpl) = tempt
                lbc_found1 = .TRUE.
              END IF
              IF (ibc_data(i,cbc_i)==iin(j) .AND. &
                  ibc_data(i,cbc_j)==jin(j) .AND. &
                  ibc_data(i,cbc_k)==kin(j) .AND. &
                  ibc_data(i,cbc_pv)==pv_head .AND. &
                  ibc_data(i,cbc_bt)==bt_neum) THEN
                dbc_data(i,1,ismpl) = pumpc
                lbc_found2 = .TRUE.
              END IF
            END DO
            IF ( .NOT. lbc_found) THEN
              WRITE(*,'(1A,3I8,1A)') 'error: no CONC Dirichlet-boundar&
                &y point found for injection [', iin(j), jin(j), &
                kin(j), '] !'
              STOP
            END IF
            IF ( .NOT. lbc_found1) THEN
              WRITE(*,'(1A,3I8,1A)') 'error: no TEMP Dirichlet-boundar&
                &y point found for injection [', iin(j), jin(j), &
                kin(j), '] !'
              STOP
            END IF
            IF ( .NOT. lbc_found2) THEN
              WRITE(*,'(1A,3I8,1A)') 'error: no HEAD Dirichlet-boundar&
                &y point found for injection [', iin(j), jin(j), &
                kin(j), '] !'
              STOP
            END IF
          END DO
        ELSE
!       reset the initial values
          DO i = 1, nbc_data
            dbc_data(i,1,ismpl) = dbc_dataold(i)
          END DO
        END IF

!VR --- special case, danger !!! ---

        RETURN
      END


text:
In 2008, Cal.NET merged with Sierra Advantage in Shingle Springs, and has changed their service and support practices drastically from their early years.  Wiki Comments and reviews of their service prior to this merge (dated 2005 through 2007) were Cal.NET/2005 through 2007 Reviews moved here.  Most of these reviews were negative.  Cal.NET also used to offer MDL service to Glacier Point Apartments Glacier Point and Saratoga West Apartments Saratoga West apartment complexes, ending this service their their 2008 changes. An angry and extensive review of that service and mostly of the company circa 2003/2004 can be read Cal.NET/Business and MDL Review 2004 here.

Cal.net currently provides dialup and DSL services, as well as Web Hosting, and other high speed data alternatives such as T1 lines.

20090212 14:11:40 nbsp If you are viewing this page and are having issues with your service, please feel free to contact me by either leaving a comment on my entry, or by emailing me at cthompson at cal dot net. I will do my best to get your issues resolved within one business day. Users/ChrisThompson

20090831 10:44:11 nbsp Greetings, i have a client with their DNS servers that are incorrectly set to 4.2.2.1 & 4.2.2.2. I just called your office # at 5307578088 but no one picked up, I had to leave a message. Hopefully someone will contact me within the hour?

I just need to know what the correct DNS servers are for Cal.net/AT&T and i can be contacted at Cybershield  5307572010 or 9165431635

Thank you in advance for any help you can provide us and the community, im sure others are misconfigured. Users/Cybershield
 Those are L3 DNS servers. Theyre actually anycasted, so they work well just about anywhere. I would hardly consider using them to be a misconfiguration. Nonetheless, Cal.net puts the IP addresses of their DNS servers on their website in a http://main.cal.net/tech_support/calnetsystem_settings.php very easy to find place, if you bothered to check.  Users/WilliamLewis, who is not affiliated with Cal.NET in any way.


text:
c
c
c     ###################################################
c     ##  COPYRIGHT (C)  1990  by  Jay William Ponder  ##
c     ##              All Rights Reserved              ##
c     ###################################################
c
c     ############################################################
c     ##                                                        ##
c     ##  subroutine extra2  --  atomwise user defined Hessian  ##
c     ##                                                        ##
c     ############################################################
c
c
c     "extra2" calculates second derivatives of any additional
c     user defined potential energy contribution for a single
c     atom at a time
c
c
      subroutine extra2 (i)
      implicit none
      include 'sizes.i'
      include 'atoms.i'
      include 'hessn.i'
      integer i
c
c
c     compute the Hessian elements for extra energy terms
c
c     do j = 1, n
c        hessx(1,j) = hessx(1,j) + ......
c        hessy(2,j) = hessy(2,j) + ......
c        hessz(3,j) = hessz(3,j) + ......
c     end do
c
      return
      end

max_stars_repo_name: htna/HCsbLib

max_issues_repo_name: htna/HCsbLib

max_forks_repo_name: htna/HCsbLib


text:
	SUBROUTINE GQSATN ( navtyp, imgnam, iret )
C************************************************************************
C* GQSATN								*
C*									*
C* This subroutine returns satellite navigation information.  The       *
C* navigation type may be up to 8 characters long.  The name may be up 	*
C* to 256 characters in length.						*
C*									*
C* GQSATN ( NAVTYP, IMGNAM, IRET )					*
C*									*
C* Output parameters:							*
C*	NAVTYP		CHAR*8		Satellite navigation type	*
C*	IMGNAM		CHAR*256	Satellite image name		*
C*	IRET		INTEGER		Return code			*
C**									*
C* Log:									*
C* M. desJardins/GSFC	10/85						*
C* M. desJardins/GSFC	 3/86	Fixed error in return code		*
C* M. desJardins/GSFC	 6/88	Documentation				*
C* L. Williams/EAi       3/94   Removed blank comments from header      *
C* S. Jacobs/NMC	 7/94	Removed AOI projection			*
C* M. Linda/GSC		12/95	Removed NPGS satellite navigation	*
C* A. Hardy/GSC          6/98   Cleaned up prolog                       *
C* S. Jacobs/NCEP	 6/00	Increased image name from 80 to 256 char*
C************************************************************************
	INCLUDE 	'FUNCCODE.PRM'
	INCLUDE 	'ERROR.PRM'
C*
	CHARACTER*(*) 	imgnam, navtyp
C
	INTEGER		isend (2), ircv (64)
C------------------------------------------------------------------------
C*	Load input parameters into buffer and write them to the mailbox.
C
	isend (1) = 2
	isend (2) = FQSATN
C
	CALL GPUT  ( isend, 2, iret )
	IF  ( iret .ne. NORMAL )  RETURN
C
	CALL GGET  ( iret, 1, ierr )
	IF  ( ierr .ne. NORMAL )  THEN
	    iret = ierr
	    RETURN
	END IF
C
	CALL GGET  ( ircv, 2,  ierr )
	IF  ( ierr .ne. NORMAL )  THEN
	    iret = ierr
	    RETURN
	END IF
	CALL ST_ITOS  ( ircv, 2, nc, navtyp, ierr )
C
	CALL GGET  ( ircv, 64, ierr )
	IF  ( ierr .ne. NORMAL )  THEN
	    iret = ierr
	    RETURN
	END IF
	CALL ST_ITOS  ( ircv, 64, nc, imgnam, ier )
C*
	RETURN
	END


text:
program ch1202
    ! The abs Generic Function
    implicit none

    complex :: c = cmplx(1.,1.)
    real :: r = 10.9
    integer :: i = -27

    print *, abs(c)
    print *, abs(r)
    print *, abs(i)
end program

max_stars_repo_name: hechengda/fortran

max_issues_repo_name: hechengda/fortran

max_forks_repo_name: hechengda/fortran


text: Ecological Management and Restoration is a major offered by the College of Agricultural and Environmental Sciences as of 2010.  The major is designed for students who are interested in managing and restoring wildland and rangeland plant communities. Courses are selected to provide an interdisciplinary background that encompasses ecology, plant biology applied plant biology, and the social sciences.


text:
c	Simplified TE model in .MEX format

c	N. L. Ricker, U. of Washington

c	STATE VARIABLES are molar holdups [kmol]
c	in the specified location unless stated otherwise.

c	 1	A in reactor
c	 2	B in reactor
c	 3	C in reactor
c	 4	D in reactor
c	 5	E in reactor
c	 6	F in reactor
c	 7	G in reactor
c	 8	H in reactor
c	 9	A in separator
c	10	B in separator
c	11	C in separator
c	12	D in separator
c	13	E in separator
c	14	F in separator
c	15	G in separator
c	16	H in separator
c	17	A in feed zone
c	18	B in feed zone
c	19	C in feed zone
c	20	D in feed zone
c	21	E in feed zone
c	22	F in feed zone
c	23	G in feed zone
c	24	H in feed zone
c	25	D in stripper bottoms
c	26	E in stripper bottoms
c	27	F in stripper bottoms
c	28	G in stripper bottoms
c	29	H in stripper bottoms
c	30	Separator internal energy [MJ/kmol]
c	31	Stripper internal energy [MJ/kmol]

c	INPUTS are flows [kmol/h] unless noted otherwise:

c	 1	Feed stream 1 (pure A)
c	 2	Feed stream 2 (pure D)
c	 3	Feed stream 3 (pure E)
c	 4	Feed stream 4 (A & C)
c	 5	Compressor valve [Percent]
c	 6	Purge (stream 9)
c	 7	Separator underflow (stream 10)
c	 8	Stripper underflow (stream 11)
c	 9	Steam flow [kg/h]
c	10	Reactor temperature [C]
c	11	Condenser valve [%]
c	12	A in stream 4 [Mole Percent]
c	13	B in stream 4 [Mole Percent]
c	14	Reaction 1 activity factor [Percent].
c	15	Reaction 2 activity factor [Percent].
c	16	Stream 10 bias flow [kmol/h].
c	17	Reactor/Sep flow parameter [Percent].
c	18	Feed/Reactor flow parameter [Percent].
c	19	NOT USED in this version, i.e., ignored.
c	20	Adjustment to VLE of D to G in separator [Percent]
c	21	Adjustment to H VLE in separator [Percent]
c	22	C bias flow to feed zone [kmol/h]
c	23	D bias flow to feed zone [kmol/h]
c	24	E bias flow to feed zone [kmol/h]
c	25	Reaction 3 activity factor [Percent]
c	26	Adjustment to VLE in reactor [Percent]
c	27	Bias correction to compressor performance curve
c	28	Bias correction to compressor power prediction
c	29	Product density correction [Percent]
c	30	# equil. stages for D in stripper
c	31	# equil. stages for E in stripper
c	32	# equil. stages for F in stripper
c	33	Condenser cooling water inlet temperature [C]
c	34	UA in condenser [kW/K]
c	35	Stripper heat loss [kW]

c	OUTPUTS are mole % in stream unless noted otherwise:

c	 1	Reactor pressure [kPa]
c	 2	Reactor liq. holdup [%]
c	 3	Separator pressure [kPa]
c	 4	Separator liq. holdup [%]
c	 5	Product liq. holdup [%]
c	 6	Feed zone pressure [kPa]
c	 7	Total feed entering reactor (stream 6) [kscmh]
c	 8	A in reactor feed (stream 6)
c	 9	B in reactor feed (stream 6)
c	10	C in reactor feed (stream 6)
c	11	D in reactor feed (stream 6)
c	12	E in reactor feed (stream 6)
c	13	F in reactor feed (stream 6)
c	14	A in purge (stream 9)
c	15	B in purge (stream 9)
c	16	C in purge (stream 9)
c	17	D in purge (stream 9)
c	18	E in purge (stream 9)
c	19	F in purge (stream 9)
c	20	G in purge (stream 9)
c	21	H in purge (stream 9)
c	22	D in product (stream 11)
c	23	E in product (stream 11)
c	24	F in product (stream 11)
c	25	G in product (stream 11)
c	26	H in product (stream 11)
c	27	Volumetric production rate [m^3/h].
c	28	Recycle rate (stream 8) [kscmh]
c	29	Compressor power [kW]
c	30	Stripper Temperature [C]
c	31	Separator Temperature [C]
c	32	Condenser coolant outlet temp [C]
c	33	Cost [cents/kmol product]
c	34	Rate of reaction 1 [kmol G produced/h]
c	35	Rate of reaction 2 [kmol H produced/h]
c	36	Rate of reaction 3 [kmol F produced/h]
c	37	Partial pressure of A in reactor [kPa]
c	38	Partial pressure of C in reactor [kPa]
c	39	Partial pressure of D in reactor [kPa]
c	40	Partial pressure of E in reactor [kPa]
c
C---------------------------------------------------------------------

      SUBROUTINE mexFunction(NLHS, PLHS, NRHS, PRHS)
      INTEGER PLHS(*), PRHS(*)
      INTEGER NLHS, NRHS
C
      INTEGER mxCreateFull, mxGetPr, mxGetPi, mexGetMatrixPtr
C
      DOUBLE PRECISION mxGetScalar
C
C THE ABOVE SUBROUTINE, ARGUMENT, AND FUNCTION DECLARATIONS ARE STANDARD
C IN ALL FORTRAN MEX FILES.
C---------------------------------------------------------------------
C

c	This is a MEX interface routine for Matlab version 4.x

c	It is designed to emulate a Simulink S-function.

c	The following will hold pointers:

	integer t,x,u,flag

c		The number of states, inputs, and outputs are specified
c		by the parameters nx, nu, and ny, respectively.
c		Also define corresponding character constants
c		lenx, lenu, to use in error messages, and double precision
c		constants for the SIZES vector.

	integer iflag,m,n,nx,nu,ny
	character*(*) lenu,lenx
	double precision dnx, dnu, dny
	parameter (nx=31, nu=35, ny=40, lenx='31', lenu='35', dnx=31.,
     +		dnu=35.,dny=40.)

	double precision 	rx(nx),rx0(nx),dxdt(nx),sizes(6),rt,ru(nu),
     +			ry(ny)

	character*80 message	! For debug messages to screen

c		The following is the appropriate definition of the SIZES
c		vector for SIMULINK system functions when there are no
c		discrete-time states and no direct feed-through from input
c		to output.

	data sizes/dnx, 0., dny, dnu, 0., 0./

	data rx0/nx*0.d0/

c		Check the number of input arguments

	if (nrhs .gt. 0 .and. nrhs .ne. 4) then
		call mexErrMsgTxt('TEest5 needs 4 input variables.')
	else if (nrhs .le. 0) then
		iflag=0
	else

c			There are 4 input arguments.  Set pointers and copy the
c			inputs to local storage.

		m=mxGetM(prhs(1))				! Get size of "t" array.
		n=mxGetN(prhs(1))
		if (m.gt.0 .and. n.gt.0) then		! If it's non-empty...
			rt=mxGetScalar(prhs(1))		! gets time value (a scalar).
		else
			rt=0.0				! Otherwise set to zero.
		end if

		m=mxGetM(prhs(4))				! FLAG
		n=mxGetN(prhs(4))
		if (m.gt.0 .and. n.gt.0) then
			if (m.eq.1 .and. n.eq.1) then
				iflag=abs(nint(mxGetScalar(prhs(4))))
			else
				call mexErrMsgTxt('FLAG must be a scalar.')
			end if
		else
			iflag=0
		end if

		m=mxGetM(prhs(2))				! X
		n=mxGetN(prhs(2))
		if (m.gt.0 .and. n.gt.0) then
		   if (min(m,n).ne.1 .or. max(m,n) .ne. nx) then
			call mexErrMsgTxt('X must be a vector, length '//lenx)
		   end if
		   x=mxGetPr(prhs(2))
		   call mxCopyPtrToReal8(x,rx,nx)		! copy x into local storage
		else if (iflag.eq.1 .or. iflag.eq.3) then
			call mexErrMsgTxt('X must be supplied.')
		end if

		m=mxGetM(prhs(3))				! U
		n=mxGetN(prhs(3))
		if (m.gt.0 .and. n.gt.0) then
			if (min(m,n).ne.1 .or. max(m,n) .ne. nu) then
				call mexErrMsgTxt
     +			('U must be a vector, length '//lenu)
			end if
			u=mxGetPr(prhs(3))
			call mxCopyPtrToReal8(u,ru,nu)		! copy u into local storage
		else if (iflag.eq.1 .or. iflag.eq.3) then
			call mexErrMsgTxt('U must be supplied.')
		end if

	end if

c		Now return outputs according to the value of "flag"

	if (iflag.eq.0) then

c			Return default initial states and size info

		plhs(1)=mxCreateFull(6,1,0)				! "sizes"
		call mxCopyReal8ToPtr(sizes,mxGetPr(plhs(1)),6)
		if (nlhs.gt.1) then
			plhs(2)=mxCreateFull(nx,1,0)			! "x0"
			call mxCopyReal8ToPtr(rx0,mxGetPr(plhs(2)),nx)
		end if
		if (nlhs.gt.2) then
			plhs(3)=mxCreateFull(0,0,0)			! "xstr" (empty)
		end if

	else if (iflag.eq.1 .or. iflag.eq.3) then

		call TEsimf(nx,nu,ny,rx,ru,ry,dxdt)    	! Main simulation call

		if (iflag.eq.1) then

c				Return dxdt.  If there are 2 LHS arguments
c				return y as the second one.  (non-standard!!)

			plhs(1)=mxCreateFull(nx,1,0)
			call mxCopyReal8ToPtr(dxdt,mxGetPr(plhs(1)),nx)  	! returns dxdt to Matlab
			if (nlhs.ge.2) then
				plhs(2)=mxCreateFull(ny,1,0)
				call mxCopyReal8ToPtr(ry,mxGetPr(plhs(2)),ny)  	! returns y to Matlab
			end if
		else
			plhs(1)=mxCreateFull(ny,1,0)
			call mxCopyReal8ToPtr(ry,mxGetPr(plhs(1)),ny)  	! returns y to Matlab
		end if
	else

c			For all other flags, return an empty matrix.

		plhs(1)=mxCreateFull(0,0,0)

	end if

	return
	end

C---------------------------------------------------------------------

	SUBROUTINE TEsimf(nx,nu,ny,x,u,y,dxdt)

c	Calculates outputs and rates of change of states.

c	This subroutine can be used independently, i.e., without
c	the MEX interface.

c	NOTE:  all non-integer variables are double precision.

c	Inputs to routine:

c	nx	number of states.  Set to integer 31.
c	nu	number of inputs.  Set to integer 35.
c	ny	number of outputs.  Set to integer 40.
c	x	current state vector, length nx.
c	u	current input vector, length nu.

c	Outputs from routine

c	y	output vector, length ny.
c	dxdt	time derivatives of states, length nx.

	IMPLICIT DOUBLE PRECISION (A-H,O-Z)

	integer nx,nu,ny,nc
	double precision cf2cm,dzero,pctGH

	double precision x(nx),u(nu),y(ny),dxdt(nx)

	parameter (nc=8,cf2cm=0.028317,TauTr=0.08333,TauTs=0.08333,
     +		TauF=0.005,TauC=0.08333,dzero=0.0d0,Rg=8.314)

	double precision 	rlmol(nc),slmol(nc),rvmol(nc),svmol(nc),
     +			plmol(nc),fvmol(nc),xlr(nc),xls(nc),RR(4),
     +			ppr(nc),pps(nc),xvr(nc),xvs(nc),xvf(nc),
     +			Fcmol(nc),CFbias,CWbias,dP,Cflow,Crecyc,
     +			F8mass,Kval(nc),Pvc(nc),Sfac(nc),F4i(nc),
     +			F5i(nc),pvapr(nc),pvaps(nc),pvapc(nc),
     +			x4(8),x11(8)

	double precision 	mwts(nc),molvol(nc),
     +			sfr(nc),gami(nc),gamr(nc)

	integer code

	data mwts/2.0,25.4,28.0,32.0,46.0,48.0,62.0,76.0/	! molec. wts.
	data molvol/3*0.0,0.1070,0.1260,0.1463,0.1013,0.1231/	! [m3/kmol]
	data sfr/6*1.0,0.07,0.04/					! stripping factors
	data gami/8*1.0/				! sep. activ. coeffs.
	data gamr/8*1.0/				! react. activ. coeffs.
	data Cpw/4.18/				! specific heat of water [kJ/kg-K]
	data Cvw/2725./				! condenser coolant coefficient  [kg/h-%]

c		Initialization

	do i=1,8
		x4(i)=0.0
		F4i(i)=0.0
	end do

	Fsteam=u(9)				! steam flow [kg/h]

c		Get disturbances from input vector. Flow disturbances are scaled so
c		they're comparable to manipulated variables (0-100% basis)

	x4(1)=u(12)/100.			! mol frac A in stream 4 (A+C feed)
	x4(2)=u(13)/100.			! mol frac B in stream 4 (A+C feed)
	x4(3)=1.0-x4(1)-x4(2)		! get C in stream 4 by difference
	R1F=u(14)/100.			! Reactivity factor, reaction 1
	R2F=u(15)/100.			! Reactivity factor, reaction 2
	F10b=u(16)				! Bias adjustment for stream 10 flow [kmol/h]
	RSflow=u(17)/100.			! Adjusts flow between reactor and separator
	FRflow=u(18)/100.			! Adjusts flow between feed zone and reactor
	fracGH=u(19)/100.			! Adjusts for impurity in product (components other
						! than G and H)
	do i=4,7
		gami(i)=u(20)/100.	! Adjusts VLE of D to G in separator
	end do
	gami(8)=u(21)/100.		! Adjusts VLE of H in separator
	Cbias=u(22)				! C flow bias at feed point
	Dbias=u(23)				! D flow bias at feed point
	Ebias=u(24)				! E flow bias at feed point
	R3F=u(25)/100.			! Reactivity factor, reaction 3
	do i=4,8
		gamr(i)=u(26)/100.	! Adjusts VLE in reactor
	end do
	CFbias=u(27)			! Bias correction to compressor flow curve
	CWbias=u(28)			! Bias correction to compressor power curve
	Pvolcf=u(29)/100.			! Product density correction
	DNeq=u(30)				! Equil. stages for D
	ENeq=u(31)				! Equil. stages E
	FNeq=u(32)				! Equil. stages F
	Tw0=u(33)				! Condenser coolant inlet temp [C]
	UAcon=u(34)*3600.			! UA in condenser [kJ/h-K]
	QLoss=u(35)*3600.			! Heat loss in stripper [kJ/h]

c		Get molar holdups in reactor, separator, and feed vapor.
c		Assumes A, B, C insoluble in liquid.

	rlsum=0.0
	slsum=0.0
	plsum=0.0
	rvsum=0.0
	svsum=0.0
	VLR=0.0
	VLS=0.0
	VLP=0.0

	do i=1,3
		rlmol(i)=0.0
		slmol(i)=0.0
		plmol(i)=0.0
	end do

	do i=4,nc
		rlmol(i)=x(i)	! gets liq. mol from state vector
		slmol(i)=x(i+8)
		plmol(i)=x(i+21)
		rlsum=rlsum+rlmol(i)		! accumulate total moles
		slsum=slsum+slmol(i)
		plsum=plsum+plmol(i)
		VLR=VLR+rlmol(i)*molvol(i)	! accumulate liquid volume [m3]
		VLS=VLS+slmol(i)*molvol(i)
		VLP=VLP+plmol(i)*molvol(i)
	end do

	fvsum=0.0
	do i=1,nc
		fvmol(i)=x(i+16)			! gets feed vapor from state vector
		fvsum=fvsum+fvmol(i)		! total feed vapor moles
	end do

c	Temperatures in reactor, separator, & stripper

	TCR=u(10)				! reactor temp [C]
	TKR=TCR+273.2			!              [K]
	enint=1000.0*x(30)*slsum	! get extensive internal energy [kJ]
	call UTfunc(enint,slmol,TCS)	! get separator temp [C]
	TKS=TCS+273.2			!                [K]
	enint=1000.0*x(31)*plsum	! get extensive internal energy [kJ]
	call UTfunc(enint,plmol,TCC)	! get stripper temp [C]

c		Volumes in reactor, separator, & feed zone

	VTR=36.8				! total reactor volume [m3]
	VTS=99.1				! total separator volume [m3]
	VTV=150.0				! total feed vapor holdup volume [m3]
	VVR=VTR-VLR				! vapor in reactor [m3]
	VVS=VTS-VLS				! vapor in separator [m3]

c		Mole fractions, partial pressures, etc.

	PTR=0.0
	PTS=0.0
	PTV=fvsum*Rg*(273.2+86.1)/VTV		! Total pressure in feed zone [kPa]
							! (assumes constant temperature of 86.1 C)

	do i=1,3							! for A, B, C
		ppr(i)=x(i)*Rg*TKR/VVR		! pp in reactor [kPa]
		pps(i)=x(i+8)*Rg*TKS/VVS	! pp in separator [kPa]
		PTR=PTR+ppr(i)			! total P in reactor [kPa]
		PTS=PTS+pps(i)			! total P in separator [kPa]
		xlr(i)=0.0
		xls(i)=0.0
		x11(i)=0.0
	end do

	call VPfunc(tcr,pvapr)						! vap. pres. in reactor [kPa]
	call VPfunc(tcs,pvaps)						! vap. pres. in separator [kPa]
	do i=4,nc								! For D,E,F,G,H
		xlr(i)=rlmol(i)/rlsum					! reactor liquid mol frac
		ppr(i)=pvapr(i)*xlr(i)*gamr(i)			! VLE
		xls(i)=slmol(i)/slsum					! Repeat for separator
		pps(i)=pvaps(i)*xls(i)*gami(i)				! Note use of VLE adjustment
		PTR=PTR+ppr(i)						! Reactor total press [kPa]
		PTS=PTS+pps(i)						! Separator total press [kPa]
		x11(i)=plmol(i)/plsum					! Product mol fractions
	end do

	wtmolr=0.0
	wtmolf=0.0
	wtmols=0.0
	do i=1,8
		xvr(i)=ppr(i)/PTR			! reactor vapor mol frac
		xvs(i)=pps(i)/PTS			! sep vapor mol frac
		xvf(i)=fvmol(i)/fvsum		! feed vapor mol frac
		wtmolr=wtmolr+xvr(i)*mwts(i)	! average mol wt of reactor vapor
		wtmols=wtmols+xvs(i)*mwts(i)	! average mol wt of separator vapor (& purge)
		wtmolf=wtmolf+xvf(i)*mwts(i)	! average mol wt of feed vapor
	end do

c		Rate Laws.  All rates in [kmol/h].

      RR(1)=VVR*DEXP(44.06-42600.0/1.987/TKR)*R1F*(PPR(1)**1.08)
     +		*(PPR(3)**0.311)*(PPR(4)**0.874)
      RR(2)=VVR*DEXP(10.27-19500.0/1.987/TKR)*R2F*(PPR(1)**1.15)
     +		*(PPR(3)**0.370)*(PPR(5)**1.00)
      RR(3)=VVR*DEXP(59.50-59500.0/1.987/TKR)*R3F*PPR(1)
     +		*(0.77*PPR(4)+PPR(5))

c		Define total flows [kmol/h].  First get those defined directly
c		as inputs.

	F1=u(1)
	F2=u(2)
	F3=u(3)
	F4=u(4)
	F9=u(6)
	F10=u(7) - F10b
	F11=u(8)

c		Now use the compressor curve and valve setting to get the recycle rate.
c		dP is pressure drop across compresor (feed zone P - separator P)
c		Cflow is compressor throughput in kg/h.
c		Crecyc is back-flow through compressor valve [kg/h]
c		F8mass is net mass flow of recycle (stream 8) [kg/h]
c		F8 is corresponding molar flow.
c		Cpower is compressor power [kW]

	dP=PTV-PTS				! force non-negative pressure drop
	Cflow=-0.05152*dP*dP - 92.52*dP + 1.17e5 + CFbias
	Crecyc=u(5) * 60.83 * sqrt(dP)
	F8mass=Cflow-Crecyc
	F8=F8mass/wtmols									! [kmol/h]
	Cpower=TKS*Rg*(-9.365e-8*dP*dP + 4.300e-4*dP - 8.667e-3 + CWbias)	! [kJ/kmol]
	Cpower=Cpower*(Cflow/wtmols/3600.0)				! [kW]

c		Stripper modeling.  Removal of D, E, & F is based on the KSB equations.
c		Product composition is calculated on basis of zero A, B, and C.
c		Note that stream 10 bias is included in liquid rate for consistency
c		with steady-state calculations in TEest6_ss.m.

	call VPfunc(Tcc,pvapc)
	do i=4,8
		Sfac(i)=(pvapc(i)/PTV)*F4/(F10+F10b)		! Ki * V / L
	end do

	sfr(4)=1.0d0 - (1.0d0 - Sfac(4))/(1.0d0 - Sfac(4)**(DNeq+1.0d0))
	sfr(5)=1.0d0 - (1.0d0 - Sfac(5))/(1.0d0 - Sfac(5)**(ENeq+1.0d0))
	sfr(6)=1.0d0 - (1.0d0 - Sfac(6))/(1.0d0 - Sfac(6)**(FNeq+1.0d0))
	sfr(7)=Sfac(7)*0.809
	sfr(8)=Sfac(8)*0.720

c		Stripper balance to get molar rates in stream 5.

	F4i(1)=x4(1)*F4
	F4i(2)=x4(2)*F4
	F4i(3)=x4(3)*F4
	F5i(1)=F4i(1)
	F5i(2)=F4i(2)
	F5i(3)=F4i(3)

	do i=4,8
		F5i(i)=sfr(i)*F10*xls(i)
	end do

c		Stripper balance to get combined feeds (streams 1+2+3+5)
c		Also account for flow biases.

	Fcmol(1)=F5i(1)+F1
	Fcmol(2)=F5i(2)
	Fcmol(3)=F5i(3)+Cbias
	Fcmol(4)=F5i(4)+F2+Dbias
	Fcmol(5)=F5i(5)+F3+Ebias
	Fcmol(6)=F5i(6)
	Fcmol(7)=F5i(7)
	Fcmol(8)=F5i(8)
	Fcomb=0.0
	do i=1,8
		Fcomb=Fcomb+Fcmol(i)
	end do

c		Pressure drop vs. flow to get flows of streams 6 & 7.

	F6=(2413.7/wtmolf)*((PTV-PTR)**0.5)*FRflow
	if (ptr .gt. ptv) then
		F6=-F6
	end if
	F7=(5722.0/wtmolr)*((PTR-PTS)**0.5)*RSflow
	if (pts .gt. ptr) then
		F7=-F7
	end if

c		Balances to get rates of change of states.  First do
c		convection terms (in - out).

	do i=1,8
		dxdt(i)=F6*xvf(i) - F7*xvr(i)					! Reactor
		dxdt(i+8)=F7*xvr(i) - (F8+F9)*xvs(i)
     +		 		- F10*xls(i)				! Separator
     		dxdt(i+16)=Fcmol(i)+F8*xvs(i)-F6*xvf(i)			! Feed zone
	end do

	dmdt1=0.0	! Accumulation of liquid mass in separator
	dmdt2=0.0	! Accumulation of liquid mass in stripper
	do i=4,8
		dxdt(i+21)=(1.0-sfr(i))*F10*xls(i) - F11*x11(i)		! Stripper
		dmdt1=dmdt1+dxdt(i+8)
		dmdt2=dmdt2+dxdt(i+21)
	end do

c		Now account for reaction.

	dxdt(1)=dxdt(1) - rr(1) - rr(2) - 0.333*rr(3)
	dxdt(2)=dxdt(2)
	dxdt(3)=dxdt(3) - rr(1) - rr(2)
	dxdt(4)=dxdt(4) - rr(1) - rr(3)
	dxdt(5)=dxdt(5) - rr(2) - 0.333*rr(3)
	dxdt(6)=dxdt(6) + rr(3)
	dxdt(7)=dxdt(7) + rr(1)
	dxdt(8)=dxdt(8) + rr(2)

c		Energy balances

	T4=45.0						! temperature of stream 4 (from base case)
	call HVfunc(T4,F4i,H4)				! enthalpy of stream 4 [kJ/h]
	call HVfunc(tcc,F5i,H5)				! enthalpy of stream 5 [kJ/h]
	call HVfunc(tcr,xvr,H7)
	H7=F7*H7						! enthalpy of reactor effluent [kJ/h]
	call HVfunc(tcs,xvs,H89)
	H89=(F8+F9)*H89					! enthalpy of vapor leaving separator
	call HLfunc(tcs,xls,H10)
	H10=F10*H10						! enthalpy of liquid leaving separator
	call HLfunc(tcc,x11,H11)
	H11=F11*H11						! enthalpy of product
	Qsteam=Fsteam*2260.				! stripper heat duty (steam) [kJ/h]
	Fw=Cvw*u(11)					! condenser coolant rate  [kg/h]

c		Combine heat balance on condenser coolant with heat transfer rate
c		equation to get coolant outlet temperature and duty in condenser.

	alph=Fw*Cpw/UAcon					! Dimensionless group
	Tw=((2*alph-1.)*Tw0 + 2*Tcs)/(2*alph+1.)	! Temp of coolant leaving condenser [C]
	Qcond=UAcon*(0.5*(Tw+Tw0) - Tcs)			! Condenser duty [kJ/h]

c		Accumulation of internal energy

	dxdt(30)=0.001*(Qcond+H7-H89-H10)			! separator heat balance [MJ/h]
	dxdt(31)=0.001*(Qsteam-QLoss+H4+H10-H5-H11)	! stripper heat balance [MJ/h]

c		Convert to intensive form [MJ/kmol]

	dxdt(30)=(dxdt(30) - x(30)*dmdt1)/slsum
	dxdt(31)=(dxdt(31) - x(31)*dmdt2)/plsum

c		Calculate the outputs.

	y(1)=PTR-101.					! Reactor pressure 	[kPa guage]
	y(2)=5.263*VLR - 12.105				! Reactor liq. 		[%]
	y(3)=PTS-101.					! Separator pressure	[kPa guage]
	y(4)=12.28*VLS - 10.53				! Separator liquid	[%]
	y(5)=22.58*VLP - 49.03				! Product liq. holdup 	[%]
	y(6)=PTV-101.					! Pressure in feed zone [kPa guage]
	y(7)=F6/44.79					! Stream 6 flowrate 	[kscmh]
	do i=1,6
		y(i+7)=xvf(i)*100.0			! Reactor feed 		[mol %]
		y(i+13)=xvs(i)*100.0			! Purge 			[mol %]
	end do

	y(20)=xvs(7)*100.0				! G in purge		[mol %]
	y(21)=xvs(8)*100.0				! H in purge		[mol %]
	y(22)=x11(4)*100.0				! D in product		[mol %]
	y(23)=x11(5)*100.0				! E in product		[mol %]
	y(24)=x11(6)*100.0				! F in product		[mol %]
	y(25)=x11(7)*100.0				! G in product		[mol %]
	y(26)=x11(8)*100.0				! H in product		[mol %]
	y(27)=F11*VLP/plsum*Pvolcf			! Production 		[m^3/h]
	y(28)=F8/44.79					! Recycle rate		[kscmh]
	y(29)=Cpower					! Compressor power	[kW]
	y(30)=Tcc						! Stripper Temperature	[C]
	y(31)=Tcs						! Separator temperature  [C]
	y(32)=Tw						! Condenser outlet temperature  [C]

	purgls=F9*(221.*xvs(1)
     +		+ 618.*xvs(3) + 2210.*xvs(4)
     +		+ 1460.*xvs(5) + 1790.*xvs(6)
     +		+ 3040.*xvs(7) + 2290.*xvs(8))	! Purge loss cost [$/h]
	prodls=F11*(2210.*x11(4) + 1460.*x11(5) + 1790.*x11(6))	! Product loss cost [$/h]

	y(33)=(5.36*Cpower + 3.18*u(9) + purgls + prodls)/F11     	! [cents/kmol]

	y(34)=RR(1)						! G production		[kmol/h]
	y(35)=RR(2)						! H production		[kmol/h]
	y(36)=RR(3)						! F production		[kmol/h]
	y(37)=PPR(1)					! A PP in reactor		[kPa]
	y(38)=PPR(3)					! C PP in reactor		[kPa]
	y(39)=PPR(4)					! D PP in reactor		[kPa]
	y(40)=PPR(5)					! F PP in reactor		[kPa]

	return
	end

C---------------------------------------------------------------------

	subroutine HLfunc(T,x,hl)

c		Calculates liquid enthalpy for TE problem [kJ].
c
c		T = temperature [C]
c		x = vector of molar compositions (length 8).
c		hl= liquid enthalpy [kJ]

c		Reference is liquid at 100 C.  Components A, B,
c		and C are ignored for liquid-phase calculations.

	double precision T,x(8),hl

	double precision mwts(8),CpL(8)

        data mwts/2.0,25.4,28.0,32.0,46.0,48.0,62.0,76.0/       ! molec. wts.
        data CpL/3*0.0d0,7.66,4.17,4.45,2.55,2.45/              ! liquid ht cap

        hl=0.0
        do i=4,8
        	hl=hl+(T-100.0)*CpL(i)*x(i)*mwts(i)
        end do
        return
        end

C---------------------------------------------------------------------

	subroutine HVfunc(T,x,hv)

c		Calculates vapor enthalpy for TE problem [kJ].
c
c		T = temperature [C]
c		x = vector of molar compositions (length 8).
c		hv= enthalpy of vapor [kJ]

c		Reference is liquid at 100 C.

	double precision T,x(8),hv

	double precision mwts(8),CpV(8),dhv(8)

        data mwts/2.0,25.4,28.0,32.0,46.0,48.0,62.0,76.0/       ! molec. wts.
        data CpV/14.6,2.04,1.05,1.85,1.87,2.02,.712,.628/       ! vapor ht cap
        data dhv/3*0.0d0,202.,372.,372.,523.,486./              ! latent heats

        hv=0.0
        do i=1,8
        	hv=hv+(dhv(i)+(T-100.0)*CpV(i))*x(i)*mwts(i)
        end do
        return
        end

C---------------------------------------------------------------------

	subroutine UTfunc(U,x,T)

c		Calculates temperature [C] for specified
c		internal energy and liquid composition.
c		For TE problem.

c		U	Internal energy of liquid [kJ]
c		x	molar composition [kmol]
c		T	calculated temperature [C]

c		The first 3 elements of x are ignored

	double precision U,x(8),T

	double precision mwts(8),CpL(8)

	data mwts/2.0,25.4,28.0,32.0,46.0,48.0,62.0,76.0/       ! molec. wts.
	data CpL/3*0.0d0,7.66,4.17,4.45,2.55,2.45/              ! liquid ht cap

c		First get heat capacity, then back-calculate
c		temperature.  Assumes heat capacity is
c		independent of temperature (which is the
c		case for the data supplied by Downs and Vogel).

	T=0.0
	do i=4,8
		T=T+x(i)*CpL(i)*mwts(i)
	end do
	T=100.0+U/T
	return
	end

C---------------------------------------------------------------------

	subroutine VPfunc(T,pvap)

c		Calculates pure-component vapor pressure [kPa]

c		T	temperature [C]

	double precision T,pvap(8)
	integer i

	double precision AVP(8),BVP(8),CVP(8)

	data AVP/3*20.72,20.81,2*21.24,21.32,22.10/ 		! A in Antoine
	data BVP/3*0.0,-1444.0,2*-2114.0,-2748.0,-3318.0/ 	! B in Antoine
	data CVP/3*273.0,259.0,2*266.0,233.0,250.0/		! C in Antoine

	do i=1,8
        	pvap(i)=1.d-3*dexp(AVP(i)+
     +			BVP(i)/(CVP(i)+T))
     	end do
     	return
     	end

max_stars_repo_name: camaramm/tennesse-eastman-challenge

max_issues_repo_name: camaramm/tennesse-eastman-challenge

max_forks_repo_name: camaramm/tennesse-eastman-challenge


subroutine funct3d

      use stel_kinds, only: dp, p4
      use name0, only: czero, cp25, cp5, c1p0, c1p5
      use name1, only: nvac, nznt, nrztd, mnmax
      use realsp, only: workb, &
                        r1, ru, rv, z1, zu, zv, rcon, zcon, &
                        ru0, zu0, rcon0, zcon0, gcon
      use rforces, only: worka, &
                         armn, brmn, crmn, azmn, bzmn, czmn, &
                         blmn, clmn
      use scalars, only: nrzt, mns, ns, neqs, iter1, iter2, &
                         hs, ohs, irst,  meven, modd, iequi, dnorm, &
                         twopi, voli, isigng
      use xstuff, only: xc, gc
      use mnarray, only: xrz3, xrz4
      use scalefac, only: scalxc, sqrts, shalf, wint
      use fsqu, only: wp, fsqr, fsqz
      use profs, only: mass, vp, pres
      use inputdat, only: gam, nvacskip
      use extfld, only: ivac, ivac2
      use magfield, only: curpol, rbtor, ctor, bsqvac, bsqsav, rbsq, dbsq
      use time, only: timer

      implicit none

      ! TODO: more elegant definition of these BLAS functions
      real(kind=dp) :: ddot
      external dcopy

      real(kind=dp) :: rmnc(mnmax)
      real(kind=dp) :: zmns(mnmax)
      real(kind=dp) :: lmns(mnmax)
      real(kind=dp) :: xm(mnmax) ! for handing over to vacuum()
      real(kind=dp) :: xn(mnmax) ! for handing over to vacuum()

      real(kind=dp) :: rax(nznt)
      real(kind=dp) :: zax(nznt)

      integer       :: lodd, l, lk
      real(kind=dp) :: bz0
      real(kind=p4) :: t4
      real(kind=dp) :: timeon, timeoff

      ! let the madness begin ...
      real(kind=dp), pointer :: guu(:)
      real(kind=dp), pointer :: guv(:)
      real(kind=dp), pointer :: gvv(:)
      real(kind=dp), pointer :: lu(:)
      real(kind=dp), pointer :: lv(:)

      guu => workb(1+13*nrztd:14*nrztd) ! rcon(1+nrztd)
      guv => workb(1+15*nrztd:16*nrztd) ! zcon(1+nrztd)
      gvv => workb(1+ 6*nrztd: 7*nrztd) ! z1
      lu  => worka(1+10*nrztd:12*nrztd) ! czmn
      lv  => worka(1+ 4*nrztd: 6*nrztd) ! crmn

      lodd = 1+nrzt

      ! EXTRAPOLATE M>2 MODES AT JS = 2
      call extrap(xc,          xc(1+  mns), &
                  xc(1+2*mns), xc(1+3*mns), &
                  xc(1+4*mns), xc(1+5*mns), &
                  xrz3, xrz4)

      ! temporary re-use of gc for scaled xc
      do l = 1,neqs
        gc(l) = xc(l) * scalxc(l)
      enddo

      ! INVERSE FOURIER TRANSFORM TO S,THETA,ZETA SPACE
      ! R, Z, AND LAMBDA ARRAYS IN FOURIER SPACE
      do l = 1, nrzt
        lu(l) = c1p0 !!! This makes lu == (1 + d(lambda)/du)
        lv(l) = czero
        lu(l+nrzt) = czero
        lv(l+nrzt) = czero
      enddo
      ! note that lv coming out of totzsp is actually -d(lambda)/dv
      call totzsp(gc, gc(1+mns), gc(1+2*mns), gc(1+3*mns), gc(1+4*mns), gc(1+5*mns), &
                  r1, ru,        rv,          z1,          zu,          zv,          &
                  lu, lv,        rcon,        zcon, &
                  worka,worka,worka)

      ! COMPUTE CONSTRAINT FORCE (GCON)
      do l = 1,nrzt
        ru0(l)  = ru  (l) + ru  (l+nrzt)*sqrts(l)
        zu0(l)  = zu  (l) + zu  (l+nrzt)*sqrts(l)
        rcon(l) = rcon(l) + rcon(l+nrzt)*sqrts(l)
        zcon(l) = zcon(l) + zcon(l+nrzt)*sqrts(l)

        ! init gcon to zero
        gcon(l) = czero
        ! azmn: temporary storage for not-yet-dealiased gcon
        ! rcon0, zcon0 are initialized to 0 in vsetup --> ok to use in first iteration
        azmn(l) = (rcon(l)-rcon0(l))*ru0(l) + (zcon(l)-zcon0(l))*zu0(l)
      enddo

      do l = 1,2*mns
        gc(l) = czero
      enddo

      ! first iteration: save rcon into rcon0, zcon into zcon0
      if (iter2.eq.1 .and. nvac.eq.0) then
        call dcopy(nrzt, rcon, 1, rcon0, 1)
        call dcopy(nrzt, zcon, 1, zcon0, 1)
      endif

      ! only put non-zero values into gcon if iter2>1
      ! --> in the first iteration, gcon stays zero
      ! --> no constraint force in first iteration
      if (iter2.gt.1) &
        call alias(gcon, azmn, worka,worka,worka, gc, gc(1+mns))

      ! COMPUTE S AND THETA DERIVATIVE OF R AND Z AND JACOBIAN ON HALF-GRID
      call jacobian(r1, ru, z1, zu,                            &
      !             zu12        ru12        zs          rs
                    armn,       azmn,       brmn,       bzmn,  &
                    azmn(lodd), armn(lodd), brmn(lodd)        )
      !             gsqrt       r12         tau
      if (irst.eq.2 .and. iequi.eq.0) then
        ! Jacobian changes sign somewhere in volume
        ! --> flux surfaces cross, need to restart with smaller timestep
        return
      end if

      ! COMPUTE PRESSURE AND VOLUME ON HALF-GRID
      ! see above for re-use of arrays; additionally:
      !             gsqrt       bsq
      call pressure(azmn(lodd), bzmn(lodd), wint)
      if (iter2.eq.1) then
        voli = (twopi**2)*hs*sum(vp(2:ns))
      end if

      ! COMPUTE COVARIANT COMPONENTS OF B, MAGNETIC PRESSURE,
      ! AND METRIC ELEMENTS ON HALF-GRID
      do lk = 1,nznt
        rax(lk) = r1(1+ns*(lk-1))
        zax(lk) = z1(1+ns*(lk-1))
      enddo

      ! backup z1 in gco (since z1 gets overwritten in bcovar for iequi=1?)
      if (iequi.eq.1) &
        call dcopy(nrzt, z1, 1, gcon, 1)

      ! see above for re-use of arrays; additionally:
      !           bsubu  bsubv  gsqrt       bsq         r12         rs    zs
      call bcovar(clmn,  blmn,  azmn(lodd), bzmn(lodd), armn(lodd), bzmn, brmn, &
                  azmn,  armn,  guu,        guv,        gvv,        brmn(lodd), lu, lv)
      !           ru12   zu12                                       phipog

      ! restore z1 from gcon (since z1 got overwritten in bcovar for iequi=1 ?)
      if (iequi.eq.1) &
        call dcopy(nrzt, gcon, 1, z1, 1)

      ! bsubv@LCFS --> <B_zeta> is being computed here
      bz0 = ddot(nznt, blmn(ns), ns, wint(ns), ns)

      ! COMPUTE VACUUM MAGNETIC PRESSURE AT PLASMA EDGE
      ! NOTE: FOR FREE BOUNDARY RUNS, THE PLASMA VOLUME CAN
      ! BE INCREASED BY DECREASING CURPOL AT FIXED PHIPS.  THE
      ! VALUE FOR CURPOL GIVEN BELOW SHOULD BE USED AS AN
      ! INITIAL GUESS, WHICH CAN BE CHANGED BY READING IN
      ! CURPOL FROM DATA STATEMENT
      if (nvac.ne.0.and.iter2.gt.1) then
        call second(t4)
        timeon = real(t4)

        ivac2=mod(iter2-iter1, nvacskip)

        if ((fsqr+fsqz).le.1.e-2) &
          ivac=ivac+1

        if (ivac.eq.1) &
          ivac2 = 0

        rbtor = curpol

        if (ivac.eq.1 .and. rbtor.eq.0.) &
          rbtor = twopi*dnorm*bz0

        if (ivac.ge.1) then
          ctor = isigng*twopi*dnorm*(                          &
                   c1p5*ddot(nznt, clmn(ns  ), ns, wint(ns), ns)   &
                 - cp5 *ddot(nznt, clmn(ns-1), ns, wint(ns), ns) )

          call convert(rmnc, zmns, lmns, xm, xn, ns,          &
                       xc,          xc(1+  mns), xc(1+2*mns), &
                       xc(1+3*mns), xc(1+4*mns), xc(1+5*mns) )

          call vacuum(rmnc, zmns, xm, xn, ctor, rbtor, bsqvac, rax, zax)
        endif

        do lk = 1, nznt
          !                   bsq
          bsqsav(lk,3) = c1p5*bzmn(ns*lk+nrzt) - cp5*bzmn(ns*lk-1+nrzt)
          rbsq(lk)     = bsqvac(lk) * ohs*(r1(ns*lk) + r1(ns*lk+nrzt))
          dbsq(lk)     = abs(bsqvac(lk)-bsqsav(lk,3))
        enddo

        if(ivac.eq.1) then
          ! copy initial plasma and vacuum field into bsqsav(:,1:2)
          call dcopy(nznt, bzmn(ns+nrzt:ns+2*nrzt), ns, bsqsav(1,1), 1)
          call dcopy(nznt, bsqvac,                   1, bsqsav(1,2), 1)
        endif

        call second(t4)
        timeoff = real(t4)
        timer(1) = timer(1) + (timeoff-timeon)

      endif ! vacuum/free-boundary contribution

      if(iequi.eq.1)then
        ! COMPUTE REMAINING COVARIANT COMPONENT OF B (BSUBS),
        ! CYLINDRICAL COMPONENTS OF B (BR, BPHI, BZ), AND
        ! AVERAGE EQUILIBRIUM PROPERTIES AT END OF RUN

        !        r12         rs    zs    ru12  zu12  bsubs
        call bss(armn(lodd), bzmn, brmn, azmn, armn, crmn(lodd), &
                 lu, lv, rcon, czmn(lodd), zcon)
        !                br    bphi        bz

        ! [bcovar] bsubu bsubv bsq
        ! [eqfor]  bu    bv    bsq         rmag  zmag
        call eqfor(clmn, blmn, bzmn(lodd), xc,   xc(1+2*mns))

        !          bsq         gsqrt       bsubu bsubv
        call wrout(bzmn(lodd), azmn(lodd), clmn, blmn,        &
                   crmn(lodd), rcon, czmn(lodd), zcon, lu, lv)
        !          bsubs       br    bphi        bz

        return
      endif ! output file quantities

      ! COMPUTE MHD FORCES ON INTEGER-MESH
      call forces(guu, guv, gvv)

      ! FOURIER-TRANSFORM MHD FORCES TO (M,N)-SPACE
      do l = 1, neqs
        gc(l) = czero
      enddo
      call tomnsp(gc, gc(1+mns), gc(1+2*mns), gc(1+3*mns), gc(1+4*mns), gc(1+5*mns), &
                  armn, brmn, crmn, azmn, bzmn, czmn, blmn, clmn, rcon, zcon, workb,workb,workb)
      do l = 1, neqs
        ! TODO: Can one maybe omit this if the odd-m force components are not scaled by sqrt(s) ?
        ! --> comment in newer VMEC versions; are the odd-m force components different there?
        ! The hnr22_kis test case converges in (3 ... 5) % less iterations if this is commented out!
        gc(l) = gc(l) * scalxc(l)
      enddo

      ! COMPUTE FORCE RESIDUALS
      call residue(gc, gc(1+2*mns), gc(1+4*mns), bz0, workb)

      return

end


! { dg-do run }
! { dg-require-effective-target tls_runtime }
!$ use omp_lib

  integer, save, allocatable :: a(:, :)
  logical :: l
!$omp threadprivate (a)
  if (allocated (a)) call abort
  l = .false.
!$omp parallel copyin (a) num_threads (4) reduction(.or.:l)
  l = l.or.allocated (a)
!$omp end parallel
  if (l.or.allocated (a)) call abort
end


text:
C$Procedure      PCKGSS ( PCK get segment summary )

      SUBROUTINE PCKGSS ( SEGID,  SEGBOD, SEGFRM, SEGTYP,
     .                    SEGBTM, SEGETM, SEGBAD, SEGEAD )

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

      CHARACTER*(*)         SEGID
      INTEGER               SEGBOD
      INTEGER               SEGFRM
      INTEGER               SEGTYP
      DOUBLE PRECISION      SEGBTM
      DOUBLE PRECISION      SEGETM
      INTEGER               SEGBAD
      INTEGER               SEGEAD

C$ Brief_I/O
C
C     Variable  I/O  Description
C     --------  ---  --------------------------------------------------
C      SEGID     O   Segment ID for the segment in a PCK file.
C      SEGBOD    O   Body for the segment in a PCK file.
C      SEGFRM    O   Reference frame for the segment in a PCK file.
C      SEGTYP    O   Ephemeris type for the segment in a PCK file.
C      SEGBTM    O   Begin time (ET) for the segment in a PCK file.
C      SEGETM    O   End time (ET) for the segment in a PCK file.
C      SEGBAD    O   Begin address in the PCK file of a PCK segment.
C      SEGEAD    O   End address in the PCK file of a PCK segment.
C
C$ Detailed_Input
C
C     None.
C
C$ Detailed_Output
C
C      SEGID    Segment ID for the current segment in a PCK file.
C
C      SEGBOD   Body for the current segment in a PCK file. This is the
C               NAIF integer code for the body.
C
C      SEGFRM   Inertial reference frame for the current segment in
C               a PCK file. this is the NAIF integer code for the
C               inertial reference frame.
C
C      SEGTYP   Data type for the current segment in a PCK file.
C               This is an integer code which represents the PCK segment
C               data type.
C
C      SEGBTM   Begin time (ET) for the current segment in a PCK file.
C
C      SEGETM   End time (ET) for the current segment in a PCK file.
C
C      SEGBAD   Beginning DAF address for the data of the current
C               segment in a PCK file.
C
C      SEGEAD   Ending DAF address for the data of the current segment
C               in a PCK file.
C
C$ Parameters
C
C     None.
C
C$ Exceptions
C
C     1) If any errors occur, they will be signalled by routines called
C        by this routine.
C
C$ Files
C
C     None.
C
C$ Particulars
C
C     This routine provides the summary of the current segment in the
C     current PCK file. A PCK file must have been opened, and a
C     search must have been initiated before this routine may called.
C
C$ Examples
C
C     None.
C
C$ Restrictions
C
C     None.
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     K.R. Gehringer (JPL)
C
C$ Version
C
C-    Beta Version 1.0.0, 22-JAN-1993 (KRG)
C
C-&

C$ Index_Entries
C
C      get the summary of the current pck segment
C
C-&

C
C     SPICELIB functions
C
      LOGICAL               RETURN
      LOGICAL               FAILED
C
C     Local parameters
C
C     Length of the double precision portion of an PCK segment summary.
C
      INTEGER               DSLEN
      PARAMETER           ( DSLEN = 2 )
C
C     Length of the integer portion of an PCK segment summary.
C
      INTEGER               ISLEN
      PARAMETER           ( ISLEN = 5 )
C
C     Length of a segment summary in an PCK file.
C
      INTEGER               SUMLEN
      PARAMETER           ( SUMLEN = DSLEN + (ISLEN + 1 ) / 2 )
C
C     Length of an PCK segment ID.
C
      INTEGER               IDLEN
      PARAMETER           ( IDLEN = 8 * SUMLEN )
C
C     Set up mnemonic names for the segment begin and end time indices.
C
      INTEGER               BTMIDX
      PARAMETER           ( BTMIDX = 1          )

      INTEGER               ETMIDX
      PARAMETER           ( ETMIDX = BTMIDX + 1 )
C
C     Set up mnemonic names for the target body, the center body, the
C     inertial reference frame, the ephemeris type, the begin address
C     of the segment, and the end address of the segment.
C
      INTEGER               BODIDX
      PARAMETER           ( BODIDX = 1          )

      INTEGER               FRMIDX
      PARAMETER           ( FRMIDX = BODIDX + 1 )

      INTEGER               TYPIDX
      PARAMETER           ( TYPIDX = FRMIDX + 1 )

      INTEGER               BADIDX
      PARAMETER           ( BADIDX = TYPIDX + 1 )

      INTEGER               EADIDX
      PARAMETER           ( EADIDX = BADIDX + 1 )
C
C     Local variables
C
      CHARACTER*(IDLEN)     TMPID

      DOUBLE PRECISION      DSUM  (DSLEN )
      DOUBLE PRECISION      SUMRY (SUMLEN)

      INTEGER               ISUM  (ISLEN)
C
C     Standard SPICE error handling.
C
      IF ( RETURN () ) THEN
         RETURN
      ELSE
         CALL CHKIN ( 'PCKGSS' )
      END IF
C
C     Get the segment ID, the segment summary, and unpack the segment
C     summary. If an error occurred during any of these operations,
C     checkout and return to the caller, as an appropriate error message
C     will already be set. We use a temporary variable to read the
C     segment ID so that we do not modify the input value unless we get
C     one.
C
      CALL DAFGN  ( TMPID                           )
      CALL DAFGS  ( SUMRY                           )
      CALL DAFUS  ( SUMRY, DSLEN, ISLEN, DSUM, ISUM )

      IF ( FAILED() ) THEN

         CALL CHKOUT ( 'PCKGSS' )
         RETURN

      END IF
C
C     Set the output values for the begin and end ET times for the
C     segment.
C
      SEGBTM = DSUM(BTMIDX)
      SEGETM = DSUM(ETMIDX)
C
C     Set the output values for: the segment ID,  the center body, the
C     inertial reference frame code, the segment data type code, the
C     begin address of the segment data, and the end address of the
C     segment data.
C
      SEGID  = TMPID
      SEGBOD = ISUM(BODIDX)
      SEGFRM = ISUM(FRMIDX)
      SEGTYP = ISUM(TYPIDX)
      SEGBAD = ISUM(BADIDX)
      SEGEAD = ISUM(EADIDX)

      CALL CHKOUT ( 'PCKGSS' )
      RETURN
      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
program arrays_10
  implicit none
  integer, pointer :: x(:)
  allocate(x(1))
  x = [ 1 ]
end program arrays_10

max_stars_repo_name: Thirumalai-Shaktivel/lfortran

max_issues_repo_name: Thirumalai-Shaktivel/lfortran

max_forks_repo_name: Thirumalai-Shaktivel/lfortran


text: Matthew Clancy (18421925) was a California pioneer who eventually settled in the Davis area.  After growing up in impoverished Ireland,he left for the United States at the age of 16 with the hope of improving his quality of life.  Having little initial success in the Eastern United States, he traveled westward, banking on the agricultural riches of California.  After some brief stints in San Francisco as a dairy worker, Clancy was at least able to find consistent work as a farmhand in Yolo County; at one point he even worked for Isaac Chiles.  At the age of 31, he began to run a farm on a lease, but he was not able to actually buy his own farmland until around the age of 54.  But at his peak ownership, he came to possess almost 400 acres near the outskirts of Davis while he continued to lease and work other lands at the same time.  He primarily grew grain and raised cattle.  With his wife, Elizabeth Rowan, and their seven daughters, Clancy finally moved into a stately Davis Alpha Chi Omega home at the corner of C and 2nd near the end of his life.  He passed away at the age of 82 and is buried in Davis Cemetery near several of his immediate family members.


!
!             Parallel Sparse BLAS  version 3.5
!   (C) Copyright 2006, 2010, 2015, 2017
!       Salvatore Filippone    Cranfield University
!       Alfredo Buttari        CNRS-IRIT, Toulouse
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
!   1. Redistributions of source code must retain the above copyright
!      notice, this list of conditions and the following disclaimer.
!   2. Redistributions in binary form must reproduce the above copyright
!      notice, this list of conditions, and the following disclaimer in the
!      documentation and/or other materials provided with the distribution.
!   3. The name of the PSBLAS group or the names of its contributors may
!      not be used to endorse or promote products derived from this
!      software without specific written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
! PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE PSBLAS GROUP OR ITS CONTRIBUTORS
! BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
! CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
! SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
! CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
!
integer function psi_exist_ovr_elem(ovr_elem, dim_list,elem_searched)
  use psb_const_mod
  !    PURPOSE:
  !    == = ====
  !
  !    If ELEM_SEARCHED exist in the list OVR_ELEM returns its position in
  !    the list, else returns -1
  !
  !
  !     INPUT
  !     == = ===
  !     OVRLAP_ELEMENT_D.: Contains for all overlap points belonging to
  !                        the current process:
  !                          1. overlap point index
  !                          2. Number of domains sharing that overlap point
  !                        the end is marked by a -1...............................
  !
  !    DIM_LIST..........: Dimension of list OVRLAP_ELEMENT_D
  !
  !    ELEM_SEARCHED.....:point's  Local index identifier to be searched.

  implicit none

  !     ....Scalars parameters....
  integer(psb_ipk_) :: dim_list,elem_searched
  !     ...array parameters....
  integer(psb_ipk_) :: ovr_elem(dim_list,*)

  !     ...local scalars....
  integer(psb_ipk_) :: i

  i=1
  do while ((i.le.dim_list).and.(ovr_elem(i,1).ne.elem_searched))
    i=i+1
  enddo
  if ((i.le.dim_list).and.(ovr_elem(i,1).eq.elem_searched)) then
    psi_exist_ovr_elem=i
  else
    psi_exist_ovr_elem=-1
  endif
end function psi_exist_ovr_elem