d! $Id: JVALUE.f,v 1.1 2009/09/16 14:06:47 bmy Exp $
      SUBROUTINE JVALUE( SA, DIRRAD1, DIFRAD1, DIRRAD2, DIFRAD2,
     &                   DIRRAD3, DIFRAD3, DIRRAD4, DIFRAD4,
     &                   DIRRAD5, DIFRAD5, DIRRAD6, DIFRAD6,
     &                   DIRRAD7, DIFRAD7 )

!-----------------------------------------------------------------------
c  Calculate the actinic flux at each level for the current SZA value.
!        quit when SZA > 98.0 deg ==> tangent height = 63 km
!             or         99.                           80 km
!-----------------------------------------------------------------------
!  Add the following input variables for CTM interface (bmy, 9/13/99)
!
!  Variable  Type    Dimensn Units   Description
!  --------  ----    ------- -----   -----------
!  SA        dble    -       -       Surface Albedo
!-----------------------------------------------------------------------
c
c     AVGF   Attenuation of beam at each level for each wavelength
c     FFF    Actinic flux at each desired level
c     WAVE   Effective wavelength of each wavelength bin
c     XQO2   Absorption cross-section of O2
c     XQO3   Absorption cross-section of O3
c
!-----------------------------------------------------------------------
#     include "define.h"

      USE CMN_FJ_MOD, ONLY : IPAR, JPAR, LPAR, JPMAX, JPPJ, JPNL
      USE JV_CMN_MOD

      IMPLICIT NONE

!=============== INPUT PARAMETERS ======================================
      REAL*8, INTENT(INOUT) :: SA
      REAL*8, INTENT(OUT) :: DIRRAD1, DIFRAD1, DIRRAD2, DIFRAD2
      REAL*8, INTENT(OUT) :: DIRRAD3, DIFRAD3, DIRRAD4, DIFRAD4
      REAL*8, INTENT(OUT) :: DIRRAD5, DIFRAD5, DIRRAD6, DIFRAD6
      REAL*8, INTENT(OUT) :: DIRRAD7, DIFRAD7

!=============== LOCAL VARIABLES =======================================
      integer j, k, IIPAR, JJPAR, I, M, N
      real*8  wave, xseco3, xseco2, ZFLUXX, FJBOTT
      real*8  AVGF(lpar),XQO3(NB),XQO2(NB)

      do J=1,jpnl
        do K=NW1,NW2
          FFF(K,J) = 0.d0
        enddo
      enddo

c---SZA check
c      write(6,1000) SZA, RFLECT, (OD(nslon,nslat,j),j=1,lpar)
      if(SZA.gt.szamax) GOTO 99     
c
!---Calculate spherical weighting functions
      CALL SPHERE
c
!---Loop over all wavelength bins
      do K=NW1,NW2
        WAVE = WL(K)
        do J=1,NB
          XQO3(J) = XSECO3(K,dble(TJ(J)))
        enddo
        do J=1,NB
          XQO2(J) = XSECO2(K,dble(TJ(J)))
        enddo
!-----------------------------------------
        CALL OPMIE(K,WAVE,XQO2,XQO3,AVGF,ZFLUXX,FJBOTT)
!-----------------------------------------
C-----------------------------------------mcz 2/20/2011 (below) (updated 4/17/2012)
        do J=1,jpnl
          !if (J .LT. 2) then
          !FFF(K,J) = a0(1) * (FFF(K,J) + FL(K)*AVGF(J))
          !else
          FFF(K,J) = FFF(K,J) + FL(K)*AVGF(J)
          !write(6,*), 'K,J,FFF', K, J, FFF(K,J)
          !endif
        enddo
      enddo

      !write(6,*), 'FL(1)', FL(1)
      !write(6,*), 'FL(2)', FL(2)
      !write(6,*), 'FL(3)', FL(3)
      !write(6,*), 'FL(4)', FL(4)
      !write(6,*), 'FL(5)', FL(5)
      !write(6,*), 'FL(6)', FL(6)
      !write(6,*), 'FL(7)', FL(7)

      !write(6,*), 'ZFLUXX', ZFLUXX
      !write(6,*), 'SA', SA

      DIRRAD1 = ZFLUXX*FL(1)  !289-298 nm
      !write(6,*), 'DIRRAD1', DIRRAD1
       !write(200,*), 'DIRRAD1', DIRRAD1
      DIRRAD2 = ZFLUXX*FL(2)  !298-307 nm
      !write(6,*), 'DIRRAD2', DIRRAD2
       !write(200,*), 'DIRRAD2', DIRRAD2
      DIRRAD3 = ZFLUXX*FL(3)  !307-312 nm
      !write(6,*), 'DIRRAD3', DIRRAD3
       !write(200,*), 'DIRRAD3', DIRRAD3
      DIRRAD4 = ZFLUXX*FL(4)  !312-320 nm
      !write(6,*), 'DIRRAD4', DIRRAD4
       !write(200,*), 'DIRRAD4', DIRRAD4
      DIRRAD5 = ZFLUXX*FL(5)  !320-345 nm
      !write(6,*), 'DIRRAD5', DIRRAD5
       !write(200,*), 'DIRRAD5', DIRRAD5
      DIRRAD6 = ZFLUXX*FL(6)  !345-412 nm
      !write(6,*), 'DIRRAD6', DIRRAD6
       !write(200,*), 'DIRRAD6', DIRRAD6
      DIRRAD7 = ZFLUXX*FL(7)  !412-850 nm
      !write(6,*), 'DIRRAD7', DIRRAD7
       !write(200,*), 'DIRRAD7', DIRRAD7

      DIFRAD1 = FJBOTT*FL(1)  !289-298 nm
      !write(6,*), 'DIFRAD1', DIFRAD1
       !write(200,*), 'DIFRAD1', DIFRAD1
      DIFRAD2 = FJBOTT*FL(2)  !298-307 nm
      !write(6,*), 'DIFRAD2', DIFRAD2
       !write(200,*), 'DIFRAD2', DIFRAD2
      DIFRAD3 = FJBOTT*FL(3)  !307-312 nm
      !write(6,*), 'DIFRAD3', DIFRAD3
       !write(200,*), 'DIFRAD3', DIFRAD3
      DIFRAD4 = FJBOTT*FL(4)  !312-320 nm
      !write(6,*), 'DIFRAD4', DIFRAD4
       !write(200,*), 'DIFRAD4', DIFRAD4
      DIFRAD5 = FJBOTT*FL(5)  !320-345 nm
      !write(6,*), 'DIFRAD5', DIFRAD5
       !write(200,*), 'DIFRAD5', DIFRAD5
      DIFRAD6 = FJBOTT*FL(6)  !345-412 nm
      !write(6,*), 'DIFRAD6', DIFRAD6
       !write(200,*), 'DIFRAD6', DIFRAD6
      DIFRAD7 = FJBOTT*FL(7)  !412-850 nm
      !write(6,*), 'DIFRAD7', DIFRAD7
       !write(200,*), 'DIFRAD7', DIFRAD7

      !write(6,*), 'ZFLUXX*FL(3)', ZFLUXX*FL(3)
      !write(6,*), 'DIRRAD3_JVALUEFILE',DIRRAD3 

c---------------------------------------mcz 2/20/2011(above) (updated 4/17/2012)

   99 continue
 1000 format('  SZA=',f6.1,' Reflectvty=',f6.3,' OD=',10(1pe10.3))
      return
      end