      PROGRAM AUGUST_8

C  include file of dimensions and numbers

      INCLUDE '../includes/VLIDORT.PARS'

C  input variables

      INCLUDE '../includes/VLIDORT_INPUTS.VARS'
c      INCLUDE '../includes/VLIDORT_L_INPUTS.VARS'
      INCLUDE '../includes/VLIDORT_BOOKKEEP.VARS'

C  result variables

      INCLUDE '../includes/VLIDORT_RESULTS.VARS'
c      INCLUDE '../includes/VLIDORT_L_RESULTS.VARS'

C  help variables

      INTEGER          STATUS_INPUTREAD
      INTEGER          STATUS_INPUTCHECK
      INTEGER          STATUS_CALCULATION
      DOUBLE PRECISION PROBLEM_RAY(6,0:2)

C  Local data and holding arrays

      INTEGER          N_LAMBDAS
      DOUBLE PRECISION LAMBDAS(6),DEPOLS(6)
      DOUBLE PRECISION LN_CUMTAU_SCAT(6,101)
      DOUBLE PRECISION LN_CUMTAU_TOTAL(6,101)
      DATA LAMBDAS / 303.239,306.340,310.090, 313.391,316.791,319.991/
      DATA DEPOLS  / 0.0318759, 0.0317431, 0.0315900,
     &               0.0314577, 0.0313350, 0.0312210 /

C  output arrays

      DOUBLE PRECISION IPRLL(6,12)
      DOUBLE PRECISION IPERP(6,12)
      DOUBLE PRECISION ITOTAL(6,12)
      DOUBLE PRECISION PDEGR(6,12)

C  when needed

      INTEGER          ndum, w, n, i, o1, l, ib, um, az, v,dir
      double precision dum,agas,difgas,depol
      double precision deltau, omega, delsca
      double precision tau_s, tau_t, tau_s_prev, tau_t_prev

C  initialize Greek matrix

       DO N = 1, MAXLAYERS
	  DO L = 0, MAXMOMENTS
	    DO O1 = 1, NSTOKES_SQ
	        GREEKMAT_TOTAL_INPUT(L,N,O1)   = ZERO
c	      L_GREEKMAT_TOTAL_INPUT(1,L,N,O1) = ZERO
	    ENDDO
	  ENDDO
       ENDDO
       do_info_write = .true.

C  Test TOMRAD
C  ===========

C  description
C     101 layers, Rayleigh + O3 absorption, 6 UV wavelengths
C     quasi-spherical, classical solution
C     1 solar angles, 6 viewing angles, 2 azimuths (12 geometries)
C     Upwelling radiance at TOA
C     WITH POLARIZATION

      CALL VLIDORT_INPUT_MASTER
     &        ( 'test_tomrad.inp','VLIDORT_LOGFILE',status_inputread)
      if ( status_inputread.ne.vlidort_success)stop'badread test_tomrad'

      DO_FDTEST = .FALSE.
      FDEPSILON = 0.0d0
      IF ( DO_UPWELLING ) DIR = 1
      IF ( DO_DNWELLING ) DIR = 2


C  number of layer, number of wavelengths

      nlayers   = 101
      n_lambdas = 6

C  Read TOMRAD  optical properties

      open(1,file = 'tomrad_input.dat',status='old')
      do i = 1, 12
        read(1,*)
      enddo
      do w = 1, n_lambdas
        do i = 1, 8
         read(1,*)
        enddo
        do n = 1, nlayers
        read(1,*)ndum,dum,dum,dum,
     &     ln_cumtau_scat(w,n),ln_cumtau_total(w,n),dum,dum
        enddo
      enddo
      close(1)

C  set up the height grid

      height_grid(0) = 85.5
      open(1,file='height.vlidortv2test',status='old')
      do n = 1, nlayers
        read(1,*)dum,dum,height_grid(nlayers+1-n),dum
      enddo
      close(1)

C  Start the wavelengths

      DO W = 1, N_LAMBDAS

C  set up the bulk optical properties

        do n = 1, nlayers
         tau_s = dexp(ln_cumtau_scat(w,n))
         tau_t = dexp(ln_cumtau_total(w,n))
         if ( n.eq.1) then
           deltau = tau_t
           omega  = tau_s / tau_t
         else
           deltau = tau_t - tau_t_prev
           delsca = tau_s - tau_s_prev
           omega = delsca / deltau
         endif
         tau_t_prev = tau_t
         tau_s_prev = tau_s
         agas = deltau - delsca
         difgas = agas/deltau
         if ( omega.gt.0.99999d0)omega = 0.99999d0
         deltau_vert_input(n) = deltau
         omega_total_input(n) = omega
c         LAYER_VARY_FLAG(n)   = .TRUE.
c         LAYER_VARY_NUMBER(n) = 1
c         L_DELTAU_VERT_INPUT(1,n) =   difgas
c         L_OMEGA_TOTAL_INPUT(1,n) = - difgas
        enddo

C  set up Rayleigh problem

        DEPOL = DEPOLS(W)
        DO O1 = 1, 6
         PROBLEM_RAY(O1,0) = 0.0D0
         PROBLEM_RAY(O1,1) = 0.0D0
         PROBLEM_RAY(O1,2) = 0.0D0
        ENDDO
        PROBLEM_RAY(2,0) = 1.0D0
        PROBLEM_RAY(4,1) = 
     &     3.0D0 * ( 1.0D0 - 2.0D0*DEPOL ) / (2.0D0 + DEPOL )
        PROBLEM_RAY(2,2) = 
     &        ( 1.0D0 - DEPOL ) / (2.0D0 + DEPOL )
        PROBLEM_RAY(3,2) = 
     &     - DSQRT(6.0D0) * ( 1.0D0 - DEPOL ) / (2.0D0 + DEPOL )
        PROBLEM_RAY(1,2) = 
     &     6.0D0 * ( 1.0D0 - DEPOL ) / (2.0D0 + DEPOL )

C  Fix
C
C        PROBLEM_RAY(5,2) = 0.1D0

C  Set up of Rayleigh scattering law

        do n = 1, nlayers
         DO L = 0, 2
          GREEKMAT_TOTAL_INPUT(L,N,1)  = PROBLEM_RAY(2,L)
          IF ( L.EQ.0)GREEKMAT_TOTAL_INPUT(L,N,1) = 1.0D0
          IF ( NSTOKES .GT. 1 ) THEN
           GREEKMAT_TOTAL_INPUT(L,N,2)  =  PROBLEM_RAY(3,L)
           GREEKMAT_TOTAL_INPUT(L,N,5)  =  PROBLEM_RAY(3,L)
           GREEKMAT_TOTAL_INPUT(L,N,6)  =  PROBLEM_RAY(1,L)
           GREEKMAT_TOTAL_INPUT(L,N,11) =  PROBLEM_RAY(6,L)
           GREEKMAT_TOTAL_INPUT(L,N,12) = - PROBLEM_RAY(5,L)
           GREEKMAT_TOTAL_INPUT(L,N,15) =   PROBLEM_RAY(5,L)
           GREEKMAT_TOTAL_INPUT(L,N,16) =   PROBLEM_RAY(4,L)
          ENDIF
         ENDDO
        enddo

C  Call to VLIDORT_MASTER

        write(*,'(/a,f10.4)')'wavelength:',LAMBDAS(W)
        CALL VLIDORT_MASTER ( STATUS_INPUTCHECK, STATUS_CALCULATION )
        CALL VLIDORT_STATUS ( STATUS_INPUTCHECK, STATUS_CALCULATION )

C  store the output (Stokes vector)

        DO IB = 1, NBEAMS
         DO UM = 1, N_USER_STREAMS
          DO AZ = 1, N_USER_RELAZMS
           V = UMOFF(IB,UM) + AZ
           IPRLL(W,V)  = PIE*HALF*(STOKES(1,V,1,DIR)+STOKES(1,V,2,DIR))           
           IPERP(W,V)  = PIE*HALF*(STOKES(1,V,1,DIR)-STOKES(1,V,2,DIR)) 
           ITOTAL(W,V) = PIE*STOKES(1,V,1,DIR)
c           PDEGR(W,V)  = PIE *DSQRT(STOKES(1,V,2,DIR)*STOKES(1,V,2,DIR)
c     &                    + STOKES(1,V,3,DIR)*STOKES(1,V,3,DIR)
c     &                    + STOKES(1,V,4,DIR)*STOKES(1,V,4,DIR)) /
c     &                ITOTAL(W,V)
           PDEGR(W,V)  = -PIE*STOKES(1,v,2,dir)/
     &                 ITOTAL(W,V)
          IF ( DABS(PDEGR(W,V)).LT.1.0D-12)PDEGR(W,V) = ZERO
          if (V.eq.1)write(*,*)PDEGR(W,V)
          ENDDO
         ENDDO
        ENDDO

C  finish wavelength loop

      ENDDO

C  output

      IF ( NSTOKES.EQ.4 ) THEN
       OPEN(1,FILE='TEST_TOMRAD_RADIANCES.out_4',STATUS='unknown')
       DO W = 1, N_LAMBDAS
        WRITE(1,*)LAMBDAS(W)
        DO IB = 1, NBEAMS
         DO AZ = 1, N_USER_RELAZMS
          DO UM = 1, N_USER_STREAMS
           V = UMOFF(IB,UM) + AZ
           WRITE(1,'(3f9.2,1p4e16.6)')
     &       BEAM_SZAS(IB),USER_ANGLES_INPUT(UM),USER_RELAZMS(AZ),
     &      IPRLL(W,V),IPERP(W,V),ITOTAL(W,V),PDEGR(W,V)
          ENDDO
         ENDDO
        ENDDO
       ENDDO
       CLOSE(1)
      ELSE
       OPEN(1,FILE='TEST_TOMRAD_RADIANCES.out_1',STATUS='unknown')
       DO W = 1, N_LAMBDAS
        WRITE(1,*)LAMBDAS(W)
        DO IB = 1, NBEAMS
         DO AZ = 1, N_USER_RELAZMS
          DO UM = 1, N_USER_STREAMS
           V = UMOFF(IB,UM) + AZ
           WRITE(1,'(3f9.2,1p4e16.6)')
     &       BEAM_SZAS(IB),USER_ANGLES_INPUT(UM),USER_RELAZMS(AZ),
     &      IPRLL(W,V),IPERP(W,V),ITOTAL(W,V),PDEGR(W,V)
          ENDDO
         ENDDO
        ENDDO
       ENDDO
       CLOSE(1)
      ENDIF

c        DO N = 1, NLAYERS
c         DO O1 = 1, NSTOKES
c           DO IB = 1, NBEAMS
c            HELPOUT(O1,IB) = PIE*ATMOSWF(1,N,1,IB,O1,drn)
c            IF ( DABS(HELPOUT(O1,IB)).LT.1.0D-10)  HELPOUT(O1,IB) = ZERO
c           ENDDO
c         ENDDO
c         write(1,'(f8.2,1p20e16.6)')
c     &        (height_grid(n)+height_grid(n-1))*0.5d0,
c     &           (HELPOUT(O1,1),O1 = 1, NSTOKES),
c     &           (HELPOUT(O1,2),O1 = 1, NSTOKES),
c     &           (HELPOUT(O1,3),O1 = 1, NSTOKES)
c        ENDDO

C  Stop

       stop
       end
