      PROGRAM MULTI
C
C  MAIN CALLING PROGRAM FOR SOLVING MULTI-LEVEL NON-LTE PROBLEMS
C  FOR DOCUMENTATION, COMMENTS ON EXPERIENCE WITH THE CODE AND
C  A DESCRIPTION OF THE METHOD SEE
C  CARLSSON 1986: 'UPPSALA ASTRONOMICAL OBSERVATORY, REPORT NO. 33'.
C  A MORE DETAILED DESCRIPTION OF THE METHOD IS IN
C  SCHARMER, CARLSSON 1985: 'JOURNAL OF COMPUTATIONAL PHYSICS',59,56
C  FOR LIST OF VARIABLES AND FILES, SEE  ROUTINE  START.
C
C  PARAMETER VARIABLES:
C  MDEP1  NUMBER OF DEPTH POINTS. MUST BE THE NUMBER ACTUALLY USED.
C  MK1    NUMBER OF LEVELS INCLUDING CONTINUUM LEVELS. -''-.
C  MDEP   MAXIMUM NUMBER OF DEPTH POINTS.
C  MK     MAXIMUM NUMBER OF LEVELS INCLUDING CONTINUUM LEVELS
C  MLINE  MAXIMUM NUMBER OF BOUND-BOUND LINES
C  MWIDE  MAXIMUM NUMBER OF BROAD LINES + CONTINUA IN DETAIL
C  MRAD   MLINE+MWIDE
C  MRFIX  MAXIMUM NUMBER OF FIXED TRANSITIONS
C  MQ     MAXIMUM NUMBER OF FREQUENCY POINTS IN ONE TRANSITION
C  MMU    MAXIMUM NUMBER OF ANGLE POINTS
C
C  ONLY MDEP1 AND MK1 MUST BE EQUAL TO THE VALUES ACTUALLY USED, NDEP1
C  AND NK1. THE OTHERS MAY BE LARGER THAN VALUES ACTUALLY USED.
C
C: MULTI  90-07-31  MODIFICATIONS: (MATS CARLSSON)
C:        CHANGED TO CONTAIN DECLARATION OF W AND E
C:        THESE PASSED IN ARGUMENT LIST INSTEAD OF IN COMMON TO
C:        AVOID NK=MK1 AND NDEP=MDEP1 REQUIREMENT
C:
      INCLUDE 'PREC'
      INCLUDE 'PARAM'
      INCLUDE 'PARAMW'
      INCLUDE 'CATOM'
      INCLUDE 'CATMOS'
      INCLUDE 'CATMO2'
      INCLUDE 'CTRAN'
      INCLUDE 'CSLINE'
      INCLUDE 'CGAUSI'
      INCLUDE 'CCONST'
      INCLUDE 'CINPUT'
      INCLUDE 'CLGMX'
      INCLUDE 'CLU'
      INCLUDE 'CMPI'
      INCLUDE 'COPCL'
      INCLUDE 'C3D'
C
      DIMENSION W(MK1,MDEP1,MK1,MDEP1),E(MK1,MDEP1)
      integer, parameter :: dp=kind(1.0d0)
      real(4), dimension(:,:,:), allocatable :: temp3d,ne3d,vz3d,rho3d
      real(4), dimension(:,:,:,:), allocatable :: n3d
      real(4), dimension(:), allocatable :: widthx,widthy,height_in
      real(4), dimension(:,:,:), allocatable :: taulg3d,cmasslg3d
      real(4), dimension(:,:,:), allocatable :: xnorm3d
      real(4), dimension(:,:,:), allocatable :: height3d
      real(4), dimension(:,:,:), allocatable :: dpnew3d
      real(4), dimension(:,:,:), allocatable :: x3d
      real(4), dimension(:,:,:), allocatable :: icont
      real(8), dimension(:), allocatable :: xl
      real(8) :: qn
      real(4), dimension(:), allocatable :: buf,buf2
      integer, dimension(2) :: ibuf
      integer, dimension(:,:), allocatable :: iredo,idxy
      real(4) :: dummy
      integer :: nu0,itag,itag_base
      character :: cdummy
      dimension dpnew(mdep)
      external string3
      character*3 string3
      character*8 cname
      character*10 cmpmyid
      character*132 file2
      integer, parameter :: mpmproc=1000
      integer :: icxy(0:mpmproc),ncxy            ! number of columns per processor and max number
      logical debug,more_to_receive(0:mpmproc)
      integer nrad_write,kr_write(mrad),nq_write(mrad)
      logical lkr_write(mrad)
c
c      debug=.false.
      debug=.true.
c
c  init MPI
c
      call mpi_init(mperr)
      call mpi_comm_size(MPI_COMM_WORLD,mpnproc,mperr)
      call mpi_comm_rank(MPI_COMM_WORLD,mpmyid,mperr)
      master=mpmyid.eq.masterid
      slave=.not.master
      if(mpnproc.gt.mpmproc) then
        call stop('main: mpnproc.gt.mpmproc')
      endif
      do i=0,mpnproc
        more_to_receive(i)=.true.
        icxy(i)=0
      enddo
c
c  read parameters from INPUT_3D file
c  new format with many runs in same input file
c  bext read moved. Old format will cause read error
c  changed to give start and stop index and step
c  if krsave is negative, the number of transitions to write is
c  read and their indices. logical array lkr_write keeps track
c
      call mopen(linput_3d,'INPUT_3D',1,'OLD')
      read(linput_3d,*) icont_opt,iredo_max
      read(linput_3d,*) dscopt
      read(linput_3d,*) ix0,ix1,ixstep
      read(linput_3d,*) iy0,iy1,iystep
      read(linput_3d,*) krsave,nusave
      if(krsave.lt.0) then
        krsave=abs(krsave)
        read(linput_3d,*) nrad_write
        read(linput_3d,*) (kr_write(i),i=1,nrad_write)
        do kr=1,mrad
          lkr_write(kr)=.false.
        enddo
        do i=1,nrad_write
          lkr_write(kr_write(i))=.true.
        enddo
      else
        nrad_write=nrad
        do kr=1,nrad
          lkr_write(kr)=.true.
        enddo
      endif
      read(linput_3d,*) nruns
c
c  calculate nx and ny
c
      nx=0
      ny=0
      do ix=ix0,ix1,ixstep
        nx=nx+1
      enddo
      do iy=iy0,iy1,iystep
        ny=ny+1
      enddo
      allocate(idxy(nx,ny),stat=ierr)
      if(ierr.ne.0) write(*,*) 'allocation of idxy failed'
      do 995 irun=1,nruns
        bext=' '
        read(linput_3d,*) bext
c
C  OPEN GLOBAL FILES ALWAYS NEEDED
C
        CALL MOPEN(LOUT,'OUT',1,'MYNEW')
c        CALL MOPEN(LTIME,'TIME',1,'MYNEW')
        CALL MOPEN(LJOBLO,'JOBLOG',1,'MYNEW')
        CALL MOPEN(LDUMS,'DUMS',1,'MYNEW')
c        CALL MOPEN(LPHI,'PHI',0,'MYNEW')
        do lphi=1,maxlu
          if(.not.lu(lphi)) goto 5
        enddo
        call stop('main: no free unit number for lphi')
    5   continue
        lu(lphi)=.true.
C
C  READ INPUT PARAMETERS
C
        CALL RINPUT
        icrsw0=icrsw
        iconv0=iconv
c
c  if icont_opt=-1 out3d file is read and only columns where icont.eq.0 are
c  calculated
c  master reads the file in order to store old values and sets the
c  array iredo to 1 (redo column) or 0
c
        if(icont_opt.eq.-1) then
          allocate(iredo(nx,ny),stat=ierr)
          if(ierr.ne.0) write(*,*) 'allocation of iredo failed'
          if(master) then
            call mopen(lout3d,'out3d.'//bext(1:lent(bext)),0,'old')
            read(lout3d) itype,isize,cname
            read(lout3d) cdummy
            read(lout3d) itype,isize,cname
            read(lout3d) nx,ny,ndep
            allocate(taulg3d(nx,ny,ndep),stat=ierr)
            if(ierr.ne.0) write(*,*) 'allocation of taulg3d failed'
            allocate(cmasslg3d(nx,ny,ndep),stat=ierr)
            if(ierr.ne.0) write(*,*) 'allocation of cmasslg3d failed'
            allocate(xnorm3d(nx,ny,ndep),stat=ierr)
            if(ierr.ne.0) write(*,*) 'allocation of xnorm3d failed'
            allocate(x3d(nx,ny,ndep),stat=ierr)
            if(ierr.ne.0) write(*,*) 'allocation of x3d failed'
            allocate(height3d(nx,ny,ndep),stat=ierr)
            if(ierr.ne.0) write(*,*) 'allocation of height3d failed'
            read(lout3d) itype,isize,cname
            read(lout3d) dummy
            read(lout3d) itype,isize,cname
            read(lout3d) taulg3d
            read(lout3d) itype,isize,cname
            read(lout3d) cmasslg3d
            read(lout3d) itype,isize,cname
            nqtot=isize/nx/ny
            allocate(icont(nqtot,nx,ny),stat=ierr)
            if(ierr.ne.0) write(*,*) 'allocation of icont failed'
            read(lout3d) icont
            read(lout3d) itype,isize,cname
            read(lout3d) xnorm3d
            read(lout3d) itype,isize,cname
            read(lout3d) x3d
            read(lout3d) itype,isize,cname
            read(lout3d) height3d
            if(iwn.ne.0) then
              read(lout3d) itype,isize,cname
              nk=isize/nx/ny/ndep
              allocate(n3d(nx,ny,ndep,nk),stat=ierr)
              if(ierr.ne.0) write(*,*) 'allocation of n3d failed'
              read(lout3d) n3d
            endif
            call close(lout3d)
            do 20 iy=1,ny
              do 10 ix=1,nx
                if(icont(1,ix,iy).eq.0) then
                  iredo(ix,iy)=1
                else
                  iredo(ix,iy)=0
                endif
   10         continue
   20       continue
          endif
          call mpi_bcast(iredo,nx*ny,MPI_INTEGER,masterid,
     &     MPI_COMM_WORLD,mperr)
          call mpi_barrier(MPI_COMM_WORLD,mperr)
        endif
c
c  reformat atmos3d file for individual read from processors
c  first write the parts of atmos3d file that individual should 
c  treat into individual file, no interpolation
c  from dscale files
c
        call mopen(latmos,'ATMOS',0,'MYNEW')
        call mopen(latmos3d,'atmos3d.'//bext(1:lent(bext)),0,'old')
        read(latmos3d) itype,isize,cname
        write(*,*) itype,isize,cname
        read(latmos3d) nx2,ny2,ndep2
        write(*,*) nx2,ny2,ndep2
        allocate(temp3d(nx2,ny2,ndep2),stat=ierr)
        if(ierr.ne.0) write(*,*) 'allocation of temp3d failed'
        allocate(ne3d(nx2,ny2,ndep2),stat=ierr)
        if(ierr.ne.0) write(*,*) 'allocation of ne3d failed'
        allocate(vz3d(nx2,ny2,ndep2),stat=ierr)
        if(ierr.ne.0) write(*,*) 'allocation of vz3d failed'
        allocate(rho3d(nx2,ny2,ndep2),stat=ierr)
        if(ierr.ne.0) write(*,*) 'allocation of rho3d failed'
        allocate(widthx(nx2),stat=ierr)
        if(ierr.ne.0) write(*,*) 'allocation of widthx failed'
        allocate(widthy(ny2),stat=ierr)
        if(ierr.ne.0) write(*,*) 'allocation of widthy failed'
        allocate(height_in(ndep2),stat=ierr)
        if(ierr.ne.0) write(*,*) 'allocation of height_in failed'
        read(latmos3d) itype,isize,cname
        write(*,*) itype,isize,cname
        read(latmos3d) widthx
        read(latmos3d) itype,isize,cname
        write(*,*) itype,isize,cname
        read(latmos3d) widthy
        read(latmos3d) itype,isize,cname
        write(*,*) itype,isize,cname
        read(latmos3d) height_in
        read(latmos3d) itype,isize,cname
        write(*,*) itype,isize,cname
        read(latmos3d) ne3d
        read(latmos3d) itype,isize,cname
        write(*,*) itype,isize,cname
        read(latmos3d) temp3d
        read(latmos3d) itype,isize,cname
        write(*,*) itype,isize,cname
        read(latmos3d) vz3d        
        read(latmos3d) itype,isize,cname
        write(*,*) itype,isize,cname
        read(latmos3d) vz3d        
        read(latmos3d) itype,isize,cname
        write(*,*) itype,isize,cname
        read(latmos3d) vz3d        
        read(latmos3d) itype,isize,cname
        write(*,*) itype,isize,cname
        read(latmos3d) rho3d        
        call close(latmos3d)
        if(dscopt(1:6).eq.'dscal2' .or.
     &   (icont_opt.eq.-1 .and. master)) then
          allocate(dpnew3d(nx,ny,ndep2),stat=ierr)
          if(ierr.ne.0) write(*,*) 'allocation of dpnew3d failed'
          call mopen(ldsc3d,'dscal2.'//bext(1:lent(bext)),0,'old')
          read(ldsc3d) itype,isize,cname
          read(ldsc3d) dpnew3d
          call close(ldsc3d)
        endif
        if(istart.eq.-1) then
          call mopen(lrstrt,'RSTRT',0,'MYNEW')
          call mopen(ln3d,'n3d.'//bext(1:lent(bext)),0,'old')
          read(ln3d) itype,isize,cname
          write(*,*) itype,isize,cname
          nk2=isize/(nx2*ny2*ndep2)
          write(*,*) 'n3d: nk=',nk2
          allocate(n3d(nx2,ny2,ndep2,nk2),stat=ierr)
          if(ierr.ne.0) write(*,*) 'allocation of n3d failed'
          read(ln3d) n3d
          call close(ln3d)
        endif
c
c ixx,iyy is point in original atmosphere
c ix,iy is point in calculated atmosphere
c idxy is number of the processor that should handle that column
c
        id=0
        ix=0
        do 90 ixx=ix0,ix1,ixstep
          ix=ix+1
          iy=0
          do 80 iyy=iy0,iy1,iystep
            iy=iy+1
            if(icont_opt.eq.-1) then
              if(iredo(ix,iy).eq.0.) goto 80
            endif
c
c  write to ATMOS_x file if my column
c  write according to DSCALE option
c
            idxy(ix,iy)=id
            icxy(id)=icxy(id)+1
            if(id.eq.mpmyid) then
              write(*,'(a,5i5)') 'id,ix,iy,ixx,iyy ',id,ix,iy,ixx,iyy
              write(ljoblo,'(a,5i5)') 'id,ix,iy,ixx,iyy ',
     &         id,ix,iy,ixx,iyy
              if(dscopt(1:lent(dscopt)).eq.'dtrho') then
                call ipol_dscale(ndep2,height_in,temp3d(ixx,iyy,:),
     &           ne3d(ixx,iyy,:),vz3d(ixx,iyy,:),rho3d(ixx,iyy,:))
                write(latmos) ndep2
                if(height(2).gt.height(1)) then
                  write(latmos) -height(1:ndep2)*1.d-5
                else
                  write(latmos) height(1:ndep2)*1.d-5
                endif
                write(latmos) temp(1:ndep2)
                write(latmos) ne(1:ndep2)
                write(latmos) vel(1:ndep2)
                write(latmos) rho(1:ndep2)
              else if(dscopt(1:6).eq.'dscal2') then
                call ipol_dscal2(ndep2,height_in,temp3d(ixx,iyy,:),
     &           ne3d(ixx,iyy,:),vz3d(ixx,iyy,:),rho3d(ixx,iyy,:),
     &           dpnew3d(ixx,iyy,:))
                write(latmos) ndep2
                if(height(2).gt.height(1)) then
                  write(latmos) -height(1:ndep2)*1.d-5
                else
                  write(latmos) height(1:ndep2)*1.d-5
                endif
                write(latmos) temp(1:ndep2)
                write(latmos) ne(1:ndep2)
                write(latmos) vel(1:ndep2)
                write(latmos) rho(1:ndep2)
              else
                write(latmos) ndep2
                if(height_in(2).gt.height_in(1)) then
                  write(latmos) -height_in*1.d-5
                else
                  write(latmos) height_in*1.d-5
                endif
                write(latmos) real(temp3d(ixx,iyy,:),dp)
                write(latmos) real(ne3d(ixx,iyy,:),dp)
                write(latmos) real(vz3d(ixx,iyy,:),dp)
                write(latmos) real(rho3d(ixx,iyy,:),dp)
              endif
              if(istart.eq.-1) then
                call ipol_n(ndep2,nk2,height_in,n3d(ixx,iyy,:,:))
                write(lrstrt) n
              endif
            endif
            id=id+1
            if(id.gt.mpnproc-1) id=0
   80     continue
   90   continue
        call close(latmos)
        deallocate(temp3d,stat=ierr)
        if(ierr.ne.0) write(*,*) 'deallocation of temp3d failed'
        deallocate(ne3d,stat=ierr)
        if(ierr.ne.0) write(*,*) 'deallocation of ne3d failed'
        deallocate(vz3d,stat=ierr)
        if(ierr.ne.0) write(*,*) 'deallocation of vz3d failed'
        deallocate(rho3d,stat=ierr)
        if(ierr.ne.0) write(*,*) 'deallocation of rho3d failed'
        deallocate(widthx,stat=ierr)
        if(ierr.ne.0) write(*,*) 'deallocation of widthx failed'
        deallocate(widthy,stat=ierr)
        if(ierr.ne.0) write(*,*) 'deallocation of widthy failed'
        deallocate(height_in,stat=ierr)
        if(ierr.ne.0) write(*,*) 'deallocation of height_in failed'
        if(istart.eq.-1) then
          call close(lrstrt)
          deallocate(n3d,stat=ierr)
          if(ierr.ne.0) write(*,*) 'deallocation of n3d failed'
        endif
        if(dscopt(1:6).eq.'dscal2') then
          if(.not.master) deallocate(dpnew3d,stat=ierr)
          if(ierr.ne.0) write(*,*) 'deallocation of dpnew3d failed'
          if(master .and. icont_opt.ne.-1) deallocate(dpnew3d,stat=ierr)
          if(ierr.ne.0) write(*,*) 'deallocation of dpnew3d failed'
        endif
        ncxy=icxy(0)
        do i=1,mpnproc
          if(icxy(i).gt.ncxy) ncxy=icxy(i)
        enddo
        ncxy=ncxy+1
        itag_base=ncxy*mpnproc
        call mpi_barrier(MPI_COMM_WORLD,mperr)
        if(debug.and.master) write(*,*) 'icxy=',icxy(0:mpnproc-1)
c
c  go through all columns, calculate and write to file
c
        WRITE(CMPMYID,'(I4)') mpmyid
        K0=1
        CALL GETWRD(CMPMYID,K0,K1,K2)
        FILE2='ATMOS'//'_'//CMPMYID(K1:K2)
        call mopen(latmos,file2,0,'old')
        if(istart.eq.-1) then 
          FILE2='RSTRT'//'_'//CMPMYID(K1:K2)
          call mopen(lrstrt,file2,0,'old')
        endif
        call mopen(loutn,'OUTN',0,'MYNEW')
        id=0
        do 300 ix=1,nx
          do 200 iy=1,ny
            if(icont_opt.eq.-1) then
              if(iredo(ix,iy).eq.0) goto 200
            endif
            ext='_'//string3(ix)//'_'//string3(iy)
            if(idxy(ix,iy).eq.mpmyid) then
              ICALL=0
              CALL MCTIME('TOTAL ',0,0,5)
C
              iredo_switch=1
              iconv=iconv0
              CALL START
              IF(NK.GT.MK1) CALL STOP('MULTI: NK.GT.MK1')
              IF(NDEP.GT.MDEP1) CALL STOP('MULTI: NDEP.GT.MDEP1')
              NK1=NK
              NDEP1=NDEP
 97           continue
              CALL ITER(W,E,NK1,NDEP1)
c
c  if not converged and cr switching:
c    increase switching steps per decade by factor of two
c    this is tried three times
c  if converged, reset icrsw to INPUT values
c
              if(iconv.ne.1 .and. icrsw.gt.0 .and. 
     &         iredo_switch.le.iredo_max) then
                icrsw=icrsw*2
                iredo_switch=iredo_switch+1
                write(ljoblo,'(a,3i5)') 'retrying: id,ix,iy',
     &           mpmyid,ix,iy
                call initia
                goto 97
              endif
              icrsw=icrsw0
              if(iconv.ne.1 .and. itmax.gt.0) then
                tau(:)=1.d-10
                cmass(:)=1.d-10
                outint(:,:,:)=0.d0
                xnorm(:)=0.d0
                xsave(:)=0.d0
                dpnew(:)=height(:)
                n(:,:)=0.d0
                write(*,'(a,3i5)') 'diverged: id,ix,iy',
     &           mpmyid,ix,iy
                write(ljoblo,'(a,3i5)') 'diverged: id,ix,iy',
     &           mpmyid,ix,iy
              else
                write(*,'(a,3i5)') 'converged: id,ix,iy',id,ix,iy
                write(ljoblo,'(a,3i5)') 'converged: id,ix,iy',id,ix,iy
              endif
              CALL FORMAL
              call dscal2(dpnew)
c
c  write results to temporary file
c
              write(loutn) ix,iy
              write(loutn) tau
              write(loutn) cmass
              write(loutn) outint
              write(loutn) xnorm
              write(loutn) xsave
              write(loutn) height
              write(loutn) dpnew
              if(iwn.ne.0) write(loutn) n(1:nk,1:ndep)
C
c  clean up temporary files
c
              CALL MCTIME('TOTAL ',0,2,5)
              write(ljoblo,'(a)') 'NORMAL END'
  100         continue
              CALL MOPEN(LDUM,'DUMI',1,'UNKNOWN')
              CLOSE(LDUM,STATUS='DELETE')
              LU(LDUM)=.FALSE.
              CALL MOPEN(LDUM,'RSTRT2',1,'UNKNOWN')
              CLOSE(LDUM,STATUS='DELETE')
              LU(LDUM)=.FALSE.
              CALL REWIND(LOPC)
              CALL REWIND(LPHI)
            endif
            id=id+1
            if(id.gt.mpnproc-1) id=0
  200     continue
  300   continue
c
c  calculations ready, synchronize
c
        call close(loutn)
c 
c  calculate total number of frequencies for packing of icont
c
        nqtot=0
        do kr=1,nrad
          if(lkr_write(kr)) nqtot=nqtot+nq(kr)+1
        enddo
        qn=qnorm*1.d5/cc
c
c  allocate output variables for master if not already allocated
c
        if(icont_opt.ne.-1 .and. master) then
          allocate(icont(nqtot,nx,ny),stat=ierr)
          if(ierr.ne.0) write(*,*) 'allocation of icont failed'
          allocate(taulg3d(nx,ny,ndep),stat=ierr)
          if(ierr.ne.0) write(*,*) 'allocation of taulg3d failed'
          allocate(cmasslg3d(nx,ny,ndep),stat=ierr)
          if(ierr.ne.0) write(*,*) 'allocation of cmasslg3d failed'
          allocate(xnorm3d(nx,ny,ndep),stat=ierr)
          if(ierr.ne.0) write(*,*) 'allocation of xnorm3d failed'
          allocate(x3d(nx,ny,ndep),stat=ierr)
          if(ierr.ne.0) write(*,*) 'allocation of x3d failed'
          allocate(height3d(nx,ny,ndep),stat=ierr)
          if(ierr.ne.0) write(*,*) 'allocation of height3d failed'
          allocate(dpnew3d(nx,ny,ndep),stat=ierr)
          if(ierr.ne.0) write(*,*) 'allocation of dpnew3d failed'
          if(iwn.ne.0) allocate(n3d(nx,ny,ndep,nk),stat=ierr)
          if(ierr.ne.0) write(*,*) 'allocation of n3d failed'
        endif
c
c  allocate buffer for transfer from slaves to master
c
        allocate(buf(ndep))
        allocate(buf2(nqtot))
        call mpi_barrier(MPI_COMM_WORLD,mperr)
        if(debug .and. master) then
          write(*,*) 'Master reformats'
          write(*,*) 'nrad,nqtot=',nrad,nqtot
        endif
c
c  master reformats output from OUTN files
c
        call mopen(loutn,'OUTN',0,'MYOLD')
        icxy(0:mpnproc)=0
  350   continue
          read(loutn,end=400) ix,iy
          read(loutn) tau
          read(loutn) cmass
          read(loutn) outint
          read(loutn) xnorm
          read(loutn) xsave
          read(loutn) height
          read(loutn) dpnew
          if(iwn.ne.0) read(loutn) n(1:nk,1:ndep)
c
c  send results to master
c
          nsend=8
          if(iwn.ne.0) nsend=nsend+nk
          if(master) then  ! master just stores data without communication
            taulg3d(ix,iy,:)=log10(tau(1:ndep))
            cmasslg3d(ix,iy,:)=log10(cmass(1:ndep))
            xnorm3d(ix,iy,:)=xnorm(1:ndep)
            x3d(ix,iy,:)=xnorm(1:ndep)*xsave(1:ndep)
            height3d(ix,iy,:)=height(1:ndep)
            dpnew3d(ix,iy,:)=dpnew(1:ndep)
            nu1=1
            do kr=1,nrad
              if(lkr_write(kr)) then
                icont(nu1:nu1+nq(kr),ix,iy)=outint(0:nq(kr),nmu,kr)
                nu1=nu1+nq(kr)+1
              endif
            enddo
            if(iwn.ne.0) then
              do i=1,nk
                n3d(ix,iy,:,i)=n(i,1:ndep)
              enddo
            endif
          else
            icxy(mpmyid)=icxy(mpmyid)+1    ! increment column counter for tag purposes
            itag=mpmyid*ncxy+icxy(mpmyid)
            if(debug) write(*,*) 'sending (from, to, tag)',
     &       mpmyid,masterid,itag
            ibuf(1)=ix
            ibuf(2)=iy
            call mpi_ssend(ibuf,2,MPI_INTEGER,
     &       masterid,itag,MPI_COMM_WORLD,mperr)
            itag=itag_base+((iy-1)*nx+(ix-1))*nsend
            buf=log10(tau(1:ndep))
            if(debug) write(*,*) 'sending (from, to, tag)',
     &       mpmyid,masterid,itag+2
            call mpi_send(buf,ndep,MPI_REAL,
     &       masterid,itag+2,MPI_COMM_WORLD,mperr)
            buf=log10(cmass(1:ndep))
            if(debug) write(*,*) 'sending (from, to, tag)',
     &       mpmyid,masterid,itag+3
            call mpi_send(buf,ndep,MPI_REAL,
     &       masterid,itag+3,MPI_COMM_WORLD,mperr)
            buf=xnorm(1:ndep)
            if(debug) write(*,*) 'sending (from, to, tag)',
     &       mpmyid,masterid,itag+4
            call mpi_send(buf,ndep,MPI_REAL,
     &       masterid,itag+4,MPI_COMM_WORLD,mperr)
            buf=xnorm(1:ndep)*xsave(1:ndep)
            if(debug) write(*,*) 'sending (from, to, tag)',
     &       mpmyid,masterid,itag+5
            call mpi_send(buf,ndep,MPI_REAL,
     &       masterid,itag+5,MPI_COMM_WORLD,mperr)
            buf=height(1:ndep)
            if(debug) write(*,*) 'sending (from, to, tag)',
     &       mpmyid,masterid,itag+6
            call mpi_send(buf,ndep,MPI_REAL,
     &       masterid,itag+6,MPI_COMM_WORLD,mperr)
            buf=dpnew(1:ndep)
            if(debug) write(*,*) 'sending (from, to, tag)',
     &       mpmyid,masterid,itag+7
            call mpi_send(buf,ndep,MPI_REAL,
     &       masterid,itag+7,MPI_COMM_WORLD,mperr)
            nu1=1
            do kr=1,nrad
              if(lkr_write(kr)) then
                buf2(nu1:nu1+nq(kr))=outint(0:nq(kr),nmu,kr)
                nu1=nu1+nq(kr)+1
              endif
            enddo
            if(debug) write(*,*) 'sending (from, to, tag)',
     &       mpmyid,masterid,itag+8
            call mpi_send(buf2,nqtot,MPI_REAL,
     &       masterid,itag+8,MPI_COMM_WORLD,mperr)
            if(iwn.ne.0) then
              do i=1,nk
                buf=n(i,1:ndep)
                if(debug) write(*,*) 'sending (from, to, tag)',
     &           mpmyid,masterid,itag+8+i
                call mpi_send(buf,ndep,MPI_REAL,masterid,
     &          itag+8+i,MPI_COMM_WORLD,mperr)
              enddo
            endif
          endif ! master

! master loops through processors and receives data

  390     continue
          if(master) then
            icxy(masterid)=icxy(masterid)+1
            do id=0,mpnproc-1
              if(id.ne.masterid .and. more_to_receive(id)) then
                itag=id*ncxy+icxy(masterid)
                if(debug) write(*,*) 'receiving (from, to, tag)',
     &           id,masterid,itag
                call mpi_recv(ibuf,2,MPI_INTEGER,
     &           id,itag,MPI_COMM_WORLD,status,mperr)
                ix=ibuf(1)
                iy=ibuf(2)
                if(ix.lt.0) then  ! slave has reached EOF
                  more_to_receive(id)=.false.
                else
                  itag=itag_base+((iy-1)*nx+(ix-1))*nsend
                  if(debug) write(*,*) 'receiving (from, to, tag)',
     &             id,masterid,itag+2
                  call mpi_recv(buf,ndep,MPI_REAL,
     &             id,itag+2,MPI_COMM_WORLD,status,mperr)
                  taulg3d(ix,iy,:)=buf
                  if(debug) write(*,*) 'receiving (from, to, tag)',
     &             id,masterid,itag+3
                  call mpi_recv(buf,ndep,MPI_REAL,
     &             id,itag+3,MPI_COMM_WORLD,status,mperr)
                  cmasslg3d(ix,iy,:)=buf
                  if(debug) write(*,*) 'receiving (from, to, tag)',
     &             id,masterid,itag+4
                  call mpi_recv(buf,ndep,MPI_REAL,
     &             id,itag+4,MPI_COMM_WORLD,status,mperr)
                  xnorm3d(ix,iy,:)=buf
                  if(debug) write(*,*) 'receiving (from, to, tag)',
     &             id,masterid,itag+5
                  call mpi_recv(buf,ndep,MPI_REAL,
     &             id,itag+5,MPI_COMM_WORLD,status,mperr)
                  x3d(ix,iy,:)=buf
                  if(debug) write(*,*) 'receiving (from, to, tag)',
     &             id,masterid,itag+6
                  call mpi_recv(buf,ndep,MPI_REAL,
     &             id,itag+6,MPI_COMM_WORLD,status,mperr)
                  height3d(ix,iy,:)=buf
                  if(debug) write(*,*) 'receiving (from, to, tag)',
     &             id,masterid,itag+7
                  call mpi_recv(buf,ndep,MPI_REAL,
     &             id,itag+7,MPI_COMM_WORLD,status,mperr)
                  dpnew3d(ix,iy,:)=buf
                  if(debug) write(*,*) 'receiving (from, to, tag)',
     &             id,masterid,itag+8,nqtot
                  call mpi_recv(buf2,nqtot,MPI_REAL,
     &             id,itag+8,MPI_COMM_WORLD,status,mperr)
                  icont(1:nqtot,ix,iy)=buf2
                  if(iwn.ne.0) then
                    do i=1,nk
                      if(debug) write(*,*) 'receiving (from, to, tag)',
     &                 id,masterid,itag+8+i
                      call mpi_recv(buf,ndep,MPI_REAL,id,
     &               itag+8+i,MPI_COMM_WORLD,status,mperr)
                      n3d(ix,iy,:,i)=buf
                    enddo  ! nk-loop
                  endif   ! iwn
                endif   ! ix ge 0
              endif   ! id.eq.master .and. more_to_receive 
            enddo   ! mpnproc-loop
          endif   ! master
        goto 350 ! go to next read
 400    continue  ! EOF
        if(slave) then   ! EOF, send signal to master
          icxy(mpmyid)=icxy(mpmyid)+1    ! increment column counter for tag purposes
          itag=mpmyid*ncxy+icxy(mpmyid)
          if(debug) write(*,*) 'sending (from, to, tag)',
     &     mpmyid,masterid,itag
          ibuf(1)=-1     ! send negative ix to signal EOF to master
          ibuf(2)=-1 
          call mpi_ssend(ibuf,2,MPI_INTEGER,
     &     masterid,itag,MPI_COMM_WORLD,mperr)
        else  ! if master, check to see that everything has been received
          do id=0,mpnproc-1
            if(id.ne.masterid .and. more_to_receive(id)) goto 390
          enddo
        endif
c 
c  all send/receive finished, synchronize
c  master writes out3d file
c
        call mpi_barrier(MPI_COMM_WORLD,mperr)
        if(master) then
c
c  fill xl with wavelengths in AA in vacuum
c
          allocate(xl(nqtot),stat=ierr)
          if(ierr.ne.0) write(*,*) 'allocation of xl failed'
          nu0=1
          kt=1
          do kr=1,nline
            if(lkr_write(kr)) then
              do nu=1,nq(kr)
                xl(nu0+nu)=alamb(kr)-alamb(kr)*q(nu,kr)*qn/
     *           (1.0d0+q(nu,kr)*qn)
              enddo
              nu0=nu0+nq(kr)+1
              if(iwide(kr)) kt=kt+1
            endif
          enddo
          do kr=nline+1,nrad
            if(lkr_write(kr)) then
              do nu=1,nq(kr)
                xl(nu0+nu)=cc*1.e8/frq(nu,kt)
              enddo
              nu0=nu0+nq(kr)+1
              kt=kt+1
            endif
          enddo
          i=1
          do kr=1,nrad
            if(lkr_write(kr)) then
              nq_write(i)=nq(kr)
              i=i+1
            endif
          enddo
          call mopen(lout3d,'out3d.'//bext(1:lent(bext)),0,'MASTERNEW')
          itype=7
          isize=80
          cname='id      '
          write(lout3d) itype,isize,cname
          write(*,*) itype,isize,cname
          write(lout3d) bext
          itype=3
          isize=5+nrad_write
          cname='dim     '
          write(lout3d) itype,isize,cname
          write(*,*) itype,isize,cname
          write(lout3d) nx,ny,ndep,mq,nrad_write,nq_write(1:nrad_write)
          itype=5
          isize=nqtot
          cname='xl      '
          write(lout3d) itype,isize,cname
          write(*,*) itype,isize,cname
          write(lout3d) xl(1:nqtot)
          itype=4
          isize=nx*ny*ndep
          cname='taulg3d '
          write(lout3d) itype,isize,cname
          write(*,*) itype,isize,cname
          write(lout3d) taulg3d
          itype=4
          isize=nx*ny*ndep
          cname='cmass3d '
          write(lout3d) itype,isize,cname
          write(*,*) itype,isize,cname
          write(lout3d) cmasslg3d
          itype=4
          isize=nx*ny*nqtot
          cname='icont   '
          write(lout3d) itype,isize,cname
          write(*,*) itype,isize,cname
          write(lout3d) icont
          itype=4
          isize=nx*ny*ndep
          cname='xnorm3d '
          write(lout3d) itype,isize,cname
          write(*,*) itype,isize,cname
          write(lout3d) xnorm3d
          itype=4
          isize=nx*ny*ndep
          cname='x3d     '
          write(lout3d) itype,isize,cname
          write(*,*) itype,isize,cname
          write(lout3d) x3d
          itype=4
          isize=nx*ny*ndep
          cname='height3d'
          write(lout3d) itype,isize,cname
          write(*,*) itype,isize,cname
          write(lout3d) height3d
          if(iwn.ne.0) then
            itype=4
            isize=nx*ny*ndep*nk
            cname='n3d     '
            write(lout3d) itype,isize,cname
            write(*,*) itype,isize,cname
            write(lout3d) n3d
          endif
          call close(lout3d)
          call mopen(ldsc3d,'dscal2.'//bext(1:lent(bext)),0,'MASTERNEW')
          itype=4
          isize=nx*ny*ndep
          cname='dpnew3d '
          write(ldsc3d) itype,isize,cname
          write(*,*) itype,isize,cname
          write(ldsc3d) dpnew3d
          call close(ldsc3d)
          deallocate(taulg3d,stat=ierr)
          if(ierr.ne.0) write(*,*) 'deallocation of taulg3d failed'
          deallocate(cmasslg3d,stat=ierr)
          if(ierr.ne.0) write(*,*) 'deallocation of cmasslg3d failed'
          deallocate(xnorm3d,stat=ierr)
          if(ierr.ne.0) write(*,*) 'deallocation of xnorm3d failed'
          deallocate(x3d,stat=ierr)
          if(ierr.ne.0) write(*,*) 'deallocation of x3d failed'
          deallocate(height3d,stat=ierr)
          if(ierr.ne.0) write(*,*) 'deallocation of height3d failed'
          deallocate(icont,stat=ierr)
          if(ierr.ne.0) write(*,*) 'deallocation of icont failed'
          deallocate(dpnew3d,stat=ierr)
          if(ierr.ne.0) write(*,*) 'deallocation of depnew3d failed'
          deallocate(xl,stat=ierr)
          if(ierr.ne.0) write(*,*) 'deallocation of xl failed'
          if(iwn.ne.0) deallocate(n3d,stat=ierr)
          if(ierr.ne.0) write(*,*) 'deallocation of n3d failed'
        endif
        close(loutn,status='delete')
        lu(loutn)=.false.
        deallocate(buf)
        deallocate(buf2)
c
c  close files and delete scratch files
c
        lopnmx=1
        call mopen(ldumi,'DUMI',1,'MYOLD')
        DO 990 I=1,MAXLU
          IF(LU(I)) THEN
            lopnmx=lopnmx+1
            IF(I.EQ.LDUMS .OR. I.EQ.LOUT .OR.
     *       i.eq.latmos .or. i.eq.LJOBLO .OR.
     *       I.EQ.LDUMI) then
              CLOSE(I,STATUS='DELETE')
              lu(i)=.false.
            ELSE if(i.ne.linput_3d) then
              CLOSE(I)
              lu(i)=.false.
            ENDIF
          ENDIF
  990   CONTINUE
        call mpi_barrier(MPI_COMM_WORLD,mperr)
  995 continue
      write(*,'(a,i4,a,i5)') 'id=',mpmyid,' max open files:',lopnmx
      call close(linput_3d)
      call mpi_finalize(mperr)
      STOP
      END
C
C************************************************************************
C
      character*3 function string3(i)
c
      if(i.lt.10) then 
        write(string3,'(A2,i1)') '00',i
      else if(i.lt.100) then
        write(string3,'(A1,i2)') '0',i
      else 
        write(string3,'(i3)') i
      endif

      return
      end
c**********************************************************************
c
      function lent(text)
c
c  finds the index of the last non-blank character of a text-string
c
      character*(*) text
c
      l=len(text)
      do 10 i=l,1,-1
        if(text(i:i).ne.' ') goto 100
   10 continue
  100 continue
      lent=i
c
      end
c
c**********************************************************************
c
      subroutine ipol_dscale(ndep2,height_in,temp_in,ne_in,vz_in,rho_in)
c
c  finds a new more optimal depth-scale and interpolates input atmosphere to it
c  use multi variables for output
c
      include 'PREC'
      include 'PARAM'
      include 'CATMOS'
      include 'CATMO2'
      include 'CCONST'
c
      integer ndep2
      real(4) height_in(ndep2),temp_in(ndep2),ne_in(ndep2),vz_in(ndep2)
      real(4) rho_in(ndep2)
      dimension aind2(mdep),f(mdep),x(mdep),thmbf(mdep),xhmbf(mdep)
      logical debug
c
      debug=.false.
      ttop=5.d4
      tbottom=1.2d4 ! not used now
      taumax=100
      grph=2.26e-24
c
c  calculate tau scale from H-bf opacity alone
c  assume all hydrogen to be neutral
c
      thmbf(1)=1.e-6
      crhmbf=2.9256e-17
      do k=1,ndep2
        xhmbf(k)=1.03526d-16*ne_in(k)*crhmbf/temp_in(k)**1.5*
     &  exp(0.754d0*ee/bk/temp_in(k))*rho_in(k)/grph
      enddo
      do k=2,ndep2
        thmbf(k)=thmbf(k-1)+0.5*(xhmbf(k)+xhmbf(k-1))*
     &   abs(height_in(k-1)-height_in(k))
      enddo
c
c  set atmospheric range from temperature of ttop and tau(hmbf)=taumax
c
      k0=1
      k1=1
      aind2(:)=0.0d0
      do k=2,ndep2
        if(temp_in(k).gt.ttop .and. k0.eq.k-1) k0=k
c        if(temp_in(k).lt.tbottom) k1=k  ! old criterion
        if(thmbf(k).lt.taumax) k1=k
      enddo
c
      do k=k0+1,k1
        tdiv=abs(log10(temp_in(k))-log10(temp_in(k-1)))/log10(1.1)
        rdiv=abs(log10(rho_in(k))-log10(rho_in(k-1)))/log10(1.1)
        taudiv=(log10(thmbf(k))-log10(thmbf(k-1)))/0.1
c        write(*,'(i4,3f10.4)') k,tdiv,rdiv,taudiv
        aind2(k)=aind2(k-1)+max(max(tdiv,rdiv),taudiv)
      enddo
c
      aind2=aind2/aind2(k1)*(ndep2-1)
c
c  make new height scale and interpolate atmosphere
c
      nn=k1-k0+1
      f(1:nn)=height_in(k0:k1)
      do k=1,ndep2
        xp=(k-1)
        call intep(k,aind2(k0:k1),f,nn,iup,xp,height(k))
c        if(k.gt.1) then
c          if(height(k).gt.height(k-1)) debug=.true.
c        endif
      enddo
      f(1:nn)=alog(temp_in(k0:k1))
      x(1:nn)=height_in(k0:k1)
      do 400 k=1,ndep2
        xp=height(k)
        call intep(k,x,f,nn,iup,xp,temp(k))
        temp(k)=exp(temp(k))
        if(temp(k).lt.1.d3) then
          do kk=2,nn
            if((xp-x(kk-1))*
     &         (xp-x(kk)) .le. 0.d0) goto 350
          enddo
  350     continue
          temp(k)=f(kk-1)+(xp-x(kk-1))/
     &     (x(kk)-x(kk-1))*(f(kk)-f(kk-1))
          temp(k)=exp(temp(k))
          write(*,*) 'linear interpolation, temp(k)=',temp(k)
        endif  
  400 continue
      f(1:nn)=alog(ne_in(k0:k1))
      do 500 k=1,ndep2
        xp=height(k)
        call intep(k,x,f,nn,iup,xp,ne(k))
        ne(k)=exp(ne(k))
  500 continue
      f(1:nn)=vz_in(k0:k1)
      do 600 k=1,ndep2
        xp=height(k)
        call intep(k,x,f,nn,iup,xp,vel(k))
  600 continue
      f(1:nn)=alog(rho_in(k0:k1))
      do 700 k=1,ndep2
        xp=height(k)
        call intep(k,x,f,nn,iup,xp,rho(k))
        rho(k)=exp(rho(k))
  700 continue
c
      if(debug) then
        write(*,'(a)') 'input variables'
        write(*,'(a)') '  k       height        temp          ne'//
     &   '          vz         rho'
        do 750 k=1,ndep2
          write(*,'(i3,1p,5e12.4)') k,height_in(k),temp_in(k),ne_in(k),
     &     vz_in(k),rho_in(k)
  750   continue
        write(*,*) ' '
        write(*,'(a)') 'interpolation array'
        write(*,'(a)') '  k       aind2'
        do 760 k=k0,k1
          write(*,'(i3,1p,e12.4)') k,aind2(k)
  760   continue
        write(*,*) ' '
        write(*,'(a)') 'interpolated variables'
        write(*,'(a)') '  k       height        temp          ne'//
     &   '          vz         rho'
        do 800 k=1,ndep2
          write(*,'(i3,1p,5e12.4)') k,height(k),temp(k),ne(k),
     &     vel(k),rho(k)
  800   continue
      endif
c
      end
c
c*****************************************************************
c
      subroutine ipol_dscal2(ndep2,height_in,temp_in,ne_in,vz_in,rho_in,
     & height2)
c
c  interpolate to height2 scale
c  use multi variables for output
c
      include 'PREC'
      include 'PARAM'
      include 'CATMOS'
      include 'CATMO2'
c
      integer ndep2
      real(4) height_in(ndep2),temp_in(ndep2),ne_in(ndep2),vz_in(ndep2)
      real(4) rho_in(ndep2),height2(ndep2)
      dimension f(mdep),x(mdep)
      logical debug
c
      debug=.false.
c
c  make new height scale and interpolate atmosphere
c  height_in and height2 in cm, converted to km when written to ATMOS_ files
c
      height(1:ndep2)=height2(1:ndep2)*1.e5
      f(1:ndep2)=alog(temp_in(1:ndep2))
      x(1:ndep2)=height_in(1:ndep2)
      do 400 k=1,ndep2
        xp=height(k)
        call intep(k,x,f,ndep2,iup,xp,temp(k))
        temp(k)=exp(temp(k))
  400 continue
      f(1:ndep2)=alog(ne_in(1:ndep2))
      do 500 k=1,ndep2
        xp=height(k)
        call intep(k,x,f,ndep2,iup,xp,ne(k))
        ne(k)=exp(ne(k))
  500 continue
      f(1:ndep2)=vz_in(1:ndep2)
      do 600 k=1,ndep2
        xp=height(k)
        call intep(k,x,f,ndep2,iup,xp,vel(k))
  600 continue
      f(1:ndep2)=alog(rho_in(1:ndep2))
      do 700 k=1,ndep2
        xp=height(k)
        call intep(k,x,f,ndep2,iup,xp,rho(k))
        rho(k)=exp(rho(k))
  700 continue
c
      if(debug) then
        write(*,'(a)') 'input variables'
        write(*,'(a)') '  k       height        temp          ne'//
     &   '          vz         rho'
        do 750 k=1,ndep2
          write(*,'(i3,1p,5e12.4)') k,height_in(k),temp_in(k),ne_in(k),
     &     vz_in(k),rho_in(k)
  750   continue
        write(*,*) ' '
        write(*,'(a)') 'interpolated variables'
        write(*,'(a)') '  k       height        temp          ne'//
     &   '          vz         rho'
        do 800 k=1,ndep2
          write(*,'(i3,1p,5e12.4)') k,height(k),temp(k),ne(k),
     &     vel(k),rho(k)
  800   continue
      endif
c
      end
c
c*****************************************************************
c
      subroutine ipol_n(ndep2,nk2,height_in,n_in)
c
c  interpolate n from height_in to height scale
c  use multi variables for output
c
      include 'PREC'
      include 'PARAM'
      include 'CATOM'
      include 'CATMOS'
      include 'CATMO2'
c
      integer ndep2,nk2
      real(4) height_in(ndep2),n_in(ndep2,nk2)
      dimension y(mdep),x(mdep)
      logical debug
c
      debug=.false.
c
c  make new height scale and interpolate atmosphere
c
      x(1:ndep2)=height_in(1:ndep2)
      do 500 i=1,nk2
        y(1:ndep2)=alog(n_in(1:ndep2,i))
        do 400 k=1,ndep2
          xp=height(k)
          call intep(k,x,y,ndep2,iup,xp,n(i,k))
          n(i,k)=exp(n(i,k))
  400   continue
 500  continue
c
      if(debug) then
        write(*,'(a)') 'input variables'
        write(*,'(a)') '  k       height           n'
        do 750 k=1,ndep2
          write(*,'(i3,1p,7e10.2)') k,height_in(k),n_in(k,:)
  750   continue
        write(*,*) ' '
        write(*,'(a)') 'interpolated variables'
        write(*,'(a)') '  k       height           n'
        do 800 k=1,ndep2
          write(*,'(i3,1p,7e10.2)') k,height(k),n(1:nk2,k)
  800   continue
      endif
c
      end
c
c*****************************************************************
c
      subroutine intep(it,x,f,n,iup,xp,p)
c
c  interpolation routine based on hermite polynomials
c  ref: publications of the dominion astrophysical observatory, xvi,6,67
c       graham hill: intep, an effective interpolation subroutine
c
      include 'PREC'
      integer it,n,iup
      real(8) lp1,lp2,l1,l2
      dimension f(n),x(n)
      save
c
      if(it.eq.1) then
        ier=1
        io=1
        iup=0
        if(x(2).lt.x(1)) iup=1
        n1=n-1
        if((xp.ge.x(n).and.iup.eq.0).or.(xp.le.x(n).and.iup.eq.1)) then
          p=f(n)
          ier=2
          return
        else if((xp.le.x(1).and.iup.eq.0).or.
     &    (xp.ge.x(1).and.iup.eq.1)) then
          p=f(1)
          ier=2
          return
        endif
      endif
    8 do 1 i=io,n
        if(xp.lt.x(i).and.iup.eq.0) goto 2
        if(xp.gt.x(i).and.iup.eq.1) goto 2
    1 continue
      p=f(n)
      ier=2
      return
    2 i=i-1
      if(i.eq.io-1) goto 4
      io=i+1
      lp1=1.d0/(x(i)-x(i+1))
      lp2=1.d0/(x(i+1)-x(i))
      if(i.eq.1) fp1=(f(2)-f(1))/(x(2)-x(1))
      if(i.eq.1) goto 3
      fp1=(f(i+1)-f(i-1))/(x(i+1)-x(i-1))
    3 if(i.ge.n1) fp2=(f(n)-f(n-1))/(x(n)-x(n-1))
      if(i.ge.n1) goto 4
      fp2=(f(i+2)-f(i))/(x(i+2)-x(i))
    4 xpi1=xp-x(i+1)
      xpi=xp-x(i)
      l1=xpi1*lp1
      l2=xpi*lp2
      p=f(i)*(1.d0-2.d0*lp1*xpi)*l1*l1+f(i+1)*(1.d0-2.d0*lp2*xpi1)*
     & l2*l2+fp2*xpi1*l2*l2+fp1*xpi*l1*l1
c
      return
      end   
   
