PROGRAM PLOTMIETAB
C         PLOTMIETAB plots the single scatter phase function from a
C       Mie table produced by cloudprp.f.  The Mie table contains the
C       Legendre polynomial expansion of the phase function for several
C       effective radii.  
C
      INTEGER MAXNTAB, MAXLEG, MAXANG
      PARAMETER (MAXNTAB=300, MAXLEG=2000, MAXANG=1801)
      INTEGER NRETAB, NANGLE, NOUTREFF, NLEG(MAXNTAB)
      INTEGER I, J, K, L
      REAL    PI, PL, PL1, PL2, SCAT, MU
      REAL    ALPHA, WAVELEN, SRETAB, ERETAB
      REAL    OUTREFF(MAXNTAB), REFF(MAXNTAB)
      REAL    LEGEN(MAXLEG+1,MAXNTAB)
      REAL    EXTINCT(MAXNTAB), ALBEDO(MAXNTAB)
      REAL    ANGLE(MAXANG), PHASE(MAXANG,MAXNTAB)
      COMPLEX RINDEX
      CHARACTER  MIEFILE*72, PLOTFILE*72, distflag*1


      WRITE (*,*) 'Mie table file name'
      READ (*,'(A)') MIEFILE

      WRITE (*,'(1X,A)') 'Number of angles (<0 for input angles) : '
      READ (*,*) NANGLE
      IF (ABS(NANGLE) .GT. MAXANG)  STOP 'Too many angles'
      IF (NANGLE .LE. 0) THEN
        NANGLE = MAX(1, IABS(NANGLE))
        WRITE (*,*) 'Input the angles (degrees)'
        DO J = 1, NANGLE
          WRITE (*,'(1X,I3,A)') J, ' : '
          READ (*,*) ANGLE(J)
        ENDDO
      ELSE
        DO J = 1, NANGLE
          ANGLE(J) = 180.0* FLOAT(J-1) /(NANGLE-1)
        ENDDO
      ENDIF


      WRITE (*,'(1X,A)') 'Number of output effective radii : '
      READ (*,*) NOUTREFF
      WRITE (*,*) 'Input the effective radii (micron) : '
      READ (*,*) (OUTREFF(K), K=1,NOUTREFF)

      WRITE (*,*) 'Plotting output file name'
      READ (*,'(A)') PLOTFILE


C         Read in the Mie table and write out some information
      CALL READ_MIE_TABLE (MIEFILE, NRETAB, SRETAB, ERETAB, REFF,
     .                     WAVELEN, RINDEX, ALPHA, 
     .                 EXTINCT, ALBEDO, MAXLEG, NLEG, LEGEN, distflag)


      PI = ACOS(-1.0)
C        Loop over each effective radius in Mie table (I), pulling out ones
C        we want to output (K)
      K = 1
      DO I = 1, NRETAB
        IF (REFF(I) .EQ. OUTREFF(K)) THEN
C             Sum the Legendre series for each angle in plot
          DO J = 1, NANGLE
            MU = COS(ANGLE(J)*PI/180.0)
            SCAT = 0.0
            PL1 = 1.0
            PL = 1.0
            DO L = 0, NLEG(I)
              IF (L .GT. 0)  PL = (2*L-1)*MU*PL1/L - (L-1)*PL2/L
              SCAT = SCAT + LEGEN(L+1,I)*PL
              PL2 = PL1
              PL1 = PL
            ENDDO
            PHASE(J,K) = SCAT
          ENDDO
          K = K + 1
        ENDIF
      ENDDO
      NOUTREFF = K-1


C         Output the phase functions
      OPEN (UNIT=1, FILE=PLOTFILE, STATUS='UNKNOWN')
      WRITE (1,'(A)') '!  Mie table phase functions'
      WRITE (1,'(A,E13.6,A)') '! ', WAVELEN, '  wavelength (micron)'
      WRITE (1,'(A,2(1X,E13.6),A)') 
     .  '! ', RINDEX, '  index of refraction'
      write (1,'(a,a1,a)')'!', distflag, ' distribution type'
      WRITE (1,'(A,F7.5,A)') 
     .  '! ', ALPHA, '  distribution shape parameter'
      WRITE (1,'(A,1X,I3,2(1X,F8.4),A)') '! ', NRETAB, SRETAB, ERETAB, 
     .  '  number, starting, ending effective radius'
      WRITE (1,'(A)') 
     .  '! Angle  cos(angle)  Phase functions for effective radii (um)'
      WRITE (1,'(A,20(6X,F6.2))') '!               ', 
     .  (OUTREFF(K), K=1,NOUTREFF)
      DO J = 1, NANGLE
          WRITE (1,'(1X,F7.2,1X,F9.6,20(E12.4))')
     .             ANGLE(J), COS(ANGLE(J)*PI/180.0),
     .             (PHASE(J,K), K=1,NOUTREFF)
      ENDDO
      CLOSE (1)

      END





      SUBROUTINE READ_MIE_TABLE (MIEFILE, NRETAB, SRETAB, ERETAB, 
     .                           REFF, WAVELEN, RINDEX, ALPHA, 
     .                 EXTINCT, ALBEDO, MAXLEG, NLEG, LEGEN, distflag)
C       Reads a table of Mie scattering properties as a function of 
C     effective radius.  
C     INTEGER NRETAB, MAXLEG, NLEG(NRETAB)
      INTEGER NRETAB, MAXLEG, NLEG(*)
      REAL    SRETAB, ERETAB, WAVELEN, ALPHA, REFF(*)
      REAL    EXTINCT(*), ALBEDO(*), LEGEN(MAXLEG+1,*)
      COMPLEX RINDEX
      CHARACTER MIEFILE*72, distflag*1
      INTEGER I, L
      REAL    MR, MI

      OPEN (UNIT=3, FILE=MIEFILE, STATUS='OLD')
      READ (3,*)
      READ (3,*) WAVELEN
      READ (3,*) MR, MI
      RINDEX = CMPLX(MR,MI)
      read (3,'(a)') distflag
      READ (3,*) ALPHA
      READ (3,*) NRETAB, SRETAB, ERETAB
      DO I = 1, NRETAB
        READ (3,*) REFF(I), EXTINCT(I), ALBEDO(I), NLEG(I)
        READ (3,*) (LEGEN(L,I), L=1,NLEG(I)+1) 
      ENDDO
      CLOSE (3)
      RETURN
      END