**==aa0001.spg  processed by SPAG 4.50J  at 14:50 on 30 Jun 1995
      IMPLICIT NONE
C*** Start of declarations inserted by SPAG
      REAL a , a1 , a2 , ABIn , ABUnd , ABUnj , AH , alcool , ALF , 
     &     ALFj , ALFly , ALFlyc , alter , AM , AN , beta , BIN , 
     &     BINmin , BINsyz , brem
      REAL BRMev , BZEro , C , C1 , C1J , C2 , C2J , CA , CAJ , CION , 
     &     cionh , CJ , CNC , colcut , collhy , COLumn , CONce , crec , 
     &     CTH , delf
      REAL den , dene , DENlim , DENz , dh , dist , DNE , dtpdti , 
     &     DYNam , dzero , E , e1 , E3 , E3J , EERg , EJ , elapse , 
     &     elohy , energy , ES
      REAL ESJ , EXIt , F , FCOol , FECOOL , FEForb , FEIii , FEIiio , 
     &     FEInf , FEPerm , FJ , foii , FOUt , frac , FTOt , func , 
     &     FWV , h , halc , HALf
      REAL HALfc , hb , HBEt , HBEtc , hbtot , hcolumn , HCOOL , hdot , 
     &     HE , HEA , heat , HENeut , HEPlus , herec , herecs , HEZi , 
     &     HEZii , hgas , HMAx , hnu
      REAL hprime , hrec , hrecs , other , PARt , PB , PCOol , PDUm , 
     &     pheat , PHIhe , PHIon , PHOT , PHOt2 , PHOt3 , PHOt4 , PM , 
     &     POT , POU , POWer , PREss
      REAL PTOt , PWR , rd , RE , RECev , RF , RF1 , RFAct , rh , 
     &     rhtest , RHY , rhyeq , RMAx , S , SI , SIG , SIGj , sigma , 
     &     SJ , SOLVX
      REAL SS , st , sume , t , ta , tap , TAU , TAUhe , TAUP , TBIn , 
     &     temper , tf3 , tf3o , tfb , tfei , tfep , time , TINt , 
     &     TLOw , total
      REAL TPLus , TU , TUFev , TWOp , TZEro , veloc , vm , vs , 
     &     vshock , w , WAVe , WDUm , WJ , x , x0 , xlumin , XMAx
      INTEGER i , ibn , ICOnt , ICX , IHE , ipho , ISTar , it , itf , 
     &        ITFin , izero , j , jcont , JDUm , JMX , k , ki , LDUm , 
     &        LL , LLJ
      INTEGER n , NBIn , NF , NI , nit , NJ , nn , nnn , no , NUM
C*** End of declarations inserted by SPAG
C  MAIN SHOCK PROGRAM : UPDATED MAY 1983
      INCLUDE 'rayspec.inc'
      COMMON /FEL   / WJ(12,220) , FJ(12,220) , E3J(12,220) , 
     &                PWR(12,220) , C1J(12,30) , C2J(12,30) , EJ(12,30)
     &                , SJ(12,30) , CJ(12,30) , LLJ(12,30) , SIGj(12,30)
     &                , ALFj(12,30) , ESJ(12,30)
      COMMON /PARAMS/ NJ(12) , ABUnj(12) , ABUnd , BINmin , BINsyz , 
     &                NBIn
      COMMON /COM   / CNC(12,30) , PTOt(12,220) , ABIn(MXBINS) , 
     &                BIN(MXBINS)
      COMMON /FLN   / FWV(320) , FTOt(320) , FOUt(320) , FCOol
      COMMON /FE    / FEInf , FEPerm , FEForb , FEIii , FEIiio
      COMMON /AUG   / CAJ(12,30)
      COMMON /HOT   / HE , HEA(2)
      COMMON /SHPAR / BZEro , TINt , HEZi , HEZii , TLOw , RF1 , SS , 
     &                SI , RMAx , RFAct , NUM , ITFin , ICOnt , IHE , 
     &                ICX , NI , AH , AN , AM , XMAx , HMAx , DENlim , 
     &                EERg , ISTar
      COMMON /DAT   / E(30) , S(30) , C(30) , WAVe(220) , E3(220) , 
     &                F(220) , LL(30) , SIG(30) , ALF(30) , ES(30)
      COMMON /RESULT/ CONce(30) , CA(30) , POWer(220) , RHY , HENeut , 
     &                HEPlus , DNE , PCOol , POU , POT , RE , TU , PM(4)
      COMMON /RATES / C1(30) , C2(30) , CTH
      COMMON /PT    / RF(500) , TAU(150) , TAUhe(150) , TPLus(150)
      COMMON /ADDER / TBIn(500) , COLumn(10,6) , PHOt2 , PHOt3 , PHOt4
      CHARACTER*80 LABel
      COMMON /LAB   / LABel
      REAL nhm
      INTEGER bb
      COMMON /CONTIN/ BRMev(MXBINS) , RECev(MXBINS) , TUFev(MXBINS)
      DIMENSION rd(10)
      COMMON /SLV   / DENz , PREss , PB , DYNam
      COMMON /TZ    / TZEro , PARt
      COMMON /HLN   / EXIt(11) , HALf , HALfc , ALFly , ALFlyc , HBEt , 
     &                HBEtc , TWOp
      COMMON /PH    / PHIon , PHIhe
      DIMENSION alcool(150) , func(150) , energy(150) , hprime(150) , 
     &          den(150)
      DIMENSION time(150) , dist(150) , h(150) , x0(150) , temper(150) , 
     &          heat(150)
      DIMENSION dh(150) , dene(150) , veloc(150)
      COMMON /LNPR  / JMX(11) , NF(11) , LDUm(5) , JDUm(5) , WDUm(5) , 
     &                PDUm(5)
      DATA NF/0 , 9 , 23 , 42 , 66 , 90 , 123 , 166 , 212 , 271 , 316/
      DATA NJ/2 , 6 , 7 , 8 , 10 , 12 , 14 , 16 , 18 , 20 , 26 , 28/
C
      OPEN (UNIT=15,FILE='[RAYMN.NEW]RDIN.DAT',TYPE='OLD')
      OPEN (UNIT=17,FILE='[RAYMN.NEW]SHPARAM.DAT',TYPE='OLD')
      OPEN (UNIT=14,FILE='[RAYMN.NEW]RDOUT.DAT',TYPE='NEW')
C      OPEN(UNIT=66,FILE='APOOL:[RAYMN]IRON.III',TYPE='NEW')
C      OPEN(UNIT=67,FILE='APOOL:[RAYMN]OXYGEN.III',TYPE='NEW')
C      OPEN(UNIT=68,FILE='APOOL:[RAYMN]NEON.III',TYPE='NEW')
C      OPEN(UNIT=69,FILE='APOOL:[RAYMN]TEMP.DEN',TYPE='NEW')
C      OPEN(UNIT=70,FILE='APOOL:[RAYMN]SULFUR.III',TYPE='NEW')
C      OPEN(UNIT=71,FILE='APOOL:[RAYMN]ARGON.III',TYPE='NEW')
      OPEN (UNIT=77,FILE='[RAYMN.NEW]PARTIAL.DAT',TYPE='NEW')
C
      CALL SETPAR(vshock,colcut)
C
C      WRITE (66,8030) LABEL
C      WRITE (67,8030) LABEL
C      WRITE (68,8030) LABEL
C      WRITE (69,8030) LABEL
C      WRITE (70,8030) LABEL
C      WRITE (71,8030) LABEL
      WRITE (77,99015) LABel
      WRITE (77,99001)
C
C  NBIN SHOULD BE LESS THAN 501 TO FIT TBIN DIMENSION
C
C  ++++++++++++++++++++++   READ ATOMIC DATA
C
      CALL ATREAD(NUM)
C
C  ++++++++++++++++++++++   READ PRESHOCK RADIATION
C
      CALL RADRD(ICOnt)
C
      DO 100 i = 1 , NBIn
         ABIn(i) = ABIn(i)*RF1
         RF(i) = 1.
 100  CONTINUE
C
C     INITIALIZE FIRST STEP
      func(ITFin) = PARt
      frac = PARt
      energy(ITFin) = EERg
      EERg = energy(ITFin)*1.6021E-12
      h(ITFin) = HMAx
      CONce(1) = 1. - PARt
      CONce(2) = PARt
      RHY = PARt/(1.-PARt)
C  INITIAL CALL TO FELINE TO SET IONIZATION FRACTIONS
      dzero = (PARt+.0003)*DENz
      CALL FELINE(TZEro,dzero,NUM,1,0,0,1,1,1)
      IF ( ISTar.EQ.0 ) THEN
         DO 150 i = 1 , NBIn
            ABIn(i) = 0.
 150     CONTINUE
      ENDIF
      jcont = 1
      IF ( IHE.NE.0 ) CNC(1,1) = HEZi
      IF ( IHE.NE.0 ) CNC(1,2) = HEZii
      IF ( IHE.NE.0 ) CNC(1,3) = 1. - HEZi - HEZii
C
C  ++++++++++++++++   BEGIN STEP LOOP
C
      DO 300 nit = 1 , ITFin
         itf = nit
         HE = 0.
         HEA(1) = 0.
         HEA(2) = 0.
         it = ITFin + 1 - nit
         IF ( it.NE.1 ) THEN
            x = SOLVX(frac,EERg,h(it))
            x0(it) = x
            veloc(it) = vshock/x
            t = (PREss-PB*x*x-DYNam/x)/(DENz*1.38E-16*(1.+frac)*x)
            temper(it) = t
            IF ( t.LE.10.**TLOw ) GOTO 400
            st = SQRT(t)
            IF ( it.NE.ITFin ) THEN
               vs = vshock*1.E-5
               vm = vs/x
               nhm = DENz/AH
               DO 160 i = 1 , NBIn
                  w = 12399./(BINmin+(i-.5)*BINsyz)
                  ta = TAUP(w,ITFin)
                  RF(i) = AMIN1(RMAx,ALOG((ta+1.)/ta)-.4/(ta+1.)**2)
 160           CONTINUE
               dtpdti = -(temper(it+1)-t)/elapse
               IF ( dtpdti.GE.1.0E-35 ) hprime(it+1) = -hdot/dtpdti
               RHY = 1.0/(CONce(1)+1.0E-35) - 1.0
               other = temper(it+1)*hprime(it+1)
C     PRINT FLOW OF PREVIOUS STEP
               k = it + 1
               PRINT 99002 , k , temper(k) , dist(k) , time(k) , 
     &               elapse , den(k) , dene(k) , h(k) , other , 
     &               hprime(k) , alcool(k) , xlumin
               PRINT 99003 , RHY
            ENDIF
            dene(it) = func(it)*x*DENz
            func(it-1) = 0.
            energy(it-1) = 0.
            sigma = 6.3E-18
            beta = 157890./t
            cionh = CION(1,1,13.6,t)
            crec = 2.06E-11*(.4288+.5*ALOG(beta)+.469*beta**(-.3333))/st
            PHIon = PHOT(1,1,13.6,t,dene(it))
            pheat = HE
            tap = TAUP(911.,ITFin)
            CALL HEXC(t,dene(it),RHY,tap,cionh,pheat)
C
C     PRINT CURRENT FLOW
            PRINT 99004 , it , t , dene(it) , x , veloc(it)
            PRINT 99005 , PHIon , pheat , HEA
            cionh = cionh + PHIon/dene(it)
            brem = 2.29E-4*st*RHY/(1.+RHY)
            halc = HCOOL(t,RHY,alter)
            hrecs = crec*dene(it)*RHY/(1.+RHY)
            RE = 0.
            TU = 0.
            ipho = 1
            IF ( ISTar.NE.0 ) ipho = 2
            CALL FELINE(t,dene(it),NUM,2,0,jcont,ipho,1,ICX)
            HE = (pheat+HEA(1)+HEA(2))/dene(it)
            heat(it) = HE
            herecs = C2J(1,2)*dene(it)*HEPlus
            elohy = .0001 + RHY
C +++++++++++++++++++++    CALL FORBIDDEN LINE ROUTINE
            CALL FSF(t,dene(it),elohy)
            CALL FPRINT(NUM)
C
C +++++++++++++++++++++    NET COOLING RATE
C
            alcool(it) = brem + halc + PCOol + FCOol + 
     &                   FECOOL(dene(it),t,2) - HE
            total = alcool(it) + HE
            PRINT 99006 , it , total , DNE , halc , FCOol , brem , TU , 
     &            RE
C
C   +++++++++++++++++++++++  BEGIN EVALUATION OF NEXT STEP
C
            IF ( it.NE.ITFin ) time(it) = time(it+1) + e1/3.156E7
            IF ( it.NE.ITFin ) dist(it) = dist(it+1) + e1*veloc(it+1)
            den(it) = x*DENz
            hdot = 1.0E-23*alcool(it)*dene(it)/AN
            xlumin = hdot*den(it)
C  USE ENTHALPY OF GAS ONLY TO CALACULATE STEP LENGTH
            hgas = .5*DYNam/(DENz*x*x) + 2.5*(1.+den(it)/dene(it))
     &             *t*1.38E-16 + EERg
            elapse = ABS((1.-10.**(-TINt))*hgas/hdot)
            IF ( ITFin-it.LE.NI ) elapse = 1.E8*SI/DENz
            e1 = elapse
            IF ( PHIon.GT.0. .AND. t.LT.20000. ) THEN
               delf = AMAX1(ABS(func(it+1)-func(it))/func(it),.01)
               e1 = SS/(PHIon*delf)
               IF ( ISTar.EQ.2 ) e1 = SS/(dene(it)*crec)
               IF ( e1.LT.elapse ) THEN
                  rhyeq = cionh/crec
                  rhtest = ABS((RHY-rhyeq)/RHY)
                  IF ( ISTar.EQ.2 .AND. rhtest.LT.0.05 )
     &                 e1 = colcut/(DENz*vshock)
                  elapse = e1
                  PRINT 99007
               ENDIF
            ENDIF
            e1 = elapse
            h(it-1) = h(it) - hdot*elapse
            dh(it) = hdot*elapse
            hrec = hrec + hrecs*elapse
            herec = herec + herecs*elapse
            HENeut = CNC(1,1)
            HEPlus = CNC(1,2)
C
C  ++++++++++++++++++++   HYDROGEN OUTPUTS
C
            rh = RHY/(1.+RHY)
            a = .5E-23*vshock*elapse*DENz*dene(it)/AN
            CALL HLINE(RHY,a,t,st,beta,alter)
            ibn = (10.2-BINmin)/BINsyz + 1
            BIN(ibn) = BIN(ibn) + EXIt(6) + EXIt(9)
            CALL TWOPH(1215.,EXIt(11),0)
C
C  ADD UP OPTICAL DEPTHS, COLUMN DENSITIES AND RADIATION
C
            a1 = DENz*vshock*elapse
            a2 = a1*dene(it)*0.5E-23/AN
            CALL ADDER(a1,a2,it,ITFin,NUM,DENz,vshock)
            IF ( ISTar.EQ.2 ) CALL LINOP(NUM,t,veloc(it),a1,it,1)
C
            collhy = (COLumn(1,1)+COLumn(1,2)+COLumn(1,3))/AH
            IF ( collhy.GT.colcut ) GOTO 400
            IF ( ISTar.EQ.2 .AND. TAU(ITFin).GE.0.1 ) GOTO 400
C      WRITE (66,421) IT,CNC(11,3)
C      WRITE (67,421) IT,CNC(4,3)
C      WRITE (68,421) IT,CNC(5,3)
C      WRITE (70,421) IT,CNC(8,3)
C      WRITE (71,421) IT,CNC(9,3)
C      WRITE (69,422) IT,TEMPER(IT),DEN(IT),DENE(IT),DIST(IT)
C
C  CALCULATE IONIZATION BALANCE FOR NEXT STEP
C
            DO 200 no = 1 , NUM
               n = NJ(no)
               nn = n + 1
               a = 10.**(ABUnj(no)-12.)/AN
               sume = 0.
               DO 170 j = 1 , nn
                  CONce(j) = CNC(no,j)
                  IF ( CONce(j).GT..0001 ) JMX(no) = MAX0(JMX(no),j)
                  CA(j) = CAJ(no,j)
                  C1(j) = C1J(no,j)
                  C2(j) = C2J(no,j)
 170           CONTINUE
               nnn = nn + 1
               DO 180 j = nnn , 30
                  CA(j) = 0.
                  C1(j) = 0.
                  C2(j) = 0.
                  CONce(j) = 0.
 180           CONTINUE
               CALL CYCLE(n,dene(it),elapse)
               CNC(no,1) = CONce(1)
               DO 190 j = 2 , nn
                  CNC(no,j) = CONce(j)
                  func(it-1) = func(it-1) + a*(j-1)*CONce(j)
                  sume = sume + EJ(no,j-1)
                  energy(it-1) = energy(it-1) + a*sume*CONce(j)
 190           CONTINUE
 200        CONTINUE
            CONce(1) = 1./(1.+RHY)
            CONce(2) = RHY*CONce(1)
            C1(1) = cionh
            C2(2) = crec
            CALL CYCLE(1,dene(it),elapse)
            func(it-1) = func(it-1) + CONce(2)/AN
            energy(it-1) = energy(it-1) + 13.595*CONce(2)/AN
            EERg = energy(it-1)*1.602E-12
            frac = func(it-1)
C
C  ADD UP IRON
            tfei = tfei + FEInf*a2
            tfep = tfep + FEPerm*a2
            tfb = tfb + FEForb*a2
            tf3 = tf3 + FEIii*a2
            tf3o = tf3o + FEIiio*a2
C
            IF ( (it/10)*10.EQ.it ) CALL LNPRT(NUM)
            RHY = 1./CONce(1) - 1.0
            hcolumn = (COLumn(1,1)+COLumn(1,2)+COLumn(1,3)+1.0)/AH
            foii = FTOt(31) + FTOt(32)
            hbtot = HBEt + HBEtc
            WRITE (77,99008)  , it , hcolumn , hbtot , t , dene(it) , 
     &                       PTOt(4,18) , FTOt(26) , foii , FTOt(36)
         ENDIF
 300  CONTINUE
 400  WRITE (77,99008) izero
      PHOt2 = 0.
      ibn = (13.6-BINmin)/BINsyz + 1
      DO 500 i = ibn , NBIn
         hnu = BINmin + (i-.5)*BINsyz
         PHOt2 = PHOt2 + TBIn(i)/(hnu*1.6E-12*DENz*vshock)
 500  CONTINUE
      PRINT 99011 , hrec , herec , TINt , PHOt2 , PHOt3 , PHOt4
      PRINT 99012 , HBEt , HBEtc
      PRINT 99013 , HALf , HALfc
      PRINT 99014 , ALFly , ALFlyc
      PRINT 99016 , ((COLumn(1,k)),k=1,6) , ((COLumn(2,k)),k=1,6)
      PRINT 99016 , ((COLumn(3,k)),k=1,6) , ((COLumn(4,k)),k=1,6)
      PRINT 99016 , ((COLumn(5,k)),k=1,6) , ((COLumn(6,k)),k=1,6)
      PRINT 99016 , ((COLumn(7,k)),k=1,6) , ((COLumn(8,k)),k=1,6)
      PRINT 99016 , ((COLumn(9,k)),k=1,6)
      PRINT 99017
      DO 600 nit = 1 , itf
         it = ITFin - nit + 1
         PRINT 99020 , it , TAU(it) , TAUhe(it) , TPLus(it) , energy(it)
     &         , den(it) , veloc(it) , dh(it) , func(it)
 600  CONTINUE
      PRINT 99019
      WRITE (35,99015) LABel
      WRITE (35,99018)
      DO 700 nit = 1 , itf
         it = ITFin - nit + 1
         PRINT 99020 , it , time(it) , dist(it) , x0(it) , temper(it) , 
     &         dene(it) , alcool(it) , h(it) , hprime(it) , heat(it)
         WRITE (35,99020) it , time(it) , dist(it) , x0(it) , temper(it)
     &                    , dene(it) , alcool(it) , h(it) , heat(it) , 
     &                    energy(it) , den(it) , veloc(it)
 700  CONTINUE
      PRINT 99021
      nn = NBIn/10
      WRITE (14,99022) vshock , DENz , ICOnt
      DO 800 i = 1 , nn
         hnu = BINmin + (10*(i-1))*BINsyz
         DO 750 k = 1 , 10
            rd(k) = TBIn(10*i-10+k)
 750     CONTINUE
         vs = vshock*1.E-7
         PRINT 99023 , i , hnu , rd
         WRITE (14,99023) i , hnu , ((rd(k)),k=1,10)
 800  CONTINUE
      CALL LNPRT(NUM)
      hb = (HBEt+HBEtc)/100.
      tfei = tfei/hb
      tfep = tfep/hb
      tfb = tfb/hb
      tf3 = tf3/hb
      tf3o = tf3o/hb
      PRINT 99024 , tfei , tfep , tfb , tf3 , tf3o
      IF ( ISTar.EQ.2 ) CALL LINOP(NUM,0.,0.,0.,ki,2)
      STOP
99001 FORMAT (
     &'  IT   HCOLUMN     HBETA        T       DENE     O III]      [O I
     &]      [O II]    [O III]')
99002 FORMAT ('0',I3,'  T=',E10.2,' D=',E10.2,' TI=',E10.2,'  E1=',
     &        E10.2,' DEN=',E10.2,' DENE=',E10.2,/,'  H=',E10.2,
     &        ' OTHER=',E10.2,' HPR=',E10.2,' ALCOOL=',E10.2,' XLUM=',
     &        E10.2)
99003 FORMAT ('   +++   RHY = ',E10.2)
99004 FORMAT ('1',' STEP NO ',I3,' T=',E10.2,'  DENE=',E10.2,'  X=',
     &        E10.2,'  VEL =',E10.2)
99005 FORMAT ('  PHION,PHEAT ',8E10.2)
99006 FORMAT ('0',' STEP ',I3,' TOTAL',E10.2,'  DNE',E10.2,' HALC',
     &        E10.2,'  FC',E10.2,/'  BREM',E10.2,'  TUPH',E10.2,'  REC',
     &        E10.2)
99007 FORMAT (' *****  SHORT STEP  *****')
99008 FORMAT (I5,11E11.3)
C
C  WRITE OUT O III, NE III, FE III AND T,D,DENE
C
99009 FORMAT (4X,I2,4X,1PE9.3)
99010 FORMAT (4X,I2,4(4X,1PE9.3))
99011 FORMAT ('1',' HREC = ',F10.2,'  HEREC = ',F10.2,'  TINT=',F10.3,
     &        '  PHOT2=',E10.3,' PHOT3=',F8.3,' PHOT4=',F8.3)
99012 FORMAT (' HBET = ',E10.3,' HBETC = ',E10.3,
     &        ' ERGS CM-2 S-1 OUT FRONT')
99013 FORMAT (1X,'HALF = ',E10.3,' HALFC = ',E10.3)
99014 FORMAT ('  LYMAN ALPHA RECOMB ',E10.3,' COLLISIONAL ',E10.3)
99015 FORMAT (A80)
99016 FORMAT ('  COLUMN',6E9.2,3X,6E9.2)
99017 FORMAT (
     &' STEP    TAU       TAUHE     TAUPLUS     ENERGY      DEN       VE
     &LOC       DH       FUNC')
99018 FORMAT (
     &'STEP   TIME       DIST       X0         TEMP        DENE       AL
     &COOL      H        HEAT      ENERGY       DEN       VEL')
99019 FORMAT ('1',
     &' STEP   TIME       DIST       XO         TEMP       DENE      ALC
     &OOL       H        HPRIME      HEAT')
99020 FORMAT (1X,I4,11E11.4)
99021 FORMAT ('1','  ++  TBINS  ++')
99022 FORMAT ('   V = ',E10.2,'   DENZ=',E10.2,' ICONT=',I3)
99023 FORMAT (1X,I3,F9.1,1P10E10.2)
99024 FORMAT (
     &  ' IRON LINES:  [FE IR] ,   FE II ,  [FE II] ,FE III,FEIII   OPT'
     &  ,/,12X,10F10.2)
      END
