code:
module crc
  use, intrinsic :: iso_c_binding, only: c_int, c_loc, c_int8_t, c_bool, c_short
  interface

    function crc14 (data, length) bind (C, name="crc14")
      use, intrinsic :: iso_c_binding, only: c_short, c_ptr, c_int
      implicit none
      integer (c_short) :: crc14
      type (c_ptr), value :: data
      integer (c_int), value :: length
    end function crc14

    function crc14_check (data, length) bind (C, name="crc16_check")
      use, intrinsic :: iso_c_binding, only: c_bool, c_ptr, c_int
      implicit none
      logical (c_bool) :: crc14_check
      type (c_ptr), value :: data
      integer (c_int), value :: length
    end function crc14_check

    function crc12 (data, length) bind (C, name="crc12")
      use, intrinsic :: iso_c_binding, only: c_short, c_ptr, c_int
      implicit none
      integer (c_short) :: crc12
      type (c_ptr), value :: data
      integer (c_int), value :: length
    end function crc12

    function crc12_check (data, length) bind (C, name="crc12_check")
      use, intrinsic :: iso_c_binding, only: c_bool, c_ptr, c_int
      implicit none
      logical (c_bool) :: crc12_check
      type (c_ptr), value :: data
      integer (c_int), value :: length
    end function crc12_check

    function crc10 (data, length) bind (C, name="crc10")
      use, intrinsic :: iso_c_binding, only: c_short, c_ptr, c_int
      implicit none
      integer (c_short) :: crc10
      type (c_ptr), value :: data
      integer (c_int), value :: length
    end function crc10

    function crc10_check (data, length) bind (C, name="crc10_check")
      use, intrinsic :: iso_c_binding, only: c_bool, c_ptr, c_int
      implicit none
      logical (c_bool) :: crc10_check
      type (c_ptr), value :: data
      integer (c_int), value :: length
    end function crc10_check

  end interface
end module crc


code:
#include <boost/crc.hpp>
#include <boost/config.hpp>

extern "C"
{
   short crc10 (unsigned char const * data, int length);
   bool crc10_check (unsigned char const * data, int length);
}

#define POLY 0x08f

#ifdef BOOST_NO_CXX11_CONSTEXPR
#define TRUNCATED_POLYNOMIAL POLY
#else
namespace
{
  unsigned long constexpr TRUNCATED_POLYNOMIAL = POLY;
}
#endif

// assumes CRC is last 16 bits of the data and is set to zero
// caller should assign the returned CRC into the message in big endian byte order
short crc10 (unsigned char const * data, int length)
{
    return boost::augmented_crc<10, TRUNCATED_POLYNOMIAL> (data, length);
}

bool crc10_check (unsigned char const * data, int length)
{
   return !boost::augmented_crc<10, TRUNCATED_POLYNOMIAL> (data, length);
}

license: mit


code:
#include <boost/crc.hpp>
#include <boost/config.hpp>

extern "C"
{
   short crc12 (unsigned char const * data, int length);
   bool crc12_check (unsigned char const * data, int length);
}

#define POLY 0xc06

#ifdef BOOST_NO_CXX11_CONSTEXPR
#define TRUNCATED_POLYNOMIAL POLY
#else
namespace
{
  unsigned long constexpr TRUNCATED_POLYNOMIAL = POLY;
}
#endif

// assumes CRC is last 16 bits of the data and is set to zero
// caller should assign the returned CRC into the message in big endian byte order
short crc12 (unsigned char const * data, int length)
{
    return boost::augmented_crc<12, TRUNCATED_POLYNOMIAL> (data, length);
}

bool crc12_check (unsigned char const * data, int length)
{
   return !boost::augmented_crc<12, TRUNCATED_POLYNOMIAL> (data, length);
}


code:
#include <boost/crc.hpp>
#include <boost/config.hpp>

extern "C"
{
   short crc14 (unsigned char const * data, int length);
   bool crc14_check (unsigned char const * data, int length);
}

#define POLY 0x2757

#ifdef BOOST_NO_CXX11_CONSTEXPR
#define TRUNCATED_POLYNOMIAL POLY
#else
namespace
{
  unsigned long constexpr TRUNCATED_POLYNOMIAL = POLY;
}
#endif

// assumes CRC is last 14 bits of the data and is set to zero
// caller should assign the returned CRC into the message in big endian byte order
short crc14 (unsigned char const * data, int length)
{
    return boost::augmented_crc<14, TRUNCATED_POLYNOMIAL> (data, length);
}

bool crc14_check (unsigned char const * data, int length)
{
   return !boost::augmented_crc<14, TRUNCATED_POLYNOMIAL> (data, length);
}


code:
real function db(x)
  db=-99.0
  if(x.gt.1.259e-10) db=10.0*log10(x)
  return
end function db


code:
subroutine decode4(dat,npts,dtx,nfreq,flip,mode4,ndepth,neme,minw,           &
     mycall,hiscall,hisgrid,decoded,nfano,deepbest,qbest,ichbest)

! Decodes JT4 data, assuming that DT and DF have already been determined.
! Input dat(npts) has already been downsampled by 2: rate = 11025/2.
! ### NB: this initial downsampling should be removed in WSJT-X, since
! it restricts the useful bandwidth to < 2.7 kHz.

  use jt4
  real dat(npts)                        !Raw data
  character decoded*22,deepmsg*22,deepbest*22
  character*12 mycall,hiscall
  character*6 hisgrid
  real*8 dt,df,phi,f0,dphi,twopi,phi1,dphi1
  complex*16 cz,cz1,c0,c1
  real*4 sym(207)

  twopi=8*atan(1.d0)
  dt=2.d0/11025             !Sample interval (2x downsampled data)
  df=11025.d0/2520.d0       !Tone separation for JT4A mode
  nsym=206
  amp=15.0
  istart=nint((dtx+0.8)/dt)              !Start index for synced FFTs
  if(istart.lt.0) istart=0
  nchips=0
  qbest=0.
  qtop=0.
  deepmsg='                      '
  ichbest=-1
  c0=0.
  k=istart
  phi=0.d0
  phi1=0.d0

  ich1=minw+1
  do ich=1,7
     if(nch(ich).le.mode4) ich2=ich
  enddo

  do ich=ich1,ich2
     nchips=min(nch(ich),70)
     nspchip=1260/nchips
     k=istart
     phi=0.d0
     phi1=0.d0
     fac2=1.e-8 * sqrt(float(mode4))
     do j=1,nsym+1
        if(flip.gt.0.0) then
           f0=nfreq + (npr(j))*mode4*df
           f1=nfreq + (2+npr(j))*mode4*df
        else
           f0=nfreq + (1-npr(j))*mode4*df
           f1=nfreq + (3-npr(j))*mode4*df
        endif
        dphi=twopi*dt*f0
        dphi1=twopi*dt*f1
        sq0=0.
        sq1=0.
        do nc=1,nchips
           phi=0.d0
           phi1=0.d0
           c0=0.
           c1=0.
           do i=1,nspchip
              k=k+1
              phi=phi+dphi
              phi1=phi1+dphi1
              cz=dcmplx(cos(phi),-sin(phi))
              cz1=dcmplx(cos(phi1),-sin(phi1))
              if(k.le.npts) then
                 c0=c0 + dat(k)*cz
                 c1=c1 + dat(k)*cz1
              endif
           enddo
           sq0=sq0 + real(c0)**2 + aimag(c0)**2
           sq1=sq1 + real(c1)**2 + aimag(c1)**2
        enddo
        sq0=fac2*sq0
        sq1=fac2*sq1
        rsym=amp*(sq1-sq0)
        if(j.ge.1) then
           rsymbol(j,ich)=rsym
           sym(j)=rsym
        endif
     enddo

     call extract4(sym,ncount,decoded)          !Do the convolutional decode
     nfano=0
     if(ncount.ge.0) then
        nfano=1
        ichbest=ich
        exit
     endif

     qual=0.                                    !Now try deep search
!     if(ndepth.ge.1) then
     if(iand(ndepth,32).eq.32) then
        call deep4(sym(2),neme,flip,mycall,hiscall,hisgrid,deepmsg,qual)
        if(qual.gt.qbest) then
           qbest=qual
           deepbest=deepmsg
           ichbest=ich
        endif
     endif
  enddo
  if(qbest.gt.qtop) then
     qtop=qbest
  endif
  qual=qbest

  return
end subroutine decode4


code:
subroutine decode65a(dd,npts,newdat,nqd,f0,nflip,mode65,ntrials,     &
     naggressive,ndepth,ntol,mycall,hiscall,hisgrid,nQSOProgress,    &
     ljt65apon, nexp_decode,  &
     bVHF,sync2,a,dt,nft,nspecial,qual,nhist,nsmo,decoded)

! Apply AFC corrections to a candidate JT65 signal, then decode it.

  use jt65_mod
  use timer_module, only: timer

  parameter (NMAX=60*12000)          !Samples per 60 s
  real*4  dd(NMAX)                   !92 MB: raw data from Linrad timf2
  complex cx(NMAX/8)                 !Data at 1378.125 sps
  complex cx1(NMAX/8)                !Data at 1378.125 sps, offset by 355.3 Hz
  complex c5x(NMAX/32)               !Data at 344.53125 Hz
  complex c5a(512)
  real s2(66,126)
  real a(5)
  logical bVHF,first,ljt65apon
  character decoded*22,decoded_best*22
  character mycall*12,hiscall*12,hisgrid*6
  character*27 cr
  data first/.true./,jjjmin/1000/,jjjmax/-1000/,cr/'(C) 2016, Joe Taylor - K1JT'/
  save

! Mix sync tone to baseband, low-pass filter, downsample to 1378.125 Hz
  call timer('filbig  ',0)
  call filbig(dd,npts,f0,newdat,cx,n5,sq0)
  if(mode65.eq.4) call filbig(dd,npts,f0+355.297852,newdat,cx1,n5,sq0)
  call timer('filbig  ',1)
! NB: cx has sample rate 12000*77125/672000 = 1378.125 Hz

! Check for a shorthand message
  if(bVHF .and. mode65.ne.101) then
     call sh65(cx,n5,mode65,ntol,xdf,nspecial,sync2)
     if(nspecial.gt.0) then
        a=0.
        a(1)=xdf
        nflip=0
     endif
  endif
  if(nflip.eq.0) go to 900

! Find best DF, drift, curvature, and DT.  Start by downsampling to 344.53125 Hz
  call timer('fil6521 ',0)
  call fil6521(cx,n5,c5x,n6)
  call timer('fil6521 ',1)

  fsample=1378.125/4.

  call timer('afc65b  ',0)
! Best fit for DF, drift, and dt. fsample = 344.53125 S/s
  dtbest=dt
  call afc65b(c5x,n6,fsample,nflip,mode65,a,ccfbest,dtbest)
  call timer('afc65b  ',1)
  dtbest=dtbest+0.003628 !Remove decimation filter and coh. integrator delay
  dt=dtbest              !Return new, improved estimate of dt
  sync2=3.7e-4*ccfbest/sq0                    !Constant is empirical
  if(mode65.eq.4) cx=cx1

! 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.
  a(3)=0
  call timer('twkfreq ',0)
  call twkfreq65(cx,n5,a)
  call timer('twkfreq ',1)

! Compute spectrum for each symbol.
  nsym=126
  nfft=512
  df=1378.125/nfft
  j=int(dtbest*1378.125)

  call timer('sh_ffts ',0)
  c5a=cmplx(0.0,0.0)
  do k=1,nsym
     do i=1,nfft
        j=j+1
        if(j.ge.1 .and. j.le.NMAX/8) then
           c5a(i)=cx(j)
        else
           c5a(i)=0.
        endif
     enddo
     call four2a(c5a,nfft,1,1,1)
     do i=1,512
        jj=i
        if(i.gt.256) jj=i-512
        s1(jj,k)=real(c5a(i))**2 + aimag(c5a(i))**2
     enddo
  enddo
  call timer('sh_ffts ',1)

  call timer('dec65b  ',0)
  qualbest=0.
  qual0=-1.e30
  minsmo=0
  maxsmo=0
  if(mode65.ge.2 .and. mode65.ne.101) then
     minsmo=nint(width/df)-1
     maxsmo=2*nint(width/df)
  endif
  nn=0
  do ismo=minsmo,maxsmo
     if(ismo.gt.0) then
        do j=1,126
           call smo121(s1(-255,j),512)
           if(j.eq.1) nn=nn+1
           if(nn.ge.4) then
              call smo121(s1(-255,j),512)
              if(j.eq.1) nn=nn+1
           endif
        enddo
     endif

     do i=1,66
        jj=i
        if(mode65.eq.2) jj=2*i-1
        if(mode65.eq.4) then
           ff=4*(i-1)*df - 355.297852
           jj=nint(ff/df)+1
        endif
        s2(i,1:126)=s1(jj,1:126)
     enddo
     nadd=ismo  !### ??? ###
     call decode65b(s2,nflip,nadd,mode65,ntrials,naggressive,ndepth,        &
          mycall,hiscall,hisgrid,nQSOProgress,ljt65apon,nexp_decode,        &
          nqd,nft,qual,     &
          nhist,decoded)

     if(nft.eq.1) then
        nsmo=ismo
        param(9)=nsmo
        nsum=1
        exit
     else if(nft.eq.2) then
        if(qual.gt.qualbest) then
           decoded_best=decoded
           qualbest=qual
           nnbest=nn
           nsmobest=ismo
        endif
     endif
     if(qual.lt.qual0) exit
     qual0=qual
  enddo

  if(nft.eq.2) then
     decoded=decoded_best
     qual=qualbest
     nsmo=nsmobest
     param(9)=nsmo
     nn=nnbest
  endif

  call timer('dec65b  ',1)

900 return
end subroutine decode65a


code:
subroutine decode65b(s2,nflip,nadd,mode65,ntrials,naggressive,ndepth,      &
     mycall,hiscall,hisgrid,nQSOProgress,ljt65apon,nexp_decode,nqd,        &
     nft,qual,         &
     nhist,decoded)

  use jt65_mod
  real s2(66,126)
  real s3(64,63)
  logical ltext,ljt65apon
  character decoded*22
  character mycall*12,hiscall*12,hisgrid*6
  save

  if(nqd.eq.-99) stop                !Silence compiler warning
  do j=1,63
     k=mdat(j)                       !Points to data symbol
     if(nflip.lt.0) k=mdat2(j)
     do i=1,64
        s3(i,j)=s2(i+2,k)
     enddo
  enddo

  call extract(s3,nadd,mode65,ntrials,naggressive,ndepth,nflip,mycall,   &
      hiscall,hisgrid,nQSOProgress,ljt65apon,nexp_decode,ncount,         &
      nhist,decoded,ltext,nft,qual)

! Suppress "birdie messages" and other garbage decodes:
  if(decoded(1:7).eq.'000AAA ') ncount=-1
  if(decoded(1:7).eq.'0L6MWK ') ncount=-1
  if(nflip.lt.0 .and. ltext) ncount=-1
  if(ncount.lt.0) then
     nft=0
     decoded='                      '
  endif

  return
end subroutine decode65b


code:
subroutine decode9w(nfqso,ntol,nsubmode,ss,id2,sync,nsnr,xdt1,f0,decoded)

! Decode a weak signal in a wide/slow JT9 submode.

  parameter (NSMAX=6827,NZMAX=60*12000)
  real ss(184,NSMAX)                     !Symbol spectra at 1/2-symbol steps
  real ccfred(NSMAX)                     !Best sync vs frequency
  real ccfblue(-9:18)                    !Sync vs time at best frequency
  real a(5)                              !Fitted Lorentzian params
  integer*2 id2(NZMAX)                   !Raw 16-bit data
  integer*1 i1SoftSymbols(207)           !Binary soft symbols
  character*22 decoded                   !Decoded message

  df=12000.0/16384.0                     !Bin spacing in ss()
  nsps=6912                              !Samples per 9-FSK symbol
  tstep=nsps*0.5/12000.0                 !Half-symbol duration
  npts=52*12000
  limit=10000                            !Fano timeout parameter

  ia=max(1,nint((nfqso-ntol)/df))        !Start frequency bin
  ib=min(NSMAX,nint((nfqso+ntol)/df))    !End frequency bin
  lag1=-int(2.5/tstep + 0.9999)          !Start lag
  lag2=int(5.0/tstep + 0.9999)           !End lag
  nhsym=184                              !Number of half-symbols

! First sync pass finds approximate Doppler spread; second pass does a
! good Lorentzian fit to determine frequency f0.
  do iter=1,2
     nadd=3
     if(iter.eq.2) nadd=2*nint(0.375*a(4)) + 1
     call sync9w(ss,nhsym,lag1,lag2,ia,ib,ccfred,ccfblue,ipk,lagpk,nadd)
     s=0.
     sq=0.
     ns=0
     do i=-9,18
        if(abs(i-lagpk).gt.3) then
           s=s+ccfblue(i)
           sq=sq+ccfblue(i)**2
           ns=ns+1
        endif
     enddo
     base=s/ns
     rms=sqrt(sq/ns - base**2)
     sync=(ccfblue(lagpk)-base)/rms
     xdt0=lagpk*tstep
     call lorentzian(ccfred(ia),ib-ia+1,a)
     f0=(ia+a(3))*df
  enddo
  ccfblue=(ccfblue-base)/rms

  call softsym9w(id2,npts,xdt0,f0,a(4)*df,nsubmode,xdt1-1.05,snrdb,i1softsymbols)
  nsnr=nint(snrdb)
  call jt9fano(i1softsymbols,limit,nlim,decoded)

!###
!  do i=-9,18
!     write(81,3081) i,ccfblue(i)
!3081 format(i3,f10.3)
!  enddo
!  do i=1,NSMAX
!     write(82,3082) i*df,ccfred(i)
!3082 format(f10.1,e12.3)
!  enddo
!###

  return
end subroutine decode9w


code:
/* Reed-Solomon decoder
 * Copyright 2002 Phil Karn, KA9Q
 * May be used under the terms of the GNU General Public License (GPL)
 */

#ifdef DEBUG
#include <stdio.h>
#endif

#include <stdlib.h>
#include <string.h>

#define	min(a,b)	((a) < (b) ? (a) : (b))

#ifdef FIXED
#include "fixed.h"
#elif defined(BIGSYM)
#include "int.h"
#else
#include "char.h"
#endif

int DECODE_RS(
#ifdef FIXED
DTYPE *data, int *eras_pos, int no_eras,int pad){
#else
void *p,DTYPE *data, int *eras_pos, int no_eras){
  struct rs *rs = (struct rs *)p;
#endif
  int deg_lambda, el, deg_omega;
  int i, j, r,k;
  DTYPE u,q,tmp,num1,num2,den,discr_r;
  DTYPE lambda[NROOTS+1], s[NROOTS];	/* Err+Eras Locator poly
					 * and syndrome poly */
  DTYPE b[NROOTS+1], t[NROOTS+1], omega[NROOTS+1];
  DTYPE root[NROOTS], reg[NROOTS+1], loc[NROOTS];
  int syn_error, count;

#ifdef FIXED
  /* Check pad parameter for validity */
  if(pad < 0 || pad >= NN)
    return -1;
#endif

  /* form the syndromes; i.e., evaluate data(x) at roots of g(x) */
  for(i=0;i<NROOTS;i++)
    s[i] = data[0];

  for(j=1;j<NN-PAD;j++){
    for(i=0;i<NROOTS;i++){
      if(s[i] == 0){
	s[i] = data[j];
      } else {
	s[i] = data[j] ^ ALPHA_TO[MODNN(INDEX_OF[s[i]] + (FCR+i)*PRIM)];
      }
    }
  }

  /* Convert syndromes to index form, checking for nonzero condition */
  syn_error = 0;
  for(i=0;i<NROOTS;i++){
    syn_error |= s[i];
    s[i] = INDEX_OF[s[i]];
  }

  if (!syn_error) {
    /* if syndrome is zero, data[] is a codeword and there are no
     * errors to correct. So return data[] unmodified
     */
    count = 0;
    goto finish;
  }
  memset(&lambda[1],0,NROOTS*sizeof(lambda[0]));
  lambda[0] = 1;

  if (no_eras > 0) {
    /* Init lambda to be the erasure locator polynomial */
    lambda[1] = ALPHA_TO[MODNN(PRIM*(NN-1-eras_pos[0]))];
    for (i = 1; i < no_eras; i++) {
      u = MODNN(PRIM*(NN-1-eras_pos[i]));
      for (j = i+1; j > 0; j--) {
	tmp = INDEX_OF[lambda[j - 1]];
	if(tmp != A0)
	  lambda[j] ^= ALPHA_TO[MODNN(u + tmp)];
      }
    }

#if DEBUG >= 1
    /* Test code that verifies the erasure locator polynomial just constructed
       Needed only for decoder debugging. */

    /* find roots of the erasure location polynomial */
    for(i=1;i<=no_eras;i++)
      reg[i] = INDEX_OF[lambda[i]];

    count = 0;
    for (i = 1,k=IPRIM-1; i <= NN; i++,k = MODNN(k+IPRIM)) {
      q = 1;
      for (j = 1; j <= no_eras; j++)
	if (reg[j] != A0) {
	  reg[j] = MODNN(reg[j] + j);
	  q ^= ALPHA_TO[reg[j]];
	}
      if (q != 0)
	continue;
      /* store root and error location number indices */
      root[count] = i;
      loc[count] = k;
      count++;
    }
    if (count != no_eras) {
      printf("count = %d no_eras = %d\n lambda(x) is WRONG\n",count,no_eras);
      count = -1;
      goto finish;
    }
#if DEBUG >= 2
    printf("\n Erasure positions as determined by roots of Eras Loc Poly:\n");
    for (i = 0; i < count; i++)
      printf("%d ", loc[i]);
    printf("\n");
#endif
#endif
  }
  for(i=0;i<NROOTS+1;i++)
    //    printf("%d  %d  %d\n",i,lambda[i],INDEX_OF[lambda[i]]);
    b[i] = INDEX_OF[lambda[i]];

  /*
   * Begin Berlekamp-Massey algorithm to determine error+erasure
   * locator polynomial
   */
  r = no_eras;
  el = no_eras;
  while (++r <= NROOTS) {	/* r is the step number */
    /* Compute discrepancy at the r-th step in poly-form */
    discr_r = 0;
    for (i = 0; i < r; i++){
      if ((lambda[i] != 0) && (s[r-i-1] != A0)) {
	discr_r ^= ALPHA_TO[MODNN(INDEX_OF[lambda[i]] + s[r-i-1])];
      }
    }
    discr_r = INDEX_OF[discr_r];	/* Index form */
    if (discr_r == A0) {
      /* 2 lines below: B(x) <-- x*B(x) */
      memmove(&b[1],b,NROOTS*sizeof(b[0]));
      b[0] = A0;
    } else {
      /* 7 lines below: T(x) <-- lambda(x) - discr_r*x*b(x) */
      t[0] = lambda[0];
      for (i = 0 ; i < NROOTS; i++) {
	if(b[i] != A0)
	  t[i+1] = lambda[i+1] ^ ALPHA_TO[MODNN(discr_r + b[i])];
	else
	  t[i+1] = lambda[i+1];
      }
      if (2 * el <= r + no_eras - 1) {
	el = r + no_eras - el;
	/*
	 * 2 lines below: B(x) <-- inv(discr_r) *
	 * lambda(x)
	 */
	for (i = 0; i <= NROOTS; i++)
	  b[i] = (lambda[i] == 0) ? A0 : MODNN(INDEX_OF[lambda[i]] - discr_r + NN);
      } else {
	/* 2 lines below: B(x) <-- x*B(x) */
	memmove(&b[1],b,NROOTS*sizeof(b[0]));
	b[0] = A0;
      }
      memcpy(lambda,t,(NROOTS+1)*sizeof(t[0]));
    }
  }

  /* Convert lambda to index form and compute deg(lambda(x)) */
  deg_lambda = 0;
  for(i=0;i<NROOTS+1;i++){
    lambda[i] = INDEX_OF[lambda[i]];
    if(lambda[i] != A0)
      deg_lambda = i;
  }
  /* Find roots of the error+erasure locator polynomial by Chien search */
  memcpy(&reg[1],&lambda[1],NROOTS*sizeof(reg[0]));
  count = 0;		/* Number of roots of lambda(x) */
  for (i = 1,k=IPRIM-1; i <= NN; i++,k = MODNN(k+IPRIM)) {
    q = 1; /* lambda[0] is always 0 */
    for (j = deg_lambda; j > 0; j--){
      if (reg[j] != A0) {
	reg[j] = MODNN(reg[j] + j);
	q ^= ALPHA_TO[reg[j]];
      }
    }
    if (q != 0)
      continue; /* Not a root */
    /* store root (index-form) and error location number */
#if DEBUG>=2
    printf("count %d root %d loc %d\n",count,i,k);
#endif
    root[count] = i;
    loc[count] = k;
    /* If we've already found max possible roots,
     * abort the search to save time
     */
    if(++count == deg_lambda)
      break;
  }
  if (deg_lambda != count) {
    /*
     * deg(lambda) unequal to number of roots => uncorrectable
     * error detected
     */
    count = -1;
    goto finish;
  }
  /*
   * Compute err+eras evaluator poly omega(x) = s(x)*lambda(x) (modulo
   * x**NROOTS). in index form. Also find deg(omega).
   */
  deg_omega = deg_lambda-1;
  for (i = 0; i <= deg_omega;i++){
    tmp = 0;
    for(j=i;j >= 0; j--){
      if ((s[i - j] != A0) && (lambda[j] != A0))
	tmp ^= ALPHA_TO[MODNN(s[i - j] + lambda[j])];
    }
    omega[i] = INDEX_OF[tmp];
  }

  /*
   * Compute error values in poly-form. num1 = omega(inv(X(l))), num2 =
   * inv(X(l))**(FCR-1) and den = lambda_pr(inv(X(l))) all in poly-form
   */
  for (j = count-1; j >=0; j--) {
    num1 = 0;
    for (i = deg_omega; i >= 0; i--) {
      if (omega[i] != A0)
	num1  ^= ALPHA_TO[MODNN(omega[i] + i * root[j])];
    }
    num2 = ALPHA_TO[MODNN(root[j] * (FCR - 1) + NN)];
    den = 0;

    /* lambda[i+1] for i even is the formal derivative lambda_pr of lambda[i] */
    for (i = min(deg_lambda,NROOTS-1) & ~1; i >= 0; i -=2) {
      if(lambda[i+1] != A0)
	den ^= ALPHA_TO[MODNN(lambda[i+1] + i * root[j])];
    }
#if DEBUG >= 1
    if (den == 0) {
      printf("\n ERROR: denominator = 0\n");
      count = -1;
      goto finish;
    }
#endif
    /* Apply error to data */
    if (num1 != 0 && loc[j] >= PAD) {
      data[loc[j]-PAD] ^= ALPHA_TO[MODNN(INDEX_OF[num1] + INDEX_OF[num2] + NN - INDEX_OF[den])];
    }
  }
 finish:
  if(eras_pos != NULL){
    for(i=0;i<count;i++)
      eras_pos[i] = loc[i];
  }
  return count;
}


code:
subroutine multimode_decoder(ss,id2,params,nfsample)

  !$ use omp_lib
  use prog_args
  use timer_module, only: timer
  use jt4_decode
  use jt65_decode
  use jt9_decode
  use ft8_decode

  include 'jt9com.f90'
  include 'timer_common.inc'

  type, extends(jt4_decoder) :: counting_jt4_decoder
     integer :: decoded
  end type counting_jt4_decoder

  type, extends(jt65_decoder) :: counting_jt65_decoder
     integer :: decoded
  end type counting_jt65_decoder

  type, extends(jt9_decoder) :: counting_jt9_decoder
     integer :: decoded
  end type counting_jt9_decoder

  type, extends(ft8_decoder) :: counting_ft8_decoder
     integer :: decoded
  end type counting_ft8_decoder

  real ss(184,NSMAX)
  logical baddata,newdat65,newdat9,single_decode,bVHF,bad0,newdat,ex
  integer*2 id2(NTMAX*12000)
  type(params_block) :: params
  real*4 dd(NTMAX*12000)
  character(len=20) :: datetime
  character(len=12) :: mycall, hiscall
  character(len=6) :: mygrid, hisgrid
  save
  type(counting_jt4_decoder) :: my_jt4
  type(counting_jt65_decoder) :: my_jt65
  type(counting_jt9_decoder) :: my_jt9
  type(counting_ft8_decoder) :: my_ft8

  !cast C character arrays to Fortran character strings
  datetime=transfer(params%datetime, datetime)
  mycall=transfer(params%mycall,mycall)
  hiscall=transfer(params%hiscall,hiscall)
  mygrid=transfer(params%mygrid,mygrid)
  hisgrid=transfer(params%hisgrid,hisgrid)

  ! initialize decode counts
  my_jt4%decoded = 0
  my_jt65%decoded = 0
  my_jt9%decoded = 0
  my_ft8%decoded = 0

  single_decode=iand(params%nexp_decode,32).ne.0
  bVHF=iand(params%nexp_decode,64).ne.0
  if(mod(params%nranera,2).eq.0) ntrials=10**(params%nranera/2)
  if(mod(params%nranera,2).eq.1) ntrials=3*10**(params%nranera/2)
  if(params%nranera.eq.0) ntrials=0

  nfail=0
10 if (params%nagain) then
     open(13,file=trim(temp_dir)//'/decoded.txt',status='unknown',            &
          position='append',iostat=ios)
  else
     open(13,file=trim(temp_dir)//'/decoded.txt',status='unknown',iostat=ios)
  endif
  if(params%nmode.eq.8) then
     inquire(file=trim(temp_dir)//'/houndcallers.txt',exist=ex)
     if(.not.ex) then
        c2fox='            '
        g2fox='    '
        nsnrfox=-99
        nfreqfox=-99
        n30z=0
        nwrap=0
        nfox=0
     endif
     open(19,file=trim(temp_dir)//'/houndcallers.txt',status='unknown')
  endif

  if(ios.ne.0) then
     nfail=nfail+1
     if(nfail.le.3) then
        call sleep_msec(10)
        go to 10
     endif
  endif

  if(params%nmode.eq.8) then
! We're in FT8 mode
     call timer('decft8  ',0)
     newdat=params%newdat
     call my_ft8%decode(ft8_decoded,id2,params%nQSOProgress,params%nfqso,    &
          params%nftx,newdat,params%nutc,params%nfa,params%nfb,              &
          params%nexp_decode,params%ndepth,logical(params%nagain),           &
          logical(params%lft8apon),logical(params%lapcqonly),params%napwid,  &
          mycall,mygrid,hiscall,hisgrid)
     call timer('decft8  ',1)
     if(nfox.gt.0) then
        n30min=minval(n30fox(1:nfox))
        n30max=maxval(n30fox(1:nfox))
     endif
     j=0
     rewind 19
     if(nfox.eq.0) then
        endfile 19
        rewind 19
     else
        do i=1,nfox
           n=n30fox(i)
           if(n30max-n30fox(i).le.4) then
              j=j+1
              c2fox(j)=c2fox(i)
              g2fox(j)=g2fox(i)
              nsnrfox(j)=nsnrfox(i)
              nfreqfox(j)=nfreqfox(i)
              n30fox(j)=n
              m=n30max-n
              if(len(trim(g2fox(j))).eq.4) then
                 call azdist(mygrid,g2fox(j),0.d0,nAz,nEl,nDmiles,nDkm, &
                      nHotAz,nHotABetter)
              else
                 nDkm=9999
              endif
              write(19,1004) c2fox(j),g2fox(j),nsnrfox(j),nfreqfox(j),nDkm,m
1004          format(a12,1x,a4,i5,i6,i7,i3)
           endif
        enddo
        nfox=j
        flush(19)
     endif
     go to 800
  endif

  rms=sqrt(dot_product(float(id2(300000:310000)),            &
       float(id2(300000:310000)))/10000.0)
  if(rms.lt.2.0) go to 800

! Zap data at start that might come from T/R switching transient?
  nadd=100
  k=0
  bad0=.false.
  do i=1,240
     sq=0.
     do n=1,nadd
        k=k+1
        sq=sq + float(id2(k))**2
     enddo
     rms=sqrt(sq/nadd)
     if(rms.gt.10000.0) then
        bad0=.true.
        kbad=k
        rmsbad=rms
     endif
  enddo
  if(bad0) then
     nz=min(NTMAX*12000,kbad+100)
!     id2(1:nz)=0                ! temporarily disabled as it can breaak the JT9 decoder, maybe others
  endif

  if(params%nmode.eq.4 .or. params%nmode.eq.65) open(14,file=trim(temp_dir)// &
       '/avemsg.txt',status='unknown')
  if(params%nmode.eq.164) open(17,file=trim(temp_dir)//'/red.dat',          &
       status='unknown')

  if(params%nmode.eq.4) then
     jz=52*nfsample
     if(params%newdat) then
        if(nfsample.eq.12000) call wav11(id2,jz,dd)
        if(nfsample.eq.11025) dd(1:jz)=id2(1:jz)
     endif
     call my_jt4%decode(jt4_decoded,dd,jz,params%nutc,params%nfqso,         &
          params%ntol,params%emedelay,params%dttol,logical(params%nagain),  &
          params%ndepth,logical(params%nclearave),params%minsync,           &
          params%minw,params%nsubmode,mycall,hiscall,         &
          hisgrid,params%nlist,params%listutc,jt4_average)
     go to 800
  endif

  npts65=52*12000
  if(params%nmode.eq.164) npts65=54*12000
  if(baddata(id2,npts65)) then
     nsynced=0
     ndecoded=0
     go to 800
  endif

  ntol65=params%ntol              !### is this OK? ###
  newdat65=params%newdat
  newdat9=params%newdat

!$call omp_set_dynamic(.true.)
!$omp parallel sections num_threads(2) copyin(/timer_private/) shared(ndecoded) if(.true.) !iif() needed on Mac

!$omp section
  if(params%nmode.eq.65 .or. params%nmode.eq.164 .or.                      &
       (params%nmode.eq.(65+9) .and. params%ntxmode.eq.65)) then
! We're in JT65 or QRA64 mode, or should do JT65 first

     if(newdat65) dd(1:npts65)=id2(1:npts65)
     nf1=params%nfa
     nf2=params%nfb
     call timer('jt65a   ',0)
     call my_jt65%decode(jt65_decoded,dd,npts65,newdat65,params%nutc,      &
          nf1,nf2,params%nfqso,ntol65,params%nsubmode,params%minsync,      &
          logical(params%nagain),params%n2pass,logical(params%nrobust),    &
          ntrials,params%naggressive,params%ndepth,params%emedelay,        &
          logical(params%nclearave),mycall,hiscall,          &
          hisgrid,params%nexp_decode,params%nQSOProgress,           &
          logical(params%ljt65apon))
     call timer('jt65a   ',1)

  else if(params%nmode.eq.9 .or. (params%nmode.eq.(65+9) .and. params%ntxmode.eq.9)) then
! We're in JT9 mode, or should do JT9 first
     call timer('decjt9  ',0)
     call my_jt9%decode(jt9_decoded,ss,id2,params%nfqso,       &
          newdat9,params%npts8,params%nfa,params%nfsplit,params%nfb,       &
          params%ntol,params%nzhsym,logical(params%nagain),params%ndepth,  &
          params%nmode,params%nsubmode,params%nexp_decode)
     call timer('decjt9  ',1)
  endif

!$omp section
  if(params%nmode.eq.(65+9)) then       !Do the other mode (we're in dual mode)
     if (params%ntxmode.eq.9) then
        if(newdat65) dd(1:npts65)=id2(1:npts65)
        nf1=params%nfa
        nf2=params%nfb
        call timer('jt65a   ',0)
        call my_jt65%decode(jt65_decoded,dd,npts65,newdat65,params%nutc,   &
             nf1,nf2,params%nfqso,ntol65,params%nsubmode,params%minsync,   &
             logical(params%nagain),params%n2pass,logical(params%nrobust), &
             ntrials,params%naggressive,params%ndepth,params%emedelay,     &
             logical(params%nclearave),mycall,hiscall,       &
             hisgrid,params%nexp_decode,params%nQSOProgress,        &
             logical(params%ljt65apon))
        call timer('jt65a   ',1)
     else
        call timer('decjt9  ',0)
        call my_jt9%decode(jt9_decoded,ss,id2,params%nfqso,                &
             newdat9,params%npts8,params%nfa,params%nfsplit,params%nfb,    &
             params%ntol,params%nzhsym,logical(params%nagain),             &
             params%ndepth,params%nmode,params%nsubmode,params%nexp_decode)
        call timer('decjt9  ',1)
     end if
  endif

!$omp end parallel sections

! JT65 is not yet producing info for nsynced, ndecoded.
800 ndecoded = my_jt4%decoded + my_jt65%decoded + my_jt9%decoded + my_ft8%decoded
  write(*,1010) nsynced,ndecoded
1010 format('<DecodeFinished>',2i4)
  call flush(6)
  close(13)
  close(19)
  if(params%nmode.eq.4 .or. params%nmode.eq.65) close(14)

  return

contains

  subroutine jt4_decoded(this,snr,dt,freq,have_sync,sync,is_deep,    &
       decoded0,qual,ich,is_average,ave)
    implicit none
    class(jt4_decoder), intent(inout) :: this
    integer, intent(in) :: snr
    real, intent(in) :: dt
    integer, intent(in) :: freq
    logical, intent(in) :: have_sync
    logical, intent(in) :: is_deep
    character(len=1), intent(in) :: sync
    character(len=22), intent(in) :: decoded0
    real, intent(in) :: qual
    integer, intent(in) :: ich
    logical, intent(in) :: is_average
    integer, intent(in) :: ave

    character*22 decoded
    character*3 cflags

    if(ich.eq.-99) stop                         !Silence compiler warning
    if (have_sync) then
       decoded=decoded0
       cflags='   '
       if(decoded.ne.'                      ') cflags='f  '
       if(is_deep) then
          cflags(1:2)='d1'
          write(cflags(3:3),'(i1)') min(int(qual),9)
          if(qual.ge.10.0) cflags(3:3)='*'
          if(qual.lt.3.0) decoded(22:22)='?'
       endif
       if(is_average) then
          write(cflags(2:2),'(i1)') min(ave,9)
          if(ave.ge.10) cflags(2:2)='*'
       endif
       write(*,1000) params%nutc,snr,dt,freq,sync,decoded,cflags
1000   format(i4.4,i4,f5.1,i5,1x,'$',a1,1x,a22,1x,a3)
    else
       write(*,1000) params%nutc,snr,dt,freq
    end if

    select type(this)
    type is (counting_jt4_decoder)
       this%decoded = this%decoded + 1
    end select
  end subroutine jt4_decoded

  subroutine jt4_average (this, used, utc, sync, dt, freq, flip)
    implicit none
    class(jt4_decoder), intent(inout) :: this
    logical, intent(in) :: used
    integer, intent(in) :: utc
    real, intent(in) :: sync
    real, intent(in) :: dt
    integer, intent(in) :: freq
    logical, intent(in) :: flip
    character(len=1) :: cused, csync

    cused = '.'
    csync = '*'
    if (used) cused = '$'
    if (flip) csync = '$'
    write(14,1000) cused,utc,sync,dt,freq,csync
1000 format(a1,i5.4,f6.1,f6.2,i6,1x,a1)
  end subroutine jt4_average

  subroutine jt65_decoded(this,sync,snr,dt,freq,drift,nflip,width,     &
       decoded0,ft,qual,nsmo,nsum,minsync)

    use jt65_decode
    implicit none

    class(jt65_decoder), intent(inout) :: this
    real, intent(in) :: sync
    integer, intent(in) :: snr
    real, intent(in) :: dt
    integer, intent(in) :: freq
    integer, intent(in) :: drift
    integer, intent(in) :: nflip
    real, intent(in) :: width
    character(len=22), intent(in) :: decoded0
    integer, intent(in) :: ft
    integer, intent(in) :: qual
    integer, intent(in) :: nsmo
    integer, intent(in) :: nsum
    integer, intent(in) :: minsync

    integer i,nap,nft
    logical is_deep,is_average
    character decoded*22,csync*2,cflags*3

    if(width.eq.-9999.0) stop              !Silence compiler warning
!$omp critical(decode_results)
    decoded=decoded0
    cflags='   '
    is_deep=ft.eq.2

    if(ft.ge.80) then                      !QRA64 mode
       nft=ft-100
       csync=': '
       if(sync-3.4.ge.float(minsync) .or. nft.ge.0) csync=':*'
       if(nft.lt.0) then
          write(*,1009) params%nutc,snr,dt,freq,csync,decoded
       else
          write(*,1009) params%nutc,snr,dt,freq,csync,decoded,nft
1009      format(i4.4,i4,f5.1,i5,1x,a2,1x,a22,i2)
       endif
       write(13,1011) params%nutc,nint(sync),snr,dt,float(freq),drift,    &
            decoded,nft
1011   format(i4.4,i4,i5,f6.2,f8.0,i4,3x,a22,' QRA64',i3)
       go to 100
    endif

    if(ft.eq.0 .and. minsync.ge.0 .and. int(sync).lt.minsync) then
       write(*,1010) params%nutc,snr,dt,freq
    else
       is_average=nsum.ge.2
       if(bVHF .and. ft.gt.0) then
          cflags='f  '
          if(is_deep) then
             cflags(1:2)='d1'
             write(cflags(3:3),'(i1)') min(qual,9)
             if(qual.ge.10) cflags(3:3)='*'
             if(qual.lt.3) decoded(22:22)='?'
          endif
          if(is_average) then
             write(cflags(2:2),'(i1)') min(nsum,9)
             if(nsum.ge.10) cflags(2:2)='*'
          endif
          nap=ishft(ft,-2)
          if(nap.ne.0) then
            write(cflags(1:3),'(a1,i1)') 'a',nap
          endif
       endif
       csync='# '
       i=0
       if(bVHF .and. nflip.ne.0 .and.                         &
            sync.ge.max(0.0,float(minsync))) then
          csync='#*'
          if(nflip.eq.-1) then
             csync='##'
             if(decoded.ne.'                      ') then
                do i=22,1,-1
                   if(decoded(i:i).ne.' ') exit
                enddo
                if(i.gt.18) i=18
                decoded(i+2:i+4)='OOO'
             endif
          endif
       endif
       write(*,1010) params%nutc,snr,dt,freq,csync,decoded,cflags
1010   format(i4.4,i4,f5.1,i5,1x,a2,1x,a22,1x,a3)
    endif
    write(13,1012) params%nutc,nint(sync),snr,dt,float(freq),drift,    &
         decoded,ft,nsum,nsmo
1012 format(i4.4,i4,i5,f6.2,f8.0,i4,3x,a22,' JT65',3i3)

100 call flush(6)

!$omp end critical(decode_results)
    select type(this)
    type is (counting_jt65_decoder)
       this%decoded = this%decoded + 1
    end select
  end subroutine jt65_decoded

  subroutine jt9_decoded (this, sync, snr, dt, freq, drift, decoded)
    use jt9_decode
    implicit none

    class(jt9_decoder), intent(inout) :: this
    real, intent(in) :: sync
    integer, intent(in) :: snr
    real, intent(in) :: dt
    real, intent(in) :: freq
    integer, intent(in) :: drift
    character(len=22), intent(in) :: decoded

    !$omp critical(decode_results)
    write(*,1000) params%nutc,snr,dt,nint(freq),decoded
1000 format(i4.4,i4,f5.1,i5,1x,'@ ',1x,a22)
    write(13,1002) params%nutc,nint(sync),snr,dt,freq,drift,decoded
1002 format(i4.4,i4,i5,f6.1,f8.0,i4,3x,a22,' JT9')
    call flush(6)
    !$omp end critical(decode_results)
    select type(this)
    type is (counting_jt9_decoder)
       this%decoded = this%decoded + 1
    end select
  end subroutine jt9_decoded

  subroutine ft8_decoded (this,sync,snr,dt,freq,decoded,nap,qual)
    use ft8_decode
    implicit none

    class(ft8_decoder), intent(inout) :: this
    real, intent(in) :: sync
    integer, intent(in) :: snr
    real, intent(in) :: dt
    real, intent(in) :: freq
    character(len=37), intent(in) :: decoded
    character c1*12,c2*12,g2*4,w*4,ctmp*12
    integer i0,i1,i2,i3,i4,i5,i6,n30,nwrap
    integer, intent(in) :: nap
    real, intent(in) :: qual
    character*2 annot
    character*37 decoded0
    logical isgrid4,first,b0,b1,b2
    data first/.true./
    save

    isgrid4(w)=(len_trim(w).eq.4 .and.                                        &
         ichar(w(1:1)).ge.ichar('A') .and. ichar(w(1:1)).le.ichar('R') .and.  &
         ichar(w(2:2)).ge.ichar('A') .and. ichar(w(2:2)).le.ichar('R') .and.  &
         ichar(w(3:3)).ge.ichar('0') .and. ichar(w(3:3)).le.ichar('9') .and.  &
         ichar(w(4:4)).ge.ichar('0') .and. ichar(w(4:4)).le.ichar('9'))

    if(first) then
       c2fox='            '
       g2fox='    '
       nsnrfox=-99
       nfreqfox=-99
       n30z=0
       nwrap=0
       nfox=0
       first=.false.
    endif

    decoded0=decoded

    annot='  '
    if(nap.ne.0) then
       write(annot,'(a1,i1)') 'a',nap
       if(qual.lt.0.17) decoded0(22:22)='?'
    endif

    i0=index(decoded0,';')
    if(i0.le.0) write(*,1000) params%nutc,snr,dt,nint(freq),decoded0(1:22),annot
1000 format(i6.6,i4,f5.1,i5,' ~ ',1x,a22,1x,a2)
    if(i0.gt.0) write(*,1001) params%nutc,snr,dt,nint(freq),decoded0
1001 format(i6.6,i4,f5.1,i5,' ~ ',1x,a37)
    write(13,1002) params%nutc,nint(sync),snr,dt,freq,0,decoded0
1002 format(i6.6,i4,i5,f6.1,f8.0,i4,3x,a37,' FT8')

    i1=index(decoded0,' ')
    i2=i1 + index(decoded0(i1+1:),' ')
    i3=i2 + index(decoded0(i2+1:),' ')
    if(i1.ge.3 .and. i2.ge.7 .and. i3.ge.10) then
       c1=decoded0(1:i1-1)//'            '
       c2=decoded0(i1+1:i2-1)
       g2=decoded0(i2+1:i3-1)
       b0=c1.eq.mycall
       if(c1(1:3).eq.'DE ' .and. index(c2,'/').ge.2) b0=.true.
       if(len(trim(c1)).ne.len(trim(mycall))) then
          i4=index(trim(c1),trim(mycall))
          i5=index(trim(mycall),trim(c1))
          if(i4.ge.1 .or. i5.ge.1) b0=.true.
       endif
       b1=i3-i2.eq.5 .and. isgrid4(g2)
       b2=i3-i2.eq.1
       if(b0 .and. (b1.or.b2) .and. nint(freq).ge.1000) then
          n=params%nutc
          n30=(3600*(n/10000) + 60*mod((n/100),100) + mod(n,100))/30
          if(n30.lt.n30z) nwrap=nwrap+5760    !New UTC day, handle the wrap
          n30z=n30
          n30=n30+nwrap
          nfox=nfox+1
          c2fox(nfox)=c2
          g2fox(nfox)=g2
          nsnrfox(nfox)=snr
          nfreqfox(nfox)=nint(freq)
          n30fox(nfox)=n30
       endif
    endif

    call flush(6)
    call flush(13)

    select type(this)
    type is (counting_ft8_decoder)
       this%decoded = this%decoded + 1
    end select

    return
  end subroutine ft8_decoded

end subroutine multimode_decoder


code:
subroutine deep4(sym0,neme,flip,mycall,hiscall,hisgrid,decoded,qual)

! Deep search routine for JT4

  use prog_args
  parameter (MAXCALLS=7000,MAXRPT=63)
  real*4 sym0(206),sym(206)
  character callsign*12,grid*4,message*22,hisgrid*6,ceme*3
  character*12 mycall,hiscall
  character mycall0*12,hiscall0*12,hisgrid0*6
  character*22 decoded
  character*22 testmsg(2*MAXCALLS + 2 + MAXRPT)
  character*15 callgrid(MAXCALLS)
  character*180 line
  character*4 rpt(MAXRPT)
  integer ncode(206)
  real*4   code(206,2*MAXCALLS + 2 + MAXRPT)
  real pp(2*MAXCALLS + 2 + MAXRPT)
  data neme0/-99/
  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 mycall0,hiscall0,hisgrid0,neme0,ntot,code,testmsg

  sym=sym0
  if(mycall.eq.mycall0 .and. hiscall.eq.hiscall0 .and.         &
       hisgrid.eq.hisgrid0 .and. neme.eq.neme0) go to 30

  open(23,file=trim(data_dir)//'/CALL3.TXT',status='unknown')
  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. 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 encode4(message,ncode)
        code(1:206,k)=2*ncode(1:206)-1
        if(n.ge.2) then
! Insert CQ message
           if(j4.lt.1) callgrid(icall)=callsign(1:j2)//' '//grid
           message='CQ '//callgrid(icall)
           k=k+1
           testmsg(k)=message
           call encode4(message,ncode)
        code(1:206,k)=2*ncode(1:206)-1
        endif
     enddo
10   continue
  enddo

20 continue
  close(23)
  ntot=k

30 mycall0=mycall
  hiscall0=hiscall
  hisgrid0=hisgrid
  neme0=neme

  sq=0.
  do j=1,206
     sq=sq + sym(j)**2
  enddo
  rms=sqrt(sq/206.0)
  sym=sym/rms

  p1=-1.e30
  p2=-1.e30
  do k=1,ntot
     pp(k)=0.
! Should re-instate the following:
!     if(k.ge.2 .and. k.le.64 .and. flip.gt.0.0) cycle
! 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
        p=0.
        do j=1,206
           p=p + code(j,k)*sym(j)
        enddo
        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. testmsg(i).ne.testmsg(ip1)) p2=pp(i)
  enddo

  qual=p1-max(1.15*p2,70.0)

! ### DO NOT REMOVE ###
  rewind 77
  if(ip1.ge.1 .and. ip1.le.2*MAXCALLS+2+MAXRPT) write(77,1001) p1,p2,ntot,   &
       rms,qual,ip1,testmsg(ip1)
1001 format(2f8.2,i8,2f8.2,i6,2x,a22)
  call flush(77)

  if(qual.gt.1.0) then
     decoded=testmsg(ip1)
  else
     decoded='                      '
     qual=0.
  endif

! 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

  return
end subroutine deep4


code:
subroutine deg2grid(dlong0,dlat,grid)

  real dlong                        !West longitude (deg)
  real dlat                         !Latitude (deg)
  character grid*6

  dlong=dlong0
  if(dlong.lt.-180.0) dlong=dlong+360.0
  if(dlong.gt.180.0) dlong=dlong-360.0

! Convert to units of 5 min of longitude, working east from 180 deg.
  nlong=int(60.0*(180.0-dlong)/5.0)
  n1=nlong/240                      !20-degree field
  n2=(nlong-240*n1)/24              !2 degree square
  n3=nlong-240*n1-24*n2             !5 minute subsquare
  grid(1:1)=char(ichar('A')+n1)
  grid(3:3)=char(ichar('0')+n2)
  grid(5:5)=char(ichar('a')+n3)

! Convert to units of 2.5 min of latitude, working north from -90 deg.
  nlat=int(60.0*(dlat+90)/2.5)
  n1=nlat/240                       !10-degree field
  n2=(nlat-240*n1)/24               !1 degree square
  n3=nlat-240*n1-24*n2              !2.5 minuts subsquare
  grid(2:2)=char(ichar('A')+n1)
  grid(4:4)=char(ichar('0')+n2)
  grid(6:6)=char(ichar('a')+n3)

  return
end subroutine deg2grid


code:
subroutine degrade_snr(d2,npts,db,bw)

! Degrade S/N by specified number of dB.

  integer*2 d2(npts)

  p0=0.
  do i=1,npts
     x=d2(i)
     p0=p0 + x*x
  enddo
  p0=p0/npts
  if(bw.gt.0.0) p0=p0*6000.0/bw
  s=sqrt(p0*(10.0**(0.1*db) - 1.0))
  fac=sqrt(p0/(p0+s*s))
  do i=1,npts
     d2(i)=nint(fac*(d2(i) + s*gran()))
  enddo

  return
end subroutine degrade_snr


code:
subroutine demod64a(s3,nadd,afac1,mrsym,mrprob,mr2sym,mr2prob,ntest,nlow)

! Demodulate the 64-bin spectra for each of 63 symbols in a frame.

! Parameters
!    nadd     number of spectra already summed
!    mrsym    most reliable symbol value
!    mr2sym   second most likely symbol value
!    mrprob   probability that mrsym was the transmitted value
!    mr2prob  probability that mr2sym was the transmitted value

  implicit real*8 (a-h,o-z)
  real*4 s3(64,63),afac1
  integer mrsym(63),mrprob(63),mr2sym(63),mr2prob(63)

  if(nadd.eq.-999) return
  afac=afac1 * float(nadd)**0.64
  scale=255.999

! Compute average spectral value
  ave=sum(s3)/(64.*63.)
  i1=1                                      !Silence warning
  i2=1

! Compute probabilities for most reliable symbol values
  do j=1,63
     s1=-1.e30
     psum=0.
     do i=1,64
        x=min(afac*s3(i,j)/ave,50.d0)
        psum=psum+s3(i,j)
        if(s3(i,j).gt.s1) then
           s1=s3(i,j)
           i1=i                              !Most reliable
        endif
     enddo
     if(psum.eq.0.0) psum=1.e-6

     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=s1/psum                              !Symbol metrics for ftrsd
     p2=s2/psum
     mrsym(j)=i1-1
     mr2sym(j)=i2-1
     mrprob(j)=scale*p1
     mr2prob(j)=scale*p2
  enddo

  nlow=0
  do j=1,63
     if(mrprob(j).le.5) nlow=nlow+1
  enddo
  ntest=sum(mrprob)

  return
end subroutine demod64a


code:
real*8 function determ(array,norder)
  implicit real*8 (a-h,o-z)
  real*8 array(10,10)

  determ=1.
  do k=1,norder
     if (array(k,k).ne.0) go to 41
     do j=k,norder
        if(array(k,j).ne.0) go to 31
     enddo
     determ=0.
     go to 60

31   do i=k,norder
        s8=array(i,j)
        array(i,j)=array(i,k)
        array(i,k)=s8
     enddo
     determ=-1.*determ
41   determ=determ*array(k,k)
     if(k.lt.norder) then
        k1=k+1
        do i=k1,norder
           do j=k1,norder
              array(i,j)=array(i,j)-array(i,k)*array(k,j)/array(k,k)
           enddo
        enddo
     end if
  enddo

60 return
end function determ


code:
subroutine downsam9(id2,npts8,nsps8,newdat,nspsd,fpk,c2)

!Downsample from id2() into c2() so as to yield nspsd samples per symbol,
!mixing from fpk down to zero frequency.  The downsample factor is 432.

  use, intrinsic :: iso_c_binding
  use FFTW3
  use timer_module, only: timer

  include 'constants.f90'
  integer(C_SIZE_T) NMAX1
  parameter (NMAX1=653184)
  parameter (NFFT1=653184,NFFT2=1512)
  type(C_PTR) :: plan                        !Pointers plan for big FFT
  integer*2 id2(0:8*npts8-1)
  logical, intent(inout) :: newdat
  real*4, pointer :: x1(:)
  complex c1(0:NFFT1/2)
  complex c2(0:NFFT2-1)
  real s(5000)
  logical first
  common/patience/npatience,nthreads
  data first/.true./
  save plan,first,c1,s,x1

  df1=12000.0/NFFT1
  npts=8*npts8
  if(npts.gt.NFFT1) npts=NFFT1  !### Fix! ###

  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
! Plan the FFTs just once

     !$omp critical(fftw) ! serialize non thread-safe FFTW3 calls
     plan=fftwf_alloc_real(NMAX1)
     call c_f_pointer(plan,x1,[NMAX1])
     x1(0:NMAX1-1) => x1        !remap bounds
     call fftwf_plan_with_nthreads(nthreads)
     plan=fftwf_plan_dft_r2c_1d(NFFT1,x1,c1,nflags)
     call fftwf_plan_with_nthreads(1)
     !$omp end critical(fftw)

     first=.false.
  endif

  if(newdat) then
     x1(0:npts-1)=id2(0:npts-1)
     x1(npts:NFFT1-1)=0.                      !Zero the rest of x1
     call timer('FFTbig9 ',0)
     call fftwf_execute_dft_r2c(plan,x1,c1)
     call timer('FFTbig9 ',1)

     nadd=int(1.0/df1)
     s=0.
     do i=1,5000
        j=int((i-1)/df1)
        do n=1,nadd
           j=j+1
           s(i)=s(i)+real(c1(j))**2 + aimag(c1(j))**2
        enddo
     enddo
     newdat=.false.
  endif

  ndown=8*nsps8/nspsd                      !Downsample factor = 432
  nh2=NFFT2/2
  nf=nint(fpk)
  i0=int(fpk/df1)

  nw=100
  ia=max(1,nf-nw)
  ib=min(5000,nf+nw)
  call pctile(s(ia),ib-ia+1,40,avenoise)

  fac=sqrt(1.0/avenoise)
  do i=0,NFFT2-1
     j=i0+i
     if(i.gt.nh2) j=j-NFFT2
     c2(i)=fac*c1(j)
  enddo
  call four2a(c2,NFFT2,1,1,1)              !FFT back to time domain

  return
end subroutine downsam9


code:
subroutine dxped_fifo(cx,gx,isnrx)

  parameter (NCALLS=268)
  character*6 xcall(NCALLS)
  character*4 xgrid(NCALLS)
  integer isnr(NCALLS)

  character cx*6,gx*4
  common/dxpfifo/nc,isnr,xcall,xgrid

  if(nc.lt.NCALLS) then
     nc=nc+1
     cx=xcall(nc)
     gx=xgrid(nc)
     isnrx=isnr(nc)
  else
     cx='      '
     gx='    '
     isnrx=0
  endif

  return
end subroutine dxped_fifo


code:
subroutine encode232(dat,nsym,symbol)

! Convolutional encoder for a K=32, r=1/2 code.

  integer*1 dat(13)                 !User data, packed 8 bits per byte
  integer*1 symbol(206)             !Channel symbols, one bit per byte
  integer*1 i1
  include 'conv232.f90'

  nstate=0
  k=0
  do j=1,nsym
     do i=7,0,-1
        i1=dat(j)
        i4=i1
        if (i4.lt.0) i4=i4+256
        nstate=ior(ishft(nstate,1),iand(ishft(i4,-i),1))
        n=iand(nstate,npoly1)
        n=ieor(n,ishft(n,-16))
        k=k+1
        symbol(k)=partab(iand(ieor(n,ishft(n,-8)),255))
        n=iand(nstate,npoly2)
        n=ieor(n,ishft(n,-16))
        k=k+1
        symbol(k)=partab(iand(ieor(n,ishft(n,-8)),255))
        if(k.ge.nsym) go to 100
     enddo
  enddo

100 continue

  return
end subroutine encode232


code:
subroutine encode4(message,ncode)

  use packjt
  parameter (MAXCALLS=7000,MAXRPT=63)
  integer ncode(206)
  character*22 message          !Message to be generated
  character*3 cok               !'   ' or 'OOO'
  integer dgen(13)
  integer*1 data0(13),symbol(216)

  call chkmsg(message,cok,nspecial,flip)
  call packmsg(message,dgen,itype,.false.) !Pack 72-bit message into 12 six-bit symbols
  call entail(dgen,data0)
  call encode232(data0,206,symbol)       !Convolutional encoding
  call interleave4(symbol,1)             !Apply JT4 interleaving
  do i=1,206
     ncode(i)=symbol(i)
  enddo

end subroutine encode4


code:
subroutine encode_msk40(message,codeword)
! Encode a 16-bit message and return a 32-bit codeword.
! The code is a (32,16) regular ldpc code with column weight 3.
! The code was generated using the PEG algorithm.
! After creating the codeword, the columns are re-ordered according to
! "colorder" to make the codeword compatible with the parity-check
! matrix stored in Radford Neal's "pchk" format.
!
integer*1 codeword(32)
integer*1 colorder(32)
integer g(16)
integer*1 gen40(16,16)
integer*1 itmp(32)
integer*1 message(16)
integer*1 pchecks(16)
logical first
data first/.true./
data g/Z'4428',Z'5a6b',Z'1b04',Z'2c12',Z'60c4',Z'1071',Z'be6a',Z'36dd', &
       Z'c580',Z'ad9a',Z'eca2',Z'7843',Z'332e',Z'a685',Z'5906',Z'1efe'/
data colorder/4,1,2,3,0,8,6,10,13,28,20,23,17,15,27,25, &
             16,12,18,19,7,21,22,11,24,5,26,14,9,29,30,31/
save first,gen40

if( first ) then ! fill the generator matrix
  gen40=0
  do i=1,16
    do j=1,16
      if( btest(g(i),16-j) ) gen40(i,j)=1
    enddo
  enddo
  first=.false.
endif

do i=1,16
  nsum=0
  do j=1,16
    nsum=nsum+message(j)*gen40(i,j)
  enddo
  pchecks(i)=mod(nsum,2)
enddo
itmp(1:16)=pchecks
itmp(17:32)=message(1:16)
codeword(colorder+1)=itmp(1:32)

return
end subroutine encode_msk40


code:
/* Reed-Solomon encoder
 * Copyright 2002, Phil Karn, KA9Q
 * May be used under the terms of the GNU General Public License (GPL)
 */
#include <string.h>

#ifdef FIXED
#include "fixed.h"
#elif defined(BIGSYM)
#include "int.h"
#else
#include "char.h"
#endif

void ENCODE_RS(
#ifdef FIXED
DTYPE *data, DTYPE *bb,int pad){
#else
void *p,DTYPE *data, DTYPE *bb){
  struct rs *rs = (struct rs *)p;
#endif
  int i, j;
  DTYPE feedback;

#ifdef FIXED
  /* Check pad parameter for validity */
  if(pad < 0 || pad >= NN)
    return;
#endif

  memset(bb,0,NROOTS*sizeof(DTYPE));

  for(i=0;i<NN-NROOTS-PAD;i++){
    feedback = INDEX_OF[data[i] ^ bb[0]];
    if(feedback != A0){      /* feedback term is non-zero */
#ifdef UNNORMALIZED
      /* This line is unnecessary when GENPOLY[NROOTS] is unity, as it must
       * always be for the polynomials constructed by init_rs()
       */
      feedback = MODNN(NN - GENPOLY[NROOTS] + feedback);
#endif
      for(j=1;j<NROOTS;j++)
	bb[j] ^= ALPHA_TO[MODNN(feedback + GENPOLY[NROOTS-j])];
    }
    /* Shift */
    memmove(&bb[0],&bb[1],sizeof(DTYPE)*(NROOTS-1));
    if(feedback != A0)
      bb[NROOTS-1] = ALPHA_TO[MODNN(feedback + GENPOLY[0])];
    else
      bb[NROOTS-1] = 0;
  }
}


code:
subroutine entail(dgen,data0)

! Move 72-bit packed data from 6-bit to 8-bit symbols and add a zero tail.
  integer dgen(13)
  integer*1 data0(13)

  i4=0
  k=0
  m=0
  do i=1,12
     n=dgen(i)
     do j=1,6
        k=k+1
        i4=i4+i4+iand(1,ishft(n,j-6))
        i4=iand(i4,255)
        if(k.eq.8) then
           m=m+1
           if(i4.gt.127) i4=i4-256
           data0(m)=i4
           k=0
        endif
     enddo
  enddo
  do m=10,13
     data0(m)=0
  enddo

  return
end subroutine entail


code:
subroutine ephem(mjd0,dut,east_long,geodetic_lat,height,nspecial,     &
     RA,Dec,Az,El,techo,dop,fspread_1GHz,vr)

  implicit real*8 (a-h,o-z)
  real*8 jd                      !Time of observationa as a Julian Date
  real*8 mjd,mjd0                !Modified Julian Date
  real*8 prec(3,3)               !Precession matrix, J2000 to Date
  real*8 rmatn(3,3)              !Nutation matrix
  real*8 rme2000(6)              !Vector from Earth center to Moon, JD2000
  real*8 rmeDate(6)              !Vector from Earth center to Moon at Date
  real*8 rmeTrue(6)              !Include nutation
  real*8 raeTrue(6)              !Vector from Earth center to Obs at Date
  real*8 rmaTrue(6)              !Vector from Obs to Moon at Date
  logical km,bary                !Set km=.true. to get km, km/s from ephemeris
  common/stcomx/km,bary,pvsun(6) !Common used in JPL subroutines
  common/librcom/xl(2),b(2)

  twopi=8.d0*atan(1.d0)          !Define some constants
  rad=360.d0/twopi
  clight=2.99792458d5
  au2km=0.1495978706910000d9
  pi=0.5d0*twopi
  pio2=0.5d0*pi
  km=.true.
  freq=1000.0d6

  do jj=1,2
     mjd=mjd0
     if(jj.eq.1) mjd=mjd - 1.d0/1440.d0
     djutc=mjd
     jd=2400000.5d0 + mjd
     djtt=mjd + sla_DTT(jd)/86400.d0
     ttjd=jd + sla_DTT(jd)/86400.d0

     if(nspecial.ne.8) then
        call pleph(ttjd,10,3,rme2000)            !RME (J2000) from JPL ephemeris

        year=2000.d0 + (jd-2451545.d0)/365.25d0
        call sla_PREC (2000.0d0, year, prec)     !Get precession matrix
        rmeDate(1:3)=matmul(prec,rme2000(1:3))   !Moon geocentric xyz at Date
        rmeDate(4:6)=matmul(prec,rme2000(4:6))   !Moon geocentric vel at Date
     else
        call sla_DMOON(djtt,rmeDate)             !No JPL ephemeris, use DMOON
        rmeDate=rmeDate*au2km
     endif

     if(nspecial.eq.7) then
        rmeTrue=rmeDate
     else
!Nutation to true equinox of Date
        call sla_NUT(djtt,rmatn)
        call sla_DMXV(rmatn,rmeDate,rmeTrue)
        call sla_DMXV(rmatn,rmeDate(4),rmeTrue(4))
     endif

! Local Apparent Sidereal Time:
     djut1=djutc + dut/86400.d0
     if(nspecial.eq.6) djut1=djutc
     xlast=sla_DRANRM(sla_GMST(djut1) + sla_EQEQX(djtt) + east_long)
     call sla_PVOBS(geodetic_lat,height,xlast,raeTrue)
     rmaTrue=rmeTrue - raeTrue*au2km

     if(nspecial.ne.2) then
! Allow for planetary aberration
        tl=499.004782D0*SQRT(rmaTrue(1)**2 + rmaTrue(2)**2 + rmaTrue(3)**2)
        rmaTrue(1:3)=rmaTrue(1:3)-tl*rmaTrue(4:6)/au2km
     endif

!Topocentric RA, Dec, dist, velocity
     call sla_DC62S(rmaTrue,RA,Dec,dist,RAdot,DECdot,vr)
     dop=-2.d0 * freq * vr/clight                    !EME doppler shift
     techo=2.d0*dist/clight                          !Echo delay time (s)
     call libration(jd,RA,Dec,xl(jj),b(jj))
  enddo

  fspread_1GHz=0.0d0
  dldt=57.2957795131*(xl(2)-xl(1))
  dbdt=57.2957795131*(b(2)-b(1))
  rate=sqrt((2*dldt)**2 + (2*dbdt)**2)
  fspread_1GHz=0.5*6741*rate

  call sla_DE2H(xlast-RA,Dec,geodetic_lat,Az,El)

  return
end subroutine ephem


code:
subroutine extract(s3,nadd,mode65,ntrials,naggressive,ndepth,nflip,     &
     mycall_12,hiscall_12,hisgrid,nQSOProgress,ljt65apon,               &
     nexp_decode,ncount,      &
     nhist,decoded,ltext,nft,qual)

! Input:
!   s3       64-point spectra for each of 63 data symbols
!   nadd     number of spectra summed into s3
!   nqd      0/1 to indicate decode attempt at QSO frequency

! Output:
!   ncount   number of symbols requiring correction (-1 for no KV decode)
!   nhist    maximum number of identical symbol values
!   decoded  decoded message (if ncount >=0)
!   ltext    true if decoded message is free text
!   nft      0=no decode; 1=FT decode; 2=hinted decode

  use prog_args                       !shm_key, exe_dir, data_dir
  use packjt
  use jt65_mod
  use timer_module, only: timer

  real s3(64,63)
  character decoded*22, apmessage*22
  character*12 mycall_12,hiscall_12
  character*6 mycall,hiscall,hisgrid
  character*6 mycall0,hiscall0,hisgrid0
  integer apsymbols(7,12),ap(12)
  integer nappasses(0:5)  ! the number of decoding passes to use for each QSO state
  integer naptypes(0:5,4) ! (nQSOProgress, decoding pass)  maximum of 4 passes for now
  integer dat4(12)
  integer mrsym(63),mr2sym(63),mrprob(63),mr2prob(63)
  integer correct(63),tmp(63)
  logical first,ltext,ljt65apon
  common/chansyms65/correct
  data first/.true./
  save

  if(mode65.eq.-99) stop                   !Silence compiler warning
  if(first) then

! aptype
!------------------------
!   1        CQ     ???    ???
!   2        MyCall ???    ???
!   3        MyCall DxCall ???
!   4        MyCall DxCall RRR
!   5        MyCall DxCall 73
!   6        MyCall DxCall DxGrid
!   7        CQ     DxCall DxGrid

     apsymbols=-1
     nappasses=(/3,4,2,3,3,4/)
     naptypes(0,1:4)=(/1,2,6,0/)  ! Tx6
     naptypes(1,1:4)=(/2,3,6,7/)  ! Tx1
     naptypes(2,1:4)=(/2,3,0,0/)  ! Tx2
     naptypes(3,1:4)=(/3,4,5,0/)  ! Tx3
     naptypes(4,1:4)=(/3,4,5,0/)  ! Tx4
     naptypes(5,1:4)=(/2,3,4,5/)  ! Tx5
     first=.false.
  endif

  mycall=mycall_12(1:6)
  hiscall=hiscall_12(1:6)
! Fill apsymbols array
  if(ljt65apon .and.                                             &
     (mycall.ne.mycall0 .or. hiscall.ne.hiscall0 .or. hisgrid.ne.hisgrid0)) then
!write(*,*) 'initializing apsymbols '
     apsymbols=-1
     mycall0=mycall
     hiscall0=hiscall
     ap=-1
     apsymbols(1,1:4)=(/62,32,32,49/) ! CQ
     if(len_trim(mycall).gt.0) then
        apmessage=mycall//" "//mycall//" RRR"
        call packmsg(apmessage,ap,itype,.false.)
        if(itype.ne.1) ap=-1
        apsymbols(2,1:4)=ap(1:4)
!write(*,*) 'mycall symbols ',ap(1:4)
        if(len_trim(hiscall).gt.0) then
           apmessage=mycall//" "//hiscall//" RRR"
           call packmsg(apmessage,ap,itype,.false.)
           if(itype.ne.1) ap=-1
           apsymbols(3,1:9)=ap(1:9)
           apsymbols(4,:)=ap
           apmessage=mycall//" "//hiscall//" 73"
           call packmsg(apmessage,ap,itype,.false.)
           if(itype.ne.1) ap=-1
           apsymbols(5,:)=ap
           if(len_trim(hisgrid(1:4)).gt.0) then
              apmessage=mycall//' '//hiscall//' '//hisgrid(1:4)
              call packmsg(apmessage,ap,itype,.false.)
              if(itype.ne.1) ap=-1
              apsymbols(6,:)=ap
              apmessage='CQ'//' '//hiscall//' '//hisgrid(1:4)
              call packmsg(apmessage,ap,itype,.false.)
              if(itype.ne.1) ap=-1
              apsymbols(7,:)=ap
           endif
        endif
     endif
  endif

  qual=0.
  nbirdie=20
  npct=50
  afac1=1.1
  nft=0
  nfail=0
  decoded='                      '
  call pctile(s3,4032,npct,base)
  s3=s3/base
  s3a=s3                                            !###

! Get most reliable and second-most-reliable symbol values, and their
! probabilities
1 call demod64a(s3,nadd,afac1,mrsym,mrprob,mr2sym,mr2prob,ntest,nlow)

  call chkhist(mrsym,nhist,ipk)       !Test for birdies and QRM
  if(nhist.ge.nbirdie) then
     nfail=nfail+1
     call pctile(s3,4032,npct,base)
     s3(ipk,1:63)=base
     if(nfail.gt.30) then
        decoded='                      '
        ncount=-1
        go to 900
     endif
     go to 1
  endif

  mrs=mrsym
  mrs2=mr2sym

  call graycode65(mrsym,63,-1)        !Remove gray code
  call interleave63(mrsym,-1)         !Remove interleaving
  call interleave63(mrprob,-1)

  call graycode65(mr2sym,63,-1)      !Remove gray code and interleaving
  call interleave63(mr2sym,-1)       !from second-most-reliable symbols
  call interleave63(mr2prob,-1)

  npass=1  ! if ap decoding is disabled
  if(ljt65apon .and. len_trim(mycall).gt.0) then
     npass=1+nappasses(nQSOProgress)
!write(*,*) 'ap is on: ',npass-1,'ap passes of types ',naptypes(nQSOProgress,:)
  endif
  do ipass=1,npass
     ap=-1
     ntype=0
     if(ipass.gt.1) then
       ntype=naptypes(nQSOProgress,ipass-1)
!write(*,*) 'ap pass, type ',ntype
       ap=apsymbols(ntype,:)
       if(count(ap.ge.0).eq.0) cycle  ! don't bother if all ap symbols are -1
!write(*,'(12i3)') ap
     endif
     ntry=0
     call timer('ftrsd   ',0)
     param=0
     call ftrsdap(mrsym,mrprob,mr2sym,mr2prob,ap,ntrials,correct,param,ntry)
     call timer('ftrsd   ',1)
     ncandidates=param(0)
     nhard=param(1)
     nsoft=param(2)
     nerased=param(3)
     rtt=0.001*param(4)
     ntotal=param(5)
     qual=0.001*param(7)
     nd0=81
     r0=0.87
     if(naggressive.eq.10) then
        nd0=83
        r0=0.90
     endif

     if(ntotal.le.nd0 .and. rtt.le.r0) then
        nft=1+ishft(ntype,2)
     endif

     if(nft.gt.0) exit
  enddo
!write(*,*) nft
  if(nft.eq.0 .and. iand(ndepth,32).eq.32) then
     qmin=2.0 - 0.1*naggressive
     call timer('hint65  ',0)
     call hint65(s3,mrs,mrs2,nadd,nflip,mycall,hiscall,hisgrid,qual,decoded)
     if(qual.ge.qmin) then
        nft=2
        ncount=0
     else
        decoded='                      '
        ntry=0
     endif
     call timer('hint65  ',1)
     go to 900
  endif

  ncount=-1
  decoded='                      '
  ltext=.false.
  if(nft.gt.0) then
! Turn the corrected symbol array into channel symbols for subtraction;
! pass it back to jt65a via common block "chansyms65".
     do i=1,12
        dat4(i)=correct(13-i)
     enddo
     do i=1,63
       tmp(i)=correct(64-i)
     enddo
     correct(1:63)=tmp(1:63)
     call interleave63(correct,63,1)
     call graycode65(correct,63,1)
     call unpackmsg(dat4,decoded,.false.,'      ')     !Unpack the user message
     ncount=0
     if(iand(dat4(10),8).ne.0) ltext=.true.
  endif
900 continue
  if(nft.eq.1 .and. nhard.lt.0) decoded='                      '

  return
end subroutine extract

subroutine getpp(workdat,p)

  use jt65_mod
  integer workdat(63)
  integer a(63)

  a(1:63)=workdat(63:1:-1)
  call interleave63(a,1)
  call graycode(a,63,1,a)

  psum=0.
  do j=1,63
     i=a(j)+1
     x=s3a(i,j)
     s3a(i,j)=0.
     psum=psum + x
     s3a(i,j)=x
  enddo
  p=psum/63.0

  return
end subroutine getpp


code:
subroutine extract4(sym0,ncount,decoded)

  use packjt
  real sym0(207)
  real sym(207)
  character decoded*22
  character*72 c72
  integer*1 symbol(207)
  integer*1 data1(13)                   !Decoded data (8-bit bytes)
  integer   data4a(9)                   !Decoded data (8-bit bytes)
  integer   data4(12)                   !Decoded data (6-bit bytes)
  integer mettab(-128:127,0:1)          !Metric table
  logical first
  data first/.true./
  save first,mettab,ndelta

  if(first) then
     call getmet4(mettab,ndelta)
     first=.false.
  endif

!### Optimize these params: ...
  amp=30.0
  limit=10000

  ave0=sum(sym0)/207.0
  sym=sym0-ave0
  sq=dot_product(sym,sym)
  rms0=sqrt(sq/206.0)
  sym=sym/rms0

  do j=1,207
     n=nint(amp*sym(j))
     if(n.lt.-127) n=-127
     if(n.gt.127) n=127
     symbol(j)=n
  enddo

  nbits=72
  ncycles=0
  ncount=-1
  decoded='                      '
  call interleave4(symbol(2),-1)          !Remove the interleaving
  call fano232(symbol(2),nbits+31,mettab,ndelta,limit,data1,     &
       ncycles,metric,ncount)
  nlim=ncycles/(nbits+31)

!### Make usage here like that in jt9fano...
  if(ncount.ge.0) then
     do i=1,9
        i4=data1(i)
        if(i4.lt.0) i4=i4+256
        data4a(i)=i4
     enddo
     write(c72,1100) (data4a(i),i=1,9)
1100 format(9b8.8)
     read(c72,1102) data4
1102 format(12b6)

     call unpackmsg(data4,decoded,.false.,'      ')
     if(decoded(1:6).eq.'000AAA') then
!        decoded='***WRONG MODE?***'
        decoded='                      '
        ncount=-1
     endif
  endif

  return
end subroutine extract4


code:
subroutine extractmessage144(decoded,msgreceived,nhashflag,recent_calls,nrecent)
  use iso_c_binding, only: c_loc,c_size_t
  use packjt
  use hashing

  character*22 msgreceived
  character*12 call1,call2
  character*12 recent_calls(nrecent)
  integer*1 decoded(80)
  integer*1, target::  i1Dec8BitBytes(10)
  integer*1 i1hashdec
  integer*4 i4Dec6BitWords(12)

! Collapse 80 decoded bits to 10 bytes. Bytes 1-9 are the message, byte 10 is the hash
  do ibyte=1,10
      itmp=0
    do ibit=1,8
      itmp=ishft(itmp,1)+iand(1,decoded((ibyte-1)*8+ibit))
    enddo
    i1Dec8BitBytes(ibyte)=itmp
  enddo

! Calculate the hash using the first 9 bytes.
  ihashdec=nhash(c_loc(i1Dec8BitBytes),int(9,c_size_t),146)
  ihashdec=2*iand(ihashdec,255)

! Compare calculated hash with received byte 10 - if they agree, keep the message.
  i1hashdec=ihashdec
  if( i1hashdec .eq. i1Dec8BitBytes(10) ) then
! Good hash --- unpack 72-bit message
    do ibyte=1,12
      itmp=0
      do ibit=1,6
        itmp=ishft(itmp,1)+iand(1,decoded((ibyte-1)*6+ibit))
      enddo
      i4Dec6BitWords(ibyte)=itmp
    enddo
    call unpackmsg144(i4Dec6BitWords,msgreceived,call1,call2)
    nhashflag=1
    if( call1(1:2) .ne. 'CQ' .and. call1(1:2) .ne. '  ' ) then
      call update_recent_calls(call1,recent_calls,nrecent)
    endif
    if( call2(1:2) .ne. '  ' ) then
      call update_recent_calls(call2,recent_calls,nrecent)
    endif
  else
    msgreceived=' '
    nhashflag=-1
  endif
  return
  end subroutine extractmessage144


code:
subroutine fano232(symbol,nbits,mettab,ndelta,maxcycles,dat,     &
     ncycles,metric,ierr)

! Sequential decoder for K=32, r=1/2 convolutional code using
! the Fano algorithm.  Translated from C routine for same purpose
! written by Phil Karn, KA9Q.

  parameter (MAXBITS=103)
  parameter (MAXBYTES=(MAXBITS+7)/8)
  integer*1 symbol(0:2*MAXBITS-1)  !Soft symbols (as unsigned i*1)
  integer*1 dat(MAXBYTES)          !Decoded user data, 8 bits per byte
  integer mettab(-128:127,0:1)     !Metric table

! These were the "node" structure in Karn's C code:
  integer nstate(0:MAXBITS)        !Encoder state of next node
  integer gamma(0:MAXBITS)         !Cumulative metric to this node
  integer metrics(0:3,0:MAXBITS)   !Metrics indexed by all possible Tx syms
  integer tm(0:1,0:MAXBITS)        !Sorted metrics for current hypotheses
  integer ii(0:MAXBITS)            !Current branch being tested

  logical noback
  include 'conv232.f90'            !Polynomials defined here

  ntail=nbits-31

! Compute all possible branch metrics for each symbol pair.
! This is the only place we actually look at the raw input symbols
  i4a=0
  i4b=0
  do np=0,nbits-1
     j=2*np
     i4a=symbol(j)
     i4b=symbol(j+1)
     metrics(0,np) = mettab(i4a,0) + mettab(i4b,0)
     metrics(1,np) = mettab(i4a,0) + mettab(i4b,1)
     metrics(2,np) = mettab(i4a,1) + mettab(i4b,0)
     metrics(3,np) = mettab(i4a,1) + mettab(i4b,1)
  enddo

  np=0
  nstate(np)=0

  n=iand(nstate(np),npoly1)                  !Compute and sort branch metrics
  n=ieor(n,ishft(n,-16))                     !from the root node
  lsym=partab(iand(ieor(n,ishft(n,-8)),255))
  n=iand(nstate(np),npoly2)
  n=ieor(n,ishft(n,-16))
  lsym=lsym+lsym+partab(iand(ieor(n,ishft(n,-8)),255))
  m0=metrics(lsym,np)
  m1=metrics(ieor(3,lsym),np)
  if(m0.gt.m1) then
     tm(0,np)=m0                             !0-branch has better metric
     tm(1,np)=m1
  else
     tm(0,np)=m1                             !1-branch is better
     tm(1,np)=m0
     nstate(np)=nstate(np) + 1               !Set low bit
  endif

  ii(np)=0                                   !Start with best branch
  gamma(np)=0
  nt=0

  do i=1,nbits*maxcycles                     !Start the Fano decoder
     ngamma=gamma(np) + tm(ii(np),np)        !Look forward
     if(ngamma.ge.nt) then
! Node is acceptable.  If first time visiting this node, tighten threshold:
        if(gamma(np).lt.(nt+ndelta)) nt=nt + ndelta * ((ngamma-nt)/ndelta)
        gamma(np+1)=ngamma                   !Move forward
        nstate(np+1)=ishft(nstate(np),1)
        np=np+1
        if(np.eq.nbits) go to 100          !We're done!

        n=iand(nstate(np),npoly1)
        n=ieor(n,ishft(n,-16))
        lsym=partab(iand(ieor(n,ishft(n,-8)),255))
        n=iand(nstate(np),npoly2)
        n=ieor(n,ishft(n,-16))
        lsym=lsym+lsym+partab(iand(ieor(n,ishft(n,-8)),255))

        if(np.ge.ntail) then
           tm(0,np)=metrics(lsym,np)      !We're in the tail, now all zeros
        else
           m0=metrics(lsym,np)
           m1=metrics(ieor(3,lsym),np)
           if(m0.gt.m1) then
              tm(0,np)=m0                 !0-branch has better metric
              tm(1,np)=m1
           else
              tm(0,np)=m1                 !1-branch is better
              tm(1,np)=m0
              nstate(np)=nstate(np) + 1   !Set low bit
           endif
        endif
        ii(np)=0                          !Start with best branch
     else
        do while(.true.)
           noback=.false.                 !Threshold violated, can't go forward
           if(np.eq.0) noback=.true.
           if(np.gt.0) then
              if(gamma(np-1).lt.nt) noback=.true.
           endif

           if(noback) then               !Can't back up, either
              nt=nt-ndelta               !Relax threshold and look forward again
              if(ii(np).ne.0) then
                 ii(np)=0
                 nstate(np)=ieor(nstate(np),1)
              endif
              exit
           endif

           np=np-1                       !Back up
           if(np.lt.ntail .and. ii(np).ne.1) then
              ii(np)=ii(np)+1            !Search the next best branch
              nstate(np)=ieor(nstate(np),1)
              exit
           endif
        enddo
     endif
  enddo
  i=nbits*maxcycles

100 metric=gamma(np)                       !Final path metric
  nbytes=(nbits+7)/8                       !Copy decoded data to user's buffer
  np=7
  do j=1,nbytes-1
     i4a=nstate(np)
     dat(j)=i4a
     np=np+8
  enddo
  dat(nbytes)=0
  ncycles=i+1
  ierr=0
  if(i.ge.maxcycles*nbits) ierr=-1

  return
end subroutine fano232


code:
subroutine fast9(id2,narg,line)

! Decoder for "fast9" modes, JT9E to JT9H.

  parameter (NMAX=30*12000,NSAVE=500)
  integer*2 id2(0:NMAX)
  integer narg(0:14)
  integer*1 i1SoftSymbols(207)
  integer*1 i1save(207,NSAVE)
  integer indx(NSAVE)
  integer*8 count0,count1,clkfreq
  real s1(720000)                      !To reserve space.  Logically s1(nq,jz)
  real s2(240,340)                     !Symbol spectra at quarter-symbol steps
  real ss2(0:8,85)                     !Folded symbol spectra
  real ss3(0:7,69)                     !Folded spectra without sync symbols
  real s(1500)
  real ccfsave(NSAVE)
  real t0save(NSAVE)
  real t1save(NSAVE)
  real freqSave(NSAVE)
  real t(6)
  character*22 msg                     !Decoded message
  character*80 line(100)
  data nsubmode0/-1/,ntot/0/
  save s1,nsubmode0,ntot

! Parameters from GUI are in narg():
  nutc=narg(0)                         !UTC
  npts=min(narg(1),NMAX)               !Number of samples in id2 (12000 Hz)
  nsubmode=narg(2)                     !0=A 1=B 2=C 3=D 4=E 5=F 6=G 7=H
  if(nsubmode.lt.4) go to 900
  newdat=narg(3)                       !1==> new data, compute symbol spectra
  minsync=narg(4)                      !Lower sync limit
  npick=narg(5)
  t0=0.001*narg(6)
  t1=0.001*narg(7)
  maxlines=narg(8)                     !Max # of decodes to return to caller
  nmode=narg(9)
  nrxfreq=narg(10)                     !Targer Rx audio frequency (Hz)
  ntol=narg(11)                        !Search range, +/- ntol (Hz)

  tmid=npts*0.5/12000.0
  line(1:100)(1:1)=char(0)
  s=0
  s2=0
  nsps=60 * 2**(7-nsubmode)            !Samples per sysbol
  nfft=2*nsps                          !FFT size
  nh=nfft/2
  nq=nfft/4
  istep=nsps/4                         !Symbol spectra at quarter-symbol steps
  jz=npts/istep
  df=12000.0/nfft                      !FFT bin width
  db1=db(2500.0/df)
  nfa=max(200,nrxfreq-ntol)            !Lower frequency limit
  nfb=min(nrxfreq+ntol,2500)           !Upper frequency limit
  nline=0
  t=0.

  if(newdat.eq.1 .or. nsubmode.ne.nsubmode0) then
     call system_clock(count0,clkfreq)
     call spec9f(id2,npts,nsps,s1,jz,nq)          !Compute symbol spectra, s1
     call system_clock(count1,clkfreq)
     t(1)=t(1)+float(count1-count0)/float(clkfreq)
  endif

  nsubmode0=nsubmode
  tmsg=nsps*85.0/12000.0
  limit=2000
  nlen0=0
  i1=0
  i2=0
  ccfsave=0.
  do ilength=1,14
     nlen=1.4142136**(ilength-1)
     if(nlen.gt.jz/340) nlen=jz/340
     if(nlen.eq.nlen0) cycle
     nlen0=nlen

     db0=db(float(nlen))
     jlen=nlen*340
     jstep=jlen/4                      !### Is this about right? ###
     if(nsubmode.ge.6) jstep=jlen/2

     do ja=1,jz-jlen,jstep
        jb=ja+jlen-1
        call system_clock(count0,clkfreq)
        call foldspec9f(s1,nq,jz,ja,jb,s2)        !Fold symbol spectra into s2
        call system_clock(count1,clkfreq)
        t(2)=t(2)+float(count1-count0)/float(clkfreq)

! Find sync; put sync'ed symbol spectra into ss2 and ss3
! Might want to do a peakup in DT and DF, then re-compute symbol spectra.

        call system_clock(count0,clkfreq)
        call sync9f(s2,nq,nfa,nfb,ss2,ss3,lagpk,ipk,ccfbest)
        call system_clock(count1,clkfreq)
        t(3)=t(3)+float(count1-count0)/float(clkfreq)

        i1=i1+1
        if(ccfbest.lt.30.0) cycle
        call system_clock(count0,clkfreq)
        call softsym9f(ss2,ss3,i1SoftSymbols)     !Compute soft symbols
        call system_clock(count1,clkfreq)
        t(4)=t(4)+float(count1-count0)/float(clkfreq)

        i2=i2+1
        ccfsave(i2)=ccfbest
        i1save(1:207,i2)=i1SoftSymbols
        t0=(ja-1)*istep/12000.0
        t1=(jb-1)*istep/12000.0
        t0save(i2)=t0
        t1save(i2)=t1
        freq=ipk*df
        freqSave(i2)=freq
     enddo
  enddo
  nsaved=i2

  ccfsave(1:nsaved)=-ccfsave(1:nsaved)
  call system_clock(count0,clkfreq)
  indx=0
  call indexx(ccfsave,nsaved,indx)
  call system_clock(count1,clkfreq)
  t(5)=t(5)+float(count1-count0)/float(clkfreq)

  ccfsave(1:nsaved)=-ccfsave(1:nsaved)

  do iter=1,2
!     do isave=1,nsaved
     do isave=1,50
        i2=indx(isave)
        if(i2.lt.1 .or. i2.gt.nsaved) cycle      !### Why needed? ###
        t0=t0save(i2)
        t1=t1save(i2)
        if(iter.eq.1 .and. t1.lt.tmid) cycle
        if(iter.eq.2 .and. t1.ge.tmid) cycle
        ccfbest=ccfsave(i2)
        i1SoftSymbols=i1save(1:207,i2)
        freq=freqSave(i2)
        call system_clock(count0,clkfreq)
        call jt9fano(i1SoftSymbols,limit,nlim,msg)      !Invoke Fano decoder
        call system_clock(count1,clkfreq)
        t(6)=t(6)+float(count1-count0)/float(clkfreq)

        i=t0*12000.0
        kz=(t1-t0)/0.02
        smax=0.
        do k=1,kz
           sq=0.
           do n=1,240
              i=i+1
              x=id2(i)
              sq=sq+x*x
           enddo
           s(k)=sq/240.
           smax=max(s(k),smax)
        enddo
        call pctile(s,kz,35,base)
        snr=smax/(1.1*base) - 1.0
        nsnr=-20
        if(snr.gt.0.0) nsnr=nint(db(snr))

!        write(72,3002) nutc,iter,isave,nlen,tmid,t0,t1,ccfbest,   &
!             nint(freq),nlim,msg
!3002    format(i6.6,i1,i4,i3,4f6.1,i5,i7,1x,a22)

        if(msg.ne.'                      ') then

! Display multiple decodes only if they differ:
           do n=1,nline
              if(index(line(n),msg).gt.1) go to 100
           enddo
!### Might want to use decoded message to get a complete estimate of S/N.
           nline=nline+1
           write(line(nline),1000) nutc,nsnr,t0,nint(freq),msg,char(0)
1000       format(i6.6,i4,f5.1,i5,1x,'@ ',1x,a22,a1)
           ntot=ntot+1
!           write(70,5001) nsaved,isave,nline,maxlines,ntot,nutc,msg
!5001       format(5i5,i7.6,1x,a22)
           if(nline.ge.maxlines) go to 900
        endif
100     continue
     enddo
  enddo

900 continue
!  write(*,6001) t,t(6)/sum(t)
!6001 format(7f10.3)

  return
end subroutine fast9


code:
subroutine fast_decode(id2,narg,ntrperiod,line,mycall_12,   &
     hiscall_12)

  parameter (NMAX=30*12000)
  integer*2 id2(NMAX)
  integer*2 id2a(NMAX)
  integer*2 id2b(NMAX)
  integer narg(0:14)
  real dat(30*12000)
  complex cdat(262145),cdat2(262145)
  real psavg(450)
  logical pick,first
  character*6 cfile6
  character*80 line(100)
  character*12 mycall_12,hiscall_12
  character*6 mycall,hiscall
  data first/.true./,nutca/0/,nutcb/0/
  save npts,cdat,cdat2,id2a,id2b,nutca,nutcb

  if(first) then
     id2a=0
     id2b=0
     first=.false.
  endif

  mycall=mycall_12(1:6)
  hiscall=hiscall_12(1:6)
  nutc=narg(0)
  ndat0=narg(1)
  nsubmode=narg(2)
  newdat=narg(3)
  minsync=narg(4)
  npick=narg(5)
  t0=0.001*narg(6)
  t1=0.001*narg(7)
  tmid=0.5*(t0+t1)
  maxlines=narg(8)
  nmode=narg(9)
  nrxfreq=narg(10)
  ntol=narg(11)
  nhashcalls=narg(12)

  line(1:100)(1:1)=char(0)
  if(t0.gt.float(ntrperiod)) go to 900
  if(t0.gt.t1) go to 900

  if(nmode.eq.102) then
     call fast9(id2,narg,line)
     go to 900
  endif

  if(newdat.eq.1) then
     cdat2=cdat
     ndat=ndat0
     call wav11(id2,ndat,dat)
     nzz=11025*ntrperiod
     if(ndat.lt.nzz) dat(ndat+1:nzz)=0.0
     ndat=min(ndat,30*11025)
     call ana932(dat,ndat,cdat,npts)          !Make downsampled analytic signal
  endif

! Now cdat() is the downsampled analytic signal.
! New sample rate = fsample = BW = 11025 * (9/32) = 3100.78125 Hz
! NB: npts, nsps, etc., are all reduced by 9/32

  write(cfile6,'(i6.6)') nutc
  nfreeze=1
  mousedf=0
  mousebutton=0
  mode4=1
  if(nsubmode.eq.1) mode4=2
  nafc=0
  ndebug=0
  t2=0.
  ia=1
  ib=npts
  pick=.false.

  if(npick.gt.0) then
     pick=.true.
     dt=1.0/11025.0 * (32.0/9.0)
     ia=t0/dt + 1.
     ib=t1/dt + 1.
     t2=t0
  endif
  jz=ib-ia+1
  line(1:100)(1:1)=char(0)
  if(npick.eq.2) then
     call iscat(cdat2(ia),jz,3,40,t2,pick,cfile6,minsync,ntol,NFreeze,    &
          MouseDF,mousebutton,mode4,nafc,ndebug,psavg,nmax,nlines,line)
  else
     call iscat(cdat(ia),jz,3,40,t2,pick,cfile6,minsync,ntol,NFreeze,     &
          MouseDF,mousebutton,mode4,nafc,ndebug,psavg,maxlines,nlines,line)
  endif

900 return
end subroutine fast_decode


code:
program fcal

! Compute Intercept (A) and Slope (B) for a series of FreqCal measurements.
  parameter(NZ=1000)
  implicit real*8 (a-h,o-z)
  real*8 fd(NZ),deltaf(NZ),r(NZ)
  character infile*50
  character line*80
  character cutc*8

  nargs=iargc()
  if(nargs.ne.1) then
     print*,'Usage:   fcal <infile>'
     print*,'Example: fcal fmtave.out'
     go to 999
  endif
  call getarg(1,infile)

  open(10,file=infile,status='old',err=997)
  open(12,file='fcal.out',status='unknown')
  open(13,file='fcal.plt',status='unknown')

  i=0
  do j=1,9999
     read(10,1000,end=10) line
1000 format(a80)
     i0=index(line,' 0 ')
     i1=index(line,' 1 ')
     if(i0.le.0 .and. i1.le.0) then
        read(line,*,err=5) f,df
        ncal=1
        i=i+1
        fd(i)=f
        deltaf(i)=df
     else if(i1.gt.0) then
        i=i+1
        read(line,*,err=5) f,df,ncal,nn,rr,cutc
        fd(i)=f
        deltaf(i)=df
        r(i)=0.d0
     endif
5    continue
  enddo

10 iz=i
  if(iz.lt.2) go to 998
  call fit(fd,deltaf,r,iz,a,b,sigmaa,sigmab,rms)

  write(*,1002)
1002 format('    Freq      DF     Meas Freq     Resid'/        &
            '   (MHz)     (Hz)      (MHz)        (Hz)'/        &
            '-----------------------------------------')
  do i=1,iz
     fm=fd(i) + 1.d-6*deltaf(i)
     calfac=1.d0 + 1.d-6*deltaf(i)/fd(i)
     write(*,1010) fd(i),deltaf(i),fm,r(i)
     write(13,1010) fd(i),deltaf(i),fm,r(i)
1010 format(f8.3,f9.3,f14.9,f9.3,2x,a6)
  enddo
  calfac=1.d0 + 1.d-6*b
  err=1.d-6*sigmab

  if(iz.ge.3) then
     write(*,1100) a,b,rms
1100 format(/'A:',f8.2,' Hz    B:',f9.4,' ppm    StdDev:',f7.3,' Hz')
  if(iz.gt.2) write(*,1110) sigmaa,sigmab
1110 format('err:',f6.2,9x,f9.4,23x,f13.9)
  else
     write(*,1120) a,b
1120 format(/'A:',f8.2,' Hz    B:',f9.4)
  endif

  write(12,1130) a,b
1130 format(f10.4)

  go to 999

997 print*,'Cannot open input file: ',infile
  go to 999
998 print*,'Input file must contain at least 2 valid measurement pairs'

999 end program fcal

subroutine fit(x,y,r,iz,a,b,sigmaa,sigmab,rms)
  implicit real*8 (a-h,o-z)
  real*8 x(iz),y(iz),r(iz)

  sx=0.d0
  sy=0.d0
  sxy=0.d0
  sx2=0.d0
  do i=1,iz
     sx=sx + x(i)
     sy=sy + y(i)
     sxy=sxy + x(i)*y(i)
     sx2=sx2 + x(i)*x(i)
  enddo
  delta=iz*sx2 - sx*sx
  a=(sx2*sy - sx*sxy)/delta
  b=(iz*sxy - sx*sy)/delta

  sq=0.d0
  do i=1,iz
     r(i)=y(i) - (a + b*x(i))
     sq=sq + r(i)**2
  enddo
  rms=0.
  sigmaa=0.
  sigmab=0.
  if(iz.ge.3) then
     rms=sqrt(sq/(iz-2))
     sigmaa=sqrt(rms*rms*sx2/delta)
     sigmab=sqrt(iz*rms*rms/delta)
  endif

  return
end subroutine fit


code:
real function fchisq(c3,npts,fsample,a)

  parameter (NMAX=85*16)
  complex c3(npts)
  complex c4(NMAX)
  real a(3)
  complex z
  data a1,a2,a3/99.,99.,99./
  include 'jt9sync.f90'
  save

  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)
     call twkfreq(c3,c4,npts,fsample,a)
  endif

! Get sync power.
  nspsd=16
  sum1=0.
  sum0=0.
  k=-1
  do i=1,85
     z=0.
     do j=1,nspsd
        k=k+1
        z=z+c4(k+1)
     enddo
     pp=real(z)**2 + aimag(z)**2
     if(isync(i).eq.1) then
        sum1=sum1+pp
     else
        sum0=sum0+pp
     endif
  enddo
  sync_4992=(sum1/16.0)/(sum0/69.0) - 1.0             !r4992
  sync_4993=sum1/10000.0                              !r4993+
!  write(80,3001) 1.e-5*sum1,1.e-5*sum0,sync_4992,sync_4993,sync
!3001 format(5f11.4)
  fchisq=-sync_4993

  return
end function fchisq


code:
real function fchisq0(y,npts,a)

  real y(npts),a(4)

!  rewind 51
  chisq = 0.
  do i=1,npts
     x=i
     z=(x-a(3))/(0.5*a(4))
     yfit=a(1)
     if(abs(z).lt.3.0) then
        d=1.0 + z*z
        yfit=a(1) + a(2) * (1.0/d - 0.1)
     endif
     chisq=chisq + (y(i) - yfit)**2
!     write(51,3001) i,y(i),yfit,y(i)-yfit
!3001 format(i5,3f10.4)
  enddo
  fchisq0=chisq

  return
end function fchisq0


code:
real function fchisq65(cx,npts,fsample,nflip,a,ccfmax,dtmax)

  use timer_module, only: timer

  parameter (NMAX=60*12000)          !Samples per 60 s
  complex cx(npts)
  real a(5)
  complex w,wstep,z
  real ss(3000)
  complex csx(0:NMAX/8)
  data twopi/6.283185307/a1,a2,a3/99.,99.,99./
  save

  call timer('fchisq65',0)
  baud=11025.0/4096.0
  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

 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)

! Mix and integrate the complex signal
     csx(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
           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)
     enddo
  endif

! Compute whole-symbol powers at 1/16-symbol steps.
  fac=1.e-4
  do i=1,nout
     j=nsps+(i-1)*nsps/16 !steps by 8 samples (1/16 of a symbol)
     k=j-nsps
     ss(i)=0.
     if(k.ge.0 .and. j.le.npts) then
        z=csx(j)-csx(k) ! difference over span of 128 pts
        ss(i)=fac*(real(z)**2 + aimag(z)**2)
     endif
  enddo

  ccfmax=0.
  call timer('ccf2    ',0)
  call ccf2(ss,nout,nflip,ccf,xlagpk)
  call timer('ccf2    ',1)
  if(ccf.gt.ccfmax) then
     ccfmax=ccf
     dtmax=xlagpk*dtstep
  endif
  fchisq65=-ccfmax
  call timer('fchisq65',1)

  return
end function fchisq65


code:
#!/bin/bash

# jt65 end-to-end simulator. Run enough cases to get a given
# number of successful decodes or a given number of errors
# For better performance, copy this script and jt65, jt65sim
# executables to a ramdisk and run from there.

mode=C
spread=4.0

for i in `seq 0 21`;
do
nruns=0
ndecodes=0
ngood=0
nbad=0
snr=$((-25+$i))
#while [ $[ $nruns - $ngood ] -lt 200 ]
#while [ $ngood -lt 200 ]
while [ $nruns -lt 1000 ]
do
  rm 000000_0001.wav
  ./jt65sim -m $mode -n 1 -d $spread -s \\$snr > /dev/null
  rm decoded.txt
  ./jt65 -m $mode -a 10 -f 1500 -n 1000 -s 000000_0001.wav > decoded.txt
  if [ -e decoded.txt ]; then
    nd=$( wc -l < decoded.txt )
    ng=$( grep "K1ABC" decoded.txt | wc -l )
    nb=$[$nd-$ng]
    ndecodes=$[$ndecodes + $nd]
    ngood=$[$ngood + $ng]
    nbad=$[$nbad + $nb]
  fi
  nruns=$[$nruns + 1]
  nmissed=$[$nruns-$ngood]
  pe=$( echo "scale=5; $nmissed/$nruns" | bc -l )
  ps=$( echo "scale=5; $ngood/$nruns" | bc -l )
  echo -en snr: $snr nruns: $nruns ndecodes: $ndecodes ngood: $ngood nbad: $nbad nmissed: $nmissed wer: $pe success: $ps"\r"
done
echo snr: $snr nruns: $nruns ndecodes: $ndecodes ngood: $ngood nbad: $nbad nmissed: $nmissed wer: $pe success: $ps
done

path: lib/fer


code:
program fer65

! End-to-end simulator for testing JT65.

! Options
!  jt65sim                             jt65
!----------------------------------------------------------------
!                                      -a aggressive
!  -d Doppler spread                   -d depth
!  -f Number of files                  -f freq
!  -m (sub)mode                        -m (sub)mode
!  -n number of generated sigs         -n ntrials
!  -t Time offset (s)                  -r robust sync
!  -p Do not seed random #s            -c mycall
!                                      -x hiscall
!                                      -g hisgrid
!                                      -X hinted-decode flags
!  -s S/N in 2500 Hz                   -s single-decode mode

  implicit real*8 (a-h,o-z)
  real*8 s(7),sq(7)
  character arg*12,cmnd*100,decoded*22,submode*1,csync*1,f1*15,f2*15
  character*12 outfile
  logical syncok

  nargs=iargc()
  if(nargs.ne.7) then
     print*,'Usage:   fer65 submode fspread snr1 snr2 Navg  DS  iters'
     print*,'Example: fer65    C      3.0   -28  -12    8    1  1000'
     go to 999
  endif

  call getarg(1,submode)
  call getarg(2,arg)
  read(arg,*) d
  call getarg(3,arg)
  read(arg,*) snr1
  call getarg(4,arg)
  read(arg,*) snr2
  call getarg(5,arg)
  read(arg,*) navg
  call getarg(6,arg)
  read(arg,*) nds
  call getarg(7,arg)
  read(arg,*) iters

  write(outfile,1001) submode,d,navg,nds
1001 format(a1,f6.2,'_',i2.2,'_',i1)
  if(outfile(2:2).eq.' ') outfile(2:2)='0'
  if(outfile(3:3).eq.' ') outfile(3:3)='0'

  ndepth=3
  if(navg.gt.1) ndepth=ndepth+16
  if(nds.ne.0) ndepth=ndepth+32

  dfmax=3
  if(submode.eq.'b' .or. submode.eq.'B') dfmax=6
  if(submode.eq.'c' .or. submode.eq.'C') dfmax=11

  ntrials=1000
  naggressive=10

  open(20,file=outfile,status='unknown')
  open(21,file='fer65.21',status='unknown')

  write(20,1000) submode,iters,ntrials,naggressive,d,ndepth,navg,nds
1000 format(/'JT65',a1,'  Iters:',i5,'  T:',i6,'  Aggr:',i3,  &
          '  Dop:',f6.2,'  Depth:',i2,'  Navg:',i3,'  DS:',i2)
  write(20,1002)
1002 format(/'  dB  nsync ngood nbad     sync       dsnr        ',     &
            'DT       Freq      Nsum     Width'/85('-'))
  flush(20)

  do isnr=0,20
     snr=snr1+isnr
     if(snr.gt.snr2) exit
     nsync=0
     ngood=0
     nbad=0
     s=0.
     sq=0.
     do iter=1,iters
        write(cmnd,1010) submode,d,snr,navg
1010    format('./jt65sim -n 1 -m ',a1,' -d',f7.2,' -s \\',f5.1,' -f',i3,' >devnull')
        call unlink('000000_????.wav')
        call system(cmnd)
        if(navg.gt.1) then
           do i=navg,2,-1
              j=2*i-1
              write(f1,1011) i
              write(f2,1011) j
1011          format('000000_',i4.4,'.wav')
              call rename(f1,f2)
           enddo
        endif
        call unlink('decoded.txt')
        call unlink('fort.13')
        isync=0
        nsnr=0
        dt=0.
        nfreq=0
        ndrift=0
        nwidth=0
        cmnd='./jt65 -m A -a 10 -c K1ABC -f 1500 -n 1000 -d  5 -s 000000_????.wav > decoded.txt'
        cmnd(11:11)=submode
        write(cmnd(47:48),'(i2)') ndepth
        call system(cmnd)
        open(13,file='fort.13',status='old',err=20)
        do i=1,navg
           read(13,1012) nutc,isync,nsnr,dt,nfreq,ndrift,nwidth,decoded,     &
             nft,nsum,nsmo
1012       format(i4,i4,i5,f6.2,i5,i4,i3,1x,a22,5x,3i3)
           if(nft.gt.0) exit
        enddo
        close(13)
        syncok=abs(dt).lt.0.2 .and. float(abs(nfreq-1500)).lt.dfmax
        csync=' '
        if(syncok) csync='*'
        write(21,1014) nutc,isync,nsnr,dt,nfreq,ndrift,nwidth,     &
             nft,nsum,nsmo,csync,decoded(1:16),nft,nsum,nsmo
1014    format(i4,i4,i5,f6.2,i5,i4,3x,4i3,1x,a1,1x,a16,i2,2i3)
        flush(21)

        if(syncok) then
           nsync=nsync+1
           s(1)=s(1) + isync
           sq(1)=sq(1) + isync*isync
           s(6)=s(6) + nwidth
           sq(6)=sq(6) + nwidth*nwidth
           if(decoded.eq.'K1ABC W9XYZ EN37      ') then
              ngood=ngood+1
              s(2)=s(2) + nsnr
              s(3)=s(3) + dt
              s(4)=s(4) + nfreq
              s(5)=s(5) + ndrift
              s(7)=s(7) + nsum

              sq(2)=sq(2) + nsnr*nsnr
              sq(3)=sq(3) + dt*dt
              sq(4)=sq(4) + nfreq*nfreq
              sq(5)=sq(5) + ndrift*ndrift
              sq(7)=sq(7) + nsum*nsum
           else if(decoded.ne.'                      ') then
              nbad=nbad+1
              print*,'Nbad:',nbad,decoded
           endif
        endif
20      continue
        fsync=float(nsync)/iter
        fgood=float(ngood)/iter
        fbad=float(nbad)/iter
        write(*,1020) nint(snr),iter,isync,nsnr,dt,nfreq,ndrift,nwidth,fsync,  &
             fgood,fbad,decoded(1:16),nft,nsum,nsmo
1020    format(i3,i5,i3,i4,f6.2,i5,i3,i3,2f6.3,f7.4,1x,a16,i2,2i3)
     enddo

     if(nsync.ge.1) then
        xsync=s(1)/nsync
        xwidth=s(6)/nsync
     endif
     esync=0.
     if(nsync.ge.2) then
        esync=sqrt(sq(1)/nsync - xsync**2)
        ewidth=sqrt(sq(6)/nsync - xwidth**2)
     endif

     if(ngood.ge.1) then
        xsnr=s(2)/ngood
        xdt=s(3)/ngood
        xfreq=s(4)/ngood
        xdrift=s(5)/ngood
        xsum=s(7)/ngood
     endif
     if(ngood.ge.2) then
        esnr=sqrt(sq(2)/ngood - xsnr**2)
        edt=sqrt(sq(3)/ngood - xdt**2)
        efreq=sqrt(sq(4)/ngood - xfreq**2)
        edrift=sqrt(sq(5)/ngood - xdrift**2)
        esum=sqrt(sq(7)/ngood - xsum**2)
     endif

     dsnr=xsnr-snr
     dfreq=xfreq-1500.0
     if(ngood.eq.0) then
        dsnr=0.
        dfreq=0.
     endif
     write(20,1100) snr,nsync,ngood,nbad,xsync,esync,dsnr,esnr,  &
          xdt,edt,dfreq,efreq,xsum,esum,xwidth,ewidth
1100 format(f5.1,2i6,i4,2f6.1,f6.1,f5.1,f6.2,f5.2,6f5.1)
     flush(20)
     if(ngood.ge.int(0.99*iters)) exit
  enddo

999 end program fer65


code:
program fersum

  character mode*5
  character infile*40
  real dop(0:9)
  real thresh(0:9,12),threshsync(0:9,12)
  data dop/0.25,0.5,1.0,2.0,4.0,8.0,16.0,32.0,64.0,128.0/

  nargs=iargc()
  if(nargs.ne.1) then
     print*,'Usage: fersum <infile>'
     go to 999
  endif
  call getarg(1,infile)
  open(10,file=infile,status='old')
  thresh=0.
  threshsync=0.

  do iblk=1,999
1    read(10,1002,end=100) mode,iters,ntot,naggr,d,navg,nds
1002 format(a5,8x,i5,4x,i6,7x,i3,6x,f6.2,17x,i3,5x,i2)
     write(33,*) iblk,mode
     if(mode.eq.'     ') go to 1
     read(10,1002)
     read(10,1002)
     read(10,1002)

     nsync0=0
     ngood0=0
     xsum0=0.
     do n=1,99
        read(10,1010,end=100) snr,nsync,ngood,nbad,xsync,esync,dsnr,esnr,   &
             xdt,edt,dfreq,efreq,xsum,esum,xwidth,ewidth
1010    format(f5.1,2i6i4,2f6.1,f6.1,f5.1,f6.2,f5.2,6f5.1)
        if(snr.eq.0.0) exit
        if(mode(5:5).eq.'A') nmode=1
        if(mode(5:5).eq.'B') nmode=2
        if(mode(5:5).eq.'C') nmode=3
        j=nint(log(d)/log(2.0) + 2.0)
        if(navg.eq.1 .and. nds.eq.0) k=nmode
        if(navg.eq.1 .and. nds.eq.1) k=nmode+3
        if(navg.gt.1 .and. nds.eq.0) k=nmode+6
        if(navg.gt.1 .and. nds.eq.1) k=nmode+9
        if(nsync0.le.iters/2 .and. nsync.ge.iters/2) then
           threshsync(j,k)=snr-float(nsync-iters/2)/(nsync-nsync0)
        endif
        if(ngood0.le.iters/2 .and. ngood.ge.iters/2) then
           threshold=snr-float(ngood-iters/2)/(ngood-ngood0)
           xsumavg=max(1.0,0.5*(xsum0+xsum))
!           write(*,1020) mode,iters,ntot,naggr,d,navg,nds,threshold,xsumavg
!1020       format(a5,i7,i7,i3,f7.2,i3,i3,f7.1,f6.1)
           thresh(j,k)=threshold
        endif
        nsync0=nsync
        ngood0=ngood
        xsum0=xsum
     enddo
  enddo

100 write(12,1100)
1100 format(' ')
  do i=0,9
     write(12,1110) dop(i),thresh(i,1:12)
1110 format(f6.2,13f6.1)
  enddo

  write(12,1110)
  do i=0,9
     write(12,1110) dop(i),threshsync(i,1:12)
  enddo

999 end program fersum