NASA Ocean Color

ocssw V2022

 module retrieval_solution_logic
  
   implicit none
  
   integer, parameter :: top_nk_asl_contour_ice = 1, &
     top_nk_asl_contour_water = 2
   
 contains
  
   ! ssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss
 subroutine find_zero_crossings (residual, crossing_vector, crossing_num)
   ! ssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss
  
   !   *** S. Platnick, 17-Oct-2005
   !   *** G. Wind 20-Oct-2005 -- removed the dynamic memory 
   !  allocation to speed things up a bit
   !
   !  residual
   !  crossing_vector
   !  crossing_num
   !
   implicit none
  
   REAL, intent(in)     :: residual(:)
   INTEGER, intent(out) :: crossing_num, crossing_vector(:)
  
   !   local variables
   ! INTEGER, allocatable, dimension(:) :: k
   integer :: k(18)
  
   INTEGER :: i,length
  
   length = size(residual)
   ! allocate (k(length))
  
   ! G.Wind 12.19.05: Replaced the where statement below with a do loop.  This 
   ! eliminates ambiguities arising from the use of nonconformable arrays. 
   k =   1
   !   where(residual<0) k=-1
  
   do i=1, length
     if ( residual(i) < 0 ) k(i) = -1
   end do  
   ! crossing index in vector is defined for the smaller index 
   !  between which the zero lies
   crossing_num=0; crossing_vector=0.
   Do i=1,length-1
     if( abs(k(i+1)-k(i)) > 1) then!{
       crossing_num = crossing_num + 1
       crossing_vector(crossing_num) = i
     end if
   End do
     
   ! deallocate (k)
     
 end subroutine find_zero_crossings
  
  
   ! ssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss
 subroutine find_local_minima (residual, local_min_vector, local_min_num)
   ! ssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss
  
   !   *** S. Platnick, 17-Oct-2005
   !   This give the local minima for the absolute residual, not the residual**2.
   !   10-18-05: modify so that minima at end points are not 
   !  considered local minima
   !    to the extent that they do not indicate a change in curvature.
   !  residual
   !  local_min_vector
   !  local_min_num
  
   implicit none
  
   REAL, intent(in)     :: residual(:)
   INTEGER, intent(out) :: local_min_num, local_min_vector(:)
  
   !   local variables
   INTEGER i,length
  
   length = size(residual)
  
   local_min_num=0; local_min_vector=0.
  
   if (length > 2) then!{
     Do i=2,length-1
       if( residual(i) < residual(i-1) .AND. residual(i) < residual(i+1) ) then!{
         local_min_num = local_min_num + 1
         local_min_vector(local_min_num) = i
       end if
     End do
   end if
  
 end subroutine find_local_minima
  
   ! ssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss
 subroutine find_local_maxima (residual, local_max_vector, local_max_num)
   ! ssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss
  
   !   *** S. Platnick, 17-Oct-2005
   !   This give the local maxima for the absolute residual, not the residual**2.
   !   10-18-05: modify so that minima at end points are not 
   !   considered local maxima
   !    to the extent that they do not indicate a change in curvature.
   !
   !  residual
   !  local_max_vector
   !  local_max_num
   !
   implicit none
  
   REAL, intent(in)     :: residual(:)
   INTEGER, intent(out) :: local_max_num, local_max_vector(:)
  
   !   local variables
   INTEGER i,length
  
   length = size(residual)
  
   local_max_num=0; local_max_vector=0.
  
   if (length > 2) then!{
     Do i=2,length-1
       if( residual(i) > residual(i-1) .AND. residual(i) > residual(i+1) ) then!{
         local_max_num = local_max_num + 1
         local_max_vector(local_max_num) = i
       end if
     End do
   end if
  
 end subroutine find_local_maxima
  
   ! sssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss
 subroutine solution_re (re, residual, crossing_vector, crossing_num,    &
       local_min_vector, local_min_num, local_max_vector, local_max_num, &
       retrieval, quality )
   ! sssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss
  
   use science_parameters 
  
   !  *** S. Platnick, 17-Oct-2005
   !      Using input from sep routines 'find_zero_crossings', 
   !          'find_local_minima', and 'find_local_maxima'
   !
   !  *** S. Platnick, 7-March-2006
   !
   !   Modification of collection 5 solution logic:
   !   1. process_zero_crossing_residual = .FALSE., i.e., re for 
   !      pixels with no zero crossings are set to fill
   !   2. solve_for_zero_crossing modified to use spline interpolation 
   !      and iterative root finder
   !   3. No longer need 'find_local_minima' and 'find_local_maxima' 
   !      routines, though left intact for now since
   !      they provide useful information that might still be needed sometime.
   !   4. No retrieval when crossing_num > 2
  
   implicit none
  
   integer, intent(in)    :: crossing_vector(:), crossing_num, &
     local_min_vector(:), local_min_num, local_max_vector(:), local_max_num
   real, intent(in)       :: residual(:), re(:)
  
   real, intent(out)      :: retrieval
   integer*1, intent(out) :: quality
  
   integer                :: residual_crossing_index, index1(1)
   real                   :: fillval_r4, resid_lo, resid_hi, slope_1, slope_2
   logical                :: process_zero_crossing_residual
  
  
   fillval_r4 = invalid_attempted_but_failed_re; residual_crossing_index = 0
  
   ! ===================================
   ! SOLUTION LOGIC QUALITY ASSIGNMENTS: 
   ! ===================================
   ! Can be used for diagnostic purposes and setting retrieval 
   ! confidence QA in core science (not presently being
   !  done as of this writing)
   !
   !   quality =  -6  => spline interpolation, iteration /< iternum-5
   !   quality =  -5  => crossing_num < 0 (should never occur, but 
   !         included for completeness)
   !   quality =  -4  => crossing_num = 0 and process_zero_crossing_residual
   !        =.FALSE.
   !   quality =  -3  => crossing_num = 0 and process_zero_crossing_residual
   !        =.TRUE.
   !   quality =  -2  => crossing_num > 2
   !   quality =  -1  => DEFAULT: this quality index should never occur, 
   !         and if it does something went awry
   !   quality =   0  => occur for truncated residual vector having fill value(s)
   !   quality =   1  => linear interpolation w/maxval(abs(residual)) > 
   !         resid_max_for_spline
   !   quality =   2  => linear interpolation when the root lies in the 
   !         last interval of the residual vector
   !   quality =   3  => spline interpolation (or no interpolation), 
   !         and iteration < iternum-5
   !
   ! If ever used in MOD06 for assigning confidence QA then possible 
   !         negative assignments are:
   !    -6->3 or 2;   -5->0;   -4->0;   -3->1 or 0;   
   !     -2->1 or 0 depending on solution logic
   
   ! --------------------------------------------------------------
   ! .TRUE. corresponds to the logic that had apparently been in 
   !  place since the beginning of the c5 deliveries
   ! 2 March 06: Changed to .FALSE., i.e., for a non-zero crossing 
   !  residual, return a fill value. Do not return with
   !   the library radius corresponding to smallest residual as the solution.
   ! -------------------------------------------------------------
  
   process_zero_crossing_residual = .false.   
  
   ! -------------------------
   ! Default retrieval output:
   ! -------------------------
   retrieval = fillval_r4; quality = -1
  
   ! ------------------------------------------
   ! Check for a proper value for crossing_num:
   ! ------------------------------------------
  
   if (crossing_num < 0) then
     quality = -5
     return
   end if
   ! ---------------------------------------------------------------------------
   ! If there are more than 2 roots for the residual vector return a fill value:
   ! ---------------------------------------------------------------------------
  
   if (crossing_num > 2) then
     quality = -2
     return
   end if
  
   ! -----------------------------------------------------------------
   ! The calling routine in C5 truncates the residual array to get 
   !  rid of contiguous fill values starting with the largest radii, and
   !  then passes the useful max radii to the re solution routines. 
   !  However, the calling routine does not look for remaining fill values
   !  that are elsewhere in the residual vector but with non-fill 
   !  neighbors. If such a residual vector remains, then return with a
   !  fill value for effective radius.
   !
   !  However, I can't determine if resdiual vector is initialized 
   !  to fillvalue_real !?
   !
   ! ----------------------------------------------------------------
   
   !  if (min(residual) == fillvalue_real) then
   !     quality = 0
   !     return
   !  end if
   
   ! ------------------------------------------------------------------------
   ! If there are NO roots for the residual vector (i.e., no zero crossings):
   ! ------------------------------------------------------------------------
   IF (crossing_num == 0) THEN
     if (process_zero_crossing_residual) then
  
       index1 = minloc(abs(residual))
       retrieval = re(index1(1))
       quality = -3
     else
       quality = -4
     end if
     return
   END IF
   
   ! --------------------------------------------------
   ! If there are 1 or 2 roots for the residual vector:
   ! --------------------------------------------------
  
   IF (crossing_num == 1 .OR. crossing_num == 2) THEN
  
   ! If there are 2 zero crossings for the residual vector 
   !  then choose the crossing having the negative slope
   !  regardless of the cloud phase. A derivative of residual 
   !  w.r.t. re < 0 corresponds to the physical situation
   !  where single scattering albedo decreases with re (and no 
   !  significant change in effective radius or extinction).
  
   If (crossing_num == 1) Then
     residual_crossing_index = crossing_vector(crossing_num)
  
   Else
     slope_1 = residual(crossing_vector(1)+1) - residual(crossing_vector(1))
     slope_2 = residual(crossing_vector(2)+1) - residual(crossing_vector(2))
     !   resid_lo = 0.5 * ( residual(crossing_vector(1)) + &
     !      residual(crossing_vector(1) + 1) )
     !   resid_hi = 0.5 * ( residual(crossing_vector(2)) + &
     !     residual(crossing_vector(2) + 1) )
  
     !  For crossing_num =2, both slopes can't be of the same sign. 
     !  If they are something's not correct, but choose larger residual.
     if ( slope_1 <= 0. ) then
       residual_crossing_index = crossing_vector(1)
     end if
     if ( slope_2 <= 0. ) then 
       residual_crossing_index = crossing_vector(2)
     end if
  
  
   End If
   
   ! it could potentially happen that slopes are > 0, then this 
   !  would be a segfault because residual_crossing_index 
   ! defaults to 0. --- G.Wind 12.15.11
   if (residual_crossing_index > 0) &
     call solve_for_zero_crossing (re, residual, residual_crossing_index, &
       fillval_r4, local_min_vector, local_min_num, local_max_vector, &
       local_max_num, retrieval, quality)
   END IF    
  
 end subroutine solution_re
  
   ! sssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss
 subroutine solve_for_zero_crossing (re, residual, &
     residual_crossing_index, fillval_r4, &
     local_min_vector, local_min_num, local_max_vector, local_max_num, &
     retrieval, quality)
   ! ssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssssss
  
   use generalauxtype
   use spline_module
  
   !  *** S. Platnick, 7-March-2006
   !
   !   Modification of collection 5 solution logic:
   !
   !   1. Use natural spline interpolation instead of quadratic(s) fits.
   !      The piece-wise quadratic approach was producing discontinuities 
   !      in the re histograms. Epecially for ice clouds where
   !      the residual vector often has changes in curvature which exacerbates 
   !      the piece-wise approach. The advantage of the 
   !      quadratic fits is that they are analytic and computationally 
   !      fast; the spline fit must use an iterative procedure for
   !      for finding the root of the residual vector.
   !
   !   2. Also, see comments in solution_re routine.
  
   implicit none
  
   real,    intent (in)   :: re(:), residual(:), fillval_r4
   integer, intent (in)   :: residual_crossing_index,  &
     local_min_vector(:), local_min_num, local_max_vector(:), local_max_num
  
   real,      intent(inout) :: retrieval
   integer*1, intent(inout) :: quality
    
   integer i, icount, iternum
   real    y2(size(re)), z, delta_resid, delta_resid_fraction, &
     resid_min_allowable, resid_max_for_spline
   real    yl(2), xl(2), re_lower, re_upper, re_iter
    
   character*15 numerical_type
  
   retrieval=fillval_r4
  
   !  -------------------------------------------------------------
   !  If the root is close enough to a library re, don't bother 
   !   interpolating for the solution
   !   (also, spline bisection search was found to fail for a residual 
   !   index of zero which did occur for test granules)
   !  --------------------------------------------------------------
  
   resid_min_allowable = 0.0001
   if (abs(residual(residual_crossing_index)) < resid_min_allowable) then
     retrieval = re(residual_crossing_index)
     quality = 3
     return
   end if
   if (abs(residual(residual_crossing_index+1)) < resid_min_allowable) then
     retrieval = re(residual_crossing_index+1)
     quality = 3
     return
   end if
  
   !  -----------------------------------------------------------------
   !  If bounding library points are so near to one another as to be 
   !    representationally equal, no need to interpolate
   !  ----------------------------------------------------------------
  
   if (real_s_equal(residual(residual_crossing_index), &
     residual(residual_crossing_index+1))) then
     retrieval = re(residual_crossing_index)
     quality = 3
     return
   end if
  
   !  -----------------------------------------------------------------
   !  Use either a linear or spline interpolation to determine the 
   !  residual vector root:
   !  ------------------------------------------------------------------
   numerical_type = 'spline'
  
   !  The following 2 Do Loops set linear interpolation for residual 
   !  roots near a local minimum or maximum. Commented out
   !   to use spline interpolation in these instances. Use of linear 
   !   interpoltion gave some re histogram discontinuities in
   !   test granules, at least for ice clouds where local min and max 
   !   seemed to occur more often.
  
   !   Do i=1,local_max_num
   !     if (residual_crossing_index == local_max_vector(i) .OR. &
   !       residual_crossing_index+1 == local_max_vector(i)) then
   !       numerical_type = 'linear'
   !       exit
   !     end if
   !   End do
   !   IF (numerical_type /= 'linear') THEN
   !     Do i=1,local_min_num
   !       if (residual_crossing_index == local_min_vector(i) .OR. &
   !         residual_crossing_index+1 == local_min_vector(i)) then
   !         numerical_type = 'linear'
   !         exit
   !       end if
   !     End do
   !   End If
  
  
   !  Use linear interpolation if the root is before end of re vector. 
   !   Especially for ice clouds, the spline appears to give undue 
   !   curvature to fit the large distance from 60 to 90 microns. For 
   !   water clouds, the accuracy in establishing the root between
   !   28 and 30 microns isn't relevant. Use of this logic did not 
   !   appear to cause any undue features in the re histogram.
  
   ! The following innocuous check on 'linear' remains for historical 
   !  reasons in case the block immediately above ever gains favor
  IF (numerical_type /= 'linear') THEN
    if (residual_crossing_index == size(re)-1) then
      numerical_type = 'linear'
      quality = 2
    end if
   End If
  
   !  Use linear interpolation if there are large discontinuities in 
   !  the residual which would give rise to spline oscillations
   !
   !   Example print output for 3.7 um residual, MYD Katrina granule: 
   !   iterational doesn't converge because the spline does not provide
   !   a root between 30 and 35 um points due to large oscillation set 
   !   up by large residual at larger radii.
   !
   !   *** number of spline iterations: icount = 30  residual_crossing_index= 6
   !   re:        5.0   10.0   15.0   20.0   25.0   30.0   35.0   
   !   40.0   45.0   50.0   55.0   60.0   90.0
   !   residual: 3.9147  2.1260  0.6555  0.0979  -0.0182  -0.0686  0.0729  
   !    0.4781  2.1382  10000.0  383.20  10000.0  10000.0
  
   resid_max_for_spline = 100.
   IF (numerical_type /= 'linear') THEN
     if ( maxval(abs(residual)) > resid_max_for_spline) then
       numerical_type = 'linear'
       quality = 1
     end if
   End If
  
   !  ------------------------------------------------------------------
   !  Implement interpolation scheme and solve for residual vector root:
   !  ------------------------------------------------------------------
  
   IF (numerical_type == 'linear') THEN
     yl(1) = residual(residual_crossing_index); yl(2) = &
       residual(residual_crossing_index+1)
     if (sign(1.,yl(1))*sign(1.,yl(2)) <= 0.) then
       call linear_interpolate_for_root (yl, &
         re(residual_crossing_index:residual_crossing_index+1), retrieval)
     end if
   END IF
  
   IF (numerical_type == 'spline') THEN
     call spline (size(re), re, residual, y2)
  
     iternum = 30; icount = 1; delta_resid_fraction=0.001
  
     delta_resid = abs( delta_resid_fraction * &
       max(residual(residual_crossing_index), &
       residual(residual_crossing_index+1)) )
     if(delta_resid < resid_min_allowable) delta_resid = resid_min_allowable
     re_lower = re(residual_crossing_index)
     re_upper = re(residual_crossing_index+1)
     re_iter = 0.5*(re_lower + re_upper)
   
     !  The interpolation finds the y (residual) for the input x (re), so 
     !  they iterate till the y they find is close to 0.  Now, if the 
     !  relation is reversed, and y is re and x is residual, you could do 
     !  the spline 1 time for x = 0.  All the min, max points could be used to 
     !  isolate the correct segment so a function wuld be analyzed.`
     DO WHILE (icount < iternum)
       z= splint(size(re), re_iter, re, residual, y2)
  
       if(abs(z) < delta_resid) then
         retrieval = re_iter
         exit
       end if
     
       if(sign(1.,z)*sign(1.,residual(residual_crossing_index)) > 0.) then
         re_lower = re_iter
       else
         re_upper = re_iter
       end if  
       re_iter = 0.5*(re_lower + re_upper)
       icount = icount + 1
     END DO
  
     ! If consistently near the max iteration number, it is possible a 
     !  higher max iteration number might be more suitable.
     if(icount < iternum-5) then
       quality =  3
     else
       quality = -6
     end if
  
   END IF
  
 end subroutine solve_for_zero_crossing
  
  
   logical function is_water_phase(library_radii)
  
   use libraryarrays
  
   real(single), intent(in) :: library_radii(:)
  
   ! if it's not ice phase, it's water
   is_water_phase = (size(library_radii) .ne. size(ice_radii)) 
  
   ! NOTE:  if at some point ice & water libs have equal sizes, 
   !  this needs reconsideration
  
   end function is_water_phase
  
 subroutine ray_corr_nearest(refl_source, As, iw, tau, re, Pc,  &
     sfr, fti1, fti0, fluxup_solar, fluxup_sensor, &
     theta, theta0, phi, location, crefl)
   !    don't know what this does specifically
   !  refl_source
   !  As
   !  iw
   !  tau
   !  re
   !  Pc
   !  sfr
   !  fti1
   !  fti0
   !  fluxup_solar
   !  fluxup_sensor
   !  theta
   !  theta0
   !  phi
   !  location  I  radius location
   !  crefl  O returned reflectance?
  
   use science_parameters
   use libraryarrays
   use mod06_run_settings
   use modis_numerical_module
  
   real, intent(in) :: tau, re, Pc, sfr, fti1, fti0, theta, theta0, phi, &
     As, refl_source
   integer, dimension(2), intent(in) :: location
   integer, intent(in) :: iw
   real, dimension(:,:,:,:), intent(in) :: fluxup_solar, fluxup_sensor
   real, intent(inout) :: crefl
   
   real, parameter :: Ps = 1013.
   real, parameter :: Cm = 0.84
   real, parameter :: taur0(2) = (/0.044, 0.025/)
  
   real :: taur, mu0, mu, x, pray, refl_s
   real :: fti1_as, fti0_as
   real :: angles(2), functionvalues(2)
   real :: fluxup_solar_tau, fluxup_sensor_tau
   integer :: lowbound, highbound
  
   lowbound = 0  ! WDR-UIV
   highbound = 0  ! WDR-UIV
  
   taur = taur0(iw) * pc / ps
   mu0 = cos(d2r*theta0)
   mu = cos(d2r*theta)
  
   x = -mu0*mu + sqrt((1.-mu0**2)*(1.-mu**2))*cos(d2r*phi)
     pray = 3.*(1.+x**2)/4.0
  
   !
   !  Interpolation for solar zenith angles and scaled optical thickness
   !
   call bisectionsearch(library_fluxsolarzenith, mu0, lowbound, highbound)
   angles(1) = library_fluxsolarzenith(lowbound)
   angles(2) = library_fluxsolarzenith(highbound)
  
   if (location(1) == 1) then 
     functionvalues = 0.
   else 
     functionvalues(1) = fluxup_solar(lowbound,location(1)-1,1,location(2))
     functionvalues(2) = fluxup_solar(highbound,location(1)-1,1,location(2))
   endif
  
   fluxup_solar_tau = linearinterpolation(angles,functionvalues,mu0)
  
   !
   !  Interpolation for sensor zenith angles and scaled optical thickness
   !
   call bisectionsearch(library_fluxsensorzenith, mu, lowbound, highbound)
   angles(1) = library_fluxsensorzenith(lowbound)
   angles(2) = library_fluxsensorzenith(highbound)
  
   if (location(1) == 1) then 
     functionvalues = 0.
   else 
     functionvalues(1) = fluxup_sensor(lowbound,location(1)-1,1,location(2))
     functionvalues(2) = fluxup_sensor(lowbound,location(1)-1,1,location(2))
   endif
   fluxup_sensor_tau = linearinterpolation(angles,functionvalues,mu)
  
   fti1_as = fluxup_sensor_tau
   fti0_as = fluxup_solar_tau
  
   if (as > 0.) then 
     fti0_as = fti0_as + fti0*as*(1-sfr)/(1 - sfr*as)
     fti1_as = fti1_as + fti1*as*(1-sfr)/(1 - sfr*as)
   endif
  
   refl_s = taur*pray/(4.*mu0*mu) + 0.5*taur* fti1_as * exp(-taur/mu )/mu0 + &
     0.5*taur* fti0_as  * exp(-taur/mu0)/mu
   crefl = (refl_source - refl_s)*exp(cm*taur*(1./mu0 + 1./mu))
  
 end subroutine ray_corr_nearest
  
 subroutine solveretrieval_nearest(xx_pt,yy_pt,Ram, Rbm, twobands, &
     radii, tau, re, lib_dist, phase_liquid, Ram_corr,quality_in,CH37_IDX, &
     CTopT,CH37_NUM, platFormName)
   !  xx_pt, IN  x, loc of point in the chunk
   !  yy_pt  IN  y loc of point in the chunk
   !  Ram  IN  non-abs reflectance
   !  Rbm   IN  absorbing band reflectance
   !  twobands  IN  size 2 indicies of non-abs and absorbing indiciew
   !  radii   IN  array of table re values
   !  tau    OUT  retrieved optical thickness
   !  re   OUT  retrieved final radius
   !  lib_dist  OUT  I believe a uncertainty measure (root sum square) for the 
   !       derived re / tau?
   !  phase_liquid   IN  switch for liquid phase if true
   !  Ram_corr    
   !  quality_in
   !  --- extra args? ---
   !  CH37_IDX
   !  CTopT,CH37_NUM
   !  platFormName
   !
   !  alternate solution exterior, just the boundary (for water, re=4, 
   !  re=30 and tau = 158.8)
   !  for ice re=5,re=60 and tau = 158.8 (top, bottom and rightmost
  
   use generalauxtype
   use core_arrays
   use science_parameters
   use nonscience_parameters
   use libraryinterpolates
   use libraryarrays
   use planck_functions
  
   integer, intent(in) ::xx_pt,yy_pt
   real, intent(in) :: Ram, Rbm
   real, dimension(:), intent(in) :: radii
   integer, dimension(:), intent(in) :: twobands
   real, intent(inout) :: tau, re, lib_dist, Ram_corr
   logical, intent(in) :: phase_liquid
     integer, intent(in) :: quality_in
     CHARACTER(len=*),INTENT(IN),OPTIONAL::platFormName
   INTEGER,INTENT(IN),OPTIONAL::CH37_NUM,CH37_IDX
   REAL, INTENT(IN), OPTIONAL::CTopT
   
   REAL :: Pc, theta, theta0, phi
   real::CTopT_lcl, intCTop,intSurf
   integer :: size_tau, size_re
  
   size_tau = number_taus + 1
   size_re = size(refliba, 2)
   
   pc = cloud_top_pressure(xx_pt,yy_pt)
   theta = sensor_zenith_angle(xx_pt,yy_pt)
   theta0= solar_zenith_angle(xx_pt,yy_pt)
   phi=relative_azimuth_angle(xx_pt,yy_pt)
   !
   !  WDR have the CTopT go to a local variable and handle presence then
   IF(PRESENT(ctopt)) THEN
     ctopt_lcl = ctopt
   else
     ctopt_lcl = 0.
   endif
   !  3.7 not in PACE, so return later
   IF(PRESENT(ch37_idx) .AND. PRESENT(ctopt) .and. PRESENT(platformname) &
     .and. PRESENT(ch37_num) .and. ctopt_lcl > 0.0)THEN     
     intctop= modis_planck(platformname, ctopt_lcl, ch37_num, 1)
     intsurf = modis_planck(platformname, surface_temperature(xx_pt,yy_pt), &
       ch37_num, 1)
        
     if( (.NOT. cox_munk))THEN
       if(phase_liquid)then
         call calc37radianceliblamb(intctop,intsurf, &
           thermal_correction_oneway(1), &
           thermal_correction_twoway(1), theta0, &
           spherical_albedo_water(:,ch37_idx,:), &
           int_fluxdownwater_sensor(:,ch37_idx,:), &
           int_fluxdownwater_solar(:,ch37_idx,:), &
           int_fluxupwater_sensor(:,ch37_idx,:))
       else
         call calc37radianceliblamb(intctop,intsurf,  &
           thermal_correction_oneway(1), thermal_correction_twoway(1), theta0, &
           spherical_albedo_ice(:,ch37_idx,:), &
           int_fluxdownice_sensor(:,ch37_idx,:), &
           int_fluxdownice_solar(:,ch37_idx,:), &
           int_fluxupice_sensor(:,ch37_idx,:))
       endif
             
     elseif(cox_munk)THEN
       if(phase_liquid)then
         call calc37radiancelibcm(intctop,intsurf, &
           thermal_correction_oneway(1), &
           thermal_correction_twoway(1), theta0, &
           int_cloud_emissivity_water(:,1,:), &
           int_surface_emissivity_water(:,1,:))
       else
         call calc37radiancelibcm(intctop,intsurf,  &
           thermal_correction_oneway(1), thermal_correction_twoway(1), &
           theta0, int_cloud_emissivity_ice(:,1,:), &
           int_surface_emissivity_ice(:,1,:))
       endif
     endif
   ENDIF   ! End 3.7 related logic
  
   if(quality_in == -4 .OR. quality_in == -3)then
     CALL asl_boundary(ram, rbm, twobands, radii, tau, re, lib_dist, pc, & 
       theta, theta0, phi, phase_liquid, ram_corr)                   
   else
     CALL asl_interior(ram, rbm, twobands, radii, tau, re, lib_dist, pc, & 
       theta, theta0, phi, phase_liquid, ram_corr)
   endif
                
   return
  
 end subroutine solveretrieval_nearest
  
 subroutine calc37radianceliblamb(intensity,intensity_g, &
     thermal_trans_1way, thermal_trans_2way, solar_zenith, &
     sfr,fti1,fti0,fri1)
  
   use generalauxtype
   use science_parameters
   use libraryarrays, only : reflibb,rad37lib
  
   implicit none
  
   real, intent(in)     ::   intensity,intensity_g,solar_zenith
  
   real, intent(in)      :: thermal_trans_2way,thermal_trans_1way
   real(single), intent(in)  :: sfr(:,:),  fti1(:,:), fti0(:,:), fri1(:,:)
 !  rfi(:,:)
  
   integer              :: total_taus, radii, wavelengths,i,j
   real                 :: absorbing_albedo
   real                 :: tc, local_trans_1way, local_trans_2way
  
   if (thermal_trans_1way < 0. .or. thermal_trans_1way > 1.) then!{
     local_trans_1way = 1.
   else!}{
     local_trans_1way = thermal_trans_1way
   endif!}
  
   if (thermal_trans_2way < 0. .or. thermal_trans_2way > 1.) then!{
     local_trans_2way = 1.
   else!}{
     local_trans_2way = thermal_trans_2way
   endif!}
  
  
   !populate rad37lib here, in reflectance units
   absorbing_albedo = reflibb(1,1)  
   
   rad37lib(1,:) = (1.0 - absorbing_albedo) * intensity_g 
   rad37lib(2:,:) = (1.0- fti1(:,:) - fri1(:,:)) * intensity + &
     fti1(:,:) * (1.0 - absorbing_albedo) * intensity_g / &
     (1.0 - absorbing_albedo * sfr(:,:))
   rad37lib(:,:) = rad37lib(:,:) * local_trans_1way * pi/ &
     ( cos(solar_zenith*d2r)*solar_constant_37)                                 
         rad37lib(:,:) = rad37lib(:,:) + local_trans_2way * reflibb(:,:)
           
 end subroutine calc37radianceliblamb
  
 subroutine calc37radiancelibcm(intensity,intensity_g,     &
     thermal_trans_1way, thermal_trans_2way, solar_zenith,  &
     cl_emis,sf_emis)
  
   use generalauxtype
   use science_parameters
   use libraryarrays, only : reflibb,rad37lib
  
   implicit none
  
   real, intent(in) :: intensity,intensity_g,solar_zenith
   real, intent(in)      :: thermal_trans_2way,thermal_trans_1way
   real(single), intent(in)     ::  cl_emis(:,:), sf_emis(:,:)
                             
   !   real,intent(out)::rad037lib(:,:),rtherm037lib(:,:)
                            
  
   integer              :: total_taus, radii, wavelengths,i,j
   real                 :: absorbing_albedo
   real                 :: tc, local_trans_1way, local_trans_2way
  
   if (thermal_trans_1way < 0. .or. thermal_trans_1way > 1.) then!{
     local_trans_1way = 1.
   else!}{
     local_trans_1way = thermal_trans_1way
   endif!}
  
   if (thermal_trans_2way < 0. .or. thermal_trans_2way > 1.) then!{
     local_trans_2way = 1.
   else!}{
     local_trans_2way = thermal_trans_2way
   endif!}
  
   !populate rad37lib here, in reflectance units
  
   rad37lib(:,:) = cl_emis(:,:) * intensity + sf_emis(:,:) * intensity_g
   rad37lib(:,:) = rad37lib(:,:) * local_trans_1way * pi/ &
     ( cos(solar_zenith*d2r)*solar_constant_37)
   rad37lib(:,:)  = rad37lib(:,:) + local_trans_2way * reflibb(:,:)
                       
 end subroutine calc37radiancelibcm
  
 subroutine asl_interior(Ram, Rbm_in, bands, radii, tau, re, lib_dist, Pc, & 
                    theta, theta0, phi, phase_liquid, Ram_corr)
   !
   !  For solutions inside the LUT with multiple solutions - derive one and
   !  essentialy an uncertainty in the form of a distance - This is for
   !alternate solution interior                   
   !
   !  Ram
   !  Rbm_in
   !  bands
   !  radii
   !  tau
   !  re
   !  lib_dist
   !  Pc
   !  theta
   !  theta0
   !  phi
   !  phase_liquid
   !  Ram_corr
   !
   use generalauxtype
   use libraryarrays
   use libraryinterpolates
   use nonscience_parameters
   use mod06_run_settings
   use science_parameters
  
   real, intent(in) :: Ram, Rbm_in, Pc, theta, theta0, phi
   real, dimension(:), intent(in) :: radii
   integer, dimension(:), intent(in) :: bands
   real, intent(inout) :: tau, re, lib_dist, Ram_corr
   logical, intent(in) :: phase_liquid
  
   
   real, dimension(:,:), allocatable ::  rel_residual,reflibBL
   integer :: size_tau, size_re, loc(2), rel_loc(2)
   logical :: re_altered
   
   real :: Ramc, crefl,Rbm
   integer :: iter, max_iter
  
   real :: albedoA, albedoB,normConst
   real :: sfr, fti1, fti0
  
   size_tau = number_taus + 1
   size_re = size(refliba, 2)
   allocate(rel_residual(size_tau, size_re),reflibbl(size_tau, size_re))
   
   reflibbl(:,:) = reflibb(:,:)
   rbm = rbm_in
       
   if( bands(2) == band_0370)then 
     normconst = pi/(cos(theta0*d2r)*solar_constant_37)
     rbm = rbm_in * normconst
     reflibbl(:,:) = rad37lib(:,:)
   endif
  
   albedoa = refliba(1,1)
   albedob = reflibbl(1,1)
  
  
 #if ASTER_INST
   if (bands(1) == band_0065 .or. bands(1) == band_0086) then 
 #else
   if (bands(1) == band_0065) then 
 #endif  
     max_iter = 3
   else 
     max_iter = 1
   endif
   
   ramc = ram
   iter = 1
   
   do while (iter <= max_iter) 
  
     ! the point has no chance, it is darker than the surface albedo
     if (ramc < albedoa) then
       tau = fillvalue_real
       re  = fillvalue_real
       lib_dist = max_cost_function
       ram_corr = ramc
       deallocate (rel_residual, reflibbl )  ! WDR add the reflibBL deallocate
       return
     endif
     ! collect the value of cost function
     rel_residual = sqrt( (refliba - ramc)**2 + (reflibbl-rbm)**2 ) / &
       sqrt(ramc**2 + rbm**2)
     loc = minloc(rel_residual)  
     lib_dist = rel_residual(loc(1), loc(2))
     re = radii(loc(2))
   
     ! don't want any 2's and such 
     if (re < 4. ) re = 4.
     ! the library taus don't have 0., so we need to explicitly count for 
     !  it. But we can't put in a 0. straight because
     !  L3 code can't handle an explicit 0. value. 
     if (loc(1) == 1) then 
       tau = 0.01
     else
       tau = library_taus(loc(1)-1)
     endif
   
     ! call rayleigh here
 #if ASTER_INST
     if (bands(1) == band_0065 .or. bands(1) == band_0086) then 
 #else
     if (bands(1) == band_0065) then 
 #endif  
       if (phase_liquid) then 
         if (loc(1) == 1) then 
           sfr = 0.
           fti1 = 1.
           fti0 = 1.
         else
           sfr = spherical_albedo_water(loc(1)-1,bands(1),loc(2))
           fti1 = int_fluxdownwater_sensor(loc(1)-1,bands(1),loc(2))
           fti0 = int_fluxdownwater_solar(loc(1)-1,bands(1),loc(2))
         endif
           
         call ray_corr_nearest (ram, albedoa, bands(1), tau, re, pc, &
           sfr, fti1, fti0, flux_up_water_solar, flux_up_water_sensor, &
           theta, theta0, phi, loc, crefl)
       else 
         if (loc(1) == 1) then 
           sfr = 0.
           fti1 = 1.
           fti0 = 1.
         else
           sfr = spherical_albedo_ice(loc(1)-1,bands(1),loc(2))
           fti1 = int_fluxdownice_sensor(loc(1)-1,bands(1),loc(2))
           fti0 = int_fluxdownice_solar(loc(1)-1,bands(1),loc(2))
         endif
         
         call ray_corr_nearest (ram, albedoa, bands(1), tau, re, pc, &
           sfr, fti1, fti0, flux_up_ice_solar, flux_up_ice_sensor, &
           theta, theta0, phi, loc, crefl)
       endif             
       ramc = crefl
   
     endif
     
     iter = iter + 1
  
   end do
   
   ram_corr = ramc
  
   if (ramc < albedoa) then
     tau = fillvalue_real
     re  = fillvalue_real
     lib_dist = max_cost_function
     deallocate (rel_residual, reflibbl )  !  WDR add this to complete deallocs
     return
   endif
   ! collect the value of cost function
   rel_residual = sqrt( (refliba - ramc)**2 + (reflibbl-rbm)**2 ) / &
     sqrt(ramc**2 + rbm**2)
   loc = minloc(rel_residual)  
   lib_dist = rel_residual(loc(1), loc(2)) 
   
   ! do the alternate retrieval  
   re = radii(loc(2))
   ! don't want any 2's and such 
   if (re < 4. ) re = 4.
   ! the library taus don't have 0., so we need to explicitly count 
   !  for it. But we can't put in a 0. straight because
   ! L3 code can't handle an explicit 0. value. 
   if (loc(1) == 1) then 
     tau = 0.01
   else
     tau = library_taus(loc(1)-1)
   endif
  
   deallocate (rel_residual,reflibbl)
  
 end subroutine asl_interior
  
  
 subroutine asl_boundary(Ram, Rbm_in, bands, radii, tau, re, lib_dist, Pc, & 
     theta, theta0, phi, phase_liquid, Ram_corr)
  
   !  For solutions outside the LUT - derive one and 
   !  essentialy an uncertainty in the form of a distance - This is for 
   !alternate solution exterior, just the boundary (for water, re=4, 
   !  re=30 and tau = 158.8)
   !  for ice re=5,re=60 and tau = 158.8 (top, bottom and rightmost)
   !
   !  Ram  IN  non-abs reflectance
   !  Rbm_in  IN  absorbing band reflectance
   !  bands  IN  size 2 indicies of non-abs and absorbing indicies
   !  radii  IN  array of table re values
   !  tau  IN/OUT  retrieved optical thickness
   !  re  IN/OUT  retrieved radius
   !  lib_dist  a distance outside the LUT, a measure of amount off
   !  Pc  IN  cloud top pressure
   !  theta  IN  I believe the sat zen angle
   ! theta0  IN  I believe the sol zen angle
   ! phi   IN  I believe the delta azimuth angle
   ! phase_liquid  IN  switch for liquid phase if true
   !  Ram_corr  OUT  may be abs band refl
  
   use generalauxtype
   use libraryarrays
   use libraryinterpolates
   use nonscience_parameters
   use mod06_run_settings
   use science_parameters
  
   real, intent(in) :: Ram, Rbm_in, Pc, theta, theta0, phi
   real, dimension(:), intent(in) :: radii
   integer, dimension(:), intent(in) :: bands
   real, intent(inout) :: tau, re, lib_dist, Ram_corr
   logical, intent(in) :: phase_liquid
   
   real, dimension(:,:), allocatable :: temp_refl,reflibBL
   real,dimension(:),allocatable::rel_residual
   integer :: size_tau, size_re, loc, rel_loc(1),ire_ct_1, &
     temp_array_size,s_rt,iL,iR
   logical :: re_altered
   
   real :: Ramc, crefl, Rbm
   integer :: iter, max_iter,iiiloc
  
   real :: albedoA, albedoB,normConst
   real :: sfr, fti1, fti0
  
   size_tau = number_taus + 1
   size_re = size(refliba, 2)
  
   ramc = ram
   rbm = rbm_in
  
   allocate(reflibbl(size_tau, size_re))
   
   reflibbl(:,:) = reflibb(:,:)
     
   if( bands(2) == band_0370)then 
     normconst = pi/(cos(theta0*d2r)*solar_constant_37)
     rbm = rbm_in * normconst
     reflibbl(:,:) = rad37lib(:,:)
   endif
   
   !  set the top contour (CT 1) for ice =1 (re =5) for water = 2(re = 4)
   !  can do it in mod06_run_settings
  
   ire_ct_1 = top_nk_asl_contour_ice 
   if(phase_liquid)ire_ct_1 = top_nk_asl_contour_water
   
   !  reflibA is the modified copied non-abs reflectance table made back in 
   !  vis non absorbing work and reflibBL is the abs table
   albedoa = refliba(1,ire_ct_1)
   albedob = reflibbl(1,ire_ct_1)
   
   s_rt = size_tau+size_re
   temp_array_size = size_tau + s_rt - (ire_ct_1 + 1)
   
   allocate(rel_residual(temp_array_size))
   allocate(temp_refl(2,temp_array_size))
   ! EXPLAIN YOUR WORK!  
   temp_refl(1,1:size_tau) = refliba(1:size_tau,ire_ct_1)
   
   temp_refl(1, size_tau+1 : s_rt-ire_ct_1) = &
     refliba(size_tau,ire_ct_1+1:size_re)  !
   
   temp_refl(1, s_rt-ire_ct_1+1 : temp_array_size) = &
     refliba(1:size_tau-1,size_re)
  
   !   assign the radiance+emission of the 3.7 LUT
   temp_refl(2, 1 : size_tau) = reflibbl(1:size_tau,ire_ct_1)
   temp_refl(2, size_tau+1 : s_rt-ire_ct_1) = &
     reflibbl(size_tau,ire_ct_1+1:size_re)  !
   temp_refl(2, s_rt-ire_ct_1+1 : temp_array_size) = &
     reflibbl(1:size_tau-1,size_re)
   
   !  abs refl is below sfc albedo: region I
   ! the point has no chance, it is darker than the surface albedo
   if (ramc < albedoa  .OR.  & 
     ( rbm < minval(reflibbl(size_tau,ire_ct_1:size_re)) .and. &
     ramc > minval(refliba(size_tau,ire_ct_1:size_re)) ) .OR. & 
     ( rbm > maxval(reflibbl(size_tau,ire_ct_1:size_re)) .and. &
     ramc > maxval(refliba(size_tau,ire_ct_1:size_re)) )) then
  
     tau = fillvalue_real
     re  = fillvalue_real
     lib_dist = max_cost_function
     ram_corr = ramc
     deallocate (rel_residual, temp_refl, reflibbl)
     return
     
   endif
   ! collect the value of cost function
   rel_residual = sqrt( (temp_refl(1,:) - ramc)**2 + &
     (temp_refl(2,:) -rbm)**2 ) / sqrt(ramc**2 + rbm**2)
   rel_loc = minloc(rel_residual)  
   lib_dist = rel_residual(rel_loc(1))
   
   !  reclaculate lib_dist with both x and y in terms of reflectances, 
   !  if band(2)=3.7
  
   ! do the alternate retrieval  
   if(rel_loc(1) <= size_tau)then
     ! contour 1
     re = radii(ire_ct_1)
     loc = rel_loc(1)
        
     if(bands(1) == band_0065)then
       if(phase_liquid) then       
         ramc = rayleigh_liq(ire_ct_1)
       else
         ramc = rayleigh_ice(ire_ct_1)
       endif
             
       CALL calcdistanceandminloc(ramc,rbm,temp_refl(1,1:size_tau), &
         temp_refl(2,1:size_tau),lib_dist,loc)
       ram_corr = ramc
     endif
              
     if (loc == 1) then 
       tau = 0.01
     else
       tau = library_taus(loc-1)
     endif
       
     if(ramc > maxval(refliba(size_tau,ire_ct_1:size_re)))then 
       if(bands(1) == band_0163  .AND. bands(2) == band_0213 )then
         tau = fillvalue_real
       else
         if((rbm - reflibbl(size_tau,ire_ct_1)) < 0.0)then
           do iiiloc = ire_ct_1+1,size_re
             il = iiiloc-1
             ir = iiiloc
             if((rbm - reflibbl(size_tau,iiiloc)) >=0.0)exit
           enddo     
           tau = max_tau_rtrieved
           re=radii(il)
  
           if(abs(reflibbl(size_tau,ir)-reflibbl(size_tau,il)) > 0.000001)then
             re = (rbm-reflibbl(size_tau,il))*(radii(ir) - &
               radii(il))/(reflibbl(size_tau,ir)-reflibbl(size_tau,il)) + &
               radii(il)
             if(re > radii(ir) .OR. re < radii(il))re = radii(il)
           endif
         endif
       endif         
     endif
       
   elseif(rel_loc(1) >  size_tau .and. rel_loc(1) <= s_rt - ire_ct_1)then
     !contour 4
     loc = rel_loc(1) - size_tau + 1
     re = radii(loc)
  
     if(loc==1 .OR. &
       ((rbm - reflibbl(size_tau,loc)) < 0.0 .AND. loc <= size_re))then
       do iiiloc = loc , size_re
         if((rbm - reflibbl(size_tau,iiiloc)) >=0.0)exit
         loc = iiiloc
       enddo
           
       il = loc
       ir = loc+1
     elseif((rbm - reflibbl(size_tau,loc)) >= 0.0 .AND. loc <= size_re)then    
       do iiiloc = loc-1 , 1,-1
         if((rbm - reflibbl(size_tau,iiiloc)) <= 0.0)exit
         loc = iiiloc
       enddo
       il = loc-1
       ir = loc
     endif
        
     if(ir <= size_re .and. il >= 1)then
       re = radii(loc)
       if(abs(reflibbl(size_tau,ir)-reflibbl(size_tau,il)) > 0.000001 &
         .and. abs(rbm - reflibbl(size_tau,loc)) > 0.0)then
         re = radii(loc) + (rbm-reflibbl(size_tau,loc))*(radii(ir) - &
           radii(il))/(reflibbl(size_tau,ir)-reflibbl(size_tau,il))  
         if(re > radii(ir) .OR. re < radii(il))re = radii(loc)            
       endif
     endif
         
     if(bands(1) == band_0065)then
       if(phase_liquid)then       
         ramc = rayleigh_liq(loc)
       else
         ramc = rayleigh_ice(loc)
       endif
       ram_corr = ramc
     endif
         
     tau = fillvalue_real
     if(bands(1) == band_0163  .AND. bands(2) == band_0213 )then
       tau = fillvalue_real
     else
       tau = max_tau_rtrieved
     endif
         
     if(ramc > 0 .and. ramc < maxval(refliba(1:size_tau,size_re)) .and. &
       ir > size_re)THEN
       tau = fillvalue_real
       re = fillvalue_real
       lib_dist = max_cost_function          
     endif
   else
     !contour 2
     re = radii(size_re)
     loc = rel_loc(1) - (s_rt-ire_ct_1) + 1
  
     if(bands(1) == band_0065)then
       if(phase_liquid)then       
         ramc = rayleigh_liq(size_re)
       else
         ramc = rayleigh_ice(size_re)
       endif
             
       CALL calcdistanceandminloc(ramc,rbm, &
         temp_refl(1,s_rt-ire_ct_1+1: temp_array_size), &
         temp_refl(2,s_rt-ire_ct_1+1: temp_array_size),lib_dist,loc)
         ram_corr = ramc            
     endif
       
     if ( loc == 1 ) then 
       tau = 0.01
     else
       tau = library_taus(loc-1 )
     endif
     
     if(ramc > maxval(refliba(1:size_tau,size_re)))then
       tau =  fillvalue_real
       re = fillvalue_real
       lib_dist = max_cost_function
  
     endif
               
   endif     
     
   if(tau > library_taus(size_tau-1)-1.0 )then
     tau =  fillvalue_real
     re = fillvalue_real
     lib_dist = max_cost_function
   endif
     
   ! don't want any 2's and such 
   ! if (re < 4. ) re = 4.
   ! the library taus don't have 0., so we need to explicitly count 
   !  for it. But we can't put in a 0. straight because
   ! L3 code can't handle an explicit 0. value. 
  
   deallocate (rel_residual,temp_refl,reflibbl)
  
 end subroutine asl_boundary
  
 subroutine calcdistanceandminloc(R1,R2,R1vec, R2vec,mindist,loc_index)
   !
   !  R1
   !  R2
   !  R1vec
   !  R2vec
   !  mindist
   !  loc_index
  
   REAL,INTENT(IN)::R1,R2
   REAL,DIMENSION(:),INTENT(IN)::R1vec,R2vec
   REAL,INTENT(OUT)::mindist
   INTEGER,INTENT(OUT)::loc_index
  
   real,dimension(:),allocatable::rel_residual
   INTEGER::rel_loc(1)
  
   allocate(rel_residual(size(r1vec)))
  
   rel_residual = sqrt( (r1vec - r1)**2 + (r2vec -r2)**2 ) / &
     sqrt(r1**2 + r2**2)
   rel_loc = minloc(rel_residual)  
   mindist = rel_residual(rel_loc(1))
   loc_index = rel_loc(1)
     
   deallocate(rel_residual)
     
   return
 end subroutine calcdistanceandminloc
  
 subroutine solve_retrieval_noswir(optical_thickness_vector, &
     library_radii, effective_radius, optical_thickness) 
   !
   !  optical_thickness_vector
   !  library_radii
   !  effective_radius
   !  optical_thickness
     
   use nonscience_parameters
   use science_parameters
   use libraryarrays, only:library_taus                 
  
   real, intent(in)   :: optical_thickness_vector(:)
   real, intent(in) :: library_radii(:)
   real, intent(in) :: effective_radius
   real, intent(out)  :: optical_thickness
  
  
   integer :: ilo, ihi
  
   ! we have no SWIR (or knocked it out on purpose) 
   ! we are fixing the value of effective radius on some value that 
   !  is generally agreed on 
   ! as being reasonable. We will return an optical thickness retrieval 
   !  that's all nice
   ! and interpolated, but the re is going to be a constant
   ! There will be no ASL retrieval because we are just staying on the 
   !  re=const curve in the library
  
   call findinterpolationrange(3, effective_radius, size(library_radii), &
     real(library_radii), ilo,  ihi)
  
   optical_thickness = lagrangeinterp(effective_radius, &
     real(library_radii(ilo:ilo+2)), optical_thickness_vector(ilo:ilo+2) )
                                       
   if (optical_thickness > 150.) optical_thickness = 150.
   if (optical_thickness < 0.) then 
     optical_thickness = fillvalue_real
   endif
  
 end subroutine solve_retrieval_noswir
  
 subroutine solveretrieval( residual, optical_thickness_vector, &
   library_radii, effective_radius, optical_thickness, debug, &
   use_nearest, quality_out)
   !
   !  residual IN  vector of fractional diff of rho(abs) measured - 
   !      rho(abs) of the rho library for the set of solved tau (by the 
   !      non-abs algorithm) at each re (make sense?)
   !  optical_thickness_vector IN  tau solved for at each re
   !  library_radii  IN  the re axis values
   !  effective_radius  OUT  solved re
   !  optical_thickness  OUT  solved tau
   !  debug  IN  switch for debugging
   !  use_nearest  OUT  apparently, this failed and solveretrieval_nearest needs
   !       to be called
   !  quality_out  OUT  many states of the re retrieval
  
   use generalauxtype
   use nonscience_parameters
   use science_parameters
   use libraryarrays, only:library_taus                    
  
   real, intent(in)   :: residual(:), optical_thickness_vector(:)
   real(single), intent(in) :: library_radii(:)
   logical, intent(in)  :: debug
   logical, intent(inout) :: use_nearest
   real, intent(out)  :: effective_radius, optical_thickness
   integer,intent(out)::quality_out
  
   real :: max_allowable_tau,local_max_allowable           !!added
   integer:: total_taus                                    !!added
  
   integer              :: number_local_minima, ilo, ihi, zero_count
   integer  :: local_minima(8)
  
   INTEGER, parameter   :: nre=18
   INTEGER :: crossing_num, crossing_vector(nre), local_max_num, &
     local_max_vector(nre), local_min_num, local_min_vector(nre), k
   integer(integer_onebyte) :: quality
  
   integer:: res_size
  
   real  re_water_outofbounds_high, re_water_outofbounds_low
   logical :: re_altered, correction_made
  
   local_minima(:) = 0
   re_water_outofbounds_high = 30.; re_water_outofbounds_low = 4.
   
   correction_made = .false.
   
   ! get effective radius solution                        !added/changed
   ! changed for water residual going in from re = 4
   !  This decides if we are doing ice or water retr by checking the 
   !  library array size  The find_zero_crossings, find_local_minima,
   !  find_local_maxima gets the indicies of the tops, bottoms and 
   !  (value before) the zero crossing of the residual and the # of these
   !  Then, solution_re
   if(is_water_phase(library_radii))then
     call find_zero_crossings (residual(top_nk_asl_contour_water:), &
       crossing_vector,  crossing_num)
     call find_local_minima   (residual(top_nk_asl_contour_water:), &
       local_min_vector, local_min_num)
     call find_local_maxima   (residual(top_nk_asl_contour_water:), &
       local_max_vector, local_max_num)
  
     call solution_re (library_radii(top_nk_asl_contour_water:), &
       residual(top_nk_asl_contour_water:) , crossing_vector, crossing_num, &
       local_min_vector, local_min_num, local_max_vector, local_max_num, &
       effective_radius, quality )
   else
     call find_zero_crossings (residual, crossing_vector,  crossing_num)
     call find_local_minima   (residual, local_min_vector, local_min_num)
     call find_local_maxima   (residual, local_max_vector, local_max_num)
  
     call solution_re (library_radii, residual, crossing_vector, crossing_num, &
       local_min_vector, local_min_num, local_max_vector, local_max_num, &
                         effective_radius, quality )
   endif
  
   quality_out = quality   !added/changed
     
   use_nearest = .false.
   if ( quality == -2 .or. quality == -3 .or. quality == -4 &
     .OR. quality == -6) then    !added/changed
  
     use_nearest = .true.
     RETURN
  
   endif
   ! get optical thickness solution
  
   ! PARTIAL RETRIEVAL LOGIC ADDENDUM FOR RETENTION OF TAU SOLUTION 
   !  FOR FAILED RE
   !  order is important here--the following block needs to be done 
   !  here before the call to findinterpolationrange
   ! for performing a partial tau retrieval after having retrieved the 
   !  maximum or minimum desireable, or failed to retrieve liquid library re
  
   ! sep: 3=8-06 
   !  For water clouds there are 4 possibilities at this point: re>=30, 
   !  re<=4, re=INVALID, or re valid
   !  For ice clouds there are 2 possibilities at this point: re=INVALID 
   !  or re valid
   !  It was decided to perform "partial" optical thickness retrievals 
   !  for any non-valid case by temporarily assigning re a 
   !  value of 10 or 30 microns for ice and water, respectively. 
   !  The logic to account for this is simpler than the previous version.
  
   re_altered = .false.
   if (is_water_phase(library_radii) .and. &  ! changed/added <= to <
     ( (effective_radius < re_water_outofbounds_low) .or. & 
     (effective_radius == invalid_attempted_but_failed_re) ) ) then!{
     effective_radius = 10.
     re_altered = .true.
   else!}{
   ! for performing a partial tau retrieval after having failed to 
   !  retrieve ice library re
     if (effective_radius == invalid_attempted_but_failed_re) then!{
       effective_radius = 30.
       re_altered = .true.
     endif!}
   endif!}
   ! END PARTIAL RETRIEVAL LOGIC ADDENDUM FOR RETENTION OF TAU 
   !   SOLUTION FOR FAILED RE 
  
   call findinterpolationrange(3, effective_radius, size(library_radii), &
     real(library_radii), ilo, ihi)
                           
   !!! added !!!  Cloud folks emphasis
   !!! correction starts here !!!!!!!!
   total_taus = size(library_taus)
   max_allowable_tau = library_taus(total_taus)
   local_max_allowable = max_allowable_tau - 1.0
  
   !  we go in if tau_vector(ilo:ilo+2) has a 158.78 
   !  or -1 (i.e, lagrange interpolation points)
   if(.not.(use_nearest) .and. &
     (maxval(optical_thickness_vector(ilo:ilo+2)) > local_max_allowable  &
     .OR. minval(optical_thickness_vector(ilo:ilo+2)) < 0))then
  
     correction_made = .true.                                    
     if(ilo <= 2)then
       if(maxval(optical_thickness_vector(1:2)) > local_max_allowable  & 
         .OR. minval(optical_thickness_vector(1:2)) < 0)then 
         optical_thickness = fillvalue_real
         use_nearest = .true.
         quality_out = -4
       else
         ilo =1
         optical_thickness = optical_thickness_vector(ilo) + &
           ( ( effective_radius - real(library_radii(ilo))) /  &
           (real(library_radii(ilo+1)-library_radii(ilo))) ) * &
           (optical_thickness_vector(ilo+1) - optical_thickness_vector(ilo) )  
       endif
  
     elseif(maxval(optical_thickness_vector(ilo-1:ilo)) > local_max_allowable  &
       .OR. minval(optical_thickness_vector(ilo-1:ilo)) < 0)then
       optical_thickness = fillvalue_real
       use_nearest = .true.
       quality_out = -4
             
     elseif(maxval(optical_thickness_vector(ilo-1:ilo)) < local_max_allowable & 
       .AND. minval(optical_thickness_vector(ilo-1:ilo)) > 0)then
       ilo = ilo-2
       if(optical_thickness_vector(ilo+3) < local_max_allowable & 
         .AND. optical_thickness_vector(ilo+3) > 0)ilo = ilo+1
  
       optical_thickness = optical_thickness_vector(ilo+1) + &
         ( ( effective_radius - real(library_radii(ilo+1))) /  &
         (real(library_radii(ilo+2)-library_radii(ilo+1))) ) * &
         (optical_thickness_vector(ilo+2) - optical_thickness_vector(ilo+1) )   
     else
       !  this else can be removed. If all above fails it should come here. 
       !  Tested 5 granuels and didn't come here
       !  but to be safe make it filled
       !  optical_thickness = fillvalue_real
       use_nearest = .true.
       quality_out = -4
     endif    
   endif
  
   !!! correction done !!!! The end of their correction
  
   if(.NOT.(correction_made))then 
     optical_thickness = lagrangeinterp(effective_radius, &
       real(library_radii(ilo:ilo+2)), &
       optical_thickness_vector(ilo:ilo+2) )
   endif
  
   if( optical_thickness < optical_thickness_vector(ilo) .AND. &
     (.NOT.(correction_made))) then!{    !!!!!!added /changed
  
     call findinterpolationrange(2, effective_radius, size(library_radii), &
       real(library_radii), ilo,  ihi)
       optical_thickness = optical_thickness_vector(ilo) + &
       ( ( effective_radius - real(library_radii(ilo))) / &
       (real(library_radii(ilo+1)-library_radii(ilo))) ) * &
       (optical_thickness_vector(ilo+1) - optical_thickness_vector(ilo) ) 
   endif!}
   
   ! KILL RE FOR PARTIAL RETRIEVAL
   ! Now that the optical thickness has been determined, set the 
   !  effective radius to fill if re_altered = .TRUE.
   if (re_altered) effective_radius = fillvalue_real
  
   ! RANGE CHECKING of optical thickness and effective radius
  
   ! OPTICAL THICKNESS 
  
   ! tau low boundary of range
   if (optical_thickness /= fillvalue_real .and. optical_thickness < 0.) then 
     optical_thickness = fillvalue_real
     effective_radius = fillvalue_real
     use_nearest = .true.
     quality_out = -4
   endif 
   
   if (optical_thickness /= fillvalue_real .and. &
     optical_thickness < 0.01 .and. optical_thickness > 0.) then!{ 
     optical_thickness = 0.01
   endif!}
   ! end tau low boundary of range
  
   ! tau high boundary of range
   !  QA for large optical thickness
   !   if (optical_thickness > 150.) then!{
   !     optical_thickness_outofbounds = 2
   !   elseif (optical_thickness > 100. .and. optical_thickness <= 150.) then!}{
   !     optical_thickness_outofbounds = 1
   !   else!}{
   ! optical_thickness_outofbounds = 0
   !   endif!}
  
   !   if ( optical_thickness > 100.) then!{
   !     optical_thickness = 100.
   !   endif!}
   ! end tau high boundary of range
  
   ! EFFECTIVE RADIUS 
  
   ! NOTE: Ideally, it would be more satisfying to assign re retrievals 
   !  outside of the library space to fill value,
   !  even here after other checks have done the same however, retrievals 
   !  can be within epsilon of boundary 
   !  (e.g., Liam's pixel with ice re = 90.00000763).  That is the intent 
   !  of this strict range enforcement here.
  
   ! re low boundary of range
   if (effective_radius /= fillvalue_real .and. &
     effective_radius < library_radii(1)) then!{ 
     effective_radius = library_radii(1)
   endif!}
   ! end re low boundary of range
  
   ! re high boundary of range
   if (effective_radius > library_radii(size(library_radii))) then!{ 
      effective_radius = library_radii(size(library_radii))
   endif!}
   ! end re high boundary of range
  
 end subroutine solveretrieval
  
 subroutine linear_interpolate_for_root(residual, rad, re)
   !
   !  residual
   !  rad
   !  re
   !
   implicit none
  
   real, intent(in)     :: residual(2), rad(2)
   real, intent(out)    :: re
  
   re = -1.0*( rad(1)*residual(2)/(residual(1)-residual(2)) + &
     rad(2)*residual(1)/(residual(2)-residual(1)) )
  
 end subroutine linear_interpolate_for_root
  
 subroutine findinterpolationrange( n1, xx, n,x, k1, k2)
   !
   !  n1
   !  xx
   !  n
   !  x
   !  k1
   !  k2
   !
   implicit none
  
   integer, intent(in)   :: n1, n
   real, intent(in)      :: xx, x(n)
  
   integer,intent(out)   :: k1, k2
   
   integer               :: i, n2
  
   if (n > n1) then!{
     n2 = 0.5*n1
     i = 1
     do while (xx > x(i) .and. i < n)
       i = i + 1
     enddo
     k1 = max(1,i-n2)
     k2 = k1 + n1 - 1
     if (k2 > n) then!{
       k2 = n
       k1 = k2 - n1 + 1
     end if
   else!}{
   !    if # of exist data set n <= n1, i
   !    use n1 = n for interp_mas_geoolation.
     k1 = 1
     k2 = n
   end if
  
 end subroutine findinterpolationrange
  
 real function lagrangeinterp( xx, x, y)
   ! in the previous interation of the COT retrieval code this was 
   !  known as 'BINTPB'
   !
   !  xx
   !  x
   !  y
   !
   implicit none
   real, intent(in)  :: x(3), y(3), xx 
   lagrangeinterp = (xx-x(2))* (xx-x(3))/ (x(1)-x(2))/ (x(1)-x(3))*y(1) + &
                    (xx-x(3))* (xx-x(1))/ (x(2)-x(3))/ (x(2)-x(1))*y(2) + &
                    (xx-x(1))* (xx-x(2))/ (x(3)-x(1))/ (x(3)-x(2))*y(3)
 end function lagrangeinterp
  
   ! sssssssssssssssssssssssssssssssssssssssssssssss
 subroutine check_for_signchange (x, signchange)
   ! sssssssssssssssssssssssssssssssssssssssssssssss
   !
   !  x
   !  signchange
   !
   !  *** Copied from collect 5 core_science code
  
   ! this routine reports on sign change for a set of 3 values
   ! signchnage == 0 -> no sign change
   ! signchange == 1 -> sign changes between x1 and x2
   ! signchange == 2 -> sign changes bweteen x2 and x3
  
   implicit none
  
   real, intent(in)      :: x(3)
   integer, intent(out)  :: signchange
  
   if ( (x(1)==0 .and. x(2)==0) .or. (x(1)==0 .and. x(3)==0) .or. &
     (x(2)==0 .and. x(3)==0) ) Then
     signchange = 0
   else!}{
     if ( sign(1.,x(1))*sign(1.,x(2))*sign(1.,x(3)) /= -1.) Then
       signchange =1
       if ( x(1) > 0 .AND. x(2) > 0 .AND. x(3) > 0) then!{
         signchange =0
       end if
  
       if ( x(1) < 0 .AND. x(2) < 0 .AND. x(3) < 0) then!{
         signchange =0
       end if
  
     else!}{
       signchange =1
     endif!}
   endif!}
  
   if (signchange == 1) then!{
     if ( sign(1.,x(1))*sign(1.,x(2)) == -1. ) signchange = 1
     if ( sign(1.,x(2))*sign(1.,x(3)) == -1. ) signchange = 2
   endif!}
  
 end subroutine check_for_signchange
  
 subroutine quad_interpolate_for_root ( radii, residual, radius_solution, &
     status)
   !
   !  radii
   !   residual
   !  radius_solution
   !  status
   !
   use generalauxtype
   implicit none
  
   real(single), intent(in)  :: radii(:), residual(:)
  
   real(single), intent(out) :: radius_solution
   integer,      intent(out) :: status
  
   real(double) :: d1, d2, d3, a0, a1, a2, root1, root2, z
  
   d1 = residual(1) / ((radii(1)-radii(2))*(radii(1)-radii(3)))
   d2 = residual(2) / ((radii(2)-radii(1))*(radii(2)-radii(3)))
   d3 = residual(3) / ((radii(3)-radii(1))*(radii(3)-radii(2)))
  
   a0 = d1*radii(2)*radii(3) + d2*radii(1)*radii(3) + d3*radii(1)*radii(2)
   a1 = -d1*(radii(2)+radii(3)) - d2*(radii(1)+radii(3)) - d3*(radii(1)+radii(2))
   a2 = d1 + d2 + d3
  
  
   z = a1**2 - 4*a0*a2
  
   if(z >= 0. .and. a2 /= 0.) then!{
     root1 = (-a1 + sqrt(z))/(2*a2)
     root2 = (-a1 - sqrt(z))/(2*a2)
   else!}{
     root1 = -999; root2 = -999
     status = 1
     return
   end if
   if (root1 < radii(1) .and. root2 < radii(1)) then!{
     radius_solution =-999.
     status = 2
   elseif (root1 > radii(3) .and. root2 > radii(3)) then!}{
     radius_solution =-999.
     status = 3
   elseif (root1 >= radii(1) .and. root1 <= radii(3)) then!}{
     radius_solution = root1
     status = 4
   elseif(root2 >= radii(1) .and. root2 <= radii(3)) then!}{
     radius_solution = root2
     status = 5
   else!}{
     radius_solution = -999.
     status = 1
   endif!}
  
 end subroutine quad_interpolate_for_root
  
 end module retrieval_solution_logic

planck_functions

Definition: planck_functions.f90:1

retrieval_solution_logic::findinterpolationrange

subroutine findinterpolationrange(n1, xx, n, x, k1, k2)

Definition: retrieval_solution_logic.f90:1636

core_arrays

Definition: core_arrays.f90:1

core_arrays::thermal_correction_oneway

real, dimension(set_number_of_bands) thermal_correction_oneway

Definition: core_arrays.f90:133

libraryarrays::rad37lib

real, dimension(:,:), allocatable rad37lib

Definition: libraryarrays.f90:77

retrieval_solution_logic::solveretrieval

subroutine solveretrieval(residual, optical_thickness_vector, library_radii, effective_radius, optical_thickness, debug, use_nearest, quality_out)

Definition: retrieval_solution_logic.f90:1361

core_arrays::cloud_top_pressure

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

Definition: core_arrays.f90:156

libraryarrays::flux_up_water_solar

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

Definition: libraryarrays.f90:59

libraryarrays

Definition: libraryarrays.f90:1

modis_numerical_module::linearinterpolation

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

Definition: modis_numerical_module.f90:143

retrieval_solution_logic::solve_retrieval_noswir

subroutine solve_retrieval_noswir(optical_thickness_vector, library_radii, effective_radius, optical_thickness)

Definition: retrieval_solution_logic.f90:1318

spline

subroutine spline(s, x, y, n, in, t, il, iu, vl, vu, e, u)

Definition: phs.f:1348

retrieval_solution_logic

Definition: retrieval_solution_logic.f90:1

retrieval_solution_logic::check_for_signchange

subroutine check_for_signchange(x, signchange)

Definition: retrieval_solution_logic.f90:1691

retrieval_solution_logic::find_local_maxima

subroutine find_local_maxima(residual, local_max_vector, local_max_num)

Definition: retrieval_solution_logic.f90:97

libraryinterpolates::int_fluxdownwater_solar

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

Definition: libraryinterpolates.f90:20

sign

#define sign(x)

Definition: misc.h:95

mod06_run_settings::band_0370

integer, parameter band_0370

Definition: mod06_run_settings.f90:26

real

#define real

Definition: DbAlgOcean.cpp:26

core_arrays::thermal_correction_twoway

real, dimension(set_number_of_bands) thermal_correction_twoway

Definition: core_arrays.f90:132

retrieval_solution_logic::find_local_minima

subroutine find_local_minima(residual, local_min_vector, local_min_num)

Definition: retrieval_solution_logic.f90:61

science_parameters

Definition: science_parameters.f90:1

mod06_run_settings::band_0086

integer, parameter band_0086

Definition: mod06_run_settings.f90:26

libraryinterpolates::int_cloud_emissivity_ice

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

Definition: libraryinterpolates.f90:24

mod06_run_settings::band_0213

integer, parameter band_0213

Definition: mod06_run_settings.f90:26

nonscience_parameters::max_tau_rtrieved

real, parameter max_tau_rtrieved

Definition: nonscience_parameters.f90:20

retrieval_solution_logic::calcdistanceandminloc

subroutine calcdistanceandminloc(R1, R2, R1vec, R2vec, mindist, loc_index)

Definition: retrieval_solution_logic.f90:1287

libraryinterpolates::int_cloud_emissivity_water

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

Definition: libraryinterpolates.f90:23

science_parameters::cox_munk

logical cox_munk

Definition: science_parameters.f90:11

libraryarrays::library_fluxsolarzenith

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

Definition: libraryarrays.f90:43

libraryinterpolates

Definition: libraryinterpolates.f90:1

libraryarrays::library_taus

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

Definition: libraryarrays.f90:40

libraryarrays::refliba

real, dimension(:,:), allocatable refliba

Definition: libraryarrays.f90:77

generalauxtype

Definition: GeneralAuxType.f90:1

retrieval_solution_logic::find_zero_crossings

subroutine find_zero_crossings(residual, crossing_vector, crossing_num)

Definition: retrieval_solution_logic.f90:12

libraryarrays::spherical_albedo_water

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

Definition: libraryarrays.f90:58

core_arrays::solar_zenith_angle

real, dimension(:,:), allocatable solar_zenith_angle

Definition: core_arrays.f90:6

modis_numerical_module

Definition: modis_numerical_module.f90:1

libraryinterpolates::int_surface_emissivity_ice

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

Definition: libraryinterpolates.f90:24

mod06_run_settings

Definition: mod06_run_settings.f90:1

retrieval_solution_logic::solve_for_zero_crossing

subroutine solve_for_zero_crossing(re, residual, residual_crossing_index, fillval_r4, local_min_vector, local_min_num, local_max_vector, local_max_num, retrieval, quality)

Definition: retrieval_solution_logic.f90:320

libraryinterpolates::int_fluxdownice_solar

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

Definition: libraryinterpolates.f90:21

libraryinterpolates::int_surface_emissivity_water

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

Definition: libraryinterpolates.f90:23

max

#define max(A, B)

Definition: main_biosmap.c:61

retrieval_solution_logic::is_water_phase

logical function is_water_phase(library_radii)

Definition: retrieval_solution_logic.f90:520

libraryarrays::spherical_albedo_ice

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

Definition: libraryarrays.f90:67

endif

endif() set(LIBS $

Definition: CMakeLists.txt:6

retrieval_solution_logic::solution_re

subroutine solution_re(re, residual, crossing_vector, crossing_num, local_min_vector, local_min_num, local_max_vector, local_max_num, retrieval, quality)

Definition: retrieval_solution_logic.f90:136

planck_functions::modis_planck

real function, public modis_planck(platform_name, TEMP, BAND, UNITS)

Definition: planck_functions.f90:510

libraryarrays::rayleigh_liq

real, dimension(:), allocatable rayleigh_liq

Definition: libraryarrays.f90:78

libraryarrays::ice_radii

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

Definition: libraryarrays.f90:63

retrieval_solution_logic::top_nk_asl_contour_ice

integer, parameter top_nk_asl_contour_ice

Definition: retrieval_solution_logic.f90:5

libraryarrays::flux_up_water_sensor

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

Definition: libraryarrays.f90:60

libraryarrays::reflibb

real, dimension(:,:), allocatable reflibb

Definition: libraryarrays.f90:77

generalauxtype::real_s_equal

logical function real_s_equal(x, y)

Definition: GeneralAuxType.f90:64

retrieval_solution_logic::ray_corr_nearest

subroutine ray_corr_nearest(refl_source, As, iw, tau, re, Pc, sfr, fti1, fti0, fluxup_solar, fluxup_sensor, theta, theta0, phi, location, crefl)

Definition: retrieval_solution_logic.f90:536

mod06_run_settings::band_0163

integer, parameter band_0163

Definition: mod06_run_settings.f90:26

retrieval_solution_logic::calc37radiancelibcm

subroutine calc37radiancelibcm(intensity, intensity_g, thermal_trans_1way, thermal_trans_2way, solar_zenith, cl_emis, sf_emis)

Definition: retrieval_solution_logic.f90:795

nonscience_parameters

Definition: nonscience_parameters.f90:4

retrieval_solution_logic::calc37radianceliblamb

subroutine calc37radianceliblamb(intensity, intensity_g, thermal_trans_1way, thermal_trans_2way, solar_zenith, sfr, fti1, fti0, fri1)

Definition: retrieval_solution_logic.f90:749

retrieval_solution_logic::asl_boundary

subroutine asl_boundary(Ram, Rbm_in, bands, radii, tau, re, lib_dist, Pc, theta, theta0, phi, phase_liquid, Ram_corr)

Definition: retrieval_solution_logic.f90:1016

nonscience_parameters::fillvalue_real

real, parameter fillvalue_real

Definition: nonscience_parameters.f90:13

#define pi

Definition: vincenty.c:23

spline_module

Definition: spline_module.f90:1

core_arrays::sensor_zenith_angle

real, dimension(:,:), allocatable sensor_zenith_angle

Definition: core_arrays.f90:113

science_parameters::d2r

real, parameter d2r

Definition: science_parameters.f90:60

core_arrays::relative_azimuth_angle

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

Definition: core_arrays.f90:118

libraryinterpolates::int_fluxupice_sensor

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

Definition: libraryinterpolates.f90:21

core_arrays::surface_temperature

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

Definition: core_arrays.f90:152

single

subroutine single

Definition: single.f:2

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::flux_up_ice_solar

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

Definition: libraryarrays.f90:68

retrieval_solution_logic::solveretrieval_nearest

subroutine solveretrieval_nearest(xx_pt, yy_pt, Ram, Rbm, twobands, radii, tau, re, lib_dist, phase_liquid, Ram_corr, quality_in, CH37_IDX, CTopT, CH37_NUM, platFormName)

Definition: retrieval_solution_logic.f90:634

retrieval_solution_logic::linear_interpolate_for_root

subroutine linear_interpolate_for_root(residual, rad, re)

Definition: retrieval_solution_logic.f90:1620

nonscience_parameters::max_cost_function

real, parameter max_cost_function

Definition: nonscience_parameters.f90:21

libraryinterpolates::int_fluxdownwater_sensor

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

Definition: libraryinterpolates.f90:20

retrieval_solution_logic::quad_interpolate_for_root

subroutine quad_interpolate_for_root(radii, residual, radius_solution, status)