PROGRAM PRINTCOF
C
C     TITLE: PRINT GALERKIN COEFFICIENTS FILE.
C     ----------------------------------------
C
C        THIS PROGRAM WILL READ TABLES FROM A FORMATFREE FILE PRODUCED BY
C     SUBROUTINE 'GNRT'.  THESE WILL THEN BE PRINTED IN A CONVENIENT FORM
C     IN DECIMAL.  THE USER WILL BE ASKED FOR TWO QUANTITIES.
C
C     INPUT
C          FNAME: THE NAME FOR THE FORMATFREE FILE OF COEFFICIENTS.
C          DEG: THE MAXIMUM DEGEE OF THE TABLES FROM FILE 'FNAME' THAT
C              ARE TO BE PRINTED.  THEY WILL BE PRINTED IN DECIMAL, WITH
C              THE INDICES OF THE ASSOCIATED SPHERICAL HARMONICS.
C              *** NOTE ***  THERE IS AN UPPER LIMIT ON 'DEG',
C              GIVEN AS 'MAXDEG' IN THE BELOW PARAMETER STATEMENT.
C
C     OUTPUT
C          FNAMED: THE NAME OF THE OUTPUT FILE OF DECIMAL COEFFICIENTS.
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      INTEGER DEG,ORDER,DEGREE,FILEIN,FILOUT
      DOUBLE PRECISION  KERMAT
      CHARACTER FNAME*50,FNAMED*50
      PARAMETER(MAXDEG=10,MAXORD=(MAXDEG+1)**2)
      DIMENSION KERMAT(MAXORD,MAXORD)
C
C     ****** INPUT/OUTPUT ****************************************************
C     FOLLOWING ARE THE I/O FILE NUMBERS BEING USED.  THESE MAY NEED TO
C     BE CHANGED FOR YOUR COMPUTER.
      DATA FILEIN/8/, FILOUT/9/
C     ************************************************************************
C
C     OBTAIN INPUT INFORMATION.
      PRINT *, ' GIVE THE NAME OF YOUR INPUT FILE OF FORMATFREE'
      PRINT *, ' COEFFICIENTS.'
      READ(*,'(A)') FNAME
      PRINT *, ' GIVE A NAME TO BE USED TO DENOTE THE OUTPUT FILE'
      PRINT *, ' OF DECIMAL COEFFICIENTS.'
      READ(*,'(A)') FNAMED
      PRINT *, ' GIVE THE MAXIMUM DEGREE OF THE COEFFICIENTS THAT YOU'
      PRINT *, ' WANT TO HAVE OUTPUT.'
      READ *, DEG
C
      L=LEN(FNAME)
      OPEN(FILEIN,FILE=FNAME(1:L),FORM='UNFORMATTED')
      L=LEN(FNAMED)
      OPEN(FILOUT,FILE=FNAMED(1:L))
C
C     MAIN LOOP.
      DO 30 K=1,DEG
        READ(FILEIN)DEGREE,NINNER,NOUTER
        WRITE(FILOUT,1000) DEGREE,NINNER,NOUTER
1000    FORMAT(//,' GALERKIN COEFFICIENTS',//,
     *     ' DEGREE=',I2,5X,'N_INNER=',I2,5X,'N_OUTER=',I2,//,
     *     ' THE INDICES OF A SPHERICAL HARMONIC ARE (N,M,L),'/,
     *     ' THE MEANING OF WHICH IS GIVEN IN FUNCTION LOCATN.',/,
     *     ' THE SUFFIXES R AND C USED BELOW ARE TO INDICATE WHETHER',/,
     *     ' THE ROW OR COLUMN IS BEING REFERENCED.',//)
        ORDER=(DEGREE+1)**2
        READ(FILEIN)((KERMAT(I,J),J=1,ORDER),I=1,ORDER)
        DO 20 MR=0,DEGREE
        DO 20 NR=MR,DEGREE
        DO 20 LR=0,1
          IF((LR .EQ. 1) .AND. (MR .EQ. 0)) GO TO 20
          LOCR=LOCATN(NR,MR,LR,DEGREE)
          WRITE(FILOUT,1001) LOCR,NR,MR,LR
1001      FORMAT(/,' ROW ',I2,3X,'(NR,MR,LR)=',3I3,/,
     *          7X,'COL   NC   MC  LC       ELEMENT')
          DO 10 MC=0,DEGREE
          DO 10 NC=MC,DEGREE
          DO 10 LC=0,1
            IF((LC .EQ. 1) .AND. (MC .EQ. 0)) GO TO 10
            LOCC=LOCATN(NC,MC,LC,DEGREE)
            WRITE(FILOUT,1002) LOCC,NC,MC,LC,KERMAT(LOCR,LOCC)
1002        FORMAT(7X,I3,3X,I2,3X,I2,3X,I1,1PD20.10)
10          CONTINUE
20        CONTINUE
30      CONTINUE
      CLOSE(FILEIN)
      CLOSE(FILOUT)
      STOP
      END

      FUNCTION LOCATN(N,M,L,NDEG)
C     ---------------
C
C     GIVES THE INDEX OF THE SPHERICAL HARMONIC BASIS FUNCTION
C     WITH INDICES (N,M,L).  FOR Q=(X,Y,Z) ON THE UNIT SPHERE, WRITE
C        Q = (SIN(THETA)*COS(PHI),SIN(THETA)*SIN(PHI),COS(THETA))
C     FOR M=0, THE SPHERICAL HARMONIC CORRESPONDING TO (N,M,L) IS
C
C        H(Q) = P (Z).
C                N
C
C     FOR M>0 AND L=0, THE SPHERICAL HARMONIC IS
C
C                M
C        H(Q) = P (Z)*COS(M*PHI)
C                N
C
C     FOR M>0 AND L=1, THE SPHERICAL HARMONIC IS
C
C                M
C        H(Q) = P (Z)*SIN(M*PHI)
C                N
C
C                              M
C     THE FUNCTIONS P (Z) AND P (Z) ARE THE LEGENDRE POLYNOMIALS AND
C                    N         N
C     THE ASSOCIATED LEGENDRE FUNCTIONS, DEFINED IN SUBROUTINE 'LEGEND'.
C
      IF(M .EQ. 0) THEN
          LOCATN=N+1
        ELSE
          LOCATN=NDEG*(2*M-1)-M*(M-3)+2*(N-M)+L
        END IF
      RETURN
      END