SUBROUTINE RADT

#include "com2d.h"
#include "comphot.h"

C     REDUCED FLUXES AT ALL 
C     WAVELENGTHS FOR ALTITUDE ZP WITHOUT RADIATIVE TRANSFER.
C
      LOGICAL LGRIDP

	DO 2000 IJ=1,L$

	CALL FILL(IJ,Z$)

C     DETERMINE DIURNAL AVERAGING COEFFICIENTS.
      IF(.NOT.LDAV)GO TO 100


C     PARAMETERS FOR TABLE LOOKUP OF DIURNAL AVERAGE FACTORS. THE
C     DECLINATION PARAMETERS DEFINE THE TABLE DECLINATIONS BETWEEN
C     WHICH INTERPOLATION IS PERFORMED FOR POSITIVE LATITUDES. FOR
C     NEGATIVE LATITUDES, USE THE CORRESPONDING NEGATIVE DECLINATIONS.
      IF(DECD.GE.0.0)THEN
        I=2
        JI=(IDA+1)/2
      ELSE
        I=(IDA+1)/2+1
        JI=IDA
      END IF
      DO 200 ID1=I,JI
      IF(DECD.GT.DECA(ID1))GO TO 220
  200 CONTINUE
  220 ID0=ID1-1
      DECC1=DECA(ID0)
      DECC2=DECA(ID1)
      R=DECC2-DECC1
      DDEC1=ABS((DECD-DECC1)/R)
      DDEC2=ABS((DECD-DECC2)/R)

C     LATITUDE PARAMETERS. REFERENCE LATITUDES ARE ALL POSITIVE. IP
C     IS THE INDEX OF THE REFERENCE LATITUDE .GE. ABS(LATP). LGRIDP
C     FLAGS WHETHER CURRENT LATITUDE MATCHES A REFERENCE LATITUDE.
      ALATP=ABS(LATP)
      DO 240 IP=1,ILATR
      IF(ALATP.LT.LATR(IP))GO TO 260
  240 CONTINUE
  260 IF(ABS(ALATP-LATR(IP-1)).LT.1.e-5)THEN
        IP=IP-1
        LGRIDP=.TRUE.
      ELSE
        LGRIDP=.FALSE.
      END IF
      IF(LATP.GE.0.)THEN
C       INDICES FOR POSITIVE LATITUDES.
        ID=ID0
        ID1=ID+1
      ELSE
C       INDICES FOR NEGATIVE LATITUDES.
        ID=IDA-ID0+1
        ID1=ID-1
      END IF

  100 CONTINUE

C     CALCULATION WITHOUT RADIATIVE TRANSFER (BEER'S LAW).
C     LOOP OVER WAVELENGTHS.
      DO 300 JL=1,IL$

	DO 1000 IK=1,Z$

         FLUXMS=FLUXMULT(JL,IJ,IK)
         TAUZ=TAUMULT(JL,IJ,IK)
C OPTICAL DEPTH AND REDUCED FLUX COMPUTED IN MULTSC

      IF(LDAV.AND.TAUZ.LT.50.0)THEN

C       DIURNAL AVERAGING PROCEDURE BASED ON HASH TABLE ITAUA OF OPTICAL
C       DEPTHS. OPTICAL DEPTHS AND SOLAR FLUX RATIOS IN LOG FORM.
C       THE COMPUTED DIURNAL AVERAGE FACTOR FAC IS ALSO IN LOG FORM.
C       TWO-PARAMETER INTERPOLATION ON DECLINATION AND OPTICAL DEPTH.
        I=10.0*TAUZ+1.0
        IT=ITAUA(I)
        IT1=IT+1
        S1=FAVG(IT,ID,IP)
        S2=FAVG(IT,ID1,IP)
        S3=FAVG(IT1,ID,IP)
        S4=FAVG(IT1,ID1,IP)
        TAU1=TAUA(IT)
        TAU2=TAUA(IT1)
        DTAU=TAU2-TAU1
        R=S1*DDEC2+S2*DDEC1
        S=S3*DDEC2+S4*DDEC1
        TAUL=1.e-20
        IF(TAUZ.NE.0.)TAUL=aLOG(TAUZ)
        FAC=(R*ABS(TAUL-TAU2)+S*ABS(TAUL-TAU1))/DTAU
c	if(ij.eq.18 .and. jl.eq.38)then
c	print *,' ij=',ij,' ik=',ik,
c     *  ' r=',r,' s=',s,' s1=',s1,' s2=',s2,' s3=',s3,
c     *  ' s4=',s4,' ddec1=',ddec1,' ddec2=',ddec2
c	print *,' fac=',fac,' taul=',taul,' tau1=',tau1,
c     *  ' tau2=',tau2,' dtau=',dtau,' it=',it,' id=',
c     *  id,' ip=',ip
c	endif

        IF(.NOT.LGRIDP)THEN
C         LATITUDE INTERPOLATION IF NOT ON REFERENCE LATITUDE.
          S11=FAVG(IT,ID,IP-1)
          S21=FAVG(IT,ID1,IP-1)
          S31=FAVG(IT1,ID,IP-1)
          S41=FAVG(IT1,ID1,IP-1)
          R=S11*DDEC2+S21*DDEC1
          S=S31*DDEC2+S41*DDEC1
          FAC1=(R*ABS(TAUL-TAU2)+S*ABS(TAUL-TAU1))/DTAU
          FAC=FAC+(ALATP-LATR(IP))*(FAC-FAC1)/(LATR(IP)-LATR(IP-1))
        END IF
      ELSE
        FAC=0.
      END IF

	FAC1=-FAC
C This is the correct expression when multiple scattering is used

c	if(ij.eq.18 .and. jl.eq.38)then
c	print *,' ij=',ij,' ik=',ik,' fac=',fac,' fac1=',fac1
c	endif

C     NOW COMPUTE EXP(FAC1), WHERE FAC=LN(DIURNAL
C     AVERAGE FACTOR).
      IF(FAC1.GT.50.)THEN
        EFAC=0.
      ELSE IF(FAC1.GT.1.e-4)THEN
        EFAC=EXP(-FAC1)
      ELSE
        EFAC=1.-FAC1
      END IF

C     HERE STORE THE REDUCED FLUX.
	F=FLUXMS*EFAC

	IF(JL.EQ.19 .AND. F.GE.1.E-20)THEN
	   O2JINT(IJ,IK)=O2JINT(IJ,IK)*EFAC
C SCALE J(O2) TO DIURNAL AVERAGE VALUES AT ABOUT 200 NM - 4/7/95
        ENDIF

	IF(JL.EQ.19 .AND. F.LT.1.E-20)THEN
	   O2JINT(IJ,IK)=0.0E0
C SCALE J(O2) TO DIURNAL AVERAGE VALUES AT ABOUT 200 NM - 4/7/95
        ENDIF

      IF(F.LT.1.e-20)F=0.
      RFLUX(JL,IJ,IK)=F

       IF(TFD(IJ) .GT. 1.E-5)RFLUX(JL,IJ,IK)=RFLUX(JL,IJ,IK)/TFD(IJ)
C  CORRECT BACK TO DAYTIME AVERAGE  4/3/87

       IF(JL.EQ.19 .AND. TFD(IJ).GT.1.E-5)THEN
	  O2JINT(IJ,IK)=O2JINT(IJ,IK)/TFD(IJ)
C  CORRECT BACK TO DAYTIME AVERAGE  4/7/95
       ENDIF

c	if(ij.eq.18 .and. jl.eq.38)print *,' ij=',ij,' ik=',ik,' jl=',
c     *  jl
c	if(ij.eq.18 .and. jl.eq.38)print *,' ldav=',ldav,' tfd=',
c     *  tfd(ij),' tauz=',tauz,' fluxms=',fluxms,' rflux=',
c     *  rflux(jl,ij,ik),' efac=',efac

1000	CONTINUE

300	CONTINUE

2000	CONTINUE

      SAVE
      RETURN
      END