text:
C     PROGRAM TO CONVERT 32 DETECTOR BT1 FILES INTO SINGLE HISTOGRAM GSAS FILES
C     Original: REM 6/93
C     Modifications: J. Stalick
C     Reorganized: B. Toby
C     Updated 9/07:  JKS
C
      INTEGER MAX_IN,MAX_OUT
      PARAMETER (MAX_IN=2000)
      PARAMETER (MAX_OUT=99999)
C     input data
      DIMENSION FINT(32,MAX_IN)
      DIMENSION SIG(32,MAX_IN)

C     normalization
      REAL PRESET,PRESET1,PNORM
C     rebinned data
      real WSINT(32,MAX_IN),sumW(32,MAX_IN),TTSTART(32)
      INTEGER NPT(32)
C     overlapped differences
      DIMENSION AVINTDIFF(32),FCORR(32)

      DIMENSION FINTB(32,MAX_IN)
      DIMENSION SIGB(32,MAX_IN)

      DIMENSION COEF(20)
C     pseudo-single detector arrays
      DIMENSION FINTV(MAX_OUT)
      DIMENSION SIGV(MAX_OUT)
      REAL CHEB(30)

      CHARACTER*80 GCARD,LINE
      CHARACTER*72 CARD
      CHARACTER*5 MONO
      CHARACTER*51 TITLE
      CHARACTER*10 TEN
      CHARACTER*4 ID
      CHARACTER*1 IREF(20)
      integer lench
      integer iargc
      logical errorflag,splineflag,interpolate,overwrite,sumflag,
     $     pipeflag,filterflag,absflag,rescale
      integer i,j,nargs,nfiles,ln,lastdet,glast,NTMP
      character*50 arg,filename,outfile,parmfile
      character*20 files(50)
      integer statcheck
C
C     FIRST SET UP THE INITIAL PARAMETERS
C
      NTERMR = 5
      NIN = 1
      NOUT = 7
      NINST = 8
      NTMP = 9
      statcheck = 0
      open(unit=9,file='stat.check',status='unknown',access='append')
C
C     get commandline options--------------------------------------------------
      nargs =  iargc()
      splineflag = .true.
      interpolate = .true.
      errorflag = .false.
      overwrite = .false.
      pipeflag = .false.
      filterflag = .false.
      sumflag = .false.
      absflag = .false.
      rescale = .true.
      lastdet = 32
      nfiles = 0
C     loop over arguments
      do I=1,nargs
         call getarg(I,arg)
         ln = lench(arg)
C     is this a flag?
         if (arg(1:2) .eq. '-l' .or. arg(1:2) .eq. '-L') then
            read(arg(3:),*,iostat=j) lastdet
            if (j .ne. 0) then
               write(0,*) 'Unable to parse: ',arg
               errorflag = .true.
            elseif (lastdet .le. 0 .or. lastdet .gt. 32) then
               write(0,*) 'Invalid detector number in: ',arg,lastdet
               errorflag = .true.
            endif
         elseif (arg(1:1) .eq. '-') then
            do j=2,ln
               if (arg(j:j) .eq. 'o' .or. arg(j:j) .eq. 'O') then
                  splineflag = .false.
               elseif (arg(j:j) .eq. 'm' .or. arg(j:j) .eq. 'M') then
                  interpolate = .false.
               elseif (arg(j:j) .eq. 'c' .or. arg(j:j) .eq. 'C') then
                  overwrite = .true.
               elseif (arg(j:j) .eq. 's' .or. arg(j:j) .eq. 'S') then
                  sumflag = .true.
               elseif (arg(j:j) .eq. 'f' .or. arg(j:j) .eq. 'F') then
                  filterflag = .true.
		  pipeflag = .true.
               elseif (arg(j:j) .eq. 'a' .or. arg(j:j) .eq. 'A') then
                  absflag = .true.
               elseif (arg(j:j) .eq. 'r' .or. arg(j:j) .eq. 'R') then
                  rescale = .false.
               else
                  write(0,'(/,5x,2a,/)')
     $                 'Invalid command line flag: -',arg(j:j)
                  nargs = 0
                  goto 100
               endif
            enddo
         else
C     remove the .raw, if present
            if (arg(ln-3:ln) .eq. '.raw') then
               ln = ln - 4
            endif
C     is this a valid file?
            open (unit=1,file=arg(1:ln)//'.raw',status='old',
     $           iostat=j)
            if (j .ne. 0) then
               write(0,'(/,5x,2a)') 'Error opening file ',
     $              arg(1:ln)//'.raw'
               errorflag = .true.
            else
               nfiles = nfiles + 1
               if (nfiles .gt. 50) then
                  write(0,*)
     $                 'Sorry, this program can only process 50 files'
                  stop
               endif
               files(nfiles) = arg(:ln)
            endif
            close(unit=1)
         endif
C     end of loop over arguments
      enddo

      if (nfiles .eq. 1 .and. sumflag) then
         write (0,*) 'Only 1 file to process, ignoring -s flag'
         sumflag = .false.
      endif
      if (nfiles .ne. 0 .and. filterflag) then
         write (0,*) 'You cannot specify -f and a input file'
         filterflag = .false.
         nfiles = 0
      elseif (filterflag) then
         nfiles = 1
      endif
 100  if (nargs .le. 0 .or. errorflag .or. nfiles .le. 0) then
         write(0,'(/,1x,a)')
     $        'Transform BT1 data to pseudo-single'//
     $        ' detector data for GSAS'
         write(0,'(/,1x,a)') 'Usage:',
     $        '  gformat [-Lxx] [-coms] file1[.raw] [file2[.raw]] ...',
     $        '  gformat -f [-Lxx] [-om] < file1.raw > file2.gsas'
         write(0,'(9x,3a)') '-c: ',
     $        ' overwrite (clobber) existing output files'
         write(0,'(9x,3a)') '-m:  move points to closest bin',
     $        '[default: interpolate]'
         write(0,'(9x,2a)') '-s: ',
     $        'sum multiple files '//
     $        '[default: process each file independently'
         write(0,'(9x,3a)') '-o: ',
     $        ' omit background spline',
     $        ' [default: add offsets to spline backgd]'
         write(0,'(9x,2a)')
     $        '-f: run as filter: read from stdin and write to stdout'
         write(0,'(9x,2a)') '-Lxx or -lxx: ',
     $        'drop detectors xx-32 [default: use all]'
         write(0,'(9x,2a,/)') '-a: ',
     $        'apply absorption correction from file abs.corr'
         write(0,'(9x,2a,/)') '-r: ',
     $        'prevent rescaling to ~counts'
         stop
      endif
C     done getting commandline options-----------------------------------------
      WRITE (0,10)
 10   FORMAT ('Transformation of 32 det BT1 data to GSAS format',/)

C     status messages
      if (splineflag) then
         write (0,*)
     $        'Offsets will be added to data to spline backgrounds'
      else
         write (0,*)
     $        'No spline correction will be applied to backgrounds'
      endif
      if (interpolate) then
         write (0,*) 'Data will be interpolated'
      else
         write (0,*) 'Data will not be interpolated'
      endif
      if (overwrite) then
         write (0,*) 'Will overwrite existing output files!'
      endif

C-----------------------------------------------------------------------
C     loop over files       ii is the file number
C-----------------------------------------------------------------------
      do ii=1,nfiles
C     open the input file
	 if (.not. filterflag) then
            filename = files(ii)(:lench(files(ii))) // '.raw'
            OPEN (UNIT=NIN,  FILE= filename, STATUS= 'OLD',
     $           iostat=j)
            if (j .ne. 0) then
               write(0,'(/,5x,2a)') 'Error opening file ',
     $              files(ii)(:lench(files(ii))) // '.raw'
               stop
            endif
            WRITE (0,'(a)') ' Reading input from '//filename
C     put the filename into the statcheck file, always
            write (9,*) 'file ',files(ii)(:lench(files(ii))) // '.raw'
         else
            filename = 'stdin'
            WRITE (0,'(a)') ' Reading input from stdin'
            nin = 5
         endif
         if (ii .eq. 1 .or. .not. sumflag) then
C     zero the arrays
            DO N=1,MAX_IN
               DO J=1,32
                  WSINT(J,N) = 0
                  sumW(J,N) = 0
               ENDDO
            ENDDO
            DO J=1,32
               npt(j) = 0
            ENDDO
         endif
C     open output file
         if ((ii .eq. 1 .or. .not. sumflag) .and. (.not. pipeflag)) then
            errorflag = .false.
            outfile = files(ii)(:lench(files(ii))) // '.gsas'
            parmfile = files(ii)(:lench(files(ii))) // '.inst'
            OPEN (UNIT=NOUT, FILE= outfile, STATUS= 'NEW', iostat=j)
            if (overwrite .and. j .ne. 0) THEN
               OPEN (UNIT=NOUT, FILE= outfile, STATUS= 'OLD', iostat=j)
               if (j .ne. 0) then
                  WRITE(0,'(3a)') '  ***ERROR***  Unable to overwrite ',
     $                 outfile(:lench(outfile)),
     $                 ' Protection problem?'
                  errorflag = .true.
               ENDIF
            elseif (j .ne. 0) then
               WRITE(0,'(3a)') '  ***ERROR***  Unable to create ',
     $              outfile(:lench(outfile)),
     $              ' Does it already exist?'
               errorflag = .true.
            endif
            if (.not. errorflag) write(0,*)
     $           'Creating ',outfile(:lench(outfile))
            if (errorflag) stop
         elseif (ii .eq. 1 .or. .not. sumflag) then
            WRITE (0,'(a)') ' Sending output to stdout'
            nout = 6
         endif
C
C
C     READ IN FIRST TWO LINES OF INPUT FILE
C
C
         READ (nin,'(a50,f8.0,2i3,1x,a5,g14.1)')
     $        TITLE,WL1,icol1,icol2,mono,preset
         if (ii .eq. 1 .or. .not. sumflag) then
            PRESET1 = PRESET
         ENDIF
         IF (sumflag .and. (PRESET .EQ. 0 .or. PRESET1 .eq. 0)) THEN
            PNORM = 1.0
            WRITE(0,*) 'No prefactor in ',
     $           files(ii)(:lench(files(ii))) // '.raw',
     $           '. Rerun proprep to sum.'
         ELSEIF (sumflag) THEN
            PNORM = PRESET1 /  PRESET
         ELSE
            PNORM = 1.0
         ENDIF
         write (0,*) 'PNORM=',pnorm

         READ (NIN,'(A72)') CARD
         ID = CARD(1:4)
!     write(0,*) 'colls = ',icol1,icol2
!     write(0,*) 'mono  = ',mono
C
C-----OUTPUT HEADER #1
C
C     create the title line
         if (ii .eq. 1 .or. .not. sumflag) then
            GCARD = ' '
            GCARD(1:51) = TITLE
            glast = 80
         endif
C     add the name of each file into the header
         ln = lench(files(ii))
         if (glast .gt. 20) gcard(glast-ln:glast) = files(ii)(:ln)
!     write (*,*) glast
!     write (*,*) gcard,ln,files(ii)(:ln)

         glast = glast - ln - 1
C
C
C     READ IN THE INTENSITIES, AND CALCULATE THE CORRECT 2THETA
C     VALUES ASSOCIATED WITH THEM
C
C
         NUMDET = 0
C-----------------------------------------------------------------------
C     loop over detectors      J is the detector number
C-----------------------------------------------------------------------
         do j=1,lastdet
C
C     READ BANK HEADER
            READ (NIN,'(F8.0,F8.0,F8.0,2F8.4,2I5)',err=412)
     $           TMIN, FSTEP, TMAX, FMULT, ZERO,IDET,IPT
            TMIN = TMIN/100.
            TMAX = TMAX/100.
            FSTEP = FSTEP/100.
            ZERO = ZERO/100.
            NUMDET = NUMDET + 1
            if (ipt .gt. MAX_IN) then
               write(0,*) 'Number of input points greater than MAX_IN',
     $              ipt,max_in
               stop
            endif

            READ (NIN,'(12F6.0)') (FINT(J,I),I=1,ipt)

C
C-----CORRECT INTENSITIES and esd's USING DETECTOR SCALES (FMULT)
C
            DO I = 1,IPT
               if (fint(j,i) .gt. 0)
     $              sig(J,i) = fmult * sqrt(fint(J,i)) * PNORM
               FINT(J,I) = FINT(J,I) * FMULT * PNORM
               IF (fint(j,i) .eq. 0) sig(J,i) = 1.
            enddo

C     1st file, 1st detector
            IF (ii .eq. 1 .and. j .eq. 1 .or. .not. sumflag) then
               IF (INTERPOLATE) THEN
C     GSTEP is output step size (output=input for interpolation, doubled otherwise)
C     this could be a user controlled option someday
                  GSTEP = FSTEP
               ELSE
!     GSTEP = 2.*FSTEP
                  GSTEP = FSTEP
               ENDIF
            ENDIF

C     1st file, 1st point (all detectors) -- save starting angle
            IF (ii .eq. 1 .or. .not. sumflag) then
C     correct 2theta
               FTHETA = TMIN - ZERO
               IF (INTERPOLATE) THEN
                  REMAIN = FTHETA - INT(FTHETA/GSTEP)*GSTEP
                  N = 1
C     ignore a difference of less than 0.01 centidegree (plus 0.00001 for roundoff)
                  IF (ABS(REMAIN) .LT. 0.00011) N = 0
C     compute the 2theta that is an even multiple of the output stepsize and
C     larger than the starting angle
                  TTSTART(J)  = (INT(FTHETA/GSTEP)+N)*GSTEP
               ELSE
C     compute a 2theta that is an even multiple of the output stepsize)
                  TTSTART(J) = (NINT(FTHETA/GSTEP))*GSTEP
               ENDIF
!     write (*,*) j,' 2theta, bin start', FTHETA, TTSTART(J)
            ENDIF
C
C-------CORRECT INTENSITIES FOR BIN--Assuming straight lines
C     successive points
C
C
            IF (INTERPOLATE) IPT = IPT - 1
            DO I = 1,IPT
C     correct 2theta
               FTHETA = TMIN - ZERO + (I-1)*FSTEP
               IF (INTERPOLATE) THEN
C
C------CALCULATE THE RIGHT THETA VALUE FOR GSAS BIN
C     THIS ROUNDS UP TO THE NEAREST MULTIPLE OF GSTEP
                  REMAIN = FTHETA - INT(FTHETA/GSTEP)*GSTEP
                  N = 1
C     ignore a difference of less than 0.01 centidegree (plus 0.00001 for roundoff)
                  IF (ABS(REMAIN) .LT. 0.00011) N = 0
C     compute a 2theta that is an even multiple of the output stepsize)
                  FTHETAC = (INT(FTHETA/GSTEP)+N)*GSTEP
                  RATIO = (FTHETAC - FTHETA)/FSTEP
                  FINTERP = FINT(J,I) +
     $                 (FINT(J,I+1)-FINT(J,I)) * RATIO
                  sigintepb = sig(J,i) +
     $                 (sig(J,i+1)-sig(J,i)) * RATIO
				  c1 = (1.-1./RATIO) * (1.-1./RATIO)
				  c2 = 1./RATIO/RATIO
				  sigj2 = sig(J,i) * sig(J,i)
				  sigj2b = sig(J,i+1) * sig(J,i+1)
                  sigintep = sqrt( c1 * sigj2 + c2 * sigj2b)

                  WRITE (0,*) 'sig is times',sigintep
                  WRITE (0,*) 'sigb is times',sigintepb

C     deal with negative intensities (ignore them)
                  IF (FINT(J,I) .lt. 0 .or. FINT(J,I+1) .lt. 0)
     $                 sigintep = 0.0
               ELSE
C     no interpolation
C------CALCULATE THE CLOSEST GSAS BIN
                  FTHETAC = (NINT(FTHETA/GSTEP))*GSTEP
                  FINTERP = FINT(J,I)
                  sigintep = sig(J,i)
C     deal with negative intensities (ignore them)
                  IF (FINT(J,I) .lt. 0) sigintep = 0.0
               ENDIF
C     compute the point number for this 2theta
               N = NINT((FTHETAC - TTSTART(J))/GSTEP) + 1
               NPT(J) = MAX(NPT(J),N)
               if (sigintep .gt. 0) then
                  s2 = 1./sigintep**2
C     Stat check#1 -- one detector to previous measurements of the same detector
                  IF (sumW(J,N) .gt. 0) then
                     diff =  FINTERP-WSINT(J,N)/sumW(J,N)
                     sigdif = sqrt (1./sumW(J,N) + 1./s2)
                     IF (abs(diff) .gt. 4.*sigdif) then
                        statcheck = statcheck+1
                        write (9,*)
     $     '  Warning: file not consistent with previous'
     $     ,' by ',nint(10*abs(diff)/sigdif)/10.,' sigma'
                        write (9,*) 'detector ',J,' angle',FTHETA
                        write (9,*) 'Intensities',
     $                       WSINT(J,N)/sumW(J,N),FINTERP
                        write (9,*) 'ESD''s      ',
     $                       sqrt(1./sumW(J,N)),sigintep
                     endif
                  ENDIF
C     weighted sum of intensities
                  WSINT(J,N) = WSINT(J,N) + FINTERP*s2
C     sum of weights
                  sumW(J,N) = sumW(J,N) + s2
               ENDIF
            ENDDO
C     End loop over all detectors
         ENDDO
 412     CLOSE(NIN)

C     process and write data unless there are more files to be summed
         IF (II .eq. nfiles .or. .not. SUMFLAG) then
C
C-----Correct for different backgrounds in different
C-----detectors. Use detector #20 as a reference baseline.
C
C-----Find the points in the previous detector that overlap the current
C-----Start with the second detector
            DO J = 1,lastdet-1
               OVERSUM1 = 0.0
               OVERWT1 = 0.0
               OVERSUM2 = 0.0
               OVERWT2 = 0.0
C     compute the point number for the first overlapped point
               N = NINT((TTSTART(J+1) - TTSTART(J))/GSTEP) + 1
               I1 = 0
C     write(*,*) NPT(J)-N+1,' Points',(NPT(J)-N+1)*GSTEP,' deg.'
               DO I=N,NPT(J)
                  I1 = I1 + 1
                  OVERSUM1 = OVERSUM1 + WSINT(J,I)
                  OVERSUM2 = OVERSUM2 + WSINT(J+1,I1)
                  OVERWT1 = OVERWT1 + sumW(J,I)
                  OVERWT2 = OVERWT2 + sumW(J+1,I1)
               ENDDO
               OVERSUM1 =OVERSUM1/OVERWT1
               OVERWT1 = 1./SQRT(OVERWT1)
               OVERSUM2 =OVERSUM2/OVERWT2
               OVERWT2 = 1./SQRT(OVERWT2)
!     write(*,'(2(a,i3,f8.2,a,f7.2))') 'overlap ',J,
!     $              OVERSUM1,' +-',OVERWT1,' compared to ',j+1,
!     $              OVERSUM2,' +-',OVERWT2
C     Is this difference statistically significant?
C     It must be greater than lets say 2 sigma
               IF (ABS(OVERSUM1-OVERSUM2) .gt.
     $              2*SQRT(OVERWT1**2 + OVERWT2**2)) THEN
                  AVINTDIFF(J) = OVERSUM1 - OVERSUM2
!     write(*,*) 'correction=',
!     $   OVERSUM1-OVERSUM2,' +-',SQRT(OVERWT1**2 + OVERWT2**2),
!     $                 ' is meaningful'
               ELSE
                  AVINTDIFF(J) = 0.0
!     write(*,*) 'correction=',
!     $   OVERSUM1-OVERSUM2,' +-',SQRT(OVERWT1**2 + OVERWT2**2),
!     $                 ' not meaningful, not used'
               ENDIF
            ENDDO
C
C---- Calculate the actual corrections that need
C---- to be applied to each detector.
C---- Detector 20 = baseline
C
C     note that Fcorr will be added to the appropriate detector
            do j=1,32
               fcorr(j) = 0.
            enddo

            IF (SPLINEFLAG) THEN
c
c---- Then first 19 detectors
c
               do  i = 19,1,-1
                  fcorr(i) = fcorr(i+1) - AVintdiff(i)
               enddo
C
C---- then remaining detectors
C
               do i = 21,lastdet
                  fcorr(i) = fcorr(i-1) + AVintdiff(i-1)
               enddo
C     tell the user what we will do
               write(0,'(2a)') 'The following background spline ',
     $              'corrections will be added to each detector'
               write(0,'(6(i4,a,f8.2))') (j,':',fcorr(j),j=1,lastdet)
            ENDIF

C
C
C---- Normalize the intensities and apply these corrections to the detectors
C
C
            do  i = 1,lastdet
               do  j = 1,NPT(I)
                  if (sumW(I,J) .gt. 0) THEN
                     fintb(I,j) = wsint(I,j)/sumW(I,J) + fcorr(i)
                     sigb(I,J) = 1./SQRT(sumW(I,J))
                  ELSE
                     fintb(I,j) = 0.0
                     sigb(I,J) = 0.0
                  ENDIF
C     write (13,445) TTSTART(I)+(J-1.)*GSTEP,fintb(i,j)
C     445       format(f8.3,f8.1)
               enddo
            enddo

C     compute range of data and zero array
            TSTART = TTSTART(1)
            TEND = TTSTART(LASTDET) + (NPT(LASTDET)-1.)*GSTEP
C     write (*,*) TTSTART(1),TTSTART(LASTDET),NPT(LASTDET)
            NPNTS = 1 + NINT( (TEND - TSTART)/GSTEP )
            IF (NPNTS .gt. MAX_OUT) THEN
               write(0,*) 'Warning: # points > MAX_OUT',
     $              npnts,MAX_OUT
               write(0,*) 'Not all data will be written'
               NPNTS = MAX_OUT
            ENDIF
            DO I=1,NPNTS
               FINTV(I) = 0.
               SIGV(I) = 0.
            ENDDO
C     Now put into pseudo-single detector array
            DO J=1,LASTDET
               DO I=1,NPT(J)
C     compute point number in master array
                  N = I + NINT((TTSTART(J) - TTSTART(1))/GSTEP)
                  IF (sigb(j,i) .gt. 0 .and. N .le. MAX_OUT) THEN
                     s2 = 1./sigb(j,i)**2
C     Stat check#2 -- current detector to previous one
                     IF (SIGV(N) .gt. 0) then
                        diff =   FINTB(j,i) - FINTV(N)/SIGV(N)
                        sigdif = sqrt (1./SIGV(N) + 1./s2)
                        IF (abs(diff) .gt. 4.*sigdif) then
                           statcheck = statcheck+1
                           write (9,*)
     $    '  Warning: detector not consistent with previous'
     $    ,' by ',nint(10*abs(diff)/sigdif)/10.,' sigma'
                           write (9,*) 'detector ',J,' angle',
     $                          GSTEP*(N-1) + TTSTART(1)
                           write (9,*) 'Intensities',
     $                          FINTV(N)/SIGV(N),FINTB(j,i)
                           write (9,*) 'ESD''s      ',
     $                          sqrt(1./SIGV(N)),sigb(j,i)
                        endif
                     ENDIF
                     FINTV(N) = FINTV(N) + FINTB(j,i) * s2
                     SIGV(N) = SIGV(N) + s2
                  ENDIF
               ENDDO
            ENDDO
C     normalize
            ICOUNT = 0
            DO N=1,NPNTS
               IF ( SIGV(N) .GT. 0) THEN
                  FINTV(N) = FINTV(N) / SIGV(N)
                  SIGV(N) = 1./SQRT(SIGV(N))
               ELSE
                  ICOUNT = ICOUNT+1
               ENDIF
            ENDDO
            IF (ICOUNT .gt. 0) WRITE(0,*) 'Warning ',ICOUNT,
     $           ' Points were not defined'
            IF (ABSflag) THEN
C     absorption correction
               WRITE (0,*) 'Applying Absorption Correction'
               OPEN(unit=99,file='abs.corr',status='old')
C     file abs.corr contains a one line header
C     and then the number of Chebyshev terms, followed by the terms
               read(99,'(A)') title
               read(99,*) NTERM,(CHEB(I),I=1,NTERM)
               close(unit=99)
               WRITE (0,*) '  from File abs.corr', title
               WRITE (0,*) NTERM,'coef: ',(CHEB(I),I=1,NTERM)
               absmin = 1e9
               absmax = 0
               DO N=1,NPNTS
C     compute 2theta in degrees
                  TTHETA = TSTART + (N-1.)*GSTEP
C     change range to +-1
                  T = -1. + TTHETA/90.
                  DY1 = 1.
                  DY2 = T
                  abscor = CHEB(1) + T*CHEB(2)
                  DO I=3,NTERM
                     DYI = 2. * T * DY2-DY1
                     DY1 = DY2
                     DY2 = DYI
                     abscor = abscor + CHEB(I)*DYI
                  ENDDO
                  absmin = min(absmin, abscor)
                  absmax = max(absmax, abscor)
                  FINTV(N) = FINTV(N) * ABSCOR
                  SIGV(N) = SIGV(N) * ABSCOR
               ENDDO
               WRITE (0,*) 'Correction ranged from ',absmin,
     $              ' to ', absmax
            ENDIF

            IF (RESCALE) THEN
               ynorm = 0
               icount = 0
               DO N=1,NPNTS
                  IF ( SIGV(N) .GT. 0) THEN
                     ynorm = ynorm + FINTV(N) / (SIGV(N)*SIGV(N))
                     ICOUNT = ICOUNT+1
                  ENDIF
               ENDDO
               IF (ICOUNT .GT. 0) THEN
                  YNORM = YNORM / ICOUNT
               ELSE
                  YNORM = 1
               ENDIF
               WRITE (0,*) 'Rescaling factor is times',ynorm
               DO N=1,NPNTS
                  IF ( SIGV(N) .GT. 0) THEN
                     FINTV(N) = ynorm * FINTV(N)
                     SIGV(N)  = ynorm * SIGV(N)
                  ENDIF
               ENDDO
            ENDIF
C
C---- OUTPUT TO GSAS FORMAT FILE (NOUT)
C---- NPNTS = NO OF GSAS POINTS
C---- CALCULATE NO OF GSAS RECORDS
C
            IREC = INT((NPNTS-1)/5.)+1

            sstart = INT(100*TSTART)
            sstep = INT(100*GSTEP)
            WRITE(NOUT,'(3A)') GCARD
            if (.not. pipeflag) then
               WRITE(LINE,'(A26,A)') 'Instrument parameter file:',
     $              parmfile(:lench(parmfile))
               WRITE(NOUT,'(3A)') LINE
            endif
            WRITE(LINE,500) NPNTS,IREC,sstart,sstep
 500        FORMAT('BANK  1 ',2i8,' CONST ',2f10.2,' 0 0 ESD',21x)
            WRITE(NOUT,'(3A)') LINE
            N = 1
            DO 530 I = 1,IREC
               WRITE (LINE,520) (FINTV(J),sigv(j),J=N,N+4)
 520           FORMAT (5(F8.0,f8.0))
               WRITE(NOUT,'(3A)') LINE
               N = N + 5
 530        CONTINUE
C
C
            if (.not. pipeflag) then
               OPEN (UNIT=NINST, FILE= parmfile, STATUS= 'unknown',
     $              err=919)
               write(0,*) 'Creating ',parmfile(:lench(parmfile))
C-----OUTPUT INSTRUMENT PARAMETER FILE (NINST)
C
C     OUTPUT FIRST 3 LINES TO INSTRUMENT PARAMETER FILE
C
               TEN = '1234567890'
               WRITE (LINE,535) TEN,TEN,TEN,TEN,TEN,TEN
 535           FORMAT (12X,6A10)
               WRITE(NINST,'(3A)') LINE
               WRITE (LINE,540) 1
 540           FORMAT ('INS   BANK  ',I5)
               WRITE(NINST,'(3A)') LINE
               WRITE (LINE,550)
 550           FORMAT ('INS   HTYPE   PNCR')
               WRITE(NINST,'(3A)') LINE
C
C
C
C-----MISCELANEOUS SET UP PARAMETERS
               DO  I = 1,20
                  IREF(I) = 'N'
               ENDDO
               DO I = 5,20
                  COEF(I) = 0.0
               ENDDO
               IF (mono .eq. 'GE311') THEN
C     Ge311 parameters
                  if (icol1 .eq. 7) then
                     COEF(1) = 246.7
                     COEF(2) = -298.6
                     COEF(3) = 157.6
                     COEF(4) = 4.0305
                  else
                     COEF(1) = 398.5
                     COEF(2) = -343.2
                     COEF(3) = 163.0
                     COEF(4) = 4.0305
                  endif
               ELSEIF  (mono .eq. 'GE733') THEN
C     Ge733 parameters
                  if (icol1 .eq. 7) then
                     COEF(1) = 36.0
                     COEF(2) = -100.7
                     COEF(3) = 106.0
                     COEF(4) = 8.0
                  else
                     COEF(1) = 51.3
                     COEF(2) = -109.0
                     COEF(3) =  105.6
                     COEF(4) = 8.0
                  endif
               ELSE
C     Cu311 parameters
                  if (icol1 .eq. 7) then
                     COEF(1) = 181.2
                     COEF(2) = -299.3
                     COEF(3) = 186.3
                     COEF(4) = 5.0
                  else
                     COEF(1) = 239.7
                     COEF(2) = -298.2
                     COEF(3) = 180.8
                     COEF(4) = 5.0
                  endif
               ENDIF
C
               WRITE (LINE,555) WL1,0.0,ZERO,0
 555           FORMAT ('INS  1 ICONS',3F10.5,I10,10X)
               WRITE(NINST,'(3A)') LINE

               WRITE (LINE,560) ID
 560           FORMAT ('INS  1I HEAD',2X,A4,2X,60X)
               WRITE(NINST,'(3A)') LINE

               WRITE (LINE,565) 0,TSTART,TMAX,1
 565           FORMAT ('INS  1I ITYP',I5,2F10.4,I10)
               WRITE(NINST,'(3A)') LINE

               IPTYP=1
               IDAMP = 0
               WRITE (LINE,575) IPTYP,IPTYP,6,0.005,IDAMP,
     $              (IREF(I),I=1,6)
 575           FORMAT ('INS  1PRCF',I1,1X,2I5,F10.5,4X,I1,20A1)
               WRITE(NINST,'(3A)') LINE

               WRITE (LINE,580) IPTYP,1,(COEF(I),I=1,4)
 580           FORMAT ('INS  1PRCF',2I1,4E15.6)
               WRITE(NINST,'(3A)') LINE

               WRITE (LINE,580) IPTYP,2,(COEF(I),I=5,6)
               WRITE(NINST,'(3A)') LINE
               IPTYP=2
               COEF(7) = COEF(4)
               COEF(4) = 0.0
               WRITE (LINE,575) IPTYP,IPTYP,12,0.005,IDAMP,
     $              (IREF(I),I=1,12)
               WRITE(NINST,'(3A)') LINE

               WRITE (LINE,580) IPTYP,1,(COEF(I),I=1,4)
               WRITE(NINST,'(3A)') LINE

               WRITE (LINE,580) IPTYP,2,(COEF(I),I=5,8)
               WRITE(NINST,'(3A)') LINE

               WRITE (LINE,580) IPTYP,3,(COEF(I),I=9,12)
               WRITE(NINST,'(3A)') LINE
               IPTYP=3
               COEF(7) = 0.04
               COEF(8) = 0.03
               WRITE (LINE,575) IPTYP,IPTYP,19,0.005,IDAMP,
     $              (IREF(I),I=1,19)
               WRITE(NINST,'(3A)') LINE

               WRITE (LINE,580) IPTYP,1,(COEF(I),I=1,4)
               WRITE(NINST,'(3A)') LINE
               WRITE (LINE,580) IPTYP,2,(COEF(I),I=5,8)
               WRITE(NINST,'(3A)') LINE
               WRITE (LINE,580) IPTYP,3,(COEF(I),I=9,12)
               WRITE(NINST,'(3A)') LINE
               WRITE (LINE,580) IPTYP,4,(COEF(I),I=13,16)
               WRITE(NINST,'(3A)') LINE
               WRITE (LINE,580) IPTYP,5,(COEF(I),I=17,19)
               WRITE(NINST,'(3A)') LINE
               CLOSE(ninst)
            endif
            if (.not. pipeflag) then
               CLOSE(nout)
            ENDIF

C     end of if block to write data
 919        CONTINUE
         ENDIF
C     End loop over all files
      ENDDO
      if (statcheck .gt. 0) then
         write(0,*) 'Note: ',statcheck,
     $        ' inconsistency warnings were written to file stat.check'
      endif
      close(9)
      STOP 'OK'
      END
C---------------------------------------------------------------------------
C---------------------------------------------------------------------------
c     Function LENCH
c
c     This function takes a character string and finds out how long the
c     "actual" string is (i.e. not including padded blanks on the right).
c
C---------------------------------------------------------------------------
      integer function lench(string)
      character*(*) string
      character*1 nul
      data nul/"\0"/
      lench=0
      if (string.eq.' '.or.string(1:1).eq.nul) return
      do lench=len(string),1,-1
         if (string(lench:lench).ne.' '.and.string(lench:lench).ne.nul)
     $        return
      enddo
      lench=0
      return
      end
C---------------------------------------------------------------------------
C     Convert strings to lower case characters
      SUBROUTINE DOWNCASE(STRING)
      CHARACTER*(*)   STRING    !String to convert
      NCH = LEN(STRING)
      DO ICH=1,NCH
         JCH = ICHAR(STRING(ICH:ICH))
         IF ( JCH.GE.ICHAR('A') .AND. JCH.LE.ICHAR('Z') )
     $        STRING(ICH:ICH) = CHAR(JCH-ICHAR('A')+ICHAR('a'))
      END DO
      RETURN
      END
C---------------------------------------------------------------------------
C     Convert strings to upper case characters
      SUBROUTINE UPCASE(STRING)
      CHARACTER*(*)   STRING    !String to convert
      NCH = LEN(STRING)
      DO ICH=1,NCH
         JCH = ICHAR(STRING(ICH:ICH))
         IF ( JCH.GE.ICHAR('a') .AND. JCH.LE.ICHAR('z') )
     $        STRING(ICH:ICH) = CHAR(JCH-ICHAR('a')+ICHAR('A'))
      END DO
      RETURN
      END

max_stars_repo_name: scattering/gformat

max_issues_repo_name: scattering/gformat

max_forks_repo_name: scattering/gformat


module forthread_types

use iso_c_binding
#include "fortypes.h"

end module forthread_types


text:
!===============================================================================
!> $Id: eb90519f67f73959d43f9476317a256ab739df87 $
!> This module handles buffering and output of data to disk.  It receives
!> call of the type "call output%buffer (id, it, data)", where id is the task
!> id, iout is the output number, and data is a pointer to a one-dimensional
!> buffer. The actual content of the buffers has no relevance for this module,
!> and is left to the procedure that packs the information into the buffer.
!> For simplicity, the buffer size (ndata) is assumed to remain constant.
!===============================================================================
MODULE buffered_io_mod
  USE io_mod
  USE patch_mod
  USE trace_mod
  implicit none
  private
  !-----------------------------------------------------------------------------
  ! The dump data type holds a pointer to the data to be dumped to disk, with id
  ! the task id, it the time slot, and buffer the packed task info.  It is
  ! assumed that the caller packed data into a temporarily allocated buffer,
  ! which is deallocated once the data has been written to disk.
  !-----------------------------------------------------------------------------
  type, public:: dump_t
    class(dump_t), pointer:: next => null()
    integer:: id, iout
    real(8):: time
    integer:: ndata
    real(4), dimension(:), pointer:: data
  end type
  !-----------------------------------------------------------------------------
  ! The output object contains procedures for opening, buffering, writing, and
  ! closing buffered output to a direct access file dumps.dat, with index data
  ! saved to a separate index.dat file.
  !-----------------------------------------------------------------------------
  type, public:: buffered_io_t
    integer:: ndata
    integer:: iout = -1
    integer(8):: record  = 0
    integer:: n_buffer   = 0
    integer:: max_buffer = 10
    class(dump_t), pointer:: head  => null()
    class(dump_t), pointer:: tail  => null()
  contains
    procedure:: open
    procedure, nopass:: output
    procedure:: buffer
    procedure:: check
    procedure:: write
    procedure:: close
    procedure:: find
    procedure:: test
  end type
  integer(8):: revision=1
  type(buffered_io_t), public:: buffered_io
CONTAINS

!===============================================================================
!> Write results to disk, using the buffered I/O in output_mod
!===============================================================================
SUBROUTINE output (patch)
  class(patch_t):: patch
  !.............................................................................
  type(header_t):: header
  integer:: n_buf, n_data, ibuf, it1, it2, iv
  real(4), dimension(:), pointer:: buffer
  real, dimension(:,:,:), pointer:: var
  real:: pt
  !-----------------------------------------------------------------------------
  ! Compute size of buffer and allocate -- it will be freed by the writer
  !-----------------------------------------------------------------------------
  call trace_begin ('patch_t%output_new')
  !$omp critical (buffered_cr)
  n_buf = product(patch%gn)
  n_data = n_header + patch%nv*n_buf
  allocate (buffer(n_data), var(patch%gn(1),patch%gn(2),patch%gn(3)))
  !-----------------------------------------------------------------------------
  ! Copy relevant patch info to sequenced header, and copy that to the buffer
  !-----------------------------------------------------------------------------
  call patch%patch_to_header (header)
  ibuf = 1
  call anonymous_copy (n_header, header, buffer(ibuf))
  ibuf = ibuf + n_header
  !-----------------------------------------------------------------------------
  ! Copy the variables to the output buffer; these are first interpolated in time
  !-----------------------------------------------------------------------------
  it1 = patch%iit(patch%nt-2)
  it2 = patch%iit(patch%nt-1)
  pt = (patch%out_next-patch%t(it1))/max(patch%t(it2)-patch%t(it1),1d-30)
  do iv=1,patch%nv
    var = patch%mem(:,:,:,iv,it1,1)*(1.-pt) + patch%mem(:,:,:,iv,it2,1)*pt
    call anonymous_copy (n_buf, var, buffer(ibuf))
    ibuf = ibuf + n_buf
  end do
  deallocate (var)
  call buffered_io%buffer (patch%id, patch%iout, patch%out_next, buffer)
  call buffered_io%check
  !$omp end critical (buffered_cr)
  call trace_end
END SUBROUTINE output

!===============================================================================
!> Open an index file and a data file.  The first record of the index file starts
!> with the data format revision number and the buffer size.
!===============================================================================
SUBROUTINE open (self, nwords, iout)
  class(buffered_io_t):: self
  integer:: nwords, iout
  logical, save:: first_time=.true., exist
  character(len=64):: filename
  !-----------------------------------------------------------------------------
  call trace_begin('buffered_io_mod::open')
  self%iout= iout
  self%ndata = nwords
  write (filename,'(a,i5.5,"/buffered")') trim(io%outputname), iout
  if (io%verbose>0) &
    print'(a)', 'OPEN: '//trim(filename)
  !---------------------------------------------------------------------------
  ! Check if the index file exists.  If it does, read the record counter, and
  ! if not, reset the counter and write it out.
  !---------------------------------------------------------------------------
  if (io%verbose>1) &
    print*, 'index file: ', trim(filename)//'.idx'
  inquire (file=trim(filename)//'.idx', exist=exist)
  if (exist) then
    open (io_unit%index, file=trim(filename)//'.idx', &
      form='unformatted', access='direct', status='unknown', recl=2*io%word_size)
    read (io_unit%index, rec=3) self%record
  else
    open (io_unit%index, file=trim(filename)//'.idx', &
      form='unformatted', access='direct', status='unknown', recl=2*io%word_size)
    self%record = 0
    write (io_unit%index, rec=3) self%record
  end if
  write (io_unit%index, rec=1) revision
  write (io_unit%index, rec=2) int(nwords,kind=8)
  !---------------------------------------------------------------------------
  ! Open the data file for writing
  !---------------------------------------------------------------------------
  if (io%verbose>1) &
    print*, ' dump file: ', trim(filename)//'.dat'
  open (io_unit%dump,  file=trim(filename)//'.dat', &
    form='unformatted', access='direct', status='unknown', recl=nwords*io%word_size)
  call trace_end
END SUBROUTINE open

!===============================================================================
!> Append a dump instance to a linked list.
!===============================================================================
SUBROUTINE buffer (self, id, iout, time, data)
  class(buffered_io_t):: self
  integer:: id, iout
  real(8):: time
  real(4), dimension(:), pointer:: data
  class(dump_t), pointer:: dump
  !-----------------------------------------------------------------------------
  call trace_begin('buffered_io_mod::buffer')
  allocate (dump)
  dump%id = id
  dump%iout = iout
  dump%time = time
  dump%ndata = size(data)
  dump%data => data
  if (associated(self%tail)) then
    self%tail%next => dump
  else
    self%head => dump
  end if
  self%tail => dump
  self%n_buffer = self%n_buffer+1
  if (io%verbose>1) &
    print*,'buffered_io_mod::buffer: id, iout, n_buffer=', id, iout, self%n_buffer
  call trace_end
  call self%check
END SUBROUTINE buffer

!===============================================================================
!> Check if the number of dumps has reached the maximum => write dumps to disk.
!> This is done as a task-list independent step by all threads.
!===============================================================================
SUBROUTINE check (self)
  class(buffered_io_t):: self
  class(dump_t), pointer:: dump
  !-----------------------------------------------------------------------------
  if (self%n_buffer >= self%max_buffer) then
    call self%write
  end if
END SUBROUTINE check

!===============================================================================
!> Write the buffers in the linked list to disk and deallocate the buffers.
!> The index file records, sequentially, the task id and dump number of each
!> dump written.  By reading the entire index file one can search very rapidly
!> for a specific dump, by treating the id+iout combination as an 8 byte integer
!===============================================================================
SUBROUTINE write (self)
  class (buffered_io_t):: self
  class (dump_t), pointer:: dump, old
  !-----------------------------------------------------------------------------
  call trace_begin('buffered_io_mod::write')
  dump => self%head
  if (io%verbose>0) &
    print'(a,i7,i6,i9,g15.6)', 'output: it, nbuf, ndata, time', &
        dump%id, self%n_buffer, dump%ndata, dump%time
  do while (associated(dump))
    !---------------------------------------------------------------------------
    ! If time slot differs, close the previous one and open the current one
    !---------------------------------------------------------------------------
    if (dump%iout /= self%iout) then
      if (self%iout >= 0) call self%close
      call self%open (dump%ndata, dump%iout)
    end if
    !---------------------------------------------------------------------------
    ! Increment the record counter, record its value, and write the data
    !---------------------------------------------------------------------------
    !$omp atomic
    self%record = self%record+1
    write (io_unit%index, rec=3) self%record
    write (io_unit%index, rec=self%record+3) dump%id, dump%iout
    flush (io_unit%index)
    write (io_unit%dump, rec=self%record) dump%data
    flush (io_unit%dump)
    !$omp atomic
    self%n_buffer = self%n_buffer-1
    if (io%verbose>1) &
      print'(a,i7,i6,i9,g15.6)', 'output: id, it, ndata, time', &
        dump%id, dump%iout, dump%ndata, dump%time
    !---------------------------------------------------------------------------
    ! Delete the buffer data, and the buffer
    !---------------------------------------------------------------------------
    old => dump
    dump => dump%next
    deallocate (old%data)
    deallocate (old)
  end do
  nullify (self%head)
  nullify (self%tail)
  call trace_end
END SUBROUTINE write

!===============================================================================
!> Close the files, writing out the total number of records in slot 3 of the
!> index file.
!===============================================================================
SUBROUTINE close (self)
  class (buffered_io_t):: self
  !-----------------------------------------------------------------------------
  if (io%do_legacy) return
  call trace_begin('buffered_io_mod::close')
  close (io_unit%index)
  close (io_unit%dump)
  call trace_end
END SUBROUTINE close

!===============================================================================
!> Find the record that contains the combination of task id and dump number iout.
!> This will typically be done in Python or IDL; this is just a demo.
!===============================================================================
FUNCTION find (self, id, iout)
  class (buffered_io_t):: self
  integer, optional:: id, iout
  integer:: find, loc(1)
  integer(8), dimension(:), pointer:: index
  integer(8):: key, ndata, nrecord, revision, i
  !-----------------------------------------------------------------------------
  call trace_begin('buffered_io_mod::find')
  !-----------------------------------------------------------------------------
  ! Read the number of records and the buffer size
  !-----------------------------------------------------------------------------
  open (io_unit%index, file=trim(io%outputname)//'index.dat', &
    form='unformatted', access='direct', status='old', recl=8)
  read (io_unit%index,rec=1) revision
  read (io_unit%index,rec=2) ndata
  read (io_unit%index,rec=3) nrecord
  print *, revision, ndata, nrecord
  close (io_unit%index)
  !-----------------------------------------------------------------------------
  ! Get the entire index in one read
  !-----------------------------------------------------------------------------
  allocate (index(nrecord))
  open (io_unit%index, file=trim(io%outputname)//'index.dat', &
    form='unformatted', access='direct', status='old', recl=8*(nrecord+3))
  read (io_unit%index,rec=1) revision, ndata, nrecord, index
  close (io_unit%index)
  !-----------------------------------------------------------------------------
  ! Find the id+iout combo
  !-----------------------------------------------------------------------------
  key = int(id,kind=8) + int(iout,kind=8)*2_8**32
  if (io%verbose>1) &
    print*, 'searching for id, iout, key =', id, iout, key
  if (io%verbose>2) then
    do i=1,nrecord
      print *, i, index(i)
    end do
  end if
  loc = minloc(abs(index-key))
  find = loc(1)
  if (io%verbose>1) &
    print*, 'record =', find
  deallocate (index)
  call trace_end
END FUNCTION find

!===============================================================================
!> Test the output module, by buffering enough dumps for a self-triggered write,
!> followed by a write triggered by a close.
!===============================================================================
SUBROUTINE test (self)
  class (buffered_io_t):: self
  real(4), dimension(:), pointer:: data
  integer, dimension(:), pointer:: p
  integer:: n=20, nv=8, id, it, ndata, record
  real(8):: time
  !-----------------------------------------------------------------------------
  call trace_begin ('buffered_io_mod::test')
  ndata = n**3*nv+n**2
  call self%open (ndata,0)
  do it=1,3
    time = it-1d0
    do id=1,9
      allocate (data(ndata))
      call self%buffer (id, it, time, data)
    end do
  end do
  call self%close
  record = self%find (id=2,iout=3)
  print *,'record =', record
  call trace_end
END SUBROUTINE test

END MODULE buffered_io_mod

max_stars_repo_name: applejwjcat/dispatch

max_issues_repo_name: applejwjcat/dispatch

max_forks_repo_name: applejwjcat/dispatch


text:
! THIS VERSION: GALAHAD 2.1 - 22/03/2007 AT 09:00 GMT.

!-*-*-*-*-*-*-*-*-*- G A L A H A D _ R P D   M O D U L E -*-*-*-*-*-*-*-*-

!  Copyright reserved, Gould/Orban/Toint, for GALAHAD productions
!  Principal author: Nick Gould

!  History -
!   originally released with GALAHAD Version 2.0. January 22nd 2006

!  For full documentation, see
!   http://galahad.rl.ac.uk/galahad-www/specs.html

   MODULE GALAHAD_RPD_double

!     ----------------------------------------------
!     |                                            |
!     | Read data for the quadratic program        |
!     |                                            |
!     |    minimize     1/2 x(T) H x + g(T) x + f  |
!     |    subject to     c_l <= A x <= c_u        |
!     |                   x_l <=  x  <= x_u        |
!     |                                            |
!     | from a problem-data file                   |
!     |                                            |
!     ----------------------------------------------

      USE GALAHAD_SMT_double, only: SMT_put
      USE GALAHAD_QPT_double

      IMPLICIT NONE

      PRIVATE
      PUBLIC :: RPD_read_problem_data, QPT_problem_type

!--------------------
!   P r e c i s i o n
!--------------------

      INTEGER, PARAMETER :: wp = KIND( 1.0D+0 )

!----------------------
!   P a r a m e t e r s
!----------------------

      REAL ( KIND = wp ), PARAMETER :: zero = 0.0_wp
      INTEGER, PARAMETER :: input_line_length = 256

!-------------------------------------------------
!  D e r i v e d   t y p e   d e f i n i t i o n s
!-------------------------------------------------

      TYPE, PUBLIC :: RPD_inform_type
        INTEGER :: status, alloc_status, io_status, line
        CHARACTER ( LEN = 1 ) :: bad_alloc
      END TYPE

   CONTAINS

!-*-*-   R P D _ R E A D _ P R O B L E M _ D A T A   S U B R O U T I N E   -*-*-

      SUBROUTINE RPD_read_problem_data( input, prob, inform )
      INTEGER, INTENT( IN ) :: input
      TYPE ( QPT_problem_type ), INTENT( INOUT ) :: prob
      TYPE ( RPD_inform_type ), INTENT( out ) :: inform

!  Read the problem-data file from unit input into the derived type prob

!  ****************************************************************************

!  The data should be input in a file on unit 5. The data is in
!  free format (blanks separate values),but must occur in the order
!  given here with NO SEPARATING BLANK LINES. Each term in "quotes"
!  denotes a required value. Any strings beyond those required on a
!  given lines will be regarded as comments and ignored.

!  "problem name"
!  "number varibales"
!  "number general linear constraints"
!  "number of nonzeros in lower traingle of H"
!  "row" "column" "value" for each entry of H (if any), one triple on each line
!  "default value for entries in g"
!  "number of non-default entries in g"
!  "index" "value" for each non-default term in g (if any), one pair per line
!  "value of f"
!  "number of nonzeros in A"
!  "row" "column" "value" for each entry of A (if any), one triple on each line
!  "default value for entries in c_l"
!  "number of non-default entries in c_l"
!  "index" "value" for each non-default term in c_l (if any), one pair per line
!  "default value for entries in c_u"
!  "number of non-default entries in c_u"
!  "index" "value" for each non-default term in c_u (if any), one pair per line
!  "default value for entries in x_l"
!  "number of non-default entries in x_l"
!  "index" "value" for each non-default term in x_l (if any), one pair per line
!  "default value for entries in x_u"
!  "number of non-default entries in x_u"
!  "index" "value" for each non-default term in x_u (if any), one pair per line
!  "default value for starting value for variables x"
!  "number of non-default starting entries in x"
!  "index" "value" for each non-default term in x (if any), one pair per line
!  "default value for starting value for Lagrange multipliers y for constraints"
!  "number of non-default starting entries in y"
!  "index" "value" for each non-default term in y (if any), one pair per line
!  "default value for starting value for dual varibales z for simple bounds"
!  "number of non-default starting entries in z"
!  "index" "value" for each non-default term in z (if any), one pair per line

!  *****************************************************************************

!  Local variables

     INTEGER :: i, j, k, A_ne, H_ne, nnzx_0, nnzy_0, nnzz_0
     INTEGER :: nnzg, nnzc_l, nnzc_u, nnzx_l, nnzx_u, smt_stat
     REAL ( KIND = wp ) :: rv, default
     CHARACTER ( LEN = 10 ) :: pname
     CHARACTER ( LEN = input_line_length ) :: input_line, blank_line

     inform%line = 0
     inform%alloc_status = 0 ; inform%bad_alloc = ' '

     DO i = 1, input_line_length
       blank_line( i : i ) = ' '
     END DO

!  Determine the problem name

     pname = '          '
     DO
       inform%line = inform%line + 1
       input_line = blank_line
       READ( input, "( A )", END = 930, ERR = 940 ) input_line
       IF ( .NOT. RPD_ignore_string( input_line ) ) EXIT
     END DO
     READ( input_line, *, IOSTAT = inform%io_status,                           &
           END = 930, ERR = 940 ) pname
     ALLOCATE( prob%name( 10 ) )
     prob%name = TRANSFER( pname, prob%name )

!  Determine the number of variables and constraints

     DO
       inform%line = inform%line + 1
       input_line = blank_line
       READ( input, "( A )", END = 930, ERR = 940 ) input_line
       IF ( .NOT. RPD_ignore_string( input_line ) ) EXIT
     END DO
     READ( input_line, *, IOSTAT = inform%io_status,                           &
           END = 930, ERR = 940 ) prob%n
     DO
       inform%line = inform%line + 1
       input_line = blank_line
       READ( input, "( A )", END = 930, ERR = 940 ) input_line
       IF ( .NOT. RPD_ignore_string( input_line ) ) EXIT
     END DO
     READ( input_line, *, IOSTAT = inform%io_status,                           &
           END = 930, ERR = 940 ) prob%m

!  Allocate suitable arrays

     ALLOCATE( prob%X( prob%n ), prob%X_l( prob%n ), prob%X_u( prob%n ),       &
               prob%G( prob%n ), prob%Z( prob%n ), STAT = inform%alloc_status )
     IF ( inform%alloc_status /= 0 ) THEN
       inform%status = - 2 ; inform%bad_alloc = 'X' ; RETURN
     END IF

     ALLOCATE( prob%C_l( prob%m ), prob%C_u( prob%m ), prob%Y( prob%m ),       &
               prob%C( prob%m ), STAT = inform%alloc_status )
     IF ( inform%alloc_status /= 0 ) THEN
       inform%status = - 2 ; inform%bad_alloc = 'Y' ; RETURN
     END IF

!  Fill component H

     DO
       inform%line = inform%line + 1
       input_line = blank_line
       READ( input, "( A )", END = 930, ERR = 940 ) input_line
       IF ( .NOT. RPD_ignore_string( input_line ) ) EXIT
     END DO

     READ( input_line, *, IOSTAT = inform%io_status,                           &
           END = 930, ERR = 940 ) prob%H%ne
     ALLOCATE( prob%H%row( prob%H%ne ), prob%H%col( prob%H%ne ),               &
               prob%H%val( prob%H%ne ), STAT = inform%alloc_status )
     IF ( inform%alloc_status /= 0 ) THEN
       inform%status = - 2 ; inform%bad_alloc = 'H' ; RETURN
     END IF

     H_ne = 0
     DO k = 1, prob%H%ne
       DO
         inform%line = inform%line + 1
         input_line = blank_line
         READ( input, "( A )", END = 930, ERR = 940 ) input_line
         IF ( .NOT. RPD_ignore_string( input_line ) ) EXIT
       END DO
       READ( input_line, *, IOSTAT = inform%io_status,                         &
             END = 930, ERR = 940 ) i, j, rv
       IF ( rv == zero ) CYCLE
       H_ne = H_ne + 1 ; prob%H%val( H_ne ) = rv
       IF ( i >= j ) THEN
         prob%H%row( H_ne ) = i
         prob%H%col( H_ne ) = j
       ELSE
         prob%H%row( H_ne ) = j
         prob%H%col( H_ne ) = i
       END IF
     END DO
     prob%H%ne = H_ne
     IF ( ALLOCATED( prob%H%type ) ) DEALLOCATE( prob%H%type )
     CALL SMT_put( prob%H%type, 'COORDINATE', smt_stat )

!  Fill component g

     DO
       inform%line = inform%line + 1
       input_line = blank_line
       READ( input, "( A )", END = 930, ERR = 940 ) input_line
       IF ( .NOT. RPD_ignore_string( input_line ) ) EXIT
     END DO
     READ( input_line, *, IOSTAT = inform%io_status,                           &
           END = 930, ERR = 940 ) default
     prob%G = default
     DO
       inform%line = inform%line + 1
       input_line = blank_line
       READ( input, "( A )", END = 930, ERR = 940 ) input_line
       IF ( .NOT. RPD_ignore_string( input_line ) ) EXIT
     END DO
     READ( input_line, *, IOSTAT = inform%io_status,                           &
           END = 930, ERR = 940 ) nnzg
     DO k = 1, nnzg
       DO
         inform%line = inform%line + 1
         input_line = blank_line
         READ( input, "( A )", END = 930, ERR = 940 ) input_line
         IF ( .NOT. RPD_ignore_string( input_line ) ) EXIT
       END DO
       READ( input_line, *, IOSTAT = inform%io_status,                         &
             END = 930, ERR = 940 ) i, rv
       prob%G( i ) = rv
     END DO

!  Fill component f

     DO
       inform%line = inform%line + 1
       input_line = blank_line
       READ( input, "( A )", END = 930, ERR = 940 ) input_line
       IF ( .NOT. RPD_ignore_string( input_line ) ) EXIT
     END DO
     READ( input_line, *, IOSTAT = inform%io_status,                           &
           END = 930, ERR = 940 ) prob%f

!  Fill component A

     DO
       inform%line = inform%line + 1
       input_line = blank_line
       READ( input, "( A )", END = 930, ERR = 940 ) input_line
       IF ( .NOT. RPD_ignore_string( input_line ) ) EXIT
     END DO
     READ( input_line, *, IOSTAT = inform%io_status,                           &
           END = 930, ERR = 940 ) prob%A%ne
     ALLOCATE( prob%A%row( prob%A%ne ), prob%A%col( prob%A%ne ),               &
               prob%A%val( prob%A%ne ), STAT = inform%alloc_status )
     IF ( inform%alloc_status /= 0 ) THEN
       inform%status = - 2 ; inform%bad_alloc = 'A' ; RETURN
     END IF

     A_ne = 0
     DO k = 1, prob%A%ne
       DO
         inform%line = inform%line + 1
         input_line = blank_line
         READ( input, "( A )", END = 930, ERR = 940 ) input_line
         IF ( .NOT. RPD_ignore_string( input_line ) ) EXIT
       END DO
       READ( input_line, *, IOSTAT = inform%io_status,                         &
             END = 930, ERR = 940 ) i, j, rv
       IF ( rv == zero ) CYCLE
       A_ne = A_ne + 1 ; prob%A%val( A_ne ) = rv
       prob%A%row( A_ne ) = i ; prob%A%col( A_ne ) = j
     END DO
     prob%A%ne = A_ne
     IF ( ALLOCATED( prob%A%type ) ) DEALLOCATE( prob%A%type )
     CALL SMT_put( prob%A%type, 'COORDINATE', smt_stat )

!  Fill component c_l

     DO
       inform%line = inform%line + 1
       input_line = blank_line
       READ( input, "( A )", END = 930, ERR = 940 ) input_line
       IF ( .NOT. RPD_ignore_string( input_line ) ) EXIT
     END DO
     READ( input_line, *, IOSTAT = inform%io_status,                           &
           END = 930, ERR = 940 ) default
     prob%C_l = default
     DO
       inform%line = inform%line + 1
       input_line = blank_line
       READ( input, "( A )", END = 930, ERR = 940 ) input_line
       IF ( .NOT. RPD_ignore_string( input_line ) ) EXIT
     END DO
     READ( input_line, *, IOSTAT = inform%io_status,                           &
           END = 930, ERR = 940 ) nnzc_l
     DO k = 1, nnzc_l
       DO
         inform%line = inform%line + 1
         input_line = blank_line
         READ( input, "( A )", END = 930, ERR = 940 ) input_line
         IF ( .NOT. RPD_ignore_string( input_line ) ) EXIT
       END DO
       READ( input_line, *, IOSTAT = inform%io_status,                         &
             END = 930, ERR = 940 ) i, rv
       prob%C_l( i ) = rv
     END DO

!  Fill component c_u

     DO
       inform%line = inform%line + 1
       input_line = blank_line
       READ( input, "( A )", END = 930, ERR = 940 ) input_line
       IF ( .NOT. RPD_ignore_string( input_line ) ) EXIT
     END DO
     READ( input_line, *, IOSTAT = inform%io_status,                           &
           END = 930, ERR = 940 ) default
     prob%C_u = default
     DO
       inform%line = inform%line + 1
       input_line = blank_line
       READ( input, "( A )", END = 930, ERR = 940 ) input_line
       IF ( .NOT. RPD_ignore_string( input_line ) ) EXIT
     END DO
     READ( input_line, *, IOSTAT = inform%io_status,                           &
           END = 930, ERR = 940 ) nnzc_u
     DO k = 1, nnzc_u
       DO
         inform%line = inform%line + 1
         input_line = blank_line
         READ( input, "( A )", END = 930, ERR = 940 ) input_line
         IF ( .NOT. RPD_ignore_string( input_line ) ) EXIT
       END DO
       READ( input_line, *, IOSTAT = inform%io_status,                         &
             END = 930, ERR = 940 ) i, rv
       prob%C_u( i ) = rv
     END DO

!  Fill component x_l

     DO
       inform%line = inform%line + 1
       input_line = blank_line
       READ( input, "( A )", END = 930, ERR = 940 ) input_line
       IF ( .NOT. RPD_ignore_string( input_line ) ) EXIT
     END DO
     READ( input_line, *, IOSTAT = inform%io_status,                           &
           END = 930, ERR = 940 ) default
     prob%X_l = default
     DO
       inform%line = inform%line + 1
       input_line = blank_line
       READ( input, "( A )", END = 930, ERR = 940 ) input_line
       IF ( .NOT. RPD_ignore_string( input_line ) ) EXIT
     END DO
     READ( input_line, *, IOSTAT = inform%io_status,                           &
           END = 930, ERR = 940 ) nnzx_l
     DO k = 1, nnzx_l
       DO
         inform%line = inform%line + 1
         input_line = blank_line
         READ( input, "( A )", END = 930, ERR = 940 ) input_line
         IF ( .NOT. RPD_ignore_string( input_line ) ) EXIT
       END DO
       READ( input_line, *, IOSTAT = inform%io_status,                         &
             END = 930, ERR = 940 ) i, rv
       prob%X_l( i ) = rv
     END DO

!  Fill component x_u

     DO
       inform%line = inform%line + 1
       input_line = blank_line
       READ( input, "( A )", END = 930, ERR = 940 ) input_line
       IF ( .NOT. RPD_ignore_string( input_line ) ) EXIT
     END DO
     READ( input_line, *, IOSTAT = inform%io_status,                           &
           END = 930, ERR = 940 ) default
     prob%X_u = default
     DO
       inform%line = inform%line + 1
       input_line = blank_line
       READ( input, "( A )", END = 930, ERR = 940 ) input_line
       IF ( .NOT. RPD_ignore_string( input_line ) ) EXIT
     END DO
     READ( input_line, *, IOSTAT = inform%io_status,                           &
           END = 930, ERR = 940 ) nnzx_u
     DO k = 1, nnzx_u
       DO
         inform%line = inform%line + 1
         input_line = blank_line
         READ( input, "( A )", END = 930, ERR = 940 ) input_line
         IF ( .NOT. RPD_ignore_string( input_line ) ) EXIT
       END DO
       READ( input_line, *, IOSTAT = inform%io_status,                         &
             END = 930, ERR = 940 ) i, rv
       prob%X_u( i ) = rv
     END DO

!  Fill component x

     DO
       inform%line = inform%line + 1
       input_line = blank_line
       READ( input, "( A )", END = 930, ERR = 940 ) input_line
       IF ( .NOT. RPD_ignore_string( input_line ) ) EXIT
     END DO
     READ( input_line, *, IOSTAT = inform%io_status,                           &
           END = 930, ERR = 940 ) default
     prob%X = default
     DO
       inform%line = inform%line + 1
       input_line = blank_line
       READ( input, "( A )", END = 930, ERR = 940 ) input_line
       IF ( .NOT. RPD_ignore_string( input_line ) ) EXIT
     END DO
     READ( input_line, *, IOSTAT = inform%io_status,                           &
           END = 930, ERR = 940 ) nnzx_0
     DO k = 1, nnzx_0
       DO
         inform%line = inform%line + 1
         input_line = blank_line
         READ( input, "( A )", END = 930, ERR = 940 ) input_line
         IF ( .NOT. RPD_ignore_string( input_line ) ) EXIT
       END DO
       READ( input_line, *, IOSTAT = inform%io_status,                         &
             END = 930, ERR = 940 ) i, rv
       prob%X( i ) = rv
     END DO

!  Fill component y

     DO
       inform%line = inform%line + 1
       input_line = blank_line
       READ( input, "( A )", END = 930, ERR = 940 ) input_line
       IF ( .NOT. RPD_ignore_string( input_line ) ) EXIT
     END DO
     READ( input_line, *, IOSTAT = inform%io_status,                           &
           END = 930, ERR = 940 ) default
     prob%Y = default
     DO
       inform%line = inform%line + 1
       input_line = blank_line
       READ( input, "( A )", END = 930, ERR = 940 ) input_line
       IF ( .NOT. RPD_ignore_string( input_line ) ) EXIT
     END DO
     READ( input_line, *, IOSTAT = inform%io_status,                           &
           END = 930, ERR = 940 ) nnzy_0
     DO k = 1, nnzy_0
       DO
         inform%line = inform%line + 1
         input_line = blank_line
         READ( input, "( A )", END = 930, ERR = 940 ) input_line
         IF ( .NOT. RPD_ignore_string( input_line ) ) EXIT
       END DO
       READ( input_line, *, IOSTAT = inform%io_status,                         &
             END = 930, ERR = 940 ) i, rv
       prob%Y( i ) = rv
     END DO

!  Fill component z

     DO
       inform%line = inform%line + 1
       input_line = blank_line
       READ( input, "( A )", END = 930, ERR = 940 ) input_line
       IF ( .NOT. RPD_ignore_string( input_line ) ) EXIT
     END DO
     READ( input_line, *, IOSTAT = inform%io_status,                           &
           END = 930, ERR = 940 ) default
     prob%Z = default
     DO
       inform%line = inform%line + 1
       input_line = blank_line
       READ( input, "( A )", END = 930, ERR = 940 ) input_line
       IF ( .NOT. RPD_ignore_string( input_line ) ) EXIT
     END DO
     READ( input_line, *, IOSTAT = inform%io_status,                           &
           END = 930, ERR = 940 ) nnzz_0
     DO k = 1, nnzz_0
       DO
         inform%line = inform%line + 1
         input_line = blank_line
         READ( input, "( A )", END = 930, ERR = 940 ) input_line
         IF ( .NOT. RPD_ignore_string( input_line ) ) EXIT
       END DO
       READ( input_line, *, IOSTAT = inform%io_status,                         &
             END = 930, ERR = 940 ) i, rv
       prob%Z( i ) = rv
     END DO

     inform%status = 0
     RETURN

!  Error returns

!  - end of file encountered

 930 CONTINUE
     inform%status = 3
     RETURN

!  - other error encountered

 940 CONTINUE
     inform%status = 4
     RETURN

!  End of RPD_read_problem_data

     END SUBROUTINE RPD_read_problem_data

!-*-*-*-*-*-   R P D _ I G N O R E _ S T R I N G   F U N C T I O N   -*-*-*-*-*-

     FUNCTION RPD_ignore_string( input_line )
     LOGICAL :: RPD_ignore_string

!  Ignore a string if it is (a) blank or (b) starts with "!", "%" or "#"

     CHARACTER ( LEN = input_line_length ), INTENT( IN ) :: input_line

!  Local variables

     INTEGER :: i, length_string

     length_string = LEN_TRIM( input_line )
     IF ( length_string <= 0 ) THEN
       RPD_ignore_string = .TRUE.
       RETURN
     END IF

     DO i = 1, length_string
       IF ( input_line( i : i ) == ' ' ) CYCLE
       IF ( input_line( i : i ) == '!' .OR. input_line( i : i ) == '#' .OR. &
            input_line( i : i ) == '%' .OR. input_line( i : i ) == '|' ) THEN
         RPD_ignore_string = .TRUE.
         RETURN
       END IF
       EXIT
     END DO
     RPD_ignore_string = .FALSE.

     RETURN

!  End of RPD_ignore_string

     END FUNCTION RPD_ignore_string

!  End of module RPD

   END MODULE GALAHAD_RPD_double

max_forks_repo_name: dynamics-of-stellar-systems/dynamite


text:
      MODULE nine0_I
      INTERFACE
!...Generated by Pacific-Sierra Research 77to90 4.3E  10:07:58  11/16/01
      SUBROUTINE nine0 (J1, J2, J3, L1, L2, L3, K1, K2, K3, AA)
      USE vast_kind_param,ONLY: DOUBLE
      INTEGER, INTENT(IN) :: J1
      INTEGER, INTENT(IN) :: J2
      INTEGER, INTENT(IN) :: J3
      INTEGER, INTENT(IN) :: L1
      INTEGER, INTENT(IN) :: L2
      INTEGER, INTENT(IN) :: L3
      INTEGER, INTENT(IN) :: K1
      INTEGER, INTENT(IN) :: K2
      INTEGER, INTENT(IN) :: K3
      REAL(DOUBLE), INTENT(OUT) :: AA
!VAST.../CONSTS/ ZERO(IN), ONE(IN)
!VAST...Calls: SIXJ
      END SUBROUTINE
      END INTERFACE
      END MODULE

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: sylas/grasp-continuum

max_forks_repo_name: sylas/grasp-continuum


text:
!******************************************************************************
!
! Licensing:
!
!   This code is distributed under the MIT license.
!
! Modified:
!
!   22 April 2019
!
! Author:
!
!   Nescirem
!
!==============================================================================
!
!>M mpi_mod
!

    module mpi_mod

        implicit none

        !use mpi
        include 'mpif.h'

        integer             :: required,provided
        integer,parameter   :: root = 0
        integer             :: istat( mpi_status_size )
        integer             :: num_p,num_slave
        integer             :: err
        integer             :: pid,myid,myleft,myright
        integer             :: p_namelen
        character(len=64)   :: p_name

    end module mpi_mod

max_stars_repo_name: nescirem/nesDockerfile

max_issues_repo_name: nescirem/nesDockerfile

max_forks_repo_name: nescirem/nesDockerfile


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

integer, PARAMETER :: AAA(3:4) = (/ -325, 400/)
integer , PARAMETER :: YYY(3:3) = (/ (abs((AAA(i))), i = 3,3) /)
integer , PARAMETER :: ZZZ = abs(aaa(3)) !!! works
integer pass

pass = 1;
if (zzz .eq. 325) then
  do i = 3, 3
    if (yyy(i) .ne. 325) then
      pass = 0
    end if
  end do
else
  pass = 0
endif

if (pass.eq.0) then
  print *, "FAIL"
else
  print *, "PASS"
endif

end

max_stars_repo_name: abrahamtovarmob/flang

max_issues_repo_name: abrahamtovarmob/flang

max_forks_repo_name: abrahamtovarmob/flang


text:
Filthy Pete and the Moustaches have two official websites:
  http://profiles.myspace.com/users/1944079
  http://www.livejournal.com/users/filthypetemusic/

Users/StanOklobdzija says this about the band: Filthy Pete and the Moustaches have invented rock n roll in order to destroy it. Born of childhood friends Jeff Robinson and Chris Caskey, aka Filthy Pete and Strummin Joe Moustache, Filthy Pete and the Moustaches combine rock and roll guitars that belong in the 50s with electronic drums that belongs in the future. Their sound is best described as the auditory mulatto of Elvis Presley and Alec Empire. In May, the band shared the podium for first place at Battle of the Bands, only to have the title stripped from them in a questionable game of roshambo. The ensuing bout of depression only further fueled the bands creative genius, leading Robinson and Caskey to churn out a flurry of new material. The result shall be an exciting year for anyone within 100 miles of their modest residence in east Davis who yearns to experience the future of music today.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
passive component C {
  enum E { X, Y } default "abc"
}

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


text:
      program type25

C     Check multi-type initialization routine

      character*20 s

      call init(s, x, i)

      print *, s, x, i

      end

      subroutine init(fs, fx, fi)

      character*(*) fs
      integer fi

      fs = "Hello World!"
      fx = 2.0
      fi = 1

      end


text:
| on linux ...

6 func: mmap | addr len prot flags fd offset
2 func: munmap

: memmap ( fd n -- b n )
	swap >r					| n r:fd
	0 over					| n 0 n r:h
	3 1 r@ 0
	mmap swap
	r> close
	;
: slurp2 ( a n -- b m )
	open/rw dup fsize
	memmap
	;

max_stars_repo_name: ronaaron/reva

max_issues_repo_name: ronaaron/reva

max_forks_repo_name: ronaaron/reva


text:
      SUBROUTINE DL7NVR(N, LIN, L)
C
C  ***  COMPUTE  LIN = L**-1,  BOTH  N X N  LOWER TRIANG. STORED   ***
C  ***  COMPACTLY BY ROWS.  LIN AND L MAY SHARE THE SAME STORAGE.  ***
C
C  ***  PARAMETERS  ***
C
      INTEGER N
      DOUBLE PRECISION L(1), LIN(1)
C     DIMENSION L(N*(N+1)/2), LIN(N*(N+1)/2)
C
C  ***  LOCAL VARIABLES  ***
C
      INTEGER I, II, IM1, JJ, J0, J1, K, K0, NP1
      DOUBLE PRECISION ONE, T, ZERO
C/6
C     DATA ONE/1.D+0/, ZERO/0.D+0/
C/7
      PARAMETER (ONE=1.D+0, ZERO=0.D+0)
C/
C
C  ***  BODY  ***
C
      NP1 = N + 1
      J0 = N*(NP1)/2
      DO 30 II = 1, N
         I = NP1 - II
         LIN(J0) = ONE/L(J0)
         IF (I .LE. 1) GO TO 999
         J1 = J0
         IM1 = I - 1
         DO 20 JJ = 1, IM1
              T = ZERO
              J0 = J1
              K0 = J1 - JJ
              DO 10 K = 1, JJ
                   T = T - L(K0)*LIN(J0)
                   J0 = J0 - 1
                   K0 = K0 + K - I
 10                CONTINUE
              LIN(J0) = T/L(K0)
 20           CONTINUE
         J0 = J0 - 1
 30      CONTINUE
 999  RETURN
C  ***  LAST CARD OF DL7NVR FOLLOWS  ***
      END

max_stars_repo_name: DanIverson/OpenVnmrJ

max_issues_repo_name: timburrow/openvnmrj-source

max_forks_repo_name: timburrow/openvnmrj-source


text:
!! R1118 parent-identifier
!     is ancestor-module-name [ : parent-submodule-name ]
!
module ancestor
  contains
end

11 submodule (ancestor) parent
  integer :: i
end submodule

SUBMODULE (ancestor:parent) subby
  integer :: j
13 END SUBMODULE subby

max_stars_repo_name: OpenFortranProject/ofp-sdf

max_issues_repo_name: OpenFortranProject/ofp-sdf

max_forks_repo_name: OpenFortranProject/ofp-sdf


! { dg-do run }
! Test of the fix to the bug triggered by NIST fm908.for.
! Contributed by Paul Thomas  <pault@gcc.gnu.org>
!
program past_eor
  character(len=82)         ::  buffer
  real                      ::  a(2), b(2), c(2), d(2), e(2)

  e = (/2.34,2.456/)

! tests 28-31 from fm908.for

  buffer = '  2.34 ,  2.456     2.34 ,  2.456     0.234E01,  2.456E00&
 &   0.234E+001, 2.456E-000'

  READ (UNIT=buffer,FMT=10) a, b, c, d
10 FORMAT (2(2(G7.5,1X),2X),2(G10.4E2,1X),1X,2(G11.7E4,1X))

  if (any (a.ne.e).or.any (b.ne.e).or.any (c.ne.e).or.any (d.ne.e)) STOP 1

end program past_eor


text:
SUBROUTINE radint(JATOM,is)
  use potnlc
  use intu
  ! The major and the minor components of the Dirac solution are u=(g/r, f/r).
  !   It turns out that f \approx (dg/dr-g/r)/2 = (du/dr)*r/2 and hence du/dr = 2*f/r
  !
  !   It then follows: <u_1|d/dr|u_2> = Int[ g_1/r  2/r f_2 r^2  dr] = 2 Int[ g_1 f_2 dr]   == 2<g_1|f_2>
  !                    <u_1|1/r |u_2> = Int[ g_1/r  1/r g_2/r r^2 dr] = Int[ g_1 g_2 /r dr] == <g_1/sqrt(r)|g_2/sqrt(r)>
  !
  !  The following components are used
  !     u   <- g   | up  <- dotg
  !     du  <- f   | dup <- dotf
  !     ul  <- go  for localized orbital
  !     dul <- fo  for localized orbital
  !ad  calculation of radial integrals:
  !ad  U   = u
  !ad  DU  = du/dr
  !ad  UP  = du/dE
  !ad  DUP = dUP/dr
  !ad
  !ad  for LO's:
  !ad  UL  = u_lo
  !ad  DUL = du_lo/dr
  !ad
  !ad  the radial functions are calculated in ATPAR and LOMAIN
  !ad
  INCLUDE 'param.inc'
  !
  IMPLICIT REAL*8 (A-H,O-Z)
  LOGICAL  lapw(0:lmax2),loor(0:lomax),lloor(0:lmax2)
  !ad
  COMMON /XA/     R(NRAD),BK(3)
  COMMON /RADFU/  U(NRAD,LMAX1),UP(NRAD,LMAX1),DU(NRAD,LMAX1),DUP(NRAD,LMAX1)
  COMMON /UHELP/  UA(NRAD), UB(NRAD), UAB(NRAD), UBA(NRAD), UAA(NRAD), UBB(NRAD)
  common /loabc/  alo(0:lomax,nloat),blo(0:lomax,nloat),clo(0:lomax,nloat),elo(0:lomax,nloat),plo(0:lomax),dplo(0:lomax),pelo(0:lomax),dpelo(0:lomax),peilo(0:lomax),pi12lo(0:lomax),pe12lo(0:lomax),UL(nrad,lomax+1),DUL(nrad,lomax+1)
  common /lolog/  nlo,nlov,nlon,lapw,ilo(0:lomax),loor,lloor
!ad..................calculation of radial integrals....................
  TWO=2.d0
  IMAX=JRI(JATOM)
  DO I=1,IMAX ! 2
     R(I)=Rnot(JATOM)*EXP((I-1)*DX(JATOM))
  ENDDO       ! 2

!!! BRISI
  !if (is.eq.1) then
  !   open(991,FILE='radials_up.dat')
  !else
  !   open(991,FILE='radials_dn.dat')
  !endif
  !do i=1,imax
  !   WRITE(991,'(7f12.7)') R(i), U(i,1), U(i,2), U(i,3), dU(i,1), dU(i,2), dU(i,3)
  !enddo
  !close(991)
!!! BRISI

  !ad
  !ad..................RUU(L)..RUUP(L)..RUPU(L)..RUPUP(L).................
  !ad
  DO L=1,LMAX     !103
     DO I=1,IMAX  !105
        SR=SQRT(R(I))
        UA(I) =U(I,L)/SR       ! g_l/sqrt(r)
        UB(I) =U(I,L+1)/SR     ! g_{l+1}/sqrt(r)
        UAB(I)=UP(I,L)/SR      ! dotg_l/sqrt(r)
        UBA(I)=UP(I,L+1)/SR    ! dotg_{l+1}/sqrt(r)
     ENDDO
     ! Kristjan : This can not be correct. I think minor and major component are here taken to be the same
     !            I think one should get correct relativistic components from previous routines.
     !            and set .FALSE. to .TRUE. Right now minor components are neglected.
      CALL RINT13(.FALSE.,UA,UA,UB,UB,RUU(L,JATOM,is),JATOM)      ! ruu  =        <u_l|1/r|u_{l+1}>
      CALL RINT13(.FALSE.,UB,UB,UAB,UAB,RUUP(L,JATOM,is),JATOM)   ! ruup =    <u_{l+1}|1/r|dotu_l>
      CALL RINT13(.FALSE.,UBA,UBA,UA,UA,RUPU(L,JATOM,is),JATOM)   ! rupu = <dotu_{l+1}|1/r|u_l>
      CALL RINT13(.FALSE.,UBA,UBA,UAB,UAB,RUPUP(L,JATOM,is),JATOM)! rupup= <dotu_{l+1}|1/r|dotu_l>
   ENDDO         !103
   !
   !ad
   !ad...............DUU1(L)..DUUP1(L)..DUPU1(L)..DUPUP1(L)................
   !ad
   DO L=1,LMAX       ! 203
      DO I=1,IMAX    ! 205
         UA(I) =U(I,L+1)       ! g_{l+1}
         UB(I) =DU(I,L)*TWO    ! 2*f_{l}
         UAB(I)=DUP(I,L)*TWO   ! 2*dotf_l
         UBA(I)=UP(I,L+1)      ! dotg_{l+1}
      ENDDO          ! 205
      CALL RINT13(.FALSE.,UA,UA,UB,UB,DUU1(L,JATOM,is),JATOM)      ! duu1  = <g_{l+1}|2*f_l>       ==    <u_{l+1}|d/dr|u_l>
      CALL RINT13(.FALSE.,UA,UA,UAB,UAB,DUUP1(L,JATOM,is),JATOM)   ! duup1 = <g_{l+1}|2*dotf_l>    ==    <u_{l+1}|d/dr|dotu_l>
      CALL RINT13(.FALSE.,UBA,UBA,UB,UB,DUPU1(L,JATOM,is),JATOM)   ! dupu1 = <dotg_{l+1}|2*f_l>    == <dotu_{l+1}|d/dr|u_l>
      CALL RINT13(.FALSE.,UBA,UBA,UAB,UAB,DUPUP1(L,JATOM,is),JATOM)! dupup1= <dotg_{l+1}|2*dotf_l> == <dotu_{l+1}|d/dr|dotu_l>
   ENDDO             ! 203
   !ad...............DUU2(L)..DUUP2(L)..DUPU2(L)..DUPUP2(L)................
   DO L=1,LMAX             ! 303
      DO I=1,IMAX          ! 305
         UA(I) =U(I,L)               ! g_l
         UB(I) =DU(I,L+1)*TWO        ! 2*f_{l+1}
         UAB(I)=DUP(I,L+1)*TWO       ! 2*dotf_{l+1}
         UBA(I)=UP(I,L)              ! dotg_l
      ENDDO                ! 305     CONTINUE
      CALL RINT13(.FALSE.,UA,UA,UB,UB,DUU2(L,JATOM,is),JATOM)      ! duu2  = <g_l|2*f_{l+1}>       == <u_l|d/dr|u_{l+1}>
      CALL RINT13(.FALSE.,UA,UA,UAB,UAB,DUUP2(L,JATOM,is),JATOM)   ! duup2 = <g_l|2*dotf_{l+1}>    == <u_l|d/dr|dotu_{l+1}>
      CALL RINT13(.FALSE.,UBA,UBA,UB,UB,DUPU2(L,JATOM,is),JATOM)   ! dupu2 = <dotg_l|2*f_{l+1}>    == <dotu_l|d/dr|u_{l+1}>
      CALL RINT13(.FALSE.,UBA,UBA,UAB,UAB,DUPUP2(L,JATOM,is),JATOM)! dupup2= <dotg_l|2*dotf_{l+1}> == <dotu_l|d/dr|dotu_{l+1}>
   ENDDO                   ! 303
   !ad..............................for LO's...............................
   DO L=0,LOMAX  ! 503
      L1=L+1      ! l1 is actually the index to the l-th component of U(l), UP(l), UL(l)
      L11=L1+1    ! l11 is the index for U(l+1), UP(l+1), UL(l+1)
      IF (.not.loor(L)) THEN
         DUUL1(L1,JATOM,is) =  0
         DUPUL1(L1,JATOM,is) = 0
         DULU2(L1,JATOM,is)  = 0
         DULUP2(L1,JATOM,is) = 0
         RUUL(L1,JATOM,is)   = 0
         RUPUL(L1,JATOM,is)  = 0
      ELSE
         DO I=1,IMAX     ! 505
            SR=SQRT(R(I))
            UA(I) =U (I,L11)/SR     ! g_{l+1}/sqrt(r)
            UB(I) =UP(I,L11)/SR     ! dotg_{l+1}/sqrt(r)
            UAB(I)=UL(I,L1) /SR     ! go_l/sqrt(r)
         ENDDO           ! 505
         !ad............DUUL1..DUPUL1..DULU2..DULUP2..RUUL..RUPUL................
         CALL RINT13(.FALSE.,U(1,L11),U(1,L11),DUL(1,L1),DUL(1,L1),DUUL1(L1,JATOM,is),JATOM)   ! <g_{1+1}|fo_l>
         DUUL1(L1,JATOM,is)=TWO*DUUL1(L1,JATOM,is)                                             ! duul1 = 2*<g_{1+1}|fo_l>     == <u_{l+1}|d/dr|ulo_l>
         CALL RINT13(.FALSE.,UP(1,L11),UP(1,L11),DUL(1,L1),DUL(1,L1),DUPUL1(L1,JATOM,is),JATOM)! <dotg_{l+1}|fo_l>
         DUPUL1(L1,JATOM,is)=TWO*DUPUL1(L1,JATOM,is)                                           ! dupul1 = 2*<dotg_{l+1}|fo_l> == <dotu_{l+1}|d/dr|ulo_l>
         CALL RINT13(.FALSE.,UL(1,L1),UL(1,L1),DU(1,L11),DU(1,L11),DULU2(L1,JATOM,is),JATOM)   ! <go_l|f_{l+1}>
         DULU2(L1,JATOM,is)=TWO*DULU2(L1,JATOM,is)                                             ! dulu2 = 2*<go_l|f_{l+1}>     == <ulo_l|d/dr|u_{l+1}>
         CALL RINT13(.FALSE.,UL(1,L1),UL(1,L1),DUP(1,L11),DUP(1,L11),DULUP2(L1,JATOM,is),JATOM)! <go_l|dotf_{l+1}>
         DULUP2(L1,JATOM,is)=TWO*DULUP2(L1,JATOM,is)                                           ! dulup2 = 2*<go_l|dotf_{l+1}> == <ulo_l|d/dr|dotu_{l+1}>
         CALL RINT13(.FALSE.,UA,UA,UAB,UAB,RUUL(L1,JATOM,is),JATOM)                            ! ruul = <g_{l+1}|1/r|go_l>    == <u_{l+1}|1/r|ulo_l>
         CALL RINT13(.FALSE.,UB,UB,UAB,UAB,RUPUL(L1,JATOM,is),JATOM)                           ! rupul= <dotg_{l+1}|1/r|go_l> == <dotu_{l+1}|1/r|ulo_l>
      END IF
      IF (L+1.le.LOMAX) THEN
         IF (.not.loor(L+1)) THEN
            DULU1(L1,JATOM,is)  = 0
            DULUP1(L1,JATOM,is) = 0
            DUUL2(L1,JATOM,is)  = 0
            DUPUL2(L1,JATOM,is) = 0
            DULUL1(L1,JATOM,is) = 0
            DULUL2(L1,JATOM,is) = 0
            RULU(L1,JATOM,is)   = 0
            RULUP(L1,JATOM,is)  = 0
            RULUL(L1,JATOM,is)  = 0
         ELSE
            DO I=1,IMAX    ! 506
               SR=SQRT(R(I))
               UA(I)=U(I,L1)/SR    ! g_l/sqrt(r)
               UB(I)=UP(I,L1)/SR   ! dotg_l/sqrt(r)
               UAB(I)=UL(I,L11)/SR ! go_{l+1}/sqrt(r)
               UBA(I)=UL(I,L1)/SR  ! go_{l}/sqrt(r)
            ENDDO          ! 506
            !ad
            !ad..........DULU1..DULUP1..DUUL2..DUPUL2..DULUP2..DULUL1..DULUL2.......
            !ad
            CALL RINT13(.FALSE.,UL(1,L11),UL(1,L11),DU(1,L1),DU(1,L1),DULU1(L1,JATOM,is),JATOM)    ! <go_{l+1}|f_l>
            DULU1(L1,JATOM,is)=TWO*DULU1(L1,JATOM,is)                                              ! dulu1 = 2*<go_{l+1}|f_l>      == <ulo_{l+1}|d/dr|u_l>
            CALL RINT13(.FALSE.,UL(1,L11),UL(1,L11),DUP(1,L1),DUP(1,L1),DULUP1(L1,JATOM,is),JATOM) ! <go_{l+1}|dotf_l>
            DULUP1(L1,JATOM,is)=TWO*DULUP1(L1,JATOM,is)                                            ! dulup1 = 2*<go_{l+1}|dotf_l>  == <ulo_{l+1}|d/dr|dotu_l>
            CALL RINT13(.FALSE.,U(1,L1),U(1,L1),DUL(1,L11),DUL(1,L11),DUUL2(L1,JATOM,is),JATOM)    ! <g_l|fo_{l+1}>
            DUUL2(L1,JATOM,is)=TWO*DUUL2(L1,JATOM,is)                                              ! duul2 = 2*<g_l|fo_{l+1}>      == <u_l|d/dr|ulo_{l+1}>
            CALL RINT13(.FALSE.,UP(1,L1),UP(1,L1),DUL(1,L11),DUL(1,L11),DUPUL2(L1,JATOM,is),JATOM) ! <dotg|fo_{l+1}>
            DUPUL2(L1,JATOM,is)=TWO*DUPUL2(L1,JATOM,is)                                            ! dupul2 = 2*<dotg_l|fo_{l+1}>  == <dotu_l|d/dr|ulo_{l+1}>
            CALL RINT13(.FALSE.,UL(1,L11),UL(1,L11),DUL(1,L1),DUL(1,L1),DULUL1(L1,JATOM,is),JATOM) ! <go_{l+1}|fo_l>
            DULUL1(L1,JATOM,is)=TWO*DULUL1(L1,JATOM,is)                                            ! dulul1 = 2*<go_{l+1}|fo_l>    == <ulo_{l+1}|d/dr|ulo_l>
            CALL RINT13(.FALSE.,UL(1,L1),UL(1,L1),DUL(1,L11),DUL(1,L11),DULUL2(L1,JATOM,is),JATOM) ! <go_l|fo_{l+1}>
            DULUL2(L1,JATOM,is)=TWO*DULUL2(L1,JATOM,is)                                            ! dulul2 = 2*<go_l|fo_{l+1}>== <ulo_l|d/dr|ulo_{l+1}>
            !ad...........................RULU..RULUP..RULUL........................
            CALL RINT13(.FALSE.,UAB,UAB,UA,UA,RULU(L1,JATOM,is),JATOM)      ! rulu = <go_{l+1}|1/r|g_l>    = <ulo_{l+1}|1/r|u_l>
            CALL RINT13(.FALSE.,UAB,UAB,UB,UB,RULUP(L1,JATOM,is),JATOM)     ! rulup = <go_{l+1}|1/r|dotg_l>= <ulo_{l+1}|1/r|dotu_l>
            CALL RINT13(.FALSE.,UAB,UAB,UBA,UBA,RULUL(L1,JATOM,is),JATOM)   ! rulul = <go_{l+1}|1/r|go_l>  = <ulo_{l+1}|1/r|ulo_l>
         END IF
      END IF
   ENDDO    ! 503

   !
   !....... END CALCULATION OF RADIAL INTEGRALS ...............
   RETURN
 END SUBROUTINE radint


text:
\ Animated hedgehog logo on OLED screen, ~430 Bytes
\ From well known words: MS

: EGEL      ( y x -- )          \ Print hedgehog at pos y x
    empty                       \ Erase custom characters
    dm 59 0 >cc  dm 52 0 >cc    \ Build four custom characters: p
    dm 58 1 >cc  dm 51 1 >cc    \ q
    dm 53 2 >cc  dm 56 2 >cc  9 2 >cc  \ r
    dm 45 3 >cc  dm 11 3 >cc    \ s
    >r  dup r@ xy" 0000000000W:`8c58050"
    dup r@ 1+  xy" 5700000000o8^^n]WjnXgV0"
    dup r@ 2 + xy" kpigbhhhge;q]^27o6n4[ii0"
    dup r@ 3 + xy" 0036grd700c0ef6hk07335of0"
    dup r@ 4 + xy" 0000000366kbkVk00b3f`46s0"
        r> 5 + xy" 00000000000000000003bf0" ;

: SHOW      ( -- )
    display-setup  graphic  &page   \ Init. and show Hedgehog
    38 00 do  80 i 2* - 1 egel  10 ms loop \ Move hedgehog to center screen
    80 ms  thin  4 6 xy" Graphic chars"
    A00 ms  &page                   \ Wait then show message
     4 0 xy" Project Forth"         \ Startup message
    28 2 xy" Works"                 \ To line 2
    1A 4 xy" Graphics"              \ To line 4
    20 6 xy" by W.O." ;             \ To line 6

\ End ;;;

max_stars_repo_name: embeddingforth/embeddinForth

max_issues_repo_name: embeddingforth/embeddinForth

max_forks_repo_name: embeddingforth/embeddinForth


MODULE subdirModuleA
END MODULE


text:
      subroutine gendenspsihalo

      include 'commonblocks'

      real psi

      open(file='denspsihalo.dat',unit=50,status='replace')

c      fcut_halo = dfhalo(psic)
      fcut_halo = 0.

      do i=1,npsi
         psi = tableE(i)
         call gethalodens(psi,denspsihalo(i))
         if(i.eq.npsi) denspsihalo(i) = 0.
         write(50,*) tableE(i),denspsihalo(i)
      enddo

      close(50)
      return

      end

      subroutine gethalodens(psi,rho)

      include 'commonblocks'

      real psi

      m = 200
      vmin = -10.
      vmax = log(sqrt(2.*(psi-psic)))
      dlogv = (vmax - vmin)/float(m-1)
      rho = 0.
      do j=1,4
         vlog = vmin + dlogv*float(j-1)
         v = exp(vlog)
         e = psi - v*v/2.
         df_lowered = dfhalo(e)-fcut_halo
         if(df_lowered.gt.0.)
     +         rho = rho + 4.*pi*dlogv*v*v*v*coef(j)*df_lowered
      enddo
      do j=5,m-4
         vlog = vmin + dlogv*float(j-1)
         v = exp(vlog)
         e = psi - v*v/2.
         df_lowered = dfhalo(e)-fcut_halo
         if(df_lowered.gt.0.)
     +         rho = rho + 4.*pi*dlogv*v*v*v*df_lowered
      enddo
      do jj=4,2,-1
         j = m-jj+1
         vlog = vmin + dlogv*float(j-1)
         v = exp(vlog)
         e = psi - v*v/2.
         df_lowered = dfhalo(e)-fcut_halo
         if(df_lowered.gt.0.)
     +         rho = rho + 4.*pi*dlogv*v*v*v*coef(jj)*df_lowered
      enddo
      return
      end

      function coef(j)

      if(j.eq.1) coef=17./48.
      if(j.eq.2) coef=59./48.
      if(j.eq.3) coef=43./48.
      if(j.eq.4) coef=49./48.

      return
      end

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
	SUBROUTINE PC_CKPM ( index, nparms, prmlst, iret )
C************************************************************************
C* PC_CKPM								*
C*									*
C* This subroutine saves the interpolation and extrapolation		*
C* flags for the PC computations.  The character flags, QCHR,		*
C* are no longer set in this subroutine.  If the parameter type		*
C* table cannot be opened or if it is too large for the buffer,		*
C* an error message will be written but no error will be returned.	*
C*									*
C* PC_CKPM ( INDEX, NPARMS, PRMLST, IRET )				*
C*									*
C* Input parameters:							*
C*	INDEX		INTEGER		Table index			*
C*	NPARMS		INTEGER		Number of parameters		*
C*	PRMLST (NPARMS)	CHAR*4		Parameter list			*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C**									*
C* Log:									*
C* M. desJardins/GSFC	 8/84						*
C* M. desJardins/GSFC	11/87	Eliminated check for character type	*
C* M. desJardins/GSFC	 9/88	Cleaned up				*
C* M. Linda/GSC		10/97	Corrected the prologue right border	*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	INCLUDE		'GMBDTA.CMN'
	INCLUDE		'pccmn.cmn'
C*
	CHARACTER*(*)	prmlst (*)
C*
	LOGICAL		found
C------------------------------------------------------------------------
C*	Read in parameter conversion table if it has not already been
C*	read.
C
	iret = 0
	IF  ( .not. pptabl )  THEN
	    CALL TB_PRMT  ( MAXPM, jptabl, pparm, pchr, pint, pext,
     +			    pang, ier )
	    pptabl = .true.
	END IF
C
C*	Loop through each parameter.
C
	DO  i = 1, nparms
	    found = .false.
	    j = 1
C
C*	    Find this parameter in the table.
C
	    DO WHILE  ( ( .not. found ) .and. ( j .le. jptabl ) )
		IF  ( prmlst (i) .eq. pparm (j) )  THEN
		    qint (i, index) = pint (j)
		    qext (i, index) = pext (j)
		    qang (i, index) = pang (j)
		    found = .true.
		ENDIF
		j = j + 1
	    ENDDO
C
C*	    If parameter was not found, use default values.
C
	    IF  ( .not. found )  THEN
		qint (i, index) = .true.
		qext (i, index) = .false.
		qang (i, index) = .false.
	    END IF
	ENDDO
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
module periodic_table

  use datatypes

  implicit none

  integer, parameter  :: n_species=94

  real(double), parameter ::  atomic_mass(n_species) = (/ &
       1.01_double,    4.00_double,  6.94_double,  9.01_double, 10.81_double, 12.01_double, &
       14.01_double,  16.00_double, 19.00_double, 20.18_double, 22.99_double, 24.31_double, &
       26.98_double,  28.09_double, 30.97_double, 32.07_double, 35.45_double, 39.95_double, &
       39.10_double,  40.08_double, 44.96_double, 47.88_double, 50.94_double, 52.00_double, &
       54.94_double,  55.85_double, 58.93_double, 58.69_double, 63.55_double, 65.39_double, &
       69.72_double,  72.61_double, 74.92_double, 78.96_double, 79.90_double, 83.80_double, &
       85.47_double,  87.62_double, 88.91_double, 91.22_double, 92.91_double, 95.94_double, &
       98.91_double, 101.07_double,102.91_double,106.42_double,107.87_double,112.41_double, &
       114.82_double,118.71_double,121.75_double,127.60_double,126.90_double,131.29_double, &
       132.91_double,137.33_double,138.91_double,140.12_double,140.91_double,144.24_double, &
       146.92_double,150.36_double,151.97_double,157.25_double,158.93_double,162.50_double, &
       164.93_double,167.26_double,168.93_double,173.04_double,174.97_double,178.49_double, &
       180.95_double,183.85_double,186.21_double,190.20_double,192.22_double,195.08_double, &
       196.97_double,200.59_double,204.38_double,207.20_double,208.98_double,208.98_double, &
       209.99_double,222.02_double,223.02_double,226.03_double,227.03_double,232.04_double, &
       231.04_double,238.03_double,237.05_double,244.06_double/)

  character(len=2), dimension(103) :: pte = &
       (/  "H ", "He", "Li", "Be", "B ", "C ", "N ", "O ", "F ", "Ne", &
         "Na", "Mg", "Al", "Si", "P ", "S ", "Cl", "Ar", "K ", "Ca", &
         "Sc", "Ti", "V ", "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn", &
         "Ga", "Ge", "As", "Se", "Br", "Kr", "Rb", "Sr", "Y ", "Zr", &
         "Nb", "Mo", "Tc", "Ru", "Rh", "Pd", "Ag", "Cd", "In", "Sn", &
         "Sb", "Te", "I ", "Xe", "Cs", "Ba", "La", "Ce", "Pr", "Nd", &
         "Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb", &
         "Lu", "Hf", "Ta", "W ", "Re", "Os", "Ir", "Pt", "Au", "Hg", &
         "Tl", "Pb", "Bi", "Po", "At", "Rn", "Fr", "Ra", "Ac", "Th", &
         "Pa", "U ", "Np", "Pu", "Am", "Cm", "Bk", "Cf", "Es", "Fm", &
         "Md", "No", "Lr"/)

end module periodic_table

max_stars_repo_name: negf/CONQUEST-release

max_issues_repo_name: negf/CONQUEST-release

max_forks_repo_name: negf/CONQUEST-release


! ******************************************************
!   WARNING: need 32 Bit version of the HDF5 library !
! ******************************************************

!>    @brief creates a HDF5 file
!>    @param[in] f_name : hdf5 file name
!>    @details
!>    create a new hdf5-file, used only for restart output.\n
      SUBROUTINE create_hdf5(f_name)
#ifndef noHDF
        USE hdf5
        use mod_hdf5_vars, only: error
#endif
        IMPLICIT NONE

!      arrayname and filename
        character (len=*) :: f_name

#ifndef noHDF
!      File identifiers
        INTEGER (hid_t) file_id
#endif

#ifndef noHDF

!      Initialize FORTRAN interface.
!aw      Call h5open_f(error)

!      Create a new file, later only open it for read and writes
!org      call h5fcreate_f(f_name, H5F_ACC_TRUNC_F, file_id, error)
        CALL h5fcreate_f(f_name,h5f_acc_trunc_f,file_id,error, &
          h5p_default_f,h5p_default_f)
        CALL h5fclose_f(file_id,error)

!      Close FORTRAN interface.
!aw      Call h5Close_f(error)

#else
        WRITE(*,*) 'error: HDF5 support was not compiled in'
        STOP
#endif
        RETURN
      END


text:
        subroutine stp_interp(s,t,p,n,s0,t0,p0)
ccc
ccc
ccc
ccc     DESCRIPTION :   Interpolate salinity and in situ temperature
ccc                     on a cast by linearly interpolating salinity
ccc                     and potential temperature
ccc
ccc     PRECISION :     Double
ccc
ccc     INPUT :         s(n)            array of cast salinities
ccc                     t(n)            array of cast in situ temperatures
ccc                     p(n)            array of cast pressures
ccc                     n               length of cast
ccc                     p0              pressure for which salinity and
ccc                                     in situ temperature are required
ccc
ccc     OUTPUT :        s0              interpolated value of salinity
ccc                     t0              interpolated value of situ temperature
ccc
ccc     UNITS :         salinities      psu (IPSS-78)
ccc                     temperatures    degrees C (IPTS-68)
ccc                     pressures       db
ccc
ccc
ccc     AUTHOR :        David Jackett
ccc
ccc     CREATED :       June 1993
ccc
ccc     REVISION :      1.1             30/6/93
ccc
ccc
ccc
        implicit double precision (a-h,o-z)

        dimension s(n),t(n),p(n)

        external indx

        data pr0/0.0/

        call indx(p,n,p0,k)

        r = (p0-p(k))/(p(k+1)-p(k))

        s0 = s(k) + r*(s(k+1)-s(k))

        thk = theta(s(k),t(k),p(k),pr0)

        th0 = thk + r*(theta(s(k+1),t(k+1),p(k+1),pr0)-thk)

        t0 = theta(s0,th0,pr0,p0)

        return
        end


text:
!------------------------------------------------------------------------------!
!  The Community Multiscale Air Quality (CMAQ) system software is in           !
!  continuous development by various groups and is based on information        !
!  from these groups: Federal Government employees, contractors working        !
!  within a United States Government contract, and non-Federal sources         !
!  including research institutions.  These groups give the Government          !
!  permission to use, prepare derivative works of, and distribute copies       !
!  of their work in the CMAQ system to the public and to permit others         !
!  to do so.  The United States Environmental Protection Agency                !
!  therefore grants similar permission to use the CMAQ system software,        !
!  but users are requested to provide copies of derivative works or            !
!  products designed to operate in the CMAQ system to the United States        !
!  Government without restrictions as to use by others.  Software              !
!  that is used with the CMAQ system but distributed under the GNU             !
!  General Public License or the GNU Lesser General Public License is          !
!  subject to their copyright restrictions.                                    !
!------------------------------------------------------------------------------!

SUBROUTINE vstamp

!-------------------------------------------------------------------------------
! Name:     Version Stamp
! Purpose:  Writes version and release identification to log file.
! Revised:  20 Sep 2001  Original version.  (T. Otte)
!           18 Oct 2001  Updated for second beta release of MCIPv2.  (T. Otte)
!           20 Nov 2001  Updated for third beta release of MCIPv2.  (T. Otte)
!           19 Mar 2002  Updated for release of MCIP Version 2.0.  (T. Otte)
!           27 Mar 2003  Updated for release of MCIP Version 2.1.  (T. Otte)
!           11 Jun 2003  Updated for release of MCIP Version 2.2.  (T. Otte)
!           18 Aug 2004  Updated for release of MCIP Version 2.3.  (T. Otte)
!           19 Aug 2005  Updated for release of MCIP Version 3.0   (T. Otte)
!           27 Feb 2006  Updated for release of MCIP Version 3.1.  (T. Otte)
!           18 Aug 2006  Updated for release of MCIP Version 3.2.  (T. Otte)
!           31 Jul 2007  Moved version descriptors to MCIPPARM.  Changed
!                        "RELEASED" to "FROZEN" in the prints to be more
!                        technically correct.  (T. Otte)
!           07 Sep 2011  Updated disclaimer.  (T. Otte)
!-------------------------------------------------------------------------------

  USE mcipparm

  IMPLICIT NONE

  WRITE (*, "(/, 1x, 78('='), /)")
  WRITE (*, "(18x, a)") 'US EPA COMMUNITY MULTISCALE AIR QUALITY MODEL'
  WRITE (*, "(20x, a, /)") 'METEOROLOGY-CHEMISTRY INTERFACE PROCESSOR'
  WRITE (*, "(27x, a, 1x, a, ' FROZEN ', a)")  &
                                     TRIM(progname), TRIM(ver), vdate
  WRITE (*, "(/, 1x, 78('='), ///)")

END SUBROUTINE vstamp

max_stars_repo_name: Simeng-unique/CMAQ-changed

max_issues_repo_name: Simeng-unique/CMAQ-changed

max_forks_repo_name: Simeng-unique/CMAQ-changed