      PROGRAM MDINTERP
      USE DFLIB
      use dfport
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (kw=99)
      DIMENSION TAU(kw),TAULOG(kw),DUMR(kw,2),DUMK(kw,2),xsc(2,2)
      COMMON /PASS1/ ND,I1,J1,nn
      COMMON /PASS2/ TEFF(02),GLOG(2),
     &RHX(kw,02,2),RHOXS(kw,2),
     &TMP(kw,02,2),TEMPE(kw,2),
     &PRS(kw,02,2),PRESS(kw,2),
     &XNE(kw,02,2),XNEDS(kw,2),
     &ABK(kw,02,2),ABKAP(kw,2),
     &ACR(kw,02,2),ACRAD(kw,2),
     &VTU(kw,02,2),VTURB(kw,2),
     &DUM(3,02,2),DUMMY(3,2)
      logical loexist
      INTEGER*2 NA
      CHARACTER ARGV*60, fnmesh*60
      CHARACTER*17 FNAME
      CHARACTER*80 LINE(30),linedir,moddir
      character*6 ziffer(2)
      character labte*5,labgl*3,labxl*3,labpm*1,labmol*1
      EXP10(X)=EXP(2.3025851E0*X)
c      open (unit=6,file='mdintout.out')
      open (unit=6,file='mdintout.out',disp='delete')
      NA=1
      CALL GETARG(NA,ARGV)
      READ(ARGV,*) TE
      ite=int(te+0.5)
      write(labte,'(i5)') ite
      if(ite.lt.10000) labte(1:1)='0'
      NA=2
      CALL GETARG(NA,ARGV)
      READ(ARGV,*) GL
      if(gl.lt.0.) then
      write(*,*) 'Negative gl not allowed!'
      stop
      end if
      igl=int(gl*100.+0.5)
      write(labgl,'(i3)') igl
      if(igl.lt.100) labgl(1:1)='0'
      if(igl.lt.10)  labgl(2:2)='0'
      NA=3
      CALL GETARG(NA,ARGV)
      READ(ARGV,*) XSLOG
      ixl=int(abs(xslog)*100.+0.5)
      write(labxl,'(i3)') ixl
      if(ixl.lt.100) labxl(1:1)='0'
      if(ixl.lt.10)  labxl(2:2)='0'
      if(xslog.ge.0.) then
        labpm='+'
      else
        labpm='-'
      end if
      NA=4
      CALL GETARG(NA,ARGV)
      READ(ARGV,*) IFMOLE
      write(labmol,'(i1)') ifmole
      write(6,'(''te:'',f7.0)') te
      write(6,'(''gl:'',f7.2)') gl
      write(6,'(''xl:'',f7.2)') xslog
      fname='t'//labte//'g'//labgl//'x'//labpm//labxl//'m'//labmol
      write(6,'(''fn:'',a17)') fname
c      if(xslog.ge.-5.0.and.xslog.lt.-4.5) then
      if(xslog.gt.-5.5.and.xslog.lt.-4.5) then
      ziffer(1)='am50k2'
      ziffer(2)='am45k2'
      xt=(xslog-(-5.0))/((-4.5)-(-5.0))
      else if (xslog.ge.-4.5.and.xslog.lt.-4.0) then
      ziffer(1)='am45k2'
      ziffer(2)='am40k2'
      xt=(xslog-(-4.5))/((-4.0)-(-4.5))
      else if (xslog.ge.-4.0.and.xslog.lt.-3.5) then
      ziffer(1)='am40k2'
      ziffer(2)='am35k2'
      xt=(xslog-(-4.0))/((-3.5)-(-4.0))
      else if (xslog.ge.-3.5.and.xslog.lt.-3.0) then
      ziffer(1)='am35k2'
      ziffer(2)='am30k2'
      xt=(xslog-(-3.5))/((-3.0)-(-3.5))
      else if (xslog.ge.-3.0.and.xslog.lt.-2.5) then
      ziffer(1)='am30k2'
      ziffer(2)='am25k2'
      xt=(xslog-(-3.0))/((-2.5)-(-3.0))
      else if (xslog.ge.-2.5.and.xslog.lt.-2.0) then
      ziffer(1)='am25k2'
      ziffer(2)='am20k2'
      xt=(xslog-(-2.5))/((-2.0)-(-2.5))
      else if (xslog.ge.-2.0.and.xslog.lt.-1.5) then
      ziffer(1)='am20k2'
      ziffer(2)='am15k2'
      xt=(xslog-(-2.0))/((-1.5)-(-2.0))
      else if (xslog.ge.-1.5.and.xslog.lt.-1.0) then
      ziffer(1)='am15k2'
      ziffer(2)='am10k2'
      xt=(xslog-(-1.5))/((-1.0)-(-1.5))
      else if (xslog.ge.-1.0.and.xslog.lt.-0.5) then
      ziffer(1)='am10k2'
      ziffer(2)='am05k2'
      xt=(xslog-(-1.0))/((-0.5)-(-1.0))
      else if (xslog.ge.-0.5.and.xslog.lt. 0.0) then
      ziffer(1)='am05k2'
      ziffer(2)='ap00k2'
      xt=(xslog-(-0.5))/(( 0.0)-(-0.5))
      else if (xslog.ge. 0.0.and.xslog.lt. 0.5) then
      ziffer(1)='ap00k2'
      ziffer(2)='ap05k2'
      xt=(xslog-( 0.0))/(( 0.5)-( 0.0))
      else if (xslog.ge. 0.5.and.xslog.lt. 1.0) then
      ziffer(1)='ap05k2'
      ziffer(2)='ap10k2'
      xt=(xslog-( 0.5))/(( 1.0)-( 0.5))
      else
      write(6,'( ''xslog='',f9.5,'' is out of range!'')') xslog
      stop
      end if
c     get MOD_DIR
c      call getenv("MOD_DIR",moddir)
      call retenv("MOD_DIR",moddir)
c      open (88,file='passdir.tmp',disp='delete')
c      read(88,'(a)') linedir
c      read(88,'(a)') moddir
c      close(88)
c      na=5
c      call getarg(na,argv)
c      read(argv,'(a)') moddir
      call lablen(moddir,80,msdir,medir)


      NA=5
      CALL GETARG(NA,ARGV)
      read(argv,'(a)',err=69,end=69) fnmesh
      inquire(file=fnmesh,exist=loexist)
      if(.not.loexist) go to 69
c      write(*,'(''fnmesh:'',(a))') fnmesh
c      write(*,*) nj1,nj2
c      call lablen(fnmesh,60,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=75,file=fnmesh)
      go to 68
   69 continue
      open (unit=75,file=moddir(msdir:medir)//'\'//'taumesh.500')
   68 continue
      rewind 75
      READ(75,*) ND
      DO J=1,ND
        READ(75,*) TAULOG(J)
      END DO
      close(75)
      DO 4 J=1,ND
    4 TAU(J)=EXP10(TAULOG(J))
c      OPEN (UNIT=70, FILE='models2.tmp')
      OPEN (UNIT=70, FILE='models2.tmp',disp='delete')
      REWIND 70
      do 8888 nn=1,2

      DO 555 IT=1,2
      DO 555 IG=1,2
      READ(70,170) TEF,GRA,WLTE
  170 FORMAT(5X,F7.0,9X,F8.5,1X,A4)
      DUM(1,IT,IG)=TEF
      DUM(2,IT,IG)=GRA
      DO 777 K=1,30
      READ(70,'(A80)') LINE(K)
      IF(LINE(K)(1:4).EQ.'READ') GO TO 779
  777 CONTINUE
  779 BACKSPACE 70
      read(line(4),'(16x,f10.5)') xsc(it,ig)
      KM=K-1
      READ(70,176) ND
  176 FORMAT(10X,I3)
      DO 256 J=1,ND
      READ(70,177) RHX(J,IT,IG),TMP(J,IT,IG),PRS(J,IT,IG),
     &XNE(J,IT,IG),ABK(J,IT,IG),ACR(J,IT,IG),VTU(J,IT,IG)
  177 FORMAT(E15.8,F9.1,5E10.3)                                        
  256 CONTINUE
      DO 29 J=1,ND
      RHX(J,IT,IG)=LOG10(RHX(J,IT,IG))
      PRS(J,IT,IG)=LOG10(PRS(J,IT,IG))
      XNE(J,IT,IG)=LOG10(XNE(J,IT,IG))
      ABK(J,IT,IG)=LOG10(ABK(J,IT,IG))
      ACR(J,IT,IG)=LOG10(ACR(J,IT,IG))
   29 CONTINUE
      READ(70,1173) PRD                                             
 1173 FORMAT(5X,E11.4)                                          
      DUM(3,IT,IG)=LOG10(PRD)
      READ(70,1175)                                                     
 1175 FORMAT(48X)                
  555 CONTINUE
C
      TEFF(1)=DUM(1,1,1)
      TEFF(2)=DUM(1,2,1)
      GLOG(1)=DUM(2,1,1)
      GLOG(2)=DUM(2,1,2)
      I1=1
      J1=1
      CALL DINPOL(TE,GL)
      WRITE(6,*)'================ NN=',nn,' ================'
      write(6,179) xsc(1,1),xsc(2,1),xsc(1,2),xsc(2,2)
  179 format(1h ,'xsc:',4f12.5)
      WRITE(6,180) TE,TEFF(I1),TEFF(I1+1)
  180 FORMAT(1H ,'  TE=',F7.0,3X,'TEFF1=',F7.0,3X,'TEFF2=',F7.0)
      WRITE(6,190) GL,GLOG(J1),GLOG(J1+1)
  190 FORMAT(1H ,'  GL=',F5.2,5X,'GLOG1=',F5.2,5X,'GLOG2=',F5.2)
      WRITE(6,*) 'T/G INTERPOLATION COMPLETED'

 8888 continue

      CLOSE (UNIT=70)

      WRITE(6,'(''MODEL OUTPUT TO FILE:'',A17)') FNAME
      OPEN (30,FILE=FNAME)
c
CCC The following three lines may be useful foe evaluation of RHOX
c   in low-gravity stars (i.e., supergiants near the Eddington limit)
c      DO 1801 J=1,ND
c 1801 DUMK(J,nn)=1./EXP10(ABKAP(J))
c      CALL INTEG(TAU,DUMK(1,nn),DUMR(1,nn),ND,EXP10(RHOXS(1)))
CCC
C
      WRITE(30,2171) TE,GL,WLTE           
 2171 FORMAT(5HTEFF ,F7.0,9H  GRAVITY,F8.5,1X,A4)         
      write(30,'(''TITLE (T/G/M)-INTERPOLATED MODEL [BY IPOL.BAT]'')') 
      DO 1850 K=2,3
      call lablen(line(k),80,ns,ne)
 1850 WRITE(30,'(A)') LINE(K)(1:ne)
      call lablen(line(4),80,ns,ne)
      xs=exp10(xslog)
      if(xslog.ge.-2.) then
        write(30,1890) xs,line(4)(26:ne)
      else if(xslog.ge.-3.) then
        write(30,1891) xs,line(4)(26:ne)
      else if(xslog.ge.-4.) then
        write(30,1892) xs,line(4)(26:ne)
      else
        write(30,1893) xs,line(4)(26:ne)
      end if
 1890 format('ABUNDANCE SCALE ',f9.5,(a))
 1891 format('ABUNDANCE SCALE ',f9.6,(a))
 1892 format('ABUNDANCE SCALE ',f9.7,(a))
 1893 format('ABUNDANCE SCALE ',f9.8,(a))
      DO 1851 K=5,KM
      call lablen(line(k),80,ns,ne)
 1851 WRITE(30,'(A)') LINE(K)(1:ne)
c
CCC Replace with these if necesary (e.g., for the case of supergiants)
c      WRITE(30,2172)
c       ND,
c     %((1.-xt)*DUMK(J,1)+xt*DUMK(J,2),
c     %(1.-xt)*TEMPE(J,1)+xt*TEMPE(J,2),
c     %EXP10((1.-xt)*PRESS(J,1)+xt*PRESS(J,2)),
c     %EXP10((1.-xt)*XNEDS(J,1)+xt*XNEDS(J,2)),
c     %EXP10((1.-xt)*ABKAP(J,1)+xt*ABKAP(J,2)),
c     %EXP10((1.-xt)*ACRAD(J,1)+xt*ACRAD(J,2)),
c     %(1.-xt)*VTURB(J,1)+xt*VTURB(J,2),J=1,ND)                                                 
CCC
c
      WRITE(30,2172)
     % ND,
     % (EXP10((1.-xt)*RHOXS(J,1)+xt*RHOXS(J,2)),
     % (1.-xt)*TEMPE(J,1)+xt*TEMPE(J,2),
     % EXP10((1.-xt)*PRESS(J,1)+xt*PRESS(J,2)),
     % EXP10((1.-xt)*XNEDS(J,1)+xt*XNEDS(J,2)),
     % EXP10((1.-xt)*ABKAP(J,1)+xt*ABKAP(J,2)),
     % EXP10((1.-xt)*ACRAD(J,1)+xt*ACRAD(J,2)),
     % (1.-xt)*VTURB(J,1)+xt*VTURB(J,2),J=1,ND)
 2172 FORMAT(10HREAD DECK6,I3,33H RHOX,T,P,XNE,ABROSS,ACCRAD,VTURB/     
     1(1PE15.8,0PF9.1,1P5E10.3))                                        
      WRITE(30,2173) EXP10((1.-xt)*DUMMY(3,1)+xt*DUMMY(3,2))                                           
 2173 FORMAT(5HPRADK,1PE11.4)                                           
      WRITE(30,2175)                                                   
 2175 FORMAT(5HBEGIN,10X,'(INTERPOLATED MODEL)')                 
C
      CLOSE (UNIT=30)
      close (6)
  999 STOP
      END
      SUBROUTINE DINPOL(TE,GL)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (kw=99)
      COMMON /PASS1/ ND,I1,J1,nn
      COMMON /PASS2/ TEFF(02),GLOG(2),
     &RHX(kw,02,2),RHOXS(kw,2),
     &TMP(kw,02,2),TEMPE(kw,2),
     &PRS(kw,02,2),PRESS(kw,2),
     &XNE(kw,02,2),XNEDS(kw,2),
     &ABK(kw,02,2),ABKAP(kw,2),
     &ACR(kw,02,2),ACRAD(kw,2),
     &VTU(kw,02,2),VTURB(kw,2),
     &DUM(3,02,2),DUMMY(3,2)
      I2=I1+1
      J2=J1+1
      DT=TEFF(I2)-TEFF(I1)
      DG=GLOG(J2)-GLOG(J1)
      P=(TE-TEFF(I1))/DT
      Q=(GL-GLOG(J1))/DG
      DO 10 J=1,ND
      T0=RHX(J,I1,J1)
      T1=-RHX(J,I1,J1)+RHX(J,I2,J1)
      T2=-RHX(J,I1,J1)+RHX(J,I1,J2)
      T3=RHX(J,I1,J1)-RHX(J,I1,J2)-RHX(J,I2,J1)+RHX(J,I2,J2)
      RHOXS(J,nn)=T0+T1*P+T2*Q+T3*P*Q
   10 CONTINUE
      DO 20 J=1,ND
      T0=TMP(J,I1,J1)
      T1=-TMP(J,I1,J1)+TMP(J,I2,J1)
      T2=-TMP(J,I1,J1)+TMP(J,I1,J2)
      T3=TMP(J,I1,J1)-TMP(J,I1,J2)-TMP(J,I2,J1)+TMP(J,I2,J2)
      TEMPE(J,nn)=T0+T1*P+T2*Q+T3*P*Q
   20 CONTINUE
      DO 30 J=1,ND
      T0=PRS(J,I1,J1)
      T1=-PRS(J,I1,J1)+PRS(J,I2,J1)
      T2=-PRS(J,I1,J1)+PRS(J,I1,J2)
      T3=PRS(J,I1,J1)-PRS(J,I1,J2)-PRS(J,I2,J1)+PRS(J,I2,J2)
      PRESS(J,nn)=T0+T1*P+T2*Q+T3*P*Q
   30 CONTINUE
      DO 40 J=1,ND
      T0=XNE(J,I1,J1)
      T1=-XNE(J,I1,J1)+XNE(J,I2,J1)
      T2=-XNE(J,I1,J1)+XNE(J,I1,J2)
      T3=XNE(J,I1,J1)-XNE(J,I1,J2)-XNE(J,I2,J1)+XNE(J,I2,J2)
      XNEDS(J,nn)=T0+T1*P+T2*Q+T3*P*Q
   40 CONTINUE
      DO 50 J=1,ND
      T0=ABK(J,I1,J1)
      T1=-ABK(J,I1,J1)+ABK(J,I2,J1)
      T2=-ABK(J,I1,J1)+ABK(J,I1,J2)
      T3=ABK(J,I1,J1)-ABK(J,I1,J2)-ABK(J,I2,J1)+ABK(J,I2,J2)
      ABKAP(J,nn)=T0+T1*P+T2*Q+T3*P*Q
   50 CONTINUE
      DO 60 J=1,ND
      T0=ACR(J,I1,J1)
      T1=-ACR(J,I1,J1)+ACR(J,I2,J1)
      T2=-ACR(J,I1,J1)+ACR(J,I1,J2)
      T3=ACR(J,I1,J1)-ACR(J,I1,J2)-ACR(J,I2,J1)+ACR(J,I2,J2)
      ACRAD(J,nn)=T0+T1*P+T2*Q+T3*P*Q
   60 CONTINUE
      DO 70 J=1,ND
      T0=VTU(J,I1,J1)
      T1=-VTU(J,I1,J1)+VTU(J,I2,J1)
      T2=-VTU(J,I1,J1)+VTU(J,I1,J2)
      T3=VTU(J,I1,J1)-VTU(J,I1,J2)-VTU(J,I2,J1)+VTU(J,I2,J2)
      VTURB(J,nn)=T0+T1*P+T2*Q+T3*P*Q
   70 CONTINUE
      DO 80 J=1,3
      T0=DUM(J,I1,J1)
      T1=-DUM(J,I1,J1)+DUM(J,I2,J1)
      T2=-DUM(J,I1,J1)+DUM(J,I1,J2)
      T3=DUM(J,I1,J1)-DUM(J,I1,J2)-DUM(J,I2,J1)+DUM(J,I2,J2)
      DUMMY(J,nn)=T0+T1*P+T2*Q+T3*P*Q
   80 CONTINUE
      RETURN
      END
      SUBROUTINE INTEG(X,F,FINT,N,START)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (kw=99)
C     SUBROUTINE INTEG(X,F,FINT,N)
      DIMENSION X(1),F(1),FINT(1)
      DIMENSION A(kw),B(kw),C(kw)
      CALL PARCOE(F,X,A,B,C,N)
      FINT(1)=START
C     FINT(1)=(A(1)+(B(1)/2.+C(1)/3.*X(1))*X(1))*X(1)
C     FINT(2)=(A(1)+(B(1)/2.+C(1)/3.*X(2))*X(2))*X(2)
C     IF(N.EQ.2)RETURN
      N1=N-1
C     DO 10 I=2,N1
      DO 10 I=1,N1
   10 FINT(I+1)=FINT(I)+(A(I)+B(I)/2.*(X(I+1)+X(I))+
     1C(I)/3.*((X(I+1)+X(I))*X(I+1)+X(I)*X(I)))*(X(I+1)-X(I))
      RETURN
      END
      SUBROUTINE PARCOE(F,X,A,B,C,N)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (kw=99)
      DIMENSION F(1),X(1),A(1),B(1),C(1)
      C(1)=0.
      B(1)=(F(2)-F(1))/(X(2)-X(1))
      A(1)=F(1)-X(1)*B(1)
      N1=N-1
      C(N)=0.
      B(N)=(F(N)-F(N1))/(X(N)-X(N1))
      A(N)=F(N)-X(N)*B(N)
      IF(N.EQ.2)RETURN
      DO 1 J=2,N1
      J1=J-1
      D=(F(J)-F(J1))/(X(J)-X(J1))
      C(J)=F(J+1)/((X(J+1)-X(J))*(X(J+1)-X(J1)))-F(J)/((X(J)-X(J1))*
     1(X(J+1)-X(J)))+F(J1)/((X(J)-X(J1))*(X(J+1)-X(J1)))
      B(J)=D-(X(J)+X(J1))*C(J)
    1 A(J)=F(J1)-X(J1)*D+X(J)*X(J1)*C(J)
      C(2)=0.
      B(2)=(F(3)-F(2))/(X(3)-X(2))
      A(2)=F(2)-X(2)*B(2)
      C(3)=0.
      B(3)=(F(4)-F(3))/(X(4)-X(3))
      A(3)=F(3)-X(3)*B(3)
      DO 2 J=2,N1
      IF(C(J).EQ.0.)GO TO 2
      J1=J+1
      WT=ABS(C(J1))/(ABS(C(J1))+ABS(C(J)))
      A(J)=A(J1)+WT*(A(J)-A(J1))
      B(J)=B(J1)+WT*(B(J)-B(J1))
      C(J)=C(J1)+WT*(C(J)-C(J1))
    2 CONTINUE
      A(N1)=A(N)
      B(N1)=B(N)
      C(N1)=C(N)
      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

      subroutine retenv(varname,string)
      use dfport
      character*7 varname
      character*80 sptool,string
      call getenv('SPTOOL',sptool)
      call lablen(sptool,80,is,ie)
      call lablen(varname,7,ivs,ive)
      if(varname(ivs:ive).eq.'SPTOOL') then
          string=sptool(is:ie)
        else if(varname(ivs:ive).eq.'SPT_DIR') then
          string=sptool(is:ie)//'\exefiles'
        else if(varname(ivs:ive).eq.'MOD_DIR') then
          string=sptool(is:ie)//'\k2models'
        else if(varname(ivs:ive).eq.'KGF_DIR') then
          string=sptool(is:ie)//'\GF10'
        else if(varname(ivs:ive).eq.'SFX_DIR') then
          string=sptool(is:ie)//'\sflxkpno'
        else if(varname(ivs:ive).eq.'DAT_DIR') then
          string=sptool(is:ie)//'\miscdata'
        else
          write(*,*) 'I do not understand '//varname(ivs:ive)
          stop
      end if
      return
      end
