CFTN7X,Q,T,S
C$CDS ON
c      INCLUDE 'CALPGM.INC'
C  supplimentary routines for CALFIT : SUBFIT.FOR
C
C   Copyright (C) 1989, California Institute of Technology
C   All rights reserved.  U. S. Government Sponsorship under
C   NASA Contract NAS7-918 is acknowledged.
C
C   Herbert M. Pickett, 24 Feb 1989
C
      SUBROUTINE GETLBL(IDPAR,PARLBL,LU,FIL,IDIV)
C
C  gets parameter label from user file
C
C     IDPAR = list of parameter ID numbers , IDPAR(0) = number of parameters
C     PARLBL = returned list of parameter labels
C     LU = unit to use for user file
C     FIL = file name to use for userfile
C     IDIV = divisor for vibrational field
C
      INTEGER LU,IDIV
      REAL*8 IDPAR(0:*)
      CHARACTER*(*) PARLBL(*),FIL
      CHARACTER*64 FNAM
      INTEGER NPAR,IBGN,KNT,I,JDIV
C      EMA IDPAR
C
      INTEGER*4 IDTMP,JDTMP
      CHARACTER*30 PARX
      NPAR=IDPAR(0)
      IBGN=1
      KNT=NPAR
C     .. initialize labels
      DO 1 I=1,NPAR
 1        PARLBL(I)=' '
      FNAM=FIL
C     .. open file of parameter names
 2    CALL OPENLU(LU,FNAM,'RL')
C     .. first line contains divisor to ignore lower digits of ID
      READ(LU,*,END=99,ERR=99) JDIV
      JDIV=MAX(JDIV,1)
C
C     DO until EOF
C
 20   READ(LU,'(I10,A)',END=99,ERR=99) JDTMP,PARX
C     .. negative ID signals that parameter name field is a new file
C        to find further parameter names
      IF ( JDTMP. LT.0 ) THEN
          CALL OPENLU(LU,FNAM,'C')
          FNAM=PARX
          GO TO 2
      ENDIF
      JDTMP=JDTMP/JDIV
      DO 10 I=IBGN,NPAR
          IDTMP=ABS(IDPAR(I))/IDIV+0.5
          IF (JDTMP.EQ.IDTMP) THEN
              IF (PARLBL(I).EQ.' ') THEN
                  PARLBL(I)=PARX
                  IF (I.EQ.IBGN) IBGN=IBGN+1
                  KNT=KNT-1
              ENDIF
          ENDIF
 10       CONTINUE
      IF (KNT.GT.0) GO TO 20
 99   CALL OPENLU(LU,FNAM,'C')
      RETURN
      END
      SUBROUTINE GETBLK(IBLK,INDX,IQNUM,NBLKPF,IPOS,NQN)
      INTEGER IBLK,INDX,NBLKPF,IPOS,NQN,IQNUM(6)
C
C  gets block (IBLK) and INDEX from quantum numbers (IQNUM)
C
C
C     NBLKPF IS THE NUMBER OF BLOCKS PER "F"
C     IPOS   IS THE POSITION IN IQNUM TO FIND "F"
C     NQN    IS THE NUMBER OF QUANTUM NUMBERS
C
      INTEGER I0,NSIZE,KQNUM(6),IQNF,IBGN,IDBLK,NCOD,IDGN,
     :        NN,KBGN,K,KDTAU,KK,NBLK
      PARAMETER (I0=0)
      INDX=0
      IBLK=0
C   ..check for indication of null level
      IF (IQNUM(1).LT.0) RETURN
      IQNF=IQNUM(IPOS)
      IBGN=IQNF
      NBLK=NBLKPF
      IF(IQNF.LT.0) THEN
          IBGN=-IBGN
          K=10-IBGN
          IF(K.GT.1) NBLK=NBLK*K
      ENDIF    
      IBGN=IBGN*NBLKPF
      IQNUM(IPOS)=IQNUM(1)
C   ..loop over blocks for given F
      DO 10 IDBLK=1,NBLK
          IBLK=IBGN+IDBLK
          CALL GETQN(IBLK,I0,KQNUM,NSIZE,NCOD)
          INDX = 1
C       ..search for match within block
C         DO WHILE(INDX.LE.NSIZE)
  30          IF (INDX.LE.NSIZE) THEN
              CALL GETQN(IBLK,INDX,KQNUM,IDGN,NCOD)
              KQNUM(IPOS)=KQNUM(1)
              NN=ABS(NCOD)
C           ..NN is the size of a wang sub-block
C           ..NCOD < 0 if oblate basis
C
              IF (NN.LE.1) THEN
                  NN=1
                  KBGN=2
              ELSE
                  KBGN=4
              ENDIF
              DO 20 K=NQN,KBGN,-1
                  IF(IQNUM(K).NE.KQNUM(K)) GO TO 15
  20              CONTINUE
C           ..all quanta tested equal
              IF (NN.EQ.1) GO TO 100
C           ..assignment could be in this sub-block
              KDTAU=IQNUM(2)-IQNUM(3)-KQNUM(2)+KQNUM(3)
              KK=KDTAU/4 
              IF (KK*4.EQ.KDTAU) THEN
C               ..symmetry is good
                  IF (NCOD.LT.0) KK=-KK
                  IF (KK.GE.0.AND.KK.LT.NN) THEN
C                   ..value is in range
                      INDX=INDX+KK
                      GO TO 100
                  ENDIF
              ENDIF
C           ..to get here, quanta were not right
  15          INDX=INDX+NN
              GO TO 30
          ENDIF
C         END DO
  10      CONTINUE
C   ..no match found
      IBLK=0
      NBLK=0
C   .. standard return
 100  IF(NBLK.GT.NBLKPF) IQNF=IQNF-(IBLK-IBGN-1)/NBLKPF
      IQNUM(IPOS)=IQNF
      IF(IQNF.LT.0) INDX=-INDX
      RETURN
      END
      SUBROUTINE DTRI(T,LDT,N,NS)
      INTEGER LDT,N,NS
      REAL*8 T(LDT,1)
C      EMA T
C
C     DTRDI COMPUTES THE DETERMINANT AND INVERSE OF A REAL
C     TRIANGULAR MATRIX.
C
C     ON ENTRY
C
C        T       REAL(LDT,N)
C                T CONTAINS THE LOWER TRIANGULAR MATRIX. THE ZERO
C                ELEMENTS OF THE MATRIX ARE NOT REFERENCED, AND
C                THE CORRESPONDING ELEMENTS OF THE ARRAY CAN BE
C                USED TO STORE OTHER INFORMATION.
C                NS ADDITIONAL ROWS USED FOR Y VECTORS
C
C        LDT     INTEGER
C                LDT IS THE LEADING DIMENSION OF THE ARRAY T.
C
C        N       INTEGER
C                N IS THE ORDER OF THE SYSTEM.
C
C        NS      INTEGER
C                NS IS THE NUMBER OF SOLUTION VECTORS
C     ON RETURN
C
C        T       INVERSE OF ORIGINAL MATRIX PLUS SOLUTIONS
C
C
C     ADAPTED FROM LINPACK.
C
C     SUBROUTINES AND FUNCTIONS
C
C     BLAS DAXPY,DSCAL
C
C     INTERNAL VARIABLES
C
      REAL*8 TEMP
      INTEGER J,K,I1,NT
      PARAMETER (I1=1)
C
      NT=NS
C     .. COMPUTE INVERSE OF LOWER TRIANGULAR
C
      DO 150 K = N ,1, -1
         TEMP = 1 / T(K,K)
         T(K,K) = TEMP
         CALL DSCAL(NT,TEMP,T(K+1,K),I1)
         NT=NT+1
         DO 130 J = 1, K - 1
            TEMP = -T(K,J)
            T(K,J) = 0
            CALL DAXPY(NT,TEMP,T(K,K),I1,T(K,J),I1)
  130       CONTINUE
  150    CONTINUE
      RETURN
      END
      SUBROUTINE FILBAK(LU,FLU,LUBAK,FBAK)
C     make backup par file with / in column 80
C     LU = input file unit number
C     FLU= input file name (.par)
C     LUBAK=backup file unit number
C     FBAK =backup file name (.bak)
C
      INTEGER LU,LUBAK,K
      CHARACTER*(*) FLU,FBAK
      CHARACTER*80 CLINE
      CALL OPENLU(LU,FLU,'R')
      CALL OPENLU(LUBAK,FBAK,'WP')
  10      READ(LU,100,END=50) CLINE
          DO 15 K=80,1,-1
              IF(CLINE(K:K).NE.' ') THEN
                  WRITE(LUBAK,100) CLINE(1:K)
                  GO TO 10
              ENDIF
  15          CONTINUE
          GO TO 10
  50  CALL OPENLU(LU,FLU,'C')
      CALL OPENLU(LUBAK,FBAK,'C')
 100  FORMAT(A)
      RETURN
      END
      SUBROUTINE PRCORR(LUFIT,NPAR,COR,NDCOR,ERR,IDPAR)
      INTEGER LUFIT,NPAR,NDCOR,J,I1,N,NDEL,I,NPARTOT
      REAL*8 COR(NDCOR,*),ERR(*),IDPAR(0:*),DDOT,VAL
C  calculates and prints correlation coefficients 
      PARAMETER (I1=1)
C      EMA VAR,ERPAR 
      NPARTOT=IDPAR(0)
      NDEL=NPAR+1
      I=0
      DO 400 J=1,NPARTOT
        IF(IDPAR(J).GE.0) THEN
          I=I+1
          VAL=1./ERR(J)
          N=NDEL-I
          CALL DSCAL(N,VAL,COR(I,I),I1)
        ENDIF  
  400   CONTINUE
      DO 480 I=NPAR,1,-1
          N=NDEL-I
          DO 470 J=1,I
  470       COR(J,I)=DDOT(N,COR(I,J),I1,COR(I,I),I1)
  480     CONTINUE
      DO 490 I=2,NPAR
          N=I-1
          CALL DCOPY(N,COR(I1,I),I1,COR(I,I1),NDCOR)
  490     CONTINUE
      WRITE (LUFIT,690) ((I,J,COR(J,I),J=1,NPAR),I=1,NPAR)
      RETURN
  690 FORMAT (8(2I3,F10.6))
      END
C$EMA //
      SUBROUTINE DNUADD(NPAR,NPARX,INITL,INDX,IFAC,EGY,
     :                  EGYDER,NSIZE,LINE,PAR,FAC)
C     subroutine to add energies and derivatives to frequency lists
C  ON INPUT:
C     NPAR   = number of parameters
C     NPARX  = number of parameters if INDX < 0
C     INITL  = first energy computed for this frequency  IF < 0
C     INDX   = index in energy list
C     IFAC   = scaling for energies
C     EGY    = energy vector
C     EGYDER = energy derivative
C     NSIZE   = column size of EGYDER
C     LINE   = line counter for frequency derivatives in the heap
C     PAR    = parameters
C     FAC    = scaling values
C
      INTEGER NPAR,NPARX,INITL,INDX,IFAC,NSIZE,LINE
      REAL*8 FAC(0:3),EGY(*),EGYDER(*),PAR(*),
     :       HEAP(0:9000),DDOT,F,ZERO(2)
      INTEGER*4 NDEL,LDNUDP,LBUFOF
      INTEGER RDWRIT,I0,I1,NPARN,NOFF,ITMP
      PARAMETER (RDWRIT=0,I0=0,I1=1)
      COMMON HEAP
      SAVE ZERO
      DATA ZERO/0.,0./
C      EMA PAR,EGY,EGYDER
C..................
      LDNUDP=LBUFOF(RDWRIT,LINE)
      NDEL=NPAR+LDNUDP
      IF (INITL.LT.0) HEAP(NDEL)=0
      F=FAC(IFAC)
C
C   ..add derivatives
      IF (INDX.LT.0) THEN
C       ..set up to ignore last few derivatives
          INDX=-INDX
          NPARN=NPARX
          NOFF=NPAR-NPARN
          ITMP=INDX+NPARN*NSIZE
C       ..add energies subtracting effect of ignored derivatives
          HEAP(NDEL)=HEAP(NDEL)+F*(EGY(INDX)-
     :      DDOT(NOFF,EGYDER(ITMP),NSIZE,PAR(NPARN+1),I1))
          IF (INITL.LT.0) CALL DCOPY(NOFF,ZERO,I0,HEAP(LDNUDP+NPARN),I1)
      ELSE
          NPARN=NPAR
C       ..add energies
          HEAP(NDEL)=HEAP(NDEL)+F*EGY(INDX)
      ENDIF
C     .. sum derivatives
      IF (INITL.LT.0) THEN
          CALL DCOPY(NPARN,EGYDER(INDX),NSIZE,HEAP(LDNUDP),I1)
          IF(IFAC.NE.0) CALL DSCAL(NPARN,F,HEAP(LDNUDP),I1)
      ELSE
          CALL DAXPY(NPARN,F,EGYDER(INDX),NSIZE,HEAP(LDNUDP),I1)
      ENDIF
      RETURN
      END
C$EMA //
      SUBROUTINE DNUGET(IFLG,NPAR,F,LINE,DVEC)
C     subroutine to get derivatives and add to DVEC
C  ON INPUT:
C     IFLG = flag to determine type of operation
C     NPAR = number of parameters
C     F    = scale factor (if IFLG <> 0 )
C     LINE   = line counter for frequency derivatives in the heap
C  RETURNS:
C     F    = frequency ( if IFLG = 0 )
C     DVEC = modified vector of derivatives and frequency
C
      INTEGER IFLG,NPAR,LINE
      REAL*8  F,DVEC(*),HEAP(0:9000)
      INTEGER*4 LDNUDP,LBUFOF
      INTEGER I1,RDFLG
      PARAMETER (I1=1,RDFLG=1)
      COMMON HEAP
C      EMA DVEC
C..................
      LDNUDP=LBUFOF(RDFLG,LINE)
      IF (IFLG.GT.0) THEN
C         .. copy scaled derivatives into DVEC
          CALL DCOPY(NPAR,HEAP(LDNUDP),I1,DVEC,I1)
          CALL DSCAL(NPAR,F,DVEC,I1)
      ELSE IF (IFLG.LT.0) THEN
C         .. sum scaled derivatives with DVEC
          CALL DAXPY(NPAR,F,HEAP(LDNUDP),I1,DVEC,I1)
      ELSE
C         .. get calculated frequency
          F=HEAP(LDNUDP+NPAR)
      ENDIF
      RETURN
      END
C
      FUNCTION LBUFOF(IFLG,IPOS)
C
C  FUNCTION TO HANDLE ADDRESSES IN HEAP , WITH POSSIBLE STORE
C
      INTEGER*4 LBUFOF,KDEL,IORG,RDORG
      LOGICAL REWRIT
      INTEGER IFLG,IPOS,BPOS,RELPOS,NEXT,NBUF,NREAD,NSAVE,
     &         MAXREC,INC,NPOS,LU,K
      SAVE NEXT,KDEL,BPOS,NBUF,IORG,MAXREC,INC,LU,REWRIT
      DATA MAXREC/0/
C..................
      IF(IFLG.LE.1.AND.IFLG.GE.-2) THEN
C       .. find offset from buffer origin, IFLG = -2, -1, 0, 1
          NPOS=IPOS
          RELPOS=NPOS-BPOS
          NSAVE=NEXT
          NREAD=1
          IF(RELPOS.GE.NEXT) THEN
              IF(RELPOS.LT.NBUF) THEN
                  K=RELPOS-NEXT+1
                  IF(K.LE.INC) THEN
C                     record nearly next in buffer
                      NREAD=K
                      NSAVE=0
                      NPOS=NEXT+BPOS
                      RDORG=KDEL*NEXT
                      NEXT=RELPOS+1
                  ENDIF
              ENDIF
              IF(IPOS.GT.MAXREC) THEN
                  MAXREC=IPOS
                  NREAD=NREAD-1
              ENDIF
          ELSE IF(RELPOS.GE.0) THEN
C             record within buffer
              NSAVE=0
              NREAD=0
          ENDIF
C         .. save old data
          IF(NSAVE.GT.0) THEN
              IF(REWRIT) CALL DSAVE(LU,KDEL,BPOS,NSAVE)
              BPOS=NPOS
              RDORG=0
              LBUFOF=0
              NEXT=1
              REWRIT=.FALSE.
          ELSE
              LBUFOF=KDEL*RELPOS
          ENDIF
          IF(ABS(IFLG).NE.1) REWRIT=.TRUE.
          IF(NREAD.GT.0) CALL DGET(LU,KDEL,RDORG,NPOS,NREAD)
          IF(IFLG.LT.0)  LBUFOF=LBUFOF+IORG
      ELSE IF(IFLG.EQ.-4) THEN
C       .. save record length
          KDEL=IPOS
C       .. compute read grouping (1024 bytes)
          INC=128/IPOS+1
C       .. return record length
          LBUFOF=KDEL
      ELSE IF(IFLG.EQ.-3) THEN
C       .. save origin
          IORG=IPOS
          LBUFOF=IORG
      ELSE IF(IFLG.EQ.-5) THEN
C       .. initialize
C             IPOS is blank common size in records
          NBUF=IPOS
          IF(MAXREC.EQ.0) THEN
              LU=-1
              NEXT=0
              BPOS=1
              REWRIT=.FALSE.
          ELSE IF(NEXT.GT.NBUF) THEN
C           .. more records now than new size
              IF(REWRIT.OR.LU.LT.0) CALL DSAVE(LU,KDEL,BPOS,NEXT)
              NEXT=NBUF
              REWRIT=.FALSE.
          ENDIF
C       .. return maximum number of lines
          LBUFOF=NBUF
      ENDIF
      RETURN
      END
C$EMA //
      SUBROUTINE DSAVE(LU,RLEN,KREC,NREC)
C
C  FUNCTION TO SAVE HEAP
C
      INTEGER LU,KREC,NREC
      INTEGER*4 RLEN,LBGN
      INTEGER IOS,I
      REAL*8 HEAP(0:9000)
      COMMON HEAP
C
C     SAVE HEAP TO DISK
C..................
      IF(LU.LT.0) THEN
          LU=RLEN
          CALL OPENLU(LU,'Heap Scratch','S')
      ENDIF
      LBGN=-1
 5    IF(NREC.LE.0) RETURN
      WRITE(LU,REC=KREC,ERR=9,IOSTAT=IOS) (HEAP(LBGN+I),I=1,RLEN)
      NREC=NREC-1
      LBGN=LBGN+RLEN
      KREC=KREC+1
      GO TO 5
C
 9    WRITE(*,*) ' Heap write ERROR, IOSTAT=',IOS
      STOP 'Aborted'
      END
C$EMA //
      SUBROUTINE DGET(LU,RLEN,IORG,KREC,NREC)
C
C  FUNCTION TO READ HEAP
C
      INTEGER LU,KREC,NREC
      INTEGER*4 RLEN,IORG
      INTEGER IOS,I
      REAL*8 HEAP(0:9000)
      COMMON HEAP
C
C     GET HEAP FROM DISC
      IORG=IORG-1
 10   READ(LU,REC=KREC,ERR=19,IOSTAT=IOS) (HEAP(IORG+I),I=1,RLEN)
      NREC=NREC-1
      IF(NREC.LE.0) RETURN
      IORG=IORG+RLEN
      KREC=KREC+1
      GO TO 10
C
 19   WRITE(*,*) ' Heap read ERROR, IOSTAT=',IOS
      STOP 'Aborted'
      END
C
C     this set of subroutines manipulates a record structure
C         of mixed type. the record elements are:
C
C     1. DNUDP = freq. derivatives      (REAL*8 VECTOR)
C     2. CFRQ  = calc. frequency        (REAL*8)
C     3. XFRQ= line frequency           (REAL*8)
C     4. XWT = weight of line           (REAL*8)
C     5. XERR= error  of line           (REAL*8)
C     6. LINKU=link to next upper state (INTEGER*2)       1
C     7. LINKL=link to next lower state (INTEGER*2)       2
C     8. IBU = upper state block number (INTEGER*2)       3
C     9. INU = upper state index number (INTEGER*2)       4
C    10. IBL = lower state block number (INTEGER*2)       5
C    11. INL = lower state index number (INTEGER*2)       6
C    12. IBLN= blend flag               (INTEGER*2)       7
C    11. IQN = quantum number input     (INTEGER*2) *12   8 ... 19
C
      SUBROUTINE HEAPNU(NPAR,NLIM)
      INTEGER   NPAR
      INTEGER*4 NLIM
C    .. set up heap information  for MAIN
C      ON INPUT:
C         NPAR = NUMBER OF PARAMETERS
C         NLIM = LENGTH OF BLANK COMMON IN REAL*8 UNITS
C      ON RETURN:
C         NLIM = POSITION OF IBU ELEMENT OF RECORD
      INTEGER*4 LBUFOF,IADR
      INTEGER   NDBLE,IDM
      PARAMETER (IDM=19)
C     NDBLE returns number of REAL*8 words for number of
C           INTEGER*2 words
C     IDM is the number of INTEGER*2 values to be saved
      INTEGER I,NDEL,IORG,IVAL
C..................
      IORG=NPAR+4
C         .. save origin
      I=-3
      IORG=LBUFOF(I,IORG)
      NDEL=IORG+NDBLE(IDM)
C         .. save record length , return size of heap sub-blocks
      I=-4
      NDEL=LBUFOF(I,NDEL)
      IADR=NLIM/NDEL
      IF(IADR.LT.2) IADR=2
C     .. fix maximum buffer size
      I=-5
      IVAL=IADR
      IF(LBUFOF(I,IVAL).LT.0) THEN
          IVAL=1
          I=LBUFOF(IVAL,IVAL)
      ENDIF
      NLIM=IADR*NDEL
      RETURN
      END
C$EMA //
      SUBROUTINE GETDBK(LINK,II,IBLK,INDX,INITL,IFAC)
      INTEGER LINK,II,IBLK,INDX,INITL,IFAC
C
C   find data for line = abs(LINK)      neg sign points to lower
C   returns II  = line number
C           IBLK= requested block number
C           INDX= requested index
C           INITL < 0 if first block
C           IFAC= code for energy conversion
C                 0 for 1., 1 for -1., 2 for c, 3 for -c
      INTEGER*4 LBUFOF,IADR
      INTEGER*2 LVEC(6)
      INTEGER   RDORG,I1,I6
      PARAMETER (I1=1,I6=6,RDORG=-1)
      REAL*8 HEAP(0:9000),LVECD(2)
      COMMON HEAP
C      EMA LVEC
      EQUIVALENCE (LVECD,LVEC)
C..................
      II=ABS(LINK)
      IADR=LBUFOF(RDORG,II)
      CALL ICOPY(I1,I6,HEAP(IADR),LVECD)
      INITL=LVEC(3)-LVEC(5)
      IF(LINK.GE.0) THEN
          LINK=LVEC(1)
          IBLK=LVEC(3)
          INDX=LVEC(4)
          IF(INITL.EQ.IBLK.OR.INITL.EQ.0) INITL=-1
          IFAC=0
      ELSE
          LINK=LVEC(2)
          IBLK=LVEC(5)
          INDX=LVEC(6)
          INITL=-INITL
          IFAC=1
      ENDIF
C     if XERR < 0 use cm-1
      IF(HEAP(IADR-1).LT.0.) IFAC=IFAC+2
      RETURN
      END
C$EMA //
      SUBROUTINE FRQDAT(II,IBLN,XFRQ,XWT,XERR,IQN)
      INTEGER   II,IBLN,IQN(12)
      REAL*8 XFRQ,XWT,XERR
C   .. gets blend code, frequency, weight, error , and quantum numbers
      INTEGER*4 LBUFOF,IADR
      INTEGER   RDORG,LBLND,NT,I
      PARAMETER (LBLND=7,RDORG=-1,NT=19)
      INTEGER*2 LVEC(NT)
      REAL*8 HEAP(0:9000),LVECD(2)
      COMMON HEAP
C      EMA LVEC
      EQUIVALENCE (LVECD,LVEC)
C..................
      IADR=LBUFOF(RDORG,II)
      CALL ICOPY(LBLND,NT,HEAP(IADR),LVECD)
      IBLN=LVEC(LBLND)
      XFRQ=HEAP(IADR-3)
      XWT= HEAP(IADR-2)
      XERR=HEAP(IADR-1)
      DO 11 I=1,12
 11       IQN(I)=LVEC(I+LBLND)
      RETURN
      END
C$EMA //
      SUBROUTINE SAVDAT(II,LKU,LKL,IBU,INU,IBL,INL,IBLN,
     :          XFRQ,XWT,XERR,IQN)
      INTEGER   II,LKU,LKL,IBU,INU,IBL,INL,IBLN,IQN(12)
      REAL*8 XFRQ,XWT,XERR
C   .. saves block numbers, indices, blend code, frequency, weight,
C   ..        error , and quantum numbers
      INTEGER*4 LBUFOF,IADR
      INTEGER   WRDORG,I1,NT,I
      PARAMETER (WRDORG=-2,I1=1,NT=19)
      INTEGER*2 LVEC(NT)
      REAL*8 HEAP(0:9000),LVECD(2)
C      EMA LVEC
      COMMON HEAP
      EQUIVALENCE (LVECD,LVEC)
C..................
      IADR=LBUFOF(WRDORG,II)
      LVEC(1)=LKU
      LVEC(2)=LKL
      LVEC(3)=IBU
      LVEC(4)=INU
      LVEC(5)=IBL
      LVEC(6)=INL
      LVEC(7)=IBLN
      DO 21 I=1,12
 21       LVEC(I+7)=IQN(I)
      CALL ICOPY(I1,NT,LVECD,HEAP(IADR))
      HEAP(IADR-3)=XFRQ
      HEAP(IADR-2)=XWT
      HEAP(IADR-1)=XERR
      RETURN
      END
C$EMA //
      SUBROUTINE SAVXWT(II,XWT,IBLN)
      INTEGER   II,IBLN
      REAL*8 XWT
C   .. renormalize weight and update blend code
      INTEGER*4 LBUFOF,IADR
      REAL*8 HEAP(0:9000),LVECD(2)
      INTEGER*2 LVEC(19)
      INTEGER   WRDORG,LBLND
      PARAMETER (WRDORG=-2,LBLND=7)
C      EMA LVEC
      COMMON HEAP
      EQUIVALENCE (LVECD,LVEC)
C..................
      IADR=LBUFOF(WRDORG,II)
      LVEC(LBLND)=IBLN
      CALL ICOPY(LBLND,LBLND,LVECD,HEAP(IADR))
      HEAP(IADR-2)=HEAP(IADR-2)/XWT
      RETURN
      END
C$EMA //
      INTEGER FUNCTION LNLINK(HEAD,PRVBLK,NBLK,IBLK,LINE,PLINK,FLG)
C     find place to insert link for line
C     return pointer to next line
C
C          PRVBLK contains pointer to last line using given block
C
      INTEGER*2 PRVBLK(*),LVEC(19)
      INTEGER   HEAD,NBLK,IBLK,LINE,PLINK,FLG,LFLG,I,IBGN,LAST,NOW
      INTEGER*4 LBUFOF,IADR
      INTEGER   WRDORG
      PARAMETER (WRDORG=-2)
      REAL*8    HEAP(0:9000),LVECD(2)
      COMMON HEAP
C      EMA LVEC
      EQUIVALENCE (LVECD,LVEC)
C..................
      IF(IBLK.LE.0) THEN
          LNLINK=0
          RETURN
      ENDIF
C     search down for non zero element of PRVBLK
      LFLG=FLG
      IBGN=IBLK
      IF(IBGN.GT.NBLK) IBGN=NBLK
      DO 20 I=IBGN,1,-1
          LAST=PRVBLK(I)
          IF(LAST.NE.0) GO TO 30
 20       CONTINUE
 30   PRVBLK(IBGN)=LINE
      IF(LAST.EQ.0) THEN
          LNLINK=HEAD
          HEAD=LINE
          IF(LFLG.EQ.0) RETURN
          LFLG=0
          LAST=LINE
          NOW=LNLINK
      ELSE IF(LAST+LINE.EQ.0.AND.LFLG.EQ.0) THEN
C       special case of in-line link
          LNLINK=PLINK
          PLINK=LINE
          RETURN
      ELSE
          NOW=LINE
      ENDIF
  40  IF(LAST.GE.0) THEN
          IBGN=1
          I=LAST
      ELSE
          IBGN=2
          I=-LAST
      ENDIF
      IADR=LBUFOF(WRDORG,I)
      CALL ICOPY(IBGN,IBGN,HEAP(IADR),LVECD)
      LNLINK=LVEC(IBGN)
      LVEC(IBGN)=NOW
      CALL ICOPY(IBGN,IBGN,LVECD,HEAP(IADR))
      IF(LFLG.EQ.0) RETURN
        LAST=LINE
      NOW=LNLINK
      LFLG=0
      GO TO 40
      END