C                                                                       
C     Program for spectral-synthesis calculation
C     (developed by Takeda based on WIDTH9 program)
C     NSET>0 ---> IFVAC=0   (wavelength in the atmosphere)
C     NSET<0 ---> IFVAC=1   (vacuum wavelength)
C
      PROGRAM KAPREPX       
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (kw=99)
      PARAMETER(KPNT=10001,KCN=100)
      EXTERNAL BLOCKE,BLOCKR,BLOCKT
      EXTERNAL FSET
      COMMON /ABROSS/ABROSS(kw),TAUROS(kw)                              
      COMMON /ABTOT/ABTOT(kw),ALPHA(kw)                                 
      COMMON /CONT/ABTOTC(kw),ALPHAC(kw),TAUNUC(kw),SNUC(kw),HNUC(kw),          
     1             JNUC(kw),JMINSC(kw),RESIDC(kw)                               
      REAL*8 JNUC,JMINSC                                                        
      COMMON /CONV/DLTDLP(kw),HEATCP(kw),DLRDLT(kw),VELSND(kw),         
     1             GRDADB(kw),HSCALE(kw),FLXCNV(kw),VCONV(kw),MIXLTH,   
     2             IFCONV,OVERWT,FLXCNV0(kw),FLXCNV1(kw)                                               
      REAL*8 MIXLTH                                                       
      COMMON /DEPART/BHYD(kw,6),BMIN(kw),NLTEON                         
      COMMON /ELEM/ABUND(99),ATMASS(99),ELEM(99)                        
      character*2 elem
      COMMON /EDENS/EDENS(kw),IFEDNS
      COMMON /FILENAME/FILENAME                                                 
      CHARACTER*60 FILENAME                                                     
      COMMON /FLUX/ FLUX,FLXERR(kw),FLXDRV(kw),FLXRAD(kw)               
      COMMON /FREQ/FREQ,FREQLG,EHVKT(kw),STIM(kw),BNU(kw)               
      COMMON /FRESET/FRESET(1563),RCOSET(1563),NULO,NUHI,NUMNU            
      COMMON /HEIGHT/HEIGHT(kw)                                         
      COMMON /IF/IFCORR,IFPRES,IFSURF,IFSCAT,TAUSCAT,IFMOL                      
      COMMON /IFOP/IFOP(20)                                             
      COMMON /IONS/XNFPH(kw,2),XNFPHE(kw,3),XNFH(kw,2),XNFHE(kw,3)                             
      COMMON /ITER/ ITER,IFPRNT(15),IFPNCH(15),NUMITS                   
      COMMON /JUNK/TITLE(74),FREQID(6),WLTE,XSCALE                      
      COMMON /MUS/ANGLE(20),SURFI(20),NMU                               
      COMMON /OPS/AHYD(kw),AH2P(kw),AHMIN(kw),SIGH(kw),AHE1(kw),        
     1       AHE2(kw),AHEMIN(kw),SIGHE(kw),ACOOL(kw),ALUKE(kw),AHOT(kw),
     2            SIGEL(kw),SIGH2(kw),AHLINE(kw),ALINES(kw),SIGLIN(kw), 
     3            AXLINE(kw),SIGXL(kw),AXCONT(kw),SIGX(kw),SHYD(kw),    
     4            SHMIN(kw),SHLINE(kw),SXLINE(kw),SXCONT(kw)            
      COMMON /OPTOT/ACONT(kw),SCONT(kw),ALINE(kw),SLINE(kw),SIGMAC(kw), 
     1              SIGMAL(kw)                                          
      COMMON /PTOTAL/PTOTAL(kw)
      COMMON /PUT/PUT,IPUT
      COMMON /PZERO/PZERO,PCON,PRADK0,PTURB0,KNU(kw),PRADK(kw),RADEN(kw)
      REAL*8 KNU
      COMMON /RAD/ ACCRAD(kw),PRAD(kw)
      COMMON /RHOX/RHOX(kw),NRHOX                                       
      COMMON /STATE/P(kw),XNE(kw),XNATOM(kw),RHO(kw)                    
      COMMON /TAUSHJ/TAUNU(kw),SNU(kw),HNU(kw),JNU(kw),JMINS(kw)        
      REAL*8 JNU,JMINS                                                    
      COMMON /STEPLG/STEPLG,TAU1LG,KRHOX
      COMMON /TEFF/TEFF,GRAV,GLOG
      COMMON /TEMP/T(kw),TKEV(kw),TK(kw),HKT(kw),TLOG(kw),ITEMP         
      COMMON /TURBPR/VTURB(kw),PTURB(kw),TRBFDG,TRBCON,TRBPOW,TRBSND,   
     1               IFTURB                                             
      COMMON /WAVEY/WBEGIN,DELTAW,IFWAVE                                
      COMMON /XABUND/XABUND(99),WTMOLE                                  
      COMMON /XNF/XNFC(kw,6),XNFN(kw,6),XNFO(kw,6),XNFNE(kw,6),                 
     1            XNFMG(kw,6),XNFSI(kw,6),XNFS(kw,6),XNFFE(kw,5)                
      COMMON /XNFP/XNFPC(kw,4),XNFPN(kw,5),XNFPO(kw,6),XNFPNE(kw,6),            
     1             XNFPAL(kw,1),XNFPMG(kw,2),XNFPSI(kw,2),XNFPCA(kw,2),         
     2             XNFPFE(kw,1),XNFPCH(kw),XNFPOH(kw)                           
      DIMENSION PART(kw,6)                                                      
C                                                                       
      COMMON /CONTIN/CONTIN                                             
      COMMON /CURVE/MINLOG,DABLOG,NABLOG,LOCONV
      LOGICAL LOCONV
      REAL*8 MINLOG                                                       
      COMMON /GAM/GLOGR,GLOGS,GLOGW                                     
      COMMON /IFPROF/IFPROF                                             
      COMMON /LINDAT/WL,GFLOG,XJ,E,XJP,EP,CODE,LABEL,LABELP,NELION,     
     1               GAMMAR,GAMMAS,GAMMAW,REF,NBLO,NBUP,ISO1,X1,ISO2,   
     2               X2,OTHER1,OTHER2,ELO,GF                            
      COMMON /LINEY/LINE0(3,kw),DOPWAV(3,kw),ADAMP(3,kw)                
      REAL*8 LINE0                                                        
      COMMON /OBS/EW,WLOBS,STAR(9)                                      
      COMMON /RESIDU/RESIDU(15),TAUONE(15),DELWAV(15)                   
      COMMON /VTS/VTS(3),NVT                                            
      COMMON /XLABEL/LAB,EXPOT                                          
      COMMON /ABDIF/ABCOR                                               
      COMMON /PLOTCG/MODECG,PDATA(9,3,2)                                
      COMMON /WDAMP/TXNXN(kw)                                           

      DIMENSION POTION(594)                                             
      DIMENSION POTA(114),POTB(114),POTC(114),POTD(114),POTE(114)       
      DIMENSION POTF(24)                                                
      DIMENSION XD(300),YD(300)                                         
      EQUIVALENCE (POTION(1),POTA(1)),(POTION(115),POTB(1))             
      EQUIVALENCE (POTION(229),POTC(1)),(POTION(343),POTD(1))           
      EQUIVALENCE (POTION(457),POTE(1)),(POTION(571),POTF(1))           
      CHARACTER LAB*28,STAR*8                                           
      CHARACTER*1 GCHAR(4),CHAR(300)                                    
      CHARACTER*80 ABLAB(2)                                             
      DATA GCHAR/'#','@','&','$'/                                       
      DIMENSION NREAD(100),STWV(100),ENWV(100),NLIN1(100),NLIN2(100)    
      DIMENSION PABCO(kw,KPNT),WVCON(KCN),CONFLX(KCN)                   
      DIMENSION WLTH(KPNT),AINT(KPNT),MMP(100),NNCON(KCN),RESIDF(KPNT)  
      DIMENSION KCSV(KPNT),ACONSV(kw,KCN),ETA0(kw)                      
      DIMENSION TAU5(kw),PRHOX(kw),PT(kw),PXNE(kw),                     
     #PTXNXN(kw),VTSAVE(3)
c     #PTXNXN(kw),VTSAVE(3),LABSYN(20)                                   
      CHARACTER*4 LABSYN(20)
      DIMENSION AHSAVE(kw,KCN),SAVEAL(kw)                               
      DIMENSION SV1(kw,KCN),SV2(kw,KCN),SV3(kw,KCN),SV4(kw,KCN)
      DIMENSION SV5(kw,KCN),SV6(kw,KCN)
      DIMENSION ICD(20),SABCO(kw,KPNT)
      LOGICAL LOSEP
      REAL*8 STWV,ENWV,SSTWV,EENWV,WLTH,DELTA,WVCON,DELC,WLDP,DMW       
      DATA QLCAL,QSYNT/4HLCAL,4HSYNT/                                   
      DIMENSION CARD(20)                                                
      DATA QSTOP,QEND,QCURV,QAVER,QVTUR,QLINE,QPROF/4HSTOP,4HEND ,      
     1 4HCURV,4HAVER,4HVTUR,4HLINE,4HPROF/                              
      data qmole/'MOLE'/
      DATA POTA/                                                        
     1   13.598,   0.   ,   0.   ,   0.   ,   0.   ,   0.   ,           
     2   24.587,  54.416,   0.   ,   0.   ,   0.   ,   0.   ,           
     3    5.392,  75.638, 122.451,   0.   ,   0.   ,   0.   ,           
     4    9.322,  18.211, 153.893, 217.713,   0.   ,   0.   ,           
     5    8.298,  25.154,  37.930, 259.368, 340.217,   0.   ,           
     6   11.260,  24.383,  47.887,  64.492, 392.077, 489.981,           
     7   14.534,  29.601,  47.448,  77.472,  97.888, 552.057,           
     8   13.618,  35.116,  54.934,  77.412, 113.896, 138.116,           
     9   17.422,  34.970,  62.707,  87.138, 114.240, 157.161,           
     T   21.564,  40.962,  63.45 ,  97.11 , 126.21 , 157.93 ,           
     1    5.139,  47.286,  71.64 ,  98.91 , 138.39 , 172.15 ,           
     2    7.646,  15.035,  80.143, 109.24 , 141.26 , 186.50 ,           
     3    5.986,  18.828,  28.447, 119.99 , 153.71 , 190.47 ,           
     4    8.151,  16.345,  33.492,  45.141, 166.77 , 205.05 ,           
     5   10.486,  19.725,  30.18 ,  51.37 ,  65.023, 220.43 ,           
     6   10.36 ,  23.33 ,  34.83 ,  47.30 ,  72.68 ,  88.049,           
     7   12.967,  23.81 ,  39.61 ,  53.46 ,  67.80 ,  97.03 ,           
     8   15.759,  27.629,  40.74 ,  59.81 ,  75.02 ,  91.007,           
     9    4.341,  31.625,  45.72 ,  60.91 ,  82.66 , 100.0  /           
      DATA POTB/                                                        
     1    6.113,  11.871,  50.908,  67.10 ,  84.41 , 108.78 ,           
     2    6.54 ,  12.80 ,  24.76 ,  73.47 ,  91.66 , 111.1  ,           
     3    6.82 ,  13.58 ,  27.491,  43.266,  99.22 , 119.36 ,           
     4    6.74 ,  14.65 ,  29.310,  46.707,  65.23 , 128.12 ,           
     5    6.766,  16.50 ,  30.96 ,  49.1  ,  69.3  ,  90.56 ,           
     6    7.435,  15.640,  33.667,  51.2  ,  72.4  ,  95.   ,           
     7    7.870,  16.18 ,  30.651,  54.8  ,  75.0  ,  99.0  ,           
     8    7.86 ,  17.06 ,  33.50 ,  51.3  ,  79.5  , 102.0  ,           
     9    7.635,  18.168,  35.17 ,  54.9  ,  75.5  , 108.0  ,           
     T    7.726,  20.292,  36.83 ,  55.2  ,  79.9  , 103.0  ,           
     1    9.394,  17.964,  39.722,  59.4  ,  82.6  , 108.0  ,           
     2    5.999,  20.51 ,  30.71 ,  64.0  ,   0.   ,   0.   ,           
     3    7.899,  15.934,  34.22 ,  45.71 ,  93.5  ,   0.   ,           
     4    9.81 ,  18.633,  28.351,  50.13 ,  62.63 , 127.6  ,           
     5    9.752,  21.19 ,  30.820,  42.944,  68.3  ,  81.70 ,           
     6   11.814,  21.8  ,  36.   ,  47.3  ,  59.7  ,  88.6  ,           
     7   13.999,  24.359,  36.95 ,  52.5  ,  64.7  ,  78.5  ,           
     8    4.177,  27.28 ,  40.0  ,  52.6  ,  71.0  ,  84.4  ,           
     9    5.695,  11.030,  43.6  ,  57.0  ,  71.6  ,  90.8  /           
      DATA POTC/                                                        
     1    6.38 ,  12.24 ,  20.52 ,  61.8  ,  77.0  ,  93.0  ,           
     2    6.84 ,  13.13 ,  22.99 ,  34.34 ,  81.50 ,   0.   ,           
     3    6.88 ,  14.32 ,  25.04 ,  38.3  ,  50.55 , 102.6  ,           
     4    7.099,  16.15 ,  27.16 ,  46.4  ,  61.2  ,  68.0  ,           
     5    7.28 ,  15.26 ,  29.54 ,   0.   ,   0.   ,   0.   ,           
     6    7.37 ,  16.76 ,  28.47 ,   0.   ,   0.   ,   0.   ,           
     7    7.46 ,  18.08 ,  31.06 ,   0.   ,   0.   ,   0.   ,           
     8    8.34 ,  19.43 ,  32.93 ,   0.   ,   0.   ,   0.   ,           
     9    7.576,  21.49 ,  34.83 ,   0.   ,   0.   ,   0.   ,           
     T    8.993,  16.908,  34.48 ,   0.   ,   0.   ,   0.   ,           
     1    5.786,  18.869,  28.03 ,  54.0  ,   0.   ,   0.   ,           
     2    7.344,  14.632,  30.502,  40.734,  72.28 ,   0.   ,           
     3    8.641,  16.53 ,  25.3  ,  44.2  ,  56.0  , 108.0  ,           
     4    9.009,  18.6  ,  27.96 ,  37.41 ,  58.75 ,  70.7  ,           
     5   10.451,  19.131,  33.0  ,   0.   ,   0.   ,   0.   ,           
     6   12.130,  21.21 ,  32.1  ,   0.   ,   0.   ,   0.   ,           
     7    3.894,  25.1  ,  35.0  ,   0.   ,   0.   ,   0.   ,           
     8    5.212,  10.004,  37.000,   0.   ,   0.   ,   0.   ,           
     9    5.577,  11.06 ,  19.177,  49.954,   0.   ,   0.   /           
      DATA POTD/                                                        
     1    5.47 ,  10.85 ,  20.197,  36.758,   0.   ,   0.   ,           
     2    5.42 ,  10.55 ,  21.624,  38.981,   0.   ,   0.   ,           
     3    5.49 ,  10.72 ,  22.14 ,  40.42 ,   0.   ,   0.   ,           
     4    5.55 ,  10.90 ,  22.42 ,  41.09 ,   0.   ,   0.   ,           
     5    5.63 ,  11.07 ,  23.45 ,  41.47 ,   0.   ,   0.   ,           
     6    5.67 ,  11.25 ,  24.71 ,  42.65 ,   0.   ,   0.   ,           
     7    6.14 ,  12.1  ,  20.38 ,  44.03 ,   0.   ,   0.   ,           
     8    5.85 ,  11.52 ,  21.98 ,  39.84 ,   0.   ,   0.   ,           
     9    5.93 ,  11.67 ,  22.83 ,  41.56 ,   0.   ,   0.   ,           
     T    6.02 ,  11.80 ,  22.84 ,  42.51 ,   0.   ,   0.   ,           
     1    6.10 ,  11.93 ,  22.74 ,  42.66 ,   0.   ,   0.   ,           
     2    6.18 ,  12.05 ,  23.68 ,  42.69 ,   0.   ,   0.   ,           
     3    6.254,  12.17 ,  25.03 ,  43.74 ,   0.   ,   0.   ,           
     4    5.426,  13.9  ,  20.960,  45.193,   0.   ,   0.   ,           
     5    7.0  ,  14.9  ,  23.3  ,  33.319,   0.   ,   0.   ,           
     6    7.89 ,  16.200,  24.0  ,   0.   ,   0.   ,   0.   ,           
     7    7.98 ,  17.70 ,  25.0  ,   0.   ,   0.   ,   0.   ,           
     8    7.88 ,  16.6  ,  26.0  ,   0.   ,   0.   ,   0.   ,           
     9    8.7  ,  17.0  ,  27.0  ,   0.   ,   0.   ,   0.   /           
      DATA POTE/                                                        
     1    9.1  ,  20.0  ,  28.0  ,   0.   ,   0.   ,   0.   ,           
     2    9.0  ,  18.563,  29.0  ,   0.   ,   0.   ,   0.   ,           
     3    9.225,  20.5  ,  30.0  ,   0.   ,   0.   ,   0.   ,           
     4   10.437,  18.756,  34.2  ,   0.   ,   0.   ,   0.   ,           
     5    6.108,  20.428,  29.83 ,   0.   ,   0.   ,   0.   ,           
     6    7.416,  15.032,  31.937,  42.32 ,  68.8  ,   0.   ,           
     7    7.289,  16.69 ,  25.56 ,  45.3  ,  56.0  ,  88.3  ,           
     8    8.42 ,  19.0  ,  27.0  ,   0.   ,   0.   ,   0.   ,           
     9    9.3  ,  20.0  ,  30.   ,   0.   ,   0.   ,   0.   ,           
     T   10.748,  20.0  ,  30.0  ,   0.   ,   0.   ,   0.   ,           
     1    4.0  ,  22.0  ,  33.0  ,   0.   ,   0.   ,   0.   ,           
     2    5.279,  10.147,  34.0  ,   0.   ,   0.   ,   0.   ,           
     3    6.9  ,  12.1  ,  20.0  ,   0.   ,   0.   ,   0.   ,           
     4    6.0  ,  11.5  ,  20.0  ,  28.8  ,   0.   ,   0.   ,           
     5    6.0  ,  12.0  ,  20.0  ,   0.   ,   0.   ,   0.   ,           
     6    6.0  ,  12.0  ,  20.0  ,   0.   ,   0.   ,   0.   ,           
     7    6.0  ,  12.0  ,  20.0  ,   0.   ,   0.   ,   0.   ,           
     8    5.800,  12.0  ,  20.0  ,   0.   ,   0.   ,   0.   ,           
     9    6.0  ,  12.0  ,  20.0  ,   0.   ,   0.   ,   0.   /           
      DATA POTF/                                                        
     1    6.0  ,  12.0  ,  20.0  ,   0.   ,   0.   ,   0.   ,           
     2    6.0  ,  12.0  ,  20.0  ,   0.   ,   0.   ,   0.   ,           
     3    6.0  ,  12.0  ,  20.0  ,   0.   ,   0.   ,   0.   ,           
     4    6.0  ,  12.0  ,  20.0  ,   0.   ,   0.   ,   0.   /
c
      character*10 molcode
      character*2 elm
      dimension iza(8)
c
      character*64 fnopen

      iter=1
      ifprnt(1)=0
      ifpnch(1)=0
      numits=1

c      write(*,*) 'Input the name of the file including fileopen data' 
c      read(*,'(A)') fnopen
c      open (unit=1,file=fnopen)
      open (unit=1,file='kapfiles.lst')
      rewind 1
      call fileopen
      ABLAB(1)='X: LOG (EPS) ;Y: LOG (W) <IN MA>'                       
      ABLAB(2)='X: LOG (EPS*GF) ;Y: LOG (W/L)'                          
      MODECG=0                                                          
      IFLCAL=0                                                          
      IFSYNT=0                                                          
      IFPROF=0                                                          
      ITEMP=0                                                           
      NABLOG=0                                                          
 2222 REWIND 13                                                         
      VTS(1)=0.                                                         
      VTS(2)=0.                                                         
      VTS(3)=0.                                                         
      NVT=1                                                             
      IFSURF=1                                                          
      MODEAV=1                                                          
c**** mole ****b
      ifmol=0
      ifpres=0
c**** mole ****e
    1 READ(5,2)CARD                                                     
      IF(CARD(1).EQ.QSTOP)CALL EXIT                                     
      WRITE(13,2)CARD                                                   
    2 FORMAT(20A4)                                                      
c**** mole ****b
      if(card(1).eq.'MOLE') then
      ifmol=1
      ifpres=1
      rewind 2
      call readmol
      end if
c**** mole ****e
      IF(CARD(1).NE.QEND)GO TO 1                                        
      REWIND 13                                                         
    3 CALL READIN(20)                                                   
C     IFOP(14) VS. IFOP(15)                                             
      IF(IFOP(15).EQ.1) THEN                                            
       IF(IFOP(14).EQ.1) THEN                                           
       WRITE(6,*)'WHY BOTH OF IFOP(14) AND IFOP(15) EQUAL TO UNITY ???' 
       CALL EXIT                                                        
       END IF                                                           
      IFOP(15)=0                                                        
      IFOP(14)=1                                                        
      WRITE(6,*)                                                        
      WRITE(6,*) '$ $ $ CHANGED TO IFOP(14)-->1 AND IFOP(15)-->0 $ $ $' 
      WRITE(6,*)                                                        
      END IF                                                            
c      IFPRES=0                                                          
      IF(NRHOX.EQ.0)GO TO 2222                                          
      ABCOR=DLOG10(ABUND(1))                                            
      ITEMP=ITEMP+1                                                     
      CALL INTEG(RHOX,ABROSS,TAUROS,NRHOX,ABROSS(1)*RHOX(1))            
      IF(IFPRES.NE.0)CALL POPS(0.D0,1,XNE)                                      
      CALL POPS(1.01D0,11,XNFPH)                                                
      CALL POPS(2.02D0,11,XNFPHE)                                               
      CALL POPS(6.03D0,11,XNFPC)                                                
      CALL POPS(7.04D0,11,XNFPN)                                                
      CALL POPS(8.05D0,11,XNFPO)                                                
      CALL POPS(10.05D0,11,XNFPNE)                                              
      CALL POPS(12.01D0,11,XNFPMG)                                              
      CALL POPS(13.00D0,11,XNFPAL)                                              
      CALL POPS(14.01D0,11,XNFPSI)                                              
      CALL POPS(20.01D0,11,XNFPCA)                                              
      CALL POPS(26.00D0,11,XNFPFE)                                              
      IF(IFMOL.EQ.0)THEN                                                        
      CALL POPS(1.01D0,12,XNFH)                                                 
      CALL POPS(2.02D0,12,XNFHE)                                                
      CALL POPS(6.05D0,12,XNFC)                                                 
      CALL POPS(7.05D0,12,XNFN)                                                 
      CALL POPS(8.05D0,12,XNFO)                                                 
      CALL POPS(10.05D0,12,XNFNE)                                               
      CALL POPS(12.05D0,12,XNFMG)                                               
      CALL POPS(14.05D0,12,XNFSI)                                               
      CALL POPS(16.05D0,12,XNFS)                                                
      CALL POPS(26.04D0,12,XNFFE)                                               
      ENDIF                                                                     
      IF(IFMOL.EQ.1)THEN                                                        
      CALL POPS(106.00D0,11,XNFPCH)                                             
      CALL POPS(108.00D0,11,XNFPOH)                                             
      CALL W(6HXNFPCH,XNFPCH,NRHOX)                                             
      CALL W(6HXNFPOH,XNFPOH,NRHOX)                                             
C     THE POPS WILL NOT RETURN NUMBER DENSITIES WHEN MOLECULES ARE ON           
C     SO WE COMPUTE NUMBER DENSITIES/PART FUNCTIONS  AND PART FUNCTIONS         
      CALL POPS(6.05D0,11,XNFC)                                                 
      CALL POPS(7.05D0,11,XNFN)                                                 
      CALL POPS(8.05D0,11,XNFO)                                                 
      CALL POPS(10.05D0,11,XNFNE)                                               
      CALL POPS(12.05D0,11,XNFMG)                                               
      CALL POPS(14.05D0,11,XNFSI)                                               
      CALL POPS(16.05D0,11,XNFS)                                                
      CALL POPS(26.04D0,11,XNFFE)                                               
      DO 444 J=1,NRHOX                                                          
      CALL PFSAHA(J,1,1,3,PART)                                                 
      XNFH(J,1)=XNFPH(J,1)*PART(J,1)                                            
      XNFH(J,2)=XNFPH(J,2)                                                      
      CALL PFSAHA(J,2,2,13,PART)                                                
      XNFHE(J,1)=XNFPHE(J,1)*PART(J,1)                                          
      XNFHE(J,2)=XNFPHE(J,2)*PART(J,2)                                          
      XNFHE(J,3)=XNFPHE(J,3)                                                    
      CALL PFSAHA(J,6,6,13,PART)                                                
      XNFC(J,1)=XNFC(J,1)*PART(J,1)                                             
      XNFC(J,2)=XNFC(J,2)*PART(J,2)                                             
      XNFC(J,3)=XNFC(J,3)*PART(J,3)                                             
      XNFC(J,4)=XNFC(J,4)*PART(J,4)                                             
      XNFC(J,5)=XNFC(J,5)*PART(J,5)                                             
      XNFC(J,6)=XNFC(J,6)*PART(J,6)                                             
      CALL PFSAHA(J,7,6,13,PART)                                                
      XNFN(J,1)=XNFN(J,1)*PART(J,1)                                             
      XNFN(J,2)=XNFN(J,2)*PART(J,2)                                             
      XNFN(J,3)=XNFN(J,3)*PART(J,3)                                             
      XNFN(J,4)=XNFN(J,4)*PART(J,4)                                             
      XNFN(J,5)=XNFN(J,5)*PART(J,5)                                             
      XNFN(J,6)=XNFN(J,6)*PART(J,6)                                             
      CALL PFSAHA(J,8,6,13,PART)                                                
      XNFO(J,1)=XNFO(J,1)*PART(J,1)                                             
      XNFO(J,2)=XNFO(J,2)*PART(J,2)                                             
      XNFO(J,3)=XNFO(J,3)*PART(J,3)                                             
      XNFO(J,4)=XNFO(J,4)*PART(J,4)                                             
      XNFO(J,5)=XNFO(J,5)*PART(J,5)                                             
      XNFO(J,6)=XNFO(J,6)*PART(J,6)                                             
      CALL PFSAHA(J,10,6,13,PART)                                               
      XNFNE(J,1)=XNFNE(J,1)*PART(J,1)                                           
      XNFNE(J,2)=XNFNE(J,2)*PART(J,2)                                           
      XNFNE(J,3)=XNFNE(J,3)*PART(J,3)                                           
      XNFNE(J,4)=XNFNE(J,4)*PART(J,4)                                           
      XNFNE(J,5)=XNFNE(J,5)*PART(J,5)                                           
      XNFNE(J,6)=XNFNE(J,6)*PART(J,6)                                           
      CALL PFSAHA(J,12,6,13,PART)                                               
      XNFMG(J,1)=XNFMG(J,1)*PART(J,1)                                           
      XNFMG(J,2)=XNFMG(J,2)*PART(J,2)                                           
      XNFMG(J,3)=XNFMG(J,3)*PART(J,3)                                           
      XNFMG(J,4)=XNFMG(J,4)*PART(J,4)                                           
      XNFMG(J,5)=XNFMG(J,5)*PART(J,5)                                           
      XNFMG(J,6)=XNFMG(J,6)*PART(J,6)                                           
      CALL PFSAHA(J,14,6,13,PART)                                               
      XNFSI(J,1)=XNFSI(J,1)*PART(J,1)                                           
      XNFSI(J,2)=XNFSI(J,2)*PART(J,2)                                           
      XNFSI(J,3)=XNFSI(J,3)*PART(J,3)                                           
      XNFSI(J,4)=XNFSI(J,4)*PART(J,4)                                           
      XNFSI(J,5)=XNFSI(J,5)*PART(J,5)                                           
      XNFSI(J,6)=XNFSI(J,6)*PART(J,6)                                           
      CALL PFSAHA(J,16,6,13,PART)                                               
      XNFS(J,1)=XNFS(J,1)*PART(J,1)                                             
      XNFS(J,2)=XNFS(J,2)*PART(J,2)                                             
      XNFS(J,3)=XNFS(J,3)*PART(J,3)                                             
      XNFS(J,4)=XNFS(J,4)*PART(J,4)                                             
      XNFS(J,5)=XNFS(J,5)*PART(J,5)                                             
      XNFS(J,6)=XNFS(J,6)*PART(J,6)                                             
      CALL PFSAHA(J,26,5,13,PART)                                               
      XNFFE(J,1)=XNFFE(J,1)*PART(J,1)                                           
      XNFFE(J,2)=XNFFE(J,2)*PART(J,2)                                           
      XNFFE(J,3)=XNFFE(J,3)*PART(J,3)                                           
      XNFFE(J,4)=XNFFE(J,4)*PART(J,4)                                           
      XNFFE(J,5)=XNFFE(J,5)*PART(J,5)                                           
  444 CONTINUE                                                                  
      ENDIF                                                                     
      WL=500.0                                                          
      FREQ=2.997925E17/WL                                               
      FREQLG=DLOG(FREQ)                                                 
      DO 33 J=1,NRHOX                                                   
      EHVKT(J)=EXP(-FREQ*HKT(J))                                        
      STIM(J)=1.0-EHVKT(J)                                              
      BNU(J)=1.47439D-47*FREQ**3*EHVKT(J)/STIM(J)                       
   33 CONTINUE                                                          
      N=1
      CALL KAPP(N,NSTEPS,STEPWT)
      DO 34 J=1,NRHOX                                                   
   34 ACONT(J)=ACONT(J)+SIGMAC(J)                                       
      CALL INTEG(RHOX,ACONT,TAU5,NRHOX,ACONT(1)*RHOX(1))                
      WRITE(6,333)                                                      
  333 FORMAT(1H1)                                                       
      IF(TRBSND.EQ.0.)GO TO 6                                           
      WRITE(6,7)                                                        
    7 FORMAT(40HREWRITE WIDTH5 TO INCLUDE SOUND VELOCITY)               
      CALL EXIT                                                         
    6 DO 5 J=1,NRHOX                                                    
      VELSND(J)=0.                                                      
    5 VTURB(J)=(TRBFDG*RHO(J)**TRBPOW+TRBSND*VELSND(J)/1.E5+TRBCON)*1.E5
    8 READ(13,9)XLABEL,NMLINE,STAR                                      
    9 FORMAT(A4,I4,9A8)                                                 
C     CURV(E OF GROWTH)                                                 
      IF(XLABEL.EQ.QCURV)GO TO 170                                      
C     AVER(AGE)                                                         
      IF(XLABEL.EQ.QAVER)GO TO 20                                       
C     VTUR(B)                                                           
      IF(XLABEL.EQ.QVTUR)GO TO 40                                       
C     END                                                               
      IF(XLABEL.EQ.QEND)GO TO 50                                        
C     LINE                                                              
      IF(XLABEL.EQ.QLINE)GO TO 60                                       
C     PROF(ILE)                                                         
      IF(XLABEL.EQ.QPROF)GO TO 80                                       
C     LCAL                                                              
      IF(XLABEL.EQ.QLCAL) GO TO 700                                     
C     SYNT(HESIS)                                                       
      IF(XLABEL.EQ.QSYNT) GO TO 900                                     
C     MOLE
      if(xlabel.eq.qmole) then
      write(6,*)
      write(6,'('' <<< molecules formation to be included >>>'')')
      write(6,'(''     ifmol = '',i2,4x,''ifpres = '',i2)') ifmol,ifpres
      write(6,*)
      go to 8
      end if
      WRITE(6,11)XLABEL                                                 
      CALL EXIT                                                         
   11 FORMAT(9H WHAT IS A4)                                             
   80 IFPROF=1                                                          
      GO TO 8                                                           
   20 IF(MODEAV.EQ.2)GO TO 30                                           
      MODEAV=2                                                          
      GO TO 8                                                           
   30 MODEAV=3                                                          
      CALL AVERAG(MODEAV,IVT,ABLG,TAULG)                                
      GO TO 8                                                           
   40 READ(13,41)NVT,(VTS(IVT),IVT=1,NVT)                               
   41 FORMAT(I5,3F5.2)                                                  
      GO TO 8                                                           
  170 READ(13,171)NABLOG,MINLOG,DABLOG,MODECG                           
C     CONVERT LOG(N/NH)+12 TO LOG(N/NT)                                 
      MINLOG=MINLOG-12.0+ABCOR                                          
  171 FORMAT(I5,2F8.2,I5)                                               
      GO TO 8                                                           
  700 IFLCAL=1                                                          
      WRITE(6,333)                                                      
      WRITE(6,1400)                                                     
      WRITE(6,736)                                                      
  736 FORMAT(1H0,25X,'%%%%%%%%%%  LCAL-MODE  %%%%%%%%%%'/)              
      DO 708 I=1,3                                                      
      VTSAVE(I)=VTS(I)                                                  
      VTS(I)=0.0                                                        
  708 CONTINUE                                                          
      NVTSAV=NVT                                                        
      NVT=1                                                             
      READ(13,720) NMSET,NUML1                                          
  720 FORMAT(2I5)                                                       
      NUML=NUML1                                                        
      IF(NUML.LE.0) NUML=1                                              
      NUML=NUML-1                                                       
      READ(13,903) (NREAD(I),I=1,NMSET)                                 
c!      ifvac=0
c!      do 4901 i=1,nmset
c!      if(nread(i).le.0) then
c!      nread(i)=iabs(nread(i))
c!      ifvac=1
c!      end if
c! 4901 continue
c!      write(6,4902) ifvac
c! 4902 format(1h ,'(((((( IFVAC=',i2,' ))))))')
      ISET=0
      REWIND 9                                                          
  740 CONTINUE                                                          
      IF(ISET.NE.0) WRITE(6,732) ISET,NUML,WL                           
  732 FORMAT(1H0,15X,'SET= ',I2,10X,'IND= ',I6,10X,'WLEND= ',F9.4/)     
      wlnend=wl
      WRITE(6,1400)                                                     
      ISET=ISET+1                                                       
      IF(ISET.LE.NMSET) GO TO 745                                       
      IFLCAL=0                                                          
      REWIND 9                                                          
      DO 789 I=1,3                                                      
  789 VTS(I)=VTSAVE(I)                                                  
      NVT=NVTSAV                                                        
      GO TO 8                                                           
  745 CONTINUE
CGFALL ==================== START ===========================
c!      READ(9,7760) WL,GFLOG,CODE,EXLOW,XJLOW,EXUP,XJUP,molcode
c! 7760 FORMAT(F11.4,F7.3,F6.2,2(F12.3,F5.1,11X),a10)
c!      if(code.le.0.) read(molcode,*) code

      read(9,*) code,wla,expot,gflog,gammar,gammas,gammaw
      wl=wla/10.

c!      if(ifvac.eq.1) then
c!      wltemp=wl
c!      wl=1.d7/dabs(exup-exlow)
c!c      write(6,1767) wltemp,wl
c!      end if
CGFALL ======================END=============================
      MMLL=NUML+1                                                       
      WRITE(6,733) ISET,MMLL,WL                                         
  733 FORMAT(1H0,15X,'SET= ',I2,10X,'IND= ',I6,10X,'WLBGN= ',F9.4/)     
      wlnbgn=wl
      BACKSPACE 9
      NCOUNT=0                                                          
  770 CONTINUE                                                          
      IF(NCOUNT.GE.NREAD(ISET)) GO TO 740                               
      NCOUNT=NCOUNT+1                                                   
      NUML=NUML+1                                                       
c!      GAMMAR=0.0                                                        
c!      GAMMAS=0.0                                                        
c!      GAMMAW=0.0                                                        
CGFALL ===================== START ===========================
c!      READ(9,7760) WL,GFLOG,CODE,EXLOW,XJLOW,EXUP,XJUP,molcode
c!      if(code.le.0.) read(molcode,*) code

      read(9,*) code,wla,expot,gflog,gammar,gammas,gammaw
      wl=wla/10.
      
c!      if(ifvac.eq.1) then
c!      wltemp=wl
c!      wl=1.d7/dabs(exup-exlow)
c!      write(6,1767) wltemp,wl
c! 1767 format(1h ,'wl value of ',f11.4,' has been changed to ',f11.4,
c!     % ' because of ifvac=1')
c!      end if
      IJEL=INT(CODE+0.5)
      KLEL=INT((CODE-FLOAT(IJEL))*100.+0.5)
C     READ(9,760) WL,GFLOG,EXLOW,EXUP,IJEL,KLEL                         
C 760 FORMAT(F10.4,F7.3,2F12.3,12X,I2,I1,24X)                           
C     CODE=FLOAT(IJEL*100+KLEL)/100.0                                   
CGFALL ===================== END ============================
c      EXPOT=min(EXLOW,EXUP)*1.23981E-4                                
      GO TO 769                                                         
   50 REWIND 13                                                         
      CLOSE(UNIT=40)
      WRITE(6,810)                                                      
  810 FORMAT(1H1,8X,'J',10X,'RHOX',15X,'T',13X,'TAUROSS',10X,           
     #'LGTAUR',9X,'TAU(5000)',8X,'LG(TXNXN)'/)                          
      DO 812 J=1,NRHOX                                                  
      ATAUR=DLOG10(TAUROS(J))                                           
      ATXN=DLOG10(TXNXN(J))                                             
      WRITE(6,813)J,RHOX(J),T(J),TAUROS(J),ATAUR,TAU5(J),ATXN           
  813 FORMAT(1H ,7X,I2,7X,1PE11.3,7X,0PF8.1,7X,1PE11.3,7X,0PF8.3,       
     #7X,1PE11.3,7X,0PF8.3)                                             
  812 CONTINUE                                                          
      GO TO 3                                                           
   60 CONTINUE                                                          
      NCOUNT=0                                                          
  230 CONTINUE                                                          
      READ(13,446) CODE,WL,EXPOT,GFLOG,GAMMAR,GAMMAS,GAMMAW,EW,LAB      
  446 FORMAT(F6.2,F9.3,5F6.3,F7.1,A28)                                  
      IF(MODECG.EQ.0) GO TO 178                                         
      READ(13,176) XMN,XMX,YMN,YMX                                      
  176 FORMAT(4F8.2)                                                     
  178 CONTINUE                                                          
      NCOUNT=NCOUNT+1                                                   
      WRITE(6,449)                                                      
  449 FORMAT(1H0,'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@',  
     1'@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@'/)  
      WRITE(6,447) CODE,WL,EXPOT,GFLOG,GAMMAR,GAMMAS,GAMMAW,EW,LAB      
  447 FORMAT(1H0,10X,'INPUT LINE DATA:',F6.2,F9.3,5F6.3,F7.1,A28/)      
      WRITE(6,1400)                                                     
 1400 FORMAT(//)                                                        
C     CONVERT A TO NM AND MA TO PM                                      
      WL=WL/10.0                                                        
      EW=EW/10.0                                                        
  769 CONTINUE                                                          
C     ESTIMATE NELION                                                   
      MDUM=CODE*100.0+0.5                                               
      MION=MOD(MDUM,100)                                                
      MELM=(MDUM-MION)/100                                              
      NELION=6*(MELM-1)+MION+1                                          
C     WRITE(6,*) 'NELION=',NELION                                       
      GF=EXP(GFLOG*2.30258509299405E0)                                  
      EPUP=EXPOT+1.23981E3/WL                                           
      ELO=EXPOT/1.23981E-4                                              
      E=ELO                                                             
      EP=EPUP/1.23981E-4                                                
      GAMMAR=EXP(GAMMAR*2.30258509299405E0)                             
      GAMMAS=EXP(GAMMAS*2.30258509299405E0)                             
      GAMMAW=EXP(GAMMAW*2.30258509299405E0)                             
C     CLASSICAL DAMPING CONSTANT                                        
      IF(GAMMAR.EQ.1.)GAMMAR=2.223E13/WL**2                             
      IF(GAMMAS.NE.1..AND.GAMMAW.NE.1.)GO TO 141                        
      IF(GAMMAS.NE.1.)GO TO 138                                         
      IF(CODE.GE.100.)GO TO 137                                         
c      ZEFF=(CODE-FLOAT(IFIX(CODE)))*100.+1.                             
      ZEFF=(CODE-DFLOAT(INT(CODE)))*100.+1.                             
      EFFNSQ=25.                                                        
      DEPUP=POTION(NELION)-EPUP                                         
      IF(DEPUP.GT.0.)EFFNSQ=13.595*ZEFF**2/DEPUP                        
      GAMMAS=1.0E-8*EFFNSQ**2*SQRT(EFFNSQ)                              
      GO TO 138                                                         
  137 GAMMAS=1.0E-5                                                     
  138 IF(GAMMAW.NE.1.)GO TO 141                                         
      IF(CODE.GE.100.)GO TO 139                                         
c      ZEFF=(CODE-FLOAT(IFIX(CODE)))*100.+1.                             
      ZEFF=(CODE-DFLOAT(INT(CODE)))*100.+1.                             
      EFFNSQ=25.                                                        
      DEPUP=POTION(NELION)-EPUP                                         
      IF(DEPUP.GT.0.)EFFNSQ=13.595*ZEFF**2/DEPUP                        
      RSQ=2.5*(EFFNSQ/ZEFF)**2                                          
      NSEQ=CODE-ZEFF+1.5                                                
C     WRITE(6,*) 'NSEQ=',NSEQ                                           
      IF(NSEQ.GT.20.AND.NSEQ.LT.29)RSQ=(45.-FLOAT(NSEQ))/ZEFF           
      GAMMAW=4.5E-9*RSQ**.4                                             
      GO TO 141                                                         
  139 GAMMAW=1.E-7/ZEFF                                                 
  141 GLOGR=DLOG10(GAMMAR)                                              
C     CALL W(6HGLOGR ,GLOGR,1)                                          
      GLOGS=DLOG10(GAMMAS)                                              
C     CALL W(6HGLOGS ,GLOGS,1)                                          
      GLOGW=DLOG10(GAMMAW)                                              
C     CALL W(6HGLOGW ,GLOGW,1)                                          
      IF(IFSYNT.NE.1) GO TO 928                                         
  900 CONTINUE                                                          
      IFSYNT=1                                                          
      READ(13,904) (LABSYN(I),I=1,20)                                   
      READ(13,901) NMSET
      if(nmset.ne.1) then
      write(6,*) 'Only nmset=1 is allowed in this version!'
      call exit
      end if
      READ(13,902) (STWV(I),ENWV(I),I=1,NMSET)                          
      READ(13,903) (MMP(I),I=1,NMSET)                                   
      READ(13,903) (NNCON(I),I=1,NMSET)                                 
      READ(13,903) (NLIN1(I),NLIN2(I),I=1,NMSET)                        
  901 FORMAT(I5)                                                        
  902 FORMAT(8F10.3)                                                    
  903 FORMAT(8I10)                                                      
  904 FORMAT(20A4)                                                      
C     NOW SPECIFY ELEMENTS TO BE SEPARATELY TREATED                     
      IF(NVT.NE.1) THEN                                                 
      WRITE(6,*) 'NVT SHOULD BE 1 IN THIS VERSION]'                     
      CALL EXIT                                                         
      END IF                                                            
      READ(13,860) NELV                                                 
      READ(13,860) (ICD(L),L=1,NELV)                                    
  860 FORMAT(20I4)                                                      
      REWIND 30                                                         
      WRITE(30,610) (LABSYN(I),I=1,20)                                  
      WRITE(30,6808) TEFF,GLOG,XSCALE,VTS(1)                             
 6808 FORMAT(F8.1,F7.3,F9.5,F6.2)                                       
      WRITE(30,6840) (DLOG10(XABUND(IE)),IE=1,99)                        
 6840 FORMAT(20F8.3)                                                    
      WRITE(30,871) (ATMASS(IE),IE=1,99)                                
  871 FORMAT(10F12.3)                                                   
      WRITE(30,872) NRHOX,'NRHOX'                                       
  872 FORMAT(I4,1X,A8)                                                  
      WRITE(30,870) (RHOX(J),J=1,NRHOX)                                 
  870 FORMAT(1P10E12.5)                                                 
      WRITE(30,873) (T(J),J=1,NRHOX)                                    
  873 FORMAT(10F12.1)                                                   
      WRITE(30,872) NELV,'NELV'                                         
      WRITE(30,860) (ICD(L),L=1,NELV)                                   
c      WRITE(30) (LABSYN(I),I=1,20)                                  
c      WRITE(30) TEFF,GLOG,XSCALE,VTS(1)                             
c      WRITE(30) (DLOG10(XABUND(IE)),IE=1,99)                        
c      WRITE(30) (ATMASS(IE),IE=1,99)                                
c      WRITE(30) NRHOX                                       
c      WRITE(30) (RHOX(J),J=1,NRHOX)                                 
c      WRITE(30) (T(J),J=1,NRHOX)                                    
c      WRITE(30) NELV                                         
c      WRITE(30) (ICD(L),L=1,NELV)                                   
C                                                                       
      REWIND 14                                                         
      WRITE(14,610) (LABSYN(I),I=1,20)                                  
  610 FORMAT(20A4)
      WRITE(14,620) NMSET,NVT
  620 FORMAT(2I5)
      WRITE(6,333)                                                      
      WRITE(6,905)(LABSYN(I),I=1,20)                                    
  905 FORMAT(1H0,20X,'SYNTHESIS LABEL: ',20A4/)                         
      ETACRT=1.00E-5                                                    
      TAUREF=0.20                                                       
      ISET=0                                                            
  909 CONTINUE                                                          
      ISET=ISET+1                                                       
      WRITE(6,449)                                                      
      WRITE(6,916) ISET,NVT,(VTS(IVT),IVT=1,NVT)                        
  916 FORMAT(1H0,15X,'SET= ',I2,15X,'NVT= ',I1,10X,'VTURB=',3F6.1/)     
      WRITE(6,969) STWV(ISET),ENWV(ISET)                                
  969 FORMAT(1H0,20X,'WV1=',F10.3,10X,'WV2=',F10.3//)                   
      SSTWV=STWV(ISET)/10.0                                             
      EENWV=ENWV(ISET)/10.0
      NPOINT=MMP(ISET)                                                  
      NCON=NNCON(ISET)                                                  
      MIDCON=(1+NCON)/2                                                 
      DELC=(EENWV-SSTWV)/FLOAT(NCON)                                    
      DO 930 L=1,NCON                                                   
  930 WVCON(L)=SSTWV+DELC/2.0+FLOAT(L-1)*DELC                           
      LCONT=0                                                           
  932 CONTINUE                                                          
      LCONT=LCONT+1                                                     
      WL=WVCON(LCONT)                                                   
  928 CONTINUE                                                          
      FREQ=2.997925E17/WL                                               
      FREQLG=DLOG(FREQ)                                                 
      DO 61 J=1,NRHOX                                                   
      EHVKT(J)=EXP(-FREQ*HKT(J))                                        
      STIM(J)=1.-EHVKT(J)                                               
      BNU(J)=1.47439D-47*FREQ**3*EHVKT(J)/STIM(J)                       
   61 CONTINUE                                                          
      IF(IFLCAL.NE.1) GO TO 709                                         
      CALL LINCEN(ABUND1,VTS,NVT)                                       
      IND=NUML                                                          
      WRITE(40) IND,CODE,IJEL,KLEL,WL,EXPOT,GFLOG,GLOGR,GLOGS,GLOGW,    
     %(LINE0(1,J),DOPWAV(1,J),ADAMP(1,J),J=1,NRHOX)                     
C     CALL TIME(ITM)
C     WRITE(6,*) 'LINE IND',IND,'  TIME=',ITM
      GO TO 770                                                         
  709 CONTINUE                                                          
C     CALCULATES CONTINUOUS OPACITIES                                   
      N=1
      CALL KAPP(N,NSTEPS,STEPWT)
      IF(IFSYNT.NE.1) GO TO 931                                         
      DO 991 J=1,NRHOX                                                  
      SAVEAL(J)=ALINE(J)                                                
      ALINE(J)=0.0                                                      
  991 CONTINUE                                                          
  931 CONTINUE                                                          
C     CALCULATES THE SOURCE FUNCTION AND SURFACE FLUX OR SURFACE INTENSI
      CALL JOSH(IFSCAT,IFSURF)                                          
      IF(IFSURF.LT.2)CONTIN=HNU(1)                                      
      IF(IFSURF.EQ.2)CONTIN=SURFI(1)                                    
C     CALL W(6HCONTIN,CONTIN,1)                                         
C     TREAT HYDROGEN LINES AS CONTINUUM                                 
      IF(IFSCAT.EQ.0.)GO TO 162                                         
      DO 63 J=1,NRHOX                                                   
      SCONT(J)=(ACONT(J)*SCONT(J)+ALINE(J)*SLINE(J))/(ACONT(J)+ALINE(J))
      ACONT(J)=ACONT(J)+ALINE(J)                                        
      SLINE(J)=BNU(J)                                                   
   63 CONTINUE                                                          
      GO TO 164                                                         
C     SOURCE FUNCTION WILL BE WEIGHTED MEAN OF CONTINUUM SOURCE         
C        FUNCTION (INCLUDING JNU) AND BNU                               
  162 DO 163 J=1,NRHOX                                                  
      SCONT(J)=SNU(J)                                                   
      SLINE(J)=BNU(J)                                                   
  163 ACONT(J)=ACONT(J)+ALINE(J)                                        
  164 CONTINUE                                                          
      IF(IFSYNT.NE.1) GO TO 988                                         
      DO 917 J=1,NRHOX                                                  
      AHSAVE(J,LCONT)=SAVEAL(J)                                         
  917 ACONSV(J,LCONT)=ACONT(J)                                          
      IF(LCONT.NE.MIDCON) GO TO 919                                     
      CALL INTEG(RHOX,ACONT,TAUNU,NRHOX,ACONT(1)*RHOX(1))               
      DO 918 J=1,NRHOX                                                  
      IF(TAUNU(J).GT.TAUREF) GO TO 973                                  
  918 CONTINUE                                                          
  973 JREF=J                                                            
      WRITE(6,3456) JREF,TAU5(JREF),DLOG10(TAU5(JREF))
 3456 FORMAT(1H ,3X,'JREF=',I3,3X,'TAU5=',1PE11.4,3X,'LOG TAU5=',0PF8.4)
  919 CONTINUE                                                          
      CONFLX(LCONT)=CONTIN                                              
      JJ=LCONT                                                          
      DO 915 J=1,NRHOX
      SV1(J,JJ)=EHVKT(J)
      SV2(J,JJ)=STIM(J)
      SV3(J,JJ)=BNU(J)
      SV4(J,JJ)=SCONT(J)
      SV5(J,JJ)=ACONT(J)
  915 SV6(J,JJ)=SIGMAC(J)
      IF(LCONT.LT.NCON) GO TO 932                                       
      DELTA=(EENWV-SSTWV)/FLOAT(NPOINT-1)                               
      DO 910 K=1,NPOINT                                                 
      WLTH(K)=SSTWV+FLOAT(K-1)*DELTA                                    
      KCSV(K)=(WLTH(K)-SSTWV)/DELC+1.000D0                              
  910 CONTINUE                                                          
C                                                                       
      WRITE(30,872) NPOINT,'NPOINT'                                     
      WRITE(30,872) NCON,'NCON'                                         
      WRITE(30,6871) SSTWV,EENWV                                         
 6871 format(2f12.4)
c      WRITE(30) NPOINT                                     
c      WRITE(30) NCON                                         
c      WRITE(30) SSTWV,EENWV                                         
      NSEP=0                                                            
C                                                                       
      KCSV(1)=1                                                         
      KCSV(NPOINT)=NCON                                                 
      N1=NLIN1(ISET)                                                    
      N2=NLIN2(ISET)                                                    
      DO 998 IVT=1,NVT                                                  
      DO 926 K=1,NPOINT                                                 
      DO 927 J=1,NRHOX                                                  
      SABCO(J,K)=0.0
  927 PABCO(J,K)=0.0                                                    
  926 CONTINUE                                                          
      REWIND 40
      DO 920 IL=1,NREAD(ISET)                                                   
      READ(40) IND,CODE,IJEL,KLEL,WL,EXPOT,GFLOG,GLOGR,GLOGS,GLOGW,     
     %(LINE0(1,J),DOPWAV(1,J),ADAMP(1,J),J=1,NRHOX)                     
      IF(IND.LT.N1.OR.IND.GT.N2) GO TO 920
C
      ijm=invcode(code)
      if(code.lt.100..and.ijm.ne.ijel) then
        write(*,*) 'ijm=',ijm,'  ijel=',ijel
        call exit
      end if
      LOSEP=.FALSE.                                                     
      DO 880 L=1,NELV                                                   
      IF(ijm.EQ.ICD(L)) THEN                                           
      LOSEP=.TRUE.                                                      
      GO TO 881                                                         
      END IF                                                            
  880 CONTINUE                                                          
      GO TO 832                                                         
  881 NSEP=NSEP+1                                                       
      write(6,9875) nsep,code,ijm,klel,wl,expot,gflog
 9875 format(1h ,'#### nsep=',i4,' #### ',f8.2,2i4,f10.4,f7.2,f7.3)
      WRITE(30,875) NSEP
  875 FORMAT('####',I4,'  NSEP  ####')                                  
      WRITE(30,6874) ijm,KLEL,WL,EXPOT,GFLOG                            
 6874 FORMAT(2I5,F12.4,F9.3,F8.3)                                        
      WRITE(30,6870) (BNU(J),J=1,NRHOX)                                  
      WRITE(30,6870) (LINE0(1,J),J=1,NRHOX)                              
      WRITE(30,6870) (DOPWAV(1,J),J=1,NRHOX)                             
      WRITE(30,6870) (ADAMP(1,J),J=1,NRHOX)
 6870 format(10e24.16)                              
c      WRITE(30) NSEP                                                
c      WRITE(30) ijm,KLEL,WL,EXPOT,GFLOG                            
c      WRITE(30) (BNU(J),J=1,NRHOX)                                  
c      WRITE(30) (LINE0(1,J),J=1,NRHOX)                              
c      WRITE(30) (DOPWAV(1,J),J=1,NRHOX)                             
c      WRITE(30) (ADAMP(1,J),J=1,NRHOX)                              
  832 CONTINUE                                                          
C                                                                       
      WLDP=WL                                                           
      KCT=(WLDP-SSTWV)/DELC+1.000D0                                     
      KCT=MAX0(1,KCT)                                                   
      KCT=MIN0(NCON,KCT)                                                
      if(code.lt.100.) then
        elm=elem(ijel)
        xabmul=xabund(ijel)
        tatmas=atmass(ijel)
      else
        call calzam(ijel,elm,tatmas,mola,iza)
        xabmul=1.
        do i=1,mola
          xabmul=xabmul*xabund(iza(i))
        end do
      end if
      DO 922 J=1,NRHOX                                                  
      ETA0(J)=xabmul*LINE0(1,J)/ACONSV(J,MIDCON)                  
      IF(ETA0(J).LT.ETACRT) GO TO 922                                   
      DOPCOR=SQRT(1.0+60.1391*tatmas*VTS(IVT)**2/T(J))            
      ETA1=ETA0(J)/DOPCOR                                               
      DMW=1.0E-3*DOPWAV(1,J)*SQRT(ETA0(J)*ADAMP(1,J)/ETACRT*.56419)     
      MK1=(WLDP-SSTWV-DMW)/DELTA                                        
      MK2=(WLDP-SSTWV+DMW)/DELTA                                        
      KMN=MAX0(MK1,1)                                                   
      KMN=MIN0(KMN,NPOINT)                                              
      KMX=MIN0(MK2+1,NPOINT)                                            
      KMX=MAX0(KMX,1)                                                   
      DO 923 K=KMN,KMX                                                  
      V=ABS(WLTH(K)-WLDP)*1.0E3/(DOPWAV(1,J)*DOPCOR)                    
      DUM=xabmul*(LINE0(1,J)/DOPCOR)*VOIGT(V,ADAMP(1,J)/DOPCOR)   
      PABCO(J,K)=PABCO(J,K)+DUM                                         
      IF(LOSEP) GO TO 923
      SABCO(J,K)=SABCO(J,K)+DUM
  923 CONTINUE                                                          
  922 CONTINUE                                                          
      IION=KLEL+1                                                       
      IF(IVT.GT.1) GO TO 920                                            
      WRITE(15,941) WLDP*10.,code,elm,IION,ETA0(JREF),EXPOT,GFLOG,        
     #GLOGR,GLOGS,GLOGW                                                 
  941 FORMAT(F9.3,2X,f5.2,2x,A2,I1,2X,0Pf10.5,2X,0P5F9.3)                   
  920 CONTINUE                                                          
C                                                                       
      WRITE(30,876)                                                     
  876 FORMAT('@@@@ CONT DATA FOLLOW @@@@')                              
      WRITE(30,872) NCON,'NCON'                                         
c***** if nsel to be read is -999 then end of line data *****
      ieflg=-999
c      write(30) ieflg
c      WRITE(30) NCON
      DO 890 II=1,NCON                                                  
      JJ=II
      DO 891 J=1,NRHOX
      EHVKT(J)=SV1(J,JJ)
      STIM(J)=SV2(J,JJ)
      BNU(J)=SV3(J,JJ)
      SCONT(J)=SV4(J,JJ)
      ACONT(J)=SV5(J,JJ)
  891 SIGMAC(J)=SV6(J,JJ)
      WRITE(30,872) JJ,'JJ'                                             
      WRITE(30,6891) (EHVKT(J),STIM(J),BNU(J),SCONT(J),ACONT(J),         
     #SIGMAC(J),J=1,NRHOX)                                              
 6891 format(6e24.16)
c      WRITE(30) JJ                                             
c      WRITE(30) (EHVKT(J),STIM(J),BNU(J),SCONT(J),ACONT(J),         
c     #SIGMAC(J),J=1,NRHOX)                                              
  890 CONTINUE                                                          
      WRITE(30,878)                                                     
  878 FORMAT('@@@@ LINE OPACITY DATA (SABCO) FOLLOW @@@@')              
C                                                                       
      LCONT=0                                                           
      DO 953 K=1,NPOINT                                                 
      LLAB=KCSV(K)                                                      
      IF(LLAB.EQ.LCONT) GO TO 955                                       
      LCONT=LLAB                                                        
      IF(LCONT.GT.NCON) LCONT=NCON                                      
      JJ=LCONT                                                          
      DO 914 J=1,NRHOX
      EHVKT(J)=SV1(J,JJ)
      STIM(J)=SV2(J,JJ)
      BNU(J)=SV3(J,JJ)
      SCONT(J)=SV4(J,JJ)
      ACONT(J)=SV5(J,JJ)
  914 SIGMAC(J)=SV6(J,JJ)
      DO 954 J=1,NRHOX                                                  
  954 SLINE(J)=BNU(J)                                                   
  955 CONTINUE                                                          
      DO 956 J=1,NRHOX                                                  
  956 ALINE(J)=PABCO(J,K)                                               
      IF(IFOP(14).EQ.0) GO TO 938                                       
      LCT=LCONT                                                         
      IF(LCT.EQ.NCON) LCT=NCON-1                                        
      DO 939 J=1,NRHOX                                                  
      DUM=AHSAVE(J,LCT)+(WLTH(K)-WVCON(LCT))*             
     #(AHSAVE(J,LCT+1)-AHSAVE(J,LCT))/(WVCON(LCT+1)-WVCON(LCT))         
      ALINE(J)=ALINE(J)+DUM
      SABCO(J,K)=SABCO(J,K)+DUM
  939 CONTINUE                                                          
  938 CONTINUE                                                          
C                                                                       
      WRITE(30,6877) K,WLTH(K)                                           
 6877 FORMAT('WLTH(',I4,') =',F12.4,'   SABCO DATA ')                     
      WRITE(30,6870) (SABCO(J,K),J=1,NRHOX)                              
c      WRITE(30) K,WLTH(K)                                           
c      WRITE(30) (SABCO(J,K),J=1,NRHOX)                              
C                                                                       
      CALL JOSH(IFSCAT,IFSURF)                                          
      AINT(K)=HNU(1)                                                    
      RESIDF(K)=AINT(K)/CONFLX(LCONT)                                   
  953 CONTINUE                                                          

      SUMRES=0.5*RESIDF(1)
      DO 959 M=2,NPOINT-1
  959 SUMRES=SUMRES+RESIDF(M)
      SUMRES=SUMRES+0.5*RESIDF(NPOINT)
      BLOCKF=1.-SUMRES/FLOAT(NPOINT-1)
      WRITE(14,630) ISET,STWV(ISET),ENWV(ISET)
      WRITE(14,640) IVT,VTS(IVT),BLOCKF
      WRITE(14,650) NPOINT,NCON
      WRITE(14,660) (WLTH(K),K=1,NPOINT)
      WRITE(14,670) (AINT(K),K=1,NPOINT)
      WRITE(14,660) (WVCON(L),L=1,NCON)
      WRITE(14,670) (CONFLX(L),L=1,NCON)
      WRITE(14,680) (KCSV(K),K=1,NPOINT)
      WRITE(14,690) (RESIDF(K),K=1,NPOINT)
  630 FORMAT(I4,2F12.3)
  640 FORMAT(I4,F5.2,F7.3)
  650 FORMAT(2I10)
  660 FORMAT(10F12.4)
  670 FORMAT(1P10E12.4)
  680 FORMAT(20I6)
  690 FORMAT(15F8.4)
      WRITE(6,564)ISET,STWV(ISET),ENWV(ISET),NPOINT,IVT,VTS(IVT),BLOCKF 
  564 FORMAT(1H /1X,'<<<<< ISET=',I3,'  STWV=',F10.3,'  ENWV=',F10.3,2X,
     \'NPOINT=',I4,'  IVT=',I1,'  VTS=',F6.2,'  BLOCKF=',F6.3,' >>>>>')
      WRITE(6,*)
      WRITE(6,977) (RESIDF(K),K=1,NPOINT)                               
  977 FORMAT(1H ,'RESID:',20F6.3)                                       
  998 CONTINUE                                                          
      IF(ISET.NE.NMSET) GO TO 909                                       
      IFSYNT=0                                                          
      XNFH2=0.                                                          
      DO 989 J=1,NRHOX                                                  
      TXNXN(J)=(XNFPH(J,1)*2.+0.41336*XNFPHE(J,1)*1.+.85*XNFH2)*        
     #(T(J)/10000.)**.3                                                 
  989 CONTINUE                                                          
      GO TO 8                                                           
  988 CONTINUE                                                          
C     CALCULATES THE LINE CENTER MASS ABSORPTION COEFFICIENT FOR UNIT AB
C     CALL W(6HLINCEN,164.,1)                                           
      CALL LINCEN(ABUND1,VTS,NVT)                                       
C     CALL W(6HLINCEN,65.,1)                                            
      DO 65 IVT=1,NVT                                                   
      WRITE(6,565) IVT,VTS(IVT)                                         
  565 FORMAT(1H0,10X,'\\\\\\\\\\\\\\\\\\\\ IVT=',I2,'  VTURB=',F4.1,    
     1' \\\\\\\\\\\\\\\\\\\\'/)                                         
C     CALCULATES THE CURVE OF GROWTH AND THE FINAL ABUNDANCE            
C     FIRST GUESS FOR ABLG                                              
      ABLG=DLOG10(ABUND1)                                               
      IF(IFPROF.EQ.1)WRITE(6,64)                                        
   64 FORMAT(/23H PROFILES FOR NEXT LINE)                               
      CALL COG(IVT,ABLG,TAULG)                                          
C     CALCULATES THE AVERAGE ABUNDANCE                                  
      CALL AVERAG(MODEAV,IVT,ABLG,TAULG)                                
   65 CONTINUE                                                          
      WRITE(6,1400)                                                     
c      IF(MODECG.EQ.0) GO TO 183                                         
c      IFNPR=1                                                           
c      DO 182 IVT=1,NVT                                                  
c      NS=(IVT-1)*NABLOG+1                                               
c      DO 184 K=1,NABLOG                                                 
c      XD(NS+K-1)=PDATA(K,IVT,1)                                         
c      YD(NS+K-1)=PDATA(K,IVT,2)                                         
c      CHAR(NS+K-1)=GCHAR(IVT)                                           
c  184 CONTINUE                                                          
c  182 CONTINUE                                                          
c      NL=NVT*NABLOG                                                     
c      CALL PLOT1(EXPOT,ABLAB(MODECG),XMN,XMX,YMN,YMX,XD,YD,CHAR,        
c     # 1,NL,0,IFNPR)                                                    
  183 CONTINUE                                                          
      IF(NCOUNT.LT.NMLINE) GO TO 230                                    
      MODECG=0                                                          
      NABLOG=0                                                          
      GO TO 8                                                           
      END                                                               

      function invcode(code)
      implicit real*8 (a-h,o-z)
      dimension cm(68)
      data cm/
     a 100.00,101.00,101.01,106.00,106.01,107.00,108.00,108.01,
     b 109.00,112.00,113.00,114.00,114.01,116.00,117.00,600.00,
     c 606.00,607.00,607.01,608.00,613.00,614.00,616.00,707.00,
     d 707.01,708.00,708.01,712.00,713.00,714.00,716.00,800.00,
     e 808.00,808.01,812.00,813.00,814.00,816.00,820.00,822.00,
     f 826.00,900.00,1216.00,1300.00,1316.00,1416.00,1600.00,1616.00,
     g 1620.00,1626.00,1700.00,2600.00,10100.00,10106.00,10107.00,
     h 10108.00,10116.00,10606.00,10607.00,10608.00,10708.00,
     i 10800.00,10808.00,60600.00,60700.00,60808.00,70700.00,
     j 80800.00/
      if(code.lt.100.) then
        invcode=int(code+0.5)
        return
      end if
      do j=1,68
        if(code.eq.cm(j)) then
          invcode=-j
          return
        end if
      end do
      write(*,*) 'Do not understand code=',code
      call exit
      end function

      subroutine calzam(ijel,elm,tatmas,mola,iza)
      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON /ELEM/ABUND(99),ATMASS(99),ELEM(99)
      character*2 elem
      dimension xcode(8),iza(8)
      DATA XCODE/1.E14,1.E12,1.E10,1.E8,1.E6,1.E4,1.E2,1.E0/
      character*2 elm
      character*1 eq(2)
      c=dfloat(ijel)
      do 11 ii=1,8
      if(c.ge.xcode(ii)) go to 12
   11 continue
      call exit
   12 continue
      ia=0
      am=0.
      do 13 i=ii,8
      id=c/xcode(i)
      am=am+atmass(id)
      ia=ia+1
      iza(ia)=id
      c=c-dfloat(id)*xcode(i)
   13 continue
      mola=ia
      tatmas=am
      if(mola.eq.1) elm=elem(iza(1))
      if(mola.eq.2) then
        do j=1,2
          if(iza(j).gt.0) then
            eq(j)=elem(iza(j))
          else
            eq(j)='-'
          end if
        end do
        elm=eq(1)//eq(2)
      end if
      if(mola.eq.3) elm='tt'
      if(mola.eq.4) elm='xx'
      return
      end

      FUNCTION EXP10(X)                                                         
      IMPLICIT REAL*8 (A-H,O-Z)
      EXP10=EXP(X*2.30258509299405D0)                                           
      RETURN                                                                    
      END                                                                       
      SUBROUTINE COG(IVT,ABLG,TAULG)                                    
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (kw=99)
      COMMON /CONTIN/CONTIN                                             
      COMMON /CURVE/MINLOG,DABLOG,NABLOG,LOCONV
      LOGICAL LOCONV
      REAL*8 MINLOG                                                       
      COMMON /GAM/GLOGR,GLOGS,GLOGW                                     
      COMMON /LINDAT/WL,GFLOG,XJ,E,XJP,EP,CODE,LABEL,LABELP,NELION,     
     1               GAMMAR,GAMMAS,GAMMAW,REF,NBLO,NBUP,ISO1,X1,ISO2,   
     2               X2,OTHER1,OTHER2,ELO,GF                            
      COMMON /OBS/EW,WLOBS,STAR(9)                                      
      COMMON /RESIDU/RESIDU(15),TAUONE(15),DELWAV(15)                   
      COMMON /TAUBAR/TAUBAR                                             
      COMMON /VTS/VTS(3),NVT                                            
      COMMON /XLABEL/LAB,EXPOT                                          
      COMMON /ABDIF/ABCOR                                               
      COMMON /PLOTCG/MODECG,PDATA(9,3,2)                                
      DIMENSION SEWLOG(9)                                               
      DIMENSION ABLOG(9),EWLOG(9),TAULOG(9)                             
      DIMENSION RESID0(9)                                               
      DIMENSION WDOLM(9)                                                
      DIMENSION PEWLOG(9),PABLOG(9)                                     
      CHARACTER STAR*8,LAB*28                                           
C                                                                       
      IF(NABLOG.EQ.0)GO TO 63                                           
      ABLOG(1)=MINLOG                                                   
      EWLOG(1)=DLOG10(WID(EXP10(ABLOG(1)),IVT))                         
      RESID0(1)=RESIDU(1)                                               
      TAULOG(1)=DLOG10(TAUBAR)                                          
      IF(NABLOG.EQ.1) GO TO 23                                          
      DO 162 I=2,NABLOG                                                 
      ABLOG(I)=ABLOG(I-1)+DABLOG                                        
      EWLOG(I)=DLOG10(WID(EXP10(ABLOG(I)),IVT))                         
      RESID0(I)=RESIDU(1)                                               
  162 TAULOG(I)=DLOG10(TAUBAR)                                          
   23 CONTINUE                                                          
      NAB=NABLOG                                                        
      OBEWLG=0.                                                         
      GO TO 76                                                          
   63 OBEWLG=DLOG10(EW)                                                 
      ABLOG(1)=ABLG                                                     
      EWLOG(1)=DLOG10(WID(EXP10(ABLOG(1)),IVT))                         
      TAULOG(1)=DLOG10(TAUBAR)                                          
      IF(EWLOG(1).GT.OBEWLG)GO TO 71                                    
      ABLOG(2)=ABLOG(1)+1.                                              
      EWLOG(2)=DLOG10(WID(EXP10(ABLOG(2)),IVT))                         
      TAULOG(2)=DLOG10(TAUBAR)                                          
      ABLOG(3)=ABLOG(2)+1.                                              
      IF(EWLOG(2).GT.OBEWLG)ABLOG(3)=ABLOG(3)-1.5                       
      GO TO 72                                                          
   71 ABLOG(2)=ABLOG(1)-1.                                              
      EWLOG(2)=DLOG10(WID(EXP10(ABLOG(2)),IVT))                         
      TAULOG(2)=DLOG10(TAUBAR)                                          
      ABLOG(3)=ABLOG(2)+.5                                              
      IF(EWLOG(2).GT.OBEWLG)ABLOG(3)=ABLOG(3)-1.5                       
   72 EWLOG(3)=DLOG10(WID(EXP10(ABLOG(3)),IVT))                         
      TAULOG(3)=DLOG10(TAUBAR)                                          
      DO 73 IAB=3,8                                                     
      NAB=IAB                                                           
      IF(ABS(OBEWLG-EWLOG(IAB)).LT..005)GO TO 76                        
      NAB1=NAB-1                                                        
      DO 79 KK=1,NAB1                                                   
      DO 79 II=1,NAB1                                                   
      IF(ABLOG(II+1).GE.ABLOG(II))GO TO 79                              
      DUMMY=ABLOG(II+1)                                                 
      ABLOG(II+1)=ABLOG(II)                                             
      ABLOG(II)=DUMMY                                                   
      DUMMY=EWLOG(II+1)                                                 
      EWLOG(II+1)=EWLOG(II)                                             
      EWLOG(II)=DUMMY                                                   
      DUMMY=TAULOG(II+1)                                                
      TAULOG(II+1)=TAULOG(II)                                           
      TAULOG(II)=DUMMY                                                  
   79 CONTINUE                                                          
      DO 80 I=2,NAB                                                     
      K=I                                                               
      IF(OBEWLG.LT.EWLOG(I))GO TO 81                                    
   80 CONTINUE                                                          
   81 ABLOG(NAB+1)=ABLOG(K-1)+(ABLOG(K)-ABLOG(K-1))/(EWLOG(K)-EWLOG(K-1)
     1)*(OBEWLG-EWLOG(K-1))                                             
      EWLOG(IAB+1)=DLOG10(WID(EXP10(ABLOG(IAB+1)),IVT))                 
   73 TAULOG(IAB+1)=DLOG10(TAUBAR)                                      
      WRITE(6,75)                                                       
   75 FORMAT(//14H NOT CONVERGED)                                       
      NAB=NAB+1                                                         
   76 CONTINUE                                                          
C     CONVERT LOG(N/NT) TO LOG(N/NH)+12                                 
C     CONVERT PM TO MA AND NM TO A                                      
      EEW=EW*10.0                                                       
      WWL=WL*10.0                                                       
      DO 85 IAB=1,NAB                                                   
      PEWLOG(IAB)=EWLOG(IAB)+1.0                                        
   85 PABLOG(IAB)=ABLOG(IAB)+12.0-ABCOR                                 
      WRITE(6,3333) CODE,WWL,EXPOT,GFLOG,GLOGR,GLOGS,GLOGW,LAB          
 3333 FORMAT(1H0,5X,'EL:',F5.2,5X,'WL:',F9.3,5X,'CHI:',F5.2,5X,         
     1 'LOGGF:',F6.3,5X,'GRLOG:',F6.2,5X,'GSLOG:',F6.2,5X,              
     2 'GWLOG:',F6.2,5X/10X,'LABEL:'A28/)                               
      IF(NABLOG.EQ.0) WRITE(6,820) EEW,PABLOG(NAB)                      
  820 FORMAT(1H0,10X,'***** FOR OBSERVED E.W. ',F6.1,' M.A.  ',         
     1 'LOG(N/NH)+12 IS ',F6.2,' *****'/)                               
      WRITE(6,77)(PABLOG(IAB),IAB=1,NAB)                                
   77 FORMAT(9X,5HVTURB,2X,5HABUND,9F9.2)                               
      WRITE(6,78)VTS(IVT),(PEWLOG(IAB),IAB=1,NAB)                       
   78 FORMAT(9X,F5.2,2X,'EWLOG',9F9.3)                                   
      ALGWL=DLOG10(WWL)                                                 
      DO 113 IAB=1,NAB                                                  
      WDOLM(IAB)=PEWLOG(IAB)-ALGWL-3.0                                  
      SEWLOG(IAB)=PEWLOG(IAB)                                           
  113 PEWLOG(IAB)=EXP10(PEWLOG(IAB))                                    
      WRITE(6,114)(PEWLOG(IAB),IAB=1,NAB)                               
  114 FORMAT(18X,'EW ',9F9.2)                                           
      WRITE(6,177)(TAULOG(IAB),IAB=1,NAB)                               
  177 FORMAT(16X,5HLGTAU,9F9.2)                                         
      IF(NABLOG.GT.0)WRITE(6,178)(RESID0(IAB),IAB=1,NAB)                
  178 FORMAT(16X,5HRESID,9F9.2)                                         
      IF(NABLOG.GT.0) WRITE(6,179) (WDOLM(IAB),IAB=1,NAB)               
  179 FORMAT(14X,'LG(W/L)',9F9.2)                                       
      ABLG=PABLOG(NAB)                                                  
      TAULG=TAULOG(NAB)                                                 
      IF(MODECG.EQ.0) RETURN                                            
      IF(MODECG.EQ.2) GO TO 182                                         
      DO 181 I=1,NABLOG                                                 
      PDATA(I,IVT,1)=PABLOG(I)                                          
  181 PDATA(I,IVT,2)=SEWLOG(I)                                          
      RETURN                                                            
  182 DO 183 I=1,NABLOG                                                 
  183 PABLOG(I)=PABLOG(I)+GFLOG                                         
      DO 184 I=1,NABLOG                                                 
      PDATA(I,IVT,1)=PABLOG(I)                                          
  184 PDATA(I,IVT,2)=WDOLM(I)                                           
      RETURN                                                            
      END                                                               
      FUNCTION WID(ABUND,IVT)                                           
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (kw=99)
      COMMON /IF/IFCORR,IFPRES,IFSURF,IFSCAT,TAUSCAT,IFMOL                      
      COMMON /MUS/ANGLE(20),SURFI(20),NMU                               
      COMMON /OPTOT/ACONT(kw),SCONT(kw),ALINE(kw),SLINE(kw),SIGMAC(kw), 
     1              SIGMAL(kw)                                          
      COMMON /RHOX/RHOX(kw),NRHOX                                       
      COMMON /ABROSS/ABROSS(kw),TAUROS(kw)                              
      COMMON /TAUSHJ/TAUNU(kw),SNU(kw),HNU(kw),JNU(kw),JMINS(kw)        
      REAL*8 JNU,JMINS                                                    
      COMMON /TURBPR/VTURB(kw),PTURB(kw),TRBFDG,TRBCON,TRBPOW,TRBSND,   
     1               IFTURB                                             
C                                                                       
      COMMON /CONTIN/CONTIN                                             
      COMMON /IFPROF/IFPROF                                             
      COMMON /LINDAT/WL,GFLOG,XJ,E,XJP,EP,CODE,LABEL,LABELP,NELION,     
     1               GAMMAR,GAMMAS,GAMMAW,REF,NBLO,NBUP,ISO1,X1,ISO2,   
     2               X2,OTHER1,OTHER2,ELO,GF                            
      COMMON /LINEY/LINE0(3,kw),DOPWAV(3,kw),ADAMP(3,kw)                
      REAL*8 LINE0                                                        
      COMMON /RESIDU/RESIDU(15),TAUONE(15),DELWAV(15)                   
      COMMON /TAUBAR/TAUBAR                                             
      COMMON /VTS/VTS(3),NVT                                            
      DIMENSION DEL(15),WEIGHT(15),ALINE0(kw)                           
      DIMENSION PDELS(15)                                               
C     IN PM                                                             
      DATA DEL/0.,.5,1.,1.5,2.,2.5,3.,3.5,4.,4.5,5.,6.25,7.5,8.75,10./  
      DATA WEIGHT/.166667,.666667,.333333,.666667,.333333,.666667,      
     1.333333,.666667,.333333,.666667,.583333,1.666667,.833333,1.666667,
     2 10.416667/                                                       
c      EXP10(X)=EXP(X*2.30258509299405E0)                                
      DO 1 J=1,NRHOX                                                    
    1 ALINE0(J)=LINE0(IVT,J)*ABUND                                      
      ASSIGN 73 TO ISWTCH                                               
C     MICROTURBULENCE IN PM                                             
      DOPVT=(VTS(IVT)*(1.E5/2.997925E10))*(WL*1000.)                    
      SCALE=SQRT(2.5**2+DOPVT**2)/2.5                                   
      DO 70 MDUMMY=1,10                                                 
      CCC=5.*SCALE                                                      
      DO 71 J=1,NRHOX                                                   
   71 ALINE(J)=ALINE0(J)*VOIGT(CCC/DOPWAV(IVT,J),ADAMP(IVT,J))          
      CALL JOSH(IFSCAT,IFSURF)                                          
      IF(IFSURF.LT.2)RESID=HNU(1)/CONTIN                                
      IF(IFSURF.EQ.2)RESID=SURFI(1)/CONTIN                              
      IF(MDUMMY.EQ.1)GO TO 75                                           
      IF(RESID.GT..93)GO TO 74                                          
      GO TO ISWTCH,(73,76)                                              
   75 IF(RESID.GE..91)GO TO 80                                          
   73 SCALE1=SCALE                                                      
      SCALE=SCALE*2.                                                    
      ASSIGN 73 TO ISWTCH                                               
      GO TO 70                                                          
   74 SCALE2=SCALE                                                      
      SCALE=(SCALE1+SCALE2)/2.                                          
      ASSIGN 76 TO ISWTCH                                               
      GO TO 70                                                          
   76 SCALE1=SCALE                                                      
      SCALE=(SCALE1+SCALE2)/2.                                          
   70 CONTINUE                                                          
   80 WIDTH=0.                                                          
      TAUSUM=0.                                                         
      DO 81 I=1,15                                                      
      CCC=DEL(I)*SCALE                                                  
      DELWAV(I)=CCC                                                     
      DO 82 J=1,NRHOX                                                   
   82 ALINE(J)=ALINE0(J)*VOIGT(CCC/DOPWAV(IVT,J),ADAMP(IVT,J))          
      CALL JOSH(IFSCAT,IFSURF)                                          
      IF(IFSURF.LT.2)RESIDU(I)=HNU(1)/CONTIN                            
      IF(IFSURF.EQ.2)RESIDU(I)=SURFI(1)/CONTIN                          
      MAX=MAP1(TAUNU,TAUROS,NRHOX,1.,TAUONE(I),1)                       
      TAUONE(I)=DLOG10(TAUONE(I))                                       
      TAUSUM=TAUSUM+WEIGHT(I)*TAUONE(I)*(1.-RESIDU(I))                  
   81 WIDTH=WIDTH+WEIGHT(I)*(1.-RESIDU(I))                              
      TAUBAR=EXP10(TAUSUM/WIDTH)                                        
      WID=WIDTH*SCALE*2.                                                
      IF(IFPROF.EQ.0)RETURN                                             
C     CONVERT PM TO MA                                                  
      DO 85 I=1,15                                                      
   85 PDELS(I)=DELWAV(I)*10.0                                           
      WRITE(6,87)PDELS,TAUONE,RESIDU                                    
   87 FORMAT(1H0,17X,3HDEL,15F7.2/16X,5HLGTAU,15F7.3/16X,5HRESID,15F7.3)
      RETURN                                                            
      END                                                               
      SUBROUTINE LINCEN(ABUND1,VTS,NVT)                                 
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (kw=99)
      COMMON /ELEM/ABUND(99),ATMASS(99),ELEM(99)                        
      character*2 elem
      COMMON /FREQ/FREQ,FREQLG,EHVKT(kw),STIM(kw),BNU(kw)
      COMMON /IONS/XNFPH(kw,2),XNFPHE(kw,3),XNFH(kw,2),XNFHE(kw,3)                             
      COMMON /IF/IFCORR,IFPRES,IFSURF,IFSCAT,TAUSCAT,IFMOL                            
      COMMON /RHOX/RHOX(kw),NRHOX                                       
      COMMON /STATE/P(kw),XNE(kw),XNATOM(kw),RHO(kw)                    
      COMMON /TEMP/T(kw),TKEV(kw),TK(kw),HKT(kw),TLOG(kw),ITEMP         
      COMMON /TURBPR/VTURB(kw),PTURB(kw),TRBFDG,TRBCON,TRBPOW,TRBSND,   
     1               IFTURB                                             
      COMMON /XABUND/XABUND(99),WTMOLE                                  
C                                                                       
      COMMON /LINDAT/WL,GFLOG,XJ,E,XJP,EP,CODE,LABEL,LABELP,NELION,     
     1               GAMMAR,GAMMAS,GAMMAW,REF,NBLO,NBUP,ISO1,X1,ISO2,   
     2               X2,OTHER1,OTHER2,ELO,GF                            
      COMMON /LINEY/LINE0(3,kw),DOPWAV(3,kw),ADAMP(3,kw)                
      COMMON /WDAMP/TXNXN(kw)                                           
      REAL*8 LINE0                                                        
      COMMON /RESIDU/RESIDU(15),TAUONE(15),DELWAV(15)                   
      DIMENSION VTS(3)                                                  
      DIMENSION TDOP(kw),ADAMP1(kw),LINE(kw),XNFPRC(kw)                 
      REAL*8 LINE                                                         
      DIMENSION XCODE(8)                                                
      DATA XCODE/1.E14,1.E12,1.E10,1.E8,1.E6,1.E4,1.E2,1.E0/            
      DATA ITEMP1/0/                                                    
      IF(ITEMP.EQ.ITEMP1)GO TO 2                                        
      XNFH2=0.                                                          
      DO 1 J=1,NRHOX                                                    
C     NUMBER DENSITIES FOR VAN DER WAALS BROADENING                     
    1 TXNXN(J)=(XNFPH(J,1)*2.+.42*XNFPHE(J,1)*1.+.85*XNFH2)*            
     1(T(J)/10000.)**.3                                                 
      ITEMP1=ITEMP                                                      
      SAVE=0.                                                           
    2 IF(CODE.EQ.SAVE)GO TO 18                                          
      AMASS=0.                                                          
      ABUND1=1.                                                         
      C=CODE                                                            
      DO 11 II=1,8                                                      
      IF(C.GE.XCODE(II))GO TO 12                                        
   11 CONTINUE                                                          
      CALL EXIT                                                         
   12 DO 13 I=II,8                                                      
      ID=C/XCODE(I)                                                     
      if(id.eq.0) go to 13
      IF(ID.LT.100)ABUND1=ABUND1*XABUND(ID)
      AMASS=AMASS+ATMASS(ID)                                            
      C=C-FLOAT(ID)*XCODE(I)
   13 CONTINUE
      SAVE=CODE                                                         
      CALL POPS(CODE,1,XNFPRC)                                          
      TWOMAS=2./AMASS/1.660E-24                                         
      DO 17 J=1,NRHOX                                                   
      XNFPRC(J)=XNFPRC(J)*(.026538/1.77245)/ABUND1/RHO(J)               
   17 TDOP(J)=TK(J)*TWOMAS                                              
   18 ELOC=ELO*2.997925E10                                              
      DO 5 J=1,NRHOX                                                    
      ADAMP1(J)=(GAMMAR+GAMMAS*XNE(J)+GAMMAW*TXNXN(J))/12.5664          
    5 LINE(J)=GF*EXP(-ELOC*HKT(J))*STIM(J)*XNFPRC(J)                    
      FREQC=FREQ/2.997925E10                                            
      WAVENU=1.E10/FREQC/FREQ                                           
      IF(VTS(1).EQ.0..AND.NVT.EQ.1)GO TO 7                              
      DO 6 IVT=1,NVT                                                    
      VT2=VTS(IVT)**2*1.E10                                             
      DO 6 J=1,NRHOX                                                    
      DOPNU=SQRT(TDOP(J)+VT2)*FREQC                                     
C     DOPWAV IS DOPLER WIDTH IN PM                                      
      DOPWAV(IVT,J)=DOPNU*WAVENU                                        
C     LINE ABSORPTION COEFFICIENT PER UNIT ABUNDANCE                    
      LINE0(IVT,J)=LINE(J)/DOPNU                                        
    6 ADAMP(IVT,J)=ADAMP1(J)/DOPNU                                      
      RETURN                                                            
    7 DO 8 J=1,NRHOX                                                    
      DOPNU=SQRT(TDOP(J)+VTURB(J)**2)*FREQC                             
      DOPWAV(1,J)=DOPNU*WAVENU                                          
      LINE0(1,J)=LINE(J)/DOPNU                                          
    8 ADAMP(1,J)=ADAMP1(J)/DOPNU                                        
      RETURN                                                            
      END                                                               
      FUNCTION H(V,A)                                                   
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (kw=99)
C     FUNCTION VOIGT(V,A)                                               
C     INTERPOLATES FROM TABLES OF HARRIS AP.J. 1948                     
      DIMENSION H0(41),H1(81),H2(41)                                    
      DATA H0 /                                                         
     1 1.0000000, 0.9900500, 0.9607890, 0.9139310, 0.8521440, 0.7788010,
     2 0.6976760, 0.6126260, 0.5272920, 0.4448580, 0.3678790, 0.2981970,
     3 0.2369280, 0.1845200, 0.1408580, 0.1053990, 0.0773050, 0.0555760,
     4 0.0391640, 0.0270520, 0.0183156, 0.0121552, 0.0079071, 0.0050418,
     5 0.0031511, 0.0019305, 0.0011592, 0.0006823, 0.0003937, 0.0002226,
     6 0.0001234, 0.0000671, 0.0000357, 0.0000186, 0.0000095, 0.0000048,
     7 0.0000024, 0.0000011, 0.0000005, 0.0000002, 0.0000001/           
      DATA H1/                                                          
     1-1.1283800,-1.1059600,-1.0404800,-0.9370300,-0.8034600,-0.6494500,
     2-0.4855200,-0.3219200,-0.1677200,-0.0301200, 0.0859400, 0.1778900,
     3 0.2453700, 0.2898100, 0.3139400, 0.3213000, 0.3157300, 0.3009400,
     4 0.2802700, 0.2564800, 0.2317260, 0.207528 , 0.1848820, 0.1643410,
     5 0.1461280, 0.1302360, 0.1165150, 0.1047390, 0.0946530, 0.0860050,
     6 0.0785650, 0.0721290, 0.0665260, 0.0616150, 0.0572810, 0.0534300,
     7 0.0499880, 0.0468940, 0.0440980, 0.0415610, 0.0392500, 0.0351950,
     8 0.0317620, 0.0288240, 0.0262880, 0.0240810, 0.0221460, 0.0204410,
     9 0.0189290, 0.0175820, 0.0163750, 0.0152910, 0.0143120, 0.0134260,
     A 0.0126200, 0.0118860, 0.0112145, 0.0105990, 0.0100332, 0.0095119,
     B 0.0090306, 0.0085852, 0.0081722, 0.0077885, 0.0074314, 0.0070985,
     C 0.0067875, 0.0064967, 0.0062243, 0.0059688, 0.0057287, 0.0055030,
     D 0.0052903, 0.0050898, 0.0049006, 0.0047217, 0.0045526, 0.0043924,
     E 0.0042405, 0.0040964, 0.0039595/                                 
      DATA H2 /                                                         
     1 1.0000000, 0.9702000, 0.8839000, 0.7494000, 0.5795000, 0.3894000,
     2 0.1953000, 0.0123000,-0.1476000,-0.2758000,-0.3679000,-0.4234000,
     3-0.4454000,-0.4392000,-0.4113000,-0.3689000,-0.3185000,-0.2657000,
     4-0.2146000,-0.1683000,-0.1282100,-0.0950500,-0.0686300,-0.0483000,
     5-0.0331500,-0.0222000,-0.0145100,-0.0092700,-0.0057800,-0.0035200,
     6-0.0021000,-0.0012200,-0.0007000,-0.0003900,-0.0002100,-0.0001100,
     7-0.0000600,-0.0000300,-0.0000100,-0.0000100, 0.0000000/           
      V0=V*10.                                                          
      N=V0                                                              
      IF(N.LT.40)GO TO 1                                                
      IF(N.LT.120)GO TO 2                                               
      VOIGT=(.56419+.846/(V**2))/(V**2)*A                               
      H=VOIGT                                                           
      RETURN                                                            
    1 V1=N                                                              
      N=N+1                                                             
      V2=V0-V1                                                          
      N1=N+1                                                            
      VOIGT=V2*(H0(N1)-H0(N)+A*(H1(N1)-H1(N)+A*(H2(N1)-H2(N))))+        
     1H0(N)+A*(H1(N)+A*H2(N))                                           
      H=VOIGT                                                           
      RETURN                                                            
    2 N=N/2+20                                                          
      V1=(N-20)*2                                                       
      N=N+1                                                             
      V2=(V0-V1)/2.                                                     
      N1=N+1                                                            
      VOIGT=A*((H1(N1)-H1(N))*V2+H1(N))                                 
      H=VOIGT                                                           
      RETURN                                                            
      END                                                               
      SUBROUTINE AVERAG(MODEAV,IVT,ABLG,TAULG)                          
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (kw=99)
      COMMON /LINDAT/WL,GFLOG,XJ,E,XJP,EP,CODE,LABEL,LABELP,NELION,     
     1               GAMMAR,GAMMAS,GAMMAW,REF,NBLO,NBUP,ISO1,X1,ISO2,   
     2               X2,OTHER1,OTHER2,ELO,GF                            
      COMMON /OBS/EW,WLOBS,STAR(9)                                      
      COMMON /VTS/VTS(3),NVT                                            
      COMMON ABUNDI(250,3),TAUBAR(250,3)                                
      DIMENSION KABUND(250,3)                                           
      EQUIVALENCE (KABUND(1,1),ABUNDI(1,1))                                 
      DIMENSION EPOT(250),WDTH(250)                                     
      DIMENSION STAR1(9)                                                
      CHARACTER*8 STAR,STAR1                                            
      DATA NLINES/0/                                                    
c      EXP10(X)=EXP(X*2.30258509299405E0)                                
      IF(MODEAV.EQ.1)RETURN                                             
      IF(MODEAV.EQ.2)GO TO 64                                           
      IF(MODEAV.EQ.3)GO TO 30                                           
      CALL EXIT                                                         
   64 IF(IVT.EQ.1)NLINES=NLINES+1                                       
      ABUNDI(NLINES,IVT)=ABLG                                           
      TAUBAR(NLINES,IVT)=TAULG                                          
      EPOT(NLINES)=ELO*1.23981E-4                                       
      WDTH(NLINES)=EW                                                   
      DO 25 I=1,9                                                       
   25 STAR1(I)=STAR(I)                                                  
      SAVE=CODE                                                         
      RETURN                                                            
   30 IF(NLINES.NE.1) GO TO 31                                          
      DEV=0.                                                            
      CODE=SAVE                                                         
      WRITE(6,33)NLINES,CODE,ABLG,DEV                                   
      GO TO 35                                                          
   31 DO 34 IVT=1,NVT                                                   
      SUMM=0.                                                           
      CSUMXY=0.                                                         
      CSUMX=0.                                                          
      CSUMXX=0.                                                         
      WSUMXY=0.                                                         
      WSUMX=0.                                                          
      WSUMXX=0.                                                         
      DO 6005 I=1,NLINES                                                
      SUMM=SUMM+ABUNDI(I,IVT)                                           
      CSUMXY=CSUMXY+ABUNDI(I,IVT)*EPOT(I)                               
      CSUMX=CSUMX+EPOT(I)                                               
      CSUMXX=CSUMXX+EPOT(I)**2                                          
      WSUMXY=WSUMXY+ABUNDI(I,IVT)*WDTH(I)                               
      WSUMX=WSUMX+WDTH(I)                                               
      WSUMXX=WSUMXX+WDTH(I)**2                                          
 6005 CONTINUE                                                          
      X=NLINES                                                          
      WB=( X*WSUMXY-WSUMX*SUMM)/(X*WSUMXX-WSUMX**2)                     
      CB=(X*CSUMXY-CSUMX*SUMM)/(X*CSUMXX-CSUMX**2)                      
      AVABN=SUMM/X                                                      
      SUM=0.                                                            
      DO 32 I=1,NLINES                                                  
   32 SUM=SUM+(AVABN-ABUNDI(I,IVT))**2                                  
      DEV=SQRT(SUM/X)                                                   
      DO 26 I=1,6                                                       
   26 STAR(I)=STAR1(I)                                                  
      CODE=SAVE                                                         
      WRITE(6,1111)VTS(IVT),STAR                                        
 1111 FORMAT(6H1VTURB,F6.3,3X,9A8)                                      
      WRITE(6,33)NLINES,CODE,AVABN,DEV                                  
   33 FORMAT(///26H ***** THE ABUNDANCE FROM I3,2X,F11.2,9H LINES IS,
     1 F7.2,3H+/-,F5.2)                                                      
      WRITE(6,6006)                                                     
 6006 FORMAT(//5X,17HLOG ABUND VS. CHI)                                 
      DO 50 I=1,NLINES                                                  
   50 KABUND(I,IVT)=(AVABN-ABUNDI(I,IVT))*50.+51.5                      
      WRITE(6,6009)CB,NLINES                                            
 6009 FORMAT(5X,7H SLOPE=,E10.3,14H PER EV USING ,I4,6H LINES)          
      CALL PLOTIT(AVABN,KABUND(1,IVT),EPOT,NLINES)                      
      WRITE(6,1111)VTS(IVT),STAR                                        
      WRITE(6,6007)                                                     
 6007 FORMAT(//5X,23HLOG ABUND VS. EQ. WIDTH)                           
      WRITE(6,6008)WB,NLINES                                            
 6008 FORMAT(5X,7H SLOPE=,E10.3,14H PER MA USING ,I4,6H LINES)          
      CALL PLOTIT(AVABN,KABUND(1,IVT),WDTH,NLINES)                      
      WRITE(6,1111)VTS(IVT),STAR                                        
      WRITE(6,1112)                                                     
 1112 FORMAT(//5X,24HLOG ABUND VS. LOG HEIGHT)                          
      CALL PLOTIT(AVABN,KABUND(1,IVT),TAUBAR(1,IVT),NLINES)             
   34 CONTINUE                                                          
   35 NLINES=0                                                          
      MODEAV=1                                                          
      WRITE(6,3333)                                                     
 3333 FORMAT(1H1)                                                       
      RETURN                                                            
      END                                                               
      SUBROUTINE PLOTIT(AVX,KX,Y,N)                                     
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (kw=99)
      DIMENSION KX(1),Y(1)                                              
      DIMENSION IPLOT(101),AXIS(11)                                     
      DATA IBLANK,IDOT,IMINUS,II,IX,I2/1H ,1H.,1H-,1HI,1HX,1H2/         
      DO 6  I=1,11                                                      
    6 AXIS(I)=AVX+.2*FLOAT(6-I)                                         
      YMAX=-1.E30                                                       
      YMIN=1.E30                                                        
      DO 1 I=1,N                                                        
      YMIN=min(YMIN,Y(I))                                             
    1 YMAX=max(YMAX,Y(I))                                             
      IF(YMIN.GT.0.)YMIN=0.                                             
      SLOPE=50./(YMAX-YMIN)                                             
      DO 3  J=1,51                                                      
      ISYMB=IBLANK                                                      
      IF(MOD(J,10).EQ.1)ISYMB=IMINUS                                    
      DO 4  K=1,101                                                     
    4 IPLOT(K)=ISYMB                                                    
      DO 44 K=2,100,2                                                   
   44 IPLOT(K)=IBLANK                                                   
      DO 40 K=1,101,10                                                  
   40 IPLOT(K)=IDOT                                                     
      IPLOT(51)=II                                                      
      DO 5  I=1,N                                                       
      IY=51.5-SLOPE*(Y(I)-YMIN)                                         
      IF(IY.NE.J)GO TO 5                                                
      L=KX(I)                                                           
      IF(L.LT.1) GO TO 5                                                
      IF(L.GT.101)  GO TO 5                                             
      IF(IPLOT(L).EQ.IX) GO TO 7                                        
      IPLOT(L)=IX                                                       
      GO TO 5                                                           
    7 IPLOT(L)=I2                                                       
    5 CONTINUE                                                          
      YY=FLOAT(51-J)/SLOPE+YMIN                                         
      WRITE(6,10)YY,IPLOT                                               
   10 FORMAT(F12.2,1X,101A1)                                             
    3 CONTINUE                                                          
      WRITE(6,20) AXIS                                                  
   20 FORMAT(5X,11F10.2)                                                
      RETURN                                                            
      END                                                               
      FUNCTION VOIGT(V,A)                                               
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (kw=99)
C     FUNCTION H(A,V)                                                   
C     FROM ERIC PEYTREMANN                                              
      LOGICAL Q                                                         
      VV=V*V                                                            
      Q=A.LT.0.2                                                        
      IF(Q.AND.V.GT.5.)GOTO 1                                           
      IF(.NOT.Q.AND.(A.GT.1.4.OR.A+V.GT.3.2))GOTO 2                     
      HO=EXP(-VV)                                                       
      H2=(1.-2.*VV)*HO                                                  
      IF(V.GT.2.4)GOTO 3                                                
      IF(V.GT.1.3)GOTO 4                                                
      H1=(.42139*VV-2.34358*V+3.28868)*VV-.15517*V-1.1247               
    5 IF(Q)GOTO 6                                                       
      HH1=H1+HO*1.12838                                                 
      HH2=H2+HH1*1.12838-HO                                             
      HH3=(1.-H2)*.37613-HH1*.66667*VV+HH2*1.12838                      
      HH4=(3.*HH3-HH1)*.37613+HO*.66667*VV*VV                           
C     H=((((HH4*A+HH3)*A+HH2)*A+HH1)*A+HO)*(((-.122727278*A+.532770573)*
      VOIGT=                                                            
     A  ((((HH4*A+HH3)*A+HH2)*A+HH1)*A+HO)*(((-.122727278*A+.532770573)*
     1A-.96284325)*A+.979895032)                                        
      RETURN                                                            
    1 VOIGT=((2.12/VV+.8463)/VV+.5642)*A/VV                             
C   1 H=((2.12/VV+.8463)/VV+.5642)*A/VV                                 
      RETURN                                                            
    2 AA=A*A                                                            
      U=(AA+VV)*1.4142                                                  
      UU=U*U                                                            
C     H=((((AA-10.*VV)*AA*3.+15.*VV*VV)/UU+3.*VV-AA)/UU+1.)*A*.79788/U  
      VOIGT=                                                            
     A  ((((AA-10.*VV)*AA*3.+15.*VV*VV)/UU+3.*VV-AA)/UU+1.)*A*.79788/U  
      RETURN                                                            
    3 H1=((-.0032783*VV+.0429913*V-.188326)*VV+.278712*V+.55415)/(VV-1.5
     1)                                                                 
      GOTO 5                                                            
    4 H1=(-.220416*VV+1.989196*V-6.61487)*VV+9.39456*V-4.4848           
      GOTO 5                                                            
    6 VOIGT=(H2*A+H1)*A+HO                                              
C   6 H=(H2*A+H1)*A+HO                                                  
      RETURN                                                            
      END
      subroutine fileopen
c
c     version for sequential files
c
c     This is a subprogram which opens as many files (with prescribed
c     unit number and filename) as one may wish.
c     !!! Note that SEQUENTIAL FILES are implicitly !!!
c     !!! assumed in this program !!!
c     The necessary file-data should be prepared in advance
c     in a text file, e.g., the one named 'fileopen.dat'.
c     In the first column, write any one of the following flags
c     specifying the data-form of the file (blank or any other
c     character is not allowed):
c     F or f  -----> FORMATTED
c     U or u  -----> UNFORMATTED
c     B or b  -----> BINARY
c     The parameter nunit0, assumed to be within 2<=|nunit0|<=99,
c     must be written in columns 2-4 with the format (i3).
c     If nunit0 is positive, the access mode is 'SEQUENTIAL'.
c     Meanwhile, negative nunit0 implies the access mode of 'APPEND'.
c     The name of the file (fname) may be written whereever you wish
c     if it isplaced within columns 5-68.
c     (i.e.,head-blanks and tail-blanks are neglected).
c     Note, however, that blanks should not be inserted within fname.
c     When all of the fname-field is kept blank, a scratch file
c     will be assigned.
c
      implicit real*8 (a-h,o-z)
      character*64 fname
      character*12 fform
      character*1 flag
      write(*,*) '<<<<<<<<<<<< files used for this job >>>>>>>>>>>>'
    1 read(1,'(a1,i3,a64)',err=998,end=999) flag,nunit0,fname
      nunit=iabs(nunit0)
      if(flag.eq.'F'.or.flag.eq.'f') then
      fform='formatted'
      go to 8
      end if
      if(flag.eq.'U'.or.flag.eq.'u') then
      fform='unformatted'
      go to 8
      end if
      if(flag.eq.'B'.or.flag.eq.'b') then
      fform='binary'
      go to 8
      end if
      write(*,*) 'I do not understand the flag --> ',flag
      stop
    8 continue
      do 10 i=1,64
      if(fname(i:i).ne.' ') go to 20
   10 continue
      open(unit=nunit,form=fform,status='scratch')
      write(*,'('' flag='',a1,2x,''nunit0='',i3,2x,''(scratch file)'')'
     %) flag,nunit0
      go to 30
   20 continue
      is=i
      do 110 i=64,1,-1
      if(fname(i:i).ne.' ') go to 120
  110 continue
  120 ie=i
      if(nunit0.gt.0) then
      open(unit=nunit,form=fform,file=fname,access='SEQUENTIAL')
      else
      open(unit=nunit,form=fform,file=fname,access='APPEND')
      end if
      write(*,'('' flag='',a1,2x,''nunit0='',i3,2x,''fname='',(a))')
     % flag,nunit0,fname(is:ie)
   30 continue
      go to 1
  998 write(*,*) '!!! read error in subroutine fileopen !!!'
      stop
  999 write(*,*)
      return
      end

      subroutine lablen(label,nlen,ns,ne)
      implicit real*8 (a-h,o-z)
      character*(*) label
      do n=nlen,1,-1
      if((label(n:n).ne.' ').and.(label(n:n).ne.'')) go to 10
      end do
      n=1
  10  ne=n
      do n=1,nlen
      if((label(n:n).ne.' ').and.(label(n:n).ne.'')) go to 20
      end do
      n=nlen
  20  ns=n
      return
      end

