text:
      SUBROUTINE INFALL(IBH,IESC,NBH2,ISUB)
*
*
*       Disruption of star by BH.
*       -------------------------
*
      INCLUDE 'common6.h'
      PARAMETER  (NMX=10,NMX3=3*NMX,NMX4=4*NMX,NMXm=NMX*(NMX-1)/2)
      REAL*8  M,M1,MASS,MC,MMIJ,XCM(3),VCM(3),CG(6),XX(3),VV(3)
      COMMON/ARCHAIN/XCH(NMX3),VCH(NMX3),WTTL,M(NMX),
     &           XCDUM(NMX3),WCDUM(NMX3),MC(NMX),
     &           XI(NMX3),VI(NMX3),MASS,RINV(NMXm),RSUM,INAME(NMX),NN
      COMMON/ARCHAIN2/ MMIJ,CMX(3),CMV(3),ENERGY,EnerGR,CHTIME
      common/TIMECOMMON/Taika,timecomparison
      COMMON/CHAINC/  XC(3,NCMAX),UC(3,NCMAX),BODYC(NCMAX),ICH,
     &                LISTC(LMAX)
      COMMON/CHREG/  TIMEC,TMAX,RMAXC,CM(10),NAMEC(NMX),NSTEP1,KZ27,KZ30
      COMMON/CLUMP/   BODYS(NCMAX,5),T0S(5),TS(5),STEPS(5),RMAXS(5),
     &                NAMES(NCMAX,5),ISYS(5)
      COMMON/CCOLL2/  QK(NMX4),PK(NMX4),RIK(NMX,NMX),SIZE(NMX),VSTAR1,
     &                ECOLL1,RCOLL,QPERI,ISTAR(NMX),ICOLL,ISYNC,NDISS1
      COMMON/POSTN2/  SEMIGR,ECCGR,DEGR,ISPIN
      COMMON/INCOND/  X4(3,NMX),XDOT4(3,NMX)
      REAL*8  G0(3),RVEC(3),VVEC(3)
*
*
*       Copy chain variables to standard form.
      LK = 0
      DO 4 L = 1,NCH
          DO 1 K = 1,3
              LK = LK + 1
              X4(K,L) = XCH(LK)
              XDOT4(K,L) = VCH(LK)
    1     CONTINUE
    4 CONTINUE
*
*       Obtain original c.m. kinetic energy.
      ZK1 = 0.0
      DO 5 K = 1,3
          ZK1 = ZK1 + XDOT(K,ICH)**2
    5 CONTINUE
      ZK1 = 0.5*BODY(ICH)*ZK1
*
*       Define new local c.m. for two-body system IBH & IESC.
      LX = IBH
      LK = 3*(IESC - 1)
      LN = 3*(LX - 1)
*       Choose half mass for ghost or whole star to be swallowed.
      IF (KZ(43).GE.2.AND.ISTAR(IESC).LT.10) THEN
          M1 = 0.1*M(IESC)
      ELSE
          M1 = M(IESC)
      END IF
*       Note use of XCM & VCM even if adopting direct escape.
      DO 10 K = 1,3
          LK = LK + 1
          LN = LN + 1
          XCM(K) = (M1*XCH(LK) + M(LX)*XCH(LN))/(M1 + M(LX))
          VCM(K) = (M1*VCH(LK) + M(LX)*VCH(LN))/(M1 + M(LX))
*       NB! XCH is not at actual pericentre (determined by A*(1-ECC)).
   10 CONTINUE
*
*       Ensure largest stellar type and reduce membership.
      ISTAR(LX) = MAX(ISTAR(IBH),ISTAR(IESC))
      NCH = NCH - 1
      NN = NCH
*
*       Check optional treatment of tidal disruption by BH.
      IF (KZ(43).GE.2.AND.ISTAR(IESC).LT.10) THEN
          NDISR = NDISR + 1
          ZMB =  BODYC(IESC) + BODYC(LX)
          M(IESC) = 0.1*M(IESC)
*       Accrete 1/10 the mass of #IESC and copy new BH variables.
          M(LX) = M(LX) + M(IESC)
          BODYC(LX) = BODYC(LX) + M(IESC)
          ZMASS = ZMASS - 9.0*M(IESC)
*       Update system masses.
          BODY(ICH) = BODY(ICH) - 9.0*M(IESC)
          MASS = MASS - 9.0*M(IESC)
          DO 20 K = 1,3
              X4(K,LX) = XCM(K)
              XDOT4(K,LX) = VCM(K)
   20     CONTINUE
*
*       Obtain external energy for ghost (internal pot in ECH).
          POT1 = 0.0
          DO 40 J = IFIRST,NTOT
              IF (J.EQ.ICH) GO TO 40
              RIJ2 = 0.0
              DO 35 K = 1,3
                  RIJ2 = RIJ2 + (X(K,ICH) - X(K,J))**2
   35         CONTINUE
              POT1 = POT1 + BODY(J)/SQRT(RIJ2)
   40     CONTINUE
*
*       Subtract ghost mass contribution to the potential energy.
          ECOLL = ECOLL - 9.0*M(IESC)*POT1
          WRITE (24,45)  TIME+TOFF, NDISR, NAMEC(IESC), ISTAR(IESC),
     &                   ECCGR, 10.0*M(IESC)*SMU, M(LX)*SMU, SEMIGR
   45     FORMAT (' DISRUPT2    T NDISR NM K* E M1 M2 SEMI ',
     &                          F8.1,I5,I7,I4,F10.6,2F6.1,1P,E10.2)
          CALL FLUSH(24)
*       Note that modified M(LX) and another component will be new KS.
*
          VR2 = 0.0
          RR2 = 0.0
          DO 46 K = 1,3
              RR2 = RR2 + (X4(K,IESC) - X4(K,LX))**2
              VI2 = VI2 + (XDOT4(K,IESC) - XDOT4(K,LX))**2
              RVEC(K) = X4(K,IESC) - X4(K,LX)
              VVEC(K) = XDOT4(K,IESC) - XDOT4(K,LX)
   46     CONTINUE
          RR = SQRT(RR2)
          SEMI = 2.0/RR - VI2/ZMB
          SEMI = 1.0/SEMI
          VA2 = ZMB/SEMI
          XFAC = SEMI/RR
          VFAC = SQRT(RR/SEMI)
          DO 47 K = 1,3
              X4(K,IESC) = XFAC*BODYC(LX)*RVEC(K)/ZMB
              X4(K,LX) = -XFAC*BODYC(IESC)*RVEC(K)/ZMB
              XDOT4(K,IESC) = VFAC*BODYC(LX)*VVEC(K)/ZMB
              XDOT4(K,LX) = -VFAC*BODYC(IESC)*VVEC(K)/ZMB
              XCM(K) = X4(K,LX)
              VCM(K) = XDOT4(K,LX)
   47     CONTINUE
          RI2 = 0.0
          VI2 = 0.0
          DO 200 K = 1,3
              RI2 = RI2 + (X4(K,IESC) - X4(K,LX))**2
              VI2 = VI2 + (XDOT4(K,IESC) - XDOT4(K,LX))**2
  200     CONTINUE
          RIJ = SQRT(RI2)
          ANEW = 2.0/RIJ - VI2/ZMB
          ANEW = 1.0/ANEW
          WRITE (6,210)  RVEC, RIJ, ANEW, VA2, VI2
  210     FORMAT (' INFALL TRANSF    RVEC R A VA2 VI2 ',1P,7E10.2)
      ELSE
*       Implement complete swallowing of compact star.
          M(LX) = M(LX) + M(IESC)
          BODYC(LX) = BODYC(LX) + M(IESC)
          NCOLL = NCOLL + 1
          WRITE (6,50)  NAMEC(IESC), ISTAR(IESC), BODYC(IESC)*SMU,
     &                  M(LX)*SMU
   50     FORMAT (' SWALLOWED STAR    NAM K* M1 M2 ',I6,I4,2F7.2)
          NDISR = NDISR + 1
          WRITE (24,45)  TIME+TOFF, NDISR, NAMEC(IESC), ISTAR(IESC),
     &                   ECCGR, BODYC(IESC)*SMU, M(LX)*SMU, SEMIGR
      END IF
*
*       Check possible reduction of dominant body index.
      IF (LX.GT.IESC) LX = LX - 1
*
*       Identify global index of coalescence/escaper body.
      I = 0
      DO 60 J = IFIRST,NTOT
          IF (NAME(J).EQ.NAMEC(IESC)) THEN
              I = J
          END IF
   60 CONTINUE
*
*       Define M(IESC) as ghost (partial or complete swallowing).
      CALL GHOST(I)
*       Ensure escape condition (R**2 < 1D+10).
      X0(1,I) = 1.0D+04
      X(1,I) = 1.0D+04
*
*       Remove chain (and clump) mass & reference name of absorbed member.
      LI = 3*(IESC - 1)
      DO 70 L = IESC,NCH
          M(L) = M(L+1)
          BODYC(L) = BODYC(L+1)
          NAMEC(L) = NAMEC(L+1)
          SIZE(L) = SIZE(L+1)
          ISTAR(L) = ISTAR(L+1)
          BODYS(L,ISUB) = BODYS(L+1,ISUB)
          NAMES(L,ISUB) = NAMES(L+1,ISUB)
          DO 65 K = 1,3
              XCH(LI+K) = XCH(LI+K+3)
              VCH(LI+K) = VCH(LI+K+3)
   65     CONTINUE
          LI = LI + 3
   70 CONTINUE
*
*       Set XCM & VCM in XCH & VCH of the current IBH.
      LK = 3*(LX - 1)
      DO 55 K = 1,3
          LK = LK + 1
          XCH(LK) = XCM(K)
          VCH(LK) = VCM(K)
   55 CONTINUE
*
*       Copy new chain coordinates & velocities to standard variables.
      LK = 0
      DO 80 L = 1,NCH
          DO 75 K = 1,3
              LK = LK + 1
              X4(K,L) = XCH(LK)
              XDOT4(K,L) = VCH(LK)
   75     CONTINUE
   80 CONTINUE
*
*       Perform c.m. check (not perfect before XX & VV correction).
      DO 85 K = 1,6
          CG(K) = 0.0
   85 CONTINUE
*
      LK = 0
      DO 95 L = 1,NCH
          DO 90 K = 1,3
              LK = LK + 1
              CG(K) = CG(K) + BODYC(L)*XCH(LK)
              CG(K+3) = CG(K+3) + BODYC(L)*VCH(LK)
   90     CONTINUE
   95 CONTINUE
*
      DO 96 K = 1,6
          CG(K) = CG(K)/BODY(ICH)
   96 CONTINUE
*
*       Adopt c.m. condition and copy to X4 & XDOT4.
      LK = 0
      DO 98 L = 1,NCH
          DO 97 K = 1,3
              LK = LK + 1
              XCH(LK) = XCH(LK) - CG(K)
              VCH(LK) = VCH(LK) - CG(K+3)
              X4(K,L) = XCH(LK)
              XDOT4(K,L) = VCH(LK)
   97     CONTINUE
   98 CONTINUE
*
      IF (KZ(30).NE.0) THEN
          WRITE (6,99)  TIME+TOFF, (CG(K),K=1,6)
   99     FORMAT (' INFALL CHECK:   T CG ',F10.2,1P,6E9.1)
      END IF
*
*       Form new c.m. kinetic energy.
      ZK2 = 0.0
      DO 120 K = 1,3
          ZK2 = ZK2 + XDOT(K,ICH)**2
  120 CONTINUE
      ZK2 = 0.5*BODY(ICH)*ZK2
*
*       Add c.m. kinetic energy change for conservation.
      ECOLL = ECOLL + (ZK1 - ZK2)
*
      RETURN
*
      END

ext: f


! { dg-do compile }
Program p
  Use ieee_arithmetic
  Use iso_c_binding
  Print *, ieee_value(x=1.0, class=ieee_negative_inf)
  Print *, c_associated(c_ptr_1=c_null_ptr)
End Program


text:
      double precision function ubtzdsq(p,j1,j2,j3,j4,j6)
C     Matrix element squared averaged over initial colors and spins
c     u(-p1)+b(p2)->Z(p3,p4)+t(p5)+d(p6)
      implicit none
      include 'constants.f'
      include 'ewcouple.f'
      double precision p(mxpart,4),ubtzd
      double complex amp(2,2)
      integer j1,j2,j3,j4,j6

      call ubtzdamp(p,j1,j2,j3,j4,j6,amp)
      ubtzd=cdabs(amp(1,1))**2+cdabs(amp(1,2))**2
     &     +cdabs(amp(2,1))**2+cdabs(amp(2,2))**2

      ubtzdsq=4d0*xn**2*aveqq*gwsq**2*esq**2*ubtzd

      return
      end

      subroutine ubtzdamp(p,j1,j2,j3,j4,j6,amp)
c--- Author: R. Rontsch (15 Feb 2013)
c--- Amplitudes taken from Campbell, Ellis, Rontsch arXiv:1302.3856 (appendix A)
c--- note: overall factor of 1/(s34-mz^2) extracted in a slightly
c--- different way for numerical stability
      implicit none
      include 'constants.f'
      include 'masses.f'
      include 'ewcouple.f'
      include 'zcouple.f'
      include 'zprods_decl.f'
      double precision p(mxpart,4),q(mxpart,4),dot
      double complex facuLl,facuRl,facdLl,facdRl
      double precision s25,s16,s34,s234,s235,s345,s346,s134,mtsq
      double complex amplower(2,2),ampmiddle(2,2),ampextra(2,2),
     &     ampupper(2,2),amp(2,2)
      double complex prW,zab2,cprop,iprZ
      integer eta,p1,p2,p3,p4,pl,pa,j,j1,j2,j3,j4,j6,jb,nu
      prW(s16)=cone/dcmplx(s16-wmass**2,zip)
      zab2(p1,p2,p3,p4)=za(p1,p2)*zb(p2,p4)+za(p1,p3)*zb(p3,p4)

      mtsq=mt**2

      jb=j2   ! momentum to make top massless
      eta=7
      do nu=1,4
      q(1,nu)=p(j1,nu)
      q(2,nu)=p(j2,nu)
      q(3,nu)=p(j3,nu)
      q(4,nu)=p(j4,nu)
      q(5,nu)=p(5,nu)-mtsq*p(jb,nu)/(2d0*dot(p,5,jb))
      q(6,nu)=p(j6,nu)
      q(eta,nu)=mtsq*p(jb,nu)/(2d0*dot(p,5,jb))
      enddo
      s16=2d0*dot(p,j1,j6)
      s25=2d0*dot(p,j2,5)+mt**2
      s345=2d0*dot(p,j1,j2)+2d0*dot(p,j2,j6)+2d0*dot(p,j1,j6)
      s134=2d0*dot(p,j1,j3)+2d0*dot(p,j1,j4)+2d0*dot(p,j3,j4)
      s234=2d0*dot(p,j2,j3)+2d0*dot(p,j2,j4)+2d0*dot(p,j3,j4)
      s346=2d0*dot(p,j6,j3)+2d0*dot(p,j6,j4)+2d0*dot(p,j3,j4)
      s34=2d0*dot(p,j3,j4)
      s235=2d0*dot(p,j1,j6)+2d0*dot(p,j1,j4)+2d0*dot(p,j6,j4)

c--- top production (j3=3,j4=4) as-is
      if (j3.eq.3) then
        call spinoru(7,q,za,zb)
      else
c--- antitop production (j3=4,j4=3) requires c.c.
        call spinoru(7,q,zb,za)
      endif

c--- Implementation of Baur-Zeppenfeld treatment of Z width
      cprop=dcmplx(1d0/dsqrt((s34-zmass**2)**2+(zmass*zwidth)**2))
      iprZ=dcmplx(s34-zmass**2)

      do j=1,2
      if (j .eq. 1) then
         pl=3
         pa=4
         facuLl=dcmplx(Qu*q1)*iprZ+dcmplx(L(2)*le)*s34
         facuRl=dcmplx(Qu*q1)*iprZ+dcmplx(R(2)*le)*s34
         facdLl=dcmplx(Qd*q1)*iprZ+dcmplx(L(1)*le)*s34
         facdRl=dcmplx(Qd*q1)*iprZ+dcmplx(R(1)*le)*s34
      else
         pl=4
         pa=3
         facuLl=dcmplx(Qu*q1)*iprZ+dcmplx(L(2)*re)*s34
         facuRl=dcmplx(Qu*q1)*iprZ+dcmplx(R(2)*re)*s34
         facdLl=dcmplx(Qd*q1)*iprZ+dcmplx(L(1)*re)*s34
         facdRl=dcmplx(Qd*q1)*iprZ+dcmplx(R(1)*re)*s34
      endif

      ampupper(j,1)=prW(s25)/s34*( facuLl/s134*za(5,6)*zb(1,pa)*
     &     zab2(pl,1,pa,2) - facdLl/s346*za(pl,6)*zb(1,2)*
     &     (zab2(5,pl,6,pa)))

      ampupper(j,2)=prW(s25)*mt/s34*( facuLl/s134*za(6,eta)*zb(1,pa)/
     &     za(5,eta)*zab2(pl,1,pa,2) + facdLl/s346*za(pl,6)*zb(1,2)
     &     /za(5,eta)*zab2(eta,pl,6,pa))

      amplower(j,1)=prW(s16)/s34*( -facuRl*mtsq/(s345-mtsq)*za(pl,6)
     &     *zb(1,2)*zb(pa,eta)/zb(5,eta) + facdLl/s234*zab2(pl,2,pa,1)*
     &     za(6,5)*zb(2,pa) - facuLl/(s345-mtsq)*zab2(6,1,2,pa)
     &     *za(pl,5)*zb(1,2) )

      amplower(j,2)=prW(s16)*mt/s34*(
     &     facuRl/(s345-mtsq)*za(pl,6)*zb(1,2)* zb(pa,5)
     &     - facdLl/s234*zab2(pl,2,pa,1)*za(6,eta)*zb(2,pa)/
     &     za(5,eta)
     &     +facuLl/(s345-mtsq)*zab2(6,1,2,pa)*za(pl,eta)*
     &     zb(1,2)/za(5,eta) )

      if (j .eq. 1) then
         ampextra(j,1)=prW(s25)*prW(s16)/(2d0*xw*s235)*za(pl,5)*zb(1,pa)
     &        *zab2(6,1,pa,2)*iprZ
         ampextra(j,2)=-prW(s25)*prW(s16)*mt/(2d0*xw*s235)*za(pl,eta)*
     &        zb(1,pa)*zab2(6,1,pa,2)/za(5,eta)*iprZ
      else
         ampextra(j,1)=czip
         ampextra(j,2)=czip
      endif
      ampmiddle(j,1)=prW(s25)*prW(s16)/s34*(-(facuLl-facdLl)*
     &     ( zab2(pl,1,6,pa)*za(6,5)*zb(1,2) + zab2(5,1,6,2)*
     &     za(pl,6)*zb(1,pa) + zab2(6,pl,pa,1)*za(pl,5)*zb(2,pa) )
     &     + mtsq/2d0*za(pl,6)*zb(1,pa)*zb(2,eta)/zb(5,eta)*
     &     (facuLl-facdLl-facuRl+facdRl) )

      ampmiddle(j,2)=prW(s25)*prW(s16)*mt/s34*( (facuLl-facdLl)*(
     &     (za(6,eta)*zb(1,2)*zab2(pl,1,6,pa) + zab2(eta,1,6,2)
     &     *za(pl,6)*zb(1,pa) + zab2(6,pl,pa,1)*za(pl,eta)*zb(2,pa))/
     &     za(5,eta)) + za(pl,6)/2d0*zb(1,pa)*zb(2,5)*
     &     (facuLl-facdLl+facuRl-facdRl))

      enddo
      amp=(amplower+ampupper+ampmiddle+ampextra)*cprop

      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
      PROGRAM xanneal
C     driver for routine anneal
      INTEGER NCITY
      PARAMETER (NCITY=10)
      INTEGER i,idum,ii,iorder(NCITY)
      REAL ran3,x(NCITY),y(NCITY)
C     create points of sale
      idum=-111
      do 11 i=1,NCITY
        x(i)=ran3(idum)
        y(i)=ran3(idum)
        iorder(i)=i
11    continue
      call anneal(x,y,iorder,NCITY)
      write(*,*) '*** System Frozen ***'
      write(*,*) 'Final path:'
      write(*,'(1x,t3,a,t13,a,t23,a)') 'city','x','y'
      do 12 i=1,NCITY
        ii=iorder(i)
        write(*,'(1x,i4,2f10.4)') ii,x(ii),y(ii)
12    continue
      END

ext: for

max_stars_repo_path: NR-Functions/Numerical Recipes- Example & Functions/Examples/xanneal.for


text:
	subroutine Vegas_Pomp(integrand,result,absacc,relacc,ndim,ncall,
     >                        maxiter,init)
! Uses kahan summation to guarantee identical results independent of number of threads
	implicit none
	include 'nf.f'
	include 'lc.f'
	include 'mxdim.f'
	include 'mxpart.f'
        include 'phasemin.f'
        include 'cutoff.f'
        include 'jetcuts.f'
        include 'breit.f'
        include 'zerowidth.f'
        include 'srdiags.f'
        include 'interference.f'
        include 'qcdcouple.f'
        include 'ewcouple.f'
        include 'masses.f'
        include 'ipsgen.f'
	include 'scale.f'
	include 'stopscales.f'
        include 'gridinfo.f'
	include 'qlfirst.f'
	include 'ptilde.f'
	include 'bitflags.f'
	include 'flags.f'
	include 'lastphot.f'
	include 'b0.f'
	include 'nodecay.f'
	include 'swapxz.f'
	include 'heavyflav.f'
	include 'nflav.f'
	include 'notag.f'
	include 'reset.f'
	double precision result, relacc, absacc,integrand
	integer ndim, ncall, maxiter, neval,init,i
	external integrand
	double precision fun, sfun, sfun2
	double precision sint,sint2,sweight
	double precision fun2, weight
	double precision r, dr, xo, xn, err,ran2
	double precision Ingrid,Incall
	integer iter, calls, dim, grid, g, c, cmax
	integer ngrid,j,jj
	parameter (ngrid = 10)
	double precision xi(ngrid, MXDIM), d(ngrid, MXDIM)
	double precision x(mxdim), imp(ngrid), tmp(ngrid - 1)
	integer pos(NDIM)
	double precision t,dt,cfun,cfun2,cd(ngrid,MXDIM)
        double precision p1ext(4),p2ext(4)

	character*255 runname
	integer nlength
	logical bin,dryrun
	common/bin/bin
	common/runname/runname
	common/nlength/nlength
	common/dryrun/dryrun
        common/pext/p1ext,p2ext
        save xi
!$omp threadprivate(/pext/)
	bin=.false.
	if(init.le.0)then
	   do i=1,ndim
	      xi(1,i)=1d0
           enddo
*       define the initial distribution of intervals
	   Ingrid=1d0/dble(ngrid)
	   do dim = 1, ndim
	      do grid = 1, ngrid
	         r = dble(grid)*Ingrid
	         xi(grid, dim) = r
	      enddo
	   enddo
	endif
	if (init.le.1)then
	   neval = 0
	   sint=0d0
	   sweight=0d0
	   sint2=0d0
	   iter=0
	endif

c--- read-in grid if necessary
        if (readin) then
           if (dryrun) then
	      open(unit=11,file=ingridfile,status='unknown')
           else
	      open(unit=11,file=runname(1:nlength)//'_'
     .             //ingridfile,status='unknown')
           endif

	   write(6,*)'*********************************************'
	   write(6,*)'* Reading in vegas grid from ',
     .     	runname(1:nlength)//'_'//ingridfile,' *'
           write(6,*)'*********************************************'
           call flush(6)
           do j=1,ndim
             read(11,203) jj,(xi(i,j),i=1,ngrid)
           enddo
           close(11)
           readin=.false.
	endif

*       iterations loop
 1	continue
	iter = iter + 1
*       initialize iteration variables
	sfun = 0d0
	sfun2 = 0d0
	d(:,:)=0d0
	cfun = 0d0
	cfun2 = 0d0
	cd(:,:)=0d0
	Incall=1d0/dble(ncall)
!$omp  parallel do
!$omp& schedule(dynamic)
!$omp& default(private)
!$omp& shared(incall,xi,ncall,ndim,sfun,sfun2,d,cfun,cfun2,cd)
!$omp& copyin(/xmin/,/taumin/,/cutoff/,/jetcuts/,/breit/,/zerowidth/)
!$omp& copyin(/srdiags/,/vsymfact/,/qcdcouple/,/ewcouple/,/masses/)
!$omp& copyin(/interference/,/ipsgen/,/scale/,/stopscales/)
!$omp& copyin(/pext/,/ColC/,/qlfirst/,/ptildes/)
!$omp& copyin(/bitflags/,/flags/,/lastphot/)
!$omp& copyin(/QCDb0/,/nodecay/,/swapxz/,/heavyflav/)
!$omp& copyin(/notag/,/nflav/,/reset/)
	do calls = 1, ncall
	   weight = Incall
!$omp critical
           do i=1,ndim
              x(i)=ran2()
           enddo
!	   call GetRandom(x)
!$omp end critical
	   do dim = 1, ndim+2
	      r = x(dim)*ngrid + 1
	      grid = int(r)
	      xo = 0
	      if( grid .gt. 1 ) xo = xi(grid - 1, dim)
	      xn = xi(grid, dim) - xo
	      x(dim) = xo + (r - grid)*xn
	      pos(dim) = grid
	      weight = weight*xn*ngrid
	   enddo
*       compute the integrand
	   fun=integrand(x,weight)
	   fun = fun*weight
	   fun2 = fun**2
!$omp critical
           t=sfun+fun
	   if (abs(sfun).ge.abs(fun)) then
	      cfun=cfun+((sfun-t)+fun)
	   else
	      cfun=cfun+((fun-t)+sfun)
	   endif
	   sfun=t
           t=sfun2+fun2
	   if (abs(sfun2).ge.abs(fun2)) then
	      cfun2=cfun2+((sfun2-t)+fun2)
	   else
	      cfun2=cfun2+((fun2-t)+sfun2)
	   endif
	   sfun2=t
	   do dim = 1, ndim
	      i=pos(dim)
	      dt=d(i,dim)
	      t=dt+fun2
	      if (abs(dt).ge.abs(fun2)) then
		 cd(i,dim)=cd(i,dim)+((dt-t)+fun2)
	      else
		 cd(i,dim)=cd(i,dim)+((fun2-t)+dt)
	      endif
	      d(i,dim)=t
	   enddo
!$omp end critical
 666	   continue
	enddo
!$omp end parallel do
	sfun=sfun+cfun
	sfun2=sfun2+cfun2
	d(:,:)=d(:,:)+cd(:,:)
	neval = neval + ncall
*       compute the integral and error values
	err = 0
	fun2 = sfun**2
	sint2 = sint2 + fun2
	r = sfun2*ncall - fun2
	if( r .ne. 0 ) then
	   weight = fun2/abs(r)*(ncall - 1)
	   sweight = sweight + weight
	   sint = sint + sfun*weight
	endif
	if( sweight .eq. 0 ) then
	   result = 0
	else
	   r = sint/sweight
	   result = r
*       if the integrand is very close to zero, it is pointless (and costly)
*       to insist on a certain relative accuracy
	   if( abs(r) .gt. absacc ) then
              r = sint2/(sint*r)
	      if( r .gt. err ) then
		 err = r
	      endif
	   endif
	endif
	err = sqrt(err/iter)

        print *, "iteration ", iter, ":",result,"+/-",result*err
        if (abs(result) .lt. 1d-9) then
          print *, "integral is zero, exiting"
          return
        endif
!	if( err .le. relacc ) then
!	   print *, "iteration ", iter, ":",result,"+/-",result*err
!	   call Outhist(iter)
!	   return
!	endif

*       redefine the grid (importance sampling)
*       - smooth the f^2 value stored for each interval
	do dim = 1, ndim
	   xo = d(1, dim)
	   xn = d(2, dim)
	   d(1, dim) = .5D0*(xo + xn)
	   x(dim) = d(1, dim)
	   do grid = 2, ngrid - 1
	      r = xo + xn
	      xo = xn
	      xn = d(grid + 1, dim)
	      d(grid, dim) = (r + xn)/3D0
	      x(dim) = x(dim) + d(grid, dim)
	   enddo
	   d(ngrid, dim) = .5D0*(xo + xn)
	   x(dim) = x(dim) + d(ngrid, dim)
	enddo

*       - compute the importance function of each interval
	do dim = 1, ndim
	   r = 0
	   do grid = 1, ngrid
	      imp(grid) = 0
	      if( d(grid, dim) .gt. 0 ) then
		 xo = x(dim)/d(grid, dim)
		 imp(grid) = ((xo - 1)/xo/log(xo))**1.5D0
	      endif
	      r = r + imp(grid)
	   enddo
	   r = r/ngrid

*       - redefine the size of each interval
	   dr = 0
	   xn = 0
	   g = 0
	   do grid = 1, ngrid - 1
	      do while( dr .lt. r )
		 g = g + 1
		 dr = dr + imp(g)
		 xo = xn
		 xn = xi(g, dim)
	      enddo
	      dr = dr - r
	      tmp(grid) = xn - (xn - xo)*dr/imp(g)
	   enddo
	   do grid = 1, ngrid - 1
	      xi(grid, dim) = tmp(grid)
	   enddo
	   xi(ngrid, dim) = 1
	enddo
c--- added to write out intermediate results
!	if ((bin) .and. (iter .lt. maxiter)) then
!	   write(6,*)'Writing out intermediate results for iteration',iter
!	   call histofin(tgral,sd,iter,maxiter)
!	endif
 28	continue
c--- write-out grid if necessary
	if (writeout) then
           open(unit=11,file=runname(1:nlength)//'_'
     .           //outgridfile,status='unknown')
	   write(6,*)'***********************************************'
	   write(6,*)'* Writing out vegas grid to ',
     .               runname(1:nlength)//'_'//outgridfile,'  *'
           write(6,*)'***********************************************'
           call flush(6)
           do j=1,ndim
             write(11,203) j,(xi(i,j),i=1,ngrid)
           enddo
           close(11)
         endif

	if( iter .ge. maxiter ) then
!	   print *,"iterations reached set maximum of ",maxiter
!	   print *, "iteration ", iter, ":",result,"+/-",result*err
	   return
	endif

	goto 1
203     FORMAT(/(5z16))

	end

************************************************************************
*       * IniRandom sets up the random-number generator to produce at most
*       * max dims-dimensional quasi-random vectors

	subroutine IniRandom(dims)
	implicit none
	integer max, dims

	integer ndim,a,b,ni(55),n(55)
	data ni /
     1  980629335, 889272121, 422278310,1042669295, 531256381,
     2  335028099,  47160432, 788808135, 660624592, 793263632,
     3  998900570, 470796980, 327436767, 287473989, 119515078,
     4  575143087, 922274831,  21914605, 923291707, 753782759,
     5  254480986, 816423843, 931542684, 993691006, 343157264,
     6  272972469, 733687879, 468941742, 444207473, 896089285,
     7  629371118, 892845902, 163581912, 861580190,  85601059,
     8  899226806, 438711780, 921057966, 794646776, 417139730,
     9  343610085, 737162282,1024718389,  65196680, 954338580,
     1  642649958, 240238978, 722544540, 281483031,1024570269,
     2  602730138, 915220349, 651571385, 405259519, 145115737 /
	common /rngdata/ ndim,a,b,n

	a=55
	b=31
	n(:)=ni(:)
	ndim = dims
	end

************************************************************************
*       * GetRandom is a subtractive Mitchell-Moore random-number generator.
*       * The algorithm is n(i) = (n(i - 24) - n(i - 55)) mod m, implemented
*       * as a circular array with n(i + 55) = n(i) and m = 2^30 in this
*       * version.  The array n has been initialized by setting n(i) = i and
*       * running the algorithm 100,000 times.  Code by Ronald Kleiss.

	subroutine GetRandom(array)
	implicit none
	double precision array(*)
	integer ndim
	integer dim, a, b, j, m, n(55)
	common /rngdata/ ndim,a,b,n
	parameter (m = 2**30)

	do dim = 1, ndim
	   a = mod(a, 55) + 1
	   b = mod(b, 55) + 1
	   j = n(b) - n(a)
	   if( j .lt. 0 ) j = j + m
	   n(a) = j
	   array(dim) = dble(j)/m
	enddo
	end


text:
program test
      character(len=80) :: string
      character(len=20) :: string2
      string = "hello world"
      string2 = string
      print *, string2, string, string2
end program test

max_stars_repo_name: clementval/fc

max_issues_repo_name: clementval/fc

max_forks_repo_name: clementval/fc


text:
!
!   VnVFPRetJ.f90
!
!
!
!
!
!
!   Created by Kurt Mitman on 14/10/14.
!

DOUBLE PRECISION FUNCTION FnVFPRetJ(bsav)

USE Params
USE Globals
USE funcs
USE Procedures

IMPLICIT NONE

DOUBLE PRECISION, INTENT(IN) :: bsav
DOUBLE PRECISION             :: cons,houseval
DOUBLE PRECISION             :: vals(2)
INTEGER                      :: inds(2)

cons=gliq-FnQb(bsav,gPh*ghouse-gmort)*bsav

if(cons .LT. cmin) then
    FnVFPRetJ = 1.d10+10.0d0*(cmin-cons)

else

    FnVFPRetJ = -1.0d0*(uown(cons,ghouse,gij)+ubequest(bsav+gPh*ghouse-gmort*(1.0d0+grm)-FnAdj(gPh*ghouse)))

endif

END FUNCTION FnVFPRetJ

max_stars_repo_name: lnsongxf/housing-boom-bust

max_issues_repo_name: floswald/housing-boom-bust

max_forks_repo_name: kurtmitman/housing-boom-bust


text:
C Copyright (c) 2007-2017 National Technology & Engineering Solutions
C of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions are
C met:
C
C     * Redistributions of source code must retain the above copyright
C       notice, this list of conditions and the following disclaimer.
C
C     * Redistributions in binary form must reproduce the above
C       copyright notice, this list of conditions and the following
C       disclaimer in the documentation and/or other materials provided
C       with the distribution.
C
C     * Neither the name of NTESS nor the names of its
C       contributors may be used to endorse or promote products derived
C       from this software without specific prior written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
C

C=======================================================================
      SUBROUTINE DBIVTT (NDB, ISEVOK, ITMP, NELBLK, NVAREL)
C=======================================================================
C
C   --*** DBIVTT *** Read element variable truth table
C   --   Modified for ExodusII format 8/26/95
C   --
C   --Parameters:
C   --   NDB    - IN  - the database number
C   --   NELBLK - IN  - the number of element blocks
C   --   NVAREL - IN  - the number of element variables; <0 if end-of-file
C   --   ISEVOK - OUT - the dynamic memory array of the element block variable
C   --                  truth table;variable i,block j exists iff ISEVOK(j,i)

      INTEGER NDB
      INTEGER NELBLK, NVAREL
      LOGICAL ISEVOK(NELBLK,*)
      INTEGER ITMP(NVAREL, NELBLK)

C     Read the element block variable truth table
C       isevok - num_elem_var cycles faster
      if (nvarel .gt. 0) then
        CALL EXGVTT(NDB, NELBLK, NVAREL, ITMP, IERR)
        DO 110 I = 1, NVAREL
          DO 100 IELB = 1, NELBLK
            ISEVOK(IELB,I) = (ITMP(I,IELB) .NE. 0)
 100      CONTINUE
 110    CONTINUE
      end if
      RETURN
      END

max_stars_repo_name: brian-kelley/seacas

max_issues_repo_name: brian-kelley/seacas

max_forks_repo_name: brian-kelley/seacas


text:
      program main

c*********************************************************************72
c
cc MAIN is the main program for BROWNIAN_MOTION_SIMULATION_PRB.
c
c  Discussion:
c
c    BROWNIAN_MOTION_SIMULATION_PRB tests the BROWNIAN_MOTION_SIMULATION library.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license.
c
c  Modified:
c
c    30 September 2012
c
c  Author:
c
c    John Burkardt
c
      implicit none

      integer k
      parameter ( k = 40 )
      integer m_max
      parameter ( m_max = 3 )
      integer n
      parameter ( n = 1001 )

      double precision d
      double precision dsq(k,n)
      character * ( 80 ) header
      integer m
      integer seed
      double precision t
      double precision x(m_max,n)

      call timestamp ( )
      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 'BROWNIAN_MOTION_SIMULATION_PRB'
      write ( *, '(a)' ) '  FORTRAN77 version'
      write ( *, '(a)' )
     &  '  Test the BROWNIAN_MOTION_SIMULATION library.'
c
c  Compute the path of a particle undergoing Brownian motion.
c
      do m = 1, 2

        d = 10.0D+00
        t = 1.0D+00
        seed = 123456789
        call brownian_motion_simulation ( m, n, d, t, seed, x )
        if ( m .eq. 1 ) then
          header = 'motion_1d'
        else if ( m .eq. 2 ) then
          header = 'motion_2d'
        end if
        call brownian_motion_display ( m, n, x, header )

      end do
c
c  Estimate the average displacement of the particle from the origin
c  as a function of time.
c
      do m = 1, 3

        d = 10.0D+00
        t = 1.0D+00
        seed = 123456789

        call brownian_displacement_simulation ( k, n, m, d, t, seed,
     &    dsq )
        if ( m .eq. 1 ) then
          header = 'displacement_1d'
        else if ( m .eq. 2 ) then
          header = 'displacement_2d'
        else if ( m .eq. 3 ) then
          header = 'displacement_3d'
        end if
        call brownian_displacement_display ( k, n, d, t, dsq, header )

      end do
c
c  Terminate.
c
      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 'BROWNIAN_MOTION_SIMULATION_PRB'
      write ( *, '(a)' ) '  Normal end of execution.'
      write ( *, '(a)' ) ''
      call timestamp ( );

      return
      end


text:
! ##################################################################################################################################
! Begin MIT license text.
! _______________________________________________________________________________________________________

! Copyright 2019 Dr William R Case, Jr (dbcase29@gmail.com)

! Permission is hereby granted, free of charge, to any person obtaining a copy of this software and
! associated documentation files (the "Software"), to deal in the Software without restriction, including
! without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to
! the following conditions:

! The above copyright notice and this permission notice shall be included in all copies or substantial
! portions of the Software and documentation.

! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
! OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
! THE SOFTWARE.
! _______________________________________________________________________________________________________

! End MIT license text.

      SUBROUTINE TDOF_PROC ( TDOF_MSG )

! TDOF table generation. TDOF is a table that has the DOF number for every DOF. The table has NDOFG rows and MTDOF columns where:

!   NDOFG = the total number of degrees of freedom in the model (6 times the number of grids plus the number of SPOINT's)
!   MTDOF = the number of different displ sets (currently 14 plus 2 USET user's sets) plus 4 cols for the internal grid/component
!           and external grid/component

! The TDOF table is sorted in grid numerical order. It's sister table, TDOFI, is TDOF sorted in internal G-set order.
! Cols 5 thru MTDOF of these tables give the degree of freedom (DOF) number for each of the 16 different displ sets/user's sets.
! An example of a TDOF table written out to an F06 file is shown below:

!                                              DEGREE OF FREEDOM TABLE SORTED ON GRID NUMBER (TDOF)

! EXTERNAL  INTERNAL                                     DOF NUMBER FOR DISPLACEMENT SET:
! GRD-COMP  GRD-COMP ---------------------------------------------------------------------------------------------------------------
! NUMBER    NUMBER         G      M      N     SA     SB     SG     SZ     SE      S      F      O      A     R      L     U1     U2

!    1022-1       1-1      1      0      1      0      0      0      0      0      0      1      0      1     0      1      0      0
!        -2        -2      2      0      2      0      0      0      0      0      0      2      0      2     0      2      0      1
!        -3        -3      3      0      3      0      0      0      0      0      0      3      0      3     0      3      0      0
!        -4        -4      4      0      4      0      0      0      0      0      0      4      0      4     0      4      0      2
!        -5        -5      5      0      5      0      0      0      0      0      0      5      0      5     0      5      0      0
!        -6        -6      6      0      6      0      0      0      0      0      0      6      0      6     0      6      0      3

!    1043-1       2-1      7      0      7      0      1      0      1      0      1      0      0      0     0      0      1      0
!        -2        -2      8      0      8      0      2      0      2      0      2      0      0      0     0      0      0      0
!        -3        -3      9      0      9      0      0      0      0      0      0      7      0      7     0      7      2      0
!        -4        -4     10      0     10      0      0      0      0      0      0      8      0      8     0      8      0      0
!        -5        -5     11      0     11      0      0      0      0      0      0      9      0      9     0      9      3      0
!        -6        -6     12      0     12      0      0      0      0      0      0     10      0     10     0     10      0      0

!                     ------ ------ ------ ------ ------ ------ ------ ------ ------ ------ ------ ------ ----- ------ ------ ------
! TOTAL NUMBER OF DOF:    12      0     12      0      2      6      2      0      2     12      0     10     0     10      3      3

      USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG, DOUBLE
      USE IOUNT1, ONLY                :  WRT_ERR, WRT_LOG, ERR, F04, F06, SC1
      USE SCONTR, ONLY                :  BLNK_SUB_NAM, FATAL_ERR, LDOFG, MTDOF, NDOFA, NDOFF, NDOFG, NDOFL, NDOFM, NDOFN, NDOFO,   &
                                         NDOFR, NDOFS, NDOFSA, NDOFSB, NDOFSE, NDOFSG, NDOFSZ, NGRID, NUM_USET_U1, NUM_USET_U2,    &
                                         SOL_NAME, WARN_ERR
      USE PARAMS, ONLY                :  EIGESTL, PRTDOF
      USE TIMDAT, ONLY                :  TSEC
      USE SUBR_BEGEND_LEVELS, ONLY    :  DOF_PROC_BEGEND
      USE DOF_TABLES, ONLY            :  TSET, TDOF, TDOFI, TDOF_ROW_START, USET
      USE DEBUG_PARAMETERS, ONLY      :  DEBUG
      USE MODEL_STUF, ONLY            :  EIG_N2, GRID, GRID_ID, GRID_SEQ, INV_GRID_SEQ

      USE TDOF_PROC_USE_IFs

      IMPLICIT NONE

      CHARACTER, PARAMETER            :: CR13 = CHAR(13)   ! This causes a carriage return simulating the "+" action in a FORMAT
      CHARACTER(LEN=LEN(BLNK_SUB_NAM)):: SUBR_NAME = 'TDOF_PROC'
      CHARACTER(LEN=*), INTENT(IN)    :: TDOF_MSG          ! Message to be printed out regarding at what point in the run the TDOF,I
!                                                            tables are printed out
      CHARACTER(  5*BYTE)             :: SET_NAME          ! A data set name for output purposes

      INTEGER(LONG)                   ::  A_SET_COL        ! Col no. in array TDOF where the  A-set is (from subr TDOF_COL_NUM)
      INTEGER(LONG)                   ::  F_SET_COL        ! Col no. in array TDOF where the  F-set is (from subr TDOF_COL_NUM)
      INTEGER(LONG)                   ::  G_SET_COL        ! Col no. in array TDOF where the  G-set is (from subr TDOF_COL_NUM)
      INTEGER(LONG)                   ::  L_SET_COL        ! Col no. in array TDOF where the  L-set is (from subr TDOF_COL_NUM)
      INTEGER(LONG)                   ::  M_SET_COL        ! Col no. in array TDOF where the  M-set is (from subr TDOF_COL_NUM)
      INTEGER(LONG)                   ::  N_SET_COL        ! Col no. in array TDOF where the  N-set is (from subr TDOF_COL_NUM)
      INTEGER(LONG)                   ::  O_SET_COL        ! Col no. in array TDOF where the  O-set is (from subr TDOF_COL_NUM)
      INTEGER(LONG)                   ::  R_SET_COL        ! Col no. in array TDOF where the  R-set is (from subr TDOF_COL_NUM)
      INTEGER(LONG)                   ::  S_SET_COL        ! Col no. in array TDOF where the  S-set is (from subr TDOF_COL_NUM)
      INTEGER(LONG)                   :: SA_SET_COL        ! Col no. in array TDOF where the SA-set is (from subr TDOF_COL_NUM)
      INTEGER(LONG)                   :: SB_SET_COL        ! Col no. in array TDOF where the SB-set is (from subr TDOF_COL_NUM)
      INTEGER(LONG)                   :: SE_SET_COL        ! Col no. in array TDOF where the SE-set is (from subr TDOF_COL_NUM)
      INTEGER(LONG)                   :: SG_SET_COL        ! Col no. in array TDOF where the SG-set is (from subr TDOF_COL_NUM)
      INTEGER(LONG)                   :: SZ_SET_COL        ! Col no. in array TDOF where the SZ-set is (from subr TDOF_COL_NUM)
      INTEGER(LONG)                   :: U1_SET_COL        ! Col no. in array TDOF where the U1-set is (from subr TDOF_COL_NUM)
      INTEGER(LONG)                   :: U2_SET_COL        ! Col no. in array TDOF where the U2-set is (from subr TDOF_COL_NUM)
      INTEGER(LONG)                   :: I,J               ! DO loop indices
      INTEGER(LONG)                   :: I_USET_U1         ! Counter for USET U1
      INTEGER(LONG)                   :: I_USET_U2         ! Counter for USET U2
      INTEGER(LONG)                   :: IGRID             ! Internal grid number
      INTEGER(LONG)                   :: IROW              ! Row number in array TDOF or TDOFI
      INTEGER(LONG)                   :: NUM_COMPS         ! Number of displ components (1 for SPOINT, 6 for physical grid)
      INTEGER(LONG), PARAMETER        :: SUBR_BEGEND = DOF_PROC_BEGEND

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9001) SUBR_NAME,TSEC
 9001    FORMAT(1X,A,' BEGN ',F10.3)
      ENDIF

! **********************************************************************************************************************************
      WRITE(SC1, * ) '     TDOF PROC'

! Call routine to calc what row in TDOF each grid's data begins

      CALL CALC_TDOF_ROW_START ( 'Y' )

! First, set NDOFG = LDOFG. It will be counted later.

      NDOFG = LDOFG

! Get column numbers for array TDOF different displ sets

      CALL TDOF_COL_NUM ( 'G ',  G_SET_COL )
      CALL TDOF_COL_NUM ( 'N ',  N_SET_COL )
      CALL TDOF_COL_NUM ( 'M ',  M_SET_COL )
      CALL TDOF_COL_NUM ( 'F ',  F_SET_COL )
      CALL TDOF_COL_NUM ( 'S ',  S_SET_COL )
      CALL TDOF_COL_NUM ( 'SA', SA_SET_COL )
      CALL TDOF_COL_NUM ( 'SB', SB_SET_COL )
      CALL TDOF_COL_NUM ( 'SG', SG_SET_COL )
      CALL TDOF_COL_NUM ( 'SZ', SZ_SET_COL )
      CALL TDOF_COL_NUM ( 'SE', SE_SET_COL )
      CALL TDOF_COL_NUM ( 'O ',  O_SET_COL )
      CALL TDOF_COL_NUM ( 'A ',  A_SET_COL )
      CALL TDOF_COL_NUM ( 'R ',  R_SET_COL )
      CALL TDOF_COL_NUM ( 'L ',  L_SET_COL )
      CALL TDOF_COL_NUM ( 'U1', U1_SET_COL )
      CALL TDOF_COL_NUM ( 'U2', U2_SET_COL )

! Set 1st 4 cols of TDOF (actual grid ID - component number, internal grid ID - component number)

      IROW = 0
      DO I = 1,NGRID
         WRITE(SC1,22345,ADVANCE='NO') '       Process col 1-4 of TDOF', I, NGRID, CR13
         CALL GET_GRID_NUM_COMPS ( GRID_ID(INV_GRID_SEQ(I)), NUM_COMPS, SUBR_NAME )
         DO J = 1,NUM_COMPS
            IROW = IROW + 1
            TDOF(IROW,1) = GRID_ID(I)
            TDOF(IROW,2) = J
            TDOF(IROW,3) = GRID_SEQ(I)
            TDOF(IROW,4) = J
         ENDDO
      ENDDO

! Calc TDOF for G-set (col 5 in TDOF). We can do this at this point since all components go in G-set.

      NDOFG = 0
      DO I=1,NGRID
         WRITE(SC1,22345,ADVANCE='NO') '       Process G -set         ', I, NGRID, CR13
         IGRID = INV_GRID_SEQ(I)
         CALL GET_GRID_NUM_COMPS ( GRID_ID(INV_GRID_SEQ(I)), NUM_COMPS, SUBR_NAME )
         DO J=1,NUM_COMPS
            IROW = TDOF_ROW_START(IGRID) + J - 1
            NDOFG = NDOFG + 1
            IF (NDOFG > LDOFG) CALL ARRAY_SIZE_ERROR_1 ( SUBR_NAME, LDOFG, 'TDOF' )
            TDOF(IROW,G_SET_COL) = NDOFG
         ENDDO
      ENDDO

! Make sure NDOFG exactly equals LDOFG (otherwise coding error)

      IF (NDOFG /= LDOFG) THEN
         FATAL_ERR = FATAL_ERR + 1
         WRITE(ERR,1311) SUBR_NAME, NDOFG, LDOFG
         WRITE(F06,1311) SUBR_NAME, NDOFG, LDOFG
         CALL OUTA_HERE ( 'Y' )
      ENDIF

! Put M-set DOF numbers into TDOF (at M_SET_COL)

      IF (NDOFM > 0) THEN
         NDOFM = 0
         DO I=1,NGRID
            WRITE(SC1,22345,ADVANCE='NO') '       Process M -set         ', I, NGRID, CR13
            IGRID = INV_GRID_SEQ(I)
            CALL GET_GRID_NUM_COMPS ( GRID_ID(INV_GRID_SEQ(I)), NUM_COMPS, SUBR_NAME )
            DO J=1,NUM_COMPS
               IF (TSET(IGRID,J) == 'M ') THEN
                  IROW = TDOF_ROW_START(IGRID) + J - 1
                  NDOFM = NDOFM + 1
                  TDOF(IROW,M_SET_COL) = NDOFM
               ENDIF
            ENDDO
         ENDDO
      ENDIF

! Put SA-set DOF numbers into TDOF (at SA_SET_COL)

      IF (NDOFSA > 0) THEN
         NDOFSA = 0
         DO I=1,NGRID
            WRITE(SC1,22345,ADVANCE='NO') '       Process SA-set         ', I, NGRID, CR13
            IGRID = INV_GRID_SEQ(I)
            CALL GET_GRID_NUM_COMPS ( GRID_ID(INV_GRID_SEQ(I)), NUM_COMPS, SUBR_NAME )
            DO J=1,NUM_COMPS
               IF (TSET(IGRID,J) == 'SA') THEN
                  IROW = TDOF_ROW_START(IGRID) + J - 1
                  NDOFSA = NDOFSA + 1
                  TDOF(IROW,SA_SET_COL) = NDOFSA
               ENDIF
            ENDDO
         ENDDO
      ENDIF

! Put SB-set DOF numbers into TDOF (at SB_SET_COL)

      IF (NDOFSB > 0) THEN
         NDOFSB = 0
         DO I=1,NGRID
            WRITE(SC1,22345,ADVANCE='NO') '       Process SB-set         ', I, NGRID, CR13
            IGRID = INV_GRID_SEQ(I)
            CALL GET_GRID_NUM_COMPS ( GRID_ID(INV_GRID_SEQ(I)), NUM_COMPS, SUBR_NAME )
            DO J=1,NUM_COMPS
               IF (TSET(IGRID,J) == 'SB') THEN
                  IROW = TDOF_ROW_START(IGRID) + J - 1
                  NDOFSB = NDOFSB + 1
                  TDOF(IROW,SB_SET_COL) = NDOFSB
               ENDIF
            ENDDO
         ENDDO
      ENDIF

! Put SG-set DOF numbers into TDOF (at SG_SET_COL)

      IF (NDOFSG > 0) THEN
         NDOFSG = 0
         DO I=1,NGRID
            WRITE(SC1,22345,ADVANCE='NO') '       Process SG-set         ', I, NGRID, CR13
            IGRID = INV_GRID_SEQ(I)
            CALL GET_GRID_NUM_COMPS ( GRID_ID(INV_GRID_SEQ(I)), NUM_COMPS, SUBR_NAME )
            DO J=1,NUM_COMPS
               IF (TSET(IGRID,J) == 'SG') THEN
                  IROW = TDOF_ROW_START(IGRID) + J - 1
                  NDOFSG = NDOFSG + 1
                  TDOF(IROW,SG_SET_COL) = NDOFSG
               ENDIF
            ENDDO
         ENDDO
      ENDIF

! Put SE-set DOF numbers into TDOF (at SE_SET_COL)

      IF (NDOFSE > 0) THEN
         NDOFSE = 0
         DO I=1,NGRID
            WRITE(SC1,22345,ADVANCE='NO') '       Process SE-set         ', I, NGRID, CR13
            IGRID = INV_GRID_SEQ(I)
            CALL GET_GRID_NUM_COMPS ( GRID_ID(INV_GRID_SEQ(I)), NUM_COMPS, SUBR_NAME )
            DO J=1,NUM_COMPS
               IF (TSET(IGRID,J) == 'SE') THEN
                  IROW = TDOF_ROW_START(IGRID) + J - 1
                  NDOFSE = NDOFSE + 1
                  TDOF(IROW,SE_SET_COL) = NDOFSE
               ENDIF
            ENDDO
         ENDDO
      ENDIF

! Put O-set DOF numbers into TDOF (at O_SET_COL)

      IF (NDOFO > 0) THEN
         NDOFO = 0
         DO I=1,NGRID
            WRITE(SC1,22345,ADVANCE='NO') '       Process O -set         ', I, NGRID, CR13
            IGRID = INV_GRID_SEQ(I)
            CALL GET_GRID_NUM_COMPS ( GRID_ID(INV_GRID_SEQ(I)), NUM_COMPS, SUBR_NAME )
            DO J=1,NUM_COMPS
               IF (TSET(IGRID,J) == 'O ') THEN
                  IROW = TDOF_ROW_START(IGRID) + J - 1
                  NDOFO = NDOFO + 1
                  TDOF(IROW,O_SET_COL) = NDOFO
               ENDIF
            ENDDO
         ENDDO
      ENDIF

! Put R-set DOF numbers into TDOF (at R_SET_COL)

      IF (NDOFR > 0) THEN
         NDOFR = 0
         DO I=1,NGRID
            IGRID = INV_GRID_SEQ(I)
            WRITE(SC1,22345,ADVANCE='NO') '       Process R -set         ', I, NGRID, CR13
            CALL GET_GRID_NUM_COMPS ( GRID_ID(INV_GRID_SEQ(I)), NUM_COMPS, SUBR_NAME )
            DO J=1,NUM_COMPS
               IF (TSET(IGRID,J) == 'R ') THEN
                  IROW = TDOF_ROW_START(IGRID) + J - 1
                  NDOFR = NDOFR + 1
                  TDOF(IROW,R_SET_COL) = NDOFR
               ENDIF
            ENDDO
         ENDDO
      ENDIF

! Remove code for L-set here. The L-set is done below also

! Calc TDOF for N-set based on G-set minus M-set = S-set + O-set + R-set + L-set

         NDOFN = 0
         DO I=1,NGRID
            WRITE(SC1,22345,ADVANCE='NO') '       Process N -set         ', I, NGRID, CR13
            IGRID = INV_GRID_SEQ(I)
            CALL GET_GRID_NUM_COMPS ( GRID_ID(INV_GRID_SEQ(I)), NUM_COMPS, SUBR_NAME )
            DO J=1,NUM_COMPS
               IROW  = TDOF_ROW_START(IGRID) + J - 1
               IF ((TDOF(IROW,G_SET_COL) > 0) .AND. (TDOF(IROW,M_SET_COL) == 0)) THEN
                  NDOFN = NDOFN + 1
                  TDOF(IROW,N_SET_COL) = NDOFN
               ELSE
                  TDOF(IROW,N_SET_COL) = 0
               ENDIF
            ENDDO
         ENDDO

! Calc DOF'S in SZ-set (all zero SPC's) based on SA + SB + SG

      IF ((NDOFSA > 0) .OR. (NDOFSB > 0) .OR. (NDOFSG > 0)) THEN
         NDOFSZ = 0
         DO I=1,NGRID
            WRITE(SC1,22345,ADVANCE='NO') '       Process SZ-set         ', I, NGRID, CR13
            IGRID = INV_GRID_SEQ(I)
            CALL GET_GRID_NUM_COMPS ( GRID_ID(INV_GRID_SEQ(I)), NUM_COMPS, SUBR_NAME )
            DO J=1,NUM_COMPS
               IROW  = TDOF_ROW_START(IGRID) + J - 1
               IF ((TDOF(IROW,SA_SET_COL) > 0) .OR. (TDOF(IROW,SB_SET_COL) > 0) .OR. (TDOF(IROW,SG_SET_COL) > 0)) THEN
                  NDOFSZ = NDOFSZ + 1
                  TDOF(IROW,SZ_SET_COL) = NDOFSZ
               ELSE
                  TDOF(IROW,SZ_SET_COL) = 0
               ENDIF
            ENDDO
         ENDDO
      ENDIF

! Calc DOF'S in S-set based on SZ + SE

      IF ((NDOFSZ > 0) .OR. (NDOFSE > 0)) THEN
         NDOFS = 0
         DO I=1,NGRID
            WRITE(SC1,22345,ADVANCE='NO') '       Process S -set         ', I, NGRID, CR13
            IGRID = INV_GRID_SEQ(I)
            CALL GET_GRID_NUM_COMPS ( GRID_ID(INV_GRID_SEQ(I)), NUM_COMPS, SUBR_NAME )
            DO J=1,NUM_COMPS
               IROW  = TDOF_ROW_START(IGRID) + J - 1
               IF ((TDOF(IROW,SZ_SET_COL) > 0) .OR. (TDOF(IROW,SE_SET_COL) > 0)) THEN
                  NDOFS = NDOFS + 1
                  TDOF(IROW,S_SET_COL) = NDOFS
               ELSE
                  TDOF(IROW,S_SET_COL) = 0
               ENDIF
            ENDDO
         ENDDO
      ENDIF

! Calc TDOF for F-set based on N-set minus S-set

         NDOFF = 0
         DO I=1,NGRID
            WRITE(SC1,22345,ADVANCE='NO') '       Process F -set         ', I, NGRID, CR13
            IGRID = INV_GRID_SEQ(I)
            CALL GET_GRID_NUM_COMPS ( GRID_ID(INV_GRID_SEQ(I)), NUM_COMPS, SUBR_NAME )
            DO J=1,NUM_COMPS
               IROW  = TDOF_ROW_START(IGRID) + J - 1
               IF ((TDOF(IROW,N_SET_COL) > 0) .AND. (TDOF(IROW,S_SET_COL) == 0)) THEN
                  NDOFF = NDOFF + 1
                  TDOF(IROW,F_SET_COL) = NDOFF
               ELSE
                  TDOF(IROW,F_SET_COL) = 0
               ENDIF
            ENDDO
         ENDDO

! Calc TDOF for A-set based on F-set minus O-set

         NDOFA = 0
         DO I=1,NGRID
            WRITE(SC1,22345,ADVANCE='NO') '       Process A -set         ', I, NGRID, CR13
            IGRID = INV_GRID_SEQ(I)
            CALL GET_GRID_NUM_COMPS ( GRID_ID(INV_GRID_SEQ(I)), NUM_COMPS, SUBR_NAME )
            DO J=1,NUM_COMPS
               IROW  = TDOF_ROW_START(IGRID) + J - 1
               IF ((TDOF(IROW,F_SET_COL) > 0) .AND. (TDOF(IROW,O_SET_COL) == 0)) THEN
                  NDOFA = NDOFA + 1
                  TDOF(IROW,A_SET_COL) = NDOFA
               ELSE
                  TDOF(IROW,A_SET_COL) = 0
               ENDIF
            ENDDO
         ENDDO

! Calc TDOF for L-set based on A-set minus R-set

      NDOFL = 0
      DO I=1,NGRID
         WRITE(SC1,22345,ADVANCE='NO') '       Process L -set         ', I, NGRID, CR13
         IGRID = INV_GRID_SEQ(I)
         CALL GET_GRID_NUM_COMPS ( GRID_ID(INV_GRID_SEQ(I)), NUM_COMPS, SUBR_NAME )
         DO J=1,NUM_COMPS
            IROW  = TDOF_ROW_START(IGRID) + J - 1
            IF ((TDOF(IROW,A_SET_COL) > 0) .AND. (TDOF(IROW,R_SET_COL) == 0)) THEN
               NDOFL = NDOFL + 1
               TDOF(IROW,L_SET_COL) = NDOFL
            ELSE
               TDOF(IROW,L_SET_COL) = 0
            ENDIF
         ENDDO
      ENDDO

! Calc TDOF for USET U1-set and put DOF numbers into TDOF at U1_SET_COL

      IF (NUM_USET_U1 > 0) THEN
         I_USET_U1 = 0
         DO I=1,NGRID
            WRITE(SC1,22345,ADVANCE='NO') '       Process U1-set         ', I, NGRID, CR13
            IGRID = INV_GRID_SEQ(I)
            CALL GET_GRID_NUM_COMPS ( GRID_ID(INV_GRID_SEQ(I)), NUM_COMPS, SUBR_NAME )
            DO J=1,NUM_COMPS
               IF (USET(IGRID,J) == 'U1') THEN
                  IROW = TDOF_ROW_START(IGRID) + J - 1
                  I_USET_U1 = I_USET_U1 + 1
                  TDOF(IROW,U1_SET_COL) = I_USET_U1
               ENDIF
            ENDDO
         ENDDO
      ENDIF

! Calc TDOF for USET U2-set and put DOF numbers into TDOF at U2_SET_COL

      IF (NUM_USET_U2 > 0) THEN
         I_USET_U2 = 0
         DO I=1,NGRID
            WRITE(SC1,22345,ADVANCE='NO') '       Process U2-set         ', I, NGRID, CR13
            IGRID = INV_GRID_SEQ(I)
            CALL GET_GRID_NUM_COMPS ( GRID_ID(INV_GRID_SEQ(I)), NUM_COMPS, SUBR_NAME )
            DO J=1,NUM_COMPS
               IF (USET(IGRID,J) == 'U2') THEN
                  IROW = TDOF_ROW_START(IGRID) + J - 1
                  I_USET_U2 = I_USET_U2 + 1
                  TDOF(IROW,U2_SET_COL) = I_USET_U2
               ENDIF
            ENDDO
         ENDDO
      ENDIF

! Sort TDOF so that G-set DOF's are in numerical order

      WRITE(SC1,12345,ADVANCE='NO') '       Setting up to get TDOFI', CR13
      DO I=1,NDOFG
         DO J=1,MTDOF
            TDOFI(I,J) = TDOF(I,J)
         ENDDO
      ENDDO

      WRITE(SC1,12345,ADVANCE='NO') '       Sort TDOF to get TDOFI ', CR13
      CALL SORT_TDOF ( SUBR_NAME, 'TDOF', NDOFG, TDOFI, G_SET_COL )

! Table TDOF is printed in the F06 file if B.D. PARAM PRTDOF = 1 or 3

      IF (PRTDOF > 0) THEN
         CALL WRITE_TDOF ( TDOF_MSG )
      ENDIF

! Make sure that an R-set has been defined if this is CB soln

      IF (SOL_NAME(1:12) == 'GEN CB MODEL') THEN
         IF (NDOFR < 6) THEN
            SET_NAME = 'R-SET'
            WRITE(ERR,1312) SET_NAME, NDOFR
            WRITE(F06,1312) SET_NAME, NDOFR
            FATAL_ERR = FATAL_ERR + 1
            CALL OUTA_HERE ( 'Y' )
         ENDIF
      ENDIF

! Make sure that if the R set exists it has at least 6 DOF's defined

      IF (NDOFR > 0) THEN
         IF (NDOFR < 6) THEN
            SET_NAME = 'R-SET'
            WRITE(ERR,1304) SET_NAME, NDOFR
            WRITE(F06,1304) SET_NAME, NDOFR
            FATAL_ERR = FATAL_ERR + 1
            CALL OUTA_HERE ( 'Y' )
         ENDIF
      ENDIF

! If this is an eigen analysis (or CB) make sure that some spec for the number of modes to be found exists

      IF ((SOL_NAME(1:5) == 'MODES') .OR. (SOL_NAME(1:12) == 'GEN CB MODEL')) THEN
         IF ((NDOFL > EIGESTL) .AND. (EIG_N2 == 0)) THEN
            FATAL_ERR = FATAL_ERR + 1
            WRITE(ERR,1313) NDOFL, EIGESTL
            WRITE(F06,1313) NDOFL, EIGESTL
            CALL OUTA_HERE ( 'Y' )
         ENDIF
      ENDIF

      WRITE(SC1,*) CR13

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9002) SUBR_NAME,TSEC
 9002    FORMAT(1X,A,' END  ',F10.3)
      ENDIF

      RETURN

! **********************************************************************************************************************************
 1304 FORMAT(' *ERROR  1304: AN ',A,' HAS BEEN DEFINED VIA SUPORT BULK DATA ENTRIES AND THEREFORE MUST HAVE AT LEAST 6 DOF IN IT.' &
                        ,/,14X,' HOWEVER, THERE WERE ONLY ',I2,' DOF DEFINED IN THIS SET')

 1308 FORMAT('  Array TDOF_ROW_START gives the row in array TDOF where the DOF data begins for the GRID listed below',//,          &
             '                  I     GRID   TDOF_ROW_START',/)

 1311 FORMAT(' *ERROR  1311: PROGRAMMING ERROR IN SUBROUTINE ',A                                                                   &
                    ,/,14X,' NDOFG = ',I12,' CALCULATED IN THIS SUBR MUST AGREE WITH LDOFG = ',I12,'  DETERMINED EARLIER')

 1312 FORMAT(' *ERROR  1312: FOR SOL = ''GEN CB MODEL'' THERE MUST BE AN ',A,'-SET WITH AT LEAST NDOFR = 6 DOF''s.'                &
                    ,/,14X,' HOWEVER ONLY ',I1,' DOF''s WERE DEFINED ON BULK DATA SUPORT ENTRIES')

 1313 format(' *ERROR  1313: FOR SOL = "MODES" OR "GEN CB MODEL" THE EIGRL ENTRY MUST HAVE THE NUMBER OF DESIRED MODES > 0 OR THE',&
                           ' PROBLEM DOF SIZE'                                                                                &
                    ,/,14X,' (NDOFL = ',I8,') MUST BE LESS THAT PARAM EIGESTL = ',I8,' (OR USE LARGER VALUE FOR PARAM EIGESTL)')

12345 FORMAT(A, A)

22345 FORMAT(A, ', I = ', I8, ' of ', I8, A)

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

      END SUBROUTINE TDOF_PROC

max_stars_repo_name: luzpaz/MYSTRAN

max_issues_repo_name: luzpaz/MYSTRAN

max_forks_repo_name: luzpaz/MYSTRAN


module dateutil

character(len=512) :: input_file
character(len=512) :: output_file
character(len=32) :: date_option
character(len=32) :: date_format
logical :: iWantMilliseconds

end module

!==================================================================================================!

program main
!----
! Utility for converting to/from calendar dates and calculating numbers of days
!----

use io, only: stdin, stdout, fileExists
use calendar, only: parse_date, date2jd, jd2date

use dateutil, only: input_file, &
                    output_file, &
                    date_option, &
                    date_format, &
                    iWantMilliseconds

implicit none

! Local variables
integer :: i, luin, luout, ios, date1(7), date2(7), ierr
character(len=512) :: input_line, arg1, arg2
character(len=32) :: string(14)
double precision :: nday, day1, day2

! Parse command line
call gcmdln()

! Check that all necessary inputs/outputs defined
if (date_option.eq.'') then
    call usage('dateutil: no date computation option defined')
endif

! Check that input file exists, if defined
if (input_file.ne.'') then
    if (.not.fileExists(input_file)) then
        call usage('dateutil: no input file found named '//trim(input_file))
    endif
endif

! Open input and output streams
if (input_file.eq.'') then
    luin = stdin
else
    luin = 11
    open(unit=luin,file=input_file,status='old')
endif

if (output_file.eq.'') then
    luout = stdout
else
    luout = 12
    open(unit=luout,file=output_file,status='unknown')
endif

! Calculate the new date or the number of days between two dates

do
    ! Read input line
    read(luin,'(A)',iostat=ios) input_line
    if (ios.ne.0) then
        exit
    endif

    ! Split it into two arguments, which depends on formatting
    ! This is an annoying block...
    if (date_format.eq.'YYYY-MM-DD' .or. date_format.eq.'YYYY-MM-DDT') then
        read(input_line,*) arg1, arg2

    elseif (date_format.eq.'YYYY MM DD') then
        if (date_option.eq.'nday') then
            read(input_line,*) (string(i),i=1,6)
            write(arg1,*) (string(i),i=1,3)
            write(arg2,*) (string(i),i=4,6)
        elseif (date_option.eq.'date') then
            read(input_line,*) (string(i),i=1,4)
            write(arg1,*) (string(i),i=1,3)
            write(arg2,*) string(4)
        endif

    elseif (date_format.eq.'YYYYMMDD') then
        read(input_line,*) arg1, arg2

    elseif (date_format.eq.'YYYY-MM-DDTHH:MM:SS') then
        read(input_line,*) arg1, arg2

    elseif (date_format.eq.'YYYY MM DD HH MM SS') then
        if (date_option.eq.'nday') then
            read(input_line,*) (string(i),i=1,12)
            write(arg1,*) (string(i),i=1,6)
            write(arg2,*) (string(i),i=7,12)
        elseif (date_option.eq.'date') then
            read(input_line,*) (string(i),i=1,7)
            write(arg1,*) (string(i),i=1,6)
            write(arg2,*) string(7)
        endif

    elseif (date_format.eq.'YYYYMMDDHHMMSS') then
        read(input_line,*) arg1, arg2

    elseif (date_format.eq.'YYYY-MM-DDTHH:MM:SS.MMM') then
        read(input_line,*) arg1, arg2

    elseif (date_format.eq.'YYYY MM DD HH MM SS MMM') then
        if (date_option.eq.'nday') then
            read(input_line,*) (string(i),i=1,14)
            write(arg1,*) (string(i),i=1,7)
            write(arg2,*) (string(i),i=8,14)
        elseif (date_option.eq.'date') then
            read(input_line,*) (string(i),i=1,8)
            write(arg1,*) (string(i),i=1,7)
            write(arg2,*) string(8)
        endif

    elseif (date_format.eq.'YYYYMMDDHHMMSSMMM') then
        read(input_line,*) arg1, arg2

    else
        call usage('dateutil: date_format "'//trim(date_format)//'" not recognized')
    endif

    ! The first argument is always a date
    call parse_date(arg1,date1,date_format,ierr)
    call date2jd(date1,day1)

    ! The second argument depends on the mode
    if (date_option.eq.'nday') then
        call parse_date(arg2,date2,date_format,ierr)
        call date2jd(date2,day2)
        write(luout,*) day2-day1
    elseif (date_option.eq.'date') then
        read(arg2,*) nday
        day2 = day1+nday
        call jd2date(day2,date2)
        if (len(trim(date_format)).le.10) then
            write(luout,1001) (date2(i),i=1,3)
        else
            if (iWantMilliseconds) then
                write(luout,1002) (date2(i),i=1,7)
            else
                if (date2(7).ge.500) then
                    date2(6) = date2(6) + 1
                endif
                if (date2(6).ge.60) then
                    date2(5) = date2(5) + 1
                    date2(6) = date2(6) - 60
                endif
                if (date2(5).ge.60) then
                    date2(4) = date2(4) + 1
                    date2(5) = date2(5) - 60
                endif
                if (date2(4).ge.24) then
                    date2(3) = date2(3) + 1
                    date2(4) = date2(4) - 24
                endif
                write(luout,1003) (date2(i),i=1,6)
            endif
        endif
        1001 format(I0.4,"-",I0.2,"-",I0.2)
        1002 format(I0.4,"-",I0.2,"-",I0.2,"T",I0.2,":",I0.2,":",I0.2,".",I0.3)
        1003 format(I0.4,"-",I0.2,"-",I0.2,"T",I0.2,":",I0.2,":",I0.2)
    endif
enddo

! Shut it down
if (input_file.ne.'') then
    close(luin)
endif
if (output_file.ne.'') then
    close(luout)
endif

end

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

subroutine gcmdln()

use dateutil, only: input_file, &
                    output_file, &
                    date_option, &
                    date_format, &
                    iWantMilliseconds

implicit none

! Local variables
integer :: i, narg
character(len=512) :: tag

! Initialize control parameters
input_file = ''
output_file = ''
date_option = ''
date_format = 'YYYY-MM-DDTHH:MM:SS'
iWantMilliseconds = .false.

! Number of arguments
narg = command_argument_count()
if (narg.eq.0) then
  call usage('')
endif

! Read arguments
i = 1
do while (i.le.narg)

    call get_command_argument(i,tag)

    if (trim(tag).eq.'-nday') then
        date_option = 'nday'
    elseif (tag(1:5).eq.'-date') then
        date_option = 'date'
    elseif (trim(tag).eq.'-f') then
        i = i + 1
        call get_command_argument(i,input_file)
    elseif (trim(tag).eq.'-o') then
        i = i + 1
        call get_command_argument(i,output_file)
    elseif (trim(tag).eq.'-format') then
        i = i + 1
        call get_command_argument(i,date_format)
    else
        call usage('dateutil: No option '//tag)
    endif

    i = i + 1

enddo

return
end subroutine

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

subroutine usage(str)

use io, only: stderr
implicit none

character(len=*) :: str

if (str.ne.'') then
    write(stderr,*) trim(str)
    write(stderr,*)
endif

write(stderr,*) 'Usage: dateutil ...options...'
write(stderr,*)
write(stderr,*) '-f IFILE      Input file (default reads stdin)'
write(stderr,*) '-o OFILE      Output file (default prints to stdout)'
write(stderr,*) '-nday         Compute number of days from DATE1 DATE2'
write(stderr,*) '-date         Compute calendar date from DATE1 NDAY'
write(stderr,*) '-format FMT   Input date format (default: YYYY-MM-DDTHH:MM:SS)'
write(stderr,*)
write(stderr,*) 'dateutil recognizes the following date formats:'
write(stderr,*) '    YYYY-MM-DD[THH:MM:SS[.MMM]]'
write(stderr,*) '    YYYY MM DD[ HH MM SS[ MMM]]'
write(stderr,*) '    YYYYMMDD[HHMMSS[MMM]]'
write(stderr,*) 'Use quotations when using second format to read spaces correctly'
write(stderr,*)

call error_exit(1)
end subroutine


!------------------------------------------------------------------------------
!
! `prms` Source File
!
! prms.f90 source function file. This file contains all of the subroutines
! responsible for defining physical constants, reading input parameters from
! sige.in, and broadcasting all of these parameters to the parallel
! processes. Also contains subroutines for creating the initial diamond cubic
! crystal structure for Si and for swapping between the fast/slow timestepping
! and thermostat regimes based on current ion energy.
!
!------------------------------------------------------------------------------

MODULE prms

  USE parallel

  IMPLICIT none

  integer              :: Natm                ! number of atoms
  integer              :: Nlj 				  ! number of lj atoms
  integer			         :: Nsg				  ! number of si-ge atoms
  integer              :: Nrand                 ! number of random points
  integer              :: impact               ! impact number
  integer              :: startimpact          ! starting impact number, only valid if which=0, if 0, start from last, else strat from that impact using restart file.
  integer              :: ions                ! write impact number to restart file
  integer, dimension(3):: Npc                 ! number of pieces in the x,y,z direction for atom division
  integer              :: Nt                  ! number timesteps
  real                 :: fwhm              !full width half max of ion distribution
  real                 :: phiz, phixy         !angles of incidence for ion
  real                 :: Ts                  ! timestep
  real                 :: Ts_i, Ts_r          ! initial and secondary timestep
  real,dimension(3)    :: Lb                  ! box sixe
  real,dimension(3)    :: iLb                 ! inverse box sixe
  real, dimension(4)   :: outz                ! places to collect the out of simulation atoms, in multiples of Lb(3)/10
  real, dimension(2)   :: atomz               ! keep out of sim atoms
  real                 :: ionz, knockz ! place to keep out-of-sim atoms, repository of ions before sim, starting position of ions
  real                 :: sidewidth           ! sidewidth to control temperature in nanometers
  real                 :: ke_limit_2            ! kinetic energy threshhold, for which the atoms are considered slow
  real                 :: ke_limit_1          ! threshold to increase timestep
  logical              :: atom_is_slow, atom_is_fast
  real                 :: restart_time
  integer, allocatable, dimension(:) :: xsides

  real                 :: A                   ! S-W prm
  real                 :: B                   ! S-W prm
  real                 :: psi                 ! S-W prm
  real                 :: qsi                 ! S-W prm
  real                 :: al                  ! S-W prm
  real                 :: lambdaSi            ! S-W prm
  real                 :: lambdaGe            ! S-W prm
  real                 :: gamma               ! S-W prm
  real                 :: epsSi               ! S-W prm
  real                 :: epsGe               ! S-W prm
  real                 :: sigmaSi             ! S-W prm
  real                 :: sigmaGe             ! S-W prm
  real                 :: massSi
  real                 :: massGe
  real                 :: massGa              ! mass of Ga atoms

  real                 :: ZSi                 ! Atomic number Si
  real                 :: ZGe                 ! Atomic number Ge
  real                 :: ZGa                 ! Atomic number Ar
  real                 :: Amol, Bmol, Cmol    ! Moliere prms
  real                 :: alpha_mol           ! Moliere prm
  real                 :: beta_mol            ! Moliere prm
  real                 :: gamma_mol           ! Moliere prm
  real,dimension(3,3)  :: aamol               ! A Moliere
  real,dimension(3)    :: Zmol                ! Z Moliere
  real                 :: ie0                 ! one over eps0 (constant for unit matching)

  real,dimension(2,2)  :: isigma              ! 1./sigma
  real,dimension(2,2)  :: eps                 ! eps
  real,dimension(2,2,2):: leps                ! lambda * eps
  real,dimension(2,2)  :: epssig              ! eps * sigma

  real                 :: Ttar1               ! Target Temperature 1
  real                 :: Qheat               ! Heat flux for #3
  real                 :: Tau                 ! Termostat prm
  real                 :: Tau_i, Tau_r
  integer              :: Tinit               ! # Berendsen on all

  integer              :: whichthrm           ! 1 - Anderson; 2 - Berendsen; 3 - Energy
  integer              :: whichtemp           ! 1 - ends; 2 - left,middle
  integer              :: freeze              ! 1 - freeze ends, 2 - no freeze ends

  integer              :: Nltemp1             ! Temp control # lat left
  integer              :: Nltemp2             ! Temp control # lat right

  real                 :: nu                  ! probability of collision for And Therm.
  real                 :: Vmax                ! Max boltzman velocity for And Therm.

  real                 :: Uo                  ! Initial velocity

  real                 :: Xint                ! location of interface
  integer              :: Nsp                ! N of monolayers in superlattice layer
  real                 :: Amp                 ! forcing amplitude
  real                 :: alpha               ! wave angle
  real                 :: Xofrac              ! postion (fraction of Lb(1))

  integer              :: atm_out             ! atom position output interval
  integer              :: atm_out_i, atm_out_r
  integer              :: abn_out             ! atom pos-vel output interval
  integer              :: res_out             ! restart file output interval
  integer              :: zen_out             ! zone energy change
  integer              :: eng_out             ! energy output
  integer              :: tmp_out             ! mean temperature output

  real                 :: listfac             ! factor rl/rc for neighborlists
  integer              :: ntlist, ntlist_i, ntlist_r              ! how often to reconstruct the neighborlist
  integer              :: atlist, atlist_i, atlist_r              ! how often updates global X
  integer              :: dslist, dslist_i, dslist_r              ! displaced atom update intercval

  integer              :: Ntable              ! number of values is Fsr table

  integer              :: ranseed             ! random number seed
  integer              :: whichic             ! which initial conditions
  integer              :: amorcrys             ! amorphous or crystalline target? 0-cryst, 1-amorphous
  real 				   :: Teps               ! epsilon of temperature

  integer, parameter :: prm_unit  =   4  ! read unit for parameters
  integer, parameter :: out_unit  =   6  ! stdout
  integer, parameter :: atm_unit  =   7  ! atom locations
  integer, parameter :: rdf_unit  =   8  ! radial distribution function
  integer, parameter :: abn_unit  =  10  ! atoms+velocities in binary
  integer, parameter :: lat_unit  =  11  ! lattice definition
  integer, parameter :: res_unit  =  12  ! restart
  integer, parameter :: wav_unit  =  13  ! perturbed position unit
  integer, parameter :: zen_unit  =  14  ! zone energy change
  integer, parameter :: eng_unit  =  15  ! energy diagnostics
  integer, parameter :: tmp_unit  =  16  ! temperature diagnostics
  integer, parameter :: max_unit  =  17  ! maxwell (?) distribution
  integer, parameter :: egk_unit  =  18  ! heat current
  integer, parameter :: fly_unit  =  19  ! y dissipation
  integer, parameter :: yen_unit  =  20  ! y energy output
  integer, parameter :: pbx_unit  =  21  ! x histories
  integer, parameter :: pby_unit  =  22  ! y histories
  integer, parameter :: bce_unit  =  23  ! bc energy
  integer, parameter :: ion_unit  =  24  ! ion trajectory
  integer, parameter :: ds_unit   =  25  ! atoms trajectory
  integer, parameter :: ran_unit  =  26  ! generated random ion positions
  integer, parameter :: snap_unit =  27  ! snapshot xyz file

  real               :: Pi                         ! Pi (duh)
  real               :: srPi                       ! SQRT(Pi)
  real               :: thrd                       ! 1./3.
  real   , parameter :: kB        =   1.381E-23    ! k-Boltzmann (J/K)
  real   , parameter :: eps0      =   8.854187E-12 ! permittivity constant (F/m = C^2/J m)
  real   , parameter :: ec        =   1.60217646E-19 ! elementary charge (C)
  real   , parameter :: ao        =   0.053E-09    ! Bohr radius

  integer, parameter :: Nmaxg=200      ! RDF num bins
  integer, parameter :: Nmaxw=200      ! M-B num bins
  real               :: eV
  real               :: dti            ! time between subsequent impacts

CONTAINS

!readprms() defines material parameters and reads input from siga.in
  SUBROUTINE readprms ()

    integer   :: ierr
    integer :: Nc1, Nc2, Nc3
    integer, dimension (8) :: seed

    if (myid.eq.0) write(out_unit,*) "READING PARAMETERS"

    Pi = 4.*ATAN(1.0)
    srPi = SQRT(Pi)
    thrd = 1./3.
    ZSi = 14.0
    ZGe = 32.0
    ZGa = 31.0
    massSi = 46.637063E-27
    massGe = 120.6049350E-27
    massGa = 1.157777417E-25

    Amol     = 0.35                     ! Moliere prm
    Bmol     = 0.55                     ! Moliere prm
    Cmol     = 0.10                     ! Moliere prm
    alpha_mol= 0.3                      ! Moliere prm
    beta_mol = 1.2                      ! Moliere prm
    gamma_mol= 6.0                      ! Moliere prm

    A     = 7.049556277                 ! S-W prm
    B     = 0.6022245584                ! S-W prm
    psi   = 4.0                         ! S-W prm
    qsi   = 0.0                         ! S-W prm
    al    = 1.8                         ! S-W prm

    lambdaSi = 21.0                     ! S-W prm
    lambdaGe = 31.0						! S-W prm
    gamma = 1.2                     	! S-W prm
    epsSi = 3.473928E-19                ! S-W prm
    epsGe = 3.085E-19                   ! S-W prm
    sigmaSi = 2.0951E-10                ! S-W prm
    sigmaGe = 2.0951E-10                ! S-W prm, made same as Silicon, since we do not have any germanium in our system, 2.1810E-10 for Ge

    nu = 0.01
    Vmax = 3500.0
    Tinit = 2000
    whichthrm = -2
    whichtemp = -1
    freeze = 0
    Qheat = 8.E-23
    Nltemp1 = 15
    Nltemp2 = 15
    listfac = 1.1   ! this is extra padding of neighborlist. prev 1.05, higher value lists more atoms as neighbors, while they will be given 0 force. Kallol

    abn_out = -1
    zen_out = -1
    Ntable = 12000  ! number of division in table, prev 10000 kallol

    if (myid.eq.0) then

        call DATE_AND_TIME(values=seed)
        open(prm_unit,file='siga.in')

        read(prm_unit,*) amorcrys
        read(prm_unit,*) whichic
        read(prm_unit,*) startimpact
        read(prm_unit,*) eV
        read(prm_unit,*) Nc1,Nc2,Nc3
        read(prm_unit,*) Npc(1), Npc(2), Npc(3)
        read(prm_unit,*) Nlj

        read(prm_unit,*) Nt
        read(prm_unit,*) Ts_i
        read(prm_unit,*) Ts_r

        read(prm_unit,*) Ttar1
        read(prm_unit,*) ke_limit_1
        read(prm_unit,*) ke_limit_2

        read(prm_unit,*) Tau_i
        read(prm_unit,*) Tau_r
        read(prm_unit,*) Teps

        read(prm_unit,*) Uo

        read(prm_unit,*) ntlist_i
        read(prm_unit,*) ntlist_r
        read(prm_unit,*) atlist_i
        read(prm_unit,*) atlist_r
        read(prm_unit,*) dslist_i
        read(prm_unit,*) dslist_r

        read(prm_unit,*) atm_out_i
        read(prm_unit,*) atm_out_r
        read(prm_unit,*) res_out
        read(prm_unit,*) eng_out
        read(prm_unit,*) tmp_out

        read(prm_unit,*) Nrand

        read(prm_unit,*) outz(1)
        read(prm_unit,*) outz(2)
        read(prm_unit,*) outz(3)
        read(prm_unit,*) outz(4)

        read(prm_unit,*) atomz(1)
        read(prm_unit,*) atomz(2)
        read(prm_unit,*) ionz
        read(prm_unit,*) knockz
        read(prm_unit,*) sidewidth
        read(prm_unit,*) dti

        read(prm_unit,*) fwhm

        read(prm_unit,*) phiz
        read(prm_unit,*) phixy

        close(prm_unit)

        if (amorcrys .eq. 0) then
            call pureepi(Nc1, Nc2,Nc3) !calculates Lb, Natm, Nsg of crystalline target
            ! disabling prev line adding new lines
            ! also enable xsides!!!!!!!!!!!!!!!!!#####################
            !Lb = 5.431073E-10
            !Lb(1) = Lb(1)*Nc1
            !Lb(2) = Lb(2)*Nc2
            !Lb(3) = Lb(3)*Nc3*10

        else
            Lb = 5.431073E-10
            Lb(1) = Lb(1)*Nc1
            Lb(2) = Lb(2)*Nc2
            Lb(3) = Lb(3)*Nc3*10
            Natm = Nsg+Nlj
        end if

        ranseed = -seed(7)*seed(8)

    end if

    call MPI_BCAST(amorcrys,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(eV,1,MPI_REAL8,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(Natm,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(Nsg,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(Nc1,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(Nc2,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(Nc3,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(Npc,3,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(Nlj,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(Nt,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(Ts_i,1,MPI_REAL8,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(Ts_r,1,MPI_REAL8,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(Lb,3,MPI_REAL8,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(Ttar1,1,MPI_REAL8,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(ke_limit_1,1,MPI_REAL8,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(ke_limit_2,1,MPI_REAL8,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(Tau_i,1,MPI_REAL8,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(Tau_r,1,MPI_REAL8,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(Teps,1,MPI_REAL8,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(Uo,1,MPI_REAL8,0,MPI_COMM_WORLD,ierr)

    call MPI_BCAST(ntlist_i,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(ntlist_r,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(atlist_i,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(atlist_r,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(dslist_i,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(dslist_r,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)

    call MPI_BCAST(atm_out_i,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(atm_out_r,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(res_out,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(eng_out,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(tmp_out,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)

    call MPI_BCAST(ranseed,1,MPI_REAL8,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(whichic,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(startimpact,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(Nrand,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(outz,4,MPI_REAL8,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(atomz,2,MPI_REAL8,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(ionz,1,MPI_REAL8,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(knockz,1,MPI_REAL8,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(sidewidth,1,MPI_REAL8,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(dti,1,MPI_REAL8,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(fwhm,1,MPI_REAL8,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(phiz,1,MPI_REAL8,0,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(phixy,1,MPI_REAL8,0,MPI_COMM_WORLD,ierr)

    call MPI_BCAST(Lb,3,MPI_REAL8,0,MPI_COMM_WORLD,ierr)

    !add these two lines in production run
    if(myid .ne. 0) allocate (xsides(Natm))
    call MPI_BCAST(xsides,Nsg,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)

    call setupSiGe
    call setupMol
    call ionrun(1)
    iLb = 1./Lb

    !if (myid.eq.Np-1) then
    if(myid .eq. -1) then

       write(out_unit,*)
       write(out_unit,*) " ---- Numerical parameters ----- "
       write(out_unit,*) "Nt      = ",Nt
       write(out_unit,*) "Ts      = ",Ts
       write(out_unit,*) "ntlist  = ",ntlist
       write(out_unit,*) "atlist  = ",atlist
       write(out_unit,*) "listfac = ",listfac
       write(out_unit,*) "ranseed = ",ranseed
       write(out_unit,*) "whichic = ",whichic
       write(out_unit,*) "outz 1-2= ",outz(1), outz(2)
       write(out_unit,*) "outz 3-4= ",outz(3), outz(4)
       write(out_unit,*) "atomz   = ", atomz(1), atomz(2)
       write(out_unit,*) "ionz    = ", ionz
       write(out_unit,*) "knockz  = ", knockz
       write(out_unit,*) "Npc =     ", Npc
       write(out_unit,*) "phiz =     ", phiz
       write(out_unit,*) "phixy =     ", phixy
       write(out_unit,*)
       write(out_unit,*) " ----- Physical parameters ----- "
       write(out_unit,*) "Natm    = ",Natm
       write(out_unit,*) "Nsg     = ",Nsg
       write(out_unit,*) "Nlj     = ",Nlj
       write(out_unit,*) "Nrand   = ",Nrand
       write(out_unit,*) "Uo      = ",Uo
       write(out_unit,*) "Lb      = ",Lb
       write(out_unit,*) "A       = ",A
       write(out_unit,*) "B       = ",B
       write(out_unit,*) "psi     = ",psi
       write(out_unit,*) "qsi     = ",qsi
       write(out_unit,*) "al      = ",al
       write(out_unit,*) "lambdaSi= ",lambdaSi
       write(out_unit,*) "lambdaGe= ",lambdaGe
       write(out_unit,*) "gamma   = ",gamma
       write(out_unit,*) "epsSi   = ",epsSi
       write(out_unit,*) "epsGe   = ",epsGe
       write(out_unit,*) "sigmaSi = ",sigmaSi
       write(out_unit,*) "sigmaGe = ",sigmaGe
       write(out_unit,*) "massSi  = ",massSi
       write(out_unit,*) "massGe  = ",massGe
       write(out_unit,*)
       write(out_unit,*) " --- Temperature parameters ---- "
       write(out_unit,*) "Ttar1   = ",Ttar1
       write(out_unit,*) "Qheat   = ",Qheat
       write(out_unit,*) "Tau     = ",Tau
       write(out_unit,*) "Tinit   = ",Tinit
       write(out_unit,*) "whichthrm = ",whichthrm
       write(out_unit,*) "whichtemp = ",whichtemp
       write(out_unit,*) "freeze  = ",freeze
       write(out_unit,*) "Nltemp1 = ",Nltemp1
       write(out_unit,*) "Nltemp2 = ",Nltemp2
       write(out_unit,*) "nu      = ",nu
       write(out_unit,*) "Vmax    = ",Vmax
       write(out_unit,*)
       write(out_unit,*) " ------ Output parameters ----- "
       write(out_unit,*) "atm_out   = ",atm_out
       write(out_unit,*) "abn_out   = ",atm_out
       write(out_unit,*) "res_out   = ",res_out
       write(out_unit,*) "zen_out   = ",zen_out
       write(out_unit,*) "eng_out   = ",eng_out
       write(out_unit,*) "tmp_out   = ",tmp_out
       write(out_unit,*)
       write(out_unit,*) " ----- Ranging parameters ----- "
       write(out_unit,*) "listfac   = ",listfac
       write(out_unit,*) "Ntlist    = ",Ntlist
       write(out_unit,*) "Ntable    = ",Ntable

    end if

  END SUBROUTINE readprms

!setupMol defines cross interaction coefficients for moliere potential
  SUBROUTINE setupMol

    Zmol(1) = ZSi
    Zmol(2) = ZGe                                                   ! previously    0. zubaer
    Zmol(3) = ZGa

    aamol = 0.
    aamol(1,3) = 0.885*ao*(SQRT(ZSi) + SQRT(ZGa))**(-2./3.)
    aamol(2,3) = 0.885*ao*(SQRT(ZGe) + SQRT(ZGa))**(-2./3.)         ! previously set  aamol(1,3)        zubaer
    aamol(3,1) = aamol(1,3)
    aamol(3,2) = aamol(2,3)
    aamol(3,3) = 0.885*ao*(SQRT(ZGa) + SQRT(ZGa))**(-2./3.)

    ie0 = 1/(4*Pi*eps0)    !  Coulomb force constant
    !  where eps0 is the permittivity of free space

  END SUBROUTINE setupMol

!setupSiGe sets up cross interaction constants for Si-Ge Stillinger-Weber potential
  SUBROUTINE setupSiGe

    eps(1,1) =  epsSi
    eps(1,2) = 0.5*(epsSi+epsGe)
    eps(2,2) =  epsGe
    eps(2,1) = 0.5*(epsSi+epsGe)

    isigma(1,1) = 1./sigmaSi
    isigma(1,2) = 2./(sigmaSi+sigmaGe)
    isigma(2,2) = 1./sigmaGe
    isigma(2,1) = 2./(sigmaSi+sigmaGe)

    leps(1,1,1) = lambdaSi*epsSi
    leps(1,1,2) = 2.*lambdaSi*epsSi/3. + lambdaGe*epsGe/3.
    leps(1,2,1) = 2.*lambdaSi*epsSi/3. + lambdaGe*epsGe/3.
    leps(2,1,1) = 2.*lambdaSi*epsSi/3. + lambdaGe*epsGe/3.
    leps(1,2,2) = lambdaSi*epsSi/3. + 2.*lambdaGe*epsGe/3.
    leps(2,2,1) = lambdaSi*epsSi/3. + 2.*lambdaGe*epsGe/3.
    leps(2,1,2) = lambdaSi*epsSi/3. + 2.*lambdaGe*epsGe/3.
    leps(2,2,2) = lambdaGe*epsGe

    epssig(1,1) = epsSi/sigmaSi
    epssig(1,2) = (epsSi+epsGe)/(sigmaSi+sigmaGe)
    epssig(2,1) = (epsSi+epsGe)/(sigmaSi+sigmaGe)
!    epssig(1,2) = 0.5*(epsSi*sigmaSi+epsGe*sigmaGe)
!    epssig(2,1) = 0.5*(epsSi*sigmaSi+epsGe*sigmaGe)
    epssig(2,2) = epsGe/sigmaGe

  END SUBROUTINE setupSiGe

!pureepi(Nc1,Nc2,Nc3) defines a perfect starting diamond cubic lattice,
!according to unit cell counts Nc1,Nc2,Nc3 in the x,y,z directions respectively
  SUBROUTINE pureepi(Nc1,Nc2,Nc3)
    ! GROWTH DIRECTION is Z
    real,dimension(8,3)                :: uc
    integer                            :: Nc1, Nc2, Nc3, m, mnew
    real, allocatable, dimension(:,:)  :: XX
    integer                            ::  l,i,j,k,cnt, ns
    real                               :: shift, maxx, maxy

    allocate(XX(8*Nc1*Nc2*Nc3,3))

    Lb = 5.431073E-10;
    !Lb(3) = 10E-10;

    ! VALUES r fr MATERIAL ONE with X MLs in [001] DIRECTION
    uc(1,1) = 0.00;         uc(1,2) = 0.00;            uc(1,3) = 0.00
    uc(2,1) = 0.00;         uc(2,2) = 0.50*Lb(2) ;     uc(2,3) = 0.50*Lb(3)
    uc(3,1) = 0.250*Lb(1);  uc(3,2) = 0.250*Lb(2) ;    uc(3,3) = 0.750*Lb(3)
    uc(4,1) = 0.250*Lb(1);  uc(4,2) = 0.750*Lb(2) ;    uc(4,3) = 0.250*Lb(3)
    uc(5,1) = 0.50*Lb(1);   uc(5,2) = 0.0*Lb(2) ;      uc(5,3) = 0.50*Lb(3)
    uc(6,1) = 0.50*Lb(1);   uc(6,2) = 0.50*Lb(2) ;     uc(6,3) = 0.0*Lb(3)
    uc(7,1) = 0.750*Lb(1);  uc(7,2) = 0.250*Lb(2) ;    uc(7,3) = 0.250*Lb(3)
    uc(8,1) = 0.750*Lb(1);  uc(8,2) = 0.750*Lb(2) ;    uc(8,3) = 0.750*Lb(3)

    m = 1
    do i = 1,Nc1
        do j = 1,Nc2
            do k = 1,Nc3
                do l = 1,8
                    XX(m,1) = uc(l,1) + Lb(1)*REAL(i-1)
                    XX(m,2) = uc(l,2) + Lb(2)*REAL(j-1)
                    XX(m,3) = uc(l,3) + Lb(3)*REAL(k-1)
                    m = m+1
                end do
            end do
        end do
    end do

    mnew = m-1
    cnt = 0
    open (2,file = 'si.dat')
    do i = 1,m-1
        write(2,"(3E20.10,I10)")XX(i,:),1
        cnt =cnt+1
    end do
    close(2)

    print *,"Si atoms",m-1
    print*, "Domain size (nm): x,y,z"
    print*,Lb(1)*Nc1*1e9
    print*,Lb(2)*Nc2*1e9
    print*,Lb(3)*Nc3*1e9

    Lb(1) = Lb(1)*Nc1
    Lb(2) = Lb(2)*Nc2
    Lb(3) = Lb(3)*Nc3*10

    Nsg = cnt
    Natm = Nsg+Nlj

    maxx = MAXVAL(XX(:,1))
    maxy = MAXVAL(XX(:,2))

!    open(2, file='si_sides.dat')
    allocate(xsides(Natm))
    xsides = 0
    do i=1,Nsg
        if(XX(i,1).lt.1.0e-10 .or. XX(i,2).lt.1.0e-10 .or. XX(i,1).gt.maxx-1.0e-10 .or. XX(i,2).gt.maxy-1.0e-10) then
            xsides(i) = i
!            write(2,"(I10)")i
        end if
    end do
!    close(2)

  END SUBROUTINE pureepi

!ionrun(tic) adjusts frequency of operations based on whether simulations is
!in fast or slow regimes (by ion energy theshold), tic = 1 or 2 respectively
!controls file write frequency, neighborlist generation frequency, core assignment frequency
!thermostat coefficient, dslist
  SUBROUTINE ionrun(tic)
    integer  :: tic

    if(myid .eq. 0) print*, 'called ionrun', tic

    if(tic .eq. 1) then
        Ts = Ts_i
        atm_out = atm_out_i
        ntlist = ntlist_i
        atlist = atlist_i
        Tau = Tau_i
        dslist = dslist_i
    elseif (tic .eq. 2) then
        Ts = Ts_r
        atm_out = atm_out_r
        ntlist = ntlist_r
        atlist = atlist_r
        Tau = Tau_r
        dslist = dslist_r
    end if

  END SUBROUTINE ionrun

END MODULE prms