diff --git a/libm65/JT65code.f90 b/libm65/JT65code.f90
deleted file mode 100644
index b232fcdd6..000000000
--- a/libm65/JT65code.f90
+++ /dev/null
@@ -1,47 +0,0 @@
-program JT65code
-
-! Provides examples of message packing, bit and symbol ordering,
-! Reed Solomon encoding, and other necessary details of the JT65
-! protocol.
-
-  character*22 msg0,msg,decoded,cok*3
-  integer dgen(12),sent(63),recd(12),era(51)
-
-  nargs=iargc()
-  if(nargs.ne.1) then
-     print*,'Usage: JT65code "message"'
-     go to 999
-  endif
-
-  call getarg(1,msg0)                     !Get message from command line
-  msg=msg0
-
-  call chkmsg(msg,cok,nspecial,flip)      !See if it includes "OOO" report
-  if(nspecial.gt.0) then                  !or is a shorthand message
-     write(*,1010) 
-1010 format('Shorthand message.')
-     go to 999
-  endif
-
-  call packmsg(msg,dgen)                  !Pack message into 72 bits
-  write(*,1020) msg0
-1020 format('Message:   ',a22)            !Echo input message
-  if(iand(dgen(10),8).ne.0) write(*,1030) !Is plain text bit set?
-1030 format('Plain text.')         
-  write(*,1040) dgen
-1040 format('Packed message, 6-bit symbols: ',12i3) !Display packed symbols
-
-  call rs_encode(dgen,sent)               !RS encode
-  call interleave63(sent,1)               !Interleave channel symbols
-  call graycode(sent,63,1)                !Apply Gray code
-  write(*,1050) sent
-1050 format('Channel symbols, including FEC:'/(i5,20i3))
-
-  call graycode(sent,63,-1)
-  call interleave63(sent,-1)
-  call rs_decode(sent,era,0,recd,nerr)
-  call unpackmsg(recd,decoded)            !Unpack the user message
-  write(*,1060) decoded,cok
-1060 format('Decoded message: ',a22,2x,a3)
-
-999 end program JT65code
diff --git a/libm65/afc65b.f b/libm65/afc65b.f
deleted file mode 100644
index d2e3a0823..000000000
--- a/libm65/afc65b.f
+++ /dev/null
@@ -1,68 +0,0 @@
-      subroutine afc65b(cx,cy,npts,fsample,nflip,ipol,xpol,a,
-     +  ccfbest,dtbest)
-
-      logical xpol
-      complex cx(npts)
-      complex cy(npts)
-      real a(5),deltaa(5)
-
-      a(1)=0.
-      a(2)=0.
-      a(3)=0.
-      a(4)=45.0*(ipol-1.0)
-      deltaa(1)=2.0
-      deltaa(2)=2.0
-      deltaa(3)=2.0
-      deltaa(4)=22.5
-      deltaa(5)=0.05
-      nterms=3
-      if(xpol) nterms=4
-      chisqr=0.
-
-C  Start the iteration
-      chisqr0=1.e6
-      do iter=1,3                               !One iteration is enough?
-         do j=1,nterms
-            chisq1=fchisq(cx,cy,npts,fsample,nflip,a,ccfmax,dtmax)
-            fn=0.
-            delta=deltaa(j)
- 10         a(j)=a(j)+delta
-            chisq2=fchisq(cx,cy,npts,fsample,nflip,a,ccfmax,dtmax)
-            if(chisq2.eq.chisq1) go to 10
-            if(chisq2.gt.chisq1) then
-               delta=-delta                      !Reverse direction
-               a(j)=a(j)+delta
-               tmp=chisq1
-               chisq1=chisq2
-               chisq2=tmp
-            endif
- 20         fn=fn+1.0
-            a(j)=a(j)+delta
-            chisq3=fchisq(cx,cy,npts,fsample,nflip,a,ccfmax,dtmax)
-            if(chisq3.lt.chisq2) then
-               chisq1=chisq2
-               chisq2=chisq3
-               go to 20
-            endif
-
-C  Find minimum of parabola defined by last three points
-            delta=delta*(1./(1.+(chisq1-chisq2)/(chisq3-chisq2))+0.5)
-            a(j)=a(j)-delta
-            deltaa(j)=deltaa(j)*fn/3.
-         enddo
-         chisqr=fchisq(cx,cy,npts,fsample,nflip,a,ccfmax,dtmax)
-         if(chisqr/chisqr0.gt.0.9999) go to 30
-         chisqr0=chisqr
-      enddo
-
- 30   ccfbest=ccfmax * (1378.125/fsample)**2
-      dtbest=dtmax
-
-      if(a(4).lt.0.0) a(4)=a(4)+180.0
-      if(a(4).ge.180.0) a(4)=a(4)-180.0
-      if(nint(a(4)).eq.180) a(4)=0.
-      ipol=nint(a(4)/45.0) + 1
-      if(ipol.gt.4) ipol=ipol-4
-
-      return
-      end
diff --git a/libm65/alignmsg.f90 b/libm65/alignmsg.f90
deleted file mode 100644
index 4608c2da3..000000000
--- a/libm65/alignmsg.f90
+++ /dev/null
@@ -1,38 +0,0 @@
-subroutine alignmsg(word0,nmin,msg,msglen,idone)
-
-  character*(*) word0
-  character*29 msg,word
-
-  word=word0//' '
-  idone=0
-
-! Test for two (or more) <space> characters
-  if(word(1:2).eq.'  ' .and. len(word).eq.2) then
-     i2=index(msg,'  ')
-     if((i2.ge.1.and.i2.lt.msglen) .or.                                 &
-          (msg(1:1).eq.' '.and.msg(msglen:msglen).eq.' ')) then
-        if(i2.eq.1) msg=msg(i2+2:msglen)           !Align on EOM
-        if(i2.ge.2) msg=msg(i2+2:msglen)//msg(1:i2-1)
-        idone=1
-     endif
-
-! Align on single <space> (as last resort)
-  else if(word(1:1).eq.' ' .and. len(word).eq.1) then
-     i3=index(msg,' ')
-     if(i3.ge.1 .and. i3.lt.msglen) msg=msg(i3+1:msglen)//msg(1:i3)
-     if(i3.eq.msglen) msg=msg(1:msglen)
-     msg=msg(1:msglen)//msg(1:msglen)
-     idone=1
-
-! Align on specified word
-  else
-     call match(word,msg(1:msglen),nstart,nmatch)
-     if(nmatch.ge.nmin) then
-        if(nstart.eq.1) msg=msg(nstart:msglen)
-        if(nstart.gt.1) msg=msg(nstart:msglen)//msg(1:nstart-1)
-        idone=1
-     endif
-  endif
-
-  return
-end subroutine alignmsg
diff --git a/libm65/astro.f b/libm65/astro.f
deleted file mode 100644
index a92b6ad12..000000000
--- a/libm65/astro.f
+++ /dev/null
@@ -1,109 +0,0 @@
-      subroutine astro(nyear,month,nday,uth,nfreq,Mygrid,
-     +     NStation,MoonDX,AzSun,ElSun,AzMoon0,ElMoon0,
-     +     ntsky,doppler00,doppler,dbMoon,RAMoon,DecMoon,HA,Dgrd,sd,
-     +     poloffset,xnr,day,lon,lat,LST)
-
-C  Computes astronomical quantities for display and tracking.
-C  NB: may want to smooth the Tsky map to 10 degrees or so.
-
-      character*6 MyGrid,HisGrid
-      real LST
-      real lat,lon
-      integer*2 nt144(180)
-
-!      common/echo/xdop(2),techo,AzMoon,ElMoon,mjd
-      real xdop(2)
-
-      data rad/57.2957795/
-      data nt144/
-     +  234, 246, 257, 267, 275, 280, 283, 286, 291, 298,
-     +  305, 313, 322, 331, 341, 351, 361, 369, 376, 381, 
-     +  383, 382, 379, 374, 370, 366, 363, 361, 363, 368, 
-     +  376, 388, 401, 415, 428, 440, 453, 467, 487, 512, 
-     +  544, 579, 607, 618, 609, 588, 563, 539, 512, 482, 
-     +  450, 422, 398, 379, 363, 349, 334, 319, 302, 282, 
-     +  262, 242, 226, 213, 205, 200, 198, 197, 196, 197, 
-     +  200, 202, 204, 205, 204, 203, 202, 201, 203, 206, 
-     +  212, 218, 223, 227, 231, 236, 240, 243, 247, 257, 
-     +  276, 301, 324, 339, 346, 344, 339, 331, 323, 316, 
-     +  312, 310, 312, 317, 327, 341, 358, 375, 392, 407, 
-     +  422, 437, 451, 466, 480, 494, 511, 530, 552, 579, 
-     +  612, 653, 702, 768, 863,1008,1232,1557,1966,2385, 
-     + 2719,2924,3018,3038,2986,2836,2570,2213,1823,1461, 
-     + 1163, 939, 783, 677, 602, 543, 494, 452, 419, 392, 
-     +  373, 360, 353, 350, 350, 350, 350, 350, 350, 348, 
-     +  344, 337, 329, 319, 307, 295, 284, 276, 272, 272, 
-     +  273, 274, 274, 271, 266, 260, 252, 245, 238, 231/
-      save
-
-      call grid2deg(MyGrid,elon,lat)
-      lon=-elon
-      call sun(nyear,month,nday,uth,lon,lat,RASun,DecSun,LST,
-     +    AzSun,ElSun,mjd,day)
-
-      freq=nfreq*1.e6
-      if(nfreq.eq.2) freq=1.8e6
-      if(nfreq.eq.4) freq=3.5e6
-
-      call MoonDop(nyear,month,nday,uth,lon,lat,RAMoon,DecMoon,
-     +  LST,HA,AzMoon,ElMoon,vr,dist)
-
-C  Compute spatial polarization offset
-      xx=sin(lat/rad)*cos(ElMoon/rad) - cos(lat/rad)*
-     +     cos(AzMoon/rad)*sin(ElMoon/rad)
-      yy=cos(lat/rad)*sin(AzMoon/rad)
-      if(NStation.eq.1) poloffset1=rad*atan2(yy,xx)
-      if(NStation.eq.2) poloffset2=rad*atan2(yy,xx)
-
-      techo=2.0 * dist/2.99792458e5                 !Echo delay time
-      doppler=-freq*vr/2.99792458e5                 !One-way Doppler
-
-      call coord(0.,0.,-1.570796,1.161639,RAMoon/rad,DecMoon/rad,el,eb)
-      longecl_half=nint(rad*el/2.0)
-      if(longecl_half.lt.1 .or. longecl_half.gt.180) longecl_half=180
-      t144=nt144(longecl_half)
-      tsky=(t144-2.7)*(144.0/nfreq)**2.6 + 2.7      !Tsky for obs freq
-
-      xdop(NStation)=doppler
-      if(NStation.eq.2) then
-         HisGrid=MyGrid
-         go to 900
-      endif
-
-      doppler00=2.0*xdop(1)
-      doppler=xdop(1)+xdop(2)
-!      if(mode.eq.3) doppler=2.0*xdop(1)
-      dBMoon=-40.0*log10(dist/356903.)
-      sd=16.23*370152.0/dist
-
-!      if(NStation.eq.1 .and. MoonDX.ne.0 .and. 
-!     +    (mode.eq.2 .or. mode.eq.5)) then
-      if(NStation.eq.1 .and. MoonDX.ne.0) then
-         poloffset=mod(poloffset2-poloffset1+720.0,180.0)
-         if(poloffset.gt.90.0) poloffset=poloffset-180.0
-         x1=abs(cos(2*poloffset/rad))
-         if(x1.lt.0.056234) x1=0.056234
-         xnr=-20.0*log10(x1)
-         if(HisGrid(1:1).lt.'A' .or. HisGrid(1:1).gt.'R') xnr=0
-      endif
-
-      tr=80.0                              !Good preamp
-      tskymin=13.0*(408.0/nfreq)**2.6      !Cold sky temperature
-      tsysmin=tskymin+tr
-      tsys=tsky+tr
-      dgrd=-10.0*log10(tsys/tsysmin) + dbMoon
- 900  AzMoon0=Azmoon
-      ElMoon0=Elmoon
-      ntsky=nint(tsky)
-
-!      auxHA = 15.0*(LST-auxra)                       !HA in degrees
-!      pi=3.14159265
-!      pio2=0.5*pi
-!      call coord(pi,pio2-lat/rad,0.0,lat/rad,auxha*pi/180.0,
-!     +  auxdec/rad,azaux,elaux)
-!      AzAux=azaux*rad
-!      ElAux=ElAux*rad
-
-      return
-
-      end
diff --git a/libm65/astro0.f90 b/libm65/astro0.f90
deleted file mode 100644
index ec234215e..000000000
--- a/libm65/astro0.f90
+++ /dev/null
@@ -1,81 +0,0 @@
-subroutine astro0(nyear,month,nday,uth8,nfreq,mygrid,hisgrid,              &
-     AzSun8,ElSun8,AzMoon8,ElMoon8,AzMoonB8,ElMoonB8,ntsky,ndop,ndop00,    &
-     dbMoon8,RAMoon8,DecMoon8,HA8,Dgrd8,sd8,poloffset8,xnr8,dfdt,dfdt0,    &
-     width1,width2,w501,w502,xlst8)
-
-  parameter (DEGS=57.2957795130823d0)
-  character*6 mygrid,hisgrid
-  real*8 AzSun8,ElSun8,AzMoon8,ElMoon8,AzMoonB8,ElMoonB8
-  real*8 dbMoon8,RAMoon8,DecMoon8,HA8,Dgrd8,xnr8,dfdt,dfdt0,dt
-  real*8 sd8,poloffset8,day8,width1,width2,w501,w502,xlst8
-  real*8 uth8
-  data uth8z/0.d0/
-  save
-
-  uth=uth8
-  call astro(nyear,month,nday,uth,nfreq,hisgrid,2,1,                 &
-       AzSun,ElSun,AzMoon,ElMoon,ntsky,doppler00,doppler,            &
-       dbMoon,RAMoon,DecMoon,HA,Dgrd,sd,poloffset,xnr,               &
-       day,xlon2,xlat2,xlst)
-  AzMoonB8=AzMoon
-  ElMoonB8=ElMoon
-  call astro(nyear,month,nday,uth,nfreq,mygrid,1,1,                  &
-       AzSun,ElSun,AzMoon,ElMoon,ntsky,doppler00,doppler,            &
-       dbMoon,RAMoon,DecMoon,HA,Dgrd,sd,poloffset,xnr,               &
-       day,xlon1,xlat1,xlst)
-
-  day8=day
-  xlst8=xlst
-  call tm2(day8,xlat1,xlon1,xl1,b1)
-  call tm2(day8,xlat2,xlon2,xl2,b2)
-  call tm2(day8+1.d0/1440.0,xlat1,xlon1,xl1a,b1a)
-  call tm2(day8+1.d0/1440.0,xlat2,xlon2,xl2a,b2a)
-  fghz=0.001*nfreq
-  dldt1=DEGS*(xl1a-xl1)
-  dbdt1=DEGS*(b1a-b1)
-  dldt2=DEGS*(xl2a-xl2)
-  dbdt2=DEGS*(b2a-b2)
-  rate1=2.0*sqrt(dldt1**2 + dbdt1**2)
-  width1=0.5*6741*fghz*rate1
-  rate2=sqrt((dldt1+dldt2)**2 + (dbdt1+dbdt2)**2)
-  width2=0.5*6741*fghz*rate2
-
-  fbend=0.7
-  a2=0.0045*log(fghz/fbend)/log(1.05)
-  if(fghz.lt.fbend) a2=0.0
-  f50=0.19 * (fghz/fbend)**a2
-  if(f50.gt.1.0) f50=1.0
-  w501=f50*width1
-  w502=f50*width2
-
-  AzSun8=AzSun
-  ElSun8=ElSun
-  AzMoon8=AzMoon
-  ElMoon8=ElMoon
-  dbMoon8=dbMoon
-  RAMoon8=RAMoon/15.0
-  DecMoon8=DecMoon
-  HA8=HA
-  Dgrd8=Dgrd
-  sd8=sd
-  poloffset8=poloffset
-  xnr8=xnr
-  ndop=nint(doppler)
-  ndop00=nint(doppler00)
-
-  if(uth8z.eq.0.d0) then
-     uth8z=uth8-1.d0/3600.d0
-     dopplerz=doppler
-     doppler00z=doppler00
-  endif
-     
-  dt=60.0*(uth8-uth8z)
-  if(dt.le.0) dt=1.d0/60.d0
-  dfdt=(doppler-dopplerz)/dt
-  dfdt0=(doppler00-doppler00z)/dt
-  uth8z=uth8
-  dopplerz=doppler
-  doppler00z=doppler00
-
-  return
-end subroutine astro0
diff --git a/libm65/astrosub.f90 b/libm65/astrosub.f90
deleted file mode 100644
index 64dddb00c..000000000
--- a/libm65/astrosub.f90
+++ /dev/null
@@ -1,14 +0,0 @@
-subroutine astrosub(nyear,month,nday,uth8,nfreq,mygrid,hisgrid,          &
-     AzSun8,ElSun8,AzMoon8,ElMoon8,AzMoonB8,ElMoonB8,ntsky,ndop,ndop00,  &
-     RAMoon8,DecMoon8,Dgrd8,poloffset8,xnr8)
-
-  implicit real*8 (a-h,o-z)
-  character*6 mygrid,hisgrid
-
-  call astro0(nyear,month,nday,uth8,nfreq,mygrid,hisgrid,                &
-     AzSun8,ElSun8,AzMoon8,ElMoon8,AzMoonB8,ElMoonB8,ntsky,ndop,ndop00,  &
-     dbMoon8,RAMoon8,DecMoon8,HA8,Dgrd8,sd8,poloffset8,xnr8,dfdt,dfdt0,  &
-     width1,width2,w501,w502,xlst8)
-
-  return  
-end subroutine astrosub
diff --git a/libm65/ccf2.f b/libm65/ccf2.f
deleted file mode 100644
index 090db5151..000000000
--- a/libm65/ccf2.f
+++ /dev/null
@@ -1,45 +0,0 @@
-      subroutine ccf2(ss,nz,nflip,ccfbest,lagpk)
-
-      parameter (LAGMAX=60)
-!      parameter (LAGMAX=200)
-      real ss(nz)
-      real ccf(-LAGMAX:LAGMAX)
-      integer npr(126)
-
-C  The JT65 pseudo-random sync pattern:
-      data npr/
-     + 1,0,0,1,1,0,0,0,1,1,1,1,1,1,0,1,0,1,0,0,
-     + 0,1,0,1,1,0,0,1,0,0,0,1,1,1,0,0,1,1,1,1,
-     + 0,1,1,0,1,1,1,1,0,0,0,1,1,0,1,0,1,0,1,1,
-     + 0,0,1,1,0,1,0,1,0,1,0,0,1,0,0,0,0,0,0,1,
-     + 1,0,0,0,0,0,0,0,1,1,0,1,0,0,1,0,1,1,0,1,
-     + 0,1,0,1,0,0,1,1,0,0,1,0,0,1,0,0,0,0,1,1,
-     + 1,1,1,1,1,1/
-      save
-
-      ccfbest=0.
-      lag1=-LAGMAX
-      lag2=LAGMAX
-      do lag=lag1,lag2
-         s0=0.
-         s1=0.
-         do i=1,126
-            j=2*(8*i + 43) + lag
-            if(j.ge.1 .and. j.le.nz-8) then
-               x=ss(j)+ss(j+8)             !Add two half-symbol contributions
-               if(npr(i).eq.0) then
-                  s0=s0 + x
-               else
-                  s1=s1 + x
-               endif
-            endif
-         enddo
-         ccf(lag)=nflip*(s1-s0)
-         if(ccf(lag).gt.ccfbest) then
-            ccfbest=ccf(lag)
-            lagpk=lag
-         endif
-      enddo
-
-      return
-      end
diff --git a/libm65/ccf65.f90 b/libm65/ccf65.f90
deleted file mode 100644
index 455b919c9..000000000
--- a/libm65/ccf65.f90
+++ /dev/null
@@ -1,118 +0,0 @@
-subroutine ccf65(ss,nhsym,ssmax,sync1,ipol1,jpz,dt1,flipk,syncshort,      &
-     snr2,ipol2,dt2)
-
-  parameter (NFFT=512,NH=NFFT/2)
-  real ss(4,322)                   !Input: half-symbol powers, 4 pol'ns
-  real s(NFFT)                     !CCF = ss*pr
-  complex cs(0:NH)                 !Complex FT of s
-  real s2(NFFT)                    !CCF = ss*pr2
-  complex cs2(0:NH)                !Complex FT of s2
-  real pr(NFFT)                    !JT65 pseudo-random sync pattern
-  complex cpr(0:NH)                !Complex FT of pr
-  real pr2(NFFT)                   !JT65 shorthand pattern
-  complex cpr2(0:NH)               !Complex FT of pr2
-  real tmp1(322)
-  real tmp2(322)
-  real ccf(-27:27,4)
-  logical first
-  integer npr(126)
-  data first/.true./
-  equivalence (s,cs),(pr,cpr),(s2,cs2),(pr2,cpr2)
-  save
-
-! The JT65 pseudo-random sync pattern:
-  data npr/                                        &
-      1,0,0,1,1,0,0,0,1,1,1,1,1,1,0,1,0,1,0,0,     &
-      0,1,0,1,1,0,0,1,0,0,0,1,1,1,0,0,1,1,1,1,     &
-      0,1,1,0,1,1,1,1,0,0,0,1,1,0,1,0,1,0,1,1,     &
-      0,0,1,1,0,1,0,1,0,1,0,0,1,0,0,0,0,0,0,1,     &
-      1,0,0,0,0,0,0,0,1,1,0,1,0,0,1,0,1,1,0,1,     &
-      0,1,0,1,0,0,1,1,0,0,1,0,0,1,0,0,0,0,1,1,     &
-      1,1,1,1,1,1/
-
-  if(first) then
-! Initialize pr, pr2; compute cpr, cpr2.
-     fac=1.0/NFFT
-     do i=1,NFFT
-        pr(i)=0.
-        pr2(i)=0.
-        k=2*mod((i-1)/8,2)-1
-        if(i.le.NH) pr2(i)=fac*k
-     enddo
-     do i=1,126
-        j=2*i
-        pr(j)=fac*(2*npr(i)-1)
-! Not sure why, but it works significantly better without the following line:
-!        pr(j-1)=pr(j)
-     enddo
-     call four2a(pr,NFFT,1,-1,0)
-     call four2a(pr2,NFFT,1,-1,0)
-     first=.false.
-  endif
-
-! Look for JT65 sync pattern and shorthand square-wave pattern.
-  ccfbest=0.
-  ccfbest2=0.
-  ipol1=1
-  ipol2=1
-  do ip=1,jpz                                  !Do jpz polarizations
-     do i=1,nhsym-1
-!        s(i)=ss(ip,i)+ss(ip,i+1)
-        s(i)=min(ssmax,ss(ip,i)+ss(ip,i+1))
-     enddo
-     s(nhsym:NFFT)=0.
-     call four2a(s,NFFT,1,-1,0)                !Real-to-complex FFT
-     do i=0,NH
-        cs2(i)=cs(i)*conjg(cpr2(i))            !Mult by complex FFT of pr2
-        cs(i)=cs(i)*conjg(cpr(i))              !Mult by complex FFT of pr
-     enddo
-     call four2a(cs,NFFT,1,1,-1)               !Complex-to-real inv-FFT
-     call four2a(cs2,NFFT,1,1,-1)              !Complex-to-real inv-FFT
-     
-     do lag=-27,27                             !Check for best JT65 sync
-        ccf(lag,ip)=s(lag+28)                  
-        if(abs(ccf(lag,ip)).gt.ccfbest) then
-           ccfbest=abs(ccf(lag,ip))
-           lagpk=lag
-           ipol1=ip
-           flipk=1.0
-           if(ccf(lag,ip).lt.0.0) flipk=-1.0
-        endif
-     enddo
-
-     do lag=-8,7                               !Check for best shorthand
-        ccf2=s2(lag+28)
-        if(ccf2.gt.ccfbest2) then
-           ccfbest2=ccf2
-           lagpk2=lag
-           ipol2=ip
-        endif
-     enddo
-     
-  enddo
-
-! Find rms level on baseline of "ccfblue", for normalization.
-  sum=0.
-  do lag=-26,26
-     if(abs(lag-lagpk).gt.1) sum=sum + ccf(lag,ipol1)
-  enddo
-  base=sum/50.0
-  sq=0.
-  do lag=-26,26
-     if(abs(lag-lagpk).gt.1) sq=sq + (ccf(lag,ipol1)-base)**2
-  enddo
-  rms=sqrt(sq/49.0)
-  sync1=ccfbest/rms - 4.0
-  dt1=2.5 + lagpk*(2048.0/11025.0)
-
-! Find base level for normalizing snr2.
-  do i=1,nhsym
-     tmp1(i)=ss(ipol2,i)
-  enddo
-  call pctile(tmp1,tmp2,nhsym,40,base)
-  snr2=0.398107*ccfbest2/base                !### empirical
-  syncshort=0.5*ccfbest2/rms - 4.0           !### better normalizer than rms?
-  dt2=2.5 + lagpk2*(2048.0/11025.0)
-
-  return
-end subroutine ccf65
diff --git a/libm65/chkhist.f b/libm65/chkhist.f
deleted file mode 100644
index 1b1a694b1..000000000
--- a/libm65/chkhist.f
+++ /dev/null
@@ -1,23 +0,0 @@
-      subroutine chkhist(mrsym,nmax,ipk)
-
-      integer mrsym(63)
-      integer hist(0:63)
-
-      do i=0,63
-         hist(i)=0
-      enddo
-      do j=1,63
-         i=mrsym(j)
-         hist(i)=hist(i)+1
-      enddo
-
-      nmax=0
-      do i=0,63
-         if(hist(i).gt.nmax) then
-            nmax=hist(i)
-            ipk=i+1
-         endif
-      enddo
-
-      return
-      end
diff --git a/libm65/chkmsg.f b/libm65/chkmsg.f
deleted file mode 100644
index 9dd7ef577..000000000
--- a/libm65/chkmsg.f
+++ /dev/null
@@ -1,32 +0,0 @@
-      subroutine chkmsg(message,cok,nspecial,flip)
-
-      character message*22,cok*3
-
-      nspecial=0
-      flip=1.0
-      cok="   "
-
-      do i=22,1,-1
-         if(message(i:i).ne.' ') go to 10
-      enddo
-      i=22
-
- 10   if(i.ge.11) then
-         if ((message(i-3:i).eq.' OOO') .or. 
-     +        (message(20:22).eq.' OO')) then
-            cok='OOO'
-            flip=-1.0
-            if(message(20:22).eq.' OO') then
-               message=message(1:19)
-            else
-               message=message(1:i-4)
-            endif
-         endif
-      endif
-
-      if(message(1:3).eq.'RO ')  nspecial=2
-      if(message(1:4).eq.'RRR ') nspecial=3
-      if(message(1:3).eq.'73 ')  nspecial=4
-
-      return
-      end
diff --git a/libm65/coord.f b/libm65/coord.f
deleted file mode 100644
index 3fb7c647d..000000000
--- a/libm65/coord.f
+++ /dev/null
@@ -1,41 +0,0 @@
-        SUBROUTINE COORD(A0,B0,AP,BP,A1,B1,A2,B2)
-
-C  Examples:
-C  1. From ha,dec to az,el:
-C       call coord(pi,pio2-lat,0.,lat,ha,dec,az,el)
-C  2. From az,el to ha,dec:
-C       call coord(pi,pio2-lat,0.,lat,az,el,ha,dec)
-C  3. From ra,dec to l,b
-C       call coord(4.635594495,-0.504691042,3.355395488,0.478220215,
-C         ra,dec,l,b)
-C  4. From l,b to ra,dec
-C       call coord(1.705981071d0,-1.050357016d0,2.146800277d0,
-C         0.478220215d0,l,b,ra,dec)
-C  5. From ra,dec to ecliptic latitude (eb) and longitude (el):
-C       call coord(0.d0,0.d0,-pio2,pio2-23.443*pi/180,ra,dec,el,eb)
-C  6. From ecliptic latitude (eb) and longitude (el) to ra,dec:
-C       call coord(0.d0,0.d0,-pio2,pio2-23.443*pi/180,el,eb,ra,dec)
-
-
-      SB0=sin(B0)
-      CB0=cos(B0)
-      SBP=sin(BP)
-      CBP=cos(BP)
-      SB1=sin(B1)
-      CB1=cos(B1)
-      SB2=SBP*SB1 + CBP*CB1*cos(AP-A1)
-      CB2=SQRT(1.e0-SB2**2)
-      B2=atan(SB2/CB2)
-      SAA=sin(AP-A1)*CB1/CB2
-      CAA=(SB1-SB2*SBP)/(CB2*CBP)
-      CBB=SB0/CBP
-      SBB=sin(AP-A0)*CB0
-      SA2=SAA*CBB-CAA*SBB
-      CA2=CAA*CBB+SAA*SBB
-      TA2O2=0.0 !Shut up compiler warnings. -db
-      IF(CA2.LE.0.e0) TA2O2=(1.e0-CA2)/SA2 
-      IF(CA2.GT.0.e0) TA2O2=SA2/(1.e0+CA2)
-      A2=2.e0*atan(TA2O2)
-      IF(A2.LT.0.e0) A2=A2+6.2831853
-      RETURN
-      END
diff --git a/libm65/dcoord.f b/libm65/dcoord.f
deleted file mode 100644
index b3d12d846..000000000
--- a/libm65/dcoord.f
+++ /dev/null
@@ -1,40 +0,0 @@
-        SUBROUTINE DCOORD(A0,B0,AP,BP,A1,B1,A2,B2)
-
-        implicit real*8 (a-h,o-z)
-C  Examples:
-C  1. From ha,dec to az,el:
-C       call coord(pi,pio2-lat,0.,lat,ha,dec,az,el)
-C  2. From az,el to ha,dec:
-C       call coord(pi,pio2-lat,0.,lat,az,el,ha,dec)
-C  3. From ra,dec to l,b
-C       call coord(4.635594495,-0.504691042,3.355395488,0.478220215,
-C         ra,dec,l,b)
-C  4. From l,b to ra,dec
-C       call coord(1.705981071d0,-1.050357016d0,2.146800277d0,
-C         0.478220215d0,l,b,ra,dec)
-C  5. From ecliptic latitude (eb) and longitude (el) to ra, dec:
-C       call coord(0.d0,0.d0,-pio2,pio2-23.443*pi/180,ra,dec,el,eb)
-
-      SB0=sin(B0)
-      CB0=cos(B0)
-      SBP=sin(BP)
-      CBP=cos(BP)
-      SB1=sin(B1)
-      CB1=cos(B1)
-      SB2=SBP*SB1 + CBP*CB1*cos(AP-A1)
-      CB2=SQRT(1.D0-SB2**2)
-      B2=atan(SB2/CB2)
-      SAA=sin(AP-A1)*CB1/CB2
-      CAA=(SB1-SB2*SBP)/(CB2*CBP)
-      CBB=SB0/CBP
-      SBB=sin(AP-A0)*CB0
-      SA2=SAA*CBB-CAA*SBB
-      CA2=CAA*CBB+SAA*SBB
-      TA2O2=0.0 !Shut up compiler warnings. -db
-      IF(CA2.LE.0.D0) TA2O2=(1.D0-CA2)/SA2 
-      IF(CA2.GT.0.D0) TA2O2=SA2/(1.D0+CA2)
-      A2=2.D0*atan(TA2O2)
-      IF(A2.LT.0.D0) A2=A2+6.2831853071795864D0
-
-      RETURN
-      END
diff --git a/libm65/decode1a.f b/libm65/decode1a.f
deleted file mode 100644
index 5102a7a7c..000000000
--- a/libm65/decode1a.f
+++ /dev/null
@@ -1,149 +0,0 @@
-      subroutine decode1a(dd,newdat,f0,nflip,mode65,nfsample,xpol,
-     +     mycall,hiscall,hisgrid,neme,ndepth,nqd,dphi,
-     +     nutc,nkhz,ndf,ipol,sync2,a,dt,pol,nkv,nhist,qual,decoded)
-
-!  Apply AFC corrections to a candidate JT65 signal, then decode it.
-
-      parameter (NMAX=60*96000)          !Samples per 60 s
-      real*4  dd(4,NMAX)                 !92 MB: raw data from Linrad timf2
-      complex cx(NMAX/64), cy(NMAX/64)   !Data at 1378.125 samples/s
-      complex c5x(NMAX/256),c5y(NMAX/256),c5tmp(NMAX/256)  !Data at 344.53125 Hz
-      complex c5a(512)
-      complex z
-      real s2(66,126)
-      real s3(64,63),sy(63)
-      real a(5)
-      logical first,xpol
-      character decoded*22
-      character mycall*12,hiscall*12,hisgrid*6
-      data first/.true./,jjjmin/1000/,jjjmax/-1000/
-      data nutc0/-999/,nkhz0/-999/
-      save
-
-!  Mix sync tone to baseband, low-pass filter, downsample to 1378.125 Hz
-      dt00=dt
-      call timer('filbig  ',0)
-      call filbig(dd,NMAX,f0,newdat,nfsample,xpol,cx,cy,n5)
-!  NB: cx, cy have sample rate 96000*77125/5376000 = 1378.125 Hz
-      call timer('filbig  ',1)
-      joff=0
-      sqa=0.
-      sqb=0.
-      do i=1,n5
-         sqa=sqa + real(cx(i))**2 + aimag(cx(i))**2
-         sqb=sqb + real(cy(i))**2 + aimag(cy(i))**2
-      enddo
-      sqa=sqa/n5
-      sqb=sqb/n5
-
-!  Find best DF, f1, f2, DT, and pol.  Start by downsampling to 344.53125 Hz
-      z=cmplx(cos(dphi),sin(dphi))
-      cy(:n5)=z*cy(:n5)                   !Adjust for cable length difference
-      call timer('fil6521 ',0)
-      call fil6521(cx,n5,c5x,n6)
-      call fil6521(cy,n5,c5y,n6)
-      call timer('fil6521 ',1)
-
-! Add some zeros at start of c5 arrays -- empirical fix for negative DT's
-! NB: might be better to add zeros to cx and cy, rather than here.
-! Q: is the DT search range big enough?
-
-      nadd=200
-      c5tmp(1:nadd)=0.
-      c5tmp(1+nadd:n6+nadd)=c5x(1:n6)
-      c5x(1:n6+nadd)=c5tmp(1:n6+nadd)
-      c5tmp(1+nadd:n6+nadd)=c5y(1:n6)
-      c5y(1:n6+nadd)=c5tmp(1:n6+nadd)
-      n6=n6+nadd
-
-      fsample=1378.125/4.
-      a(5)=dt00
-      i0=nint((a(5)+0.5)*fsample) - 2 + 200
-      if(i0.lt.1) then
-         write(13,*) 'i0 too small in decode1a:',i0,f0
-         flush(13)
-         i0=1
-      endif
-      nz=n6+1-i0
-
-! We're looking only at sync tone here... so why not downsample by another
-! factor of 1/8, say?  Should be a significant execution speed-up.
-      call timer('afc65b  ',0)
-! Best fit for DF, f1, f2, pol
-      call afc65b(c5x(i0),c5y(i0),nz,fsample,nflip,ipol,xpol,a,
-     +     ccfbest,dtbest)
-      call timer('afc65b  ',1)
-
-      pol=a(4)/57.2957795
-      aa=cos(pol)
-      bb=sin(pol)
-      sq0=aa*aa*sqa + bb*bb*sqb
-      sync2=3.7*ccfbest/sq0
-
-!  Apply AFC corrections to the time-domain signal
-!  Now we are back to using the 1378.125 Hz sample rate, enough to 
-!  accommodate the full JT65C bandwidth.
-
-      call timer('twkfreq ',0)
-      call twkfreq(cx,cy,n5,a)
-      call timer('twkfreq ',1)
-
-!  Compute spectrum at best polarization for each half symbol.
-!  Adding or subtracting a small number (e.g., 5) to j may make it decode.\
-!  NB: might want to try computing full-symbol spectra (nfft=512, even for
-!  submodes B and C).
-
-      nsym=126
-      nfft=512/mode65
-      j=(dt00+dtbest+2.685)*1378.125 + joff
-      if(j.lt.0) j=0
-
-      call timer('sh_ffts ',0)
-
-! Perhaps should try full-symbol-length FFTs even in B, C sub-modes?
-! (Tried this, found no significant difference in decodes.)
-
-      do k=1,nsym
-         do n=1,mode65
-            do i=1,nfft
-               j=j+1
-               c5a(i)=aa*cx(j) + bb*cy(j)
-            enddo
-            call four2a(c5a,nfft,1,1,1)
-            if(n.eq.1) then
-               do i=1,66
-                  s2(i,k)=real(c5a(i))**2 + aimag(c5a(i))**2
-               enddo
-            else
-               do i=1,66
-                  s2(i,k)=s2(i,k) + real(c5a(i))**2 + aimag(c5a(i))**2
-               enddo
-            endif
-         enddo
-      enddo
-
-      call timer('sh_ffts ',1)
-
-      flip=nflip
-      call timer('dec65b  ',0)
-      call decode65b(s2,flip,mycall,hiscall,hisgrid,mode65,neme,ndepth,
-     +    nqd,nkv,nhist,qual,decoded,s3,sy)
-      dt=dt00 + dtbest
-      call timer('dec65b  ',1)
-
-      if(nqd.eq.1 .and. nkv.eq.0) then
-         if(nutc.ne.nutc0) syncbest=0.
-         if(sync2.gt.syncbest) then
-            if(nutc.eq.nutc0) nsave=nsave-1
-            if(nkhz.ne.nkhz0) nsave=0
-            nkhz0=nkhz
-            nsave=min(32,nsave+1)
-            npol=nint(57.296*pol)
-            call s3avg(nsave,mode65,nutc,ndf,dt+0.8,npol,s3,nkv,decoded)
-            syncbest=sync2
-         endif
-         nutc0=nutc
-      endif
-
-      return
-      end
diff --git a/libm65/decode65b.f b/libm65/decode65b.f
deleted file mode 100644
index b94c05534..000000000
--- a/libm65/decode65b.f
+++ /dev/null
@@ -1,49 +0,0 @@
-      subroutine decode65b(s2,flip,mycall,hiscall,hisgrid,mode65,
-     +  neme,ndepth,nqd,nkv,nhist,qual,decoded,s3,sy)
-
-      real s2(66,126)
-      real s3(64,63),sy(63)
-      logical first,ltext
-      character decoded*22,deepmsg*22
-      character mycall*12,hiscall*12,hisgrid*6
-      common/prcom/pr(126),mdat(126),mref(126,2),mdat2(126),mref2(126,2)
-      data first/.true./
-      save
-
-      if(first) call setup65
-      first=.false.
-
-      do j=1,63
-         k=mdat(j)                       !Points to data symbol
-         if(flip.lt.0.0) k=mdat2(j)
-         do i=1,64
-            s3(i,j)=s2(i+2,k)
-         enddo
-         k=mdat2(j)                       !Points to data symbol
-         if(flip.lt.0.0) k=mdat(j)
-         sy(j)=s2(1,k)
-      enddo
-
-      nadd=mode65
-      call extract(s3,nadd,ncount,nhist,decoded,ltext)     !Extract the message
-C  Suppress "birdie messages" and other garbage decodes:
-      if(decoded(1:7).eq.'000AAA ') ncount=-1
-      if(decoded(1:7).eq.'0L6MWK ') ncount=-1
-      if(flip.lt.0.0 .and. ltext) ncount=-1
-      nkv=1
-      if(ncount.lt.0) then 
-         nkv=0
-         decoded='                      '
-      endif
-
-      qual=0.
-      if(ndepth.ge.1 .and. (nqd.eq.1 .or. flip.eq.1.0)) then
-         call deep65(s3,mode65,neme,flip,mycall,hiscall,
-     +       hisgrid,deepmsg,qual)
-         if(nqd.ne.1 .and. qual.lt.10.0) qual=0.0
-         if(ndepth.lt.2 .and. qual.lt.6.0) qual=0.0
-      endif
-      if(nkv.eq.0 .and. qual.ge.1.0) decoded=deepmsg
-
-      return
-      end
diff --git a/libm65/decodemsk.f90 b/libm65/decodemsk.f90
deleted file mode 100644
index 9a08edd82..000000000
--- a/libm65/decodemsk.f90
+++ /dev/null
@@ -1,39 +0,0 @@
-subroutine decodemsk(cdat,npts,cw,i1,nchar,s2,msg)
-
-! DF snd sync have been established, now decode the message
-
-  complex cdat(npts)
-  complex cw(168,0:63)                  !Complex waveforms for codewords
-  real s2(0:63,400)
-  character msg*400
-  complex z,zmax
-  character cc*64
-!                    1         2         3         4         5         6
-!          0123456789012345678901234567890123456789012345678901234567890123
-  data cc/'0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ./?-                 _     @'/
-
-  msg=' '
-  do j=1,nchar                         !Find best match for each character
-     ia=i1 + (j-1)*168
-     smax=0.
-     do k=0,40
-        kk=k
-        if(k.eq.40) kk=44
-        z=0.
-        do i=1,168
-           z=z + cdat(ia+i)*conjg(cw(i,kk))
-        enddo
-        ss=abs(z)
-        s2(k,j)=ss
-        if(ss.gt.smax) then
-           smax=ss
-           zmax=z
-           kpk=kk
-        endif
-     enddo
-     msg(j:j)=cc(kpk+1:kpk+1)
-     if(kpk.eq.44) msg(j:j)=' '
-  enddo
-
-  return
-end subroutine decodemsk
diff --git a/libm65/deep65.f90 b/libm65/deep65.f90
deleted file mode 100644
index 1ddfe844b..000000000
--- a/libm65/deep65.f90
+++ /dev/null
@@ -1,169 +0,0 @@
-subroutine deep65(s3,mode65,neme,flip,mycall,hiscall,hisgrid,decoded,qual)
-
-  parameter (MAXCALLS=7000,MAXRPT=63)
-  real s3(64,63)
-  character callsign*12,grid*4,message*22,hisgrid*6,c*1,ceme*3
-  character*12 mycall,hiscall
-  character*22 decoded
-  character*22 testmsg(2*MAXCALLS + 2 + MAXRPT)
-  character*15 callgrid(MAXCALLS)
-  character*180 line
-  character*4 rpt(MAXRPT)
-  integer ncode(63,2*MAXCALLS + 2 + MAXRPT)
-  real pp(2*MAXCALLS + 2 + MAXRPT)
-  common/mrscom/ mrs(63),mrs2(63)
-  common/c3com/ mcall3a
-  data rpt/'-01','-02','-03','-04','-05',          &
-           '-06','-07','-08','-09','-10',          &
-           '-11','-12','-13','-14','-15',          &
-           '-16','-17','-18','-19','-20',          &
-           '-21','-22','-23','-24','-25',          &
-           '-26','-27','-28','-29','-30',          &
-           'R-01','R-02','R-03','R-04','R-05',     &
-           'R-06','R-07','R-08','R-09','R-10',     &
-           'R-11','R-12','R-13','R-14','R-15',     &
-           'R-16','R-17','R-18','R-19','R-20',     &
-           'R-21','R-22','R-23','R-24','R-25',     &
-           'R-26','R-27','R-28','R-29','R-30',     &
-           'RO','RRR','73'/
-  save
-
-  if(mcall3a.eq.0) go to 30
-
-  call timer('deep65a ',0)
-  mcall3a=0
-  rewind 23
-  k=0
-  icall=0
-  do n=1,MAXCALLS
-     if(n.eq.1) then
-        callsign=hiscall
-        do i=4,12
-           if(ichar(callsign(i:i)).eq.0) callsign(i:i)=' '
-        enddo
-        grid=hisgrid(1:4)
-        if(ichar(grid(3:3)).eq.0) grid(3:3)=' '
-        if(ichar(grid(4:4)).eq.0) grid(4:4)=' '
-     else
-        read(23,1002,end=20) line
-1002    format (A80)
-        if(line(1:4).eq.'ZZZZ') go to 20
-        if(line(1:2).eq.'//') go to 10
-        i1=index(line,',')
-        if(i1.lt.4) go to 10
-        i2=index(line(i1+1:),',')
-        if(i2.lt.5) go to 10
-        i2=i2+i1
-        i3=index(line(i2+1:),',')
-        if(i3.lt.1) i3=index(line(i2+1:),' ')
-        i3=i2+i3
-        callsign=line(1:i1-1)
-        grid=line(i1+1:i2-1)
-        ceme=line(i2+1:i3-1)
-        if(neme.eq.1 .and. ceme.ne.'EME') go to 10
-     endif
-
-     icall=icall+1
-     j1=index(mycall,' ') - 1
-     if(j1.le.-1) j1=12
-     if(j1.lt.3) j1=6
-     j2=index(callsign,' ') - 1
-     if(j2.le.-1) j2=12
-     if(j2.lt.3) j2=6
-     j3=index(mycall,'/')                 ! j3>0 means compound mycall
-     j4=index(callsign,'/')               ! j4>0 means compound hiscall
-     callgrid(icall)=callsign(1:j2)
-
-     mz=1
-! Allow MyCall + HisCall + rpt (?)
-     if(n.eq.1 .and. j3.lt.1 .and. j4.lt.1 .and.                       &
-          flip.gt.0.0 .and. callsign(1:6).ne.'      ') mz=MAXRPT+1
-     do m=1,mz
-        if(m.gt.1) grid=rpt(m-1)
-        if(j3.lt.1 .and.j4.lt.1) callgrid(icall)=callsign(1:j2)//' '//grid
-        message=mycall(1:j1)//' '//callgrid(icall)
-        k=k+1
-        testmsg(k)=message
-        call encode65(message,ncode(1,k))
-        
-! Insert CQ message
-        if(j4.lt.1) callgrid(icall)=callsign(1:j2)//' '//grid
-        message='CQ '//callgrid(icall)
-        k=k+1
-        testmsg(k)=message
-        call encode65(message,ncode(1,k))
-     enddo
-10   continue
-  enddo
-
-20 continue
-  ntot=k
-  call timer('deep65a ',1)
-
-30 continue
-  call timer('deep65b ',0)
-  ref0=0.
-  do j=1,63
-     ref0=ref0 + s3(mrs(j),j)
-  enddo
-
-  p1=-1.e30
-  p2=-1.e30
-  do k=1,ntot
-     pp(k)=0.
-! Test all messages if flip=+1; skip the CQ messages if flip=-1.
-     if(flip.gt.0.0 .or. testmsg(k)(1:3).ne.'CQ ') then
-        sum=0.
-        ref=ref0
-        do j=1,63
-           i=ncode(j,k)+1
-           sum=sum + s3(i,j)
-           if(i.eq.mrs(j)) ref=ref - s3(i,j) + s3(mrs2(j),j)
-        enddo
-        p=sum/ref
-        pp(k)=p
-        if(p.gt.p1) then
-           p1=p
-           ip1=k
-        endif
-     endif
-  enddo
-
-  do i=1,ntot
-     if(pp(i).gt.p2 .and. pp(i).ne.p1) p2=pp(i)
-  enddo
-
-! ### DO NOT REMOVE ### 
-  rewind 77
-  write(77,*) p1,p2
-! ### Works OK without it (in both Windows and Linux) if compiled 
-! ### without optimization.  However, in Windows this is a colossal 
-! ### pain because of the way F2PY wants to run the compile step.
-
-  if(mode65.eq.1) bias=max(1.12*p2,0.335)
-  if(mode65.eq.2) bias=max(1.08*p2,0.405)
-  if(mode65.ge.4) bias=max(1.04*p2,0.505)
-
-  if(p2.eq.p1 .and. p1.ne.-1.e30) stop 'Error in deep65'
-  qual=100.0*(p1-bias)
-
-  decoded='                      '
-  c=' '
-
-  if(qual.gt.1.0) then
-     if(qual.lt.6.0) c='?'
-     decoded=testmsg(ip1)
-  else
-     qual=0.
-  endif
-  decoded(22:22)=c
-
-! Make sure everything is upper case.
-  do i=1,22
-     if(decoded(i:i).ge.'a' .and. decoded(i:i).le.'z')                &
-          decoded(i:i)=char(ichar(decoded(i:i))-32)
-  enddo
-  call timer('deep65b ',1)
-
-  return
-end subroutine deep65
diff --git a/libm65/demod64a.f b/libm65/demod64a.f
deleted file mode 100644
index b5fbd6087..000000000
--- a/libm65/demod64a.f
+++ /dev/null
@@ -1,73 +0,0 @@
-      subroutine demod64a(s3,nadd,mrsym,mrprob,
-     +  mr2sym,mr2prob,ntest,nlow)
-
-C  Demodulate the 64-bin spectra for each of 63 symbols in a frame.
-
-C  Parameters
-C     nadd     number of spectra already summed
-C     mrsym    most reliable symbol value
-C     mr2sym   second most likely symbol value
-C     mrprob   probability that mrsym was the transmitted value
-C     mr2prob  probability that mr2sym was the transmitted value
-
-      implicit real*8 (a-h,o-z)
-      real*4 s3(64,63)
-      real*8 fs(64)
-      integer mrsym(63),mrprob(63),mr2sym(63),mr2prob(63)
-      common/mrscom/ mrs(63),mrs2(63)
-
-      afac=1.1 * float(nadd)**0.64
-      scale=255.999
-
-C  Compute average spectral value
-      sum=0.
-      do j=1,63
-         do i=1,64
-            sum=sum+s3(i,j)
-         enddo
-      enddo
-      ave=sum/(64.*63.)
-      i1=1                                      !Silence warning
-      i2=1
-
-C  Compute probabilities for most reliable symbol values
-      do j=1,63
-         s1=-1.e30
-         fsum=0.
-         do i=1,64
-            x=min(afac*s3(i,j)/ave,50.d0)
-            fs(i)=exp(x)
-            fsum=fsum+fs(i)
-            if(s3(i,j).gt.s1) then
-               s1=s3(i,j)
-               i1=i                              !Most reliable
-            endif
-         enddo
-
-         s2=-1.e30
-         do i=1,64
-            if(i.ne.i1 .and. s3(i,j).gt.s2) then
-               s2=s3(i,j)
-               i2=i                              !Second most reliable
-            endif
-         enddo
-         p1=fs(i1)/fsum                          !Normalized probabilities
-         p2=fs(i2)/fsum
-         mrsym(j)=i1-1
-         mr2sym(j)=i2-1
-         mrprob(j)=scale*p1
-         mr2prob(j)=scale*p2
-         mrs(j)=i1
-         mrs2(j)=i2
-      enddo
-
-      sum=0.
-      nlow=0
-      do j=1,63
-         sum=sum+mrprob(j)
-         if(mrprob(j).le.5) nlow=nlow+1
-      enddo
-      ntest=sum/63
-
-      return
-      end
diff --git a/libm65/display.f90 b/libm65/display.f90
deleted file mode 100644
index 22f4a6933..000000000
--- a/libm65/display.f90
+++ /dev/null
@@ -1,169 +0,0 @@
-subroutine display(nkeep,ftol)
-
-  parameter (MAXLINES=400,MX=400)
-  integer indx(MAXLINES),indx2(MX)
-  character*81 line(MAXLINES),line2(MX),line3(MAXLINES)
-  character out*50,cfreq0*3,cqlive*52
-  character*6 callsign,callsign0
-  character*12 freqcall(100)
-  real freqkHz(MAXLINES)
-  integer utc(MAXLINES),utc2(MX),utcz
-  real*8 f0
-
-  rewind 26
-
-  do i=1,MAXLINES
-     read(26,1010,end=10) line(i)
-1010 format(a80)
-     read(line(i),1020) f0,ndf,nh,nm
-1020 format(f8.3,i5,25x,i3,i2)
-     utc(i)=60*nh + nm
-     freqkHz(i)=1000.d0*(f0-144.d0) + 0.001d0*ndf
-  enddo
-
-10 nz=i-1
-  utcz=utc(nz)
-  nz=nz-1
-  if(nz.lt.1) go to 999
-  nquad=max(nkeep/4,3)
-  do i=1,nz
-     nage=utcz-utc(i)
-     if(nage.lt.0) nage=nage+1440
-     iage=nage/nquad
-     write(line(i)(80:81),1021) iage
-1021 format(i2)
-  enddo
-
-  nage=utcz-utc(1)
-  if(nage.lt.0) nage=nage+1440
-  if(nage.gt.nkeep) then
-     do i=1,nz
-        nage=utcz-utc(i)
-        if(nage.lt.0) nage=nage+1440
-        if(nage.le.nkeep) go to 20
-     enddo
-20   i0=i
-     nz=nz-i0+1
-     rewind 26
-     if(nz.lt.1) go to 999
-     do i=1,nz
-        j=i+i0-1
-        line(i)=line(j)
-        utc(i)=utc(j)
-        freqkHz(i)=freqkHz(j)
-        write(26,1010) line(i)
-     enddo
-  endif
-
-  call flush(26)
-  call indexx(nz,freqkHz,indx)
-
-  nstart=1
-  k3=0
-  k=1
-  m=indx(1)
-  if(m.lt.1 .or. m.gt.MAXLINES) then
-     print*,'Error in display.f90: ',nz,m
-     m=1
-  endif
-  line2(1)=line(m)
-  utc2(1)=utc(m)
-  do i=2,nz
-     j0=indx(i-1)
-     j=indx(i)
-     if(freqkHz(j)-freqkHz(j0).gt.2.0*ftol) then
-        if(nstart.eq.0) then
-           k=k+1
-           line2(k)=""
-           utc2(k)=-1
-        endif
-        kz=k
-        if(nstart.eq.1) then
-           call indexx(kz,utc2,indx2)
-           k3=0
-           do k=1,kz
-              k3=min(k3+1,400)
-              line3(k3)=line2(indx2(k))
-           enddo
-           nstart=0
-        else
-           call indexx(kz,utc2,indx2)
-           do k=1,kz
-              k3=min(k3+1,400)
-              line3(k3)=line2(indx2(k))
-           enddo
-        endif
-        k=0
-     endif
-     if(i.eq.nz) then
-        k=k+1
-        line2(k)=""
-        utc2(k)=-1
-     endif
-     k=k+1
-     line2(k)=line(j)
-     utc2(k)=utc(j)
-     j0=j
-  enddo
-  kz=k
-  call indexx(kz,utc2,indx2)
-  do k=1,kz
-     k3=min(k3+1,400)
-     line3(k3)=line2(indx2(k))
-  enddo
-
-  rewind 19
-  rewind 20
-  cfreq0='   '
-  nc=0
-  callsign0='      '
-  do k=1,k3
-     out=line3(k)(6:13)//line3(k)(28:31)//line3(k)(39:43)//       &
-          line3(k)(35:38)//line3(k)(44:67)//line3(k)(77:81)
-     if(out(1:3).ne.'   ') then
-        cfreq0=out(1:3)
-        if(iw.lt.MAXLINES-1) iw=iw+1
-        cqlive=line3(k)(6:13)//line3(k)(28:31)//line3(k)(39:43)//       &
-             line3(k)(23:27)//line3(k)(35:38)//line3(k)(44:67)//        &
-             line3(k)(80:81)
-        if(index(cqlive,' CQ ').gt.0 .or. index(cqlive,' QRZ ').gt.0 .or.   &
-           index(cqlive,' QRT ').gt.0 .or. index(cqlive,' CQV ').gt.0 .or.  &
-           index(cqlive,' CQH ').gt.0) write(19,1029) cqlive
-1029    format(a52)
-        write(*,1030) out
-1030    format('@',a50)
-        i1=index(out(24:),' ')
-        callsign=out(i1+24:)
-        i2=index(callsign,' ')
-        if(i2.gt.1) callsign(i2:)='      '
-        if(callsign.ne.'      ' .and. callsign.ne.callsign0) then
-           len=i2-1
-           if(len.lt.0) len=6
-           if(len.ge.4) then                        !Omit short "callsigns"
-              nc=nc+1
-              freqcall(nc)=cfreq0//' '//callsign//line3(k)(80:81)
-              callsign0=callsign
-           endif
-        endif
-        if(callsign.ne.'      ' .and. callsign.eq.callsign0) then
-           freqcall(nc)=cfreq0//' '//callsign//line3(k)(80:81)
-        endif
-     endif
-  enddo
-  flush(19)
-  nc=nc+1
-  freqcall(nc)='            '
-  nc=nc+1
-  freqcall(nc)='            '
-  freqcall(nc+1)='            '
-  freqcall(nc+2)='            '
-  iz=(nc+2)/3
-
-  do i=1,nc
-  write(*,1042) freqcall(i)
-1042 format('&',a12)
-  enddo
-
-999  continue
-  return
-end subroutine display
diff --git a/libm65/dot.f b/libm65/dot.f
deleted file mode 100644
index 8e2d826bf..000000000
--- a/libm65/dot.f
+++ /dev/null
@@ -1,11 +0,0 @@
-      real*8 function dot(x,y)
-
-      real*8 x(3),y(3)
-
-      dot=0.d0
-      do i=1,3
-         dot=dot+x(i)*y(i)
-      enddo
-
-      return
-      end
diff --git a/libm65/dpol.f90 b/libm65/dpol.f90
deleted file mode 100644
index 152d69349..000000000
--- a/libm65/dpol.f90
+++ /dev/null
@@ -1,41 +0,0 @@
-real function dpol(mygrid,hisgrid)
-
-! Compute spatial polartzation offset in degrees for the present 
-! time, between two specified grid locators.
-
-  character*6 MyGrid,HisGrid
-  real lat,lon,LST
-  character cdate*8,ctime2*10,czone*5,fnamedate*6
-  integer  it(8)
-  data rad/57.2957795/
-
-  call date_and_time(cdate,ctime2,czone,it)
-  nyear=it(1)
-  month=it(2)
-  nday=it(3)
-  nh=it(5)-it(4)/60
-  nm=it(6)
-  ns=it(7)
-  uth=nh + nm/60.0 + ns/3600.0
-
-  call grid2deg(MyGrid,lon,lat)
-  call MoonDop(nyear,month,nday,uth,-lon,lat,RAMoon,DecMoon,         &
-       LST,HA,AzMoon,ElMoon,vr,dist)
-  xx=sin(lat/rad)*cos(ElMoon/rad) - cos(lat/rad)*                    &
-       cos(AzMoon/rad)*sin(ElMoon/rad)
-  yy=cos(lat/rad)*sin(AzMoon/rad)
-  poloffset1=rad*atan2(yy,xx)
-
-  call grid2deg(hisGrid,lon,lat)
-  call MoonDop(nyear,month,nday,uth,-lon,lat,RAMoon,DecMoon,         &
-       LST,HA,AzMoon,ElMoon,vr,dist)
-  xx=sin(lat/rad)*cos(ElMoon/rad) - cos(lat/rad)*                    &
-       cos(AzMoon/rad)*sin(ElMoon/rad)
-  yy=cos(lat/rad)*sin(AzMoon/rad)
-  poloffset2=rad*atan2(yy,xx)
-
-  dpol=mod(poloffset2-poloffset1+720.0,180.0)
-  if(dpol.gt.90.0) dpol=dpol-180.0
-
-  return
-end function dpol
diff --git a/libm65/encode65.f b/libm65/encode65.f
deleted file mode 100644
index 670c2e583..000000000
--- a/libm65/encode65.f
+++ /dev/null
@@ -1,13 +0,0 @@
-      subroutine encode65(message,sent)
-
-      character message*22
-      integer dgen(12)
-      integer sent(63)
-
-      call packmsg(message,dgen)
-      call rs_encode(dgen,sent)
-      call interleave63(sent,1)
-      call graycode(sent,63,1)
-
-      return
-      end
diff --git a/libm65/extract.F b/libm65/extract.F
deleted file mode 100644
index 04313c241..000000000
--- a/libm65/extract.F
+++ /dev/null
@@ -1,109 +0,0 @@
-      subroutine extract(s3,nadd,ncount,nhist,decoded,ltext)
-
-      real s3(64,63)
-      real tmp(4032)
-      character decoded*22
-      integer era(51),dat4(12),indx(64)
-      integer mrsym(63),mr2sym(63),mrprob(63),mr2prob(63)
-      logical first,ltext
-      data first/.true./,nsec1/0/
-      save
-
-      nfail=0
- 1    continue
-!      call timer('demod64a',0)
-      call demod64a(s3,nadd,mrsym,mrprob,mr2sym,mr2prob,ntest,nlow)
-!      call timer('demod64a',1)
-      if(ntest.lt.50 .or. nlow.gt.20) then
-         ncount=-999                         !Flag bad data
-         go to 900
-      endif
-      call chkhist(mrsym,nhist,ipk)
-
-      if(nhist.ge.20) then
-         nfail=nfail+1
-         call pctile(s3,tmp,4032,50,base)     ! ### or, use ave from demod64a
-         do j=1,63
-            s3(ipk,j)=base
-         enddo
-         if(nfail.gt.30) then
-            decoded='                      '
-            ncount=-1
-            go to 900
-         endif
-         go to 1
-      endif
-
-      call graycode(mrsym,63,-1)
-      call interleave63(mrsym,-1)
-      call interleave63(mrprob,-1)
-
-      ndec=1
-      nemax=30                              !Was 200 (30)
-      maxe=8
-      xlambda=13.0                          !Was 12
-
-      if(ndec.eq.1) then
-         call graycode(mr2sym,63,-1)
-         call interleave63(mr2sym,-1)
-         call interleave63(mr2prob,-1)
-
-         nsec1=nsec1+1
-         write(22,rec=1) nsec1,xlambda,maxe,200,
-     +        mrsym,mrprob,mr2sym,mr2prob
-         call flush(22)
-!         call timer('kvasd   ',0)
-#ifdef UNIX
-         iret=system('./kvasd -q > dev_null')
-#else
-         iret=system('kvasd -q > dev_null')
-#endif
-!         call timer('kvasd   ',1)
-         if(iret.ne.0) then
-            if(first) write(*,1000) iret
- 1000       format('Error in KV decoder, or no KV decoder present.'/
-     +         'Return code:',i8,'.  Will use BM algorithm.')
-            ndec=0
-            first=.false.
-            go to 20
-         endif
-
-         read(22,rec=2) nsec2,ncount,dat4
-         j=nsec2                !Silence compiler warning
-         decoded='                      '
-         ltext=.false.
-         if(ncount.ge.0) then
-            call unpackmsg(dat4,decoded)     !Unpack the user message
-            if(iand(dat4(10),8).ne.0) ltext=.true.
-            do i=2,12
-               if(dat4(i).ne.dat4(1)) go to 20
-            enddo
-            write(13,*) 'Bad decode?',nhist,nfail,ipk,
-     +      ' ',dat4,decoded
-            ncount=-1           !Suppress supposedly bogus decodes
-            decoded='                      '
-         endif
-      endif
- 20   if(ndec.eq.0) then
-         call indexx(63,mrprob,indx)
-         do i=1,nemax
-            j=indx(i)
-            if(mrprob(j).gt.120) then
-               ne2=i-1
-               go to 2
-            endif
-            era(i)=j-1
-         enddo
-         ne2=nemax
- 2       decoded='                      '
-         do nerase=0,ne2,2
-            call rs_decode(mrsym,era,nerase,dat4,ncount)
-            if(ncount.ge.0) then
-               call unpackmsg(dat4,decoded)
-               go to 900
-            endif
-         enddo
-      endif
-
- 900  return
-      end
diff --git a/libm65/fchisq.f b/libm65/fchisq.f
deleted file mode 100644
index f4c17d79c..000000000
--- a/libm65/fchisq.f
+++ /dev/null
@@ -1,76 +0,0 @@
-      real function fchisq(cx,cy,npts,fsample,nflip,a,ccfmax,dtmax)
-
-      parameter (NMAX=60*96000)          !Samples per 60 s
-      complex cx(npts),cy(npts)
-      real a(5)
-      complex w,wstep,za,zb,z
-      real ss(2600)
-      complex csx(0:NMAX/64),csy(0:NMAX/64)
-      data twopi/6.283185307/a1,a2,a3/99.,99.,99./
-      save
-
-      call timer('fchisq  ',0)
-      baud=11025.0/4096.0
-      if(a(1).ne.a1 .or. a(2).ne.a2 .or. a(3).ne.a3) then
-         a1=a(1)
-         a2=a(2)
-         a3=a(3)
-
-C  Mix and integrate the complex X and Y signals
-         csx(0)=0.
-         csy(0)=0.
-         w=1.0
-         x0=0.5*(npts+1)
-         s=2.0/npts
-         do i=1,npts
-            x=s*(i-x0)
-            if(mod(i,100).eq.1) then
-               p2=1.5*x*x - 0.5
-!               p3=2.5*(x**3) - 1.5*x
-!               p4=4.375*(x**4) - 3.75*(x**2) + 0.375
-               dphi=(a(1) + x*a(2) + p2*a(3)) * (twopi/fsample)
-               wstep=cmplx(cos(dphi),sin(dphi))
-            endif
-            w=w*wstep
-            csx(i)=csx(i-1) + w*cx(i)
-            csy(i)=csy(i-1) + w*cy(i)
-         enddo
-      endif
-
-C  Compute 1/2-symbol powers at 1/16-symbol steps.
-      fac=1.e-4
-      pol=a(4)/57.2957795
-      aa=cos(pol)
-      bb=sin(pol)
-      nsps=nint(fsample/baud)                  !Samples per symbol
-      nsph=nsps/2                              !Samples per half-symbol
-
-      ndiv=16                                  !Output ss() steps per symbol
-      nout=ndiv*npts/nsps
-      dtstep=1.0/(ndiv*baud)                   !Time per output step
-
-      do i=1,nout
-         j=i*nsps/ndiv
-         k=j-nsph
-         ss(i)=0.
-         if(k.ge.1) then
-            za=csx(j)-csx(k)
-            zb=csy(j)-csy(k)
-            z=aa*za + bb*zb
-            ss(i)=fac*(real(z)**2 + aimag(z)**2)
-         endif
-      enddo
-
-      ccfmax=0.
-      call timer('ccf2    ',0)
-      call ccf2(ss,nout,nflip,ccf,lagpk)
-      call timer('ccf2    ',1)
-      if(ccf.gt.ccfmax) then
-         ccfmax=ccf
-         dtmax=lagpk*dtstep
-      endif
-      fchisq=-ccfmax
-
-      call timer('fchisq  ',1)
-      return
-      end
diff --git a/libm65/fil6521.f b/libm65/fil6521.f
deleted file mode 100644
index 32feb8f44..000000000
--- a/libm65/fil6521.f
+++ /dev/null
@@ -1,44 +0,0 @@
-      subroutine fil6521(c1,n1,c2,n2)
-
-C  FIR lowpass filter designed using ScopeFIR
-
-C                   Pass #1   Pass #2  
-C-----------------------------------------------
-C  fsample    (Hz)  1378.125   Input sample rate
-C  Ntaps            21         Number of filter taps
-C  fc         (Hz)  40         Cutoff frequency
-C  fstop      (Hz)  172.266    Lower limit of stopband
-C  Ripple     (dB)  0.1        Ripple in passband
-C  Stop Atten (dB)  38         Stopband attenuation
-C  fout       (Hz)  344.531    Output sample rate
-
-      parameter (NTAPS=21)
-      parameter (NH=NTAPS/2)
-      parameter (NDOWN=4)                !Downsample ratio = 1/4
-      complex c1(n1)
-      complex c2(n1/NDOWN)
-
-C  Filter coefficients:
-      real a(-NH:NH)
-      data a/
-     +  -0.011958606980,-0.013888627387,-0.015601306443,-0.010602249570,
-     +   0.003804023436, 0.028320058273, 0.060903935217, 0.096841904411,
-     +   0.129639871228, 0.152644580853, 0.160917511283, 0.152644580853,
-     +   0.129639871228, 0.096841904411, 0.060903935217, 0.028320058273,
-     +   0.003804023436,-0.010602249570,-0.015601306443,-0.013888627387,
-     +  -0.011958606980/
-
-      n2=(n1-NTAPS+NDOWN)/NDOWN
-      k0=NH-NDOWN+1
-
-C  Loop over all output samples
-      do i=1,n2
-         c2(i)=0.
-         k=k0 + NDOWN*i
-         do j=-NH,NH
-            c2(i)=c2(i) + c1(j+k)*a(j)
-         enddo
-      enddo
-
-      return
-      end
diff --git a/libm65/filbig.f b/libm65/filbig.f
deleted file mode 100644
index 60206ddce..000000000
--- a/libm65/filbig.f
+++ /dev/null
@@ -1,134 +0,0 @@
-      subroutine filbig(dd,nmax,f0,newdat,nfsample,xpol,c4a,c4b,n4)
-
-C  Filter and downsample complex data stored in array dd(4,nmax).  
-C  Output is downsampled from 96000 Hz to 1375.125 Hz.
-
-      parameter (MAXFFT1=5376000,MAXFFT2=77175)
-      real*4  dd(4,nmax)                         !Input data
-      complex ca(MAXFFT1),cb(MAXFFT1)            !FFTs of input
-      complex c4a(MAXFFT2),c4b(MAXFFT2)          !Output data
-      real*8 df
-      real halfpulse(8)                 !Impulse response of filter (one sided)
-      complex cfilt(MAXFFT2)                     !Filter (complex; imag = 0)
-      real rfilt(MAXFFT2)                        !Filter (real)
-      integer*8 plan1,plan2,plan3,plan4,plan5
-      logical first,xpol
-      include 'fftw3.f'
-      equivalence (rfilt,cfilt)
-      data first/.true./,npatience/1/
-      data halfpulse/114.97547150,36.57879257,-20.93789101,
-     +  5.89886379,1.59355187,-2.49138308,0.60910773,-0.04248129/
-      save
-
-      nfft1=MAXFFT1
-      nfft2=MAXFFT2
-      if(nfsample.eq.95238) then
-         nfft1=5120000
-         nfft2=74088
-      endif
-      if(nmax.lt.0) go to 900
-      if(first) then
-         nflags=FFTW_ESTIMATE
-         if(npatience.eq.1) nflags=FFTW_ESTIMATE_PATIENT
-         if(npatience.eq.2) nflags=FFTW_MEASURE
-         if(npatience.eq.3) nflags=FFTW_PATIENT
-         if(npatience.eq.4) nflags=FFTW_EXHAUSTIVE
-C  Plan the FFTs just once
-         call timer('FFTplans ',0)
-         call sfftw_plan_dft_1d(plan1,nfft1,ca,ca,
-     +        FFTW_BACKWARD,nflags)
-         call sfftw_plan_dft_1d(plan2,nfft1,cb,cb,
-     +        FFTW_BACKWARD,nflags)
-         call sfftw_plan_dft_1d(plan3,nfft2,c4a,c4a,
-     +        FFTW_FORWARD,nflags)
-         call sfftw_plan_dft_1d(plan4,nfft2,c4b,c4b,
-     +        FFTW_FORWARD,nflags)
-         call sfftw_plan_dft_1d(plan5,nfft2,cfilt,cfilt,
-     +        FFTW_BACKWARD,nflags)
-         call timer('FFTplans ',1)
-
-C  Convert impulse response to filter function
-         do i=1,nfft2
-            cfilt(i)=0.
-         enddo
-         fac=0.00625/nfft1
-         cfilt(1)=fac*halfpulse(1)
-         do i=2,8
-            cfilt(i)=fac*halfpulse(i)
-            cfilt(nfft2+2-i)=fac*halfpulse(i)
-         enddo
-         call timer('FFTfilt ',0)
-         call sfftw_execute(plan5)
-         call timer('FFTfilt ',1)
-
-         base=cfilt(nfft2/2+1)
-         do i=1,nfft2
-            rfilt(i)=real(cfilt(i))-base
-         enddo
-
-         df=96000.d0/nfft1
-         if(nfsample.eq.95238) df=95238.1d0/nfft1
-         first=.false.
-      endif
-
-C  When new data comes along, we need to compute a new "big FFT"
-C  If we just have a new f0, continue with the existing ca and cb.
-
-      if(newdat.ne.0) then
-         nz=min(nmax,nfft1)
-         do i=1,nz
-            ca(i)=cmplx(dd(1,i),dd(2,i))
-            if(xpol) cb(i)=cmplx(dd(3,i),dd(4,i))
-         enddo
-
-         if(nmax.lt.nfft1) then
-            do i=nmax+1,nfft1
-               ca(i)=0.
-               if(xpol) cb(i)=0.
-            enddo
-         endif
-         call timer('FFTbig  ',0)
-         call sfftw_execute(plan1)
-         if(xpol) call sfftw_execute(plan2)
-         call timer('FFTbig  ',1)
-         newdat=0
-      endif
-
-C  NB: f0 is the frequency at which we want our filter centered.
-C      i0 is the bin number in ca and cb closest to f0.
-
-      i0=nint(f0/df) + 1
-      nh=nfft2/2
-      do i=1,nh                                !Copy data into c4a and c4b,
-         j=i0+i-1                              !and apply the filter function
-         if(j.ge.1 .and. j.le.nfft1) then
-            c4a(i)=rfilt(i)*ca(j)
-            if(xpol) c4b(i)=rfilt(i)*cb(j)
-         else
-            c4a(i)=0.
-            if(xpol) c4b(i)=0.
-         endif
-      enddo
-      do i=nh+1,nfft2
-         j=i0+i-1-nfft2
-         if(j.lt.1) j=j+nfft1                  !nfft1 was nfft2
-         c4a(i)=rfilt(i)*ca(j)
-         if(xpol) c4b(i)=rfilt(i)*cb(j)
-      enddo
-
-C  Do the short reverse transform, to go back to time domain.
-      call timer('FFTsmall',0)
-      call sfftw_execute(plan3)
-      if(xpol) call sfftw_execute(plan4)
-      call timer('FFTsmall',1)
-      n4=min(nmax/64,nfft2)
-      go to 999
-
- 900  call sfftw_destroy_plan(plan1)
-      call sfftw_destroy_plan(plan2)
-      call sfftw_destroy_plan(plan3)
-      call sfftw_destroy_plan(plan4)
-      call sfftw_destroy_plan(plan5)
-
- 999  return
-      end
diff --git a/libm65/foldmsk.f90 b/libm65/foldmsk.f90
deleted file mode 100644
index f992957fd..000000000
--- a/libm65/foldmsk.f90
+++ /dev/null
@@ -1,50 +0,0 @@
-subroutine foldmsk(s2,msglen,nchar,mycall,msg,msg29)
-
-! Fold the 2-d "goodness of fit" array s2 modulo message length, 
-! then decode the folded message.
-
-  real s2(0:63,400)
-  real fs2(0:63,29)
-  integer nfs2(29)
-  character mycall*12
-  character msg*400,msg29*29
-  character cc*64
-!                    1         2         3         4         5         6
-!          0123456789012345678901234567890123456789012345678901234567890123
-  data cc/'0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ./?-                 _     @'/
-
-  fs2=0.
-  nfs2=0
-  do j=1,nchar                           !Fold s2 into fs2, modulo msglen
-     jj=mod(j-1,msglen)+1
-     nfs2(jj)=nfs2(jj)+1
-     do i=0,40
-        fs2(i,jj)=fs2(i,jj) + s2(i,j)
-     enddo
-  enddo
-
-  msg=' '
-  do j=1,msglen
-     smax=0.
-     do k=0,40
-        if(fs2(k,j).gt.smax) then
-           smax=fs2(k,j)
-           kpk=k
-        endif
-     enddo
-     if(kpk.eq.40) kpk=57
-     msg(j:j)=cc(kpk+1:kpk+1)
-     if(kpk.eq.57) msg(j:j)=' '
-  enddo
-
-  msg29=msg(1:msglen)
-
-  call alignmsg('  ',2,msg29,msglen,idone)
-  if(idone.eq.0) call alignmsg('CQ',  3,msg29,msglen,idone)
-  if(idone.eq.0) call alignmsg('QRZ', 3,msg29,msglen,idone)
-  if(idone.eq.0) call alignmsg(mycall,4,msg29,msglen,idone)
-  if(idone.eq.0) call alignmsg(' ',   1,msg29,msglen,idone)
-  msg29=adjustl(msg29)
-
-  return
-end subroutine foldmsk
diff --git a/libm65/gen65.f90 b/libm65/gen65.f90
deleted file mode 100644
index 92c45e7fe..000000000
--- a/libm65/gen65.f90
+++ /dev/null
@@ -1,93 +0,0 @@
-subroutine gen65(message,mode65,samfac,nsendingsh,msgsent,iwave,nwave)
-
-! Encodes a JT65 message into a wavefile.  
-! Executes in 17 ms on opti-745.
-
-  parameter (NMAX=60*11025)     !Max length of wave file
-  character*22 message          !Message to be generated
-  character*22 msgsent          !Message as it will be received
-  character*3 cok               !'   ' or 'OOO'
-  real*8 dt,phi,f,f0,dfgen,dphi,twopi,samfac
-  integer*2 iwave(NMAX)         !Generated wave file
-  integer dgen(12)
-  integer sent(63)
-  logical first
-  integer nprc(126)
-  real pr(126)
-  data nprc/1,0,0,1,1,0,0,0,1,1,1,1,1,1,0,1,0,1,0,0,  &
-            0,1,0,1,1,0,0,1,0,0,0,1,1,1,0,0,1,1,1,1,  &
-            0,1,1,0,1,1,1,1,0,0,0,1,1,0,1,0,1,0,1,1,  &
-            0,0,1,1,0,1,0,1,0,1,0,0,1,0,0,0,0,0,0,1,  &
-            1,0,0,0,0,0,0,0,1,1,0,1,0,0,1,0,1,1,0,1,  &
-            0,1,0,1,0,0,1,1,0,0,1,0,0,1,0,0,0,0,1,1,  &
-            1,1,1,1,1,1/
-  data twopi/6.283185307179586476d0/,first/.true./
-  save
-
-  if(first) then
-     do i=1,126
-        pr(i)=2*nprc(i)-1
-     enddo
-     first=.false.
-  endif
-
-  call chkmsg(message,cok,nspecial,flip)
-  if(nspecial.eq.0) then
-     call packmsg(message,dgen)          !Pack message into 72 bits
-     nsendingsh=0
-     if(iand(dgen(10),8).ne.0) nsendingsh=-1    !Plain text flag
-
-     call rs_encode(dgen,sent)
-     call interleave63(sent,1)           !Apply interleaving
-     call graycode(sent,63,1)            !Apply Gray code
-     nsym=126                            !Symbols per transmission
-     nsps=4096
-  else
-     nsym=32
-     nsps=16384
-     nsendingsh=1                         !Flag for shorthand message
-  endif
-  if(mode65.eq.0) go to 900
-
-! Set up necessary constants
-  dt=1.d0/(samfac*11025.d0)
-  f0=118*11025.d0/1024
-  dfgen=mode65*11025.d0/4096.d0
-  phi=0.d0
-  i=0
-  k=0
-  do j=1,nsym
-     f=f0
-     if(nspecial.ne.0 .and. mod(j,2).eq.0) f=f0+10*nspecial*dfgen
-     if(nspecial.eq.0 .and. flip*pr(j).lt.0.0) then
-        k=k+1
-        f=f0+(sent(k)+2)*dfgen
-     endif
-     dphi=twopi*dt*f
-     do ii=1,nsps
-        phi=phi+dphi
-        if(phi.gt.twopi) phi=phi-twopi
-        xphi=phi
-        i=i+1
-        iwave(i)=32767.0*sin(xphi)
-     enddo
-  enddo
-
-  iwave(nsym*nsps+1:)=0
-  nwave=nsym*nsps + 5512
-  call unpackmsg(dgen,msgsent)
-  if(flip.lt.0.0) then
-     do i=22,1,-1
-        if(msgsent(i:i).ne.' ') goto 10
-     enddo
-10   msgsent=msgsent(1:i)//' OOO'
-  endif
-
-  if(nsendingsh.eq.1) then
-     if(nspecial.eq.2) msgsent='RO'
-     if(nspecial.eq.3) msgsent='RRR'
-     if(nspecial.eq.4) msgsent='73'
-  endif
-
-900 return
-end subroutine gen65
diff --git a/libm65/genjtms3.f90 b/libm65/genjtms3.f90
deleted file mode 100644
index b5c889cb9..000000000
--- a/libm65/genjtms3.f90
+++ /dev/null
@@ -1,158 +0,0 @@
-subroutine genjtms3(msg28,iwave,nwave)
-!subroutine genjtms3(msg28,iwave,cwave,isrch,nwave)
-
-! Generate a JTMS3 wavefile.
-
-  parameter (NMAX=30*48000)     !Max length of wave file
-  integer*2 iwave(NMAX)         !Generated wave file
-  complex cwave(NMAX)           !Alternative for searchms
-  character*28 msg28            !User message
-  character*29 msg
-  character cc*64
-  integer sentsym(203)          !Transmitted symbols (0/1)
-  real sentsam(4872)            !Transmitted waveform
-  real*8 dt,phi,f0,dphi,pi,twopi,samfac
-  real p(0:420)
-  real carrier(4872)
-  real dat(4872),bb(4872),wave(4872)
-  complex cdat(0:2436)
-  logical first
-  integer np(9)
-  data np/5,7,9,11,13,17,19,23,29/  !Permissible message lengths
-!                   1         2         3         4         5         6
-!          0123456789012345678901234567890123456789012345678901234567890123
-  data cc/'0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ./?-                 _     @'/
-  data samfac/1.d0/,first/.true./
-  equivalence (dat,cdat)
-  save
-
-  sinc(x)=sin(pi*x)/(pi*x)
-
-  if(first) then
-     pi=4.d0*atan(1.d0)
-     twopi=2.d0*pi
-     x=0.
-     dx=1.0/24.0
-     width=3.0
-     do i=1,420                                !Generate the BPSK pulse shape
-        x=x+dx
-        fac=0.0
-        if(x/width.lt.0.5*pi) then
-           fac=(cos(x/width))**2
-           ipz=i
-        endif
-        p(i)=fac*sinc(x)
-     enddo
-     p(0)=1.0
-
-     f0=10000.d0/7.d0
-     dt=1.d0/48000.d0
-     dphi=twopi*f0*dt
-     phi=0.d0
-     do i=1,4872                             !Generate the carrier
-        phi=phi+dphi
-        if(phi.gt.twopi)phi=phi-twopi
-        xphi=phi
-        carrier(i)=sin(xphi)
-     enddo
-
-     first=.false.
-  endif
-
-  msg=msg28//' '                               !Extend to 29 characters
-  do i=28,1,-1                                 !Find user's message length
-     if(msg(i:i).ne.' ') go to 1
-  enddo
-1 iz=i+1                                       !Add one for space at EOM
-  msglen=iz
-  if(isrch.ne.0) go to 3
-  do i=1,9
-     if(np(i).ge.iz) go to 2
-  enddo
-  i=8
-2 msglen=np(i)
-
-! Convert message to a bit sequence, 7 bits per character (6 + even parity)
-3  sentsym=0
-  k=0
-  do j=1,msglen
-     if(msg(j:j).eq.' ') then
-        i=58
-        go to 5
-     else
-        do i=1,64
-           if(msg(j:j).eq.cc(i:i)) go to 5
-        enddo
-     endif
-5    m=0
-     do n=5,0,-1                            !Each character gets 6 bits
-        k=k+1
-        sentsym(k)=iand(1,ishft(i-1,-n))
-        m=m+sentsym(k)
-     enddo
-     k=k+1
-     sentsym(k)=iand(m,1)                   !Insert bit for even parity
-  enddo
-  nsym=7*msglen                             !# symbols in message
-  nsam=24*nsym                              !# samples in message
-
-  bb(1:nsam)=0.
-  do j=1,nsym
-     fac=1.0
-     if(sentsym(j).eq.0) fac=-1.0
-     k0=24*j - 23
-     do i=0,ipz
-        k=k0+i
-        if(k.gt.nsam) k=k-nsam
-        bb(k)=bb(k) + fac*p(i)
-        if(i.gt.0) then
-           k=k0-i
-           if(k.lt.1) k=k+nsam
-           bb(k)=bb(k) + fac*p(i)
-        endif
-     enddo
-  enddo
-
-  sq=0.
-  wmax=0.
-  do i=1,nsam
-     wave(i)=carrier(i)*bb(i)
-     sq=sq + wave(i)**2
-     wmax=max(wmax,abs(wave(i)))
-!     write(15,3002) i*dt,bb(i),wave(i)
-!3002 format(f12.6,2f12.3)
-  enddo
-
-  rms=sqrt(sq/nsam)
-!  print*,rms,wmax,wmax/rms
-
-  fac=32767.0/wmax
-  iwave(1:nsam)=fac*wave(1:nsam)
-
-!  nwave=nsam
-
-  nblk=30*48000/nsam
-  do n=2,nblk
-     i0=(n-1)*nsam
-     iwave(i0+1:i0+nsam)=iwave(1:nsam)
-  enddo
-  nwave=i0+nsam
-
-! Compute the spectrum
-!  nfft=nsam
-!  df=48000.0/nfft
-!  ib=4000.0/df
-!  fac=10.0/nfft
-!  dat(1:nfft)=fac*bb(1:nfft)
-!  call four2a(dat,nfft,1,-1,0)
-!  do i=0,ib
-!     sq=real(cdat(i))**2 + aimag(cdat(i))**2
-!     write(14,3010) i*df,sq,10.0*log10(sq)
-!3010 format(3f12.3)
-!  enddo
-
-!  if(isrch.eq.0) iwave(k+1:)=0
-!  nwave=k
-
-  return
-end subroutine genjtms3
diff --git a/libm65/geocentric.f b/libm65/geocentric.f
deleted file mode 100644
index 1767f9cf8..000000000
--- a/libm65/geocentric.f
+++ /dev/null
@@ -1,17 +0,0 @@
-      subroutine geocentric(alat,elev,hlt,erad)
-
-      implicit real*8 (a-h,o-z)
-
-C  IAU 1976 flattening f, equatorial radius a
-      f = 1.d0/298.257d0
-      a = 6378140.d0
-      c = 1.d0/sqrt(1.d0 + (-2.d0 + f)*f*sin(alat)*sin(alat))
-      arcf = (a*c + elev)*cos(alat)
-      arsf = (a*(1.d0 - f)*(1.d0 - f)*c + elev)*sin(alat)
-      hlt = datan2(arsf,arcf)
-      erad = sqrt(arcf*arcf + arsf*arsf)
-      erad = 0.001d0*erad
-
-      return
-      end
-
diff --git a/libm65/getdphi.f90 b/libm65/getdphi.f90
deleted file mode 100644
index e55d7120d..000000000
--- a/libm65/getdphi.f90
+++ /dev/null
@@ -1,18 +0,0 @@
-subroutine getdphi(qphi)
-
-  real qphi(12)
-
-  s=0.
-  c=0.
-  do i=1,12
-     th=i*30/57.2957795
-     s=s+qphi(i)*sin(th)
-     c=c+qphi(i)*cos(th)
-  enddo
-
-  dphi=57.2957795*atan2(s,c)
-  write(*,1010) nint(dphi)
-1010 format('!Best-fit Dphi =',i4,' deg')
-
-  return
-  end
diff --git a/libm65/hipass.f90 b/libm65/hipass.f90
deleted file mode 100644
index 9a24fa89a..000000000
--- a/libm65/hipass.f90
+++ /dev/null
@@ -1,23 +0,0 @@
-subroutine hipass(y,npts,nwidth)
-
-!  Hipass filter for time-domain data.  Removes an RC-type running 
-!  mean (time constant nwidth) from array y(1:npts).  
-
-  real y(npts)
-
-  c1=1.0/nwidth
-  c2=1.0-c1
-  s=0.
-  do i=1,nwidth                      !Get initial average
-     s=s+y(i)
-  enddo
-  ave=c1*s
-
-  do i=1,npts                        !Do the filtering
-     y0=y(i)
-     y(i)=y0-ave                     !Remove the mean
-     ave=c1*y0 + c2*ave              !Update the mean
-  enddo
-
-return
-end subroutine hipass
diff --git a/libm65/interleave63.f b/libm65/interleave63.f
deleted file mode 100644
index f6793023a..000000000
--- a/libm65/interleave63.f
+++ /dev/null
@@ -1,25 +0,0 @@
-      subroutine interleave63(d1,idir)
-
-C  Interleave (idir=1) or de-interleave (idir=-1) the array d1.
-
-      integer d1(0:6,0:8)
-      integer d2(0:8,0:6)
-
-      if(idir.ge.0) then
-         do i=0,6
-            do j=0,8
-               d2(j,i)=d1(i,j)
-            enddo
-         enddo
-         call move(d2,d1,63)
-      else
-         call move(d1,d2,63)
-         do i=0,6
-            do j=0,8
-               d1(i,j)=d2(j,i)
-            enddo
-         enddo
-      endif
-         
-      return
-      end
diff --git a/libm65/iqcal.f90 b/libm65/iqcal.f90
deleted file mode 100644
index c0c4fce2d..000000000
--- a/libm65/iqcal.f90
+++ /dev/null
@@ -1,30 +0,0 @@
-subroutine iqcal(nn,c,nfft,gain,phase,zsum,ipk,reject)
-
-  complex c(0:nfft-1)
-  complex z,zsum,zave
-
-  if(nn.eq.0) then
-     zsum=0.
-  endif
-  nn=nn+1
-  smax=0.
-  ipk=1
-  do i=1,nfft-1                       !Find strongest signal
-     s=real(c(i))**2 + aimag(c(i))**2
-     if(s.gt.smax) then
-        smax=s
-        ipk=i
-     endif
-  enddo
-  pimage=real(c(nfft-ipk))**2 + aimag(c(nfft-ipk))**2
-  p=smax + pimage
-  z=c(ipk)*c(nfft-ipk)/p              !Synchronous detection of image
-  zsum=zsum+z
-  zave=zsum/nn
-  tmp=sqrt(1.0 - (2.0*real(zave))**2)
-  phase=asin(2.0*aimag(zave)/tmp)    !Estimate phase
-  gain=tmp/(1.0-2.0*real(zave))      !Estimate gain
-  reject=10.0*log10(pimage/smax)
-
-  return
-end subroutine iqcal
diff --git a/libm65/iqfix.f90 b/libm65/iqfix.f90
deleted file mode 100644
index 336f7928a..000000000
--- a/libm65/iqfix.f90
+++ /dev/null
@@ -1,29 +0,0 @@
-subroutine iqfix(c,nfft,gain,phase)
-
-  complex c(0:nfft-1)
-  complex z,h,u,v
-  real*8 sq1,sq2
-
-  nh=nfft/2
-  h=gain*cmplx(cos(phase),sin(phase))
-
-  do i=1,nh-1
-     u=c(i)
-     v=c(nfft-i)
-     x=real(u)  + real(v)  - (aimag(u) + aimag(v))*aimag(h) +         &
-          (real(u) - real(v))*real(h)
-     y=aimag(u) - aimag(v) + (aimag(u) + aimag(v))*real(h)  +         &
-          (real(u) - real(v))*aimag(h)
-     c(i)=0.5*cmplx(x,y)
-     z=u
-     u=v
-     v=z
-     x=real(u)  + real(v)  - (aimag(u) + aimag(v))*aimag(h) +         &
-          (real(u) - real(v))*real(h)
-     y=aimag(u) - aimag(v) + (aimag(u) + aimag(v))*real(h)  +         &
-          (real(u) - real(v))*aimag(h)
-     c(nfft-i)=0.5*cmplx(x,y)
-  enddo
-
-  return
-end subroutine iqfix
diff --git a/libm65/jtmsk.f90 b/libm65/jtmsk.f90
deleted file mode 100644
index ab3e048f2..000000000
--- a/libm65/jtmsk.f90
+++ /dev/null
@@ -1,72 +0,0 @@
-subroutine jtmsk(dat,npts,cfile6,tpk,mswidth,ndb,nrpt,Nfreeze,       &
-     ntol,MouseDF,pick,mycall)
-
-! Decode a JTMSK ping
-
-  parameter (NZ=2*48000)
-  real dat(npts)                        !Raw data
-  complex cdat(NZ)                      !Analytic form of signal
-  character*6 cfile6                    !FileID
-  logical pick
-  character*12 mycall
-  real s(NZ)                            !Power spectrum
-  real s2(0:63,400)
-  real r(60000)
-  complex cw(168,0:63)                   !Complex waveforms for all codewords
-  complex cwb(168)                       !Complex waveform for <space>
-  logical first
-  character msg*400,msg29*29
-  data first/.true./
-  save first,cw,cwb
-  save cdat                             !Fix its address, for four2
-
-  if(first) call setupmsk(cw,cwb)       !Calculate waveforms for codewords
-  first=.false.
-
-  nsps=24                               !Samples per symbol
-  f0=1000.0                             !Nominal frequency for bit=0
-  n=log(float(npts))/log(2.0) + 1.0
-  nfft1=2**n                            !FFT length
-  call analytic(dat,npts,nfft1,s,cdat)  !Convert to analytic signal
-
-  call mskdf(cdat,npts,tpk,nfft1,f0,nfreeze,mousedf,ntol,     &
-       dfx,snrsq2)                      !Get DF
-  sq2lim=7.0
-  if(pick) sq2lim=5.0
-  if(snrsq2.lt.sq2lim) go to 900        !Reject non-JTMS signals
-
-  call tweak1(cdat,npts,-dfx,cdat)      !Mix to standard frequency
-
-  call syncmsk(cdat,npts,cwb,r,i1)      !Get character sync
-
-  call lenmsk(r,npts,msglen)            !Find message length
-  
-  s2=0.
-  nchar=(npts-168+1-i1)/168
-  if(nchar.gt.400) nchar=400
-
-  call decodemsk(cdat,npts,cw,i1,nchar,s2,msg)   !Decode the message
-
-  ia=1
-  if(nchar.ge.40) ia=min(nchar/3,nchar-28)    
-  ib=min(ia+28,nchar)                   !Can better limits ia, ib be found?
-  msg29=adjustl(msg(ia:ib))
-  msg=adjustl(msg)
-  ib=min(nchar,45)
-  ndf=nint(dfx)
-  nchk=max(20,nint(1.5*msglen))
-
-  if(msglen.eq.0 .or. nchar.lt.nchk .or. pick) then
-     write(*,1110) cfile6,tpk,mswidth,ndb,nrpt,ndf,msg(1:45)
-1110 format(a6,f5.1,i5,i3,1x,i2.2,i5,5x,a45)
-  endif
-
-  if(msglen.gt.0 .and. nchar.ge.nchk) then
-     call foldmsk(s2,msglen,nchar,mycall,msg,msg29)   !Decode folded message
-     write(*,1120) cfile6,tpk,mswidth,ndb,nrpt,ndf,msg29
-1120 format(a6,f5.1,i5,i3,1x,i2.2,i5,5x,a29,11x,'*')
-   endif
-
-900 continue
-  return
-end subroutine jtmsk
diff --git a/libm65/lenmsk.f90 b/libm65/lenmsk.f90
deleted file mode 100644
index 4f4ea09c2..000000000
--- a/libm65/lenmsk.f90
+++ /dev/null
@@ -1,56 +0,0 @@
-subroutine lenmsk(r,npts,msglen)
-
-! Determine length of the user message in a JTMS ping.
-
-  real r(60000)
-  real acf(4872)
-  integer np(9)
-  data np/5,7,9,11,13,17,19,23,29/      !Permissible message lengths
-  save acf                              !Why necessary?  (But don't remove!)
-
-  msglen=0                              !Use ACF to find msg length
-  if(npts.ge.8*168) then
-     r=r-sum(r(1:npts))/npts
-     acfmax=0.
-     acf0=dot_product(r(1:npts),r(1:npts))
-     kz=min(nint(0.75*npts),29*168)
-     do k=8,kz
-        fac=float(npts)/(npts-k)
-        acf(k)=fac*dot_product(r(1:npts),r(1+k:npts+k))/acf0
-     enddo
-     call hipass(acf(8),kz-7,50)
-
-     do k=8,kz                          !Find acfmax, kpk
-        if(acf(k).gt.acfmax) then
-           acfmax=acf(k)
-           kpk=k
-        endif
-     enddo
-
-     sumsq=0.
-     n=0
-     do k=8,kz                          !Find rms, skipping around kpk
-        if(abs(k-kpk).gt.10) then
-           sumsq=sumsq+acf(k)**2
-           n=n+1
-        endif
-     enddo
-     rms=sqrt(sumsq/n)
-     acf=acf/rms                        !Normalize the acf
-
-     amax2=0.
-     acflim=3.5
-     do i=1,9
-        k=168*np(i)                     !Check only the permitted lengths
-        if(k.gt.kz) go to 10
-        if(acf(k).gt.acflim .and. acf(k).gt.amax2) then  
-           amax2=acf(k)                 !Save best value >3.5 sigma
-           msglen=np(i)                 !Save message length
-           kpk2=k
-        endif
-     enddo
-10   continue
-  endif
-
-  return
-end subroutine lenmsk
diff --git a/libm65/m65.f90 b/libm65/m65.f90
deleted file mode 100644
index 27d9b373c..000000000
--- a/libm65/m65.f90
+++ /dev/null
@@ -1,130 +0,0 @@
-program m65
-
-! Decoder for map65.  Can run stand-alone, reading data from *.tf2 files;
-! or as the back end of map65, with data placed in a shared memory region.
-
-  parameter (NSMAX=60*96000)
-  parameter (NFFT=32768)
-  integer*2 i2(4,87)
-  real*8 hsym
-  real*4 ssz5a(NFFT)
-  logical*1 lstrong(0:1023)
-  common/tracer/limtrace,lu
-  real*8 fc0,fcenter
-  character*80 arg,infile
-  character mycall*12,hiscall*12,mygrid*6,hisgrid*6,datetime*20
-  common/datcom/dd(4,5760000),ss(4,322,NFFT),savg(4,NFFT),fc0,nutc0,junk(34)
-  common/npar/fcenter,nutc,idphi,mousedf,mousefqso,nagain,                &
-       ndepth,ndiskdat,neme,newdat,nfa,nfb,nfcal,nfshift,                 &
-       mcall3,nkeep,ntol,nxant,nrxlog,nfsample,nxpol,mode65,              &
-       mycall,mygrid,hiscall,hisgrid,datetime
-
-  nargs=iargc()
-  if(nargs.lt.1) then
-     print*,'Usage: m65 [95238] file1 [file2 ...]'
-     print*,'       Reads data from *.tf2 files.'
-     print*,''
-     print*,'       m65 -s'
-     print*,'       Gets data from shared memory region.'
-     go to 999
-  endif
-  call getarg(1,arg)
-  if(arg(1:2).eq.'-s') then
-     call m65a
-     call ftnquit
-     go to 999
-  endif
-  nfsample=96000
-  nxpol=1
-  mode65=2
-  ifile1=1
-  if(arg.eq.'95238') then
-     nfsample=95238
-     call getarg(2,arg)
-     ifile1=2
-  endif
-
-  limtrace=0
-  lu=12
-  nfa=100
-  nfb=162
-  nfshift=6
-  ndepth=2
-  nfcal=344
-  idphi=-50
-  ntol=500
-  nkeep=10
-
-  call ftninit('.')
-
-  do ifile=ifile1,nargs
-     call getarg(ifile,infile)
-     open(10,file=infile,access='stream',status='old',err=998)
-     i1=index(infile,'.tf2')
-     read(infile(i1-4:i1-1),*,err=1) nutc0
-     go to 2
-1    nutc0=0
-2    hsym=2048.d0*96000.d0/11025.d0          !Samples per half symbol
-     nhsym0=-999
-     k=0
-     fcenter=144.125d0
-     mousedf=0
-     mousefqso=125
-     newdat=1
-     mycall='K1JT'
-
-     if(ifile.eq.ifile1) call timer('m65     ',0)
-     do irec=1,9999999
-        call timer('read_tf2',0)
-        read(10) i2
-        call timer('read_tf2',1)
-        
-        call timer('float   ',0)
-        do i=1,87
-           k=k+1
-           dd(1,k)=i2(1,i)
-           dd(2,k)=i2(2,i)
-           dd(3,k)=i2(3,i)
-           dd(4,k)=i2(4,i)
-        enddo
-        call timer('float   ',1)
-        nhsym=(k-2048)/hsym
-        if(nhsym.ge.1 .and. nhsym.ne.nhsym0) then
-           ndiskdat=1
-           nb=0
-! Emit signal readyForFFT
-           call timer('symspec ',0)
-           fgreen=-13.0
-           iqadjust=1
-           iqapply=1
-           nbslider=100
-           gainx=0.9962
-           gainy=1.0265
-           phasex=0.01426
-           phasey=-0.01195
-           call symspec(k,nxpol,ndiskdat,nb,nbslider,idphi,nfsample,fgreen,  &
-                iqadjust,iqapply,gainx,gainy,phasex,phasey,rejectx,rejecty,  &
-                pxdb,pydb,ssz5a,nkhz,ihsym,nzap,slimit,lstrong)
-           call timer('symspec ',1)
-           nhsym0=nhsym
-           if(ihsym.ge.278) go to 10
-        endif
-     enddo
-
-10   continue
-     if(iqadjust.ne.0) write(*,3002) rejectx,rejecty
-3002 format('Image rejection:',2f7.1,' dB')
-     nutc=nutc0
-     nstandalone=1
-     call decode0(dd,ss,savg,nstandalone,nfsample)
-  enddo
-
-  call timer('m65     ',1)
-  call timer('m65     ',101)
-  call ftnquit
-  go to 999
-
-998 print*,'Cannot open file:'
-  print*,infile
-
-999 end program m65
diff --git a/libm65/m65a.F90 b/libm65/m65a.F90
deleted file mode 100644
index b5fc4b11f..000000000
--- a/libm65/m65a.F90
+++ /dev/null
@@ -1,97 +0,0 @@
-subroutine m65a
-
-! NB: this interface block is required by g95, but must be omitted
-!     for gfortran.  (????)
-
-#ifndef UNIX
-  interface
-     function address_m65()
-     end function address_m65
-  end interface
-#endif
-  
-  integer*1 attach_m65,lock_m65,unlock_m65
-  integer size_m65
-  integer*1, pointer :: address_m65,p_m65
-  character*80 cwd
-  logical fileExists
-  common/tracer/limtrace,lu
-
-  call getcwd(cwd)
-  call ftninit(trim(cwd))
-  limtrace=0
-  lu=12
-  i1=attach_m65()
-
-10 inquire(file=trim(cwd)//'/.lock',exist=fileExists)
-  if(fileExists) then
-     call sleep_msec(100)
-     go to 10
-  endif
-
-  inquire(file=trim(cwd)//'/.quit',exist=fileExists)
-  if(fileExists) then
-     call ftnquit
-     i=detach_m65()
-     go to 999
-  endif
-  
-  nbytes=size_m65()
-  if(nbytes.le.0) then
-     print*,'m65a: Shared memory mem_m65 does not exist.' 
-     print*,'Program m65a should be started automatically from within map65.'
-     go to 999
-  endif
-  p_m65=>address_m65()
-  call m65b(p_m65,nbytes)
-
-  write(*,1010) 
-1010 format('<m65aFinished>')
-  flush(6)
-
-100 inquire(file=trim(cwd)//'/.lock',exist=fileExists)
-  if(fileExists) go to 10
-  call sleep_msec(100)
-  go to 100
-
-999 return
-end subroutine m65a
-
-subroutine m65b(m65com,nbytes)
-  integer*1 m65com(0:nbytes-1)
-  kss=4*4*60*96000
-  ksavg=kss+4*4*322*32768
-  kfcenter=ksavg+4*4*32768
- call m65c(m65com(0),m65com(kss),m65com(ksavg),m65com(kfcenter))
-  return
-end subroutine m65b
-
-subroutine m65c(dd,ss,savg,nparams0)
-  integer*1 detach_m65
-  real*4 dd(4,5760000),ss(4,322,32768),savg(4,32768)
-  real*8 fcenter
-  integer nparams0(37),nparams(37)
-  character*12 mycall,hiscall
-  character*6 mygrid,hisgrid
-  character*20 datetime
-  common/npar/fcenter,nutc,idphi,mousedf,mousefqso,nagain,              &
-       ndepth,ndiskdat,neme,newdat,nfa,nfb,nfcal,nfshift,               &
-       mcall3,nkeep,ntol,nxant,nrxlog,nfsample,nxpol,mode65,            &
-       mycall,mygrid,hiscall,hisgrid,datetime
-  equivalence (nparams,fcenter)
-  
-  nparams=nparams0                     !Copy parameters into common/npar/
-  npatience=1
-  if(iand(nrxlog,1).ne.0) then
-     write(21,1000) datetime(:17)
-1000 format(/'UTC Date: 'a17/78('-'))
-     flush(21)
-  endif
-  if(iand(nrxlog,2).ne.0) rewind 21
-  if(iand(nrxlog,4).ne.0) rewind 26
-
-  nstandalone=0
-  if(sum(nparams).ne.0) call decode0(dd,ss,savg,nstandalone)
-
-  return
-end subroutine m65c
diff --git a/libm65/makepings.f90 b/libm65/makepings.f90
deleted file mode 100644
index eb96fa93e..000000000
--- a/libm65/makepings.f90
+++ /dev/null
@@ -1,29 +0,0 @@
-subroutine makepings(iwave,nwave)
-
-  parameter (NFSAMPLE=48000)
-  integer*2 iwave(nwave)
-  real*8 t
-
-  iping0=-999
-  dt=1.0/NFSAMPLE
-  do i=1,nwave
-     iping=i/(3*NFSAMPLE)
-     if(iping.ne.iping0) then
-        ip=mod(iping,3)
-        w=0.015*(4-ip)
-        ig=(iping-1)/3
-        amp=sqrt((3.0-ig)/3.0)
-        t0=dt*(iping+0.5)*(3*NFSAMPLE)
-        iping0=iping
-     endif
-     t=(i*dt-t0)/w
-     if(t.lt.0.d0 .and. t.lt.10.0) then
-        fac=0.
-     else
-        fac=2.718*t*dexp(-t)
-     endif
-     iwave(i)=nint(fac*amp*iwave(i))
-  enddo
-
-  return
-end subroutine makepings
diff --git a/libm65/map65a.f90 b/libm65/map65a.f90
deleted file mode 100644
index 282b65836..000000000
--- a/libm65/map65a.f90
+++ /dev/null
@@ -1,438 +0,0 @@
-subroutine map65a(dd,ss,savg,newdat,nutc,fcenter,ntol,idphi,nfa,nfb,        &
-     mousedf,mousefqso,nagain,ndecdone,ndiskdat,nfshift,ndphi,              &
-     nfcal,nkeep,mcall3b,nsave,nxant,rmsdd,mycall,mygrid,                    &
-     neme,ndepth,hiscall,hisgrid,nhsym,nfsample,nxpol,mode65)
-
-!  Processes timf2 data from Linrad to find and decode JT65 signals.
-
-  parameter (MAXMSG=1000)            !Size of decoded message list
-  parameter (NSMAX=60*96000)
-  parameter (NFFT=32768)
-  real dd(4,NSMAX)
-  real*4 ss(4,322,NFFT),savg(4,NFFT)
-  real tavg(-50:50)                  !Temp for finding local base level
-  real base(4)                       !Local basel level at 4 pol'ns
-  real tmp (200)                     !Temp storage for pctile sorting
-  real sig(MAXMSG,30)                !Parameters of detected signals
-  real a(5)
-  real*8 fcenter
-  character*22 msg(MAXMSG)
-  character*3 shmsg0(4)
-  character mycall*12,hiscall*12,mygrid*6,hisgrid*6,grid*6,cp*1
-  integer indx(MAXMSG),nsiz(MAXMSG)
-  logical done(MAXMSG)
-  logical xpol
-  character decoded*22,blank*22
-  real short(3,NFFT)                 !SNR dt ipol for potential shorthands
-  real qphi(12)
-  common/c3com/ mcall3a
-  common/testcom/ifreq
-  
-  data blank/'                      '/
-  data shmsg0/'ATT','RO ','RRR','73 '/
-  data nfile/0/,nutc0/-999/,nid/0/,ip000/1/,ip001/1/,mousefqso0/-999/
-  save
-
-  mcall3a=mcall3b
-  mousefqso0=mousefqso
-  xpol=(nxpol.ne.0)
-  if(.not.xpol) ndphi=0
-
-!### Should use AppDir! ###
-!  open(23,file='release/CALL3.TXT',status='unknown')
-  open(23,file='CALL3.TXT',status='unknown')
-
-  if(nutc.ne.nutc0) nfile=nfile+1
-  nutc0=nutc
-  df=96000.0/NFFT                     !df = 96000/NFFT = 2.930 Hz
-  if(nfsample.eq.95238) df=95238.1/NFFT
-  ftol=0.010                          !Frequency tolerance (kHz)
-  dphi=idphi/57.2957795
-  foffset=0.001*(1270 + nfcal)              !Offset from sync tone, plus CAL
-  fqso=mousefqso + foffset - 0.5*(nfa+nfb) + nfshift !fqso at baseband (khz)
-  iloop=0
-
-2  if(ndphi.eq.1) dphi=30*iloop/57.2957795
-
-  do nqd=1,0,-1
-     if(nqd.eq.1) then                     !Quick decode, at fQSO
-        fa=1000.0*(fqso+0.001*mousedf) - ntol
-        fb=1000.0*(fqso+0.001*mousedf) + ntol + 4*53.8330078
-     else                                  !Wideband decode at all freqs
-        fa=-1000*0.5*(nfb-nfa) + 1000*nfshift
-        fb= 1000*0.5*(nfb-nfa) + 1000*nfshift
-     endif
-     ia=nint(fa/df) + 16385
-     ib=nint(fb/df) + 16385
-     ia=max(51,ia)
-     ib=min(32768-51,ib)
-
-     km=0
-     nkm=1
-     nz=n/8
-     freq0=-999.
-     sync10=-999.
-     fshort0=-999.
-     syncshort0=-999.
-     ntry=0
-     short=0.                                 !Zero the whole short array
-     jpz=1
-     if(xpol) jpz=4
-
-     do i=ia,ib                               !Search over freq range
-        freq=0.001*(i-16385)*df
-!  Find the local base level for each polarization; update every 10 bins.
-        if(mod(i-ia,10).eq.0) then
-           do jp=1,jpz
-              do ii=-50,50
-                 iii=i+ii
-                 if(iii.ge.1 .and. iii.le.32768) then
-                    tavg(ii)=savg(jp,iii)
-                 else
-                    write(13,*) ,'Error in iii:',iii,ia,ib,fa,fb
-                    flush(13)
-                    go to 999
-                 endif
-              enddo
-              call pctile(tavg,tmp,101,50,base(jp))
-           enddo
-        endif
-
-!  Find max signal at this frequency
-        smax=0.
-        do jp=1,jpz
-           if(savg(jp,i)/base(jp).gt.smax) then
-              smax=savg(jp,i)/base(jp)
-              jpmax=jp
-           endif
-        enddo
-
-        if(smax.gt.1.1) then
-
-!  Look for JT65 sync patterns and shorthand square-wave patterns.
-           call timer('ccf65   ',0)
-!           ssmax=4.0*(rmsdd/22.5)**2
-           ssmax=savg(jpmax,i)
-           call ccf65(ss(1,1,i),nhsym,ssmax,sync1,ipol,jpz,dt,flipk,     &
-                syncshort,snr2,ipol2,dt2)
-           call timer('ccf65   ',1)
-
-! ########################### Search for Shorthand Messages #################
-!  Is there a shorthand tone above threshold?
-           thresh0=1.0
-!  Use lower thresh0 at fQSO
-           if(nqd.eq.1 .and. ntol.le.100) thresh0=0.
-           if(syncshort.gt.thresh0) then
-! ### Do shorthand AFC here (or maybe after finding a pair?) ###
-              short(1,i)=syncshort
-              short(2,i)=dt2
-              short(3,i)=ipol2
-           
-!  Check to see if lower tone of shorthand pair was found.
-              do j=2,4
-                 i0=i-nint(j*mode65*10.0*(11025.0/4096.0)/df)
-!  Should this be i0 +/- 1, or just i0?
-!  Should we also insist that difference in DT be either 1.5 or -1.5 s?
-                 if(short(1,i0).gt.thresh0) then
-                    fshort=0.001*(i0-16385)*df
-                    noffset=0
-                    if(nqd.eq.1) noffset=nint(1000.0*(fshort-fqso)-mousedf)
-                    if(abs(noffset).le.ntol) then
-!  Keep only the best candidate within ftol.
-!### NB: sync2 was not defined here!
-!                       sync2=syncshort                   !### try this ???
-                       if(fshort-fshort0.le.ftol .and. syncshort.gt.syncshort0 &
-                            .and. nkm.eq.2) km=km-1
-                       if(fshort-fshort0.gt.ftol .or.                     &
-                            syncshort.gt.syncshort0) then
-                          if(km.lt.MAXMSG) km=km+1
-                          sig(km,1)=nfile
-                          sig(km,2)=nutc
-                          sig(km,3)=fshort + 0.5*(nfa+nfb)
-                          sig(km,4)=syncshort
-                          sig(km,5)=dt2
-                          sig(km,6)=45*(ipol2-1)/57.2957795
-                          sig(km,7)=0
-                          sig(km,8)=snr2
-                          sig(km,9)=0
-                          sig(km,10)=0
-!                           sig(km,11)=rms0
-                          sig(km,12)=savg(ipol2,i)
-                          sig(km,13)=0
-                          sig(km,14)=0
-                          sig(km,15)=0
-                          sig(km,16)=0
-!                           sig(km,17)=0
-                          sig(km,18)=0
-                          msg(km)=shmsg0(j)
-                          fshort0=fshort
-                          syncshort0=syncshort
-                          nkm=2
-                       endif
-                    endif
-                 endif
-              enddo
-           endif
-
-! ########################### Search for Normal Messages ###########
-!  Is sync1 above threshold?
-           thresh1=1.0
-!  Use lower thresh1 at fQSO
-           if(nqd.eq.1 .and. ntol.le.100) thresh1=0.
-           noffset=0
-           if(nqd.eq.1) noffset=nint(1000.0*(freq-fqso)-mousedf)
-
-           if(sync1.gt.thresh1 .and. abs(noffset).le.ntol) then
-!  Keep only the best candidate within ftol.
-!  (Am I deleting any good decodes by doing this?)
-              if(freq-freq0.le.ftol .and. sync1.gt.sync10 .and.       &
-                   nkm.eq.1) km=km-1
-              if(freq-freq0.gt.ftol .or. sync1.gt.sync10) then
-                 nflip=nint(flipk)
-                 f00=(i-1)*df          !Freq of detected sync tone (0-96000 Hz)
-                 ntry=ntry+1
-                 if((nqd.eq.1 .and. ntry.ge.40) .or.                  &
-                          (nqd.eq.0 .and. ntry.ge.400)) then
-! Too many calls to decode1a!
-                    write(*,*) '! Signal too strong, or suspect data?  Decoding aborted.'
-                    write(13,*) 'Signal too strong, or suspect data?  Decoding aborted.'
-                    call flush(13)
-                    go to 999
-                 endif
-                 call timer('decode1a',0)
-                 ifreq=i
-                 ikHz=nint(freq+0.5*(nfa+nfb)-foffset)-nfshift
-                 idf=nint(1000.0*(freq+0.5*(nfa+nfb)-foffset-(ikHz+nfshift)))
-                 call decode1a(dd,newdat,f00,nflip,mode65,nfsample,xpol,  &
-                      mycall,hiscall,hisgrid,neme,ndepth,nqd,dphi,        &
-                      nutc,ikHz,idf,ipol,sync2,a,dt,pol,nkv,nhist,qual,decoded)
-                 dt=dt+0.8                           !### empirical tweak
-                 call timer('decode1a',1)
-
-                 if(km.lt.MAXMSG) km=km+1
-                 sig(km,1)=nfile
-                 sig(km,2)=nutc
-                 sig(km,3)=freq + 0.5*(nfa+nfb)
-                 sig(km,4)=sync1
-                 sig(km,5)=dt
-                 sig(km,6)=pol
-                 sig(km,7)=flipk
-                 sig(km,8)=sync2
-                 sig(km,9)=nkv
-                 sig(km,10)=qual
-!                 sig(km,11)=idphi
-                 sig(km,12)=savg(ipol,i)
-                 sig(km,13)=a(1)
-                 sig(km,14)=a(2)
-                 sig(km,15)=a(3)
-                 sig(km,16)=a(4)
-!                     sig(km,17)=a(5)
-                 sig(km,18)=nhist
-                 msg(km)=decoded
-                 freq0=freq
-                 sync10=sync1
-                 nkm=1
-              endif
-           endif
-        endif
-!70      continue
-     enddo
-
-     if(nqd.eq.1) then
-        nwrite=0
-        do k=1,km
-           decoded=msg(k)
-           if(decoded.ne.'                      ') then
-              nutc=sig(k,2)
-              freq=sig(k,3)
-              sync1=sig(k,4)
-              dt=sig(k,5)
-              npol=nint(57.2957795*sig(k,6))
-              flip=sig(k,7)
-              sync2=sig(k,8)
-              nkv=sig(k,9)
-              nqual=sig(k,10)
-!              idphi=nint(sig(k,11))
-              if(flip.lt.0.0) then
-                 do i=22,1,-1
-                    if(decoded(i:i).ne.' ') go to 8
-                 enddo
-                 stop 'Error in message format'
-8                if(i.le.18) decoded(i+2:i+4)='OOO'
-              endif
-              nkHz=nint(freq-foffset)-nfshift
-              mhz=fcenter                         ! ... +fadd ???
-              f0=mhz+0.001*nkHz
-              ndf=nint(1000.0*(freq-foffset-(nkHz+nfshift)))
-              nsync1=sync1
-              nsync2=nint(10.0*log10(sync2)) - 40 !### empirical ###
-              if(decoded(1:4).eq.'RO  ' .or. decoded(1:4).eq.'RRR  ' .or.  &
-                   decoded(1:4).eq.'73  ') nsync2=nsync2-6
-              nwrite=nwrite+1
-              if(nxant.ne.0) then
-                 npol=npol-45
-                 if(npol.lt.0) npol=npol+180
-              endif
-
-!  If Tx station's grid is in decoded message, compute optimum TxPol
-              i1=index(decoded,' ')
-              i2=index(decoded(i1+1:),' ') + i1
-              grid='      '
-              if(i2.ge.8 .and. i2.le.18) grid=decoded(i2+1:i2+4)//'mm'
-              ntxpol=0
-              cp=' '
-              if(xpol) then
-                 if(grid(1:1).ge.'A' .and. grid(1:1).le.'R' .and.           &
-                      grid(2:2).ge.'A' .and. grid(2:2).le.'R' .and.         &
-                      grid(3:3).ge.'0' .and. grid(3:3).le.'9' .and.         &
-                      grid(4:4).ge.'0' .and. grid(4:4).le.'9') then                 
-                    ntxpol=mod(npol-nint(2.0*dpol(mygrid,grid))+720,180)
-                    if(nxant.eq.0) then
-                       cp='H'
-                       if(ntxpol.gt.45 .and. ntxpol.le.135) cp='V'
-                    else
-                       cp='/'
-                       if(ntxpol.ge.90 .and. ntxpol.lt.180) cp='\\'
-                    endif
-                 endif
-              endif
-              
-              if(ndphi.eq.0) then
-                 write(*,1010) nkHz,ndf,npol,nutc,dt,nsync2,    &
-                      decoded,nkv,nqual,ntxpol,cp
-1010             format('!',i3,i5,i4,i5.4,f5.1,i4,2x,a22,i4,i5,i5,1x,a1)
-              else
-                 if(iloop.ge.1) qphi(iloop)=sig(k,10)
-                 write(*,1010) nkHz,ndf,npol,nutc,dt,nsync2,    &
-                      decoded,nkv,nqual,30*iloop
-                 write(27,1011) 30*iloop,nkHz,ndf,npol,nutc,  &
-                      dt,sync2,nkv,nqual,decoded
-1011             format(i3,i4,i5,i4,i5.4,f5.1,f7.1,i3,i5,2x,a22)
-              endif
-           endif
-        enddo
-
-        if(nwrite.eq.0) then
-           write(*,1012) mousefqso,nutc
-1012       format('!',i3,9x,i5.4,' ')
-        endif
-   
-     endif
-     if(ndphi.eq.1 .and.iloop.lt.12) then
-        iloop=iloop+1
-        go to 2
-     endif
-     
-     if(ndphi.eq.1 .and.iloop.eq.12) call getdphi(qphi)
-     if(nagain.eq.1) go to 999
-  enddo
-
-!  Trim the list and produce a sorted index and sizes of groups.
-!  (Should trimlist remove all but best SNR for given UTC and message content?)
-  call trimlist(sig,km,ftol,indx,nsiz,nz)
-
-  do i=1,km
-     done(i)=.false.
-  enddo
-  j=0
-  ilatest=-1
-  do n=1,nz
-     ifile0=0
-     do m=1,nsiz(n)
-        i=indx(j+m)
-        ifile=sig(i,1)
-        if(ifile.gt.ifile0 .and.msg(i).ne.blank) then
-           ilatest=i
-           ifile0=ifile
-        endif
-     enddo
-     i=ilatest
-
-     if(i.ge.1) then
-        if(.not.done(i)) then
-           done(i)=.true.
-           nutc=sig(i,2)
-           freq=sig(i,3)
-           sync1=sig(i,4)
-           dt=sig(i,5)
-           npol=nint(57.2957795*sig(i,6))
-           flip=sig(i,7)
-           sync2=sig(i,8)
-           nkv=sig(i,9)
-           nqual=min(sig(i,10),10.0)
-!                  rms0=sig(i,11)
-           do k=1,5
-              a(k)=sig(i,12+k)
-           enddo
-           nhist=sig(i,18)
-           decoded=msg(i)
-           
-           if(flip.lt.0.0) then
-              do i=22,1,-1
-                 if(decoded(i:i).ne.' ') go to 10
-              enddo
-              stop 'Error in message format'
-10            if(i.le.18) decoded(i+2:i+4)='OOO'
-           endif
-           mhz=fcenter                             !... +fadd ???
-           nkHz=nint(freq-foffset)-nfshift
-           f0=mhz+0.001*nkHz
-           ndf=nint(1000.0*(freq-foffset-(nkHz+nfshift)))
-           ndf0=nint(a(1))
-           ndf1=nint(a(2))
-           ndf2=nint(a(3))
-           nsync1=sync1
-           nsync2=nint(10.0*log10(sync2)) - 40 !### empirical ###
-           if(decoded(1:4).eq.'RO  ' .or. decoded(1:4).eq.'RRR  ' .or.  &
-                decoded(1:4).eq.'73  ') nsync2=nsync2-6
-           if(nxant.ne.0) then
-              npol=npol-45
-              if(npol.lt.0) npol=npol+180
-           endif
-
-!  If Tx station's grid is in decoded message, compute optimum TxPol
-           i1=index(decoded,' ')
-           i2=index(decoded(i1+1:),' ') + i1
-           grid='      '
-           if(i2.ge.8 .and. i2.le.18) grid=decoded(i2+1:i2+4)//'mm'
-           ntxpol=0
-           cp=' '
-           if(xpol) then
-              if(grid(1:1).ge.'A' .and. grid(1:1).le.'R' .and.           &
-                   grid(2:2).ge.'A' .and. grid(2:2).le.'R' .and.         &
-                   grid(3:3).ge.'0' .and. grid(3:3).le.'9' .and.         &
-                   grid(4:4).ge.'0' .and. grid(4:4).le.'9') then                 
-                 ntxpol=mod(npol-nint(2.0*dpol(mygrid,grid))+720,180)
-                 if(nxant.eq.0) then
-                    cp='H'
-                    if(ntxpol.gt.45 .and. ntxpol.le.135) cp='V'
-                 else
-                    cp='/'
-                    if(ntxpol.ge.90 .and. ntxpol.lt.180) cp='\\'
-                 endif
-              endif
-           endif
-           write(26,1014) f0,ndf,ndf0,ndf1,ndf2,dt,npol,nsync1,       &
-                nsync2,nutc,decoded,nkv,nqual,nhist,cp
-           write(21,1014) f0,ndf,ndf0,ndf1,ndf2,dt,npol,nsync1,       &
-                nsync2,nutc,decoded,nkv,nqual,nhist
-1014       format(f8.3,i5,3i3,f5.1,i4,i3,i4,i5.4,2x,a22,3i3,1x,a1)
-
-        endif
-     endif
-     j=j+nsiz(n)
-  enddo
-  write(26,1015) nutc
-1015 format(39x,i4.4)
-  call flush(21)
-  call flush(26)
-  call display(nkeep,ftol)
-  ndecdone=2
-
-999 close(23)
-  ndphi=0
-  nagain=0
-  mcall3b=mcall3a
-
-  return
-end subroutine map65a
diff --git a/libm65/match.f90 b/libm65/match.f90
deleted file mode 100644
index e9570dabc..000000000
--- a/libm65/match.f90
+++ /dev/null
@@ -1,28 +0,0 @@
-subroutine match(s1,s2,nstart,nmatch)
-
-  character*(*) s1,s2
-  character s1a*29
-
-  nstart=-1
-  nmatch=0
-  n1=len_trim(s1)+1
-  n2=len(s2)
-  s1a=s1//' '
-  if(n2.ge.n1) then
-     do j=1,n2
-        n=0
-        do i=1,n1
-           k=j+i-1
-           if(k.gt.n2) k=k-n2
-           if(s2(k:k).eq.s1a(i:i)) n=n+1
-        enddo
-        if(n.gt.nmatch) then
-           nmatch=n
-           nstart=j
-        endif
-     enddo
-  endif
-
-  return
-end subroutine match
-        
diff --git a/libm65/moon2.f b/libm65/moon2.f
deleted file mode 100644
index 8144b675f..000000000
--- a/libm65/moon2.f
+++ /dev/null
@@ -1,167 +0,0 @@
-      subroutine moon2(y,m,Day,UT,lon,lat,RA,Dec,topRA,topDec,
-     +  LST,HA,Az,El,dist)
-
-      implicit none
-
-      integer y                           !Year
-      integer m                           !Month
-      integer Day                         !Day
-      real*8 UT                           !UTC in hours
-      real*8 RA,Dec                       !RA and Dec of moon
-
-C  NB: Double caps are single caps in the writeup.
-
-      real*8 NN                           !Longitude of ascending node
-      real*8 i                            !Inclination to the ecliptic
-      real*8 w                            !Argument of perigee
-      real*8 a                            !Semi-major axis
-      real*8 e                            !Eccentricity
-      real*8 MM                           !Mean anomaly
-
-      real*8 v                            !True anomaly
-      real*8 EE                           !Eccentric anomaly
-      real*8 ecl                          !Obliquity of the ecliptic
-
-      real*8 d                            !Ephemeris time argument in days
-      real*8 r                            !Distance to sun, AU
-      real*8 xv,yv                        !x and y coords in ecliptic
-      real*8 lonecl,latecl                !Ecliptic long and lat of moon
-      real*8 xg,yg,zg                     !Ecliptic rectangular coords
-      real*8 Ms                           !Mean anomaly of sun
-      real*8 ws                           !Argument of perihelion of sun
-      real*8 Ls                           !Mean longitude of sun (Ns=0)
-      real*8 Lm                           !Mean longitude of moon
-      real*8 DD                           !Mean elongation of moon
-      real*8 FF                           !Argument of latitude for moon
-      real*8 xe,ye,ze                     !Equatorial geocentric coords of moon
-      real*8 mpar                         !Parallax of moon (r_E / d)
-      real*8 lat,lon                      !Station coordinates on earth
-      real*8 gclat                        !Geocentric latitude
-      real*8 rho                          !Earth radius factor
-      real*8 GMST0,LST,HA
-      real*8 g
-      real*8 topRA,topDec                 !Topocentric coordinates of Moon
-      real*8 Az,El
-      real*8 dist
-
-      real*8 rad,twopi,pi,pio2
-      data rad/57.2957795131d0/,twopi/6.283185307d0/
-
-      d=367*y - 7*(y+(m+9)/12)/4 + 275*m/9 + Day - 730530 + UT/24.d0
-      ecl = 23.4393d0 - 3.563d-7 * d
-
-C  Orbital elements for Moon:  
-      NN = 125.1228d0 - 0.0529538083d0 * d
-      i = 5.1454d0
-      w = mod(318.0634d0 + 0.1643573223d0 * d + 360000.d0,360.d0)
-      a = 60.2666d0
-      e = 0.054900d0
-      MM = mod(115.3654d0 + 13.0649929509d0 * d + 360000.d0,360.d0)
-
-      EE = MM + e*rad*sin(MM/rad) * (1.d0 + e*cos(MM/rad))
-      EE = EE - (EE - e*rad*sin(EE/rad)-MM) / (1.d0 - e*cos(EE/rad))
-      EE = EE - (EE - e*rad*sin(EE/rad)-MM) / (1.d0 - e*cos(EE/rad))
-
-      xv = a * (cos(EE/rad) - e)
-      yv = a * (sqrt(1.d0-e*e) * sin(EE/rad))
-
-      v = mod(rad*atan2(yv,xv)+720.d0,360.d0)
-      r = sqrt(xv*xv + yv*yv)
-
-C  Get geocentric position in ecliptic rectangular coordinates:
-
-      xg = r * (cos(NN/rad)*cos((v+w)/rad) -
-     +          sin(NN/rad)*sin((v+w)/rad)*cos(i/rad))
-      yg = r * (sin(NN/rad)*cos((v+w)/rad) + 
-     +          cos(NN/rad)*sin((v+w)/rad)*cos(i/rad))
-      zg = r * (sin((v+w)/rad)*sin(i/rad))
-
-C  Ecliptic longitude and latitude of moon:
-      lonecl = mod(rad*atan2(yg/rad,xg/rad)+720.d0,360.d0)
-      latecl = rad*atan2(zg/rad,sqrt(xg*xg + yg*yg)/rad)
-
-C  Now include orbital perturbations:
-      Ms = mod(356.0470d0 + 0.9856002585d0 * d + 3600000.d0,360.d0)
-      ws = 282.9404d0 + 4.70935d-5*d
-      Ls = mod(Ms + ws + 720.d0,360.d0)
-      Lm = mod(MM + w + NN+720.d0,360.d0)
-      DD = mod(Lm - Ls + 360.d0,360.d0)
-      FF = mod(Lm - NN + 360.d0,360.d0)
-
-      lonecl = lonecl 
-     +  -1.274d0 * sin((MM-2.d0*DD)/rad)
-     +  +0.658d0 * sin(2.d0*DD/rad)
-     +  -0.186d0 * sin(Ms/rad)
-     +  -0.059d0 * sin((2.d0*MM-2.d0*DD)/rad)
-     +  -0.057d0 * sin((MM-2.d0*DD+Ms)/rad)
-     +  +0.053d0 * sin((MM+2.d0*DD)/rad)
-     +  +0.046d0 * sin((2.d0*DD-Ms)/rad)
-     +  +0.041d0 * sin((MM-Ms)/rad)
-     +  -0.035d0 * sin(DD/rad)
-     +  -0.031d0 * sin((MM+Ms)/rad)
-     +  -0.015d0 * sin((2.d0*FF-2.d0*DD)/rad)
-     +  +0.011d0 * sin((MM-4.d0*DD)/rad)
-
-      latecl = latecl
-     +  -0.173d0 * sin((FF-2.d0*DD)/rad)
-     +  -0.055d0 * sin((MM-FF-2.d0*DD)/rad)
-     +  -0.046d0 * sin((MM+FF-2.d0*DD)/rad)
-     +  +0.033d0 * sin((FF+2.d0*DD)/rad)
-     +  +0.017d0 * sin((2.d0*MM+FF)/rad)
-
-      r = 60.36298d0
-     +  - 3.27746d0*cos(MM/rad)
-     +  - 0.57994d0*cos((MM-2.d0*DD)/rad)
-     +  - 0.46357d0*cos(2.d0*DD/rad)
-     +  - 0.08904d0*cos(2.d0*MM/rad)
-     +  + 0.03865d0*cos((2.d0*MM-2.d0*DD)/rad)
-     +  - 0.03237d0*cos((2.d0*DD-Ms)/rad)
-     +  - 0.02688d0*cos((MM+2.d0*DD)/rad)
-     +  - 0.02358d0*cos((MM-2.d0*DD+Ms)/rad)
-     +  - 0.02030d0*cos((MM-Ms)/rad)
-     +  + 0.01719d0*cos(DD/rad)
-     +  + 0.01671d0*cos((MM+Ms)/rad)
-
-      dist=r*6378.140d0
-
-C  Geocentric coordinates:
-C  Rectangular ecliptic coordinates of the moon:
-
-      xg = r * cos(lonecl/rad)*cos(latecl/rad)
-      yg = r * sin(lonecl/rad)*cos(latecl/rad)
-      zg = r *                 sin(latecl/rad)
-
-C  Rectangular equatorial coordinates of the moon:
-      xe = xg
-      ye = yg*cos(ecl/rad) - zg*sin(ecl/rad)
-      ze = yg*sin(ecl/rad) + zg*cos(ecl/rad)
-   
-C  Right Ascension, Declination:
-      RA = mod(rad*atan2(ye,xe)+360.d0,360.d0)
-      Dec = rad*atan2(ze,sqrt(xe*xe + ye*ye))
-
-C  Now convert to topocentric system:
-      mpar=rad*asin(1.d0/r)
-C      alt_topoc = alt_geoc - mpar*cos(alt_geoc)
-      gclat = lat - 0.1924d0*sin(2.d0*lat/rad)
-      rho = 0.99883d0 + 0.00167d0*cos(2.d0*lat/rad)
-      GMST0 = (Ls + 180.d0)/15.d0
-      LST = mod(GMST0+UT+lon/15.d0+48.d0,24.d0)    !LST in hours
-      HA = 15.d0*LST - RA                          !HA in degrees
-      g = rad*atan(tan(gclat/rad)/cos(HA/rad))
-      topRA = RA - mpar*rho*cos(gclat/rad)*sin(HA/rad)/cos(Dec/rad)
-      topDec = Dec - mpar*rho*sin(gclat/rad)*sin((g-Dec)/rad)/sin(g/rad)
-
-      HA = 15.d0*LST - topRA                       !HA in degrees
-      if(HA.gt.180.d0) HA=HA-360.d0
-      if(HA.lt.-180.d0) HA=HA+360.d0
-
-      pi=0.5d0*twopi
-      pio2=0.5d0*pi
-      call dcoord(pi,pio2-lat/rad,0.d0,lat/rad,ha*twopi/360,
-     +  topDec/rad,az,el)
-      Az=az*rad
-      El=El*rad
-
-      return
-      end
diff --git a/libm65/moondop.f b/libm65/moondop.f
deleted file mode 100644
index ebac631c4..000000000
--- a/libm65/moondop.f
+++ /dev/null
@@ -1,73 +0,0 @@
-      subroutine MoonDop(nyear,month,nday,uth4,lon4,lat4,RAMoon4,
-     +  DecMoon4,LST4,HA4,AzMoon4,ElMoon4,vr4,dist4)
-
-      implicit real*8 (a-h,o-z)
-      real*4 uth4                    !UT in hours
-      real*4 lon4                    !West longitude, degrees
-      real*4 lat4                    !Latitude, degrees
-      real*4 RAMoon4                 !Topocentric RA of moon, hours
-      real*4 DecMoon4                !Topocentric Dec of Moon, degrees
-      real*4 LST4                    !Locat sidereal time, hours
-      real*4 HA4                     !Local Hour angle, degrees
-      real*4 AzMoon4                 !Topocentric Azimuth of moon, degrees
-      real*4 ElMoon4                 !Topocentric Elevation of moon, degrees
-      real*4 vr4                     !Radial velocity of moon wrt obs, km/s
-      real*4 dist4                   !Echo time, seconds
-
-      real*8 LST
-      real*8 RME(6)                  !Vector from Earth center to Moon
-      real*8 RAE(6)                  !Vector from Earth center to Obs
-      real*8 RMA(6)                  !Vector from Obs to Moon
-      real*8 pvsun(6)
-      real*8 rme0(6)
-      logical km,bary
-
-      data rad/57.2957795130823d0/,twopi/6.28310530717959d0/
-
-      km=.true.
-      dlat=lat4/rad
-      dlong1=lon4/rad
-      elev1=200.d0
-      call geocentric(dlat,elev1,dlat1,erad1)
-
-      dt=100.d0                       !For numerical derivative, in seconds
-      UT=uth4
-
-C  NB: geodetic latitude used here, but geocentric latitude used when 
-C  determining Earth-rotation contribution to Doppler.
-
-      call moon2(nyear,month,nDay,UT-dt/3600.d0,dlong1*rad,dlat*rad,
-     +  RA,Dec,topRA,topDec,LST,HA,Az0,El0,dist)
-      call toxyz(RA/rad,Dec/rad,dist,rme0)      !Convert to rectangular coords
-
-      call moon2(nyear,month,nDay,UT,dlong1*rad,dlat*rad,
-     +  RA,Dec,topRA,topDec,LST,HA,Az,El,dist)
-      call toxyz(RA/rad,Dec/rad,dist,rme)       !Convert to rectangular coords
-
-      phi=LST*twopi/24.d0
-      call toxyz(phi,dlat1,erad1,rae)           !Gencentric numbers used here!
-      radps=twopi/(86400.d0/1.002737909d0)
-      rae(4)=-rae(2)*radps                      !Vel of Obs wrt Earth center
-      rae(5)=rae(1)*radps
-      rae(6)=0.d0
-
-      do i=1,3
-         rme(i+3)=(rme(i)-rme0(i))/dt
-         rma(i)=rme(i)-rae(i)
-         rma(i+3)=rme(i+3)-rae(i+3)
-      enddo
-
-      call fromxyz(rma,alpha1,delta1,dtopo0)     !Get topocentric coords
-      vr=dot(rma(4),rma)/dtopo0
-
-      RAMoon4=topRA
-      DecMoon4=topDec
-      LST4=LST
-      HA4=HA
-      AzMoon4=Az
-      ElMoon4=El
-      vr4=vr
-      dist4=dist
-
-      return
-      end
diff --git a/libm65/mskdf.f90 b/libm65/mskdf.f90
deleted file mode 100644
index 2420052ae..000000000
--- a/libm65/mskdf.f90
+++ /dev/null
@@ -1,75 +0,0 @@
-subroutine mskdf(cdat,npts,t2,nfft1,f0,nfreeze,mousedf,ntol,dfx,snrsq2)
-
-! Determine DF for a JTMS signal.  Also find ferr, a measure of
-! frequency differerence between 1st and 2nd harmonic.  
-! (Should be 0.000)
-
-  parameter (NZ=128*1024)
-  complex cdat(npts)
-  integer ntol
-  real sq(NZ)
-  real ccf(-2600:2600)                  !Correct limits?
-  real tmp(NZ)
-  complex c(NZ)
-  data nsps/8/
-  save c
-
-  df1=48000.0/nfft1
-  nh=nfft1/2
-  fac=1.0/(nfft1**2)
-
-  do i=1,npts
-     c(i)=fac*cdat(i)**2
-  enddo
-  c(npts+1:nfft1)=0.
-  call four2a(c,nfft1,1,-1,1)
-
-! In the "doubled-frequencies" spectrum of squared cdat:
-  fa=2.0*(f0-400)
-  fb=2.0*(f0+400)
-  j0=nint(2.0*f0/df1)
-  ja=nint(fa/df1)
-  jb=nint(fb/df1)
-  jd=nfft1/nsps
-
-  do j=1,nh+1
-     sq(j)=real(c(j))**2 + aimag(c(j))**2
-!     if(j*df1.lt.6000.0) write(54,3009) j*df1,sq(j),db(sq(j))
-!3009 format(3f12.3)
-  enddo
-
-  ccf=0.
-  do j=ja,jb
-     ccf(j-j0-1)=sq(j)+sq(j+jd)
-  enddo
-
-  call pctile(ccf(ja-j0-1),tmp,jb-ja+1,50,base)
-  ccf=ccf/base
-
-  if(NFreeze.gt.0) then
-     fa=2.0*(f0+MouseDF-ntol)
-     fb=2.0*(f0+MouseDF+ntol)
-  endif  
-  ja=nint(fa/df1)
-  jb=nint(fb/df1)
-
-!  rewind 51
-  smax=0.
-  do j=ja,jb
-     k=j-j0-1
-     if(ccf(k).gt.smax) then
-        smax=ccf(k)
-        jpk=j
-     endif
-     f=0.5*k*df1
-!     write(51,3002) f,ccf(k)
-!3002 format(2f12.3)
-  enddo
-!  call flush(51)
-
-  fpk=(jpk-1)*df1  
-  dfx=0.5*fpk-f0
-  snrsq2=smax
-
-  return
-end subroutine mskdf
diff --git a/libm65/ping.f90 b/libm65/ping.f90
deleted file mode 100644
index 73e8ce23a..000000000
--- a/libm65/ping.f90
+++ /dev/null
@@ -1,42 +0,0 @@
-subroutine ping(s,nz,dtbuf,slim,wmin,pingdat,nping)
-
-! Detect pings and make note of their start time, duration, and peak strength.
-
-  real*4 s(nz)
-  real*4 pingdat(3,100)
-  logical inside
-
-  nping=0
-  peak=0.
-  inside=.false.
-  snrlim=10.0**(0.1*slim) - 1.0
-  sdown=10.0*log10(0.25*snrlim+1.0)
-
-  i0=0                                 !Silence compiler warnings.
-  tstart=0.0
-  do i=2,nz
-     if(s(i).ge.slim .and. .not.inside) then
-        i0=i
-        tstart=i0*dtbuf
-        inside=.true.
-        peak=0.
-     endif
-     if(inside .and. s(i).gt.peak) then
-        peak=s(i)
-     endif
-     if(inside .and. (s(i).lt.sdown .or. i.eq.nz)) then
-        if(i.gt.i0) then
-           if(dtbuf*(i-i0).ge.wmin) then
-              if(nping.le.99) nping=nping+1
-              pingdat(1,nping)=tstart
-              pingdat(2,nping)=dtbuf*(i-i0)
-              pingdat(3,nping)=peak
-           endif
-           inside=.false.
-           peak=0.
-        endif
-     endif
-  enddo
-
-  return
-end subroutine ping
diff --git a/libm65/recvpkt.f90 b/libm65/recvpkt.f90
deleted file mode 100644
index f55288d88..000000000
--- a/libm65/recvpkt.f90
+++ /dev/null
@@ -1,70 +0,0 @@
-subroutine recvpkt(nsam,nblock2,userx_no,k,buf4,buf8,buf16)
-
-! Reformat timf2 data from Linrad and stuff data into r*4 array dd().
-
-  parameter (NSMAX=60*96000)          !Total sample intervals per minute
-  parameter (NFFT=32768)
-  integer*1 userx_no
-  real*4 d4,buf4(*)                   !(348)
-  real*8 d8,buf8(*)                   !(174)
-  complex*16 c16,buf16(*)             !(87)
-  integer*2 jd(4),kd(2),nblock2
-  real*4 xd(4),yd(2)
-  real*8 fcenter
-  common/datcom/dd(4,5760000),ss(4,322,NFFT),savg(4,NFFT),fcenter,nutc,junk(34)
-  equivalence (kd,d4)
-  equivalence (jd,d8,yd)
-  equivalence (xd,c16)
-
-  if(nsam.eq.-1) then
-! Move data from the UDP packet buffer into array dd().
-     if(userx_no.eq.-1) then
-        do i=1,174                    !One RF channel, r*4 data
-           k=k+1
-           d8=buf8(i)
-           dd(1,k)=yd(1)
-           dd(2,k)=yd(2)
-        enddo
-     else if(userx_no.eq.1) then
-        do i=1,348                    !One RF channel, i*2 data
-           k=k+1
-           d4=buf4(i)
-           dd(1,k)=kd(1)
-           dd(2,k)=kd(2)
-        enddo
-     else if(userx_no.eq.-2) then
-        do i=1,87                    !Two RF channels, r*4 data
-           k=k+1
-           c16=buf16(i)
-           dd(1,k)=xd(1)
-           dd(2,k)=xd(2)
-           dd(3,k)=xd(3)
-           dd(4,k)=xd(4)
-        enddo
-     else if(userx_no.eq.2) then
-        do i=1,174                    !Two RF channels, i*2 data
-           k=k+1
-           d8=buf8(i)
-           dd(1,k)=jd(1)
-           dd(2,k)=jd(2)
-           dd(3,k)=jd(3)
-           dd(4,k)=jd(4)
-        enddo
-     endif
-  else
-     if(userx_no.eq.1) then
-        do i=1,nsam                    !One RF channel, r*4 data
-           k=k+1
-           d4=buf4(i)
-           dd(1,k)=kd(1)
-           dd(2,k)=kd(2)
-
-           k=k+1
-           dd(1,k)=kd(1)
-           dd(2,k)=kd(2)
-        enddo
-     endif
-  endif
-
-  return
-end subroutine recvpkt
diff --git a/libm65/rfile3a.f90 b/libm65/rfile3a.f90
deleted file mode 100644
index 7e2266513..000000000
--- a/libm65/rfile3a.f90
+++ /dev/null
@@ -1,14 +0,0 @@
-subroutine rfile3a(infile,ibuf,n,fcenter,ierr)
-
-  character*(*) infile
-  integer*8 ibuf(n)
-  real*8 fcenter
-
-  open(10,file=infile,access='stream',status='old',err=998)
-  read(10,end=998) (ibuf(i),i=1,n/8),fcenter
-  ierr=0
-  go to 999
-998 ierr=1002
-999 close(10)
-  return
-end subroutine rfile3a
diff --git a/libm65/rtping.f90 b/libm65/rtping.f90
deleted file mode 100644
index 886c9fd12..000000000
--- a/libm65/rtping.f90
+++ /dev/null
@@ -1,92 +0,0 @@
-subroutine rtping(dat,k,cfile6,MinSigdB,MouseDF,ntol,mycall)
-
-! Called from datasink() every 2048 sample intervals (approx 43 ms).
-! Detects pings (signal level above MinSigdB).  When a ping ends, its
-! MSK signal is synchronized and decoded.
-
-  parameter (NSMAX=30*48000)
-  parameter (NZMAX=703)          !703 = NSMAX/2048
-  real dat(NSMAX)                !Raw audio data
-  character*6 cfile6             !Time hhmmss at start of this Rx interval
-  character*12 mycall
-  real sig(NZMAX)                !Sq-law detected signal, sampled at 43 ms
-  real sigdb(NZMAX)              !Signal in dB, sampled at 43 ms
-  real tmp(NZMAX)
-  logical inside,pingFound
-  data k0/9999999/
-  save
-
-  if(k.lt.k0) then
-     inside=.false.
-     pingFound=.false.
-     j0=0
-     t1=0.
-     width=0.
-     peak=0.
-     tpk=0.
-     dt=1.0/48000.0
-     kstep=2048
-     sdt=dt*kstep
-     wmin=0.043
-  endif
-  k0=k
-
-  slim=MinSigdB
-  snrlim=10.0**(0.1*slim) - 1.0
-  sdown=10.0*log10(0.25*snrlim+1.0)
-
-! Find signal power
-  j=k/kstep
-  sig(j)=dot_product(dat(k-kstep+1:k),dat(k-kstep+1:k))/kstep
-  if(j.lt.20) return
-
-! Determine baseline noise level
-  if(mod(j,20).eq.0) call pctile (sig,tmp,j,50,base)
-  sigdb(j)=db(sig(j)/base)                             ! (S+N)/N in dB
-
-!  write(72,3001) j*sdt,base,sig(j),sigdb(j)
-!3001 format(f10.3,3f12.6)
-
-  if(sigdb(j).ge.slim .and. .not.inside) then
-     j0=j                                       !Mark the start of a ping
-     t1=j0*sdt
-     inside=.true.
-     peak=0.
-  endif
-
-  if(inside .and. sigdb(j).gt.peak) then
-     peak=sigdb(j)                              !Save peak strength
-     tpk=j*sdt                                  ! and time of peak
-  endif
-  
-  if(inside .and. (sigdb(j).lt.sdown .or. j.eq.NZMAX)) then
-     width=(j-j0)*sdt                           !Save ping width
-     if(width.ge.wmin) pingFound=.true.
-  endif
-  if(.not.pingFound) return
-
-! A ping was found!  Assemble a signal report:
-  nwidth=0
-  if(width.ge.0.04) nwidth=1
-  if(width.ge.0.12) nwidth=2
-  if(width.gt.1.00) nwidth=3
-  nstrength=6
-  if(peak.ge.11.0) nstrength=7
-  if(peak.ge.17.0) nstrength=8
-  if(peak.ge.23.0) nstrength=9
-  nrpt=10*nwidth + nstrength
-
-  mswidth=10*nint(100.0*width)
-  i1=(t1-0.02)/dt
-  if(i1.lt.1) i1=1
-  iz=nint((width+0.02)/dt) + 1
-  iz=min(iz,k+1-i1,2*48000)                   !Length of ping in samples
-
-  call jtmsk(dat(i1),iz,cfile6,tpk,mswidth,int(peak),nrpt,      &
-       nfreeze,DFTolerance,MouseDF,pick,mycall)
-
-  pingFound=.false.
-  inside=.false.
-
-  return
-end subroutine rtping
diff --git a/libm65/s3avg.f90 b/libm65/s3avg.f90
deleted file mode 100644
index c50927bf3..000000000
--- a/libm65/s3avg.f90
+++ /dev/null
@@ -1,42 +0,0 @@
-subroutine s3avg(nsave,mode65,nutc,ndf,xdt,npol,s3,nkv,decoded)
-
-  real s3(64,63),s3b(64,63)
-  real s3a(64,63,32)
-  integer iutc(32),idf(32),ipol(32)
-  real dt(32)
-  character*22 decoded
-  logical ltext
-  save
-  
-  n=nsave
-  iutc(n)=nutc
-  idf(n)=ndf
-  ipol(n)=npol
-  dt(n)=xdt
-  s3a(1:64,1:63,n)=s3
-
-  s3b=0.
-  nsum=0
-  idfdiff=100
-  dtdiff=0.2
-  do i=1,n
-     if(mod(iutc(i),2).ne.mod(nutc,2)) cycle
-     if(abs(ndf-idf(i)).gt.idfdiff) cycle
-     if(abs(xdt-dt(i)).gt.dtdiff) cycle
-     s3b=s3b + s3a(1:64,1:63,i)
-     nsum=nsum+1
-  enddo
-
-  decoded='                      '
-  if(nsum.ge.2) then
-     nadd=mode65*nsum
-     call extract(s3b,nadd,ncount,nhist,decoded,ltext)     !Extract the message
-     nkv=nsum
-     if(ncount.lt.0) then 
-        nkv=0
-        decoded='                      '
-     endif
-  endif
-
-  return
-end subroutine s3avg
diff --git a/libm65/setup65.f b/libm65/setup65.f
deleted file mode 100644
index 25b821e5a..000000000
--- a/libm65/setup65.f
+++ /dev/null
@@ -1,96 +0,0 @@
-      subroutine setup65
-
-C  Defines arrays related to the JT65 pseudo-random synchronizing pattern.
-C  Executed at program start.
-
-      integer nprc(126)
-      common/prcom/pr(126),mdat(126),mref(126,2),mdat2(126),mref2(126,2)
-
-C  JT65
-      data nprc/
-     + 1,0,0,1,1,0,0,0,1,1,1,1,1,1,0,1,0,1,0,0,
-     + 0,1,0,1,1,0,0,1,0,0,0,1,1,1,0,0,1,1,1,1,
-     + 0,1,1,0,1,1,1,1,0,0,0,1,1,0,1,0,1,0,1,1,
-     + 0,0,1,1,0,1,0,1,0,1,0,0,1,0,0,0,0,0,0,1,
-     + 1,0,0,0,0,0,0,0,1,1,0,1,0,0,1,0,1,1,0,1,
-     + 0,1,0,1,0,0,1,1,0,0,1,0,0,1,0,0,0,0,1,1,
-     + 1,1,1,1,1,1/
-      data mr2/0/                !Silence compiler warning
-
-C  Put the appropriate pseudo-random sequence into pr
-      nsym=126
-      do i=1,nsym
-         pr(i)=2*nprc(i)-1
-      enddo
-
-C  Determine locations of data and reference symbols
-      k=0
-      mr1=0
-      do i=1,nsym
-         if(pr(i).lt.0.0) then
-            k=k+1
-            mdat(k)=i
-         else
-            mr2=i
-            if(mr1.eq.0) mr1=i
-         endif
-      enddo
-      nsig=k
-
-C  Determine the reference symbols for each data symbol.
-      do k=1,nsig
-         m=mdat(k)
-         mref(k,1)=mr1
-         do n=1,10                     !Get ref symbol before data
-            if((m-n).gt.0) then
-               if (pr(m-n).gt.0.0) go to 10
-            endif
-         enddo
-         go to 12
- 10      mref(k,1)=m-n
- 12      mref(k,2)=mr2
-         do n=1,10                     !Get ref symbol after data
-            if((m+n).le.nsym) then
-               if (pr(m+n).gt.0.0) go to 20
-            endif
-         enddo
-         go to 22
- 20      mref(k,2)=m+n
- 22   enddo
-
-C  Now do it all again, using opposite logic on pr(i)
-      k=0
-      mr1=0
-      do i=1,nsym
-         if(pr(i).gt.0.0) then
-            k=k+1
-            mdat2(k)=i
-         else
-            mr2=i
-            if(mr1.eq.0) mr1=i
-         endif
-      enddo
-      nsig=k
-
-      do k=1,nsig
-         m=mdat2(k)
-         mref2(k,1)=mr1
-         do n=1,10
-            if((m-n).gt.0) then
-               if (pr(m-n).lt.0.0) go to 110
-            endif
-         enddo
-         go to 112
- 110      mref2(k,1)=m-n
- 112      mref2(k,2)=mr2
-         do n=1,10
-            if((m+n).le.nsym) then
-               if (pr(m+n).lt.0.0) go to 120
-            endif
-         enddo
-         go to 122
- 120     mref2(k,2)=m+n
- 122  enddo
-
-      return
-      end
diff --git a/libm65/setupmsk.f90 b/libm65/setupmsk.f90
deleted file mode 100644
index c61d2d92a..000000000
--- a/libm65/setupmsk.f90
+++ /dev/null
@@ -1,51 +0,0 @@
-subroutine setupmsk(cw,cwb)
-
-! Calculate the JTMS character waveforms.
-
-  complex cw(168,0:63)
-  complex cwb(168)
-  integer nb(7)
-!  real*8 twopi,dt,f0,f1
-  character cc*64
-!                    1         2         3         4         5         6
-!          0123456789012345678901234567890123456789012345678901234567890123
-  data cc/'0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ./?-    _                  @'/
-
-  nsps=24
-  twopi=8.d0*atan(1.d0)
-  dt=1.d0/48000.d0                     !Sample interval
-  f0=1000.d0
-  f1=2000.d0
-  dphi0=twopi*dt*f0
-  dphi1=twopi*dt*f1
-
-  do i=0,63
-     k=0
-     m=0
-     do n=5,0,-1                          !Each character gets 6+1 bits
-        k=k+1
-        nb(k)=iand(1,ishft(i,-n))
-        m=m+nb(k)
-     enddo
-     k=k+1
-     nb(k) = 1 - iand(m,1)                !Insert odd parity bit
-
-     phi=0.
-     j=0
-     do k=1,7                             !Generate the waveform
-        if(nb(k).eq.0) then
-           dphi=dphi0
-        else
-           dphi=dphi1
-        endif
-        do ii=1,nsps
-           j=j+1
-           phi=phi+dphi
-           cw(j,i)=cmplx(cos(phi),sin(phi))
-        enddo
-     enddo
-  enddo
-  cwb=cw(1:168,44)
-
-  return
-end subroutine setupmsk
diff --git a/libm65/specjtms.f90 b/libm65/specjtms.f90
deleted file mode 100644
index 60a0dce9a..000000000
--- a/libm65/specjtms.f90
+++ /dev/null
@@ -1,106 +0,0 @@
-subroutine specjtms(k,px,pxsmo,spk0,f0)
-
-! Compute noise level and 2D spectra, for GUI display.
-
-  parameter (NSMAX=30*48000)
-  parameter (MAXFFT=8192)
-  integer*2 id
-  real x(MAXFFT)
-  complex cx(MAXFFT),cx2(MAXFFT)
-  logical first
-  common/mscom/id(1440000),s1(215),s2(215),kin,ndiskdat,kline,nutc
-  common/jt8com/nn(1800*12000),ss(184,32768),savg(32768)
-  data first/.true./
-  save
- 
-  if(first) then
-     pi=4.0*atan(1.0)
-     twopi=2.0*pi
-     kstep=2048
-     first=.false.
-     sqsmo=0.
-  endif
-
-  if(k.lt.kstep .or. k.gt.NSMAX) return
-
-  t=k/48000.0
-  nfft=4096
-  df=48000.0/nfft
-  nh=nfft/2
-  ib=k
-  ia=k-kstep+1
-  i0=k-nfft+1
-  sq=0.
-  do i=ia,ib
-     d=id(i)
-     sq=sq + d*d
-  enddo
-  sq=sq/2048.0
-  sqsmo=0.95*sqsmo + 0.05*sq
-  rms=sqrt(sq)
-  px=db(sq) - 23.0
-  pxsmo=db(sqsmo) - 23.0
-
-  if(i0.gt.0) then
-     do i=i0,ia-1
-        d=id(i)
-        sq=sq + d*d
-     enddo
-  endif
-  sq0=sq
-!  write(13,1010) t,rms,sq,px,pxsmo
-!1010 format(5f12.3)
-  if(k.lt.nfft) return
-
-
-  x(1:nfft)=id(i0:ib)
-  fac=0.002/nfft
-  cx=fac*x
-  call four2a(cx,nfft,1,-1,1)                    !Forward c2c FFT
-
-  iz=nint(2500.0/df)
-
-  do i=1,iz                                      !Save spectrum for plot
-     s1(i)=real(cx(i+1))**2 + aimag(cx(i+1))**2
-  enddo
-
-  cx(1)=0.5*cx(1)
-  cx(nh+2:nfft)=0.
-  call four2a(cx,nfft,1,1,1)                     !Inverse c2c FFT
-
-  cx2(1:nfft)=cx(1:nfft)*cx(1:nfft)
-  cx2(nfft+1:)=0.0
-
-  nfft=8192
-  df=48000.0/nfft
-  
-  call four2a(cx2,nfft,1,-1,1)                   !Forward c2c FFT of cx2
-
-!  j0=nint(2.0*1428.57/df)
-  j0=nint(2.0*1500.0/df)
-  ja=j0-107
-  jb=j0+107
-  do j=ja,jb
-     s2(j-ja+1)=1.e-4*(real(cx2(j))**2 + aimag(cx2(j))**2)
-  enddo
-
-  spk0=0.
-  fac=(5e8/sq0)**2
-  s2=fac*s2
-  do j=1,215
-     if(s2(j).gt.spk0) then
-        spk0=s2(j)
-        f=(j+ja-1)*df
-        f0=0.5*(f-3000.0)
-        phi0=0.5*atan2(aimag(cx2(j)),real(cx2(j)))
-     endif
-  enddo
-
-  spk0=0.5*db(spk0) - 7.0
-  kline=k/2048
-!  write(14,3001) k/2048,spk0,f0,phi0
-!3001 format(i8,3f12.3)
-!  call flush(14)
-
-  return
-end subroutine specjtms
diff --git a/libm65/sun.f b/libm65/sun.f
deleted file mode 100644
index a3f326e7f..000000000
--- a/libm65/sun.f
+++ /dev/null
@@ -1,88 +0,0 @@
-      subroutine sun(y,m,DD,UT,lon,lat,RA,Dec,LST,Az,El,mjd,day)
-
-      implicit none
-
-      integer y                         !Year
-      integer m                         !Month
-      integer DD                        !Day
-      integer mjd                       !Modified Julian Date
-      real UT                           !UTC in hours
-      real RA,Dec                       !RA and Dec of sun
-
-C  NB: Double caps here are single caps in the writeup.
-
-C  Orbital elements of the Sun (also N=0, i=0, a=1):
-      real w                            !Argument of perihelion
-      real e                            !Eccentricity
-      real MM                           !Mean anomaly
-      real Ls                           !Mean longitude
-
-C  Other standard variables:
-      real v                            !True anomaly
-      real EE                           !Eccentric anomaly
-      real ecl                          !Obliquity of the ecliptic
-      real d                            !Ephemeris time argument in days
-      real r                            !Distance to sun, AU
-      real xv,yv                        !x and y coords in ecliptic
-      real lonsun                       !Ecliptic long and lat of sun
-C Ecliptic coords of sun (geocentric)
-      real xs,ys
-C Equatorial coords of sun (geocentric)
-      real xe,ye,ze
-      real lon,lat
-      real GMST0,LST,HA
-      real xx,yy,zz
-      real xhor,yhor,zhor
-      real Az,El
-
-      real day
-      real rad
-      data rad/57.2957795/
-
-C  Time in days, with Jan 0, 2000 equal to 0.0:
-      d=367*y - 7*(y+(m+9)/12)/4 + 275*m/9 + DD - 730530 + UT/24.0
-      mjd=d + 51543
-      ecl = 23.4393 - 3.563e-7 * d
-
-C  Compute updated orbital elements for Sun:
-      w = 282.9404 + 4.70935e-5 * d
-      e = 0.016709 - 1.151e-9 * d
-      MM = mod(356.0470d0 + 0.9856002585d0 * d + 360000.d0,360.d0)
-      Ls = mod(w+MM+720.0,360.0)
-
-      EE = MM + e*rad*sin(MM/rad) * (1.0 + e*cos(M/rad))
-      EE = EE - (EE - e*rad*sin(EE/rad)-MM) / (1.0 - e*cos(EE/rad))
-
-      xv = cos(EE/rad) - e
-      yv = sqrt(1.0-e*e) * sin(EE/rad)
-      v = rad*atan2(yv,xv)
-      r = sqrt(xv*xv + yv*yv)
-      lonsun = mod(v + w + 720.0,360.0)
-C  Ecliptic coordinates of sun (rectangular):
-      xs = r * cos(lonsun/rad)
-      ys = r * sin(lonsun/rad)
-
-C  Equatorial coordinates of sun (rectangular):
-      xe = xs
-      ye = ys * cos(ecl/rad)
-      ze = ys * sin(ecl/rad)
-
-C  RA and Dec in degrees:
-      RA = rad*atan2(ye,xe)
-      Dec = rad*atan2(ze,sqrt(xe*xe + ye*ye))
-
-      GMST0 = (Ls + 180.0)/15.0
-      LST = mod(GMST0+UT+lon/15.0+48.0,24.0)    !LST in hours
-      HA = 15.0*LST - RA                        !HA in degrees
-      xx = cos(HA/rad)*cos(Dec/rad)
-      yy = sin(HA/rad)*cos(Dec/rad)
-      zz =             sin(Dec/rad)
-      xhor = xx*sin(lat/rad) - zz*cos(lat/rad)
-      yhor = yy
-      zhor = xx*cos(lat/rad) + zz*sin(lat/rad)
-      Az = mod(rad*atan2(yhor,xhor) + 180.0 + 360.0,360.0)
-      El = rad*asin(zhor)
-      day=d-1.5
-
-      return
-      end
diff --git a/libm65/sync9.f90 b/libm65/sync9.f90
index b989c1035..3c85438bc 100644
--- a/libm65/sync9.f90
+++ b/libm65/sync9.f90
@@ -34,8 +34,6 @@ subroutine sync9(ss,tstep,f0a,df3,lagpk,fpk)
   enddo
 
   fpk=f0a + (npk-1)*df3
-  write(*,1010) lagpk,npk,fpk
-1010 format('lagpk:',i4,'   npk:',i6,'   fpk:',f8.2)
 
   do lag=-lagmax,lagmax
      sum=0.
@@ -43,10 +41,7 @@ subroutine sync9(ss,tstep,f0a,df3,lagpk,fpk)
         k=ii(i) + lag
         if(k.ge.1) sum=sum + ss(k,npk)
      enddo
-!     write(73,3000) lag,sum
-!3000 format(i3,f12.3)
   enddo
-  flush(73)
 
   return
 end subroutine sync9
diff --git a/libm65/syncmsk.f90 b/libm65/syncmsk.f90
deleted file mode 100644
index 52b51d178..000000000
--- a/libm65/syncmsk.f90
+++ /dev/null
@@ -1,38 +0,0 @@
-subroutine syncmsk(cdat,npts,cwb,r,i1)
-
-! Establish character sync within a JTMS ping.
-
-  complex cdat(npts)                    !Analytic signal
-  complex cwb(168)                      !Complex waveform for 'space'
-  real r(60000)
-  real tmp(60000)
-  integer hist(168),hmax(1)
-  complex z
-
-  r=0.
-  jz=npts-168+1
-  do j=1,jz
-     z=0.
-     ss=0.
-     do i=1,168
-        ss=ss + abs(cdat(i+j-1))          !Total power
-        z=z + cdat(i+j-1)*conjg(cwb(i))   !Signal matching <space>
-     enddo
-     r(j)=abs(z)/ss                       !Goodness-of-fit to <space>
-!     write(52,3001) j/168.0,r(j),cdat(j)
-!3001 format(4f12.3)
-  enddo
-
-  ncut=99.0*float(jz-10)/float(jz)
-  call pctile(r,tmp,jz,ncut,rlim)
-  hist=0
-  do j=1,jz
-     k=mod(j-1,168)+1
-     if(r(j).gt.rlim) hist(k)=hist(k)+1
-  enddo
-
-  hmax=maxloc(hist)
-  i1=hmax(1)
-
-  return
-end subroutine syncmsk
diff --git a/libm65/tm2.f90 b/libm65/tm2.f90
deleted file mode 100644
index bd8361187..000000000
--- a/libm65/tm2.f90
+++ /dev/null
@@ -1,14 +0,0 @@
-subroutine tm2(day,xlat4,xlon4,xl4,b4)
-
-  implicit real*8 (a-h,o-z)
-  parameter (RADS=0.0174532925199433d0)
-
-  real*4 day4,xlat4,xlon4,xl4,b4
-
-  glat=xlat4*RADS
-  glong=xlon4*RADS
-  call tmoonsub(day,glat,glong,el,rv,xl,b,pax)
-  xl4=xl
-  b4=b
-
-end subroutine tm2
diff --git a/libm65/tmoonsub.c b/libm65/tmoonsub.c
deleted file mode 100644
index 29ac28b49..000000000
--- a/libm65/tmoonsub.c
+++ /dev/null
@@ -1,514 +0,0 @@
-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>
-
-#define RADS 0.0174532925199433
-#define DEGS 57.2957795130823
-#define TPI 6.28318530717959
-#define PI 3.1415927
-
-/* ratio of     earth radius to astronomical unit */
-#define ER_OVER_AU 0.0000426352325194252
-
-/* all prototypes here */
-
-double getcoord(int coord);
-void getargs(int argc, char *argv[], int *y, int *m, double *tz, double *glong, double *glat);
-double range(double y);
-double rangerad(double y);
-double days(int y, int m, int dn, double hour);
-double days_(int *y, int *m, int *dn, double *hour);
-void moonpos(double, double *, double *, double *);
-void sunpos(double , double *, double *, double *);
-double moontransit(int y, int m, int d, double timezone, double glat, double glong, int *nt);
-double atan22(double y, double x);
-double epsilon(double d);
-void equatorial(double d, double *lon, double *lat, double *r);
-void ecliptic(double d, double *lon, double *lat, double *r);
-double gst(double d);
-void topo(double lst, double glat, double *alp, double *dec, double *r);
-double alt(double glat, double ha, double dec);
-void libration(double day, double lambda, double beta, double alpha, double *l, double *b, double *p);
-void illumination(double day, double lra, double ldec, double dr, double sra, double sdec, double *pabl, double *ill);
-int daysinmonth(int y, int m);
-int isleap(int y);
-void tmoonsub_(double *day, double *glat, double *glong, double *moonalt, 
-   double *mrv, double *l, double *b, double *paxis);
-
-static const char
-*usage = "  Usage: tmoon date[yyyymm] timz[+/-h.hh] long[+/-dddmm] lat[+/-ddmm]\n"
-            "example: tmoon 200009 0 -00155 5230\n";
-
-/*
-  getargs() gets the arguments from the command line, does some basic error
-  checking, and converts arguments into numerical form. Arguments are passed
-  back in pointers. Error messages print to stderr so re-direction of output
-  to file won't leave users blind. Error checking prints list of all errors
-  in a command line before quitting.
-*/
-void getargs(int argc, char *argv[], int *y, int *m, double *tz,
-             double *glong, double *glat) {
-
-  int date, latitude, longitude;
-  int mflag = 0, yflag = 0, longflag = 0, latflag = 0, tzflag = 0;
-  int longminflag = 0, latminflag = 0, dflag = 0;
-
-  /* if not right number of arguments, then print example command line */
-
-  if (argc !=5) {
-    fprintf(stderr, usage);
-    exit(EXIT_FAILURE);
-  }
-
-  date = atoi(argv[1]);
-  *y = date / 100;
-  *m = date - *y * 100;
-  *tz = (double) atof(argv[2]);
-  longitude = atoi(argv[3]);
-  latitude = atoi(argv[4]);
-  *glong = RADS * getcoord(longitude);
-  *glat = RADS * getcoord(latitude);
-
-  /* set a flag for each error found */
-
-  if (*m > 12 || *m < 1) mflag = 1;
-  if (*y > 2500) yflag = 1;
-  if (date < 150001) dflag = 1;
-  if (fabs((float) *glong) > 180 * RADS) longflag = 1;
-  if (abs(longitude) % 100 > 59) longminflag = 1;
-  if (fabs((float) *glat) > 90 * RADS) latflag = 1;
-  if (abs(latitude) % 100 > 59) latminflag = 1;
-  if (fabs((float) *tz) > 12) tzflag = 1;
-
-  /* print all the errors found */
-  
-  if (dflag == 1) {
-    fprintf(stderr, "date: dates must be in form yyyymm, gregorian, and later than 1500 AD\n");
-  }
-  if (yflag == 1) {
-    fprintf(stderr, "date: too far in future - accurate from 1500 to 2500\n");
-  }
-  if (mflag == 1) {
-    fprintf(stderr, "date: month must be in range 0 to 12, eg - August 2000 is entered as 200008\n");
-  }
-  if (tzflag == 1) {
-    fprintf(stderr, "timz: must be in range +/- 12 hours, eg -6 for Chicago\n");
-  }
-  if (longflag == 1) {
-    fprintf(stderr, "long: must be in range +/- 180 degrees\n");
-  }
-  if (longminflag == 1) {
-    fprintf(stderr, "long: last two digits are arcmin - max 59\n");
-  }
-  if (latflag == 1) {
-    fprintf(stderr, " lat: must be in range +/- 90 degrees\n");
-  }
-  if (latminflag == 1) {
-    fprintf(stderr, " lat: last two digits are arcmin - max 59\n");
-  }
-
-  /* quits if one or more flags set */
-
-  if (dflag + mflag + yflag + longflag + latflag + tzflag + longminflag + latminflag > 0) {
-    exit(EXIT_FAILURE);
-  }
-  
-}
-
-/*
-   returns coordinates in decimal degrees given the
-   coord as a ddmm value stored in an integer.
-*/
-double getcoord(int coord) {
-  int west = 1;
-  double glg, deg;
-  if (coord < 0) west = -1;
-  glg = fabs((double) coord/100);
-  deg = floor(glg);
-  glg = west* (deg + (glg - deg)*100 / 60);
-  return(glg);
-}
-
-/*
-  days() takes the year, month, day in the month and decimal hours
-  in the day and returns the number of days since J2000.0.
-  Assumes Gregorian calendar.
-*/
-double days(int y, int m, int d, double h) {
-  int a, b;
-  double day;
-  
-  /*
-    The lines below work from 1900 march to feb 2100
-    a = 367 * y - 7 * (y + (m + 9) / 12) / 4 + 275 * m / 9 + d;
-    day = (double)a - 730531.5 + hour / 24;
-  */
-
-  /*  These lines work for any Gregorian date since 0 AD */
-  if (m ==1 || m==2) {
-    m +=12;
-    y -= 1;
-  }
-  a = y / 100;
-  b = 2 - a + a/4;
-  day = floor(365.25*(y + 4716)) + floor(30.6001*(m + 1))
-    + d + b - 1524.5 - 2451545 + h/24;
-  return(day);
-}
-double days_(int *y0, int *m0, int *d0, double *h0) 
-{
-  return days(*y0,*m0,*d0,*h0);
-}
-
-/*
-Returns 1 if y a leap year, and 0 otherwise, according
-to the Gregorian calendar
-*/
-int isleap(int y) {
-  int a = 0;
-  if(y % 4 == 0) a = 1;
-  if(y % 100 == 0) a = 0;
-  if(y % 400 == 0) a = 1;
-  return(a);
-}
-
-/*
-Given the year and the month, function returns the
-number of days in the month. Valid for Gregorian
-calendar.
-*/
-int daysinmonth(int y, int m) {
-  int b = 31;
-  if(m == 2) {
-    if(isleap(y) == 1) b= 29;
-    else b = 28;
-  }
-  if(m == 4 || m == 6 || m == 9 || m == 11) b = 30;
-  return(b);
-}
-
-/*
-moonpos() takes days from J2000.0 and returns ecliptic coordinates
-of moon in the pointers. Note call by reference.
-This function is within a couple of arcminutes most of the time,
-and is truncated from the Meeus Ch45 series, themselves truncations of
-ELP-2000. Returns moon distance in earth radii.
-Terms have been written out explicitly rather than using the
-table based method as only a small number of terms is
-retained.
-*/
-void moonpos(double d, double *lambda, double *beta, double *rvec) {
-  double dl, dB, dR, L, D, M, M1, F, e, lm, bm, rm, t;
-
-  t = d / 36525;
-
-  L = range(218.3164591  + 481267.88134236  * t) * RADS;
-  D = range(297.8502042  + 445267.1115168  * t) * RADS;
-  M = range(357.5291092  + 35999.0502909  * t) * RADS;
-  M1 = range(134.9634114  + 477198.8676313  * t - .008997 * t * t) * RADS;
-  F = range(93.27209929999999  + 483202.0175273  * t - .0034029*t*t)*RADS;
-  e = 1 - .002516 * t;
-
-  dl =      6288774 * sin(M1);
-  dl +=     1274027 * sin(2 * D - M1);
-  dl +=      658314 * sin(2 * D);
-  dl +=      213618 * sin(2 * M1);
-  dl -=  e * 185116 * sin(M);
-  dl -=      114332 * sin(2 * F) ;
-  dl +=       58793 * sin(2 * D - 2 * M1);
-  dl +=   e * 57066 * sin(2 * D - M - M1) ;
-  dl +=       53322 * sin(2 * D + M1);
-  dl +=   e * 45758 * sin(2 * D - M);
-  dl -=   e * 40923 * sin(M - M1);
-  dl -=       34720 * sin(D) ;
-  dl -=   e * 30383 * sin(M + M1) ;
-  dl +=       15327 * sin(2 * D - 2 * F) ;
-  dl -=       12528 * sin(M1 + 2 * F);
-  dl +=       10980 * sin(M1 - 2 * F);
-  lm = rangerad(L + dl / 1000000 * RADS);
-
-  dB =   5128122 * sin(F);
-  dB +=   280602 * sin(M1 + F);
-  dB +=   277693 * sin(M1 - F);
-  dB +=   173237 * sin(2 * D - F);
-  dB +=    55413 * sin(2 * D - M1 + F);
-  dB +=    46271 * sin(2 * D - M1 - F);
-  dB +=    32573 * sin(2 * D + F);
-  dB +=    17198 * sin(2 * M1 + F);
-  dB +=     9266 * sin(2 * D + M1 - F);
-  dB +=     8822 * sin(2 * M1 - F);
-  dB += e * 8216 * sin(2 * D - M - F);
-  dB +=     4324 * sin(2 * D - 2 * M1 - F);
-  bm = dB / 1000000 * RADS;
-
-  dR =    -20905355 * cos(M1);
-  dR -=     3699111 * cos(2 * D - M1);
-  dR -=     2955968 * cos(2 * D);
-  dR -=      569925 * cos(2 * M1);
-  dR +=   e * 48888 * cos(M);
-  dR -=        3149 * cos(2 * F);
-  dR +=      246158 * cos(2 * D - 2 * M1);
-  dR -=  e * 152138 * cos(2 * D - M - M1) ;
-  dR -=      170733 * cos(2 * D + M1);
-  dR -=  e * 204586 * cos(2 * D - M);
-  dR -=  e * 129620 * cos(M - M1);
-  dR +=      108743 * cos(D);
-  dR +=  e * 104755 * cos(M + M1);
-  dR +=       79661 * cos(M1 - 2 * F);
-  rm = 385000.56  + dR / 1000;
-
-  *lambda = lm;
-  *beta = bm;
-  /* distance to Moon must be in Earth radii */
-  *rvec = rm / 6378.14;
-}
-
-/*
-topomoon() takes the local siderial time, the geographical
-latitude of the observer, and pointers to the geocentric
-equatorial coordinates. The function overwrites the geocentric
-coordinates with topocentric coordinates on a simple spherical
-earth model (no polar flattening). Expects Moon-Earth distance in
-Earth radii.    Formulas scavenged from Astronomical Almanac 'low
-precision formulae for Moon position' page D46.
-*/
-
-void topo(double lst, double glat, double *alp, double *dec, double *r) {
-  double x, y, z, r1;
-  x = *r * cos(*dec) * cos(*alp) - cos(glat) * cos(lst);
-  y = *r * cos(*dec) * sin(*alp) - cos(glat) * sin(lst);
-  z = *r * sin(*dec)  - sin(glat);
-  r1 = sqrt(x*x + y*y + z*z);
-  *alp = atan22(y, x);
-  *dec = asin(z / r1);
-  *r = r1;
-}
-
-/*
-moontransit() takes date, the time zone and geographic longitude
-of observer and returns the time (decimal hours) of lunar transit
-on that day if there is one, and sets the notransit flag if there
-isn't. See Explanatory Supplement to Astronomical Almanac
-section 9.32 and 9.31 for the method.
-*/
-
-double moontransit(int y, int m, int d, double tz, double glat, double glong, int *notransit) {
-  double hm, ht, ht1, lon, lat, rv, dnew, lst;
-  int itcount;
-
-  ht1 = 180 * RADS;
-  ht = 0;
-  itcount = 0;
-  *notransit = 0;
-  do {
-    ht = ht1;
-    itcount++;
-    dnew = days(y, m, d, ht * DEGS/15) - tz/24;
-    lst = gst(dnew) + glong;
-    /* find the topocentric Moon ra (hence hour angle) and dec */
-    moonpos(dnew, &lon, &lat, &rv);
-    equatorial(dnew, &lon, &lat, &rv);
-    topo(lst, glat, &lon, &lat, &rv);
-    hm = rangerad(lst -  lon);
-    ht1 = rangerad(ht - hm);
-    /* if no convergence, then no transit on that day */
-    if (itcount > 30) {
-      *notransit = 1;
-      break;
-    }
-  }
-  while (fabs(ht - ht1) > 0.04 * RADS);
-  return(ht1);
-}
-
-/*
-  Calculates the selenographic coordinates of either the sub Earth point
-  (optical libration) or the sub-solar point (selen. coords of centre of
-  bright hemisphere).  Based on Meeus chapter 51 but neglects physical
-  libration and nutation, with some simplification of the formulas.
-*/
-void libration(double day, double lambda, double beta, double alpha, double *l, double *b, double *p) {
-  double i, f, omega, w, y, x, a, t, eps;
-  t = day / 36525;
-  i = 1.54242 * RADS;
-  eps = epsilon(day);
-  f = range(93.2720993 + 483202.0175273 * t - .0034029 * t * t) * RADS;
-  omega = range(125.044555 - 1934.1361849 * t + .0020762 * t * t) * RADS;
-  w = lambda - omega;
-  y = sin(w) * cos(beta) * cos(i) - sin(beta) * sin(i);
-  x = cos(w) * cos(beta);
-  a = atan22(y, x);
-  *l = a - f;
-
-  /*  kludge to catch cases of 'round the back' angles  */
-  if (*l < -90 * RADS) *l += TPI;
-  if (*l > 90 * RADS)  *l -= TPI;
-  *b = asin(-sin(w) * cos(beta) * sin(i) - sin(beta) * cos(i));
-
-  /*  pa pole axis - not used for Sun stuff */
-  x = sin(i) * sin(omega);
-  y = sin(i) * cos(omega) * cos(eps) - cos(i) * sin(eps);
-  w = atan22(x, y);
-  *p = rangerad(asin(sqrt(x*x + y*y) * cos(alpha - w) / cos(*b)));
-}
-
-/*
-  Takes: days since J2000.0, eq coords Moon, ratio of moon to sun distance,
-  eq coords Sun
-  Returns: position angle of bright limb wrt NCP, percentage illumination
-  of Sun
-*/
-void illumination(double day, double lra, double ldec, double dr, double sra, double sdec, double *pabl, double *ill) {
-  double x, y, phi, i;
-  y = cos(sdec) * sin(sra - lra);
-  x = sin(sdec) * cos(ldec) - cos(sdec) * sin(ldec) * cos (sra - lra);
-  *pabl = atan22(y, x);
-  phi = acos(sin(sdec) * sin(ldec) + cos(sdec) * cos(ldec) * cos(sra-lra));
-  i = atan22(sin(phi) , (dr - cos(phi)));
-  *ill = 0.5*(1 + cos(i));
-}
-
-/*
-sunpos() takes days from J2000.0 and returns ecliptic longitude
-of Sun in the pointers. Latitude is zero at this level of precision,
-but pointer left in for consistency in number of arguments.
-This function is within 0.01 degree (1 arcmin) almost all the time
-for a century either side of J2000.0. This is from the 'low precision
-fomulas for the Sun' from C24 of Astronomical Alamanac
-*/
-void sunpos(double d, double *lambda, double *beta, double *rvec) {
-  double L, g, ls, bs, rs;
-
-  L = range(280.461 + .9856474 * d) * RADS;
-  g = range(357.528 + .9856003 * d) * RADS;
-  ls = L + (1.915 * sin(g) + .02 * sin(2 * g)) * RADS;
-  bs = 0;
-  rs = 1.00014 - .01671 * cos(g) - .00014 * cos(2 * g);
-  *lambda = ls;
-  *beta = bs;
-  *rvec = rs;
-}
-
-/*
-this routine returns the altitude given the days since J2000.0
-the hour angle and declination of the object and the latitude
-of the observer. Used to find the Sun's altitude to put a letter
-code on the transit time, and to find the Moon's altitude at
-transit just to make sure that the Moon is visible.
-*/
-double alt(double glat, double ha, double dec) {
-  return(asin(sin(dec) * sin(glat) + cos(dec) * cos(glat) * cos(ha)));
-}
-
-/* returns an angle in degrees in the range 0 to 360 */
-double range(double x) {
-  double a, b;
-  b = x / 360;
-  a = 360 * (b - floor(b));
-  if (a < 0)
-    a = 360 + a;
-  return(a);
-}
-
-/* returns an angle in rads in the range 0 to two pi */
-double rangerad(double x) {
-  double a, b;
-  b = x / TPI;
-  a = TPI * (b - floor(b));
-  if (a < 0)
-    a = TPI + a;
-  return(a);
-}
-
-/*
-gets the atan2 function returning angles in the right
-order and  range
-*/
-double atan22(double y, double x) {
-  double a;
-
-  a = atan2(y, x);
-  if (a < 0) a += TPI;
-  return(a);
-}
-
-/*
-returns mean obliquity of ecliptic in radians given days since
-J2000.0.
-*/
-double epsilon(double d) {
-  double t = d/ 36525;
-  return((23.4392911111111 - (t* (46.8150 + 0.00059*t)/3600)) *RADS);
-}
-
-/*
-replaces ecliptic coordinates with equatorial coordinates
-note: call by reference destroys original values
-R is unchanged.
-*/
-void equatorial(double d, double *lon, double *lat, double *r) {
-  double  eps, ceps, seps, l, b;
-
-  l = *lon;
-  b = * lat;
-  eps = epsilon(d);
-  ceps = cos(eps);
-  seps = sin(eps);
-  *lon = atan22(sin(l)*ceps - tan(b)*seps, cos(l));
-  *lat = asin(sin(b)*ceps + cos(b)*seps*sin(l));
-}
-
-/*
-replaces equatorial coordinates with ecliptic ones. Inverse
-of above, but used to find topocentric ecliptic coords.
-*/
-void ecliptic(double d, double *lon, double *lat, double *r) {
-  double  eps, ceps, seps, alp, dec;
-  alp = *lon;
-  dec = *lat;
-  eps = epsilon(d);
-  ceps = cos(eps);
-  seps = sin(eps);
-  *lon = atan22(sin(alp)*ceps + tan(dec)*seps, cos(alp));
-  *lat = asin(sin(dec)*ceps - cos(dec)*seps*sin(alp));
-}
-
-/*
-returns the siderial time at greenwich meridian as
-an angle in radians given the days since J2000.0
-*/
-double gst( double d) {
-  double t = d / 36525;
-  double theta;
-  theta = range(280.46061837 + 360.98564736629 * d + 0.000387933 * t * t);
-  return(theta * RADS);
-}
-
-void tmoonsub_(double *day, double *glat, double *glong, double *moonalt, 
-   double *mrv, double *l, double *b, double *paxis)
-{
-  double mlambda, mbeta;
-  double malpha, mdelta;
-  double lst, mhr;
-  double tlambda, tbeta, trv;
-
-  lst = gst(*day) + *glong;
-      
-  /* find Moon topocentric coordinates for libration calculations */
-
-  moonpos(*day, &mlambda, &mbeta, mrv);
-  malpha = mlambda;
-  mdelta = mbeta;
-  equatorial(*day, &malpha, &mdelta, mrv);
-  topo(lst, *glat, &malpha, &mdelta, mrv);
-  mhr = rangerad(lst - malpha);
-  *moonalt = alt(*glat, mhr, mdelta);
-      
-  /* Optical libration and Position angle of the Pole */
-
-  tlambda = malpha;
-  tbeta = mdelta;
-  trv = *mrv;
-  ecliptic(*day, &tlambda, &tbeta, &trv);
-  libration(*day, tlambda, tbeta, malpha,  l, b, paxis);
-}
diff --git a/libm65/toxyz.f b/libm65/toxyz.f
deleted file mode 100644
index 9f75d5de1..000000000
--- a/libm65/toxyz.f
+++ /dev/null
@@ -1,25 +0,0 @@
-      subroutine toxyz(alpha,delta,r,vec)
-
-      implicit real*8 (a-h,o-z)
-      real*8 vec(3)
-
-      vec(1)=r*cos(delta)*cos(alpha)
-      vec(2)=r*cos(delta)*sin(alpha)
-      vec(3)=r*sin(delta)
-
-      return
-      end
-
-      subroutine fromxyz(vec,alpha,delta,r)
-
-      implicit real*8 (a-h,o-z)
-      real*8 vec(3)
-      data twopi/6.283185307d0/
-
-      r=sqrt(vec(1)**2 + vec(2)**2 + vec(3)**2)
-      alpha=atan2(vec(2),vec(1))
-      if(alpha.lt.0.d0) alpha=alpha+twopi
-      delta=asin(vec(3)/r)
-
-      return
-      end
diff --git a/libm65/tweak1.f90 b/libm65/tweak1.f90
deleted file mode 100644
index f6799dae4..000000000
--- a/libm65/tweak1.f90
+++ /dev/null
@@ -1,24 +0,0 @@
-subroutine tweak1(ca,jz,f0,cb)
-
-! Shift frequency of analytic signal ca, with output to cb
-
-  complex ca(jz),cb(jz)
-  real*8 twopi
-  complex*16 w,wstep
-  data twopi/0.d0/
-  save twopi
-
-  if(twopi.eq.0.d0) twopi=8.d0*atan(1.d0)
-  w=1.d0
-  dphi=twopi*f0/11025.d0
-  wstep=cmplx(cos(dphi),sin(dphi))
-  x0=0.5*(jz+1)
-  s=2.0/jz
-  do i=1,jz
-     x=s*(i-x0)
-     w=w*wstep
-     cb(i)=w*ca(i)
-  enddo
-
-  return
-end subroutine tweak1
diff --git a/libm65/twkfreq.f b/libm65/twkfreq.f
deleted file mode 100644
index 364e17c02..000000000
--- a/libm65/twkfreq.f
+++ /dev/null
@@ -1,29 +0,0 @@
-      subroutine twkfreq(c4aa,c4bb,n5,a)
-
-      complex c4aa(n5)
-      complex c4bb(n5)
-      real a(5)
-      complex w,wstep
-      data twopi/6.283185307/
-
-C  Apply AFC corrections to the c4aa and c4bb data
-      w=1.0
-      wstep=1.0
-      x0=0.5*(n5+1)
-      s=2.0/n5
-      do i=1,n5
-         x=s*(i-x0)
-         if(mod(i,1000).eq.1) then
-            p2=1.5*x*x - 0.5
-!            p3=2.5*(x**3) - 1.5*x
-!            p4=4.375*(x**4) - 3.75*(x**2) + 0.375
-            dphi=(a(1) + x*a(2) + p2*a(3)) * (twopi/1378.125)
-            wstep=cmplx(cos(dphi),sin(dphi))
-         endif
-         w=w*wstep
-         c4aa(i)=w*c4aa(i)
-         c4bb(i)=w*c4bb(i)
-      enddo
-
-      return
-      end
diff --git a/libm65/wrapkarn.c b/libm65/wrapkarn.c
deleted file mode 100644
index 9e0a51caf..000000000
--- a/libm65/wrapkarn.c
+++ /dev/null
@@ -1,70 +0,0 @@
-#include <math.h>
-#include <stdio.h>
-#include <float.h>
-#include <limits.h>
-#include <stdlib.h>
-#include "rs.h"
-
-static void *rs;
-static int first=1;
-
-void rs_encode_(int *dgen, int *sent)
-// Encode JT65 data dgen[12], producing sent[63].
-{
-  int dat1[12];
-  int b[51];
-  int i;
-
-  if(first) {
-    // Initialize the JT65 codec
-    rs=init_rs_int(6,0x43,3,1,51,0);
-    first=0;
-  }
-
-  // Reverse data order for the Karn codec.
-  for(i=0; i<12; i++) {
-    dat1[i]=dgen[11-i];
-  }
-  // Compute the parity symbols
-  encode_rs_int(rs,dat1,b);
-
-  // Move parity symbols and data into sent[] array, in reverse order.
-  for (i = 0; i < 51; i++) sent[50-i] = b[i];
-  for (i = 0; i < 12; i++) sent[i+51] = dat1[11-i];
-}
-
-void rs_decode_(int *recd0, int *era0, int *numera0, int *decoded, int *nerr)
-// Decode JT65 received data recd0[63], producing decoded[12].
-// Erasures are indicated in era0[numera].  The number of corrected
-// errors is *nerr.  If the data are uncorrectable, *nerr=-1 is returned.
-{
-  int numera;
-  int i;
-  int era_pos[50];
-  int recd[63];
-
-  if(first) {
-    rs=init_rs_int(6,0x43,3,1,51,0);
-    first=0;
-  }
-
-  numera=*numera0;
-  for(i=0; i<12; i++) recd[i]=recd0[62-i];
-  for(i=0; i<51; i++) recd[12+i]=recd0[50-i];
-  if(numera) 
-    for(i=0; i<numera; i++) era_pos[i]=era0[i];
-  *nerr=decode_rs_int(rs,recd,era_pos,numera);
-  for(i=0; i<12; i++) decoded[i]=recd[11-i];
-}
-
-
-void rs_encode__(int *dgen, int *sent)
-{
-  rs_encode_(dgen, sent);
-}
-
-void rs_decode__(int *recd0, int *era0, int *numera0, int *decoded, int *nerr)
-{
-  rs_decode_(recd0, era0, numera0, decoded, nerr);
-}
-
diff --git a/wsjtx.pro b/wsjtx.pro
index c88a54c89..1b9987872 100644
--- a/wsjtx.pro
+++ b/wsjtx.pro
@@ -54,7 +54,7 @@ RC_FILE = wsjtx.rc
 unix {
 INCLUDEPATH += $$quote(/usr/include/qwt-qt4)
 LIBS += -lfftw3f /usr/lib/libgfortran.so.3
-LIBS += ../wsjtx/libm65/libm65.a
+LIBS += ../wsjtx/libm65/libjt9.a
 LIBS += /usr/lib/libqwt-qt4.so
 LIBS += -lportaudio
 #LIBS +- -lusb
@@ -62,7 +62,7 @@ LIBS += -lportaudio
 
 win32 {
 INCLUDEPATH += c:/qwt-6.0.1/include
-LIBS += ../wsjtx/libm65/libm65.a
+LIBS += ../wsjtx/libm65/libjt9.a
 LIBS += ../wsjtx/libfftw3f_win.a
 LIBS += ../QtSupport/palir-02.dll
 LIBS += libwsock32