NASA Ocean Color

ocssw V2022

 module atmospheric_correction_module
  
     use specific_ancillary
  
 implicit none
  
 private
  
 public :: atmospheric_correction
  
 !... scale in p, tpw and miu dimension:
  
 INTEGER, PARAMETER :: NumberOfPressureHeights = 10
 REAL, DIMENSION(NumberOfPressureHeights), PARAMETER :: &
 heightScale = (/ 100., 200., 300., 400., 500., 600., 700., 800., 900., 1000. /)
  
  
 ! changed by G. Wind 5.6.04 to match the new axis of transmittance table
 INTEGER, PARAMETER :: NumberOfPrecipitableWater = 53
 REAL, DIMENSION(NumberOfPrecipitableWater), PARAMETER :: &
  pwScale = (/ 0.0, 0.02, 0.04, 0.06, 0.08, 0.1, 0.12, 0.14, 0.16, 0.18, 0.2, &
                     0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6,1.8, &
                    2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, &
                    4.2, 4.4, 4.6, 4.8, 5.0, 5.2, 5.4, 5.6, 5.8, 6.0, 6.2, &
                    6.4, 6.6, 6.8, 7.0, 7.2, 7.4, 7.6, 7.8, 8.0, 8.2, 8.4, &
                    8.6 /)
  
 INTEGER, PARAMETER :: NumberOfMiu = 20 
 REAL, DIMENSION(NumberOfMiu), PARAMETER :: &
 miuScale = (/ 0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40, 0.45, &
               0.50, 0.55, 0.60, 0.65, 0.70, 0.75, 0.80, 0.85, 0.90, &
               0.95, 1.00 /) 
  
  
 contains
  
 #if AMS_INST | ASTER_INST | AVIRIS_INST | RSP_INST | EPIC_INST | SSFR_INST
 subroutine fascode_vals(xpoint, ypoint, P, model, emission, trans2way, trans1way, trans2way_mean, &
                         trans2way_low, trans2way_high, trans1way_low, trans1way_high) 
     use global_model_grids
     
     real, intent(in) :: P
     type(ancillary_type), intent(in) :: model
     real, intent(inout) :: emission, trans2way, trans1way
     real, intent(inout) :: trans2way_mean, trans1way_low, trans1way_high, trans2way_low, trans2way_high
     integer, intent(in) :: xpoint, ypoint
  
 end subroutine fascode_vals
  
 #else
 subroutine fascode_vals(xpoint, ypoint, P, model, emission, trans2way, trans1way, trans2way_mean, &
                         trans2way_low, trans2way_high, trans1way_low, trans1way_high) 
     
     use science_parameters, only: model_levels, delta_pc
     use rtm_support, only: rtm_rad_atm_clr, rtm_trans_2way, rtm_trans_atm_clr, &
                     rtm_trans_2way_mean, rtm_trans_2way_low, rtm_trans_atm_clr_low, &
                     rtm_trans_2way_high, rtm_trans_atm_clr_high
     use core_arrays, only : cloud_height_method
     use global_model_grids
  
     real, intent(in) :: P
     type(ancillary_type), intent(in) :: model
     real, intent(inout) :: emission, trans2way, trans1way
     real, intent(inout) :: trans2way_mean, trans1way_low, trans1way_high, trans2way_low, trans2way_high
     integer, intent(in) :: xpoint, ypoint
  
     integer :: i, n, ilow, ihi, istart, iend
     
     real :: elow, ehi, t2low, t2hi, t1low, t1hi
  
     istart = model%trop_level
     iend = model%surface_level
  
     do i=istart, iend
         if (p <= model%pressure_profile(i)) exit
     end do
     
     if (i > iend) i = iend
     
     emission = rtm_rad_atm_clr(2, i)
     trans2way = rtm_trans_2way(2,i) ! 2 way correction for 3.7um, to be used later
     trans1way = rtm_trans_atm_clr(2,i)  ! 1 way correction for 3.7um, to be user later
  
     if (cloud_height_method(xpoint, ypoint) == 6) then 
         trans2way_mean = rtm_trans_2way_mean(i) ! 2 way correction for 3.7um, to be used later
  
         trans2way_low = rtm_trans_2way_low(2,i)
         trans2way_high = rtm_trans_2way_high(2,i)
  
         trans1way_low = rtm_trans_atm_clr_low(2,i)  
         trans1way_high = rtm_trans_atm_clr_high(2,i)    
  
     else
  
         do ilow=istart, iend
             if ((p-delta_pc) <= model%pressure_profile(ilow)) exit
         end do
         do ihi=istart, iend
             if ((p+delta_pc) <= model%pressure_profile(ihi)) exit
         end do
         
         t2low = rtm_trans_2way(2,ilow) ! 2 way correction for 3.7um, to be used later
         t1low = rtm_trans_atm_clr(2,ilow)  ! 1 way correction for 3.7um, to be user later
         
         t2hi = rtm_trans_2way(2,ihi) ! 2 way correction for 3.7um, to be used later
         t1hi = rtm_trans_atm_clr(2,ihi)  ! 1 way correction for 3.7um, to be user later
  
         trans2way_mean = ( abs(t2low - trans2way) + abs(t2hi - trans2way) ) / 2.0
  
         trans2way_low = t2low
         trans2way_high = t2hi
  
         trans1way_low = t1low
         trans1way_high = t1hi
  
         
     endif
  
  
  
 end subroutine fascode_vals
 #endif
  
 subroutine atmospheric_correction(xpoint,ypoint, iteration, meas_out, model, debug,status)
  
    use generalauxtype
    use science_parameters
    use core_arrays
    use libraryarrays
    use libraryinterpolates
    use mod06_run_settings
    use nonscience_parameters
    use global_model_grids
    use specific_ancillary
    use ch_xfr, only: c2_cmp_there
    implicit none
  
    logical, intent(in)  :: debug
    integer, intent(in)  :: xpoint,ypoint, iteration
    integer, intent(out) :: status
    real, dimension(:), intent(inout) :: meas_out
    type(ancillary_type), intent(in) :: model
  
    real                 :: oneway_correction(2), twoway_correction(trans_nband), &
                            oneway_standarddev(2),twoway_standarddev(trans_nband), &
                            twoway_correction_lower(trans_nband), twoway_correction_upper(trans_nband)
    real                 :: cos_solarzenith, cos_viewingangle
    integer              :: errorlevel
    integer, parameter   :: oneway_correction_index = 2, twoway_correction_index =trans_nband
    real                 :: angle
     real :: emission37, trans2way37, trans1way37
     real :: trans2way37_mean, trans2way37_low, trans2way37_high, trans1way37_low, trans1way37_high
     real :: o3_trans_low, o3_trans_high, o3_trans_550, o3_trans_470
     
     real :: sdev_2way(trans_nband), sdev_1way(2), oneway_correction_lower(2), oneway_correction_upper(2)
     
    status = 0
  
     call setup_atm_corr
  
  
    cos_solarzenith = cos(d2r*solar_zenith_angle(xpoint,ypoint))
   
    cos_viewingangle = cos(d2r*sensor_zenith_angle(xpoint,ypoint))
  
    angle = 1./cos_solarzenith + 1./cos_viewingangle
  
 !  currently the transmittance table breaks the MODIS band ordering by
 !  choosing index the wavelength without regard to the SDS ordering implicit
 !  in the data
 !  bottom line; in the transmittance data .94 site between .86 and 1.2
 !               in the measurements .94 sits betweem 2.1 and 3.7 (because .94 is in the 1km dataset)
  
    if ( cloud_top_pressure(xpoint,ypoint)     >= 0. .and.  &
         cloud_top_temperature(xpoint, ypoint) >= 0. ) then
  
      call gettransmittancedata(oneway_correction_index, &
                            twoway_correction_index, &
                            transmit_correction_table, &
                            transmit_stddev_table, &
                            cloud_top_pressure(xpoint,ypoint), &
                            abovecloud_watervapor(xpoint,ypoint), &
                            cos_solarzenith, &
                            cos_viewingangle, &
                            oneway_correction, &
                            twoway_correction, &
                            oneway_standarddev, &
                            twoway_standarddev, &
                            errorlevel)
     if (have_fascode) then 
      call gettransmittance_simple( &
                            transmit_correction_table, &
                            cloud_top_pressure(xpoint,ypoint), &
                            abovecloud_watervapor(xpoint,ypoint)*(1. - watervapor_error), &
                            cos_solarzenith, &
                            cos_viewingangle, &
                            twoway_correction_lower, &
                            errorlevel)
      call gettransmittance_simple( &
                            transmit_correction_table, &
                            cloud_top_pressure(xpoint,ypoint), &
                            abovecloud_watervapor(xpoint,ypoint)*(1. + watervapor_error), &
                            cos_solarzenith, &
                            cos_viewingangle, &
                            twoway_correction_upper, &
                            errorlevel)
     else
      call gettransmittancedata(oneway_correction_index, &
                            twoway_correction_index, &
                            transmit_correction_table, &
                            transmit_stddev_table, &
                            cloud_top_pressure(xpoint,ypoint), &
                            abovecloud_watervapor(xpoint,ypoint)*(1. - watervapor_error), &
                            cos_solarzenith, &
                            cos_viewingangle, &
                            oneway_correction_lower, &
                            twoway_correction_lower, &
                            sdev_1way, &
                            sdev_2way, &
                            errorlevel)
      call gettransmittancedata(oneway_correction_index, &
                            twoway_correction_index, &
                            transmit_correction_table, &
                            transmit_stddev_table, &
                            cloud_top_pressure(xpoint,ypoint), &
                            abovecloud_watervapor(xpoint,ypoint)*(1. + watervapor_error), &
                            cos_solarzenith, &
                            cos_viewingangle, &
                            oneway_correction_upper, &
                            twoway_correction_upper, &
                            sdev_1way, &
                            sdev_2way, &
                            errorlevel)  
     endif
  
  
 ! the values will never be fill, because there has been a change in the gettransmittance routines 
 ! used above. G.Wind 2.2.05
  
  
 ! we need to re-order the mean
  
      meandelta_trans(1:trans_nband) = ( abs(twoway_correction_lower - twoway_correction) + &
                                             abs(twoway_correction_upper - twoway_correction) ) /2.
      transmittance_twoway(1:trans_nband) = twoway_correction
      transmittance_stddev(1:trans_nband) = twoway_standarddev
  
 ! we are not using the table, so transmittance for 3.7um is from FASCODE and there is no standard deviation as it's not 
 ! statistical. 
     if (do_37) then
         if (have_fascode) then  
             call fascode_vals(xpoint, ypoint, cloud_top_pressure(xpoint, ypoint), model, emission37, trans2way37, trans1way37, &
                         trans2way37_mean, trans2way37_low, trans2way37_high, trans1way37_low, trans1way37_high )
                                                 
             meandelta_trans(band_0370) = trans2way37_mean
             transmittance_stddev(band_0370) = 0. 
         else
             emission37 = 0.
         endif
  
     endif
  
    else
  
 !    these values are set to have a benign effect on subsequent computations
 !    status is set an an early return performed.  It is expected that subsequest
 !    retrieval computation on this scene will be terminated  
      oneway_correction = 1.
      twoway_correction = 1.
  
      transmittance_twoway(:) = 1.
      transmittance_stddev(:) = 0.
      meandelta_trans(:) = 0.
  
      status = 1
      
      emission37 = 0.
      trans2way37 = 1.
      trans1way37 = 1.
      trans2way37_mean = 1.
      
      return
    endif
 !  Beer's law for ozone transmittance from M.D. King
    ozone_transmittance = exp(-1.0 * column_ozone(xpoint, ypoint)* ozone_absorp_coeff * angle)
    
    o3_trans_low = exp(-1.0 * 0.8*column_ozone(xpoint, ypoint)* ozone_absorp_coeff * angle)
    o3_trans_high = exp(-1.0 * 1.2*column_ozone(xpoint, ypoint)* ozone_absorp_coeff * angle)
    mean_delta_ozone = (abs(o3_trans_low - ozone_transmittance) + abs(o3_trans_high - ozone_transmittance))/2.0
     
 ! the only thing that channels 3 and 4 have absorption from is ozone
 ! used for discriminating very uniform stratocumulus from aerosol when CSR=2 for the GMAO FFNET code
  
     if (have_band(1)) then
         o3_trans_470 = exp(-1.0 * column_ozone(xpoint, ypoint)* o3_coef_470 * angle)
         if (meas_out(band_0047) /= fillvalue_real)&
             meas_out(band_0047) = meas_out(band_0047) / o3_trans_470
     endif
     if (have_band(2)) then
         o3_trans_550 = exp(-1.0 * column_ozone(xpoint, ypoint)* o3_coef_550 * angle)
         if (meas_out(band_0055) /= fillvalue_real)&
             meas_out(band_0055) = meas_out(band_0055) / o3_trans_550
     endif
  
     if (have_band(3)) then
        if (meas_out(band_0065) /= fillvalue_real) &
             meas_out(band_0065) = meas_out(band_0065)/ &
                           (ozone_transmittance*twoway_correction(band_0065)) ! .68
        meandelta_trans(band_0065) = meandelta_trans(band_0065) * ozone_transmittance
        transmittance_twoway(band_0065) = transmittance_twoway(band_0065) * ozone_transmittance
        transmittance_stddev(band_0065 ) = transmittance_stddev(band_0065) * ozone_transmittance
     endif
  
     if (have_band(4)) then
        if (meas_out(band_0086) /= fillvalue_real) &
             meas_out(band_0086) = meas_out(band_0086)/&
                             twoway_correction(band_0086) ! .86
     endif
     ! WDR orig, out: if (have_band(5)) then
     if (c2_cmp_there(band_0935) == 1 ) then
         if (meas_out(band_0935) /= fillvalue_real) &
             meas_out(band_0935) = meas_out(band_0935)/&
                                  twoway_correction(band_0935) ! .94
     endif
     if (have_band(6)) then
        if (meas_out(band_0124) /= fillvalue_real) &
            meas_out(band_0124) = meas_out(band_0124)/&
                           twoway_correction(band_0124) ! 1.2
     endif
     if (have_band(7)) then
        if (meas_out(band_0163) /= fillvalue_real) &
             meas_out(band_0163) = meas_out(band_0163)/&
                                    twoway_correction(band_0163) ! 1.6
     endif
     if (have_band(8)) then
        if (meas_out(band_0213) /= fillvalue_real) &
             meas_out(band_0213) = meas_out(band_0213)/&
                                    twoway_correction(band_0213) ! 2.1
     endif
  
    if (have_band(9)) then
       if (meas_out(band_0226) /= fillvalue_real) &
            meas_out(band_0226) = meas_out(band_0226)/&
                                 twoway_correction(band_0226) ! 2.2
     endif
  
     
    thermal_correction_twoway(1:2) = 1.
    thermal_correction_oneway(1:2) = 1.
  
 !  3.7 transmittance corrections cannot be applied until we have seperated the visible/thermal contribution to 
 !  the total radiance
     if (do_37) then 
 ! subtract the emission due to atmosphere above cloud from measured radiance
 #ifndef SIM_NORAD
         meas_out(band_0370) = meas_out(band_0370) - emission37
 #endif
  
         if (have_fascode) then 
             transmittance_twoway(band_0370) = trans2way37
         else
 #if !VIIRS_INST && !AHI_INST 
             transmittance_twoway(band_0370) = transmittance_twoway(band_0370)
             trans2way37 = twoway_correction(band_0370)
             trans1way37 = oneway_correction(1)
             trans1way37_low = oneway_correction_lower(1)
             trans1way37_high = oneway_correction_upper(1)
             trans2way37_low = twoway_correction_lower(band_0370)
             trans2way37_high = twoway_correction_upper(band_0370)
 #endif
         endif
     else
         trans2way37 = 1.
         trans1way37 = 1.
         trans1way37_low = 1.
         trans1way37_high = 1.
         trans2way37_low = 1.
         trans2way37_high = 1.
     endif
     
       thermal_correction_twoway(1) = trans2way37 ! 2 way correction for 3.7um, to be used later
       thermal_correction_oneway(1) = trans1way37  ! 1 way correction for 3.7um, to be user later
  
       thermal_correction_oneway_low(1) = trans1way37_low 
       thermal_correction_oneway_high(1) = trans1way37_high
       thermal_correction_twoway_low(1) = trans2way37_low
       thermal_correction_twoway_high(1) = trans2way37_high
  
  
 end subroutine atmospheric_correction
  
 SUBROUTINE gettransmittancedata(KDIM_1WAY,   &
                                 KDIM_2WAY,   &
                                 BigTauTable, &
                                 BigSdevTable,&
                                 p,           &
                                 tpw,         &
                                 miu0,        &
                                 miu1,        &
                                 tau1way,     &
                                 tau2way,     &
                                 sdev1way,    &
                                 sdev2way,    &
                                 errorLevel)
 !................................................................................
 !
 !................................................................................
 ! !F90
 !
 ! !DESCRIPTION:     
 !     given pressure, integrated precipitable water amount, miu0 and miu1, this
 !     program returns: 
 !     (1) one way transmittance for 3.7 and 11 um bands.
 !     (2) two way transmittance for 0.67, 0.86, 1.2, 1.6, 2.1 and 3.7 um bands.
 ! 
 !
 ! !INPUT PARAMETERS:
 !
 !       KDIM_1WAY = delcared band-dimension for 1-way tau array; >= 2
 !       KDIM_2WAY = declared band-dimension for 2-way tau array; >= 6
 !
 !       p(XDIM, YDIM) = pressure height in mb.
 !     tpw(XDIM, YDIM) = precipitable water (g/cm2) from level p to TOA.
 !    miu0(XDIM, YDIM) = cosine of solar zenith angle.
 !    miu1(XDIM, YDIM) = cosine of viewing zenith at level p.
 !
 !
 ! !OUTPUT PARAMETERS:
 !
 !      tau1way(XDIM, YDIM, 2) = contains 1-way spectral transmittance information.
 !      tau2way(XDIM, YDIM, 6) = contains 2-way spectral transmittance information.
 !      sdev1way(XDIM, YDIM, 2) = contains 1-way standard deviations.
 !      sdev2way(XDIM, YDIM, 6) = contains 2-way standard deviations.
 ! 
 !      *** Note ***: Fill value for tau1way and tau2way is -1.
 !
 !      errorLevel = return error status
 !                   0 : none. 
 !                   1 : second nearest neighbour maybe used (warning).
 !                   2 : cannot find required files (fatal)
 ! 
 !
 !
 !
 ! !REQUIRED EXTERNAL PROGRAMS: 
 !                 netcdf.f90 and ModisAuxType.f90
 !................................................................................
 ! !Revision History:
 !
 !
 ! Revision 2.5 2004/06/30    MAG
 ! fixed confusion between transmittance table and standard deviation table
 !
 ! Revision 2.5 2004/05/06    wind
 !   added standard deviation processing in order to facilitate uncertainty calculations
 ! Revision 2.0  2002/10/25           wind,mag
 ! add bilinear smoothing to fields
 !
 ! Revision 1.5  2001/05/01 18:50:08  jyli
 ! read entire transmittance table once to speed up process.
 !
 ! Revision 1.4  2001/04/18 21:56:35  jyli
 ! return both 1-way (for longwave) and 2-way (for shortwave) tau.
 !
 ! Revision 1.3  2001/04/06 19:11:57  jyli
 ! first delivery version which includs MODIS run environment stuff.
 !
 ! Revision 1.1  2001/04/06 18:46:11  jyli
 ! Initial revision
 !
 !................................................................................
 ! !Team-unique Header:
 !             Cloud Retrieval Group, NASA/GSFC
 !
 ! !PROGRAMMER:
 !
 !             Mark Gray (L3 Com)
 !             Climate and Radiation Branch
 !             NASA Goddard Space Flight Center
 !             Greenbelt, Maryland 20771, U.S.A.
 !
 !             Gala Wind (L3 Com)
 !             Climate and Radiation Branch
 !             NASA Goddard Space Flight Center
 !             Greenbelt, Maryland 20771, U.S.A.
 !
 !             Jason Li (Emergent-IT)
 !             Climate and Radiation Branch
 !             NASA Goddard Space Flight Center
 !             Greenbelt, Maryland 20771, U.S.A.
 !*******************************************************************************
 ! !END
 use modis_numerical_module, only: linearinterpolation
 use mod06_run_settings
 use specific_ancillary
 use ch_xfr, only : c2_cmp_there
  
 IMPLICIT NONE
  
  
 !... i/o variable delcarations:
  
 INTEGER, intent(in) :: KDIM_1WAY, KDIM_2WAY
 REAL , intent(in)            :: p, tpw, miu0, miu1
 REAL, intent(out),DIMENSION(KDIM_1WAY) :: tau1way, sdev1way
 REAL, intent(out),DIMENSION(KDIM_2WAY) :: tau2way, sdev2way
 INTEGER, intent(out)                               :: errorLevel
 REAL, intent(in), DIMENSION(:,:,:,:)   :: bigTauTable, bigSdevTable
  
  
 INTEGER :: i, pIndex, miu1Index, miu2Index, pwIndex, set_fill
 INTEGER, DIMENSION(4) :: start, counts, mindex
  
 INTEGER           :: NumberOfPixels_1km, NumberOfLines_1km
  
 REAL            :: miu, tauValue
 real            :: tau11, tau12, tau21, tau22, tau_mid1, tau_mid2
 real            :: sdev11, sdev12, sdev21, sdev22, sdev_mid1, sdev_mid2
 Integer         :: first_pwIndex, first_miuindex, second_pwIndex, second_miuindex
 real, dimension(8) :: temp_trans, temp_sdev
 !................................................................................
  
  
 ! initialize values
 temp_trans = -999.
 temp_sdev = -999.
  
  
 errorlevel = 0
  
    miu = (miu0*miu1)/(miu0+miu1)
  
 !
 !... find array indicies along p, pw and miu dimension:
 !
    pindex = nint( (p - heightscale(1)) / 100.0 ) + 1
    set_fill = 0 
    if (p< 0. .and. tpw < 0.) set_fill = 1
  
    if(pindex < 1) then
       pindex = 1
       errorlevel = 1
    elseif(pindex > numberofpressureheights) then
       pindex = numberofpressureheights
       errorlevel = 1
    endif
  
    miu1index = nint( (miu1 - miuscale(1)) / 0.05 ) + 1
    if (miu1index < 1) then
        miu1index = 1
        errorlevel = 1
    elseif(miu1index > numberofmiu) then
        miu1index = numberofmiu
        errorlevel = 1
    endif
  
    miu2index = nint( (miu - miuscale(1)) / 0.05 ) + 1
    if (miu2index < 1) then
        miu2index = 1
        errorlevel = 1
    elseif(miu2index > numberofmiu) then
        miu2index = numberofmiu
        errorlevel = 1
    endif
  
  
 ! changed by G. Wind 5.6.04 to match the new bins in transmittance table
    if(tpw < 0.0) then
        pwindex = 1
        errorlevel = 1
    elseif(tpw < 0.2) then
        pwindex = nint( (tpw- pwscale(1)) / 0.02 ) + 1
    else
        pwindex = nint( (tpw - pwscale(11)) / 0.20 ) + 11
        if(pwindex > numberofprecipitablewater) then
           pwindex = numberofprecipitablewater
           errorlevel = 1
        endif
    endif
  
    !... 
    !... 1-way case (for longwave 3.7 and 11 micron bands):
    !... 
  
    mindex(1) = pindex
    mindex(2) = pwindex
    mindex(3) = miu1index
    mindex(4) = 1
  
    tauvalue = bigtautable(mindex(1),mindex(2),mindex(3),mindex(4))
  
 ! we need to step back until we find a value that is non-fill, the second neighbor will not be
 ! enough in many cases
  
    if (tauvalue < 0.0) then  ! do I need 2nd nearest neighbour?
 !       errorLevel = 1
 !       mindex(2) = pwIndex - 1  
       
        do i=mindex(2), 1, -1
           tauvalue =  bigtautable(mindex(1),i,mindex(3),mindex(4))
           if (tauvalue >= 0.) exit
        end do
  
 ! now we have a good value of trans, non-fill, so we set the index accordingly
        mindex(2) = i
  
    endif
  
 !  grab the 4 points for bilinear interpolation: pw1,u1 pw2,u1 pw1,u2, pw2,u2
  
    if (pwscale(mindex(2)) >= tpw ) then
       second_pwindex = mindex(2)
       if (mindex(2) == 1) then 
          first_pwindex = mindex(2)
       else
          first_pwindex = mindex(2) - 1
          if (bigtautable(mindex(1), first_pwindex, mindex(3), mindex(4)) < 0.) &
             first_pwindex = mindex(2)
       endif
  
    else
       first_pwindex = mindex(2)
       if (mindex(2) == numberofprecipitablewater) then 
          second_pwindex = mindex(2) 
       else 
          second_pwindex = mindex(2) + 1
          if (bigtautable(mindex(1),second_pwindex,mindex(3),mindex(4)) < 0.0) &
             second_pwindex = mindex(2)
       endif
    endif
  
    if (miuscale(miu1index) >= miu1 ) then
       first_miuindex = miu1index -1
       second_miuindex = miu1index
    else
       first_miuindex = miu1index
       second_miuindex = miu1index + 1
    endif
  
    if( c2_cmp_there(band_0370) == 1) then
      do i = 1, numberoflongwavelengths
         mindex(4) = bandindexmaplw(i)
   
         tau11 = bigtautable(pindex, first_pwindex, first_miuindex, mindex(4))
         tau12 = bigtautable(pindex, first_pwindex, second_miuindex, mindex(4))
   
         tau21 = bigtautable(pindex, second_pwindex, first_miuindex, mindex(4))
         tau22 = bigtautable(pindex, second_pwindex, second_miuindex, mindex(4))
   
         sdev11 = bigsdevtable(pindex, first_pwindex, first_miuindex, mindex(4))
         sdev12 = bigsdevtable(pindex, first_pwindex, second_miuindex, mindex(4))
   
         sdev21 = bigsdevtable(pindex, second_pwindex, first_miuindex, mindex(4))
         sdev22 = bigsdevtable(pindex, second_pwindex, second_miuindex, mindex(4))
   
   
   !     first we interpolate along the PW lines
         tau_mid1 = linearinterpolation( (/pwscale(first_pwindex), pwscale(second_pwindex) /), &
                                         (/tau11, tau21 /), &
                                         tpw)
         tau_mid2 = linearinterpolation( (/pwscale(first_pwindex), pwscale(second_pwindex) /), &
                                         (/tau12, tau22 /), &
                                         tpw)
   
         sdev_mid1 =linearinterpolation( (/pwscale(first_pwindex),  pwscale(second_pwindex)/), &
                                         (/sdev11, sdev21/), &
                                         tpw)
         sdev_mid2 =linearinterpolation( (/pwscale(first_pwindex),  pwscale(second_pwindex)/), &
                                         (/sdev12, sdev22/), &
                                         tpw)
   
   
   !     now we interpolate along the MIU line through tau_mid1 and tau_mid2
   !     (bilinear interpolation)
         tau1way(i) = linearinterpolation( (/miuscale(first_miuindex), miuscale(second_miuindex)/), &
                                            (/tau_mid1, tau_mid2/), &
                                            miu1)
         sdev1way(i) = linearinterpolation( (/miuscale(first_miuindex), miuscale(second_miuindex)/), &
                                            (/sdev_mid1, sdev_mid2/), &
                                            miu1) 
         if (set_fill == 1 )then
            tau1way(i) = -999.
            sdev1way(i) = -999.
         endif
      enddo
   endif
  
    !
    !... 2-way case (for shortwave channels upto 3.7 micron band):
    !
  
    mindex(2) = pwindex
    mindex(3) = miu2index
    mindex(4) = 1
  
 ! the pw_index values had already been determined and they don't change with changing miu, 
 ! so pray tell why should we recompute them. This code had been deleted. G.Wind 2.2.05
  
    if (miuscale(miu2index) >= miu ) then
       first_miuindex = miu2index -1
       second_miuindex = miu2index
    else
       first_miuindex = miu2index
       second_miuindex = miu2index + 1
    endif
  
    do i = 1, numberofshortwavelengths
       mindex(4) = bandindexmapsw(i)
  
       tau11 = bigtautable(pindex, first_pwindex, first_miuindex, mindex(4))
       tau12 = bigtautable(pindex, first_pwindex, second_miuindex, mindex(4))
  
       tau21 = bigtautable(pindex, second_pwindex, first_miuindex, mindex(4))
       tau22 = bigtautable(pindex, second_pwindex, second_miuindex, mindex(4))
  
       sdev11 = bigsdevtable(pindex, first_pwindex, first_miuindex, mindex(4))
       sdev12 = bigsdevtable(pindex, first_pwindex, second_miuindex, mindex(4))
  
       sdev21 = bigsdevtable(pindex, second_pwindex, first_miuindex, mindex(4))
       sdev22 = bigsdevtable(pindex, second_pwindex, second_miuindex, mindex(4))
  
 !     first we interpolate along the PW lines
       tau_mid1 = linearinterpolation( (/pwscale(first_pwindex), pwscale(second_pwindex) /), &
                                       (/tau11, tau21 /), &
                                       tpw)
       tau_mid2 = linearinterpolation( (/pwscale(first_pwindex), pwscale(second_pwindex) /), &
                                       (/tau12, tau22 /), &
                                       tpw)
 !     then we repeat for the standard deviations
       sdev_mid1 =linearinterpolation( (/pwscale(first_pwindex),  pwscale(second_pwindex)/), &
                                       (/sdev11, sdev21/), &
                                       tpw)
       sdev_mid2 =linearinterpolation( (/pwscale(first_pwindex),  pwscale(second_pwindex)/), &
                                       (/sdev12, sdev22/), &
                                       tpw)
 !     now we interpolate along the MIU line through tau_mid1 and tau_mid2
 !     nothing but rather standard bilinear interpolation
 !     and that is our answer.
       temp_trans(i) = linearinterpolation( (/miuscale(first_miuindex),  miuscale(second_miuindex)/), &
                                         (/tau_mid1, tau_mid2/), &
                                         miu)      
  
 !     and these are the final standard deviations
       temp_sdev(i) = linearinterpolation( (/miuscale(first_miuindex), miuscale(second_miuindex)/), &
                                          (/sdev_mid1, sdev_mid2/), &
                                          miu )
       if (set_fill == 1 )then
          temp_trans(i) = -999.
          temp_sdev(i) = -999.
       endif
          
    end do
  
  
 ! remap the bands, so they can be indexed by the main instrument band index
     call remap_bands(tau2way, temp_trans, sdev2way, temp_sdev)
  
  
 9999 RETURN
  
 END SUBROUTINE gettransmittancedata
 SUBROUTINE gettransmittance_simple( &
                                    BigTauTable, &
                                    p,           &
                                    tpw,         &
                                    miu0,        &
                                    miu1,        &
                                    tau2way,     &
                                    errorLevel)
 !................................................................................
 !
 !................................................................................
 ! !F90
 !
 ! !DESCRIPTION:     
 !     given pressure, integrated precipitable water amount, miu0 and miu1, this
 !     program returns: 
 !     (1) one way transmittance for 3.7 and 11 um bands.
 !     (2) two way transmittance for 0.67, 0.86, 1.2, 1.6, 2.1 and 3.7 um bands.
 ! 
 !
 ! !INPUT PARAMETERS:
 !
 !       KDIM_2WAY(xdim, ydim, n) = declared band-dimension for 2-way tau array; >= n bands
 !
 !       p(XDIM, YDIM) = pressure height in mb.
 !     tpw(XDIM, YDIM) = precipitable water (g/cm2) from level p to TOA.
 !    miu0(XDIM, YDIM) = cosine of solar zenith angle.
 !    miu1(XDIM, YDIM) = cosine of viewing zenith at level p.
 !
 !
 ! !OUTPUT PARAMETERS:
 !
 !      tau2way(XDIM, YDIM, n) = contains 2-way spectral transmittance information.
 ! 
 !      *** Note ***: Fill value for tau1way and tau2way is -1.
 !
 !      errorLevel = return error status
 !                   0 : none. 
 !                   1 : second nearest neighbour maybe used (warning).
 !                   2 : cannot find required files (fatal)
 ! 
 !................................................................................
 ! !Revision History:
 !................................................................................
 ! !Team-unique Header:
 !             Cloud Retrieval Group, NASA/GSFC
 !
 ! !PROGRAMMER:
 !
 !             Mark Gray (L3 Com)
 !             Climate and Radiation Branch
 !             NASA Goddard Space Flight Center
 !             Greenbelt, Maryland 20771, U.S.A.
 !
 !             Gala Wind (L3 Com)
 !             Climate and Radiation Branch
 !             NASA Goddard Space Flight Center
 !             Greenbelt, Maryland 20771, U.S.A.
 !
 !             Jason Li (Emergent-IT)
 !             Climate and Radiation Branch
 !             NASA Goddard Space Flight Center
 !             Greenbelt, Maryland 20771, U.S.A.
 !*******************************************************************************
 ! !END
  
 use modis_numerical_module, only: linearinterpolation
 use mod06_run_settings
 use specific_ancillary
  
 IMPLICIT NONE
  
  
  
 !... i/o variable delcarations:
  
 REAL , intent(in)            :: p, tpw, miu0, miu1
 REAL, intent(out),DIMENSION(:) :: tau2way 
 INTEGER, intent(out)                               :: errorLevel
 REAL, intent(in), DIMENSION(:,:,:,:)   :: bigTauTable
  
  
  
 INTEGER :: i, pIndex, miu1Index, miu2Index, pwIndex, set_fill
 INTEGER, DIMENSION(4) :: start, counts, mindex
  
 INTEGER           :: NumberOfPixels_1km, NumberOfLines_1km
  
 REAL            :: miu, tauValue
 real            :: tau11, tau12, tau21, tau22, tau_mid1, tau_mid2
 real            :: sdev11, sdev12, sdev21, sdev22, sdev_mid1, sdev_mid2
 Integer         :: first_pwIndex, first_miuindex, second_pwIndex, second_miuindex
 real, dimension(8) :: temp_trans, temp_sdev, sdev2way
 !................................................................................
  
  
 ! initialize values
 temp_trans = -999.
  
  
  
 errorlevel = 0
  
    miu = (miu0*miu1)/(miu0+miu1)
  
 !
 !... find array indicies along p, pw and miu dimension:
 !
    pindex = nint( (p - heightscale(1)) / 100.0 ) + 1
    set_fill = 0 
    if (p< 0. .and. tpw < 0.) set_fill = 1
  
    if(pindex < 1) then
       pindex = 1
       errorlevel = 1
    elseif(pindex > numberofpressureheights) then
       pindex = numberofpressureheights
       errorlevel = 1
    endif
  
    miu1index = nint( (miu1 - miuscale(1)) / 0.05 ) + 1
    if (miu1index < 1) then
        miu1index = 1
        errorlevel = 1
    elseif(miu1index > numberofmiu) then
        miu1index = numberofmiu
        errorlevel = 1
    endif
  
    miu2index = nint( (miu - miuscale(1)) / 0.05 ) + 1
    if (miu2index < 1) then
        miu2index = 1
        errorlevel = 1
    elseif(miu2index > numberofmiu) then
        miu2index = numberofmiu
        errorlevel = 1
    endif
  
  
 ! changed by G. Wind 5.6.04 to match the new bins in transmittance table\
    pwindex = 1
    if(tpw < 0.0) then
        errorlevel = 1
    elseif(tpw < 0.2) then
        pwindex = nint( (tpw- pwscale(1)) / 0.02 ) + 1
    else
        pwindex = nint( (tpw - pwscale(11)) / 0.20 ) + 11
        if(pwindex > numberofprecipitablewater) then
           pwindex = numberofprecipitablewater
           errorlevel = 1
        endif
    endif
  
    !
    !... 2-way case (for shortwave channels upto 3.7 micron band):
    !
  
    mindex(1) = pindex
    mindex(2) = pwindex
    mindex(3) = miu2index
    mindex(4) = 1
       
    if (pwindex  < 1) print*, "HERE" , mindex, tpw, p
       
    tauvalue =  bigtautable(mindex(1),mindex(2),mindex(3),mindex(4))
   
 ! we need to step back until we find a value that is non-fill, the second neighbor will not be
 ! enough in many cases
  
    if (tauvalue < 0.0) then  ! do I need 2nd nearest neighbour?
 !       errorLevel = 1
 !       mindex(2) = pwIndex - 1  
       
        do i=mindex(2), 1, -1
           tauvalue =  bigtautable(mindex(1),i,mindex(3),mindex(4))
           if (tauvalue >= 0.) exit
        end do
  
 ! now we have a good value of trans, non-fill, so we set the index accordingly
        mindex(2) = i
  
    endif
  
  
 !  grab the 4 points for bilinear interpolation: pw1,u1 pw2,u1 pw1,u2, pw2,u2
    if (pwscale(mindex(2)) >= tpw ) then
       second_pwindex = mindex(2)
       if (mindex(2) == 1) then 
          first_pwindex = mindex(2)
       else
          first_pwindex = mindex(2) - 1
          if (bigtautable(mindex(1), first_pwindex, mindex(3), mindex(4)) < 0.) &
             first_pwindex = mindex(2)
       endif
  
    else
       first_pwindex = mindex(2)
       if (mindex(2) == numberofprecipitablewater) then 
          second_pwindex = mindex(2) 
       else 
          second_pwindex = mindex(2) + 1
          if (bigtautable(mindex(1),second_pwindex,mindex(3),mindex(4)) < 0.0) &
             second_pwindex = mindex(2)
       endif
    endif
  
  
    if (miuscale(miu2index) >= miu ) then
       first_miuindex = miu2index -1
       second_miuindex = miu2index
    else
       first_miuindex = miu2index
       second_miuindex = miu2index + 1
    endif
  
    do i = 1, numberofshortwavelengths
       mindex(4) = bandindexmapsw(i)
  
       tau11 = bigtautable(pindex, first_pwindex, first_miuindex, mindex(4))
       tau12 = bigtautable(pindex, first_pwindex, second_miuindex, mindex(4))
  
       tau21 = bigtautable(pindex, second_pwindex, first_miuindex, mindex(4))
       tau22 = bigtautable(pindex, second_pwindex, second_miuindex, mindex(4))
  
 !     first we interpolate along the PW lines
       tau_mid1 = linearinterpolation( (/pwscale(first_pwindex), pwscale(second_pwindex) /), &
                                       (/tau11, tau21 /), &
                                       tpw)
       tau_mid2 = linearinterpolation( (/pwscale(first_pwindex), pwscale(second_pwindex) /), &
                                       (/tau12, tau22 /), &
                                       tpw)
  
 !     now we interpolate along the MIU line through tau_mid1 and tau_mid2
 !     nothing but rather standard bilinear interpolation
 !     and that is our answer.
       temp_trans(i) = linearinterpolation( (/miuscale(first_miuindex),  miuscale(second_miuindex)/), &
                                         (/tau_mid1, tau_mid2/), &
                                         miu)
  
       if (set_fill == 1 )then
          temp_trans(i) = -999.
       endif
  
    end do
  
  
 ! remap the bands, so they can be indexed by the main instrument band index
     temp_sdev(1) = -10000. ! we don't do standard deviation in this subroutine
     call remap_bands(tau2way, temp_trans, sdev2way, temp_sdev)
  
  
 9999 RETURN
  
 END SUBROUTINE gettransmittance_simple
  
  
  
  
 end module atmospheric_correction_module

ch_xfr

Definition: ch_xfr.f90:1

rtm_support::rtm_trans_atm_clr

real, dimension(nchan, model_levels), public rtm_trans_atm_clr

Definition: rtm_support.f90:47

specific_ancillary::setup_atm_corr

subroutine setup_atm_corr

Definition: specific_ancillary.f90:27

mod06_run_settings::band_0124

integer, parameter band_0124

Definition: mod06_run_settings.f90:26

core_arrays

Definition: core_arrays.f90:1

specific_ancillary::bandindexmapsw

integer, dimension(8) bandindexmapsw

Definition: specific_ancillary.f90:7

core_arrays::ozone_transmittance

real ozone_transmittance

Definition: core_arrays.f90:129

rtm_support::rtm_trans_atm_clr_low

real, dimension(nchan, model_levels), public rtm_trans_atm_clr_low

Definition: rtm_support.f90:58

core_arrays::thermal_correction_oneway

real, dimension(set_number_of_bands) thermal_correction_oneway

Definition: core_arrays.f90:133

core_arrays::thermal_correction_oneway_high

real, dimension(2) thermal_correction_oneway_high

Definition: core_arrays.f90:198

core_arrays::meandelta_trans

real, dimension(set_number_of_bands) meandelta_trans

Definition: core_arrays.f90:128

specific_ancillary::bandindexmaplw

integer, dimension(numberoflongwavelengths) bandindexmaplw

Definition: specific_ancillary.f90:8

core_arrays::thermal_correction_twoway_high

real, dimension(2) thermal_correction_twoway_high

Definition: core_arrays.f90:199

specific_ancillary::ozone_absorp_coeff

real, parameter ozone_absorp_coeff

Definition: specific_ancillary.f90:14

science_parameters::watervapor_error

real, parameter watervapor_error

Definition: science_parameters.f90:20

specific_ancillary::do_37

logical do_37

Definition: specific_ancillary.f90:19

core_arrays::cloud_top_pressure

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

Definition: core_arrays.f90:156

specific_ancillary::numberoflongwavelengths

integer, parameter numberoflongwavelengths

Definition: specific_ancillary.f90:6

mod06_run_settings::band_0047

integer, parameter band_0047

Definition: mod06_run_settings.f90:26

specific_ancillary::numberofshortwavelengths

integer numberofshortwavelengths

Definition: specific_ancillary.f90:5

libraryarrays

Definition: libraryarrays.f90:1

modis_numerical_module::linearinterpolation

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

Definition: modis_numerical_module.f90:143

specific_ancillary::trans_nband

integer, parameter trans_nband

Definition: specific_ancillary.f90:10

global_model_grids

Definition: global_model_grids.f90:1

core_arrays::transmittance_twoway

real, dimension(set_number_of_bands) transmittance_twoway

Definition: core_arrays.f90:130

mod06_run_settings::band_0370

integer, parameter band_0370

Definition: mod06_run_settings.f90:26

core_arrays::thermal_correction_twoway

real, dimension(set_number_of_bands) thermal_correction_twoway

Definition: core_arrays.f90:132

science_parameters

Definition: science_parameters.f90:1

mod06_run_settings::band_0226

integer, parameter band_0226

Definition: mod06_run_settings.f90:42

mod06_run_settings::band_0086

integer, parameter band_0086

Definition: mod06_run_settings.f90:26

specific_ancillary::have_fascode

logical have_fascode

Definition: specific_ancillary.f90:19

mod06_run_settings::band_0213

integer, parameter band_0213

Definition: mod06_run_settings.f90:26

rtm_support::rtm_trans_atm_clr_high

real, dimension(nchan, model_levels), public rtm_trans_atm_clr_high

Definition: rtm_support.f90:63

libraryinterpolates

Definition: libraryinterpolates.f90:1

rtm_support::rtm_rad_atm_clr

real, dimension(nchan, model_levels), public rtm_rad_atm_clr

Definition: rtm_support.f90:48

science_parameters::delta_pc

real, parameter delta_pc

Definition: science_parameters.f90:33

generalauxtype

Definition: GeneralAuxType.f90:1

core_arrays::solar_zenith_angle

real, dimension(:,:), allocatable solar_zenith_angle

Definition: core_arrays.f90:6

rtm_support::rtm_trans_2way_high

real, dimension(nchan, model_levels), public rtm_trans_2way_high

Definition: rtm_support.f90:62

modis_numerical_module

Definition: modis_numerical_module.f90:1

core_arrays::column_ozone

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

Definition: core_arrays.f90:161

mod06_run_settings

Definition: mod06_run_settings.f90:1

core_arrays::mean_delta_ozone

real mean_delta_ozone

Definition: core_arrays.f90:129

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

Definition: PrintDebug.hpp:38

rtm_support::rtm_trans_2way_low

real, dimension(nchan, model_levels), public rtm_trans_2way_low

Definition: rtm_support.f90:57

specific_ancillary::o3_coef_470

real, parameter o3_coef_470

Definition: specific_ancillary.f90:15

libraryarrays::transmit_stddev_table

real, dimension(:,:,:,:), allocatable transmit_stddev_table

Definition: libraryarrays.f90:38

specific_ancillary::o3_coef_550

real, parameter o3_coef_550

Definition: specific_ancillary.f90:16

mod06_run_settings::band_0163

integer, parameter band_0163

Definition: mod06_run_settings.f90:26

nonscience_parameters

Definition: nonscience_parameters.f90:4

core_arrays::cloud_top_temperature

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

Definition: core_arrays.f90:155

nonscience_parameters::fillvalue_real

real, parameter fillvalue_real

Definition: nonscience_parameters.f90:13

atmospheric_correction_module::atmospheric_correction

subroutine, public atmospheric_correction(xpoint, ypoint, iteration, meas_out, model, debug, status)

Definition: atmospheric_correction.f90:125

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

science_parameters::model_levels

integer, parameter model_levels

Definition: science_parameters.f90:36

global_model_grids::ancillary_type

Definition: global_model_grids.f90:14

mod06_run_settings::band_0065

integer, parameter band_0065

Definition: mod06_run_settings.f90:26

atmospheric_correction_module::gettransmittancedata

subroutine gettransmittancedata(KDIM_1WAY, KDIM_2WAY, BigTauTable, BigSdevTable, p, tpw, miu0, miu1, tau1way, tau2way, sdev1way, sdev2way, errorLevel)

Definition: atmospheric_correction.f90:405

mod06_run_settings::band_0935

integer, parameter band_0935

Definition: mod06_run_settings.f90:26

mod06_run_settings::band_0055

integer, parameter band_0055

Definition: mod06_run_settings.f90:26

core_arrays::abovecloud_watervapor

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

Definition: core_arrays.f90:160

specific_ancillary

Definition: specific_ancillary.f90:1