!!$---------------------------------------------------------------------
SUBROUTINE AGNPREP(MJTOT,DUM3,FnuSTAR,MJnu1,Inu1,LAMBDA,WT3,WT4,&
     CTHETA,Z,NK0,NRAD,NWAV,NQ,MXTH,SFLAG)
!!$---------------------------------------------------------------------
!!$  This subroutine adds the "AGN mode" to the 2-Dust code.
!!$  
!!$  In collaboration with R. Szczerba, Z. Loska, & B. Czerny at the
!!$  N. Copernicus Astronomical Institute of Poland.  
!!$                                         June 2003, Toshiya Ueta, ROB
!!$---------------------------------------------------------------------
!!$  AGN is represented here as an artificial central point source that
!!$  consists of the following two components;
!!$   (1) the disk component
!!$       this has the accretion disk spectrum (log(F_dsk)) whosse 
!!$       radiation is anisotropic and decreases with the cosine of the
!!$       inclination angle
!!$   (2) the spherical component
!!$       this has a power law spectrum (log(F_sph))
!!$
!!$  This subroutine needs "sourceflux.dat" with the following format
!!$   line 1-2: blank
!!$   line 3  : Rdsk, Rsph in AU
!!$   line 4-8: blank
!!$   line 9- : log(lam[um]), log(I_dsk[cgs]), log(I_sph[cgs])
!!$---------------------------------------------------------------------
  USE DEFVARS
!!$
  IMPLICIT NONE
!!$
  INTEGER, INTENT(IN)                      :: NRAD,NWAV,NQ,MXTH,SFLAG
  REAL(PRC), DIMENSION(NWAV), INTENT(IN)   :: LAMBDA,WT3,WT4
  REAL(PRC), DIMENSION(0:NQ), INTENT(IN)   :: CTHETA
  REAL(PRC), DIMENSION(NQ), INTENT(IN)     :: Z
  INTEGER, DIMENSION(0:NRAD+1), INTENT(IN) :: NK0
!!$
  REAL(PRC), INTENT(INOUT)                   :: MJTOT,DUM3
  REAL(PRC), DIMENSION(MXTH,NWAV,0:NRAD+1,0:2*NQ+1), INTENT(INOUT) :: &
       Inu1
  REAL(PRC), DIMENSION(NWAV,0:NRAD+1,0:2*NQ+1), INTENT(INOUT)      :: &
       FnuSTAR,MJnu1
!!$  
  CHARACTER(30) :: DUMMY
!!$
  INTEGER :: IOFLAG,IERROR,I,II,M,K,L,NPTS,NQ2
!!$
  REAL(PRC) :: DUM1,DUM2,F1,F2,RDUM,RDSK,RSPH,LS,DOMEGA,AA,BB
!!$
  REAL(PRC), ALLOCATABLE, DIMENSION(:,:)     :: DUMARR1,DUMARR2,DUMARR3
!!$
  NQ2 = 2*NQ+1
!!$
  OPEN(UNIT=30,FILE='datafiles.dat',STATUS='OLD',IOSTAT=IERROR)
  READ(30,*) IOFLAG
  CLOSE(30)
  IF (IOFLAG == 0) THEN
     WRITE(*,'(" Invoking the AGN mode... ")')
     WRITE(*,'("   ... reading sourceflux.dat ")')
     WRITE(*,'(" ")')
  END IF
  ALLOCATE(DUMARR1(3,NDUM),STAT=IERROR)
  DUMARR1 = 0.0_PRC
  OPEN(UNIT=30,FILE='sourceflux.dat',STATUS='OLD',IOSTAT=IERROR)
  READ(30,'(A30)') DUMMY
  READ(30,'(A30)') DUMMY
  READ(30,*)       RDSK,RSPH ! in AU
  RDSK  = RDSK * AU ! in cm
  RSPH  = RSPH * AU ! in cm
  READ(30,'(A30)') DUMMY
  READ(30,'(A30)') DUMMY
  READ(30,'(A30)') DUMMY
  READ(30,'(A30)') DUMMY
  READ(30,'(A30)') DUMMY
  I = 0
  DO II=1,NDUM
     READ(30,*,IOSTAT=IERROR) DUM1,DUM2,F1
     IF (IERROR < 0.0_PRC) EXIT
     I = I + 1
     DUMARR1(1,II) = 10.0_PRC**DUM1 ! lambda in um
     F2 = LOG10(DUMARR1(1,II)*DUMARR1(1,II)/CSP)
     IF (DUM2 /= 0.0_PRC) THEN
        DUMARR1(2,II) = DUM2 + F2 ! I_disk in log(I_nu)
     ELSE
        DUMARR1(2,II) = 0.0_PRC
     END IF
     IF (F1 /= 0.0_PRC) THEN
        DUMARR1(3,II) = F1   + F2 ! I_sph  in log(I_nu)
     ELSE
        DUMARR1(3,II) = 0.0_PRC
     END IF
  END DO
  CLOSE(UNIT=30)
  NPTS = I
!!$
!!$  Copy only the filled part of the DUMARR1 into DUMARR3
!!$   and invert the order to ascending order for interpolation
!!$
  ALLOCATE(DUMARR3(4,NPTS), STAT=IERROR)
  DUMARR3 = 0.0_PRC
  DO M=1,3
     DO II=1,NPTS
        DUMARR3(M,II) = DUMARR1(M,NPTS-II+1)
     END DO
  END DO
  DEALLOCATE(DUMARR1,STAT=IERROR) ! no longer needed
  
!!$  DO II=1,NPTS
!!$     write(*,*) DUMARR3(1,ii),DUMARR3(2,ii),DUMARR3(3,ii)
!!$     write(51,*) DUMARR3(1,ii),DUMARR3(2,ii),DUMARR3(3,ii)
!!$  END DO
  
!!$---------------------------------------------------------------------  
!!$  Interpolate log(l_dsk[cgs]) and log(l_sph[cgs]) values
!!$  at the wavelengths specified in the specgrid.dat
!!$---------------------------------------------------------------------  
  ALLOCATE(DUMARR2(3,NWAV),STAT=IERROR)
  DUMARR2 = 0.0_PRC
  RDUM = 1.0E+30_PRC
  DO M = 2,3
     CALL SPLINE(DUMARR3(1,:),DUMARR3(M,:),NPTS,RDUM,RDUM,&
          DUMARR3(4,:))
     DO II=1,NWAV
        DUMARR2(1,II) = LAMBDA(II)
        CALL SPLINT(DUMARR3(1,:),DUMARR3(M,:),DUMARR3(4,:),NPTS,&
             LAMBDA(II),DUMARR2(M,II))
     END DO
  END DO
  
  OPEN(UNIT=30,FILE='sourcefluxout.dat',STATUS='unknown',IOSTAT=IERROR)
  WRITE(30,'(2ES22.13)') RDSK,RSPH ! in cm
  DO II = 1,NWAV
!!$        write(*,'(3ES22.13)') DUMARR2(1,II),DUMARR2(2,II),DUMARR2(3,II)
!!$        write(52,'(3ES22.13)') DUMARR2(1,II),DUMARR2(2,II),DUMARR2(3,II)
     DUM1 = DUMARR2(2,II)
     DUM2 = DUMARR2(3,II)
     IF (DUM1 /= 0.0_PRC) THEN
        DUMARR2(2,II) = (10.0_PRC**DUMARR2(2,II))* &
             (RDSK/RSPH)*(RDSK/RSPH)
     ELSE
        DUMARR2(2,II) = 0.0_PRC
     END IF
     IF (DUM2/= 0.0_PRC) THEN
        DUMARR2(3,II) = 10.0_PRC**DUMARR2(3,II)
     ELSE
        DUMARR2(3,II) = 0.0_PRC
     END IF
!!$        WRITE(*,'(3ES22.13)') DUMARR2(1,II),DUMARR2(2,II),DUMARR2(3,II)
        WRITE(30,'(3ES22.13)') DUMARR2(1,II),DUMARR2(2,II),DUMARR2(3,II)
  END DO
  CLOSE(UNIT=30)
  DEALLOCATE(DUMARR3,STAT=IERROR)
!!$---------------------------------------------------------------------  
!!$  Now input lat-dependent source rad field is set up
!!$    DUMARR2(1,NWAV) = lambda in um
!!$    DUMARR2(2,NWAV) = I_nu,dsk scaled to RSPH in cgs
!!$    DUMARR2(3,NWAV) = I_nu,sph in cgs
!!$  and these values along with R_dsk and R_sph are written in a file
!!$  called "sourcefluxout.dat".      
!!$---------------------------------------------------------------------  
!!$  Now integrate the source rad field to get luminosity
!!$---------------------------------------------------------------------  
  DUM1   = STARRADIUS/RMIN
  DOMEGA = PI * DUM1 * DUM1
  DUM3   = 0.0_PRC
  DO II = 1,NWAV
!!$     First set FnuSTAR and MJnu1 along the axis
!!$       Input intensity is at the surface of the source, and
!!$       thus need to dilute by "DOMEGA"
     LS = (CTHETA(0) * DUMARR2(2,II) + DUMARR2(3,II)) * DOMEGA 
     FnuSTAR(II,0,0)   = LS
     FnuSTAR(II,0,NQ2) = FnuSTAR(II,0,0)
     MJnu1(II,0,0)     = LS/FOURPI ! meanintensity
     MJnu1(II,0,NQ2)   = MJnu1(II,0,0)
!!$     Then set FnuSTAR and MJnu1 for THETA(NQ2) angles
!!$       Input intensity is at the surface of the source, and
!!$       thus need to be diluted by "DOMEGA"
     BB = 0.0_PRC
     AA = 0.0_PRC
     DO M=1,NQ
        LS = CTHETA(M) * DUMARR2(2,II) + DUMARR2(3,II)
!!$        write(*,*) II,M,LS,CTHETA(M),DUMARR2(2,II),DUMARR2(3,II)
        BB = BB + LS * Z(M)  
        LS = LS * DOMEGA
        FnuSTAR(II,0,M)     = LS
        FnuSTAR(II,0,NQ2-M) = FnuSTAR(II,0,M)
        LS                  = LS / FOURPI
        AA                  = AA + LS * Z(M)
        MJnu1(II,0,M)       = LS
        MJnu1(II,0,NQ2-M)   = MJnu1(II,0,M)
!!$     Initially anisotropic radiation field not supported yet.
!!$     So, make the initial field isotropic
        IF (SFLAG == 1) THEN
           DO L=0,NRAD+1
              DO K=1,NK0(L)
                 Inu1(K,II,0,M)     = MJnu1(II,0,M)
                 Inu1(K,II,0,NQ2-M) = Inu1(K,II,0,M)
              END DO
           END DO
        END IF
     ENDDO    
!!$     Below, DUM1 is B_nu (theta, phi) being integrated over the 
!!$     entire wavelength range.  The factor "4 PI" consists of "2 PI"
!!$     for the phi angle integration (azimuthal symmetry) and "2" for 
!!$     the fact that the radiation field is symmetric w.r.t. the 
!!$     equatorial plane (we calculated above for the half sphere).
!!$
!!$     MJTOT will be needed to estimate how much radiation is due to
!!$     star and how much is due to dust.  We do exactly the same 
!!$     integration for MJTOT as DUM1 but then MJTOT is "mean"intensity
!!$     averaged over "4 PI" radians.  So, no factor "4 PI" here.
     DUM3  = DUM3  + WT4(II) * BB * FOURPI ! BB for half-sphere
     MJTOT = MJTOT + WT3(II) * AA
  END DO
  DEALLOCATE(DUMARR2,STAT=IERROR)  
END SUBROUTINE AGNPREP
