!********************************************************
!  These subroutines are in a module so I can pass
!  allocated arrays back and forth
!********************************************************

MODULE RADIANT2S_DRIVER

CONTAINS

!********************************************************
!  RADIANT_DRIVER SUBROUTINE
!********************************************************

SUBROUTINE RADIANT_DRIVER(SZA,TSURF,VIEW_ANGLE,Z_LVL,P_LVL,T_LVL,&
  TAU_CLD,OMEGA_CLD,ASPAR_CLD,NUMLAY,O3_LVL,CO2_LVL,H2O_LVL,&
  SIGS_RAY_IN,KABS_H2O,KABS_CO2,ALBEDO,ITPTOT)

!THIS PROGRAM SERVES AS A DRIVER PROGRAM FOR RADIANT 2.0

!PROGRAMMER: MATT CHRISTI

!DATE LAST MODIFIED: 8/09/04

!INPUTS*************************************************************************
!  HEIGHT, TEMPERATURE, AND MASS MIXING RATIOS OF GASES IN CENTER
!  OF LAYERS (in km, K, and unitless)
!  PRESSURE AT CENTER OF LAYERS (in Pa)
!*******************************************************************************


!DATA DICTIONARY****************************************************************
!
! ALBEDO       = SURFACE ALBEDO WHEN SURF=1
! B_IO         = SUBROUTINE BUILD_LAYER I/O FLAG
!		   (0=OFF, 1=ON)
! CHI_AERO     = MATRIX CONTAINING COEFFICIENTS FOR LEGENDRE POLYNOMIAL
!                EXPANSION OF AEROSOL PHASE FUNCTION FOR ANY LAYER (IF USED)
! CHI_CLD      = MATRIX CONTAINING COEFFICIENTS FOR LEGENDRE POLYNOMIAL
!                EXPANSION OF CLOUD PHASE FUNCTION FOR ANY LAYER (IF USED)
! CHI_SURF     = VECTOR CONTAINING COEFFICIENTS FOR LEGENDRE POLYNOMIAL
!                EXPANSION OF SURFACE BIDIRECTIONAL REFLECTIVITY (IF USED)
! DELTA_M_AERO = VECTOR OF FLAGS TO INDICATE WHETHER OR NOT DELTA_M SCALING IS
!                USED IN SUBROUTINE LOCAL FOR POSSIBLE AEROSOL PHASE FUNCTION
!		   IN EACH LAYER
!		   (0=NO, 1=YES)
! DELTA_M_CLD  = VECTOR OF FLAGS TO INDICATE WHETHER OR NOT DELTA_M SCALING IS
!                USED IN SUBROUTINE LOCAL FOR POSSIBLE CLOUD PHASE FUNCTION
!		   IN EACH LAYER
!		   (0=NO, 1=YES)
! DELTAZ       = VECTOR CONTAINING THICKNESSES OF LAYERS UNDER COMPUTATION
!	          (IN km)
! DG	        = FLAG TO INDICATE WHICH VERSION OF GAUSS QUADRATURE WILL BE
!	          USED IN SUBROUTINE LOCAL IF QUAD = 0 IS SELECTED
!	          (0=NORMAL GAUSS, 1=DOUBLE GAUSS)
! DISPLAY_IN   = I/O FLAG TO PERMIT DISPLAY OF USER INPUT TO SCREEN
!		   (0=OFF, 1=ON)
! DISPLAY_OUT  = I/O FLAG TO PERMIT DISPLAY OF RADIANCE OUTPUT TO SCREEN
!		   (0=OFF, 1=ON)
! EPS	        = AMOUNT BY WHICH OPTICAL PROPERTIES ARE PERTURBED BY FOR THE
!		   COMPUTATION OF JACOBIAN ELEMENTS BY FINITE DIFFERENCE WHEN
!		   THE PERTURBATION SCHEME IS IN USE
! FIB	        = FLUX INTO THE BOTTOM BOUNDARY
! FIT          = FLUX INTO THE TOP BOUNDARY
! FLUX_IO      = FLUX TEST I/O FLAG FOR SUBROUTINE RAD
!		   (0=OFF, 1=ON)
! FLUXES       = RADIANT FLUX OUTPUT I/O FLAG
!		   (0=OFF, 1=ON)
! FOB          = FLUX OUT OF THE BOTTOM BOUNDARY
! FOT          = FLUX OUT OF THE TOP BOUNDARY
! FSUN         = FLUX AT THE TOP OF THE ATMOSPHERE AT A GIVEN WAVELENGTH
! G_AERO       = VECTOR OF AEROSOL EFFECTIVE ASYMMETRY PARAMETERS
! G1_AERO      = VECTOR OF AEROSOL ASYMMETRY PARAMETERS FOR FORWARD SCATTERING
! G2_AERO      = VECTOR OF AEROSOL ASYMMETRY PARAMETERS FOR BACKWARD SCATTERING
! G_CLD        = VECTOR OF CLOUD EFFECTIVE ASYMMETRY PARAMETERS
! G1_CLD       = VECTOR OF CLOUD ASYMMETRY PARAMETERS FOR FORWARD SCATTERING
! G2_CLD       = VECTOR OF CLOUD ASYMMETRY PARAMETERS FOR BACKWARD SCATTERING
! INV_IO       = SUBROUTINE INVERT I/O FLAG (FOR MATRIX INVERSE)
!		   (0=OFF, 1=ON)
! IBMTOT       = VECTOR OF TOTAL DOWNWELLING RADIANCES AT THE BOTTOM OF THE
!		   ATMOSPHERE FROM ALL FOURIER COMPONENTS COMPUTED
! IBPTOT       = VECTOR OF TOTAL UPWELLING RADIANCE AT THE BOTTOM OF THE
!		   ATMOSPHERE FROM ALL FOURIER COMPONENTS COMPUTED
! ITM	        = VECTOR OF TOTAL DOWNWELLING DIFFUSE RADIANCES FROM THE TOP
!		   BOUNDARY
! ITMS	        = VECTOR OF DOWNWELLING DIFFUSE SOLAR RADIANCES FROM THE TOP
!		   BOUNDARY
! ITMT	        = VECTOR OF DOWNWELLING DIFFUSE THERMAL RADIANCES FROM THE TOP
!		   BOUNDARY
! ITPTOT       = VECTOR OF TOTAL UPWELLING RADIANCES AT THE TOP BOUNDARY
!		   FROM ALL FOURIER COMPONENTS COMPUTED
! KFLAG	 = FLAG TO INDICATE WHETHER OR NOT ADDITIONAL COMPUTATIONS WILL
!		   BE PERFORMED INSIDE RADIANT TO PREPARE FOR SUBSEQUENT CALLS
!		   TO RADIANT WHERE SAVED MATRICES WILL BE USED
!		   (0=NO, 1=YES)
!		   NOTE: KFLAG USED IN CONJUNCTION WITH WVN_FLAG
! L_IO         = SUBROUTINE LOCAL I/O FLAG
!		   (0=OFF, 1=ON)
! LAYCHG       = LAYER NUMBER OF CHANGED LAYER WHEN IN LAYER-SAVING MODE
! NUMLAY       = NUMBER OF LAYERS IN THE MODEL ATMOSPHERE
! OMEGA_AERO   = VECTOR OF SINGLE SCATTER ALBEDO FOR AEROSOL
! OMEGA_CLD    = VECTOR OF SINGLE SCATTER ALBEDO FOR CLOUD
! P_IO         = SUBROUTINE RADIANT I/O FLAG
!		   (0=OFF, 1=ON)
! PERT_FLAG    = FLAG TO INDICATE WHETHER OR NOT THE PERTURBATION SCHEME WILL
!		   BE USED TO BUILD LAYERS ON SUBSEQUENT CALLS TO RADIANT DURING
!		   LAYER SAVING MODE
!		   (0=NO, 1=YES)
! PF_AERO      = VECTOR OF FLAGS TO DETERMINE HOW COEFFICIENTS OF AEROSOL PHASE
!	          FUNCTION EXPANSION ARE OBTAINED
!		   (1=HENYEY-GREENSTEIN - SINGLE OR DOUBLE DEPENDING ON INPUT,
!		    2=DETERMINED BY CHI_AERO)
! PF_CLD       = VECTOR OF FLAGS TO DETERMINE HOW COEFFICIENTS OF CLOUD PHASE
!	          FUNCTION EXPANSION ARE OBTAINED
!		   (1=HENYEY-GREENSTEIN - SINGLE OR DOUBLE DEPENDING ON INPUT,
!		    2=DETERMINED BY CHI_CLD)
! PHI	        = AZIMUTH ANGLE FOR WHICH THE RADIANCE VECTOR IS COMPUTED
!		   (in degrees)
!		   NOTE: REFERENCE AZIMUTH ASSUMED TO BE 0
! PLANCK_FLAG  = FLAG TO CHOOSE WHETHER SPECTRAL RADIANCES OR RADIANCES WILL
!		   BE USED FOR THERMAL SOURCES
!		   (0=SPECTRAL RADIANCES (in mW/(m^2 cm^-1 ster))
!		    1=RADIANCES          (in mW/(m^2 ster))     )
! QUAD	        = FLAG TO INDICATE WHICH QUADRATURE SCHEME WILL BE USED
!		   (0=GAUSS, 1=LOBATTO)
! R_IO         = SUBROUTINE RAD I/O FLAG
!		   (0=OFF, 1=ON)
! SRC_FLAG     = FLAG INDICATING TYPE OF SOURCES INCLUDED IN COMPUTING
!	          RADIANCES
!		   (1=SOLAR ONLY, 2=SOLAR & THERMAL, 3=THERMAL ONLY)
! STREAMS      = NUMBER OF STREAMS TO COMPUTE RADIANCE FIELD
! SURF         = FLAG INDICATING TYPE OF REFLECTING SURFACE
!		   (0=NONE, 1=LAMBERTIAN, 2=DETERMINED BY CHI_SURF)
! SZA	        = THE SOLAR ZENITH ANGLE (in degrees from zenith)
! TAU_AERO     = VECTOR OF OPTICAL DEPTHS FOR AEROSOL
! TAU_CLD      = VECTOR OF OPTICAL DEPTHS FOR CLOUD
! TAU_GAS      = VECTOR OF OPTICAL DEPTHS FOR ATMOSPHERIC GAS
! TAU_RAY      = VECTOR OF OPTICAL DEPTHS FOR RAYLEIGH SCATTERING
! TEMISS       = EMISSIVITY OF TOP BOUNDARY FOR THERMAL SOURCES
! TLEV         = VECTOR (I.E. VERTICAL PROFILE) OF LEVEL-BASED ATMOSPHERIC
!		   TEMPERATURES (in K)
! TSURF	 = TEMPERATURE OF THE SURFACE (in K)
! TTOP	        = TEMPERATURE OF THE TOP BOUNDARY (in K)
! VIEW_ANGLE   = THE USER-DEFINED VIEW ANGLE (in degrees from zenith)
! VIEW_FLAG    = FLAG INDICATING USER-DEFINED VIEW ANGLE IN USE
!		   (0=OFF, 1=ON)
! WVN          = WAVENUMBER FOR WHICH SPECTRAL RADIANCE IS COMPUTED FOR PLANCK
!	          FUNCTION (in cm^-1)
! WVNLO        = LOW WAVENUMBER IN SPECTRAL INTERVAL FOR WHICH RADIANCE IS
!	          COMPUTED FROM INTEGRAL OF PLANCK FUNCTION (in cm^-1)
! WVNHI        = HIGH WAVENUMBER IN SPECTRAL INTERVAL FOR WHICH RADIANCE IS
!	          COMPUTED FROM INTEGRAL OF PLANCK FUNCTION (in cm^-1)
! WVN_FLAG     = FLAG TO INDICATE WHETHER OR NOT SAVED MATRICES WILL BE
!	          USED TO COMPUTE RADIANCES ON CURRENT CALL TO SUBROUTINE
!	          RADIANT
!		   ('NEW'=NO, 'OLD'=YES)
!		   NOTE: WVN_FLAG USED IN CONJUNCTION WITH KFLAG
! WARNING_FLAG = RADIANT WARNING MESSAGES I/0 FLAG
!		   (0=OFF, 1=ON)
!
!*******************************************************************************

	USE RADIANT2S

	IMPLICIT NONE
!INPUT VARIABLES (DRIVER)
	INTEGER :: &
	  DISPLAY_IN,DISPLAY_OUT,STREAMS
!INPUT VARIABLES (RADIANT)
	INTEGER :: &
	  B_IO,DG,FLUX_IO,FLUXES,INV_IO,KFLAG,L_IO,N,&
	  NUMLAY,P_IO,PERT_FLAG,PLANCK_FLAG,QUAD,R_IO,SRC_FLAG,SURF,&
	  VIEW_FLAG,WARNING_FLAG,LAYCHG
	INTEGER, DIMENSION(:), ALLOCATABLE :: &
	  DELTA_M_AERO,DELTA_M_CLD,PF_AERO,PF_CLD
	DOUBLE PRECISION :: &
	  ALBEDO,EPS,FSUN,FOURIER_TOL,PHI,TEMISS,SZA,TSURF,TTOP,&
	  VIEW_ANGLE,WVN,WVNLO,WVNHI
	DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: &
	  G_AERO,G_CLD,G1_AERO,G1_CLD,G2_AERO,G2_CLD,&
	  ITMS,ITMT,OMEGA_AERO,OMEGA_CLD,TAU_AERO,TAU_CLD,TAU_GAS,&
	  TAU_RAY,TLEV,CHI_SURF
	DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: &         
         CHI_AERO,CHI_CLD
	CHARACTER(LEN=3) :: &
	  WVN_FLAG
	LOGICAL :: &
	  AZIMUTHAL_RAD_ONLY  
!OUTPUT VARIABLES (RADIANT)
	DOUBLE PRECISION :: &
	  FIT,FOT,FOB,FIB
	DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: &
	  ITM,IBMTOT,IBPTOT,ITPTOT,MU
!INTERNAL VARIABLES (DRIVER)
	INTEGER :: &
	  L,DUMMY,I,J,NUMEXP,START,WVN_INDEX,X_SIZE
        DOUBLE PRECISION :: &
	  DUMVAR,GRAVITY,M_CO2,M_DRY_AIR,M_H2O,&
	  R_DRY_AIR,TAU_GAS_SAVE,TSLOPE
        DOUBLE PRECISION :: &
	  PI=3.1415926535897932D0
	DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: &
	  DELTAZ,F,RHO_AIR_LVL,TAU_CO2,TAU_H2O
	DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: &
	  K
	CHARACTER(LEN=5) :: &
	  STRING2
	CHARACTER(LEN=17) :: &
	  FORM
	CHARACTER(LEN=30) :: &         
         PHASE_INFILE
!MODEL ATMOSPHERE INPUT PARAMETERS
	DOUBLE PRECISION :: &
	  PSURF
	INTEGER, DIMENSION(:), ALLOCATABLE :: &
	  LVL
	DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: &
	  ASPAR_AERO,ASPAR_CLD,CO2_LVL,H2O_LVL,&
	  O3_LVL,P_LVL,T_LVL,Z_LVL
	DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: &
	  KABS_CO2,KABS_H2O,SIGS_RAY_IN
         
!START PROGRAM
	print*,'entering radiant 2s driver'

!DEFINE DRIVER INPUT
!Flag for displaying user input to the screen, 0=no,1=yes
	DISPLAY_IN = 0

!DEFINE SOME RADIANT INPUTS (constants)
    STREAMS = 16   ! was 8 SJC
    QUAD = 1
    DG = 0
    VIEW_FLAG = 1  !user defined view angle
    PHI = 0.0D0    !azimuth angle
    EPS = 1.0D-5   !FOR RADIANTN CHECK
    AZIMUTHAL_RAD_ONLY = .FALSE.!.TRUE.
    FOURIER_TOL = 1.0D-8
    SRC_FLAG = 1   !solar sources only
    FSUN = 1.0D0   !flux at toa at this wavelength
    PLANCK_FLAG = 2  !radiant flag
    WVN = 0.0D0    !wavelength for planck function
    WVNLO = 0.0D0  !radiant variable
    WVNHI = 1.0D0  !radiant variable
    TTOP = 0.0D0   !temperature of the top boundary
    TEMISS = 0.0D0 !emissivity of top boundary
    PERT_FLAG = 0  !use perterbation scheme 0=no,1=yes
    FLUXES = 0     !flux output flag 0=off,1=on
    WARNING_FLAG = 0  !radiant warning messages 0=off,1=on
    P_IO = 0       !radiant I/O flag 0=off,1=on
    R_IO = 0       !rad I/O flag 0=off,1=on
    FLUX_IO = 0    !flux test I/O flag 0=off,1=on
    L_IO = 0       !local I/O flag 0=off,1=on
    B_IO = 0       !build_layer I/O flag 0=off,1=on
    INV_IO = 0     !invert I/O flag 0=off,1=on
    SURF = 1.0     !type of reflecting surface 1=lambertian
    WVN_FLAG = 'NEW'  !use saved matrices
    KFLAG = 1      !use saved matrices
    LAYCHG = 0     !layer saving mode

!ALLOCATE MEMORY FOR RADIANT ARRAYS
    N = STREAMS/2 + VIEW_FLAG

!*** numlay is the number of layer midpoints
!*** numlay+1 gives the number of boundary layers

!allocate aerosol arrays (numlay)  
    ALLOCATE( PF_AERO(NUMLAY),DELTA_M_AERO(NUMLAY),&
    CHI_AERO(0:(4*N),NUMLAY),&
    G1_AERO(NUMLAY),G2_AERO(NUMLAY),G_AERO(NUMLAY),&
    OMEGA_AERO(NUMLAY),TAU_AERO(NUMLAY),ASPAR_AERO(NUMLAY))

!allocate arrays (numlay)  
    ALLOCATE ( PF_CLD(NUMLAY),DELTA_M_CLD(NUMLAY),&
    CHI_CLD(0:(4*N),NUMLAY),&
    G1_CLD(NUMLAY),G2_CLD(NUMLAY),G_CLD(NUMLAY),&
    TAU_GAS(NUMLAY),TAU_RAY(NUMLAY),&
    LVL(NUMLAY),RHO_AIR_LVL(NUMLAY),&
    TAU_CO2(NUMLAY),TAU_H2O(NUMLAY),DELTAZ(NUMLAY))

!boundary values, not midpoints
    ALLOCATE( TLEV(NUMLAY+1) )

!allocate arrays (N)
    ALLOCATE(CHI_SURF(0:(4*N)),ITMS(N),ITMT(N),ITM(N),&
    IBMTOT(N),IBPTOT(N),MU(N))!,ITPTOT(N), )

!  Only allocate if it hasn't already been allocated
    if (.not. allocated(ITPTOT)) allocate(ITPTOT(N))

!This sets the aerosol parameter profiles to 0, we aren't considering aerosols
!DEFINE INPUT AEROSOL PARAMETERS (OPTICAL DEPTH, SINGLE SCATTER
!ALBEDO, ASSYMETRY PARAMETER).  
    TAU_AERO   = 0.0D0
    OMEGA_AERO = 0.0D0
    ASPAR_AERO = 0.0D0

!VECTOR OF DOWNWELLING DIFFUSE SOLAR RADIANCES FROM THE TOP BOUNDARY
    ITMS = 0.0D0  !radiant variable
!VECTOR OF DOWNWELLING DIFFUSE THERMAL RADIANCES FROM THE TOP BOUNDARY    
    ITMT = 0.0D0  !radiant variable
!VECTOR CONTAINING COEFFICIENTS FOR LEGENDRE POLYNOMIAL
!EXPANSION OF SURFACE BIDIRECTIONAL REFLECTIVITY (IF USED)    
    CHI_SURF = 0.0D0

    !THICKNESSES OF LAYERS (in km)
    !also checks heights to confirm they are toa to surface
    DO I=1,NUMLAY-1
      DELTAZ(I) = Z_LVL(I) - Z_LVL(I+1)
      IF (DELTAZ(I) < 0.0D0) THEN
        WRITE(*,*) 'DELTAZ(I) < 0 FOR I = ',I
        STOP
      END IF
    END DO
    DELTAZ(NUMLAY)=DELTAZ(NUMLAY-1)   !assign last element

    !AIR DENSITY AT CENTER OF LAYERS (in kg/m^3)
    R_DRY_AIR = 287.05D0  !in J/(kg*K)
    DO I=1,NUMLAY
      RHO_AIR_LVL(I) = P_LVL(I)/(R_DRY_AIR*T_LVL(I))
    END DO

    !*** should this be given values instead of set to zero?
    TLEV = 0.0D0  !vertical profile of level based atmospheric temps (K)
    !print*,'tlev',(tlev(i),i=1,NUMLAY+1)

    !  Inititalize gas variables
    DO I=1,NUMLAY
      !RAYLEIGH OPTICAL DEPTH
      TAU_RAY(I) = SIGS_RAY_IN(I)*DELTAZ(I)
        
      !GAS OPTICAL DEPTHS
      TAU_CO2(I) = KABS_CO2(I)*CO2_LVL(I)* &
                   RHO_AIR_LVL(I)*1.0D2*DELTAZ(I)
    
      TAU_H2O(I) = KABS_H2O(I)*H2O_LVL(I)* &
                   RHO_AIR_LVL(I)*1.0D2*DELTAZ(I)

      TAU_GAS(I) = TAU_CO2(I) + TAU_H2O(I)		
    END DO  !end of loop through layers

    !  Inititalize cloud variables
    DO I=1,NUMLAY
      PF_CLD(I)      = 1
      G1_CLD(I)      = ASPAR_CLD(I)
      G2_CLD(I)      = 0.0D0
      G_CLD(I)       = ASPAR_CLD(I)                        
	    
      CHI_CLD(:,I)   = 0.0D0
      DELTA_M_CLD(I) = 0.0
    END DO  !end of loop through numlay

    !  Initializing aerosol variables
    DO I=1,NUMLAY
      PF_AERO(I)      = 1
      G1_AERO(I)      = ASPAR_AERO(I)
      G2_AERO(I)      = 0.0D0
      G_AERO(I)       = ASPAR_AERO(I)
	    
      CHI_AERO(:,I)   = 0.0D0
      DELTA_M_AERO(I) = 1
    END DO  !end of loop through numlay
         

!************************************************************
!  This is the call to radiant that we want
!************************************************************

    !DISPLAY INPUTS IF DESIRED
    IF (DISPLAY_IN == 1) THEN
      CALL DISPLAY_INPUT(STREAMS,N,NUMLAY,VIEW_FLAG,VIEW_ANGLE,&
      SRC_FLAG,FSUN,SZA,ITMS,ITMT,TTOP,TEMISS,TLEV,TSURF,&
      PLANCK_FLAG,WVN,WVNLO,WVNHI,PHI,QUAD,DG,&
      TAU_RAY,TAU_GAS,TAU_AERO,TAU_CLD,OMEGA_AERO,OMEGA_CLD,&
      G1_CLD,G2_CLD,G_CLD,G1_AERO,G2_AERO,G_AERO,PF_CLD,&
      CHI_CLD,PF_AERO,CHI_AERO,DELTA_M_CLD,DELTA_M_AERO,&
      SURF,ALBEDO,CHI_SURF,WVN_FLAG,KFLAG,AZIMUTHAL_RAD_ONLY,&
      FOURIER_TOL,FLUXES,WARNING_FLAG,&
      P_IO,R_IO,FLUX_IO,L_IO,B_IO,INV_IO)
    END IF

    print*,'before radiant'
    !COMPUTE RADIANCE OF BASE STATE AT THE GIVEN WAVENUMBER
    CALL RADIANT(N,NUMLAY,VIEW_FLAG,VIEW_ANGLE,SRC_FLAG,FSUN,SZA,&
      ITMS,ITMT,TTOP,TEMISS,TLEV,TSURF,PLANCK_FLAG,WVN,WVNLO,WVNHI,&
      PHI,QUAD,DG,TAU_RAY,TAU_GAS,TAU_AERO,TAU_CLD,OMEGA_AERO,&
      OMEGA_CLD,G1_CLD,G2_CLD,G_CLD,G1_AERO,G2_AERO,G_AERO,&
      PF_CLD,CHI_CLD,PF_AERO,CHI_AERO,DELTA_M_CLD,DELTA_M_AERO,&
      SURF,ALBEDO,CHI_SURF,WVN_FLAG,KFLAG,LAYCHG,PERT_FLAG,EPS,&
      AZIMUTHAL_RAD_ONLY,FOURIER_TOL,FLUXES,WARNING_FLAG,&
      P_IO,R_IO,FLUX_IO,L_IO,B_IO,INV_IO,ITM,IBMTOT,IBPTOT,ITPTOT,&
      FIT,FOT,FOB,FIB,MU)
    print*,'itptot',itptot(9),mu(9)  ! final answer we need
              

    !print*,'itptot',itptot(9),mu(9)
    !write(*,*) 'top radiance',ITPTOT
    !write(*,*)  ITPTOT
      !PRINT*,MU,'MU'
    !OPEN OUTPUT FILE
    !OPEN(UNIT=17, FILE='forward_output.dat')
    !WRITE(17,*) 'TOP RADIANCE'
    !WRITE(17,*) ITPTOT(1), MU(1)
    !WRITE(17,*) ITPTOT(2), MU(2)
    !WRITE(17,*) ITPTOT(3), MU(3)
    !WRITE(17,*) ITPTOT(4), MU(4)
    !WRITE(17,*) ITPTOT(5), MU(5)
    !WRITE(17,*) ITPTOT(6), MU(6)
    !WRITE(17,*) ITPTOT(7), MU(7)
    !WRITE(17,*) ITPTOT(8), MU(8)
    !WRITE(17,*) ITPTOT(9), MU(9)
    !close(17)

!DEALLOCATE ARRAYS TO PREPARE FOR NEXT TEST   

	DEALLOCATE ( PF_AERO,PF_CLD,DELTA_M_AERO,DELTA_M_CLD,&
         CHI_AERO,CHI_CLD,OMEGA_AERO,&
         G1_AERO,G2_AERO,G_AERO,G1_CLD,G2_CLD,G_CLD,&
         TAU_AERO,TAU_GAS,TAU_RAY,&
         CHI_SURF,TLEV,ITMS,ITMT,ITM,IBMTOT,IBPTOT)
        
	DEALLOCATE ( LVL,RHO_AIR_LVL,DELTAZ,ASPAR_AERO,TAU_CO2,TAU_H2O)

    !  Don't deallocate these, they are allocated in the main
    !     Z_LVL,T_LVL,P_LVL,CO2_LVL,H2O_LVL,O3_LVL,ASPAR_CLD,KABS_CO2,KABS_H2O,SIGS_RAY_IN &
    !     TAU_CLD,OMEGA_CLD,ITPTOT

	!PRINT*,'radiant2s_driver_example2.f90 DONE'

	END subroutine RADIANT_DRIVER

!*******************************************************************************
!*******************************************************************************
	SUBROUTINE DISPLAY_INPUT(STREAMS,N,NUMLAY,VIEW_FLAG,VIEW_ANGLE,&
	       SRC_FLAG,FSUN,SZA,ITMS,ITMT,TTOP,TEMISS,TLEV,TSURF,&
		PLANCK_FLAG,WVN,WVNLO,WVNHI,PHI,QUAD,DG,&
		TAU_RAY,TAU_GAS,TAU_AERO,TAU_CLD,OMEGA_AERO,OMEGA_CLD,&
		G1_CLD,G2_CLD,G_CLD,G1_AERO,G2_AERO,G_AERO,PF_CLD,&
		CHI_CLD,PF_AERO,CHI_AERO,DELTA_M_CLD,DELTA_M_AERO,&
		SURF,ALBEDO,CHI_SURF,WVN_FLAG,KFLAG,AZIMUTHAL_RAD_ONLY,&
	       FOURIER_TOL,FLUXES,WARNING_FLAG,&
		P_IO,R_IO,FLUX_IO,L_IO,B_IO,INV_IO)

!INPUT : STREAMS,N,NUMLAY,VIEW_FLAG,VIEW_ANGLE,
!        SRC_FLAG,FSUN,SZA,ITMS,ITMT,TTOP,TEMISS,TLEV,TSURF,
!	  PLANCK_FLAG,WVN,WVNLO,WVNHI,PHI,QUAD,DG,
!	  TAU_RAY,TAU_GAS,TAU_AERO,TAU_CLD,OMEGA_AERO,OMEGA_CLD,
!	  G1_CLD,G2_CLD,G_CLD,G1_AERO,G2_AERO,G_AERO,PF_CLD,
!	  CHI_CLD,PF_AERO,CHI_AERO,DELTA_M_CLD,DELTA_M_AERO,
!	  SURF,ALBEDO,CHI_SURF,WVN_FLAG,KFLAG,FLUXES,WARNING_FLAG,
!	  P_IO,R_IO,FLUX_IO,L_IO,B_IO,INV_IO
!OUTPUT: NONE

!THIS PROGRAM DISPLAYS THE INPUTS FOR RADIANT

!PROGRAMMER: MATT CHRISTI

	IMPLICIT NONE
!INPUT VARIABLES (DRIVER)
	INTEGER :: &
	  STREAMS
!INPUT VARIABLES (RADIANT)
	INTEGER :: &
	  B_IO,DG,FLUX_IO,FLUXES,INV_IO,KFLAG,L_IO,N,&
	  NUMLAY,P_IO,PLANCK_FLAG,QUAD,R_IO,SRC_FLAG,&
	  SURF,TEST,VIEW_FLAG,WARNING_FLAG
	INTEGER, DIMENSION(NUMLAY) :: &
	  DELTA_M_AERO,DELTA_M_CLD,PF_AERO,PF_CLD
	DOUBLE PRECISION :: &
	  ALBEDO,BSURF,BTOP,FSUN,PHI,FOURIER_TOL,&
	  SZA,TEMISS,TSURF,TTOP,VIEW_ANGLE,WVN,WVNLO,WVNHI
	DOUBLE PRECISION, DIMENSION(N) :: &
	  ITMS,ITMT
	DOUBLE PRECISION, DIMENSION(NUMLAY) :: &
	  G_AERO,G_CLD,G1_AERO,G1_CLD,G2_AERO,G2_CLD,&
	  TAU_RAY,TAU_GAS,TAU_AERO,TAU_CLD,OMEGA_AERO,OMEGA_CLD
	DOUBLE PRECISION, DIMENSION(NUMLAY+1) :: &
	  TLEV
  	DOUBLE PRECISION, DIMENSION(0:(4*N)) :: &
         CHI_SURF       
	DOUBLE PRECISION, DIMENSION(0:(4*N),NUMLAY) :: &         
         CHI_AERO,CHI_CLD         
	CHARACTER (LEN=3) :: &
         WVN_FLAG
	LOGICAL :: &
	  AZIMUTHAL_RAD_ONLY
!INTERNAL VARIABLES
 	INTEGER :: &
	  I,MOMENT
	DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: &
	  TAU
	CHARACTER (LEN=3) :: &
	  FORMAT_STRING
	CHARACTER (LEN=4) :: &
	  FORM1
	CHARACTER (LEN=20) :: &
	  FORM2                     

!START PROGRAM
!	PRINT*
!	PRINT*,'*** THE INPUTS FOR THE ATMOSPHERIC SCENE ARE: ***'

!	PRINT*
!	PRINT*,'# OF STREAMS TO COMPUTE RADIANCE FIELD = ',STREAMS
	IF (QUAD == 0) THEN
	  IF (DG == 0) THEN
	    PRINT*,'USING GAUSS'
	  ELSE
	    PRINT*,'USING DOUBLE GAUSS'
	  END IF
	ELSE IF (QUAD == 1) THEN
	  PRINT*,'USING LOBATTO'
	END IF
!	PRINT*,'NUMLAY = ',NUMLAY

!	PRINT*
!	PRINT*,'VIEW_FLAG = ',VIEW_FLAG
	IF (VIEW_FLAG == 1) THEN
	  PRINT*,'VIEW_ANGLE = ',VIEW_ANGLE
	END IF
	PRINT*
	PRINT*,'PHI = ',PHI
	PRINT*,'AZIMUTHAL_RAD_ONLY = ',AZIMUTHAL_RAD_ONLY
	PRINT*,'FOURIER_TOL = ',FOURIER_TOL

	PRINT*
	PRINT*,'SRC_FLAG = ',SRC_FLAG
	IF (SRC_FLAG == 1) THEN
	  PRINT*,'ONLY SOLAR SOURCES IN USE'
	ELSE IF (SRC_FLAG == 2) THEN
	  PRINT*,'SOLAR AND THERMAL SOURCES IN USE'
	ELSE IF (SRC_FLAG == 3) THEN
	  PRINT*,'ONLY THERMAL SOURCES IN USE'
	END IF

	IF ((SRC_FLAG == 1).OR.(SRC_FLAG == 2)) THEN
	  PRINT*
	  PRINT*,'FSUN = ',FSUN
	  PRINT*,'SZA = ',SZA

	  PRINT*
	  PRINT*,'ITMS:'
	  DO I=N,1,-1
	    WRITE(*,10) N-I+1,ITMS(I)
	  END DO
	END IF

	IF ((SRC_FLAG == 2).OR.(SRC_FLAG == 3)) THEN
	  IF (TTOP == 0.0D0) THEN
	    PRINT*
	    PRINT*,'ITMT:'
	    DO I=N,1,-1
	      WRITE(*,10) N-I+1,ITMT(I)
	    END DO
	  ELSE
	    PRINT*
	    PRINT*,'TTOP = ',TTOP
	    PRINT*
	    PRINT*,'TEMISS = ',TEMISS
	  END IF

	  PRINT*
	  WRITE(*,40) 'LEVEL','  TLEV  '
	  DO I=1,NUMLAY+1
	    WRITE(*,10) I,TLEV(I)
	  END DO
	  PRINT*
	  PRINT*,'TSURF = ',TSURF

	  PRINT*
	  PRINT*,'PLANCK_FLAG = ',PLANCK_FLAG
	  IF (PLANCK_FLAG == 1) THEN
	    PRINT*,'WVN = ',WVN
	  ELSE IF (PLANCK_FLAG == 2) THEN
	    PRINT*,'WVNLO = ',WVNLO
	    PRINT*,'WVNHI = ',WVNHI
	  END IF
	END IF

	PRINT*
	WRITE(*,40) 'LAYER','TAU_RAY '
	DO I=1,NUMLAY
	  WRITE(*,10) I,TAU_RAY(I)
	END DO

	PRINT*
	WRITE(*,40) 'LAYER','TAU_GAS '
	DO I=1,NUMLAY
	  WRITE(*,10) I,TAU_GAS(I)
	END DO

       !PREPARE FORMATTING FOR CLOUD & AEROSOL DATA
       FORM1 = '(I1)'
       WRITE(FORMAT_STRING,FORM1) 5  
       FORM2 = '(' // FORMAT_STRING // '(1X,E15.9E2,1X))'       
             
	PRINT*
	WRITE(*,37) 'LAYER','OMEGA_CLD',' TAU_CLD ','PF_CLD ', &
		     'DELTA_M_CLD ','G1_CLD ','G2_CLD ','G_CLD  '
	DO I=1,NUMLAY
	  WRITE(*,33) I,OMEGA_CLD(I),TAU_CLD(I),PF_CLD(I), &
		       DELTA_M_CLD(I),G1_CLD(I),G2_CLD(I),G_CLD(I)           
	END DO
      
	PRINT*
	WRITE(*,*) 'CHI_CLD ='
	DO I=1,NUMLAY
         WRITE(*,*) 'LAYER',I
	  WRITE(*,FORM2) (CHI_CLD(MOMENT,I),MOMENT=0,4*N)            
	END DO		     

	PRINT*
	WRITE(*,37) 'LAYER','OMEGA_AERO','TAU_AERO ','PF_AERO', &
		     'DELTA_M_AERO','G1_AERO','G2_AERO','G_AERO '
	DO I=1,NUMLAY
	  WRITE(*,33) I,OMEGA_AERO(I),TAU_AERO(I),PF_AERO(I), &
		       DELTA_M_AERO(I),G1_AERO(I),G2_AERO(I),G_AERO(I)
	END DO
       
	PRINT*
	WRITE(*,*) 'CHI_AERO ='
	DO I=1,NUMLAY
         WRITE(*,*) 'LAYER',I       
	  WRITE(*,FORM2) (CHI_AERO(MOMENT,I),MOMENT=0,4*N)            
	END DO		        

33	FORMAT(2X,I3,2X,2(2X,E15.9E2,1X),4X,I2,10X,I1,8X,3(2X,E15.9E2,1X))
37	FORMAT(1X,A5,1X,2(5X,A9,4X),2X,A7,2X,1(1X,A12,1X),3(6X,A7,5X), &
	       1X,1(1X,A12,1X))       
       
	PRINT*
	IF (SURF == 0) THEN
	  PRINT*,'SURFACE = NONE'
	ELSE IF (SURF == 1) THEN
	  PRINT*,'SURFACE = LAMBERTIAN'
	  PRINT*,'ALBEDO = ',ALBEDO
	ELSE
	  PRINT*,'SURFACE = NONLAMBERTIAN'
	  PRINT*
	  WRITE(*,*) 'CHI_SURF ='
	  WRITE(*,FORM2) (CHI_SURF(MOMENT),MOMENT=0,4*N)
	END IF

	PRINT*
	PRINT*,'WVN_FLAG = ',WVN_FLAG
	PRINT*,'KFLAG = ',KFLAG

	IF ((WVN_FLAG == 'NEW').AND.(KFLAG == 0)) THEN
	  PRINT*,'WVN_FLAG AND KFLAG SET TO COMPUTE RADIANCES ONLY'
	ELSE IF ((WVN_FLAG == 'NEW').AND.(KFLAG == 1)) THEN
	  PRINT*,'WVN_FLAG AND KFLAG SET TO COMPUTE RADIANCES AND ' // &
	         'SAVE LAYERS FOR NEXT CALL'
	ELSE IF ((WVN_FLAG == 'OLD').AND.(KFLAG == 1)) THEN
	  PRINT*,'WVN_FLAG AND KFLAG SET TO COMPUTE RADIANCES ' // &
		  'USING SAVED LAYERS'
	END IF

	PRINT*
	PRINT*,'FLUXES       = ',FLUXES
	PRINT*,'WARNING_FLAG = ',WARNING_FLAG
	PRINT*,'P_IO         = ',P_IO
	PRINT*,'R_IO         = ',R_IO
	PRINT*,'FLUX_IO      = ',FLUX_IO
	PRINT*,'L_IO         = ',L_IO
	PRINT*,'B_IO         = ',B_IO
	PRINT*,'INV_IO       = ',INV_IO

10     FORMAT(2X,I3,5X,E15.9E2)
20     FORMAT(1X,I3,3X,3(3X,E15.9E2))
30     FORMAT(1X,I3,3X,2(3X,E15.9E2))
32     FORMAT(2X,I3,5X,I1,7X,A7,9X,I1)

40     FORMAT(1X,A5,7X,A8)
45     FORMAT(2X,A5,6X,A10)
65     FORMAT(2X,A5,2(5X,A11,2X))

	PRINT*
	PRINT*,'DERIVED LAYER OPTICAL DEPTHS FROM INPUT PARAMETERS:'
	PRINT*,'(FOR INFORMATION ONLY)'
	WRITE(*,40) 'LAYER','  TAU   '
	ALLOCATE ( TAU(NUMLAY) )
	DO I=1,NUMLAY
	  TAU(I) = TAU_CLD(I) + TAU_AERO(I) + TAU_GAS(I) + TAU_RAY(I)
	  WRITE(*,10) I,TAU(I)
	END DO
	DEALLOCATE (TAU)

	PRINT*
	PRINT*,'*** END OF INPUTS ***'
	  
	RETURN
	END SUBROUTINE DISPLAY_INPUT
	
!*******************************************************************************
!*******************************************************************************
	FUNCTION STRING2(NUM)

!INPUT : NUM
!OUTPUT: STRING

!THIS FUNCTION CONVERTS AN INTEGER WITH UP TO 5 DIGITS TO A CHARACTER STRING
!(WITH POSSIBLE TRAILING BLANKS - USE THIS FUNCTION WITH INTRINSIC FUNCTION
!"TRIM" TO ELIMINATE ANY TRAILING BLANKS)

!PROGRAMMER: MATT CHRISTI

!DATA DICTIONARY********************************************************
!
!
!***********************************************************************

!INTRINSIC SUBPROGRAMS USED ********************************************
!	DBLE,MOD,INT,ADJUSTL
!***********************************************************************

!EXTERNAL SUBPROGRAMS USED *********************************************
!	NONE
!***********************************************************************

	IMPLICIT NONE
!INPUT VARIABLES
	INTEGER :: &
         NUM
!OUTPUT VARIABLES
	CHARACTER(LEN=5) :: &
         STRING2
!INTERNAL VARIABLES
	INTEGER :: &
         I,N1,DIGIT
	CHARACTER(LEN=1), DIMENSION(0:10) :: &
         NUM_LET
	CHARACTER(LEN=1), DIMENSION(5) :: &
         STR
	DATA NUM_LET/'0','1','2','3','4','5','6','7','8','9',' '/


!START PROGRAM
!	PRINT*
!	PRINT*,'ENTERING STRING'

!CONVERT DIGITS OF INTEGER TO STRING
	N1 = NUM
	DO I=5,1,-1
	  DIGIT = INT(MOD(DBLE(N1),10.0D0))
	  STR(I) = NUM_LET(DIGIT)
	  N1 = (N1 - DIGIT)/10
	END DO

!CHANGE LEADING ZEROS OF INTEGER PORTION TO BLANKS
	I = 0
	DO
	  I = I + 1
	  IF ((STR(I) == '0').AND.(I /= 5)) THEN
	    STR(I) = ' '
	  ELSE
	    EXIT
	  END IF
	END DO

!PREPARE OUTPUT
	STRING2 = STR(1)//STR(2)//STR(3)//STR(4)//STR(5)
	STRING2 = ADJUSTL(STRING2)

!END PROGRAM
!	PRINT*
!	PRINT*,'LEAVING STRING'

	RETURN
	END FUNCTION STRING2

!**********************************************************************
!**********************************************************************

END MODULE radiant2s_driver