      program mprep
      USE DFLIB
      use dfport
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (kw=99)
      character*80 line,linedir,moddir,datdir
      character argv*60, fnmesh*60
      integer*2 na
      logical loexist
CT****
      EXTERNAL BLOCKE,BLOCKR,BLOCKT
      EXTERNAL FSET
CT****
C     STANDARD OPTICAL DEPTH IS TAU(5000)
      PARAMETER (MD=99)
      DIMENSION TAUL1(MD),TAUL2(MD),VT1(MD),VT2(MD),RHOL1(MD),RHOL2(MD)
      DIMENSION RHOXL1(MD),T1(MD),PL1(MD),XNEL1(MD),ABSL1(MD),ACCRL1(MD)
      DIMENSION RHOXL2(MD),T2(MD),PL2(MD),XNEL2(MD),ABSL2(MD),ACCRL2(MD)
      DIMENSION DUMMY(MD),HP1(MD),HP2(MD),TAUD(MD),GMACC2(MD),PRDM(MD)
      DIMENSION AKAP5(MD),TAU5(MD)
      DIMENSION ABS5L1(MD),ABS5L2(MD),TAURL1(MD),TAURL2(MD)
      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)
      CHARACTER*3 A,B
c      character*80 fname
c      write(*,*)
c     % 'Please input the name of the file including fileopen data:'
c      read(*,'(a80)') fname
c      open (unit=1,file=fname)
c      call fileopen
c      open (unit=6,file='mprepout.out')

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

      open (unit=6,file='mprepout.out',disp='delete')
      na=1
      call getarg(na,argv)
      read(argv,*) ifmole

c      open (88,file='passdir.tmp')
c      read(88,'(a)') linedir
c      read(88,'(a)') moddir

c     get SPT_DIR
c      call getenv("SPT_DIR",linedir)
      call retenv("SPT_DIR",linedir)
c      na=2
c      call getarg(na,argv)
c      read(argv,'(a)') linedir
      call lablen(linedir,80,nsdir,nedir)

c     get MOD_DIR
c      call getenv("MOD_DIR",moddir)
      call retenv("MOD_DIR",moddir)
c      na=3
c      call getarg(na,argv)
c      read(argv,'(a)') moddir
      call lablen(moddir,80,msdir,medir)

c     get DAT_DIR
      call retenv("DAT_DIR",datdir)
      call lablen(datdir,80,ndsdir,ndedir)

      na=2
      call getarg(na,argv)
      read(argv,'(a)',err=69,end=69) fnmesh
      inquire(file=fnmesh,exist=loexist)
      if(.not.loexist) go to 69
c      call lablen(fnmesh,nj1,nj2)
c      write(*,'(''fnmesh:'',(a))') fnmesh
c      write(*,*) nj1,nj2
c      if((nj1.eq.1.and.nj2.eq.1.).and.
c     & (fnmesh(1:1).eq.'0'.or.fnmesh(1:1).eq.'1')) go to 69
      open (unit=65,file=fnmesh)
      go to 68
   69 continue
      open (unit=65,file=moddir(msdir:medir)//'\'//'taumesh.500')
   68 continue

c      inquire(unit=5,opened=loopen)
c      if(lopen) close(5)
c      open (unit=5,file='models1.tmp')
      open (unit=5,file='models1.tmp',disp='delete')
      open (unit=19,file='models2.tmp')
      MC=0
      IF(IFMOLE.EQ.1) THEN
      IFPRES=1
      IFMOL=1
      inquire(unit=2,opened=loopen)
      if(lopen) close(2)
      open (unit=2, file=datdir(ndsdir:ndedir)//'\'//'molec.dat')
      REWIND 2
      CALL READMOL
      close (2)
      ELSE
      IFMOL=0
      END IF
      REWIND 65
      READ(65,*) N2
      DO J=1,N2
        READ(65,*) TAUL2(J)
      END DO
      close (65)
      write(6,*)
      write(6,*) '================ OUTPUT FROM MPREP ================'
      write(6,7) N2
    7 FORMAT(1X,'DEPTH POINTS ARE: ',I2,' WITH TAULOG(5000) SET OF:')
      write(6,8) (TAUL2(J),J=1,N2)
    8 FORMAT(1X,10F8.4)
      IF(IFMOL.EQ.1) write(6,*) 'MOLECULE FORMATION EFFECT INCLUDED...'
    1 CALL READIN(2)
      ITEMP=ITEMP+1
      MC=MC+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
      FREQ=2.997925E17/500.
      FREQLG=DLOG(FREQ)
      DO 30 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)
   30 CONTINUE
      N=0
      CALL KAPP(N,NSTEPS,STEPWT)
      DO 31 J=1,NRHOX
   31 AKAP5(J)=ACONT(J)+SIGMAC(J)
      CALL INTEG(RHOX,AKAP5,TAU5,NRHOX,AKAP5(1)*RHOX(1))
      N1=NRHOX-1
      DO 50 J=1,N1
      J1=J+1
      TAURL1(J)=DLOG10(TAUROS(J1))
      TAUL1(J)=DLOG10(TAU5(J1))
      RHOXL1(J)=DLOG10(RHOX(J1))
      T1(J)=T(J1)
      PL1(J)=DLOG10(P(J1))
      XNEL1(J)=DLOG10(XNE(J1))
      ABSL1(J)=DLOG10(ABROSS(J1))
      ABS5L1(J)=DLOG10(AKAP5(J1))
      ACCRL1(J)=DLOG10(ACCRAD(J1))
      VT1(J)=VTURB(J1)
      RHOL1(J)=DLOG10(RHO(J1))
   50 CONTINUE
cc      N1=NRHOX
cc      DO 50 J=1,N1
cc      TAURL1(J)=DLOG10(TAUROS(J))
cc      TAUL1(J)=DLOG10(TAU5(J))
cc      RHOXL1(J)=DLOG10(RHOX(J))
cc      T1(J)=T(J)
cc      PL1(J)=DLOG10(P(J))
cc      XNEL1(J)=DLOG10(XNE(J))
cc      ABSL1(J)=DLOG10(ABROSS(J))
cc      ABS5L1(J)=DLOG10(AKAP5(J))
cc      ACCRL1(J)=DLOG10(ACCRAD(J))
cc      VT1(J)=VTURB(J)
cc      RHOL1(J)=DLOG10(RHO(J))
cc   50 CONTINUE
      DO 20 J=1,N1
      DUMMY(J)=-1./(EXP10(RHOL1(J))*EXP10(ABS5L1(J)))
      TAUD(J)=EXP10(TAUL1(J))
   20 CONTINUE
      CALL INTEG(TAUD,DUMMY,HP1,N1,0.)
      IDUM=MAP1(TAUL1,HP1,N1,0.,HREF,1)
      DO 21 J=1,N1
   21 HP1(J)=HP1(J)-HREF
      IDUM=MAP1(TAUL1,TAURL1,N1,TAUL2,TAURL2,N2)
      IDUM=MAP1(TAUL1,T1,N1,TAUL2,T2,N2)
      IDUM=MAP1(TAUL1,PL1,N1,TAUL2,PL2,N2)
      IDUM=MAP1(TAUL1,XNEL1,N1,TAUL2,XNEL2,N2)
      IDUM=MAP1(TAUL1,ABSL1,N1,TAUL2,ABSL2,N2)
      IDUM=MAP1(TAUL1,ABS5L1,N1,TAUL2,ABS5L2,N2)
      IDUM=MAP1(TAUL1,ACCRL1,N1,TAUL2,ACCRL2,N2)
      IDUM=MAP1(TAUL1,VT1,N1,TAUL2,VT2,N2)
      IDUM=MAP1(TAUL1,RHOL1,N1,TAUL2,RHOL2,N2)
      DO 40 J=1,N2
      DUMMY(J)=-1./(EXP10(RHOL2(J))*EXP10(ABS5L2(J)))
      TAUD(J)=EXP10(TAUL2(J))
   40 CONTINUE
      CALL INTEG(TAUD,DUMMY,HP2,N2,0.)
      IDUM=MAP1(TAUL2,HP2,N2,0.,HREF,1)
      DO 41 J=1,N2
   41 HP2(J)=HP2(J)-HREF
      DO 59 J=1,N2
      DUMMY(J)=1./(GRAV-EXP10(ACCRL2(J)))
      PRDM(J)=EXP10(PL2(J))
   59 CONTINUE
      CALL INTEG(PRDM,DUMMY,RHOXL2,N2,EXP10(TAUL2(1)-ABS5L2(1)))
      DO 51 J=1,N2
   51 RHOXL2(J)=DLOG10(RHOXL2(J))
      write(6,*)
      write(6,'(''<<<<< MODEL'',I2,'' >>>>>'')') MC
      write(6,*) 'STANDARD OPTICAL DEPTH ... TAU(5000)'
      write(6,500) NRHOX,N1
  500 FORMAT(1H ,'PREVIOUS MODEL-MESH WITH NRHOX=',I2,2X,'N1=',I2)
      write(6,501)
  501 FORMAT(1H ,' J',2X,'  RHOXL ',2X,'  TAUL5 ',2X,'   T   ',2X,
     #'   PL   ',2X,'  XNEL  ',2X,'  ABSL  ',2X,'  ACCRL ',2X,
     #'    VT   ',2X,'   RHOL  ',2X,'     HP     ',2X,
     #'  TAULR  ',2X,'  ABSL5  ')
      DO 505 J=1,N1
      J1=J+1
  505 write(6,510)J1,RHOXL1(J),TAUL1(J),T1(J),PL1(J),XNEL1(J),ABSL1(J),
     #ACCRL1(J),VT1(J),RHOL1(J),HP1(J),TAURL1(J),ABS5L1(J)
  510 FORMAT(1H ,I2,2X,0PF8.4,2X,0PF8.4,2X,0PF7.1,4(2X,0PF8.4),
     #2X,1PE9.3,2X,0PF9.4,2X,1PE12.5,2(2X,0PF9.4))
      write(6,600) N2
  600 FORMAT(1H ,' ADOPTED MODEL-MESH WITH N2=',I3)
      write(6,501)
      DO 605 J=1,N2
  605 write(6,510) J,RHOXL2(J),TAUL2(J),T2(J),PL2(J),XNEL2(J),ABSL2(J),
     #ACCRL2(J),VT2(J),RHOL2(J),HP2(J),TAURL2(J),ABS5L2(J)
C
C     OUTPUT MODEL DATA TO FILE 11                          
      IF(IFCONV.EQ.0) THEN
      A='OFF'
      ELSE
      A='ON '
      END IF
      IF(IFTURB.EQ.0) THEN
      B='OFF'
      ELSE
      B='ON '
      END IF
      WRITE(19,1171) TEFF,GLOG,WLTE,TITLE,IFOP,A,MIXLTH,B,TRBFDG,               
     1TRBPOW,TRBSND,TRBCON
 1171 FORMAT(5HTEFF ,F7.0,9H  GRAVITY,F8.5,1X,A4/6HTITLE ,74A1/                 
     1 13H OPACITY IFOP,20I2/                                                   
     2 12H CONVECTION ,A3,F6.2,12H TURBULENCE , A3,4F6.2)                       
      XSL=LOG10(XSCALE)
      IF(XSL.GE.-2.) THEN
      WRITE(19,1181) XSCALE,(IZ,ABUND(IZ),IZ=1,2)
      ELSE IF (XSL.GE.-3.) THEN
      WRITE(19,1182) XSCALE,(IZ,ABUND(IZ),IZ=1,2)
      ELSE IF (XSL.GE.-4.) THEN
      WRITE(19,1183) XSCALE,(IZ,ABUND(IZ),IZ=1,2)
      ELSE
      WRITE(19,1184) XSCALE,(IZ,ABUND(IZ),IZ=1,2)
      END IF
 1181 FORMAT(16HABUNDANCE SCALE ,F9.5,17H ABUNDANCE CHANGE,2(I2,F8.5))                                    
 1182 FORMAT(16HABUNDANCE SCALE ,F9.6,17H ABUNDANCE CHANGE,2(I2,F8.5))                                        
 1183 FORMAT(16HABUNDANCE SCALE ,F9.7,17H ABUNDANCE CHANGE,2(I2,F8.5))                                        
 1184 FORMAT(16HABUNDANCE SCALE ,F9.8,17H ABUNDANCE CHANGE,2(I2,F8.5))                                        
      WRITE(19,1190) (IZ,ABUND(IZ),IZ=3,99)                        
 1190 FORMAT((17H ABUNDANCE CHANGE,6(I3,F7.2)))    
      GRAV=10.**GLOG                                                            
      FLUX=5.6697E-5/12.5664*TEFF**4                                            
      IF(WLTE.EQ.4HLTE )NLTEON=0                                                
      IF(WLTE.EQ.4HNLTE)NLTEON=1                                                
      IF(A.EQ.'OFF')IFCONV=0                                                    
      IF(A.EQ.'ON ')IFCONV=1                                                    
      IF(B.EQ.'OFF')IFTURB=0                                                    
      IF(B.EQ.'ON ')IFTURB=1                                                    
      WRITE(19,1172)N2,(EXP10(RHOXL2(J)),T2(J),EXP10(PL2(J)),        
     %EXP10(XNEL2(J)),EXP10(ABSL2(J)),EXP10(ACCRL2(J)),        
     %VT2(J),J=1,N2)                                                       
 1172 FORMAT(10HREAD DECK6,I3,33H RHOX,T,P,XNE,ABROSS,ACCRAD,VTURB/             
     1(1PE15.8,0PF9.1,1P5E10.3))                                                
      WRITE(19,1173)PRADK0                                                      
 1173 FORMAT(5HPRADK,1PE11.4)                                                   
      CALL INTEG(RHOX,ACCRAD,PRAD,NRHOX,ACCRAD(1)*RHOX(1))                      
      DO 1177 J=1,NRHOX                                                         
 1177 PRADK(J)=PRAD(J)+PRADK0                                                   
      IF(NLTEON.EQ.1)THEN                                                       
      write(6,*) 'NLTE MODEL ARE NOT ALLOWED IN THIS VERSION!'
      STOP
      ENDIF                                                                     
      WRITE(19,1175)                                                     
 1175 FORMAT(5HBEGIN,20X,10HITERATION ,3X,10H COMPLETED )                 
C
      GO TO 1
      END
      FUNCTION EXP10(X)
      IMPLICIT REAL*8 (A-H,O-Z)
      EXP10=EXP(X*2.30258509299405D0)
      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.'')) exit
      end do
      ne=n
      do n=1,nlen
      if((label(n:n).ne.' ').and.(label(n:n).ne.'')) exit
      end do
      ns=n
      return
      end
