NASA Ocean Color

ocssw V2022

 module get_retrieval_uncertainty
  
  implicit none
  private
  
   logical :: print_info
   
   real :: new_water_radii(25), new_ice_radii(25)
  
  public ::  getuncertainties, getradiibounds, init_half_radii
  
  contains
  
 subroutine init_half_radii
  
     use libraryarrays
  
  
     integer :: i, n
     
     n = number_waterradii
     new_water_radii(1) = water_radii(1)
     do i=2, n
         new_water_radii(i) = water_radii(i) - (water_radii(i) - water_radii(i-1)) / 2.0 
     end do
     new_water_radii(n+1) = water_radii(n)
     
     n = number_iceradii
     new_ice_radii(1) = ice_radii(1)
     do i=2, n
         new_ice_radii(i) = ice_radii(i) - (ice_radii(i)-ice_radii(i-1)) / 2.0
     end do
     new_ice_radii(n+1) = ice_radii(n)
     
  
 end subroutine init_half_radii
  
  
  
  
 subroutine getuncertainties  (thickness,               &
                               re,                      &
                               water_path,              &
                               phase,                   &
                               R1R2wavelengthIdx,       &
                               meas_error, &
                               surface_albedo,          &
                               transmittance_w1,        &
                               transmittance_w2,        &
                               delta_transmittance_w1,  &
                               delta_transmittance_w2,  &
                               transmittance_stddev_w1, &
                               transmittance_stddev_w2, &
                               emission_pw, &
                               emission_Tc, &
                               sigma_R37_pw, &
                               uTau,uRe,uWaterPath, xpoint, ypoint)
 !-----------------------------------------------------------------------
 ! !f90
 !
 ! !Description:
 !     This subroutine computes the final uncertainties in thickness, re and water vapor
 !
 ! !Input Parameters:
 !     NONE
 !
 ! !Output Parameters:
 !     NONE
 !
 ! !Revision History:
 !     2004/06/01 bwind: Brad Wind
 !     Revision 1.0 Initial Revision
 !
 ! !Team-Unique Header:
 !    Developed by the  Cloud Retrieval Group, NASA GSFC, Greenbelt, Maryland, USA.
 !
 ! !References and Credits:
 !    Written by:
 !    Bradley P. Wind
 !    L3 GSI
 !    Code 913, NASA/GSFC
 !    Greenbelt, MD 20771
 !    bwind@climate.gsfc.nasa.gov
 ! 
 !    Mark Gray
 !    L3-GSI
 !    Climate and Radiation Branch, Code 913
 !    NASA/GSFC
 !    Greenbelt MD 20771
 !    gray@climate.gsfc.nasa.gov
 !
 ! !Design Notes:
 !   NONE
 !
 ! !END
 !
 !-----------------------------------------------------------------------
  
 ! 
 ! S. Platnick, 2-23-05: 
 !
 !  Wrote variance calculations in terms of matrix formulation
 !  Added covariance(tau,re) to water path uncertainty
 !
  
   use generalauxtype
   use science_parameters
   use mod06_run_settings
   use ch_xfr, only: c2_sensor_id, oci_id, ocis_id
 #if NOSWIR
   use core_arrays, only: mean_delta_ozone, ozone_transmittance, optical_thickness_37_final, &
     optical_thickness_1621_final
 #else
   use core_arrays, only: mean_delta_ozone, ozone_transmittance
 #endif
   use nonscience_parameters, only : fillvalue_real
  
   implicit none
  
     integer, intent(in) :: xpoint, ypoint
   real, intent(in) :: thickness, re, water_path,            &
                               surface_albedo(2),     &
                               transmittance_w1,transmittance_w2, &
                               delta_transmittance_w1,delta_transmittance_w2, &
                               transmittance_stddev_w1, transmittance_stddev_w2, meas_error(2), &
                               emission_pw(:), emission_tc(:), sigma_r37_pw(:)
  
    
   integer, intent(in)       :: r1r2wavelengthidx(2)
   character(10), intent(in) :: phase
   real, intent(out) :: utau, ure,uwaterpath
  
   real :: meancloudtopreflw1, meancloudtopreflw2
   real :: partialderivtauwrtr1atconstr2, partialderivtauwrtr2atconstr1 
   real :: partialderivrewrtr2atconstr1 , partialderivrewrtr1atconstr2
   real :: meandeltar1trans, meandeltar2trans
   real :: deltarcloudtopmean(2), density, sdev_refl(2)
  
   real :: correlation_trans, correlation_wp
   real :: emis_pw_val, emis_tc_val, sigma_r37_val
  
 ! declare matrices
   real, dimension(2,2) :: s, s_transpose, retrieval_covariance, dummy,  &
                                   refl_cov_total, refl_cov_sfc_albedo, &
                                   refl_cov_trans, refl_cov_trans_pw, refl_cov_trans_table, &
                                   refl_cov_calibration, refl_cov_cm_sdev, refl_cov_ws, &
                                   refl_cov_emission_pw, refl_cov_emission_tc, refl_cov_trans_o3
  
     integer :: band37_index
  
 #ifdef AIRBORNE
     band37_index = band_0370 + 3
 #else
  
 #ifdef AHI_INST 
     band37_index = band_0370
 #else
     band37_index = band_0370 - 1
 #endif
  
 #endif
     
     
     print_info = .false.
  
 ! -----------------------------------------------------------------------
 ! Calculate change in cloud top reflectance due to sfc albedo uncertainty
 ! -----------------------------------------------------------------------
     sdev_refl = 0.
  
     if (cox_munk) then
  
  
         call getctrefl_windspeeddiff(re, thickness, phase,   &
                             r1r2wavelengthidx,  &  
                             deltarcloudtopmean, meancloudtopreflw1, meancloudtopreflw2)
     else
     
         call getctref_albedodiff(re, thickness, phase,     &
                           r1r2wavelengthidx,  &
                           surface_albedo,     &
                           deltarcloudtopmean, meancloudtopreflw1, meancloudtopreflw2)
  
     endif
     
     call get_refl_windvector(re, thickness, phase, r1r2wavelengthidx, sdev_refl)
         
 ! -----------------------------------------------------------------------
 ! Calculate sensitivity derivatives and generate the sensitivity matrix S    
 ! -----------------------------------------------------------------------
  
  
     call sensitivitypartialderivatives( r1r2wavelengthidx, &
                                     re, thickness, phase,  &
                                     surface_albedo, &
                                     partialderivtauwrtr1atconstr2, &
                                     partialderivtauwrtr2atconstr1, &
                                     partialderivrewrtr1atconstr2,  &
                                     partialderivrewrtr2atconstr1 )
  
                                     
     s(1,1) = partialderivtauwrtr1atconstr2
     s(1,2) = partialderivtauwrtr2atconstr1  
     s(2,1) = partialderivrewrtr1atconstr2
     s(2,2) = partialderivrewrtr2atconstr1
  
     s_transpose = transpose(s)
     
  
 ! ----------------------------------------------------------------------
 ! Calculate the uncertainties attributable to each source of uncertainty 
 ! ----------------------------------------------------------------------
  
 ! (1) Reflectance covariance matrix due to Surface albedo uncertainty
  
     refl_cov_cm_sdev(1,1) = sdev_refl(1)**2
 #if NOSWIR
     refl_cov_cm_sdev(2,2) = 0.0
 #else
     refl_cov_cm_sdev(2,2) = sdev_refl(2)**2
 #endif
     refl_cov_cm_sdev(1,2) = 0.
     refl_cov_cm_sdev(2,1) = refl_cov_cm_sdev(1,2)
  
     if (cox_munk) then 
 ! wind speed uncertainty
         refl_cov_ws(1,1) = deltarcloudtopmean(1)**2
 #if NOSWIR
         refl_cov_ws(2,2) = 0.0
         
 #else
         refl_cov_ws(2,2) = deltarcloudtopmean(2)**2
 #endif
         refl_cov_ws(1,2) = 0.
         refl_cov_ws(2,1) = refl_cov_ws(1,2)
 ! wind vector uncertainty       
  
         refl_cov_sfc_albedo = refl_cov_ws + refl_cov_cm_sdev
  
     else
         refl_cov_sfc_albedo(1,1) = deltarcloudtopmean(1)**2
 #if NOSWIR
         refl_cov_sfc_albedo(2,2) = 0.0
 #else
         refl_cov_sfc_albedo(2,2) = deltarcloudtopmean(2)**2
 #endif
         refl_cov_sfc_albedo(1,2) = 0.
         refl_cov_sfc_albedo(2,1) = refl_cov_sfc_albedo(1,2)
  
         refl_cov_sfc_albedo = refl_cov_sfc_albedo + refl_cov_cm_sdev
     endif
  
 ! (2) Reflectance covariance matrix due to atmospheric transmittance uncertainty
  
 !  2a. Calculate covariance matrix due to percipitable water (PW) first:
  
 !  deltaR due to PW uncertainty only:
    meandeltar1trans =  meancloudtopreflw1 * abs(delta_transmittance_w1)/transmittance_w1
 #if NOSWIR
    meandeltar2trans =  0.0
 #else
    meandeltar2trans =  meancloudtopreflw2 * abs(delta_transmittance_w2)/transmittance_w2
 #endif
  
    correlation_trans = 1.0
    refl_cov_trans_pw(1,1) = meandeltar1trans**2
    refl_cov_trans_pw(2,2) = meandeltar2trans**2
  
    if( ( c2_sensor_id /= oci_id ) .AND. ( c2_sensor_id /= ocis_id ) )then
      if (r1r2wavelengthidx(2) == band37_index) then 
   
        ! the sigma_r37_val is already squared. 
   
          call get_emission_values(re,  sigma_r37_pw, phase, sigma_r37_val)
          refl_cov_trans_pw(1,2) = abs(meandeltar1trans * sqrt(sigma_r37_val + &
          refl_cov_trans_pw(2,2) ))*correlation_trans
          refl_cov_trans_pw(2,1) = refl_cov_trans_pw(1,2)
         
          refl_cov_trans_pw(2,2) = refl_cov_trans_pw(2,2) + sigma_r37_val
       else
         refl_cov_trans_pw(1,2) = abs(meandeltar1trans*meandeltar2trans)* &
          correlation_trans
         refl_cov_trans_pw(2,1) = refl_cov_trans_pw(1,2)
     
       endif
    endif
  
                 
 !  2b. Calculate covariance matrix due to table uncertainty next:
  
 !  deltaR due to table uncertainty only:
    meandeltar1trans =  meancloudtopreflw1 * transmittance_stddev_w1/transmittance_w1
 #if NOSWIR
    meandeltar2trans =  0.0
 #else
    meandeltar2trans =  meancloudtopreflw2 * transmittance_stddev_w2/transmittance_w2
 #endif
  
    refl_cov_trans_table(1,1) = meandeltar1trans**2
 #if NOSWIR
    refl_cov_trans_table(2,2) = 0.0
 #else
    refl_cov_trans_table(2,2) = meandeltar2trans**2
 #endif
    refl_cov_trans_table(1,2) = 0.
    refl_cov_trans_table(2,1) = refl_cov_trans_table(1,2)
  
 ! 2c. uncertainty due to ozone
  
     refl_cov_trans_o3 = 0.
     if (r1r2wavelengthidx(1) == band_0065 .and. mean_delta_ozone /= fillvalue_real &
             .and. ozone_transmittance /= fillvalue_real ) then 
         refl_cov_trans_o3(1,1) = (meancloudtopreflw1 * abs(mean_delta_ozone)/ozone_transmittance)**2
     endif
  
  
 !  2d. Combined covariance matrix:
  
    refl_cov_trans = refl_cov_trans_pw + refl_cov_trans_table + refl_cov_trans_o3
  
 ! (3) Reflectance covariance matrix due to calibration uncertainty
  
    refl_cov_calibration(1,1) = (meas_error(1) * meancloudtopreflw1)**2
 #if NOSWIR
    refl_cov_calibration(2,2) = 0.0
 #else
    refl_cov_calibration(2,2) = (meas_error(2) * meancloudtopreflw2)**2 
 #endif
 ! delta F0 = 0.42 taken from difference between Fontenla and Thekaekara (band averaged Aqua and Terra)
  
    if( ( c2_sensor_id /= oci_id ) .and. ( c2_sensor_id /= ocis_id ) )then
      if (r1r2wavelengthidx(2) == band37_index) then 
        refl_cov_calibration(2,2) = refl_cov_calibration(2,2) + &
          (( 0.42 * meancloudtopreflw2 ) / solar_constant_37)**2
      endif
    endif
    
    refl_cov_calibration(1,2) = 0.
    refl_cov_calibration(2,1) = refl_cov_calibration(1,2)
  
 ! Total tau & re covariance matrix
  
    refl_cov_total = refl_cov_sfc_albedo + refl_cov_trans + refl_cov_calibration
  
 ! ------------------------------------------------------------------------------
 ! Calculate retrieval covariance matrix = S * retrieval_covariance * S_transpose
 ! ------------------------------------------------------------------------------
  
    dummy = matmul(refl_cov_total, s_transpose)
    retrieval_covariance = matmul(s, dummy)
 #if NOSWIR
     ! Use the standard deviation of all possible CER values (calculated from the LUT CER grid),
     ! uniformly distributed (i.e., each is equally possible). It is important to have a somewhat
     ! reasonable CER uncertainty for calculating the water path uncertainty.
     if ( ( phase == 'water') .and. ( c2_sensor_id /= oci_id ) .and &
      ( c2_sensor_id /= ocis_id ) ) then
         retrieval_covariance(1,1) = retrieval_covariance(1,1) + &
             (optical_thickness_37_final(xpoint,ypoint))**2  ! Add uncertainty due to unknown CER
         retrieval_covariance(2,2) = 64.0    ! stddev ~8.0, assuming 14µm CER
     else
         retrieval_covariance(1,1) = retrieval_covariance(1,1) + &
             (optical_thickness_1621_final(xpoint,ypoint))**2    ! Add uncertainty due to unknown CER
         retrieval_covariance(2,2) = 289.0   ! stddev ~17.0, assuming 30µm CER
     end if
 #endif
  
  
 ! ---------------------
 ! Tau, re uncertainties
 ! ---------------------
  
    utau = sqrt( retrieval_covariance(1,1) ) 
    ure  = sqrt( retrieval_covariance(2,2) ) 
  
  
 ! ----------------------
 ! Water path uncertainty
 ! ----------------------
  
  
    if ((100.*utau/thickness) >= 200. .OR. (100.*ure/re) >= 200.  ) then
         uwaterpath = 200. 
  
  
    else
  
       if (phase == 'water') then
          density = liquid_water_density
       else
          density = ice_water_density
       endif
       
       uwaterpath = (re**2)*(utau**2) + (thickness**2)*(ure**2) + (utau**2)*(ure**2) + &
                     2*retrieval_covariance(1,2)*utau*ure + retrieval_covariance(1,2)**2
       if (uwaterpath < 0.) then
         uwaterpath = 200.
       else
          uwaterpath= sqrt( density**2 * uwaterpath )
       end if
    endif
  
  
 ! ----------------------
 ! Relative uncertainties
 ! ----------------------
  
    utau = 100.*utau/thickness
    ure  = 100.*ure/re
    uwaterpath = 100.*uwaterpath/water_path
  
 ! ------------------------------------------------------
 ! Limit answer to a maximum of 200% relative uncertainty
 ! ------------------------------------------------------
  
 !  sep: Was told by mag that max value should be 200 (or 255) 
 !       because of the scaling value used in the filespec
  
    if (utau > 200.  ) utau = 200.
    if (ure > 200. ) ure = 200.
    if (uwaterpath > 200.) uwaterpath = 200.
    
    
  
 end subroutine getuncertainties
  
 subroutine get_emission_values(re, sigma37,  phase,  sigma37_val)
  
     use libraryarrays, only : water_radii, ice_radii, number_waterradii, number_iceradii
     use modis_numerical_module, only : bisectionsearch, linearinterpolation
     
     implicit none
  
     character*(*), intent(in) :: phase
     real, intent(in) :: sigma37(:), re
     real, intent(inout) :: sigma37_val
     
     integer :: i, n, ilow, ihigh
  
    ilow = 0 ! WDR-UIV
    ihigh = 0 ! WDR-UIV 
  
    ilow = 0 ! WDR-UIV
    ihigh = 0 ! WDR-UIV
     
     if (phase == 'water') then 
         call bisectionsearch(water_radii, re, ilow, ihigh)
         sigma37_val = linearinterpolation( (/water_radii(ilow), water_radii(ihigh) /), &
                         (/ sigma37(ilow), sigma37(ihigh) /), re)
     else 
         call bisectionsearch(ice_radii, re, ilow, ihigh)    
         sigma37_val = linearinterpolation( (/ice_radii(ilow), ice_radii(ihigh) /), &
                         (/ sigma37(ilow), sigma37(ihigh) /), re)
     endif
  
  
 end subroutine get_emission_values
  
  
  
 subroutine getctrefl_windspeeddiff(&
                               radius,      &
                              optical_thickness, &
                              phase,             &
                              wave_index,        &
                              reflectancediff, meancloudtopreflw1, meancloudtopreflw2)
  
  
   use generalauxtype
   use libraryinterpolates
   use libraryarrays
   use science_parameters
   use mod06_run_settings
   use interpolate_libraries, only: interpolate_wind_speed
   implicit none
  
   real, intent(in)  :: radius
   integer, intent(in)       :: wave_index(2)
   real, intent(in)  :: optical_thickness
   character(10),intent(in)  :: phase
   real, intent(out) :: reflectancediff(2), meancloudtopreflw1, meancloudtopreflw2
  
   integer :: radii_size , radius_upperbound, radius_lowerbound, radius_index, i,j, vectorsize
   real :: reflectance_lower(2), reflectance_middle(2), reflectance_upper(2)
   real, dimension(:), allocatable :: opticalthick_vector  !(6)
   
   real, dimension(:,:,:), allocatable :: temp_refl
   integer :: status
   
   integer :: dummy
   
   integer :: TOTAL_POINTS
   
      integer :: start_time, end_time, cmax, crate
  
  
   
   total_points = size(library_taus) + 1
   
   
   allocate(opticalthick_vector(total_points))
   
   
   reflectance_upper = 0.
   reflectance_lower = 0.
   reflectance_middle = 0.
  
  
   if (phase == 'water') then
     call getclosestradius(number_waterradii,  &
                    water_radii,radius,  &
                    i)
  
   else
     call getclosestradius(number_iceradii,  &
                         ice_radii,radius,  &
                         i)
   endif
  
  
     opticalthick_vector(1) = 0.
     opticalthick_vector(2:total_points) = library_taus(1:(total_points-1))
  
  
     if (phase == 'water') then 
         
         
         do j =1, 2
             call nonasymptotic_calcu_cox_munk(opticalthick_vector, &
                                               optical_thickness, &
                                               int_reflectance_water(:,wave_index(j),i), &
                                               reflectance_middle(j))
         end do
     
         
         allocate(temp_refl(total_points, number_wavelengths, number_waterradii))
  
         call interpolate_wind_speed(lastinterp_wind_speed*(1.0+wind_speed_error), &
                                         int_reflectance_water_wspeed, &
                                         temp_refl(:,:,:) )
  
  
         
  
         do j =1, 2
             call nonasymptotic_calcu_cox_munk(opticalthick_vector, &
                                               optical_thickness, &
                                               temp_refl(:,wave_index(j),i), &
                                               reflectance_upper(j))
         end do
  
         call interpolate_wind_speed(lastinterp_wind_speed*(1.0-wind_speed_error), &
                                         int_reflectance_water_wspeed, &
                                         temp_refl(:,:,:) )
  
  
         do j =1, 2
             call nonasymptotic_calcu_cox_munk(opticalthick_vector, &
                                               optical_thickness, &
                                               temp_refl(:,wave_index(j),i), &
                                               reflectance_lower(j))
         end do
  
  
         deallocate(temp_refl)
  
         do j=1, 2
  
             reflectancediff(j) = (abs(reflectance_lower(j) - reflectance_middle(j)) +  &
                          abs(reflectance_upper(j) - reflectance_middle(j)))/2.     
             if(j .eq. 1) then
                 meancloudtopreflw1 = reflectance_middle(j)
             else
                 meancloudtopreflw2 = reflectance_middle(j)
             endif
  
  
         end do
  
     else
     
         do j =1, 2
             call nonasymptotic_calcu_cox_munk(opticalthick_vector, &
                                               optical_thickness, &
                                               int_reflectance_ice(:,wave_index(j),i), &
                                               reflectance_middle(j))
         end do
  
     
         allocate(temp_refl(total_points, number_wavelengths, number_iceradii))
  
         call interpolate_wind_speed(lastinterp_wind_speed*(1.0+wind_speed_error), &
                                         int_reflectance_ice_wspeed, &
                                         temp_refl(:,:,:) )
  
  
         do j =1, 2
             call nonasymptotic_calcu_cox_munk(opticalthick_vector, &
                                               optical_thickness, &
                                               temp_refl(:,wave_index(j),i), &
                                               reflectance_upper(j))
         end do
  
         call interpolate_wind_speed(lastinterp_wind_speed*(1.0-wind_speed_error), &
                                         int_reflectance_ice_wspeed, &
                                         temp_refl(:,:,:) )
  
  
         do j =1, 2
             call nonasymptotic_calcu_cox_munk(opticalthick_vector, &
                                               optical_thickness, &
                                               temp_refl(:,wave_index(j),i), &
                                               reflectance_lower(j))
         end do
  
  
         deallocate(temp_refl)
  
         do j=1, 2
  
             reflectancediff(j) = (abs(reflectance_lower(j) - reflectance_middle(j)) +  &
                          abs(reflectance_upper(j) - reflectance_middle(j)))/2.     
             if(j .eq. 1) then
                 meancloudtopreflw1 = reflectance_middle(j)
             else
                 meancloudtopreflw2 = reflectance_middle(j)
             endif
  
  
         end do
     
     
     endif
  
    
    deallocate(opticalthick_vector)
  
  
  
 end subroutine getctrefl_windspeeddiff
  
 subroutine get_refl_windvector(radius, optical_thickness, phase, wave_index, sdev_refl)
  
     
   use generalauxtype
   use libraryinterpolates
   use libraryarrays
   use science_parameters
   use mod06_run_settings
   implicit none
  
   real, intent(in)  :: radius
   integer, intent(in)       :: wave_index(2)
   real, intent(in)  :: optical_thickness
   character(10),intent(in)  :: phase
   real, intent(out) :: sdev_refl(2)
  
   integer :: radii_size , radius_upperbound, radius_lowerbound, radius_index, i,j, vectorsize
   real, dimension(:), allocatable :: opticalthick_vector  !(6)
   
   integer :: status
   
   integer :: dummy
   
   integer :: TOTAL_POINTS
   
   total_points = size(library_taus) + 1
   
   allocate(opticalthick_vector(total_points))
   
   
  
   if (phase == 'water') then
     call getclosestradius(number_waterradii,  &
                    water_radii,radius,  &
                    i)
  
   else
     call getclosestradius(number_iceradii,  &
                         ice_radii,radius,  &
                         i)
   endif
  
  
     opticalthick_vector(1) = 0.
     opticalthick_vector(2:total_points) = library_taus(1:(total_points-1))
  
  
     if (phase == 'water') then 
         
         do j =1, 2
             call nonasymptotic_calcu_cox_munk(opticalthick_vector, &
                                               optical_thickness, &
                                               int_reflectance_water_sdev(:,wave_index(j),i), &
                                               sdev_refl(j))
         end do
  
     else
     
         do j =1, 2
             call nonasymptotic_calcu_cox_munk(opticalthick_vector, &
                                               optical_thickness, &
                                               int_reflectance_ice_sdev(:,wave_index(j),i), &
                                               sdev_refl(j))
         end do
     
     endif
  
  
     deallocate(opticalthick_vector)
  
 end subroutine get_refl_windvector
  
  
  
 subroutine getctref_albedodiff(&
                               radius,      &
                              optical_thickness, &
                              phase,             &
                              wave_index,        &
                              surface_albedo,    &
                              reflectancediff, meancloudtopreflw1, meancloudtopreflw2)
   use generalauxtype
   use libraryinterpolates
   use libraryarrays
   use science_parameters
   use mod06_run_settings
   implicit none
  
   real, intent(in)  :: radius
   integer, intent(in)       :: wave_index(2)
   real, intent(in)  :: optical_thickness, surface_albedo(2)
   character(10),intent(in)  :: phase
   real, intent(out) :: reflectancediff(2), meancloudtopreflw1, meancloudtopreflw2
  
   integer :: radii_size , radius_upperbound, radius_lowerbound, radius_index, i,j, vectorsize
   real :: reflectance_lower, reflectance_middle, reflectance_upper
   real, dimension(:), allocatable :: opticalthick_vector  !(6)
   integer :: dummy
   
   integer :: TOTAL_POINTS
   
   total_points = size(library_taus) + 1
   
   
   allocate(opticalthick_vector(total_points))
   
   
   reflectance_upper = 0.
   reflectance_lower = 0.
   reflectance_middle = 0.
  
  
   if (phase == 'water') then
     call getclosestradius(number_waterradii,  &
                    water_radii,radius,  &
                    i)
  
   else
     call getclosestradius(number_iceradii,  &
                         ice_radii,radius,  &
                         i)
                         
   endif
  
  
     opticalthick_vector(1) = 0.
     opticalthick_vector(2:total_points) = library_taus(1:(total_points-1))
  
   if (phase == 'water') then
    do j = 1,2
  
  
       call nonasymptotic_albedovar(opticalthick_vector, &
                               optical_thickness,                         &
                               spherical_albedo_water(:,wave_index(j),i),&
                               int_reflectance_water(:,wave_index(j),i), &
                               int_fluxdownwater_sensor(:,wave_index(j),i), &
                               int_fluxdownwater_solar(:,wave_index(j),i), & 
                               surface_albedo(j),                              &
                               reflectance_lower, &
                               reflectance_middle, &
                               reflectance_upper)
  
  
         reflectancediff(j) = (abs(reflectance_lower - reflectance_middle) +  &
                          abs(reflectance_upper - reflectance_middle))/2.     
         if(j .eq. 1) then
             meancloudtopreflw1 = reflectance_middle
         else
             meancloudtopreflw2 = reflectance_middle
         endif
     enddo 
   else 
    do j = 1,2
  
              call nonasymptotic_albedovar(opticalthick_vector, &
                               optical_thickness,                         &
                               spherical_albedo_ice(:,wave_index(j),i),&
                               int_reflectance_ice(:,wave_index(j),i), &
                               int_fluxdownice_sensor(:,wave_index(j),i), &
                               int_fluxdownice_solar(:,wave_index(j),i), &
                               surface_albedo(j),                              &
                               reflectance_lower, &
                               reflectance_middle, &
                               reflectance_upper)
  
             reflectancediff(j) = (abs(reflectance_lower - reflectance_middle) + &
                            abs(reflectance_upper - reflectance_middle))/2. 
      if(j .eq. 1) then
         meancloudtopreflw1 = reflectance_middle
      else
         meancloudtopreflw2 = reflectance_middle
      endif
     enddo
    endif
    
    
    deallocate(opticalthick_vector)
  
 end subroutine getctref_albedodiff
  
  
 subroutine getclosestradius(numberradii, radiidata, radius, index)
  
   integer, intent(in):: numberradii
   real, intent(in)   :: radiidata(numberradii), radius
   integer, intent(out):: index
   integer :: index1(1)
  
   index1 = minloc(abs(radiidata-radius))
   index=index1(1)
   
 end subroutine getclosestradius
  
  
  
 subroutine nonasymptotic_calcu_cox_munk(tau_vector,tc, rf, rf_calculated)
   use generalauxtype
   use modis_numerical_module
   use mod06_run_settings
  
   implicit none
  
   real, intent(in)    :: tau_vector(:),tc
   real, intent(in)    :: rf(:)
   real, intent(out)   :: rf_calculated
  
  
    integer  :: lowbound, highbound, taubounds(2)
    real     :: iftau(2), refvals(2)
  
    lowbound = 0 ! WDR-UIV
    highbound = 0 ! WDR-UIV
  
       call bisectionsearch(tau_vector, tc, lowbound, highbound)
       taubounds(1) = lowbound
       taubounds(2) = highbound
       iftau(1)    = tau_vector(taubounds(1))
       iftau(2)    = tau_vector(taubounds(2))
  
       refvals(1) = rf(lowbound)
       refvals(2) = rf(highbound)
  
       rf_calculated = linearinterpolation(iftau,refvals,tc)
       
 end subroutine nonasymptotic_calcu_cox_munk
  
  
  
 subroutine nonasymptotic_calcu(tau_vector,tc, &
                             sfr, rf, ft1,ft0, &
                             albedo,rf_calculated)
   use generalauxtype
   use modis_numerical_module
   use mod06_run_settings
   use science_parameters
  
   implicit none
  
   real, intent(in)    :: tau_vector(:),tc,albedo
   real, intent(in)    :: rf(:), ft0(:), ft1(:), sfr(:)
   real, intent(out)   :: rf_calculated
  
   real, dimension(:), allocatable :: reflectance_vector
  
    integer  :: lowbound, highbound, taubounds(2)
    real     :: iftau(2), refvals(2)
  
   integer :: TOTAL_POINTS
  
   lowbound = 0 ! WDR-UIV
   highbound = 0 ! WDR-UIV
   
   total_points = size(tau_vector)
  
   allocate(reflectance_vector(total_points))
  
     reflectance_vector(1) =  albedo
  
     reflectance_vector(2:total_points) = rf(1:(total_points-1)) + &
         ft1(1:(total_points-1)) * ft0(1:(total_points-1)) * albedo /( 1. - albedo*sfr(1:(total_points-1)))
  
  
       call bisectionsearch(tau_vector, tc, lowbound, highbound)
       taubounds(1) = lowbound
       taubounds(2) = highbound
       iftau(1)    = tau_vector(taubounds(1))
       iftau(2)    = tau_vector(taubounds(2))
  
       refvals(1) = reflectance_vector(lowbound)
       refvals(2) = reflectance_vector(highbound)
  
       rf_calculated = linearinterpolation(iftau,refvals,tc)
  
     deallocate(reflectance_vector)
  
 end subroutine nonasymptotic_calcu
  
 subroutine nonasymptotic_albedovar(tau_vector,tc, &
                             sfr, rf, ft1,ft0, &
                             albedo,rf_lower, rf_middle, rf_upper)
   use generalauxtype
   use science_parameters
   implicit none
  
   real, intent(in)    :: tau_vector(:),tc,albedo
   real, intent(in)    :: rf(:), ft0(:), ft1(:), sfr(:)
   real, intent(out)   :: rf_middle, rf_lower, rf_upper
   real :: albedoUpper
  
  
   call nonasymptotic_calcu(tau_vector,tc, &
                             sfr, rf, ft1,ft0, &
                             albedo,rf_middle)
  
   call nonasymptotic_calcu(tau_vector,tc, &
                             sfr, rf, ft1,ft0, &
                             (albedo * (1.0 - albedo_error)),rf_lower)
  
   albedoupper = albedo * (1.0 + albedo_error)
   if(albedoupper .gt. 0.99) albedoupper = 0.99
   call nonasymptotic_calcu(tau_vector,tc, &
                             sfr, rf, ft1,ft0, &
                             albedoupper,rf_upper)
  
 end subroutine nonasymptotic_albedovar
  
  
   subroutine sensitivitypartialderivatives( &
                                            R1R2wavelengthIdx,&
                                            re, tau,  &
                                            phase, surface_albedo, &
                                            partialDerivTauWrtR1AtConstR2,&
                                            partialDerivTauWrtR2AtConstR1,&
                                            partialDerivReWrtR1AtConstR2, &
                                            partialDerivReWrtR2AtConstR1)
 !-----------------------------------------------------------------------
 ! !f90
 !
 ! !Description:
 !   This subroutine computes the partial derivatives of tau and re with 
 !   respect to reflectance in one band as the other one is held constant
 !
 !
 ! !Input Parameters:
 !
 ! !Output Parameters:
 !   NONE
 !
 ! !Revision History:
 !   2004/05/28 bwind: Brad Wind
 !   Revision 1.0 Initial Revision
 !
 ! !Team-Unique Header:
 !   Developed by the  Cloud Retrieval Group, NASA GSFC, Greenbelt, Maryland, USA.
 !
 ! !References and Credits:
 !   Written by:
 !   Bradley P. Wind
 !   L3 GSI
 !   Code 913, NASA/GSFC
 !   Greenbelt, MD 20771
 !   bwind@climate.gsfc.nasa.gov
 ! 
 !   Mark Gray
 !   L3-GSI
 !   Climate and Radiation Branch, Code 913
 !   NASA/GSFC
 !   Greenbelt MD 20771
 !   gray@climate.gsfc.nasa.gov
 !
 ! !Design Notes:
 !   NONE
 !
 ! !END
 !
 !-----------------------------------------------------------------------
   use generalauxtype
   use science_parameters
   use libraryarrays
   implicit none
  
   integer, intent(in)      :: R1R2wavelengthIdx(2)
   real, intent(in) :: re, tau
   real, intent(in) :: surface_albedo(2)
   character(10), intent(in):: phase
   real, intent(out) :: partialDerivTauWrtR1AtConstR2, &
                                partialDerivTauWrtR2AtConstR1, &
                                partialDerivReWrtR1AtConstR2,  &
                                partialDerivReWrtR2AtConstR1
  
   real :: partialDerivR1wrtTauAtConstRe, &
                   partialDerivR2wrtTauAtConstRe, &
                   partialDerivR1wrtReAtConstTau, &
                   partialDerivR2wrtReAtConstTau
  
   integer :: i
   integer :: wrtParam
  
   real :: partialDerivativesTau, partialDerivativesRe
   real :: upperLimit, lowerLimit,abscissaintervalstartingpoint
   real :: halfAbscissaInterval
   real, parameter :: halfAbscissaIntervalTau = 1., halfabscissaintervalre = 1.
   real, parameter :: tauUpperLimit = 200.
   real, parameter :: tauLowerLimit = 0.
  
   integer, parameter :: tauParamPosition = 4
   integer, parameter :: reParamPosition = 5
   real :: small_number
  
 ! for tau, determine abscissa interval taking into account the 
 ! extrema and potentially rough transition point(s)
  
   halfabscissainterval = halfabscissaintervaltau
   abscissaintervalstartingpoint = tau
   wrtparam = tauparamposition
  
  
 #if NOSWIR
  
     upperlimit=tauupperlimit
     lowerlimit=taulowerlimit
  
     i = 1
     call reflpartialderivative( &
         upperlimit, lowerlimit, &
         halfabscissainterval, &
         abscissaintervalstartingpoint, &
         re, tau,phase,  &
         surface_albedo(i), &
         wrtparam,                                &
         r1r2wavelengthidx(i), &
         partialderivativestau)
  
     small_number = epsilon(partialderivativestau)
  
     if(abs(partialderivativestau) .lt. small_number) partialderivativestau = &
         small_number
  
     partialderivtauwrtr1atconstr2 = 1. / partialderivativestau
     partialderivtauwrtr2atconstr1 = 0.0
     partialderivrewrtr1atconstr2 = 0.0
     partialderivrewrtr2atconstr1 = 0.0
  
 #else
  
 ! loop over wavelengths
   do i =1, 2
  
       upperlimit=tauupperlimit
       lowerlimit=taulowerlimit
  
 !   Compute derivative for dR1/dT )re, and dR2/dT )re 
     call reflpartialderivative( &
                                upperlimit, lowerlimit, &
                                halfabscissainterval, & 
                                abscissaintervalstartingpoint, &
                                re, tau,phase,  &
                                surface_albedo(i), &
                                wrtparam,                                &
                                r1r2wavelengthidx(i), &
                                partialderivativestau)
     if (i .eq. 1) then
 !     dR1/dT | re
       partialderivr1wrttauatconstre = partialderivativestau
     else
 !     dR2/dT | re
       partialderivr2wrttauatconstre = partialderivativestau
     endif
   end do
  
 ! for re, determine abscissa interval taking into account the extrema
   halfabscissainterval = halfabscissaintervalre
   abscissaintervalstartingpoint = re
  
  
  
 ! for re must determine which phase 
   if(phase == 'water') then
     upperlimit = water_radii(number_waterradii)
     lowerlimit = water_radii(1)
   else
     upperlimit = ice_radii(number_iceradii)
     lowerlimit = ice_radii(1)
   endif
  
 ! Compute derivative for dR1/dre )T, and dR2/dre)T
   wrtparam = reparamposition
  
   do i = 1, 2
    call reflpartialderivative(  &
                               upperlimit, lowerlimit,  &
                               halfabscissainterval,    &
                               abscissaintervalstartingpoint, &
                               re, tau,phase,           &
                               surface_albedo(i),       &
                               wrtparam,                &
                               r1r2wavelengthidx(i),    &
                               partialderivativesre)
  
     if (i == 1 ) then
 !     dR1 / dre | T
       partialderivr1wrtreatconsttau = partialderivativesre
     else
 !     dR2 / dre | T
       partialderivr2wrtreatconsttau = partialderivativesre
     endif
  
   enddo
  
  
 ! head-off divide by 0
     small_number = epsilon(partialderivr1wrttauatconstre)
  
   if(abs(partialderivr1wrttauatconstre) .lt. small_number) partialderivr1wrttauatconstre = &
      small_number
   if(abs(partialderivr2wrttauatconstre) .lt. small_number) partialderivr2wrttauatconstre = &
       small_number
   if(abs(partialderivr1wrtreatconsttau) .lt.  small_number) partialderivr1wrtreatconsttau = &
       small_number
   if(abs(partialderivr2wrtreatconsttau) .lt.  small_number) partialderivr2wrtreatconsttau = &
       small_number
  
 ! Compute the derivatives
   partialderivtauwrtr1atconstr2 = 1. / &
                                   (partialderivr1wrttauatconstre - &
                                    partialderivr1wrtreatconsttau * &
                                    (partialderivr2wrttauatconstre/partialderivr2wrtreatconsttau) )
     
   partialderivtauwrtr2atconstr1 = 1. / &
                                    (partialderivr2wrttauatconstre - &
                                     partialderivr2wrtreatconsttau * &
                                     (partialderivr1wrttauatconstre/partialderivr1wrtreatconsttau) )
    
   partialderivrewrtr1atconstr2 = 1. / &
                                   (partialderivr1wrtreatconsttau - &
                                    partialderivr1wrttauatconstre * &
                                    (partialderivr2wrtreatconsttau/partialderivr2wrttauatconstre) )
     
   partialderivrewrtr2atconstr1 = 1. / &
                                   (partialderivr2wrtreatconsttau - &
                                    partialderivr2wrttauatconstre * &
                                    (partialderivr1wrtreatconsttau/partialderivr1wrttauatconstre) )
  
 #endif
  
 if(debugprn) then
 print*, 'partialDerivTauWrtR1AtConstR2 = ', partialderivtauwrtr1atconstr2
 print*, 'partialDerivTauWrtR2AtConstR1 = ', partialderivtauwrtr2atconstr1
 print*, 'partialDerivReWrtR1AtConstR2 = ', partialderivrewrtr1atconstr2
 print*, 'partialDerivReWrtR2AtConstR1 = ', partialderivrewrtr2atconstr1
 endif
  
 end subroutine sensitivitypartialderivatives
   subroutine reflpartialderivative(upperLimit, lowerLimit, &
                                    halfAbscissaInterval, abscissaIntervalStartingPoint, &
                                    re, tau, phase, &
                                    surface_albedo, &
                                    wrtParam, wavelengthIdx, &
                                    partialDerivative)
 !-----------------------------------------------------------------------
 ! !f90
 !
 ! !Description:
 !     This subroutine computes a partial derivative of tau or re using the
 !     central difference method. Special cases such as crossing to asymptotic regime
 !     and being at libraries' ends are handled also. 
 !
 ! !Input Parameters:
 !     wrtParam -- integer -- flag controlling whether or not to compute the derivatives
 !
 ! !Output Parameters:
 !     partialDerivatives -- real*4, dimension(:) -- array of partial derivatives
 !
 ! !Revision History:
 !     2004/05/28 bwind: Brad Wind
 !     Revision 1.0 Initial Revision
 !
 ! !Team-Unique Header:
 !    Developed by the  Cloud Retrieval Group, NASA GSFC, Greenbelt, Maryland, USA.
 !
 ! !References and Credits:
 !    Written by:
 !    Bradley P. Wind
 !    L3 GSI
 !    Code 913, NASA/GSFC
 !    Greenbelt, MD 20771
 !    bwind@climate.gsfc.nasa.gov
 !
 !    Mark Gray
 !    L3-GSI
 !    Climate and Radiation Branch, Code 913
 !    NASA/GSFC
 !    Greenbelt MD 20771
 !    gray@climate.gsfc.nasa.gov
 !
 ! !Design Notes:
 !      NONE
 !
 ! !END
 !
 !-----------------------------------------------------------------------
   use generalauxtype
   use science_parameters
   implicit none
  
   integer, intent(in)       :: wrtParam, wavelengthIdx
   real, intent(in)  :: re, tau
   real, intent(in)  :: surface_albedo
   real, intent(in)  :: upperLimit, lowerLimit
   character(10),intent(in)  :: phase
   real, intent(out) :: partialDerivative
  
   real :: Rb, Ra
   real :: tmpTau, tmpRe
   
   real :: halfAbscissaInterval, abscissaIntervalStartingPoint
   real :: abscissaIntervalLower, abscissaIntervalUpper
   integer, parameter :: tauParamPosition = 4
   integer, parameter :: reParamPosition = 5
  
   integer :: idxJ
   call abscissainterval(upperlimit, lowerlimit, &
                         halfabscissainterval, abscissaintervalstartingpoint, &
                         abscissaintervallower, abscissaintervalupper)
     
  
  
  
   rb = 0.
   ra = 0.
  
 ! get reflectances
   select case (wrtparam)
     case(tauparamposition)
 !     in tau case
 !     ...abscissaIntervalUpper/Lower are upper and lower optical thicknesses
       call getctref_constradius(      &
                                 abscissaintervalupper,      &
                                 abscissaintervallower,      &
                                 re, phase,                  &
                                 surface_albedo,             &
                                 wavelengthidx, rb, ra)
         partialderivative = (rb - ra) 
         partialderivative = partialderivative / (abscissaintervalupper - abscissaintervallower)
  
     case(reparamposition)
 !     in re case
 !     ...abscissaIntervalUpper/Lower are upper and lower effective radii
  
         call getctref_consttau( &
                                abscissaintervalupper,      &
                                abscissaintervallower,      &
                                re, tau, phase,                  &
                                surface_albedo,             &
                                wavelengthidx, partialderivative)      
  
  
     case default
   end select
  
  
   ! perform the division
  
   end subroutine reflpartialderivative
  
  
 subroutine abscissainterval(upperLimit, lowerLimit, &
                             halfAbscissaInterval, abscissaIntervalStartingPoint, &
                             abscissaIntervalLower, abscissaIntervalUpper)
  
 !-----------------------------------------------------------------------
 ! !f90
 !
 ! !Description:
 !   This subroutine figures out the interval for use in the
 !   central difference method. Special cases such as crossing to asymptotic regime
 !   and being at libraries' ends are handled also. 
 !
 ! !Input Parameters:
 !   NONE
 !
 ! !Output Parameters:
 !   NONE
 !
 ! !Revision History:
 !   2004/05/28 bwind: Brad Wind
 !   Revision 1.0 Initial Revision
 !
 ! !Team-Unique Header:
 !   Developed by the  Cloud Retrieval Group, NASA GSFC, Greenbelt, Maryland, USA.
 !
 ! !References and Credits:
 !   Written by:
 !   Bradley P. Wind
 !   L3 GSI
 !   Code 913, NASA/GSFC
 !   Greenbelt, MD 20771
 !   bwind@climate.gsfc.nasa.gov
 !
 !   Mark Gray
 !   L3-GSI
 !   Climate and Radiation Branch, Code 913
 !   NASA/GSFC
 !   Greenbelt MD 20771
 !   gray@climate.gsfc.nasa.gov
 !
 ! !Design Notes:
 !   NONE
 !
 ! !END
 !
 !-----------------------------------------------------------------------\
   use science_parameters
  
   implicit none
  
   real, intent(in) :: upperLimit, lowerLimit, &
                               halfAbscissaInterval, abscissaIntervalStartingPoint
   real, intent(out) :: abscissaIntervalLower, abscissaIntervalUpper
  
   ! Reminder: halfAbscissaIntervalTau = 1., halfAbscissaIntervalRe = 1.
    
   if((abscissaintervalstartingpoint+halfabscissainterval) > upperlimit) then 
       abscissaintervalupper = abscissaintervalstartingpoint
   else
      abscissaintervalupper = (abscissaintervalstartingpoint+halfabscissainterval)
   end if
    
   if((abscissaintervalstartingpoint-halfabscissainterval) < lowerlimit) then 
      abscissaintervallower = abscissaintervalstartingpoint 
   else
      abscissaintervallower = (abscissaintervalstartingpoint-halfabscissainterval)
   end if
  
 end subroutine abscissainterval
 subroutine getctref_constradius( &
                                 optical_thickness_upper, &
                                 optical_thickness_lower, &
                                 radius,                  &
                                 phase,                   &
                                 surface_albedo,          &
                                 wave_index,              &
                                 reflectance_upper,       &
                                 reflectance_lower)
  
   use libraryinterpolates
   use libraryarrays
   use modis_numerical_module
   use mod06_run_settings
   use science_parameters
   implicit none
  
   real, intent(in)  :: optical_thickness_upper, optical_thickness_lower
   integer, intent(in)       :: wave_index
   real, intent(in)  :: radius, surface_albedo
   character(10),intent(in)  :: phase
   real, intent(out) :: reflectance_upper, reflectance_lower
  
   integer :: radii_size , radius_upperbound, radius_lowerbound, radius_index, i,j, points
   real, dimension(:), allocatable :: opticalthick_vector!(6)
   real :: reflectance_vector_upper(3), reflectance_vector_lower(3), yy2(3)
  
     integer :: TOTAL_POINTS
     
     total_points = size(library_taus) + 1
  
  
     allocate(opticalthick_vector(total_points))
  
  
      opticalthick_vector(1) = 0.
       opticalthick_vector(2:total_points) = library_taus(1:(total_points-1))
  
  
     if (phase == 'water') then
        call getclosestradius(number_waterradii,  &
                              water_radii,radius,  &
                              i)
  
         if (cox_munk) then 
  
                 call nonasymptotic_calcu_cox_munk(opticalthick_vector, &
                                               optical_thickness_lower, &
                                               int_reflectance_water(:,wave_index,i), &
                                               reflectance_lower)
  
                 call nonasymptotic_calcu_cox_munk(opticalthick_vector, &
                                               optical_thickness_upper, &
                                               int_reflectance_water(:,wave_index,i), &
                                               reflectance_upper)
         
         
         else 
             call nonasymptotic_calcu(opticalthick_vector, optical_thickness_lower,  &
                                   spherical_albedo_water(:,wave_index,i), int_reflectance_water(:,wave_index,i), &
                                   int_fluxdownwater_sensor(:,wave_index,i), int_fluxdownwater_solar(:,wave_index,i), & 
                                   surface_albedo, reflectance_lower)
             call nonasymptotic_calcu(opticalthick_vector, optical_thickness_upper, &
                                   spherical_albedo_water(:,wave_index,i), int_reflectance_water(:,wave_index,i), &
                                   int_fluxdownwater_sensor(:,wave_index,i), int_fluxdownwater_solar(:,wave_index,i), & 
                                   surface_albedo, reflectance_upper)
         endif
  
      else
        call getclosestradius(number_iceradii,  &
                              ice_radii,radius,  &
                              i)
                              
          if (cox_munk) then 
  
                 call nonasymptotic_calcu_cox_munk(opticalthick_vector, &
                                               optical_thickness_lower, &
                                               int_reflectance_ice(:,wave_index,i), &
                                               reflectance_lower)
  
                 call nonasymptotic_calcu_cox_munk(opticalthick_vector, &
                                               optical_thickness_upper, &
                                               int_reflectance_ice(:,wave_index,i), &
                                               reflectance_upper)
                                          
          else 
             call nonasymptotic_calcu(opticalthick_vector, optical_thickness_lower,  &
                                   spherical_albedo_ice(:,wave_index,i), int_reflectance_ice(:,wave_index,i), &
                                   int_fluxdownice_sensor(:,wave_index,i), int_fluxdownice_solar(:,wave_index,i), & 
                                   surface_albedo, reflectance_lower)
             call nonasymptotic_calcu(opticalthick_vector, optical_thickness_upper,  &
                                   spherical_albedo_ice(:,wave_index,i), int_reflectance_ice(:,wave_index,i), &
                                   int_fluxdownice_sensor(:,wave_index,i), int_fluxdownice_solar(:,wave_index,i), & 
                                     surface_albedo, reflectance_upper)
          endif
      endif
  
   
   deallocate(opticalthick_vector)
   
 end subroutine getctref_constradius
 subroutine getctref_consttau(&
                              radius_upper,      &
                              radius_lower,      &
                              radius, optical_thickness, &
                              phase,             &
                              surface_albedo,    &
                              wave_index,        &
                              partial_derivative)
   use generalauxtype
   use libraryinterpolates
   use libraryarrays
   use modis_numerical_module
   use mod06_run_settings
   use science_parameters
   implicit none
  
   real, intent(inout)  :: radius_upper, radius_lower
   integer, intent(in)       :: wave_index
   real, intent(in)  :: optical_thickness, radius, surface_albedo
   character(10),intent(in)  :: phase
   real, intent(out) :: partial_derivative
   
   real :: reflectance_upper, reflectance_lower
  
   integer :: radii_size , radius_upperbound, radius_lowerbound, radius_index, i,j, vectorsize
   real, dimension(:), allocatable :: opticalthick_vector, opticalthick_vector2!(6)
   real,dimension(20):: reflectance_vector
   
   integer :: dummy, il, jl, iu, ju, k
   real :: my_radius, myrefl(90), myres(90), myx(3), myy(3), half_rad(2)
   integer :: TOTAL_POINTS
   real :: mymid1, mymid2, reflectance_middle, delre
   real :: new_derivatives(25)
   integer :: istart, jstart
   
   total_points = size(library_taus) + 1
   
   
   allocate(opticalthick_vector(total_points), opticalthick_vector2(total_points))
   
  
      opticalthick_vector(1) = 0.
       opticalthick_vector(2:total_points) = library_taus(1:(total_points-1))
      opticalthick_vector2 = opticalthick_vector
      if (phase == 'water') then
  
         call getradiibounds(number_waterradii, water_radii, radius, j, i)
  
 !       radius_lower = water_radii(i)
 !       radius_upper = water_radii(j)
 !       call nonasymptotic_calcu(opticalthick_vector, optical_thickness,  &
 !                                spherical_albedo_water(:,wave_index,i), int_reflectance_water(:,wave_index,i), &
 !                                int_fluxdownwater_sensor(:,wave_index,i), int_fluxdownwater_solar(:,wave_index,i), & 
 !                                surface_albedo, reflectance_lower)
  
 !       call nonasymptotic_calcu(opticalthick_vector2, optical_thickness,  &
 !                                spherical_albedo_water(:,wave_index,j), int_reflectance_water(:,wave_index,j), &
 !                                int_fluxdownwater_sensor(:,wave_index,j), int_fluxdownwater_solar(:,wave_index,j), & 
 !                                surface_albedo, reflectance_upper)
  
         new_derivatives = 0.
         
         i = i-1
         j = j+1
         if (i<1) i = 1
         if (j>number_waterradii) j = number_waterradii
         
         if (cox_munk) then 
         
             do il = i, j-1
                 call nonasymptotic_calcu_cox_munk(opticalthick_vector, &
                                               optical_thickness, &
                                               int_reflectance_water(:,wave_index,il), &
                                               reflectance_lower)
  
                 call nonasymptotic_calcu_cox_munk(opticalthick_vector, &
                                               optical_thickness, &
                                               int_reflectance_water(:,wave_index,il+1), &
                                               reflectance_upper)
  
                 new_derivatives(il+1) = (reflectance_upper - reflectance_lower) / &
                                             (water_radii(il+1) - water_radii(il))
             end do
         
         else 
         
             do il = i, j-1
                 call nonasymptotic_calcu(opticalthick_vector, optical_thickness,  &
                                 spherical_albedo_water(:,wave_index,il), int_reflectance_water(:,wave_index,il), &
                                 int_fluxdownwater_sensor(:,wave_index,il), int_fluxdownwater_solar(:,wave_index,il), & 
                                 surface_albedo, reflectance_lower)
  
                 call nonasymptotic_calcu(opticalthick_vector2, optical_thickness,  &
                                 spherical_albedo_water(:,wave_index,il+1), int_reflectance_water(:,wave_index,il+1), &
                                 int_fluxdownwater_sensor(:,wave_index,il+1), int_fluxdownwater_solar(:,wave_index,il+1), & 
                                 surface_albedo, reflectance_upper)
                                 
                 new_derivatives(il+1) = (reflectance_upper - reflectance_lower) / &
                                             (water_radii(il+1) - water_radii(il))
                                                                             
             end do 
         endif
  
  
         new_derivatives(1) = new_derivatives(2)
         new_derivatives(number_waterradii+1) = new_derivatives(number_waterradii)
  
         call getradiibounds(number_waterradii+1, new_water_radii, radius, iu, il)
             
         partial_derivative = linearinterpolation( (/ new_water_radii(il), new_water_radii(iu) /), &
                                         (/ new_derivatives(il), new_derivatives(iu) /), &
                                            radius)
  
  
      else
      
         call getradiibounds(number_iceradii, ice_radii, radius, j, i)
  
  
         new_derivatives = 0.
         
         i  = i-1
         j  = j+1
         if (i <1) i = 1
         if (j>number_iceradii) j = number_iceradii
         
         if (cox_munk) then 
             do il = i, j-1
                 call nonasymptotic_calcu_cox_munk(opticalthick_vector, &
                                               optical_thickness, &
                                               int_reflectance_ice(:,wave_index,il), &
                                               reflectance_lower)
  
                 call nonasymptotic_calcu_cox_munk(opticalthick_vector, &
                                               optical_thickness, &
                                               int_reflectance_ice(:,wave_index,il+1), &
                                               reflectance_upper)
                                 
                 new_derivatives(il+1) = (reflectance_upper - reflectance_lower) / &
                                             (ice_radii(il+1) - ice_radii(il))
                                                                             
             end do 
         
         else 
             do il = i, j-1
                 call nonasymptotic_calcu(opticalthick_vector, optical_thickness,  &
                                 spherical_albedo_ice(:,wave_index,il), int_reflectance_ice(:,wave_index,il), &
                                 int_fluxdownice_sensor(:,wave_index,il), int_fluxdownice_solar(:,wave_index,il), & 
                                 surface_albedo, reflectance_lower)
  
                 call nonasymptotic_calcu(opticalthick_vector, optical_thickness ,  &
                                 spherical_albedo_ice(:,wave_index,il+1), int_reflectance_ice(:,wave_index,il+1), &
                                 int_fluxdownice_sensor(:,wave_index,il+1), int_fluxdownice_solar(:,wave_index,il+1), & 
                                 surface_albedo, reflectance_upper)
                                 
                 new_derivatives(il+1) = (reflectance_upper - reflectance_lower) / &
                                             (ice_radii(il+1) - ice_radii(il))
                                                                             
             end do 
         endif
  
         new_derivatives(1) = new_derivatives(2)
         new_derivatives(number_iceradii+1) = new_derivatives(number_iceradii)
  
  
         call getradiibounds(number_iceradii+1, new_ice_radii, radius, iu, il)
                     
         partial_derivative = linearinterpolation( (/ new_ice_radii(il), new_ice_radii(iu) /), &
                                         (/ new_derivatives(il), new_derivatives(iu) /), &
                                            radius)
  
      endif
  
  
   deallocate(opticalthick_vector, opticalthick_vector2)
  
  
 end subroutine getctref_consttau
  
 subroutine getradiibounds(radiisize, radiidata, radius, upperbound, lowerbound)
   use generalauxtype
   implicit none
  
   integer,      intent(in)  :: radiisize
   real, intent(in)  :: radiidata(:), radius
   integer,      intent(out) :: upperbound, lowerbound
  
   integer :: i
  
   upperbound = -1
   lowerbound = -1
   do i = 1, radiisize - 1
     if (radius >= radiidata(i).and. radius <= radiidata(i+1)) then
       upperbound = i+1
       lowerbound = i
      exit
     endif
   enddo
  
  
 end subroutine getradiibounds
  
 end module get_retrieval_uncertainty 

ch_xfr

Definition: ch_xfr.f90:1

libraryarrays::number_wavelengths

integer number_wavelengths

Definition: libraryarrays.f90:5

core_arrays

Definition: core_arrays.f90:1

core_arrays::ozone_transmittance

real ozone_transmittance

Definition: core_arrays.f90:129

get_retrieval_uncertainty::nonasymptotic_calcu

subroutine nonasymptotic_calcu(tau_vector, tc, sfr, rf, ft1, ft0, albedo, rf_calculated)

Definition: get_retrieval_uncertainty.f90:873

ch_xfr::ocis_id

integer ocis_id

Definition: ch_xfr.f90:51

get_retrieval_uncertainty::getctrefl_windspeeddiff

subroutine getctrefl_windspeeddiff(radius, optical_thickness, phase, wave_index, reflectancediff, meancloudtopreflw1, meancloudtopreflw2)

Definition: get_retrieval_uncertainty.f90:466

core_arrays::optical_thickness_37_final

real(single), dimension(:,:), allocatable optical_thickness_37_final

Definition: core_arrays.f90:41

get_retrieval_uncertainty::nonasymptotic_albedovar

subroutine nonasymptotic_albedovar(tau_vector, tc, sfr, rf, ft1, ft0, albedo, rf_lower, rf_middle, rf_upper)

Definition: get_retrieval_uncertainty.f90:922

libraryarrays

Definition: libraryarrays.f90:1

modis_numerical_module::linearinterpolation

real function, public linearinterpolation(X, Y, XX)

Definition: modis_numerical_module.f90:143

libraryinterpolates::int_reflectance_ice

real(single), dimension(:,:,:), allocatable int_reflectance_ice

Definition: libraryinterpolates.f90:4

libraryinterpolates::int_fluxdownwater_solar

real(single), dimension(:,:,:), allocatable int_fluxdownwater_solar

Definition: libraryinterpolates.f90:20

mod06_run_settings::band_0370

integer, parameter band_0370

Definition: mod06_run_settings.f90:26

science_parameters

Definition: science_parameters.f90:1

get_retrieval_uncertainty::getctref_constradius

subroutine getctref_constradius(optical_thickness_upper, optical_thickness_lower, radius, phase, surface_albedo, wave_index, reflectance_upper, reflectance_lower)

Definition: get_retrieval_uncertainty.f90:1366

get_retrieval_uncertainty::getctref_consttau

subroutine getctref_consttau(radius_upper, radius_lower, radius, optical_thickness, phase, surface_albedo, wave_index, partial_derivative)

Definition: get_retrieval_uncertainty.f90:1466

get_retrieval_uncertainty::reflpartialderivative

subroutine reflpartialderivative(upperLimit, lowerLimit, halfAbscissaInterval, abscissaIntervalStartingPoint, re, tau, phase, surface_albedo, wrtParam, wavelengthIdx, partialDerivative)

Definition: get_retrieval_uncertainty.f90:1181

libraryarrays::water_radii

real(single), dimension(:), allocatable water_radii

Definition: libraryarrays.f90:54

science_parameters::cox_munk

logical cox_munk

Definition: science_parameters.f90:11

libraryinterpolates

Definition: libraryinterpolates.f90:1

libraryarrays::library_taus

real(single), dimension(:), allocatable library_taus

Definition: libraryarrays.f90:40

libraryinterpolates::int_reflectance_ice_sdev

real(single), dimension(:,:,:), allocatable int_reflectance_ice_sdev

Definition: libraryinterpolates.f90:5

get_retrieval_uncertainty

Definition: get_retrieval_uncertainty.f90:1

generalauxtype

Definition: GeneralAuxType.f90:1

libraryinterpolates::int_reflectance_water_sdev

real(single), dimension(:,:,:), allocatable int_reflectance_water_sdev

Definition: libraryinterpolates.f90:5

science_parameters::ice_water_density

real, parameter ice_water_density

Definition: science_parameters.f90:17

libraryarrays::spherical_albedo_water

real(single), dimension(:,:,:), allocatable spherical_albedo_water

Definition: libraryarrays.f90:58

get_retrieval_uncertainty::getuncertainties

subroutine, public getuncertainties(thickness, re, water_path, phase, R1R2wavelengthIdx, meas_error, surface_albedo, transmittance_w1, transmittance_w2, delta_transmittance_w1, delta_transmittance_w2, transmittance_stddev_w1, transmittance_stddev_w2, emission_pw, emission_Tc, sigma_R37_pw, uTau, uRe, uWaterPath, xpoint, ypoint)

Definition: get_retrieval_uncertainty.f90:58

modis_numerical_module

Definition: modis_numerical_module.f90:1

interpolate_libraries::interpolate_wind_speed

subroutine, public interpolate_wind_speed(wspeed, data_in, data_out)

Definition: interpolate_libraries.f90:1490

mod06_run_settings

Definition: mod06_run_settings.f90:1

core_arrays::mean_delta_ozone

real mean_delta_ozone

Definition: core_arrays.f90:129

libraryinterpolates::int_fluxdownice_solar

real(single), dimension(:,:,:), allocatable int_fluxdownice_solar

Definition: libraryinterpolates.f90:21

void print(std::ostream &stream, const char *format)

Definition: PrintDebug.hpp:38

libraryarrays::spherical_albedo_ice

real(single), dimension(:,:,:), allocatable spherical_albedo_ice

Definition: libraryarrays.f90:67

get_retrieval_uncertainty::getradiibounds

subroutine, public getradiibounds(radiisize, radiidata, radius, upperbound, lowerbound)

Definition: get_retrieval_uncertainty.f90:1638

science_parameters::albedo_error

real, parameter albedo_error

Definition: science_parameters.f90:21

and

===========================================================================V5.0.48(Terra) 03/20/2015 Changes shown below are differences from MOD_PR02 V5.0.46(Terra)============================================================================Changes noted for V6.1.20(Terra) below were also instituted for this version.============================================================================V6.1.20(Terra) 03/12/2015 Changes shown below are differences from MOD_PR02 V6.1.18(Terra)============================================================================Changes from v6.1.18 which may affect scientific output:A situation can occur in which a scan which contains sector rotated data has a telemetry value indicating the completeness of the sector rotation. This issue is caused by the timing of the instrument command to perform the sector rotation and the recording of the telemetry point that reports the status of sector rotation. In this case a scan is considered valid by L1B and pass through the calibration - reporting extremely high radiances. Operationally the TEB calibration uses a 40 scan average coefficient, so the 20 scans(one mirror side) after the sector rotation are contaminated with anomalously high radiance values. A similar timing issue appeared before the sector rotation was fixed in V6.1.2. Our analysis indicates the ‘SET_FR_ENC_DELTA’ telemetry correlates well with the sector rotation encoder position. The use of this telemetry point to determine scans that are sector rotated should fix the anomaly occured before and after the sector rotation(usually due to the lunar roll maneuver). The fix related to the sector rotation in V6.1.2 is removed in this version.============================================================================V6.1.18(Terra) 10/01/2014 Changes shown below are differences from MOD_PR02 V6.1.16(Terra)============================================================================Added doi attributes to NRT(Near-Real-Time) product.============================================================================V6.1.16(Terra) 01/27/2014 Changes shown below are differences from MOD_PR02 V6.1.14(Terra)============================================================================Migrate to SDP Toolkit 5.2.17============================================================================V6.1.14(Terra) 06/26/2012 Changes shown below are differences from MOD_PR02 V6.1.12(Terra)============================================================================Added the doi metadata to L1B product============================================================================V6.1.12(Terra) 04/25/2011 Changes shown below are differences from MOD_PR02 V6.1.8(Terra)============================================================================1. The algorithm to calculate uncertainties for reflective solar bands(RSB) is updated. The current uncertainty in L1B code includes 9 terms from prelaunch analysis. The new algorithm regroups them with the new added contributions into 5 terms:u1:the common term(AOI and time independent) and

Definition: HISTORY.txt:126

get_retrieval_uncertainty::abscissainterval

subroutine abscissainterval(upperLimit, lowerLimit, halfAbscissaInterval, abscissaIntervalStartingPoint, abscissaIntervalLower, abscissaIntervalUpper)

Definition: get_retrieval_uncertainty.f90:1291

libraryarrays::ice_radii

real(single), dimension(:), allocatable ice_radii

Definition: libraryarrays.f90:63

get_retrieval_uncertainty::get_emission_values

subroutine get_emission_values(re, sigma37, phase, sigma37_val)

Definition: get_retrieval_uncertainty.f90:427

science_parameters::wind_speed_error

real, parameter wind_speed_error

Definition: science_parameters.f90:23

ch_xfr::c2_sensor_id

integer c2_sensor_id

Definition: ch_xfr.f90:50

interpolate_libraries

Definition: interpolate_libraries.f90:1

nonscience_parameters

Definition: nonscience_parameters.f90:4

nonscience_parameters::fillvalue_real

real, parameter fillvalue_real

Definition: nonscience_parameters.f90:13

get_retrieval_uncertainty::getclosestradius

subroutine getclosestradius(numberradii, radiidata, radius, index)

Definition: get_retrieval_uncertainty.f90:824

libraryinterpolates::int_reflectance_water

real(single), dimension(:,:,:), allocatable int_reflectance_water

Definition: libraryinterpolates.f90:4

science_parameters::solar_constant_37

real solar_constant_37

Definition: science_parameters.f90:35

mod06_run_settings::band_0065

integer, parameter band_0065

Definition: mod06_run_settings.f90:26

libraryarrays::number_waterradii

integer number_waterradii

Definition: libraryarrays.f90:5

core_arrays::optical_thickness_1621_final

real(single), dimension(:,:), allocatable optical_thickness_1621_final

Definition: core_arrays.f90:42

#define re

Definition: l1_czcs.c:695

phase

subroutine phase

Definition: phase.f:2

libraryinterpolates::int_fluxdownwater_sensor

real(single), dimension(:,:,:), allocatable int_fluxdownwater_sensor

Definition: libraryinterpolates.f90:20

science_parameters::lastinterp_wind_speed

real lastinterp_wind_speed

Definition: science_parameters.f90:47

get_retrieval_uncertainty::getctref_albedodiff

subroutine getctref_albedodiff(radius, optical_thickness, phase, wave_index, surface_albedo, reflectancediff, meancloudtopreflw1, meancloudtopreflw2)

Definition: get_retrieval_uncertainty.f90:722

libraryarrays::number_iceradii

integer number_iceradii

Definition: libraryarrays.f90:5

science_parameters::liquid_water_density

real, parameter liquid_water_density

Definition: science_parameters.f90:17

ch_xfr::oci_id

integer oci_id

Definition: ch_xfr.f90:52

modis_numerical_module::bisectionsearch

subroutine, public bisectionsearch(xx, x, jlo, jhi)

Definition: modis_numerical_module.f90:51

libraryinterpolates::int_fluxdownice_sensor

real(single), dimension(:,:,:), allocatable int_fluxdownice_sensor

Definition: libraryinterpolates.f90:21

get_retrieval_uncertainty::get_refl_windvector

subroutine get_refl_windvector(radius, optical_thickness, phase, wave_index, sdev_refl)

Definition: get_retrieval_uncertainty.f90:642

libraryinterpolates::int_reflectance_water_wspeed

real(single), dimension(:,:,:,:), allocatable int_reflectance_water_wspeed

Definition: libraryinterpolates.f90:7

radius

void radius(double A)

Definition: proj_report.c:132

abs

#define abs(a)

Definition: misc.h:90

libraryinterpolates::int_reflectance_ice_wspeed

real(single), dimension(:,:,:,:), allocatable int_reflectance_ice_wspeed

Definition: libraryinterpolates.f90:7

get_retrieval_uncertainty::sensitivitypartialderivatives

subroutine sensitivitypartialderivatives(R1R2wavelengthIdx, re, tau, phase, surface_albedo, partialDerivTauWrtR1AtConstR2, partialDerivTauWrtR2AtConstR1, partialDerivReWrtR1AtConstR2, partialDerivReWrtR2AtConstR1)

Definition: get_retrieval_uncertainty.f90:957

generalauxtype::debugprn

logical debugprn

Definition: GeneralAuxType.f90:36

get_retrieval_uncertainty::init_half_radii

subroutine, public init_half_radii

Definition: get_retrieval_uncertainty.f90:15

transpose

void transpose(float *in[], float *out[])

Definition: get_zenaz.c:97