      character*80 filein
      character*80 fileout
      logical exists
      character*80 string
c
      parameter (maxl=36)
      parameter (maxd=1)
      dimension vecmod((maxl+1)**2,maxd)
      dimension ylmcof((maxl+1)**2)
      dimension wk1(maxl+1)
      dimension wk2(maxl+1)
      dimension wk3(maxl+1)
      parameter (maxk=16)
      parameter (maxuknbig=maxk*(maxl+1)**2)
      dimension hrvmodel(maxuknbig)
      dimension maskk(maxk)
      dimension maskl(maxl+1)
      dimension chebyshev(maxk)
      parameter (r0=6371.)
      parameter (rmoho=6371.0-24.4)
      parameter (r670=6371.-670.)
      parameter (rcmb=3480.0)
      logical gotmodel
      character*128 model
      character*128 modeldef
      logical upper
      logical lower
      data modeldef/' '/
      save
c
      write(6,"('This program evaluates a mantle model expressed')")
      write(6,"('in terms of spherical harmonics and')")
      write(6,"('Chebyshev polynomials.')")
      write(6,"(' ')")
      write(6,"('The user provides an input file with latitude and')")
      write(6,"('longitude for a point on each line, and the')")
      write(6,"('program creates an output file with the model')")
      write(6,"('value added.')")
      write(6,"(' ')")
      write(6,"('The model value is a percentage perturbation w.r.t.')")
      write(6,"('the value of PREM at the give depth.')")
      write(6,"(' ')")
      write(6,"('Type depth (km) at which to evaluate ')")
      write(6,"('the model (greater than 24.4 km)')")
      read(5,*) depth
c
c---- loop on depths
c
      modeldef='S20A.SH'
      gotmodel=.false.
      do while(.not.gotmodel)
        write(6,"('change earth model [',a,']:')") 
     #      modeldef(1:lnblnk(modeldef))
        read(5,"(a)") string
        if(lnblnk(string).ne.0) then
        else
          string=modeldef
        endif
        inquire(file=string,exist=exists)
        if(exists) then
             call rdhetmod(1,string,maxl,maxk,
     #            hrvmodel,kmax,lmax,maskk,maskl,0)
            if(ierror.eq.0) then
              model=string
              modeldef=model
              gotmodel=.true.
            endif
        else
          gotmodel=.false.
        endif
      enddo
      write(6,"('maximum angular order of model:',2i5)") lmax
c
c---- evaluate model at this depth
c
      upper=.false.
      lower=.false.
      radius=r0-depth
      if(radius.gt.rcmb.and.radius.lt.r670) then
        lower=.true.
      else if(radius.ge.r670.and.radius.lt.rmoho) then
        upper=.true.
      endif
      do i=1,(lmax+1)**2
        vecmod(i,1)=0.
      enddo
      if(upper) then
              u=(radius+radius-rmoho-r670)/(rmoho-r670)
              call chebyfun(u,13,chebyshev)
              do k=3,10
                j1=(k-1)*((maxl+1)**2)+1
                do j=1,(lmax+1)**2
                  vecmod(j,1)=
     #              vecmod(j,1)+100.*chebyshev(k-2)*hrvmodel(j1)
                j1=j1+1
                enddo
              enddo
      else if(lower) then
              u=(radius+radius-r670-rcmb)/(r670-rcmb)
c              write(6,"('lower mantle:',2f10.3)") radius,u
              call chebyfun(u,13,chebyshev)
              do k=11,16
                j1=(k-1)*((maxl+1)**2)+1
                do j=1,(lmax+1)**2
                  vecmod(j,1)=
     #              vecmod(j,1)+100.*chebyshev(k-10)*hrvmodel(j1)
                  j1=j1+1
                enddo
              enddo
      else
              write(6,"('neither upper nor lower mantle?',f10.2)") radius
              stop
      endif
c
      write(6,"('type name of file with latitudes and longitudes')")
      read(5,"(a)") filein
      open(12,file=filein,iostat=ios,status='old')
      write(fileout,"(a,'.',i4.4)") filein(1:lnblnk(filein)),nint(depth)
      write(6,"('output file:',a)") fileout(1:lnblnk(fileout))
      open(14,file=fileout)
      if(ios.eq.0) then
      else
        write(6,"('input file does not exist')")
        stop
      endif
      do while (ios.eq.0)
        read(12,*,iostat=ios) xlat,xlon
        if(ios.eq.0) then
            call ylm(xlat,xlon,lmax,ylmcof,wk1,wk2,wk3)
            value=0.
            do i=1,(lmax+1)**2
                value=value+ylmcof(i)*vecmod(i,1)
            enddo
            write(14,"(4f10.3)") xlat,xlon,value
        endif
      enddo
      close(14)
      close(12)
      end
c
 
      subroutine rdhetmod(lu,file,maxl,maxr,x,nstr,lmax,mask,lask,inorm)
c
c---- reads in a heterogeneous model description
c
      dimension x((maxl+1)**2,1),mask(1),lask(0:1)
      character*(*) file
      character*132 abuf
      character*80 form
c
      sqth=sqrt(.5)
c
      open(lu,file=file,status='old')
      read(lu,"(a132)") abuf
      llen=132
      do 30 i=1,132
      if(abuf(llen:llen).ne.' ') goto 40
   30 llen=llen-1
   40 read(abuf,"(10x,i5)") lmax
      write(6,*) lmax
c
      if(lmax.gt.maxl) then
        write(6,"('gethetmod not ready for this l',2i5)") maxl,lmax
        call exit(3)
      endif
c
      num=llen-15-1-(lmax+1)-1-3-1
c
      if(lmax.gt.maxl) then
        write(6,"('gethetmod not ready for this number of '
     #     'radial functions',2i5)") maxr,num
        call exit(3)
      endif

      if(llen.gt.15) then
        idef=0
        write(form,"('(16x,',i3,'i1,1x,i3,1x,80i1)')") lmax+1
        write(*,"(a30,2i8)") form(1:30),lmax,num
        read(abuf,form) (lask(i),i=0,lmax),nstr,(mask(i),i=1,num)
        if(num.ne.nstr) pause 'error 1 in heread'
      else
        idef=1
        do 50 i=0,lmax
        lask(i)=1
   50   continue
      endif

      leny=0
      do 60 l=0,lmax
      if(lask(l).ne.0) leny=leny+2*l+1
   60 continue

      ip=0
  100 read(lu,"(a132)",end=99) abuf
      backspace lu
      ip=1+ip
      if(idef.ne.0) mask(ip)=1

cc      ind=1+(ip-1)*leny
      ind=1
      do 704 l=0,lmax
      if(lask(l).ne.0) then
        ind1=ind+2*l
        read(lu,"(11e12.4)")(x(i,ip),i=ind,ind1)
        if(inorm.eq.1) x(ind,ip)=x(ind,ip)/sqth
        ind=ind1+1
      endif
  704 continue
      goto 100

   99 if(idef.ne.0) then
      nstr=ip
      nstruc=ip
      else
      nstruc=0
      do 67 i=1,nstr
      if(mask(i).ne.0) nstruc=1+nstruc
   67 continue
      if(nstruc.ne.ip) pause 'error 3 in heread'
      endif

      lenf=0
      do 22 i=1,len(file)
      ii=i
      if(file(i:i).eq.' '.or.file(i:i).eq.'<0>') goto 23
   22 continue
   23 ii=ii-1
      write(*,"('model from',1x,80a1)") (file(i:i),i=1,ii)
      write(*,"(i4,' parameters:  mask ',80i1)") nstr,(mask(i),i=1,nstr)
      write(*,"(4x,' lmax =',i4,':  mask ',80i1)") lmax,
     1(lask(i),i=0,lmax)
      close(lu)
      return
      end
      subroutine chebyfun(u,kmax,f)
      dimension chebycoeff(0:13),f(0:kmax)
      data  chebycoeff /
     . 0.70710678118655,1.2247448713916,1.0350983390135,1.0145993123918,
     . 1.00803225754840,1.0050890913907,1.0035149493262,1.0025740068320,
     . 1.00196657023780,1.0015515913133,1.0012554932754,1.0010368069141,
     . 1.00087070107920,1.0007415648034 /
       
      if(kmax.gt.13)then
         write(*,"(' kmax exceeds the limit in chebyfun')")
         stop
      endif
       
      f(0)=1.0
      f(1)=u
      twou=2.0*u
       
      do k=2,kmax
         f(k) = twou*f(k-1)-f(k-2)
      enddo
       
      do k=0,kmax
         f(k)=f(k)*chebycoeff(k)
      enddo
       
      return
      end
      SUBROUTINE YLM(XLAT,XLON,LMAX,Y,WK1,WK2,WK3)
C
      COMPLEX TEMP,FAC,DFAC
      DIMENSION WK1(1),WK2(1),WK3(1),Y(1)
C
C     WK1,WK2,WK3 SHOULD BE DIMENSIONED AT LEAST (LMAX+1)*4
C
      DATA RADIAN/57.2957795/
C
      THETA=(90.-XLAT)/RADIAN
      PHI=XLON/RADIAN
C
      IND=0
      LM1=LMAX+1
      DO 10 IL1=1,LM1
      L=IL1-1
      CALL LEGNDR(THETA,L,L,WK1,WK2,WK3)
C
      FAC=(1.,0.)
      DFAC=CEXP(CMPLX(0.,PHI))
C
      DO 20 IM=1,IL1
      TEMP=FAC*CMPLX(WK1(IM),0.)
      IND=IND+1
      Y(IND)=REAL(TEMP)
      IF(IM.EQ.1) GOTO 20
      IND=IND+1
      Y(IND)=AIMAG(TEMP)
   20 FAC=FAC*DFAC
C
   10 CONTINUE
      RETURN
      END
c
      SUBROUTINE LEGNDR(THETA,L,M,X,XP,XCOSEC)
      DIMENSION X(2),XP(2),XCOSEC(2)
      DOUBLE PRECISION SMALL,SUM,COMPAR,CT,ST,FCT,COT,FPI,X1,X2,X3,
     1F1,F2,XM,TH,DFLOAT
      DATA FPI/12.56637062D0/
      DFLOAT(I)=FLOAT(I)
      SUM=0.D0
      LP1=L+1
      TH=THETA
      CT=DCOS(TH)
      ST=DSIN(TH)
      MP1=M+1
      FCT=DSQRT(DFLOAT(2*L+1)/FPI)
      SFL3=SQRT(FLOAT(L*(L+1)))
      COMPAR=DFLOAT(2*L+1)/FPI
      DSFL3=SFL3
      SMALL=1.D-16*COMPAR
      DO 1 I=1,MP1
      X(I)=0.
      XCOSEC(I)=0.
    1 XP(I)=0.
      IF(L.GT.1.AND.ABS(THETA).GT.1.E-5) GO TO 3
      X(1)=FCT
      IF(L.EQ.0) RETURN
      X(1)=CT*FCT
      X(2)=-ST*FCT/DSFL3
      XP(1)=-ST*FCT
      XP(2)=-.5D0*CT*FCT*DSFL3
      IF(ABS(THETA).LT.1.E-5) XCOSEC(2)=XP(2)
      IF(ABS(THETA).GE.1.E-5) XCOSEC(2)=X(2)/ST
      RETURN
    3 X1=1.D0
      X2=CT
      DO 4 I=2,L
      X3=(DFLOAT(2*I-1)*CT*X2-DFLOAT(I-1)*X1)/DFLOAT(I)
      X1=X2
    4 X2=X3
      COT=CT/ST
      COSEC=1./ST
      X3=X2*FCT
      X2=DFLOAT(L)*(X1-CT*X2)*FCT/ST
      X(1)=X3
      X(2)=X2
      SUM=X3*X3
      XP(1)=-X2
      XP(2)=DFLOAT(L*(L+1))*X3-COT*X2
      X(2)=-X(2)/SFL3
      XCOSEC(2)=X(2)*COSEC
      XP(2)=-XP(2)/SFL3
      SUM=SUM+2.D0*X(2)*X(2)
      IF(SUM-COMPAR.GT.SMALL) RETURN
      X1=X3
      X2=-X2/DSQRT(DFLOAT(L*(L+1)))
      DO 5 I=3,MP1
      K=I-1
      F1=DSQRT(DFLOAT((L+I-1)*(L-I+2)))
      F2=DSQRT(DFLOAT((L+I-2)*(L-I+3)))
      XM=K
      X3=-(2.D0*COT*(XM-1.D0)*X2+F2*X1)/F1
      SUM=SUM+2.D0*X3*X3
      IF(SUM-COMPAR.GT.SMALL.AND.I.NE.LP1) RETURN
      X(I)=X3
      XCOSEC(I)=X(I)*COSEC
      X1=X2
      XP(I)=-(F1*X2+XM*COT*X3)
    5 X2=X3
      RETURN
      end
