interpolate_libraries.f90

Go to the documentation of this file.

 module interpolate_libraries
 
    implicit none
    
    private
    
    public :: libraryinterpolate, scatangle, interpolate_wind_speed, lib_clean, lib_init ! GetReflForGivenWindSpeedIce, GetReflForGivenWindSpeedWater
    
!   REAL*8                  ::  CONST_PI,   RAD_PER_DEG , DEG_PER_RAD   
    real, parameter, dimension(3) :: wind_speed = (/3.0, 7.0, 15.0/)
    integer :: NUM_OF_TAUS, NUM_OF_WAVELENGTHS_RFTF
    
    real :: RLphase
    real, dimension(:), allocatable :: aeroPhase
    real, dimension(:,:), allocatable :: fOmega_water, fOmega_ice
    real, dimension(:,:,:), allocatable :: tauC_source_water, tauC_source_ice
    
    REAL,ALLOCATABLE::multiScatRefl_water_lamb(:,:,:)
    REAL,ALLOCATABLE::multiScatRefl_ice_lamb(:,:,:)
 
    REAL,ALLOCATABLE::multiScatRefl_water_CM(:,:,:,:)
    REAL,ALLOCATABLE::multiScatRefl_ice_CM(:,:,:,:)
 
! single-scattering arrays  
    real, allocatable :: fPhase_lamb(:,:), tauC(:,:)
    REAL, ALLOCATABLE :: gsquare(:), aero_denom(:)
 
    integer :: lib_call_counter
 
        REAL,ALLOCATABLE::phaseFunVals_water(:,:), phaseFunVals_ice(:,:), ssRefl_water(:,:,:), ssRefl_ice(:,:,:)
        real, dimension(:), allocatable :: LUT_TAU_VALUES
 
    real:: COS_SCAT, COS_VZA, COS_SZA
 
 
! ** extra SS stuff
 
    real:: prev_ratio 
    
    REAL, ALLOCATABLE ::  fPhase_liq(:,:),tauLR_liq(:,:),EXP0_liq(:,:),EXP1_liq(:,:)
    REAL, ALLOCATABLE ::  fPhase_ice(:,:),tauLR_ice(:,:),EXP0_ice(:,:),EXP1_ice(:,:)
    REAL, dimension(:), allocatable ::ocean_aero_Omega, ocean_aero_phase
 
! ** end extra SS stuff 
 
 
 contains
 
 
 subroutine lib_init
    
    use libraryarrays
 
    integer :: nwl, iwl, itau, ntau, i
    real, dimension(:,:), allocatable :: fqe_water, fqe_ice
 
    lib_call_counter = 0
    
    nwl = number_wavelengths
    ntau = number_taus
    allocate(aerophase(nwl))
 
    allocate(fqe_water(number_wavelengths, number_waterradii), fomega_water(number_wavelengths, number_waterradii))
    allocate(fqe_ice(number_wavelengths, number_iceradii), fomega_ice(number_wavelengths, number_iceradii))
 
! ** extra SS stuff 
 
    ALLOCATE(fphase_liq(1:nwl,1:number_waterradii),taulr_liq(1:nwl,1:number_waterradii), &
                exp0_liq(1:nwl,1:number_waterradii),exp1_liq(1:nwl,1:number_waterradii))
    ALLOCATE(fphase_ice(1:nwl,1:number_iceradii),taulr_ice(1:nwl,1:number_iceradii), &
                exp0_ice(1:nwl,1:number_iceradii),exp1_ice(1:nwl,1:number_iceradii))
 
    allocate(ocean_aero_omega(nwl), ocean_aero_phase(nwl))
   ! WDR make re-entrant
    if( .not. allocated(aero_denom) ) allocate(aero_denom(nwl))
 
    do i=1, nwl
        if (aerosol_tau(i) > 0.) then 
            ocean_aero_omega(i) = ((aerosol_tau(i) -0.1) * 1.0 + 0.1 * aerosol_ssa(i))/aerosol_tau(i)
        else
            ocean_aero_omega(i) = 0. 
        endif
    end do
 
    aero_denom = (aerosol_tau -0.1) + 0.1 * aerosol_ssa 
 
    prev_ratio = -999.
 
 
! ** end extra SS stuff 
 
 
 
    fqe_water(1,:) = 1.0
    do iwl = 2,nwl
        fqe_water(iwl,:) = extinction_water(iwl,:)/extinction_water(1,:)
    enddo
    
    fomega_water = truncation_factor_water * singlescattering_water              !ssalbedo * ftrunc= (omega*f)
 
    fqe_ice(1,:) = 1.0
    do iwl = 2,nwl
        fqe_ice(iwl,:) = extinction_ice(iwl,:)/extinction_ice(1,:)
    enddo
    
    fomega_ice = truncation_factor_ice * singlescattering_ice              !ssalbedo * ftrunc= (omega*f)
 
    allocate(tauc_source_water(number_taus+1, number_wavelengths, number_waterradii))
    allocate(tauc_source_ice(number_taus+1, number_wavelengths, number_iceradii))
 
    tauc_source_water(1,:,:) = 0.0
        tauc_source_ice(1,:,:) = 0.0
 
    do itau = 2, number_taus+1
        tauc_source_water(itau,:,:) = (library_taus(itau-1) * fqe_water) * (1.0 - fomega_water)   !Qe sacling, tau scaling, phase func truncation
        tauc_source_ice(itau,:,:) = (library_taus(itau-1) * fqe_ice) * (1.0 - fomega_ice)   !Qe sacling, tau scaling, phase func truncation
    end do
 
    deallocate(fqe_water, fqe_ice)
 
!  ice_int(ntau+1,number_wavelengths,number_iceradii) )     
 
    allocate( multiscatrefl_water_lamb(ntau+1,number_wavelengths,number_waterradii), &
              multiscatrefl_ice_lamb(ntau+1,number_wavelengths,number_iceradii), &
              multiscatrefl_water_cm(ntau+1,number_wavelengths,number_waterradii, 3), &
              multiscatrefl_ice_cm(ntau+1,number_wavelengths,number_iceradii, 3) )
        
    allocate( fphase_lamb(nwl, number_waterradii))
            
    allocate( gsquare(nwl))
        
 
 
    gsquare = aerosol_asym*aerosol_asym
    aerophase = 0.0
 
      num_of_taus = number_taus + 1
      num_of_wavelengths_rftf = 2
     
 
      allocate(lut_tau_values(num_of_taus))
       lut_tau_values(1) = 0.
       lut_tau_values(2:num_of_taus) = library_taus(1:number_taus)
 
        ALLOCATE(ssrefl_water(num_of_taus,number_wavelengths,number_waterradii), &
                 ssrefl_ice(num_of_taus,number_wavelengths,number_iceradii), &
                 phasefunvals_water(number_wavelengths,number_waterradii),  &
                 phasefunvals_ice(number_wavelengths,number_iceradii))
 
 
 end subroutine lib_init
 
 
 subroutine lib_clean
 
    deallocate(aerophase)
    deallocate(fomega_water, fomega_ice)
 
    deallocate(tauc_source_water, tauc_source_ice)
    
    deallocate(multiscatrefl_water_lamb, multiscatrefl_ice_lamb)
    deallocate(multiscatrefl_water_cm, multiscatrefl_ice_cm)
    
    deallocate(fphase_lamb)
    deallocate(gsquare)
    deallocate(lut_tau_values) 
        DEALLOCATE(phasefunvals_water, phasefunvals_ice, ssrefl_water, ssrefl_ice)
 
! ** extra SS stuff 
    deallocate(fphase_liq, taulr_liq, exp0_liq, exp1_liq)
    deallocate(fphase_ice, taulr_ice, exp0_ice, exp1_ice)
    deallocate(ocean_aero_omega, ocean_aero_phase)
! ** end extra SS stuff 
 
 
 end subroutine lib_clean
 
 
 
 subroutine libraryinterpolate(local_solarzenith, &
                              local_sensorzenith, &
                              local_relativeazimuth, &
                              local_scatangle, &
                              local_wind_speed, &
                              wind_speed_only, interp_MS, interp_SS, &
                              debug, &
                              status, i, j)
  !
  !  so, this interpolates library values to the current conditions of
  !  1st 5 inputs
  !    What it interpolates?
   use core_arrays
   use generalauxtype
   use libraryarrays
   use libraryinterpolates
   use science_parameters
   
   implicit none
 
   real, intent(in)   :: local_solarzenith, local_sensorzenith, local_relativeazimuth, local_wind_speed
   real, intent(in) :: local_scatangle
   logical,intent(in)   :: debug
   integer, intent(in) ::  i, j
   integer,intent(out)  :: status
   logical, intent(in) :: wind_speed_only, interp_ms, interp_ss
 
    real :: dtheta
    integer :: ianghi, ianglow
 
   status = 0
   
   
!       CONST_PI    =   4.d0 * ATAN(1.d0) !value of PI
!   RAD_PER_DEG = CONST_PI/180.d0   !radians per degree 
!   DEG_PER_RAD = 180.d0/CONST_PI       !degrees per radian
 
    cos_vza = cos(local_sensorzenith * d2r)
    cos_sza = cos(local_solarzenith * d2r)
    cos_scat  = cos(local_scatangle * d2r) 
!   NUM_OF_TAUS = number_taus + 1
!   NUM_OF_WAVELENGTHS_RFTF = 2
 
! For Lambertian libraries: land, also ocean when 0.86um saturates.
    if (.not. cox_munk) then 
    
        int_reflectance_ice = 0.
        int_reflectance_water = 0.
 
        call getrefl(cos_sza, cos_vza, local_relativeazimuth, local_scatangle, &
                    int_reflectance_water(:,:,:), int_reflectance_ice(:,:,:), interp_ms, interp_ss, status)
      if( status == 5 ) return  ! outside table
 
        if (interp_ms) then 
 
            int_reflectance_water_sdev = 0.
            int_reflectance_ice_sdev = 0.
 
         call getsdevrefllamb(cos_sza, cos_vza, local_relativeazimuth,&
                        int_reflectance_water_sdev(:,:,:), &
                        int_reflectance_ice_sdev(:,:,:), &
                        status)
         if( status == 5 ) return  ! out of table space
            call setup_emissivity_flux(cos_sza, library_fluxsolarzenith, ianghi, ianglow, dtheta)           
 
            call interpolatefluxes(cos_sza, library_fluxsolarzenith, flux_down_ice_solar, int_fluxdownice_solar, &
                        ianghi, ianglow, dtheta, status)
            call interpolatefluxes(cos_sza, library_fluxsolarzenith, flux_down_water_solar, int_fluxdownwater_solar, &
                        ianghi, ianglow, dtheta, status)
        
            call setup_emissivity_flux(cos_vza, library_fluxsensorzenith, ianghi, ianglow, dtheta)          
 
            call interpolatefluxes(cos_vza, library_fluxsensorzenith, flux_down_ice_sensor, int_fluxdownice_sensor, &
                        ianghi, ianglow, dtheta, status)
            call interpolatefluxes(cos_vza, library_fluxsensorzenith, flux_up_ice_sensor, int_fluxupice_sensor, &
                        ianghi, ianglow, dtheta, status)
    
            call interpolatefluxes(cos_vza, library_fluxsensorzenith, flux_down_water_sensor, int_fluxdownwater_sensor, &
                        ianghi, ianglow, dtheta, status)
            call interpolatefluxes(cos_vza, library_fluxsensorzenith, flux_up_water_sensor, int_fluxupwater_sensor, &
                        ianghi, ianglow, dtheta, status)
 
        endif
 
    else
 
        int_reflectance_ice = 0.
        int_reflectance_water = 0.
 
        call getreflforgivenwindspeed(cos_sza, cos_vza, local_relativeazimuth, &
          local_scatangle, cos_scat, local_wind_speed, &
          int_reflectance_water(:,:,:), int_reflectance_ice(:,:,:), &
          wind_speed_only, interp_ms, interp_ss, status)
      if( status == 5 ) return
 
        if (interp_ms) then                                 
                                        
            int_reflectance_water_sdev = 0.
            int_reflectance_ice_sdev = 0.
                            
        call getsdevreflforgivenwindspeed(cos_sza, cos_vza, &
           local_relativeazimuth, local_wind_speed, &
           int_refl_water_sdev_wspeed(:,:,:,:), &
           int_reflectance_water_sdev(:,:,:), &
           int_refl_ice_sdev_wspeed(:,:,:,:), &
           int_reflectance_ice_sdev(:,:,:), &          
           wind_speed_only, status)
         if( status == 5 ) return  !  point beyond table range
 
            call  setup_emissivity_flux(cos_vza, library_sensor_zenith, ianghi, ianglow, dtheta) 
 
            call interpemissforagivenwspeed(cos_vza, local_wind_speed, surface_emissivity_water,&
                            int_surface_emis_water_wspeed, int_surface_emissivity_water, &
                            wind_speed_only, ianghi, ianglow, dtheta, status)
            call interpemissforagivenwspeed(cos_vza, local_wind_speed, surface_emissivity_ice, &
                            int_surface_emis_ice_wspeed, int_surface_emissivity_ice, &
                            wind_speed_only, ianghi, ianglow, dtheta,  status)
            call interpemissforagivenwspeed(cos_vza, local_wind_speed, cloud_emissivity_water, &
                            int_cloud_emis_water_wspeed, int_cloud_emissivity_water, &
                            wind_speed_only,  ianghi, ianglow, dtheta, status)
            call interpemissforagivenwspeed(cos_vza, local_wind_speed, cloud_emissivity_ice, &
                            int_cloud_emis_ice_wspeed, int_cloud_emissivity_ice, &
                            wind_speed_only,  ianghi, ianglow, dtheta, status)
 
! These quantities are not used for Collection 6 calculations
            call interpemissforagivenwspeed(cos_vza, local_wind_speed, surface_emissivity_water_sdev,&
                            int_surface_emis_water_sdev_wspeed, int_surface_emissivity_water_sdev, &
                            wind_speed_only,  ianghi, ianglow, dtheta, status)
            call interpemissforagivenwspeed(cos_vza, local_wind_speed, surface_emissivity_ice_sdev, &
                            int_surface_emis_ice_sdev_wspeed, int_surface_emissivity_ice_sdev, &
                            wind_speed_only,  ianghi, ianglow, dtheta, status)
 
            call interpemissforagivenwspeed(cos_vza, local_wind_speed, cloud_emissivity_water_sdev, &
                        int_cloud_emis_water_sdev_wspeed, int_cloud_emissivity_water_sdev, &
                        wind_speed_only,  ianghi, ianglow, dtheta, status)
            call interpemissforagivenwspeed(cos_vza, local_wind_speed, cloud_emissivity_ice_sdev, &
                        int_cloud_emis_ice_sdev_wspeed, int_cloud_emissivity_ice_sdev, &
                        wind_speed_only,  ianghi, ianglow, dtheta, status)
        endif
 
    endif
 
 
 
 
end subroutine libraryinterpolate
 
REAL FUNCTION scatangle(solarAng,viewAng,relAzm)
!................................................................................
!
!................................................................................
! !F90
!
! !DESCRIPTION:     
!     This function computes the scattering angle for a given sun-satellite geometry! 
!
! !INPUT PARAMETERS:
!   solarAng : solar angle in degrees (scalar)
!   viewAng  : view angle in degrees (scalar)
!   relAzm   : rel azimuth (from the principle plane of the SUN) (scalar)
!
!
! !OUTPUT PARAMETERS:
!   scatAngle : scattering angle in degrees (scalar)
! 
!................................................................................
! !Revision History:
!       Revision 1.0  2009/06/03 GW
!       Initial integration of the standalone version. No changes from the standalone.
!       
!................................................................................
! !Team-unique Header:
!             Cloud Retrieval Group, NASA/GSFC
!
! !PROGRAMMER:
!
!             Nandana Amarasinghe (SSAI)
!             Climate and Radiation Branch
!             NASA Goddard Space Flight Center
!             Greenbelt, Maryland 20771, U.S.A.
!
!             Gala Wind (SSAI)
!             Climate and Radiation Branch
!             NASA Goddard Space Flight Center
!             Greenbelt, Maryland 20771, U.S.A.
!
!*******************************************************************************
! !END
    use science_parameters, only:d2r
 
    IMPLICIT NONE
    !CALLING ARGUMENTS
    REAL, INTENT(IN):: solarang,viewang,relazm
 
    !LOCAL VARIABLES
    REAL::solarmu,viewmu !,pi,torad,todeg
 
!   pi      = acos(-1.d0 )
!   torad   = pi / 180.d0
!   todeg   = 180.d0 / pi
 
    solarmu = cos(solarang * d2r)
    viewmu = cos(viewang * d2r)
 
    scatangle= -solarmu*viewmu + sqrt( ( 1.- solarmu**2 )* &
                      ( 1.- viewmu**2 ) )*cos( relazm*d2r)
 
 
    IF(scatangle < -1.0)scatangle = -1.0 !guard against round off error 
    IF(scatangle > 1.0)scatangle = 1.0
        
                     
    scatangle = acos(scatangle) / d2r !* todeg
 
END FUNCTION scatangle
 
 
SUBROUTINE getphasefunctionvalues(scatAngle, scatAngleArray, num_angles, phaseFunArray, phaseNormConst, & 
                                                                     phaseFunVals, ierror)
!................................................................................
!
!................................................................................
! !F90
!
! !DESCRIPTION:     
!   This program interpolates the phase fun value array for a given scattering angle an returns
!   the corresponding values for particular lambda and re (see the inputs and outputs)!
! !INPUT PARAMETERS:
    !   scatAngle       :   scattering angle in degrees (scalar)
    !   scatAngleArray  :   1-D array of sacattering angles for which phase fun vals are precomputed
    !   phaseFunArray   :   3-D array of phase fun vals: 1st DIM = num of angles, 
    !                       2nd DIM = num of wavelengths, 3rd DIM = num of REs
    !   phaseNormConst  :   Normalization factor for the phase function. If it is normalized to 1
    !                       then it is OK. But, this number (= pmom(0)) from the dfit routines  has
    !                       been saved and included in the library file in adition to d_fit coeffs.
!
!
! !OUTPUT PARAMETERS:
    !   phaseFunVals    :   2-D array of interpolated phase fun values at scatAngle
    !                       1st dim = num of wavelengths, 2nd dim = num of rE's
    !   ierror          :   integer value
    !                       ierror  = 0     success, 
    !                               = 3     invalid data
!................................................................................
! !Revision History:
!       Revision 1.0  2009/06/03 GW
!       Initial integration of the standalone version. No changes from the standalone.
!       
!................................................................................
! !Team-unique Header:
!             Cloud Retrieval Group, NASA/GSFC
!
! !PROGRAMMER:
!
!             Nandana Amarasinghe (SSAI)
!             Climate and Radiation Branch
!             NASA Goddard Space Flight Center
!             Greenbelt, Maryland 20771, U.S.A.
!
!             Gala Wind (SSAI)
!             Climate and Radiation Branch
!             NASA Goddard Space Flight Center
!             Greenbelt, Maryland 20771, U.S.A.
!
!*******************************************************************************
! !END
 
 use libraryarrays, only: number_wavelengths
 use core_arrays, only: platform_name
    IMPLICIT NONE
    !Arguments
    REAL,INTENT(IN)::scatAngle, scatAngleArray(:),phaseFunArray(:,:,:),phaseNormConst(:,:)
    REAL,INTENT(OUT)::phaseFunVals(:,:)
    integer, intent(in):: num_angles
    INTEGER,INTENT(OUT)::ierror
 
    !local variables
    INTEGER:: iAng, iAngLow, iAngHi, nwl
    REAL::dtheta
    CHARACTER(150)::tempstr
 
    !start execution
    
    if (platform_name == "AVIRIS" .or. platform_name == "RSP" .or. platform_name == "AHI" .or. platform_name == "EPIC" .or. &
        platform_name == "SSFR") then 
        nwl = number_wavelengths
    else
        nwl = number_wavelengths-1
    endif
        
 
    ierror = 0 !start with no error
    ianglow = 1
    ianghi = num_angles
    IF(scatangle > scatanglearray(ianghi))THEN
        ianglow = ianghi-1 !no need to do the search; do extrapolation
        
    ELSEIF(scatangle < scatanglearray(ianglow))THEN
        ianghi = ianglow+1
    ELSE
        DO                                  !do a bisection search
            IF((ianghi - ianglow) == 1)EXIT
            iang = (ianghi + ianglow)/2
            IF(scatanglearray(iang) < scatangle)THEN
            ianglow = iang
            ELSE
                ianghi = iang
            ENDIF
        ENDDO
    ENDIF
 
 
    dtheta = scatanglearray(ianghi) - scatanglearray(ianglow)
    !print*,iAngHi,iAngLow,dtheta,scatAngleArray(iAngHi)-scatAngle,scatAngle - scatAngleArray(iAngLow)
!   IF(dtheta ==0)THEN
!       ierror = 3          !invalid data read
!       RETURN
!   ENDIF
    
    phasefunvals(1:nwl,:) = ( phasefunarray(ianglow,1:nwl,:)  * (scatanglearray(ianghi)-scatangle) + &
                    phasefunarray(ianghi,1:nwl,:)  * (scatangle - scatanglearray(ianglow) ))/dtheta 
 
 
    phasefunvals(1:nwl,:) = phasefunvals(1:nwl,:)/phasenormconst(1:nwl,:)
 
    RETURN
 
END SUBROUTINE getphasefunctionvalues
 
 
subroutine get_aero_params(cos_scatAngle, aeroG)
    use libraryarrays, only: number_wavelengths
 use core_arrays, only: platform_name
        
    real, intent(in) :: aeroG(:), cos_scatAngle
        
    INTEGER:: nwl
    
    if (platform_name == "AVIRIS" .or. platform_name == "RSP" .or. platform_name == "AHI" .or. platform_name == "EPIC" .or. &
        platform_name == "SSFR" ) then 
        nwl = number_wavelengths
    else
        nwl = number_wavelengths-2
    endif
        
    rlphase = 0.75 * (1 + cos_scatangle*cos_scatangle)
 
! it is generally advisable to avoid a real power whenever possible in order to make the code run faster
! GW: 10.25.13
    aerophase(1:nwl) = (1.0 - gsquare(1:nwl))/ &
                            (sqrt(1 + gsquare(1:nwl) &
                            - 2.0 * aerog(1:nwl) * cos_scatangle)**3)
 
end subroutine get_aero_params
 
 
SUBROUTINE interpolatefluxes(solarOrViewAng, solarAngMuArray, inFluxArray, outFluxArray, iAngHi, iAngLow, dtheta, ierror)
!................................................................................
! !F90
!
! !DESCRIPTION:     
    !This program interpolates(linear) the RF or TF array for a given solar/oview angle and returns
    !the corresponding values for particular tau, lambda and re (see the inputs and outputs)
! !INPUT PARAMETERS:
    !   solarOrViewAng      :   input angle in degrees (scalar)
    !   solarAngleMuArray   :   1-D array of solar angles for which RF/TF vals are precomputed
    !   inFluxArray         :   4-D array of flux vals: 1st DIM = num of angles, 2nd DIM = ntaus
    !                           3rdd DIM = num of wavelengths, 4th DIM = num of REs
!
!
! !OUTPUT PARAMETERS:
    !   outFluxArray        : 3-D array of interpolated RF or TF values at the input angle
    !                         1st dim = ntaus, 2nd dim = num of wavelengths, 3rdd dim = num of rE's
    !   ierror          :   integer value
    !                       ierror  = 0     success, 
    !                               = 3     invalid data
!................................................................................
! !Revision History:
!       Revision 1.0  2009/06/03 GW
!       Initial integration of the standalone version. No changes from the standalone.
!       
!................................................................................
! !Team-unique Header:
!             Cloud Retrieval Group, NASA/GSFC
!
! !PROGRAMMER:
!
!             Nandana Amarasinghe (SSAI)
!             Climate and Radiation Branch
!             NASA Goddard Space Flight Center
!             Greenbelt, Maryland 20771, U.S.A.
!
!             Gala Wind (SSAI)
!             Climate and Radiation Branch
!             NASA Goddard Space Flight Center
!             Greenbelt, Maryland 20771, U.S.A.
!
!*******************************************************************************
! !END
    IMPLICIT NONE
    !Arguments
    REAL,INTENT(IN):: solarOrViewAng, solarAngMuArray(:),inFluxArray(:,:,:,:)
    REAL,INTENT(OUT)::outFluxArray(:,:,:)
    INTEGER,INTENT(OUT)::ierror
    integer, intent(in) :: iAngLow, iAngHi
    real, intent(in) :: dtheta
 
    !local variables
    real*8::umu
 
    real :: diff_one_way, diff_other_way
 
    !start esecution
    ierror = 0
    
    umu  = solarorviewang 
 
    diff_one_way = solarangmuarray(ianghi)- umu
    diff_other_way = umu - solarangmuarray(ianglow)
 
    outfluxarray(:,:,:) = 0.
    outfluxarray(:,:,:) =  ( influxarray(ianglow,:,:,:)  * diff_one_way + &
                                influxarray(ianghi,:,:,:)  * diff_other_way ) / dtheta
 
   RETURN
 
END SUBROUTINE interpolatefluxes
 
    SUBROUTINE calcsingscatpartofreflectance(phaseFunVals_water, phaseFunVals_ice, &
                            theta, theta0, ssRefl_water, ssRefl_ice, ierror)
 
!................................................................................
! !F90
!
! !DESCRIPTION:     
    !Calculates the single scattering part of the reflectance according to Nakajima/Tanaka method
    !implemented in DISORT Version2. We had to implement this way, because of the way single 
    !scattering part is takent out in DISORT. Please refer to the DISORT2 documentation
    !Formula Used here: refl = omega * PC * {1.0 - EXP[-tauc *(1/mu + 1/mu0)]}
    !                   refl = refl/(4.0 * (mu + mu0))
    !                   where tauc = tau * [Qe(lambda,re)/Qe(CH1,re)] * (1.0 - ftrunc * omega)
    !                          PC  = Ptrue / (1 - ftrunc * omega)
    !                   also remember refl_lUT = (I_disort) * pi/( mu0* solarflux)
    !                   so that multiplication is already taken care of here, no need to do it
    !                   outside of this routine
! !INPUT PARAMETERS:
    !   PhaseFunVals : phase function calculated for all the wavelengths and effective radii
    !                  DIMS[1:nlambda,1:nradii] (look at the program getPhaseFunValues.f90)
    !   ssAlbedo     : single scattering albedo DIMS[1:nlambda,1:nradii]
    !   extCoeff     : extinction coefficient  DIMS[1:nlambda,1:nradii]
    !   dfitF        : phase function truncation factor in dfit method, DIMS[1:nlambda,1:nradii]
    !   tauVals      : opticat thickness values in the LUT, DIMS[1:ntau]
    !   theta        : view angle in degrees
    !   theta0       : solar angle in degrees
!
!
! !OUTPUT PARAMETERS:
    !   ssRefl       : single scattering part of the reflectance , here no need to input the beam
    !                 strength, since we calculate the reflectance.
    !   ierror          :   integer value
    !                       ierror  = 0     success, 
    !                               = 2     memory allocation error
    !                               = 3     invalid data
!................................................................................
! !Revision History:
!       Revision 1.0  2009/06/03 GW
!       Initial integration of the standalone version. No changes from the standalone.
!       
!................................................................................
! !Team-unique Header:
!             Cloud Retrieval Group, NASA/GSFC
!
! !PROGRAMMER:
!
!             Nandana Amarasinghe (SSAI)
!             Climate and Radiation Branch
!             NASA Goddard Space Flight Center
!             Greenbelt, Maryland 20771, U.S.A.
!
!             Gala Wind (SSAI)
!             Climate and Radiation Branch
!             NASA Goddard Space Flight Center
!             Greenbelt, Maryland 20771, U.S.A.
!
!*******************************************************************************
! !END
    use libraryarrays
 use core_arrays, only: platform_name
 
    REAL, INTENT(IN)  :: phaseFunVals_water(:,:), phaseFunVals_ice(:,:)
    REAL, INTENT(IN)  :: theta, theta0
    REAL, INTENT(OUT) :: ssRefl_water(:,:,:), ssRefl_ice(:,:,:)
    INtEGER,INTENT(OUT)::ierror
    
    
    INTEGER:: ntau, nwl, nre, itau,iwl, allocateStatus, ii,jj
    REAL ::  MUPLUSMU0,ratio, mytau
    
    ntau = number_taus !size(tauVals)
 
    if( platform_name == "AVIRIS" .or. platform_name == "RSP" .or. platform_name == "AHI" .or. platform_name == "EPIC" .or. &
        platform_name == "SSFR" ) then 
        nwl = number_wavelengths + 1 ! we have to make sure not to skip the last channel
    else
        nwl = number_wavelengths
    endif
                
 
    muplusmu0 = theta + theta0
    ratio = muplusmu0/(theta * theta0)
    
 
    fphase_lamb = phasefunvals_water/ (1.0 - fomega_water)   !add something here to safeguard div by zero
    
    ssrefl_water = 0.0
    ssrefl_ice = 0.0
    
    do itau = 1, ntau
    
        do ii=1, nwl-1
            do jj=1, number_waterradii
 
                mytau = tauc_source_water(itau+1,ii,jj)*ratio
 
                if (mytau == 0) then 
                    exp1_liq(ii,jj) = 1.
                else if (mytau > 10.) then 
                    exp1_liq(ii,jj) = 0.
                else  
                    exp1_liq(ii,jj) = exp(-mytau)
 
!                   EXP1_liq(ii,jj) = 1. + (mytau ) + (mytau )**2/2. + (mytau )**3/6. + &
!                                                               (mytau )**4/24. + (mytau )**5/120. + mytau**6/720. 
!                   EXP1_liq(ii,jj) = 1. / EXP1_liq(ii,jj)
                endif
    
            ssrefl_water(itau,ii,jj) = (singlescattering_water(ii,jj)* fphase_lamb(ii,jj) * (1 - exp1_liq(ii,jj))) / (4.0*muplusmu0) 
            
            end do
        end do
    
    
    enddo
    
    fphase_lamb(:, 1:number_iceradii) = phasefunvals_ice/ (1.0 - fomega_ice)   !add something here to safeguard div by zero
    
    do itau = 1, ntau
 
 
                do ii=1, nwl-1
                        do jj=1, number_iceradii
 
 
                            mytau = tauc_source_ice(itau+1,ii,jj)*ratio
 
                                if (mytau == 0) then
                                        exp1_ice(ii,jj) = 1.
                                else if (mytau > 10.) then
                                        exp1_ice(ii,jj) = 0.
                                else 
                                        exp1_ice(ii,jj) = exp(-mytau)
                                endif
 
                        ssrefl_ice(itau,ii,jj) = (singlescattering_ice(ii,jj)* fphase_lamb(ii,jj) * (1 - exp1_ice(ii,jj))) / (4.0*muplusmu0) 
 
                        end do
                end do
 
 
    enddo
    
    END SUBROUTINE calcsingscatpartofreflectance    
 
 
    SUBROUTINE single_scattering_calcs_ocean(phaseFunVals_liq, phaseFunVals_ice, ssAlbedo_liq, ssAlbedo_ice, &
                                                RLphase, aeroPhase, RLTau,aeroTau, aeroOmega,&
                                                  theta, theta0,ssRefl_liq, ssRefl_ice)
 
    use libraryarrays, only: number_wavelengths, number_waterradii, number_iceradii
 use core_arrays, only: platform_name
 
    REAL, INTENT(IN)  :: phaseFunVals_liq(:,:), ssAlbedo_liq(:,:),  &
                         phaseFunVals_ice(:,:), ssAlbedo_ice(:,:), &
                         RLTau(:), aeroTau(:), &
                         aeroOmega(:), RLphase, aeroPhase(:)
    REAL, INTENT(IN)  :: theta, theta0
    REAL, INTENT(OUT) :: ssRefl_liq(:,:,:), ssRefl_ice(:,:,:)
    
    
    INTEGER:: ntau, nwl, itau, iwl
    REAL :: MUPLUSMU0,ratio
    
    real :: ocean_mul(20)
    logical:: changed_ratio
    integer :: ii, jj
 
    ntau = num_of_taus
    
    if (platform_name == "AVIRIS" .or. platform_name == "RSP" .or. platform_name == "AHI" .or. platform_name == "EPIC" .or. &
        platform_name == "SSFR") then
        nwl = number_wavelengths+1 ! capture the last channel, it is not to be skipped
    else
        nwl = number_wavelengths
    endif
                                
    muplusmu0 = theta + theta0
    ratio = muplusmu0/(theta * theta0)
 
    fphase_liq = phasefunvals_liq/ (1.0 - fomega_water)   !add something here to safeguard div by zero
    fphase_ice = phasefunvals_ice/ (1.0 - fomega_ice)
    
 
#if AVIRIS_INST | RSP_INST | AHI_INST | EPIC_INST | SSFR_INST
    do iwl = 1,  number_wavelengths
#else
    do iwl = 1,  number_wavelengths-2
#endif
        ocean_aero_phase(iwl) = ( (aerotau(iwl) -0.1) * rlphase  + 0.1 * aeroomega(iwl) * aerophase(iwl) ) / aero_denom(iwl) 
        ocean_mul(iwl) = ocean_aero_phase(iwl) * ocean_aero_omega(iwl)
    end do
 
    exp1_liq = 0.
    exp1_ice = 0.
 
    DO itau = 1, ntau
 
        taulr_liq =  tauc_source_water(itau,:,:) * ratio
        taulr_ice =  tauc_source_ice(itau,:,:) * ratio
 
 
        do ii=1, nwl-1
            do jj=1, number_waterradii
 
                if (taulr_liq(ii,jj) == 0) then 
                    exp1_liq(ii,jj) = 1.
                else if (taulr_liq(ii,jj) > 10.) then 
                    exp1_liq(ii,jj) = 0.
                else 
                    exp1_liq(ii,jj) = exp(-taulr_liq(ii,jj))
                endif
 
            end do
        end do
 
                do ii=1, nwl-1
                        do jj=1, number_iceradii
 
                                if (taulr_ice(ii,jj) == 0) then
                                        exp1_ice(ii,jj) = 1.
                                else if (taulr_ice(ii,jj) > 10.) then
                                        exp1_ice(ii,jj) = 0.
                                else
                                        exp1_ice(ii,jj) = exp(-taulr_ice(ii,jj))
                                endif
 
                        end do
                end do
 
    
        ssrefl_liq(itau,:,:) = ssalbedo_liq * fphase_liq * (1. - exp1_liq)
        ssrefl_ice(itau,:,:) = ssalbedo_ice * fphase_ice * (1. - exp1_ice)
        
        exp0_liq = exp1_liq
        exp0_ice = exp1_ice
 
#if AVIRIS_INST | RSP_INST | AHI_INST | EPIC_INST | SSFR_INST
        DO iwl =1,  number_wavelengths
#else
        DO iwl = 1,  number_wavelengths-2
#endif
                       taulr_liq(iwl,:) =  (tauc_source_water(itau,iwl,:) + rltau(iwl))*ratio       
                       taulr_ice(iwl,:) =  (tauc_source_ice(itau,iwl,:) + rltau(iwl))*ratio
 
            ii = iwl
                        do jj=1, number_waterradii
 
                                if (taulr_liq(ii,jj) == 0) then
                                        exp1_liq(ii,jj) = 1.
                                else if (taulr_liq(ii,jj) > 10.) then
                                        exp1_liq(ii,jj) = 0.
                                else
                                        exp1_liq(ii,jj) = exp(-taulr_liq(ii,jj))
                                endif
 
                        end do
 
            ii = iwl
                        do jj=1, number_iceradii
 
                                if (taulr_ice(ii,jj) == 0) then
                                        exp1_ice(ii,jj) = 1.
                                else if (taulr_ice(ii,jj) > 10.) then
                                        exp1_ice(ii,jj) = 0.
                                else
                                        exp1_ice(ii,jj) = exp(-taulr_ice(ii,jj))
                                endif
 
                        end do
 
            
            ssrefl_liq(itau,iwl,:) = ssrefl_liq(itau,iwl,:) +  & 
                                        rlphase * 1.0 * (exp0_liq(iwl,:) - exp1_liq(iwl,:))
            ssrefl_ice(itau,iwl,:) = ssrefl_ice(itau,iwl,:) +  & 
                                        rlphase * 1.0 * (exp0_ice(iwl,:) - exp1_ice(iwl,:))
            
            exp0_liq(iwl,:) = exp1_liq(iwl,:)
            exp0_ice(iwl,:) = exp1_ice(iwl,:)
            
                       taulr_liq(iwl,:) =  (tauc_source_water(itau,iwl,:) + rltau(iwl) + aerotau(iwl)) * ratio
                       taulr_ice(iwl,:) =  (tauc_source_ice(itau,iwl,:) + rltau(iwl) + aerotau(iwl)) * ratio
            
            !coming in aeroTau is the (rayleightau + 0.1) from DISORT run
 
                       ii = iwl
                        do jj=1, number_waterradii
 
                                if (taulr_liq(ii,jj) == 0) then
                                        exp1_liq(ii,jj) = 1.
                                else if (taulr_liq(ii,jj) > 10.) then
                                        exp1_liq(ii,jj) = 0.
                                else
                                        exp1_liq(ii,jj) = exp(-taulr_liq(ii,jj))
                                endif
 
                        end do
 
                        ii = iwl
                        do jj=1, number_iceradii
                       
                                if (taulr_ice(ii,jj) == 0) then
                                        exp1_ice(ii,jj) = 1.
                                else if (taulr_ice(ii,jj) > 10.) then
                                        exp1_ice(ii,jj) = 0.
                                else
                                        exp1_ice(ii,jj) = exp(-taulr_ice(ii,jj))
                                endif
 
                        end do
 
 
 
            ssrefl_liq(itau,iwl,:) = ssrefl_liq(itau,iwl,:) + ocean_mul(iwl) * &
                                                            (exp0_liq(iwl,:) - exp1_liq(iwl,:))
            ssrefl_ice(itau,iwl,:) = ssrefl_ice(itau,iwl,:) + ocean_mul(iwl) * &
                                                            (exp0_ice(iwl,:) - exp1_ice(iwl,:))
        ENDDO
                 
    ENDDO
 
    
    ssrefl_liq = ssrefl_liq/(4.0 * muplusmu0)
    ssrefl_ice = ssrefl_ice/(4.0 * muplusmu0)
        
 
 
    END SUBROUTINE single_scattering_calcs_ocean
 
 
 
 
 
    !----------------------SUBROUTINE CalcSingScatPartOfOceanReflectance---------------------------------
    !Calculates the single scattering part of the reflectance according to Nakajima/Tanaka method
    !implemented in DISORT Version2. We had to implement this way, because of the way single 
    !scattering part is takent out in DISORT. Please refer to the DISORT2 documentation
    !Formula Used here: refl = omega * PC * {1.0 - EXP[-tauc *(1/mu + 1/mu0)]}
    !                   refl = refl/(4.0 * (mu + mu0))
    !                   where tauc = tau * [Qe(lambda,re)/Qe(CH1,re)] * (1.0 - ftrunc * omega)
    !                          PC  = Ptrue / (1 - ftrunc * omega)
    !                   also remember refl_lUT = (I_disort) * pi/( mu0* solarflux)
    !                   so that multiplication is already taken care of here, no need to do it
    !                   outside of this routine
    !INPUTS:
    !   PhaseFunVals : phase function calculated for all the wavelengths and effective radii
    !                  DIMS[1:nlambda,1:nradii] (look at the program getPhaseFunValues.f90)
    !   RLphase      : Rayleigh phase func value (fn of the scattering angle)
    !   aeroPhse     : mixed phase function of rayleigh layer and the aerosol layer
    !   ssAlbedo     : single scattering albedo DIMS[1:nlambda,1:nradii]
    !   extCoeff     : extinction coefficient  DIMS[1:nlambda,1:nradii]
    !   dfitF        : phase function truncation factor in dfit method, DIMS[1:nlambda,1:nradii]
    !   tauVals      : opticat thickness values in the LUT, DIMS[1:ntau]
    !   RLTau        : Rayleigh optical thickness (combined several layers upto 8km)
    !   aeroTau      : aerosol optical thickness(0.1) + rayleigh layer(lowest)
    !   aeroOmega    : ss albedo of aerosol (not the mixed, mixing will be done in SS calc routine)
    !   theta        : view angle in degrees
    !   theta0       : solar angle in degrees
    !OUTPUTS:
    !   ssRefl       : single scattering part of the reflectance , here no need to input the beam
    !                 strength, since we calculate the reflectance.
    !   ierror          :   integer value
    !                       ierror  = 0     success, 
    !                               = 2     memory allocation error
    !                               = 3     invalid data
    !-----------------------------------------------------------------------------------------------    
    
    subroutine setup_emissivity_flux(angle, angle_array, idx_hi, idx_lo, dtheta) 
    
        real, intent(in) :: angle, angle_array(:)
        real, intent(out) :: dtheta
        integer, intent(out) :: idx_hi, idx_lo
    
    
        real :: UMU
        integer :: iAng
    
        umu  = angle !COS(solarOrViewAng * RAD_PER_DEG)
        idx_lo = 1
        idx_hi = size(angle_array)
        IF(umu < angle_array(idx_lo))THEN
            idx_hi = idx_lo+1 !no need to do the search; do extrapolation
            if (idx_hi > size(angle_array)) idx_hi = idx_lo ! we don't have enough points
        ELSE
            DO                                  !do a bisection search
                IF((idx_hi - idx_lo) == 1)EXIT
                iang = (idx_hi + idx_lo)/2
                IF(angle_array(iang) < umu)THEN
                    idx_lo = iang
                ELSE
                idx_hi = iang
                ENDIF
            ENDDO
        ENDIF
 
 
    dtheta = angle_array(idx_hi) - angle_array(idx_lo)
        
    end subroutine setup_emissivity_flux
    
    
    !--------------------SUBROUTINE InterpCloudEmissForAGivenWSpeed ---------------------------------
    !This program interpolates(linear) the cloud emissivity arrays for a given wind speed and
    !returns the corresponding values for particular tau, lambda and re(see the inputs and outputs)
    !INPUTS:
    !   iphase              :   integer  1 =Water Cloud, 2 =Ice Cloud
    !   solarOrViewAng      :   input angle in degrees (scalar)
    !   wspeed              :   real scalar, wind speed in m/s
    !OUTPUTS:
    !   outFluxArray        : 3-D array ofinterpolated cloud emissivity values at the input ang
    !                         1st dim = ntaus, 2nd dim = num of wavelengths, 3rdd dim = num of rE's
    !   ierror              :   integer value
    !                           ierror  = 0     success, 
    !                               = 3     invalid data
    !-----------------------------------------------------------------------------------------------
    !                       
    SUBROUTINE interpemissforagivenwspeed(solarOrViewAng,wspeed, inEmissArray, outEmissArray_ws, outEmissArray, &
                                                wspeed_only, iAngHi, iAngLow, dtheta,ierror)
 
   use libraryarrays
 
    IMPLICIT NONE
    !Arguments
    REAL,INTENT(IN):: solarOrViewAng,wspeed, dtheta
    real, intent(in) :: inEmissArray(:,:,:,:,:)
    REAL,INTENT(OUT)::outEmissArray(:,:,:)
    real, intent(inout) :: outEmissArray_ws(:,:,:,:)
    INTEGER,INTENT(OUT)::ierror
    logical, intent(in) :: wspeed_only
    integer, intent(in) :: iAngHi, iAngLow
 
 
    real :: diff_one_way, diff_other_way
 
 
        diff_one_way = library_sensor_zenith(ianghi)- solarorviewang
        diff_other_way = solarorviewang - library_sensor_zenith(ianglow)
 
        outemissarray_ws(:,:,:,1) = (inemissarray(ianglow,:,:,:,1)  * diff_one_way + &
                                     inemissarray(ianghi,:,:,:,1) * diff_other_way ) / dtheta
        outemissarray_ws(:,:,:,2) = (inemissarray(ianglow,:,:,:,2)  * diff_one_way + &
                                     inemissarray(ianghi,:,:,:,2) * diff_other_way) / dtheta
        outemissarray_ws(:,:,:,3) = (inemissarray(ianglow,:,:,:,3)  * diff_one_way + &
                                     inemissarray(ianghi,:,:,:,3) * diff_other_way) / dtheta
 
    outemissarray = 0. 
    call interpolate_wind_speed(wspeed, outemissarray_ws, outemissarray)
 
    RETURN
 
    END SUBROUTINE interpemissforagivenwspeed
    
 
!   SUBROUTINE TriLinearInterpolationRefl(solarAng,viewAng,azmAng, &
!                                               solarMuArray, viewMuArray, relAzmArray, &
!                                               multiScatReflLUT_water, multiScatRefl_water, &
!                                               multiScatReflLUT_ice, multiScatRefl_ice,ierror)
!!................................................................................
!! !F90
!!
!! !DESCRIPTION:     
!   !perform 3-D linear interpolation on multiple scattering part of the reflectance. Instead of 
!   !indexing, a search is done to find the the indices of the cube that enclose the interpolation
!   !point.
!! !INPUT PARAMETERS:
!   !   solarAng        :   solar zenith angle in degrees (scalar)
!   !   viewAng         :   sensor zenith angle in degrees (scalar)
!   !   azmAng          :   rel azm ang in degrees (scalar)]
!   !   solarMuArray    :   array of solar zenith angles (from the LUT) 1-D array
!   !   viewMuArray     :   array of sensor zenith angles (from the LUT) 1D array
!   !   relAzmArray     :   array of rel azimuth angles (from the LUT) 1-D array
!   !   multiScatReflLUT:   multi scat part of the reflectance (from LUT) 6D array
!   !                       DIMS[1:nview, 1:nsolar,1:nazm, 1:ntau, 1:nlambda, 1:nre]
!!
!!
!! !OUTPUT PARAMETERS:
!   !   multiScatRefl   :   interpolated MS reflectance for the input geometry 3D array
!   !                       DIMS[1:ntau, 1:nlambda, 1:nre]
!   !   ierror          :   integer value
!   !                       ierror  = 0     success, 
!   !                               = 1     I/O error
!   !                               = 2     memory allocation error
!   !                               = 3     invalid data    
!!................................................................................
!! !Revision History:
!!      Revision 1.0  2009/06/03 GW
!!      Initial integration of the standalone version. No changes from the standalone.
!!      
!!................................................................................
!! !Team-unique Header:
!!             Cloud Retrieval Group, NASA/GSFC
!!
!! !PROGRAMMER:
!!
!!             Nandana Amarasinghe (SSAI)
!!             Climate and Radiation Branch
!!             NASA Goddard Space Flight Center
!!             Greenbelt, Maryland 20771, U.S.A.
!!
!!             Gala Wind (SSAI)
!!             Climate and Radiation Branch
!!             NASA Goddard Space Flight Center
!!             Greenbelt, Maryland 20771, U.S.A.
!!
!!*******************************************************************************
!! !END
!   use libraryarrays
!
!    IMPLICIT NONE
!
!   REAL,INTENT(IN):: solarAng, viewAng, azmAng
!   REAL,INTENT(IN):: solarMuArray(:), viewMuArray(:), relAzmArray(:)
!   REAL,INTENT(IN):: multiScatReflLUT_water(:,:,:,:,:,:)
!   REAL,INTENT(IN):: multiScatReflLUT_ice(:,:,:,:,:,:)
!
!   REAL,INTENT(OUT)::multiScatRefl_water(:,:,:)
!   REAL,INTENT(OUT)::multiScatRefl_ice(:,:,:)
!   INTEGER,INTENT(OUT)::ierror
!
!
!   REAL*8::prod11,prod12,prod21,prod22,vol1,vol2,vol3,vol4,vol5,vol6,vol7,vol8,volume
!   REAL*4::solarmu,viewmu
!   REAL*8::dmu0,dmu,dmu0_1,dmu0_2,dmu_1,dmu_2,dazm,dazm_1,dazm_2
!   INTEGER:: ishi,islow, ivhi,ivlow, iazhi,iazlow,n1,n2,n3
!
!   n1 = number_solarzenith !size(solarMuArray)
!   n2 = number_sensorzenith !size(viewMuArray)
!   n3 = number_relazimuth !size(relAzmArray)
!   
!   
!   ierror = 0   !no error
!   solarmu  = solarAng !COS(solarAng * RAD_PER_DEG)
!   viewmu   = viewAng !COS(viewAng  * RAD_PER_DEG)
!   
!   CALL BisectionSearch(solarmu,solarMuArray(1:n1),n1,islow,ishi)
!   
!   CALL BisectionSearch(viewmu,viewMuArray(1:n2),n2,ivlow,ivhi)
!   
!   CALL BisectionSearch(azmAng,relAzmArray(1:n3),n3,iazlow,iazhi)
!   
!   dmu = solarMuArray(ishi) - solarMuArray(islow)
!   dmu0  = viewMuArray(ivhi)  - viewMuArray(ivlow)
!   dazm = relAzmArray(iazhi) - relAzmArray(iazlow)
!   
!   dmu_1 = solarMu - solarMuArray(islow) ; dmu_2 = solarMuArray(ishi) - solarMu
!   dmu0_1  = viewMu  - viewMuArray(ivlow)  ; dmu0_2  = viewMuArray(ivhi)  - viewMu
!   dazm_1 = azmAng  -  relAzmArray(iazlow); dazm_2 = relAzmArray(iazhi) - azmAng
!       
!   volume = dmu0 * dmu * dazm
!   IF(volume == 0.0)THEN
!       ierror = 3      !invalid data 
!       RETURN
!   ENDIF
!   !print*,ivhi,ishi,iazlow
!   !print*,dmu0_1,dmu0_2,dmu_1,dmu_2,dazm_1,dazm_2
!
!   prod22 = dmu_2 * dmu0_2
!   prod12 = dmu_1 * dmu0_2
!   prod11 = dmu_1 * dmu0_1
!   prod21 = dmu_2 * dmu0_1
!
!   vol8 =  abs(prod22 * dazm_2)
!   vol7 =  abs(prod22 * dazm_1)
!   vol6 =  abs(prod12 * dazm_1)
!   vol5 =  abs(prod12 * dazm_2)
!
!   vol1 =  abs(prod11 * dazm_1)
!   vol2 =  abs(prod21 * dazm_1)
!   vol3 =  abs(prod21 * dazm_2)
!   vol4 =  abs(prod11 * dazm_2)
!   !volume = vol1 + vol2 + vol3 + vol4 + vol5 + vol6 + vol7 + vol8
!       
!   
!   multiScatRefl_water(:,:,:) = ( multiScatReflLUT_water(ivlow,islow,iazlow,:,:,:) * vol8  + &
!                            multiScatReflLUT_water(ivhi,islow,iazlow,:,:,:)  * vol3  + &
!                                multiScatReflLUT_water(ivhi,islow,iazhi,:,:,:)   * vol2  + &
!                                multiScatReflLUT_water(ivlow,islow,iazhi,:,:,:)  * vol7  + &
!                            multiScatReflLUT_water(ivlow,ishi,iazlow,:,:,:)  * vol5  + &
!                            multiScatReflLUT_water(ivhi,ishi,iazlow,:,:,:)   * vol4  + &
!                                multiScatReflLUT_water(ivhi,ishi,iazhi,:,:,:)    * vol1  + &
!                                multiScatReflLUT_water(ivlow,ishi,iazhi,:,:,:)   * vol6 )/volume
!
!   multiScatRefl_ice(:,:,:) = ( multiScatReflLUT_ice(ivlow,islow,iazlow,:,:,:) * vol8  + &
!                            multiScatReflLUT_ice(ivhi,islow,iazlow,:,:,:)  * vol3  + &
!                                multiScatReflLUT_ice(ivhi,islow,iazhi,:,:,:)   * vol2  + &
!                                multiScatReflLUT_ice(ivlow,islow,iazhi,:,:,:)  * vol7  + &
!                            multiScatReflLUT_ice(ivlow,ishi,iazlow,:,:,:)  * vol5  + &
!                            multiScatReflLUT_ice(ivhi,ishi,iazlow,:,:,:)   * vol4  + &
!                                multiScatReflLUT_ice(ivhi,ishi,iazhi,:,:,:)    * vol1  + &
!                                multiScatReflLUT_ice(ivlow,ishi,iazhi,:,:,:)   * vol6 )/volume
!
!
! 
!   RETURN
!                  
!
!   CONTAINS
!
!       SUBROUTINE BisectionSearch(theta,thetaArray,ntheta,iAngLow, iAngHi)
!   
!       INTEGER,INTENT(IN)::ntheta
!       REAL,INTENT(IN)::theta,thetaArray(1:ntheta)
!       INTEGER,INTENT(OUT)::iAngLow, iAngHi
!   
!       INTEGER:: iAng
!   
!
!       iAngLow = 1
!       iAngHi = ntheta
!
!       if (ntheta == 1) then 
!           iAngLow = 1
!           iAngHi = 1
!           return
!       endif
!
!       DO
!           IF((iAngHi - iAngLow) == 1)EXIT
!           iAng = (iAngHi + iAngLow)/2
!           IF(thetaArray(iAng) < theta)THEN
!               iAngLow = iAng
!           ELSE
!               iAngHi = iAng
!           ENDIF
!       ENDDO
!
!       if (iAngLow < 1) iAngLow = 1
!       if (ianghi < 1) ianghi = 1
!       if (iAngHi > ntheta) iAngHi = ntheta
!       if (ianglow > ntheta) ianglow = ntheta
!
!   
!       RETURN
!   
!       END SUBROUTINE BisectionSearch
!   
!   END SUBROUTINE TriLinearInterpolationRefl
 
 
    SUBROUTINE getrefl(solarAng,viewAng,azmAng,in_scat, refl_water, refl_ice, interp_MS, interp_SS, ierror)
!................................................................................
! !F90
!
! !DESCRIPTION:     
    !This subroutine outouts the Reflectance table (R(tau,lambda, re)) for a given sun-satellite
    !geometry 
! !INPUT PARAMETERS:
    !   solarAng    :   solar zenith in degrees (scalar)
    !   viewAng     :   viewing angle in degrees (scalar)
    !   azmAng      :   relative azimuth between sat and solar planes
!
!
! !OUTPUT PARAMETERS:
    !   refl    :   reflectance as a function of tau,lambda,re : 3-D Array
    !                   DIM-1 = ntau, DIM-2 = nlambda, DIM-3 = nre
    !   ierror          :   integer value
    !                       ierror  = 0     success, 
    !                               = 1     I/O error
    !                               = 2     memory allocation error
    !                               = 3     invalid data    
!................................................................................
! !Revision History:
!       Revision 1.0  2009/06/03 GW
!       Initial integration of the standalone version. Changed from the standalone
!       delivered version to not apply the albedo inside this routine. That is done in a
!       different place later in corescience. 
!       
!................................................................................
! !Team-unique Header:
!             Cloud Retrieval Group, NASA/GSFC
!
! !PROGRAMMER:
!
!             Nandana Amarasinghe (SSAI)
!             Climate and Radiation Branch
!             NASA Goddard Space Flight Center
!             Greenbelt, Maryland 20771, U.S.A.
!
!             Gala Wind (SSAI)
!             Climate and Radiation Branch
!             NASA Goddard Space Flight Center
!             Greenbelt, Maryland 20771, U.S.A.
!
!*******************************************************************************
! !END
 
    use libraryarrays
    use generalauxtype
   ! WDR to get the scan line
   use ch_xfr, only : c2_scan
    !
    REAL, INTENT(IN) :: solarAng,viewAng,azmAng, in_scat
    REAL,INTENT(OUT) :: refl_water(:,:,:), refl_ice(:,:,:)
    INTEGER,INTENT(OUT)::ierror
    logical, intent(in) :: interp_MS, interp_SS
    
    INTEGER::allocateStatus
    REAL :: scat_angle
   INTEGER::do_mgr
 
   !WDR TEMP set some loop controls for refl table output
   !INTEGER::ksen,klam,ktau,kre
 
   scat_angle = in_scat ! scatAngle(solarAng,viewAng,azmAng)
 
   if (interp_ms) then
 
   !interpolate multi scattering part of the reflection
     multiscatrefl_water_lamb = 0
     multiscatrefl_ice_lamb = 0
     call mng_ms_get_lib( solarang, viewang, azmang, 0, 0, c2_scan, &
       multiscatrefl_water_lamb, multiscatrefl_ice_lamb, ierror )
     if( ierror /= 0 ) return
    endif                                    
 
    if (interp_ss) then
 
        !get the phase function values
        CALL getphasefunctionvalues(scat_angle, phase_angles_water, number_phase_angles_water, phase_funcs_water, &
                                     phase_fun_norm_constant_water, phasefunvals_water,ierror)                                          
        CALL getphasefunctionvalues(scat_angle, phase_angles_ice, number_phase_angles_ice, phase_funcs_ice, &
                                     phase_fun_norm_constant_ice, phasefunvals_ice,ierror)                                          
 
 
!       IF(ierror /= 0)RETURN      !if there is any error no point of continuing
    
        !calculate the single scattering part of the reflection
        CALL calcsingscatpartofreflectance(phasefunvals_water, phasefunvals_ice, &
                     viewang, solarang, ssrefl_water, ssrefl_ice, ierror)
 
    endif
 
 
!   IF(ierror /= 0)RETURN      !if there is any error no point of continuing                        
                                    
 
    refl_water = ssrefl_water + multiscatrefl_water_lamb  !total reflection with asurf = 0.0    
    refl_ice = ssrefl_ice + multiscatrefl_ice_lamb  !total reflection   with asurf = 0.0    
 
    RETURN
    END SUBROUTINE getrefl
 
 
 
    !--------------------------SUBROUTINE GetReflForGivenWindSpeedWater-----------------------------
    !This subroutine outputs the Reflectance table (R(tau,lambda, re)) for a given sun-satellite
    !geometry and a wind speed
    !INPUTS :
    !   solarAng    :   solar zenith in degrees (scalar)
    !   viewAng     :   viewing angle in degrees (scalar)
    !   azmAng      :   relative azimuth between sat and solar planes
    !   wspeed      :   wind speed m/s
    !OUTPUTS :
    !   reflAsurf   :   reflectance as a function of tau,lambda,re : 3-D Array
    !                   DIM-1 = ntau, DIM-2 = nlambda, DIM-3 = nre
    !   ierror          :   integer value
    !                       ierror  = 0     success, 
    !                               = 1     I/O error
    !                               = 2     memory allocation error
    !                               = 3     invalid data    
    !ROUTINES CAlLED :
    !   ReadPhaseFunctionDataHDF
    !   ReadWaterLibraryHDF 
    !   getAllPhaseFunctionValues
    !   TriLinearInterpolationReflSearch
    !   CalcSingScatPartOfReflectance`
    !-----------------------------------------------------------------------------------------------
    SUBROUTINE getreflforgivenwindspeed(solarAng,viewAng,azmAng, in_scat, cos_scat, wspeed, reflAsurf_water, &
                                            reflAsurf_ice, wind_speed_only, interp_MS, interp_SS, ierror)
 
 
    use libraryarrays
    use libraryinterpolates
   use ch_xfr, only : c2_scan
    !
    REAL, INTENT(IN) :: solarAng,viewAng,azmAng,wspeed, in_scat, cos_scat
    REAL,INTENT(OUT) :: reflAsurf_water(:,:,:), reflAsurf_ice(:,:,:)
    INTEGER,INTENT(OUT)::ierror
    logical, intent(in) :: wind_speed_only, interp_MS, interp_SS
    
    INTEGER::allocateStatus,iuhi,iulow
    REAL :: scat_angle,deltau
   INTEGER::do_mgr
    
    scat_angle = in_scat !scatAngle(solarAng,viewAng,azmAng)
 
    if (interp_ss) then
 
        !get the phase function values
        call get_aero_params(cos_scat, aerosol_asym)
                                                                
        CALL getphasefunctionvalues(scat_angle, phase_angles_water, number_phase_angles_water, phase_funcs_water, &
                                     phase_fun_norm_constant_water, phasefunvals_water,ierror)      
 
        CALL getphasefunctionvalues(scat_angle, phase_angles_ice, number_phase_angles_ice, phase_funcs_ice, &
                                     phase_fun_norm_constant_ice, phasefunvals_ice,ierror)                                          
 
        !calculate the single scattering part of the reflection
        call single_scattering_calcs_ocean(phasefunvals_water, phasefunvals_ice, singlescattering_water, singlescattering_ice, &
                rlphase, aerophase, rayleigh_tau, aerosol_tau, aerosol_ssa, &
                    solarang,viewang, ssrefl_water, ssrefl_ice)
 
    endif
 
 
                                     
    if (interp_ms) then 
 
!     ! WDR wire up the sfc typ 1 for 3 m/s
       multiscatrefl_water_cm(:,:,:,1) = 0
       multiscatrefl_ice_cm(:,:,:,1) = 0
       call mng_ms_get_lib( solarang, viewang, azmang, 1, 0, c2_scan, &
          multiscatrefl_water_cm(:,:,:,1), multiscatrefl_ice_cm(:,:,:,1), &
          ierror )
       if( ierror /= 0 ) return
 
!
!     ! WDR wire up the sfc typ 1 for 7 m/s
        multiscatrefl_water_cm(:,:,:,2) = 0
        multiscatrefl_ice_cm(:,:,:,2) = 0
        call mng_ms_get_lib( solarang, viewang, azmang, 2, 0, c2_scan, &
          multiscatrefl_water_cm(:,:,:,2), multiscatrefl_ice_cm(:,:,:,2), &
          ierror )
       if( ierror /= 0 ) return
 
!
!     ! WDR wire up the sfc typ 1 for 15 m/s
       multiscatrefl_water_cm(:,:,:,3) = 0
       multiscatrefl_ice_cm(:,:,:,3) = 0
       call mng_ms_get_lib( solarang, viewang, azmang, 3, 0, c2_scan, &
         multiscatrefl_water_cm(:,:,:,3), multiscatrefl_ice_cm(:,:,:,3), &
         ierror )
       if( ierror /= 0 ) return
    endif
 
 
    int_reflectance_water_wspeed(:,:,:,1) = multiscatrefl_water_cm(:,:,:,1) + ssrefl_water(:,:,:)                   
    int_reflectance_water_wspeed(:,:,:,2) = multiscatrefl_water_cm(:,:,:,2) + ssrefl_water(:,:,:)
    int_reflectance_water_wspeed(:,:,:,3) = multiscatrefl_water_cm(:,:,:,3) + ssrefl_water(:,:,:)                   
 
    int_reflectance_ice_wspeed(:,:,:,1) = multiscatrefl_ice_cm(:,:,:,1) + ssrefl_ice(:,:,:)                 
    int_reflectance_ice_wspeed(:,:,:,2) = multiscatrefl_ice_cm(:,:,:,2) + ssrefl_ice(:,:,:)                 
    int_reflectance_ice_wspeed(:,:,:,3) = multiscatrefl_ice_cm(:,:,:,3) + ssrefl_ice(:,:,:)                 
 
                                                            
                                                                
    call interpolate_wind_speed(wspeed, int_reflectance_water_wspeed, reflasurf_water)
    call interpolate_wind_speed(wspeed, int_reflectance_ice_wspeed, reflasurf_ice)
 
    END SUBROUTINE getreflforgivenwindspeed
 
 
    
    subroutine interpolate_wind_speed(wspeed, data_in, data_out)
    
      use ch_xfr, only : c2_scan
        real, intent(in) :: wspeed
        real, dimension(:,:,:,:), intent(in) :: data_in
        real, dimension(:,:,:), intent(inout) :: data_out
        
        real :: deltau
        integer :: iuhi, iulow
 
        IF( (wspeed - wind_speed(1)) <=  0.2)THEN
    
            data_out(:,:,:) = data_in(:,:,:,1)
    
        
        ELSEIF( (wspeed - wind_speed(3)) >=  -0.2)THEN
        
            data_out(:,:,:) = data_in(:,:,:,3) 
            
 
        ELSEIF(abs(wspeed - wind_speed(2)) <= 0.2)THEN
 
 
            data_out(:,:,:) = data_in(:,:,:,2) 
 
        ELSE
        
            IF(wspeed - wind_speed(2) < 0.0)THEN
                deltau = wind_speed(2) - wind_speed(1)
                iuhi  = 2
                iulow = 1
 
                data_out(:,:,:) = ( data_in(:,:,:,1)  * (wind_speed(iuhi)- wspeed) + &
                            data_in(:,:,:,2)  * (wspeed - wind_speed(iulow) ))/deltau
 
                
            ELSE
 
            deltau = wind_speed(3) - wind_speed(2)
            iuhi = 3
            iulow = 2
 
                data_out(:,:,:) = ( data_in(:,:,:,2)  * (wind_speed(iuhi)- wspeed) + &
                            data_in(:,:,:,3)  * (wspeed - wind_speed(iulow) ))/deltau
 
 
 
            ENDIF
 
       
        ENDIF
 
 
    end subroutine interpolate_wind_speed
    
    
 
    !--------------------------SUBROUTINE GetReflForGivenWindSpeedIce-------------------------------
    !This subroutine outputs the Reflectance table (R(tau,lambda, re)) for a given sun-satellite
    !geometry and a wind speed(m/s)
    !INPUTS :
    !   solarAng    :   solar zenith in degrees (scalar)
    !   viewAng     :   viewing angle in degrees (scalar)
    !   azmAng      :   relative azimuth between sat and solar planes
    !   wspeed      :   wind speed m/s
    !OUTPUTS :
    !   reflAsurf   :   reflectance as a function of tau,lambda,re : 3-D Array
    !                   DIM-1 = ntau, DIM-2 = nlambda, DIM-3 = nre
    !   ierror          :   integer value
    !                       ierror  = 0     success, 
    !                               = 1     I/O error
    !                               = 2     memory allocation error
    !                               = 3     invalid data    
    !ROUTINES CAlLED :
    !   ReadPhaseFunctionDataHDF
    !   ReadIceLibraryHDF   
    !   getAllPhaseFunctionValues
    !   TriLinearInterpolationReflSearch
    !   CalcSingScatPartOfReflectance
    !-----------------------------------------------------------------------------------------------
    SUBROUTINE getsdevreflforgivenwindspeed(solarAng,viewAng,azmAng,wspeed, &
     inrefl_ws_water, reflAsurf_water, inrefl_ws_ice, reflAsurf_ice, &
     wind_speed_only, ierror)
 
    use libraryarrays
    use ch_xfr, only : c2_scan
    !
    REAL, INTENT(IN) :: solarAng,viewAng,azmAng,wspeed
    real, intent(inout) :: inrefl_ws_water(:,:,:,:), inrefl_ws_ice(:,:,:,:)
    REAL,INTENT(OUT) :: reflAsurf_water(:,:,:), reflAsurf_ice(:,:,:)
    INTEGER,INTENT(OUT) :: ierror
    logical, intent(in) :: wind_speed_only
    
    INTEGER :: allocateStatus,iuhi,iulow
    REAL :: deltau
    integer :: num_radii
   INTEGER::do_mgr
    
    if (.not. wind_speed_only) then 
 
!
!     ! WDR wire up the sfc typ 1 for 3 m/s
       inrefl_ws_water(:,:,:,1) = 0
       inrefl_ws_ice(:,:,:,1) = 0
       call mng_ms_get_lib( solarang, viewang, azmang, 1, 1, c2_scan, &
         inrefl_ws_water(:,:,:,1), inrefl_ws_ice(:,:,:,1), ierror )
       if( ierror /= 0 ) return
 
!
!     ! WDR wire up the sfc typ 2 for 7 m/s
       inrefl_ws_water(:,:,:,2) = 0
       inrefl_ws_ice(:,:,:,2) = 0
       call mng_ms_get_lib( solarang, viewang, azmang, 2, 1, c2_scan, &
         inrefl_ws_water(:,:,:,2), inrefl_ws_ice(:,:,:,2), ierror )
       if( ierror /= 0 ) return
 
!
!     ! WDR wire up the sfc typ 3 for 15 m/s
       inrefl_ws_water(:,:,:,3) = 0
       inrefl_ws_ice(:,:,:,3) = 0
       call mng_ms_get_lib( solarang, viewang, azmang, 3, 1, c2_scan, &
         inrefl_ws_water(:,:,:,3), inrefl_ws_ice(:,:,:,3), ierror )
       if( ierror /= 0 ) return
 
     !deallocate( wtrsd_int, icesd_int )
    
    endif
    
    
    call interpolate_wind_speed(wspeed, inrefl_ws_water, reflasurf_water)
    call interpolate_wind_speed(wspeed, inrefl_ws_ice, reflasurf_ice)
 
 
    RETURN
 
    END SUBROUTINE getsdevreflforgivenwindspeed
    !-----------------------------------------------------------------------------------------------
   SUBROUTINE getsdevrefllamb(solarAng,viewAng,azmAng, reflAsurf_water, &
     reflAsurf_ice, ierror)
 
                                        
    use libraryarrays
   ! WDR access the scan line
   use ch_xfr, only : c2_scan
    !
    REAL, INTENT(IN) :: solarAng,viewAng,azmAng
    REAL,INTENT(OUT) :: reflAsurf_water(:,:,:), reflAsurf_ice(:,:,:)
    INTEGER,INTENT(OUT) :: ierror
    
    INTEGER :: allocateStatus,iuhi,iulow
    REAL :: deltau
    integer :: num_radii
   INTEGER::do_mgr
 
     reflasurf_water = 0
     reflasurf_ice = 0
     call mng_ms_get_lib( solarang, viewang, azmang, 0, 1, c2_scan, &
       reflasurf_water, reflasurf_ice, ierror )
     if( ierror /= 0 ) return
 
!   deallocate( wtrsd_int, icesd_int )
 
    END SUBROUTINE getsdevrefllamb
 
 end module interpolate_libraries