text:
Subroutine crnd(irun, px, py, pz, srt, i1, i2, iblock, signn, sig, sigk, xsk1, xsk2, xsk3, xsk4, xsk5, nt, ipert1)
  Parameter (maxstr=150001, maxr=1, amn=0.939457, amp=0.93828, ap1=0.13496, aka=0.498, aphi=1.020, ap2=0.13957, am0=1.232, pi=3.1415926, cutoff=1.8966, avmass=0.9383)
  Parameter (mx=4, my=4, mz=8, mpx=4, mpy=4, mpz=10, mpzp=10)
  Parameter (xmd=1.8756, npdmax=10000)
  Common /aa/r(3, maxstr)
  Common /bb/p(3, maxstr)
  Common /cc/e(maxstr)
  Common /ee/id(maxstr), lb(maxstr)
  Common /ff/f(-mx:mx, -my:my, -mz:mz, -mpx:mpx, -mpy:mpy, -mpz:mpzp)
  Common /gg/dx, dy, dz, dpx, dpy, dpz
  Common /input/nstar, ndirct, dir
  Common /nn/nnn
  Common /bg/betax, betay, betaz, gamma
  Common /run/num
  Common /pa/rpion(3, maxstr, maxr)
  Common /pb/ppion(3, maxstr, maxr)
  Common /pc/epion(maxstr, maxr)
  Common /pd/lpion(maxstr, maxr)
  Common /input1/masspr, massta, iseed, iavoid, dt
  Common /leadng/lb1, px1, py1, pz1, em1, e1, xfnl, yfnl, zfnl, tfnl, px1n, py1n, pz1n, dp1n
  Common /rndf77/nseed
  Common /dpert/dpertt(maxstr, maxr), dpertp(maxstr), dplast(maxstr), dpdcy(maxstr), dpdpi(maxstr, maxr), dpt(maxstr, maxr), dpp1(maxstr, maxr), dppion(maxstr, maxr)
  Common /dpi/em2, lb2
  Common /para8/idpert, npertd, idxsec
  Dimension ppd(3, npdmax), lbpd(npdmax)
  Save
  n12 = 0
  m12 = 0
  iblock = 0
  ntag = 0
  em1 = e(i1)
  em2 = e(i2)
  pr = sqrt(px**2+py**2+pz**2)
  c2 = pz/pr
  x1 = ranart(nseed)
  ianti = 0
  If (lb(i1)<0 .And. lb(i2)<0) ianti = 1
  Call sbbdm(srt, sdprod, ianti, lbm, xmm, pfinal)
  If (idpert==1 .And. ipert1==1) Then
     If (srt<2.012) Return
     If ((iabs(lb(i1))==1 .Or. iabs(lb(i1))==2) .And. (iabs(lb(i2))>=6 .And. iabs(lb(i2))<=13)) Then
        Goto 108
     Else If ((iabs(lb(i2))==1 .Or. iabs(lb(i2))==2) .And. (iabs(lb(i1))>=6 .And. iabs(lb(i1))<=13)) Then
        Goto 108
     Else
        Return
     End If
  End If
  If (x1<=signn/sig) Then
     as = (3.65*(srt-1.8766))**6
     a = 6.0*as/(1.0+as)
     ta = -2.0*pr**2
     x = ranart(nseed)
     t1 = sngl(dlog(dble(1.-x)*dexp(dble(a)*dble(ta))+dble(x)))/a
     c1 = 1.0 - t1/ta
     t1 = 2.0*pi*ranart(nseed)
     iblock = 1
     Goto 107
  Else
     If (srt<2.04) Return
     If (((iabs(lb(i1))==2 .Or. iabs(lb(i2))==2) .And. (lb(i1)*lb(i2))==20) .Or. (lb(i1)*lb(i2))==13) Then
        If (x1<=((signn+sdprod)/sig)) Goto 108
     End If
     prf = sqrt(0.25*srt**2-avmass**2)
     If (em1>1.) Then
        deltam = em1
     Else
        deltam = em2
     End If
     renom = deltam*prf**2/denom(srt, 1.)/pr
     renomn = deltam*prf**2/denom(srt, 2.)/pr
     renom1 = deltam*prf**2/denom(srt, -1.)/pr
     If ((iabs(lb(i1))==2) .And. (iabs(lb(i2))==6)) renom = 0.
     If ((iabs(lb(i2))==2) .And. (iabs(lb(i1))==6)) renom = 0.
     If ((iabs(lb(i1))==1) .And. (iabs(lb(i2))==9)) renom = 0.
     If ((iabs(lb(i2))==1) .And. (iabs(lb(i1))==9)) renom = 0.
     Call m1535(iabs(lb(i1)), iabs(lb(i2)), srt, x1535)
     x1440 = (3./4.)*sigma(srt, 2, 0, 1)
     If (((iabs(lb(i1))==2) .And. (iabs(lb(i2))==6)) .Or. ((iabs(lb(i2))==2) .And. (iabs(lb(i1))==6)) .Or. ((iabs(lb(i1))==1) .And. (iabs(lb(i2))==9)) .Or. ((iabs(lb(i2))==1) .And. (iabs(lb(i1))==9))) Then
        If (x1<=((signn+sdprod)/sig)) Goto 108
        If ((sigk+signn+sdprod)/sig>=x1) Goto 306
     End If
     If (lb(i1)*lb(i2)==18 .And. (iabs(lb(i1))==2 .Or. iabs(lb(i2))==2)) Then
        signd = sigma(srt, 1, 1, 0) + 0.5*sigma(srt, 1, 1, 1)
        sigdn = 0.25*signd*renom
        If (x1<=((signn+sdprod)/sig)) Goto 108
        If (x1>(signn+sigdn+x1440+x1535+sigk+sdprod)/sig) Return
        If (sigk/(sigk+sigdn+x1440+x1535)>ranart(nseed)) Goto 306
        If (ranart(nseed)<sigdn/(sigdn+x1440+x1535)) Then
           m12 = 3
           Goto 206
        Else
           If (ranart(nseed)<x1440/(x1440+x1535)) Then
              m12 = 37
           Else
              Return
           End If
           Goto 204
        End If
     End If
     If (lb(i1)*lb(i2)==6 .And. ((iabs(lb(i1))==1) .Or. (iabs(lb(i2))==1))) Then
        signd = sigma(srt, 1, 1, 0) + 0.5*sigma(srt, 1, 1, 1)
        sigdn = 0.25*signd*renom
        If (x1<=((signn+sdprod)/sig)) Goto 108
        If (x1>(signn+sigdn+x1440+x1535+sigk+sdprod)/sig) Return
        If (sigk/(sigk+sigdn+x1440+x1535)>ranart(nseed)) Goto 306
        If (ranart(nseed)<sigdn/(sigdn+x1440+x1535)) Then
           m12 = 6
           Goto 206
        Else
           If (ranart(nseed)<x1440/(x1440+x1535)) Then
              m12 = 47
           Else
              Return
           End If
           Goto 204
        End If
     End If
     If (lb(i1)*lb(i2)==8 .And. (iabs(lb(i1))==1 .Or. iabs(lb(i2))==1)) Then
        signd = 1.5*sigma(srt, 1, 1, 1)
        sigdn = 0.25*signd*renom
        If (x1<=((signn+sdprod)/sig)) Goto 108
        If (x1>(signn+sigdn+x1440+x1535+sigk+sdprod)/sig) Return
        If (sigk/(sigk+sigdn+x1440+x1535)>ranart(nseed)) Goto 306
        If (ranart(nseed)<sigdn/(sigdn+x1440+x1535)) Then
           m12 = 4
           Goto 206
        Else
           If (ranart(nseed)<x1440/(x1440+x1535)) Then
              m12 = 39
           Else
              m12 = 40
           End If
           Goto 204
        End If
     End If
     If (lb(i1)*lb(i2)==14 .And. (iabs(lb(i1))==2 .Or. iabs(lb(i2))==2)) Then
        signd = 1.5*sigma(srt, 1, 1, 1)
        sigdn = 0.25*signd*renom
        If (x1<=((signn+sdprod)/sig)) Goto 108
        If (x1>(signn+sigdn+x1440+x1535+sigk+sdprod)/sig) Return
        If (sigk/(sigk+sigdn+x1440+x1535)>ranart(nseed)) Goto 306
        If (ranart(nseed)<sigdn/(sigdn+x1440+x1535)) Then
           m12 = 5
           Goto 206
        Else
           If (ranart(nseed)<x1440/(x1440+x1535)) Then
              m12 = 48
           Else
              m12 = 49
           End If
           Goto 204
        End If
     End If
     If (lb(i1)*lb(i2)==16 .And. (iabs(lb(i1))==2 .Or. iabs(lb(i2))==2)) Then
        signd = 0.5*sigma(srt, 1, 1, 1) + 0.25*sigma(srt, 1, 1, 0)
        sigdn = 0.5*signd*renom
        If (x1<=((signn+sdprod)/sig)) Goto 108
        If (x1>(signn+sigdn+2.*x1440+2.*x1535+sigk+sdprod)/sig) Return
        If (sigk/(sigk+sigdn+2*x1440+2*x1535)>ranart(nseed)) Goto 306
        If (ranart(nseed)<sigdn/(sigdn+2.*x1440+2.*x1535)) Then
           m12 = 1
           Goto 206
        Else
           If (ranart(nseed)<x1440/(x1440+x1535)) Then
              m12 = 41
              If (ranart(nseed)<=0.5) m12 = 43
           Else
              m12 = 42
              If (ranart(nseed)<=0.5) m12 = 44
           End If
           Goto 204
        End If
     End If
     If (lb(i1)*lb(i2)==7) Then
        signd = 0.5*sigma(srt, 1, 1, 1) + 0.25*sigma(srt, 1, 1, 0)
        sigdn = 0.5*signd*renom
        If (x1<=((signn+sdprod)/sig)) Goto 108
        If (x1>(signn+sigdn+2.*x1440+2.*x1535+sigk+sdprod)/sig) Return
        If (sigk/(sigk+sigdn+2*x1440+2*x1535)>ranart(nseed)) Goto 306
        If (ranart(nseed)<sigdn/(sigdn+2.*x1440+2.*x1535)) Then
           m12 = 2
           Goto 206
        Else
           If (ranart(nseed)<x1440/(x1440+x1535)) Then
              m12 = 50
              If (ranart(nseed)<=0.5) m12 = 51
           Else
              m12 = 52
              If (ranart(nseed)<=0.5) m12 = 53
           End If
           Goto 204
        End If
     End If
     If (lb(i1)*lb(i2)==10 .And. (iabs(lb(i1))==1 .Or. iabs(lb(i2))==1)) Then
        signd = (3./4.)*sigma(srt, 2, 0, 1)
        sigdn = signd*renomn
        If (x1<=((signn+sdprod)/sig)) Goto 108
        If (x1>(signn+sigdn+x1535+sigk+sdprod)/sig) Return
        If (sigk/(sigk+sigdn+x1535)>ranart(nseed)) Goto 306
        If (ranart(nseed)<sigdn/(sigdn+x1535)) Then
           m12 = 7
           Goto 206
        Else
           m12 = 54
           If (ranart(nseed)<=0.5) m12 = 55
        End If
        Goto 204
     End If
     If (lb(i1)*lb(i2)==22 .And. (iabs(lb(i1))==2 .Or. iabs(lb(i2))==2)) Then
        signd = (3./4.)*sigma(srt, 2, 0, 1)
        sigdn = signd*renomn
        If (x1<=((signn+sdprod)/sig)) Goto 108
        If (x1>(signn+sigdn+x1535+sigk+sdprod)/sig) Return
        If (sigk/(sigk+sigdn+x1535)>ranart(nseed)) Goto 306
        If (ranart(nseed)<sigdn/(sigdn+x1535)) Then
           m12 = 8
           Goto 206
        Else
           m12 = 56
           If (ranart(nseed)<=0.5) m12 = 57
        End If
        Goto 204
     End If
     If ((iabs(lb(i1))==12) .Or. (iabs(lb(i1))==13) .Or. (iabs(lb(i2))==12) .Or. (iabs(lb(i2))==13)) Then
        signd = x1535
        sigdn = signd*renom1
        If (x1<=((signn+sdprod)/sig)) Goto 108
        If (x1>(signn+sigdn+sigk+sdprod)/sig) Return
        If (sigk/(sigk+sigdn)>ranart(nseed)) Goto 306
        If (lb(i1)*lb(i2)==24) m12 = 10
        If (lb(i1)*lb(i2)==12) m12 = 12
        If (lb(i1)*lb(i2)==26) m12 = 11
        If (lb(i1)*lb(i2)==13) m12 = 9
        Goto 206
     End If
204  Continue
     dmax = srt - avmass - 0.005
     dmin = 1.078
     If ((m12==37) .Or. (m12==39) .Or. (m12==41) .Or. (m12==43) .Or. (m12==46) .Or. (m12==48) .Or. (m12==50) .Or. (m12==51)) Then
        If (dmax<1.44) Then
           fm = fns(dmax, srt, 0.)
        Else
           xdmass = 1.44
           fm = fns(xdmass, srt, 1.)
        End If
        If (fm==0.) fm = 1.E-09
        ntry2 = 0
11      dm = ranart(nseed)*(dmax-dmin) + dmin
        ntry2 = ntry2 + 1
        If ((ranart(nseed)>fns(dm,srt,1.)/fm) .And. (ntry2<=10)) Goto 11
        If (dm>2.14) Goto 11
        Goto 13
     Else
        If (dmax<1.535) Then
           fm = fd5(dmax, srt, 0.)
        Else
           xdmass = 1.535
           fm = fd5(xdmass, srt, 1.)
        End If
        If (fm==0.) fm = 1.E-09
        ntry1 = 0
12      dm = ranart(nseed)*(dmax-dmin) + dmin
        ntry1 = ntry1 + 1
        If ((ranart(nseed)>fd5(dm,srt,1.)/fm) .And. (ntry1<=10)) Goto 12
        If (dm>1.84) Goto 12
     End If
13   Continue
     prf = 0.
     pf2 = ((srt**2-dm**2+avmass**2)/(2.*srt))**2 - avmass**2
     If (pf2>0.) prf = sqrt(pf2)
     If (m12==37) Then
        If (iabs(lb(i1))==9) Then
           lb(i1) = 1
           e(i1) = amp
           lb(i2) = 11
           e(i2) = dm
        Else
           lb(i2) = 1
           e(i2) = amp
           lb(i1) = 11
           e(i1) = dm
        End If
        Goto 207
     End If
     If (m12==38) Then
        If (iabs(lb(i1))==9) Then
           lb(i1) = 1
           e(i1) = amp
           lb(i2) = 13
           e(i2) = dm
        Else
           lb(i2) = 1
           e(i2) = amp
           lb(i1) = 13
           e(i1) = dm
        End If
        Goto 207
     End If
     If (m12==39) Then
        If (iabs(lb(i1))==8) Then
           lb(i1) = 1
           e(i1) = amp
           lb(i2) = 11
           e(i2) = dm
        Else
           lb(i2) = 1
           e(i2) = amp
           lb(i1) = 11
           e(i1) = dm
        End If
        Goto 207
     End If
     If (m12==40) Then
        If (iabs(lb(i1))==8) Then
           lb(i1) = 1
           e(i1) = amp
           lb(i2) = 13
           e(i2) = dm
        Else
           lb(i2) = 1
           e(i2) = amp
           lb(i1) = 13
           e(i1) = dm
        End If
        Goto 207
     End If
     If (m12==41) Then
        If (iabs(lb(i1))==8) Then
           lb(i1) = 2
           e(i1) = amn
           lb(i2) = 11
           e(i2) = dm
        Else
           lb(i2) = 2
           e(i2) = amn
           lb(i1) = 11
           e(i1) = dm
        End If
        Goto 207
     End If
     If (m12==42) Then
        If (iabs(lb(i1))==8) Then
           lb(i1) = 2
           e(i1) = amn
           lb(i2) = 13
           e(i2) = dm
        Else
           lb(i2) = 2
           e(i2) = amn
           lb(i1) = 13
           e(i1) = dm
        End If
        Goto 207
     End If
     If (m12==43) Then
        If (iabs(lb(i1))==8) Then
           lb(i1) = 1
           e(i1) = amp
           lb(i2) = 10
           e(i2) = dm
        Else
           lb(i2) = 1
           e(i2) = amp
           lb(i1) = 10
           e(i1) = dm
        End If
        Goto 207
     End If
     If (m12==44) Then
        If (iabs(lb(i1))==8) Then
           lb(i1) = 1
           e(i1) = amp
           lb(i2) = 12
           e(i2) = dm
        Else
           lb(i2) = 1
           e(i2) = amp
           lb(i1) = 12
           e(i1) = dm
        End If
        Goto 207
     End If
     If (m12==46) Then
        If (iabs(lb(i1))==6) Then
           lb(i1) = 2
           e(i1) = amn
           lb(i2) = 10
           e(i2) = dm
        Else
           lb(i2) = 2
           e(i2) = amn
           lb(i1) = 10
           e(i1) = dm
        End If
        Goto 207
     End If
     If (m12==47) Then
        If (iabs(lb(i1))==6) Then
           lb(i1) = 2
           e(i1) = amn
           lb(i2) = 12
           e(i2) = dm
        Else
           lb(i2) = 2
           e(i2) = amn
           lb(i1) = 12
           e(i1) = dm
        End If
        Goto 207
     End If
     If (m12==48) Then
        If (iabs(lb(i1))==7) Then
           lb(i1) = 2
           e(i1) = amn
           lb(i2) = 11
           e(i2) = dm
        Else
           lb(i2) = 2
           e(i2) = amn
           lb(i1) = 11
           e(i1) = dm
        End If
        Goto 207
     End If
     If (m12==49) Then
        If (iabs(lb(i1))==7) Then
           lb(i1) = 2
           e(i1) = amn
           lb(i2) = 12
           e(i2) = dm
        Else
           lb(i2) = 2
           e(i2) = amn
           lb(i1) = 12
           e(i1) = dm
        End If
        Goto 207
     End If
     If (m12==50) Then
        If (iabs(lb(i1))==7) Then
           lb(i1) = 1
           e(i1) = amp
           lb(i2) = 10
           e(i2) = dm
        Else
           lb(i2) = 1
           e(i2) = amp
           lb(i1) = 10
           e(i1) = dm
        End If
        Goto 207
     End If
     If (m12==51) Then
        If (iabs(lb(i1))==7) Then
           lb(i1) = 2
           e(i1) = amn
           lb(i2) = 11
           e(i2) = dm
        Else
           lb(i2) = 2
           e(i2) = amn
           lb(i1) = 11
           e(i1) = dm
        End If
        Goto 207
     End If
     If (m12==52) Then
        If (iabs(lb(i1))==7) Then
           lb(i1) = 1
           e(i1) = amp
           lb(i2) = 12
           e(i2) = dm
        Else
           lb(i2) = 1
           e(i2) = amp
           lb(i1) = 12
           e(i1) = dm
        End If
        Goto 207
     End If
     If (m12==53) Then
        If (iabs(lb(i1))==7) Then
           lb(i1) = 2
           e(i1) = amn
           lb(i2) = 13
           e(i2) = dm
        Else
           lb(i2) = 2
           e(i2) = amn
           lb(i1) = 13
           e(i1) = dm
        End If
        Goto 207
     End If
     If (m12==54) Then
        If (iabs(lb(i1))==10) Then
           lb(i1) = 2
           e(i1) = amn
           lb(i2) = 13
           e(i2) = dm
        Else
           lb(i2) = 2
           e(i2) = amn
           lb(i1) = 13
           e(i1) = dm
        End If
        Goto 207
     End If
     If (m12==55) Then
        If (iabs(lb(i1))==10) Then
           lb(i1) = 1
           e(i1) = amp
           lb(i2) = 12
           e(i2) = dm
        Else
           lb(i2) = 1
           e(i2) = amp
           lb(i1) = 12
           e(i1) = dm
        End If
        Goto 207
     End If
     If (m12==56) Then
        If (iabs(lb(i1))==11) Then
           lb(i1) = 2
           e(i1) = amn
           lb(i2) = 13
           e(i2) = dm
        Else
           lb(i2) = 2
           e(i2) = amn
           lb(i1) = 13
           e(i1) = dm
        End If
        Goto 207
     End If
     If (m12==57) Then
        If (iabs(lb(i1))==11) Then
           lb(i1) = 1
           e(i1) = amp
           lb(i2) = 12
           e(i2) = dm
        Else
           lb(i2) = 1
           e(i2) = amp
           lb(i1) = 12
           e(i1) = dm
        End If
     End If
     Goto 207
206  If (m12==1) Then
        If (iabs(lb(i1))==8) Then
           lb(i2) = 2
           lb(i1) = 1
           e(i1) = amp
        Else
           lb(i1) = 2
           lb(i2) = 1
           e(i2) = amp
        End If
        Goto 207
     End If
     If (m12==2) Then
        If (iabs(lb(i1))==7) Then
           lb(i2) = 1
           lb(i1) = 2
           e(i1) = amn
        Else
           lb(i1) = 1
           lb(i2) = 2
           e(i2) = amn
        End If
        Goto 207
     End If
     If (m12==3) Then
        lb(i1) = 1
        lb(i2) = 1
        e(i1) = amp
        e(i2) = amp
        Goto 207
     End If
     If (m12==4) Then
        lb(i1) = 1
        lb(i2) = 1
        e(i1) = amp
        e(i2) = amp
        Goto 207
     End If
     If (m12==5) Then
        lb(i1) = 2
        lb(i2) = 2
        e(i1) = amn
        e(i2) = amn
        Goto 207
     End If
     If (m12==6) Then
        lb(i1) = 2
        lb(i2) = 2
        e(i1) = amn
        e(i2) = amn
        Goto 207
     End If
     If (m12==7) Then
        If (iabs(lb(i1))==1) Then
           lb(i1) = 1
           lb(i2) = 2
           e(i1) = amp
           e(i2) = amn
        Else
           lb(i1) = 2
           lb(i2) = 1
           e(i1) = amn
           e(i2) = amp
        End If
        Goto 207
     End If
     If (m12==8) Then
        If (iabs(lb(i1))==2) Then
           lb(i1) = 2
           lb(i2) = 1
           e(i1) = amn
           e(i2) = amp
        Else
           lb(i1) = 1
           lb(i2) = 2
           e(i1) = amp
           e(i2) = amn
        End If
        Goto 207
     End If
     If (m12==9) Then
        lb(i1) = 1
        lb(i2) = 1
        e(i1) = amp
        e(i2) = amp
        Goto 207
     End If
     If (m12==12) Then
        lb(i1) = 2
        lb(i2) = 1
        e(i1) = amn
        e(i2) = amp
        Goto 207
     End If
     If (m12==11) Then
        lb(i1) = 2
        lb(i2) = 1
        e(i1) = amn
        e(i2) = amp
        Goto 207
     End If
     If (m12==12) Then
        lb(i1) = 1
        lb(i2) = 2
        e(i1) = amp
        e(i2) = amn
     End If
207  pr = prf
     c1 = 1.0 - 2.0*ranart(nseed)
     If (srt<=2.14) c1 = 1.0 - 2.0*ranart(nseed)
     If (srt>2.14 .And. srt<=2.4) c1 = ang(srt, iseed)
     If (srt>2.4) Then
        xptr = 0.33*pr
        cc1 = ptr(xptr, iseed)
        scheck = pr**2 - cc1**2
        If (scheck<0) Then
           Write (99, *) 'scheck4: ', scheck
           scheck = 0.
        End If
        c1 = sqrt(scheck)/pr
     End If
     t1 = 2.0*pi*ranart(nseed)
     iblock = 3
  End If
  If (ianti==1 .And. lb(i1)>=1 .And. lb(i2)>=1) Then
     lb(i1) = -lb(i1)
     lb(i2) = -lb(i2)
  End If
107 If (px==0.0 .And. py==0.0) Then
     t2 = 0.0
  Else
     t2 = atan2(py, px)
  End If
  scheck = 1.0 - c1**2
  If (scheck<0) Then
     Write (99, *) 'scheck5: ', scheck
     scheck = 0.
  End If
  s1 = sqrt(scheck)
  scheck = 1.0 - c2**2
  If (scheck<0) Then
     Write (99, *) 'scheck6: ', scheck
     scheck = 0.
  End If
  s2 = sqrt(scheck)
  ct1 = cos(t1)
  st1 = sin(t1)
  ct2 = cos(t2)
  st2 = sin(t2)
  pz = pr*(c1*c2-s1*s2*ct1)
  ss = c2*s1*ct1 + s2*c1
  px = pr*(ss*ct2-s1*st1*st2)
  py = pr*(ss*st2+s1*st1*ct2)
  Return
306 Continue
  If (xsk5/sigk>ranart(nseed)) Then
     pz1 = p(3, i1)
     pz2 = p(3, i2)
     lb(i1) = 1 + int(2*ranart(nseed))
     lb(i2) = 1 + int(2*ranart(nseed))
     nnn = nnn + 1
     lpion(nnn, irun) = 29
     epion(nnn, irun) = aphi
     iblock = 222
     Goto 208
  End If
  iblock = 11
  If (ianti==1) iblock = -11
  pz1 = p(3, i1)
  pz2 = p(3, i2)
  nnn = nnn + 1
  lpion(nnn, irun) = 23
  epion(nnn, irun) = aka
  If (srt<=2.63) Then
     ic = 1
     lb(i1) = 1 + int(2*ranart(nseed))
     lb(i2) = 14
     Goto 208
  End If
  If (srt<=2.74 .And. srt>2.63) Then
     If (xsk1/(xsk1+xsk2)>ranart(nseed)) Then
        ic = 1
        lb(i1) = 1 + int(2*ranart(nseed))
        lb(i2) = 14
     Else
        lb(i1) = 1 + int(2*ranart(nseed))
        lb(i2) = 15 + int(3*ranart(nseed))
        ic = 2
     End If
     Goto 208
  End If
  If (srt<=2.77 .And. srt>2.74) Then
     If (xsk1/(xsk1+xsk2+xsk3)>ranart(nseed)) Then
        ic = 1
        lb(i1) = 1 + int(2*ranart(nseed))
        lb(i2) = 14
        Goto 208
     Else
        If (xsk2/(xsk2+xsk3)>ranart(nseed)) Then
           ic = 2
           lb(i1) = 1 + int(2*ranart(nseed))
           lb(i2) = 15 + int(3*ranart(nseed))
        Else
           ic = 3
           lb(i1) = 6 + int(4*ranart(nseed))
           lb(i2) = 14
        End If
        Goto 208
     End If
  End If
  If (srt>2.77) Then
     If (xsk1/(xsk1+xsk2+xsk3+xsk4)>ranart(nseed)) Then
        ic = 1
        lb(i1) = 1 + int(2*ranart(nseed))
        lb(i2) = 14
        Goto 208
     Else
        If (xsk3/(xsk2+xsk3+xsk4)>ranart(nseed)) Then
           ic = 3
           lb(i1) = 6 + int(4*ranart(nseed))
           lb(i2) = 14
           Goto 208
        Else
           If (xsk2/(xsk2+xsk4)>ranart(nseed)) Then
              lb(i1) = 1 + int(2*ranart(nseed))
              lb(i2) = 15 + int(3*ranart(nseed))
              ic = 2
           Else
              ic = 4
              lb(i1) = 6 + int(4*ranart(nseed))
              lb(i2) = 15 + int(3*ranart(nseed))
           End If
           Goto 208
        End If
     End If
  End If
208 Continue
  If (ianti==1 .And. lb(i1)>=1 .And. lb(i2)>=1) Then
     lb(i1) = -lb(i1)
     lb(i2) = -lb(i2)
     If (lpion(nnn,irun)==23) lpion(nnn, irun) = 21
  End If
  lbi1 = lb(i1)
  lbi2 = lb(i2)
  ntry1 = 0
128 Call bbkaon(ic, srt, px3, py3, pz3, dm3, px4, py4, pz4, dm4, ppx, ppy, ppz, icou1)
  ntry1 = ntry1 + 1
  If ((icou1<0) .And. (ntry1<=20)) Goto 128
  Call rotate(px, py, pz, px3, py3, pz3)
  Call rotate(px, py, pz, px4, py4, pz4)
  Call rotate(px, py, pz, ppx, ppy, ppz)
  e1cm = sqrt(dm3**2+px3**2+py3**2+pz3**2)
  p1beta = px3*betax + py3*betay + pz3*betaz
  transf = gamma*(gamma*p1beta/(gamma+1)+e1cm)
  pt1i1 = betax*transf + px3
  pt2i1 = betay*transf + py3
  pt3i1 = betaz*transf + pz3
  eti1 = dm3
  e2cm = sqrt(dm4**2+px4**2+py4**2+pz4**2)
  p2beta = px4*betax + py4*betay + pz4*betaz
  transf = gamma*(gamma*p2beta/(gamma+1.)+e2cm)
  pt1i2 = betax*transf + px4
  pt2i2 = betay*transf + py4
  pt3i2 = betaz*transf + pz4
  eti2 = dm4
  epcm = sqrt(aka**2+ppx**2+ppy**2+ppz**2)
  ppbeta = ppx*betax + ppy*betay + ppz*betaz
  transf = gamma*(gamma*ppbeta/(gamma+1.)+epcm)
  ppion(1, nnn, irun) = betax*transf + ppx
  ppion(2, nnn, irun) = betay*transf + ppy
  ppion(3, nnn, irun) = betaz*transf + ppz
  dppion(nnn, irun) = dpertp(i1)*dpertp(i2)
  rpion(1, nnn, irun) = r(1, i1)
  rpion(2, nnn, irun) = r(2, i1)
  rpion(3, nnn, irun) = r(3, i1)
  p(1, i1) = pt1i1
  p(2, i1) = pt2i1
  p(3, i1) = pt3i1
  e(i1) = eti1
  lb(i1) = lbi1
  p(1, i2) = pt1i2
  p(2, i2) = pt2i2
  p(3, i2) = pt3i2
  e(i2) = eti2
  lb(i2) = lbi2
  px1 = p(1, i1)
  py1 = p(2, i1)
  pz1 = p(3, i1)
  em1 = e(i1)
  id(i1) = 2
  id(i2) = 2
  id1 = id(i1)
  If (lpion(nnn,irun)/=29) iblock = 11
  lb1 = lb(i1)
  lb2 = lb(i2)
  am1 = em1
  am2 = em2
  e1 = sqrt(em1**2+px1**2+py1**2+pz1**2)
  Return
108 Continue
  If (idpert==1 .And. ipert1==1 .And. npertd>=1) Then
     ndloop = npertd
  Else If (idpert==2 .And. npertd>=1) Then
     ndloop = npertd + 1
  Else
     ndloop = 1
  End If
  dprob1 = sdprod/sig/float(npertd)
  Do idloop = 1, ndloop
     Call bbdangle(pxd, pyd, pzd, nt, ipert1, ianti, idloop, pfinal, dprob1, lbm)
     Call rotate(px, py, pz, pxd, pyd, pzd)
     xmass = xmd
     e1dcm = sqrt(xmass**2+pxd**2+pyd**2+pzd**2)
     p1dbeta = pxd*betax + pyd*betay + pzd*betaz
     transf = gamma*(gamma*p1dbeta/(gamma+1.)+e1dcm)
     pxi1 = betax*transf + pxd
     pyi1 = betay*transf + pyd
     pzi1 = betaz*transf + pzd
     If (ianti==0) Then
        lbd = 42
     Else
        lbd = -42
     End If
     If (idpert==1 .And. ipert1==1 .And. npertd>=1) Then
        nnn = nnn + 1
        ppion(1, nnn, irun) = pxi1
        ppion(2, nnn, irun) = pyi1
        ppion(3, nnn, irun) = pzi1
        epion(nnn, irun) = xmd
        lpion(nnn, irun) = lbd
        rpion(1, nnn, irun) = r(1, i1)
        rpion(2, nnn, irun) = r(2, i1)
        rpion(3, nnn, irun) = r(3, i1)
        dppion(nnn, irun) = sdprod/sig/float(npertd)
     Else If (idpert==2 .And. idloop<=npertd) Then
        ppd(1, idloop) = pxi1
        ppd(2, idloop) = pyi1
        ppd(3, idloop) = pzi1
        lbpd(idloop) = lbd
     Else
        e(i1) = xmm
        e2picm = sqrt(xmm**2+pxd**2+pyd**2+pzd**2)
        p2pibeta = -pxd*betax - pyd*betay - pzd*betaz
        transf = gamma*(gamma*p2pibeta/(gamma+1.)+e2picm)
        pxi2 = betax*transf - pxd
        pyi2 = betay*transf - pyd
        pzi2 = betaz*transf - pzd
        p(1, i1) = pxi2
        p(2, i1) = pyi2
        p(3, i1) = pzi2
        lb(i1) = lbm
        px1 = p(1, i1)
        py1 = p(2, i1)
        pz1 = p(3, i1)
        em1 = e(i1)
        id(i1) = 2
        id1 = id(i1)
        e1 = sqrt(em1**2+px1**2+py1**2+pz1**2)
        lb1 = lb(i1)
        p(1, i2) = pxi1
        p(2, i2) = pyi1
        p(3, i2) = pzi1
        lb(i2) = lbd
        lb2 = lb(i2)
        e(i2) = xmd
        eti2 = e(i2)
        id(i2) = 2
        If (idpert==2 .And. idloop==ndloop) Then
           Do ipertd = 1, npertd
              nnn = nnn + 1
              ppion(1, nnn, irun) = ppd(1, ipertd)
              ppion(2, nnn, irun) = ppd(2, ipertd)
              ppion(3, nnn, irun) = ppd(3, ipertd)
              epion(nnn, irun) = xmd
              lpion(nnn, irun) = lbpd(ipertd)
              rpion(1, nnn, irun) = r(1, i1)
              rpion(2, nnn, irun) = r(2, i1)
              rpion(3, nnn, irun) = r(3, i1)
              dppion(nnn, irun) = 1./float(npertd)
           End Do
        End If
     End If
  End Do
  iblock = 501
  Return
End Subroutine crnd

max_stars_repo_name: xiaohaijin/AMPT

max_issues_repo_name: xiaohaijin/AMPT

max_forks_repo_name: xiaohaijin/AMPT


text:
! { dg-do compile }
! { dg-options "-fwhole-file" }
!
! Tests the fix PR40011 comment 17 in which the explicit interface was
! being ignored and the missing argument was not correctly handled, which
! led to an ICE.
!
! Contributed by Dominique d'Humieres  <dominiq@lps.ens.fr
!
          Implicit None
          call sub(1,2)
          call sub(1,2,3)

          contains

          subroutine sub(i,j,k)
          Implicit None
          Integer, Intent( In )           :: i
          Integer, Intent( In )           :: j
          Integer, Intent( In ), Optional :: k
          intrinsic present
          write(*,*)' 3 presence flag ',present(k)
          write(*,*)' 1st arg ',i
          write(*,*)' 2nd arg ',j
          if (present(k)) then
            write(*,*)' 3rd arg ',k
          else
            write(*,*)' 3rd arg is absent'
          endif
          return
          end subroutine

          end

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
program AOC

  integer :: ios, u = 42, result, i, j, k, m, m_t, m_b, m_l, m_r
  integer, dimension(100,100) :: map
  character(len = 100) :: line

  open(unit=u, file='day_9.txt', iostat=ios)

  do i = 1, 100
     read(u, *, iostat=ios) line
     do j = 1, 100
        map(j, i) = ichar(line(j:j)) - 48
     end do
  end do

  do i = 1, 100
     do j = 1, 100
        m = map(j, i)
        if (j == 1 .or. m < map(j - 1, i)) then
           if (j == 100 .or. m < map(j + 1, i)) then
              if (i == 1 .or. m < map(j, i - 1)) then
                 if (i == 100 .or. m < map(j, i + 1)) then
                    result = result + 1 + m
                 end if
              end if
           end if
        end if
     end do
  end do

  print*, "Answer = ", result

  close(u)
end program AOC

max_stars_repo_name: jbouwman/musical-potato

max_issues_repo_name: jbouwman/musical-potato

max_forks_repo_name: jbouwman/musical-potato


!! Print leading zeros for real, integer in Fortran
use, intrinsic :: iso_fortran_env, only: stdout=>output_unit, dp=>real64

implicit none

real(dp):: x
integer :: i

x = 1234.56_dp
i = 123456

write(stdout,'(a5,f12.7)') 'hello',x
write(stdout,'(i0.12)')  i

print *,'there should be leading zeros and no blanks for the real and integer above.'

end program


text:
*DECK ACOSH
      FUNCTION ACOSH (X)
C***BEGIN PROLOGUE  ACOSH
C***PURPOSE  Compute the arc hyperbolic cosine.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C4C
C***TYPE      SINGLE PRECISION (ACOSH-S, DACOSH-D, CACOSH-C)
C***KEYWORDS  ACOSH, ARC HYPERBOLIC COSINE, ELEMENTARY FUNCTIONS, FNLIB,
C             INVERSE HYPERBOLIC COSINE
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C ACOSH(X) computes the arc hyperbolic cosine of X.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  R1MACH, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770401  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   900326  Removed duplicate information from DESCRIPTION section.
C           (WRB)
C***END PROLOGUE  ACOSH
      SAVE ALN2,XMAX
      DATA ALN2 / 0.6931471805 5994530942E0/
      DATA XMAX /0./
C***FIRST EXECUTABLE STATEMENT  ACOSH
      IF (XMAX.EQ.0.) XMAX = 1.0/SQRT(R1MACH(3))
C
      IF (X .LT. 1.0) CALL XERMSG ('SLATEC', 'ACOSH', 'X LESS THAN 1',
     +   1, 2)
C
      IF (X.LT.XMAX) ACOSH = LOG (X + SQRT(X*X-1.0))
      IF (X.GE.XMAX) ACOSH = ALN2 + LOG(X)
C
      RETURN
      END

ext: f


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2021 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine gencycsymelemcfd(cs,islav,nslav,imast,nmast,inomat,
     &  nk,co,ne,ipkon,lakon,kon,nkon,ielmat,mi,vold,ielslav,ielmast,
     &  inoslav,inomast,iponoel,inoel)
!
!     creates an additional layer of elements on the slave and the
!     master surface of a cyclic symmetric structure
!
      implicit none
!
      character*8 lakon(*)
!
      integer imast(*),i,j,k,nmast,islav(*),nslav,
     &  nk,ne,ipkon(*),kon(*),nkon,
     &  mi(*),index,nodes,nodem,nope,
     &  ielmat(mi(3),*),indexe,inoslav(*),inomast(*),ielslav(*),
     &  ielmast(*),iponoel(*),inoel(3,*),nodeslav,nodemast,ielemslav,
     &  ielemmast,inomat(*)
!
      real*8 cs(17,*),phibasis,phi,xa(3),xn(3),dd,cphi,sphi,dphi,
     &  c(3,3),xp(3),xq(3),al,co(3,*),vold(0:mi(2),*)
!
!     mcs=1 is assumed
!     phi is the angle of the cyclic basis sector
!
      phibasis=8.d0*datan(1.d0)/cs(1,1)
!
!     xn is an normed vector along the axis
!     xa is a point on the axis
!
      do i=1,3
         xa(i)=cs(5+i,1)
         xn(i)=cs(8+i,1)-xa(i)
      enddo
      dd=dsqrt(xn(1)*xn(1)+xn(2)*xn(2)+xn(3)*xn(3))
      do i=1,3
         xn(i)=xn(i)/dd
      enddo
!
!     new nodes on the master side
!
!     phi is the thickness of the layer in radians
!
      phi=phibasis
      cphi=dcos(phi)
      sphi=dsin(phi)
      dphi=1.d0-cphi
!
!     rotation matrix
!
      c(1,1)=cphi+dphi*xn(1)*xn(1)
      c(1,2)=     dphi*xn(1)*xn(2)-sphi*xn(3)
      c(1,3)=     dphi*xn(1)*xn(3)+sphi*xn(2)
      c(2,1)=     dphi*xn(2)*xn(1)+sphi*xn(3)
      c(2,2)=cphi+dphi*xn(2)*xn(2)
      c(2,3)=     dphi*xn(2)*xn(3)-sphi*xn(1)
      c(3,1)=     dphi*xn(3)*xn(1)-sphi*xn(2)
      c(3,2)=     dphi*xn(3)*xn(2)+sphi*xn(1)
      c(3,3)=cphi+dphi*xn(3)*xn(3)
!
!     new nodes and elements on the master side
!
      do i=1,nmast
!
!        corresponding slave node
!
         nodeslav=inoslav(imast(i))
!
!        loop over all elements to which nodeslav belongs
!
         index=iponoel(nodeslav)
         if(index.eq.0) cycle
!
         do
            ielemslav=inoel(1,index)
            if(ielmast(ielemslav).eq.0) then
               read(lakon(ielemslav)(4:4),'(i1)') nope
               indexe=ipkon(ielemslav)
               do k=1,nope
                  nodes=kon(indexe+k)
                  if(inomast(nodes).eq.0) then
!
!                    generate new node on the master side
!
                     nk=nk+1
                     inomast(nodes)=nk
!
!                    calculating the coordinates of the node
!
                     do j=1,3
                        xp(j)=co(j,nodes)
                     enddo
!
                     al=(xp(1)-xa(1))*xn(1)+
     &                    (xp(2)-xa(2))*xn(2)+
     &                    (xp(3)-xa(3))*xn(3)
!
!                    xq is the orthogonal projection of xp on the axis
!
                     do j=1,3
                        xq(j)=xa(j)+al*xn(j)
                     enddo
!
                     do j=1,3
                        co(j,nk)=xq(j)+c(j,1)*(xp(1)-xq(1))+
     &                       c(j,2)*(xp(2)-xq(2))+
     &                       c(j,3)*(xp(3)-xq(3))
                     enddo
!
!                    initial conditions (global rectangular system)
!
                     do j=0,4
                        vold(j,nk)=vold(j,nodes)
                     enddo
!
!                    assign material to node
!
                     inomat(nk)=inomat(nodes)
!
                  endif
               enddo
!
!              create a new element on the master side
!
               ne=ne+1
               ielmast(ielemslav)=ne
               ipkon(ne)=nkon
               lakon(ne)=lakon(ielemslav)
               ielmat(1,ne)=ielmat(1,ielemslav)
               do j=1,nope
                  kon(nkon+j)=inomast(kon(indexe+j))
               enddo
               nkon=nkon+nope
            endif
!
            index=inoel(3,index)
            if(index.eq.0) exit
         enddo
      enddo
!
!     new nodes on the slave side
!
      phi=-phibasis
      cphi=dcos(phi)
      sphi=dsin(phi)
      dphi=1.d0-cphi
!
!     rotation matrix
!
      c(1,1)=cphi+dphi*xn(1)*xn(1)
      c(1,2)=     dphi*xn(1)*xn(2)-sphi*xn(3)
      c(1,3)=     dphi*xn(1)*xn(3)+sphi*xn(2)
      c(2,1)=     dphi*xn(2)*xn(1)+sphi*xn(3)
      c(2,2)=cphi+dphi*xn(2)*xn(2)
      c(2,3)=     dphi*xn(2)*xn(3)-sphi*xn(1)
      c(3,1)=     dphi*xn(3)*xn(1)-sphi*xn(2)
      c(3,2)=     dphi*xn(3)*xn(2)+sphi*xn(1)
      c(3,3)=cphi+dphi*xn(3)*xn(3)
!
!     new nodes and elements on the slave side
!
      do i=1,nslav
!
!        corresponding slave node
!
         nodemast=inomast(islav(i))
!
!        loop over all elements to which nodeslav belongs
!
         index=iponoel(nodemast)
         if(index.eq.0) cycle
!
         do
            ielemmast=inoel(1,index)
            if(ielslav(ielemmast).eq.0) then
               read(lakon(ielemmast)(4:4),'(i1)') nope
               indexe=ipkon(ielemmast)
               do k=1,nope
                  nodem=kon(indexe+k)
                  if(inoslav(nodem).eq.0) then
!
!                    generate new node on the master side
!
                     nk=nk+1
                     inoslav(nodem)=nk
!
!                    calculating the coordinates of the node
!
                     do j=1,3
                        xp(j)=co(j,nodem)
                     enddo
!
                     al=(xp(1)-xa(1))*xn(1)+
     &                    (xp(2)-xa(2))*xn(2)+
     &                    (xp(3)-xa(3))*xn(3)
!
!                    xq is the orthogonal projection of xp on the axis
!
                     do j=1,3
                        xq(j)=xa(j)+al*xn(j)
                     enddo
!
                     do j=1,3
                        co(j,nk)=xq(j)+c(j,1)*(xp(1)-xq(1))+
     &                       c(j,2)*(xp(2)-xq(2))+
     &                       c(j,3)*(xp(3)-xq(3))
                     enddo
!
!                    initial conditions (global rectangular system)
!
                     do j=0,4
                        vold(j,nk)=vold(j,nodem)
                     enddo
!
!                    assign material to node
!
                     inomat(nk)=inomat(nodem)
!
                  endif
               enddo
!
!              create a new element on the slave side
!
               ne=ne+1
               ielslav(ielemmast)=ne
               ipkon(ne)=nkon
               lakon(ne)=lakon(ielemmast)
               ielmat(1,ne)=ielmat(1,ielemmast)
               do j=1,nope
                  kon(nkon+j)=inoslav(kon(indexe+j))
               enddo
               nkon=nkon+nope
            endif
!
            index=inoel(3,index)
            if(index.eq.0) exit
         enddo
      enddo
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
*
*     ------------------------------------------------------------------
*	B N D I N V
*     ------------------------------------------------------------------
*
        SUBROUTINE BNDINV(A,EL,N,DETERM,EPSIL,ITEST,NSIZE)
C
C       DOUBLE PRECISION MATRIX INVERSION SUBROUTINE
C       FROM "DLYTAP".
C
C*      DOUBLE PRECISION E,F
C*      DOUBLE PRECISION A,EL,D,DSQRT,C,S,DETERP
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        DIMENSION A(NSIZE,1),EL(NSIZE,1)
        IF(N.LT.2)GO TO 140
        ISL2=0
        K000FX=2
        IF(ISL2.EQ.0)INDSNL=2
        IF(ISL2.EQ.1)INDSNL=1
C       CALL SLITET(2,INDSNL)
C       CALL OVERFL(K000FX)
C       CALL DVCHK(K000FX)
C
C       SET EL = IDENTITY MATRIX
        DO 30 I=1,N
        DO 10 J=1,N
 10     EL(I,J)=0.0D0
 30     EL(I,I)=1.0D0
C
C       TRIANGULARIZE A, FORM EL
C
        N1=N-1
        M=2
        DO 50 J=1,N1
        DO 45 I=M,N
        IF(A(I,J).EQ.0.0D0)GO TO 45
        D=DSQRT(A(J,J)*A(J,J)+A(I,J)*A(I,J))
        C=A(J,J)/D
        S=A(I,J)/D
 38     DO 39 K=J,N
        D=C*A(J,K)+S*A(I,K)
        A(I,K)=C*A(I,K)-S*A(J,K)
        A(J,K)=D
 39     CONTINUE
        DO 40 K=1,N
        D=C*EL(J,K)+S*EL(I,K)
        EL(I,K)=C*EL(I,K)-S*EL(J,K)
        EL(J,K)=D
 40     CONTINUE
 45     CONTINUE
 50     M=M+1
C       CALL OVERFL(K000FX)
C       GO TO (140,51),K000FX
C
C       CALCULATE THE DETERMINANT
 51     DETERP=A(1,1)
        DO 52 I=2,N
 52     DETERP=DETERP*A(I,I)
        DETERM=DETERP
C       CALL OVERFL(K000FX)
C       GO TO (140,520,520),K000FX
C
C       IS MATRIX SINGULAR
 520    F=A(1,1)
        E=A(1,1)
        DO 58 I=2,N
        IF(DABS(F).LT.DABS(A(I,I)))F=A(I,I)
        IF(DABS(E).GT.DABS(A(I,I)))E=A(I,I)
 58     CONTINUE
        EPSILP=EPSIL
        IF(EPSILP.LE.0)EPSILP=1.0E-8
        RAT=E/F
        IF(ABS(RAT).LT.EPSILP)GO TO 130
C
C       INVERT TRIANGULAR MATRIX
        J=N
        DO 100 J1=1,N
C       CALL SLITE(2)
        I=J
        ISL2=1
        DO 90 I1=1,J
C       CALL SLITET(2,K000FX)
        IF(ISL2.EQ.0)K000FX=2
        IF(ISL2.EQ.1)K000FX=1
        IF(ISL2.EQ.1)ISL2=0
        GO TO (70,75),K000FX
 70     A(I,J)=1.0D0/A(I,I)
        GO TO 90
 75     KS=I+1
        D=0.0D0
        DO 80 K=KS,J
 80     D=D+A(I,K)*A(K,J)
        A(I,J)=-D/A(I,I)
 90     I=I-1
 100    J=J-1
C       CALL OVERFL(K000FX)
C       GO TO (140,103,103),K000FX

C103    CALL DVCHK(K000FX)
C       GO TO (140,105),K000FX
C
C       PREMULTIPLY EL BY INVERTED TRIANGULAR MATRIX
 105    M=1
        DO 120 I=1,N
        DO 118 J=1,N
        D=0.0D0
        DO 107 K=M,N
 107    D=D+A(I,K)*EL(K,J)
        EL(I,J)=D
 118    CONTINUE
 120    M=M+1
C       CALL OVERFL(K000FX)
C       GO TO (140,123,123),K000FX
C
C       RECOPY EL TO A
 123    DO 124 I=1,N
        DO 124 J=1,N
 124    A(I,J)=EL(I,J)
        ITEST=0
C126    IF(INDSNL.EQ.1)CALL SLITE(2)
 126    IF(INDSNL.EQ.1)ISL2=1
        RETURN
C
 130    ITEST=1
        GO TO 126
 140    ITEST=-1
        GO TO 126
        END


text:
        !COMPILER-GENERATED INTERFACE MODULE: Tue Oct 10 17:23:02 2017
        MODULE GEIDER_DISC__genmod
          INTERFACE
            SUBROUTINE GEIDER_DISC
            END SUBROUTINE GEIDER_DISC
          END INTERFACE
        END MODULE GEIDER_DISC__genmod

max_stars_repo_name: BingzhangChen/NPZDFeCONT

max_issues_repo_name: BingzhangChen/NPZDFeCONT

max_forks_repo_name: BingzhangChen/NPZDFeCONT


text:
SUBROUTINE DCSYMBDL(A,M,N,NN,EPS,Z,W,L,LI,LJ,IER)
!**********************************************************************
!  GAUSS METHOD FOR A SYMMETRIC BAND MATRIX WITH A SHORT WIDTH.       *
!  THIS ROUTINE IS FOR A VECTOR COMPUTER.                             *
!                                                                     *
!  PARAMETERS:                                                        *
!   ON ENTRY:                                                         *
!     A      THE ARRAY WITH DIMENSION (M+1)*N WHICH CONTAINS          *
!            THE LOWER BAND MATRIX IN ROW-WISE.                       *
!     M      THE HALF BAND WIDTH OF THE MATRIX A, EXCLUDING DIADONALS.*
!     N      THE ORDER OF THE MATRIX A.                               *
!     NN     THE ORDER OF WORKING ARRAYS L, LI, AND LJ.               *
!     EPS    THE TOLERANCE FOR PIVOTAL ELEMENTS.                      *
!   ON RETURN:                                                        *
!     A      THE LOWER TRIANGULAR MATRIX L AFTER DECOMPOSITION.       *
!     IER    ERROR CODE. IF IER=0, NORMAL RETURN.                     *
!   OTHERS: WORKING PARAMETERS.                                       *
!                                                                     *
!  COPYRIGHT:       FUMIKO NAGAHORI    1 SEP. 1991      VER. 1        *
!    modified: Kensuke Konishi 2018 in Paris
!**********************************************************************
    use parameters
    implicit none
    INTEGER:: N,M,MM,NN
    INTEGER:: L(NN),LI(NN),LJ(NN),IER
    double precision::EPS
    COMPLEX(dp):: A((M+1)*N),Z(M+1),W(M+1)
    INTEGER:: I,J,K,IJ,KK,NK,NKK,NKI
    COMPLEX(dp):: PIV

    IER = 0
    IJ = 0
    DO I=1,M
        IJ = IJ + 1
        LI(IJ) = I + 1
        LJ(IJ) = 2
        L(IJ) = M * (I-1) + 1
        DO J=I+1,M
            IJ = IJ + 1
            L(IJ) = L(IJ-1) + M + 1
            LI(IJ) = LI(IJ-1) + 1
            LJ(IJ) = LJ(IJ-1) + 1
        enddo
    enddo
    MM = (M+1) * M / 2
    DO K=1,N-M
        NK = (M+1) * (K-1) + 1
        IF (cdabs(A(NK+M)) < EPS) THEN
            WRITE(*,*) '(SUBR. SYMBDL) SINGULAR AT STEP = ', K
            IER = 1
            RETURN
        ENDIF
        PIV = 1 / A(NK+M)
        DO J=2,M+1
            Z(J) = -A(NK+M*J)
            W(J) = A(NK+M*J) * PIV
            A(NK+M*J) = W(J)
        enddo

        KK = NK + M + M
        !VORTION VEC
        A(KK+L(1:mm)) = A(KK+L(1:mm)) + W(LJ(1:mm)) * Z(LI(1:mm))
    enddo

    DO K=N-M+1,N-1
        NK = (M+1) * (K-1) + 1
        NKK = (M+1) * K - 1
        IF (cdabs(A(NK+M)) < EPS) THEN
            WRITE(*,*) '(SUBR. SYMBDL) SINGULAR AT STEP = ', K
            IER = 1
            RETURN
        ENDIF
        PIV = 1 / A(NK+M)
        DO  J=2,N-K+1
            Z(J) = -A(NK+M*J)
            W(J) = A(NK+M*J) * PIV
            A(NK+M*J) = W(J)
            NKI = NKK + M * (J-1)
            A(NKI+2:nki+j) = A(NKI+2:nki+j) + W(J) * Z(2:j)
        enddo
    enddo
    RETURN
END SUBROUTINE DCSYMBDL

SUBROUTINE DCSBDLV(A,B,M,N,NP,Z)
!**********************************************************************
!  GAUSS METHOD FOR A SYMMETRIC BAND MATRIX WITH A SHORT WIDTH.       *
!  THIS ROUTINE IS FOR A VECTOR COMPUTER.                             *
!                                                                     *
!  PARAMETERS:                                                        *
!   ON ENTRY:                                                         *
!     A      THE ARRAY WITH DIMENSION (M+1),N WHICH CONTAINS          *
!            THE LEFT HAND SIDE LOWER BAND MATRIX.                    *
!     M      THE HALF BAND WIDTH OF THE MATRIX A, EXCLUDING DIADONALS.*
!     N      THE ORDER OF THE MATRIX A.                               *
!     EPS    THE TOLERANCE FOR PIVOTAL ELEMENTS.                      *
!     B      RIGHT HAND SIDE VECTOR                                   *
!   ON RETURN:                                                        *
!     B      SOLUTION.                                                *
!     IER    ERROR CODE. IF IER=0, NORMAL RETURN.                     *
!   OTHERS: WORKING PARAMETERS.                                       *
!                                                                     *
!  COPYRIGHT:       FUMIKO NAGAHORI    1 SEP. 1991      VER. 1        *
! modified: Kensuke Konishi 2018 in Paris
!**********************************************************************
    use parameters
    implicit none
    INTEGER:: M,N
    integer,intent(in):: NP
    COMPLEX(dp):: A(M+1,N),B(N),Z(N)
    INTEGER:: MM,J,K,I1
    COMPLEX(dp):: SUM

    ! FORWARD SUBSTITUTION
    MM = M + 1
    IF (MM < 3) THEN
        Z(NP) = B(NP)
        DO J=NP+1,N
            Z(J) = B(J) - A(1,J) * Z(J-1)
        enddo
        B(N) = Z(N) / A(M+1,N)
    ELSE
        Z(NP) = B(NP)
        Z(NP+1) = B(NP+1) - A(MM-1,NP+1) * Z(NP)
        DO J=NP+2,N
            IF (J > NP-1+MM) THEN
                I1 = 1
            ELSE
                I1 = NP-1+MM - J + 1
            ENDIF
            SUM = 0
            DO K=I1,MM-1
                SUM = SUM + A(K,J) * Z(J-MM+K)
            enddo
            Z(J) = B(J) - SUM
        enddo

        Z(n-1:n) = Z(n-1:n) / A(M+1,n-1:n)

        B(N) = Z(N)
        B(N-1) = Z(N-1) - A(MM-1,N) * Z(N)
    ENDIF
    RETURN
END SUBROUTINE DCSBDLV

SUBROUTINE DCSBDLV0(A,B,M,N,Z)
!**********************************************************************
!  GAUSS METHOD FOR A SYMMETRIC BAND MATRIX WITH A SHORT WIDTH.       *
!  THIS ROUTINE IS FOR A VECTOR COMPUTER.                             *
!                                                                     *
!  PARAMETERS:                                                        *
!   ON ENTRY:                                                         *
!     A      THE ARRAY WITH DIMENSION (M+1),N WHICH CONTAINS          *
!            THE LEFT HAND SIDE LOWER BAND MATRIX.                    *
!     M      THE HALF BAND WIDTH OF THE MATRIX A, EXCLUDING DIADONALS.*
!     N      THE ORDER OF THE MATRIX A.                               *
!     EPS    THE TOLERANCE FOR PIVOTAL ELEMENTS.                      *
!     B      RIGHT HAND SIDE VECTOR                                   *
!   ON RETURN:                                                        *
!     B      SOLUTION.                                                *
!     IER    ERROR CODE. IF IER=0, NORMAL RETURN.                     *
!   OTHERS: WORKING PARAMETERS.                                       *
!                                                                     *
!  COPYRIGHT:       FUMIKO NAGAHORI    1 SEP. 1991      VER. 1        *
!**********************************************************************
    use parameters
    IMPLICIT NONE
    integer:: M,N
    complex(dp):: A(M+1,N),B(N),Z(N)
    integer:: MM,J,K,I1,J1
    complex(dp):: SUM

    ! FORWARD SUBSTITUTION
    MM = M + 1
    IF (MM < 3) THEN
        Z(1) = B(1)
        DO J=2,N
            Z(J) = B(J) - A(1,J) * Z(J-1)
        enddo
        DO J=1,N
            Z(J) = Z(J) / A(M+1,J)
        enddo
        B(N) = Z(N)
        DO J=1,N-1
            B(N-J) = Z(N-J) - A(1,N-J+1) * B(N-J+1)
        enddo
    ELSE
        Z(1) = B(1)
        Z(2) = B(2) - A(M-1,2) * Z(1)
        DO J=3,N
            IF (J > MM) THEN
                I1 = 1
            ELSE
                I1 = MM - J + 1
            ENDIF
            SUM = 0
            DO K=I1,MM-1
                SUM = SUM + A(K,J) * Z(J-MM+K)
            enddo
            Z(J) = B(J) - SUM
        enddo

        Z(1:n) = Z(1:n) / A(M+1,1:n)

        B(N) = Z(N)
        B(N-1) = Z(N-1) - A(MM-1,N) * Z(N)
        DO J=3,N
            J1 = N - J + 1
            I1 = 1
            IF (J < MM) I1 = MM - J + 1
            SUM = 0
            DO K=I1,MM-1
                SUM = SUM + A(K,MM-K+J1) * B(MM-K+J1)
            enddo
            B(J1) = Z(J1) - SUM
        enddo
    ENDIF
    RETURN
END

max_stars_repo_name: seismobassoon/dsmpy

max_issues_repo_name: seismobassoon/dsmpy

max_forks_repo_name: seismobassoon/dsmpy


text:
	SUBROUTINE YGTORB(ISTAT)
*
*  Module number:
*
*  Module name: YGTORB
*
*  Keyphrase:
*  ----------
*       Read shp header keywords for orbital and positional information
*
*  Description:
*  ------------
*       This routine reads keywords from the shp that give
*	positional and orbital information about HST and
*       places the values in a common block.
*
*  FORTRAN name: ygtorb.for
*
*  Keywords of accessed files and tables:
*  --------------------------------------
*  Name                         I/O     Description / Comments
*       <rootname>.SHH          I       FOS data file
*
*  Subroutines Called:
*  -------------------
*  CDBS:
*       ymsput
*  SDAS:
*       uhdgs*
*
*  History:
*  --------
*  Version      Date        Author          Description
*       1       May 91  S. Hulbert      Designed and coded
*     1.1       Jun 91  S. Hulbert      Fixed bug in reading attitiude keywords
*-------------------------------------------------------------------------------
        INTEGER ISTAT
C                                    --->ERROR status
C
C     UMSPUT DESTINATIONS -- CB, DAO, 4-SEP-87
C
      INTEGER STDOUT
      PARAMETER (STDOUT = 1)
      INTEGER STDERR
      PARAMETER (STDERR = 2)
C
C Common block containing input/output file descriptions
C
C       IDS - file id numbers
C       GCOUNT - group counts
C       NAXIS - naxis
C       NAXIS1 - first dimensions
C       NAXIS2 - second dimensions
C       FILL - Fill values
C
        INTEGER IDS(30),NAXIS(30),NAXIS1(30),NAXIS2(30),GCOUNT(30)
        REAL FILL(30)
        COMMON /IOCOM/IDS,GCOUNT,NAXIS,NAXIS1,NAXIS2,FILL

C
C Common block containing orbital/positional parameters
C
        DOUBLE PRECISION EPCHTIME,SDMEANAN,CIRVELOC,COSINCLI
        DOUBLE PRECISION ECCENTRY,ECCENTX2,ECBDX4D3,ESQDX5D2,ECBDX3
        DOUBLE PRECISION FDMEANAN,RCASCNRV,ARGPERIG,MEANANOM,RCARGPER
        DOUBLE PRECISION RASCASCN,SINEINCL,SEMILREC
        DOUBLE PRECISION PSANGLV3,RTASCNV1,DECLNV1
        COMMON /ORBPAR/ EPCHTIME,SDMEANAN,CIRVELOC,COSINCLI,
     $  ECCENTRY,ECCENTX2,ECBDX4D3,ESQDX5D2,ECBDX3,
     $  FDMEANAN,RCASCNRV,ARGPERIG,MEANANOM,RCARGPER,
     $  RASCASCN,SINEINCL,SEMILREC,
     $  PSANGLV3,RTASCNV1,DECLNV1
C
C Local variables
C
        INTEGER I, ID, ISTATS(20)
        CHARACTER*8 NAMES(20)
        CHARACTER*80 CONTXT
        DATA NAMES/'EPCHTIME','SDMEANAN','CIRVELOC','COSINCLI',
     $  'ECCENTRY','ECCENTX2','ECBDX4D3','ESQDX5D2','ECBDX3',
     $  'FDMEANAN','RCASCNRV','ARGPERIG','MEANANOM','RCARGPER',
     $  'RASCASCN','SINEINCL','SEMILREC',
     $  'PA_V3','RA_V1','DEC_V1'/
C----------------------------------------------------------------
        ID = IDS(6)
        CALL UHDGSD(ID,NAMES(1),EPCHTIME,ISTATS(1))
        CALL UHDGSD(ID,NAMES(2),SDMEANAN,ISTATS(2))
        CALL UHDGSD(ID,NAMES(3),CIRVELOC,ISTATS(3))
        CALL UHDGSD(ID,NAMES(4),COSINCLI,ISTATS(4))
        CALL UHDGSD(ID,NAMES(5),ECCENTRY,ISTATS(5))
        CALL UHDGSD(ID,NAMES(6),ECCENTX2,ISTATS(6))
        CALL UHDGSD(ID,NAMES(7),ECBDX4D3,ISTATS(7))
        CALL UHDGSD(ID,NAMES(8),ESQDX5D2,ISTATS(8))
        CALL UHDGSD(ID,NAMES(9),ECBDX3,ISTATS(9))
        CALL UHDGSD(ID,NAMES(10),FDMEANAN,ISTATS(10))
        CALL UHDGSD(ID,NAMES(11),RCASCNRV,ISTATS(11))
        CALL UHDGSD(ID,NAMES(12),ARGPERIG,ISTATS(12))
        CALL UHDGSD(ID,NAMES(13),MEANANOM,ISTATS(13))
        CALL UHDGSD(ID,NAMES(14),RCARGPER,ISTATS(14))
        CALL UHDGSD(ID,NAMES(15),RASCASCN,ISTATS(15))
        CALL UHDGSD(ID,NAMES(16),SINEINCL,ISTATS(16))
        CALL UHDGSD(ID,NAMES(17),SEMILREC,ISTATS(17))
        CALL UHDGSD(ID,NAMES(18),PSANGLV3,ISTATS(18))
        CALL UHDGSD(ID,NAMES(19),RTASCNV1,ISTATS(19))
        CALL UHDGSD(ID,NAMES(20),DECLNV1,ISTATS(20))
C
        DO 100 I=1,20
            IF(ISTATS(I).NE.0)THEN
                CONTXT='ERROR reading keyword '//NAMES(I)//
     *                  ' from .shh file'
                GO TO 999
            ENDIF
100     CONTINUE
C
        ISTAT=0
        GO TO 1000
999     CALL YMSPUT(CONTXT,STDOUT+STDERR,0,ISTAT)
        ISTAT=1
1000    RETURN
        END

max_stars_repo_name: iraf-community/stsdas


text:
!
! Test_CRTM_TauCoeff
!
! Program to test the CRTM_TauCoeff module and dependencies.
!
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, CIMSS/SSEC 04-Nov-2006
!                       paul.vandelst@ssec.wisc.edu
!

PROGRAM Test_CRTM_TauCoeff

  ! -----------------
  ! Environment setup
  ! -----------------
  ! Module usage
  USE Message_Handler, ONLY: SUCCESS, Display_Message
  USE CRTM_TauCoeff  , ONLY: CRTM_Load_TauCoeff, CRTM_Destroy_TauCoeff
  ! Disable all implicit typing
  IMPLICIT NONE

  ! ----------
  ! Parameters
  ! ----------
  CHARACTER(*), PARAMETER :: PROGRAM_NAME   = 'Test_CRTM_TauCoeff'
  CHARACTER(*), PARAMETER :: PROGRAM_RCS_ID = &
  INTEGER, PARAMETER :: NSENSORS=4
  CHARACTER(*), PARAMETER, DIMENSION(NSENSORS) :: SENSORID=&
    (/ 'amsua_n17', &
       'hirs3_n17', &
       'ssmis_f16', &
       'imgr_g11 ' /)

  ! ---------
  ! Variables
  ! ---------
  CHARACTER(256) :: Message
  INTEGER :: Error_Status
  INTEGER :: n

  ! Loop over sensors
  DO n = 1, NSENSORS
    ! Output info
    WRITE(*,'(/5x, "Reading TauCoeff data for:", 10(1x,a,:,","))') SENSORID(1:n)
    ! Message
    WRITE(Message,'("n=",i0)') n
    ! Load the TauCoeff data
    Error_Status = CRTM_Load_TauCoeff( SensorID=SENSORID(1:n) )
    IF ( Error_Status/=SUCCESS ) THEN
      CALL Display_Message( PROGRAM_NAME, &
                            'Error loading TauCoeff data for '//TRIM(Message), &
                            Error_Status )
      STOP
    END IF
    ! Destroy the TauCoeff data
    Error_Status = CRTM_Destroy_TauCoeff()
    IF ( Error_Status/=SUCCESS ) THEN
      CALL Display_Message( PROGRAM_NAME, &
                            'Error destroying TauCoeff data for '//TRIM(Message), &
                            Error_Status )
      STOP
    END IF
  END DO

  ! Default datafile read
  WRITE(*,'(/5x, "Reading default TauCoeff data file")')
  ! Load the TauCoeff data
  Error_Status = CRTM_Load_TauCoeff()
  IF ( Error_Status/=SUCCESS ) THEN
    CALL Display_Message( PROGRAM_NAME, &
                          'Error loading default TauCoeff data', &
                          Error_Status )
    STOP
  END IF
  ! Destroy the TauCoeff data
  Error_Status = CRTM_Destroy_TauCoeff()
  IF ( Error_Status/=SUCCESS ) THEN
    CALL Display_Message( PROGRAM_NAME, &
                          'Error destroying default TauCoeff data', &
                          Error_Status )
    STOP
  END IF

END PROGRAM Test_CRTM_TauCoeff

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm


      subroutine dr_db_read

      use dr_module
      use input_file_module
      use constituent_mass_module
      !use organic_mineral_mass_module
      use maximum_data_module

      character (len=80) :: titldum, header
      character (len=16) :: namedum
      integer :: eof, imax
      logical :: i_exist              !none       |check to determine if file exists

      eof = 0

      !read all delivery ratio data
      inquire (file=in_delr%del_ratio, exist=i_exist)
      if (i_exist .or. in_delr%del_ratio /= "null") then
        do
          open (107,file=in_delr%del_ratio)
          read (107,*,iostat=eof) titldum
          if (eof < 0) exit
          read (107,*,iostat=eof) header
          if (eof < 0) exit
          imax = 0
          do while (eof == 0)
            read (107,*,iostat=eof) titldum
            if (eof < 0) exit
            imax = imax + 1
          end do

          db_mx%dr = imax

          allocate (dr_db(imax))
          rewind (107)
          read (107,*,iostat=eof) titldum
          if (eof < 0) exit
          read (107,*,iostat=eof) header
          if (eof < 0) exit

        do ii = 1, imax
          read (107,*,iostat=eof) dr_db(ii)
          if (eof < 0) exit
        end do
          close (107)
          exit
        end do
      end if

      ! read delivery ratio data for all constituent types
      call dr_read_om
      if (cs_db%num_pests > 0) call dr_read_pest
      if (cs_db%num_paths > 0) call dr_path_read
      if (cs_db%num_metals > 0) call dr_read_hmet
      if (cs_db%num_salts > 0) call dr_read_salt

      return
      end subroutine dr_db_read


text:
subroutine  ALBEDO (SOLSIN)
  use simsphere_mod, only: albg, albf, wgg, wmax, frveg, xlai, albdoe, sigf, eq
  implicit none

  real :: SOLSIN
  logical :: ALGFLG = .false. ,ALFFLG = .false.
  integer :: ALBFLG = 0

!      INCLUDE 'modvars.h'

! **  Modified 1/8/90 to allow user to specify Albedoes.

! **  If the albedoes are ommitted in START then they are calculated,
! **  otherwise the albedo of the ground and the foliage input are used
! **  to calculate the weighted albedo (ALBDOE).
! **  Note: ALBF depends on solar angle (SOLSIN); ALBG on WGG.

  IF (ALBFLG == 0) THEN
    IF (eq(ALBG,0.0)) ALGFLG = .TRUE.
    IF (eq(ALBF,0.0)) ALFFLG = .TRUE.
    ALBFLG = 1
  ENDIF

  IF ( ALGFLG ) THEN
    ALBG = 0.25 - 0.20 * WGG / WMAX     !  Toby's Value
!	ALBG = 0.20 - 0.15 * WGG / WMAX     !  Fudged Dim
!	ALBG = 0.30 - 0.20 * WGG / WMAX     !  Fudged Bright
  ENDIF

  IF ( ALFFLG ) THEN
!	ALBF = 0.025 /( 0.1 + 0.1 * SOLSIN ) + ( 1 - SOLSIN**2 ) * 0.1	! T
    ALBF = 0.032 /( 0.1 + 0.1 * SOLSIN ) + ( 1 - SOLSIN**2 ) * 0.1 ! A
  ENDIF

  SIGF = 1 - EXP(-0.4 * XLAI )
  ALBDOE = (SIGF*ALBF + (1 - SIGF)*ALBG)*FRVEG + (1-FRVEG) * ALBG

  return
end

max_stars_repo_name: derekmorr/simsphere

max_issues_repo_name: derekmorr/simsphere

max_forks_repo_name: derekmorr/simsphere


text:
       SUBROUTINE INEDIT(BUF1,LBUF1,IBUF,ISCR)
c
C                             PROCESS INPUT (READS MAGNETIC TAPE)
C                             ------- -----
C
C  DOCUMENT DATE:
C
C     PROCESS INPUT READS SEISMIC DATA FROM MAGNETIC TAPE.  THE DATA MUST BE IN
C  SEGY FORMAT.  INPUT IS TRACE ORIENTED RATHER THAN SHOT OR RP ORIENTED AND
C  FEW LOGICAL BOUNDARIES EXIST.  PROCESS INPUT READS A SINGLE TRACE FROM TAPE
C  AND IMMEDIATELY PASSES IT ON TO THE NEXT PROCESS IN THE PROC LIST.
C      THE DATA MAY BE READ FROM TAPES USING ONE OF THREE DIFFERENT CONTROLS.
C  FIELD TAPES MAY BE CONTROLLED BY EITHER SHOT NUMBERS OR GMT TIME.  TAPES THAT
C  HAVE BEEN PROCESSED THROUGH PROCESS GATHER MUST BE CONTROLLED BY RP NUMBERS
C  STACKED TAPES MAY BE READ USING EITHER RP NUMBERS OR BY GMT TIME.
C      INPUT REEL CHANGES ARE HANDLED BY PROCESS INPUT THROUGH OPERATOR
C  INTEVENTION WHEN EITHER END OF TAPE IS REACHED (A FILE MARK) OR WHEN THE
C  USER SPECIFIES A NEW REEL NUMBER VIA THE IREELN PARAMETER.  OPERATOR REEL
C  CHANGES ARE ACCOMPLISHED BY THE COMPUTER OPERATOR SPECIFYING THE NEW TAPE
C  UNIT NUMBER VIA FILE 'IN'.  THE FILE 'IN' MUST BE IN THE SAME DIRECTORY
c  THAT THE JOB IS RUNNING IN.  A -1 UNIT NUMBER TERMINATES THE JOB.
c      SIOSEIS allows the user to break" the SEGY standard by putting more
c  than one file on a tape (some Exabyte tapes hold 5GB of data!).  In order
c  to read multiple files, INPUT parameter NFILES MUST be used.  Remember
c  that tapes are terminated by a file mark; in order to read multiple
c  tapes NFILES must be used.
C      SIOSEIS ASSUMES THAT SHOT NUMBERS INCREASE ON THE TAPE unless parameter
c  order is used.  IF THE SHOT NUMBERS DO NOT INCREASE, OR IF THE SHOTS ARE
c  NOT TO BE PROCESSED IN THE ORDER THEY ARE ON TAPE, MULTIPLE PARAMETER LIST
c  MAY BE USED.  E.G. IF THE SHOTS ON TAPE ARE ORDERED 10001 TO 11000 FOLLOWED
c  BY SHOTS 500 TO 600, USE EITHER FIS 99999     OR USE TWO LISTS -
c                                  FIS 10001 LIS 11000 END
C                                  FIS 500 LIS 600 END
C
C  THE PARAMETER DICTIONARY
C  --- --------- ----------
C
C  REQUIRED PARAMETERS
C  -------- ----------
C  IUNIT  - THE TAPE UNIT NUMBER OF THE FIRST INPUT REEL.  Tape units may be
c           changed by changing iunit in a new fno/lno list or by the above
c           method (file in) when a file mark is detected on the SEGY input tape.
c           Unix tape assignment is /dev/nrstIUNIT. e.g. iunit 0 is /dev/nrst0.
c           VMS tape assignments are:
c               0<= LUN < 10  tape drive MSA, 10<= LUN < 20  tape drive MSB
c              20<= LUN < 30  tape drive MTA, 30<= LUN < 40  tape drive MTB
c              40<= LUN < 50  tape drives MUA, 50<= LUN < 60  tape drives MUB
c              60<= lun < 70  tape drives XTA, 70<= lun < 80  tape drives XTB
C           Preset 0.        E.G. IUNIT 1
c  DEVICE - The UNIX tape device name to read from.  When given, the
c           Fortran magnetic tape device independent driver is used to
c           read the specified device.  The "no rewind" device should be used.
c           Preset = none       e.g.  device  /dev/nrst8
C
C  PARAMETERS FOR PROCESSING BY SHOT
C  ---------- --- ---------- -- ----
C  FIS    - THE FIRST INPUT SHOT NUMBER OF A SET OF SHOTS (OF A FIS-LIS-END LIST).
C           A NEGATIVE FIS MEANS THAT THE FIRST SHOT OF THE SET IS IN A REVERSE
C           OR BACKWARDS DIRECTION ON TAPE FROM THE PREVIOUS SHOT INPUT.
c           The shots will be read from tape in the order specified by:
c           FIS to LIS in increments of SINC.  E.G. fis 1 lis 5 means
c           that the shots must be ordered 1, 2, 3, 4, 5 on tape.  If a
c           shot within FIS - LIS is missing, the program will continue
c           to search for the missing shot and will stop at the end of
c           tape without finding the shot.
C           A FIS OF 99999 WILL CAUSE THE PROGRAM TO READ ALL SHOTS ON THE
C           INPUT TAPE REGARDLESS OF SHOT NUMBER (THE SHOT NUMBERS DO NOT
C           HAVE TO BE CONSECUTIVE).  USE OF FIS 99999 CAUSES LIS TO BE
C           IGNORED, THUS THE ONLY WAY THE JOB CAN BE STOPPED IS AT THE
C           REEL CHANGE WHEN THE OPERATOR RESPONDS WITH A -1 UNIT NUMBER.
C           PRESET = 99999.  E.G. FIS 41
C  LIS    - THE LAST INPUT SHOT NUMBER OF A SET OF SHOTS (OF A LIST).
C           PRESET = 32767.  E.G. LIS 201
C  NIS    - THE NUMBER OF INPUT SHOTS TO READ WHEN FIS = 99999.
C           PRESET = 32767  E.G. NIS 5
C  SINC   - THE SHOT INCREMENT BETWEEN THE FIRST SHOT (IABS(FIS)) AND THE LAST
C           SHOT (LIS) OF A LIST.  ALTERNATE SHOTS MAY BE PROCESSED BY USING
C           SINC 2.  A NEGATIVE SINC MUST BE USED IF SHOT LIS IS BEFORE
C           SHOT FIS ON TAPE.
C           PRESET = 1.  E.G. SINC 2
C  FTR    - THE FIRST TRACE OF EACH SHOT TO PROCESS.  A TRACE WITH A TRACE
C           NUMBER ON TAPE LESS THAN FTR WILL BE IGNORED.  An FTR of 99999
c           indicates that all traces will be read regardless of trace
c           number or trace order.
C           PRESET = 1.  E.G. FTR 23
C  LTR    - THE LAST TRACE OF EACH SHOT TO PROCESS.  A TRACE WITH A TRACE
C           NUMBER ON TAPE LARGER THAN LTR WILL BE IGNORED.
C           PRESET = TAPE HEADER.  E.G. LTR 23
C  TRINC  - THE INCREMENT BETWEEN TRACES FTR AND LTR.  ALTERNATE TRACES MAY
C           BE PROCESSED BY USING TRINC 2.
C           PRESET = 1.  E.G. TRINC 24
C
C  PARAMETERS FOR PROCESSING BY RP (CDP GATHERS)
C  ---------- --- ---------- -- -- ---- --------
C  FRP    - THE FIRST INPUT RP NUMBER OF A SET OF RPS (OF A FRP-LRP-END LIST).
C           A NEGATIVE FRP MEANS THAT THE FIRST RP OF THE SET IS IN A REVERSE
C           OR BACKWARDS DIRECTION ON TAPE FROM THE PREVIOUS RP INPUT.
C           PRESET = 99999.  E.G. FRP 41
C  LRP    - THE LAST INPUT RP NUMBER OF A SET OF RPS (OF A LIST).
C           PRESET = 32767.  E.G. LRP 201
C  RPINC  - THE RP INCREMENT BETWEEN THE FIRST RP (IABS(FRP)) AND THE LAST
C           RP (LRP) OF A LIST.  ALTERNATE RPS MAY BE PROCESSED BY USING
C           RPINC 2.  A NEGATIVE RPINC MUST BE USED IF RP LRP IS IN FRONT OF
C           RP FRP
C            PRESET = 1.  E.G. rpinc 2
C  FTR    - THE FIRST TRACE OF EACH RP TO PROCESS.  A TRACE WITH A TRACE
C           NUMBER ON TAPE LESS THAN FTR WILL BE IGNORED.  An FTR of 99999
c           indicates that all traces will be read regardless of trace
c           number or trace order.
C            PRESET = 1.  E.G. FTR 23
C  LTR    - THE LAST TRACE OF EACH RP TO PROCESS.  A TRACE WITH A TRACE
C           NUMBER ON TAPE LARGER THAN LTR WILL BE IGNORED.
C            PRESET = TAPE HEADER.  E.G. LTR 23
C  TRINC  - THE INCREMENT BETWEEN TRACES FTR AND LTR.  ALTERNATE TRACES MAY
C           BE PROCESSED BY USING TRINC 2.
C            PRESET = 1.  E.G. TRINC 24
C
C  PARAMETERS FOR PROCESSING BY TIME (GMT)
C  ---------- --- ---------- -- ---- -----
C  ALL SHOTS PRIOR TO FDAY, FGMT ARE IGNORED.  ALL SHOTS ARE THEN PROCESSED
C  AS THEY ARE READ FROM TAPE UNTIL LDAY, LGMT IS EXCEEDED.
C
C  FDAY   - THE FIRST JULIAN DAY OF THE DATA TO BE READ FROM TAPE.
C           PRESET=0.  E.G. FDAY 364
C  LDAY   - THE LAST JULIAN DAY OF THE DATA TO BE PROCESSED.
C           PRESET=FDAY
C  FGMT   - THE FIRST TIME OF DAY FDAY OF THE DATA TO BE READ FROM TAPE.
C           GMT IS EXPRESSED AS HHMM OR HOURS AND MINUTES OF THE 24 HOUR
C           CLOCK (0000-2359).
C           PRESET=0.  E.G. FGMT 1400
C  LGMT   - THE LAST TIME OF LDAY TO BE READ FROM TAPE.
C           PRESET=0.  E.G. 0605
C  FTR    - THE FIRST TRACE OF EACH SHOT TO PROCESS.  A TRACE WITH A TRACE
C           NUMBER ON TAPE LESS THAN FTR WILL BE IGNORED.  An FTR of 99999
c           indicates that all traces will be read regardless of trace
c           number or trace order.
C           PRESET = 1.  E.G. FTR 2
C  LTR    - THE LAST TRACE OF EACH SHOT TO PROCESS.  A TRACE WITH A TRACE
C           NUMBER ON TAPE LARGER THAN LTR WILL BE IGNORED.
C           PRESET = ALL TRACES.  E.G. LTR 1
C  TRINC  - THE INCREMENT BETWEEN TRACES FTR AND LTR.  ALTERNATE TRACES MAY
C           BE PROCESSED BY USING TRINC 2.
C           PRESET = 1.  E.G. TRINC 24
C
C  USEFUL OPTIONAL PARAMETERS
C  ------ -------- ----------
c  FNO    - First shot/rp number. FNO may be used instead of FIS or FRP.
c           Ordinarily SIOSEIS insists that the users know whether the
c           input tape is sorted by shot or rp and will not run if the
c           tape is described incorrectly.  The use of FNO and LNO eases
c           this!  As with FIS and LIS, FNO 99999 reads all traces from
c           tape regardless of their numbers and will stop reading only
c           when a negative tape unit number is given in file in at the
c           end of a tape (see tape change).
c           Preset = none          e.g.   fno 101
c  LNO    - Last shot/rp to read.
c           Preset = 32768          e.g.  lno 200
c  NOINC  - The increment between FNO and LNO.  The shots/rps on tape
c           will be read in the order FNO to LNO in increments of NOINC.
c           IF a shot/rp is missing, PROCESS INPUT will continue searching
c           until it finds it.  E.G.  FNO 1 LNO 5 but the shot/rps on
c           tape are 1, 2, 4, 5  then SIOSEIS will search the entire
c           tape for shot 3, stopping when it finds end of tape, and
c           never reading shots 4 or 5.  Use FNO 99999 when this occurs.
C  RENUM  - RENUMBERS EVERY SHOT/RP SO THAT THE NUMBERS ARE MONOTONICALLY
C           INCREASING BY 1, STARTING WITH THE NUMBER SPECIFIED BY RENUM.
C           RENUMBERING THE SHOT/RP NUMBERS IS HELPFUL WHEN THE INPUT SHOT/RP
C           NUMBERS ARE NONUNIQUE AND SOME PARAMETER NEEDS TO BE SPATIALLY
C           VARIED.  THE RENUMBERING IS DONE AFTER THE SHOT/RP IS READ.
C           DEFAULT = not given. E.G. RENUM 1
C  C      - COMMENT CARD IMAGES TO REPLACE IN THE TAPE HEADER.  THE COMMENT
C           CARD MUST START WITH THE LETTER C AND MUST BE FOLLOWED IMMEDIATELY
C           BY A 2 DIGIT NUMBER, FOLLOWED BY A BLANK.  THE NUMBER IS THE CARD
C           NUMBER WITHIN THE HEADER TO REPLACE.  THE COMMENT ITSELF MUST BE
C           ENCLOSED IN SINGLE QUOTES.
C           PRESET=NONE. EXAMPLE C15 'THIS IS AN EXAMPLE OF A COMMENT'
c  STIME  - Start time of the data to process.  Data prior to STIME will be
c           discarded and the deep water delay will be set to STIME.  If
c           STIME is less than the deep water delay, STIME will be ignored.
c           Preset = delay of each trace     e.g.   stime 3.0
C  SECS   - THE NUMBER OF SECONDS OF DATA TO PROCESS.  SECS IS HELD CONSTANT
C           FOR ALL TRACES WITHIN THE JOB.  THE TOTAL TIME LENGTH OF EACH TRACE
C           IS THE SUM OF THE DELAY AND SECS.  If the data is not SECS
c           long, zeroes are NOT added or padded and the trace will be
c           less tha SECS long.
C           PRESET= ALL THE DATA ON THE TRACE.
c  SET    - The start and end times of the data to read from disk.  SET
c           is a pair of times in seconds.  The use of SET causes the
c           deep water delay and the number of samples to be changed.
c           If either set is outside of the data, the data is padded
c           with zeroes.  The data will always be SET(2) - SET(1) long.
c           Preset = none     e.g.   set 2.0 3.0
C  IREELN - THE INPUT REELN NUMBER.  THIS IS USEFUL WHEN THE INPUT TAPE SHOULD
C           BE CHANGED AND THE PREVIOUS TAPE IS NOT AT THE END OF TAPE (2 FILE
C           MARKS).  E.G.  FIS 1 LIS 10 IREELN 1 END
C                          FIS 11 LIS 20 IREELN 2 END
C           CAUSES A REEL CHANGE TO OCCUR AFTER SHOT 10 IS READ.  THE NORMAL
C           REELN CHANGE PROCEDURES SHOULD BE FOLLOWED (USING FILE IN).
C  IFMT   - THE INPUT TAPE FORMAT.  USED TO OVERRIDE THE FORMAT INDICATOR ON THE
C           TAPE ITSELF.
C         =0, USE THE VALUE FROM THE TAPE.
C         =1, IBM 360 32 BIT HEX BASED FLOATING POINT.
C         =2, 32 BIT INTEGER.
C         =3, 16 BIT INTEGER.
C         =4, 16 BIT FIXED POINT WITH GAIN CODES. (not really SEGY
C         >4, HOST BIT FLOATING POINT
C           PRESET=0
C  SI     - THE INPUT SAMPLE INTERVAL IN SECONDS.  USED TO OVERRIDE THE
C           SAMPLE INTERVAL CONTAINED IN THE TRACE HEADER.
C           PRESET=TAPE
C  DELAY  - THE INPUT DEEP WATER DELAY IN SECONDS.  USED TO OVERRIDE THE
C           DELAY CONTAINED IN THE TRACE HEADER.
C           PRESET=TAPE
C  NTRCS  - THE NUMBER OF TRACES PER INPUT RECORD (SHOT OR RP).  USED TO
C           OVERRIDE THE NUMBER OF TRACES CONTAINED IN THE BINARY TAPE HEADER.
C           PRESET=TAPE HEADER.
C  NTRGAT >0, NUMBER OF TRACES PER GATHER.  NTRGAT CONVERTS A SHOT SORTED TAPE
C             INTO A RP SORTED TAPE, OR CONVERTS A FOREIGN RP SORTED TAPE INTO
C             A SIOSEIS RP SORTED TAPE.  EACH GATHER WILL CONTAIN NTRGAT TRACES.
C             THE RP NUMBERS WILL BE THE SAME AS THE SHOT NUMBERS.  EVERY
C             NTRGAT TRACES WILL BE FLAGGED AS THE END OF GATHER, SO PROCESS
C             STACK WILL STACK NTRGAT TRACES.  THIS PARAMETER IS USEFUL FOR
C             PROCESSING TAPES THAT HAVE BEEN GATHERED OUTSIDE OF SIOSEIS.  THIS
C             IS ALSO USEFUL FOR CONVERTING SINGLE CHANNEL SHOTS INTO "ZERO-
C             OFFSET" STACKED DATA (NEEDED BY MIGRATION IN DIGICON'S DISCO SYSTEM).
c          <0, CONVERTS TAPES SORTED BY RP (IF THE RP HAS A NON ZERO TRACE
C              NUMBER) TO A TAPE SORTED BY SHOT (RP TRACE NUMBER (SEGY TRACE
C              HEADER WORD 7) TO ZERO.
C           PRESET=0.  E.G.  NTRGAT 1
c  FORGAT - Foreign gather tape switch.  The use of FORGAT indicates that the
c           input gather tape was not generated by SIOSEIS and thus does not
c           have the end-of-gather convention used by SIOSEIS (a -1 in SEGY
c           header word 51).  FORGAT is similar to NTRGAT but allows each
c           gather to have a different number of traces.  FORGAT causes process
c           INPUT to wait for the next trace to be read from tape and sets the
c           end-of-gather flag to -1 if the next gather number is different,
c           or the next trace number is greater than LTR.
c           End of gather will also be set if the parameter LTR is used and
c           the trace on tape is equal to LTR.
c           The value of forgat indicates the number of rps to concatenate
c           into a single gather which is terminated by the end-of-gather
c           flag.  The penalty for this is significantly slower throughput on
c           non-Unix computers (computers that allow asyncronus or overlapped I/O).
c           LDGO gather tapes start with the largest trace number first, which
c           breaks the SIOSEIS monotonicly increasing assumption; LDGO gather
c           tapes may be read using forgat 1 and ftr 99999, in which case
c           SIOSEIS will use all traces within the gather.
c           Preset = 0       e.g. forgat 1
C  DECIMF - DECIMATION FACTOR.  THE DATA WILL BE DECIMATED OR RESAMPLED BY A
C           FACTOR OF DECIMF, THUS REDUCING THE NUMBER OF SAMPLES IN THE TRACE.
C           **  NOTE **  NO ANTI-ALIAS FILTER IS APPLIED PRIOR TO DECIMATION.
C           PRESET=1.  E.G. DECIMF 2   TAKES EVERY SECOND SAMPLE
c  NFSKIP - The number of files to skip before reading the SEGY tape header.
c           Some places (LDGO) sometimes put multiple lines or logical tapes
c           on a single physical tape, a violation of the SEGY standard!
c           Preset = 0.   e.g.   nfskip 1        skips 1 file
c  REWIND - Rewinds the input tape PRIOR to reading the first shot (or
c           doing tape positioning with NFSKIP).  REWIND is set to 1,
c           or TRUE, after the first tape so that subsequent tapes
c           will be rewound when using the "operator" tape change
c           method (file "in").  Likewise, REWIND is set to 1 (TRUE)
c           on every INPUT list (a list is terminated with "END").
c         = 1, (TRUE), rewind.
c         = 0, (FALSE), NO rewind.
c           DEFAULT = 1,      e.g.  rewind 0    # don't rewind
c  NRSKIP - The number of tape records (traces) to skip before reading
c           the first trace.  This is useful when 2 or more SEGY files
c           are concatenated into 1 file on tape.
c         >0, the record skip is done AFTER the SEGY tape header is read.
c         <0, the record skip is done BEFORE the SEGY tape header is read.
c           Preset = 0.   e.g.   nrskip 3        skips 3 records
c  ORDER  - SIOSEIS normally assumes that the shot/rp numbers are in assending
c           order on tape.  In order to reverse this assumption, set order
c           to -1, which will then force process INPUT to assume that the
c           shot/rp numbers are in decreasing order on tape.  The parameter
c           sinc/rpinc/noinc usually needs to be a negative number when
c           order -1 is given.  e.g. If the order of rps on tape is 9100,
c           9099, 9098, 9097,...  the following parameters would work:
c           frp 9100 lrp 9097 rpinc -1
c           Preset = 0      e.g. order -1
c  NFILES - The number of files to use.  The tape files must be
c           consecutive on tape.  A filemark on tape is the end of file
c           indicator.  Two consecutive file marks are considered end of
c           information and will cause input to rewind the tape and
c           expect a tape change using the normal SIOSEIS tape change
c           procedures (by creating a file named in).  When used with
c           parameter DEVICE, the "no rewind" device MUST be used.
c           PRESET = 1       e.g.  nfiles 2
C
C  WRITTEN AND COPYRIGHTED (C) BY:
C  PAUL HENKART, SCRIPPS INSTITUTION OF OCEANOGRAPHY, JANUARY 1980
C  ALL RIGHTS ARE RESERVED BY THE AUTHOR.  PERMISSION TO COPY OR REPRODUCE THIS
C  SUBROUTINE, BY COMPUTER OR OTHER MEANS, MAY BE OBTAINED ONLY FROM THE AUTHOR.
C
C
C  ARGUMENTS:
C  BUF1   - AN ARRAY THAT WILL RECEIVE THE FIRST TRACE OF THE JOB WITH THE
C           TRACE HEADER.
C  LBUF1  - THE SAME ARRAY BUT IN LONG INTEGER SINCE ARGUMENTS CAN'T BE
C           EQUIVALENCED.
C  IBUF   - THE SAME ARRAY BUT IN SHORT INTEGER SINCE ARGUMENTS CAN'T BE
C           EQUIVALENCED.
C  ISCR   - AN ARRAY TO BE USED AS SCRATCH - IE. THE CONTENTS WILL BE DESTROYED.
C
c   mod 1 Mar 91 to add nrskip parameter
c   mod 5 Mar 91 to add stime parameter
c   mod 24 July 91 to add nfiles parameter
c   mod 27 Aug 91 to add device
c   mod 17 Nov 91 change nfiles preset to 99999
c   mod 11 Feb 92 to allow fis/lis on rp sorted tapes and clean up
c           fis/lrp. frp/lrp, fno/lno presets and error messages.
c   mod 5 Oct 92 to add parameter set
c   mod 22 Oct 92 to allow renum to change from list to list.
c   mod 24 Nov 92 change nfiles preset to 1
c   mod 10 Jan 94 - reject nfiles to be used with the device parameter.
c   mod 9 June 94 - Don't reject nfiles and device if norewind device.
c   mod 22 Jun 94 - set iunit = 1 when DEVICE is given
c   mod 20 Mar 95 - Add no =   and tr =
c   mod 29 Mar 95 - Change fis/frp/fno preset to 99999
c   mod 29 Mar 95 - Add parameters allno and alltr
c   mod 15 aug 95 - Allow iunit up to 40
c   mod 13 Nov 95 - Add parameter REWIND.
c   mod 7 Dec 95 - Do some editing on NO and TR parameters.
c   mod 18 Mar 96 - allno 0 and alltr 0 didn't work
c   mod 1 Oct 97 - Use upcase for  allno, alltr, rewind etc
c   mod mar 99 - Change allno to no if there's an increment between shots etc
c              - Likewise for trinc - set alltr to 0 if trinc is given
c   mod 31 Jul 00 - mar 99 change was bad because sinc preset is 1 not 0!
c   mod 16 May 01 - Allow iunit less than 100 (was 40).
c   mod 31 Jul 02 - Add warning if DEVICE is not 'bn'
c   mod 8 May 03 - The above isn't true on SGI
c   mod 26 Jan 04 - Change nfiles preset to 0 so that nfiles terminates
c                   without the tape change question and file in.
c   mod 17 May 10 - Comment out creating real mute times in header words 47 & 48
c
      PARAMETER (NPARS = 44 )                                           ! THE NUMBER OF PARAMETERS AVAILABLE TO THE USER
      INTEGER FIS,SINC,FRP,RPINC,FTR,TRINC,RENUM,NIS,fno,forgat
      INTEGER allno, alltr, rewind, rewindi
      CHARACTER*3 BYPASS
      CHARACTER*4 LTYPE
      COMMON/INPUT/ IPARUN ,NLISTS
      COMMON /porder/ num, iorder(40)
      LOGICAL FIRST,CFIRST
      COMMON /READT/ILUN,NUMHDR,NUMDAT,IUNHDR,IREELM,INTRCS,JFMT,NSKIP,
     *    SSECS,idummy,ISRCF,IDTYPE,nfskp,jform,itxsi,itxdel,
     &    nfktrc, norigtr, nrskp, nfile, rewindi
      COMMON /SIOAP/ IASGND,IRELSE,IN,IOUT,NEXTAD,LAPSIZ
      INTEGER FDAY,LDAY,FGMT,LGMT,DECIMF,order
      CHARACTER*6 NAMES(NPARS)
      CHARACTER*1 TYPE(NPARS)
      DIMENSION LENGTH(NPARS)
      CHARACTER*80 TOKEN
c****  device is not easy.  Put it in common so rdtrc has access to it.  Use a
c****  separate common because of the word boundary problem with characters.
c****  Write it to disk too, so that it can change from list to list
      COMMON /inputdev/inputdev
      CHARACTER*80 inputdev
      DIMENSION VALS(NPARS),LVALS(NPARS)
      COMMON /EDITS/ IERROR,IWARN,IRUN,NOW,ICOMPT
      CHARACTER*80 CHEADR(40),COMMNT(40)
      DIMENSION BUF1(1000),LBUF1(1000), set(2)
      INTEGER*2 IBUF(1000),ISCR(1000)
C
C
      EQUIVALENCE (IREELN,LVALS(1)),
     2            (IUNIT,LVALS(2)),
     3            (BYPASS,LVALS(3)),
     4            (FIS,LVALS(4)),
     5            (LIS,LVALS(5)),
     6            (SINC,LVALS(6)),
     7            (FTR,LVALS(7)),
     8            (LTR,LVALS(8)),
     9            (TRINC,LVALS(9)),
     *            (FRP,LVALS(10)),
     1            (LRP,LVALS(11)),
     2            (RPINC,LVALS(12)),
     3            (LPRINT,LVALS(13)),
     4            (SECS,VALS(14)),
     5            (IFMT,LVALS(15)),
     6            (SI,VALS(16)),
     7            (DELAY,VALS(17)),
     8            (NTRCS,LVALS(18)),
     9            (NTRGAT,LVALS(19)),
     *            (RENUM,LVALS(20))
      EQUIVALENCE (FDAY,LVALS(21)),                                     ! F77 ALLOW A MAX OF 19 CONTINUATION LINES!!
     2            (LDAY,LVALS(22)),
     3            (FGMT,LVALS(23)),
     4            (LGMT,LVALS(24)),
     5            (DECIMF,LVALS(25)),
     6            (NIS,LVALS(26)),
     7            (nfskip,lvals(27)),
     8            (fno,lvals(28)),
     9            (lno,lvals(29)),
     *            (forgat,lvals(30)),
     1            (noinc,lvals(31)),
     2            (order, lvals(32)),
     3            (nrskip, lvals(33)),
     4            (stime, vals(34)),
     5            (nfiles,lvals(35)),
     6            (device,lvals(36)),
     7            (allno,lvals(42)),
     8            (alltr,lvals(43)),
     9            (rewind,lvals(44))
      DATA NAMES/'IREELN','IUNIT ','BYPASS','FIS   ','LIS   ','SINC  ',
     1  'FTR   ','LTR   ','TRINC ','FRP   ','LRP   ','RPINC ','LPRINT',
     2  'SECS  ','IFMT  ','SI    ','DELAY ','NTRCS ','NTRGAT','RENUM ',
     3  'FDAY  ','LDAY  ','FGMT  ','LGMT  ','DECIMF','NIS   ','NFSKIP',
     4  'FNO   ','LNO   ','FORGAT','NOINC ','ORDER ','NRSKIP','STIME ',
     5  'NFILES','DEVICE','SET   ','NO    ','NO=   ','TR    ','TR=   ',
     6  'ALLNO ','ALLTR ','REWIND' /
      DATA LENGTH/6,5,6,3,3,4,3,3,5,3,3,5,6,4,4,2,5,5,6,5,4,4,4,4,6,3,
     *    6,3,3,6,5,5,6,5,6,6,3,2,3,2,3,5,5,6/
      DATA TYPE/'L','L','A',10*'L','F','L','F','F',16*'L','F','L','A',
     &        'F','A','F','A','F',2*'A','A'/
      DATA FIRST /.TRUE./, CFIRST /.TRUE./
      DATA COMMNT/40*' '/
C****
C****   PRESETS
C****
      IUNIT = 0
      BYPASS='NO'
      FIS=-32767
      LIS = 32767
      SINC=1
      FTR=1
      LTR=-32767
      TRINC=1
      FRP=-32767
      LRP = 32767
      RPINC=1
      LPRINT=0
      SECS=0.
      IFMT=0
      SI=0.
      DELAY=-32.
      NTRCS=0
      NTRGAT=0
      RENUM=-32767
      DECIMF=1
      FDAY=0
      LDAY=0
      FGMT=-1                                                           ! ALLOW A GMT OF 0000
      LGMT=-1
      NIS=-32767
      TIME1=0.
      TIME2=0.
      IREELN=0
      NUMDAT=0                                                          ! SOMEONE HAS TO SET IT!
      nskip = 0
      nfskip = 0
      nrskip = 0
      fno = -32767
      lno = 32767
      forgat = 0
      noinc = 1
      order = 1
      stime = -1
      nfiles = 0
      inputdev = ' '
      iread = 0                                                         ! is input in the procs list?
      set(1) = 0.
      set(2) = 0.
      notype = 0
      noindex = 0
      itrtype = 0
      itrindex = 0
      allno = 1
      alltr = 1
      rewind = 1
      DO 10 i = 1, num
         IF( iorder(i) .EQ. 2 ) iread = 1
   10 CONTINUE
C****
C****   SET UP A DISC FILE TO HOLD THE PARAMETERS
C****
      CALL GETFIL(2,IPARUN,token,ISTAT)                                 ! GET A FREE PRIMOS FILE UNIT NUMBER
      IF(ICOMPT.EQ.1) IPARUN=IPARUN+12                                  ! PRIME COMPUTERS HAVE FILE UNIT IDIOSYNCRASIES
      OPEN( UNIT=iparun, STATUS='SCRATCH', FORM='UNFORMATTED' )
C****
C****   THE CURRENT COMMAND LINE IN THE SYSTEM BUFFER MAY HAVE THE PARAMETERS.
C****   GET A PARAMETER LIST FROM THE USER.
C****
      NTOKES=1
  100 CONTINUE
      ns = 0
      CALL GETOKE(TOKEN,NCHARS)                                         ! GET A TOKEN FROM THE USER PARAMTER LINE
      CALL UPCASE(TOKEN,NCHARS)
      IF( NCHARS .LE. 0 ) THEN
          IF(NOW.EQ.1) PRINT 140
  140        FORMAT(' <  ENTER PARAMETERS  >')
          CALL RDLINE                                                   ! GET ANOTHER USER PARAMETER LINE
          NTOKES=0
          GO TO 100
      ENDIF
      NTOKES=NTOKES+1
  150 DO 190 I=1,NPARS                                                  ! SEE IF IT IS A PARAMETER NAME
      LEN=LENGTH(I)                                                     ! GET THE LEGAL PARAMETER NAME LENGTH
      IPARAM=I                                                          ! SAVE THE INDEX
      IF(TOKEN(1:NCHARS).EQ.NAMES(I)(1:LEN).AND.NCHARS.EQ.LEN) GO TO 200
  190 CONTINUE                                                          ! STILL LOOKING FOR THE NAME
      IF(TOKEN(1:NCHARS).EQ.'END'.AND.NCHARS.EQ.3) GO TO 1000           ! END OF PARAM LIST?
      IF( TOKEN(1:1) .EQ. 'C' .AND .NCHARS .EQ. 3 ) THEN                ! IS IT A COMMENT CARD
          CFIRST = .FALSE.
          READ( TOKEN, '(1X,I2)')  ITEMP
          COMMNT(ITEMP)(1:4)=TOKEN(1:NCHARS)
          CALL GETOKE(TOKEN,NCHARS)                                     ! GET THE COMMENT (WITHIN THE QUOTES)
          COMMNT(ITEMP) (5:80) = TOKEN(1:NCHARS)
          GO TO 100
      ENDIF
      PRINT 191, TOKEN(1:NCHARS)
  191 FORMAT(' ***  ERROR  *** INPUT DOES NOT HAVE A PARAMETER ',
     *  'NAMED ',A)
      IERROR=IERROR+1
      GO TO 100
****
C****    FOUND THE PARAMETER NAME, NOW FIND THE VALUE
C****
  200 CONTINUE
      NPARAM=IPARAM
  210 CONTINUE                                                          !  NOW FIND THE VALUE
      CALL GETOKE(TOKEN,NCHARS)
      NTOKES=NTOKES+1
      IF( NCHARS .LE. 0 ) THEN                                          ! END OF LINE?
          IF( NOW .EQ. 1 ) PRINT 140                                    ! THIS ALLOWS A PARAMETER TO BE ON A DIFFERENT LINE FROM THE NAME
          CALL RDLINE                                                   ! GET ANOTHER LINE
          NTOKES=0
          GO TO 210
      ENDIF
      IF( TYPE(NPARAM) .EQ. 'A' ) THEN
          IF(NPARAM.EQ.3) BYPASS=TOKEN(1:NCHARS)
          IF( names(nparam) .EQ. 'DEVICE') THEN
              inputdev = token(1:nchars)
              iunit = 1
c****  sgi needs    nrns  (no-rewind, no-swap)
              IF( token(nchars-3:nchars) .EQ. 'nrns' ) GOTO 100
c****  sun needs   bn
              DO i = 1, nchars-1
                 IF( token(i:i+1) .EQ. 'bn' ) GOTO 100
              ENDDO
              PRINT *,
     &    ' ***  WARNING  ***  DEVICE should be BSD and No Rewind (bn).'
              iwarn = iwarn + 1
              GOTO 100
          ENDIF
          IF( names(nparam) .EQ. 'NO' ) THEN
              CALL upcase( token, nchars)
              IF( token(1:1) .EQ. 'I' ) THEN
                  notype = 1
              ELSEIF( token(1:1) .EQ. 'L' ) THEN
                  notype = 2
              ELSEIF( token(1:1) .EQ. 'R' ) THEN
                  notype = 3
              ELSE
                  PRINT *,' ***  ERROR  ***  NO value must start with ',
     &                  ' I or L or R.'
                  ierror = ierror + 1
              ENDIF
              token(1:nchars) = token(2:nchars)
              nchars = nchars - 1
              CALL dcode( token, nchars, areal, istat )
              IF( istat .NE. 2 ) THEN
                  ierror = ierror + 1
                  GOTO 100
              ENDIF
              noindex = NINT(areal)
          ENDIF
          IF( names(nparam) .EQ. 'TR' ) THEN
              CALL upcase( token, nchars)
              IF( token(1:1) .EQ. 'I' ) THEN
                  itrtype = 1
              ELSEIF( token(1:1) .EQ. 'L' ) THEN
                  itrtype = 2
              ELSEIF( token(1:1) .EQ. 'R' ) THEN
                  itrtype = 3
               ELSE
                  PRINT *,' ***  ERROR  ***  TR value must start with ',
     &                  ' I or L or R.'
                  ierror = ierror + 1
              ENDIF
              token(1:nchars) = token(2:nchars)
              nchars = nchars - 1
              CALL dcode( token, nchars, areal, istat )
              IF( istat .NE. 2 ) THEN
                  ierror = ierror + 1
                  GOTO 100
              ENDIF
              itrindex = NINT(areal)
          ENDIF
          IF( names(nparam) .EQ. 'ALLNO' ) THEN
              CALL upcase( token, nchars)
              IF( token(1:1) .EQ. 'Y' ) allno = 1
              IF( token(1:1) .EQ. 'N' ) allno = 0
          ENDIF
          IF( names(nparam) .EQ. 'ALLTR' ) THEN
              CALL upcase( token, nchars)
              IF( token(1:1) .EQ. 'Y' ) alltr = 1
              IF( token(1:1) .EQ. 'N' ) alltr = 0
          ENDIF
          IF( names(nparam) .EQ. 'REWIND' ) THEN
              CALL upcase( token, nchars)
              IF( token(1:nchars) .EQ. 'YES' ) rewind = 1
              IF( token(1:nchars) .EQ. 'NO' ) rewind = 0
              IF( token(1:nchars) .EQ. 'ON' ) rewind = 1
              IF( token(1:nchars) .EQ. 'OFF' ) rewind = 0
              IF( token(1:nchars) .EQ. '1' ) rewind = 1
              IF( token(1:nchars) .EQ. '0' ) rewind = 0
          ENDIF
	     GOTO 100
      ENDIF
      CALL UPCASE(TOKEN,NCHARS)
      CALL DCODE(TOKEN,NCHARS,AREAL,ISTAT)                              ! TRY AND DECODE IT
      IF( ISTAT .NE. 2 ) THEN                                           ! =2 MEANS IT IS A NUMERIC
          IERROR=IERROR+1                                               ! DCODE PRINTED AN ERROR
          GO TO 100
      ENDIF
      ns = ns + 1
      IF( TYPE(NPARAM) .EQ. 'F' ) THEN
          VALS(NPARAM) = AREAL                                          !  FLOATING POINT VALUES
          IF( names(nparam) .EQ. 'SET' ) THEN
              set(ns) = areal
              IF( ns .EQ. 1 ) GOTO 200
          ENDIF
          GOTO 100
      ENDIF
      IF( type(nparam) .EQ. 'L' ) THEN                                  !  32 BIT INTEGER VALUES
          LVALS(NPARAM) = AREAL
          GO TO 100
      ENDIF
C****
C****   FINISHED A LIST, NOW DO THE ERROR AND VALIDITY CHECKS
C****
 1000 CONTINUE
      IF( icompt .NE. 4 ) THEN                                          ! VMS has real screwy tape numbers!
          IF( iunit .LT. 0 .OR. iunit .GT. 100 ) THEN
              PRINT *,' ***  ERROR  ***  Illegal IUNIT.'
              ierror = ierror + 1
          ENDIF
      ELSE
          IF( iunit .LT. 0 .OR. iunit .GT. 100 ) THEN
              PRINT *,' ***  ERROR  ***  Illegal IUNIT.'
              ierror = ierror + 1
          ENDIF
      ENDIF
      IF( noindex .LT. 0 .OR. noindex .GT. 120 ) THEN
          PRINT *,' ***  ERROR  ***  NO index is too big.', noindex
          ierror = ierror + 1
      ENDIF
      IF( itrindex .LT. 0 .OR. itrindex .GT. 120 ) THEN
          PRINT *,' ***  ERROR  ***  TR index is too big.', itrindex
          ierror = ierror + 1
      ENDIF
      IF(IRUN.EQ.0) GO TO 1040                                          ! don't use the tape unless OK
      IF( iread .EQ. 0 ) GOTO 1040
      IF( ierror .NE. 0 ) GOTO 1040
      IF( FIRST ) THEN
          FIRST=.FALSE.
          ILUN=IUNIT
          IF( LPRINT .EQ. 1 ) PRINT *,' INPUT PARAMTERS:'
          IN=0                                                          !  SET THE AP ADDRESS TO 0
          SSECS=SECS
          JFMT=IFMT
          ISRCF=DECIMF
          nfskp = nfskip
          nrskp = nrskip
          nfile = nfiles
          IF(FRP.EQ.99999.AND.NIS.GT.0) LRP=NIS
          IF(FIS.EQ.99999.AND.NIS.GT.0) LIS=NIS                         ! PRESET LIS TO NIS
          IF( iread .EQ. 0 ) GOTO 1040
          rewindi = rewind
          CALL RDTRC(BUF1,BUF1,BUF1,ISTAT,0)                            ! GET THE FIRST TRACE INTO BUF1
          CALL PODISC(IUNHDR,1,0)                                       ! REWIND THE HEADER FILE
          CALL RDDISC(IUNHDR,ISCR,800,ISTAT)                            ! GET BY THE EBCDIC HEADER (VMS CAN NOT PODISC(1,800,)T POSITION!)
          CALL RDDISC(IUNHDR,ISCR,100,ISTAT)                            ! GET THE BINARY HEADER
      ENDIF
      ISI=SI*1000000.                                                   ! CONVERT TO MICROSECONDS (AS IN THE SEGY HEADER)
      IDELAY=DELAY*1000                                                 ! CONVERT TO MILLISECONDS (AS IN THE SEGY HEADER)
      SR=REAL(IBUF(59))/1000000.                                       !  SET UP THE TRACE HEADER!!!
      BUF1(46)=IBUF(55)/1000.                                           ! THE DELAY IN SECONDS
c      BUF1(47)=IBUF(56)/1000.                                           ! THE START MUTE TIME IN SECONDS
c      BUF1(48)=IBUF(57)/1000.                                           ! THE END MUTE TIME IN SECONDS
      BUF1(49)=SR                                                       ! THE SAMPLE INTERVAL IN SECONDS
      IF(ISI.NE.0) IBUF(59)=ISI                                         ! OVERRIDE THE HEADER
      IF(IDELAY.GT.-32000) IBUF(55)=IDELAY                              ! ALLOW A NEW DELAY OF ZERO!
C****
C****   DO THE PRESETS FROM THE TAPE HEADER
C****
      IF(NTRCS.NE.0) INTRCS=NTRCS                                       ! THIS TELL RDTRC TO EXPECT NTRCS PER SHOT
      IF(NTRGAT.GT.0) ISCR(15)=2                                        ! SET THE BINARY TAPE HEADER SORT FLAG TO RP
      IF( NTRGAT.LT.0) THEN
          ISCR(15)=1
          lbuf1(7) = 0
      ENDIF
      IF( fno .NE. -32767 ) GOTO 1040
      IF( FIS .NE. -32767 .AND. LBUF1(7) .NE. 0 ) THEN                  ! FIS GIVEN ON RP SORTED TAPE
          PRINT *,' ***  WARNING  ***  FIS given on an rp sorted tape.'
          iwarn = iwarn + 1
      ENDIF
      IF( FRP .NE. -32767 .AND. LBUF1(7) .EQ. 0 ) THEN                  ! FRP GIVEN ON SHOT SORTED TAPE
          PRINT *,' ***  ERROR  ***  FRP GIVEN ON A SHOT SORTED TAPE'
          IERROR=IERROR+1
      ENDIF
      IF( FIS .EQ. -32767 .AND. frp .EQ. -32767 ) THEN
          IF( lbuf1(7) .EQ. 0 ) THEN
              fis = 99999                                               ! was lbuf1(3)
          ELSE
              frp = 99999                                               ! was lbuf1(6)
          ENDIF
      ENDIF
c      IF( LTR .EQ. -32767 .AND. alltr .EQ. 1 ) LTR = INTRCS               ! PRESET THE NUMBER OF TRACES PER SHOT
 1040 BUF1(49)=SR
      IF(FDAY+LDAY+FGMT+LGMT.LE.0) GO TO 1050                           ! ARE THE PARAMETERS BY TIME?
      IF(LDAY.EQ.0) LDAY=FDAY                                           !  PRESET LDAY TO FDAY
      IF( LDAY .LT. FDAY ) THEN
          PRINT *,' ***  ERROR  ***  BAD DAY OR GMT PARAMETERS.'
          IERROR=IERROR+1
      ENDIF
      IF(FDAY.NE.0.AND.LDAY.NE.0.AND.FGMT.NE.-1.AND.LGMT.NE.-1)
     *    GO TO 1043
          PRINT *,' ***  ERROR  ***  BAD DAY OR GMT PARAMETERS.'
          IERROR=IERROR+1
 1043 CONTINUE
      TIME1=FDAY*10000.+FGMT                                            ! PACK THE DAY AND GMT INTO 1 WORD
      TIME2=LDAY*10000.+LGMT
      IF( TIME2 .LT. TIME1) THEN
          PRINT *,' ***  ERROR  ***  BAD DAY OR GMT PARAMETERS.'
          IERROR=IERROR+1
      ENDIF
      FGMT=-1
      LGMT=-1
 1050 CONTINUE
      IF(SINC.NE.1.AND.TIME1+TIME2.NE.0.) THEN
           PRINT *,' ***  WARNING  ***  SINC WILL BE IGNORED WHEN ',
     *          'PROCESSING BY GMT.'
           IWARN=IWARN+1
      ENDIF
      IF( order .NE. 1 ) order = -1
      IF( stime .LT. 0. .AND. stime .NE. -1. ) THEN
          PRINT *,' ***  ERROR  ***  STIME must be non-negative.'
          ierror = ierror + 1
      ENDIF
C****
C****
C****
      IF(FNO.EQ.99999.AND.LNO.EQ.32767) LNO=999999
      IF(FIS.EQ.99999.AND.LIS.EQ.32767) LIS=999999
      IF(FIS.EQ.-32767.OR.LIS.EQ.32767) GO TO 1150
      IF(IABS(FIS).LE.LIS.OR.FIS.EQ.99999) GO TO 1150
      PRINT *,' ***  WARNING  ***  FIS SHOULD NORMALLY BE SMALLER THAN',
     *    ' LIS'
      IWARN=IWARN+1
 1150 CONTINUE
      IF(FRP.EQ.99999.AND.LRP.EQ.32767) LRP=999999
      IF(FRP.EQ.-32767.OR.LRP.EQ.32767) GO TO 1160
      IF(IABS(FRP).LE.LRP.OR.FRP.EQ.99999) GO TO 1160
      PRINT *,' ***  WARNING  *** FRP SHOULD NORMALLY BE SMALLER THAN',
     *   ' LRP.'
      IWARN=IWARN+1
 1160 CONTINUE
C****
C****    FINISHED EDITING, NOW WRITE THE LIST TO DISC FOR THE EXECUTION PHASE
C****
      NLISTS=NLISTS+1                                                   ! INCREMENT THE LISTS COUNTER
      IF( fis .NE. -32767 ) THEN
          istart = fis
          iend = lis
          ltype = 'SHOT'
      ENDIF
      IF( frp .NE. -32767 ) THEN
          istart = frp
          iend = lrp
          ltype = 'RP'
      ENDIF
      IF( fno .NE. -32767 ) THEN
          istart = fno
          iend = lno
          ltype = 'NONE'
      ENDIF
      IF( rpinc .NE. 1 ) sinc = rpinc
      IF( noinc .NE. 1 ) sinc = noinc
      IF( sinc .NE. 1 ) allno = 0
      IF( trinc .NE. 1 ) alltr = 0
      rewindi = rewind
      WRITE(IPARUN) IREELN,IUNIT,LTYPE,istart,iend,SINC,FTR,LTR,TRINC,
     *  ISI,IDELAY,NTRGAT,TIME1,TIME2,DECIMF,lprint,forgat,order,stime,
     &  inputdev, set, renum, notype, noindex, itrtype, itrindex,
     &  allno, alltr, rewind
      IF( IAND(LPRINT,1) .EQ. 1 ) THEN
          PRINT *, NLISTS,IREELN,IUNIT,LTYPE,istart,iend,SINC,FTR,LTR,
     *       TRINC,ISI,IDELAY,NTRGAT,TIME1,TIME2,DECIMF,forgat,order
          PRINT *, stime, set, renum, notype, noindex, itrtype,
     &        itrindex, allno, alltr, rewind
	  PRINT *, inputdev
      ENDIF
c****
c****  set the defaults
c****
      rewind = 1
      renum = -32767
 2030 CALL GETOKE(TOKEN,NCHARS)                                         ! GET THE NEXT TOKEN
      CALL UPCASE(TOKEN,NCHARS)
      NTOKES=NTOKES+1
      IF(NCHARS.GT.0) GO TO 2040                                        ! WAS IT THE END OF A LINE?
      IF(NOW.EQ.1) PRINT 140
      CALL RDLINE                                                       ! GET ANOTHER LINE
      NTOKES=0
      GO TO 2030
 2040 IF(TOKEN(1:NCHARS).NE.'END'.OR.NCHARS.NE.3) GO TO 150
C****
C****  PUT THE TAPE HEADER CARD IMAGES IN THE TAPE HEADER FILE
C****
      IF( .NOT. CFIRST ) THEN                                           ! ARE THERE ANY CARDS FOR THE TAPE HEADER?
          IF( IRUN .NE. 0) THEN
              CALL PODISC(IUNHDR,1,0)                                   ! REWIND THE HEADER FILE
              DO 2100 I=1,40
                 CALL RDDISC(IUNHDR,CHEADR(I),20,ISTAT)
                 IF(COMMNT(I).NE.' ') cheadr(i) = COMMNT(I)             ! TRANSLATE AND MOVE IT
 2100         CONTINUE
              CALL PODISC(IUNHDR,1,0)
              DO 2110 I=1,40                                            ! REWRITE THE HEADER TO DISC
 2110         CALL WRDISC(IUNHDR,CHEADR(I),20)
          ENDIF
      ENDIF
      RETURN
      END


text:
! dummy routine
      subroutine mc23ad( )
      end

max_forks_repo_name: dynamics-of-stellar-systems/dynamite


text: #redirect Naomi Amaha

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search