text:
C     LAST UPDATE 16/03/89
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      SUBROUTINE CPK
      IMPLICIT NONE
C
C Purpose: Determine the major peak parameters, integral, width, position
C          skewness etc.
C
      INCLUDE 'COMMON.FOR'
C
C Calls   5: APLOT  , DAWRT  , ERRMSG , GETHDR , GETVAL , OPNFIL
C            PEKDES
C Called by: OTOKO
C
C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
C Local variables:
C
      REAL VALUE(10),SUM,SIGI,SPMAX,RMEAN,SKEW1,SKEW2
      INTEGER ICLO,JOP,IRC,IFRAME,IHFMAX,IFRMAX,IMEM,NCHAN,NFRAME
      INTEGER KREC,ISPEC,LSPEC,MEM,IFFR,ILFR,INCR,IFINC,NVAL,NPOINT
      INTEGER L1,L2,L3,L4,IMAX,IMED,IMODE,IFWHM,IACF,IACL,I,J,K
      CHARACTER*80 HEAD1,HEAD2
      CHARACTER*13 HFNAM,OFNAM
C
C ICLO   : Close file indicator
C JOP    : Open file indicator
C IRC    : Return code
C IFRAME : Total nos. of frames
C IFRMAX : Nos. of frames in sequence
C IHFMAX : Nos. of header file in sequence
C KREC   : Output file record
C ISPEC  : First header file of sequence
C LSPEC  : Last header file in sequence
C MEM    : Positional or calibration data indicator
C IFFR   : First frame in sequence
C ILFR   : Last frame in sequence
C INCR   : Header file increment
C IFINC  : Frame increment
C HFNAM  : Header filename
C OFNAM  : Output header filename
C HEAD1  : Output header record 1
C HEAD2  : Output header record 2
C IMEM   : Output memory dataset
C NCHAN  : Nos. of data points in spectrum
C NFRAME : Nos. of time frames
C
      DATA  IMEM/1/ , KREC/12/
C
C-----------------------------------------------------------------------
10    ICLO=0
      JOP=0
      NPOINT=0
C
      CALL GETHDR (ITERM,IPRINT,IUNIT,HFNAM,ISPEC,LSPEC,INCR,MEM,IFFR,
     1             ILFR,IFINC,IHFMAX,IFRMAX,NCHAN,IRC)
      IF (IRC.NE.0) GOTO 999
      IFRAME=IHFMAX+IFRMAX-1
C
20    WRITE (IPRINT,1000)
      CALL GETVAL (ITERM,VALUE,NVAL,IRC)
      IF (IRC.EQ.1) GOTO 10
      IF (IRC.EQ.2) GOTO 20
      IF (NVAL.LT.4) THEN
         CALL ERRMSG ('Error: Insufficient input')
         GOTO 20
      ELSE
         L1=INT(VALUE(1))
         L2=INT(VALUE(2))
         L3=INT(VALUE(3))
         L4=INT(VALUE(4))
         IF ((L1.LT.1.OR.L1.GT.NCHAN).OR.
     1       (L2.LT.1.OR.L2.GT.NCHAN).OR.
     2       (L3.LT.1.OR.L3.GT.NCHAN).OR.
     3       (L4.LT.1.OR.L4.GT.NCHAN)) THEN
            CALL ERRMSG ('Error: Value outside of range')
            GOTO 20
         ENDIF
      ENDIF
      DO 60 I=1,IHFMAX
C
         CALL OPNFIL (JUNIT,HFNAM,ISPEC,MEM,IFFR,ILFR,NCHAN,NFRAME,IRC)
         IF (IRC.EQ.0) THEN
C
            ISPEC=ISPEC+INCR
            DO 50 J=1,IFRMAX
               READ (JUNIT,REC=IFFR) (SP1(K),K=1,NCHAN)
               IFFR=IFFR+IFINC
               NPOINT=NPOINT+1
C
               SUM=0.0
               DO 40 K=L2,L3
                  SUM=SUM+SP1(K)
40             CONTINUE
               CALL PEKDES (SP1,NCHAN,L1,L4,IACF,IACL,SPMAX,IMAX,RMEAN,
     1                      IMED,IMODE,SIGI,IFWHM,SKEW1,SKEW2,IRC)
               SP2(NPOINT)=SUM
               SP3(NPOINT)=SIGI
               SP4(NPOINT)=SPMAX
               SP5(NPOINT)=IMAX
               SP6(NPOINT)=RMEAN
               SP7(NPOINT)=IMED
               SP8(NPOINT)=IMODE
               SP9(NPOINT)=IFWHM
               SP10(NPOINT)=SKEW1
               SP11(NPOINT)=SKEW2
               SP12(NPOINT)=IACF
               SPW(NPOINT)=IACL
50          CONTINUE
            CLOSE (UNIT=JUNIT)
         ENDIF
60    CONTINUE
      IF (IFRAME.NE.1) THEN
         CALL APLOT (ITERM,IPRINT,SP2,IFRAME)
      ENDIF
      CALL OUTFIL (ITERM,IPRINT,OFNAM,HEAD1,HEAD2,IRC)
      IF (IRC.EQ.0) THEN
         CALL DAWRT (KUNIT,OFNAM,IMEM,IFRAME,12,HEAD1,HEAD2,SP2,
     &                  1,JOP,ICLO,IRC)
         IF (IRC.NE.0) GOTO 10
         CALL DAWRT (KUNIT,OFNAM,IMEM,IFRAME,12,HEAD1,HEAD2,SP3,
     &                  2,JOP,ICLO,IRC)
         IF (IRC.NE.0) GOTO 10
         CALL DAWRT (KUNIT,OFNAM,IMEM,IFRAME,12,HEAD1,HEAD2,SP4,
     &                  3,JOP,ICLO,IRC)
         IF (IRC.NE.0) GOTO 10
         CALL DAWRT (KUNIT,OFNAM,IMEM,IFRAME,12,HEAD1,HEAD2,SP5,
     &                  4,JOP,ICLO,IRC)
         IF (IRC.NE.0) GOTO 10
         CALL DAWRT (KUNIT,OFNAM,IMEM,IFRAME,12,HEAD1,HEAD2,SP6,
     &                  5,JOP,ICLO,IRC)
         IF (IRC.NE.0) GOTO 10
         CALL DAWRT (KUNIT,OFNAM,IMEM,IFRAME,12,HEAD1,HEAD2,SP7,
     &                  6,JOP,ICLO,IRC)
         IF (IRC.NE.0) GOTO 10
         CALL DAWRT (KUNIT,OFNAM,IMEM,IFRAME,12,HEAD1,HEAD2,SP8,
     &                  7,JOP,ICLO,IRC)
         IF (IRC.NE.0) GOTO 10
         CALL DAWRT (KUNIT,OFNAM,IMEM,IFRAME,12,HEAD1,HEAD2,SP9,
     &                  8,JOP,ICLO,IRC)
         IF (IRC.NE.0) GOTO 10
         CALL DAWRT (KUNIT,OFNAM,IMEM,IFRAME,12,HEAD1,HEAD2,SP10,
     &                  9,JOP,ICLO,IRC)
         IF (IRC.NE.0) GOTO 10
         CALL DAWRT (KUNIT,OFNAM,IMEM,IFRAME,12,HEAD1,HEAD2,SP11,
     &                  10,JOP,ICLO,IRC)
         IF (IRC.NE.0) GOTO 10
         CALL DAWRT (KUNIT,OFNAM,IMEM,IFRAME,12,HEAD1,HEAD2,SP12,
     &                  11,JOP,ICLO,IRC)
         IF (IRC.NE.0) GOTO 10
         ICLO=1
         CALL DAWRT (KUNIT,OFNAM,IMEM,IFRAME,12,HEAD1,HEAD2,SPW,
     &                  12,JOP,ICLO,IRC)
         IF (IRC.NE.0) GOTO 10
      ENDIF
      GOTO 10
999   RETURN
C
1000  FORMAT (' Enter range of points on either side of peak,',/
     1        ' low_min, low_max, high_min, high_max: ',$)
1010  FORMAT (' Frame',I4,' Pos of max',I5,' Max',G12.5,' Mean',G12.5,/,
     1        ' Median',I5,' Mode',I5,' Std Dev.',G12.5,' FWHM',I5,/,
     2        ' Skew 1',G12.5,' Skew 2',G12.5,' First channel',I5,/,
     3        ' Last channel',I5,' Integral',G12.5)
      END


text:
locate type T at "T.fpp"
locate constant a at "a.fpp"
locate constant b at "b.fpp"
locate constant c at "c.fpp"
locate constant d at "d.fpp"

active component C {

   param P: T default a id b set opcode c save opcode d

}

ext: fpp

max_stars_repo_name: kevin-f-ortega/fpp

max_issues_repo_name: kevin-f-ortega/fpp

max_forks_repo_name: kevin-f-ortega/fpp


program tarefaa
    print *, "Digite o valor inteiro de N:"
    read (*,*) N

    ! Número de "Andarilhos" M
    M = 1000

    ! loop do passo
    do i = 1, N
        ! zera a soma
        soma = 0e0

        ! loop do andarilho
        do j = 1, M
            ! soma os passos de cada andarilho
            soma = soma + rand()**i
        end do

        ! Escreve a media no terminal
        write (*,'(A,I0,A,1F5.3)') "<x^",i,">: ", soma/M
    end do

end program tarefaa


text:
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      PROGRAM PROJ8
C     Project 8: Monte Carlo solution of the H2 molecule
C  COMPUTATIONAL PHYSICS (FORTRAN VERSION)
C  by Steven E. Koonin and Dawn C. Meredith
C  Copyright 1989, Addison-Wesley Publishing Company Inc.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      CALL INIT           !display header screen, setup parameters
5     CONTINUE            !main loop/ execute once for each set of param
        CALL PARAM        !get input from screen
        CALL ARCHON       !calculate the eigenvalue for this value of S
      GOTO 5
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE ARCHON
C calculates the electronic eigenvalue or energy auto-correlation
C for a given separation of the protons
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'PARAM.P8'
      INCLUDE 'IO.ALL'
C Local variables:
C     energy has two indices
C     first index is the level: sweep, group, or total
C     second index is the value: quantity, quant**2, or sigma**2
      REAL ENERGY(3,3)           !energy
      REAL CONFIG(NCOORD)        !configuration
      REAL W                     !weight for single variational config
      REAL WEIGHT(MAXENS)        !weight of ensemble members
      REAL ENSMBL(NCOORD,MAXENS) !ensemble of configurations
      REAL ESAVE(MAXCRR)         !array of local energies for corr
      REAL EPSILN                !local energy of CONFIG
      REAL ACCPT                 !acceptance ratio
      INTEGER ITHERM             !thermalization index
      INTEGER ISWP,ISMPL         !sweep and sample index
      INTEGER IQUANT             !quantity index
      INTEGER IGRP               !group index
      INTEGER NLINES             !number of lines printed to terminal
      INTEGER MORE               !how many more groups
      INTEGER SWEEP,GROUP,TOTAL  !which level of calculation
      INTEGER VALUE,SQUARE,SIGSQ !which quantity
      INTEGER CORR,EPS           !what is being calculated?
      INTEGER PIMC,VARY          !which method?
C Functions:
      REAL ELOCAL                !local energy
      INTEGER GETINT             !get integer data from screen
      DATA SWEEP,GROUP,TOTAL/1,2,3/
      DATA VALUE,SQUARE,SIGSQ/1,2,3/
      DATA EPS,CORR /1,2/
      DATA VARY,PIMC /1,2/
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     output summary of parameters
      IF (TTERM) CALL PRMOUT(OUNIT,NLINES)
      IF (TFILE) CALL PRMOUT(TUNIT,NLINES)
      IF (GFILE) CALL PRMOUT(GUNIT,NLINES)
C
C     generate initial configuration or ensemble of configurations
      IF (METHOD .EQ. VARY) THEN
         CALL INTCFG(CONFIG,W)
      ELSE IF (METHOD .EQ. PIMC) THEN
         CALL INTENS(ENSMBL,WEIGHT,CONFIG)
      END IF
C
C     take thermalization steps
      DO 10 ITHERM=1,NTHERM
        IF (ITHERM .EQ. 1) WRITE (OUNIT,*) ' Thermalizing...'
        IF (ITHERM .EQ. NTHERM) WRITE (OUNIT,*) ' '
        IF (METHOD .EQ. VARY) THEN
           CALL METROP(CONFIG,W,ACCPT)
        ELSE IF (METHOD .EQ. PIMC) THEN
           CALL TSTEP(ENSMBL,WEIGHT,EPSILN)
        END IF
10    CONTINUE
C
      DO 11 IQUANT=1,3                       !zero total sums
         ENERGY(TOTAL,IQUANT)=0.
11    CONTINUE
      ACCPT=0                                !zero acceptance
      MORE=NGROUP                            !initial number of groups
C
15    CONTINUE                               !allow for more groups
        DO 20 IGRP=NGROUP-MORE+1,NGROUP      !loop over groups
C
           DO 21 IQUANT=1,3                  !zero group sums
             ENERGY(GROUP,IQUANT)=0.
21         CONTINUE
C
           DO 30 ISWP=1,NFREQ*NSMPL          !loop over sweeps
C
             IF (METHOD .EQ. VARY) THEN      !take a Metrop step
                CALL METROP(CONFIG,W,ACCPT)
             ELSE IF (METHOD .EQ. PIMC) THEN !or a time step
                CALL TSTEP(ENSMBL,WEIGHT,EPSILN)
             END IF
C
             IF (MOD(ISWP,NFREQ) .EQ. 0) THEN !sometimes save the energy
                ISMPL=ISWP/NFREQ
                IF (METHOD .EQ. VARY) EPSILN=ELOCAL(CONFIG)
                ENERGY(GROUP,VALUE)=ENERGY(GROUP,VALUE)+EPSILN
                ENERGY(GROUP,SQUARE)=ENERGY(GROUP,SQUARE)+EPSILN**2
                IF (.NOT. TERSE) THEN
                  IF (TTERM) WRITE (OUNIT,100) IGRP,ISMPL,NSMPL,EPSILN
                  IF (TFILE) WRITE (TUNIT,100) IGRP,ISMPL,NSMPL,EPSILN
100               FORMAT (5X,' Group ',I4, ',  sample ',I4,' of ',I4,5X,
     +                   'Energy =',F9.4)
                END IF
             END IF
C
30        CONTINUE                            !this group is done
          IF (CALC .EQ. CORR) THEN            !save energy for corr
              ESAVE(IGRP)=ENERGY(GROUP,VALUE)
          ELSE                                !or calc averages
              CALL AVERAG(ENERGY,ACCPT,IGRP)
          END IF
C
20      CONTINUE
        IF (CALC .EQ. CORR) CALL CRLTNS(ESAVE)   !calc corr
C
C       allow for more groups, taking care not to exceed array bounds
        MORE=GETINT(10,0,1000,'How many more groups?')
        IF ((CALC .EQ. CORR) .AND. (NGROUP+MORE .GT. MAXCRR)) THEN
            WRITE (OUNIT,200) MAXCRR-NGROUP
200         FORMAT(' You will run out of storage space for '
     +      'corr if you do more than ',I3,' more groups')
            MORE=GETINT(MAXCRR-NGROUP,0,MAXCRR-NGROUP,
     +                    'How many more groups?')
        END IF
      IF (MORE .GT. 0) THEN
           NGROUP=NGROUP+MORE
           NLINES=0
           IF (TTERM) CALL CLEAR
           GOTO 15
      END IF
C
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE TSTEP(ENSMBL,WEIGHT,EPSILN)
C take a time step using the Path Integral Monte Carlo method
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'PARAM.P8'
C Input/Output variables
      REAL WEIGHT(MAXENS)       !weight of ensemble members (I/O)
      REAL ENSMBL(NCOORD,MAXENS)!ensemble of configurations (I/O)
      REAL EPSILN               !local energy of CONFIG (output)
C Local variables:
      REAL CONFIG(NCOORD)       !configuration
      REAL W                    !weight for single config
      REAL EBAR,WBAR            !ensmble average local energy and weight
      INTEGER IENSEM            !ensemble index
      INTEGER ICOORD            !coordinate index
      REAL NORM                 !normalization of weights
      REAL SHIFT(NCOORD)        !array containing drift vector
C Functions:
      REAL GAUSS                !Gaussian random number
      REAL ELOCAL               !local energy of the configuration
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      EBAR=0.                    !zero sums
      WBAR=0.
      DO 10 IENSEM=1,NENSEM                        !loop over ensemble
         DO 20 ICOORD=1,NCOORD
            CONFIG(ICOORD)=ENSMBL(ICOORD,IENSEM)   !get a configuration
20       CONTINUE
         CALL DRIFT(CONFIG,SHIFT)                  !calc shifts
         DO 30 ICOORD=1,NCOORD
            CONFIG(ICOORD)=CONFIG(ICOORD)+         !shift configuration
     +                     GAUSS(DSEED)*SQHBDT+SHIFT(ICOORD)
30       CONTINUE
C
         EPSILN=ELOCAL(CONFIG)                     !calculate energy
         WEIGHT(IENSEM)=WEIGHT(IENSEM)*EXP(-EPSILN*DT)  !calc weight
         EBAR=EBAR+WEIGHT(IENSEM)*EPSILN           !update sums
         WBAR=WBAR+WEIGHT(IENSEM)
C
         DO 40 ICOORD=1,NCOORD
            ENSMBL(ICOORD,IENSEM)=CONFIG(ICOORD)   !save configuration
40       CONTINUE
10    CONTINUE
C
      EPSILN=EBAR/WBAR                    !weighted average energy
      NORM=NENSEM/WBAR
      DO 50 IENSEM=1,NENSEM               !renormalize weights
         WEIGHT(IENSEM)=NORM*WEIGHT(IENSEM)
50    CONTINUE
C
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE METROP(CONFIG,W,ACCPT)
C take a Metropolis step
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'PARAM.P8'
C Input/Output variables:
      REAL CONFIG(NCOORD)        !configuration
      REAL W                     !weight for single config
      REAL ACCPT                 !acceptance ratio
C Local variables:
      INTEGER ICOORD             !coordinate index
      REAL CSAVE(NCOORD)         !temp storage for last config
      REAL WTRY                  !weight for trial config
C Function:
      REAL PHI                   !total wave function
      REAL RANNOS                !uniform random number
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      DO 10 ICOORD=1,NCOORD
       CSAVE(ICOORD)=CONFIG(ICOORD)        !save previous values
       CONFIG(ICOORD)=CONFIG(ICOORD)+DELTA*(RANNOS(DSEED)-.5)!trial step
10    CONTINUE
      WTRY=PHI(CONFIG)**2                  !trial weight
C
      IF (WTRY/W .GT. RANNOS(DSEED)) THEN  !sometimes accept the step
         W=WTRY                            !save new weight
         ACCPT=ACCPT+1                     !update accpt ratio
      ELSE
        DO 20 ICOORD=1,NCOORD
          CONFIG(ICOORD)=CSAVE(ICOORD)     !or else restore old config
20      CONTINUE
      END IF
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      REAL FUNCTION ELOCAL(CONFIG)
C calculate the local energy for CONFIG
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variable:
      INCLUDE 'PARAM.P8'
C Input variables:
      REAL CONFIG(NCOORD)        !configuration
C Local variables:
      REAL TPOP,VPOP             !kinetic and potential contributions
      REAL EECORR                !elec-elec correlation
      REAL CROSS,CROSS1,CROSS2   !cross terms
      REAL ONEE1,ONEE2           !one electron terms
      REAL X1,X2,Y1,Y2,Z1,Z2     !coordinates of 2 electrons
      REAL R1L,R1R,R2L,R2R,R12   !relative distances
      REAL CHI1,CHI2,F           !parts of the wave function
      REAL DOTR1L,DOTR2L,DOTR1R,DOTR2R !dot products with R12
      REAL R12DR1,SR12Z          !temp vars for dot products
      REAL CHI,FDCHI,SDCHI,LAPCHI !atomic orbitals
      REAL FEE,FDFEE,SDFEE,LAPFEE !elec-elec correlations
      REAL DIST                   !Euclidean distance
      REAL R,X,Y,Z                !dummy variables
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     define functions
      CHI(R)=EXP(-R/A)                        !atomic orbital
      FDCHI(R)=-CHI(R)/A                      !its first derivative,
      SDCHI(R)=CHI(R)/A/A                     !second derivative,
      LAPCHI(R)=SDCHI(R)+2*FDCHI(R)/R         !and Laplacian
C
      FEE(R)=EXP(R/(ALPHA*(1+BETA*R)))        !elec-elec correlation
      FDFEE(R)=FEE(R)/(ALPHA*(1.+BETA*R)**2)  !its first,second deriv,
      SDFEE(R)=FDFEE(R)**2/FEE(R)-2.*BETA*FEE(R)/ALPHA/(1+BETA*R)**3
      LAPFEE(R)=SDFEE(R)+2*FDFEE(R)/R         !and Laplacian
C
      DIST(X,Y,Z)=SQRT(X**2+Y**2+Z**2)        !Euclidean distance
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     get coordinates and radii
      CALL RADII(X1,X2,Y1,Y2,Z1,Z2,R1L,R2L,R1R,R2R,R12,CONFIG)
C
C     calculate dot products with R12
      R12DR1=X1*(X1-X2)+Y1*(Y1-Y2)+Z1*(Z1-Z2) !convenient starting place
      SR12Z=S*(Z1-Z2)/2                !useful constant
      DOTR1L=R12DR1+SR12Z              !dot products with R12
      DOTR1R=R12DR1-SR12Z
      DOTR2L=DOTR1L-R12**2
      DOTR2R=DOTR1R-R12**2
      DOTR1L=DOTR1L/R12/R1L            !dot products of unit vectors
      DOTR2L=DOTR2L/R12/R2L
      DOTR1R=DOTR1R/R12/R1R
      DOTR2R=DOTR2R/R12/R2R
C
      CHI1=CHI(R1R)+CHI(R1L)           !pieces of the total wave function
      CHI2=CHI(R2R)+CHI(R2L)
      F=FEE(R12)
C
      EECORR=2*LAPFEE(R12)/F           !correlation contribution
      ONEE1=(LAPCHI(R1L)+LAPCHI(R1R))/CHI1               !electron one
      ONEE2=(LAPCHI(R2L)+LAPCHI(R2R))/CHI2               !electron two
      CROSS1=(FDCHI(R1L)*DOTR1L+FDCHI(R1R)*DOTR1R)/CHI1  !cross terms
      CROSS2=(FDCHI(R2L)*DOTR2L+FDCHI(R2R)*DOTR2R)/CHI2
      CROSS=2*FDFEE(R12)*(CROSS1-CROSS2)/F
C
      TPOP=-HBM*(EECORR+ONEE1+ONEE2+CROSS)/2                !kinetic
      VPOP=-E2*(1./R1L + 1./R1R + 1./R2L + 1./R2R - 1./R12) !potential
      ELOCAL=TPOP+VPOP                                      !total
C
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE DRIFT(CONFIG,SHIFT)
C calculate the drift vector (SHIFT) for CONFIG
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'PARAM.P8'
C Input/Output variables:
      REAL CONFIG(NCOORD)        !configuration (input)
      REAL SHIFT(NCOORD)         !array containing drift vector (output)
C Local variables:
      INTEGER ICOORD              !coordinate index
      REAL X1,X2,Y1,Y2,Z1,Z2      !coordinates of 2 electrons
      REAL R1L,R1R,R2L,R2R,R12    !relative distances
      REAL CHI1,CHI2,F            !parts of the wave function
      REAL CHI,FDCHI,SDCHI,LAPCHI !atomic orbital
      REAL FEE,FDFEE,SDFEE,LAPFEE !elec-elec correlations
      REAL R                      !dummy variables
      REAL FACTA,FACTB,FACTE      !useful factors
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     define functions
      CHI(R)=EXP(-R/A)                        !atomic orbital
      FDCHI(R)=-CHI(R)/A                      !its first derivative,
      SDCHI(R)=CHI(R)/A/A                     !second derivative,
      LAPCHI(R)=SDCHI(R)+2*FDCHI(R)/R         !and Laplacian
C
      FEE(R)=EXP(R/(ALPHA*(1+BETA*R)))        !elec-elec correlation
      FDFEE(R)=FEE(R)/(ALPHA*(1.+BETA*R)**2)  !its first, second deriv,
      SDFEE(R)=FDFEE(R)**2/FEE(R)-2.*BETA*FEE(R)/ALPHA/(1+BETA*R)**3
      LAPFEE(R)=SDFEE(R)+2*FDFEE(R)/R         !and Laplacian
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     get coordinates and radii
      CALL RADII(X1,X2,Y1,Y2,Z1,Z2,R1L,R2L,R1R,R2R,R12,CONFIG)
C
      CHI1=CHI(R1R)+CHI(R1L)           !pieces of the total wave function
      CHI2=CHI(R2R)+CHI(R2L)
      F=FEE(R12)
C
      FACTA=HBMDT*(FDCHI(R1L)/R1L+FDCHI(R1R)/R1R)/CHI1  !useful factors
      FACTB=HBMDT*(FDCHI(R1L)/R1L-FDCHI(R1R)/R1R)/CHI1
      FACTE=HBMDT*FDFEE(R12)/F/R12
C
      SHIFT(1)=FACTA*X1+FACTE*(X1-X2)  !shift for electron one
      SHIFT(2)=FACTA*Y1+FACTE*(Y1-Y2)
      SHIFT(3)=FACTA*Z1+FACTE*(Z1-Z2)+FACTB*S2
C
      FACTA=HBMDT*(FDCHI(R2L)/R2L+FDCHI(R2R)/R2R)/CHI2
      FACTB=HBMDT*(FDCHI(R2L)/R2L-FDCHI(R2R)/R2R)/CHI2
C
      SHIFT(4)=FACTA*X2-FACTE*(X1-X2)  !shift for electron two
      SHIFT(5)=FACTA*Y2-FACTE*(Y1-Y2)
      SHIFT(6)=FACTA*Z2-FACTE*(Z1-Z2)+FACTB*S2
C
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      REAL FUNCTION GAUSS(DSEED)
C returns a Gaussian random number with zero mean and unit variance
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      INTEGER IGAUSS                 !sum index
      DOUBLE PRECISION DSEED         !random number seed
      REAL RANNOS                    !uniform random number
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      GAUSS=0.                       !sum 12 uniform random numbers
      DO 10 IGAUSS=1,12
         GAUSS=GAUSS+RANNOS(DSEED)
10    CONTINUE
      GAUSS=GAUSS-6.                 !subtract six so that mean=0
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      REAL FUNCTION PHI(CONFIG)
C calculates the total variational wave function for CONFIG
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'PARAM.P8'
C Input variables:
      REAL CONFIG(NCOORD)        !configuration
C Local variables:
      REAL X1,X2,Y1,Y2,Z1,Z2     !coordinates of 2 electrons
      REAL R1L,R1R,R2L,R2R,R12   !relative distances
      REAL CHI1R,CHI1L,CHI2R,CHI2L,F  !parts of the wave function
      REAL CHI,FEE               !terms in the wave function
      REAL R                     !radius
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      CHI(R)=EXP(-R/A)                   !atomic orbital
      FEE(R)=EXP(R/(ALPHA*(1+BETA*R)))   !electron-electron correlation
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     calculate the radii
      CALL RADII(X1,X2,Y1,Y2,Z1,Z2,R1L,R2L,R1R,R2R,R12,CONFIG)
C
      CHI1R=CHI(R1R)          !pieces of the total wave function
      CHI1L=CHI(R1L)
      CHI2R=CHI(R2R)
      CHI2L=CHI(R2L)
      F=FEE(R12)
      PHI=(CHI1L +CHI1R)*(CHI2L+CHI2R)*F   !the whole thing
C
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE INTENS(ENSMBL,WEIGHT,CONFIG)
C generate the ENSMBL at t=0 for PIMC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'PARAM.P8'
C Output variables:
      REAL CONFIG(NCOORD)        !configuration
      REAL WEIGHT(MAXENS)        !weight of ensemble members
      REAL ENSMBL(NCOORD,MAXENS) !ensemble of configurations
C Local variables:
      INTEGER ISTEP              !step index
      INTEGER ICOORD             !coordinate index
      REAL W                     !weight for single config
      REAL ACCPT                 !acceptance ratio
      INTEGER IENSEM             !ensemble index
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      CALL INTCFG(CONFIG,W)      !generate a single intial configuration
C
      DO 10 ISTEP=1,20                 !do 20 thermalization steps
         CALL METROP(CONFIG,W,ACCPT)   !using Metropolis algorithm
10    CONTINUE
C
      DO 30 ISTEP=1,10*NENSEM          !generate the ensemble
         CALL METROP(CONFIG,W,ACCPT)   !take a Metrop step
         IF (MOD(ISTEP,10) .EQ. 0) THEN
             IENSEM=ISTEP/10           !save every 10th config
             DO 20 ICOORD=1,NCOORD
                ENSMBL(ICOORD,IENSEM)=CONFIG(ICOORD)
20           CONTINUE
             WEIGHT(IENSEM)=1.         !set all weights=1
         END IF
30    CONTINUE
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE INTCFG(CONFIG,W)
C generate a configuration (CONFIG) and calculate its weight (W)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variable:
      INCLUDE 'PARAM.P8'
C Output variables:
      REAL CONFIG(NCOORD)        !configuration
      REAL W                     !weight for single config
C Local variables:
      INTEGER ICOORD             !coordinate index
C Function:
      REAL PHI                   !total wave function
      REAL RANNOS                !uniform random number
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      DO 10 ICOORD=1,NCOORD      !pick configuration at random
         CONFIG(ICOORD)=A*(RANNOS(DSEED)-.5)
10    CONTINUE
      CONFIG(3)=CONFIG(3)+S2    !center elec 1. at right
      CONFIG(6)=CONFIG(6)-S2    !center elec 2. at left
      W=PHI(CONFIG)**2          !weight=phi**2
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE CRLTNS(ESAVE)
C calculate the energy auto-correlations
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'PARAM .P8'
      INCLUDE 'IO.ALL'
C Input variables:
      REAL ESAVE(MAXCRR)           !array of local energies for corr.
C Local variables:
      REAL EI,EIK,ESQI,ESQIK,EIEK  !sums
      INTEGER I,K                  !index of ESAVE
      REAL ECORR(0:MAXCRR)         !energy auto-correlations
      INTEGER NI                   !number of energies in sum
      INTEGER SCREEN               !send to terminal
      INTEGER PAPER                !make a hardcopy
      INTEGER FILE                 !send to a file
      DATA SCREEN,PAPER,FILE/1,2,3/
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      DO 10 K=0,NCORR              !loop over correlation lengths
         EI=0.                     !zero sums
         EIK=0.
         ESQI=0.
         ESQIK=0.
         EIEK=0.
         NI=NGROUP-K
         DO 20 I=1,NI
            EI=EI+ESAVE(I)         !calculate sums
            EIK=EIK+ESAVE(I+K)
            ESQI=ESQI+ESAVE(I)**2
            ESQIK=ESQIK+ESAVE(I+K)**2
            EIEK=EIEK+ESAVE(I)*ESAVE(I+K)
20       CONTINUE
         EI=EI/NI                  !calculate averages
         EIK=EIK/NI
         ESQI=ESQI/NI
         ESQIK=ESQIK/NI
         EIEK=EIEK/NI
         ECORR(K)=(EIEK-EI*EIK)/(SQRT(ESQI-EI**2))/(SQRT(ESQIK-EIK**2))
10    CONTINUE
C
      IF (GTERM) THEN                   !display results
          CALL PAUSE ('to see the energy auto-correlations...',1)
          CALL GRFOUT(SCREEN,ECORR)
      ELSE IF (TTERM) THEN
          CALL PAUSE ('to see the energy auto-correlations...',1)
          CALL CRROUT(OUNIT,ECORR)
      END IF
      IF (TFILE) CALL CRROUT(TUNIT,ECORR)
      IF (GHRDCP) CALL GRFOUT(PAPER,ECORR)
      IF (GFILE) CALL GRFOUT(FILE,ECORR)
C
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE AVERAG(ENERGY,ACCPT,IGRP)
C calculate group averages, add to totals, print out
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'PARAM.P8'
      INCLUDE 'IO.ALL'
C Input variables:
C     energy has two indices
C     first index is the level: sweep, group, or total
C     second index is the value: quantity, quant**2, or sigma**2
      REAL ENERGY(3,3)           !energy
      INTEGER IGRP               !group index
      REAL ACCPT                 !acceptance ratio
C Local variables:
      REAL EVALUE                !current average energy
      REAL SIG1,SIG2             !uncertainties in energy
      REAL U                     !total pot energy of the system
      INTEGER NLINES             !number of lines printed to terminal
      INTEGER SWEEP,GROUP,TOTAL  !which level of calculation
      INTEGER VALUE,SQUARE,SIGSQ !which quantity
      DATA SWEEP,GROUP,TOTAL/1,2,3/
      DATA VALUE,SQUARE,SIGSQ/1,2,3/
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     calculate group averages and uncertainties
      ENERGY(GROUP,VALUE)=ENERGY(GROUP,VALUE)/NSMPL
      ENERGY(GROUP,SQUARE)=ENERGY(GROUP,SQUARE)/NSMPL
      ENERGY(GROUP,SIGSQ)=
     +  (ENERGY(GROUP,SQUARE)-ENERGY(GROUP,VALUE)**2)/NSMPL
      IF (ENERGY(GROUP,SIGSQ) .LT. 0.) ENERGY(GROUP,SIGSQ)=0.
C
C     add to totals
      ENERGY(TOTAL,VALUE)=ENERGY(TOTAL,VALUE)+ENERGY(GROUP,VALUE)
      ENERGY(TOTAL,SQUARE)=ENERGY(TOTAL,SQUARE)+ENERGY(GROUP,SQUARE)
      ENERGY(TOTAL,SIGSQ)=ENERGY(TOTAL,SIGSQ)+ENERGY(GROUP,SIGSQ)
C
C     calculate current grand averages
      EVALUE=ENERGY(TOTAL,VALUE)/IGRP
      SIG1=(ENERGY(TOTAL,SQUARE)/IGRP-EVALUE**2)/IGRP/NSMPL
      IF (SIG1 .LT. 0.) SIG1=0.
      SIG1=SQRT(SIG1)
      SIG2=SQRT(ENERGY(TOTAL,SIGSQ))/IGRP
C
C     calculate total energy of the system
      IF (S .GT. .01) THEN
          U=EVALUE+E2/S+E2/ABOHR
      ELSE
          U=0.
      END IF
C
      IF (TTERM) CALL TXTOUT(IGRP,ENERGY,EVALUE,SIG1,SIG2,U,ACCPT,OUNIT)
      IF (TFILE) CALL TXTOUT(IGRP,ENERGY,EVALUE,SIG1,SIG2,U,ACCPT,TUNIT)
C
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE RADII(X1,X2,Y1,Y2,Z1,Z2,R1L,R2L,R1R,R2R,R12,CONFIG)
C calculates cartesian coordinates and radii given CONFIG
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variable:
      INCLUDE 'PARAM.P8'
C Input variables:
      REAL CONFIG(NCOORD)        !configuration
C Output variables:
      REAL X1,X2,Y1,Y2,Z1,Z2     !coordinates of 2 electrons
      REAL R1L,R1R,R2L,R2R,R12   !relative distances
      REAL DIST                  !Euclidean distance
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      DIST(X1,Y1,Z1)=SQRT(X1**2+Y1**2+Z1**2)        !Euclidean distance
C
      X1=CONFIG(1)            !give the CONFIG elements their real names
      X2=CONFIG(4)
      Y1=CONFIG(2)
      Y2=CONFIG(5)
      Z1=CONFIG(3)
      Z2=CONFIG(6)
C
      R1L=DIST(X1,Y1,Z1+S2)   !calculate separations
      R1R=DIST(X1,Y1,Z1-S2)
      R2L=DIST(X2,Y2,Z2+S2)
      R2R=DIST(X2,Y2,Z2-S2)
      R12=DIST(X1-X2,Y1-Y2,Z1-Z2)
C
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE INIT
C initializes constants, displays header screen,
C initializes menu arrays for input parameters
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'IO.ALL'
      INCLUDE 'MENU.ALL'
      INCLUDE 'PARAM.P8'
C Local parameters:
      CHARACTER*80 DESCRP          !program description
      DIMENSION DESCRP(20)
      INTEGER NHEAD,NTEXT,NGRAPH   !number of lines for each description
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     get environment parameters
      CALL SETUP
C
C     display header screen
      DESCRP(1)= 'PROJECT 8'
      DESCRP(2)= 'Monte Carlo solution of the H2 molecule'
      NHEAD=2
C
C     text output description
      DESCRP(3)= 'electronic eigenvalue and its uncertainty'
      DESCRP(4)= 'or energy auto-correlation'
      NTEXT=2
C
C     graphics output description
      DESCRP(5)= 'energy auto-correlation vs. correlation length'
      NGRAPH=1
C
      CALL HEADER(DESCRP,NHEAD,NTEXT,NGRAPH)
C
C     calculate constants
      HBM=7.6359               !hbar*2/(mass)
      E2=14.409                !charge of the electron
      ABOHR=HBM/E2
C
C     setup menu arrays, beginning with constant part
      CALL MENU
C
      MTYPE(12)=TITLE
      MPRMPT(12)= 'PHYSICAL PARAMETERS'
      MLOLIM(12)=0.
      MHILIM(12)=1.
C
      MTYPE(13)=FLOAT
      MPRMPT(13)='Enter the interproton separation S (Angstroms)'
      MTAG(13)='Inter proton separation (Angstroms)'
      MLOLIM(13)=0.
      MHILIM(13)=10.
      MREALS(13)=0.
C
      MTYPE(14)=FLOAT
      MPRMPT(14)=
     + 'Enter value for variational parameter Beta (Angstroms**-1)'
      MTAG(14)='variational parameter Beta (Angstroms**-1)'
      MLOLIM(14)=0.
      MHILIM(14)=10.
      MREALS(14)=.25
C
      MTYPE(15)=SKIP
      MREALS(15)=35.
C
      MTYPE(37)=TITLE
      MPRMPT(37)= 'NUMERICAL PARAMETERS'
      MLOLIM(37)=1.
      MHILIM(37)=1.
C
      MTYPE(38)=TITLE
      MPRMPT(38)='Methods of calculation:'
      MLOLIM(38)=0.
      MHILIM(38)=0.
C
      MTYPE(39)=MTITLE
      MPRMPT(39)='1) Variational'
      MLOLIM(39)=0.
      MHILIM(39)=0.
C
      MTYPE(40)=MTITLE
      MPRMPT(40)='2) Path Integral Monte Carlo'
      MLOLIM(40)=0.
      MHILIM(40)=1.
C
      MTYPE(41)=MCHOIC
      MPRMPT(41)='Make a menu choice and press return'
      MTAG(41)='44 42'
      MLOLIM(41)=1.
      MHILIM(41)=2.
      MINTS(41)=1
      MREALS(41)=1.
C
      MTYPE(42)=NUM
      MPRMPT(42)= 'Enter size of the ensemble'
      MTAG(42)= 'Ensemble size'
      MLOLIM(42)=1.
      MHILIM(42)=MAXENS
      MINTS(42)=20.
C
      MTYPE(43)=FLOAT
      MPRMPT(43)='Enter time step (units of 1E-16 sec/hbar)'
      MTAG(43)='Time step (units of 1E-16 sec/hbar)'
      MLOLIM(43)=0.
      MHILIM(43)=10.
      MREALS(43)=.01
C
      MTYPE(44)=FLOAT
      MPRMPT(44)= 'Enter step size for sampling PHI (Angstroms)'
      MTAG(44)= 'Sampling step size (Angstroms)'
      MLOLIM(44)=.01
      MHILIM(44)=10.
      MREALS(44)=.4
C
      MTYPE(45)=NUM
      MPRMPT(45)= 'Number of thermalization sweeps'
      MTAG(45)= 'Thermalization sweeps'
      MLOLIM(45)=0
      MHILIM(45)=1000
      MINTS(45)=20
C
      MTYPE(46)=TITLE
      MPRMPT(46)='Quantity to calculate:'
      MLOLIM(46)=1.
      MHILIM(46)=0.
C
      MTYPE(47)=MTITLE
      MPRMPT(47)='1) Energy'
      MLOLIM(47)=0.
      MHILIM(47)=0.
C
      MTYPE(48)=MTITLE
      MPRMPT(48)='2) Correlations'
      MLOLIM(48)=0.
      MHILIM(48)=1.
C
      MTYPE(49)=MCHOIC
      MPRMPT(49)='Make a menu choice and press return'
      MTAG(49)='50 53'
      MLOLIM(49)=1.
      MHILIM(49)=2.
      MINTS(49)=1
      MREALS(49)=1.
C
      MTYPE(50)=NUM
      MPRMPT(50)=  'Enter sampling frequency (to avoid correlations)'
      MTAG(50)= 'Sampling frequency'
      MLOLIM(50)=1
      MHILIM(50)=100
      MINTS(50)=6
C
      MTYPE(51)=NUM
      MPRMPT(51)= 'Enter number of samples in a group'
      MTAG(51)= 'Group sample size'
      MLOLIM(51)=1
      MHILIM(51)=1000
      MINTS(51)=10
C
      MTYPE(52)=SKIP
      MREALS(52)=54.
C
      MTYPE(53)=NUM
      MPRMPT(53)= 'Enter maximum correlation length'
      MTAG(53)= 'Maximum correlation length'
      MLOLIM(53)=1
      MHILIM(53)=100.
      MINTS(53)=40
C
      MTYPE(54)=NUM
      MPRMPT(54)= 'Enter number of groups'
      MTAG(54)= 'Number of groups'
      MLOLIM(54)=1
      MHILIM(54)=1000
      MINTS(54)=10
C
      MTYPE(55)=NUM
      MPRMPT(55)= 'Integer random number seed for init fluctuations'
      MTAG(55)= 'Random number seed'
      MLOLIM(55)=1000.
      MHILIM(55)=99999.
      MINTS(55)=34767
C
      MTYPE(56)=SKIP
      MREALS(56)=60.
C
      MSTRNG(MINTS(75))= 'proj8.txt'
C
      MTYPE(76)=BOOLEN
      MPRMPT(76)='Do you want the SHORT version of the output?'
      MTAG(76)='Short version of output'
      MINTS(76)=0
C
      MTYPE(77)=SKIP
      MREALS(77)=80.
C
      MSTRNG(MINTS(86))= 'proj8.grf'
C
      MTYPE(87)=SKIP
      MREALS(87)=90.
C
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE PARAM
C gets parameters from screen
C ends program on request
C closes old files
C maps menu variables to program variables
C opens new files
C calculates all derivative parameters
C performs checks on parameters
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'MENU.ALL'
      INCLUDE 'IO.ALL'
      INCLUDE 'PARAM.P8'
C Local variables:
      REAL AOLD           !temp variable to search for A
      INTEGER CORR,EPS    !what is being calculated?
      INTEGER PIMC,VARY   !which method?
C gives the map between menu indices and parameters
      INTEGER IS,IBETA,IMETHD,INENSM,IDT,IDELTA,ITHERM,ICALC,
     +        INFREQ,INSMPL,INCORR,IGROUP,IDSEED,ITERSE
      PARAMETER (IS     = 13)
      PARAMETER (IBETA  = 14)
      PARAMETER (IMETHD = 41)
      PARAMETER (INENSM = 42)
      PARAMETER (IDT    = 43)
      PARAMETER (IDELTA = 44)
      PARAMETER (ITHERM = 45)
      PARAMETER (ICALC  = 49)
      PARAMETER (INFREQ = 50)
      PARAMETER (INSMPL = 51)
      PARAMETER (INCORR = 53)
      PARAMETER (IGROUP = 54)
      PARAMETER (IDSEED = 55)
      PARAMETER (ITERSE = 76)
C Functions:
      LOGICAL LOGCVT      !converts 1 and 0 to true and false
      INTEGER GETINT      !get integer from screen
      DATA VARY,PIMC /1,2/
      DATA EPS,CORR /1,2/
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     get input from terminal
      CALL CLEAR
      CALL ASK(1,ISTOP)
C
C     stop program if requested
      IF (MREALS(IMAIN) .EQ. STOP) CALL FINISH
C
C     close files if necessary
      IF (TNAME .NE. MSTRNG(MINTS(ITNAME)))
     +     CALL FLCLOS(TNAME,TUNIT)
      IF (GNAME .NE. MSTRNG(MINTS(IGNAME)))
     +     CALL FLCLOS(GNAME,GUNIT)
C
C     set new parameter values
C     physical and numerical
      S=MREALS(IS)
      BETA=MREALS(IBETA)
      METHOD=MINTS(IMETHD)
      NENSEM=MINTS(INENSM)
      DT=MREALS(IDT)
      DELTA=MREALS(IDELTA)
      NTHERM=MINTS(ITHERM)
      CALC=MINTS(ICALC)
      NFREQ=MINTS(INFREQ)
      NSMPL=MINTS(INSMPL)
      NCORR=MINTS(INCORR)
      NGROUP=MINTS(IGROUP)
      DSEED=DBLE(MINTS(IDSEED))
C
C     text output
      TTERM=LOGCVT(MINTS(ITTERM))
      TFILE=LOGCVT(MINTS(ITFILE))
      TNAME=MSTRNG(MINTS(ITNAME))
      TERSE=LOGCVT(MINTS(ITERSE))
C
C     graphics output
      GTERM=LOGCVT(MINTS(IGTERM))
      GHRDCP=LOGCVT(MINTS(IGHRD))
      GFILE=LOGCVT(MINTS(IGFILE))
      GNAME=MSTRNG(MINTS(IGNAME))
C
C     open files
      IF (TFILE) CALL FLOPEN(TNAME,TUNIT)
      IF (GFILE) CALL FLOPEN(GNAME,GUNIT)
      !files may have been renamed
      MSTRNG(MINTS(ITNAME))=TNAME
      MSTRNG(MINTS(IGNAME))=GNAME
C
C     check parameters for correlations, fix NFREQ, NSMPL
      IF (CALC .EQ. CORR) THEN
         NFREQ=1            !fixed for correlations
         NSMPL=1
         IF ((NGROUP .GT. MAXCRR) .OR. ((NGROUP-NCORR) .LE. 20)) THEN
           WRITE (OUNIT,*) ' '
           WRITE (OUNIT,20)
           WRITE (OUNIT,30) NGROUP,NCORR+20,MAXCRR
20         FORMAT (5X,' For reasonable values of the correlations ')
30         FORMAT (5X,' NGROUP (',I4,') must be between NCORR+20 (',
     +             I4,') and MAXCRR (',I4,')')
           WRITE (OUNIT,*) ' '
           NCORR=GETINT(NCORR,1,100,'Reenter NCORR')
           NGROUP=GETINT(NCORR+100,NCORR+20,MAXCRR,'Re-enter NGROUP')
           MINTS(INCORR)=NCORR
           MINTS(IGROUP)=NGROUP
         END IF
      END IF
C
      CALL CLEAR
C
C     calculate derivative parameters
      A=ABOHR
      AOLD=0.
10    IF (ABS(A-AOLD) .GT. 1.E-6) THEN
         AOLD=A
         A=ABOHR/(1+EXP(-S/AOLD))
      GOTO 10
      END IF
      S2=S/2
      HBMDT=HBM*DT
      SQHBDT=SQRT(HBMDT)
      ALPHA=2*ABOHR
      IF (METHOD .EQ. PIMC) DELTA=1.5*A
C
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE PRMOUT(MUNIT,NLINES)
C write out parameter summary to MUNIT
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'IO.ALL'
      INCLUDE 'PARAM.P8'
C Input variables:
      INTEGER MUNIT                   !fortran unit number
      INTEGER NLINES                  !number of lines sent to terminal
C Local variables:
      INTEGER CORR,EPS                !what is being calculated?
      INTEGER PIMC,VARY               !which method?
      DATA EPS,CORR /1,2/
      DATA VARY,PIMC /1,2/
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      IF (MUNIT .EQ. OUNIT) THEN
         CALL CLEAR
      ELSE
         WRITE (MUNIT,*) ' '
         WRITE (MUNIT,*) ' '
      END IF
C
      WRITE (MUNIT,5)
      WRITE (MUNIT,7) S
      WRITE (MUNIT,8) BETA
      WRITE (MUNIT,9) A
      WRITE (MUNIT,*) ' '
      IF (METHOD .EQ. PIMC) THEN
         WRITE (MUNIT,10) NENSEM, DT
      ELSE
         WRITE (MUNIT,11)
      END IF
      IF (CALC .EQ. CORR) WRITE (MUNIT,12) NCORR
      WRITE (MUNIT,13)DELTA
      WRITE (MUNIT,15) NTHERM
      WRITE (MUNIT,20) NFREQ,NSMPL
      WRITE (MUNIT,*) ' '
C
      NLINES=11
C
5     FORMAT (' Output from project 8:',
     +        ' Monte Carlo solution of the H2 molecule')
7     FORMAT (' Proton separation (Angstroms) = ',F7.4)
8     FORMAT (' Variational parameter Beta (Angstroms**-1) = ',F7.4)
9     FORMAT (' Wave function parameter A (Angstroms) = ',F7.4)

10    FORMAT (' Path Integral Monte Carlo with ensemble size = ',I4,
     +        ' and time step = ',1PE12.5)
11    FORMAT (' Variational Monte Carlo method')
12    FORMAT (' correlations will be calculated up to K = ', I4)
13    FORMAT (' Metropolis step in coordinate space (Angstroms)=',F7.4)
15    FORMAT (' number of thermalization sweeps =',I4)
20    FORMAT (' sweep frequency = ',I4,' group size =',I4)
C
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE TXTOUT(IGRP,ENERGY,EVALUE,SIG1,SIG2,U,ACCPT,MUNIT)
C write out results to MUNIT
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'PARAM.P8'
      INCLUDE 'IO.ALL'
C Input variables:
C     energy has two indices
C     first index is the level: sweep, group, or total
C     second index is the value: quantity, quant**2, or sigma**2
      REAL ENERGY(3,3)           !energy
      INTEGER IGRP               !group index
      REAL EVALUE                !current average energy
      REAL SIG1,SIG2             !uncertainties in energy
      REAL U                     !total energy of the system at this S
      REAL ACCPT                 !acceptance ratio
      INTEGER MUNIT              !unit to write to
C Local variables:
      INTEGER SWEEP,GROUP,TOTAL  !which level of calculation
      INTEGER VALUE,SQUARE,SIGSQ !which quantity
      INTEGER PIMC,VARY          !which method?
      DATA SWEEP,GROUP,TOTAL/1,2,3/
      DATA VALUE,SQUARE,SIGSQ/1,2,3/
      DATA VARY,PIMC /1,2/
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      WRITE (MUNIT,10) IGRP,NGROUP,
     +        ENERGY(GROUP,VALUE),SQRT(ENERGY(GROUP,SIGSQ))
      IF (METHOD .EQ. VARY) THEN
        WRITE (MUNIT,20) EVALUE,SIG1,SIG2,U,ACCPT/IGRP/NFREQ/NSMPL
      ELSE
        WRITE (MUNIT,30) EVALUE,SIG1,SIG2,U
      END IF
      IF (MUNIT .EQ. TUNIT) WRITE (MUNIT,*) ' '
C
      IF ((MUNIT .EQ. OUNIT) .AND. (.NOT. TERSE))
     +      CALL PAUSE('to continue...',1)

10    FORMAT (2X,'Group ', I4,' of ', I4,5X,'Eigenvalue = ',F9.4,
     +       ' +- ',F8.4)
20    FORMAT (2X,'Grand average E =',F9.4,'+-',F8.4,'/',F8.4,
     +        '   U=',F9.4,'  acceptance=',F6.4)
30    FORMAT (2X,'Grand average E =',F9.4,'+-',F8.4,'/',F8.4,
     +        '   U=',F9.4)
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE GRFOUT(DEVICE,ECORR)
C outputs energy auto-correlation vs. correlation length
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables
      INCLUDE 'IO.ALL'
      INCLUDE 'PARAM.P8'
      INCLUDE 'GRFDAT.ALL'
C Input variables:
      REAL ECORR(0:MAXCRR)         !energy auto-correlations
      INTEGER DEVICE               !which device is being used?
C Local variables
      REAL K(0:MAXCRR)             !correlation length
      INTEGER IK                   !correlation length index
      CHARACTER*9 CB,CS,CG         !Beta, S, NGROUP as character data
      INTEGER SCREEN               !send to terminal
      INTEGER PAPER                !make a hardcopy
      INTEGER FILE                 !send to a file
      INTEGER LB,LS,LG             !true length of character data
      DATA SCREEN,PAPER,FILE/1,2,3/
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     messages for the impatient
      IF (DEVICE .NE. SCREEN) WRITE (OUNIT,100)
C
C     calculate parameters for graphing
      IF (DEVICE .NE. FILE) THEN
          NPLOT=1                      !how many plots?
          IPLOT=1
C
          YMIN=-1.                     !limits on plot
          YMAX=1.
          XMIN=0.
          XMAX=NCORR
          X0VAL=0.
          Y0VAL=XMIN
C
          NPOINT=NCORR+1
C
          ILINE=1                      !line and symbol styles
          ISYM=1
          IFREQ=1
          NXTICK=5
          NYTICK=5
C
          CALL CONVRT(BETA,CB,LB)             !titles and labels
          CALL CONVRT(S,CS,LS)
          CALL ICNVRT(NGROUP,CG,LG)
          INFO='NGROUP = '//CG(1:LG)
          TITLE = 'H2 molecule, S='//CS(1:LS)//',  Beta='//CB(1:LB)
          LABEL(1)= 'Correlation length'
          LABEL(2)= 'Energy auto-correlation'
C
          CALL GTDEV(DEVICE)                  !device nomination
          IF (DEVICE .EQ. SCREEN) CALL GMODE  !change to graphics mode
          CALL LNLNAX                         !draw axes
      END IF
C
      DO 10 IK=0,NCORR                        !fill array of corr length
         K(IK)=REAL(IK)
10    CONTINUE
C
C     output results
      IF (DEVICE .EQ. FILE) THEN
          WRITE (GUNIT,*) ' '
          WRITE (GUNIT,25) NGROUP
          WRITE (GUNIT,70) (K(IK),ECORR(IK),IK=0,NCORR)
      ELSE
          CALL XYPLOT (K,ECORR)
      END IF
C
C     end graphing session
      IF (DEVICE .NE. FILE) CALL GPAGE(DEVICE)!close graphics package
      IF (DEVICE .EQ. SCREEN) CALL TMODE      !switch to text mode
C
70    FORMAT (2(5X,E11.3))
25    FORMAT (6X,'corr length',5X,
     +    'energy auto-correlation for NGROUP=',I5)
100   FORMAT (/,' Patience, please; output going to a file.')

      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE CRROUT(MUNIT,ECORR)
C write out correlations to MUNIT
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C Global variables:
      INCLUDE 'PARAM .P8'
      INCLUDE 'IO.ALL'
C Input variables:
      REAL ECORR(0:MAXCRR)         !energy auto-correlations
      INTEGER MUNIT                !unit to write to
C Local variables:
      INTEGER K                    !correlation length
      INTEGER NLINES               !number of lines written to screen
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      IF (MUNIT .EQ. OUNIT) THEN
         CALL CLEAR
      ELSE
         WRITE (MUNIT,*) ' '
      END IF
C
      NLINES=1
      WRITE (MUNIT,30) NGROUP
30    FORMAT(' Correlations with NGROUP = ',I5)
      DO 10 K=0,NCORR
         NLINES=NLINES+1
         WRITE (MUNIT,20) K,ECORR(K)
         IF ((MUNIT .EQ. OUNIT) .AND. (MOD(NLINES,TRMLIN-3).EQ. 0)) THEN
            CALL PAUSE('to continue...',0)
            NLINES=0
         END IF
10    CONTINUE
      IF (MUNIT .NE. OUNIT) WRITE (MUNIT,*) ' '
20    FORMAT (5X,' Correlation length = ', I3, 5X,
     +        'Energy auto-correlation = ', F12.5)
      RETURN
      END

ext: for


text:
	SUBROUTINE SFMPRM  ( sfparm, parms, nparm, colors,
     +			     prmlst, chrflg, ncprm, prcons, wndflg,
     +			     icolor, ccvals, icclrs, numccc, icrprm,
     +			     iaddcl, endflg, iret )
C************************************************************************
C* SFMPRM								*
C*									*
C* This subroutine extracts and sets up parameters to process.		*
C*									*
C* SFMPRM ( SFPARM, PARMS, NPARM, COLORS, PRMLST, CHRFLG, NCPRM,        *
C*          PRCONS, WNDFLG, ICOLOR, CCVALS, ICCLRS, NUMCCC, ICRPRM,     *
C*          IADDCL, ENDFLG, IRET )			                *
C*									*
C* Input parameters:							*
C*	SFPARM		CHAR*		Input parameter string		*
C*	PARMS (*)	CHAR*		Data set parameters		*
C*	NPARM		INTEGER		Number of parameters in dataset	*
C*	COLORS		CHAR*		COLORS input			*
C*									*
C* Output parameters:							*
C*	PRMLST (NCPRM)	CHAR*		Output parameters 		*
C*	CHRFLG (NCPRM)	LOGICAL		Character flags			*
C*	NCPRM		INTEGER		Number of valid parameters	*
C*	PRCONS (NCPRM)	CHAR*		Parameters and conditions	*
C*	WNDFLG 		LOGICAL		Flag for existence of wind parm *
C*      ICOLOR (*)      INTEGER         Colors array                    *
C*      CCVALS (*)      REAL            Values for color coding         *
C*      ICCLRS (*)      INTEGER         Colors for color coding         *
C*      NUMCCC (*)      INTEGER         Number of colors for coding     *
C*      ICRPRM (*)      INTEGER         Pointers to reference parms     *
C*	IADDCL		INTEGER		Number of added reference parms *
C*	ENDFLG		CHAR*		Value range end point flag	*
C*					  L = Lower data range		*
C*					  U = Upper data range		*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C**									*
C* Log:									*
C* I. Graffman/RDS	 7/87	GEMPAK4					*
C* M. desJardins/GSFC	 6/88	Reorganized; eliminated check for 	*
C*				same parm string since it didn't work	*
C* M. desJardins/GSFC	11/89	Added conditions on parameters		*
C* S. Schotz/GSC	 4/90   Added weather/cloud symbols , plot at	*
C*				station option, and cleaned up 		*
C* S. Schotz/GSC	 5/90   Removed calls to SFMCOF and SFMCLR	*
C* S. Schotz/GSC	 8/90	Add parameter-condition string output	*
C* K. Brill/NMC		10/91	Changes for the cloud symbols		*
C* M. desJardins/NMC	11/91	Add symbol plotting & change winds	*
C* K. Brill/NMC		11/91	Added flag for winds			*
C* M. desJardins/NMC	11/91	Remove DSET; eliminate "ONE"; clean up	*
C* K. Brill/NMC		03/92	Added DARR				*
C* P. Bruehl/Unidata	 1/94	Added SKYM or SKYK as SKYC+BRBM/K	*
C* S. Maxwell/GSC   	 3/97	Added call to IN_MRKR            	*
C* D. Kidwell/NCEP  	 2/98	Added processing of color coding info,  *
C*				removed NEWFIL from calling sequence    *
C* D. Kidwell/NCEP  	 5/98	Added DAWV                              *
C* A. Hardy/GSC          2/99   Increased MXPRPL 11 -> 25               *
C* S. Jacobs/NCEP	 3/99	Added endflg; Changed call to SNMCLR	*
C* A. Hardy/GSC         10/99   Added ICSY and TBSY                     *
C* G. Grosshans/NCEP	10/99	Added BRBS				*
C* D. Kidwell/NCEP  	 3/02	Added DASH                              *
C* D. Kidwell/NCEP  	 9/02	Added TWSY, TSKC and BRGK               *
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	REAL		ccvals (*)
	CHARACTER*(*)	sfparm, parms (*), colors, prmlst (*),
     + 			prcons(*), endflg
	INTEGER		icolor (*), icclrs (*), numccc (*), icrprm (*)
	LOGICAL		chrflg (*), wndflg
C*
	CHARACTER	prmcnd (MMPARM)*24, p*1, condtn*12,
     +			prmnew (MMPARM)*4
	LOGICAL		cmpflg (MMPARM), levflg (MMPARM), skyflg,
     +			skywnd
C
C*      MXPRPL is the maximum number of parameters to be plotted
C
	INTEGER         MXPRPL
	PARAMETER	( MXPRPL = 25 )
C------------------------------------------------------------------------
	iret = 0
	DO  i = 1, MXPRPL
	    prcons (i) = ' '
	END DO
	wndflg = .false.
	iaddcl = 0
C
C*	Break the input string into a list of parameters.
C
	CALL IN_PRMC  ( MXPRPL, sfparm, prmlst, prmcnd, ncprm, ier )
C
C*	Exit if there are no parameters to compute.
C
	IF  ( ncprm .eq. 0 )  RETURN
C
C*	Process SKYK or SKYM as SKYC + BRBK/BRBM
C*	SKYK or SKYM must be in position 1 in order to plot
C*	sky cover symbol with wind barbs.  Otherwise, numerical
C*	value is plotted.
C
	skywnd = .true.
	IF ( prmlst(1) .eq. 'SKYK' ) THEN
	    ncprm = ncprm + 1
	    prmlst(ncprm) = 'BRBK'
	    prmcnd(ncprm) = prmcnd(1)
	    prmlst(1) = 'SKYC'
	  ELSE IF ( prmlst(1) .eq. 'SKYM' ) THEN
	    ncprm = ncprm + 1
	    prmlst(ncprm) = 'BRBM'
	    prmcnd(ncprm) = prmcnd(1)
	    prmlst(1) = 'SKYC'
	  ELSE
	    skywnd = .false.
	END IF
C
C*	Translate input colors.
C
	CALL SFMCLR ( ncprm, prmlst, colors, icolor, ccvals, icclrs,
     +		      numccc, icrprm, inew, prmnew, endflg, ier )
C
C*	Make sure sky condition and barb plot in same color if SKYK
C*	or SKYM was specified.
C
	IF (skywnd) icolor (ncprm) = icolor (1)
C
C*	Add new reference parameters to the parm list.
C
	iaddcl = inew - ncprm
	DO i = ncprm + 1, inew
	    prmlst (i) = prmnew (i)
	    prmcnd (i) = ' '
	    icolor (i) = 0
	END DO
	ncprm  = inew
C
C*	Determine calculable parameters.
C
	CALL PC_DFLS  ( ncprm, prmlst, chrflg, cmpflg, levflg, np, ier )
C
C*	Check list of calculable parameters.
C
	knt = 0
	DO  i = 1, ncprm
	    IF  ( cmpflg (i) .or. ( prmlst (i) .eq. 'MARK' ) )  THEN
		knt = knt + 1
	      ELSE IF  ( ( prmlst (i) .eq. 'BLNK' ) .or.
     +			 ( prmlst (i) .eq. 'SPAC' ) )  THEN
	      ELSE
		CALL ER_WMSG  ( 'SFMAP', -4, prmlst (i), ier )
		prmlst (i) = 'BLNK'
C
C*              If this is a reference parameter for color coding,
C*		ignore color coding request and set color to 1.
C
		DO j = 1, ncprm - iaddcl
		    IF ( icolor (j) .lt. 0 ) THEN
			IF ( icrprm ( IABS ( icolor (j) ) ) .eq. i )
     +			     icolor (j) = 1
		    END IF
		END DO
	    END IF
	END DO
C
C*	Return if there are no real parameters to compute.
C
	IF  ( knt .eq. 0 )  THEN
	    ncprm = 0
	    RETURN
	END IF
C
C*      Check for weather/cloud symbol parameters and decode their
C*      inputs.
C
	DO  i = 1, ncprm
	    IF ( prmlst (i) .eq. 'WSYM' .or.
     +		 prmlst (i) .eq. 'TWSY' ) THEN
                CALL IN_WSYM  ( prmcnd (i), condtn, ier )
                prmcnd (i) = condtn
C*
              ELSE IF ( prmlst (i) .eq. 'SKYC' .or.
     +		        prmlst (i) .eq. 'TSKC' ) THEN
                CALL IN_SKYC  ( prmcnd (i), condtn, ier )
                prmcnd (i) = condtn
C*
	      ELSE IF ( prmlst (i) .eq. 'CSYL' .or.
     +		        prmlst (i) .eq. 'CSYM' .or.
     +		        prmlst (i) .eq. 'CSYH' .or.
     +		        prmlst (i) .eq. 'CSYT' ) THEN
	        CALL IN_CLDT  ( prmcnd (i), condtn, ier )
	        prmcnd (i) = condtn
C*
	      ELSE IF ( prmlst (i) .eq. 'PWTH' ) THEN
		CALL IN_PWTH  ( prmcnd (i), condtn, ier )
		prmcnd (i) = condtn
C*
	      ELSE IF ( prmlst (i) .eq. 'PTND' ) THEN
		CALL IN_PTND  ( prmcnd (i), condtn, ier )
		prmcnd (i) = condtn
C*
	      ELSE IF ( prmlst (i) .eq. 'MARK' ) THEN
		CALL IN_MRKR  ( prmcnd (i), condtn, ier )
		prmcnd (i) = condtn
C*
              ELSE IF ( prmlst (i) .eq. 'ICSY' ) THEN
                CALL IN_ICNG  ( prmcnd (i), condtn, ier )
                prmcnd (i) = condtn
C*
              ELSE IF ( prmlst (i) .eq. 'TBSY' ) THEN
                CALL IN_TURB  ( prmcnd (i), condtn, ier )
                prmcnd (i) = condtn
C*
	      ELSE IF ( ( prmlst (i) .eq. 'ARRW' ) .or.
     +			( prmlst (i) .eq. 'ARRM' ) .or.
     +			( prmlst (i) .eq. 'ARRK' ) .or.
     +			( prmlst (i) .eq. 'DARR' ) .or.
     +			( prmlst (i) .eq. 'DAWV' ) .or.
     +			( prmlst (i) .eq. 'DASH' ) .or.
     +			( prmlst (i) .eq. 'BARB' ) .or.
     +			( prmlst (i) .eq. 'BRBS' ) .or.
     +			( prmlst (i) .eq. 'BRBM' ) .or.
     +			( prmlst (i) .eq. 'BRBK' ) .or.
     +			( prmlst (i) .eq. 'BRGK' ) )  THEN
		wndflg = .true.
		p = prmlst (i) (1:1)
		IF  ( prmlst (1) .eq. 'SKYC' )  THEN
		    skyflg = .true.
		  ELSE
		    skyflg = .false.
		END IF
		CALL IN_PWND  ( prmcnd (i), p, skyflg, condtn, ier )
		prmcnd (i) = condtn
		IF ( prmlst (i) .eq. 'BARB' ) prmlst (i) = 'BRBM'
            END IF
        END DO
C
C*	Set conditions on parameters.
C
	CALL PC_SLCD  ( ncprm, prmcnd, ier )
C
C*	Build parameter-condition string
C
	DO  i = 1, ncprm
	    IF  ( prmcnd (i) .ne. ' ' ) THEN
		CALL ST_LSTR (prmlst (i), len, ier )
		prcons (i) = prmlst (i) (:len) // prmcnd (i)
              ELSE
		prcons (i) = prmlst (i)
            END IF
	END DO
C*
	RETURN
	END

ext: f

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


    SUBROUTINE STORE_USER(IRECL)
!=======================================================================
!     LATEST CORRECTION BY
!
!     PURPOSE
!           Store user defined data on binary file
!
!     PROGRAMMED BY:  bjorn.melhus@akersolutions.com
!     CREATED......:  10.05.2021
!=======================================================================
    use cfast
    use clist
    use cfuni
    IMPLICIT NONE
    INTEGER islave,ivar,nvar,ipar,irecl
    CHARACTER*13 CFILE
    CHARACTER*14, allocatable :: ctabvar(:)   ! Table of strings to be updated in the slave files ('WaveTp', 'WaveTP', etc...)
    Logical, allocatable :: ltabvar(:)   ! Table of strings constant?
    CHARACTER*14, allocatable :: cmatpar(:,:) ! Contains set of values to be updated in the slave file, one slave file for each set
    CHARACTER*4 CEXT_SLAVE,CNUM
    CHARACTER*136 CNEW_SLAVE
!
    IF (.not.LANAL) return                   ! Things will be stored if an analyses is performed...
!
!.......................................................................................
!   Should've made a count on the number of variables, i.e. 'WaveTp', 'WaveTP', etc...
!   Never mind, just count them the stupid way..:
    nvar=0
    DO ISLAVE=1,SLAVE%NSLAVE          ! lOOP over slaves...
       DO IVAR=1,SLAVE%NVAR(ISLAVE)   ! lOOP over number of variables for given slave
            nvar=nvar+1
       ENDDO
    ENDDO
!.......................................................................................
!
    allocate (ctabvar(nvar))
    allocate (ltabvar(nvar))
    allocate (CMATPAR(npara,nvar))
!
!   Write the user input to some temporary tables
!   and make a dump to binary fles for later access
    nvar=0
    DO ISLAVE=1,SLAVE%NSLAVE    ! lOOP over slaves...
       DO IVAR=1,SLAVE%NVAR(ISLAVE)
           nvar=nvar+1
           ctabvar(nvar)=SLAVE%CVAR(IVAR,ISLAVE)
           ltabvar(nvar)=SLAVE%CONST(IVAR,ISLAVE)
           DO IPAR=1,NPARA
               cmatpar(ipar,nvar)=SLAVE%CPAR(IPAR,IVAR,ISLAVE)
           ENDDO
       ENDDO
    ENDDO
!--------------------------------------------------------------------
!
!    Store user data on binary file
     cfile=cana//'_USER.BIN'
!    If file exists it will be overwritten and the needed record length
!    (each parameter is 14 characters long) will be:
!    Note: the record length is for intel fortran given in 4 byte words...
!          also..: the first record will contain the record length of the 3 binary files..
!                  so at least 4*3=12 bytes.
!          Correction: The ID of the run is stored, which contains 64 characters...
!          New correction: The maestro file is 128 characters.
     irecl=nvar*4
     if (irecl.lt.128) irecl=128
!
     CALL OPEN_BIN(IBIU,'UNKNOWN',Cfile,IRECL)
!
!     Store number of parameters (NPARA) i.e. number of openfast analyses and
!     number of user defined variables (NVAR), e.g. WaveHs, WaveTp etc.
!     on the first record.
!
     WRITE(IBIU,REC=1) NPARA,NVAR
!
!    Store the user ID for the run if analyses in turned on..
      WRITE(IBIU,REC=2) CIDENT
!
      WRITE(IBIU,REC=3) FAST%CMAESTRO(1)     ! Openfast maestro
!
!    Note: Leave a couple of spare records here, just in case. Max length of irecl as defined above.
!
!      WRITE(IBIU,REC=4)  < whatever >
!      WRITE(IBIU,REC=5)  < whatever >
!      WRITE(IBIU,REC=6)  < whatever >  ! insert whatever is needed and make a retreival routine pf<something>.f90

     WRITE(IBIU,REC=7) (ctabvar(ivar),ivar=1,nvar)      ! The variables defined by user e.g. WaveHs, WaveTp etc.:
     WRITE(IBIU,REC=8) (ltabvar(ivar),ivar=1,nvar)      ! if .true. then constant for all openfast analyses
!
!   And then the user defined parameters for each openfast analysis:
    DO IPAR=1,NPARA
        WRITE(IBIU,REC=8+IPAR) (cmatpar(ipar,ivar),ivar=1,nvar)
    ENDDO
!
    DEALLOCATE (CTABVAR)
    DEALLOCATE (ltabvar)
    DEALLOCATE (CMATPAR)
!
    RETURN
    END


text:
C    @(#)sen_bph2.f	20.8 8/30/00
      subroutine sen_bph2 (mode)
c
c     Mode = 0: Partition phase shifters subsystem
c            1: Restore original phase shifter subsystems
c
      include 'ipfinc/parametr.inc'

      include 'ipfinc/addtbx.inc'
      include 'ipfinc/alpha.inc'
      include 'ipfinc/alpha2.inc'
      include 'ipfinc/amtrx.inc'
      include 'ipfinc/area.inc'
      include 'ipfinc/beta.inc'
      include 'ipfinc/blank.inc'
      include 'ipfinc/bus.inc'
      include 'ipfinc/branch.inc'
      include 'ipfinc/dc.inc'
      include 'ipfinc/dcsln.inc'
      include 'ipfinc/ecvar.inc'
      include 'ipfinc/gendrp.inc'
      include 'ipfinc/ikk.inc'
      include 'ipfinc/intbus.inc'
      include 'ipfinc/lfiles.inc'
      include 'ipfinc/ntotzt.inc'
      include 'ipfinc/optim1.inc'
      include 'ipfinc/prt.inc'
      include 'ipfinc/slnopt.inc'
      include 'ipfinc/slnphs.inc'
      include 'ipfinc/tran.inc'
      include 'ipfinc/ordsta.inc'
      include 'ipfinc/phase.inc'

      common /scratch/ numdiff, ixref(MAXBUS), kdiff(2,MAXBUS),
     &                 fdiff(6,MAXBUS), cdiff(2,MAXBUS),
     &                 lphtx(2,MAXPHS), phtx(6,MAXPHS)
      character cdiff * 1

      double precision aim, bim, angle
      complex * 16 y_16(2,2)
      integer kerr, k1, k2, komp_ykm, sect, p, error
      character id * 1
      logical found
      external komp_ykm, swap_ykm

      kerr = 0

      if (mode. eq. 0) then
c
c        Mode = 0: Partition phase shifters subsystem
c
c        Link TXTIE with new bus numbers
c
         do jt = 1, max0 (jphno, ntxtie)
            lphtx(1,jt) = 0
            lphtx(2,jt) = 0
         enddo
         do jt = 1, ntxtie
            do j = 1, 6
               phtx(j,jt) = 0.0
            enddo
         enddo

         icount = iflag(ntot+1)
         ntopt = 0
         do jt = 1,ntxtie
            k1 = txtie(1,jt)
            k2 = txtie(2,jt)
            it = 1
            found = .false.
            do while (it .le. jphno .and. .not. found)
               k = jphid(1, it)
               m = jphid(2, it)
               if (ordvlt .eq. 1) then
                  kt = k
                  mt = m
               else
                  kt = inp2opt(k)
                  mt = inp2opt(m)
               endif
               if (kt .eq. k1 .and. mt .eq. k2) then
                  lphtx(1,it) = jt
                  lphtx(2,jt) = it
                  found = .true.
               else if (kt .eq. k2 .and. mt .eq. k1) then
                  lphtx(1,it) = -jt
                  lphtx(2,jt) = -it
                  found = .true.
               else
                  it = it + 1
               endif
            enddo
            if (k1 .lt. k2) then
               ntopt = ntopt + 1
               txtie(8,jt) = ntopt
            endif
            ltc = txtie(7,jt)
            if (ltc .gt. 0) then
               txtie(4,jt) = tap(ltc)
            else if (ltc .lt. 0) then
               txtie(4,jt) = -tap(-ltc)
            endif
         enddo

         do jt = 1,ntxtie
            k1 = txtie(1,jt)
            k2 = txtie(2,jt)
            if (k1 .lt. k2) then
            else
               found = .false.
               i = 1
               do while (i .le. ntxtie .and. .not. found)
                  if (txtie(1,i) .eq. k2 .and. txtie(2,i) .eq. k1) then
                     found = .true.
                     txtie(8,jt) = -txtie(8,i)
                  else
                     i = i + 1
                  endif
               enddo
            endif
         enddo

         num_nodes = ntot
c
c        First pass - low to high
c
         do jt = 1,ntxtie
            k1 = txtie(1,jt)
            k2 = txtie(2,jt)
            if (k1 .lt. k2) then
               lt = iabs(lphtx(2,jt))
               if (lt .gt. 0) then
                 call getchr(1, id, jphid(3,lt))
                 sect = jphid(4,lt)
                 p = numbrn (opt2inp(k1), opt2inp(k2), id, sect)
               else
                 p = numbrn (opt2inp(k1), opt2inp(k2), '*', 0)
               endif
               if (brtype(p) .eq. 4) p = brnch_nxt(p)
               call pieqiv (p, y_16, error)
               gkm = dreal(y_16(1,2))
               bkm = dimag(y_16(1,2))
               gmk = dreal(y_16(1,1))
               bmk = dimag(y_16(1,1))
               angle = txtie(4,jt)
               it = txtie(8,jt)
               itx = iabs(it)

               ls = km(k1)
               lf = kmlen(k1) + ls - 1
               l = ls
               found = .false.
               do while (l .le. lf .and. .not. found)
                  if (ikmu(l) .eq. k2) then
                     found = .true.
                  else
                     l = l + 1
                  endif
               enddo
               if (.not. found) call erexit()
C
C              Add residual admittance Y_res = 0.01 - j0.05 to
C              prevent singularity.
C
               gkmu(l) = gkmu(l) - gkm - 0.01d0
               bkmu(l) = bkmu(l) - bkm + 0.05d0
               gkku(k1) = gkku(k1) - gmk + 0.01d0
               bkku(k1) = bkku(k1) - bmk - 0.05d0
               pk = pnetu(k1)
               qk = qnetu(k1)
               call nrpqv (k1, pk, dpk, qk, dqk, vk)
               phtx(1,itx) = dpk
               phtx(2,itx) = dqk
c
c              Insert jflag(*) entity for node k1
c
               if (icount .eq. 0) then
                  do i = 1, num_nodes+1
                     iflag(i) = 1
                  enddo
               endif
               do i = k1+1, num_nodes+1
                  iflag(i) = iflag(i) + 1
               enddo
               next = iflag(k1)
               do i = icount, next, -1
                  jflag(1,i+1) = jflag(1,i)
                  jflag(2,i+1) = jflag(2,i)
               enddo
               jflag(1,next) = 11
               jflag(2,next) = jt
               icount = icount + 1
c
c              Add new pseudo-bus
c
               kt = ntot + itx
               num_nodes = kt
               e(kt) = e(k1) * cos(-angle) - f(k1) * sin(-angle)
               f(kt) = e(k1) * sin(-angle) + f(k1) * cos(-angle)
               pnetu(kt) = 0.0
               qnetu(kt) = 0.0
               kvolt(kt) = 0
               gkku(kt) = gmk
               bkku(kt) = bmk
               yptr = yptr + 1
               km(kt) = yptr
               kmlen(kt) = 1
               ikmu(yptr) = k2
               gkmu(yptr) = -gmk
               bkmu(yptr) = -bmk
               pk = pnetu(kt)
               qk = qnetu(kt)
               call nrpqv (kt, pk, dpk, qk, dqk, vk)
               phtx(3,itx) = dpk
               phtx(4,itx) = dqk
c
c              Insert jflag(*) entity for node kt
c
               next = iflag(kt)
               do i = icount, next+1, -1
                  jflag(1,i+1) = jflag(1,i)
                  jflag(2,i+1) = jflag(2,i)
               enddo
               jflag(1,next) = 12
               jflag(2,next) = jt
               icount = icount + 1
               iflag(kt+1) = icount
C
C              Compute initial flow Pkm in branch
C
               mt = ntot + ntopt + itx
               aim = gmk*e(k1) - bmk*f(k1) + gkm*e(k2) - bkm*f(k2)
               bim = gmk*f(k1) + bmk*e(k1) + gkm*f(k2) + bkm*e(k2)
               e(mt) = e(k1)*aim + f(k1)*bim
               f(mt) = -e(k1)*bim + f(k1)*aim
               phtx(1,itx) = phtx(1,itx) - e(mt)
               phtx(2,itx) = phtx(2,itx) - f(mt)
               phtx(3,itx) = phtx(3,itx) + e(mt)
               phtx(4,itx) = phtx(4,itx) + f(mt)
               write (*, 10000) intbus(k1), intbas(k1), intbus(k2),
     &           intbas(k2), e(mt)*bmva, f(mt)*bmva
10000          format (' Phase shifter ', a8, f6.1, 1x, a8, f6.1,
     &           ' S = ', 2f8.1)

            endif
         enddo
c
c        Second pass - high to low
c
         do jt = 1,ntxtie
            k1 = txtie(1,jt)
            k2 = txtie(2,jt)
            if (k1 .gt. k2) then
               lt = iabs(lphtx(2,jt))
               if (lt .gt. 0) then
                 call getchr(1, id, jphid(3,lt))
                 sect = jphid(4,lt)
                 p = numbrn (opt2inp(k1), opt2inp(k2), id, sect)
               else
                 p = numbrn (opt2inp(k1), opt2inp(k2), '*', 0)
               endif
               if (brtype(p) .eq. 4) p = brnch_nxt(p)
               call pieqiv (p, y_16, error)
               gkm = dreal(y_16(1,2))
               bkm = dimag(y_16(1,2))
               gmk = dreal(y_16(1,1))
               bmk = dimag(y_16(1,1))
               angle = txtie(4,jt)
               it = txtie(8,jt)
               itx = iabs(it)

               ls = km(k1)
               lf = kmlen(k1) + ls - 1
               l = ls
               found = .false.
               do while (l .le. lf .and. .not. found)
                  if (ikmu(l) .eq. k2) then
                     found = .true.
                  else
                     l = l + 1
                  endif
               enddo
               if (.not. found) call erexit()
C
C              Replace branch with symmetric admittance, renumber and
C              resort Y-matrix.
C
               ikmu(l) = ntot + itx
               gkmu(l) = gkmu(l) - gkm - gmk
               bkmu(l) = bkmu(l) - bkm - bmk
               call qiksrt(ls, lf, komp_ykm, swap_ykm)
               pk = pnetu(k1)
               qk = qnetu(k1)
               call nrpqv (k1, pk, dpk, qk, dqk, vk)
               phtx(5,itx) = dpk
               phtx(6,itx) = dqk

            endif
         enddo
c
c        Check modelling accuracy
c
         do jt = 1,ntxtie
            k1 = txtie(1,jt)
            k2 = txtie(2,jt)
            it = txtie(8,jt)
            itx = iabs (it)
            mt = ntot + itx
            if (k1 .lt. k2) then
               if (abs(phtx(1,itx)) .gt. 1.0e-2 .or.
     &             abs(phtx(2,itx)) .gt. 1.0e-2) then
                  write (*, 10010) intbus(k1), intbas(k1), phtx(1,itx),
     &               phtx(2,itx)
10010             format (' Phase shift error - terminal 1 ',
     &               a8, f7.1, ' error = ', 2e11.3)
               endif
               if (abs(phtx(3,itx)) .gt. 1.0e-2 .or.
     &             abs(phtx(4,itx)) .gt. 1.0e-2) then
                  write (*, 10020) intbus(k2), intbas(k2), phtx(3,itx),
     &               phtx(4,itx)
10020             format (' Phase shift error - terminal 2 ',
     &               a8, f7.1, ' error = ', 2e11.3)
               endif
               if (abs(phtx(5,itx)) .gt. 1.0e-2 .or.
     &             abs(phtx(6,itx)) .gt. 1.0e-2) then
                  intbus(mt) = intbus(k1)(1:7) // '&'
                  intbas(mt) = intbas(k1)
                  write (*, 10030) intbus(mt), intbas(mt), phtx(5,itx),
     &               phtx(6,itx)
10030             format (' Phase shift error - terminal 0 ',
     &               a8, f7.1, ' error = ', 2e11.3)
               endif
            endif
         enddo

      else
c
c        Mode = 1: Restore original phase shifter subsystems
c
         do jt = 1,ntxtie
            k1 = txtie(1,jt)
            k2 = txtie(2,jt)
            gkm=txtie(5,jt)
            bkm=txtie(6,jt)
            gmk=txtie(9,jt)
            bmk=txtie(10,jt)
            angle = txtie(4,jt)
            it = txtie(8,jt)

            ls = km(k1)
            lf = kmlen(k1) + ls - 1
            l = ls
            found = .false.
            if (k1 .le. k2) then
               kxx = k2
            else
               kxx = ntot + iabs (it)
            endif
            do while (l .le. lf .and. .not. found)
               if (ikmu(l) .eq. kxx) then
                  found = .true.
               else
                  l = l + 1
               endif
            enddo
            if (.not. found) call erexit()
            if (k1 .lt. k2) then
C
C              Restore admittance Y_res = 0.01 - j0.05 to prevent
C              singularity.
C
               gkmu(l) = gkmu(l) + gkm + 0.01d0
               bkmu(l) = bkmu(l) + bkm - 0.05d0
               gkku(k1) = gkku(k1) + gmk - 0.01d0
               bkku(k1) = bkku(k1) + bmk + 0.05d0
               pk = pnetu(kt)
               qk = qnetu(kt)
               call nrpqv (k1, pk, dpk, qk, dqk, vk)

            else
C
C              Restore branch with assymmetric admittance, renumber and
C              resort Y-matrix.
C
               ikmu(l) = k2
               gkmu(l) = gkmu(l) + gkm + gmk
               bkmu(l) = bkmu(l) + bkm + bmk
               call qiksrt(ls, lf, komp_ykm, swap_ykm)
               pk = pnetu(k2)
               qk = qnetu(k2)
               call nrpqv (k2, pk, dpk, qk, dqk, vk)
            endif
         enddo
      endif
      return
      end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


! From Advance Fortran Programming
! https://www.youtube.com/watch?v=W-JCEoV0wng&list=PLNmACol6lYY4wPUNCK03LOwVTBbca7Ix-&index=3

! Program for demonstrating statistical functions
program stats
    use stats_funcs
    implicit none
    ! Simple printing of strings
    character :: name1*10 = " Fortran " ! 10 units long
    character :: name2*10 = "C++", name6*10
    character ::name3*20, name4*20, name5*10 = " String s "

    10 format("Name ", i1, " is:",A10," Name ", i1," is:",A10, "_Check") ! i1 = integer 1 value, A10 = character of 10 units long, 10 = label 10
    20 format("Name 3 is:",A20," Name 4 is:",A20)

    print 10, 1, name1, 2, name2
    ! Concatenation of strings
    name3 = name1//name2 ! String concatenation
    name4 = trim(name1)//trim(name2)
    name6 = adjustl(name5) ! Adjusting the length of the string by removing trailing spaces from front and back. Only trailing spaces from the formatting is present

    print 20,name3, name4
    print 10,5,name5,6,name6
    print *, "length of string is:" , len(" string s")
    print *, name1              ! Printing part of the string
    print *,name1(:6)
    print *,name1(6:)
    print *,name1(4:7)
end program stats


text:
program fortran_hip
  use hipfort
  use hipfort_check

  implicit none

  interface
     ! dim3(320), dim3(256), 0, 0
     subroutine launch(grid,block,shmem,stream,out,a,b,N) bind(c)
       use iso_c_binding
       use hipfort_types
       implicit none
       type(c_ptr),value :: a, b, out
       integer(c_int), value :: N, shmem
       type(dim3) :: grid, block
       type(c_ptr),value :: stream
     end subroutine
  end interface

  integer(c_int), parameter :: N = 1000000

  real(8),allocatable,dimension(:) :: a,b,out
  real(8),pointer,dimension(:) :: da => null(), db => null(),dout => null()

  !type(dim3) :: grid  = dim3(320,1,1)
  !type(dim3) :: block = dim3(256,1,1)

  integer :: i
  type(hipDeviceProp_t),target :: props
  !
  call hipCheck(hipGetDeviceProperties(props,0))
  write(*,"(a)",advance="no") "-- Running test 'vecadd' (Fortran 2008 interfaces)"
  write(*,"(a)",advance="no") "- device: "
  i=1
  do while ( iachar(props%name(i)) .ne. 0 ) ! print till end char
    write(*,"(a)",advance="no") props%name(i)
    i = i+1
  end do
  write(*,"(a)",advance="no") " - "

  ! Allocate host memory
  allocate(a(N),b(N),out(N))

  ! Initialize host arrays
  a(:) = 1.0
  b(:) = 1.0

  ! Allocate array space on the device
  call hipCheck(hipMalloc(da,N))
  call hipCheck(hipMalloc(db,N))
  call hipCheck(hipMalloc(dout,N))

  ! Transfer data from host to device memory
  call hipCheck(hipMemcpy(da, a, N, hipMemcpyHostToDevice))
  call hipCheck(hipMemcpy(db, b, N, hipMemcpyHostToDevice))

  ! launch kernel
  call launch(dim3(320),dim3(256),0,c_null_ptr,c_loc(dout),c_loc(da),c_loc(db),N)
  !call launch(grid,block,0,c_null_ptr,c_loc(dout),c_loc(da),c_loc(db),N)
  call hipCheck(hipDeviceSynchronize())

  ! Transfer data back to host memory
  call hipCheck(hipMemcpy(out, dout, N, hipMemcpyDeviceToHost))

  if ( sum(out) .eq. N*2.0 ) then
     print *, "PASSED!"
  else
     print *, "FAILED!"
  endif

  call hipCheck(hipFree(da))
  call hipCheck(hipFree(db))
  call hipCheck(hipFree(dout))

  ! Deallocate host memory
  deallocate(a,b,out)

end program fortran_hip

max_stars_repo_name: arvindcheru/hipfort

max_issues_repo_name: arvindcheru/hipfort

max_forks_repo_name: arvindcheru/hipfort


text: Octus Energy Inc. is a company based in California that helps clients save money by becoming more energy efficient.  The company has two offices, one located in Davis and the other in Carlsbad, CA.  On July 22 of 2010 the company website reported that Octus was awarded a contract from the City of Davis to develop a more energyefficient outdoor lighting system.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
*DECK ATANH
      FUNCTION ATANH (X)
C***BEGIN PROLOGUE  ATANH
C***PURPOSE  Compute the arc hyperbolic tangent.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C4C
C***TYPE      SINGLE PRECISION (ATANH-S, DATANH-D, CATANH-C)
C***KEYWORDS  ARC HYPERBOLIC TANGENT, ATANH, ELEMENTARY FUNCTIONS,
C             FNLIB, INVERSE HYPERBOLIC TANGENT
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C ATANH(X) computes the arc hyperbolic tangent of X.
C
C Series for ATNH       on the interval  0.          to  2.50000D-01
C                                        with weighted error   6.70E-18
C                                         log weighted error  17.17
C                               significant figures required  16.01
C                                    decimal places required  17.76
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  CSEVL, INITS, R1MACH, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770401  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   900326  Removed duplicate information from DESCRIPTION section.
C           (WRB)
C***END PROLOGUE  ATANH
      DIMENSION ATNHCS(15)
      LOGICAL FIRST
      SAVE ATNHCS, NTERMS, DXREL, SQEPS, FIRST
      DATA ATNHCS( 1) /    .0943951023 93195492E0 /
      DATA ATNHCS( 2) /    .0491984370 55786159E0 /
      DATA ATNHCS( 3) /    .0021025935 22455432E0 /
      DATA ATNHCS( 4) /    .0001073554 44977611E0 /
      DATA ATNHCS( 5) /    .0000059782 67249293E0 /
      DATA ATNHCS( 6) /    .0000003505 06203088E0 /
      DATA ATNHCS( 7) /    .0000000212 63743437E0 /
      DATA ATNHCS( 8) /    .0000000013 21694535E0 /
      DATA ATNHCS( 9) /    .0000000000 83658755E0 /
      DATA ATNHCS(10) /    .0000000000 05370503E0 /
      DATA ATNHCS(11) /    .0000000000 00348665E0 /
      DATA ATNHCS(12) /    .0000000000 00022845E0 /
      DATA ATNHCS(13) /    .0000000000 00001508E0 /
      DATA ATNHCS(14) /    .0000000000 00000100E0 /
      DATA ATNHCS(15) /    .0000000000 00000006E0 /
      DATA FIRST /.TRUE./
C***FIRST EXECUTABLE STATEMENT  ATANH
      IF (FIRST) THEN
         NTERMS = INITS (ATNHCS, 15, 0.1*R1MACH(3))
         DXREL = SQRT (R1MACH(4))
         SQEPS = SQRT (3.0*R1MACH(3))
      ENDIF
      FIRST = .FALSE.
C
      Y = ABS(X)
      IF (Y .GE. 1.0) CALL XERMSG ('SLATEC', 'ATANH', 'ABS(X) GE 1', 2,
     +   2)
C
      IF (1.0-Y .LT. DXREL) CALL XERMSG ('SLATEC', 'ATANH',
     +   'ANSWER LT HALF PRECISION BECAUSE ABS(X) TOO NEAR 1', 1, 1)
C
      ATANH = X
      IF (Y.GT.SQEPS .AND. Y.LE.0.5) ATANH = X*(1.0 + CSEVL (8.*X*X-1.,
     1  ATNHCS, NTERMS))
      IF (Y.GT.0.5) ATANH = 0.5*LOG((1.0+X)/(1.0-X))
C
      RETURN
      END


text:
      subroutine viscf
c
c     ******************************************************************
c     *                                                                *
c     *   computes viscosity coefficients                              *
c     *                                                                *
c     ******************************************************************
c
      use dims
c
c     ******************************************************************
c
      use flo_var
      use mesh_var
c
c     ******************************************************************
c
      use flo_param
      use solv_param
      use mg_param
c
c     ******************************************************************
c
      real, dimension(ie,je)            :: u,v,astr,rev0
c
c     ******************************************************************
c
c     following are some useful constants
c
      pi        = 4.*atan(1.)
c     ckr       = (.062/(2.*pi)**4)
c     ckr       = .01915
      ckr       = .0256
c     cwk       = .225
      cwk       = 0.
c     scf       = re/(sqrt(gamma)*rm)
      scf       = (scal*re/chord)/(sqrt(gamma)*rm)
c
c     compute the molecular viscosity
c
      do j=1,je
      do i=1,ie
         tt       = p(i,j)/w(i,j,1)*t0
         rlv(i,j) = 1.461e-06*tt*sqrt(tt)/((tt+110.3)*rmu0)
      end do
      end do
c
c     for laminar flows we are done.
c     for turbulent flows we are also done on the coarser grids.
c
      if (kvis.le.1.or.mode.ne.0) return
c
c     if we are using the baldwin and lomax model call turbbl and return
c
      aturb     = 1.
      if (ncyc.gt.25) aturb = .5
      if (kturb.eq.1) then
         do j=1,je
         do i=1,ie
            rev0(i,j) = rev(i,j)
         end do
         end do
c        call turbbl
         call turb2
         do j=1,je
         do i=1,ie
            rev(i,j) = aturb*rev(i,j)  +(1.  -aturb)*rev0(i,j)
         end do
         end do
         return
      end if
c
c     else start the rng algebraic model
c
      do j=1,je
      do i=1,ie
         u(i,j)   = w(i,j,2)/w(i,j,1)
         v(i,j)   = w(i,j,3)/w(i,j,1)
      end do
      end do

      do i=itl+1,itu
         u(i,1)   = -u(i,2)
         v(i,1)   = -v(i,2)
      end do
c
c
c
      do j=1,jl
      do i=1,il
         dx13      = xc(i,j,1)   - xc(i+1,j+1,1)
         dy13      = xc(i,j,2)   - xc(i+1,j+1,2)
         dx24      = xc(i+1,j,1) - xc(i,j+1,1)
         dy24      = xc(i+1,j,2) - xc(i,j+1,2)
         du13      = u(i,j) - u(i+1,j+1)
         dv13      = v(i,j) - v(i+1,j+1)
         du24      = u(i+1,j) - u(i,j+1)
         dv24      = v(i+1,j) - v(i,j+1)
         ua        = .25*(u(i,j) + u(i+1,j+1) + u(i+1,j) + u(i,j+1))
         va        = .25*(v(i,j) + v(i+1,j+1) + v(i+1,j) + v(i,j+1))
         dsij      = 1./(dx13*dy24 - dx24*dy13)
         dvdx      =  dsij * (dv13*dy24 - dv24*dy13)
         dudy      = -dsij * (du13*dx24 - du24*dx13)
         dudx      =  dsij * (du13*dy24 - du24*dy13)
         dvdy      = -dsij * (dv13*dx24 - dv24*dx13)
         astr(i,j) = (dudy+dvdx)**2
     .               +2.*(dudx**2  +dvdy**2  -((dudx+dvdy)**2)/3.)
      end do
      end do
c
c
c
      call delt
c
c
c
      do 30 j=2,jl
      do 20 i=2,il

      xbi       = .5*(x(i-1,1,1)  +x(i,1,1))
      ybi       = .5*(x(i-1,1,2)  +x(i,1,2))
      astra     = .25*(astr(i-1,j-1)  +astr(i-1,j)
     .                +astr(i,j-1)    +astr(i,j))
      if (i.ge.itl.and.i.le.itu+1) then
         a3        = 1./(.225*abs(ynot(i)))
         ysci      = sqrt((xc(i,j,1)  -xbi)**2  +(xc(i,j,2)  -ybi)**2)
         ysc       = w(i,2,1)/(ysci*w(i,j,1))
         csc       = 1./(ysc+a3)**2
      else
         csc       = (cwk*ynot(i))**2
      end if
c
c     set some parameters
c
      rnul      = rlv(i,j)/w(i,j,1)
      rnut0     = rev(i,j)/w(i,j,1)
      rnul3     = rnul**3
      a11       = ckr*(csc*csc*scf*scf)/rnul**2
      a2        = 75.
      a1        = a11*(astra)
      rnut0     = rnul+rnut0
c
c     solve for the eddy viscosity
c
      if (dim(rnut0*a1,a2).eq.0.) then
         rev(i,j)  = 0.
         go to 20
      else
         rnut      = sqrt(a1)
      end if

      k      = 0
      fac    = a2 - 1.

   11 den    = 1./(4.*rnut*rnut*rnut + fac)
      rnut1  = rnut - (rnut**4+rnut*fac  -rnut0*rnut0*a1)*den

      if (abs((rnut1  -rnut)).le.1.e-3) then
         rev(i,j) = w(i,j,1)*dim(rnut1,rnul)
         go to 20
      else
         k      = k  +1
         if (k.gt.200) then
            write (6,*) ' iteration not converged ',i,j
            write (6,*) ' rnut = ',rnut,' rnut1 =',rnut1
            rev(i,j)  = w(i,j,1)*dim(rnut1,rnul)
            go to 20
         end if
         rnut   = rnut1
         go to 11
      end if

   20 continue
   30 continue
c
c     adjust the near wake
c
      ii        = ie
      do i=2,itl+1
      ii        = ii  -1
      do j=2,jl
         pex       = -(xc(i,2,1)  -xc(itl+1,2,1))/(20.*dsti(itl+1))
         rev(i,j)  = rev(i,j)  +(rev(itl+1,j)  -rev(i,j))*exp(pex)
         pex       = -(xc(ii,2,1)  -xc(itu,2,1))/(20.*dsti(itu))
         rev(ii,j) = rev(ii,j)  +(rev(itu,j)  -rev(ii,j))*exp(pex)
      end do
      end do

      do i=2,il
         ii        = ib  -i
         rev(i,je) = rev(i,jl)
         rev(i,1)  = rev(ii,2)
      end do

      do i=itl+1,itu
         if (xc(i,2,1).le.xtran) then
            do j=1,jl
               rev(i,j)  = 0
            end do
         end if
         rev(i,1)  = -rev(i,2)
      end do

      do j=1,je
         rev(1,j)  = rev(2,j)
         rev(ie,j) = rev(il,j)
      end do

      return

      end

max_stars_repo_name: AlexT-L/RANS

max_issues_repo_name: AlexT-L/RANS

max_forks_repo_name: AlexT-L/RANS


text:
deconv(f, g) = -8, -9, -3, -1, -6, 7
deconv(h, g) = -3, -6, -1, 8, -6, 3, -1, -9, -9, 3, -2, 5, 2, -2, -7, -1

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData