code:
C++++++++++++++++++++++++
C
      SUBROUTINE FSIZER1(NRECL,KSIZE,NRFILE,NAMFIL)
C
C++++++++++++++++++++++++
C
C   Version 1.0 uses the INQUIRE statement to find out the the record length
C   of the direct access file before opening it.  This procedure is non-standard,
C   but seems to work for VAX machines.
C
C  THE SUBROUTINE ALSO SETS THE VALUES OF  NRECL, NRFILE, AND NAMFIL.

C  *****************************************************************
C  *****************************************************************
C
C  THE PARAMETERS NAMFIL, NRECL, AND NRFILE ARE TO BE SET BY THE USER
C
C  *****************************************************************

C  NAMFIL IS THE EXTERNAL NAME OF THE BINARY EPHEMERIS FILE

      CHARACTER*256  NAMFIL

c      NAMFIL='JPLEPH'

C  *****************************************************************

C  NRECL=1 IF "RECL" IN THE OPEN STATEMENT IS THE RECORD LENGTH IN S.P. WORDS
C  NRECL=4 IF "RECL" IN THE OPEN STATEMENT IS THE RECORD LENGTH IN BYTES
C  (for a VAX, it is probably 1)
C
      NRECL=4

C  *****************************************************************

C  NRFILE IS THE INTERNAL UNIT NUMBER USED FOR THE EPHEMERIS FILE

c      NRFILE=12

C  *****************************************************************

C   FIND THE RECORD SIZE USING THE INQUIRE STATEMENT

c      IRECSZ=0

      INQUIRE(FILE=NAMFIL,RECL=IRECSZ)

C IF 'INQUIRE' DOES NOT WORK, USUALLY IRECSZ WILL BE LEFT AT 0

      IF(IRECSZ .LE. 0) write(*,*)
     . ' INQUIRE STATEMENT PROBABLY DID NOT WORK'

      KSIZE=IRECSZ/NRECL
      if(nrfile.eq.-99) stop              !silence compiler warning

      RETURN

      END
C++++++++++++++++++++++++
C
      SUBROUTINE FSIZER2(NRECL,KSIZE,NRFILE,NAMFIL)
C
C++++++++++++++++++++++++
C  THIS SUBROUTINE OPENS THE FILE, 'NAMFIL', WITH A PHONY RECORD LENGTH, READS
C  THE FIRST RECORD, AND USES THE INFO TO COMPUTE KSIZE, THE NUMBER OF SINGLE
C  PRECISION WORDS IN A RECORD.
C
C  THE SUBROUTINE ALSO SETS THE VALUES OF  NRECL, NRFILE, AND NAMFIL.

      IMPLICIT DOUBLE PRECISION(A-H,O-Z)

      SAVE

      INTEGER OLDMAX
      PARAMETER (OLDMAX = 400)
      INTEGER NMAX
      PARAMETER (NMAX = 1000)
      CHARACTER*6 TTL(14,3),CNAM(NMAX)
      CHARACTER*256 NAMFIL,jpleph_file_name
      DIMENSION SS(3)
      INTEGER IPT(3,13)
      common/jplcom/jpleph_file_name

C  *****************************************************************
C  *****************************************************************
C
C  THE PARAMETERS NRECL, NRFILE, AND NAMFIL ARE TO BE SET BY THE USER
C
C  *****************************************************************

C  NRECL=1 IF "RECL" IN THE OPEN STATEMENT IS THE RECORD LENGTH IN S.P. WORDS
C  NRECL=4 IF "RECL" IN THE OPEN STATEMENT IS THE RECORD LENGTH IN BYTES
C  (for UNIX, it is probably 4)
C
      NRECL=4

C  NRFILE IS THE INTERNAL UNIT NUMBER USED FOR THE EPHEMERIS FILE

      NRFILE=12

C  NAMFIL IS THE EXTERNAL NAME OF THE BINARY EPHEMERIS FILE

!      NAMFIL='JPLEPH'
      NAMFIL=jpleph_file_name

C  *****************************************************************
C  *****************************************************************

C  **  OPEN THE DIRECT-ACCESS FILE AND GET THE POINTERS IN ORDER TO
C  **  DETERMINE THE SIZE OF THE EPHEMERIS RECORD

      MRECL=NRECL*1000

        OPEN(NRFILE,
     *       FILE=NAMFIL,
     *       ACCESS='DIRECT',
     *       FORM='UNFORMATTED',
     *       RECL=MRECL,
     *       STATUS='OLD')

      READ(NRFILE,REC=1)TTL,(CNAM(K),K=1,OLDMAX),SS,NCON,AU,EMRAT,
     & ((IPT(I,J),I=1,3),J=1,12),NUMDE,(IPT(I,13),I=1,3)

      CLOSE(NRFILE)

C  FIND THE NUMBER OF EPHEMERIS COEFFICIENTS FROM THE POINTERS

      KMX = 0
      KHI = 0

      DO I = 1,13
         IF (IPT(1,I) .GT. KMX) THEN
            KMX = IPT(1,I)
            KHI = I
         ENDIF
      ENDDO

      ND = 3
      IF (KHI .EQ. 12) ND=2

      KSIZE = 2*(IPT(1,KHI)+ND*IPT(2,KHI)*IPT(3,KHI)-1)

      RETURN

      END
C++++++++++++++++++++++++
C
      SUBROUTINE FSIZER3(NRECL,KSIZE,NRFILE,NAMFIL)
C
C++++++++++++++++++++++++
C
C  THE SUBROUTINE SETS THE VALUES OF  NRECL, KSIZE, NRFILE, AND NAMFIL.

      SAVE
      CHARACTER*256 NAMFIL,jpleph_file_name
      common/jplcom/jpleph_file_name

C  *****************************************************************
C  *****************************************************************
C
C  THE PARAMETERS NRECL, NRFILE, AND NAMFIL ARE TO BE SET BY THE USER

C  *****************************************************************

C  NRECL=1 IF "RECL" IN THE OPEN STATEMENT IS THE RECORD LENGTH IN S.P. WORDS
C  NRECL=4 IF "RECL" IN THE OPEN STATEMENT IS THE RECORD LENGTH IN BYTES

       NRECL=4

C  *****************************************************************

C  NRFILE IS THE INTERNAL UNIT NUMBER USED FOR THE EPHEMERIS FILE (DEFAULT: 12)

      NRFILE=12

C  *****************************************************************

C  NAMFIL IS THE EXTERNAL NAME OF THE BINARY EPHEMERIS FILE

!      NAMFIL='JPLEPH'
      NAMFIL=jpleph_file_name

C  *****************************************************************

C  KSIZE must be set by the user according to the ephemeris to be read

C  For  de200, set KSIZE to 1652
C  For  de405, set KSIZE to 2036
C  For  de406, set KSIZE to 1456
C  For  de414, set KSIZE to 2036
C  For  de418, set KSIZE to 2036
C  For  de421, set KSIZE to 2036
C  For  de422, set KSIZE to 2036
C  For  de423, set KSIZE to 2036
C  For  de424, set KSIZE to 2036
C  For  de430, set KSIZE to 2036

      KSIZE = 2036

C  *******************************************************************

      RETURN

      END
C++++++++++++++++++++++++++
C
      SUBROUTINE PLEPH ( ET, NTARG, NCENT, RRD )
C
C++++++++++++++++++++++++++
C  NOTE : Over the years, different versions of PLEPH have had a fifth argument:
C  sometimes, an error return statement number; sometimes, a logical denoting
C  whether or not the requested date is covered by the ephemeris.  We apologize
C  for this inconsistency; in this present version, we use only the four necessary
C  arguments and do the testing outside of the subroutine.
C
C     THIS SUBROUTINE READS THE JPL PLANETARY EPHEMERIS
C     AND GIVES THE POSITION AND VELOCITY OF THE POINT 'NTARG'
C     WITH RESPECT TO 'NCENT'.
C
C     CALLING SEQUENCE PARAMETERS:
C
C       ET = D.P. JULIAN EPHEMERIS DATE AT WHICH INTERPOLATION
C            IS WANTED.
C
C       ** NOTE THE ENTRY DPLEPH FOR A DOUBLY-DIMENSIONED TIME **
C          THE REASON FOR THIS OPTION IS DISCUSSED IN THE
C          SUBROUTINE STATE
C
C     NTARG = INTEGER NUMBER OF 'TARGET' POINT.
C
C     NCENT = INTEGER NUMBER OF CENTER POINT.
C
C            THE NUMBERING CONVENTION FOR 'NTARG' AND 'NCENT' IS:
C
C                1 = MERCURY           8 = NEPTUNE
C                2 = VENUS             9 = PLUTO
C                3 = EARTH            10 = MOON
C                4 = MARS             11 = SUN
C                5 = JUPITER          12 = SOLAR-SYSTEM BARYCENTER
C                6 = SATURN           13 = EARTH-MOON BARYCENTER
C                7 = URANUS           14 = NUTATIONS (LONGITUDE AND OBLIQ)
C                            15 = LIBRATIONS, IF ON EPH FILE
C
C             (IF NUTATIONS ARE WANTED, SET NTARG = 14. FOR LIBRATIONS,
C              SET NTARG = 15. SET NCENT=0.)
C
C      RRD = OUTPUT 6-WORD D.P. ARRAY CONTAINING POSITION AND VELOCITY
C            OF POINT 'NTARG' RELATIVE TO 'NCENT'. THE UNITS ARE AU AND
C            AU/DAY. FOR LIBRATIONS THE UNITS ARE RADIANS AND RADIANS
C            PER DAY. IN THE CASE OF NUTATIONS THE FIRST FOUR WORDS OF
C            RRD WILL BE SET TO NUTATIONS AND RATES, HAVING UNITS OF
C            RADIANS AND RADIANS/DAY.
C
C            The option is available to have the units in km and km/sec.
C            For this, set km=.true. in the STCOMX common block.
C

      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      INTEGER NMAX
      PARAMETER (NMAX = 1000)

      DIMENSION RRD(6),ET2Z(2),ET2(2),PV(6,13)
      DIMENSION PVST(6,11),PNUT(4)
      DIMENSION SS(3),CVAL(NMAX),PVSUN(6),ZIPS(2)
      DATA ZIPS/2*0.d0/

      LOGICAL BSAVE,KM,BARY
      LOGICAL FIRST
      DATA FIRST/.TRUE./

      INTEGER LIST(12),IPT(39),DENUM
      COMMON/EPHHDR/CVAL,SS,AU,EMRAT,DENUM,NCON,IPT
      COMMON/STCOMX/KM,BARY,PVSUN

C     INITIALIZE ET2 FOR 'STATE' AND SET UP COMPONENT COUNT
C
      ET2(1)=ET
      ET2(2)=0.D0
      GO TO 11

C     ENTRY POINT 'DPLEPH' FOR DOUBLY-DIMENSIONED TIME ARGUMENT
C          (SEE THE DISCUSSION IN THE SUBROUTINE STATE)

      ENTRY DPLEPH(ET2Z,NTARG,NCENT,RRD)

      ET2(1)=ET2Z(1)
      ET2(2)=ET2Z(2)

  11  IF(FIRST) CALL STATE(ZIPS,LIST,PVST,PNUT)
      FIRST=.FALSE.

      IF(NTARG .EQ. NCENT) RETURN

      DO I=1,12
        LIST(I)=0
      ENDDO

C     CHECK FOR NUTATION CALL

      IF(NTARG.NE.14) GO TO 97
        IF(IPT(35).GT.0) THEN
          LIST(11)=2
          CALL STATE(ET2,LIST,PVST,PNUT)
          DO I=1,4
            RRD(I)=PNUT(I)
          ENDDO
          RRD(5) = 0.d0
          RRD(6) = 0.d0
          RETURN
        ELSE
          DO I=1,4
            RRD(I)=0.d0
          ENDDO
          WRITE(6,297)
  297     FORMAT(' *****  NO NUTATIONS ON THE EPHEMERIS FILE  *****')
          STOP
        ENDIF

C     CHECK FOR LIBRATIONS

  97  CONTINUE
      DO I=1,6
        RRD(I)=0.d0
      ENDDO

      IF(NTARG.NE.15) GO TO 98
        IF(IPT(38).GT.0) THEN
          LIST(12)=2
          CALL STATE(ET2,LIST,PVST,PNUT)
          DO I=1,6
            RRD(I)=PVST(I,11)
          ENDDO
          RETURN
        ELSE
          WRITE(6,298)
  298     FORMAT(' *****  NO LIBRATIONS ON THE EPHEMERIS FILE  *****')
          STOP
        ENDIF

C       FORCE BARYCENTRIC OUTPUT BY 'STATE'

  98  BSAVE=BARY
      BARY=.TRUE.

C       SET UP PROPER ENTRIES IN 'LIST' ARRAY FOR STATE CALL

      DO I=1,2
        K=NTARG
        IF(I .EQ. 2) K=NCENT
        IF(K .LE. 10) LIST(K)=2
        IF(K .EQ. 10) LIST(3)=2
        IF(K .EQ. 3) LIST(10)=2
        IF(K .EQ. 13) LIST(3)=2
      ENDDO

C       MAKE CALL TO STATE

      CALL STATE(ET2,LIST,PVST,PNUT)

      DO I=1,10
        DO J = 1,6
          PV(J,I) = PVST(J,I)
        ENDDO
      ENDDO

      IF(NTARG .EQ. 11 .OR. NCENT .EQ. 11) THEN
      DO I=1,6
        PV(I,11)=PVSUN(I)
      ENDDO
      ENDIF

      IF(NTARG .EQ. 12 .OR. NCENT .EQ. 12) THEN
        DO I=1,6
          PV(I,12)=0.D0
        ENDDO
      ENDIF

      IF(NTARG .EQ. 13 .OR. NCENT .EQ. 13) THEN
        DO I=1,6
          PV(I,13) = PVST(I,3)
        ENDDO
      ENDIF

      IF(NTARG*NCENT .EQ. 30 .AND. NTARG+NCENT .EQ. 13) THEN
        DO I=1,6
          PV(I,3)=0.D0
        ENDDO
        GO TO 99
      ENDIF

      IF(LIST(3) .EQ. 2) THEN
        DO I=1,6
          PV(I,3)=PVST(I,3)-PVST(I,10)/(1.D0+EMRAT)
        ENDDO
      ENDIF

      IF(LIST(10) .EQ. 2) THEN
        DO I=1,6
          PV(I,10) = PV(I,3)+PVST(I,10)
        ENDDO
      ENDIF

  99  DO I=1,6
        RRD(I)=PV(I,NTARG)-PV(I,NCENT)
      ENDDO

      BARY=BSAVE

      RETURN
      END
C+++++++++++++++++++++++++++++++++
C
      SUBROUTINE INTERP(BUF,T,NCF,NCM,NA,IFL,PV)
C
C+++++++++++++++++++++++++++++++++
C
C     THIS SUBROUTINE DIFFERENTIATES AND INTERPOLATES A
C     SET OF CHEBYSHEV COEFFICIENTS TO GIVE POSITION AND VELOCITY
C
C     CALLING SEQUENCE PARAMETERS:
C
C       INPUT:
C
C         BUF   1ST LOCATION OF ARRAY OF D.P. CHEBYSHEV COEFFICIENTS OF POSITION
C
C           T   T(1) IS DP FRACTIONAL TIME IN INTERVAL COVERED BY
C               COEFFICIENTS AT WHICH INTERPOLATION IS WANTED
C               (0 .LE. T(1) .LE. 1).  T(2) IS DP LENGTH OF WHOLE
C               INTERVAL IN INPUT TIME UNITS.
C
C         NCF   # OF COEFFICIENTS PER COMPONENT
C
C         NCM   # OF COMPONENTS PER SET OF COEFFICIENTS
C
C          NA   # OF SETS OF COEFFICIENTS IN FULL ARRAY
C               (I.E., # OF SUB-INTERVALS IN FULL INTERVAL)
C
C          IFL  INTEGER FLAG: =1 FOR POSITIONS ONLY
C                             =2 FOR POS AND VEL
C
C
C       OUTPUT:
C
C         PV   INTERPOLATED QUANTITIES REQUESTED.  DIMENSION
C               EXPECTED IS PV(NCM,IFL), DP.
C
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      SAVE
C
      DOUBLE PRECISION BUF(NCF,NCM,*),T(2),PV(NCM,*),PC(18),VC(18)

C
      DATA NP/2/
      DATA NV/3/
      DATA TWOT/0.D0/
      DATA PC(1),PC(2)/1.D0,0.D0/
      DATA VC(2)/1.D0/
C
C       ENTRY POINT. GET CORRECT SUB-INTERVAL NUMBER FOR THIS SET
C       OF COEFFICIENTS AND THEN GET NORMALIZED CHEBYSHEV TIME
C       WITHIN THAT SUBINTERVAL.
C
      DNA=DBLE(NA)
      DT1=DINT(T(1))
      TEMP=DNA*T(1)
      L=IDINT(TEMP-DT1)+1

C         TC IS THE NORMALIZED CHEBYSHEV TIME (-1 .LE. TC .LE. 1)

      TC=2.D0*(DMOD(TEMP,1.D0)+DT1)-1.D0

C       CHECK TO SEE WHETHER CHEBYSHEV TIME HAS CHANGED,
C       AND COMPUTE NEW POLYNOMIAL VALUES IF IT HAS.
C       (THE ELEMENT PC(2) IS THE VALUE OF T1(TC) AND HENCE
C       CONTAINS THE VALUE OF TC ON THE PREVIOUS CALL.)

      IF(TC.NE.PC(2)) THEN
        NP=2
        NV=3
        PC(2)=TC
        TWOT=TC+TC
      ENDIF
C
C       BE SURE THAT AT LEAST 'NCF' POLYNOMIALS HAVE BEEN EVALUATED
C       AND ARE STORED IN THE ARRAY 'PC'.
C
      IF(NP.LT.NCF) THEN
        DO 1 I=NP+1,NCF
        PC(I)=TWOT*PC(I-1)-PC(I-2)
    1   CONTINUE
        NP=NCF
      ENDIF
C
C       INTERPOLATE TO GET POSITION FOR EACH COMPONENT
C
      DO 2 I=1,NCM
      PV(I,1)=0.D0
      DO 3 J=NCF,1,-1
      PV(I,1)=PV(I,1)+PC(J)*BUF(J,I,L)
    3 CONTINUE
    2 CONTINUE
      IF(IFL.LE.1) RETURN
C
C       IF VELOCITY INTERPOLATION IS WANTED, BE SURE ENOUGH
C       DERIVATIVE POLYNOMIALS HAVE BEEN GENERATED AND STORED.
C
      VFAC=(DNA+DNA)/T(2)
      VC(3)=TWOT+TWOT
      IF(NV.LT.NCF) THEN
        DO 4 I=NV+1,NCF
        VC(I)=TWOT*VC(I-1)+PC(I-1)+PC(I-1)-VC(I-2)
    4   CONTINUE
        NV=NCF
      ENDIF
C
C       INTERPOLATE TO GET VELOCITY FOR EACH COMPONENT
C
      DO 5 I=1,NCM
      PV(I,2)=0.D0
      DO 6 J=NCF,2,-1
      PV(I,2)=PV(I,2)+VC(J)*BUF(J,I,L)
    6 CONTINUE
      PV(I,2)=PV(I,2)*VFAC
    5 CONTINUE
C
      RETURN
C
      END

C+++++++++++++++++++++++++
C
      SUBROUTINE SPLIT(TT,FR)
C
C+++++++++++++++++++++++++
C
C     THIS SUBROUTINE BREAKS A D.P. NUMBER INTO A D.P. INTEGER
C     AND A D.P. FRACTIONAL PART.
C
C     CALLING SEQUENCE PARAMETERS:
C
C       TT = D.P. INPUT NUMBER
C
C       FR = D.P. 2-WORD OUTPUT ARRAY.
C            FR(1) CONTAINS INTEGER PART
C            FR(2) CONTAINS FRACTIONAL PART
C
C            FOR NEGATIVE INPUT NUMBERS, FR(1) CONTAINS THE NEXT
C            MORE NEGATIVE INTEGER; FR(2) CONTAINS A POSITIVE FRACTION.
C
C       CALLING SEQUENCE DECLARATIONS
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      DIMENSION FR(2)

C       MAIN ENTRY -- GET INTEGER AND FRACTIONAL PARTS

      FR(1)=DINT(TT)
      FR(2)=TT-FR(1)

      IF(TT.GE.0.D0 .OR. FR(2).EQ.0.D0) RETURN

C       MAKE ADJUSTMENTS FOR NEGATIVE INPUT NUMBER

      FR(1)=FR(1)-1.D0
      FR(2)=FR(2)+1.D0

      RETURN

      END

C++++++++++++++++++++++++++++++++
C
      SUBROUTINE STATE(ET2,LIST,PV,PNUT)
C
C++++++++++++++++++++++++++++++++
C
C THIS SUBROUTINE READS AND INTERPOLATES THE JPL PLANETARY EPHEMERIS FILE
C
C     CALLING SEQUENCE PARAMETERS:
C
C     INPUT:
C
C         ET2   DP 2-WORD JULIAN EPHEMERIS EPOCH AT WHICH INTERPOLATION
C               IS WANTED.  ANY COMBINATION OF ET2(1)+ET2(2) WHICH FALLS
C               WITHIN THE TIME SPAN ON THE FILE IS A PERMISSIBLE EPOCH.
C
C                A. FOR EASE IN PROGRAMMING, THE USER MAY PUT THE
C                   ENTIRE EPOCH IN ET2(1) AND SET ET2(2)=0.
C
C                B. FOR MAXIMUM INTERPOLATION ACCURACY, SET ET2(1) =
C                   THE MOST RECENT MIDNIGHT AT OR BEFORE INTERPOLATION
C                   EPOCH AND SET ET2(2) = FRACTIONAL PART OF A DAY
C                   ELAPSED BETWEEN ET2(1) AND EPOCH.
C
C                C. AS AN ALTERNATIVE, IT MAY PROVE CONVENIENT TO SET
C                   ET2(1) = SOME FIXED EPOCH, SUCH AS START OF INTEGRATION,
C                   AND ET2(2) = ELAPSED INTERVAL BETWEEN THEN AND EPOCH.
C
C        LIST   12-WORD INTEGER ARRAY SPECIFYING WHAT INTERPOLATION
C               IS WANTED FOR EACH OF THE BODIES ON THE FILE.
C
C                         LIST(I)=0, NO INTERPOLATION FOR BODY I
C                                =1, POSITION ONLY
C                                =2, POSITION AND VELOCITY
C
C               THE DESIGNATION OF THE ASTRONOMICAL BODIES BY I IS:
C
C                         I = 1: MERCURY
C                           = 2: VENUS
C                           = 3: EARTH-MOON BARYCENTER
C                           = 4: MARS
C                           = 5: JUPITER
C                           = 6: SATURN
C                           = 7: URANUS
C                           = 8: NEPTUNE
C                           = 9: PLUTO
C                           =10: GEOCENTRIC MOON
C                           =11: NUTATIONS IN LONGITUDE AND OBLIQUITY
C                           =12: LUNAR LIBRATIONS (IF ON FILE)
C
C     OUTPUT:
C
C          PV   DP 6 X 11 ARRAY THAT WILL CONTAIN REQUESTED INTERPOLATED
C               QUANTITIES (OTHER THAN NUTATION, STOERD IN PNUT).
C               THE BODY SPECIFIED BY LIST(I) WILL HAVE ITS
C               STATE IN THE ARRAY STARTING AT PV(1,I).
C               (ON ANY GIVEN CALL, ONLY THOSE WORDS IN 'PV' WHICH ARE
C                AFFECTED BY THE  FIRST 10 'LIST' ENTRIES, AND BY LIST(12)
C                IF LIBRATIONS ARE ON THE FILE, ARE SET.
C                THE REST OF THE 'PV' ARRAYIS UNTOUCHED.)
C               THE ORDER OF COMPONENTS STARTING IN PV(1,I) IS: X,Y,Z,DX,DY,DZ.
C
C               ALL OUTPUT VECTORS ARE REFERENCED TO THE EARTH MEAN
C               EQUATOR AND EQUINOX OF J2000 IF THE DE NUMBER IS 200 OR
C               GREATER; OF B1950 IF THE DE NUMBER IS LESS THAN 200.
C
C               THE MOON STATE IS ALWAYS GEOCENTRIC; THE OTHER NINE STATES
C               ARE EITHER HELIOCENTRIC OR SOLAR-SYSTEM BARYCENTRIC,
C               DEPENDING ON THE SETTING OF COMMON FLAGS (SEE BELOW).
C
C               LUNAR LIBRATIONS, IF ON FILE, ARE PUT INTO PV(K,11) IF
C               LIST(12) IS 1 OR 2.
C
C         NUT   DP 4-WORD ARRAY THAT WILL CONTAIN NUTATIONS AND RATES,
C               DEPENDING ON THE SETTING OF LIST(11).  THE ORDER OF
C               QUANTITIES IN NUT IS:
C
C                        D PSI  (NUTATION IN LONGITUDE)
C                        D EPSILON (NUTATION IN OBLIQUITY)
C                        D PSI DOT
C                        D EPSILON DOT
C
C           *   STATEMENT # FOR ERROR RETURN, IN CASE OF EPOCH OUT OF
C               RANGE OR I/O ERRORS.
C
C     COMMON AREA STCOMX:
C
C          KM   LOGICAL FLAG DEFINING PHYSICAL UNITS OF THE OUTPUT
C               STATES. KM = .TRUE., KM AND KM/SEC
C                          = .FALSE., AU AND AU/DAY
C               DEFAULT VALUE = .FALSE.  (KM DETERMINES TIME UNIT
C               FOR NUTATIONS AND LIBRATIONS.  ANGLE UNIT IS ALWAYS RADIANS.)
C
C        BARY   LOGICAL FLAG DEFINING OUTPUT CENTER.
C               ONLY THE 9 PLANETS ARE AFFECTED.
C                        BARY = .TRUE. =\ CENTER IS SOLAR-SYSTEM BARYCENTER
C                             = .FALSE. =\ CENTER IS SUN
C               DEFAULT VALUE = .FALSE.
C
C       PVSUN   DP 6-WORD ARRAY CONTAINING THE BARYCENTRIC POSITION AND
C               VELOCITY OF THE SUN.
C
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      SAVE

      INTEGER OLDMAX
      PARAMETER ( OLDMAX = 400)
      INTEGER NMAX
      PARAMETER ( NMAX = 1000)

      DIMENSION ET2(2),PV(6,11),PNUT(4),T(2),PJD(4),BUF(1500),
     . SS(3),CVAL(NMAX),PVSUN(6)

      INTEGER LIST(12),IPT(3,13)

      LOGICAL FIRST
      DATA FIRST/.TRUE./

      CHARACTER*6 TTL(14,3),CNAM(NMAX)
      CHARACTER*256 NAMFIL

      LOGICAL KM,BARY

      COMMON/EPHHDR/CVAL,SS,AU,EMRAT,NUMDE,NCON,IPT
      COMMON/CHRHDR/CNAM,TTL
      COMMON/STCOMX/KM,BARY,PVSUN

C
C       ENTRY POINT - 1ST TIME IN, GET POINTER DATA, ETC., FROM EPH FILE
C
      IF(FIRST) THEN
        FIRST=.FALSE.

C ************************************************************************
C ************************************************************************

C THE USER MUST SELECT ONE OF THE FOLLOWING BY DELETING THE 'C' IN COLUMN 1

C ************************************************************************

C        CALL FSIZER1(NRECL,KSIZE,NRFILE,NAMFIL)
C        CALL FSIZER2(NRECL,KSIZE,NRFILE,NAMFIL)
        CALL FSIZER3(NRECL,KSIZE,NRFILE,NAMFIL)

      IF(NRECL .EQ. 0) WRITE(*,*)'  ***** FSIZER IS NOT WORKING *****'

C ************************************************************************
C ************************************************************************

      IRECSZ=NRECL*KSIZE
      NCOEFFS=KSIZE/2

        OPEN(NRFILE,
     *       FILE=NAMFIL,
     *       ACCESS='DIRECT',
     *       FORM='UNFORMATTED',
     *       RECL=IRECSZ,
     *       STATUS='OLD')

      READ(NRFILE,REC=1)TTL,(CNAM(K),K=1,OLDMAX),SS,NCON,AU,EMRAT,
     & ((IPT(I,J),I=1,3),J=1,12),NUMDE,(IPT(I,13),I=1,3)
     & ,(CNAM(L),L=OLDMAX+1,NCON)

      IF(NCON .LE. OLDMAX)THEN
        READ(NRFILE,REC=2)(CVAL(I),I=1,OLDMAX)
      ELSE
        READ(NRFILE,REC=2)(CVAL(I),I=1,NCON)
      ENDIF

      NRL=0

      ENDIF

C       ********** MAIN ENTRY POINT **********

      IF(ET2(1) .EQ. 0.D0) RETURN

      S=ET2(1)-.5D0
      CALL SPLIT(S,PJD(1))
      CALL SPLIT(ET2(2),PJD(3))
      PJD(1)=PJD(1)+PJD(3)+.5D0
      PJD(2)=PJD(2)+PJD(4)
      CALL SPLIT(PJD(2),PJD(3))
      PJD(1)=PJD(1)+PJD(3)

C       ERROR RETURN FOR EPOCH OUT OF RANGE

      IF(PJD(1)+PJD(4).LT.SS(1) .OR. PJD(1)+PJD(4).GT.SS(2)) GO TO 98

C       CALCULATE RECORD # AND RELATIVE TIME IN INTERVAL

      NR=IDINT((PJD(1)-SS(1))/SS(3))+3
      IF(PJD(1).EQ.SS(2)) NR=NR-1

        tmp1 = DBLE(NR-3)*SS(3) + SS(1)
        tmp2 = PJD(1) - tmp1
        T(1) = (tmp2 + PJD(4))/SS(3)

C       READ CORRECT RECORD IF NOT IN CORE

      IF(NR.NE.NRL) THEN
        NRL=NR
        READ(NRFILE,REC=NR,ERR=99)(BUF(K),K=1,NCOEFFS)
      ENDIF

      IF(KM) THEN
      T(2)=SS(3)*86400.D0
      AUFAC=1.D0
      ELSE
      T(2)=SS(3)
      AUFAC=1.D0/AU
      ENDIF

C   INTERPOLATE SSBARY SUN

      CALL INTERP(BUF(IPT(1,11)),T,IPT(2,11),3,IPT(3,11),2,PVSUN)

      DO I=1,6
      PVSUN(I)=PVSUN(I)*AUFAC
      ENDDO

C   CHECK AND INTERPOLATE WHICHEVER BODIES ARE REQUESTED

      DO 4 I=1,10
      IF(LIST(I).EQ.0) GO TO 4

      CALL INTERP(BUF(IPT(1,I)),T,IPT(2,I),3,IPT(3,I),
     & LIST(I),PV(1,I))

      DO J=1,6
       IF(I.LE.9 .AND. .NOT.BARY) THEN
       PV(J,I)=PV(J,I)*AUFAC-PVSUN(J)
       ELSE
       PV(J,I)=PV(J,I)*AUFAC
       ENDIF
      ENDDO

   4  CONTINUE

C       DO NUTATIONS IF REQUESTED (AND IF ON FILE)

      IF(LIST(11).GT.0 .AND. IPT(2,12).GT.0)
     * CALL INTERP(BUF(IPT(1,12)),T,IPT(2,12),2,IPT(3,12),
     * LIST(11),PNUT)

C       GET LIBRATIONS IF REQUESTED (AND IF ON FILE)

      IF(LIST(12).GT.0 .AND. IPT(2,13).GT.0)
     * CALL INTERP(BUF(IPT(1,13)),T,IPT(2,13),3,IPT(3,13),
     * LIST(12),PV(1,11))

      RETURN

  98  WRITE(*,198)ET2(1)+ET2(2),SS(1),SS(2)
 198  FORMAT(' ***  Requested JED,',f12.2,
     * ' not within ephemeris limits,',2f12.2,'  ***')

      STOP

   99 WRITE(*,'(2F12.2,A80)')ET2,'ERROR RETURN IN STATE'

      STOP

      END
C+++++++++++++++++++++++++++++
C
      SUBROUTINE CONST(NAM,VAL,SSS,N)
C
C+++++++++++++++++++++++++++++
C
C     THIS ENTRY OBTAINS THE CONSTANTS FROM THE EPHEMERIS FILE
C
C     CALLING SEQEUNCE PARAMETERS (ALL OUTPUT):
C
C       NAM = CHARACTER*6 ARRAY OF CONSTANT NAMES
C
C       VAL = D.P. ARRAY OF VALUES OF CONSTANTS
C
C       SSS = D.P. JD START, JD STOP, STEP OF EPHEMERIS
C
C         N = INTEGER NUMBER OF ENTRIES IN 'NAM' AND 'VAL' ARRAYS
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

      SAVE

      INTEGER NMAX
      PARAMETER (NMAX = 1000)

      CHARACTER*6 NAM(*),TTL(14,3),CNAM(NMAX)

      DOUBLE PRECISION VAL(*),SSS(3),SS(3),CVAL(NMAX),ZIPS(2)
      DOUBLE PRECISION PVST(6,11),PNUT(4)
      DATA ZIPS/2*0.d0/

      INTEGER IPT(3,13),DENUM,LIST(12)
      logical first
      data first/.true./

      COMMON/EPHHDR/CVAL,SS,AU,EMRAT,DENUM,NCON,IPT
      COMMON/CHRHDR/CNAM,TTL

C  CALL STATE TO INITIALIZE THE EPHEMERIS AND READ IN THE CONSTANTS

      IF(FIRST) CALL STATE(ZIPS,LIST,PVST,PNUT)
      first=.false.

      N=NCON

      DO I=1,3
        SSS(I)=SS(I)
      ENDDO

      DO I=1,N
        NAM(I)=CNAM(I)
        VAL(I)=CVAL(I)
      ENDDO

      RETURN

      END


code:
program jt49sim

! Generate simulated data for testing JT4 and JT9

  use wavhdr
  use packjt
  use jt4
  parameter (NMAX=60*12000)              ! = 648,000
  parameter (NFFT=10*65536,NH=NFFT/2)
  type(hdr) h                            !Header for .wav file
  integer*2 iwave(NMAX)                  !Generated waveform
  integer*4 itone(206)                     !Channel symbols (values 0-8)
  real*4 xnoise(NMAX)                    !Generated random noise
  real*4 dat(NMAX)                       !Generated real data
  complex cdat(NMAX)                     !Generated complex waveform
  complex cspread(0:NFFT-1)              !Complex amplitude for Rayleigh fading
  complex z
  real*8 f0,dt,twopi,phi,dphi,baud,fsample,freq,dnsps
  character message*22,fname*11,csubmode*2,arg*12
  character msgsent*22

  nargs=iargc()
  if(nargs.ne. 7) then
     print *, 'Usage:   jt49sim        "msg"     nA-nE Nsigs fDop DT Nfiles SNR'
     print *, 'Example  jt49sim "K1ABC W9XYZ EN37" 4G   10   0.2 0.0   1     0'
     print *, 'Example  jt49sim "K1ABC W9XYZ EN37" 9A    1   0.0 0.0   1    -20'
     print *, 'Use msg=@nnnn to generate a tone at nnnn Hz:'
     print *, 'Example  jt49sim "@1500"            9A    1  10.0 0.0   1    -20'
     go to 999
  endif
  call getarg(1,message)
  call fmtmsg(message, iz)
  call getarg(2,csubmode)
  imode=ichar(csubmode(1:1)) - ichar('0')
  nsubmode=ichar(csubmode(2:2)) - ichar('A')
  if(imode.ne.4 .and. imode.ne.9) go to 999
  if(nsubmode.lt.0 .or. nsubmode.gt.7) go to 999
  call getarg(3,arg)
  read(arg,*) nsigs
  call getarg(4,arg)
  read(arg,*) fspread
  call getarg(5,arg)
  read(arg,*) xdt
  call getarg(6,arg)
  read(arg,*) nfiles
  call getarg(7,arg)
  read(arg,*) snrdb

  rms=100.
  fsample=12000.d0                   !Sample rate (Hz)
  dt=1.d0/fsample                    !Sample interval (s)
  twopi=8.d0*atan(1.d0)
  npts=60*12000                      !Total samples in .wav file
  h=default_header(12000,npts)
  dfsig=2000.0/nsigs                 !Freq spacing between sigs in file (Hz)
  ichk=0

  if(imode.eq.4) then
     nsym=206                           !Number of channel symbols (JT4)
     dnsps=12000.d0/4.375d0
     baud=12000.d0/dnsps                !Keying rate = 1.7361111111
  else if(imode.eq.9) then
     nsym=85                            !Number of channel symbols (JT9)
     dnsps=6912.d0                      !Samples per symbol
     baud=12000.d0/dnsps                !Keying rate = 1.736...
  endif

  write(*,1000)
1000 format('File  Sig    Freq  Mode   S/N   DT   Dop    Message'/60('-'))

  do ifile=1,nfiles                  !Loop over requested number of files
     write(fname,1002) ifile         !Output filename
1002 format('000000_',i4.4)
     open(10,file=fname//'.wav',access='stream',status='unknown')
     xnoise=0.
     cdat=0.
     if(snrdb.lt.90) then
        do i=1,npts
           xnoise(i)=gran()          !Generate gaussian noise
        enddo
     endif

     do isig=1,nsigs                 !Generate requested number of sigs
        if(mod(nsigs,2).eq.0) f0=1500.0 + dfsig*(isig-0.5-nsigs/2)
        if(mod(nsigs,2).eq.1) f0=1500.0 + dfsig*(isig-(nsigs+1)/2)
        if(nsigs.eq.1) f0=1000.0
        xsnr=snrdb
        if(snrdb.eq.0.0) xsnr=-20 - isig

        if(imode.eq.4) call gen4(message,ichk,msgsent,itone,itype)
        if(imode.eq.9) call gen9(message,ichk,msgsent,itone,itype)

        bandwidth_ratio=2500.0/6000.0
        sig=sqrt(2*bandwidth_ratio)*10.0**(0.05*xsnr)
        if(xsnr.gt.90.0) sig=1.0
        write(*,1020) ifile,isig,f0,csubmode,xsnr,xdt,fspread,message
1020    format(i4,i4,f10.3,2x,a2,2x,f5.1,f6.2,f6.1,1x,a22)

        phi=0.d0
        dphi=0.d0
        k=(xdt+1.0)*12000                   !Start audio at t = xdt + 1.0 s
        isym0=-99

        do i=1,npts                         !Add this signal into cdat()
           isym=i/dnsps + 1
           if(isym.gt.nsym) exit
           if(isym.ne.isym0) then
              if(message(1:1).eq.'@') then
                 read(message(2:),*) freq
              else
                 if(imode.eq.4) freq=f0 + itone(isym)*baud*nch(1+nsubmode) !JT4
                 if(imode.eq.9) freq=f0 + itone(isym)*baud*(2**nsubmode)   !JT9
              endif
              dphi=twopi*freq*dt
              isym0=isym
           endif
           phi=phi + dphi
           if(phi.gt.twopi) phi=phi-twopi
           xphi=phi
           z=cmplx(cos(xphi),sin(xphi))
           k=k+1
           if(k.ge.1) cdat(k)=cdat(k) + sig*z
        enddo
     enddo

     if(fspread.ne.0) then                  !Apply specified Doppler spread
        df=12000.0/nfft
        twopi=8*atan(1.0)
        cspread(0)=1.0
        cspread(NH)=0.
        b=6.0                                     !Lorenzian 3/28 onward
        do i=1,NH
           f=i*df
           x=b*f/fspread
           z=0.
           a=0.
           if(x.lt.3.0) then                          !Cutoff beyond x=3
              a=sqrt(1.111/(1.0+x*x)-0.1)             !Lorentzian
              call random_number(r1)
              phi1=twopi*r1
              z=a*cmplx(cos(phi1),sin(phi1))
           endif
           cspread(i)=z
           z=0.
           if(x.lt.50.0) then
              call random_number(r2)
              phi2=twopi*r2
              z=a*cmplx(cos(phi2),sin(phi2))
           endif
           cspread(NFFT-i)=z
        enddo

        do i=0,NFFT-1
           f=i*df
           if(i.gt.NH) f=(i-nfft)*df
           s=real(cspread(i))**2 + aimag(cspread(i))**2
!          write(13,3000) i,f,s,cspread(i)
!3000      format(i5,f10.3,3f12.6)
        enddo
!        s=real(cspread(0))**2 + aimag(cspread(0))**2
!        write(13,3000) 1024,0.0,s,cspread(0)

        call four2a(cspread,NFFT,1,1,1)             !Transform to time domain

        sum=0.
        do i=0,NFFT-1
           p=real(cspread(i))**2 + aimag(cspread(i))**2
           sum=sum+p
        enddo
        avep=sum/NFFT
        fac=sqrt(1.0/avep)
        cspread=fac*cspread                   !Normalize to constant avg power
        cdat(1:NFFT)=cspread*cdat(1:NFFT)     !Apply Rayleigh fading

!        do i=0,NFFT-1
!           p=real(cspread(i))**2 + aimag(cspread(i))**2
!           write(14,3010) i,p,cspread(i)
!3010       format(i8,3f12.6)
!        enddo

     endif

     dat=aimag(cdat) + xnoise                 !Add the generated noise
     fac=32767.0/nsigs
     if(snrdb.ge.90.0) iwave(1:npts)=nint(fac*dat(1:npts))
     if(snrdb.lt.90.0) iwave(1:npts)=nint(rms*dat(1:npts))
     write(10) h,iwave(1:npts)                !Save the .wav file
     close(10)
  enddo

999 end program jt49sim

license: mit


code:
module jt4_decode
  type :: jt4_decoder
     procedure(jt4_decode_callback), pointer :: decode_callback => null ()
     procedure(jt4_average_callback), pointer :: average_callback => null ()
   contains
     procedure :: decode
     procedure, private :: wsjt4, avg4
  end type jt4_decoder

! Callback function to be called with each decode
  abstract interface
     subroutine jt4_decode_callback (this, snr, dt, freq, have_sync,     &
          sync, is_deep, decoded, qual, ich, is_average, ave)
       import jt4_decoder
       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) :: decoded
       real, intent(in) :: qual
       integer, intent(in) :: ich
       logical, intent(in) :: is_average
       integer, intent(in) :: ave
     end subroutine jt4_decode_callback
  end interface

! Callback function to be called with each average result
  abstract interface
     subroutine jt4_average_callback (this, used, utc, sync, dt, freq, flip)
       import jt4_decoder
       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
     end subroutine jt4_average_callback
  end interface

contains

  subroutine decode(this,decode_callback,dd,jz,nutc,nfqso,ntol0,emedelay,     &
       dttol,nagain,ndepth,nclearave,minsync,minw,nsubmode,mycall,hiscall,    &
       hisgrid,nlist0,listutc0,average_callback)

    use jt4
    use timer_module, only: timer

    class(jt4_decoder), intent(inout) :: this
    procedure(jt4_decode_callback) :: decode_callback
    integer, intent(in) :: jz,nutc,nfqso,ntol0,ndepth,minsync,minw,nsubmode,  &
         nlist0,listutc0(10)
    real, intent(in) :: dd(jz),emedelay,dttol
    logical, intent(in) :: nagain, nclearave
    character(len=12), intent(in) :: mycall,hiscall
    character(len=6), intent(in) :: hisgrid
    procedure(jt4_average_callback), optional :: average_callback

    real*4 dat(30*11025)
    character*6 cfile6

    this%decode_callback => decode_callback
    if (present (average_callback)) then
       this%average_callback => average_callback
    end if
    mode4=nch(nsubmode+1)
    ntol=ntol0
    neme=0
    lumsg=6                         !### temp ? ###
    ndiag=1
    nlist=nlist0
    listutc=listutc0

    ! Lowpass filter and decimate by 2
    call timer('lpf1    ',0)
    call lpf1(dd,jz,dat,jz2)
    call timer('lpf1    ',1)

    write(cfile6(1:4),1000) nutc
1000 format(i4.4)
    cfile6(5:6)='  '

    call timer('wsjt4   ',0)
    call this%wsjt4(dat,jz2,nutc,NClearAve,minsync,ntol,emedelay,dttol,mode4, &
         minw,mycall,hiscall,hisgrid,nfqso,NAgain,ndepth,neme)
    call timer('wsjt4   ',1)

    return
  end subroutine decode

  subroutine wsjt4(this,dat,npts,nutc,NClearAve,minsync,ntol,emedelay,dttol,  &
       mode4,minw,mycall,hiscall,hisgrid,nfqso,NAgain,ndepth,neme)

! Orchestrates the process of decoding JT4 messages.  Note that JT4
! always operates as if in "Single Decode" mode; it looks for only one
! decodable signal in the FTol range.

    use jt4
    use timer_module, only: timer

    class(jt4_decoder), intent(inout) :: this
    integer, intent(in) :: npts,nutc,minsync,ntol,mode4,minw,       &
         nfqso,ndepth,neme
    logical, intent(in) :: NAgain,NClearAve
    character(len=12), intent(in) :: mycall,hiscall
    character(len=6), intent(in) :: hisgrid
    real, intent(in) :: dat(npts) !Raw data

    real ccfblue(-5:540)                             !CCF in time
    real ccfred(-224:224)                            !CCF in frequency
    real ps0(450)

!    real z(458,65)
    logical first,prtavg
    character decoded*22,special*5
    character*22 avemsg,deepmsg,deepave,blank,deepmsg0,deepave1
    character csync*1
    data first/.true./,nutc0/-999/,nfreq0/-999999/
    save

    if(first) then
       nsave=0
       first=.false.
       blank='                      '
       ccfblue=0.
       ccfred=0.
! Silence compiler warnings
       if(dttol.eq.-99.0 .and. emedelay.eq.-99.0 .and. nagain) stop
    endif

    zz=0.
!    syncmin=3.0 + minsync
    syncmin=1.0+minsync
    naggressive=0
    if(ndepth.ge.2) naggressive=1
    nq1=3
    nq2=6
    if(naggressive.eq.1) nq1=1
    if(NClearAve) then
       nsave=0
       iutc=-1
       nfsave=0.
       listutc=0
       ppsave=0.
       rsymbol=0.
       dtsave=0.
       syncsave=0.
       nfanoave=0
       ndeepave=0
    endif

! Attempt to synchronize: look for sync pattern, get DF and DT.
    call timer('sync4   ',0)
    mousedf=nint(nfqso + 1.5*4.375*mode4 - 1270.46)
    call sync4(dat,npts,ntol,1,MouseDF,4,mode4,minw+1,dtx,dfx,    &
         snrx,snrsync,ccfblue,ccfred,flip,width,ps0)
    sync=snrsync
    dtxz=dtx-0.8
    nfreqz=dfx + 1270.46 - 1.5*4.375*mode4
    call timer('sync4   ',1)

    snrx=db(sync) - 26.
    nsnr=nint(snrx)
    if(sync.lt.syncmin) then
       if (associated (this%decode_callback)) then
          call this%decode_callback(nsnr,dtxz,nfreqz,.false.,csync,      &
               .false.,decoded,0.,ich,.false.,0)
       end if
       go to 990
    endif

! We have achieved sync
    decoded=blank
    deepmsg=blank
    special='     '
    nsync=sync
    nsnrlim=-33
    csync='*'
    if(flip.lt.0.0) csync='#'
    qbest=0.
    qabest=0.
    prtavg=.false.

    do idt=-2,2
       dtx=dtxz + 0.03*idt
       nfreq=nfreqz + 2*idf

! Attempt a single-sequence decode, including deep4 if Fano fails.
       call timer('decode4 ',0)
       call decode4(dat,npts,dtx,nfreq,flip,mode4,ndepth,neme,minw,           &
            mycall,hiscall,hisgrid,decoded,nfano,deepmsg,qual,ich)
       call timer('decode4 ',1)

       if(nfano.gt.0) then
! Fano succeeded: report the message and return              !Fano OK
          if (associated (this%decode_callback)) then
             call this%decode_callback(nsnr,dtx,nfreq,.true.,csync,      &
                  .false.,decoded,99.,ich,.false.,0)
          end if
          nsave=0
          go to 990

       else                                                  !Fano failed
          if(qual.gt.qbest) then
             dtx0=dtx
             nfreq0=nfreq
             deepmsg0=deepmsg
             ich0=ich
             qbest=qual
          endif
       endif

       if(idt.ne.0) cycle
! Single-sequence Fano decode failed, so try for an average Fano decode:
       qave=0.
! If we're doing averaging, call avg4
       if(iand(ndepth,16).eq.16 .and. (.not.prtavg)) then
          if(nutc.ne.nutc0 .or. abs(nfreq-nfreq0).gt.ntol) then
! This is a new minute or a new frequency, so call avg4.
             nutc0=nutc                                   !Try decoding average
             nfreq0=nfreq
             nsave=nsave+1
             nsave=mod(nsave-1,64)+1
             call timer('avg4    ',0)
             call this%avg4(nutc,sync,dtx,flip,nfreq,mode4,ntol,ndepth,neme,  &
                  mycall,hiscall,hisgrid,nfanoave,avemsg,qave,deepave,ich,    &
                  ndeepave)
             call timer('avg4    ',1)
          endif

          if(nfanoave.gt.0) then
! Fano succeeded: report the message                       AVG FANO OK
             if (associated (this%decode_callback)) then
                call this%decode_callback(nsnr,dtx,nfreq,.true.,csync,   &
                     .false.,avemsg,99.,ich,.true.,nfanoave)
             end if
             prtavg=.true.
             cycle
          else
             if(qave.gt.qabest) then
                dtx1=dtx
                nfreq1=nfreq
                deepave1=deepave
                ich1=ich
                qabest=qave
             endif
          endif
       endif
    enddo

    dtx=dtx0
    nfreq=nfreq0
    deepmsg=deepmsg0
    ich=ich0
    qual=qbest

    if (associated (this%decode_callback)) then
       if(int(qual).ge.nq1) then
          call this%decode_callback(nsnr,dtx,nfreqz,.true.,csync,.true., &
               deepmsg,qual,ich,.false.,0)
       else
          call this%decode_callback(nsnr,dtxz,nfreqz,.true.,csync,       &
               .false.,blank,0.,ich,.false.,0)
       endif
    end if

    dtx=dtx1
    nfreq=nfreq1
    deepave=deepave1
    ich=ich1
    qave=qabest

    if (associated (this%decode_callback) .and. ndeepave.ge.2) then
       if(int(qave).ge.nq1) then
          call this%decode_callback(nsnr,dtx,nfreq,.true.,csync,.true.,  &
               deepave,qave,ich,.true.,ndeepave)
       endif
    end if

990 return
  end subroutine wsjt4

  subroutine avg4(this,nutc,snrsync,dtxx,flip,nfreq,mode4,ntol,ndepth,neme,   &
       mycall,hiscall,hisgrid,nfanoave,avemsg,qave,deepave,ichbest,ndeepave)

! Decodes averaged JT4 data

    use jt4
    class(jt4_decoder), intent(inout) :: this

    character*22 avemsg,deepave,deepbest
    character mycall*12,hiscall*12,hisgrid*6
    character*1 csync,cused(64)
    real sym(207,7)
    integer iused(64)
    logical first
    data first/.true./
    save

    if(first) then
       iutc=-1
       nfsave=0
       dtdiff=0.2
       first=.false.
    endif

    do i=1,64
       if(nutc.eq.iutc(i) .and. abs(nfreq-nfsave(i)).le.ntol) go to 10
    enddo

! Save data for message averaging
    iutc(nsave)=nutc
    syncsave(nsave)=snrsync
    dtsave(nsave)=dtxx
    nfsave(nsave)=nfreq
    flipsave(nsave)=flip
    ppsave(1:207,1:7,nsave)=rsymbol(1:207,1:7)

10  sym=0.
    syncsum=0.
    dtsum=0.
    nfsum=0
    nsum=0

    do i=1,64
       cused(i)='.'
       if(iutc(i).lt.0) cycle
       if(mod(iutc(i),2).ne.mod(nutc,2)) cycle  !Use only same sequence
       if(abs(dtxx-dtsave(i)).gt.dtdiff) cycle  !DT must match
       if(abs(nfreq-nfsave(i)).gt.ntol) cycle   !Freq must match
       if(flip.ne.flipsave(i)) cycle            !Sync (*/#) must match
       sym(1:207,1:7)=sym(1:207,1:7) +  ppsave(1:207,1:7,i)
       syncsum=syncsum + syncsave(i)
       dtsum=dtsum + dtsave(i)
       nfsum=nfsum + nfsave(i)
       cused(i)='$'
       nsum=nsum+1
       iused(nsum)=i
    enddo
    if(nsum.lt.64) iused(nsum+1)=0

    syncave=0.
    dtave=0.
    fave=0.
    if(nsum.gt.0) then
       sym=sym/nsum
       syncave=syncsum/nsum
       dtave=dtsum/nsum
       fave=float(nfsum)/nsum
    endif

    do i=1,nsave
       csync='*'
       if(flipsave(i).lt.0.0) csync='#'
       if (associated (this%average_callback)) then
          call this%average_callback(cused(i) .eq. '$',iutc(i),               &
               syncsave(i),dtsave(i),nfsave(i),flipsave(i).lt.0.)
       end if
    enddo

    sqt=0.
    sqf=0.
    do j=1,64
       i=iused(j)
       if(i.eq.0) exit
       csync='*'
       if(flipsave(i).lt.0.0) csync='#'
       sqt=sqt + (dtsave(i)-dtave)**2
       sqf=sqf + (nfsave(i)-fave)**2
    enddo
    rmst=0.
    rmsf=0.
    if(nsum.ge.2) then
       rmst=sqrt(sqt/(nsum-1))
       rmsf=sqrt(sqf/(nsum-1))
    endif
    kbest=ich1
    do k=ich1,ich2
       call extract4(sym(1,k),ncount,avemsg)     !Do the Fano decode
       nfanoave=0
       if(ncount.ge.0) then
          ichbest=k
          nfanoave=nsum
          go to 900
       endif
       if(nch(k).ge.mode4) exit
    enddo

    deepave='                      '
    qave=0.

! Possibly should pass nadd=nused, also ?
    if(iand(ndepth,32).eq.32) then
       flipx=1.0                     !Normal flip not relevant for ave msg
       qbest=0.
       do k=ich1,ich2
          call deep4(sym(2,k),neme,flipx,mycall,hiscall,hisgrid,deepave,qave)
          if(qave.gt.qbest) then
             qbest=qave
             deepbest=deepave
             kbest=k
             ndeepave=nsum
          endif
          if(nch(k).ge.mode4) exit
       enddo

       deepave=deepbest
       qave=qbest
       ichbest=kbest
    endif

900 return
  end subroutine avg4
end module jt4_decode


code: ! The contents of this file have been migrated to lib/jt4_decode.f90


code:
program jt4code

! Provides examples of message packing, bit and symbol ordering,
! convolutional encoding, and other necessary details of the JT4
! protocol.

  use jt4
  use packjt
  character*22 msg,decoded,bad*1,msgtype*13
  integer i4tone(206)
  include 'testmsg.f90'

  nargs=iargc()
  if(nargs.ne.1) then
     print*,'Usage: jt4code "message"'
     print*,'       jt4code -t'
     go to 999
  endif

  call getarg(1,msg)
  nmsg=1
  if(msg(1:2).eq."-t") nmsg=NTEST

  write(*,1010)
1010 format("     Message                 Decoded                Err? Type"/   &
            74("-"))
  do imsg=1,nmsg
     if(nmsg.gt.1) msg=testmsg(imsg)
     call fmtmsg(msg,iz)                !To upper case, collapse multiple blanks
     ichk=0
     call gen4(msg,ichk,decoded,i4tone,itype)

     msgtype=""
     if(itype.eq.1) msgtype="Std Msg"
     if(itype.eq.2) msgtype="Type 1 prefix"
     if(itype.eq.3) msgtype="Type 1 suffix"
     if(itype.eq.4) msgtype="Type 2 prefix"
     if(itype.eq.5) msgtype="Type 2 suffix"
     if(itype.eq.6) msgtype="Free text"

     bad=" "
     if(decoded.ne.msg) bad="*"
     write(*,1020) imsg,msg,decoded,bad,itype,msgtype
1020 format(i2,'.',2x,a22,2x,a22,3x,a1,i3,": ",a13)
  enddo

  if(nmsg.eq.1) write(*,1030) i4tone
1030 format(/'Channel symbols'/(30i2))

999 end program jt4code


code:
program jt4sim

! Generate simulated jt4 data for testing the decoder using new fading models
! Based on qra64sim with appropriate changes for JT4
! Includes option to generate 1000Hz single tone with "ST" as the message

  use wavhdr
  use packjt
  use jt4
  parameter (NMAX=60*12000)              ! = 648,000
  parameter (NFFT=10*65536,NH=NFFT/2)
  type(hdr) h                            !Header for .wav file
  integer*2 iwave(NMAX)                  !Generated waveform
  integer*4 itone(206)                     !Channel symbols (values 0-8)
  real*4 xnoise(NMAX)                    !Generated random noise
  real*4 dat(NMAX)                       !Generated real data
  complex cdat(NMAX)                     !Generated complex waveform
  complex cspread(0:NFFT-1)              !Complex amplitude for Rayleigh fading
  complex z
  real*8 f0,dt,twopi,phi,dphi,baud,fsample,freq
  character message*22,fname*11,csubmode*1,arg*12
  character msgsent*22

  nargs=iargc()
  if(nargs.ne. 7) then
     print *, 'Usage:   jt4sim        "msg"     A-E Nsigs fDop DT Nfiles SNR'
     print *, 'Example  jt4sim "K1ABC W9XYZ EN37" A   10   0.2 0.0   1     0'
     print *, 'Example  jt4sim "ST"               A   10   0.2 0.0   1     0'
     go to 999
  endif
  call getarg(1,message)
  call fmtmsg(message, iz)
  call getarg(2,csubmode)
  mode4=ichar(csubmode) - ichar('A') + 1
  call getarg(3,arg)
  read(arg,*) nsigs
  call getarg(4,arg)
  read(arg,*) fspread
  call getarg(5,arg)
  read(arg,*) xdt
  call getarg(6,arg)
  read(arg,*) nfiles
  call getarg(7,arg)
  read(arg,*) snrdb

  if(mode4.ge.8) nsigs=1 ! temporary - will need sorting out
  rms=100.
  fsample=12000.d0                   !Sample rate (Hz)
  dt=1.d0/fsample                    !Sample interval (s)
  twopi=8.d0*atan(1.d0)
  npts=60*12000                      !Total samples in .wav file
  nsps=12000.d0/4.375d0
  baud=12000.d0/nsps                 !Keying rate = 1.7361111111
  nsym=206                           !Number of channel symbols
  h=default_header(12000,npts)
  dfsig=2000.0/nsigs                 !Freq spacing between sigs in file (Hz)
  ichk=0

  write(*,1000)
1000 format('File  Sig    Freq  A-E   S/N   DT   Dop    Message'/60('-'))

  do ifile=1,nfiles                  !Loop over requested number of files
     write(fname,1002) ifile         !Output filename
1002 format('000000_',i4.4)
     open(10,file=fname//'.wav',access='stream',status='unknown')
     xnoise=0.
     cdat=0.
     if(snrdb.lt.90) then
        do i=1,npts
           xnoise(i)=gran()          !Generate gaussian noise
        enddo
     endif

     do isig=1,nsigs                 !Generate requested number of sigs
        if(mod(nsigs,2).eq.0) f0=1500.0 + dfsig*(isig-0.5-nsigs/2)
        if(mod(nsigs,2).eq.1) f0=1500.0 + dfsig*(isig-(nsigs+1)/2)
        if(nsigs.eq.1) f0=1000.0
        xsnr=snrdb
        if(snrdb.eq.0.0) xsnr=-20 - isig

        call gen4(message,ichk,msgsent,itone,itype)

        bandwidth_ratio=2500.0/6000.0
        sig=sqrt(2*bandwidth_ratio)*10.0**(0.05*xsnr)
        if(xsnr.gt.90.0) sig=1.0
        write(*,1020) ifile,isig,f0,csubmode,xsnr,xdt,fspread,message
1020    format(i4,i4,f10.3,2x,a1,2x,f5.1,f6.2,f6.1,1x,a22)

        phi=0.d0
        dphi=0.d0
        k=(xdt+1.0)*12000                   !Start audio at t = xdt + 1.0 s
        isym0=-99

        do i=1,npts                         !Add this signal into cdat()
           isym=i/nsps + 1
           if(isym.gt.nsym) exit
           if(isym.ne.isym0) then
           if(message(1:2).eq.'ST') then    !Finds ST in message string for single tone
              freq=1000                     !1000Hz for single tone
           else
              freq=f0 + itone(isym)*baud*nch(mode4)
           endif
              dphi=twopi*freq*dt
              isym0=isym

           endif
           phi=phi + dphi
           if(phi.gt.twopi) phi=phi-twopi
           xphi=phi
           z=cmplx(cos(xphi),sin(xphi))
           k=k+1
           if(k.ge.1) cdat(k)=cdat(k) + sig*z
        enddo
     enddo

     if(fspread.ne.0) then                  !Apply specified Doppler spread
        df=12000.0/nfft
        twopi=8*atan(1.0)
        cspread(0)=1.0
        cspread(NH)=0.
        b=6.0                                     !Lorenzian 3/28 onward
        do i=1,NH
           f=i*df
           x=b*f/fspread
           z=0.
           a=0.
           if(x.lt.3.0) then                          !Cutoff beyond x=3
              a=sqrt(1.111/(1.0+x*x)-0.1)             !Lorentzian
              call random_number(r1)
              phi1=twopi*r1
              z=a*cmplx(cos(phi1),sin(phi1))
           endif
           cspread(i)=z
           z=0.
           if(x.lt.50.0) then
              call random_number(r2)
              phi2=twopi*r2
              z=a*cmplx(cos(phi2),sin(phi2))
           endif
           cspread(NFFT-i)=z
        enddo

        do i=0,NFFT-1
           f=i*df
           if(i.gt.NH) f=(i-nfft)*df
           s=real(cspread(i))**2 + aimag(cspread(i))**2
!          write(13,3000) i,f,s,cspread(i)
!3000      format(i5,f10.3,3f12.6)
        enddo
!        s=real(cspread(0))**2 + aimag(cspread(0))**2
!        write(13,3000) 1024,0.0,s,cspread(0)

        call four2a(cspread,NFFT,1,1,1)             !Transform to time domain

        sum=0.
        do i=0,NFFT-1
           p=real(cspread(i))**2 + aimag(cspread(i))**2
           sum=sum+p
        enddo
        avep=sum/NFFT
        fac=sqrt(1.0/avep)
        cspread=fac*cspread                   !Normalize to constant avg power
        cdat(1:NFFT)=cspread*cdat(1:NFFT)     !Apply Rayleigh fading

!        do i=0,NFFT-1
!           p=real(cspread(i))**2 + aimag(cspread(i))**2
!           write(14,3010) i,p,cspread(i)
!3010       format(i8,3f12.6)
!        enddo

     endif

     dat=aimag(cdat) + xnoise                 !Add the generated noise
     fac=32767.0/nsigs
     if(snrdb.ge.90.0) iwave(1:npts)=nint(fac*dat(1:npts))
     if(snrdb.lt.90.0) iwave(1:npts)=nint(rms*dat(1:npts))
     write(10) h,iwave(1:npts)                !Save the .wav file
     close(10)
  enddo

999 end program jt4sim


code:
program jt65

  ! Test the JT65 decoder for WSJT-X

  use options
  use timer_module, only: timer
  use timer_impl, only: init_timer
  use jt65_test
  use readwav

  character c,mode
  logical :: display_help=.false.,nrobust=.false.,single_decode=.false., ljt65apon=.false.
  type(wav_header) :: wav
  integer*2 id2(NZMAX)
  real*4 dd(NZMAX)
  character*80 infile
  character(len=500) optarg
  character*12 mycall,hiscall
  character*6 hisgrid

  type (option) :: long_options(12) = [ &
       option ('aggressive',.true.,'a','aggressiveness [0-10], default AGGR=0','AGGR'), &
       option ('depth',.true.,'d','depth=5 hinted decoding, default DEPTH=0','DEPTH'),  &
       option ('freq',.true.,'f','signal frequency, default FREQ=1270','FREQ'),         &
       option ('help',.false.,'h','Display this help message',''),                      &
       option ('mode',.true.,'m','Mode A, B, C. Default is A.','MODE'),                 &
       option ('ntrials',.true.,'n','number of trials, default TRIALS=10000','TRIALS'), &
       option ('robust-sync',.false.,'r','robust sync',''),                             &
       option ('my-call',.true.,'c','my callsign',''),                                  &
       option ('his-call',.true.,'x','his callsign',''),                                &
       option ('his-grid',.true.,'g','his grid locator',''),                            &
       option ('experience-decoding',.true.,'X'                                         &
               ,'experience decoding options (1..n), default FLAGS=0','FLAGS'),         &
       option ('single-signal-mode',.false.,'s','decode at signal frequency only','') ]

  naggressive=10
  nfqso=1500
  ntrials=100000
  nexp_decode=0
  ntol=20
  nsubmode=0
  nlow=200
  nhigh=4000
  n2pass=1
  ndepth=1
  nQSOProgress=6

  do
     call getopt('a:d:f:hm:n:rc:x:g:X:s',long_options,c,optarg,narglen,nstat,noffset,nremain,.true.)
     if( nstat .ne. 0 ) then
        exit
     end if
     select case (c)
     case ('a')
        read (optarg(:narglen), *) naggressive
     case ('d')
        read (optarg(:narglen), *) ndepth
     case ('f')
        read (optarg(:narglen), *) nfqso
     case ('h')
        display_help = .true.
     case ('m')
        read (optarg(:narglen), *) mode
        if( mode .eq. 'b' .or. mode .eq. 'B' ) then
          nsubmode=1
        endif
        if( mode .eq. 'c' .or. mode .eq. 'C' ) then
          nsubmode=2
        endif
     case ('n')
        read (optarg(:narglen), *) ntrials
     case ('r')
        nrobust=.true.
     case ('c')
        read (optarg(:narglen), *) mycall
     case ('x')
        read (optarg(:narglen), *) hiscall
     case ('g')
        read (optarg(:narglen), *) hisgrid
     case ('X')
        read (optarg(:narglen), *) nexp_decode
     case ('s')
        single_decode=.true.
        ntol=100
        nlow=nfqso-ntol
        nhigh=nfqso+ntol
        n2pass=1
     end select
  end do

  if(single_decode) nexp_decode=ior(nexp_decode,32)
  if(display_help .or. nstat.lt.0 .or. nremain.lt.1) then
     print *, ''
     print *, 'Usage: jt65 [OPTIONS] file1 [file2 ...]'
     print *, ''
     print *, '       JT65 decode pre-recorded .WAV file(s)'
     print *, ''
     print *, 'OPTIONS:'
     print *, ''
     do i = 1, size (long_options)
        call long_options(i) % print (6)
     end do
     go to 999
  endif

  call init_timer ('timer.out')
  call timer('jt65    ',0)

  ndecoded=0
  do ifile=noffset+1,noffset+nremain
     nfa=nlow
     nfb=nhigh
     minsync=0
     call get_command_argument(ifile,optarg,narglen)
     infile=optarg(:narglen)
     call timer('read    ',0)
     call wav%read (infile)
     i1=index(infile,'.wav')
     if( i1 .eq. 0 ) i1=index(infile,'.WAV')
     read(infile(i1-4:i1-1),*,err=998) nutc
     npts=52*12000
     read(unit=wav%lun) id2(1:npts)
     close(unit=wav%lun)
     call timer('read    ',1)
     dd(1:npts)=id2(1:npts)
     dd(npts+1:)=0.
     call test(dd,nutc,nfa,nfb,nfqso,ntol,nsubmode, &
          n2pass,nrobust,ntrials,naggressive,ndepth, &
          mycall,hiscall,hisgrid,nexp_decode,nQSOProgress,ljt65apon)
!     if(nft.gt.0) exit
  enddo

  call timer('jt65    ',1)
  call timer('jt65    ',101)
  !  call four2a(a,-1,1,1,1)                  !Free the memory used for plans
  !  call filbig(a,-1,1,0.0,0,0,0,0,0)        ! (ditto)
  go to 999

998 print*,'Cannot read from file:'
  print*,infile

999 continue
end program jt65


code:
module jt65_decode

  integer, parameter :: NSZ=3413, NZMAX=60*12000

  type :: jt65_decoder
     procedure(jt65_decode_callback), pointer :: callback => null()
   contains
     procedure :: decode
  end type jt65_decoder

! Callback function to be called with each decode
  abstract interface
     subroutine jt65_decode_callback(this,sync,snr,dt,freq,drift,     &
          nflip,width,decoded,ft,qual,nsmo,nsum,minsync)

       import jt65_decoder
       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) :: decoded
       integer, intent(in) :: ft
       integer, intent(in) :: qual
       integer, intent(in) :: nsmo
       integer, intent(in) :: nsum
       integer, intent(in) :: minsync

     end subroutine jt65_decode_callback
  end interface

contains

  subroutine decode(this,callback,dd0,npts,newdat,nutc,nf1,nf2,nfqso,     &
       ntol,nsubmode,minsync,nagain,n2pass,nrobust,ntrials,naggressive,   &
       ndepth,emedelay,clearave,mycall,hiscall,hisgrid,nexp_decode,       &
       nQSOProgress,ljt65apon)

!  Process dd0() data to find and decode JT65 signals.

    use jt65_mod
    use timer_module, only: timer

    include 'constants.f90'

    class(jt65_decoder), intent(inout) :: this
    procedure(jt65_decode_callback) :: callback
    real, intent(in) :: dd0(NZMAX),emedelay
    integer, intent(in) :: npts, nutc, nf1, nf2, nfqso, ntol     &
         , nsubmode, minsync, n2pass, ntrials, naggressive, ndepth      &
         , nexp_decode, nQSOProgress
    logical, intent(in) :: newdat, nagain, nrobust, clearave, ljt65apon
    character(len=12), intent(in) :: mycall, hiscall
    character(len=6), intent(in) :: hisgrid

    real dd(NZMAX)
    real ss(552,NSZ)
    real savg(NSZ)
    real a(5)
    character*22 decoded,decoded0,avemsg,deepave
    type candidate
       real freq
       real dt
       real sync
       real flip
    end type candidate
    type(candidate) ca(300)
    type accepted_decode
       real freq
       real dt
       real sync
       character*22 decoded
    end type accepted_decode
    type(accepted_decode) dec(50)
    logical :: first_time,prtavg,single_decode,bVHF

    integer h0(0:11),d0(0:11)
    real r0(0:11)
    common/decstats/ntry65a,ntry65b,n65a,n65b,num9,numfano
    common/steve/thresh0
    common/sync/ss

!            0  1  2  3  4  5  6  7  8  9 10 11
    data h0/41,42,43,43,44,45,46,47,48,48,49,49/
    data d0/71,72,73,74,76,77,78,80,81,82,83,83/

!             0    1    2    3    4    5    6    7    8    9   10   11
    data r0/0.70,0.72,0.74,0.76,0.78,0.80,0.82,0.84,0.86,0.88,0.90,0.90/
    data nutc0/-999/,nfreq0/-999/,nsave/0/
    save

    this%callback => callback
    first_time=newdat
    dd=dd0
    ndecoded=0

    if(nsubmode.ge.100) then
! This is QRA64 mode
       mode64=2**(nsubmode-100)
!###
!       open(60,file='qra64_data.bin',access='stream',position='append')
!       write(60) dd,npts,nutc,nf1,nf2,nfqso,ntol,mode64,minsync,ndepth,   &
!            mycall,hiscall,hisgrid
!       close(60)
!###
       call qra64a(dd,npts,nutc,nf1,nf2,nfqso,ntol,mode64,minsync,ndepth,   &
            emedelay,mycall,hiscall,hisgrid,sync,nsnr,dtx,nfreq,decoded,nft)
       if (associated(this%callback)) then
          ndrift=0
          nflip=1
          width=1.0
          nsmo=0
          nqual=0
          call this%callback(sync,nsnr,dtx,nfreq,ndrift,  &
               nflip,width,decoded,nft,nqual,nsmo,1,minsync)
       end if
       go to 900
    endif

    single_decode=iand(nexp_decode,32).ne.0 .or. nagain
    bVHF=iand(nexp_decode,64).ne.0

    if( bVHF ) then
      nvec=ntrials
      npass=1
      if(n2pass.gt.1) npass=2
    else
      nvec=1000
      if(ndepth.eq.1) then
         npass=2
         nvec=100
      elseif(ndepth.eq.2) then
         npass=2
         nvec=1000
      else
         npass=4
         nvec=1000
      endif
    endif
    do ipass=1,npass
       first_time=.true.
       if(ipass.eq.1) then                        !First-pass parameters
          thresh0=2.5
          nsubtract=1
          nrob=0
       elseif( ipass.eq.2 ) then                  !Second-pass parameters
          thresh0=2.0
          nsubtract=1
          nrob=0
       elseif( ipass.eq.3 ) then
          thresh0=2.0
          nsubtract=1
          nrob=0
       elseif( ipass.eq.4 ) then
          thresh0=2.0
          nsubtract=0
          nrob=1
       endif
       if(npass.eq.1) then
         nsubtract=0
         thresh0=2.0
       endif

       call timer('symsp65 ',0)
       ss=0.
       call symspec65(dd,npts,nqsym,savg)    !Get normalized symbol spectra
       call timer('symsp65 ',1)
       nfa=nf1
       nfb=nf2

!### Q: should either of the next two uses of "single_decode" be "bVHF" instead?
       if(single_decode .or. (bVHF .and. ntol.lt.1000)) then
          nfa=max(200,nfqso-ntol)
          nfb=min(4000,nfqso+ntol)
          thresh0=1.0
       endif
       df=12000.0/8192.0                     !df = 1.465 Hz
       if(bVHF) then
          ia=max(1,nint(nfa/df)-ntol)
          ib=min(NSZ,nint(nfb/df)+ntol)
          nz=ib-ia+1
          call lorentzian(savg(ia),nz,a)
          baseline=a(1)
          amp=a(2)
          f0=(a(3)+ia-1)*df
          width=a(4)*df
       endif

       ncand=0
       call timer('sync65  ',0)
       call sync65(nfa,nfb,naggressive,ntol,nqsym,ca,ncand,nrob,bVHF)
       call timer('sync65  ',1)

! If a candidate was found within +/- ntol of nfqso, move it into ca(1).
       call fqso_first(nfqso,ntol,ca,ncand)
       if(single_decode) then
          if(ncand.eq.0) ncand=1
          if(abs(ca(1)%freq - f0).gt.width) width=2*df    !### ??? ###
       endif

       mode65=2**nsubmode
       nflip=1
       nqd=0
       decoded='                      '
       decoded0=""
       freq0=0.
       prtavg=.false.
       if(.not.nagain) nsum=0
       if(clearave) then
          nsum=0
          nsave=0
       endif

       if(bVHF) then
! Be sure to search for shorthand message at nfqso +/- ntol
          if(ncand.lt.300) ncand=ncand+1
          ca(ncand)%sync=5.0
          ca(ncand)%dt=2.5
          ca(ncand)%freq=nfqso
       endif
       do icand=1,ncand
          sync1=ca(icand)%sync
          dtx=ca(icand)%dt
          freq=ca(icand)%freq
          if(bVHF) then
             flip=ca(icand)%flip
             nflip=flip
          endif
          if(sync1.lt.float(minsync)) nflip=0
          if(ipass.eq.1) ntry65a=ntry65a + 1
          if(ipass.eq.2) ntry65b=ntry65b + 1
          call timer('decod65a',0)
          nft=0
          nspecial=0
          call decode65a(dd,npts,first_time,nqd,freq,nflip,mode65,nvec,     &
               naggressive,ndepth,ntol,mycall,hiscall,hisgrid,nQSOProgress, &
               ljt65apon,nexp_decode,bVHF,sync2,a,dtx,nft,nspecial,qual,     &
               nhist,nsmo,decoded)
          if(nspecial.eq.2) decoded='RO'
          if(nspecial.eq.3) decoded='RRR'
          if(nspecial.eq.4) decoded='73'
          call timer('decod65a',1)
          if(sync1.lt.float(minsync) .and.                                  &
               decoded.eq.'                      ') nflip=0
          if(nft.ne.0) nsum=1

          nhard_min=param(1)
          nrtt1000=param(4)
          ntotal_min=param(5)
          nsmo=param(9)

          nfreq=nint(freq+a(1))
          ndrift=nint(2.0*a(2))
          if(bVHF) then
            xtmp=10**((sync1+16.0)/10.0) ! sync comes to us in dB
            s2db=1.1*db(xtmp)+1.4*(dB(width)-4.3)-52.0
!             s2db=sync1 - 30.0 + db(width/3.3)       !### VHF/UHF/microwave
             if(nspecial.gt.0) s2db=sync2
          else
             s2db=10.0*log10(sync2) - 35             !### Empirical (HF)
          endif
          nsnr=nint(s2db)
          if(nsnr.lt.-30) nsnr=-30
          if(nsnr.gt.-1) nsnr=-1
          nftt=0

!********* DOES THIS STILL WORK WHEN NFT INCLUDES # OF AP SYMBOLS USED??
          if(nft.ne.1 .and. iand(ndepth,16).eq.16 .and. (.not.prtavg)) then
! Single-sequence FT decode failed, so try for an average FT decode.
             if(nutc.ne.nutc0 .or. abs(nfreq-nfreq0).gt.ntol) then
! This is a new minute or a new frequency, so call avg65.
                nutc0=nutc
                nfreq0=nfreq
                nsave=nsave+1
                nsave=mod(nsave-1,64)+1
                call avg65(nutc,nsave,sync1,dtx,nflip,nfreq,mode65,ntol,     &
                     ndepth,nagain,ntrials,naggressive,clearave,neme,mycall, &
                     hiscall,hisgrid,nftt,avemsg,qave,deepave,nsum,ndeepave, &
                     nQSOProgress,ljt65apon)
                nsmo=param(9)
                nqave=qave

                if (associated(this%callback) .and. nsum.ge.2) then
                   call this%callback(sync1,nsnr,dtx-1.0,nfreq,ndrift,  &
                        nflip,width,avemsg,nftt,nqave,nsmo,nsum,minsync)
                   prtavg=.true.
                end if

             endif
          endif

          if(nftt.eq.1) then
!             nft=1
             decoded=avemsg
             go to 5
          endif
          n=naggressive
          rtt=0.001*nrtt1000
          if(nft.lt.2 .and. minsync.ge.0 .and. nspecial.eq.0) then
             if(nhard_min.gt.50) cycle
             if(nhard_min.gt.h0(n)) cycle
             if(ntotal_min.gt.d0(n)) cycle
             if(rtt.gt.r0(n)) cycle
          endif

5         continue
          if(decoded.eq.decoded0 .and. abs(freq-freq0).lt. 3.0 .and.    &
               minsync.ge.0) cycle                  !Don't display dupes
          if(decoded.ne.'                      ' .or. minsync.lt.0) then
             if(nsubtract.eq.1) then
                call timer('subtr65 ',0)
                call subtract65(dd,npts,freq,dtx)
                call timer('subtr65 ',1)
             endif

             ndupe=0 ! de-dedupe
             do i=1, ndecoded
                if(decoded==dec(i)%decoded) then
                   ndupe=1
                   exit
                endif
             enddo
             if(ndupe.ne.1 .and. sync1.ge.float(minsync)) then
                if(ipass.eq.1) n65a=n65a + 1
                if(ipass.eq.2) n65b=n65b + 1
                if(ndecoded.lt.50) ndecoded=ndecoded+1
                dec(ndecoded)%freq=freq+a(1)
                dec(ndecoded)%dt=dtx
                dec(ndecoded)%sync=sync2
                dec(ndecoded)%decoded=decoded
                nqual=min(qual,9999.0)

                if (associated(this%callback)) then
                   call this%callback(sync1,nsnr,dtx-1.0,nfreq,ndrift,  &
                        nflip,width,decoded,nft,nqual,nsmo,1,minsync)
                end if
             endif
             decoded0=decoded
             freq0=freq
             if(decoded0.eq.'                      ') decoded0='*'
          endif
       enddo                                 !Candidate loop
       if(ipass.eq.2 .and. ndecoded.lt.1) exit
    enddo                                    !Multiple-pass loop
900 return
  end subroutine decode

  subroutine avg65(nutc,nsave,snrsync,dtxx,nflip,nfreq,mode65,ntol,ndepth,    &
       nagain, ntrials,naggressive,clearave,neme,mycall,hiscall,hisgrid,nftt, &
       avemsg,qave,deepave,nsum,ndeepave,nQSOProgress,ljt65apon)

! Decodes averaged JT65 data

    use jt65_mod
    parameter (MAXAVE=64)
    character*22 avemsg,deepave,deepbest
    character mycall*12,hiscall*12,hisgrid*6
    character*1 csync,cused(64)
    logical nagain
    integer iused(64)
! Accumulated data for message averaging
    integer iutc(MAXAVE)
    integer nfsave(MAXAVE)
    integer nflipsave(MAXAVE)
    real s1b(-255:256,126)
    real s1save(-255:256,126,MAXAVE)
    real s2(66,126)
    real s3save(64,63,MAXAVE)
    real s3b(64,63)
    real s3c(64,63)
    real dtsave(MAXAVE)
    real syncsave(MAXAVE)
    logical first,clearave,ljt65apon
    data first/.true./
    save

    if(first .or. clearave) then
       iutc=-1
       nfsave=0
       dtdiff=0.2
       first=.false.
       s3save=0.
       s1save=0.
       nsave=1           !### ???
! Silence compiler warnings
       if(nagain .and. ndeepave.eq.-99 .and. neme.eq.-99) stop
    endif

    do i=1,64
       if(iutc(i).lt.0) exit
       if(nutc.eq.iutc(i) .and. abs(nfreq-nfsave(i)).le.ntol) go to 10
    enddo

! Save data for message averaging
    iutc(nsave)=nutc
    syncsave(nsave)=snrsync
    dtsave(nsave)=dtxx
    nfsave(nsave)=nfreq
    nflipsave(nsave)=nflip
    s1save(-255:256,1:126,nsave)=s1
    s3save(1:64,1:63,nsave)=s3a

10  syncsum=0.
    dtsum=0.
    nfsum=0
    nsum=0
    s1b=0.
    s3b=0.
    s3c=0.

    do i=1,MAXAVE                               !Consider all saved spectra
       cused(i)='.'
       if(iutc(i).lt.0) cycle
       if(mod(iutc(i),2).ne.mod(nutc,2)) cycle  !Use only same (odd/even) seq
       if(abs(dtxx-dtsave(i)).gt.dtdiff) cycle  !DT must match
       if(abs(nfreq-nfsave(i)).gt.ntol) cycle   !Freq must match
       if(nflip.ne.nflipsave(i)) cycle          !Sync type (*/#) must match
       s3b=s3b + s3save(1:64,1:63,i)
       s1b=s1b + s1save(-255:256,1:126,i)
       syncsum=syncsum + syncsave(i)
       dtsum=dtsum + dtsave(i)
       nfsum=nfsum + nfsave(i)
       cused(i)='$'
       nsum=nsum+1
       iused(nsum)=i
    enddo
    if(nsum.lt.64) iused(nsum+1)=0

    syncave=0.
    dtave=0.
    fave=0.
    if(nsum.gt.0) then
       syncave=syncsum/nsum
       dtave=dtsum/nsum
       fave=float(nfsum)/nsum
    endif

    do i=1,nsave
       csync='*'
       if(nflipsave(i).lt.0.0) csync='#'
       write(14,1000) cused(i),iutc(i),syncsave(i),dtsave(i)-1.0,nfsave(i),csync
1000   format(a1,i5.4,f6.1,f6.2,i6,1x,a1)
    enddo
    if(nsum.lt.2) go to 900

    nftt=0
    df=1378.125/512.0

! Do the smoothing loop
    qualbest=0.
    minsmo=0
    maxsmo=0
    if(mode65.ge.2) then
       minsmo=nint(width/df)
       maxsmo=2*minsmo
    endif
    nn=0
    do ismo=minsmo,maxsmo
       if(ismo.gt.0) then
          do j=1,126
             call smo121(s1b(-255,j),512)
             if(j.eq.1) nn=nn+1
             if(nn.ge.4) then
                call smo121(s1b(-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)=s1b(jj,1:126)
       enddo

       do j=1,63
          k=mdat(j)                       !Points to data symbol
          if(nflip.lt.0) k=mdat2(j)
          do i=1,64
             s3c(i,j)=4.e-5*s2(i+2,k)
          enddo
       enddo

       nadd=nsum*ismo
       call extract(s3c,nadd,mode65,ntrials,naggressive,ndepth,nflip,mycall, &
            hiscall,hisgrid,nQSOProgress,ljt65apon,nexp_decode,ncount,nhist, &
            avemsg,ltext,nftt,qual)
       if(nftt.eq.1) then
          nsmo=ismo
          param(9)=nsmo
          go to 900
       else if(nftt.eq.2) then
          if(qual.gt.qualbest) then
             deepbest=avemsg
             qualbest=qual
             nnbest=nn
             nsmobest=ismo
             nfttbest=nftt
          endif
       endif
    enddo

    if(nfttbest.eq.2) then
       avemsg=deepbest       !### ???
       deepave=deepbest
       qave=qualbest
       nsmo=nsmobest
       param(9)=nsmo
       nftt=nfttbest
    endif
900 continue

    return
  end subroutine avg65

end module jt65_decode


code:
module jt65_mod

  integer param(0:9)
  integer mrs(63)
  integer mrs2(63)
  integer mdat(126),mref(126,2),mdat2(126),mref2(126,2)    !From prcom

  real s1(-255:256,126)
  real s3a(64,63)
  real pr(126)
  real width

end module jt65_mod


code:
module jt65_test

  ! Test the JT65 decoder for WSJT-X

  implicit none

  public :: test
  integer, parameter, public :: NZMAX=60*12000
  integer, public :: nft

contains

  subroutine test (dd,nutc,nflow,nfhigh,nfqso,ntol,nsubmode,n2pass,nrobust     &
       ,ntrials,naggressive,ndepth,mycall,hiscall,hisgrid,nexp_decode,         &
        nQSOProgress,ljt65apon)
    use timer_module, only: timer
    use jt65_decode
    implicit none

    include 'constants.f90'
    real, intent(in) :: dd(NZMAX)
    integer, intent(in) :: nutc, nflow, nfhigh, nfqso, ntol, nsubmode, n2pass  &
         , ntrials, naggressive, ndepth, nexp_decode, nQSOProgress
    logical, intent(in) :: nrobust,ljt65apon
    character(len=12), intent(in) :: mycall, hiscall
    character(len=6), intent(in) :: hisgrid
    type(jt65_decoder) :: my_decoder
    logical nclearave                          !### Should be a dummy arg?
    nclearave=.false.

    call timer('jt65a   ',0)
    call my_decoder%decode(my_callback,dd,npts=52*12000,newdat=.true.,     &
         nutc=nutc,nf1=nflow,nf2=nfhigh,nfqso=nfqso,ntol=ntol,             &
         nsubmode=nsubmode, minsync=-1,nagain=.false.,n2pass=n2pass,       &
         nrobust=nrobust,ntrials=ntrials,naggressive=naggressive,          &
         ndepth=ndepth,emedelay=0.0,clearave=nclearave,mycall=mycall,      &
         hiscall=hiscall,hisgrid=hisgrid,nexp_decode=nexp_decode,          &
         nQSOProgress=nQSOProgress,ljt65apon=ljt65apon)
    call timer('jt65a   ',1)
  end subroutine test

  subroutine my_callback (this,sync,snr,dt,freq,drift,nflip,width,     &
       decoded,ft,qual,smo,sum,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) :: decoded
    integer, intent(in) :: ft
    integer, intent(in) :: qual
    integer, intent(in) :: smo
    integer, intent(in) :: sum
    integer, intent(in) :: minsync

    integer nwidth
    real t

    if(minsync+nflip+qual.eq.-9999) stop            !Silence compiler warning
    t=max(0.0,width*width-7.2)
    nwidth=max(nint(sqrt(t)),2)
!### deal with nflip here! ###
!### also single_decode, csync, etc... ###
    write(*,1012) nint(sync),snr,dt,freq,drift,nwidth,         &
         decoded,ft,sum,smo
1012 format(i4,i5,f6.2,i5,i4,i3,1x,a22,' JT65',3i3)
    nft=ft
    call flush(6)

  end subroutine my_callback

end module jt65_test


code:
program JT65code

! Provides examples of message packing, bit and symbol ordering,
! Reed Solomon encoding, and other necessary details of the JT65
! protocol.

  use packjt
  character*22 msg,msgchk,msg0,msg1,decoded,cok*3,bad*1,msgtype*10,expected
  integer dgen(12),sent(63),tmp(63),recd(12),era(51)
  include 'testmsg.f90'

  nargs=iargc()
  if(nargs.ne.1) then
     print*,'Usage: jt65code "message"'
     print*,'       jt65code -t'
     go to 999
  endif

  call getarg(1,msg)                     !Get message from command line
  msgchk=msg
  call fmtmsg(msgchk,iz)
  nmsg=1
  if(msg(1:2).eq."-t") then
     if (NTEST+5 > MAXTEST) then
        write(*,*) "NTEST exceed MAXTEST"
     endif
     testmsg(NTEST+1)="KA1ABC WB9XYZ EN34 OOO"
     testmsg(NTEST+2)="KA1ABC WB9XYZ OOO"
     testmsg(NTEST+3)="RO"
     testmsg(NTEST+4)="RRR"
     testmsg(NTEST+5)="73"
     testmsgchk(NTEST+1)="KA1ABC WB9XYZ EN34 OOO"
     testmsgchk(NTEST+2)="KA1ABC WB9XYZ OOO"
     testmsgchk(NTEST+3)="RO"
     testmsgchk(NTEST+4)="RRR"
     testmsgchk(NTEST+5)="73"
     nmsg=NTEST+5
  endif

  write(*,1010)
1010 format("    Message                Decoded              Err? Type          Expected"/   &
            76("-"))

  do imsg=1,nmsg
     if(nmsg.gt.1) then
        msg=testmsg(imsg)
        msgchk=testmsgchk(imsg)
     endif

     call fmtmsg(msg,iz)                    !To upper, collapse mult blanks
     msg0=msg                               !Input message
     call chkmsg(msg,cok,nspecial,flip)     !See if it includes "OOO" report
     msg1=msg                               !Message without "OOO"

     if(nspecial.gt.0) then                  !or is a shorthand message
        if(nspecial.eq.2) decoded="RO"
        if(nspecial.eq.3) decoded="RRR"
        if(nspecial.eq.4) decoded="73"
        itype=-1
        msgtype="Shorthand"
        go to 10
     endif

     call packmsg(msg1,dgen,itype,.false.) !Pack message into 12 six-bit bytes
     msgtype=""
     if(itype.eq.1) msgtype="Std Msg"
     if(itype.eq.2) msgtype="Type 1 pfx"
     if(itype.eq.3) msgtype="Type 1 sfx"
     if(itype.eq.4) msgtype="Type 2 pfx"
     if(itype.eq.5) msgtype="Type 2 sfx"
     if(itype.eq.6) msgtype="Free text"

     call rs_encode(dgen,sent)               !RS encode
     call interleave63(sent,1)               !Interleave channel symbols
     call graycode(sent,63,1,sent)           !Apply Gray code

     call graycode(sent,63,-1,tmp)           !Remove Gray code
     call interleave63(tmp,-1)               !Remove interleaving
     call rs_decode(tmp,era,0,recd,nerr)     !Decode the message
     call unpackmsg(recd,decoded,.false.,'      ') !Unpack the user message
     if(cok.eq."OOO") decoded(20:22)=cok
     call fmtmsg(decoded,iz)

10   bad=" "
     if(decoded.ne.msgchk) bad="*"
     expected = 'EXACT'
     if (msg0.ne.msgchk) expected = 'TRUNCATED'
     if (nmsg.eq.1) expected = 'UNKNOWN'
     write(*,1020) imsg,msg0,decoded,bad,itype,msgtype,expected
1020 format(i2,'.',1x,a22,1x,a22,1x,a1,i3,":",a10,2x,a22)
  enddo

  if(nmsg.eq.1 .and. nspecial.eq.0) then
     write(*,1030) dgen
1030 format(/'Packed message, 6-bit symbols ',12i3) !Display packed symbols

     write(*,1040) sent
1040 format(/'Information-carrying channel symbols'/(i5,20i3))
  endif

999 end program JT65code


code:
program jt65sim

! Generate simulated JT65 data for testing WSJT-X

  use wavhdr
  use packjt
  use options
  parameter (NMAX=54*12000)              ! = 648,000 @12kHz
  parameter (NFFT=10*65536,NH=NFFT/2)
  type(hdr) h                            !Header for .wav file
  integer*2 iwave(NMAX)                  !Generated waveform
  integer*4 itone(126)                   !Channel symbols (values 0-65)
  integer dgen(12)                       !Twelve 6-bit data symbols
  integer sent(63)                       !RS(63,12) codeword
  real*4 xnoise(NMAX)                    !Generated random noise
  real*4 dat(NMAX)                       !Generated real data
  complex cdat(NMAX)                     !Generated complex waveform
  complex cspread(0:NFFT-1)              !Complex amplitude for Rayleigh fading
  complex z
  real*8 f0,dt,twopi,phi,dphi,baud,fsample,freq,sps
  character msg*22,fname*11,csubmode*1,c,optarg*500,numbuf*32
!  character call1*5,call2*5
  logical :: display_help=.false.,seed_prngs=.true.
  type (option) :: long_options(12) = [ &
    option ('help',.false.,'h','Display this help message',''),                                &
    option ('sub-mode',.true.,'m','sub mode, default MODE=A','MODE'),                          &
    option ('num-sigs',.true.,'n','number of signals per file, default SIGNALS=10','SIGNALS'), &
    option ('f0',.true.,'F','base frequency offset, default F0=1500.0','F0'), &
    option ('doppler-spread',.true.,'d','Doppler spread, default SPREAD=0.0','SPREAD'),        &
    option ('time-offset',.true.,'t','Time delta, default SECONDS=0.0','SECONDS'),             &
    option ('num-files',.true.,'f','Number of files to generate, default FILES=1','FILES'),    &
    option ('no-prng-seed',.false.,'p','Do not seed PRNGs (use for reproducible tests)',''),   &
    option ('strength',.true.,'s','S/N in dB (2500Hz reference b/w), default SNR=0','SNR'),    &
    option ('11025',.false.,'S','Generate at 11025Hz sample rate, default 12000Hz',''),        &
    option ('gain-offset',.true.,'G','Gain offset in dB, default GAIN=0dB','GAIN'),            &
    option ('message',.true.,'M','Message text','Message') ]

  integer nprc(126)                                      !Sync pattern
  data nprc/1,0,0,1,1,0,0,0,1,1,1,1,1,1,0,1,0,1,0,0,  &
            0,1,0,1,1,0,0,1,0,0,0,1,1,1,0,0,1,1,1,1,  &
            0,1,1,0,1,1,1,1,0,0,0,1,1,0,1,0,1,0,1,1,  &
            0,0,1,1,0,1,0,1,0,1,0,0,1,0,0,0,0,0,0,1,  &
            1,0,0,0,0,0,0,0,1,1,0,1,0,0,1,0,1,1,0,1,  &
            0,1,0,1,0,0,1,1,0,0,1,0,0,1,0,0,0,0,1,1,  &
            1,1,1,1,1,1/

! Default parameters:
  csubmode='A'
  mode65=1
  nsigs=10
  bf0=1500.
  fspread=0.
  xdt=0.
  snrdb=0.
  nfiles=1
  nsample_rate=12000
  gain_offset=0.
  msg="K1ABC W9XYZ EN37"

  do
     call getopt('hm:n:F:d:t:f:ps:SG:M:',long_options,c,optarg,narglen,nstat,noffset,nremain,.true.)
     if( nstat .ne. 0 ) then
        exit
     end if
     select case (c)
     case ('h')
        display_help = .true.
     case ('m')
        read (optarg(:narglen), *) csubmode
        if(csubmode.eq.'A') mode65=1
        if(csubmode.eq.'B') mode65=2
        if(csubmode.eq.'C') mode65=4
     case ('n')
        read (optarg(:narglen), *,err=10) nsigs
     case ('F')
        read (optarg(:narglen), *,err=10) bf0
     case ('d')
        read (optarg(:narglen), *,err=10) fspread
     case ('t')
        read (optarg(:narglen), *) numbuf
        if (numbuf(1:1) == '\') then !'\'
           read (numbuf(2:), *,err=10) xdt
        else
           read (numbuf, *,err=10) xdt
        end if
     case ('f')
        read (optarg(:narglen), *,err=10) nfiles
     case ('p')
        seed_prngs=.false.
     case ('s')
        read (optarg(:narglen), *) numbuf
        if (numbuf(1:1) == '\') then !'\'
           read (numbuf(2:), *,err=10) snrdb
        else
           read (numbuf, *,err=10) snrdb
        end if
     case ('S')
        nsample_rate=11025
     case ('G')
        read (optarg(:narglen), *) numbuf
        if (numbuf(1:1) == '\') then !'\'
           read (numbuf(2:), *, err=10) gain_offset
        else
           read (numbuf, *, err=10) gain_offset
        end if
     case ('M')
        read (optarg(:narglen), '(A)',err=10) msg
        write(*,*) msg
     end select
     cycle
10   display_help=.true.
     print *, 'Optional argument format error for option -', c
  end do

  if(display_help .or. nstat.lt.0 .or. nremain.ge.1) then
     print *, ''
     print *, 'Usage: jt65sim [OPTIONS]'
     print *, ''
     print *, '       Generate one or more simulated JT65 signals in .WAV file(s)'
     print *, ''
     print *, 'Example: jt65sim -m B -n 10 -d 0.2 -s \\-24.5 -t 0.0 -f 4'
     print *, ''
     print *, 'OPTIONS: NB Use \ (\\ on *nix shells) to escape -ve arguments'
     print *, ''
     do i = 1, size (long_options)
       call long_options(i) % print (6)
     end do
     go to 999
  endif

  if (seed_prngs) then
     call init_random_seed()    ! seed Fortran RANDOM_NUMBER generator
     call sgran()               ! see C rand generator (used in gran)
  end if

  rms=100. * 10. ** (gain_offset / 20.)

  fsample=nsample_rate               !Sample rate (Hz)
  dt=1.d0/fsample                    !Sample interval (s)
  twopi=8.d0*atan(1.d0)
  npts=54*nsample_rate               !Total samples in .wav file
  baud=11025.d0/4096.d0              !Keying rate
  sps=real(nsample_rate)/baud        !Samples per symbol, at fsample=NSAMPLE_RATE Hz
  nsym=126                           !Number of channel symbols
  h=default_header(nsample_rate,npts)
  dfsig=2000.0/nsigs                 !Freq spacing between sigs in file (Hz)

  do ifile=1,nfiles                  !Loop over requested number of files
     write(fname,1002) ifile         !Output filename
1002 format('000000_',i4.4)
     open(10,file=fname//'.wav',access='stream',status='unknown')

     xnoise=0.
     cdat=0.
     if(snrdb.lt.90) then
        do i=1,npts
           xnoise(i)=gran()          !Generate gaussian noise
        enddo
     endif

     do isig=1,nsigs                        !Generate requested number of sigs
        if(mod(nsigs,2).eq.0) f0=bf0 + dfsig*(isig-0.5-nsigs/2)
        if(mod(nsigs,2).eq.1) f0=bf0 + dfsig*(isig-(nsigs+1)/2)
        xsnr=snrdb
        if(snrdb.eq.0.0) xsnr=-19 - isig
        if(csubmode.eq.'B' .and. snrdb.eq.0.0) xsnr=-21 - isig
        if(csubmode.eq.'C' .and. snrdb.eq.0.0) xsnr=-21 - isig

!###
!        call1="K1ABC"
!        ic3=65+mod(isig-1,26)
!        ic2=65+mod((isig-1)/26,26)
!        ic1=65
!        call2="W9"//char(ic1)//char(ic2)//char(ic3)
!        write(msg,1010) call1,call2,nint(xsnr)
!1010    format(a5,1x,a5,1x,i3.2)
!###
        call packmsg(msg,dgen,itype,.false.) !Pack message into 12 six-bit bytes
        call rs_encode(dgen,sent)           !Encode using RS(63,12)
        call interleave63(sent,1)           !Interleave channel symbols
        call graycode65(sent,63,1)          !Apply Gray code

        k=0
        do j=1,nsym                         !Insert sync and data into itone()
           if(nprc(j).eq.0) then
              k=k+1
              itone(j)=sent(k)+2
           else
              itone(j)=0
           endif
        enddo

        bandwidth_ratio=2500.0/(fsample/2.0)
        sig=sqrt(2*bandwidth_ratio)*10.0**(0.05*xsnr)
        if(xsnr.gt.90.0) sig=1.0
        write(*,1020) ifile,isig,f0,csubmode,xsnr,xdt,fspread,msg
1020    format(i4,i4,f10.3,2x,a1,2x,f5.1,f6.2,f5.1,1x,a22)

        phi=0.d0
        dphi=0.d0
        k=nsample_rate + xdt*nsample_rate     !Start audio at t = xdt + 1.0 s
        isym0=-99
        do i=1,npts                         !Add this signal into cdat()
           isym=floor(i/sps)+1
           if(isym.gt.nsym) exit
           if(isym.ne.isym0) then
              freq=f0 + itone(isym)*baud*mode65
              dphi=twopi*freq*dt
              isym0=isym
           endif
           phi=phi + dphi
           if(phi.gt.twopi) phi=phi-twopi
           xphi=phi
           z=cmplx(cos(xphi),sin(xphi))
           k=k+1
           if(k.ge.1) cdat(k)=cdat(k) + sig*z
        enddo
     enddo

     if(fspread.ne.0) then                  !Apply specified Doppler spread
        df=real(nsample_rate)/nfft
        twopi=8*atan(1.0)
        cspread(0)=1.0
        cspread(NH)=0.

! The following options were added 3/15/2016 to make the half-power tone
! widths equal to the requested Doppler spread.  (Previously we effectively
! used b=1.0 and Gaussian shape, which made the tones 1.665 times wider.)
!        b=2.0*sqrt(log(2.0))                     !Gaussian (before 3/15/2016)
!        b=2.0                                    !Lorenzian 3/15 - 3/27
        b=6.0                                     !Lorenzian 3/28 onward

        do i=1,NH
           f=i*df
           x=b*f/fspread
           z=0.
           a=0.
           if(x.lt.3.0) then                          !Cutoff beyond x=3
!              a=sqrt(exp(-x*x))                      !Gaussian
              a=sqrt(1.111/(1.0+x*x)-0.1)             !Lorentzian
              call random_number(r1)
              phi1=twopi*r1
              z=a*cmplx(cos(phi1),sin(phi1))
           endif
           cspread(i)=z
           z=0.
           if(x.lt.50.0) then
              call random_number(r2)
              phi2=twopi*r2
              z=a*cmplx(cos(phi2),sin(phi2))
           endif
           cspread(NFFT-i)=z
        enddo

        do i=0,NFFT-1
           f=i*df
           if(i.gt.NH) f=(i-nfft)*df
           s=real(cspread(i))**2 + aimag(cspread(i))**2
!          write(13,3000) i,f,s,cspread(i)
!3000      format(i5,f10.3,3f12.6)
        enddo
!        s=real(cspread(0))**2 + aimag(cspread(0))**2
!        write(13,3000) 1024,0.0,s,cspread(0)

        call four2a(cspread,NFFT,1,1,1)             !Transform to time domain

        sum=0.
        do i=0,NFFT-1
           p=real(cspread(i))**2 + aimag(cspread(i))**2
           sum=sum+p
        enddo
        avep=sum/NFFT
        fac=sqrt(1.0/avep)
        cspread=fac*cspread                   !Normalize to constant avg power
        cdat(1:npts)=cspread(1:npts)*cdat(1:npts) !Apply Rayleigh fading

!        do i=0,NFFT-1
!           p=real(cspread(i))**2 + aimag(cspread(i))**2
!           write(14,3010) i,p,cspread(i)
!3010       format(i8,3f12.6)
!        enddo

     endif

     dat=aimag(cdat) + xnoise                 !Add the generated noise
     if(snrdb.lt.90.0) then
       dat=rms*dat(1:npts)
     else
       datpk=maxval(abs(dat(1:npts)))
       fac=32766.9/datpk
       dat(1:npts)=fac*dat(1:npts)
     endif
     if(any(abs(dat(1:npts)).gt.32767.0)) print*,"Warning - data will be clipped."
     iwave(1:npts)=nint(dat(1:npts))
     write(10) h,iwave(1:npts)                !Save the .wav file
     close(10)
  enddo

999 end program jt65sim


code:
program jt9

! Decoder for JT9.  Can run stand-alone, reading data from *.wav files;
! or as the back end of wsjt-x, with data placed in a shared memory region.

  use options
  use prog_args
  use, intrinsic :: iso_c_binding
  use FFTW3
  use timer_module, only: timer
  use timer_impl, only: init_timer, fini_timer
  use readwav

  include 'jt9com.f90'

  integer(C_INT) iret
  type(wav_header) wav
  real*4 s(NSMAX)
  character c
  character(len=500) optarg, infile
  character wisfile*80
!### ndepth was defined as 60001.  Why???
  integer :: arglen,stat,offset,remain,mode=0,flow=200,fsplit=2700,          &
       fhigh=4000,nrxfreq=1500,ntrperiod=1,ndepth=1,nexp_decode=0
  logical :: read_files = .true., tx9 = .false., display_help = .false.
  type (option) :: long_options(21) = [ &
    option ('help', .false., 'h', 'Display this help message', ''),          &
    option ('shmem',.true.,'s','Use shared memory for sample data','KEY'),   &
    option ('tr-period', .true., 'p', 'Tx/Rx period, default MINUTES=1',     &
        'MINUTES'),                                                          &
    option ('executable-path', .true., 'e',                                  &
        'Location of subordinate executables (KVASD) default PATH="."',      &
        'PATH'),                                                             &
    option ('data-path', .true., 'a',                                        &
        'Location of writeable data files, default PATH="."', 'PATH'),       &
    option ('temp-path', .true., 't',                                        &
        'Temporary files path, default PATH="."', 'PATH'),                   &
    option ('lowest', .true., 'L',                                           &
        'Lowest frequency decoded (JT65), default HERTZ=200', 'HERTZ'),      &
    option ('highest', .true., 'H',                                          &
        'Highest frequency decoded, default HERTZ=4007', 'HERTZ'),           &
    option ('split', .true., 'S',                                            &
        'Lowest JT9 frequency decoded, default HERTZ=2700', 'HERTZ'),        &
    option ('rx-frequency', .true., 'f',                                     &
        'Receive frequency offset, default HERTZ=1500', 'HERTZ'),            &
    option ('patience', .true., 'w',                                         &
        'FFTW3 planing patience (0-4), default PATIENCE=1', 'PATIENCE'),     &
    option ('fft-threads', .true., 'm',                                      &
        'Number of threads to process large FFTs, default THREADS=1',        &
        'THREADS'),                                                          &
    !option ('jt65', .false., '6', 'JT65 mode', ''),                          &
    !option ('jt9', .false., '9', 'JT9 mode', ''),                            &
    option ('js8', .false., '8', 'JS8 mode', ''),                            &
    !option ('jt4', .false., '4', 'JT4 mode', ''),                            &
    !option ('qra64', .false., 'q', 'QRA64 mode', ''),                        &
    option ('sub-mode', .true., 'b', 'Sub mode, default SUBMODE=A', 'A'),    &
    option ('depth', .true., 'd',                                            &
        'Decoding depth (1-3), default DEPTH=1', 'DEPTH'),               &
    option ('tx-jt9', .false., 'T', 'Tx mode is JT9', ''),                   &
    option ('my-call', .true., 'c', 'my callsign', 'CALL'),                  &
    option ('my-grid', .true., 'G', 'my grid locator', 'GRID'),              &
    option ('his-call', .true., 'x', 'his callsign', 'CALL'),                &
    option ('his-grid', .true., 'g', 'his grid locator', 'GRID'),            &
    option ('experience-decode', .true., 'X',                                &
        'experience based decoding flags (1..n), default FLAGS=0',           &
        'FLAGS') ]

  type(dec_data), allocatable :: shared_data
  character(len=20) :: datetime=''
  character(len=12) :: mycall='K1ABC', hiscall='W9XYZ'
  character(len=6) :: mygrid='', hisgrid='EN37'
  common/patience/npatience,nthreads
  common/decstats/ntry65a,ntry65b,n65a,n65b,num9,numfano
  data npatience/1/,nthreads/1/

  nsubmode = 0

  do
     call getopt('hs:e:a:b:r:m:p:d:f:w:t:9864qTL:S:H:c:G:x:g:X:',      &
          long_options,c,optarg,arglen,stat,offset,remain,.true.)
     if (stat .ne. 0) then
        exit
     end if
     select case (c)
        case ('h')
           display_help = .true.
        case ('s')
           read_files = .false.
           shm_key = optarg(:arglen)
        case ('e')
           exe_dir = optarg(:arglen)
        case ('a')
           data_dir = optarg(:arglen)
        case ('b')
           nsubmode = ichar (optarg(:1)) - ichar ('A')
        case ('t')
           temp_dir = optarg(:arglen)
        case ('m')
           read (optarg(:arglen), *) nthreads
        case ('p')
           read (optarg(:arglen), *) ntrperiod
        case ('d')
           read (optarg(:arglen), *) ndepth
        case ('f')
           read (optarg(:arglen), *) nrxfreq
        case ('L')
           read (optarg(:arglen), *) flow
        case ('S')
           read (optarg(:arglen), *) fsplit
        case ('H')
           read (optarg(:arglen), *) fhigh
        !case ('q')
        !   mode = 164
        !case ('4')
        !   mode = 4
        !case ('6')
        !   if (mode.lt.65) mode = mode + 65
        !case ('9')
        !   if (mode.lt.9.or.mode.eq.65) mode = mode + 9
        case ('8')
           mode = 8
        case ('T')
           tx9 = .true.
        case ('w')
           read (optarg(:arglen), *) npatience
        case ('c')
           read (optarg(:arglen), *) mycall
        case ('G')
           read (optarg(:arglen), *) mygrid
        case ('x')
           read (optarg(:arglen), *) hiscall
        case ('g')
           read (optarg(:arglen), *) hisgrid
        case ('X')
           read (optarg(:arglen), *) nexp_decode
     end select
  end do

  if (display_help .or. stat .lt. 0                      &
       .or. (.not. read_files .and. remain .gt. 0)       &
       .or. (read_files .and. remain .lt. 1)) then

     print *, 'Usage: js8 [OPTIONS] file1 [file2 ...]'
     print *, '       Reads data from *.wav files.'
     print *, ''
     print *, '       js8 -s <key> [-w patience] [-m threads] [-e path] [-a path] [-t path]'
     print *, '       Gets data from shared memory region with key==<key>'
     print *, ''
     print *, 'OPTIONS:'
     print *, ''
     do i = 1, size (long_options)
       call long_options(i) % print (6)
     end do
     go to 999
  endif

  iret=fftwf_init_threads()            !Initialize FFTW threading

! Default to 1 thread, but use nthreads for the big ones
  call fftwf_plan_with_nthreads(1)

! Import FFTW wisdom, if available
  wisfile=trim(data_dir)//'/jt9_wisdom.dat'// C_NULL_CHAR
  iret=fftwf_import_wisdom_from_filename(wisfile)

  ntry65a=0
  ntry65b=0
  n65a=0
  n65b=0
  num9=0
  numfano=0

  if (.not. read_files) then
     call jt9a()          !We're running under control of WSJT-X
     go to 999
  endif

  allocate(shared_data)
  nflatten=0

  do iarg = offset + 1, offset + remain
     call get_command_argument (iarg, optarg, arglen)
     infile = optarg(:arglen)
     call wav%read (infile)
     nfsample=wav%audio_format%sample_rate
     i1=index(infile,'.wav')
     if(i1.lt.1) i1=index(infile,'.WAV')
     if(infile(i1-5:i1-5).eq.'_') then
        read(infile(i1-4:i1-1),*,err=1) nutc
     else
        read(infile(i1-6:i1-1),*,err=1) nutc
     endif
     go to 2
1    nutc=0
2    nsps=0
     if(ntrperiod.eq.1)  then
        nsps=6912
        shared_data%params%nzhsym=181
     else if(ntrperiod.eq.2)  then
        nsps=15360
        shared_data%params%nzhsym=178
     else if(ntrperiod.eq.5)  then
        nsps=40960
        shared_data%params%nzhsym=172
     else if(ntrperiod.eq.10) then
        nsps=82944
        shared_data%params%nzhsym=171
     else if(ntrperiod.eq.30) then
        nsps=252000
        shared_data%params%nzhsym=167
     endif
     if(nsps.eq.0) stop 'Error: bad TRperiod'

     kstep=nsps/2
     k=0
     nhsym0=-999
     npts=(60*ntrperiod-6)*12000
     if(iarg .eq. offset + 1) then
        call init_timer (trim(data_dir)//'/timer.out')
        call timer('jt9     ',0)
     endif

     shared_data%id2=0          !??? Why is this necessary ???

     do iblk=1,npts/kstep
        k=iblk*kstep
        if(mode.eq.8 .and. k.gt.179712) exit
        call timer('read_wav',0)
        read(unit=wav%lun,end=3) shared_data%id2(k-kstep+1:k)
        go to 4
3       call timer('read_wav',1)
        print*,'EOF on input file ',infile
        exit
4       call timer('read_wav',1)
        nhsym=(k-2048)/kstep
        if(nhsym.ge.1 .and. nhsym.ne.nhsym0) then
           if(mode.eq.9 .or. mode.eq.74) then
! Compute rough symbol spectra for the JT9 decoder
              ingain=0
              call timer('symspec ',0)
              nminw=1
              call symspec(shared_data,k,ntrperiod,nsps,ingain,nminw,pxdb,  &
                   s,df3,ihsym,npts8,pxdbmax)
              call timer('symspec ',1)
           endif
           nhsym0=nhsym
           if(nhsym.ge.181) exit
        endif
     enddo
     close(unit=wav%lun)
     shared_data%params%nutc=nutc
     shared_data%params%ndiskdat=.true.
     shared_data%params%ntr=60
     shared_data%params%nfqso=nrxfreq
     shared_data%params%newdat=.true.
     shared_data%params%npts8=74736
     shared_data%params%nfa=flow
     shared_data%params%nfsplit=fsplit
     shared_data%params%nfb=fhigh
     shared_data%params%ntol=20
     shared_data%params%kin=64800
     shared_data%params%nzhsym=181
     shared_data%params%ndepth=ndepth
     shared_data%params%lft8apon=.true.
     shared_data%params%ljt65apon=.true.
     shared_data%params%napwid=75
     shared_data%params%dttol=3.

!     shared_data%params%minsync=0       !### TEST ONLY
!     shared_data%params%nfqso=1500     !### TEST ONLY
!     mycall="G3WDG       "              !### TEST ONLY
!     hiscall="VK7MO       "             !### TEST ONLY
!     hisgrid="QE37        "             !### TEST ONLY
     if(mode.eq.164 .and. nsubmode.lt.100) nsubmode=nsubmode+100

     shared_data%params%naggressive=0
     shared_data%params%n2pass=2
!     shared_data%params%nranera=8                      !### ntrials=10000
     shared_data%params%nranera=6                      !### ntrials=3000
     shared_data%params%nrobust=.false.
     shared_data%params%nexp_decode=nexp_decode
     shared_data%params%mycall=transfer(mycall,shared_data%params%mycall)
     shared_data%params%mygrid=transfer(mygrid,shared_data%params%mygrid)
     shared_data%params%hiscall=transfer(hiscall,shared_data%params%hiscall)
     shared_data%params%hisgrid=transfer(hisgrid,shared_data%params%hisgrid)
     if (tx9) then
        shared_data%params%ntxmode=9
     else
        shared_data%params%ntxmode=65
     end if
     if (mode.eq.0) then
        shared_data%params%nmode=65+9
     else
        shared_data%params%nmode=mode
     end if
     shared_data%params%nsubmode=nsubmode
     datetime="2013-Apr-16 15:13" !### Temp
     shared_data%params%datetime=transfer(datetime,shared_data%params%datetime)
     if(mode.eq.9 .and. fsplit.ne.2700) shared_data%params%nfa=fsplit
     call multimode_decoder(shared_data%ss,shared_data%id2,shared_data%params,nfsample)
  enddo

  call timer('jt9     ',1)
  call timer('jt9     ',101)

999 continue
! Output decoder statistics
  call fini_timer ()
!  open (unit=12, file=trim(data_dir)//'/timer.out', status='unknown', position='append')
!  write(12,1100) n65a,ntry65a,n65b,ntry65b,numfano,num9
!1100 format(58('-')/'   JT65_1  Tries_1  JT65_2 Tries_2    JT9   Tries'/  &
!            58('-')/6i8)

! Save wisdom and free memory
  iret=fftwf_export_wisdom_to_filename(wisfile)
  call four2a(a,-1,1,1,1)
  call filbig(a,-1,1,0.0,0,0,0,0,0)        !used for FFT plans
  call fftwf_cleanup_threads()
  call fftwf_cleanup()
end program jt9