        Program PROCESS_Land

C-----------------------------------------------------------------------
C!F77
C
C!DESCRIPTION:
C
C   This program derives aerosol optical depth (AOD) at 10km x 10km
C   from cloud and water screened 500 meter and 250 meter MODIS measured
C   radiances over land. AOD are derived at 0.553 micron, by inverting
C   radiance data at 0.466,0.644 and 2.113 micron channels. Dark pixels
C   are identified by the 2.113 mid-shortwave-IR (SWIR) channel, and
C   sorted so that the brightest (50%) and darkest (20%) pixels are discarded
C   to reduce cloud and other contaminations. Surface reflectance in the 
C   the visible channels (0.466 and 0.644 micron) are estimated from
C   the 2.113 channel, using relationships parameterized from vegetation
C   indices in the mid-IR (MVI) channels of 1.24 and 2.113 microns), and
C   also as a function of the solar/surface/satellite scattering angles. 
C   Note, that the 2.113 SWIR is not assumed transparent to aerosol, so that
C   the inversion of the two visible and one SWIR channels results in
C   three parameters: 1) spectral AOD normalized to 0.553 micron
C   2) The spectral fine mode AOD (AODfine = AOD * non-dust weighting)
C   and 3) Surface reflectance at 2.113 micron.
C 
C   The aerosol models have been derived from the AERONET data base and
C   are considered "dynamic" models of the optical depth. 
C   and are fixed as a function of location and season. 
C 
c   The program assumes that the input data are in reflectance units
C   of PI*L/(F0*cos(theta0)) for all visible and mid-IR channels. 
C   The input data are on MODIS spatial resolutions where the pixel size 
C   of the 0.66 micron is double the resolution of the 0.466, 1.24 and 2.11 
C   micron channels
C
C   Note: the program uses input cloud mask file (cloudy=0,clear=1),
C   which is determined by the thresholds of visible reflectance and
C   the visible reflectance ratios.  It requires that for 0.47 micron
C   channel 75% of the subpixels (250 m resolution) to be cloud free
C   and 100% of the subpixels to be water free. For 0.66 micron channel
C   it requires 100% subpixels to be cloud and water free. The selection
C   of 20th to 50th percentile will eliminate the residual clouds.
C   A valid retrieval needs at least 12 remaining pixels.
C
C!INPUT PARAMETERS:
C
C   HANDLE_LUTs          Logical unit numbers to open multi-wavelength lookup tables
C   IMONTH               Calendar month
C   IDATA                Index of 10x10 km box
C   NUMSQ                Total number of 10x10 km boxes
C   MTHET0               Solar zenith angle
C   MTHET                viewing angle
C   MDPHI                Relative azimuth angle
C   LAT_CENTER           Center latitude of 10x10 km box
C   LON_CENTER           Center longitude of 10x10 km box
C
C!OUTPUT PARAMETERS:
C
C!REVISION HISTORY:
C $Log: Process_land_V6.f,v $
c 02/15/2006 levy
c Initial revision
c
c!TEAM-UNIQUE HEADER:
c
c Developed by MODIS Aerosol/Solar Water Vapor Retrieval Team
c GSFC, Greenbelt, MD
c
C
C!DESIGN NOTES:
C
C The following ERROR CODE indicates no retrieval possibly done:
C   Shana!!
C      1     Solar, satellite zenith angle and relative azimuth angles
C            are out of bound of look-up tables
C      2     Threshold for 2.13 um is not met in the entire grid box
C      3     The thresold for number of cloudfree points is not met.
C            (NUMCLDFREE see the parameter) or there is water.
C      4     The computed reflectance for 0.47 or 0.66 channels is out
C            of bounds from look-up table
C      5     Threshold for 2.13 um is not met in option 2 and 3
C      6     Aerosol type can not be determined
C
C The values of aerosol optical thickness corresponding to error code
C Error code       Value of aerosol optical thickness
C
C      1              -3.0
C      2              -4.0
C      3              -5.0
C      4              -6.0
C      5              -7.0
C      6              -8.0
C
C  Aerosol types (FTABLE)
C      1              Continental
C      2              Generic : SSA ~ 0.9
C      3              Smoke   : SSA ~ 0.85
C      4              Urban   : SSA ~ 0.95
C      5              Dust
C
C Procedure used in retrieving aerosol optical thickness:
C
C      0 = no dark targets possibly met by the following thresholds
C      1 = threshold used for wave 2.13 um 0.01 - 0.25
C      2 = threshold used for wave 2.13 um 0.25 - 0.40
C
C INTERNALS:
C
C   Subroutines:
C      RNLOOKUP      -  READS THE LOOK-UP TABLES
C      ERROR         -  SETS THE ERROR CODES
C      OUTPUT        -  WRITES THE OUTPUT TO HDF FILE
C      AEROSOL_MAP   -  Reads map of aerosol model as function of season and place
C      INTANGLE_NL_RAY  -  Interpolates LUT Rayleigh parameters to measured geometry
C      Process_land_Rob -  The main code for doing the radiance to aerosol inversion
C      COMPUTE_SCATTANGLE_LAND - Computes scattering angle
C
C
C   Variables:
C     MTHET0                         solar zenith angle (degree)
C     MTHET                          satellite viewangle angle (degree)
C     MDPHI                          relative azimuth in angle (degree)
C     MHGHT                          Topographic altitude (km)
C     yint644,yint466,slope466,slope644  
C                 characteristics of VIS/IR surface reflectance relationship
C
C!END
C----------------------------------------------------------------------
C  Parameters for the resolution of modis data

      IMPLICIT NONE
      SAVE

      INCLUDE 'mod04.inc'
      INCLUDE 'mod04_land.inc'

      INTEGER IFINISH,IMONTH,IPR,IPROCE
      REAL WEIGHT,SCAT_ANGLE_LAND

c L2 Reflectance data
      REAL REFW466_L,REFW550_L 
      REAL REFW644_L,REFW124_L,REFW164_L,REFW866_L,REFW212_L
      REAL SDW466_L,SDW550_L
      REAL SDW644_L,SDW124_L,SDW164_L,SDW866_L,SDW212_L

      REAL SAVERATIO
      REAL THR213MIN, THR213MAX

      REAL MTHET0,MTHET,MDPHI,MPHI,MPHI0,MHGHT
      REAL LATM,LONM,LAT_CENTER,LON_CENTER

      integer ij,ik ,Qcontrol_special_land

       REAL G_factor, DEGRAD


c  Parameters Direct from the new lookup table (NL0 = New Land Initialize)

      REAL SBAR_NL0(NLTHET0,NLTAU,NLWAV,NLTABLE)
      REAL OPTH_NL0(NLTAU,NLWAV,NLTABLE)
      REAL INT_NL0(NLPHI,NLTHE,NLTHET0,NLTAU,NLWAV,NLTABLE)
      REAL FdT_NL0(NLTHE,NLTHET0,NLTAU,NLWAV,NLTABLE)
      REAL Fd_NL0(NLTHET0,NLTAU,NLWAV,NLTABLE)
      REAL T_NL0(NLTHE,NLTHET0,NLTAU,NLWAV,NLTABLE)
      REAL MASSCOEF_NL0(NLTAU,NLWAV,NLTABLE)
      REAL EXTNORM_NL0(NLTAU,NLWAV,NLTABLE)

C  For Determining aerosol type from aerosol map
      INTEGER nlon, nlat
      PARAMETER (nlon = 360, nlat = 180)
      INTEGER AEROSOL(nlon,nlat)

c retrieved products
      REAL ERR644
      REAL TAULAND55,RHOSFC213
      REAL MASSCON,ANGSTROM
      REAL RHOSFC(NLWAV)
      REAL ETA
      REAL AOD(NLWAV),AODF(NLWAV), AODC(NLWAV)
      INTEGER AVE_COUNT
      INTEGER FTABLE_NL
      INTEGER ETA_FLAG
      INTEGER IAER,IERROR,IDATA

c dummy
      INTEGER IWAV 

c Rayleigh reflectance interpolated for angle

      REAL REF_RAY_NL(NLWAV)


C     End of Declarations
C ------------------------------------------------------- 
C
c testing Inputs 2001124.1525
c  x = 8, y = 34
      OPEN(87, FILE='land_retrieval.csv', STATUS='UNKNOWN')
      WRITE(87,*) 'AER_TYPE,AOD553,ETA,AODFine,Angstrom,FitErr'


      OPEN (21, file = './land_input.dat',
     +   status = 'old', form = 'formatted')
      READ (21,*) 
      READ (21,*) REFW466_L, REFW550_L, REFW644_L, 
     +   REFW866_L, REFW124_L, REFW164_L, REFW212_L
      READ (21,*) 
      READ (21,*) SDW466_L, SDW550_L, SDW644_L, 
     +   SDW866_L, SDW124_L, SDW164_L, SDW212_L
      READ (21,*) 
      READ (21,*) MTHET0, MPHI0, MTHET, MPHI  
      READ (21,*)
      READ (21,*) LON_center, LAT_center, MHGHT
      READ (21,*)
      READ (21,*) IMONTH
      CLOSE (21)
 
      MDPHI = abs(MPHI0 - MPHI -180.0)
      IF(MDPHI.GT.360.0) MDPHI=amod(MDPHI,360.0)
      IF(MDPHI.GT.180.0) MDPHI=360.0-MDPHI

c Start program 

      HANDLE_LUT466 = 10
      HANDLE_LUT553 = 11
      HANDLE_LUT644 = 12
      HANDLE_LUT213 = 13
      HANDLE_LUTMAP = 14

      OPEN (HANDLE_LUT466, file = './lookup_land_w0487.n1',
     +   status = 'old', form = 'formatted')
      OPEN (HANDLE_LUT553, file = './lookup_land_w0551.n1',
     +   status = 'old', form = 'formatted')
      OPEN (HANDLE_LUT644, file = './lookup_land_w0638.n1',
     +   status = 'old', form = 'formatted')
      OPEN (HANDLE_LUT213, file = './lookup_land_w2257.n1',
     +   status = 'old', form = 'formatted')
      OPEN (HANDLE_LUTMAP, file = './aerosol_land_map.v3',
     +   status = 'old', form = 'formatted')



C Subroutine RNLOOKUP reads the new look-up tables for land

      CALL RNLOOKUP(
     +    HANDLE_LUT466,HANDLE_LUT553,HANDLE_LUT644,HANDLE_LUT213,
     +    INT_NL0,Fd_NL0,T_NL0,OPTH_NL0,
     +    SBAR_NL0,MASSCOEF_NL0,EXTNORM_NL0)
C
C  Determine fine mode aeorsol location map
C
       CALL AEROSOL_MAP(HANDLE_LUTMAP,IMONTH, AEROSOL)

      CLOSE (HANDLE_LUT466)
      CLOSE (HANDLE_LUT553)
      CLOSE (HANDLE_LUT644)
      CLOSE (HANDLE_LUT213)
      CLOSE (HANDLE_LUTMAP)



c     Initialize retrievals

      ETA = -9999
      ERR644 = -9999
      MASSCON = -9999
      ANGSTROM = -9999
      DO IWAV = 1, NLWAV
        RHOSFC(NLWAV) = -9999
        AOD(NLWAV) = -9999
        AODF(NLWAV) = -9999
        AODC(NLWAV) = -9999
      ENDDO
      AVE_COUNT = -9999
      FTABLE_NL = -9999       

      IFINISH=0
      Qcontrol_special_land=0


C
C Checking the angle bounds
C

       IF(MTHET0 .GE. MINMTHET0 .AND. MTHET0 .LE. MAXMTHET0.AND.
     *    MTHET  .GE. MINMTHET  .AND. MTHET  .LE. MAXMTHET .AND.
     *    MDPHI  .GE. MINMPHI   .AND. MDPHI  .LE. MAXMPHI) THEN
c
        CALL  COMPUTE_SCATTANGLE_LAND(MTHET0,MTHET,MDPHI,IDATA,
     *  SCAT_ANGLE_LAND)

        CALL INTANGLE_RAY_NL(MTHET0,MTHET,MDPHI,
     *	  INT_NL0,REF_RAY_NL)

          
       LATM=Lat_center
       LONM=Lon_center
       IPROCE = 0

C  Select criterion for detecting dark targets using 2.1 micron channel
C  (reflectance at 2.1 micron between 0.01 and THR213MAX) if the second
C  criterion fails (i.e., IFINISH=0)
C

             THR213MIN = 0.01
             THR213MAX = 0.25
             G_factor=1./COS(DTR*MTHET)
     *                +1./SQRT(COS(DTR*MTHET0))
             G_factor=0.5*G_factor

             IF (REFW212_L .GE. THR213MIN
     *        .AND. REFW212_L .LE. THR213MAX) THEN
                IPROCE = 1
             ELSE IF (REFW212_L .GE. THR213MAX
     *        .AND. REFW212_L .LE. THR213MAX*G_factor
     *        .AND. REFW212_L .LE. 40.0) THEN
                THR213MIN= 0.25
                THR213MAX= 0.25*G_factor
                IPROCE = 2
             ELSE
                IPROCE = 0
             ENDIF
             IF (IPROCE .GE. 1) IFINISH = 1
             PRINT*, 'IPROCE ', IPROCE

C Now do Rob's new inversion routine
         If(IFINISH.EQ.1) Then
              IERROR=0
         CALL PROCESS_Land_Inver(IPROCE,
     *   MTHET0,MTHET,MDPHI,SCAT_ANGLE_LAND,IMONTH,LATM,LONM,MHGHT,
     *   REFW466_L,REFW550_L,REFW644_L,REFW866_L,
     *   REFW124_L,REFW164_L,REFW212_L,
     *   SDW466_L,SDW550_L,SDW644_L,SDW866_L,
     *   SDW124_L,SDW164_L,SDW212_L,
     *   INT_NL0,Fd_NL0,T_NL0,OPTH_NL0,SBAR_NL0,AEROSOL,REF_RAY_NL,
     *   ETA,ETA_FLAG,AOD,ERR644,RHOSFC,AVE_COUNT,FTABLE_NL,
     *   MASSCOEF_NL0,MASSCON,
     *   EXTNORM_NL0,AODF,AODC,ANGSTROM)
     
C Set IAER     
            IAER = FTABLE_NL
            
     
C     If 2.13 Um not met procedure and error is set     
         ELSE
           IERROR=4
           IPROCE=0
           IAER =-999
C ENDIF for Ifinish     
           ENDIF
        

          
 

C
C                        **Following else for 'ifang'

      ELSE
C                        **Write the output if the angles are out of
C                          bounds from the lookup table for different
C                          options with error indicators.
C
        IERROR=1
       
            
  
C
C Endif for checking the angle bounds
C
      ENDIF
  
           
   1  FORMAT(132A1)
      RETURN
      END


C***********************************************************************
      SUBROUTINE COMPUTE_SCATTANGLE_LAND(MTHET0,MTHET,MDPHI,IDATA,
     *SCAT_ANGLE_LAND)
C-----------------------------------------------------------------------
C!F77
C
C!DESCRIPTION:
C              This subroutine computes scattering angle from MODIS
C              geometry.
C
C!INPUT PARAMETERS:
C
C       MEHET0            Solar zenith angle
C       MTHET             Satellite Viewing angle
C       MDPHI             Relative azimuth angle
C       IDATA             Index of 10x10 km box
C
C!OUTPUT PARAMETERS:
C
C       SCAT_ANGLE_LAND   Scattering angle
C
C!REVISION HISTORY:
C 02/02/98 achu
C identified input/output arguments in INPUT / OUTPUT PARAMETERS sections
C in prolog.
C
C 01/28/98 fhliang
C fixed prolog.
C
C!TEAM-UNIQUE HEADER:
C
C Developed by MODIS Aerosol Team at NASA GSFC, Greenbelt, MD
C
C!DESIGN NOTES:
C
C!END
C-----------------------------------------------------------------------
C
      IMPLICIT NONE
      INCLUDE 'mod04.inc'
      REAL SCAT_ANGLE_LAND
      REAL MTHET0,MTHET,MDPHI
      INTEGER IDATA
      SAVE

C      DTR=ACOS(-1.)/180.
C      RTD=180./ACOS(-1.)
      SCAT_ANGLE_LAND=0.0

C
C Scattering angle of MODIS geometry
C

       SCAT_ANGLE_LAND= -COS(MTHET0*DTR)*COS(MTHET*DTR)
     *                  +SIN(MTHET0*DTR)*SIN(MTHET*DTR)
     *                  *COS(MDPHI*DTR)

       SCAT_ANGLE_LAND= ACOS(SCAT_ANGLE_LAND)*RTD

       RETURN
       END

C*****************************************************************
C ROBS NEW ROUTINES!
C*********************************************************************
       SUBROUTINE PROCESS_Land_Inver(IPROCE,
     *    MTHET0,MTHET,MDPHI,SCAT_ANGLE_LAND,IMONTH,LATM,LONM,MHGHT,
     *    REFW466,REFW553,REFW644,REFW866,REFW123,REFW163,REFW212,
     *    SDW466,SDW553,SDW644,SDW866,SDW123,SDW163,SDW212,
     *    INT_NL0,Fd_NL0,T_NL0,OPTH_NL0,SBAR_NL0,AEROSOL,REF_RAY_NL,
     *    ETA,ETA_FLAG,AOD,ERR644,RHOSFC,AVE_COUNT,FTABLE_NL,
     *    MASSCOEF_NL0,MASSCON,
     *    EXTNORM_NL0,AODF,AODC,ANGSTROM)

C-----------------------------------------------------------------------
C!F77
C
C!DESCRIPTION:
C   The subroutine derives non-cloudy aerosol optical thickness from MODIS
C   measured radiances over land, using the 0.47 (blue), 0.66 (red) 
C   and 2.13 (IR) micron channels.
C   Input data are cloud screened and averaged reflectance at at 10 x 10 km. 
C   Inversion:
C   A) Select fine mode model based on geography (dynamic aerosol models)
C   B) Mix fine mode and dust mode with selected discrete mixing ratios
C      best matches the observed spectral reflectance.
C
C   The program assumes that the input data are in reflectance units
C   of PI*L/(F0*cos(theta0)) for all channels.
C   It also assumes that the surface reflectance in the two visible channels
C   are functions of the surface reflectance in the IR. 
C   Nominally, red is one half of IR, and blue is one-half of red, respectively. 
C!INPUT PARAMETERS:
C
C   MTHET0                         solar zenith angle (degree)
C   MTHET                          satellite viewangle angle (degree)
C   MDPHI                          relative azimuth in angle (degree)
C   MHGHT                          Topographic altitude (km)
C   IMONTH                         Calendar month
C   LATM                           latitude of 10x10 km box
C   LONM                           longitude of 10x10 km box
C   REFW466,REFW644,REFW866,REFW212   averaged cloud screened reflectance
C   SDW470,SDW659,SDW866,SDW212       stddev cloud screened reflectance
C   INT_NL0			      reflectance from lookup table
C   Fd_NL0                           flux down from lookup table
C   T_NL0                           transmission from lookup table
C   OPTH_NL0                        optical thickness from lookup table
C   SBAR_NL0                        sbar from lookup table
C   MASSCOEF_NL0                    Mass Concentration from lookup table
C   EXTNORM_NL0                     Normalized Extinction Coeficients from LUT
C
C!OUTPUT PARAMETERS:
C
C  ETA                Fine mode (non dust) fraction
C  AOD                AOD at 4 wavelengths 
C  AODF/C             Fine and Coarse mode AOD at 4 wavelengths 
C  RHOSFC             Surface reflectance at 4 wavelengths
C  AVE_COUNT          Number of observations with error < 3%
C  FTABLE_NL          Which Fine model
C  MASSCON            Mass Concentration
C  ETA_FLAG           Whether ETA between 0.0 and 1.0
C  ANGSTROM           Angstrom Exponenent 0.47/0.66
C   
C  Aerosol types (FTABLE)
C      1              Continental
C      2              Generic : SSA ~ 0.9
C      3              Smoke   : SSA ~ 0.85
C      4              Urban   : SSA ~ 0.95
C      5              Dust
C
C INTERNALS:
C
C   Subroutines:
C      RNLOOKUP       -  READS THE LOOK-UP TABLES
C      SELECT_NLOOKUP -  SELECT LOOK-UP TABLE
C      INTANGLE_NL   -  Interpolates LUT parameters to measured geometry
C      TOA_SIMULATE  -  Computes simulated TOA reflectance for each TAU index
C      RETRIEVAL1    -  Derives AOT, ETA and 2.1 Surface Reflectance for
C                       IPROCE=1 (0.01 < rho2p1 < 0.25)
C      RETRIEVAL2    -  Derives AOT for Continental model, IPROCE=2(0.25 < rho2p1)
C      SHLSRT_NL     -  Sorts an array into ascending order
C      INTERP_EXTRAP -  Linear interpolation/extrapolation
C      PERCENTILE    -  Computes average value/standard deviation of selected
C                       pixels given percentiles
C      ERROR         -  SETS THE ERROR CODES
C      OUTPUT        -  WRITES THE OUTPUT TO HDF FILE
C      AEROSOL_MAP   -  Reads map of aerosol model as function of season and place
C      INT_ELEV      -  Interpolates lookup table at different wavelengths to 
C                       simulate elevation 
C
C
C!END
C----------------------------------------------------------------------
      IMPLICIT NONE
      INCLUDE 'mod04.inc'
      INCLUDE 'mod04_land.inc'

      INTEGER NUMSQDIM,NUMXBOX,NUMYBOX
      PARAMETER (NUMSQDIM=400)
      INTEGER RTN,START_DATA,END_DATA
      INTEGER IDATA,ISCAN

      INTEGER IFINISH,IMONTH,IPR,IPROCE
      REAL SCAT_ANGLE_LAND

      REAL MTHET0,MTHET,MPHI0, MPHI, MDPHI,MHGHT
      REAL LATM,LONM,LAT_CENTER,LON_CENTER
C      REAL WATER_VAPOR(NUMSQDIM),OZONE_COL(NUMSQDIM)
      CHARACTER * 5 IOFLAG

C*********************************************************************
C   Parameters used from look-up tables
C ********************************************************************
C   wav_nl=wavelengths(4)(.466,.533,.644,2.13 micrometers)       4
C   nlthet0=number of theta0(solar zenith angle)              9
C   thet0_nl=(0,12,24,36,48,54,60,66,72   degrees)
C   nlthe =number of the(observation zenith angle)           13
C   the_nl=(0,6,.....72 degrees every 6 degrees)
C   nlphi =number of phi(observation azimuth angle)          16
C   phi_nl=(0,12,24,36,48,60,72,84,96,
C        108,120,132,144,156,168,180)
C   nltau =number of opth(optical thickness)                 7 
C   opth_nl=(0.0,0.25,0.50,0.75,1.0,2.0,3.0,5.0)
C   lhght = number of hght(observation heights)              1
C   hght=(80.0 km)
C**********************************************************************

C  
C  For Determining aerosol type from aerosol map
      INTEGER nlon, nlat
      PARAMETER (nlon = 360, nlat = 180)
      INTEGER AEROSOL(nlon,nlat)
      INTEGER ilat, ilon
      INTEGER FTABLE,FTABLE_NL
      REAL LONM1
      REAL dlat, dlon
      PARAMETER (dlat = 0.5, dlon = 0.5)
      CHARACTER*5 AEROSOL_TYPEs(5), AEROSOL_TYPE




C     Data directly read from the entire set of LUTs

      REAL SBAR_NL0(NLTHET0,NLTAU,NLWAV,NLTABLE)
      REAL OPTH_NL0(NLTAU,NLWAV,NLTABLE)
      REAL INT_NL0(NLPHI,NLTHE,NLTHET0,NLTAU,NLWAV,NLTABLE)
      REAL Fd_NL0(NLTHET0,NLTAU,NLWAV,NLTABLE)
      REAL T_NL0(NLTHE,NLTHET0,NLTAU,NLWAV,NLTABLE)
      REAL MASSCOEF_NL0(NLTAU,NLWAV,NLTABLE)
      REAL EXTNORM_NL0(NLTAU,NLWAV,NLTABLE)

      integer ij,ik,set_counter_land,qcontrol_special 

C     Data for specific LUTs (specified by aerosol model)

      REAL INT_NL1(NLPHI,NLTHE,NLTHET0,NLTAU,NLWAV,NLSIZE)
      REAL Fd_NL1(NLTHET0,NLTAU,NLWAV,NLSIZE)
      REAL SBAR_NL1(NLTHET0,NLTAU,NLWAV,NLSIZE)
      REAL T_NL1(NLTHE,NLTHET0,NLTAU,NLWAV,NLSIZE)

C     Interpolated lookup table data

      REAL INT_NL(NLTAU,NLWAV,NLSIZE)
      REAL Fd_NL(NLTAU,NLWAV,NLSIZE)
      REAL T_NL(NLTAU,NLWAV,NLSIZE)
      REAL FdT_NL(NLTAU,NLWAV,NLSIZE)
      REAL SBAR_NL(NLTAU,NLWAV,NLSIZE)
      REAL OPTH_NL(NLTAU,NLWAV,NLSIZE)
      REAL EQWAV_NL(NLWAV)
      REAL MASSCOEF_NL(NLTAU,NLWAV,NLSIZE)
      REAL EXTNORM_NL(NLTAU,NLWAV,NLSIZE)

c     Rayleigh data
      REAL REF_RAY_NL(NLWAV)


C     Simulated Reflectancea
      REAL  RHO_S212(NLTAU,NLSIZE), RHOSTAR(NLTAU,NLWAV,NLSIZE)
      REAL yint_466,slope_466,yint_644,slope_644

C     Results

      REAL ERR644
      REAL TAULAND55,RHOSFC212
      REAL RHOSFC(NLWAV)
      REAL ETA
      INTEGER ETA_FLAG
      REAL AOD(NLWAV),AODF(NLWAV),AODC(NLWAV)
      INTEGER AVE_COUNT
      REAL MASSCON
      REAL ANGSTROM


C     Dummy Indices
     
      INTEGER NSOLUTION,ISIZE,IETA,JETA,IWAV
      INTEGER IPHI,ITHE,ITAU,ITHET0,J
      INTEGER NUMDATA

C     Geometry Indices

      INTEGER KSZAM1,KSZAP1,KTHEM1,KTHEP1,KPHIM1,KPHIP1
      REAL REFW466,REFW553,REFW644,REFW866,REFW123,REFW163,REFW212
      REAL SDW466,SDW553,SDW644,SDW866,SDW123,SDW163,SDW212

c     INPUT Surface reflectance and TAU
c     Match corresponding aerosol model with aerosol map index

      SAVE

      DATA AEROSOL_TYPEs/'CONTN','BACKG','SMOKE','URBAN','DUST'/
     

      

c     Dummy set some stuff

C Pre-processing

      REFW466L = REFW466
      REFW644L = REFW644
      REFW212L = REFW212
      REFW124L = REFW123
      REFW866L = REFW866
      REFW164L = REFW163
      REFW553L = REFW553

C   Following If statement checks for missing data in a box 10*10
C  in any one of wavelengths.
C  If there is missing data in any one wavelength
C  arrays are filled with _fillValues.

       IF (REFW466L.LT.99999 .AND. REFW644L .LT. 99999 .AND.
     *    REFW866L.LT.99999 .AND. REFW212L .LT. 99999) THEN


C Determine expected fine model (from date and location)

      ilon = 181 + NINT(LONM + dlon)
      ilat = 91 - NINT(LATM + dlat)

      FTABLE = AEROSOL(ilon,ilat) + 2
      AEROSOL_TYPE = AEROSOL_TYPEs(FTABLE)

C  Select LUT corresponding to season and location
C  Obtain fine mode LUT (choice of tables #2, #3 or #4)
C  Coarse mode LUT = #5
C  Continental mode LUT = #1

      IF (IPROCE .EQ. 1) THEN
        MTABLE_NL(1) = FTABLE
        MTABLE_NL(2) = DTABLE_NL
      ELSE IF (IPROCE .EQ. 2) THEN
        MTABLE_NL(1) = 1
        MTABLE_NL(2) = 1
      ENDIF
      FTABLE_NL = MTABLE_NL(1)


      DO iSIZE = 1, NLSIZE

         CALL SELECT_NLOOKUP(INT_NL0,Fd_NL0,T_NL0,OPTH_NL0,
     *      SBAR_NL0,MASSCOEF_NL0,EXTNORM_NL0,
     *      ISIZE,INT_NL1,Fd_NL1,T_NL1,SBAR_NL1,
     *      OPTH_NL,MASSCOEF_NL,EXTNORM_NL,MTABLE_NL(ISIZE))

      ENDDO


C Interpolate lookup tables to measured angles. 

	CALL  SET_index_inter_NL(MTHET0,MTHET,MDPHI,
     *        KSZAM1,KSZAP1,KTHEM1,KTHEP1,KPHIM1,KPHIP1)
 
C Interpolate lookup table to measured geometry. 
     
        CALL INTANGLE_NL(MTHET0,MTHET,MDPHI,
     *	  INT_NL1,Fd_NL1,T_NL1,SBAR_NL1,
     *    INT_NL,Fd_NL,T_NL,FdT_NL,SBAR_NL,
     *    KSZAM1,KSZAP1,KTHEM1,KTHEP1,KPHIM1,KPHIP1)

C Interpolate lookup table to target elevation by interpolating wavelengths

        CALL INT_ELEV(EQWAV_NL,INT_NL,Fd_NL,T_NL,FdT_NL,OPTH_NL,
     *        SBAR_NL,REF_RAY_NL,MHGHT)

C
C  Simulate TOA reflectance
C

           CALL TOA_SIMULATE(
     *        MTHET0,SCAT_ANGLE_LAND,
     *        INT_NL,FdT_NL,OPTH_NL,SBAR_NL,RHO_S212,RHOSTAR,
     *        yint_466,slope_466,yint_644,slope_644)

C	 
C DO RETRIEVAL!!!!!!!!!!!
C     Retrieval#1 for IPROCE=1, Retrieval#2 for IPROCE=2	 
C

           ETA_FLAG = 0
           IF (IPROCE .EQ. 1) THEN	
             CALL RETRIEVAL1 (
     *       yint_466,slope_466,yint_644,slope_644,
     *       OPTH_NL,MASSCOEF_NL,EXTNORM_NL,
     *       RHO_S212,RHOSTAR,
     *       AVE_COUNT,ETA,ETA_FLAG,AOD,AODF,AODC,
     *       RHOSFC212,ERR644,RHOSFC,MASSCON)
           ELSE IF (IPROCE .EQ. 2) THEN
            CALL RETRIEVAL2 (
     *       yint_466,slope_466,yint_644,slope_644,
     *       OPTH_NL,MASSCOEF_NL,RHO_S212,RHOSTAR,
     *       AVE_COUNT,ETA,AOD,
     *       RHOSFC212,ERR644,RHOSFC,MASSCON)
           ENDIF

C Calculate Angstrom Exponent

           IF (AOD(iwave_644) .GT. 0) THEN
             ANGSTROM = -1.0 * ALOG(AOD(iwave_466)/AOD(iwave_644))/
     *           ALOG(0.466/0.644)
           ENDIF
        
           WRITE(87,87) FTABLE_NL,AOD(iwave_553),ETA,AODF(iwave_553),ANGSTROM,ERR644


          ENDIF
          CLOSE (87) 
      RETURN		

 87   FORMAT(i3,",",5(f7.4,","))
 
      END

C*****************************************************************
      SUBROUTINE SELECT_NLOOKUP(
     *      INT_NL0,Fd_NL0,T_NL0,OPTH_NL0,SBAR_NL0,MASSCOEF_NL0,
     *      EXTNORM_NL0,
     *      ISIZE,INT_NL1,Fd_NL1,T_NL1,SBAR_NL1,
     *      OPTH_NL,MASSCOEF_NL,EXTNORM_NL,ITAB)
C-----------------------------------------------------------------------
C!F77
C
C!DESCRIPTION:
C              This subroutine selects look-up table (using ITABLE)
C	       Used to select both the "fine" and the "coarse/dust" land model
C
C!INPUT PARAMETERS:
C
C          INT_NL0      Reflectance
C          Fd_NL0       Flux Down
C          T_NL0        Transmssion factor
C          OPTH_NL0     Optical depth
C          SBAR_NL0     Sbar
C          MASSCOEF_NL0  MASSCON coeficient
C          EXTNORM_NL0  Normalized extinction coeficient
C          ITABLE     Index for desired lookup table
C          ISIZE      Index indicating "fine" or "coarse"
C
C!OUTPUT PARAMETERS:
C
C          INT_NL1        Reflectance
C          Fd_NL1       Flux Down
C          T_NL1        Transmssion factor
C          SBAR_NL1       Sbar
C          OPTH_NL       Optical depth
C          MASSCOEF_NL     Mass Conencentration coeficient
C          EXTNORM_NL     Normalized extinction coeficient
C
C!REVISION HISTORY:
C 5/1/98 DAC created
C
C!TEAM-UNIQUE HEADER:
C
C Developed by MODIS aerosol team at Goddard Space Flight Center,
C Greenbelt, Maryland.
C
C!END
C-----------------------------------------------------------------
C
      IMPLICIT NONE
      INCLUDE 'mod04_land.inc'

      INTEGER IWAV,ITAB,ITAU,ITHET0,ITHE,IPHI
      INTEGER ISIZE

      REAL OPTH_NL0(NLTAU,NLWAV,NLTABLE)
      REAL SBAR_NL0(NLTHET0,NLTAU,NLWAV,NLTABLE)
      REAL INT_NL0(NLPHI,NLTHE,NLTHET0,NLTAU,NLWAV,NLTABLE)
      REAL Fd_NL0(NLTHET0,NLTAU,NLWAV,NLTABLE)
      REAL T_NL0(NLTHE,NLTHET0,NLTAU,NLWAV,NLTABLE)
      REAL MASSCOEF_NL0(NLTAU,NLWAV,NLTABLE)
      REAL EXTNORM_NL0(NLTAU,NLWAV,NLTABLE)

      REAL INT_NL1(NLPHI,NLTHE,NLTHET0,NLTAU,NLWAV,NLSIZE)
      REAL Fd_NL1(NLTHET0,NLTAU,NLWAV,NLSIZE)
      REAL SBAR_NL1(NLTHET0,NLTAU,NLWAV,NLSIZE)
      REAL T_NL1(NLTHE,NLTHET0,NLTAU,NLWAV,NLSIZE)
      REAL OPTH_NL(NLTAU,NLWAV,NLSIZE)
      REAL MASSCOEF_NL(NLTAU,NLWAV,NLSIZE)
      REAL EXTNORM_NL(NLTAU,NLWAV,NLSIZE)
      SAVE

C
C Select look-up table variable for given ITABLE
C

      DO 10 IWAV=1,NLWAV
        DO 20 ITAU=1,NLTAU
          OPTH_NL(ITAU,IWAV,ISIZE) = OPTH_NL0(ITAU,IWAV,ITAB)
          MASSCOEF_NL(ITAU,IWAV,ISIZE) = MASSCOEF_NL0(ITAU,IWAV,ITAB)
          EXTNORM_NL(ITAU,IWAV,ISIZE) = EXTNORM_NL0(ITAU,IWAV,ITAB)
          DO 30 ITHET0=1,NLTHET0
            Fd_NL1(ITHET0,ITAU,IWAV,ISIZE) =
     +         Fd_NL0(ITHET0,ITAU,IWAV,ITAB) 
            SBAR_NL1(ITHET0,ITAU,IWAV,ISIZE) =
     +         SBAR_NL0(ITHET0,ITAU,IWAV,ITAB) 
            DO 40 ITHE=1,NLTHE
               T_NL1(ITHE,ITHET0,ITAU,IWAV,ISIZE) = 
     +            T_NL0(ITHE,ITHET0,ITAU,IWAV,ITAB) 
              DO 50 IPHI=1,NLPHI
                  INT_NL1(IPHI,ITHE,ITHET0,ITAU,IWAV,ISIZE) = 
     +              INT_NL0(IPHI,ITHE,ITHET0,ITAU,IWAV,ITAB) 
50            CONTINUE
40          CONTINUE
30        CONTINUE
20      CONTINUE
10     CONTINUE

      RETURN
      END


C*****************************************************************
       SUBROUTINE SET_index_inter_NL(MTHET0,MTHET,MDPHI,
     *            KSZAM1,KSZAP1,KTHEM1,KTHEP1,KPHIM1,KPHIP1)

C-----------------------------------------------------------------------
C!F77
C
C!DESCRIPTION: This subroutine sets the index for ending and starting
C             positions from the lookup table (based on input geometry)
C
C!INPUT PARAMETERS:
C              MTHET0      Measured solar Zenith angle
C              MTHET      Measured view  Zenith angle
C              MDPHI      Measured Azimuthal Angle
C
C              THET0      array of  solar Zenith angle in look_up table
C              THE        array of  satllite Zenith angle in look_up table
C              PHI        array of azimuth angle in look_up table
C
C!OUTPUT PARAMETERS:
C               KSZAM1     Starting Index for solar zenith angle
C               KSZAP1     Ending Index for solar zenith angle
C               KTHEM1     Starting Index for view angle
C               KTHEP1     Ending   Index for  view angle
C               KPHIM1     Starting Index for  azimuth angle
C               KPHIP1     Ending   Index for  azimuth angle
C
C!REVISION HISTORY:
C
C!TEAM-UNIQUE HEADER: - NOT YET DEFINED
C
C!REFERENCES AND CREDITS Shana Mattoo
C
C!DESIGN NOTES:
C
C!END
C----------------------------------------------------------------------

       IMPLICIT NONE
       SAVE

       INCLUDE 'mod04_land.inc'

       REAL MTHET0,MTHET,MDPHI,DEL
       INTEGER KSZAM1,KSZAP1,KTHEM1,KTHEP1,KPHIM1,KPHIP1
       INTEGER I

C      Initialize
 
       KSZAM1=0
       KSZAP1=0
       KTHEM1=0
       KTHEP1=0
       KPHIM1=0
       KPHIP1=0

C      Find THET0 indices
       DO I = 1,NLTHET0-1
          IF ((MTHET0 .GE. THET0_NL(I)) .AND. 
     *       (MTHET0 .LE. THET0_NL(I+1))) THEN
	    KSZAM1=I
	    KSZAP1=I+1
          ENDIF
       ENDDO

C      Find THE indices
       DO I = 1,NLTHE-1
          IF ((MTHET .GE. THE_NL(I)) .AND. 
     *       (MTHET .LE. THE_NL(I+1))) THEN
	    KTHEM1=I
	    KTHEP1=I+1
          ENDIF
       ENDDO

C      Find PHI indices
       DO I = 1,NLPHI-1
          IF ((MDPHI .GE. PHI_NL(I)) .AND. 
     *       (MDPHI .LE. PHI_NL(I+1))) THEN
	    KPHIM1=I
	    KPHIP1=I+1
          ENDIF
       ENDDO


       RETURN
       END


			
C*********************************************************************
      SUBROUTINE   INTANGLE_NL(MTHET0,MTHET,MDPHI,
     *	  INT_NL1,Fd_NL1,T_NL1,SBAR_NL1,
     *    INT_NL,Fd_NL,T_NL,FdT_NL,SBAR_NL,
     *    KSZAM1,KSZAP1,KTHEM1,KTHEP1,KPHIM1,KPHIP1)
C----------------------------------------------------------------------
C!F77
C
C!DESCRIPTION:  Subroutine INTANGLE_NL interpolates the lookup
C               reflectances to the measured geometry.
C
C!INPUT PARAMETERS:
C	  INT_NL1		  radiance
C	  Fd_NL1		  flux down
C	  T_NL1		  transmission
C	  SBAR_NL1		  sbar
C         KSZAM1      Starting Index for solar zenith angle
C         KSZAP1      Ending Index for solar zenith angle
C         KTHEM1      Starting Index for view angle
C         KTHEP1      Ending   Index for  view angle
C         KPHIM1      Starting Index for  azimuth angle
C         KPHIP1      Ending   Index for  azimuth angle
C         MTHET0        Solar zenith angle.
C         MTHET        View angle.
C         MDPHI        Azimuth angle.
C
C         THET0      array of  solar Zenith angle in look_up table
C         THE        array of  satllite Zenith angle in look_up table
C         PHI        array of azimuth angle in look_up table
C
C!OUTPUT PARAMETERS 
C	  INT_NL		  interpolated radiance
C	  Fd_NL		  interpolated flux down
C	  T_NL		  interpolated transmission
C	  SBAR_NL         interpolated  sbar
C	  FdT_NL	  two way transmission transmission
C         REF_RAY        Reflectance for rayleigh only
C
C!REVISION HISTORY:
C
C!TEAM-UNIQUE HEADER:
C
C!REFERENCES AND CREDITS:
C
C!DESIGN NOTES:
C
C!END
C-----------------------------------------------------------------------

      IMPLICIT NONE
      SAVE

      INCLUDE 'mod04_land.inc'

c     Inputs
      REAL  MTHET0,MTHET,MDPHI
      INTEGER KSZAM1,KSZAP1,KTHEM1,KTHEP1,KPHIM1,KPHIP1
      
      REAL INT_NL1(NLPHI,NLTHE,NLTHET0,NLTAU,NLWAV,NLSIZE)
      REAL Fd_NL1(NLTHET0,NLTAU,NLWAV,NLSIZE)
      REAL SBAR_NL1(NLTHET0,NLTAU,NLWAV,NLSIZE)
      REAL T_NL1(NLTHE,NLTHET0,NLTAU,NLWAV,NLSIZE)

c     Outputs
      REAL INT_NL(NLTAU,NLWAV,NLSIZE),Fd_NL(NLTAU,NLWAV,NLSIZE)
      REAL T_NL(NLTAU,NLWAV,NLSIZE),FdT_NL(NLTAU,NLWAV,NLSIZE)
      REAL SBAR_NL(NLTAU,NLWAV,NLSIZE)

c     Dummy
      CHARACTER*132 LINE
      CHARACTER*45  LINE1
      INTEGER ICASE,IJ,IPHI,ISIZE,ITAU,ITH,ITH0,IWAV,IETA,LOPT
      REAL  X(2),Y(2),XX1(2),YY1(2),XX2(2),YY2(2),Y1
      REAL  WW2(2),WW1(2),W1
      REAL  VV2(2),V1
      REAL  UU2(2),U1
      INTEGER LL,MM,NN

C LOOP IS AROUND THET0,NEW FILE FOR EACH THETA0
C

        Y1=0.0
        DO IJ = 1,2
           X(IJ)=0.0
           Y(IJ)=0.0
           XX1(IJ)=0.0
           XX2(IJ)=0.0
           YY1(IJ)=0.0
           YY2(IJ)=0.0
           WW1(IJ)=0.0
           WW2(IJ)=0.0
           VV2(IJ)=0.0
           UU2(IJ)=0.0
        ENDDO
        DO 5 ISIZE= 1,NLSIZE
          DO 15 IWAV=1,NLWAV
            DO 25  ITAU  = 1,NLTAU
              LL=0
              y1 = 0.0
              DO  40 ITH0  = KSZAM1,KSZAP1
                MM=0
                DO  50  ITH  = KTHEM1,KTHEP1
                  NN=0
                  DO 60  IPHI  = KPHIM1,KPHIP1
                    NN=NN+1
                    X(NN)=PHI_NL(IPHI)
                    Y(NN)=INT_NL1(IPHI,ITH,ITH0,ITAU,IWAV,ISIZE)
60                CONTINUE

                  CALL INTERP_EXTRAP(NN,MDPHI,X,Y,Y1,1)
                  MM=MM+1
                  XX1(MM)=THE_NL(ITH)
                  YY1(MM)=Y1
                  WW1(MM)=T_NL1(ITH,ITH0,ITAU,IWAV,ISIZE)

50              CONTINUE
                y1=0.0
                w1=0.0 
                CALL INTERP_EXTRAP(MM,MTHET,XX1,YY1,Y1,1)
                CALL INTERP_EXTRAP(MM,MTHET,XX1,WW1,W1,1)
                LL=LL+1
                XX2(LL)=THET0_NL(ITH0)
                YY2(LL)=Y1
                WW2(LL)=W1
                VV2(LL)=Fd_NL1(ITH0,ITAU,IWAV,ISIZE)
                UU2(LL)=SBAR_NL1(ITH0,ITAU,IWAV,ISIZE)

40            CONTINUE
              y1=0.0
              w1=0.0
              v1=0.0
              u1=0.0
              CALL INTERP_EXTRAP(LL,MTHET0,XX2,YY2,Y1,1)
              CALL INTERP_EXTRAP(LL,MTHET0,XX2,WW2,W1,1)
              CALL INTERP_EXTRAP(LL,MTHET0,XX2,VV2,V1,1)
              CALL INTERP_EXTRAP(LL,MTHET0,XX2,UU2,U1,1)

              INT_NL(ITAU,IWAV,ISIZE) = Y1
              T_NL(ITAU,IWAV,ISIZE) = W1
              Fd_NL(ITAU,IWAV,ISIZE) = V1
              SBAR_NL(ITAU,IWAV,ISIZE) = U1
              FdT_NL(ITAU,IWAV,ISIZE) = W1*V1

25          CONTINUE
15        CONTINUE
 5      CONTINUE


        RETURN
        END



C*********************************************************************
      SUBROUTINE  TOA_SIMULATE(
     *    MTHET0,SCAT_ANGLE_LAND,
     *    INT_NL,FdT_NL,OPTH_NL,SBAR_NL,RHO_S212,RHOSTAR,
     *    yint_466,slope_466,yint_644,slope_644)

C-----------------------------------------------------------------------
C!F77
C
C
C!DESCRIPTION:
C               This subroutine simulates TOA reflectance
C
C
C! INPUTS
C!OUTPUTS (INTO INCLUDE FILE)
C	       RHOSTAR	  array of TOA reflectance
C-----------------------------------------------------------------------

        IMPLICIT NONE
        Save
        INCLUDE 'mod04_land.inc'
       
c     Inputs

      REAL MTHET0,SCAT_ANGLE_LAND
      REAL INT_NL(NLTAU,NLWAV,NLSIZE),FdT_NL(NLTAU,NLWAV,NLSIZE),
     *    SBAR_NL(NLTAU,NLWAV,NLSIZE),OPTH_NL(NLTAU,NLWAV,NLSIZE)
      REAL REF_RAY_NL(NLWAV)

c     for estimating surface reflectance
      REAL yint_466, slope_466
      REAL yint_644, slope_644
      REAL RHO_S466, RHO_S644
      REAL MVI

c     Outputs

      REAL  RHO_S212(NLTAU,NLSIZE), RHOSTAR(NLTAU,NLWAV,NLSIZE)

c     Dummy

        INTEGER IETA,ITAU,IWAV,ISIZE,I


C 	At each tau index, calculate surface reflectance at 2.1 and
C           apparent reflectance at all wavelengths
       DO ISIZE = 1, NLSIZE
          DO ITAU = 1,NLTAU
	    RHO_S212(ITAU,ISIZE) = 0.0
	    RHO_S212(ITAU,ISIZE) = 
     *         (INT_NL(ITAU,iwave_212,ISIZE)-REFW212L) / 
     *         ((SBAR_NL(ITAU,iwave_212,ISIZE) * 
     *          (INT_NL(ITAU,iwave_212,ISIZE)-REFW212L)) -  
     *          (FdT_NL(ITAU,iwave_212,ISIZE))) 

            IF (RHO_S212(ITAU,ISIZE) .GT. 1.0) THEN
              RHO_S212(ITAU,ISIZE) = -9.999
            ENDIF 


c         Estiamte surface reflectance
c         Use Red/IR and Blue/Red surface correlations

c         Use MVI (function of 1.24 and 2.1)

c MVI dependence
            MVI = (REFW124L - REFW212L)/(REFW124L + REFW212L)
            IF (MVI .GT. 0.75) THEN
               yint_644 = 0.0
               slope_644 = 0.58
               yint_466 = 0.005
               slope_466 = 0.49
            ELSE IF (MVI .LT. 0.25) THEN
               yint_644 = 0.0
               slope_644 = 0.48
               yint_466 = 0.005
               slope_466 = 0.49
            ELSE !(MVI around 0.5)
               yint_644 = 0.0
               slope_644 = 0.48 + (0.10/0.50)*(MVI-0.25)
               yint_466 = 0.005
               slope_466 = 0.49
            ENDIF

c Scattering angle dependence
            yint_644 = yint_644 - 0.00025*SCAT_ANGLE_LAND + 0.033
            slope_644 = slope_644 +  0.002*SCAT_ANGLE_LAND - 0.27

            RHO_S644 = slope_644*RHO_S212(ITAU,ISIZE) + yint_644
            RHO_S466 = slope_466*RHO_S644 + yint_466
 

C          Compute model differentiated apparent reflectance
           RHOSTAR(ITAU,iwave_466,ISIZE) = 0.0
           RHOSTAR(ITAU,iwave_466,ISIZE) = 
     *       INT_NL(ITAU,iwave_466,ISIZE) +
     *       ((FdT_NL(ITAU,iwave_466,ISIZE) * RHO_S466) / 
     *       (1 - (SBAR_NL(ITAU,iwave_466,ISIZE) * RHO_S466)))

           RHOSTAR(ITAU,iwave_553,ISIZE) = 0.0

           RHOSTAR(ITAU,iwave_644,ISIZE) = 0.0
           RHOSTAR(ITAU,iwave_644,ISIZE) = 
     *       INT_NL(ITAU,iwave_644,ISIZE) +
     *       ((FdT_NL(ITAU,iwave_644,ISIZE) * RHO_S644) / 
     *       (1 - (SBAR_NL(ITAU,iwave_644,ISIZE) * RHO_S644)))

           RHOSTAR(ITAU,iwave_212,ISIZE) = 0.0
           RHOSTAR(ITAU,iwave_212,ISIZE) = 
     *          INT_NL(ITAU,iwave_212,ISIZE) +
     *          ((FdT_NL(ITAU,iwave_212,ISIZE) *
     *           (RHO_S212(ITAU,ISIZE))) / 
     *           (1 - (SBAR_NL(ITAU,iwave_212,ISIZE) *
     *            (RHO_S212(ITAU,ISIZE)))))

          ENDDO

        ENDDO
        RETURN

      END


C*********************************************************************

      SUBROUTINE RETRIEVAL1(
     *   yint_466,slope_466,yint_644,slope_644,
     *   OPTH_NL,MASSCOEF_NL,EXTNORM_NL,RHO_S212,RHOSTAR,
     *   AVE_COUNT,ETA,ETA_FLAG,AOD,AODF,AODC,
     *   RHOSFC212,ERR644,RHOSFC,MASSCON)
C-----------------------------------------------------------------------
C!F77
C
C
C!DESCRIPTION:
C               This subroutine retrieves optical thickness, surface
C	 	reflectance and error parameters by comparing
C	        MODIS observations and LUT data
C
C
C! INPUTS
C              REFW466L,REFW644L,REFW212L   reflectance
C              OPTH_NL      LUT based optical thickness
C              MASSCOEF_NL  LUT based mass concentration coeficients
C              EXTNORM_NL   LUT based normalized extinction coeficients
C              RHO_S212   simulated surface reflectance
C              RHOSTAR    simulated TOA reflectance
C              yint644,yint466,slope466,slope644  
C                 characteristics of VIS/IR surface reflectance relationship
C! OUTPUTS
C	       ERR644	  retrieval fitting errors 
C	       AOD        retrieved Optical depths 
C	       AODF/C     Fine and coarse mode optical depth 
C	       RHOSFC     retrieved surface reflectance 
C	       ETA        retrieved fine mode weighting
C	       MASSCON    retrieved mass concentration coeficient
C              ETA_FLAG   0 if eta within 0.0 - 1.0
C
C-----------------------------------------------------------------------

        IMPLICIT NONE
        INCLUDE 'mod04_land.inc'

c     INPUTS
      REAL OPTH_NL(NLTAU,NLWAV,NLSIZE)
      REAL MASSCOEF_NL(NLTAU,NLWAV,NLSIZE)
      REAL EXTNORM_NL(NLTAU,NLWAV,NLSIZE)
      REAL RHO_S212(NLTAU,NLSIZE), RHOSTAR(NLTAU,NLWAV,NLSIZE)
      REAL yint_466,slope_466,yint_644,slope_644

c     Intermediate
      REAL ETA_TEMP
      REAL AOD553_TEMP(NLETA)
      REAL RHOSFC212_TEMP(NLETA),ERR644_TEMP(NLETA)
      REAL ERR644_SORT(NLETA)
      REAL NLETA_SORT(NLETA)
      REAL RHOSFC212_SORT(NLETA)
      REAL AOD553_SORT(NLETA)


c     OUTput TAU, ETA and SFCReflectance
      REAL SFC212OUT,ETAOUT,ERR644OUT,AOD553OUT
      REAL SFC212BEST,ETABEST,ERR644BEST
      REAL SFC212TOT,ETATOT,ERR644TOT
      REAL SFC212AVE,ETAAVE,ERR644AVE
      REAL AOD553AVE,AOD553BEST,AOD553TOT
      INTEGER AVE_COUNT

     
      REAL AODF553, AODC553, AOD553, EXTNORM
      REAL RHOSTAR_TOT(NLTAU,NLWAV)
      REAL ERRWAVE(NLWAV)
      REAL RHO_S212_TOT(NLTAU)
      REAL RHOSTAR_TOT1(NLWAV,NLETA)

c     OUTPUT
      REAL ETA,ERR644,RHOSFC(NLWAV),AOD(NLWAV)
      REAL RHOSFC212
      REAL MASSCON,MASSCONF,MASSCONC
      REAL AODF(NLWAV), AODC(NLWAV)
      INTEGER ETA_FLAG

c     Dummy

      REAL XMIN,XS(100,100),YS(100,100),X(100),Y(100),Denom
      REAL y1
      INTEGER IETA,JETA,ITAU,IWAV,ISIZE,I,NSOLUTION

      SAVE	

c     Initialize

      ETA = -9999
      ERR644 = -9999
      MASSCON = -9999
      DO IWAV = 1, NLWAV
        RHOSFC(NLWAV) = -9999
        AOD(NLWAV) = -9999
      ENDDO
   

C Now loop through ETA values (kind of like ocean)
      NSOLUTION = 0
      DO 199 IETA = 1, NLETA

C Average, best, QCONTROL, FILLVALUE, OUTPUT, etc
        ETA_TEMP = -0.2 + (IETA * 0.1)

C 	Compute total apparent reflectance

        DO ITAU = 1, NLTAU
          DO IWAV = 1, NLWAV
	     RHOSTAR_TOT(ITAU,IWAV) = 
     *          ETA_TEMP * RHOSTAR(ITAU,IWAV,1) + 
     *          ((1-ETA_TEMP) * RHOSTAR(ITAU,IWAV,2))
          ENDDO
C 	Compute total surface reflectance
          RHO_S212_TOT(ITAU) = 
     *          ETA_TEMP * RHO_S212(ITAU,1) + 
     *          ((1-ETA_TEMP) * RHO_S212(ITAU,2))


        ENDDO


C                   ********COMPUTE FOR TAU FOR  USING
C                           RADAINCE VALUE OF 0.466 UM AND
C                           TAU VALUES OF WAV553
C

	 DO ITAU = 1, NLTAU
	     X(ITAU) = RHOSTAR_TOT(ITAU,iwave_466)
	     Y(ITAU) = OPTH_NL(ITAU,iwave_553,1)
          ENDDO
          Y1=0
          CALL INTERP_EXTRAP(NLTAU,REFW466L,X,Y,Y1,1)
          AOD553_TEMP(IETA) = y1

C
C                 ******* FOR TAUX55 COMPUTE appararent Reflectance FOR
C                     all wavelengths
C
         DO IWAV = 1,NLWAV
            DO ITAU = 1,NLTAU
               X(ITAU)=OPTH_NL(ITAU,iwave_553,1)
               Y(ITAU)=RHOSTAR_TOT(ITAU,IWAV)
	    ENDDO
            y1=0
            CALL INTERP_EXTRAP(NLTAU,AOD553_TEMP(IETA),X,Y,Y1,1)
            RHOSTAR_TOT1(IWAV,IETA)=Y1
         ENDDO

C
C                 ******* FOR TAUX55 COMPUTE Surface Reflectance at 2.1
C
          DO ITAU = 1,NLTAU
               X(ITAU)=OPTH_NL(ITAU,iwave_553,1)
               Y(ITAU)=RHO_S212_TOT(ITAU)
	  ENDDO
          y1=0
          CALL INTERP_EXTRAP(NLTAU,AOD553_TEMP(IETA),X,Y,Y1,1)
          RHOSFC212_TEMP(IETA)=Y1

C   	Errors

	 ERRWAVE(1) = ABS(REFW466L - RHOSTAR_TOT1(1,IETA)) / REFW466L
	 ERRWAVE(3) = ABS(REFW644L - RHOSTAR_TOT1(3,IETA)) / REFW644L
	 ERRWAVE(4) = ABS(REFW212L - RHOSTAR_TOT1(4,IETA)) / REFW212L

C
C                 ******* Determine error at 644nm  
C
C

 	 ERR644_TEMP(IETA) = ERRWAVE(iwave_644)


 199      ENDDO

C Sort results in ascending error order
	  DO IETA = 1, NLETA 
	     NLETA_SORT(IETA) = IETA
             ERR644_SORT(IETA) = ERR644_TEMP(IETA)
             RHOSFC212_SORT(IETA) = RHOSFC212_TEMP(IETA)
             AOD553_SORT(IETA) = AOD553_TEMP(IETA)
	  ENDDO
	  CALL SHLSRT_NL(NLETA,ERR644_SORT,NLETA_SORT)
	  DO IETA = 1, NLETA 
             JETA = NLETA_SORT(IETA)
             RHOSFC212_SORT(IETA) = RHOSFC212_SORT(JETA)
             AOD553_SORT(IETA) = AOD553_SORT(JETA)
	  ENDDO

C Best solution


          ETABEST = (-0.2 + (NLETA_SORT(1) * 0.1))
          ERR644BEST = ERR644_SORT(1)
          SFC212BEST = RHOSFC212_SORT(1)
          AOD553BEST = AOD553_SORT(1)
          
C Average solution
        
          AOD553TOT = 0.0
          AOD553AVE = -9999
          ETATOT = 0.0
          ERR644TOT = 0.0
          SFC212TOT = 0.0
          ETAAVE = -9999
          ERR644AVE = -9999
          SFC212AVE = -9999
          AVE_COUNT = 0
  
          DO IETA = 1, NLETA
             IF (ERR644_SORT(IETA) .LE. 0.03 
     *            .AND. AVE_COUNT .LT. 3) THEN
                ERR644OUT = ERR644_SORT(IETA)
                ETAOUT = (-0.2 + (NLETA_SORT(IETA) * 0.1))
                SFC212OUT = RHOSFC212_SORT(IETA)
                AOD553OUT = AOD553_SORT(IETA)
                ETATOT = ETATOT + ETAOUT 
                ERR644TOT = ERR644TOT + ERR644OUT 
                SFC212TOT = SFC212TOT + SFC212OUT
                AOD553TOT = AOD553TOT + AOD553OUT 
                AVE_COUNT = AVE_COUNT + 1
             ENDIF
          ENDDO
          
          IF (AVE_COUNT .GT. 0) THEN
             ETAAVE = ETATOT / AVE_COUNT
             ERR644AVE = ERR644TOT / AVE_COUNT
             SFC212AVE = SFC212TOT / AVE_COUNT
             AOD553AVE = AOD553TOT / AVE_COUNT
          ENDIF

c          IF (ETAAVE .NE. -9999) THEN
c             ETA = ETAAVE
c             ERR644 = ERR644AVE
c             RHOSFC212 = SFC212AVE
c             AOD553 = AOD553AVE
c             RHOSFC(iwave_212) = RHOSFC212
c             RHOSFC(iwave_466) = yint_466 + slope_466*RHOSFC212
c             RHOSFC(iwave_644) = yint_644 + slope_644*RHOSFC212
c          ELSE
          AVE_COUNT = -1
          ETA = ETABEST
          ERR644 = ERR644BEST
          RHOSFC212 = SFC212BEST
          AOD553 = AOD553BEST
          RHOSFC(iwave_212) = RHOSFC212
          RHOSFC(iwave_644) = yint_644 + slope_644*RHOSFC212 
          RHOSFC(iwave_466) = yint_466 + slope_466*RHOSFC(iwave_644) 
c          ENDIF

c set  out of bounds Eta to 0 or 1

       IF( ETA .LT. 0) THEN
           ETA=0
           ETA_FLAG = -1
       ENDIF
       IF( ETA .GT. 1) THEN
           ETA=1
           ETA_FLAG = -1
       ENDIF


c   Compute fine and coarse AOD for mass concentration
           AODF553 = AOD553 * ETA
           AODC553 = AOD553 * (1.-ETA)
C
C                 ******* FOR TAUX55 COMPUTE Mass Concentration
C

C   Fine
          DO ITAU = 1,NLTAU
               X(ITAU)=OPTH_NL(ITAU,iwave_553,1)
               Y(ITAU)=MASSCOEF_NL(ITAU,iwave_553,1)
	  ENDDO
          y1=0
          CALL INTERP_EXTRAP(NLTAU,AODF553,X,Y,Y1,1)
          MASSCONF=Y1*AODF553

C  Coarse 
          DO ITAU = 1,NLTAU
               X(ITAU)=OPTH_NL(ITAU,iwave_553,2)
               Y(ITAU)=MASSCOEF_NL(ITAU,iwave_553,2)
	  ENDDO
          y1=0
          CALL INTERP_EXTRAP(NLTAU,AODC553,X,Y,Y1,1)
          MASSCONC=Y1*AODC553

          MASSCON = MASSCONF + MASSCONC

C Compute AOD for fine mode at all wavelengths
         DO IWAV = 1,NLWAV
            DO ITAU = 1,NLTAU
               X(ITAU)=OPTH_NL(ITAU,iwave_553,1)
               Y(ITAU)=EXTNORM_NL(ITAU,IWAV,1)
            ENDDO
            y1=0
            CALL INTERP_EXTRAP(NLTAU,AODF553,X,Y,Y1,1)
            EXTNORM=Y1
            AODF(IWAV) = AODF553 * EXTNORM
         ENDDO

C Compute AOD for coarse mode at all wavelengths
C Compute total AOD at all wavelengths
         DO IWAV = 1,NLWAV
            DO ITAU = 1,NLTAU
               X(ITAU)=OPTH_NL(ITAU,iwave_553,2)
               Y(ITAU)=EXTNORM_NL(ITAU,IWAV,2)
            ENDDO
            y1=0
            CALL INTERP_EXTRAP(NLTAU,AODC553,X,Y,Y1,1)
            EXTNORM=Y1
            AODC(IWAV) = AODC553 * EXTNORM
            AOD(IWAV) = AODF(IWAV) + AODC(IWAV)
         ENDDO

         PRINT*, 'AOD ', (AOD(IWAV),IWAV=1,4)
         PRINT*, 'AODF ', (AODF(IWAV),IWAV=1,4)
         PRINT*, 'AODC ', (AODC(IWAV),IWAV=1,4)

      RETURN
      END

C*********************************************************************

      SUBROUTINE RETRIEVAL2(
     *   yint_466,slope_466,yint_644,slope_644,
     *   OPTH_NL,MASSCOEF_NL,RHO_S212,RHOSTAR,
     *   AVE_COUNT,ETA,AOD,RHOSFC212,ERR644,RHOSFC,MASSCON)
C-----------------------------------------------------------------------
C!F77
C
C
C!DESCRIPTION:
C               This subroutine retrieves optical thickness, surface
C	 	reflectance and error parameters by comparing
C	        MODIS observations and LUT data. The "continental"
C               aerosol model is assumed.
C
C
C! INPUTS
C              REFW466L,REFW644L,REFW212L   reflectance
C              OPTH       array of LUT optical thickness
C              RHO_S212   simulated surface reflectance
C              RHOSTAR    simulated TOA reflectance
C              yint644,yint466,slope466,slope644  
C                 characteristics of VIS/IR surface reflectance relationship
C! OUTPUTS
C	       ERR644	  retrieval fitting errors 
C	       AOD        retrieved optical depths 
C	       RHOSFC     retrieved surface reflectance 
C	       ETA        retrieved fine mode weighting
C
C-----------------------------------------------------------------------

        IMPLICIT NONE
        INCLUDE 'mod04_land.inc'

c     INPUTS
      REAL OPTH_NL(NLTAU,NLWAV,NLSIZE)
      REAL MASSCOEF_NL(NLTAU,NLWAV,NLSIZE)
      REAL RHO_S212(NLTAU,NLSIZE), RHOSTAR(NLTAU,NLWAV,NLSIZE)
      REAL yint_466,slope_466,yint_644,slope_644

c     Intermediate

      INTEGER AVE_COUNT

     
      REAL RHOSTAR_TOT(NLTAU,NLWAV)
      REAL ERRWAVE(NLWAV)
      REAL RHO_S212_TOT(NLTAU)
      REAL RHOSTAR_TOT1(NLWAV)

c     OUTPUT
      REAL ETA,ERR644,RHOSFC(NLWAV),AOD(NLWAV)
      REAL RHOSFC212
      REAL MASSCON

c     Dummy

      REAL XMIN,XS(100,100),YS(100,100),X(100),Y(100),Denom
      REAL y1
      INTEGER JETA,ITAU,IWAV,ISIZE,I,NSOLUTION

      SAVE	

c     Initialize

      ETA = -9999
      ERR644 = -9999
      MASSCON = -9999
      DO IWAV = 1, NLWAV
        RHOSFC(NLWAV) = -9999
        AOD(NLWAV) = -9999
      ENDDO

C Now loop through ETA values (kind of like ocean)
      NSOLUTION = 0
C Average, best, QCONTROL, FILLVALUE, OUTPUT, etc

C 	Compute total apparent reflectance

        DO ITAU = 1, NLTAU
          DO IWAV = 1, NLWAV
	     RHOSTAR_TOT(ITAU,IWAV) = RHOSTAR(ITAU,IWAV,1)
          ENDDO
C 	Compute total surface reflectance
          RHO_S212_TOT(ITAU) =  RHO_S212(ITAU,1)
        ENDDO


C       Interpolate everything based on measured reflectance at 466.

C                   ********COMPUTE FOR TAU FOR  USING
C                           RADAINCE VALUE OF 0.466 UM AND
C                           TAU VALUES OF WAV553
C

        DO IWAV = 1, NLWAV
	  DO ITAU = 1, NLTAU
	     X(ITAU) = RHOSTAR_TOT(ITAU,iwave_466)
	     Y(ITAU) = OPTH_NL(ITAU,IWAV,1)
	  ENDDO
          Y1=0
          CALL INTERP_EXTRAP(NLTAU,REFW466L,X,Y,Y1,1)
          AOD(IWAV)=y1
        ENDDO

C
C                 ******* FOR TAUX55 COMPUTE appararent Reflectance FOR
C                     all wavelengths
C
         DO IWAV = 1,NLWAV
            DO ITAU = 1,NLTAU
               X(ITAU)=OPTH_NL(ITAU,iwave_553,1)
               Y(ITAU)=RHOSTAR_TOT(ITAU,IWAV)
	    ENDDO
            y1=0
            CALL INTERP_EXTRAP(NLTAU,AOD(iwave_553),X,Y,Y1,1)
            RHOSTAR_TOT1(IWAV)=Y1
         ENDDO

C
C                 ******* FOR TAUX55 COMPUTE Surface Reflectance at 2.1
C
          DO ITAU = 1,NLTAU
               X(ITAU)=OPTH_NL(ITAU,iwave_553,1)
               Y(ITAU)=RHO_S212_TOT(ITAU)
	  ENDDO
          y1=0
          CALL INTERP_EXTRAP(NLTAU,AOD(iwave_553),X,Y,Y1,1)
          RHOSFC212=Y1

C   	Errors

	 ERRWAVE(1) = ABS(REFW466L - RHOSTAR_TOT1(1)) / REFW466L
	 ERRWAVE(3) = ABS(REFW644L - RHOSTAR_TOT1(3)) / REFW644L
	 ERRWAVE(4) = ABS(REFW212L - RHOSTAR_TOT1(4)) / REFW212L

C
C                 ******* Determine error at 644nm  
C
C

 	 ERR644 = ERRWAVE(iwave_644)

         AVE_COUNT = -1
  
         RHOSFC(iwave_212) = RHOSFC212
         RHOSFC(iwave_553) = -9999
         RHOSFC(iwave_644) = yint_644 + slope_644*RHOSFC212 
         RHOSFC(iwave_466) = yint_466 + slope_466*RHOSFC(iwave_644)

C
C                 ******* FOR TAUX55 COMPUTE Mass Concentration  
C

C   Continental only 
          DO ITAU = 1,NLTAU
               X(ITAU)=OPTH_NL(ITAU,iwave_553,1)
               Y(ITAU)=MASSCOEF_NL(ITAU,iwave_553,1)
          ENDDO
          y1=0
          CALL INTERP_EXTRAP(NLTAU,AOD(iwave_553),X,Y,Y1,1)
          MASSCON=Y1*AOD(iwave_553)


      RETURN
      END

C*********************************************************************

       subroutine SHLSRT_NL(n,arr1,arr2)

C-----------------------------------------------------------------------
C !F77
C
C !DESCRIPTION: From numerical recepies Sorts arr1ay 'arr1' of
C              length 'n' into ascending order ('arr1' and arr2
C              is replaced on OUTPUT by its sorted rearr1angement)
C
C !INPUT PARAMETERS:
C                 N         length of ARR1
C                ARR1       Unsorted array ARR1
C                ARR2       Unsorted array ARR2
C
C !OUTPUT PARAMETERS:
C                ARR1       Sorted array ARR1
C                ARR2       Sorted array ARR2
C
C
C !TEAM-UNIQUE HEADER:
C
C !END
C-----------------------------------------------------------------------

      implicit none

      integer i,j,k,l,lognb2,m,n,nn
      real aln2i,tiny
      parameter (aln2i = 1.4426950, tiny = 1.e-5)
      real  t,q,arr1(n),arr2(n)

*/  Modified by JC Guu  01/14/97
*/  One logical variable declared.

      logical flag

      lognb2 = int(alog(float(n))*aln2i+tiny)
      m = n
      do 12 nn = 1,lognb2
        m = m/2
        k = n-m
        do 11 j = 1,k
          i = j

*/  Modified by JC Guu  01/14/97
*/  Change a "continue ... if() go to" construct
*/  to a "DO WHILE" construct.

C  3       continue

          flag=.true.

          do while(flag.eqv..true.)
          l = i+m
          flag=.false.
          if(arr1(l).lt.arr1(i)) then
            t = arr1(i)
            arr1(i) = arr1(l)
            arr1(l) = t
            q = arr2(i)
            arr2(i) = arr2(l)
            arr2(l) = q
c            write(6,*)'nn,m,k,j,i,l,t,arr1(i),arr1(l)',
c     * nn,m,k,j,i,l,t,arr1(i),arr1(l)
            i = i-m
            if(i.ge.1) flag=.true.
          endif
          end do

C            if(i.ge.1) go to 3
C          endif

11      continue
12    continue

      return

      end

C*********************************************************************

      SUBROUTINE INTERP_EXTRAP(M,X1,X,Y,Y1,INDEX)

C----------------------------------------------------------------------
C !F77
C
C !DESCRIPTION: This subroutine is a general purpose routine and
C              interpolates linearly.  Value of y1 is interpolated
C              or extrapolated for x1
C
C !INPUT PARAMETERS:I,M,X1,X,Y
C
C !OUTPUT PARAMETERS:Y1,LOPT
C
C !REVISION HISTORY:
C $Log: Process_ocean_V2.f,v $
c
C
C !TEAM-UNIQUE HEADER:
C
C !END
C----------------------------------------------------------------------

      IMPLICIT NONE

      INTEGER IL,LL,LOPT,M,INDEX
      REAL PINTEN,PPHI,SINTEN,SPHI
      REAL  X(M),Y(M),Y1,X1,Diff

      Y1=0.0
      LOPT=0
      LL=M-1
        DO 230 IL=1,LL
C        Extrapolation on lower bound
       IF(X1 .LE.X(1))THEN
          PPHI=X(1)
          SPHI=X(2)
          PINTEN=Y(1)
          SINTEN=Y(2)
           Diff=(SPHI-PPHI)
           if(Diff . Le. 0.0) Diff=1
          Y1=PINTEN+((SINTEN-PINTEN)*((X1-PPHI)/Diff))
               LOPT=1

*/  Modified by JC Guu  01/09/97
*/  "GO TO 290" is changed to RETURN

           RETURN

C        Extrapolation on upper bound
       ELSEIF(X1 .GE.X(LL+1)) THEN
           PPHI=X(LL)
           SPHI=X(LL+1)
         PINTEN=Y(LL)
         SINTEN=Y(LL+1)
         Diff=(SPHI-PPHI)
           if(Diff .Le. 0.0) Diff=1
          Y1=PINTEN+((SINTEN-PINTEN)*((X1-PPHI)/Diff))
              LOPT=1

*/  Modified by JC Guu  01/09/97
*/  "GO TO 290" is changed to RETURN

           RETURN

C       interpolation
       ELSEIF (X1 .GE.X(IL) .AND.X1.LE. X(IL+1)) THEN
         PPHI=X(IL)
         SPHI=X(IL+1)
         PINTEN=Y(IL)
         SINTEN=Y(IL+1)
         Diff=(SPHI-PPHI)
           if(Diff .Le. 0.0) Diff=1
          Y1=PINTEN+((SINTEN-PINTEN)*((X1-PPHI)/Diff))
          LOPT=1

*/  Modified by JC Guu  01/09/97
*/  "GO TO 290" is changed to RETURN
*/  and two lines are commented out.

           RETURN
C        ELSE
C        GO TO 230

         ENDIF
  230    CONTINUE

  290    RETURN
          END



C***********************************************************************
      SUBROUTINE RNLOOKUP(
     +    HANDLE_LUT466,HANDLE_LUT553,HANDLE_LUT644,HANDLE_LUT213,
     +    INT_NL0,Fd_NL0,T_NL0,OPTH_NL0,
     +    SBAR_NL0,MASSCOEF_NL0,EXTNORM_NL0)

C-----------------------------------------------------------------------
C!F77
C
C
C!DESCRIPTION:
C               This subroutine reads 4 lookup tables
C
C!INPUT PARAMETERS:
C
C!OUTPUT PARAMETERS:
C
C          INT_NL0           Reflectance
C          Fd_NL0           Flux down
C          T_NL0           Transmission factor
C          OPTH_NL0          Optical depth
C          SBAR_NL0          Sbar
C          MASSCOEF_NL0      Mass Concentration coeficient
C          EXTNORM_NL0       Normalized extinction coeficient
C
C!OUTPUT PARAMETERS (INTO INCLUDE FILE)
C          THET0           Array of LUT Solar Zenith Angles
C          THE             Array of LUT Satellite View angles
C          PHI             Array of LUT Azimuth angle differences
C          WAV             Array of LUT wavelengths
C
C!DESIGN NOTES:
C
C  Read from look-up table.
C   wav=wavelengths(2)(.466,.533,.644 and 2.13 micrometers)  4
C   lthet0=number of theta0(solar zenith angle)              9
C   thet0=(0,12,24,36,48,54,60,66,72   degrees)
C   lthe =number of the(observation zenith angle)           16
C   the=(0,6,.....72 degrees every 6 degrees)
C   lphi =number of phi(observation azimuth angle)          31
C   phi=(0,6,12,18,24,30,36,42,48,54,60,66,72,78,84,90,96,
C        102,108,114,120,126,132,138,144,150,156,162,168,
C        174,180)
C   ltau =number of opth(optical thickness)                  7
C   opth=(0.0,0.25,0.50,1.0,2.0,3.0,5.0)
C   lhght = number of hght(observation heights)              1
C   hght=(80.0 km)
C
C!REVISION HISTORY:
C 02/08/2006 Rob Levy/Shana Mattoo
C!TEAM-UNIQUE HEADER:
C
C Developed by MODIS aerosol team at NASA Goddard Space Flight
C Center, Greenbelt, Maryland
C
C!END
C-------------------------------------------------------------------------
C
        IMPLICIT NONE 
        INCLUDE 'mod04_land.inc'
        INTEGER HANDLE_LUT466,HANDLE_LUT553
        INTEGER HANDLE_LUT644,HANDLE_LUT213
        INTEGER I,IPHI,ITAU,ITAB,NFILE
        INTEGER ITHE,ITHET0,IWAV
        REAL OPTH_NL0(NLTAU,NLWAV,NLTABLE)
        REAL MASSCOEF_NL0(NLTAU,NLWAV,NLTABLE)
        REAL EXTNORM_NL0(NLTAU,NLWAV,NLTABLE)
        REAL SSA_NL0(NLTAU,NLWAV,NLTABLE)
        REAL QEXT_NL0(NLTAU,NLWAV,NLTABLE)
        REAL BEXT_NL0(NLTAU,NLWAV,NLTABLE)
        REAL VEXT_NL0(NLTAU,NLWAV,NLTABLE)
        REAL MEXT_NL0(NLTAU,NLWAV,NLTABLE)
        REAL SBAR_NL0(NLTHET0,NLTAU,NLWAV,NLTABLE)
        REAL INT_NL0(NLPHI,NLTHE,NLTHET0,NLTAU,NLWAV,NLTABLE)
        REAL Fd_NL0(NLTHET0,NLTAU,NLWAV,NLTABLE)
        REAL T_NL0(NLTHE,NLTHET0,NLTAU,NLWAV,NLTABLE)
        CHARACTER*1 LINE(132)
        REAL  Omega0(NLTAU,NLWAV,NLTABLE),ROD(NLWAV),GOD(NLWAV)
        SAVE

C Loops through each wavelength and aerosol model

      DO 10 IWAV = 1,NLWAV
      DO 20 ITAB=1,NLTABLE

        IF (IWAV .EQ. 1) NFILE = HANDLE_LUT466        
        IF (IWAV .EQ. 2) NFILE = HANDLE_LUT553
        IF (IWAV .EQ. 3) NFILE = HANDLE_LUT644
        IF (IWAV .EQ. 4) NFILE = HANDLE_LUT213
C
C Reads observation zenith angles(the),and observation azimuth angles
C from look-up tables
C
        READ(NFILE,1111)(THE_NL(I),I=1,NLTHE)
        READ(NFILE,1000)(PHI_NL(I),I=1,NLPHI)

C
C Reads wavelength(wav),optical thickness(opth),solar zenith angle(thet0),
C reflectance ofatmosphere(sbarw) from look-up tables.
C

        DO 100 ITAU=1,NLTAU
          READ(NFILE,*)
          READ(NFILE,1001) SSA_NL0(ITAU,IWAV,ITAB),
     +        QEXT_NL0(ITAU,IWAV,ITAB),
     +        BEXT_NL0(ITAU,IWAV,ITAB),
     +        VEXT_NL0(ITAU,IWAV,ITAB),
     +        MEXT_NL0(ITAU,IWAV,ITAB),
     +        MASSCOEF_NL0(ITAU,IWAV,ITAB)
          READ(NFILE,1002) WAV_NL(IWAV),OPTH_NL0(ITAU,IWAV,ITAB),
     +        ROD(IWAV),GOD(IWAV)
 
          DO 101 ITHET0=1,NLTHET0
             READ(NFILE,1003) 
     +         THET0_NL(ITHET0), SBAR_NL0(ITHET0,ITAU,IWAV,ITAB), 
     +         Fd_NL0(ITHET0,ITAU,IWAV,ITAB)

C
C Reads transmission as a function of observation zenith angle,
C optical thickness
C

             READ(NFILE,1004)
     +          (T_NL0(ITHE,ITHET0,ITAU,IWAV,ITAB),ITHE=1,NLTHE)       
             READ(NFILE,*)
C
C Reads atmospheric radiance (int) as a function of solar zenith angle,
C optical thickness, height, observation zenith angle and azimuth angle
C from look-up table.
C
              DO 103 ITHE=1,NLTHE
                 READ(NFILE,1005)
     +      (INT_NL0(IPHI,ITHE,ITHET0,ITAU,IWAV,ITAB),IPHI=1,NLPHI)
103           CONTINUE
101       CONTINUE
100     CONTINUE

c  Set extinction parameters for "AOD = 0.0" case

        QEXT_NL0(1,IWAV,ITAB) = QEXT_NL0(2,IWAV,ITAB)
        BEXT_NL0(1,IWAV,ITAB) = BEXT_NL0(2,IWAV,ITAB)
        VEXT_NL0(1,IWAV,ITAB) = VEXT_NL0(2,IWAV,ITAB)
        MEXT_NL0(1,IWAV,ITAB) = MEXT_NL0(2,IWAV,ITAB)
        MASSCOEF_NL0(1,IWAV,ITAB) = MASSCOEF_NL0(2,IWAV,ITAB)

20    CONTINUE
10    CONTINUE

       DO ITAB=1,NLTABLE
       DO ITAU=1,NLTAU
       DO IWAV=1,NLWAV
       EXTNORM_NL0(ITAU,IWAV,ITAB) =
     *     QEXT_NL0(ITAU,IWAV,ITAB) / QEXT_NL0(ITAU,iwave_553,ITAB)
       ENDDO
       ENDDO
       ENDDO

1000  FORMAT(5X,16F8.3)
1111  FORMAT(5X,15F8.3)
1001  FORMAT(6(5x,f9.4))
1002  FORMAT(10x,F9.4,3(4x,f9.4))
1003  FORMAT(5x,f8.3, 5x, f8.5, 2(5x,e11.4))
1004  FORMAT(3X,15e11.5)
1005  FORMAT(16e11.5)

      RETURN
      END


C*********************************************************************
      SUBROUTINE INT_ELEV(EQWAV_NL,INT_NL,Fd_NL,T_NL,FdT_NL,OPTH_NL,
     *        SBAR_NL,REF_RAY_NL,MHGHT)
C----------------------------------------------------------------------
C!F77
C
C!DESCRIPTION:  Subroutine INTELEV interpolates the lookup
C               reflectances to the target elevation.
C               Basically it is fudge:
C               Interpolate between wavelengths to simulate
C               elevation by a longer wavelength
C
C!INPUT PARAMETERS:
C         WAV_NL           wavelengths
C	  INT_NL		  radiance
C	  Fd_NL		  flux down
C         T_NL            transmission
C         FdT_NL          2 way transmission factor
C	  SBAR_NL	  sbar
C         OPTH_NL          optical depth
C         MGHTH           elevation
C
C!OUTPUT PARAMETERS 
C	  INT_NL		  interpolated radiance
C	  Fd_NL		          interpolated flux down
C         T_NL                    interpolated transmission
C	  FdT_NL		  interpolated transmission
C	  SBAR_NL	  interpolated sbar
C         OPTH_NL          interpolated optical depth
C         REF_RAY        Reflectance for rayleigh only
C         EQWAV_NL        Interpolated wavelengths
C
C!REVISION HISTORY:
C 02/08/2006 Shana Mattoo/Rob Levy
C Initial revision
C
C!TEAM-UNIQUE HEADER:
C
C Developed by MODIS Aerosol Retrieval Team at NASA GSFC, Greenbelt, MD
C
C!REFERENCES AND CREDITS
C
C!DESIGN NOTES:
C
C!END
C-----------------------------------------------------------------------

      IMPLICIT NONE
      SAVE

      INCLUDE 'mod04_land.inc'

c     Inputs
      REAL  MHGHT
      
c     Inputs and Outputs
      REAL INT_NL(NLTAU,NLWAV,NLSIZE),FdT_NL(NLTAU,NLWAV,NLSIZE),
     *    Fd_NL(NLTAU,NLWAV,NLSIZE), T_NL(NLTAU,NLWAV,NLSIZE),
     *    SBAR_NL(NLTAU,NLWAV,NLSIZE),OPTH_NL(NLTAU,NLWAV,NLSIZE)
      REAL REF_RAY_NL(NLWAV)
      REAL INT_NL9(NLTAU,NLWAV,NLSIZE),FdT_NL9(NLTAU,NLWAV,NLSIZE),
     *    Fd_NL9(NLTAU,NLWAV,NLSIZE), T_NL9(NLTAU,NLWAV,NLSIZE),
     *    SBAR_NL9(NLTAU,NLWAV,NLSIZE),OPTH_NL9(NLTAU,NLWAV,NLSIZE)

      REAL ROD_1013(NLWAV)
      REAL p, p0, expfactor
      PARAMETER (p0 = 1013.0)
      REAL ROD_PRES(NLWAV)
      REAL EQWAV_NL(NLWAV)

      REAL lambda0, lambda1, lambda2, diff0, exp0, exp1, exp2
      REAL tau0, tau1, tau2

c     Dummy
      INTEGER ICASE,IJ,ISIZE,ITAU,IWAV,IWAV1,IWAV2,JWAV
      REAL  X(2),Y(2),W(2),V(2),Z(2),T(2),U(2)
      REAL  Y1,W1,V1,Z1,T1,U1
      INTEGER LL


C     Estimate surface pressure (hypsometric EQ)
      p = p0 * exp(-(MHGHT/7.5))

C     Estimate ROD at nominal wavelenghts at p0 and at pres
      DO IWAV = 1, 3
         expfactor = -4.15 + (0.2 * WAV_NL(IWAV))
         ROD_1013(IWAV) = 0.0088 * WAV_NL(IWAV)**expfactor
         ROD_PRES(IWAV) = 0.0088*(p/p0) * WAV_NL(IWAV)**expfactor

C    Estimate equivalent wavelengths for ROD at pressure
         lambda1 = 0.1
         lambda2 = 2.0
         diff0 = 99.
         DO WHILE (diff0 .GT. 0.00001)
           lambda0 = (lambda1 + lambda2) / 2.
           exp0 = -4.15 + 0.2*lambda0
           exp1 = -4.15 + 0.2*lambda1
           exp2 = -4.15 + 0.2*lambda2
           tau0 = 0.0088*lambda0**exp0
           tau1 = 0.0088*lambda1**exp1
           tau2 = 0.0088*lambda2**exp2
           IF (tau1 .GT. ROD_PRES(IWAV) .AND. 
     *          tau2 .LT. ROD_PRES(IWAV)) THEN
              IF (tau0 .GT. ROD_PRES(IWAV)) THEN
                lambda1 = (lambda1 + lambda2)/2.
              ELSE 
                lambda2 = (lambda1 + lambda2)/2.
              ENDIF
           ENDIF
           diff0 = ABS(tau0 - ROD_PRES(IWAV))
         ENDDO
         EQWAV_NL(IWAV) = lambda0 
       ENDDO

 
C     INterpolate lookup tables to equiv Waves (let's start only with
C      1st two wavelengths until we derive 0.86 table)


        DO IJ = 1,2
           X(IJ)=0.0
           Y(IJ)=0.0
           V(IJ)=0.0
           W(IJ)=0.0
           Z(IJ)=0.0
           U(IJ)=0.0
           T(IJ)=0.0
        ENDDO
        DO 5 ISIZE= 1,NLSIZE
          DO 15 IWAV=1,3
            DO 25  ITAU  = 1,NLTAU
              LL=0
              Y1=0.0
              W1=0.0
              Z1=0.0
              V1=0.0
              IF (IWAV .EQ. 3) THEN
               IWAV1 = IWAV-1
               IWAV2 = IWAV
              ELSE
               IWAV1 = IWAV
               IWAV2 = IWAV+1
              ENDIF
              DO 60  JWAV = IWAV1, IWAV2
                LL=LL+1

c    Interpolate on log log scale
 
                X(LL)=ALOG(WAV_NL(JWAV))
                Y(LL)=ALOG(INT_NL(ITAU,JWAV,ISIZE))
                W(LL)=ALOG(FdT_NL(ITAU,JWAV,ISIZE))
                Z(LL)=ALOG(SBAR_NL(ITAU,JWAV,ISIZE))
                V(LL)=OPTH_NL(ITAU,JWAV,ISIZE)
                U(LL)=ALOG(Fd_NL(ITAU,JWAV,ISIZE))
                T(LL)=ALOG(T_NL(ITAU,JWAV,ISIZE))
                IF (OPTH_NL(ITAU,JWAV,ISIZE) .GT. 0.) THEN
                  V(LL)=ALOG(OPTH_NL(ITAU,JWAV,ISIZE))
                ENDIF
60            CONTINUE

              CALL INTERP_EXTRAP(LL,ALOG(EQWAV_NL(IWAV)),X,Y,Y1,1)
              CALL INTERP_EXTRAP(LL,ALOG(EQWAV_NL(IWAV)),X,W,W1,1)
              CALL INTERP_EXTRAP(LL,ALOG(EQWAV_NL(IWAV)),X,Z,Z1,1)
              CALL INTERP_EXTRAP(LL,ALOG(EQWAV_NL(IWAV)),X,V,V1,1)
              CALL INTERP_EXTRAP(LL,ALOG(EQWAV_NL(IWAV)),X,U,U1,1)
              CALL INTERP_EXTRAP(LL,ALOG(EQWAV_NL(IWAV)),X,T,T1,1)

              INT_NL9(ITAU,IWAV,ISIZE) = EXP(Y1)
              FdT_NL9(ITAU,IWAV,ISIZE) = EXP(W1)
              Fd_NL9(ITAU,IWAV,ISIZE) = EXP(U1)
              T_NL9(ITAU,IWAV,ISIZE) = EXP(T1)
              SBAR_NL9(ITAU,IWAV,ISIZE) = EXP(Z1)
              OPTH_NL9(ITAU,IWAV,ISIZE) = EXP(V1)
              IF (V1 .EQ. 0.) THEN
                OPTH_NL9(ITAU,IWAV,ISIZE) = V1
              ENDIF


25          CONTINUE
15        CONTINUE
 5      CONTINUE


        DO ISIZE= 1,NLSIZE
          DO IWAV=1,3
            DO ITAU  = 1,NLTAU
              INT_NL(ITAU,IWAV,ISIZE) = INT_NL9(ITAU,IWAV,ISIZE)
              FdT_NL(ITAU,IWAV,ISIZE) = FdT_NL9(ITAU,IWAV,ISIZE)
              Fd_NL(ITAU,IWAV,ISIZE) = Fd_NL9(ITAU,IWAV,ISIZE)
              T_NL(ITAU,IWAV,ISIZE) = T_NL9(ITAU,IWAV,ISIZE)
              SBAR_NL(ITAU,IWAV,ISIZE) = SBAR_NL9(ITAU,IWAV,ISIZE)
              OPTH_NL(ITAU,IWAV,ISIZE) = OPTH_NL9(ITAU,IWAV,ISIZE)
            ENDDO
            REF_RAY_NL(IWAV) = INT_NL(1,IWAV,ISIZE)
          ENDDO
        ENDDO


        RETURN
        END


c*************************************************
        SUBROUTINE AEROSOL_MAP(HANDLE_LUTMAP,IMONTH,AEROSOL)
C----------------------------------------------------
C!F77
C
C!DESCRIPTION:  Subroutine AEROSOL_MAP reads
C               the aerosol map (1 degree resolution)
C               to determine expected aerosol type at
C               a given location and season
C
C!INPUT PARAMETERS:
C         HANDLE_LUTMAP   Logical Unit # for file
C	  IMONTH          Integer month 1-12
C
C!OUTPUT PARAMETERS 
C	  AEROSOL	  360x180 map of aerosol type for appropriate season
C
C!REVISION HISTORY:
C 02/08/2006 Shana Mattoo/Rob Levy
C Initial revision
C
C!TEAM-UNIQUE HEADER:
C
C Developed by MODIS Aerosol Retrieval Team at NASA GSFC, Greenbelt, MD
C
C!REFERENCES AND CREDITS
C
C!DESIGN NOTES:
C
C!END
C----------------------------------------------------

        IMPLICIT none

        INTEGER IMONTH, IMONTH1,iseason
        INTEGER HANDLE_LUTMAP
	INTEGER nlon, nlat
	PARAMETER (nlon = 360, nlat = 180)
 	INTEGER AEROSOL(nlon,nlat)
 	INTEGER AEROSOL_all(nlon,nlat,4)
        INTEGER ilat, ilon
        CHARACTER*3 MMMs(4), MMM
        DATA MMMs/'DJF','MAM','JJA','SON'/

        DO iseason = 1, 4
         READ(HANDLE_LUTMAP,*) 
         DO ilat = 1, nlat
            READ(HANDLE_LUTMAP,99) 
     +      (AEROSOL_all(ilon,ilat,iseason), ilon=1, nlon)
         ENDDO
        ENDDO

        IMONTH1 = IMONTH
        IF (IMONTH .EQ. 12) THEN 
          IMONTH1 = 0
        ENDIF
        iseason = IMONTH1/3 + 1
        MMM = MMMs(iseason)
        PRINT*, 'MMM ', MMM
        DO ilon = 1, nlon
          DO ilat = 1, nlat
            AEROSOL(ilon,ilat) = AEROSOL_all(ilon,ilat,iseason)
          ENDDO
        ENDDO
        
 99     FORMAT (i1,359(1x,i1))
  
        RETURN 
        END

C*********************************************************************
      SUBROUTINE INTANGLE_RAY_NL(MTHET0,MTHET,MDPHI,
     *	     INT_NL0,REF_RAY_NL)		
C----------------------------------------------------------------------
C!F77
C
C!DESCRIPTION:  Subroutine INTANGLE_NL interpolates the lookup
C               Rayleigh reflectances to the measured geometry.
C
C!INPUT PARAMETERS:
C	  INT_NL0	radiance from LUT
C         MTHET0        Solar zenith angle.
C         MTHET        View angle.
C         MDPHI        Azimuth angle.
C
C         THET0      array of  solar Zenith angle in look_up table
C         THE        array of  satllite Zenith angle in look_up table
C         PHI        array of azimuth angle in look_up table
C
C!OUTPUT PARAMETERS 
C         REF_RAY_NL0      Reflectance for rayleigh only
C
C!REVISION HISTORY:
C
C!TEAM-UNIQUE HEADER:
C
C!REFERENCES AND CREDITS:
C
C!DESIGN NOTES:
C
C!END
C-----------------------------------------------------------------------

      IMPLICIT NONE
      SAVE

      INCLUDE 'mod04_land.inc'

c     Inputs
      REAL  MTHET0,MTHET,MDPHI
      INTEGER KSZAM1,KSZAP1,KTHEM1,KTHEP1,KPHIM1,KPHIP1
      
      REAL INT_NL0(NLPHI,NLTHE,NLTHET0,NLTAU,NLWAV,NLTABLE)

c     Outputs
      REAL REF_RAY_NL(NLWAV)

c     Dummy
      CHARACTER*132 LINE
      CHARACTER*45  LINE1
      INTEGER ICASE,IJ,IPHI,ITABLE,ITAU,ITH,ITH0,IWAV,IETA,LOPT
      REAL  X(2),Y(2),XX1(2),YY1(2),XX2(2),YY2(2),Y1
      INTEGER LL,MM,NN

C      Initialize
 
       KSZAM1=0
       KSZAP1=0
       KTHEM1=0
       KTHEP1=0
       KPHIM1=0
       KPHIP1=0

C      Find THET0 indices
       DO ITH0 = 1,NLTHET0-1
          IF ((MTHET0 .GE. THET0_NL(ITH0)) .AND. 
     *       (MTHET0 .LE. THET0_NL(ITH0+1))) THEN
	    KSZAM1=ITH0
	    KSZAP1=ITH0+1
          ENDIF
       ENDDO

C      Find THE indices
       DO ITH = 1,NLTHE-1
          IF ((MTHET .GE. THE_NL(ITH)) .AND. 
     *       (MTHET .LE. THE_NL(ITH+1))) THEN
	    KTHEM1=ITH
	    KTHEP1=ITH+1
          ENDIF
       ENDDO

C      Find PHI indices
       DO IPHI = 1,NLPHI-1
          IF ((MDPHI .GE. PHI_NL(IPHI)) .AND. 
     *       (MDPHI .LE. PHI_NL(IPHI+1))) THEN
	    KPHIM1=IPHI
	    KPHIP1=IPHI+1
          ENDIF
       ENDDO

C
C Initialize
C
        Y1=0.0
        DO IJ = 1,2
           X(IJ)=0.0
           Y(IJ)=0.0
           XX1(IJ)=0.0
           XX2(IJ)=0.0
           YY1(IJ)=0.0
           YY2(IJ)=0.0
        ENDDO

C
C Interpolate Rayleigh
C


          ITABLE = 1
          ITAU = 1
          DO 15 IWAV=1,NLWAV
              LL=0
              y1 = 0.0
              DO  40 ITH0  = KSZAM1,KSZAP1
                MM=0
                DO  50  ITH  = KTHEM1,KTHEP1
                  NN=0
                  DO 60  IPHI  = KPHIM1,KPHIP1
                    NN=NN+1
                    X(NN)=PHI_NL(IPHI)
                    Y(NN)=INT_NL0(IPHI,ITH,ITH0,ITAU,IWAV,ITABLE)
60                CONTINUE

                  CALL INTERP_EXTRAP(NN,MDPHI,X,Y,Y1,1)
                  MM=MM+1
                  XX1(MM)=THE_NL(ITH)
                  YY1(MM)=Y1

50              CONTINUE
                y1=0.0
                CALL INTERP_EXTRAP(MM,MTHET,XX1,YY1,Y1,1)
                LL=LL+1
                XX2(LL)=THET0_NL(ITH0)
                YY2(LL)=Y1

40            CONTINUE
              y1=0.0
              CALL INTERP_EXTRAP(LL,MTHET0,XX2,YY2,Y1,1)

            REF_RAY_NL(IWAV) = Y1
15        CONTINUE


        RETURN
        END


