Return to inv2d.f CVS log

Up to [Development] / JSOC / proj / globalhs / apps

now write output splittings to file splittings.out instead of stdout

      program inv2d
c 070701 Changed output format to accommodate l=1000

      implicit double precision (a-h,o-z)
      include 'include1.f'
      include 'include2.f'

      parameter (maxpoints=maxr1*maxt1,nblock1=10)

      double precision split(maxnlm),spliti(maxl)
      double precision sigma(maxnlm),sigmai(maxl)
      double precision sigma1(maxnlm),ss1(maxnlm),ss2(maxnlm)
      double precision rmesh1(maxr1),tmesh1(maxt1)
      double precision rweights(maxr1),tweights(maxt1)
      double precision fr(maxr1),ft(maxr1),fwr(maxr1),fwt(maxt1)
      character*80 cefunc,outfile,splitfile
      real dummy(2),dtime,etime,time1
      real mur,mut
      real f1(maxmodes,maxr1),f2(maxmodes,maxr1)
      real g1(maxlm,maxt1),g2(maxlm,maxt1)
      double precision f1x(maxmodes),f2x(maxmodes)
      double precision g1x(maxlm),g2x(maxlm)
      double precision f1i(maxr1),f2i(maxr1),g1i(maxt1),g2i(maxt1)
      double precision f1o(maxt1),f2o(maxt1)
      double precision a(0:maxs1-1,maxpoints)
      double precision ax(0:maxs1-1),bx(nblock,nblock)
      integer iax(maxpoints)
      double precision b(maxpoints,maxpoints),right(maxpoints)
      double precision c1(maxr1,maxr1),c2(maxt1,maxt1)
      double precision d(maxpoints,maxpoints),help(3*maxpoints)
      integer iwork(maxpoints)
c     double precision d1(maxpoints,maxpoints),cov(maxpoints,maxpoints)
      double precision var(maxpoints),dib(maxpoints)
      double precision dibx(maxpoints,nblock1)
c     double precision di(maxpoints)
      double precision coeff(maxnlm)
      double precision omega(maxr1,maxt1)
      double precision fro(maxr1,maxt1),fto(maxr1,maxt1)
      integer jav(maxpoints),kav(maxpoints)
      real avkern(maxtav,maxrav)
      integer nnlm1(maxnlm),llist(maxnlm)
      real fnlm(maxnlm)
      logical found
      integer ijk(maxr1,maxt1)
      logical ibig(maxpoints)
c     integer lsplit(maxnlm)
c     integer msplit(maxnlm),isplit(maxnlm),lmsplit(maxnlm)
      integer mindex(maxnlm)
      real x1(maxtav),x2(maxrav),yy2(maxtav,maxrav),mjta(4)
      real xr(10)
      integer ixl(maxmodes),ixh(maxmodes)
      

c     nmode,lmode: identification of modes
c     freq: frequency of mode (muHz)
c     radmesh: mesh in radius for kernels and eigenfunctions
c     nexp: number of expansion coefficients
c     cefunc: control file for reading of eigenfunctions
c     iomega: 0 to read omega from input
c             1 to read omega from file
c     omegafile: file from which omega is read
c     nrin, ntin: number of points where omega is given
c     rin,tin,omegain: points where omega is given and the values there
c                  omega is assumed constant outside the given interval

      open (20,file='status.log',form='formatted')
      write (20,*) 1
      close (20)
c     read parameters etc.
      call clnfile(5,10)
      write (6,*) 'ierr, err'
      read (10,*) ierr,err
      write (6,*) 'name of file to read splittings from'
      read (10,'(A80)') splitfile
      write (6,*) 'control file for reading of eigenfunctions'
      read (10,'(A80)') cefunc
      write (6,*) 'number of points to skip in r for setting up f'
      read (10,*) nskipr
      write (6,*) 'number of points in x for setting up the plm''s'
      read (10,*) nsetp
      write (6,*) 'number of points to skip in x for setting up g'
      read (10,*) nskipt
      write (6,*) 'irdf, iwrf, irdg, iwrg, irdb, iwrb'
      read (10,*) irdf,iwrf,irdg,iwrg,irdb,iwrb
      write (6,*) '0 for multiple centerpoints, 1 for 1'
      read (10,*) icenter
      write (6,*) 'mur, mut'
      read (10,*) mur,mut
      nreg=2
      write (6,*) 'icovar, icoeff, iavkern, nsr, nst, nsr2, nt2'
      read (10,*) icovar, icoeff, iavkern, nsr, nst, nsr2, nt2
      if (nsr.eq.0) then
        nav=nst
        do 105,i=1,nav
          read (10,*) jav(i),kav(i)
 105    continue
      end if
c     write (6,*) 'name of output file'
c     read (10,'(A80)') outfile
      close(10)

      pi=4.0D0*atan(1.0D0)
      pi2=pi/2.0D0
      write (6,*) 'start reading splittings'
      time1=dtime(dummy)
      open (11,file=splitfile,status='old')
      i=1
 110  continue
      if (ierr.eq.0) then
        read (11,*,end=150) lnlm(i),nnlm1(i),fnlm(i),
     c                      mnlm(i),itnlm(i),it1nlm(i),split(i)
        sigma(i)=err
      else
        read (11,*,end=150) lnlm(i),nnlm1(i),fnlm(i),
     c                      mnlm(i),itnlm(i),it1nlm(i),
     c                      split(i),sigma(i)
      end if
      sigma1(i)=1.0/sigma(i)
      ss1(i)=split(i)*sigma1(i)
      ss2(i)=split(i)*sigma1(i)**2
      i=i+1
      goto 110
 150  continue
      close(11)
      write (6,*) 'splittings have been read, time=',
     c            dtime(dummy),dummy(1),dummy(2)
      nnlm=i-1
      call setilm
      write (6,*) dtime(dummy),dummy(1),dummy(2)
      lmax1=-1
      do 10,i=1,nnlm
        if (lnlm(i).gt.lmax1) lmax1=lnlm(i)
 10   continue
      i1=0
      do 100,l=0,lmax1
        nl=0
        do 20,i=1,nnlm
          if (lnlm(i).eq.l) then
            nl=nl+1
            llist(nl)=i
          end if
 20     continue
        i2=i1+1
        do 60,i=1,nl
          il=llist(i)
          found=.false.
          do 30,ii=i2,i1
            if (nnlm1(il).eq.nmode(ii)) then
              found=.true.
              ifound=ii
            end if
  30      continue
          if (found) then
            iln(il)=ifound
          else
            i1=i1+1
            lmode(i1)=l
            nmode(i1)=nnlm1(il)
            freq(i1)=fnlm(il)
            iln(il)=i1
          end if
 60     continue
 100  continue
      nmodes=i1
      write (6,*) dtime(dummy),dummy(1),dummy(2)
      write (6,*) '#nl=',nmodes
      write (6,*) '#nlm=',nnlm
c     open (11,file=splitfile,form='unformatted',status='old')
c     i=1
c     i1=0
c110  read (11,end=150) lmode(i),nmode(i),freq(i),
c    c                  (spliti(m),m=1,lmode(i))
c     m0index(i)=i1
c     do 130,m=1,lmode(i)
c       i1=i1+1
c       split(i1)=spliti(m)
c       sigma(i1)=1.0E0
c       sigma1(i1)=1.0D0/sigma(i1)
c       ss1(i1)=split(i1)*sigma1(i1)
c       lsplit(i1)=lmode(i)
c       msplit(i1)=m
c       isplit(i1)=i
c130  continue
c     i=i+1
c     goto 110
c150  nmodes=i-1
c     nnlm=i1
c     close (11)
      write (6,*) 'start reading eigenfunctions'
      time1=dtime(dummy)
      call readefun(cefunc)
c     do i=1,nrad
c       write (6,*) dradmesh(i)
c     end do
      write (6,*) 'eigenfunctions have been read, time=',
     c            dtime(dummy),dummy(1),dummy(2)
      if (mod(nrad-1,nskipr).ne.0) then
        write (6,*) 'skip to get inversion mesh in radius does not'
        write (6,*) 'divide into number of kernel meshpoints-1'
        stop
      end if
      if (mod(nrad-1,nsr2).ne.0) then
        write (6,*) 'skip to get av. kernel mesh in radius does not'
        write (6,*) 'divide into number of kernel meshpoints-1'
        stop
      end if
      nrad1=(nrad-1)/nskipr+1
      nr2=(nrad-1)/nsr2+1
      do 160,i=1,nrad1
        rmesh1(i)=radmesh(nskipr*(i-1)+1)
 160  continue
      call dsetint(1,nrad1,rmesh1,rweights)
      if (irdf.eq.0) then
        write (6,*) 'start setting up f matrix'
        time1=dtime(dummy)
        call setf(nskipr,f1,f2)
        write (6,*) 'f matrix has been set up, time=',
     c              dtime(dummy),dummy(1),dummy(2)
      else
        write (6,*) 'start reading f matrix'
        time1=dtime(dummy)
        open (21,file=tmpdir//'savef',form='unformatted',
     c        status='old')
        read (21) ((f1(i,j),i=1,nmodes),j=1,nrad1)
        read (21) ((f2(i,j),i=1,nmodes),j=1,nrad1)
        close(21)
        write (6,*) 'f matrix has been read, time=',
     c              dtime(dummy),dummy(1),dummy(2)
      end if
      if (iwrf.ne.0) then
        write (6,*) 'start writing f matrix'
        time1=dtime(dummy)
        open (21,file=tmpdir//'savef',form='unformatted',
     c        status='unknown')
        write (21) ((f1(i,j),i=1,nmodes),j=1,nrad1)
        write (21) ((f2(i,j),i=1,nmodes),j=1,nrad1)
        close(21)
        write (6,*) 'f matrix has been written, time=',
     c              dtime(dummy),dummy(1),dummy(2)
      end if
      if (mod(nsetp-1,nskipt).ne.0) then
        write (6,*) 'skip to get av. kernel mesh in latitude does not'
        write (6,*) 'divide into original number of meshpoints-1'
        stop
      end if
      nsetp1=(nsetp-1)/nskipt+1
      if (irdg.eq.0) then
        write (6,*) 'start setting up g matrix'
        time1=dtime(dummy)
        call setg(0,nsetp,nsetp1,g1,g2,tmesh1)
        write (6,*) 'g matrix has been set up, time=',
     c              dtime(dummy),dummy(1),dummy(2)
      else
        write (6,*) 'start reading g matrix'
        time1=dtime(dummy)
        open (21,file=tmpdir//'saveg',form='unformatted',
     c        status='old')
        read (21) (tmesh1(i),i=1,nsetp1)
c       read (21) nlm,(lplm(i),mplm(i),i=1,nlm)
cc      read (21) lmax,((iplm(i,j),i=0,lmax),j=0,lmax)
        read (21) ((g1(i,j),i=1,nlm),j=1,nsetp1)
        read (21) ((g2(i,j),i=1,nlm),j=1,nsetp1)
        close(21)
        write (6,*) 'g matrix has been read, time=',
     c              dtime(dummy),dummy(1),dummy(2)
      end if
      if (iwrg.ne.0) then
        write (6,*) 'start writing g matrix'
        time1=dtime(dummy)
        open (21,file=tmpdir//'saveg',form='unformatted',
     c        status='unknown')
        write (21) (tmesh1(i),i=1,nsetp1)
c       write (21) nlm,(lplm(i),mplm(i),i=1,nlm)
cc      write (21) lmax,((iplm(i,j),i=0,lmax),j=0,lmax)
        write (21) ((g1(i,j),i=1,nlm),j=1,nsetp1)
        write (21) ((g2(i,j),i=1,nlm),j=1,nsetp1)
        close(21)
        write (6,*) 'g matrix has been written, time=',
     c              dtime(dummy),dummy(1),dummy(2)
      end if
      call dsetint(1,nsetp1,tmesh1,tweights)
      npoints=nrad1*nsetp1
      nident=icenter*(nsetp1-1)
      npoi1=npoints-nident
      do 500,i2=1,npoints
        k2=mod(i2-1,nsetp1)+1
        j2=(i2-1)/nsetp1+1
        ijk(j2,k2)=max(1,i2-nident)
 500  continue
      do 600,i=1,nnlm
c       lmsplit(i)=iplm(msplit(i),lsplit(i))
 600  continue
      open (24,file='rmesh.orig',form='formatted')
      do 610,i=1,nrad
        write (24,*) radmesh(i)
 610  continue
      close (24)
      if (iavkern.ne.0) then
        open (24,file='avkern.log',form='formatted')
      end if
      open (25,file='rot.2d',form='formatted')
      if (icovar.eq.1) then
        open (26,file='err.2d',form='formatted')
      end if
      open (27,file='avkern.2d',form='unformatted')
      if (icoeff.eq.1) then
        open (28,file='coeff.2d',form='unformatted')
      end if
      xr(1)=-9.5
      xr(2)=-8.5
      xr(3)=-3.5
      xr(4)=-3.0
      xr(5)=-2.5
      xr(6)=-2.0
c     do 5000,imr=-8,-4,1
c     mur=10.**imr
c     do 5000,imr=1,6
c     mur=10.**xr(imr)
c     do 5000,imt=-7,-1,1
c     mut=10.**imt
c     write (6,*) mur,mut
      if (irdb.eq.0) then
        write (6,*) 'start setting up b matrix'
        time1=dtime(dummy)
        do 950,i=1,npoi1
          do 940,j=1,npoi1
            b(j,i)=0.0D0
 940      continue
 950    continue
	nb1=0
	nb2=0
	nb2a=0
	t2=0
	t2a=0
        do 1100,is1=1,nnlm,maxs1
          write (6,*) is1
	  nbig=0
	  nbig2=0
          is0=min(maxs1-1,nnlm-is1)
          ix=0
          do 1020,i2=1,npoints
	    ibig(i2)=.false.
            k2=mod(i2-1,nsetp1)+1
            j2=(i2-1)/nsetp1+1
            nbig1=0
            do 1010,isx=0,is0
              is=is1+isx
              i=iln(is)
              ilmi=ilm(is)
              ax(isx)=sigma1(is)*
     c          (f1(i,j2)*g1(ilmi,k2)+f2(i,j2)*g2(ilmi,k2))
	      if (abs(ax(isx)).gt.1e-15) then
		nbig1=nbig1+1
	      end if
 1010       continue
            if (nbig1.ne.0) then
	      ibig(i2)=.true.
              nbig=nbig+nbig1
              nbig2=nbig2+1
              ix=ix+1
              iax(ix)=i2
              do 1015,isx=0,is0
                a(isx,ix)=ax(isx)
 1015         continue
            end if
 1020     continue
          nx=ix
	  nb1=nb1+nbig
	  nb2=nb2+nbig2
	  nb2a=nb2a+npoints
	  t2=t2+nbig2**2*(is0+1.)
	  t2a=t2a+npoints**2*(is0+1.)
          do 1080,i1x=1,nx,nblock
            n1x=min(i1x+nblock-1,nx)-i1x+1
            do 1070,i2x=i1x,nx,nblock
              n2x=min(i2x+nblock-1,nx)-i2x+1
              call dgemm('t','n',n2x,n1x,is0+1,1.0D0,a(0,i2x),maxs1,
     c                   a(0,i1x),maxs1,0.0D0,bx,nblock)
              do 1060,i1=i1x,min(i1x+nblock-1,nx)
                i1a=max(1,iax(i1)-nident)
                do 1050,i2=i2x,min(i2x+nblock-1,nx)
                  i2a=max(1,iax(i2)-nident)
                  if (i2a.ge.i1a) then
c                   b(i2a,i1a)=b(i2a,i1a)+ddot(is0+1,a(0,i2),1,a(0,i1),1)
                    b(i2a,i1a)=b(i2a,i1a)+bx(i2-i2x+1,i1-i1x+1)
                  end if
 1050           continue
 1060         continue
 1070       continue
 1080     continue
 1100   continue
	write (6,*) (nb1+0.0)/(nnlm*npoints),(nb2+0.0)/nb2a,t2/t2a
        do 1290,i2=1,npoi1
          do 1280,i1=1,i2-1
            b(i1,i2)=b(i2,i1)
 1280     continue
 1290   continue
        write (6,*) 'b matrix has been set up, time=',
     c              dtime(dummy),dummy(1),dummy(2)
      else
        write (6,*) 'start reading b matrix'
        time1=dtime(dummy)
        open (21,file=tmpdir//'saveb',form='unformatted',
     c        status='old')
        read (21) ((b(j,i),j=1,npoi1),i=1,npoi1)
        close(21)
c       do i=1,npoi1
c         do j=1,npoi1
c           b(j,i)=real(b(j,i))
c         end do
c       end do
        write (6,*) 'b matrix has been read, time=',
     c              dtime(dummy),dummy(1),dummy(2)
      end if
c     write (6,*) (b(i,i),i=1,npoi1)
      if (iwrb.ne.0) then
        write (6,*) 'start writing b matrix'
        time1=dtime(dummy)
        open (21,file=tmpdir//'saveb',form='unformatted',
     c        status='unknown')
        write (21) ((b(j,i),j=1,npoi1),i=1,npoi1)
        close(21)
        write (6,*) 'b matrix has been written, time=',
     c              dtime(dummy),dummy(1),dummy(2)
c     to prevent recalculating b if doing trade-off curves
        irdb=1
      end if
      write (6,*) 'start adding regularization'
      time1=dtime(dummy)
      do 1310,i=1,npoi1
        do 1300,j=1,npoi1
          d(j,i)=b(j,i)
 1300   continue
 1310 continue
c     c1 is is the regularization in the 1-d radial case
c     times the radius
      if (nreg.eq.2) then
c       set up weighting functions
        do 1400,i=1,nrad1
c         fr(i)=1.0
          fr(i)=rmesh1(i)
          if (rmesh1(i).eq.0) then
            ft(i)=0.0
          else
            ft(i)=1.0D0/rmesh1(i)**4
          end if
c just a try
c         ft(i)=1.0
 1400   continue
        do 1405,i=1,nrad1
          fwr(i)=mur*rweights(i)*fr(i)
c         fwr(i)=mur*rweights(i)*rmesh1(i)
 1405   continue
        do 1430,j=1,nrad1
          do 1420,k=j,min(nrad1,j+2)
            ckj=0.0D0
            do 1410,i=max(2,j-1,k-1),min(nrad1-1,j+1,k+1)
              xm1=rmesh1(i-1)-rmesh1(i)
              xp1=rmesh1(i+1)-rmesh1(i)
              den=xm1*xp1*(xm1-xp1)
              if (i.eq.j-1) fji=-xm1
              if (i.eq.j)   fji=xm1-xp1
              if (i.eq.j+1) fji=xp1
              if (i.eq.k-1) fki=-xm1
              if (i.eq.k)   fki=xm1-xp1
              if (i.eq.k+1) fki=xp1
              ckj=ckj+fji*fki*fwr(i)/den**2
 1410       continue
c factor of 4=2^2 since second derivative is 2*a where y=a x^2+b x +c
            c1(k,j)=4.*ckj
            c1(j,k)=4.*ckj
 1420     continue
 1430   continue
        do 1500,i=1,nsetp1
          fwt(i)=mut*tweights(i)
 1500   continue
        do 1530,j=1,nsetp1
          do 1520,k=j,min(nsetp1,j+2)
            ckj=0.0D0
            do 1510,i=max(1,j-1,k-1),min(nsetp1,j+1,k+1)
              t=tmesh1(i)
              if (i.gt.1) tm1=tmesh1(i-1)
              if (i.eq.1) tm1=2*tmesh1(1)-tmesh1(2)
              if (i.lt.nsetp1) tp1=tmesh1(i+1)
              if (i.eq.nsetp1) tp1=2*tmesh1(nsetp1)-tmesh1(nsetp1-1)
              xm1=tm1-t
              xp1=tp1-t
              den=xm1*xp1*(xm1-xp1)
              if (i.eq.j-1) fji=-xm1
              if (i.eq.j)   fji=xm1-xp1
              if (i.eq.j+1) fji=xp1
              if (i.eq.1.and.j.eq.2) fji=xp1-xm1
              if (i.eq.nsetp1.and.j.eq.(nsetp1-1)) fji=xp1-xm1
              if (i.eq.k-1) fki=-xm1
              if (i.eq.k)   fki=xm1-xp1
              if (i.eq.k+1) fki=xp1
              if (i.eq.1.and.k.eq.2) fki=xp1-xm1
              if (i.eq.nsetp1.and.k.eq.(nsetp1-1)) fki=xp1-xm1
              ckj=ckj+fji*fki*fwt(i)/den**2
c             write (6,*) i,j,k,ckj
 1510       continue
c factor of 4=2^2 since second derivative is 2*a where y=a x^2+b x +c
            c2(k,j)=4.*ckj
            c2(j,k)=4.*ckj
 1520     continue
 1530   continue
      end if
      write (6,*) dtime(dummy),dummy(1),dummy(2)
c     first the regularization in theta is applied
      do 1650,j=2,nrad1
        wr3i=rweights(j)*ft(j)
c       wr3i=rweights(j)/rmesh1(j)**4
        do 1640,k1=1,nsetp1
          i1=ijk(j,k1)
          do 1630,k2=max(1,k1-2),min(nsetp1,k1+2)
            i2=ijk(j,k2)
            d(i2,i1)=d(i2,i1)+c2(k2,k1)*wr3i
 1630     continue
 1640   continue
 1650 continue
      write (6,*) dtime(dummy),dummy(1),dummy(2)
c     then the regularization in radius is added
      do 1750,k=1,nsetp1
        w=tweights(k)
        do 1740,j1=1,nrad1
          i1=ijk(j1,k)
          do 1730,j2=1,nrad1
            i2=ijk(j2,k)
            d(i2,i1)=d(i2,i1)+c1(j2,j1)*w
 1730     continue
 1740   continue
 1750 continue
      write (6,*) 'regularization has been added, time=',
     c            dtime(dummy),dummy(1),dummy(2)
      write (6,*) 'start setting up right hand side'
      time1=dtime(dummy)
      do 1810,i=1,npoi1
        right(i)=0.0D0
 1810 continue
      do 1850,i1=1,npoints
c       i1a=max(1,i1-nident)
        k1=mod(i1-1,nsetp1)+1
        j1=(i1-1)/nsetp1+1
        do 1820,i=1,nmodes
          f1x(i)=f1(i,j1)
          f2x(i)=f2(i,j1)
 1820   continue
        do 1830,i=1,nlm
          g1x(i)=g1(i,k1)
          g2x(i)=g2(i,k1)
 1830   continue
        r=0.0D0
        do 1840,is=1,nnlm
          i=iln(is)
          ilmi=ilm(is)
          r=r+ss2(is)*(f1x(i)*g1x(ilmi)+f2x(i)*g2x(ilmi))
 1840   continue
c       right(i1a)=right(i1a)+r
        help(i1)=r
 1850 continue
      do 1860,i1=1,npoints
        i1a=max(1,i1-nident)
        right(i1a)=right(i1a)+help(i1)
 1860 continue
      write (6,*) 'right hand side has been set up, time=',
     c            dtime(dummy),dummy(1),dummy(2)
c     open (21,file=tmpdir//'saved',form='unformatted',
c    c      status='unknown')
c     do i=1,npoi1
c       write (21) (d(j,i),j=1,npoi1)
c     end do
c     write (21) (right(i),i=1,npoi1)
c     close(21)
c     stop
      write (6,*) 'start decomposing d matrix'
      time1=dtime(dummy)
c     call dpoco(d,maxpoints,npoi1,rcond,help,info)
      anorm=dlansy('1','U',npoi1,d,maxpoints,help)
      call dpotrf('U',npoi1,d,maxpoints,info)
      write (6,*) 'd matrix has been decomposed, time=',
     c            dtime(dummy),dummy(1),dummy(2),info
      write (6,*) 'start finding condition number'
      call dpocon('U',npoi1,d,maxpoints,anorm,rcond,help,iwork,info1)
      write (6,*) 'condition number has been found, time=',
     c            dtime(dummy),dummy(1),dummy(2)
      write (6,*) rcond,info,info1,anorm
      write (6,*) 'start finding solution'
      time1=dtime(dummy)
c     call dposl(d,maxpoints,npoi1,right)
      call dpotrs('U',npoi1,1,d,maxpoints,right,maxpoints,info)
      write (6,*) 'solution has been found, time=',
     c            dtime(dummy),dummy(1),dummy(2)
c     write (21) (right(i),i=1,npoi1)
c     close(21)
c     stop
      write (6,*) 'rcond,info',rcond,info
      do 1875,j=1,nrad1
        do 1870,k=1,nsetp1
          omega(j,k)=right(ijk(j,k))
 1870   continue
 1875 continue
      do 1877,j=1,nrad1
        write (25,'(101f12.4)') (right(ijk(j,k)),k=1,nsetp1)
 1877 continue
      write (6,*) 'start calculating splittings'
      time1=dtime(dummy)
c set ixl to the first mode with a given (l,n)
c set ixh to the last mode with a given (l,n)
      ilast=-1
      ix=0
      do 1880,i1=1,nnlm
        i=iln(i1)
        if (i.ne.ilast) then
          ix=ix+1
          ixl(ix)=i1
          ilast=i
        end if
        ixh(ix)=i1
 1880 continue
      nx=ix
      chisq=0.0
c*$* assert do (serial)
      open(10,file='splittings.out',form='formatted')
      do 1897,ix=1,nx
c loop first over each ln
        i=iln(ixl(ix))
        do 1891,j=1,nrad1
          f1i(j)=f1(i,j)
          f2i(j)=f2(i,j)
 1891   continue
        call dgemv('t',nrad1,nsetp1,1.0D0,omega,maxr1,
     c             f1i,1,0.0D0,f1o,1)
        call dgemv('t',nrad1,nsetp1,1.0D0,omega,maxr1,
     c             f2i,1,0.0D0,f2o,1)
c now do the individual nlm's
        do 1895,i1=ixl(ix),ixh(ix)
          index=ilm(i1)
c         do 1893,k=1,nsetp1
c           g1i(k)=g1(index,k)
c           g2i(k)=g2(index,k)
c1893     continue
c         s=ddot(nsetp1,f1o,1,g1i,1)+ddot(nsetp1,f2o,1,g2i,1)
          s=0.0D0
          do 1893,k=1,nsetp1
            s=s+g1(index,k)*f1o(k)+g2(index,k)*f2o(k)
 1893     continue
          l=lmode(i)
          n=nmode(i)
          write (10,'(i4,i4,f8.2,i4,2i3,3f12.4)')
c         write (6,'(i4,i4,f8.2,i4,2i3,3f12.4)') 
     c          l,n,freq(i),mnlm(i1),itnlm(i1),it1nlm(i1),
     c          s,split(i1),sigma(i1)
          chisq=chisq+((s-split(i1))*sigma1(i1))**2
 1895   continue
 1897 continue
      write (6,*) 'splittings have been calculated, time=',
     c            dtime(dummy),dummy(1),dummy(2)
      write (6,*) 'start calculating wiggle measures'
      time1=dtime(dummy)
      do 1907,k=1,nsetp1
        do 1905,j=1,nrad1
          fro(j,k)=0.0
          fto(j,k)=0.0
c         omega(j,k)=sin(tmesh1(k))**2*rmesh1(j)**2
c         omega(j,k)=rmesh1(j)**3
 1905   continue
 1907 continue
      do 1930,j=1,nrad1
        do 1920,i=max(2,j-1),min(nrad1-1,j+1)
          xm1=rmesh1(i-1)-rmesh1(i)
          xp1=rmesh1(i+1)-rmesh1(i)
          den=xm1*xp1*(xm1-xp1)
          if (i.eq.j-1) fji=-xm1
          if (i.eq.j)   fji=xm1-xp1
          if (i.eq.j+1) fji=xp1
          do 1910,k=1,nsetp1
            fro(i,k)=fro(i,k)+omega(j,k)*2*fji/den
 1910     continue
 1920   continue
 1930 continue
      do 1960,k=1,nsetp1
        do 1950,i=max(1,k-1),min(nsetp1,k+1)
          t=tmesh1(i)
          if (i.gt.1) tm1=tmesh1(i-1)
          if (i.eq.1) tm1=2*tmesh1(1)-tmesh1(2)
          if (i.lt.nsetp1) tp1=tmesh1(i+1)
          if (i.eq.nsetp1) tp1=2*tmesh1(nsetp1)-tmesh1(nsetp1-1)
          xm1=tm1-t
          xp1=tp1-t
          den=xm1*xp1*(xm1-xp1)
          if (i.eq.k-1) fki=-xm1
          if (i.eq.k)   fki=xm1-xp1
          if (i.eq.k+1) fki=xp1
          if (i.eq.1.and.k.eq.2) fki=xp1-xm1
          if (i.eq.nsetp1.and.k.eq.(nsetp1-1)) fki=xp1-xm1
          do 1940,j=1,nrad1
            fto(j,i)=fto(j,i)+omega(j,k)*2*fki/den
 1940     continue
 1950   continue
 1960 continue
c     calculate wigr and wigt
      wigr=0.0
      wigt=0.0
      do 1980,k=1,nsetp1
        wt=tweights(k)
        do 1970,j=1,nrad1
          wigr=wigr+wt*rweights(j)*fr(j)*fro(j,k)**2
          wigt=wigt+wt*rweights(j)*ft(j)*fto(j,k)**2
 1970   continue
 1980 continue
      write (6,*) 'wiggle measures have been computed, time=',
     c            dtime(dummy),dummy(1),dummy(2)
      write (6,*) 'chisq=',chisq
      write (6,*) 'rms=',sqrt(chisq/nnlm)
      write (6,*) 'wigr, wigt=',wigr,wigt
      write (6,*) '%',mur,mut,chisq,sqrt(chisq/nnlm),wigr,wigt,
     c            info,rcond
c     end do
c     end do
c     stop
      if ((icovar.eq.1).or.(iavkern.ge.1)) then
        write (6,*) 'start inverting d matrix'
        time1=dtime(dummy)
c       call dpodi(d,maxpoints,npoi1,det,1)
        call dpotri('U',npoi1,d,maxpoints,info)
c       open (21,file=tmpdir//'saved1',form='unformatted',
c    c        status='unknown',access='direct',recl=8*npoi1)
        do 2095,i=1,npoi1
c         do 2092,j=1,npoi1
c           di(j)=0.0D0
c2092     continue
c         di(i)=1.0D0
c         call dposl(d,maxpoints,npoi1,di)
c         write (21,rec=i) (di(j),j=1,npoi1)
c         write (21,rec=i) (d(j,i),j=1,i),(d(i,j),j=i+1,npoi1)
          do 2092,j=i+1,npoi1
            d(j,i)=d(i,j)
 2092     continue
 2095   continue
c       close(21)
        write (6,*) 'd matrix has been inverted, time=',
     c              dtime(dummy),dummy(1),dummy(2)
      end if
      if (icovar.ge.1) then
        write (6,*) 'start finding variances'
c the following loop could be blocked
c       do 2150,i=1,npoi1
c         call dgemv('n',npoi1,npoi1,1.0D0,b,maxpoints,d(1,i),1,
c    c               0.0D0,dib,1)
c         var(i)=ddot(npoi1,d(1,i),1,dib,1)
c2150   continue
        do 2150,i=1,npoi1,nblock1
          ix=min(npoi1,i+nblock1-1)
          nx=ix+1-i
          call dgemm('n','n',npoi1,nx,npoi1,1.0D0,b,maxpoints,
     c               d(1,i),maxpoints,0.0D0,dibx,maxpoints)
          do 2140,j=i,ix
            var(j)=ddot1(npoi1,d(1,j),dibx(1,j-i+1))
 2140     continue
 2150   continue
        do 2190,j=1,nrad1
          write (26,'(100f11.5)') (sqrt(var(ijk(j,k))),k=1,nsetp1)
 2190   continue
c       close(21)
        write (6,*) 'variances have been found, time=',
     c              dtime(dummy),dummy(1),dummy(2)
c     else
c       write (6,*) 'start finding covariance matrix'
c       time1=dtime(dummy)
c       do 2230,i=1,npoi1
c         do 2220,j=1,npoi1
c           d(j,i)=ddot1(npoi1,b(1,j),d1(1,i))
c2220     continue
c2230   continue
c       do 2250,i=1,npoi1
c         do 2240,j=i,npoi1
c           cov(j,i)=ddot1(npoi1,d1(1,j),d(1,i))
c2240     continue
c2250   continue
c       do 2270,i=1,npoi1
c         do 2260,j=1,i-1
c           cov(j,i)=cov(i,j)
c2260     continue
c2270   continue
c       do 2280,i=1,nrad1
c         write (6,*) rmesh1(i)
c2280   continue
c       do 2290,j=1,nrad1
c         write (6,'(100f9.5)') (sqrt(cov(ijk(j,k),ijk(j,k))),
c    c      k=1,nsetp1)
c2290   continue
c       write (6,*) 'covariance matrix has been found, time=',
c    c              dtime(dummy),dummy(1),dummy(2)
      end if
      if (iavkern.eq.0) then
        write (6,*) 'skipping calculation of averaging kernels'
      else
        write (6,*) 'start finding averaging kernels'
        time1=dtime(dummy)
        write (6,*) nt2,nr2
c       open (21,file=tmpdir//'saved1',form='unformatted',
c    c        status='old',access='direct',recl=8*npoi1)
c moved up
c       open (27,file='avkern.2d',form='unformatted')
c       open (28,file='coeff',form='unformatted')
        do 2300,j1=1,nr2
          x2(j1)=radmesh(nsr2*(j1-1)+1)
 2300   continue
        pi=4.0D0*atan(1.0D0)
        pi2=pi/2.0D0
        do 2301,k1=1,nt2
          x1(k1)=pi2*(k1-1)/(nt2-1.)
 2301   end do
        write (6,*) nrad1,nsr,nsetp1,nst
        if (nsr.ne.0) then
          i=0
          do 2303,j=1,nrad1,nsr
            do 2302,k=1,nsetp1,nst
              i=i+1
              jav(i)=j
              kav(i)=k
 2302       continue
 2303     continue
          nav=i
        end if
c       do 2380,j=1,nrad1,nsr
c         do 2370,k=1,nsetp1,nst
        do 2380,iav=1,nav
          j=jav(iav)
          k=kav(iav)
            ipoi=ijk(j,k)
c           read (21,rec=ipoi) (di(j1),j1=1,npoi1)
            do 2305,is=1,nnlm
              coeff(is)=0.0D0
 2305       continue
            do 2330,i1=1,npoints
              i1a=max(1,i1-nident)
              k1=mod(i1-1,nsetp1)+1
              j1=(i1-1)/nsetp1+1
              do 2310,i=1,nmodes
                f1x(i)=f1(i,j1)
                f2x(i)=f2(i,j1)
 2310         continue
              do 2315,i=1,nlm
                g1x(i)=g1(i,k1)
                g2x(i)=g2(i,k1)
 2315         continue
c             d1i=d1(i1a,ipoi)
c             d1i=di(i1a)
              d1i=d(i1a,ipoi)
              do 2320,is=1,nnlm
                i=iln(is)
                ilmi=ilm(is)
                coeff(is)=coeff(is)+d1i*sigma1(is)*
     c            (f1x(i)*g1x(ilmi)+f2x(i)*g2x(ilmi))
 2320         continue
 2330       continue
c           s1=0.0D0
c           s2=0.0D0
c           s3=0.0D0
            vari=0.0
            do 2350,is=1,nnlm
              coeff(is)=coeff(is)*sigma1(is)
c             s1=s1+coeff(is)*split(is)
c             s2=s2+mnlm(is)*coeff(is)
c             s3=s3+coeff(is)
              vari=vari+(coeff(is)*sigma(is))**2
 2350       continue
c           write (6,*) '#',j,k
c           write (6,*) s1,right(ipoi),s1-right(ipoi),s2,s3
c           call set2davk(nnlm,isplit,msplit,coeff,avkern,nsr2,nt2)
            call set2davk(coeff,avkern,nsr2,nt2)
            if (iavkern.eq.2) then
              write (27) ((avkern(k1,j1),k1=1,nt2),j1=1,nr2)
            end if
c           call fitmn(x1,x2,nt2,nr2,avkern,yy2,mjta,4,
c    c                 real(tmesh1(k)),real(rmesh1(j)))
            write (24,*) mur,mut,real(rmesh1(j)),real(tmesh1(k)),
     c                  right(ipoi),sqrt(vari),
     c                  mjta(3),mjta(1),mjta(4),mjta(2)
            if (icoeff.eq.1) then
              write (28) (real(coeff(is)),is=1,nnlm)
            end if
c           end if
c2370     continue
 2380   continue
c       close(27)
c       close(28)
c       close(21)
        write (6,*) 'averaging kernels set up, time=',
     c              dtime(dummy),dummy(1),dummy(2)
      end if
 5000 continue
      if (iavkern.ne.0) then
        close (24)
      end if
      close (25)
      if (icovar.eq.26) then
        close (26)
      end if
      close (27)
      if (icoeff.eq.28) then
        close (28)
      end if
      open (20,file='status.log',form='formatted')
      write (20,*) 0
      close (20)

      write (6,*) 'total time used:',etime(dummy)
      write (6,*) 'user time:      ',dummy(1)
      write (6,*) 'system time:    ',dummy(2)
      end