C------------------------------------------------------------------------
C	* * *  THE FOLLOWING IS A COLLECTION OF MISC ROUTINES  * * *
C------------------------------------------------------------------------

    	FUNCTION IFMISC_DECODE( TEXT, NF, GF, KF )
C
C  PROGRAMMER:  R.S. ZHOU
C		ALFRED UNIVERSITY
C		NOV-1986
C
C	THIS ROUTINE IMPLEMENTS A FREE FORMAT INPUT PROCESSOR.
C
C	TEXT --- INPUT TEXTSTRING IN CHARACTER FORMAT TO BE DECODED.
C	NF ----- NUMBER OF REAL*4 DECIMALS TO BE DECODED FROM THE "TEXT".
C	GF ----- REAL*4 ARRAY TO STORE THE DECODED DECIMALS.
C	KF ----- INTEGER ARRAY TO FLAG WHICH ELEMENT OF THE "GF"
C	         ARRAY GETS AN INPUT: 0=NO, 1=YES.
C
C	IFMISC_DECODE:
C
C		 RETURN WITH ZERO IF ERROR OCCURS IN DECODING OR ELSE
C		 THE POSITION OF THE LAST CHARACTER IN THE "TEXT"
C	         THAT DOESN'T GET DECODED. IN OTHER WORDS, IT FLAGS
C	         HOW MANY CHARACTERS HAVE BEEN USED, COUNTED FROM
C	         THE LEFT TO RIGHT IN THE "TEXT", TO DECODE THE "NF"
C	         REAL*4 DECIMALS.
C
C	COMMENTS:
C
C	THE FIELD SEPARATOR CAN BE EITHER A SPACE OR SPACES OR COMMA.
C	TO SKIP A FIELD, YOU HAVE TO USE THE COMMAS AND THE INPUT
C	PROCESSOR RETURNS WITH ITS "KF" = 0.
C
C	INTEGER FIELDS ARE DECODED AND STORED IN THE "GF" AS DECIMAL
C	NUMBERS. THEREFORE DECIMAL NUMBERS CAN BE ENTERED AS INTEGERS
C	AND VICE VERSA.
C
C	AN ERROR RESULTS IF THE INPUT PROCESSOR ENCOUNTERS TRAILING
C	TEXT CHARACTERS IN THE "TEXT" BEFORE THE "NF"TH DATA FIELD
C	GOT DECODED FROM THE "TEXT".
C	IF THERE IS NO TRAILING TEXT CHARACTERS IN THE "TEXT", THE
C	"TEXT" CAN HAVE FEWER THAN "NF" NUMBER OF DATA FIELDS AND
C	THE REST OF "GF" ARE FILLED WITH VALUES 0.0s AND THEIR "KF"s
C	WITH 0s.
C
    	CHARACTER SYM*1, BEL*1, STRN*40, TEXT*(*)
     	INTEGER*2 JMS,NLS,I1,I2,ICI, NCS(50),NSP(50)
    	DIMENSION KF(NF), GF(NF)

    	NFS = NF
    	DO I=1, NFS
    	   KF(I) = 0
    	   GF(I) = 0.0
    	ENDDO
    	NLS = 0
    	NCH = IFMISC_LENG( TEXT )
    	IF( NCH.EQ.0 ) GOTO 999
    	JMS = 1
    	NCS(JMS) = 0
    	NSP(JMS) = 0
    	BEL = ','
    	ICI = 0

    	DO 9 I=1, NCH
    	SYM = TEXT(I:I)
     	IF( SYM.EQ.'!' ) GOTO 2
    	IF( SYM.NE.',' ) GOTO 1
    	    IF( BEL.NE.' ' .OR. ICI.EQ.0 ) THEN
    	        JMS = JMS + 1
    		IF( ICI.EQ.0 ) THEN
    		    NSP(JMS) = I - NCS(JMS-1)
    		ELSE
    		    NSP(JMS) = 0
    		ENDIF
    	    ELSE
    	        NSP(JMS) = I - NCS(JMS)
    	    ENDIF
    	    NCS(JMS) = I
    	    NLS = I + 1
    	    ICI = 0
    	GOTO 6
1	IF( SYM.NE.' ' .OR. BEL.EQ.' ' .OR. BEL.EQ.',' ) GOTO 5
    	    JMS = JMS + 1
    	    NCS(JMS) = I
    	    NSP(JMS) = 0
    	GOTO 8
5	IF( SYM.EQ.' ' ) GOTO 8
    	NLS = I
6	IF( JMS.GT.NFS .OR. JMS.EQ.50 ) GOTO 10
8	BEL = SYM
    	IF( SYM.NE.' ' .AND. SYM.NE.',' ) ICI = 1
9	CONTINUE

2	NLS = I
    	JMS = JMS + 1
    	NCS(JMS) = I
    	NSP(JMS) = 0

10	DO I=1,JMS-1
    	   I1 = NCS(I) + 1
    	   I2 = NCS(I+1) - 1 - NSP(I+1)
    	   NCH = I2 - I1 + 1
     	   IF( NCH.GT.39 ) GOTO 888
    	   IF( NCH.GT.0 ) THEN
     	       STRN = TEXT(I1:I2)
    	       ICI  = INDEX( STRN(1:NCH), '.' )
    	       IF( ICI.EQ.0 ) THEN
    		   NCH = NCH + 1
    		   STRN(NCH:NCH) = '.'
    	       ENDIF
    	       DECODE( NCH, 20, STRN, ERR=888 ) GF(I)
    	       KF(I) = 1
    	   ENDIF
    	ENDDO
20	FORMAT(G40.30)
    	GOTO 999
888	NLS = 0
999	IFMISC_DECODE = NLS
     	RETURN
    	END

    	FUNCTION IFMISC_STRIP( TEXT )
C
C	THIS ROUTINE RETURNS THE LENGTH OF CHARACTER STRING "TEXT"
C	STRIPPED OFF LEADING AND TRAILING BLANKS.
C
    	CHARACTER TEXT*(*)

    	N = IFMISC_LENG( TEXT )
     	DO  J=1, N		!CONVERT "INVISIBLE" ASCII-CHARS TO SPACES
     	    IF( ICHAR(TEXT(J:J)).LT.32 ) TEXT(J:J) = ' '
     	ENDDO
	DO  K = N, 1, -1	!REMOVE TRAILING BLANKS
    	    IF( TEXT(K:K).NE.' ' ) GOTO 1
    	ENDDO
	K = 0
    	GOTO 9
  1	DO  J=1, K		!REMOVE LEADING BLANKS
     	    IF( TEXT(J:J).NE.' ' ) GOTO 2
     	ENDDO
     	K = 0
     	GOTO 9
  2	IF( J.GT.1 ) THEN	!REPACK THE STRING
     	    TEXT = TEXT(J:)
     	    K = K-J+1
     	ENDIF
  9	IFMISC_STRIP = K
    	RETURN
	END

     	FUNCTION IFMISC_UPPER( TEXT )
     	CHARACTER TEXT*(*)

     	K = IFMISC_LENG( TEXT )
     	DO  J=1, K
     	    N = ICHAR(TEXT(J:J))
     	    IF( N.GT.96 .AND. N.LT.123 ) N = N - 32
     	    TEXT(J:J) = CHAR(N)
     	ENDDO
     	IFMISC_UPPER = K
     	RETURN
     	END

     	FUNCTION IFMISC_LENG( TEXT )
     	CHARACTER STRN*132, TEXT*(*)

     	STRN = TEXT
     	DO K=132,1,-1
     	   IF( STRN(K:K).NE.' ' ) GOTO 1
     	ENDDO
        K = 0
  1	IFMISC_LENG = K
     	RETURN
     	END

    	FUNCTION IFMISC_NEXT( TEXT )
C
C	THIS ROUTINE RETURNS THE LENGTH OF CHARACTER STRING "TEXT"
C	AFTER THE FIRST SPACE IN "TEXT".
C
    	CHARACTER TEXT*(*)

     	K = IFMISC_LENG( TEXT )
     	J = INDEX( TEXT(1:K), ' ' )
     	IF( J.EQ.0 ) THEN
     	    TEXT = ' '
     	    K = 0
     	ELSE
     	    DO  I=J+1, K
     		IF( TEXT(I:I).NE.' ' ) GOTO 1
     	    ENDDO
  1	    TEXT = TEXT(I:)
     	    K = K-I+1
     	ENDIF
  9	IFMISC_NEXT = K
    	RETURN
	END

     	SUBROUTINE IFMISC_DEFEXT( FNAME, EXT )
C
C  PROGRAMMER:  R.S. ZHOU
C		ALFRED UNIVERSITY, MAY-1987
C
C	THIS ROUTINE ATTACHES THE DEFAULT FILE EXTENSION
C	"EXT" TO A FILE "FNAME" MISSING ITS FILE EXTENSION.
C
     	CHARACTER FNAME*(*), EXT*(*)

     	N = INDEX( FNAME, ']' ) + 1
     	N = INDEX( FNAME(N:), '.' )
     	IF( N.GT.0 ) RETURN		!FILE EXTENSION EXISTS
     	N = INDEX( FNAME, ' ' )
     	FNAME(N:) = EXT

     	RETURN
     	END

    	SUBROUTINE IFMISC_NAMEXT( FNAME, EXT )
C
C  PROGRAMMER: R.S. ZHOU
C
C  	THIS ROUTINE ATTACHES FILE EXTENSION "EXT" TO "FNAME".
C
    	CHARACTER FNAME*(*), EXT*(*)

    	N = INDEX( FNAME, ']' )
    	M = INDEX( FNAME(N+1:), '.' )
    	IF( M.EQ.0 ) THEN
     	    M = INDEX( FNAME, ' ' )
     	ELSE
     	    M = N + M
     	ENDIF
    	FNAME(M:) = EXT

    	RETURN
    	END

    	FUNCTION IFMISC_DIGITS( IT )
C
C	THIS FUNCTION RETURNS WITH THE NUMBER OF DIGITS OF "IT".
C
    	ITT = IABS(IT)
    	IF( ITT.EQ.0 ) THEN
    	    IFMISC_DIGITS = 1
    	ELSE
    	    IFMISC_DIGITS = INT( LOG10(FLOAT(ITT)) ) + 1
    	    IF( IT.LT.0 ) IFMISC_DIGITS = IFMISC_DIGITS + 1
    	ENDIF

    	RETURN
    	END

     	SUBROUTINE IFMISC_SPAWN( COMM, MASK )
C
C  PROGRAMMER:  R.S. ZHOU
C		ALFRED UNIVERSITY, MAY-1987
C
C	THIS ROUTINE SPAWNS A SYSTEM SUBPROCESS USING VMS'S
C	RUN-TIME LIBRARY (RTL) ROUTINE LIB$SPAWN.
C
C	COMM --- DCL COMMAND-STRING.
C	MASK --- FLAGS MASK (SEE RTL MANAUL FOR DETAILS).
C
     	CHARACTER COMM*(*), STRN*132

     	STRN = COMM
     	NCHS = IFMISC_NEXT( STRN )
     	IF( NCHS.EQ.0 ) THEN
     	    CALL IFMAC_KEYEXP(0)
     	    CALL IFINIT_LINE(-24)
     	    CALL IFANSI_VIDEO(47)
     	    WRITE(6,1)
  1	    FORMAT(' *** Type LOGOUT to get back to program ***')
     	    CALL LIB$SPAWN()
     	    CALL IFMAC_KEYEXP(1)
     	ELSE
     	    CALL LIB$SPAWN( STRN,,, MASK )
     	ENDIF

     	RETURN
     	END

    	FUNCTION IFMISC_TYPE( PNAME, KEYWD, ILINE )
C
C  PROGRAMMER:  R.S. ZHOU
C		ALFRED UNIVERSITY, MAY 1986
C
C	THIS ROUTINE DISPLAYS AN ASCII FILE.
C
C	FNAME ---- ASCII FILE IN CHARACTER FORMAT TO BE TYPED.
C	KEYWD ---- KEYWORD TO SEARCH FOR THE DESIRED PARAGRAPH.
C		   KEYWORD MUST BE LEAD BY A LINE FEED <LF>.
C		   IF KEYWORD IS BLANK, THE WHOLE FILE IS DISPLAYED.
C	ILINE ---- ILINE>0: DISPLAYED IN "ILINE" LINES PER PAGE.
C		   ILINE=0: DISPLAYED WITHOUT PAUSE.
C		   ILINE<0: DISPLAYED AND RESET DIALOG AREA LINE AS NEEDED.
C
C	IFMISC_TYPE=0: OK READ.
C	IFMISC_TYPE=1: CAN NOT OPEN/FIND THE FILE.
C	IFMISC_TYPE=2: NO INFORMATION FOR THE INPUT KEYWORD.
C
    	CHARACTER KEY*4, FNAME*32, PNAME*(*), KEYW*4, KEYWD*(*)
	BYTE LF
     	CHARACTER TEXT*80
     	LOGICAL YES

     	IERR  = 0
	FNAME = PNAME
    	INQUIRE( FILE=FNAME, EXIST=YES )
     	IF( .NOT.YES ) GOTO 888
     	OPEN( UNIT=1, FILE=FNAME, STATUS='OLD', FORM='FORMATTED',
     +        CARRIAGECONTROL='LIST', READONLY, ERR=888 )

     	CALL IFANSI_VIDEO(47)
	ILN  = ILINE
     	KEYW = KEYWD
     	NCH  = IFMISC_STRIP(KEYW)
     	CALL IFMISC_UPPER(KEYW)
	IF( NCH.EQ.0 ) THEN
     	    CALL IFANSI_MARGIN(0,0)
	    NPS = 0
	    IF( ILN.GT.0 ) THEN
     		CALL IFINIT_LINE(24)
     	    ELSE
     		CALL IFANSI_CLEAR(0)
     	    ENDIF
1	    READ(1,2,END=111) LF, TEXT
	    IF( LF.EQ.10 ) GOTO 1
     	    NK = IFMISC_LENG( TEXT )
     	    IF( NK.GT.0 ) THEN
     		WRITE(6,3) LF,TEXT(1:NK)
     	    ELSE
     		WRITE(6,4) LF
     	    ENDIF
	    IF( ILN.EQ.0 ) GOTO 1
	    NPS = NPS + 1
	    IF( ILN.GT.0 .AND. ILN.EQ.NPS ) THEN
	        WRITE(6,5)
	        CALL IFQIOW_KEY( NPS )
     		IF( NPS.EQ.26 .OR. NPS.EQ.13 ) GOTO 999
		NPS = 0
     		CALL IFANSI_CLEAR(6)
     		CALL IFANSI_CURSOR(-2,0)
	    ENDIF
	    GOTO 1
2	FORMAT(A1,A)
	ELSE
3	FORMAT(x,A1,A<NK>)
4	FORMAT(x,A1)
10	    READ(1,2,END=777) LF,KEY
	    IF( LF.NE.10 ) GOTO 10
	    IF( KEYW(1:NCH).NE.KEY(1:NCH) ) GOTO 10
	    NPS = 0
     	    CALL IFANSI_CLEAR(0)
	    IF( ILN.GT.0 ) CALL IFINIT_LINE(ILN+2)
20	    READ(1,2,END=111) LF, TEXT
	    IF( LF.EQ.10 ) GOTO 111
     	    NK = IFMISC_LENG( TEXT )
     	    IF( NK.GT.0 ) THEN
     		WRITE(6,3) LF, TEXT(1:NK)
     	    ELSE
     		WRITE(6,4) LF
     	    ENDIF
	    IF( ILN.EQ.0 ) GOTO 20
	    NPS = NPS + 1
	    IF( ILN.GT.0 .AND. ILN.EQ.NPS ) THEN
	        WRITE(6,5)
	        CALL IFQIOW_KEY( NPS )
     		IF( NPS.EQ.26 .OR. NPS.EQ.13 ) GOTO 999
		NPS = 0
     		CALL IFANSI_CLEAR(6)
     		CALL IFANSI_CURSOR(-2,0)
	    ENDIF
	    GOTO 20
	ENDIF
111	IF( NPS.GT.0 ) THEN
	    IF( ILN.LT.0 ) CALL IFINIT_LINE( ILN-NPS )
	    IF( ILN.GT.0 ) THEN
		WRITE(6,5)
	        CALL IFQIOW_KEY( NPS )
	    ENDIF
	ENDIF
	GOTO 999
 5	FORMAT(/'$Press SPACEBAR to continue...')

777	IERR = 2
	GOTO 999
888	IERR = 1
999    	CLOSE(UNIT=1)
     	IFMISC_TYPE = IERR
    	RETURN
    	END

     	FUNCTION IFMISC_HELP(FILE,KEYW,IL,IC)
C
C	THIS ROUTINE DISPLAYS ONE-LINE ON-LINE HELP.
C
C	 FILE: FILE TO BE OPENED AND READ.
C	 KEYW: KEY-WORD TO SEARCH FOR THE LINE.
C	IL,IC: CURSOR POSTION TO PRINT THE LINE.
C
C	IFMISC_HELP=0: OK READ.
C	IFMISC_HELP=1: CAN NOT OPEN/FIND THE FILE.
C	IFMISC_HELP=2: NO INFORMATION FOR THE INPUT KEYWORD.
C
     	CHARACTER HELP*80, KEY*4, KEYW*(*), PNAME*32, FILE*(*)
     	LOGICAL   YES

     	IERR  = 0
     	PNAME = FILE
     	INQUIRE( FILE=PNAME, EXIST=YES )
     	IF( .NOT.YES ) GOTO 888
     	OPEN(UNIT=1,FILE=PNAME,STATUS='OLD',FORM='FORMATTED',
     +       CARRIAGECONTROL='LIST',READONLY,ERR=888)

     	KEY = KEYW
     	NCH = IFMISC_STRIP( KEY )
     	CALL IFMISC_UPPER( KEY )
1	READ(1,2,END=777) HELP
2	FORMAT(A)
     	IF( HELP.EQ.' ' ) GOTO 1
     	IF( HELP(1:NCH).NE.KEY(1:NCH) ) GOTO 1
     	HELP = HELP(NCH+2:)
     	NCH = IFMISC_STRIP( HELP )
     	IF( NCH.GT.80-IC+1 ) NCH = 80-IC+1
     	CALL IFANSI_CURSOR(IL,IC)
     	CALL IFANSI_CLEAR(4)
     	CALL IFANSI_TEXT( HELP, NCH, IL, IC )
     	GOTO 999

777	CALL IFANSI_VIDEO(70)
     	HELP = ' No help found for keyword> '
     	HELP(29:) = KEY
     	NCH = 28+NCH
     	CALL IFANSI_TEXT( HELP, NCH, IL, IC )
     	IERR = 2
     	GOTO 999
888	CALL IFANSI_VIDEO(70)
     	CALL IFANSI_TEXT(' Help file is not available.', 28,IL,IC)
     	IERR = 1
999	CLOSE(UNIT=1)
     	IFMISC_HELP = IERR
     	RETURN
     	END

     	SUBROUTINE IFMISC_RANGE( N, A, B, R )
C
C	THIS ROUTINE FINDS MINIMUM AND MAXIMUM VALUES
C	IN ARRAY A, B.
C
     	DIMENSION A(1), B(1), R(1)

     	NPTS = N
     	AMIN = A(1)
     	AMAX = AMIN
     	BMIN = B(1)
     	BMAX = BMIN
     	DO  I=2, NPTS
     	    AA = A(I)
     	    IF( AA.LT.AMIN ) THEN
     		AMIN = AA
     	    ELSE IF( AA.GT.AMAX ) THEN
     		AMAX = AA
     	    ENDIF
     	    AA = B(I)
     	    IF( AA.LT.BMIN ) THEN
     		BMIN = AA
     	    ELSE IF( AA.GT.BMAX ) THEN
     		BMAX = AA
     	    ENDIF
     	ENDDO
     	R(1) = AMIN
     	R(2) = AMAX
     	R(3) = BMIN
     	R(4) = BMAX
     	RETURN
     	END

     	SUBROUTINE IFMISC_SORT( N, A )
C
C	MODIFIED SINGLETON SORT ROUTINE (CACM 12, 185 (1969)).
C	N IS NUMBER OF POINTS IN THE ARRAYS.
C
     	DIMENSION A(1)
     	INTEGER I, II, IJ, IL(16), IU(16), J, K, L, M
     	REAL S, T

     	II = 1
     	I = II
     	J = N
     	M = 1
  10	IF (I .GE. J) GO TO 100
C
C  ORDER ENDS AND MIDPOINT OF SEGMENT BOUNDED BY I AND J
C
  20	K = I
     	IJ = (I + J)/2
     	T = A(IJ)
     	IF (A(I) .LE. T) GO TO 30
     	A(IJ) = A(I)
     	A(I) = T
     	T = A(IJ)

  30	L = J
     	IF (A(J) .GE. T) GO TO 50
     	A(IJ) = A(J)
     	A(J) = T
     	T = A(IJ)

     	IF (A(I) .LE. T) GO TO 50
     	A(IJ) = A(I)
     	A(I) = T
     	T = A(IJ)
     	GO TO 50
C
C  INTERCHANGE ELEMENT .LE. T WITH ELEMENT .GE. T
C
  40	S = A(L)
     	A(L) = A(K)
     	A(K) = S
C
C  SCAN DOWN SEQUENCE BETWEEN K AND L FOR ELEMENT .LE. T
C
  50	L = L - 1
     	IF (A(L) .GT. T) GO TO 50
C
C  NOW SCAN UP FOR ELEMENT .GE. T
C
  60	K = K + 1
     	IF (A(K) .LT. T) GO TO 60
     	IF (K - L) 40, 70, 80
C
C  AVOID UNNECESSARY INTERCHANGE
C
  70	K = K + 1
     	L = L - 1
C
C  STACK ENDPOINTS OF LONGER SEQUENCE
C
  80	IF (L - I .LE. J - K) GO TO 90
     	IL(M) = I
     	IU(M) = L
     	M = M + 1
     	I = K
     	GO TO 110

  90	IL(M) = K
     	IU(M) = J
     	M = M + 1
     	J = L
     	GO TO 110
C
C  UNSTACK A SAVED SEQUENCE
C
  100	M = M - 1
     	IF (M .EQ. 0) RETURN
     	I = IL(M)
     	J = IU(M)
C
C  SELECT SORTING METHOD
C
  110	IF (J - I .GE. 11) GO TO 20
     	IF (I .EQ. II) GO TO 10
C
C  SHORT SEQUENCES (EXCEPT THE FIRST) ARE ORDERED BY A SINKING SORT
C
     	I = I - 1
  120	I = I + 1
     	IF (I .EQ. J) GO TO 100
     	IF (A(I) .LE. A(I+1)) GO TO 120
     	T = A(I+1)
     	K = I
  130	A(K+1) = A(K)
     	K = K - 1
     	IF (A(K) .GT. T) GO TO 130
     	A(K+1) = T
     	GO TO 120
     	END

     	SUBROUTINE IFMISC_SORT_ONE( N, A, B )
C
C	MODIFIED SINGLETON SORT ROUTINE (CACM 12, 185 (1969)).
C	B IS PASSIVE ARRAY TO BE SORTED ON A.
C	N IS NUMBER OF POINTS IN THE ARRAYS.
C
     	DIMENSION A(1), B(1)
     	INTEGER I, II, IJ, IL(16), IU(16), J, K, L, M
     	REAL S, T, U

     	II = 1
     	I = II
     	J = N
     	M = 1
  10	IF (I .GE. J) GO TO 100
C
C  ORDER ENDS AND MIDPOINT OF SEGMENT BOUNDED BY I AND J
C
  20	K = I
     	IJ = (I + J)/2
     	T = A(IJ)
     	IF (A(I) .LE. T) GO TO 30
     	A(IJ) = A(I)
     	A(I) = T
     	T = A(IJ)
     	U = B(I)
     	B(I) = B(IJ)
     	B(IJ) = U

  30	L = J
     	IF (A(J) .GE. T) GO TO 50
     	A(IJ) = A(J)
     	A(J) = T
     	T = A(IJ)
     	U = B(IJ)
     	B(IJ) = B(J)
     	B(J) = U

     	IF (A(I) .LE. T) GO TO 50
     	A(IJ) = A(I)
     	A(I) = T
     	T = A(IJ)
     	U = B(IJ)
     	B(IJ) = B(I)
     	B(I) = U
     	GO TO 50
C
C  INTERCHANGE ELEMENT .LE. T WITH ELEMENT .GE. T
C
  40	S = A(L)
     	A(L) = A(K)
     	A(K) = S
     	U = B(L)
     	B(L) = B(K)
     	B(K) = U
C
C  SCAN DOWN SEQUENCE BETWEEN K AND L FOR ELEMENT .LE. T
C
  50	L = L - 1
     	IF (A(L) .GT. T) GO TO 50
C
C  NOW SCAN UP FOR ELEMENT .GE. T
C
  60	K = K + 1
     	IF (A(K) .LT. T) GO TO 60
     	IF (K - L) 40, 70, 80
C
C  AVOID UNNECESSARY INTERCHANGE
C
  70	K = K + 1
     	L = L - 1
C
C  STACK ENDPOINTS OF LONGER SEQUENCE
C
  80	IF (L - I .LE. J - K) GO TO 90
     	IL(M) = I
     	IU(M) = L
     	M = M + 1
     	I = K
     	GO TO 110

  90	IL(M) = K
     	IU(M) = J
     	M = M + 1
     	J = L
     	GO TO 110
C
C  UNSTACK A SAVED SEQUENCE
C
  100	M = M - 1
     	IF (M .EQ. 0) RETURN
     	I = IL(M)
     	J = IU(M)
C
C  SELECT SORTING METHOD
C
  110	IF (J - I .GE. 11) GO TO 20
     	IF (I .EQ. II) GO TO 10
C
C  SHORT SEQUENCES (EXCEPT THE FIRST) ARE ORDERED BY A SINKING SORT
C
     	I = I - 1
  120	I = I + 1
     	IF (I .EQ. J) GO TO 100
     	IF (A(I) .LE. A(I+1)) GO TO 120
     	T = A(I+1)
     	U = B(I+1)
     	K = I
  130	A(K+1) = A(K)
     	B(K+1) = B(K)
     	K = K - 1
     	IF (A(K) .GT. T) GO TO 130
     	A(K+1) = T
     	B(K+1) = U
     	GO TO 120
     	END

     	SUBROUTINE IFMISC_SORT_TWO( N, A, B, C )
C
C	MODIFIED SINGLETON SORT ROUTINE (CACM 12, 185 (1969)).
C	B AND C ARE PASSIVE ARRAYS TO BE SORTED ON A.
C	N IS NUMBER OF POINTS IN THE ARRAYS.
C
     	DIMENSION A(1), B(1), C(1)
     	INTEGER I, II, IJ, IL(16), IU(16), J, K, L, M
     	REAL S, T, U, V

     	II = 1
     	I = II
     	J = N
     	M = 1
  10	IF (I .GE. J) GO TO 100
C
C  ORDER ENDS AND MIDPOINT OF SEGMENT BOUNDED BY I AND J
C
  20	K = I
     	IJ = (I + J)/2
     	T = A(IJ)
     	IF (A(I) .LE. T) GO TO 30
     	A(IJ) = A(I)
     	A(I) = T
     	T = A(IJ)
     	U = B(I)
     	B(I) = B(IJ)
     	B(IJ) = U
     	U = C(I)
     	C(I) = C(IJ)
     	C(IJ) = U

  30	L = J
     	IF (A(J) .GE. T) GO TO 50
     	A(IJ) = A(J)
     	A(J) = T
     	T = A(IJ)
     	U = B(IJ)
     	B(IJ) = B(J)
     	B(J) = U
     	U = C(IJ)
     	C(IJ) = C(J)
     	C(J) = U

     	IF (A(I) .LE. T) GO TO 50
     	A(IJ) = A(I)
     	A(I) = T
     	T = A(IJ)
     	U = B(IJ)
     	B(IJ) = B(I)
     	B(I) = U
     	U = C(IJ)
     	C(IJ) = C(I)
     	C(I) = U
     	GO TO 50
C
C  INTERCHANGE ELEMENT .LE. T WITH ELEMENT .GE. T
C
  40	S = A(L)
     	A(L) = A(K)
     	A(K) = S
     	U = B(L)
     	B(L) = B(K)
     	B(K) = U
     	U = C(L)
     	C(L) = C(K)
     	C(K) = U
C
C  SCAN DOWN SEQUENCE BETWEEN K AND L FOR ELEMENT .LE. T
C
  50	L = L - 1
     	IF (A(L) .GT. T) GO TO 50
C
C  NOW SCAN UP FOR ELEMENT .GE. T
C
  60	K = K + 1
     	IF (A(K) .LT. T) GO TO 60
     	IF (K - L) 40, 70, 80
C
C  AVOID UNNECESSARY INTERCHANGE
C
  70	K = K + 1
     	L = L - 1
C
C  STACK ENDPOINTS OF LONGER SEQUENCE
C
  80	IF (L - I .LE. J - K) GO TO 90
     	IL(M) = I
     	IU(M) = L
     	M = M + 1
     	I = K
     	GO TO 110

  90	IL(M) = K
     	IU(M) = J
     	M = M + 1
     	J = L
     	GO TO 110
C
C  UNSTACK A SAVED SEQUENCE
C
  100	M = M - 1
     	IF (M .EQ. 0) RETURN
     	I = IL(M)
     	J = IU(M)
C
C  SELECT SORTING METHOD
C
  110	IF (J - I .GE. 11) GO TO 20
     	IF (I .EQ. II) GO TO 10
C
C  SHORT SEQUENCES (EXCEPT THE FIRST) ARE ORDERED BY A SINKING SORT
C
     	I = I - 1
  120	I = I + 1
     	IF (I .EQ. J) GO TO 100
     	IF (A(I) .LE. A(I+1)) GO TO 120
     	T = A(I+1)
     	U = B(I+1)
     	V = C(I+1)
     	K = I
  130	A(K+1) = A(K)
     	B(K+1) = B(K)
     	C(K+1) = C(K)
     	K = K - 1
     	IF (A(K) .GT. T) GO TO 130
     	A(K+1) = T
     	B(K+1) = U
     	C(K+1) = V
     	GO TO 120
     	END

     	SUBROUTINE IFMISC_FILTER( NPTS, YDAT, NDEG, NFIL )
C
C  PROGRAMMER: R.S. ZHOU	JUNE-1986
C
C	THIS ROUTINE IMPLEMENTS THE SAVITZKY-GOLAY SIMPLIFIED
C	LEAST SQUARE ALGORITHM FOR DATA SMOOTHING.
C
C	NPTS --- NUMBER OF DATA POINTS IN THE ARRAY.
C	YDAT --- ARRAY TO BE FILTERED.
C	NDEG --- POLYNOMIAL DEGREE (2,3,4,5).
C	NFIL --- FILTER LENGTH IN NUMBER OF DATA POINTS (5,7,9).
C
    	DIMENSION YDAT(1)
    	DIMENSION IC2(10,3), IC4(10,2)

    	DATA IC2 /
     + -3,  12,  17,  12,  -3,  35,   0,   0,   0,   0,
     + -2,   3,   6,   7,   6,   3,  -2,  21,   0,   0,
     +-21,  14,  39,  54,  59,  54,  39,  14, -21, 231 /

    	DATA IC4 /
     +  5,  -30,  75, 131,  75, -30,   5, 231,   0,   0,
     + 15,  -55,  30, 135, 179, 135,  30, -55,  15, 429 /

    	NF = NFIL
    	MF = NF/2 - 1
     	MM = NPTS
     	MID= 0

    	IF( NDEG.LT.4 ) THEN
     	    DO  J=1, MM-NF+1
     		SUM = 0.0
     		DO  K=1, NF
     		    SUM = SUM + YDAT(J+K-1)*IC2(K,MF)
     		ENDDO
     		MID = MID + 1
     		YDAT(MID) = SUM/IC2(K,MF)
     	    ENDDO
     	ELSE
     	    IF( NF.LT.7 ) NF = 7
    	    MF = MF - 1
     	    DO  J=1, MM-NF+1
     		SUM = 0.0
     		DO  K=1, NF
     		    SUM = SUM + YDAT(J+K-1)*IC4(K,MF)
     		ENDDO
     		MID = MID + 1
     		YDAT(MID) = SUM/IC4(K,MF)
     	    ENDDO
     	ENDIF

     	RETURN
    	END

     	SUBROUTINE IFMISC_DERIV( NPTS, YDAT, NDEG, NFIL, NDEV )
C
C  PROGRAMMER: R.S. ZHOU	JUNE-1986
C
C	THIS ROUTINE IMPLEMENTS THE SAVITZKY-GOLAY SIMPLIFIED
C	LEAST SQUARE ALGORITHM FOR DATA DIFFERENTIATION.
C
C	NPTS --- NUMBER OF DATA POINTS IN THE ARRAY.
C	YDAT --- INPUT DATA ARRAY TO BE FILTRATED.
C	NDEG --- POLYNOMIAL DEGREE (2,3,4).
C	NFIL --- FILTER LENGTH IN NUMBER OF DATA POINTS (5,7,9).
C	NDEV --- DERIVATIVE ORDER (1,2).
C
    	DIMENSION YDEV(10), YDAT(1)
    	DIMENSION ID12(10,3), ID14(10,3), ID22(10,3), ID24(10,3)

    	DATA ID12 /
     + -2,   -1,    0,    1,    2,   10,    0,    0,    0,    0,
     + -3,   -2,   -1,    0,    1,    2,    3,   28,    0,    0,
     + -4,   -3,   -2,   -1,    0,    1,    2,    3,    4,   60 /

    	DATA ID14 /
     +  1,   -8,    0,    8,   -1,   12,    0,    0,    0,    0,
     + 22,  -67,  -58,    0,   58,   67,  -22,  252,    0,    0,
     + 86, -142, -193, -126,    0,  129,  193,  142,  -86, 1188 /

    	DATA ID22 /
     +  2,   -1,   -2,   -1,   -2,    7,    0,    0,    0,    0,
     +  5,    0,   -3,   -4,   -3,    0,    5,   42,    0,    0,
     + 28,    7,   -8,  -17,  -20,  -17,   -8,    7,   28,  462 /

    	DATA ID24 /
     +    -3,   48,  -90,   48,   -3,    3,    0,    0,    0,    0,
     +  -117,  603, -171, -630, -171,  603, -117,   99,    0,    0,
     + -4158,12243, 4983,-6963,-12210,-6963,4983,12243,-4158, 4719 /

    	NO = NDEV
    	ND = NDEG
    	NF = NFIL
     	MM = NPTS
    	MF = NF/2 - 1
    	MID= 0
     	JK = 0
    	IF( NO.EQ.1 ) THEN
    	    IF( ND.LT.3 ) THEN
     		DO  J=1, MM-NF+1
     		    SUM = 0.0
     		    DO  K=1, NF
     			SUM = SUM + YDAT(J+K-1)*ID12(K,MF)
     		    ENDDO
     		    JK = JK + 1
     		    YDEV(JK) = SUM/ID12(K,MF)
     		    IF( JK.EQ.NF ) THEN
     			DO  I=1, NF
     			    MID = MID + 1
     			    YDAT(MID) = YDEV(I)
     			ENDDO
     			JK = 0
     		    ENDIF
     		ENDDO
    	    ELSE
     		DO  J=1, MM-NF+1
     		    SUM = 0.0
     		    DO  K=1, NF
     			SUM = SUM + YDAT(J+K-1)*ID14(K,MF)
     		    ENDDO
     		    JK = JK + 1
     		    YDEV(JK) = SUM/ID14(K,MF)
     		    IF( JK.EQ.NF ) THEN
     			DO  I=1, NF
     			    MID = MID + 1
     			    YDAT(MID) = YDEV(I)
     			ENDDO
     			JK = 0
     		    ENDIF
     		    MID = MID + 1
     		ENDDO
    	    ENDIF
     	ELSE
    	    IF( ND.LT.4 ) THEN
     		DO  J=1, MM-NF+1
     		    SUM = 0.0
     		    DO  K=1, NF
     			SUM = SUM + YDAT(J+K-1)*ID22(K,MF)
     		    ENDDO
     		    JK = JK + 1
     		    YDEV(JK) = SUM/ID22(K,MF)
     		    IF( JK.EQ.NF ) THEN
     			DO  I=1, NF
     			    MID = MID + 1
     			    YDAT(MID) = YDEV(I)
     			ENDDO
     			JK = 0
     		    ENDIF
     		ENDDO
     	    ELSE
     		DO  J=1, MM-NF+1
     		    SUM = 0.0
     		    DO  K=1, NF
     			SUM = SUM + YDAT(J+K-1)*ID24(K,MF)
     		    ENDDO
     		    JK = JK + 1
     		    YDEV(JK) = SUM/ID24(K,MF)
     		    IF( JK.EQ.NF ) THEN
     			DO  I=1, NF
     			    MID = MID + 1
     			    YDAT(MID) = YDEV(I)
     			ENDDO
     			JK = 0
     		    ENDIF
     		ENDDO
    	    ENDIF
     	ENDIF

     	IF( JK.GT.0 ) THEN
     	    DO  I=1, JK
     		MID = MID + 1
     		YDAT(MID) = YDEV(I)
     	    ENDDO
     	ENDIF

     	RETURN
    	END

     	SUBROUTINE IFMISC_MAP( MC, NC, NR, CL, AM, FF, XS, YS )
C
C	THIS ROUTINE CONTOURS A MAP.
C
C	MC --------- MAXIMUM NUBER OF COLUMNS IN 2-D ARRAY AM(I,J).
C	NC --------- NUMBER OF COLUMNS (X) TO BE PLOTTED.
C	NR --------- NUMBER OF    ROWS (Y) TO BE PLOTTED.
C	CL --------- CONTOUR LEVEL.
C	AM(I,J) ---- ORDINATE AT POINT X(I), Y(J)
C	             X(I),Y(J) ARE ASSUMED TO BE EQUALLY SPACED.
C
C	FF --------- PATTERN OF AREA-FILL IF FILL IS DESIRED.
C	XS,YS ------ WIDTH AND HEIGHT OF THE GRAPH IN INCHES.
C
     	DIMENSION         AM(1)
     	INTEGER           REC, DIM, FF, SPS
     	COMMON/MAP_XYBND/ MT, NT, XINC, YINC, XSIZ, YSIZ
     	COMMON/MAP_CAVIN/ DIM, IX, IY, IDX, IDY, ISS, SPS,
     +                    NP, N, CV, IS, IS0, IX0, IY0, DCP, IFF, KFF,
     +                    INX(8), INY(8), REC(1000), X(1024), Y(1024)
     	DATA              INX/ -1, -1,  0,  1,  1,  1,  0, -1 /
     	DATA              INY/  0,  1,  1,  1,  0, -1, -1, -1 /
C
C  INITIALIZE...
C
     	DIM  = IABS(MC)
     	SPS  = 0
     	IF( MC.LT.0 ) SPS = 1		!2-D SPLINE-SMOOTH ON
     	NT   = NC
     	MT   = NR
     	XSIZ = XS
     	YSIZ = YS
     	XINC = XSIZ/(NT-1)
     	YINC = YSIZ/(MT-1)
     	IFF  = FF
     	KFF  = 0
     	IF( IFF.NE.0 ) KFF = -1
     	CV   = CL
     	NP   = 0
     	ISS  = 0
     	MT1  = MT-1
     	NT1  = NT-1
C
C  CONTOUR THE MAP...
C
     	DO 110 I=1,MT1
     	IF ((AM(I)  -CV) .GE. 0.0) GO TO 110
     	IF ((AM(I+1)-CV) .LT. 0.0) GO TO 110
     	IX0=I+1
     	IX=I+1
     	IY0=1
     	IY=1
     	IS0=1
     	IS=1
     	IDX=-1
     	IDY=0
     	CALL IFMISC_MAP_TRACE( AM(1) )
  110	CONTINUE

     	J=MT-DIM
     	DO 120 I=1,NT1
     	J=J+DIM
     	IF ((AM(J)    -CV) .GE. 0.0) GO TO 120
     	IF ((AM(J+DIM)-CV) .LT. 0.0) GO TO 120
     	IX0=MT
     	IX=MT
     	IY0=I+1
     	IY=I+1
     	IDX=0
     	IDY=-1
     	IS0=7
     	IS=7
     	CALL IFMISC_MAP_TRACE( AM(1) )
  120	CONTINUE

     	J=MT+NT1*DIM+1
     	DO 130 I=1,MT1
     	J=J-1
     	IF ((AM(J)  -CV) .GE. 0.0) GO TO 130
     	IF ((AM(J-1)-CV) .LT. 0.0) GO TO 130
     	IX0=MT-I
     	IX=MT-I
     	IY0=NT
     	IY=NT
     	IDX=1
     	IDY=0
     	IS0=5
     	IS=5
     	CALL IFMISC_MAP_TRACE( AM(1) )
  130	CONTINUE

     	J=NT*DIM+1
     	DO 140 I=1,NT1
     	J=J-DIM
     	IF ((AM(J)    -CV) .GE. 0.0) GO TO 140
     	IF ((AM(J-DIM)-CV) .LT. 0.0) GO TO 140
     	IX0=1
     	IX=1
     	IY0=NT-I
     	IY=NT-I
     	IDX=0
     	IDY=1
     	IS0=3
     	IS=3
     	CALL IFMISC_MAP_TRACE( AM(1) )
  140	CONTINUE

     	ISS=1
     	L=0
     	DO 170 J=2,NT1
     	L=L+DIM
     	DO 160 I=1,MT1
     	L=L+1
     	IF ((AM(L)  -CV) .GE. 0.0) GO TO 160
     	IF ((AM(L+1)-CV) .LT. 0.0) GO TO 160
     	K=L+1
     	DO 150 ID=1,NP
     	IF (REC(ID) .EQ. K) GO TO 160
  150	CONTINUE

     	IX0=I+1
     	IX=I+1
     	IY0=J
     	IY=J
     	IDX=-1
     	IDY=0
     	IS0=1
     	IS=1
     	CALL IFMISC_MAP_TRACE( AM(1) )
  160	CONTINUE
  170	L=L-MT1
     	RETURN
     	END

     	SUBROUTINE IFMISC_MAP_TRACE( AM )

     	DIMENSION         AM(1)
     	INTEGER           REC, DIM, SPS
     	COMMON/MAP_XYBND/ MT, NT, XINC, YINC, XSIZ, YSIZ
     	COMMON/MAP_CAVIN/ DIM, IX, IY, IDX, IDY, ISS, SPS,
     +                    NP, N, CV, IS, IS0, IX0, IY0, DCP, IFF, KFF,
     +                    INX(8), INY(8), REC(1000), X(1024), Y(1024)

     	N=0
     	JY=DIM*(IY-1)+IX
     	MY=DIM*IDY+IDX+JY
  100	N=N+1
     	IF (N .GT.1020) GO TO 310
     	IF (IDX .GT. 0) GO TO 120
     	IF (IDX .LT. 0) GO TO 110
     	X(N)=FLOAT(IY-1)+FLOAT(IDY)*(AM(JY)-CV)/(AM(JY)-AM(DIM*IDY+JY))
     	Y(N)=FLOAT(IX-1)
     	GO TO 130
  110	NP=NP+1
     	REC(NP)=JY
  120	Y(N)=FLOAT(IX-1)+FLOAT(IDX)*(AM(JY)-CV)/(AM(JY)-AM(JY+IDX))
     	X(N)=FLOAT(IY-1)
  130	IS=IS+1
  140	IF (IS .LE. 8) GO TO 150
     	IS=IS-8
  150	IDX=INX(IS)
     	IDY=INY(IS)
     	IX2=IX+IDX
     	IY2=IY+IDY
     	IR=IDX*IDY
     	IF (ISS .EQ. 0) GO TO 160
     	IF ((IS.NE.IS0) .OR. (IY.NE.IY0) .OR. (IX.NE.IX0)) GO TO 170
     	N=N+1
     	X(N)=X(1)
     	Y(N)=Y(1)
     	GO TO 310
  160	IF( (IX2.EQ.0).OR.(IX2.GT.MT).OR.(IY2.EQ.0).OR.(IY2.GT.NT) )
     +	GO TO 310
  170	MY=DIM*IDY+IDX+JY
     	IF (IR .GT. 0) GO TO 200
     	IF (IR .LT. 0) GO TO 190
     	IF (CV .GT. AM(MY)) GO TO 100
     	IX=IX2
     	IY=IY2
  180	IS=IS+5
     	JY=MY
     	GO TO 140
  190	KY=JY+IDX
     	LY=MY-IDX
     	GO TO 210
  200	KY=MY-IDY
     	LY=JY+IDY
  210	DCP=(AM(JY)+AM(KY)+AM(LY)+AM(MY))*.25
     	IF (CV .LE. DCP) GO TO 220
     	CALL IFMISC_MAP_GETPT(JY,AM(1))
     	GO TO 130
  220	IF (IR .GE. 0) GO TO 230
     	IX=IX2
     	IDX=-IDX
     	CALL IFMISC_MAP_GETPT(KY,AM(1))
     	IY=IY2
     	IDY=-IDY
     	GO TO 240
  230	IY=IY2
     	IDY=-IDY
     	CALL IFMISC_MAP_GETPT(KY,AM(1))
     	IX=IX2
     	IDX=-IDX
  240	IF (CV .LE. AM(MY)) GO TO 180
     	CALL IFMISC_MAP_GETPT(MY,AM(1))
     	IF (IR .GE. 0) GO TO 250
     	IX=IX+IDX
     	IDX=-IDX
     	GO TO 260
  250	IY=IY+IDY
     	IDY=-IDY
  260	IF (CV .GT. AM(LY)) GO TO 290
     	IS=IS-1
     	JY=LY
     	GO TO 150
  290	CALL IFMISC_MAP_GETPT(LY,AM(1))
     	IF (IR .GE. 0) GO TO 300
     	IY=IY+IDY
     	GO TO 130
  300	IX=IX+IDX
     	GO TO 130
  310	DO 320 I=1, N
     	X(I) = X(I)*XINC
  320	Y(I) = Y(I)*YINC
C
C  PLOT THE CONTOUR LINE...
C  IF AREA-FILL IS ON, CLOSE AN OPEN LOOP BY ADDING A CORNER POINT...
C
     	IF( KFF.NE.0 .AND. ABS(X(1)-X(N)).GT.0.01
     +.AND. ABS(Y(1)-Y(N)).GT.0.01 ) THEN
     	    IF( X(1).EQ.0.0 .OR. X(N).EQ.0.0 ) THEN
     		X(N+1) = 0.0
     	    ELSE
     		X(N+1) = XSIZ
     	    ENDIF
     	    IF( Y(1).EQ.0.0 .OR. Y(N).EQ.0.0 ) THEN
     		Y(N+1) = 0.0
     	    ELSE
     		Y(N+1) = YSIZ
     	    ENDIF
     	    N = N + 1
     	ENDIF
     	X(N+1) = 0.0
     	Y(N+1) = 0.0
     	X(N+2) = 1.0
     	Y(N+2) = 1.0
     	IF( SPS.NE.0 ) N = -N
     	CALL IFLINE( IFF, 0.0, KFF, N, X, Y )
     	RETURN
     	END

     	SUBROUTINE IFMISC_MAP_GETPT(J,AM)

     	DIMENSION         AM(1)
     	INTEGER           REC, DIM, SPS
     	COMMON/MAP_CAVIN/ DIM, IX, IY, IDX, IDY, ISS, SPS,
     +                    NP, N, CV, IS, IS0, IX0, IY0, DCP, IFF, KFF,
     +                    INX(8), INY(8), REC(1000), X(1024), Y(1024)
     	N=N+1
     	B=AM(J)-DCP
     	IF (B.NE.0) GO TO 120
     	V=.5
     	GO TO 130
  120	V=.5*(AM(J)-CV)/B
  130	Y(N)=FLOAT(IX-1)+FLOAT(IDX)*V
     	X(N)=FLOAT(IY-1)+FLOAT(IDY)*V

     	RETURN
     	END

     	SUBROUTINE IFMISC_MAP_EXPAND(JS,MC,NC,NR,AM)
C
C	THIS ROUTINE EXPANDS A CONTOUR MAP BY A SPLINE FUNCTION.
C
C	JS ---- GRID EXPANSION FACTOR.
C	NR ---- NUMBER OF ROWS.
C	NC ---- NUMBER OF COLUMNS.
C	AM ---- THE MAP ARRAY.
C
     	DIMENSION          AM(MC,MC)
     	COMMON/MAP_EXPAND/ Q(101), Z(101), R(101)
C
C----INITIALIZE
C
     	ISP = JS
     	N   = NC
     	M   = NR
    4	MM = M + ISP*(M-1)
     	NN = N + ISP*(N-1)
     	IF( MM.GT.MC .OR. NN.GT.MC ) THEN
     	    ISP = ISP - 1
     	    GO TO 4
     	ENDIF
     	SIP = ISP + 1
C
C----INTERPOLATE FOR ROWS
C
    8	IF(ISP.EQ.0) RETURN
     	LN=N-1
     	IM=M+1
     	DO 30 IT=1,M
     	I=IM-IT
     	DO 10 J=1,N
     	Q(J) = AM(I,J)
   10	CONTINUE
     	CALL IFMISC_MAP_SPLINE(N)
     	II = I+ISP*(I-1)
     	DO 20 J=1,LN
     	JJ = J+ISP*(J-1)
     	AM(II,JJ) = Q(J)
     	DO 20 K=1,ISP
     	RAT = FLOAT(K)/SIP
     	QAT = 1. - RAT
     	JK=JJ+K
     	AM(II,JK) = QAT*(Z(J)*QAT**2+R(J))+RAT*(Z(J+1)*RAT**2+R(J+1))
   20	CONTINUE
     	AM(II,NN) = Q(N)
   30	CONTINUE
C
C----INTERPOLATE FOR COLUMNS
C
     	LM=M-1
     	DO 60 JJ=1,NN
     	DO 40 I=1,M
     	II = I+ISP*(I-1)
     	Q(I) = AM(II,JJ)
   40	CONTINUE
     	CALL IFMISC_MAP_SPLINE(M)
     	DO 50 I=1,LM
     	II = I+ISP*(I-1)
     	DO 50 K=1,ISP
     	RAT = FLOAT(K)/SIP
     	QAT = 1. - RAT
     	IK=II+K
     	AM(IK,JJ) = QAT*(Z(I)*QAT**2+R(I))+RAT*(Z(I+1)*RAT**2+R(I+1))
   50	CONTINUE
   60	CONTINUE

     	NC = NN
     	NR = MM
     	RETURN
     	END

     	SUBROUTINE IFMISC_MAP_SPLINE(M)

     	COMMON/MAP_EXPAND/ Y(101), Z(101), R(101)
     	DIMENSION          B(101), E(101), A(101,3)
     	EQUIVALENCE        (E,Z)

     	MM=M-1
     	DO 2 K=1,MM
    2	E(K) = Y(K+1)-Y(K)
     	DO 3 K=2,MM
    3	B(K)=E(K)-E(K-1)
     	A(1,2)=-2.
     	A(1,3)=1.
     	A(2,2)=1.
     	A(2,3)=0.
     	DO 4 K=3,MM
     	A(K,2)=2./3. - A(K-1,3)/6.
     	B(K)=B(K)-B(K-1)/6.
     	A(K,3)=1./(6.*A(K,2))
    4	B(K)=B(K)/A(K,2)
     	A(M,1)=2. + A(M-2,3)
     	A(M,2) = -A(M,1)*A(M-1,3) -1.
     	B(M)=B(M-2)-A(M,1)*B(M-1)
     	Z(M)=B(M)/A(M,2)
     	MN=M-2
     	DO 6 I=1,MN
     	K=M-I
    6	Z(K)=B(K)-A(K,3)*Z(K+1)
     	Z(1)=-A(1,2)*Z(2)-A(1,3)*Z(3)
     	DO 7 K=1,M
     	Z(K) = Z(K)/6.
     	R(K) = Y(K) - Z(K)
    7	CONTINUE
     	RETURN
     	END

