*-----------------------------------------------------------------------
* scribe.f - read a Portable Gray Map, write a MIDI music file
* by Andy Allinger, 2019-2022, released to the public domain
*
*    Permission  to  use, copy, modify, and distribute this software and
*    its documentation  for  any  purpose  and  without  fee  is  hereby
*    granted,  without any conditions or restrictions.  This software is
*    provided "as is" without express or implied warranty.
*
*_______________________OPTIONS_________________________________________
*   --tune [+|-]x                   performance is sharp of standard A=440
*                                     tuning by x semitones (negative for flat)
*   --patch #                       select General Midi patch/program #
*   --text hhh                      a comment of any type
*   --copyright hhh                 legal ownership of the file
*   --seqname hhh                   a title or description
*-----------------------------------------------------------------------
      PROGRAM SCRIBE
       IMPLICIT NONE

*                Constants
       INTEGER DUR,     ! max file length [seconds]
     $         ELIM,    ! length of event arrays
     $         LLIM,    ! max number of pixels
     $         PBR,     ! pitch bend range [cents]
     $         PXPS,    ! Pixels per second
     $         VERT,    ! height of image [pixels]
     $         WLIM     ! max width of image [pixels]

       PARAMETER (DUR = 5 * 60,        ! five minutes
     $            ELIM = 144000,
     $            PBR = 200,           ! two semitones
     $            PXPS = 240,
     $            VERT = 640)
       PARAMETER (WLIM = DUR * PXPS)
       PARAMETER (LLIM = WLIM * VERT)  ! bound of input array

*                     Variables
       INTEGER EK(ELIM), ET(ELIM),     ! event kinds, times
     $         EP1(ELIM), EP2(ELIM),   ! event parameters
     $         ETOT,                   ! # events
     $         I,                      ! counter
     $         HEIGHT,                 ! height of image
     $         IERR,                   ! error code
     $         L,                      ! a string length
     $         MAXBRT,                 ! value of max intensity
     $         NARG,                   ! number of command line arguments
     $         OFFSET,                 ! position in text buffer
     $         PATCH,                  ! instrument to play
     $         WIDTH                   ! width of image

       REAL A(LLIM),              ! bitmap array
     $      RAXBRT,               ! maxbrt as REAL
     $      TRNPOS                ! amount to adjust tuning

       LOGICAL LEXIST    ! file exists

       CHARACTER*256 IFIL, OFIL,   ! filenames
     $               OPT           ! command line option

       CHARACTER*1024 TEXT

*            External functions
       INTEGER LTRIM, IARGC

*-----------------------------------------------------------------------
*        Begin.  Read command line.
*-----------------------------------------------------------------------
       IERR = 0
       ETOT = 0
       OFFSET = 0
       TRNPOS = 0.

       NARG = IARGC()
       IF (NARG < 2) STOP 'files not specified'
       IF (NARG > 99) STOP 'too many arguments!'
       CALL GETARG (NARG-1, IFIL)      ! read file name from command line
       CALL GETARG (NARG, OFIL)

*         Check files for existence
       INQUIRE (FILE=IFIL, EXIST=LEXIST, IOSTAT=IERR)
       IF (IERR .NE. 0) STOP 'trouble inquiring of input file!'
       IF (.NOT. LEXIST) STOP 'file not found'
       INQUIRE (FILE=OFIL, EXIST=LEXIST, IOSTAT=IERR)
       IF (IERR .NE. 0) STOP 'trouble inquiring of output file!'
       IF (LEXIST) THEN
         L = LTRIM(OFIL)
         PRINT *, 'File:  ', OFIL(1:L), ' already exists.  ',
     $              'Please choose another name.'
         STOP
       END IF

*          Options processing loop
       DO I = 1, NARG-2
         CALL GETARG (I, OPT)
         CALL UCASE (OPT)

*             Text events.  Beware the nonstandard code numbers.
         IF (INDEX(OPT, '--TEXT') > 0) THEN
           CALL TXTEVT (OPT,TEXT,I,ET,EK,EP1,EP2,ELIM,ETOT,OFFSET,61)
         ELSE IF (INDEX(OPT, '--COPYRIGHT') > 0) THEN
           CALL TXTEVT (OPT,TEXT,I,ET,EK,EP1,EP2,ELIM,ETOT,OFFSET,62)
         ELSE IF (INDEX(OPT, '--SEQNAME') > 0) THEN
           CALL TXTEVT (OPT,TEXT,I,ET,EK,EP1,EP2,ELIM,ETOT,OFFSET,63)

*             Patch change
         ELSE IF (INDEX(OPT, '--PATCH') > 0) THEN
           CALL GETARG (I+1, OPT)
           READ (UNIT=OPT, FMT=*, IOSTAT=IERR) PATCH
           IF (IERR .NE. 0) PRINT *, 'trouble reading patch change!'
           IF (PATCH < 0 .OR. PATCH > 127)
     $          STOP 'patch numbers are from 0 to 127'
           ETOT = ETOT + 1
           ET(ETOT) = 0
           EK(ETOT) = 55
           EP1(ETOT) = PATCH

         ELSE IF (INDEX(OPT, '--TUNE') > 0) THEN
           CALL GETARG (I+1, OPT)
           READ (UNIT=OPT, FMT=*, IOSTAT=IERR) TRNPOS
           IF (IERR .NE. 0) PRINT *, 'trouble reading retuning option!'
         END IF
       END DO

*-----------------------------------------------------------------------
*         Read in a standard .PGM file
*-----------------------------------------------------------------------
       CALL RPGM (IFIL, A, LLIM, WIDTH, HEIGHT, MAXBRT, IERR)
       IF (IERR .NE. 0) THEN
         PRINT *, 'trouble reading from file: ', IFIL
         STOP
       END IF
       IF (HEIGHT .NE. VERT) THEN
         PRINT *, 'scribe: height must be ', VERT, ' pixels!'
         STOP
       END IF

       CALL FLIPUD (A, WIDTH, HEIGHT)      ! turn upside-down
       RAXBRT = FLOAT(MAXBRT)

       ETOT = ETOT + 1                    ! set pitch bend range
       ET(ETOT) = 0
       EK(ETOT) = 76
       EP1(ETOT) = PBR

*-----------------------------------------------------------------------
*          Interpret
*-----------------------------------------------------------------------
       CALL CLEF (A, WIDTH, HEIGHT, TRNPOS, RAXBRT,
     $             ET, EK, EP1, EP2, ELIM, ETOT, IERR)
       IF (IERR .NE. 0) PRINT *, 'clef returns #', IERR
       IF (ETOT .GE. ELIM) PRINT *, 'scribe:  out of memory'

*-----------------------------------------------------------------------
*          Write out
*-----------------------------------------------------------------------
       CALL WMIDI (OFIL, TEXT, ETOT, ET, EK, EP1, EP2, IERR)
       IF (IERR .NE. 0) PRINT *, 'wmidi returns #', IERR

       STOP 'program complete'
      END  ! of scribe


*-----------------------------------------------------------------------
*         Make a text event
*-----------------------------------------------------------------------
      SUBROUTINE TXTEVT (OPT,TEXT,I,ET,EK,EP1,EP2,ELIM,ETOT,OFFSET,ID)
       IMPLICIT NONE
       CHARACTER*(*) OPT
       CHARACTER*(*) TEXT
       INTEGER I, ELIM, ETOT, OFFSET, ID
       INTEGER ET(ELIM), EK(ELIM), EP1(ELIM), EP2(ELIM)

       INTEGER BUFL
       PARAMETER (BUFL = 1024)
       INTEGER L          ! local variable
       INTEGER LTRIM      ! external function

       CALL GETARG (I+1, OPT)
       L = LTRIM(OPT)
       IF (OFFSET + L > BUFL) STOP 'too much text'
       ETOT = ETOT + 1
       ET(ETOT) = 0
       EK(ETOT) = ID
       EP1(ETOT) = OFFSET + 1
       EP2(ETOT) = L
       TEXT(OFFSET+1:OFFSET+L) = OPT(1:L)
       OFFSET = OFFSET + L
       RETURN
      END  ! of txtevt


*-----------------------------------------------------------------------
* flipud - flip over the image so that high frequencies have high indices
*
*___Name___________Type______In/Out___Description_______________________
*   A(HORZ,VERT)   Real      In       Bitmap array
*   HORZ           Integer   Out      Width of image
*   VERT           Integer   Out      Height of image
*-----------------------------------------------------------------------
      SUBROUTINE FLIPUD (A, HORZ, VERT)
       IMPLICIT NONE
       INTEGER HORZ, VERT
       REAL A(HORZ,VERT)

       INTEGER I, J, K
       REAL SWAP

*          Begin.
       J = VERT
       DO I = 1, VERT/2
         DO K = 1, HORZ
           SWAP = A(K,I)
           A(K,I) = A(K,J)
           A(K,J) = SWAP
         END DO
         J = J - 1
       END DO

       RETURN
      END  ! of FLIPUD


*-----------------------------------------------------------------------
* clef - Separate pixels into dark or light.  Make midi events.
*
*___Name______________Type______In/Out___Description____________________
*   A(WIDTH,HEIGHT)   Real      In       Bitmap of pixel intensities
*   WIDTH             Integer   In       Width of image
*   HEIGHT            Integer   In       Height of image
*   TRNPOS            Real      In       Amount to adjust tuning
*   RAXBRT            Real      In       Maxval from .PPM file
*   ET(ELIM)          Integer   Out      Event times
*   EK(ELIM)          Integer   Out      Event Kinds (refer to WMIDI)
*   EP1(ELIM)         Integer   Out      Event parameter #1
*   EP2(ELIM)         Integer   Out      Event parameter #2
*   ELIM              Integer   In       Length of arrays
*   ETOT              Integer   Out      # events
*   IERR              Integer   Out      Error code 0=no errors
*-----------------------------------------------------------------------
      SUBROUTINE CLEF (A, WIDTH, HEIGHT, TRNPOS, RAXBRT,
     $               ET, EK, EP1, EP2, ELIM, ETOT, IERR)
       IMPLICIT NONE
       INTEGER WIDTH, HEIGHT, ELIM
       INTEGER ET(ELIM), EK(ELIM), EP1(ELIM), EP2(ELIM), ETOT, IERR
       REAL A(WIDTH,HEIGHT), TRNPOS, RAXBRT

*       Constants - must agree with elsewhere, including program spectr
       INTEGER DUR, KLIM, PXPS, VERT
       REAL ADPARM, ATOL, BF, BOGUS, FMARK, FMNNLO, FMNNHI,
     $  FSPAC, MDPARM, PARM, PBR, PI, PI2,
     $  PIP, RE, RESPAR, SIG, USPPX, VCRIT

       PARAMETER (
     $            ADPARM = 1.,          ! difference to split notes
     $            MDPARM = 3.,          ! max difference within a note
     $            BOGUS = -999.,        ! invalid note number
     $            DUR = 5 * 60,         ! five minutes
     $            FMNNLO = 52.5,        ! note # at least frequency
     $            FMNNHI = 84.5,        ! note # at greatest frequency
     $            PBR = 200.,           ! pitch bend range in cents
     $            PI  = 3.141592654,
     $            PI2 = 9.869604401,
     $            PIP = 0.035,          ! ignore events briefer than PIP [s]
     $            PXPS = 240,           ! pixels per second
     $            RE = 1.E-6,           ! relative error for root-finding
     $            RESPAR = 1. / 1024.,  ! resolution parameter
     $            SIG = 2.0,            ! for Lanczos interpolation
     $            VERT = 640,           ! height of image
     $            VCRIT = 0.66)         ! threshhold variance
       PARAMETER (BF = 8192. * 100. / PBR,  ! pitch bend units per semitone
     $            FSPAC = (FMNNHI - FMNNLO) / VERT,   ! semitones per pixel
     $            KLIM = PXPS * DUR,
     $            PARM = 1. + 1. / RESPAR,
     $            USPPX = 1 000 000. / PXPS)     ! microseconds per pixel
       PARAMETER (ATOL = 0.1 / (BF * FSPAC),     ! absolute error tolerance
     $            FMARK = FMNNLO - 0.5 * FSPAC)  ! frequency of zeroth pixel

*          Local variables
       INTEGER
     $         BACK(0:1,KLIM),  ! back pointers
     $         BOUND(KLIM),     ! index of pixel of a note onset
     $         EX, EXOLD,       ! expression values
     $         I, J, L,         ! counters
     $         JMAX,            ! index of maximum brightness
     $         K,               ! # segments
     $         LENGTH(KLIM),    ! workspace for GLOM
     $         N, NOLD,         ! MIDI Note Number as integer
     $         PB, PBOLD,       ! pitch bend values
     $         RIND(KLIM),      ! workspace for GLOM
     $         SMOORD,          ! order of smoothing
     $         T,               ! time in microseconds
     $         V                ! midi velocity

       REAL
     $      AMPL(VERT),     ! pixel intensity in a column
     $      AVE(KLIM),      ! average pitch each note
     $      B, C,           ! bounds of search interval
     $      BASE,           ! controls diminuendo
     $      BRT,            ! brightness
     $      COREX,          ! factor in expression formula
     $      DW(KLIM),       ! workspace for GLOM
     $      FMEAN, FVAR,    ! mean and variance of frequency
     $      FREQ(VERT),     ! Midi note number in a column
     $      GREAT(KLIM),    ! workspace for GLOM
     $      LEAST(KLIM)     ! workspace for GLOM

       REAL POS,            ! fractional index of max brightness
     $      RARG(4),        ! arguments for kernel function
     $      SRRP,           ! square root of resolution parameter
     $      SUMX(KLIM),     ! workspace for GLOM
     $      SA, SAP, SAP2,  ! sums of amplitude
     $      TRANS,          ! transition penalty
     $      WORK(VERT),     ! workspace for QWM
     $      X(KLIM),        ! midi note number each time step
     $      Y(KLIM),        ! total intensity, each time step
     $      YM              ! maximum value of amplitude

       LOGICAL SOUND        ! note is sounding at a pixel

*             External functions
       REAL KERNIN
       EXTERNAL KERNIN
       LOGICAL SAFDIV

*-----------------------------------------------------------------------
*             Begin.
*-----------------------------------------------------------------------
       IF (WIDTH < 1 .OR. HEIGHT < 1) THEN
         IERR = 1
         RETURN
       ELSE IF (HEIGHT .NE. VERT) THEN
         IERR = 2
         RETURN
       END IF
       DO J = 1, VERT
         FREQ(J) = FMARK + J * FSPAC + TRNPOS
       END DO

*-----------------------------------------------------------------------
*             Pitch and amplitude each time step.
*-----------------------------------------------------------------------
       DO 10 I = 1, WIDTH
         K = 0
         JMAX = 0
         YM = 0.
         DO J = 1, VERT
           AMPL(J) = A(I,J) / RAXBRT
           IF (A(I,J) > 0.5) K = K + 1
           IF (AMPL(J) > YM) THEN
             JMAX = J
             YM = AMPL(J)
           END IF
         END DO

         IF (K .EQ. 0) THEN       ! silence
           X(I) = BOGUS
           Y(I) = 0.
           GO TO 10
         ELSE IF (K .EQ. 1) THEN  ! one pixel
           X(I) = FREQ(JMAX)
           Y(I) = YM
           GO TO 10
         END IF

         SA = 0.                  ! weighted variance
         SAP = 0.
         DO J = 1, VERT
           SA = SA + AMPL(J)
           SAP = SAP + AMPL(J) * FLOAT(J)
         END DO
         FMEAN = SAP / SA
         SAP2 = 0.
         DO J = 1, VERT
           SAP2 = SAP2 + AMPL(J) * (FLOAT(J) - FMEAN)**2
         END DO
         FVAR = SAP2 / SA

*-----------------------------------------------------------------------
*             De-interpolate in units of pixels and linear light.
*-----------------------------------------------------------------------
         IF (FVAR < VCRIT) THEN   ! attempt if variance is small
           YM = 0.
           JMAX = 2
           DO J = 2, VERT
             SA = AMPL(J-1) + AMPL(J)
             IF (SA > YM) THEN
               YM = SA
               JMAX = J
             END IF
           END DO

           RARG(1) = FLOAT(JMAX-1)
           RARG(2) = FLOAT(JMAX)
           RARG(3) = RESPAR * (PARM**AMPL(JMAX-1) - 1.)
           RARG(4) = RESPAR * (PARM**AMPL(JMAX) - 1.)
           B = RARG(1) + ATOL
           C = RARG(2) - ATOL
           POS = B
           CALL ZERORA (KERNIN, POS, C, RE, ATOL, RARG, IERR)
           IF (IERR .EQ. 1 .OR. IERR .EQ. 2) THEN
             IF (RARG(3) > RARG(4)) THEN
               C = RARG(1)
               YM = RARG(3)
             ELSE
               C = RARG(2)
               YM = RARG(4)
             END IF
             B = PI2 * (C - POS)**2
             C = SIG * SIN(PI * (C - POS)) * SIN(PI * (C - POS) / SIG)
             IF (SAFDIV(B, C)) THEN
               BRT = B / C
             ELSE
               BRT = 1.
             END IF
             BRT = BRT * YM
             X(I) = FMARK + POS * FSPAC + TRNPOS  ! pixel -> MIDI note
             Y(I) = LOG(1. + BRT / RESPAR) / LOG(PARM)
             GO TO 10
           END IF
         END IF

*-----------------------------------------------------------------------
*             Fallback routine in units of Midi note
*-----------------------------------------------------------------------
         DO J = 1, VERT
           AMPL(J) = RESPAR * (PARM**AMPL(J) - 1.)
         END DO
         CALL QWM (FREQ, AMPL, VERT, POS, BRT, WORK, IERR)
         IF (IERR .NE. 0) THEN
           PRINT *, 'qwm returns #', IERR
           RETURN
         END IF
         BRT = LOG(1. + BRT / RESPAR) / LOG(PARM)
         BRT = MIN(MAX(0., BRT), 1.)
         X(I) = POS
         Y(I) = BRT
  10   CONTINUE

*-----------------------------------------------------------------------
* Choose sounding or silent by dynamic programming.
*-----------------------------------------------------------------------
       DO I = 1, WIDTH                           ! Limit to 0...1
         Y(I) = MIN(MAX(0., Y(I)), 1.)
         GREAT(I) = (Y(I) * RAXBRT + 1.) / (RAXBRT + 1.)  ! can't take log(0)
       END DO
       SRRP = SQRT(RESPAR)
       TRANS = LOG(RAXBRT) * FLOAT(PXPS) * PIP
       CALL DYRB (GREAT, WIDTH, SRRP, TRANS, BOUND, BACK, IERR)
       IF (IERR .NE. 0) PRINT *, 'dyrb returns #', IERR

       DO K = 1, WIDTH
         IF (BOUND(K) .EQ. 0) THEN
           X(K) = BOGUS                ! force silence
           Y(K) = 0.
         END IF
       END DO

       DO K = WIDTH-1, 2, -1
         IF (BOUND(K) .EQ. 1 .AND. X(K) < 0.) THEN
           X(K) = X(K+1)               ! fill in
           Y(K) = Y(K+1)
         END IF
       END DO

       DO K = 2, WIDTH-1, +1
         IF (BOUND(K) .EQ. 1 .AND. X(K) < 0.) THEN
           X(K) = X(K-1)
           Y(K) = Y(K-1)
         END IF
       END DO

*-----------------------------------------------------------------------
*             Merge contiguous frequencies into notes.
*-----------------------------------------------------------------------
       SOUND = .FALSE.                           ! mark note beginnings
       K = 0
       DO J = 1, WIDTH
         IF (.NOT. SOUND .AND. X(J) .GE. 0.) THEN
           SOUND = .TRUE.
           K = K + 1
           BOUND(K) = J
         ELSE IF (SOUND .AND. (X(J) < 0. .OR. J .EQ. WIDTH)) THEN
           SOUND = .FALSE.
           LENGTH(K) = J - BOUND(K)
         END IF
       END DO

*             Smoothing pass
       SMOORD = INT(PXPS * PIP)
       DO J = 1, K
         CALL SMOOTH (X(BOUND(J)), LENGTH(J), SMOORD, LEAST, GREAT, AVE)
         CALL SMOOTH (Y(BOUND(J)), LENGTH(J), SMOORD, LEAST, GREAT, AVE)
       END DO

       K = 1
       CALL GLOM (X, WIDTH, ADPARM, MDPARM, K, BOUND, LENGTH, AVE,
     $        LEAST, GREAT, SUMX, DW, BACK(0,1), BACK(0,KLIM/2+1), RIND)

*             max velocity per note in SUMX
       I = 0
       DO J = 1, K
         AVE(J) = 0.5 * (LEAST(J) + GREAT(J))  ! split the difference
         SUMX(J) = 0.
         L = LENGTH(J)
         IF (L .LE. 0) PRINT *, 'clef: warning, zero length ', L
         DO T = 1, L
           I = I + 1
           IF (I > WIDTH) PRINT *, 'clef: out of bounds at velocity'
           IF (Y(I) > SUMX(J)) SUMX(J) = Y(I)
         END DO
         IF (AVE(J) < 0.) SUMX(J) = BOGUS
       END DO

*-----------------------------------------------------------------------
*             Make events
*-----------------------------------------------------------------------
       BASE = SQRT(PARM)
       COREX = SQRT(1. + RESPAR)
       EX = 127          ! default max for expression
       PB = 8192         ! default center for pitch bend
       J = 0             ! not yet to first segment
       L = 0             ! zero offset
       N = NINT(BOGUS)
       NOLD = N
       V = 0
       DO I = 1, WIDTH
         NOLD = N
         EXOLD = EX
         PBOLD = PB
         T = NINT(I * USPPX)
         IF (I .GE. L) THEN  ! advance to next segment?
           IF (J .GE. K) THEN
             IERR = 3
             RETURN
           END IF
           J = J + 1
           L = BOUND(J) + LENGTH(J)
           N = NINT(AVE(J))                      ! per-segment pitch
           N = MIN(N, 127)

*-----------------------------------------------------------------------
* Quote from the MIDI specification:
*
*     "Preferably, application of velocity to volume should be an exponential
*      function. This is the suggested default action;"
*
* Program SPECTR sets the image intensity to the logarithm of sound amplitude.
* Here velocity is set directly proportional to image intensity.  Thus,
* an exponential interpretation of velocity restores the original volume.
*-----------------------------------------------------------------------
           IF (SUMX(J) > 0.) THEN
             V = NINT(SUMX(J) * 127.)  ! per-segment velocity
             V = MIN(MAX(1, V), 127)
           ELSE
             V = 0
           END IF
         END IF

*-----------------------------------------------------------------------
*             Note Off
*-----------------------------------------------------------------------
         IF (NOLD .GE. 0 .AND. N .NE. NOLD) THEN
           IF (ETOT .GE. ELIM) RETURN
           ETOT = ETOT + 1
           ET(ETOT) = T
           EK(ETOT) = 0
           EP1(ETOT) = NOLD
           EP2(ETOT) = 0
         END IF

         IF (V > 0) THEN                         ! sounding ?
           IF (N < 0) PRINT *, 'clef: positive v, negative n !!'

*-----------------------------------------------------------------------
*             Expression
* The MIDI spec describes the possibility of using volume only, or volume
* and expression.  The equation for volume only is:
*
*      L = 40 log (V/127)
*
* where L is amplitude in decibels, and V is the MIDI volume.
* Program SCRIBE uses expression only.  In program SPECTR, intensity was:
*
*            log(A + p A_max) - log(p A_max)
*      Y = ------------------------------------ R
*           log(A_max + p A_max) - log(p A_max)
*
* where p is on the order of the precision, and R is the maximum brightness.
* A decibel of amplitude, defined so as to avoid taking a logarithm of zero:
*
*                    A + p A_max
*      dBA = 20 log(-------------)
*                       A_max
*
* After some algebra, it appears:
*
*      E = 127 (1 + p)^(1/2) (1/p + 1)^( (1/2) (Y/R - V/127))
*
*-----------------------------------------------------------------------
           EX = NINT(127. * COREX * BASE**(Y(I) - FLOAT(V) / 127.) )
           EX = MIN(MAX(0, EX), 127)
           IF (EX .NE. EXOLD) THEN
             IF (ETOT .GE. ELIM) RETURN
             ETOT = ETOT + 1
             ET(ETOT) = T
             EK(ETOT) = 12
             EP1(ETOT) = EX
           END IF

*-----------------------------------------------------------------------
*             Pitch bend
*-----------------------------------------------------------------------
           PB = 8192 + NINT( (X(I) - FLOAT(N)) * BF )
           PB = MIN(MAX(0, PB), 16383)
           IF (PB .NE. PBOLD) THEN
             IF (ETOT .GE. ELIM) RETURN
             ETOT = ETOT + 1
             ET(ETOT) = T
             EK(ETOT) = 57
             EP1(ETOT) = PB
           END IF
         END IF

*-----------------------------------------------------------------------
*             Note On
*-----------------------------------------------------------------------
         IF (N .GE. 0 .AND. N .NE. NOLD) THEN
           IF (ETOT .GE. ELIM) RETURN
           ETOT = ETOT + 1
           ET(ETOT) = T
           EK(ETOT) = 1
           EP1(ETOT) = N
           EP2(ETOT) = V
         END IF
       END DO  ! next time step

*-----------------------------------------------------------------------
*             Turn off any note left sounding
*-----------------------------------------------------------------------
       IF (NOLD .GE. 0) THEN
         IF (ETOT .GE. ELIM) RETURN
         ETOT = ETOT + 1
         ET(ETOT) = NINT((WIDTH+1) * USPPX)
         EK(ETOT) = 0
         EP1(ETOT) = NOLD
         EP2(ETOT) = 0
       END IF
       RETURN
      END  ! of clef


*-----------------------------------------------------------------------
* kernin - invert the Lanczos kernel function
*
*___Name______Type___In/Out___Description_______________________________
*   X         Real   In       A position in pixels of the pre-interpolated data
*   RARG(4)   Real   In       RARG(1) Position X1 after interpolation
*                             RARG(2) Position X2 after interpolation
*                             RARG(3) Intensity Y1 at position X1
*                             RARG(4) Intensity Y2 at position X2
*   KERNIN    Real   Out      Residual
*-----------------------------------------------------------------------
      FUNCTION KERNIN (X, RARG)
       IMPLICIT NONE
       REAL X, RARG(4), KERNIN

       REAL PI, SIG
       PARAMETER (PI = 3.141592654,
     $            SIG = 2.0)
       REAL A, B, DX, Q1, Q2
       LOGICAL SAFDIV

*-----------------------------------------------------------------------
*             Begin
*-----------------------------------------------------------------------
       DX = RARG(1) - X
       A = DX**2
       B = SIN(PI * DX) * SIN(PI * DX / SIG)
       IF (SAFDIV(A, B)) THEN
         Q1 = A / B
       ELSE
         Q1 = 1.
       END IF
       Q1 = Q1 * RARG(3)
       DX = RARG(2) - X
       A = DX**2
       B = SIN(PI * DX) * SIN(PI * DX / SIG)
       IF (SAFDIV(A, B)) THEN
         Q2 = A / B
       ELSE
         Q2 = 1.
       END IF
       Q2 = Q2 * RARG(4)
       KERNIN = Q1 - Q2

       RETURN
      END  ! of kernin


*-----------------------------------------------------------------------
* qwm - quadratic weighted mean.  X is assumed sorted and equally spaced.
*
*___Name______Type______In/Out____Description___________________________
*   X(N)      Real      In        Data array
*   Y(N)      Real      In        Importance weights
*   N         Integer   In        Length of arrays
*   AVE       Real      Out       Average
*   BRT       Real      Out       Brightness
*   WORK(N)   Real      Neither   Cumulative Y
*   IERR      Integer   Out       Error code
*-----------------------------------------------------------------------
      SUBROUTINE QWM (X, Y, N, AVE, BRT, WORK, IERR)
       IMPLICIT NONE
       INTEGER N, IERR
       REAL X(N), Y(N), AVE, BRT, WORK(N)

*             Constant
       REAL PI
       PARAMETER (PI = 3.141592654)

*             Local variables
       INTEGER I, J, K
       REAL A, B,
     $      G2,              ! Cauchy distibution parameter
     $      H,               ! temp scalar
     $      Q,               ! a fraction of total weight
     $      WJM1,            ! W(J-1)
     $      XQ(3),           ! quartile abscissas
     $      XSPAC            ! spacing of X values
       DOUBLE PRECISION SWXY, SWY, SY

       LOGICAL DAFDIV, SAFDIV        ! external functions

*-----------------------------------------------------------------------
*             Begin.
*-----------------------------------------------------------------------
       IF (N < 1) THEN
         IERR = -1
         RETURN
       END IF
       IERR = 0
       IF (N .EQ. 1) THEN
         AVE = X(1)
         BRT = Y(1)
         RETURN
       END IF

       SY = 0.D0                  ! sum of intensities
       DO J = 1, N
         IF (Y(J) < 0.) Y(J) = 0.
         SY = SY + Y(J)
         WORK(J) = SNGL(SY)
       END DO
       IF (SY .LE. 0.) THEN
         IERR = 1
         RETURN
       END IF

*             Weighted quartiles
       XSPAC = (X(N) - X(1)) / FLOAT(N - 1)
       DO K = 1, 3
         Q = FLOAT(K) / 4. * SNGL(SY)
         CALL RBSGT (WORK, N, Q, J)
         IF (J > 1) THEN
           WJM1 = WORK(J-1)
         ELSE
           WJM1 = 0.
         END IF
         XQ(K) = X(J) + ((Q - WJM1) / Y(J) - 0.5) * XSPAC
       END DO

*-----------------------------------------------------------------------
*             Quadratic weighted mean
*-----------------------------------------------------------------------
       AVE = XQ(2)
       G2 = 0.5 * (XQ(3) - XQ(1))                ! semi-interquartile range
       G2 = G2**2

       SWY = 0.D0
       SWXY = 0.D0
       DO J = 1, N
         H = (X(J) - AVE)**2 + G2
         IF (SAFDIV(Y(J), H)) THEN
           H = Y(J) / H
         ELSE
           IERR = 2
           RETURN
         END IF
         SWY = SWY + H
         SWXY = SWXY + H * (X(J) - AVE)
       END DO

       IF (DAFDIV(SWXY, SWY)) THEN
         AVE = AVE + SNGL(SWXY / SWY)            ! weighted mean
       ELSE
         IERR = 3
         RETURN
       END IF

*-----------------------------------------------------------------------
*             Maximum weight within 2g of the median
*-----------------------------------------------------------------------
       A = XQ(2) - (XQ(3) - XQ(1))
       B = XQ(2) + (XQ(3) - XQ(1))
       CALL RBSGT (X, N, A, J)
       CALL RBSLT (X, N, B, K)
       BRT = 0.
       DO I = J, K
         IF (Y(I) > BRT) BRT = Y(I)
       END DO

       RETURN
      END  ! of qwm


*-----------------------------------------------------------------------
* rpgm - read a portable gray map
*
*___Name______Type____________In/Out___Description______________________
*   IFIL      Character*(*)   In       Filename
*   A(LLIM)   Real            In       Bitmap array
*   LLIM      Integer         In       Max length of A
*   HORZ      Integer         Out      Width of image
*   VERT      Integer         Out      Height of image
*   MAXBRT    Integer         Out      Maximum intensity
*   IERR      Integer         Out      Error code
*-----------------------------------------------------------------------
      SUBROUTINE RPGM (IFIL, A, LLIM, HORZ, VERT, MAXBRT, IERR)
       IMPLICIT NONE
       INTEGER LLIM, HORZ, VERT, MAXBRT, IERR
       REAL A(LLIM)
       CHARACTER*(*) IFIL

*             Local variables
       INTEGER*1 B
       INTEGER H, I, J, K,     ! counters
     $         HLIM,           ! 0 or 1 sets 1 byte or 2 bytes per pixel
     $         IOU,            ! i/o unit number
     $         L,              ! index into A
     $         N,              ! length of file in bytes, input
     $         PIXL,           ! total intensity for pixel
     $         POS             ! offset into file
       LOGICAL COMENT,         ! this is a comment line?
     $         PREVW           ! previous character was whitespace?
       CHARACTER*2 FTYP
       CHARACTER*32 DECLAR, HEADER

*            External function
       INTEGER IOUNIT

*-----------------------------------------------------------------------
*             Begin.
*-----------------------------------------------------------------------
       IERR = 0
       CALL LFIL (IFIL, N, IERR)
       IF (IERR .NE. 0) THEN
         PRINT *, 'trouble getting length of file: ', IFIL
         RETURN
       END IF
       IOU = IOUNIT ()
       OPEN (UNIT=IOU, FILE=IFIL, STATUS='OLD', ACCESS='DIRECT',
     $   RECL=1, FORM='UNFORMATTED', IOSTAT=IERR)
       IF (IERR .NE. 0) THEN
         PRINT *, 'trouble opening file: ', IFIL
         STOP
       END IF

*-----------------------------------------------------------------------
*             Interpret header.
*-----------------------------------------------------------------------
       POS = 0             ! position in file
       J = 0               ! # white space intervals
       L = 0               ! position in header
       HEADER = ' '
       PREVW = .FALSE.
       COMENT = .FALSE.
       DO 10 WHILE (J < 4)
         POS = POS + 1
         READ (UNIT=IOU, REC=POS, IOSTAT=IERR) B
         IF (IERR .NE. 0) PRINT *, 'rpgm: trouble in while loop!'
         IF (B .EQ. 35) COMENT = .TRUE.
         IF (B .EQ. 10) COMENT = .FALSE.
         IF (COMENT) GO TO 10
         IF ((B .GE. 9 .AND. B .LE. 13) .OR. B .EQ. 32) THEN
           IF (.NOT. PREVW) THEN
             J = J + 1
             L = L + 1
             HEADER(L:L) = ' '
             PREVW = .TRUE.
           END IF
         ELSE
           L = L + 1
           PREVW = .FALSE.
           HEADER(L:L) = CHAR(B)
         END IF
         IF (L .GE. 32) THEN
           PRINT *, 'header too long for PNM file format'
           IERR = 10
           RETURN
         END IF
  10   CONTINUE

       READ (UNIT=HEADER, FMT=*, IOSTAT=IERR)
     $      DECLAR, HORZ, VERT, MAXBRT
       IF (IERR .NE. 0) PRINT *, 'rpgm: trouble reading header!'
       FTYP = DECLAR(1:2)
       IF (      FTYP .NE. 'P2'
     $     .AND. FTYP .NE. 'P3'
     $     .AND. FTYP .NE. 'P5'
     $     .AND. FTYP .NE. 'P6') THEN
         IERR = 1
         RETURN
       ELSE IF (HORZ * VERT > LLIM) THEN
         PRINT *, 'ppm file too big!'
         IERR = 2
         RETURN
       ELSE IF (MAXBRT > 65535
     $          .AND. FTYP .NE. 'P2'
     $          .AND. FTYP .NE. 'P3') THEN
         PRINT *, 'maxbrt too big!'
         IERR = 3
         RETURN
       END IF

       IF (MAXBRT > 255) THEN
         HLIM = 1
       ELSE
         HLIM = 0
       END IF
       IF (FTYP .EQ. 'P3') GO TO 30
       IF (FTYP .EQ. 'P5') GO TO 50
       IF (FTYP .EQ. 'P6') GO TO 60

*-----------------------------------------------------------------------
*             Format P2
*-----------------------------------------------------------------------
       J = 0
       L = 0
       PREVW = .FALSE.
       COMENT = .FALSE.
       HEADER = ' '
       DO 20 I = POS+1, N
         READ (UNIT=IOU, REC=I, IOSTAT=IERR) B
         IF (IERR .NE. 0) THEN
           PRINT *, 'rpgm: trouble reading P2 data'
           RETURN
         END IF
         IF (B .EQ. 35) COMENT = .TRUE.
         IF (B .EQ. 10) COMENT = .FALSE.
         IF (COMENT) GO TO 20
         IF ((B .GE. 9 .AND. B .LE. 13) .OR. B .EQ. 32) THEN
           IF (.NOT. PREVW) THEN
             READ (UNIT=HEADER, FMT=*, IOSTAT=IERR) H
             IF (IERR .NE. 0) THEN
               PRINT *, 'rpgm: trouble internal I/O, format P2'
               RETURN
             END IF
             IF (H > MAXBRT) THEN
               PRINT *, 'rpgm: pixel too bright!'
               IERR = 4
               RETURN
             END IF
             L = L + 1
             IF (L > LLIM) THEN
               PRINT *, 'rpgm: out of memory!'
               IERR = 5
               RETURN
             END IF
             A(L) = FLOAT(H)
             PREVW = .TRUE.
             HEADER = ' '
             J = 0
           END IF
         ELSE
           J = J + 1
           PREVW = .FALSE.
           HEADER(J:J) = CHAR(B)
         END IF
  20   CONTINUE
       GO TO 90

*-----------------------------------------------------------------------
*             Format P3
*-----------------------------------------------------------------------
  30   J = 0
       K = 0
       L = 0
       PIXL = 0
       PREVW = .FALSE.
       COMENT = .FALSE.
       HEADER = ' '
       DO 40 I = POS+1, N
         READ (UNIT=IOU, REC=I, IOSTAT=IERR) B
         IF (IERR .NE. 0) THEN
           PRINT *, 'rpgm: trouble reading P3 data'
           RETURN
         END IF
         IF (B .EQ. 35) COMENT = .TRUE.
         IF (B .EQ. 10) COMENT = .FALSE.
         IF (COMENT) GO TO 40
         IF ((B .GE. 9 .AND. B .LE. 13) .OR. B .EQ. 32) THEN
           IF (.NOT. PREVW) THEN
             READ (UNIT=HEADER, FMT=*, IOSTAT=IERR) H
             IF (IERR .NE. 0) THEN
               PRINT *, 'rpgm: trouble internal I/O, format P3'
               RETURN
             END IF
             IF (H > MAXBRT) THEN
               PRINT *, 'rpgm: color too bright'
               IERR = 6
               RETURN
             END IF
             K = K + 1
             PIXL = PIXL + H
             IF (K .EQ. 3) THEN
               L = L + 1
               IF (L > LLIM) THEN
                 PRINT *, 'rpgm:  out of memory!'
                 IERR = 7
                 RETURN
               END IF
               A(L) = PIXL
               K = 0
               PIXL = 0
             END IF
             PREVW = .TRUE.
             HEADER = ' '
             J = 0
           END IF
         ELSE
           J = J + 1
           PREVW = .FALSE.
           HEADER(J:J) = CHAR(B)
         END IF
  40   CONTINUE
       MAXBRT = MAXBRT * 3
       GO TO 90

*-----------------------------------------------------------------------
*             Format P5
*-----------------------------------------------------------------------
  50   L = 0
       N = 0
       IF (VERT * HORZ > LLIM) THEN
         PRINT *,'rpgm: image to big for allocated memory!'
         IERR = 8
         RETURN
       END IF
       DO J = 1, VERT
         DO I = 1, HORZ
           DO H = HLIM, 0, -1
             POS = POS + 1
             READ (UNIT=IOU, REC=POS, IOSTAT=IERR) B
             IF (IERR .NE. 0) PRINT *, 'rpgm: reading P5 at ', I,J,H
             CALL BYTGET (B, N, H)
           END DO
           L = L + 1
           A(L) = FLOAT(N)
         END DO
       END DO
       GO TO 90


*-----------------------------------------------------------------------
*             Format P6
*-----------------------------------------------------------------------
  60   L = 0
       N = 0
       IF (VERT * HORZ > LLIM) THEN
         PRINT *,'rpgm: image to big for allocated memory!'
         IERR = 9
         RETURN
       END IF
       DO J = 1, VERT
         DO I = 1, HORZ
           PIXL = 0
           DO K = 1, 3
             DO H = HLIM, 0, -1
               POS = POS + 1
               READ (UNIT=IOU, REC=POS, IOSTAT=IERR) B
               IF (IERR .NE. 0) PRINT *, 'rpgm: reading P6 ', I,J,K,H
               CALL BYTGET (B, N, H)
             END DO
             PIXL = PIXL + N
           END DO
           L = L + 1
           A(L) = FLOAT(PIXL)
         END DO
       END DO
       MAXBRT = MAXBRT * 3

*-----------------------------------------------------------------------
*             Done.
*-----------------------------------------------------------------------
  90   CLOSE (UNIT=IOU, IOSTAT=IERR)
       IF (IERR .NE. 0) PRINT *, 'trouble closing output file!'
       RETURN
      END  ! of rpgm


*-----------------------------------------------------------------------
* wmidi - write a Standard MIDI file.  This subroutine was based on
* the much longer version appearing in project MIXMIDI, which supports
* the entire MIDI 1.0 specification.
*
*                        SUPPORTED EVENT KINDS
*
*___Kind___Description______Par._1________Par._2______________MIDI_#____
*    0     Note Off         Note #        Velocity            8
*    1     Note On          Note #        Velocity            9
*   12     Expression       7 bits                            11-11
*   55     Patch change     Patch #                           12
*   57     Pitch Bend       14 bits                           14
*   61     Text             Offset        Length              15-15-1
*   62     Copyright        Offset        Length              15-15-2
*   63     Sequence Name    Offset        Length              15-15-3
*   76     Bend Range       Cents                             RPN 0
*
*                              I/O LIST
*
*___Name________Type____________In/Out____Description___________________
*   FNAME       Character*(*)   In        Output filename
*   TEXT        Character*(*)   In        Text of name, copyright, etc.
*   ETOT        Integer         In        Count of all events
*   ET(ETOT)    Integer         In        Event times [microseconds]
*   EK(ETOT)    Integer         In        Event kind
*   EP1(ETOT)   Integer         In        Event parameter 1
*   EP2(ETOT)   Integer         In        Event parameter 2
*   IFAULT      Integer         Out       Error code 0=no errors
*-----------------------------------------------------------------------
      SUBROUTINE WMIDI (FNAME, TEXT, ETOT, ET, EK, EP1, EP2, IFAULT)
       IMPLICIT NONE
       INTEGER ETOT
       INTEGER ET(ETOT), EK(ETOT), EP1(ETOT), EP2(ETOT), IFAULT
       CHARACTER*(*) FNAME, TEXT

*           Local variables
       INTEGER*1
     $           B, C, D,             ! single bytes
     $           VARB(4)              ! a delta-time

       INTEGER
     $           G,                    ! count events
     $           H,                    ! misc. counter
     $           I,                    ! file position
     $           IOU,                  ! a file I/O Unit
     $           KIND,                 ! what kind of event
     $           L,                    ! position to write track length at
     $           M,                    ! temporary integer
     $           N,                    ! length of track data
     $           P,                    ! number of data bytes in an event
     $           S, SPREV              ! ticks elapsed so far

       REAL
     $      TPUS                       ! ticks per microsecond

*            External functions
       INTEGER
     $           IOUNIT                ! get a valid I/O Unit

*-----------------------------------------------------------------------
* Constants.  Timebase of 120 translates 240 pix/sec to ticks at 1:1.
*-----------------------------------------------------------------------
       INTEGER*1
     $           MTHD(8), MTRK(4),   ! file and track headings
     $           NULL
       INTEGER BPM, CHANEL, TIMDIV
       PARAMETER (BPM = 120,         ! always 120 beats per minute
     $            CHANEL = 0,        ! channel number always 0
     $            TIMDIV = 120)      ! timing resolution in ticks per Q.N.
       SAVE MTHD, MTRK, NULL
       DATA MTHD / 77, 84, 104, 100, 0, 0, 0, 6 /
       DATA MTRK / 77, 84, 114, 107 /
       DATA NULL / 0 /
       DATA B, C, D / 0, 0, 0 /  ! avoid compiler warnings

*-----------------------------------------------------------------------
*        Begin.
*-----------------------------------------------------------------------
       IFAULT = 0
       TPUS = TIMDIV * BPM / 60. / 1 000 000.    ! Convert times to ticks
       DO I = 1, ETOT
         ET(I) = NINT(FLOAT(ET(I)) * TPUS)
       END DO

*-----------------------------------------------------------------------
*               File I/O
*-----------------------------------------------------------------------
       IOU = IOUNIT()
       OPEN (UNIT=IOU, FILE=FNAME, ACCESS='DIRECT', RECL=1,
     $     FORM='UNFORMATTED', STATUS='NEW', IOSTAT=IFAULT)
         IF (IFAULT .NE. 0) PRINT *, 'trouble opening file ', FNAME

*                     Write the header
       DO I = 1, 8
         WRITE (UNIT=IOU, REC=I, ERR=90) MTHD(I)
       END DO

*                     File format / type 0
       B = 0
       C = 0
       WRITE (UNIT=IOU, REC=9, ERR=90) B
       WRITE (UNIT=IOU, REC=10, ERR=90) C

*                      Number of tracks
       B = 0
       C = 1
       WRITE (UNIT=IOU, REC=11, ERR=90) B
       WRITE (UNIT=IOU, REC=12, ERR=90) C

*                       Time division
       CALL BYTPUT (B, TIMDIV, 1)
       CALL BYTPUT (C, TIMDIV, 0)
       WRITE (UNIT=IOU, REC=13, ERR=90) B
       WRITE (UNIT=IOU, REC=14, ERR=90) C

       I = 14      ! at end of file header
       G = 0       ! no events processed

       DO H = 1, 4             ! write 'MTrk'
         I = I + 1
         WRITE (UNIT=IOU, REC=I, ERR=90) MTRK(H)
       END DO

*           Remember file position for track length
       L = I
       I = I + 4              ! skip ahead

*-----------------------------------------------------------------------
*            Process events
*-----------------------------------------------------------------------
       S = 0                      ! at beginning
       DO G = 1, ETOT          ! for all events in track
         KIND = EK(G)

*-----------------------------------------------------------------------
*                 Write delta-time
*-----------------------------------------------------------------------
         SPREV = S
         S = ET(G)
         M = S - SPREV
         IF (M < 0) PRINT *, 'wmidi: negative delta-time!'
         CALL WVARB (M, VARB, P)           ! convert to variable bytes
         IF (P .EQ. -1) PRINT *, 'cant convert delta time of ', M
         DO H = 1, P
           I = I + 1
           WRITE (UNIT=IOU, REC=I, ERR=90) VARB(H)
         END DO

*-----------------------------------------------------------------------
*               NoteOff - zero velocity NoteOn
*-----------------------------------------------------------------------
         IF (KIND .EQ. 0) THEN
           M = IOR(144, CHANEL)
           CALL BYTPUT (B, M, 0)                ! status byte
           CALL BYTPUT (C, EP1(G), 0)           ! note #
           WRITE (UNIT=IOU, REC=I+1, ERR=90) B
           WRITE (UNIT=IOU, REC=I+2, ERR=90) C
           WRITE (UNIT=IOU, REC=I+3, ERR=90) NULL
           I = I + 3

*-----------------------------------------------------------------------
*             NoteOn
*-----------------------------------------------------------------------
         ELSE IF (KIND .EQ. 1) THEN
           M = IOR(144, CHANEL)
           CALL BYTPUT (B, M, 0)                ! status byte
           CALL BYTPUT (C, EP1(G), 0)           ! note #
           CALL BYTPUT (D, EP2(G), 0)           ! velocity
           WRITE (UNIT=IOU, REC=I+1, ERR=90) B
           WRITE (UNIT=IOU, REC=I+2, ERR=90) C
           WRITE (UNIT=IOU, REC=I+3, ERR=90) D
           I = I + 3

*-----------------------------------------------------------------------
*             Expression
*-----------------------------------------------------------------------
         ELSE IF (KIND .EQ. 12) THEN              ! expression
           M = IOR(176, CHANEL)
           CALL BYTPUT (B, M, 0)
           C = 11
           CALL BYTPUT (D, EP1(G), 0)
           IF (D < 0) THEN
             PRINT *, 'WMIDI: negative expression!'
             IFAULT = 3
           END IF
           WRITE (UNIT=IOU, REC=I+1, ERR=90) B
           WRITE (UNIT=IOU, REC=I+2, ERR=90) C
           WRITE (UNIT=IOU, REC=I+3, ERR=90) D
           I = I + 3

*-----------------------------------------------------------------------
*             Patch change
*-----------------------------------------------------------------------
           ELSE IF (KIND .EQ. 55) THEN
             M = IOR(192, CHANEL)
             CALL BYTPUT (B, M, 0)
             CALL BYTPUT (C, EP1(G), 0)
             WRITE (UNIT=IOU, REC=I+1, ERR=90) B
             WRITE (UNIT=IOU, REC=I+2, ERR=90) C
             I = I + 2

*-----------------------------------------------------------------------
*             Pitch bend
*-----------------------------------------------------------------------
         ELSE IF (KIND .EQ. 57) THEN
           M = IOR(224, CHANEL)
           CALL BYTPUT (B, M, 0)
           M = IBITS(EP1(G), 0, 7)
           CALL BYTPUT (C, M, 0)
           M = IBITS(EP1(G), 7, 7)
           CALL BYTPUT (D, M, 0)
           WRITE (UNIT=IOU, REC=I+1, ERR=90) B
           WRITE (UNIT=IOU, REC=I+2, ERR=90) C
           WRITE (UNIT=IOU, REC=I+3, ERR=90) D
           I = I + 3

*-----------------------------------------------------------------------
*             Text-type events
*-----------------------------------------------------------------------
         ELSE IF (KIND .GE. 61 .AND. KIND .LE. 63) THEN
           B = -1             ! meta
           M = KIND - 60
           CALL BYTPUT (C, M, 0)
           WRITE (UNIT=IOU, REC=I+1, ERR=90) B
           WRITE (UNIT=IOU, REC=I+2, ERR=90) C
           I = I + 2
           M = EP2(G)
           CALL WVARB(M, VARB, P)
           IF (P .EQ. -1) PRINT *, 'cant convert text length ', m
           DO H = 1, P                               ! write length
             I = I + 1
             WRITE (UNIT=IOU, REC=I, ERR=90) VARB(H)
           END DO
           P = EP1(G)                              ! pointer
           DO H = P, P+M-1                         ! write string
             I = I + 1
             B = ICHAR(TEXT(H:H))
             WRITE (UNIT=IOU, REC=I, ERR=90) B
           END DO

*-----------------------------------------------------------------------
*             Registered parameter - pitch bend range
*-----------------------------------------------------------------------
         ELSE IF (KIND .EQ. 76) THEN
           M = IOR(176, CHANEL)
           CALL BYTPUT (B, M, 0)
           C = 101                            ! coarse adjust
           D = 0
           WRITE (UNIT=IOU, REC=I+1, ERR=90) B
           WRITE (UNIT=IOU, REC=I+2, ERR=90) C
           WRITE (UNIT=IOU, REC=I+3, ERR=90) D
           C = 100                            ! fine adjust
           D = 0
           WRITE (UNIT=IOU, REC=I+4, ERR=90) NULL   ! delta-time
           WRITE (UNIT=IOU, REC=I+5, ERR=90) B
           WRITE (UNIT=IOU, REC=I+6, ERR=90) C
           WRITE (UNIT=IOU, REC=I+7, ERR=90) D
           C = 6                                    ! slider
           M = EP1(G) / 100
           CALL BYTPUT (D, M, 0)
           WRITE (UNIT=IOU, REC=I+8, ERR=90) NULL
           WRITE (UNIT=IOU, REC=I+9, ERR=90) B
           WRITE (UNIT=IOU, REC=I+10, ERR=90) C
           WRITE (UNIT=IOU, REC=I+11, ERR=90) D
           C = 38
           M = MOD(EP1(G), 100)
           CALL BYTPUT (D, M, 0)
           WRITE (UNIT=IOU, REC=I+12, ERR=90) NULL
           WRITE (UNIT=IOU, REC=I+13, ERR=90) B
           WRITE (UNIT=IOU, REC=I+14, ERR=90) C
           WRITE (UNIT=IOU, REC=I+15, ERR=90) D
           I = I + 15

         ELSE
           PRINT *, 'unrecognized event kind ', KIND
           IFAULT = -7
         END IF  ! kind of event
       END DO  ! track events

*-----------------------------------------------------------------------
*             End of track
*-----------------------------------------------------------------------
       B = -1
       C = 47
       WRITE (UNIT=IOU, REC=I+1, ERR=90) NULL
       WRITE (UNIT=IOU, REC=I+2, ERR=90) B
       WRITE (UNIT=IOU, REC=I+3, ERR=90) C
       WRITE (UNIT=IOU, REC=I+4, ERR=90) NULL
       I = I + 4

*                 Length of track events
       N = I - L - 4
       CALL BYTPUT (B, N, 3)
       WRITE (UNIT=IOU, REC=L+1, ERR=90) B
       CALL BYTPUT (B, N, 2)
       WRITE (UNIT=IOU, REC=L+2, ERR=90) B
       CALL BYTPUT (B, N, 1)
       WRITE (UNIT=IOU, REC=L+3, ERR=90) B
       CALL BYTPUT (B, N, 0)
       WRITE (UNIT=IOU, REC=L+4, ERR=90) B

*                     done
       CLOSE (UNIT=IOU, IOSTAT=IFAULT)
       IF (IFAULT .NE. 0) PRINT *, 'WMIDI: trouble closing file'
       RETURN      ! success

*-----------------------------------------------------------------------
  90   CONTINUE         ! fatal file I/O errors
*-----------------------------------------------------------------------
       PRINT *, 'trouble writing to file ', FNAME, ' position ', I
       STOP 'unrecoverable error'

      END  ! of wmidi


*-----------------------------------------------------------------------
* Convert an integer to MIDI variable-byte format
*-----------------------------------------------------------------------
      SUBROUTINE WVARB (N, VARB, P)
       IMPLICIT NONE
       INTEGER N, P
       INTEGER*1 VARB(4)

*           1 byte size
       IF (N .LE. 127) THEN
         VARB(1) = IAND(N, 127)
         P = 1
         RETURN

*           2 byte size
       ELSE IF (N .LE. 16383) THEN
         VARB(1) = IBSET(ISHFT(IAND(N, 16256), -7), 7)
         VARB(2) = IAND(N, 127)
         P = 2
         RETURN

*          3 byte size
       ELSE IF (N .LE. 2097151) THEN
         VARB(1) = IBSET(ISHFT(IAND(N, 2080768), -14), 7)
         VARB(2) = IBSET(ISHFT(IAND(N, 16256), -7), 7)
         VARB(3) = IAND(N, 127)
         P = 3
         RETURN

*          4 byte size
       ELSE IF (N .LE. 268435455) THEN
         VARB(1) = IBSET(ISHFT(IAND(N, 266338304), -21), 7)
         VARB(2) = IBSET(ISHFT(IAND(N, 2080768), -14), 7)
         VARB(3) = IBSET(ISHFT(IAND(N, 16256), -7), 7)
         VARB(4) = IAND(N, 127)
         P = 4
         RETURN

*          Out of bounds
       ELSE
         PRINT *, 'number too large for variable-format'
       END IF
       P = -1
       RETURN
      END  ! of wvarb
*------------------------ End of file scribe.f -----------------------
