atrem_app_refl_plus_gas_removal_for_l2gen3.f90

Go to the documentation of this file.

!------------------------------------------------------------------------------|
! The September 2012 version of atrem code, 'atrem_for_PRISM.f' (atrem for JPL |
!     PRISM Imaging Spectrometer, was obtained through modifying a 2011        |
!     version of 'atrem_f90_cubeio.f'.                                         |
! The March 2013 version of atrem code, 'atrem_for_PRISM_2013.f', was updated  |
!     from the September 2012 version of atrem code. The upgrades include:     |
!    a) replacing the MODTRAN 3.5 solar irradiance curve with a new solar curve|
!       constructed from the solar irradiance curve of Thuillier et al. 2004   |
!       ATLAS 3 (<=644.7 nm) and that of Kurucz 2005 (> 644.7 nm) built in     |
!       MODTRAN 5.2. This newly contructed solar curve is more suited for      |
!       modeling imaging imaging spectrometer data below 450 nm and with a     |
!       spectral resolution finer than 5 nm;                                   |
!    b) replacing high spectral resolution gas absorption tables with new      |
!       tables calculated with the HITRAN2008+ line database and a fast        |
!       line-by-line code originally developed by Bill Ridgway. After such     |
!       upgrades, the retrieved surface reflectances near 1.45 and 1.95 micron |
!       are improved significantly;                                            |
!    c) modifying volume mixing ratios of atmospheric carbon dioxide;          |
!    d) applying a scaling factor for improved modeling of atmospheric         |
!       oxygen band absorption centered near 1.265 micron.                     |
!------------------------------------------------------------------------------|
! The 2001 version of ATREM code uses the new 'cubeio.f90' for I/O operations. |
!     The output data cubes have no headers, unlike the previous generations   |
!     of ATREM codes containing 512 bytes of headers in output data cubes.     |
!------------------------------------------------------------------------------|
!********************************************************************************
!*                                         *
!* Name:           atrem_app_refl_plus_gas_removal_2013         April,  2013    *
!*                     (ATmosphere REMoval Program)                        *
!*                                         *
!* Author: Bo-Cai Gao, Remote Sensing Division, Code 7230,                      *
!*                     Naval Research Laboratory, Washington, DC 20375 USA      *
!*                                         *
!* Purpose: To derive surface reflectances from spectral imaging data collected *
!*             by the Airborne Visible Infrared Imaging Spectrometer (AVIRIS),  *
!*             HYDICE, HSI, PHILLS, HICO, NIS, PRISM, and possibly other        *
!*             airborne and spaceborne imaging spectrometers, and to derive a   *
!*             column water vapor image for each scene.                         *
!*                                         *
!* General Principles: In order to derive surface reflectances from image data  *
!*             cube, a thorough compensation for the atmospheric absorption and *
!*             scattering is required. The spatial and temporal variations of   *
!*             atmospheric water vapor amounts pose difficulties in removing    *
!*             water vapor absorption features in data cube using standard      *
!*             atmospheric models. In this algorithm, the amount of water vapor *
!*             on a pixel-by-pixel basis is derived from data cube themselves   *
!*             using the 0.94- and the 1.14-um water vapor bands and a          *
!*             three-channel ratioing technique. The derived water vapor        *
!*             values are then used for modeling water vapor absorption effects *
!*             in the entire 0.4-2.5 um region. The absorption effects of       *
!*             well mixed atmospheric gases, such as CO2, N2O, CH4, and O2,     *
!*             are modeled using standard atmospheric models. A line-by-line    *
!*             code developed by W. Ridgway at NASA Goddard Space Flight Center *
!*             is used to generate a database of gaseous absorption             *
!*             coefficients at high spectral resolution. This database and      *
!*             ozone cross sections and NO2 cross sections are used in          *
!*             calculating transmittances of eight atmospheric gases (H2O, CO2  *
!*             O3, N2O, CO, CH4, O2, and NO2). The scattering effect is modeled *
!*             using a modified version of 6S code (Vermote et al. 1996).       *
!*                                         *
!* Algorithm:  1. An input parameter file is read in and global data are        *
!*       initialized.                              *
!*        2. The solar zenith angle is derived based on the flight time    *
!*       and latitude and longitude of the scene.              *
!*        3. A table consisting of transmittance spectra at the solar      *
!*       and observational geometry but with varying amounts of water  *
!*           vapor is created.                             *
!*           The transmittance spectra for water vapor (H2O), carbon       *
!*       dioxide (CO2), ozone (O3), nitrous oxide (N2O), carbon        *
!*           monoxide (CO), methane (CH4), oxygen (O2), and nitrogen       *
!*                dioxide (NO2) in the 0.4-2.5 micron region are calculated.    *
!*             4. Aerosol and molecular scattering terms are calculated using   *
!*       a modified version of the 6S code.                *
!*             5. The image data cube (2D spatial and 1D spectral) are accessed *
!*           one spectral slice at a time from an image file in any storage*
!*           order, Band SeQuential (BSQ), Band Interleaved by Pixel (BIP),*
!*                and Band Interleaved by Line (BIL). Each measured radiance    *
!*           spectrum is divided by the solar irradiance curve above the   *
!*       atmosphere to get apparent reflectances.              *
!*             6. The integrated water vapor amount for each apparent           *
!*       reflectance spectrum is derived from the 0.94- and the 1.14-um*
!*           water vapor absorption bands using a 3-channel ratioing       *
!*       technique and a look-up table procedure.              *
!*        7. The surface reflectances are derived using another look-up    *
!*       table procedure. Routines written in C language are used to   *
!*                perform all the input/output operations.             *
!*                                         *
!*  Limitations:                                           *
!*    1.  This algorithm works best for areas where the surface elevation  *
!*             does not vary by more than 1 km.                    *
!*    2.  The algorithm works well only for measurements taken on clear    *
!*             days. For hazy days with large aerosol optical depths, the       *
!*             algorithm does not work nearly as well.                          *
!*    3.  The algorithm does not work well over dark areas such as rivers  *
!*        and lakes. It does not perform corrections for "atmospheric      *
!*             adjacency" effects.                         *
!*    4.  The algorithm assumes horizontal surfaces having Lambertian      *
!*             reflectances.                                               *
!*         5.  Additional work is needed to port this algorithm to different    *
!*             computing systems. This algorithm was developed on an SGI        *
!*             workstation and used FORTRAN statements to read several binary   *
!*             files. The binary file format on different computers can be      *
!*             different. The record lengths defined on an SGI machine are      *
!*             different from most other computers. Some of the Dec Alpha       *
!*             workstations use 8 bytes (64 bits) to store a C-pointer, instead *
!*             of 4 bytes (32 bits). All these factors need to be taken into    *
!*             considerations when porting this algorithm to other computers.   *
!*                                         *
!*  User Input:                                        *
!*    1.  Name - Imaging spectrometer name, e.g., AVIRIS, HYDICE, HSI, PHYLLS.  *
!*    2.  Plane altitude in km, above the sea level.                            *
!*    3.  Date - The month, day, and year that the data were collected.         *
!*    4.  Time - The hour, minute, second (GMT) that the data were         *
!*                 collected.                                                   *
!*    5.  Latitude - The mean latitude of the measurement area in degrees,      *
!*           minutes, seconds, and hemisphere.                                     *
!*    6.  Longitude - The mean longitude of the measurement area in             *
!*           degrees, minutes, seconds, and hemisphere.                            *
!*    7.  View zenith angle - the angle  in degrees, minutes, and seconds       *
!*                            between the viewing direction and nadir           *
!*    8.  View azimuth angle in degree, minutes and seconds. The view azimuth   *
!*                            angle is defined as the angle between local North *
!*                            the view direction projected on the ground        *
!*                            (clockwise count). This angle can vary from 0     *
!*                            to 360 deg. Ths definition here is consistent with*
!*                            navigator's convention, but not consistent with   *
!*                            astronomer's convention. Here is an example: for  *
!*                            a downward looking instrument pointed off-nadir   *
!*                            to the West, the view azimuthal angle is 90 deg.  *
!*                            Conversely, if the instrument is pointed toward   *
!*                            East, the view azimuthal angle is 270 deg.        *
!*                                         *
!*    9.  Resolution of input spectra in nm (~ 10 nm for AVIRIS) or "0"         *
!*        if the FWHM values (in micron) are present in the wavelength file     *
!*   10.  Filename of spectrometer wavelength table - 2 or 3 column ASCII data. *
!*           The first column is the channel number, the second column is       *
!*           the wavelength and the optional third column is the FWHM value for *
!*           each channel.                                     *
!*   11.  Channel ratio parameters - Center positions and widths of window and  *
!*           water vapor channels used in channel ratios for deriving              *
!*      the amount of water vapor from imaging data. Each of the window    *
!*           and water vapor channels typically consists of several narrow      *
!*           imaging spectrometer channels.                                *
!*   12.  Atmospheric model - Temperature, pressure, water vapor volume mixing  *
!*           ratio profiles. The model can be either a predefined model or a    *
!*           user defined model.                               *
!*   13.  Gases - An indication of which of the 8 gases (H2O, CO2, O3,          *
!*      N2O, CO, CH4, O2, and NO2) should be included in the spectral      *
!*      calculations.                                              *
!*   14. Ozone - Column amount of O3 - changes with season and                  *
!*      latitude.  Typical value is .34 atm-cm.                            *
!*   15. Aerosol model number and visibility when measurements were taken.      *
!*   16. Aerosol optical depth at 550 nm (Optional, only if visibility V        *
!*           is set to 0)                              *
!*   17. Surface elevation - The average surface elevation of an imaging scene  *
!*           in km.                                                             *
!*   18. Filename of input image file                          *
!*   19. Dimensions of input image.                        *
!*   20. Filename of output image file in same storage order as input file      *
!*   21. Resolution of output surface reflectance spectra in micron             *
!*   22. Scale factor for output reflectance values                             *
!*   23. Filename of output water vapor image.                     *
!*   24. Filename of output atmospheric transmittance look-up table             *
!*                                         *
!*  Output:                                        *
!*     a. Output to user specified files:                      *
!*        1. Surface reflectance cube data - retrieved from the image data cube *
!*        2. Water vapor image - an image of the derived column water vapor     *
!*           amounts at each pixel.                        *
!*        3. Transmittance Look-up table - consisting of channel ratio values   *
!*           for the .94- and 1.14-um channels corresponding to 60 water vapor  *
!*           amounts, and the associated atmospheric transmittance spectra.     *
!*     b. Output to standard output:                                       *
!*        1. Debugging information - if the source is compiled with the         *
!*           debug flag set (-d_lines for the ULTRIX compiler),                    *
!*           then debug information is written out.                        *
!*        2. Error messages - if any of the user input is invalid, or there is  *
!*           an I/O error, a message is written out, and the program halts.     *
!*        3. Progress indicator - a message, which shows the number of the      *
!*           spectral slices that have been processed, to give the user an      *
!*           indication of the progress of the program execution.               *
!*                                         *
!*  Special Considerations:                                *
!*     a. Make sure that the imaging spectrometer's channel positions specified *
!*           in the input wavelength table are correct. The incorrect channel   *
!*           positions will introduce sharp features in derived surface         *
!*           reflectance spectra.                                               *
!*     b. The output reflectance cube file has the same size as the input       *
!*           data cube file. Make sure there is enough space in the file        *
!*           system for the output cube before running this program.           *
!*                                         *
!*  Change History:                                *
!*     ATREM 1.2 - uses full-width half-max values to smooth solar irradiance   *
!*                 curve and a new scaling factor for methane              *
!*     ATREM 1.3 - used new solar irradiance curve                 *
!*      - supports 1992 AVIRIS data                    *
!*     ATREM 1.3.1 - allows 0 and above for output resolution                   *
!*     ATREM 2.0 - allows variable scale factors for each spectrometer          *
!*                 (code submitted by Roger Clark, USGS)                        *
!*               - user can input output scale factor for reflectance values    *
!*               - user can input radiance file header size in bytes            *
!*               - new solar irradiance curve (Neckel and Labs plus ATMOS)      *
!*                 used (Green and Gao, 1993).                                  *
!*     ATREM 3.0 - replaced the 5S code by the 6S code for modeling atmospheric *
!*                 scattering effects and for modeling measurements from        *
!*                 low-altitude aircrafts. ATREM 3.0, like previous versions    *
!*                 of ATREM, only models nadir viewing geometry.                *
!*               - changed algorithms for calculating atmospheric gaseous       *
!*                 transmittances for the upward surface-sensor path            *
!*      - users need to specify aircraft altitude (km) in input file   *
!*     ATREM 4.0 - replaced the band model by a line-by-line-based algorithm    *
!*                 for calculating atmospheric gaseous transmittances. This     *
!*                 allows the modeling of imaging spectrometer data at          *
!*                 spectral resolutions better than 10 nm.                      *
!*               - increased the buffer size to 1024x1024 in order to handle    *
!*                 images larger than the AVIRIS' size (614x512).               *
!*               - modified the algorithm to allow off-nadir pointing geometry. *
!*               - users need to specify view zenith angle and view azimuth     *
!*                 angle. The view azimuth angle is defined as the angle        *
!*                 between local North and the view direction projected on the  *
!*                 ground (clockwise count). This angle can vary from 0 to      *
!*                 360 deg.                                                     *
!*               - replaced the solar irradiance curve (Neckel and Labs plus    *
!*                 ATMOS) by the high spectral resolution solar irradiance      *
!*                 curve from MODTRAN 3.5 released in December of 1996.         *
!*     ATREM 4.1 - included atmospheric NO2 in transmittance calculations.      *
!*                                         *
!*  Acknowledgments:                                   *
!*           This work was partially supported over several years by grants     *
!*           from NASA Jet Propulsion Laboratory, California Institute of       *
!*           Technology, from NASA Headquarters, and from the Office of Naval   *
!*           Research. Special thanks goes to Kathy Heidebrecht at the Center   *
!*           for the Study of Earth from Space, University of Colorado at       *
!*           Boulder, Colorado, for supporting the development of the earlier   *
!*           versions of ATREM codes. Special thanks also goes to William L.    *
!*           Ridgway for providing the line-by-line gaseous absorption database *
!*           used in the present version of ATREM.                              *
!*                                                     *
!*  References:                                        *
!*        Gao, B.-C., K. Heidebrecht, and A. F. H. Goetz, Derivation of scaled  *
!*           surface reflectances from AVIRIS data, Remote Sens. Env., 44,      *
!*           165-178, 1993.                            *
!*        Gao, B.-C., and C. O. Davis, Development of an operational algorithm  *
!*           for removing atmospheric effects from HYDICE and HSI data,         *
!*           in SPIE'96 Conference Proceedings, Vol. 2819, 45-55, 1996.         *
!*        Gao, B.-C., and A. F. H. Goetz, Column atmospheric water vapor and    *
!*           vegetation liquid water retrievals from airborne imaging          *
!*           spectrometer data, J. Geophys. Res., 95, 3549-3564, 1990.         *
!*        Goetz, A. F. H., and M. Herring, The high resolution imaging         *
!*           spectrometer (HIRIS) for Eos, IEEE Trans. Geosci. Remote Sens., 27,*
!*           136-144, 1989.                            *
!*        Goetz, A. F. H., G. Vane, J. Solomon, and B. N. Rock, Imaging        *
!*           spectrometry for Earth remote sensing, Science, 228, 1147-1153,1985*
!*        Green, R. O., and B.-C. Gao, A Proposed Update to the Solar Irradiance*
!*           Spectrum used in LOWTRAN and MODTRAN, in Summaries of the Fourth   *
!*           Annual JPL Airborne Geoscience Workshop, October 25-29, (Editor,   *
!*           R. O. Green), JPL Publ. 93-26, Vol. 1, pp. 81-84, Jet Propul. Lab, *
!*           Pasadena, Calif., 1993.                           *
!*        Kneizys, F. X., E. P. Shettle, L. W. Abreu, J. H. Chetwynd, G. P.     *
!*           Anderson, W. O. Gallery, J. E. A. Selby, and S. A. Clough, Users   *
!*           guide to LOWTRAN7, AFGL-TR-8-0177, Air Force Geophys. Lab.,        *
!*           Bedford, Mass., 1988.                         *
!*        Iqbal, M., An Introduction To Solar Radiation, Academic, San Diego,   *
!*           Calif., 1983.                             *
!*        Malkmus, W., Random Lorentz band model with exponential-tailed S line *
!*           intensity distribution function, J. Opt. Soc. Am., 57, 323-329,1967*
!*        Press, W. H., B. P. Flannery, S. A. Teukolsky, and W. T.  Vetterling, *
!*           Numerical Recipes-The ART of Scientific Computing, Cambridge       *
!*           University Press, 1986.                           *
!*        Rothman, L. S., et al., The HITRAN 2008 molecular spectroscopic       *
!*           database, JQSRT, 110, 533-572, 2009.                               *
!*        Solomon, S., R. W. Portmann, R. W. Sanders, J. S. Daniel, W. Madsen,  *
!*           B. Bartram, and E. G. Dutton, On the role of nitrogen dioxide in   *
!*           the absorption of solar radiation, J. Geophys. Res., 104,          *
!*           12047-12058, 1999.                                                 *
!*        Tanre, D., C. Deroo, P. Duhaut, M. Herman, J. J. Morcrette, J. Perbos,*
!*           and P. Y. Deschamps, Description of a computer code to simulate    *
!*           the satellite signal in the solar spectrum: the 5S code, Int.      *
!*           J. Remote Sens., 11, 659-668, 1990.                       *
!*        Tanre, D., C. Deroo, P. Duhaut, M. Herman, J. J. Morcrette, J. Perbos,*
!*           and P. Y. Deschamps, Simulation of the satellite signal in the     *
!*           solar spectrum (5S), Users' Guide (U. S. T. De Lille, 59655        *
!*           Villeneu d'Ascq, France: Laboratoire d'Optique Atmospherique),     *
!*      1986.                                  *
!*        Thuillier, G., et al., Solar irradiance reference spectra for two     *
!*           solar active levels, Adv. Space Res., 34, 256-261, 2004.           *
!*        Vane, G., R. O. Green, T. G. Chrien, H. T. Enmark, E. G. Hansen, and  *
!*           W. M. Porter, The Airborne Visible/Infrared Imaging Spectrometer,  *
!*           Remote Sens. Env., 44, 127-143, 1993.                              *
!*        Vane, G. (Ed), Airborne visible/infrared imaging spectrometer        *
!*      (AVIRIS), JPL Publ. 87-38, Jet Propul. Lab, Pasadena, Calif., 1987.*
!*        Vermote, E., D. Tanre, J. L. Deuze, M. Herman, and J. J. Morcrette,   *
!*           Second simulation of the satellite signal in the solar spectrum    *
!*           (6S), 6S User's Guide Version 1, NASA-GSFC, Greenbelt, Maryland,   *
!*           134 pages, 1994.                                                   *
!*                                         *
!********************************************************************************
!
!
!********************************************************************************
!*                                         *
!*  Name: GET_INPUT                                *
!*  Purpose: Reads and verifies all user provided input.                   *
!*  Parameters: none                                   *
!*  Algorithm: Data is read in from standard input. The user is not prompted    *
!*             interactively, so data should be redirected from an input file.  *
!*             The data is validated after it is read in.  If the data is       *
!*             invalid, a message is printed and the program stops.            *
!*  Globals used:   TPVMR - two dimensional array containing 6 predefined       *
!*                atmospheric models.  The models contain values for   *
!*            the number of atmospheric layer boundaries, altitude,*
!*            temperature, pressure, and water vapor volume mixing *
!*            ratio.                           *
!*  Global output:  IH2OVP, ICO2, IO3, IN2O, ICO, ICH4, IO2 - set to 0 to       *
!*                         indicate that the gas should NOT be used in the      *
!*                         calculations, and set to 1 if the gas should be used *
!*         H(), T(), P(), VMR(), NB, NL, MODEL - altitude (km),        *
!*            temperature (K), pressure (atm), water vapor volume  *
!*            mixing ratio (ppm), number of atmospheric layer      *
!*            boundaries, number of atmospheric layers (NB-1),     *
!*                         and model number for the selected atmospheric model. *
!*         VRTO3 - column amount of O3.                    *
!*         WAVOBS() - wavelengths for all channels.                    *
!*         FWHM() -   resolutions for each channel.                    *
!*                  NOBS - number of AVIRIS wavelengths.                        *
!*                  HSURF - the mean surface elevation of the imaging scene     *
!*                  DLT, DLT2 - resolution, in units of nm, of input            *
!*                         spectra and resolution of output surface reflectance *
!*                         spectra. If DLT2>DLT, output spectra are smoothed    *
!*                         using a gaussian function.                      *
!*             WNDOW1, WNDOW2, WP94C - center wavelength positions of two  *
!*            broad window channels and one broad .94-um water     *
!*            vapor channel.                       *
!*         WNDOW3, WNDOW4, W1P14C - center wavelength positions of two *
!*            broad window channels and one broad 1.14-um water    *
!*            vapor channel.                       *
!*         NB1, NB2, NBP94, NB3, NB4, NB1P14 - number of individual    *
!*            narrow channels that form the corresponding broad    *
!*            window and water vapor absorption channels.          *
!*         IMN, IDY, IYR, IH, IM, IS - month, day, year, hour, minute, *
!*                         and second of measurement.                  *
!*         XLATD, XLATM, XLATS, LATHEM - degrees, minutes, seconds, and*
!*                         hemisphere of latitude of measured area.             *
!*         XLONGD, XLONGM, XLONGS, LNGHEM - degrees, minutes, seconds, *
!*                         and hemisphere of latitude of measured area.         *
!*                  NAME_INSTRU: Imaging spectrometer name, e.g., AVIRIS, HYDICE*
!*                  XPSS:  = HSURF, an interface for using 6S.                  *
!*                  XPPP:  plane height (km, above the sea level).              *
!*                                                     *
!*  Return Codes: none                                 *
!*  Special Considerations: None                               *
!*                                         *
!********************************************************************************
      
      SUBROUTINE get_input
 
      use cubeio
 
      include 'COMMONS_INC.f'
 
      INTEGER              LUN_IN, LUN_OUT, LUN_VAP, I_RET
      COMMON /inout_units/ lun_in, lun_out, lun_vap
 
 
!C  Common variables
      dimension h(25), t(25), p(25), vmr(25)
      dimension wavobs(1024),fwhm(1024)
      dimension tpvmr(7,81)
      CHARACTER*1 LATHEM, LNGHEM
 
      CHARACTER (LEN = 1000) :: FINAV,FOCUB,FOH2O
      INTEGER SORDER,HDREC
 
!C  Local variables
      INTEGER ANS
      CHARACTER (LEN = 1000) :: FINPWV,FTPVMR      
      LOGICAL GOOD_DATA
      CHARACTER (LEN = 1000) :: FOUT1
      INTEGER DIMS(4)
      INTEGER ST_ORDER
      
      COMMON /getinput1/ ih2ovp,ico2,io3,in2o,ico,ich4,io2,ino2
      COMMON /getinput3/ h,t,p,vmr,nb,nl,model,iaer,v,taer55,vrto3,sno2
      COMMON /getinput4/ wavobs,fwhm
      COMMON /getinput5/ nobs,hsurf,dlt,dlt2
      COMMON /getinput6/ wndow1,wndow2,wp94c,wndow3,wndow4,w1p14c
      COMMON /getinput7/ nb1,nb2,nbp94,nb3,nb4,nb1p14
      COMMON /getinput8/ imn,idy,iyr,ih,im,is
      COMMON /getinput9/ xlatd,xlatm,xlats,lathem
      COMMON /getinput10/xlongd,xlongm,xlongs,lnghem
      COMMON /getinput11/hdrec,nsamps,nlines,nbands,sorder
      COMMON /getinput12/scalef
      COMMON /tpvmr_init1/ tpvmr
 
!C Commons for use with the C programs for cube I/O
      COMMON /outcube/ focub
      COMMON /incube/ finav
      COMMON /outh2ovap/ foh2o
 
!C Parameters for names of imaging spectrometers
      CHARACTER (LEN = 80) :: NAME_INSTRU, NAMES(10)
      COMMON /getinput13/ name_instru, names
 
      REAL XPSS, XPPP
      COMMON /getinput14/ xpss, xppp
 
      REAL XVIEWD, XVIEWM, XVIEWS
      REAL XAZMUD, XAZMUM, XAZMUS
      COMMON /getinput15/ xviewd,xviewm,xviews, xazmud,xazmum,xazmus
 
      good_data = .true.    ! initialize flag for good data
 
!C  Get the name of imaging spectrometer, e.g., AVIRIS, HYDICE
      READ(5,627) name_instru
 627  FORMAT(a10)
!C
!C
!C***Temp code for assigning names of imaging spectrometers. Based on these
!C    names, different scale factors should be used for different instrument
!C    when converting measured radiances to standard radiance units.
!C*** The coding here may need to be moved to the file "COMMONS_INC"
      names(1) = 'AVIRIS'
      names(2) = 'HYDICE'
      names(3) = 'HSI'
      names(4) = 'TRWIS-III'
      names(5) = 'PHYLLS'
      names(6) = 'Hyperion'
      names(7) = 'HICO'
      names(8) = 'NIS'
      names(9) = 'PRISM'
      names(10)= 'OTHERS?'
!C***End of temp coding ------------
 
      print*, 'Instrument Name: ', name_instru
 
!C Get the plane altitude (in km, above sea level). Must > HSURF.
      READ(*,*) xppp
 
!C--      print*, 'XPPP = ', XPPP
 
!C  Get date and time of measurements.  Format: MM DD YYYY HH MM SS
      READ (5,*) imn,idy,iyr,ih,im,is
      IF((imn.LT.1).OR.(imn.GT.12).OR.(idy.LT.1).OR.(idy.GT.31).OR. &
        (iyr.LT.1987).OR.(iyr.GT.2020)) THEN
        WRITE(*,*)'Invalid date:',imn,idy,iyr
        WRITE(*,*)'Format is MM DD YYYY where 0<MM<12, 0<DD<32,  &
      YYYY>1986.'
        good_data = .false. 
      ENDIF
      IF((ih.LT.0).OR.(ih.GT.24).OR.(im.LT.0).OR.(im.GT.60).OR. &
        (is.LT.0).OR.(is.GT.60)) THEN
        WRITE(*,*)'Invalid time:',ih,im,is
        WRITE(*,*)'Format is HH MM SS where 0<HH<24, 0<MM<60, SS<60.'
        good_data = .false. 
      ENDIF
 
!C  Get latitude of measurements.  Format: degrees minutes seconds
      READ (5,*) xlatd,xlatm,xlats
      IF((xlatd.GT.90).OR.(xlatm.GT.60).OR.(xlats.GT.60)) THEN
        WRITE(*,*)'Invalid latitude:',xlatd,xlatm,xlats
        WRITE(*,*)'Format: degrees minutes seconds'
        WRITE(*,*)'Valid values are degrees < 90, minutes < 60, &
      seconds<60.'
        good_data = .false. 
      ENDIF
 
!C Get hemisphere of latitude.  Format: "N" or "S"
      READ (5,127) lathem
  127 FORMAT(a1)
      IF((lathem.NE.'N').AND.(lathem.NE.'S')) THEN
        WRITE(*,*)'Invalid hemisphere value:',lathem
        WRITE(*,*)'Valid values are "N" or "S".'
        good_data = .false. 
      ENDIF
 
!C  Get longitude of measurements.  Format: degrees minutes seconds
      READ (5,*) xlongd,xlongm,xlongs
      IF((xlongd.GT.180).OR.(xlongm.GT.60).OR.(xlongs.GT.60)) THEN
        WRITE(*,*)'Invalid longitude:',xlongd,xlongm,xlongs
        WRITE(*,*)'Format: degrees minutes seconds'
        WRITE(*,*)'Valid values are degrees < 180 minutes < 60, &
      seconds<60.'
        good_data = .false. 
      ENDIF
 
!C Get hemisphere of longitude.  Format: "E" or "W"
      READ (5,127) lnghem
      IF((lnghem.NE.'E').AND.(lnghem.NE.'W')) THEN
        WRITE(*,*)'Invalid hemisphere:',lnghem
        WRITE(*,*)'Valid values are "E" or "W".'
        good_data = .false. 
      ENDIF
 
!C  Get view zenith angle.  Format: degrees minutes seconds
      READ (5,*) xviewd,xviewm,xviews
      IF((xviewd.GT.90).OR.(xviewm.GT.60).OR.(xviews.GT.60)) THEN
        WRITE(*,*)'Invalid latitude:',xviewd,xviewm,xviews
        WRITE(*,*)'Format: degrees minutes seconds'
        WRITE(*,*)'Valid values are degrees < 90, minutes < 60 & seconds<60.'
        good_data = .false. 
      ENDIF
 
!C  Get relative azimuth angle between the solar ray and view direction with
!C      both projected on the ground.  Format: degrees minutes seconds
      READ (5,*) xazmud,xazmum,xazmus
      IF((xazmud.GT.360.).OR.(xazmum.GT.60).OR.(xazmus.GT.60)) THEN
        WRITE(*,*)'Invalid view azimuth ang.:',xazmud,xazmum,xazmus
        WRITE(*,*)'Format: degrees minutes seconds'
        WRITE(*,*)'Valid values are degrees < 360, minutes < 60, & seconds<60.'
        good_data = .false. 
      ENDIF
 
!C Get spectral resolutions (FWHM) in micron. If the value is 0, that means
!C read the FWHM data from the wavelength file.  Otherwise, a non-zero value
!C will be used as a the constant FWHM value.
      READ(*,*)dlt
      IF (dlt.GT.0.0) THEN
        IF((dlt.LT.0.0005).OR.(dlt.GT.0.1)) THEN
          WRITE(*,*)'Invalid spectral resolution value:',dlt
          WRITE(*,*)'Valid values are 0.0005-0.1 micron.'
          good_data = .false. 
        ELSE
!C Initialize the FWHM array
!C***          DLTNEW=DLT/1000.
          DO 521 i=1,1024
            fwhm(i)=dlt
  521     CONTINUE
        ENDIF
      ENDIF
 
!C  Get imaging spectrometer wavelength file name
      READ(*,85) finpwv
  85  FORMAT(a1000)
      OPEN(11,file=finpwv,status='OLD',iostat=istat)
      IF(istat.NE.0) THEN
        WRITE(*,*)'wavelength file did not open successfully:',finpwv
        good_data = .false. 
      ELSE
 
!C  Read wavelength data.  For 1991 and before, the file contains the channel
!C  number and the wavelength for each channel.  For 1992, the file also contains
!C  a third column with the resolution (FWHM) in micron for each channel.
 
        nobs = 1024           !initially set to max allowable # of observations
        DO 510 i=1,nobs
          IF (dlt.GT.0.) THEN
            READ(11,*,END=520) X,WAVOBS(I)
          ELSE
            READ(11,*,END=520) X,WAVOBS(I), FWHM(I)
          ENDIF
 510    CONTINUE
 520    nobs=i-1              !set to actual # of observations in input file
        CLOSE(11)
      ENDIF
 
!C Determine if channel ratio parameters should be read.  Format: 0=no 1=yes
      READ(*,*) ans
      IF((ans.NE.0).AND.(ans.NE.1)) THEN
    WRITE(*,*)'Invalid value to indicate whether the channel ratio parameters should be read:',ans
        WRITE(*,*)'Valid values: 0=no 1=yes.'
        good_data = .false. 
      ELSE
        IF (ans .EQ. 1) THEN
!C  Get ranges on curve for the two atmospheric windows surrounding the .94-um
!C    water vapor absorption feature and the center point of the .94-um water
!C    vapor absorption feature.  Enter:
!C         1. the midpoint of first window (0.6-2.5)
!C         2. number of points to average for first window (1-10)
!C         3. the midpoint of second window (0.6-2.5)
!C         4. number of points to average for second window (1-10)
!C         5. the midpoint of water vapor absorption feature (0.6-2.5)
!C         6. the number of points to average absorption feature (1-30)
 
  114     READ(*,*)wndow1, nb1, wndow2, nb2, wp94c, nbp94
 
          IF((wndow1.LT.0.6).OR.(wndow1.GT.2.5)) THEN
        WRITE(*,*)'Invalid wavelength position for first atmospheric window in the .94-um region:',wndow1
            WRITE(*,*)'Valid values are 0.6-2.5.'
            good_data = .false. 
          ENDIF
 
          IF((nb1.LT.1).OR.(nb1.GT.50)) THEN
        WRITE(*,*)'Invalid number of channels for first wavelength position in the .94-um region:',nb1
            WRITE(*,*)'Valid values are 1-50.'
            good_data = .false. 
          ENDIF
 
          IF((wndow2.LT.0.6).OR.(wndow2.GT.2.5)) THEN
        WRITE(*,*)'Invalid wavelength position for second atmospheric window in the .94-um region:',wndow2
            WRITE(*,*)'Valid values are 0.6-2.5.'
            good_data = .false. 
          ENDIF
 
          IF((nb2.LT.1).OR.(nb2.GT.50)) THEN
        WRITE(*,*)'Invalid number of channels for second wavelength position in the .94-um region:',nb2
            WRITE(*,*)'Valid values are 1-50.'
            good_data = .false. 
          ENDIF
 
          IF((wp94c.LE.wndow1).OR.(wp94c.GE.wndow2)) THEN
        WRITE(*,*)'Invalid water vapor wavelength position for the .94-um region:',wp94c
            WRITE(*,*)'Valid range is:',wndow1,' < value < ',wndow2,'.'
            good_data = .false. 
          ENDIF
 
          IF((nbp94.LT.1).OR.(nbp94.GT.90)) THEN
        WRITE(*,*)'Invalid number of channels for water vapor wavelength position in the .94-um region:',nbp94
            WRITE(*,*)'Valid values are 1-90.'
            good_data = .false. 
          ENDIF
 
!C  Get ranges on curve for the two atmospheric windows surrounding the 1.14-um
!C    water vapor absorption feature and the center point of the 1.14-um water
!C    vapor absorption feature.  Enter:
!C         1. the midpoint of third window (0.6-2.5)
!C         2. number of points to average for third window (1-10)
!C         3. the midpoint of fourth window (0.6-2.5)
!C         4. number of points to average for fourth window (1-10)
!C         5. the midpoint of 1.14-um absorption feature (0.6-2.5)
!C         6. the number of points to average for the absorption feature (1-30)
 
 
          READ(*,*)wndow3, nb3, wndow4, nb4, w1p14c, nb1p14
 
          IF((wndow3.LT.0.6).OR.(wndow3.GT.2.5)) THEN
        WRITE(*,*)'Invalid wavelength position for first atmospheric window in the 1.14-um region:',wndow3
            WRITE(*,*)'Valid values are 0.6-2.5'
            good_data = .false. 
          ENDIF
 
          IF((nb3.LT.1).OR.(nb3.GT.50)) THEN
        WRITE(*,*)'Invalid number of channels for first wavelength position in the 1.14-um region:',nb3
            WRITE(*,*)'Valid values are 1-50.'
            good_data = .false. 
          ENDIF
 
          IF((wndow4.LT.0.6).OR.(wndow4.GT.2.5)) THEN
        WRITE(*,*)'Invalid wavelength position for second atmospheric window in the 1.14-um region:',wndow4
            WRITE(*,*)'Valid values are 0.6-2.5'
            good_data = .false. 
          ENDIF
 
          IF((nb4.LT.1).OR.(nb4.GT.50)) THEN
        WRITE(*,*)'Invalid number of channels for second wavelength position in the 1.14-um region:',nb4
            WRITE(*,*)'Valid values are 1-50.'
            good_data = .false. 
          ENDIF
 
          IF((w1p14c.LE.wndow3).OR.(w1p14c.GE.wndow4)) THEN
        WRITE(*,*)'Invalid water vapor wavelength position for 1.14-um:',w1p14c
            WRITE(*,*)'Valid range is:',wndow3,' < value <', wndow4,'.'
            good_data = .false. 
          ENDIF
 
          IF((nb1p14.LT.1).OR.(nb1p14.GT.110)) THEN
        WRITE(*,*)'Invalid number of channels for water vapor wavelength position in the 1.14-um region:',nb1p14
            WRITE(*,*)'Valid values are 1-110.'
            good_data = .false. 
          ENDIF
        ELSE
 
!C Initialize default atmospheric window and water vapor absorption regions
!C     for 3-channel ratioing calculations.
          wndow1 = 0.865
          nb1    = 3
          wndow2 = 1.030
          nb2    = 3
          wp94c  = 0.940
          nbp94  = 5
          wndow3 = 1.050
          nb3    = 3
          wndow4 = 1.235
          nb4    = 3
          w1p14c = 1.1375
          nb1p14 = 7
        ENDIF
      ENDIF  
 
!C  Determine which atmospheric model to use
!C     1 = tropical
!C     2 = mid latitude summer
!C     3 = mid latitude winter
!C     4 = subarctic summer
!C     5 = subarctic winter
!C     6 = US standard 1962
!C     7 = user defined model
      READ(*,*)model 
      IF((model.LT.1).OR.(model.GT.7)) THEN
    WRITE(*,*)'Invalid atmospheric model value:',model
        WRITE(*,*)'Valid values are 1-7.'
        good_data = .false. 
        model = 1       ! init model to a valid value for use below
      ENDIF
      
      IF (model.EQ.7) THEN
!C Get file name of atmospheric model information
        READ(*,101) ftpvmr
  101   FORMAT(a1000)
        OPEN(20,file=ftpvmr,status='OLD',iostat=istat)
        IF(istat.NE.0) THEN
          WRITE(*,*)'Atmospheric model file did not open  successfully:',ftpvmr
          good_data = .false. 
        ELSE
!C Read in number of boundaries, temperature (K), pressure (mb), and water
!C vapor volume mixing ratio (VMR, in units of parts per million (ppm)) profiles
          READ(20,*) nb
          IF((nb.LE.0).OR.(nb.GT.25)) THEN
        WRITE(*,*)'Invalid number of boundaries:',nb
            WRITE(*,*)'Valid values are 1-25.'
            good_data = .false. 
          ELSE
            DO 300 i=1,nb
              READ(20,*) h(i),p(i),t(i),vmr(i)
  300       CONTINUE
            CLOSE(20,iostat=istat)
          ENDIF
        ENDIF
      ELSE
!C       Initialize NB, H, P, T, VMR from predefined atmospheric model array
        nb = tpvmr(model,1)
        DO 120 i=1,nb
          h(i) = tpvmr(model,2+(4*(i-1)))
          p(i) = tpvmr(model,3+(4*(i-1)))
          t(i) = tpvmr(model,4+(4*(i-1)))
          vmr(i) = tpvmr(model,5+(4*(i-1)))
  120   CONTINUE
      ENDIF
      nl=nb-1
 
        DO i     = nb+1, 25
          h(i)   = 1000.
          p(i)   = 0.0
          t(i)   = 300.
          vmr(i) = 0.0
        END DO
!
!C Determine if various gases should be included in the calculations.  Format:
!C 1=yes or 0=no.  The order should be:
!C                     1. water vapor
!C                     2. carbon dioxide
!C                     3. ozone
!C                     4. nitrous oxide
!C                     5. carbon monoxide
!C                     6. methane
!C                     7. oxygen
!C                     8. nitrogen dioxide
 
      READ(*,*)ih2ovp, ico2, io3, in2o, ico, ich4, io2, ino2
      IF((ih2ovp.NE.0).AND.(ih2ovp.NE.1)) THEN
    WRITE(*,*)'Invalid selection for H2O vapor:',ih2ovp
        WRITE(*,*)'Valid values: 0=no, 1=yes.'
        good_data = .false. 
      ENDIF
 
      IF((ico2.NE.0).AND.(ico2.NE.1)) THEN
    WRITE(*,*)'Invalid selection for CO2:',ico2
        WRITE(*,*)'Valid values: 0=no, 1=yes.'
        good_data = .false. 
      ENDIF
 
      IF((io3.NE.0).AND.(io3.NE.1)) THEN
    WRITE(*,*)'Invalid selection for O3:',io3
        WRITE(*,*)'Valid values: 0=no, 1=yes.'
        good_data = .false. 
      ENDIF
 
      IF((in2o.NE.0).AND.(in2o.NE.1)) THEN
    WRITE(*,*)'Invalid selection for N2O:',in2o
        WRITE(*,*)'Valid values: 0=no, 1=yes.'
        good_data = .false. 
      ENDIF
 
      IF((ico.NE.0).AND.(ico.NE.1)) THEN
    WRITE(*,*)'Invalid selection for CO:',ico
        WRITE(*,*)'Valid values: 0=no, 1=yes.'
        good_data = .false. 
      ENDIF
 
      IF((ich4.NE.0).AND.(ich4.NE.1)) THEN
    WRITE(*,*)'Invalid selection for CH4:',ich4
        WRITE(*,*)'Valid values: 0=no, 1=yes.'
        good_data = .false. 
      ENDIF
 
      IF((io2.NE.0).AND.(io2.NE.1)) THEN
    WRITE(*,*)'Invalid selection for O2:',io2
        WRITE(*,*)'Valid values: 0=no, 1=yes.'
        good_data = .false. 
      ENDIF
 
      IF((ino2.NE.0).AND.(ino2.NE.1)) THEN
    WRITE(*,*)'Invalid selection for NO2:',ino2
        WRITE(*,*)'Valid values: 0=no, 1=yes.'
        good_data = .false. 
      ENDIF
 
!C Read in the column ozone amount and scaling factor for NO2 column amount.
!C    The typical ozone amount is: 0.28-0.55 (atm-cm). The built-in NO2
!C    column amount is 5.0E+15 molecules/cm^2.
      READ(*,*)vrto3, sno2
      IF((vrto3.LT.0.1).OR.(vrto3.GT.0.6)) THEN
    WRITE(*,*)'Invalid vertical column amount of ozone:',vrto3
        WRITE(*,*)'Valid values are 0.1-0.6.'
        WRITE(*,*)'Reasonable values are 0.28-0.55.'
        good_data = .false. 
      ENDIF
 
!C Get aerosol model, visibility, or aerosol optical depth at 0.55 um (TAER55)
!C for SSSSS(). If V > 0, no TAER55 is needed in the input file. If V = 0,
!C TAER55 is needed.
!C IAER: 0=no aerosol, set V = -1
!C IAER: 1=continental aerosol, 2=maritime aerosol, 3=urban aerosol, set V>5 km.
 
      READ(*,*) iaer,v
      IF((iaer.LT.0).OR.(iaer.GT.6)) THEN
    WRITE(*,*)'Invalid aerosol model value:',iaer
        WRITE(*,*)'Valid values are 0-3.'
        good_data = .false. 
      ENDIF
 
      IF(iaer.EQ.0) THEN
        v = -1
        taer55=0.0
       ELSE
        IF((v.LT.0).OR.(v.GT.300)) THEN
      WRITE(*,*)'Invalid visibility value:',v
          WRITE(*,*)'Value must be greater than 0 and less than 300.'
          good_data = .false. 
       ENDIF
      ENDIF
 
      IF(v.EQ.0) THEN
        READ(*,*) taer55
        IF((taer55.GT.10.).OR.(taer55.LE.0)) THEN
      WRITE(*,*)'Invalid aerosol optical depth at 550 nm: ',taer55
          WRITE(*,*)'Valid values are between 0 and 10.'
          good_data = .false. 
        ENDIF
      ENDIF
 
!C Get the mean surface elevation.  Range: 0.0 - H(NB-1) km
      READ(*,*)hsurf 
!C--
      xpss = hsurf
!C--
      IF((hsurf.LT.0).OR.(hsurf.GT.h(nb-1))) THEN
    WRITE(*,*)'Invalid surface elevation value:',hsurf
        WRITE(*,*)'Value must be less than the maximum elevation in the atmospheric model.'
        good_data = .false. 
      ENDIF
 
!C Test to see if plane altitude is greater than bottom surface elevation
      IF(xppp .LT. hsurf) THEN
        WRITE(*,*)'Invalid plane altitude:',xppp
        WRITE(*,*)'Plane altitude must be greater than the bottom surface elevation.'
        good_data = .false.
      ENDIF
 
!C  Read in file name of data cube and get data dimensions in DIMS, ST_ORDER
!C  If no SIPS header on data, then default dimensions are header size = 0,
!C  number of samples = 614, number of lines = 512, number of bands = 224,
!C  storage order = BIL
      READ(*,85) finav
!C------      CALL OPENINFILE(IRET,DIMS,ST_ORDER)
      lun_in  = 47
!      CALL OPENINFILE(LUN_IN, I_RET)
 
!      IF(I_RET.NE.0) THEN
!        WRITE(*,*)'Input image file did not open
!     & successfully:',FINAV
!        GOOD_DATA = .FALSE.
!      ENDIF
 
!C Determine if cube dimensions should be read in.  If not, the dimensions
!C in the image file header will be used.
      READ(*,*) ans
      IF((ans.NE.0).AND.(ans.NE.1)) THEN
    WRITE(*,*)'Invalid value to indicate whether the cube dimensions should be read:',ans
        WRITE(*,*)'Valid values: 0=no 1=yes.'
        good_data = .false. 
      ELSE
        IF (ans .EQ. 1) THEN
!C         Read the header size in bytes, samples, lines, bands, and
!C           storage order (0=BSQ, 1=BIP, 2=BIL)'
          READ(*,*) hdrec, nsamps, nlines, nbands, sorder
          IF ((hdrec.LT.0).OR.(nsamps.LE.0).OR. (nlines.LE.0).OR.(nbands.LE.0).OR. &
             (sorder.LE.0).OR.(sorder.GT.2)) THEN
        WRITE(*,*)'Invalid cube parameters:',hdrec,nsamps,nlines, &
      nbands,sorder
        WRITE(*,*)'Values must be greater than or equal to 1 and the storage order must be less than 3.'
      WRITE(*,*)'Storage Order (0=BSQ, 1=BIP, 2=BIL)'
      WRITE(*,*)' 0 = BSQ is not supported in this version of ATREM'
            good_data = .false. 
          ENDIF
        ELSE
      hdrec  = dims(1)
          nsamps = dims(2)
          nlines = dims(3)
          nbands = dims(4)
      sorder = st_order
        ENDIF  
      ENDIF  
 
!C Get the output reflectance cube file name.  The output cube always has a 512
!C byte header record.
      READ(*,85) focub
!C---      CALL OPENOUTFILE(IRET,512,NSAMPS,NLINES,NBANDS,SORDER)
      lun_out = 48
!      CALL OPENOUTFILE(LUN_OUT, I_RET)
 
!      IF(I_RET.NE.0) THEN
!        WRITE(*,*)'Output image file did not open
!     & successfully:',FOCUB
!        GOOD_DATA = .FALSE.
!      ENDIF
 
!C Get the resolution of output spectra.  Format: 0-100 nm.
      READ(*,*)dlt2
      IF((dlt2.LT.0.).OR.(dlt2.GT.100.)) THEN
    WRITE(*,*)'Invalid output resolution value:',dlt2
        WRITE(*,*)'Valid values are 0-100 nm.'
        good_data = .false. 
      ENDIF
 
!C Get the scaling factor for output reflectance values.  Range: 1. to 32000.
      READ(*,*)scalef
      IF((scalef.LT.1.).OR.(scalef.GT.32000.)) THEN
    WRITE(*,*)'Invalid output resolution value:',scalef
        WRITE(*,*)'Valid values are 1-32000 nm.'
        good_data = .false. 
      ENDIF
 
!C Get the output water vapor file name.  The water vapor output has a 512-byte
!C header, is always one band, and is in BSQ format (BSQ=0).
      READ(*,85) foh2o
!C---      CALL OPENVAPFILE(IRET,512,NSAMPS,NLINES,1,0)
      lun_vap = 49
!      CALL OPENVAPFILE(LUN_VAP, I_RET)
 
!      IF(I_RET.NE.0) THEN
!        WRITE(*,*)'Output water vapor file did not open
!     & successfully:',FOH2O
!        GOOD_DATA = .FALSE.
!      ENDIF
 
!C  Get the file name for output table consisting of atmospheric
!C  transmission spectra at the solar and observational geometry but
!C  with varying amounts of water vapor.
 
!COMMENTED out RJH
!      READ(*,85) FOUT1
!      OPEN(35,FILE=FOUT1,STATUS='UNKNOWN',IOSTAT=ISTAT)
!      IF(ISTAT.NE.0) THEN
!        WRITE(*,*)'Output file did not open successfully:',FOUT1
!        GOOD_DATA = .FALSE.
!      ENDIF
!
!      IF(GOOD_DATA .NEQV. .TRUE.) THEN
!        WRITE(*,*)'**** ERROR: Invalid input detected. ****'
!        STOP
!      ENDIF
 
!C--D     WRITE(*,*)'IH2OVP=',IH2OVP,' ICO2=',ICO2,' IO3=',IO3,' IN2O=',IN2O
!C--D     WRITE(*,*)'ICO=',ICO,' ICH4=',ICH4,' IO2=',IO2
!C--D     WRITE(*,*)'H         T        P        VMR'
!C--D     DO 998 I=1,25
!C--D 998   WRITE(*,*)H(I),' ',T(I),' ',P(I),' ',VMR(I)
!C--D     WRITE(*,*)'NB=',NB,' NL=',NL,' MODEL=',MODEL
!C--D     WRITE(*,*)' IAER=',IAER,' V=',V,' VRTO3=',VRTO3
!C--D     DO 997 I=1,NOBS
!C--D 997   WRITE(*,*)'I=',I,' WAVOBS(I)=',WAVOBS(I),' FWHM(I)=',FWHM(I)
!C--D     WRITE(*,*)'NOBS=',NOBS,' HSURF=',HSURF,' DLT=',DLT,' DLT2=',DLT2
!C--D     WRITE(*,*)'SCALEF=',SCALEF
!C--D     WRITE(*,*)'WNDOW1=',WNDOW1,' WNDOW2=',WNDOW2,' WP94C=',WP94C
!C--D     WRITE(*,*)'WNDOW3=',WNDOW3,' WNDOW4=',WNDOW4,' W1P14C=',W1P14C
!C--D     WRITE(*,*)'NB1=',NB1,' NB2=',NB2,' NBP94=',NBP94
!C--D     WRITE(*,*)'NB3=',NB3,' NB4=',NB4,' NB1P14=',NB1P14
!C--D     WRITE(*,*)'IMN=',IMN,' IDY=',IDY,' IYR=',IYR
!C--D     WRITE(*,*)'IH=',IH,' IM=',IM,' IS=',IS
!C--D     WRITE(*,*)'XLATD=',XLATD,' XLATM=',XLATM,' XLATS=',XLATS
!C--D     WRITE(*,147)LATHEM
!C--D 147 FORMAT('LATHEM=',A1)
!C--D     WRITE(*,*)'XLONGD=',XLONGD,' XLONGM=',XLONGM,' XLONGS=',XLONGS
!C--D     WRITE(*,148)LNGHEM
!C--D 148 FORMAT('LNGHEM=',A1)
!C--D     WRITE(*,*)'HDREC=',HDREC,' NSAMPS=',NSAMPS,' NLINES=',NLINES
!C--D     WRITE(*,*)'NBANDS=',NBANDS,' SORDER=',SORDER
 
      RETURN
      END
 
!********************************************************************************
!*                                                 *
!*  Name: MODEL_ADJ                                *
!*  Purpose: resets the bottom boundary of the input model if the surface       *
!*           elevation is greater than 0, and calculate the column water vapor  *
!*           amount in the selected model.                     *
!*  Parameters: none                                   *
!*  Algorithm: If the surface elevation > 0, the bottom layer temperature and   *
!*             water vapor volume mixing ratio are obtained through linear      *
!*             interpolation, while the bottom layer pressure is obtained       *
!*             through exponential interpolation.                  *
!*  Globals used:  H, T, P, VMR, NB, NL, HSURF - these values are adjusted if   *
!*                      HSURF > 0                          *
!*  Global output:  CLMVAP - Column water vapor amount in unit of cm.           *
!*                       Q - Number of molecules above the surface at one       *
!*                           atmosphere in units of molecules/cm**2        *
!*  Return Codes: none                                 *
!*  Special Considerations: none                               *
!*                                         *
!********************************************************************************
      SUBROUTINE model_adj
 
      include 'COMMONS_INC.f'
 
!C  Common variables
      dimension h(25), t(25), p(25), vmr(25)
 
      COMMON /getinput3/ h,t,p,vmr,nb,nl,model,iaer,v,taer55,vrto3,sno2
      COMMON /getinput5/ nobs,hsurf,dlt,dlt2
      COMMON /model_adj1/ clmvap,q
 
      dimension hp(25), tp(25), pp(25), vmrp(25)
      COMMON /model_adj2/ hp, tp, pp, vmrp
      COMMON /model_adj3/ k_plane, dvap_plane, dvap_layer, &
                         dp_plane, dp_layer, clmvapp
      COMMON /model_adj4/ k_surf
 
      REAL XPSS, XPPP
      COMMON /getinput14/ xpss, xppp
 
!C     H = Layer boundary height, SL=Slant path length, DSLODH=DSL/DH
 
      re=6380.0
!C     Q=# of molecules above the surface at one atmosphere (molecules/cm**2)
      q=2.152e25
 
!C--        print*, 'original model atmosphere ...'
!C--      DO I = 1, 25
!C--        print*, 'H,T,P,VMR = ', H(I), T(I), P(I), VMR(I)
!C--      END DO
 
!C--- Convert the atmospheric pressure from "mb" to "atm.":
      DO i = 1, nb
         p(i) = p(i) / 1013.
      END DO
!C--- End of conversion
 
!C     Convert the VMR from the ppm unit in the model to absolute unit.
      DO 310 i=1,nb
  310   vmr(i)=vmr(i)*1.0e-06
 
!C     Do special processing if surface altitude is greater than 0
!C---      IF(HSURF.NE.0.0) THEN
!C       Reset H(I) to a smaller value if HSURF.EQ.H(I) to avoid possible
!C         problems in table searching using LOCATE
        DO 7455 i=1,nb
 7455     IF(hsurf.EQ.h(i)) hsurf=h(i)+0.0001
 
!C       Determine index of H() such that H(index) < HSURF < H(index+1)
        CALL locate(h,nb,hsurf,k)
!C
!C K_SURF is an index relative to the original model atmosphere (0 - 100 km)
        k_surf = k
!C---        print*, 'K_SURF = ', K_SURF
!C
        IF(k.EQ.0) THEN
          WRITE(6,5237) 
 5237     FORMAT(2x,'***WARNING: Surface elevation smaller then lowest boundary of the model atmosphere.')
 
          k_surf = 1
          GOTO 5255
        ENDIF
 
        IF(k.GT.0) THEN
          dhk=h(k+1)-h(k)
          dhs=hsurf -h(k)
 
!C         linear interpolation for surface temperature (TSURF) and VMR (  VMRS)
          tsurf=t(k)+(dhs/dhk)*(t(k+1)-t(k))
          vmrs =vmr(k)+(dhs/dhk)*(vmr(k+1)-vmr(k))
 
!C         exponential interpolation for surface pressure (PSURF)
          psurf=p(k)*exp(-alog(p(k)/p(k+1))*dhs/dhk)
          h(1)=hsurf
          p(1)=psurf
          t(1)=tsurf
          vmr(1)=vmrs
 
          nb=nb-k+1
          nl=nb-1
!C---          print*,'NL = ', NL, ' NB = ', NB, 'in Model_Adj'
!C
          DO 5240 i=2,nb
            h(i)=h(i+k-1)
            p(i)=p(i+k-1)
            t(i)=t(i+k-1)
            vmr(i)=vmr(i+k-1)
 5240     CONTINUE
!C  Zero out pressures and VMRS of top atmospheric layers.
          DO 5245 i=nb+1,25
            h(i)=1000.
            p(i)=0.0
            t(i)=300.
            vmr(i)=0.0
 5245     CONTINUE
        ENDIF
!C---      ENDIF
 
 5255 CONTINUE
 
      amtvrt=0.0
 
      DO 350 i=1,nl
        damtvt=q*(p(i)-p(i+1))*(vmr(i)+vmr(i+1))/2.0
        amtvrt=amtvrt+damtvt
 350  CONTINUE
 
      clmvap=amtvrt/3.34e+22
 
      WRITE(*,*)'Column vapor amount in model atmosphere from ground'
      WRITE(*,*)'       to space = ', clmvap, ' cm'
 
!C---      WRITE(*,*)'In MODEL_ADJ, NB = ',NB,' NL = ',NL
!C---        print*, 'After adjusting for elevated surface ...'
!C---      DO I = 1, 25
!C---        print*, 'H,T,P,VMR = ', H(I), T(I), P(I), VMR(I)
!C---      END DO
!
!C
!C
!C Setting the upward atmospheric path's T, P, and VMR profiles:
!C
!C  1st duplicate the entire atmospheric profiles from the downward path
!C      to the upward path
!C
      DO i = 1, 25
         hp(i)    = h(i)
         tp(i)    = t(i)
         pp(i)    = p(i)
         vmrp(i)  = vmr(i)
      END DO
 
 
      hplane = xppp
!C   Set the highest plane altitude to the upper bound of model atmosphere
!C---      IF(HPLANE.GT.100.0) HPLANE = 100. - 0.0001
      IF(hplane.GE.100.0) hplane = 100. - 0.0001
!C
!C  Do special processing if the plane height (HPLANE) is greater than HP(1)
      IF(hplane.GT.hp(1)) THEN
!C         Reset Plane altitude HPLANE (= XPPP) to a larger value if
!C             HPLANE.EQ.HP(I) to avoid possible problems in table
!C             searching using LOCATE
        DO 7456 i=1,25
 7456     IF(hplane.EQ.hp(i)) hplane=hp(i)-0.0001
 
!C       Determine index of HP() such that HP(index) < HPLANE < H(index+1)
        CALL locate(hp,nb,hplane,kk)
 
        IF(kk.EQ.0) THEN
          WRITE(6,5239) 
 5239     FORMAT(2x,'***WARNING: Plane altitude less then lowest boundary of the model atmosphere.')
          GOTO 5256
        ENDIF
 
        IF(kk.GT.0) THEN
          dhkk = hp(kk+1) - hp(kk)
          dhss = hplane   - hp(kk)
 
!C         linear interpolation for plane temperature (TPLANE) and VMR (  VMRSP)
          tplane = tp(kk)   + (dhss/dhkk)*(tp(kk+1)-tp(kk))
          vmrsp  = vmrp(kk) + (dhss/dhkk)*(vmrp(kk+1)-vmrp(kk))
 
!C         exponential interpolation for plane pressure (PPLANE)
          pplane     = pp(kk)*exp(-alog(pp(kk)/pp(kk+1))*dhss/dhkk)
          hp(kk+1)   = hplane
          pp(kk+1)   = pplane
          tp(kk+1)   = tplane
          vmrp(kk+1) = vmrsp
 
!C  Zero out pressures and VMRP of top atmospheric layers.
          IF(kk.LT.24) THEN
           DO i=kk+2,25
              hp(i)=1000.
              pp(i)=0.0
              tp(i)=300.
              vmrp(i)=0.0
           END DO 
          END IF
 
        ENDIF
      ENDIF
 
 5256 CONTINUE
 
      amtvrtp=0.0
 
      DO 357 i=1,kk
        damtvtp=q*(pp(i)-pp(i+1))*(vmrp(i)+vmrp(i+1))/2.0
        amtvrtp=amtvrtp+damtvtp
 357  CONTINUE
 
      clmvapp=amtvrtp/3.34e+22
 
      WRITE(*,*)'Column vapor below plane (CLMVAPP) = ', &
      clmvapp, ' cm'
 
!C---      WRITE(*,*)'In MODEL_ADJ, NB = ',NB,' KK = ', KK
!C---        print*, 'After further adjusting for plane height ...'
!C---      DO I = 1, 25
!C---        print*, 'HP,TP,PP,VMRP = ', HP(I), TP(I), PP(I), VMRP(I)
!C---      END DO
 
!C--- Indices and parameters for the plane layer
      k_plane = kk
 
      dvap_plane = q*(pp(k_plane) - pp(k_plane+1))* &
       (vmrp(k_plane) + vmrp(k_plane+1))/2.0 / 3.34e+22
 
      dvap_layer = q*(p(k_plane) - p(k_plane+1))*  &
       (vmr(k_plane) + vmr(k_plane+1))/2.0 / 3.34e+22
 
      dp_plane = pp(k_plane) - pp(k_plane+1)
      dp_layer = p(k_plane)  - p(k_plane+1)
 
!C---      print*, 'K_PLANE, DVAP_PLANE, DVAP_LAYER = ',
!C---     &         K_PLANE, DVAP_PLANE, DVAP_LAYER
!C---      print*, 'DP_PLANE, DP_LAYER = ', DP_PLANE, DP_LAYER
 
 
      RETURN
      END
!********************************************************************************
!*                                                 *
!*  Name: GEOMETRY                                 *
!*  Purpose: Calculates the solar and the observational geometric factors.      *
!*  Parameters: none                                   *
!*  Algorithm: The solar geometry was obtained based on the latitude, longitude,*
!*             GMT time using programs written by W. Mankin at National Center  *
!*             for Atmospheric Research in Boulder, Colorado. The          *
!*             geometric factors for CO2, O3, N2O, CO, CH4, and O2 are based    *
!*             only on the solar and observational angles. Sixty artificial     *
!*             geometric factors for H2O are set up to produce a transmittance  *
!*             table for different atmospheric water vapor amounts.        *
!*  Globals used:  VRTO3    - Column O3 amount in units of atm-cm              *
!*      IMN,IDY,IYR,IH,IM,IS - time and date of data measurements          *
!*      XLATD,XLATM,XLATS,LATHEM    - Latitude of measured area            *
!*      XLONGD,XLONGM,XLONGS,LNGHEM - Longitude of measured area               *
!*      CLMVAP - Column water vapor in unit of cm in the model atmosphere       *
!*  Global output:                                     *
!*      SOLZNI,SOLAZ,OBSZNI,OBSPHI,IDAY - Solar zenith angle, solar azimuth     *
!*            angle, observational zenith angle, observational azimuth angle,   *
!*            and the day number in the year                       *
!*      GCO2,GO3,GN2O,GCO,GCH4,GO2,SSH2O,GGEOM,TOTLO3 - Geometric factors for   *
!*            different gases, and total O3 amount in the sun-surface ray path. *
!*            The geometric factor is defined as: if the vertical column amount *
!*            of the gas is equal 1, then GGAS is equal to the total amount of  *
!*            the gas in the combined Sun-surface-sensor ray path.         *
!*  Return Codes: none                                 *
!*  Special Considerations: none                               *
!*                                         *
!********************************************************************************
 
      SUBROUTINE geometry
 
      include 'COMMONS_INC.f'
 
      dimension vapvrt(60), vap_slant(60)
      dimension md(12)
      dimension ssh2o(60)
      dimension h(25), t(25), p(25), vmr(25)
      CHARACTER*1 LATHEM,LNGHEM
 
!C The VAPVRT array contains 60 column water vapor values used for generating
!C  a table of transmittance spectra with different amount of water vapor.
!C  The values in VAPVRT is designed so that there is approximately 2% change
!C  in the .94-um H2O channel ratio for column water vapor in the range .4-6 cm.
 
      DATA vapvrt/.00, .02, .06,  .11, .16,  .21,  .26,  .31,  .36, .40, &
                 .43, .46, .50,  .54, .58,  .62,  .66,  .70,  .75, .80, &
                 .86, .92, .98, 1.06,1.14, 1.22,  1.3,  1.4,  1.5, 1.6, &
                 1.7, 1.8, 1.9, 2.05, 2.2, 2.35, 2.55, 2.75, 2.95, 3.2, &
                 3.5, 3.8, 4.1,  4.4, 4.7,  5.0,  5.3,  5.6,  6.0, 6.4, &
                 7.0, 7.7, 8.5,  9.4,10.4, 11.6, 13.0, 15.0, 25.0, 50./
 
      DATA md/0,31,59,90,120,151,181,212,243,273,304,334/
 
      COMMON /getinput3/ h,t,p,vmr,nb,nl,model,iaer,v,taer55,vrto3,sno2
 
      COMMON /getinput8/ imn,idy,iyr,ih,im,is
      COMMON /getinput9/ xlatd,xlatm,xlats,lathem
      COMMON /getinput10/xlongd,xlongm,xlongs,lnghem
      COMMON /model_adj1/ clmvap,q
      COMMON /geometry1/ solzni,solaz,obszni,obsphi,iday
      COMMON /geometry2/ gco2,go3,gn2o,gco,gch4,go2,ssh2o,totlo3,ggeom
 
      COMMON /model_adj3/ k_plane, dvap_plane, dvap_layer, &
                          dp_plane, dp_layer, clmvapp
 
      dimension g_vap(25), g_other(25)
      COMMON /geometry3/ g_vap, g_other, g_vap_equiv, vap_slant_mdl
 
      REAL XPSS, XPPP
      COMMON /getinput14/ xpss, xppp
 
      REAL XVIEWD, XVIEWM, XVIEWS
      REAL XAZMUD, XAZMUM, XAZMUS
      COMMON /getinput15/ xviewd,xviewm,xviews, xazmud,xazmum,xazmus
 
      hplane = xppp
!C
!C VAP_SLANT is a new array for containing two-way total vapor amounts
!C     Here the number "2" can be changed to other numbers, e.g., 2.5,
!C     without any major effects in derived water vapor values from
!C     imaging spectrometer data.
!C
 
      DO i = 1, 60
        vap_slant(i) = vapvrt(i) * 2.0
      END DO
 
      xh=ih
      xm=im
      xs=is
      tt=6.28318*((xh)/24+xm/1440+xs/86400)
 
      xlat = xlatd + xlatm/60.0 + xlats/3600.0
      xlong = xlongd + xlongm/60.0 + xlongs/3600.0
 
!C northern hemisphere is positive, and western hemisphere is positive
!C (NOTE: this is not standard, but has been internally coded in this program)
      IF(lathem.NE.'N') xlat = -xlat
      IF(lnghem.NE.'W') xlong = -xlong
 
325   xlatr = xlat/57.2958
      xlongr = xlong/57.2958
 
      CALL suncor(idy,imn,iyr,tt,dec,haz)
      CALL hazel(haz+tt-xlongr,dec,solaz,el,xlatr)
!C---Note: DEC, SOLAZ,and EL DEC are in radiant units
 
      IF (el .LE. 0.) THEN
        WRITE(*,*)'ERROR: Sun is below the horizon!!!'
        WRITE(*,*)'Check input date, time, latitude and longitude.'
        stop
      ENDIF
 
!C converting observational zenith angle and relative azimuth angle into degrees
      obszni = xviewd + xviewm/60.0 + xviews/3600.0
      obsphi = xazmud + xazmum/60.0 + xazmus/3600.0
 
      write(91,*) 'OBSZNI = ',obszni, ' OBSPHI = ',obsphi, ' degrees'
 
      obszni = obszni / 57.2958
      obsphi = obsphi / 57.2958
 
      solzni = 90.0/57.2958 - el
 
      write(91,*) 'solzni=',solzni
      write(91,*) 'KPLANE=',k_plane
 
      ggeom  = 1./cos(solzni) + 1./cos(obszni)
 
      gco2   = ggeom
 
      go3    = ggeom
      IF(hplane.LT.27.) go3 = ggeom - 1./cos(obszni)
 
      gn2o   = ggeom
      gco    = ggeom
      gch4   = ggeom
      go2    = ggeom
 
      write(91,*) 'GGEOM, GCO2, GO3, GN2O, GCO, GCH4, GO2 = '
      write(91,*)  ggeom, gco2, go3, gn2o, gco, gch4, go2
 
      totlo3 = go3 * vrto3
 
      WRITE(91,*) 'TOTLO3 = ', totlo3
 
!C Initialize newly created geometrical factors for each atmospheric
!C   layers (here G_VAP and G_OTHER are true geometrical factors that
!C   account for the actual Sun-surface-plane path lengths in the
!C   model atmosphere):
!C
!C---For layers below the plane layer---
      DO i = 1, k_plane - 1
         g_vap(i)   = ggeom
         g_other(i) = ggeom
      END DO
!C
!C---For layers above the plane layer
      DO i = k_plane + 1, 25
         g_vap(i)   = ggeom - 1./cos(obszni)
         g_other(i) = ggeom - 1./cos(obszni)
      END DO
!C
!C---Special treatment for the plane layer to take account the
!C     "shorter" upward path length
      g_vap(k_plane)   = ggeom - 1./cos(obszni)  &
                         + dvap_plane/dvap_layer/cos(obszni)
      g_other(k_plane) = ggeom - 1./cos(obszni)  &
                         + dp_plane/dp_layer/cos(obszni)
 
     write(91,*) ' G_VAP, G_OTHER, I ='
      DO i = 1, 25
         write(91,*) g_vap(i), g_other(i), i
      END DO
 
!C Calculate the water vapor SCALING factor relative to the total amount
!C    of water vapor in the model atmosphere in the L-shaped
!C    Sun-surface-plane ray path.
!C
      vap_slant_mdl = clmvap/cos(solzni) + clmvapp/cos(obszni)
      write(91,*) 'VAP_SLANT_MDL =', vap_slant_mdl, ' cm'
!C
!C The "equivalent" geometrical factor corresponding to the total
!C    slant vapor amount of VAP_SLANT_MDL':
!C
      g_vap_equiv = vap_slant_mdl / clmvap
       write(91,*) 'G_VAP_EQUIV = ', g_vap_equiv
 
      DO 310 i=1,60
        ssh2o(i) = vap_slant(i) / vap_slant_mdl
        write(91,*) 'SSH2O(I), I = ', ssh2o(i), i
  310 CONTINUE
 
!C Calculate the number of days that have passed in this year.  Take leap year
!C into account.
      iday = md(imn) + idy
      lpyr = iyr - (4 * (iyr/4))
      IF((lpyr.EQ.0).AND.(iday.GT.59).AND.(imn.NE.2)) iday = iday + 1
 
!C--D     WRITE(*,*)'SOLZNI=',SOLZNI,' SOLAZ=',SOLAZ,' OBSZNI=',OBSZNI
!C--D     WRITE(*,*)'OBSPHI=',OBSPHI,' IDAY = ',IDAY
!C--D     WRITE(*,*)'GCO2=',GCO2,' GO3=',GO3,' GN2O=',GN2O,' GCO=',GCO
!C--D     WRITE(*,*)'GCH4=',GCH4,' GO2=',GO2,' TOTLO3=',TOTLO3
!C--D     WRITE(*,*)'GGEOM=',GGEOM
!C--D     DO 311 I=1,60
!C--D 311   WRITE(*,*)'I=',I,' SSH2O(I)',SSH2O(I)
 
      RETURN
      END
 
!********************************************************************************
!*                                         *
!*  Name: INIT_SPECCAL                                 *
!*  Purpose: initialize global data for spectrum calculations.             *
!*  Parameters: none                                   *
!*  Algorithm: initialize data.                            *
!*  Globals used: AH2O, APH2O, BH2O, BPH2O, SODLT, SOGAM, O3CF - Band model     *
!*                             parameters for spectral calculations.            *
!*                WNDOW1, WNDOW2, WP94C, WNDOW3, WNDOW4, W1P14C - Center        *
!*                             positions of window and water vapor absorption   *
!*                             channels used in 3-channel ratio calculations.   *
!*                NB1, NB2, NBP94, NB3, NB4, NB1P14 - Number of narrow channels *
!*                             to form broader window and absorption channels.  *
!*  Global output:                                 *
!*     IH2OLQ,RLQAMT,NGASTT,NH2O,VSTART,VEND - Flag for including liquid           *
!*                     water, liquid water amount (cm), total number of gases   *
!*                     (typically 8), number of water vapor values, starting    *
!*                     and ending wavelengths in internal calculations.         *
!*     NO3PT,NCV,NCVHAF,NCVTOT,VMIN,ISTART,IEND  - Number of O3 abs. coef.     *
!*                     points, parameters for gaussian function and spectral    *
!*                     calculations.                           *
!*     ISTCAL,IEDCAL,DP,PM,TM,VMRM - Parameters for spectral calculations       *
!*     IST1,IED1,IST2,IED2,ISTP94,IEDP94 - 3-channel ratioing parameters    for    *
!*                     the 0.94-um water vapor band.                   *
!*     IST3,IED3,IST4,IED4,IST1P14,IED1P14 - 3-channel ratioing parameters for  *
!*                     the 1.14-um water vapor band.                   *
!*     WT1,WT2,WT3,WT4,JA - Relative weights for the four window channels       *
!*                     used in channel-ratioing calculations. JA is a          *
!*                     output parameter from a table searching routine.        *
!*     NCV2,NCVHF2,NCVTT2,ISTRT2,IEND2,FINST2 - Parameters for smoothing           *
!*                     output reflectance spectra.                         *
!*     NATOT,NBTOT,NCTOT,NDTOT - Number of channels for the four AVIRIS'           *
!*                     grating spectrometers (A, B, C, and D).             *
!*  Return Codes: None.                                *
!*  Special Considerations:  Some parameters may need to be fine-tuned.        *
!*                                         *
!********************************************************************************
!*  Notes about water vapor VMRS and related quantities:                   *
!*                                         *
!*     VAPVRT(60)   - a table containing 60 column vapor values (in unit of cm) *
!*                                         *
!*     VAP_SLANT(I) = VAPVRT(I) * 2.0, VAP_SLANT is a new table for containing  *
!*                    two-way total vapor amounts. Here the number "2" can be   *
!*                    changed to other numbers, e.g., 2.5, without major        *
!*                    effects on retrieved water vapor values.                  *
!*                                         *
!*     G_VAP(I = 1,..., NL) = true vapor geometric factor for each layer in     *
!*                    the model atmosphere (after adjusting for the elevated    *
!*                    surface.                                                  *
!*                                         *
!*     VMRM(I) = VMRM(I)*G_VAP(I). The VMRS are multiplied by the geometrical   *
!*                    factor. We can calculate the vapor transmittance on the   *
!*                    Sun-surface-sensor path by assuming a vertical path in    *
!*                    the model atmosphere with geometric-factor-adjusted VMRS. *
!*                                         *
!*     CLMVAP  = vertical column amount from ground to space in model atmosphere*
!*     CLMVAPP = vertical column amount from ground to aircraft or satellite    *
!*                    sensor in model atmosphere                                *
!*     Q       = 2.152E25 = # of molecules above the surface at one atmosphere  *
!*                    (in unit of  molecules/cm**2)                *
!*                                         *
!*     VAP_SLANT_MDL= CLMVAP/COS(SOLZNI) + CLMVAPP/COS(OBSZNI) = total amount   *
!*                    of water vapor in the model atmosphere in the L-shaped    *
!*                    Sun-surface-plane ray path.                               *
!*                                         *
!*     G_VAP_EQUIV  = VAP_SLANT_MDL / CLMVAP = the "equivalent" geometrical     *
!*                    factor corresponding to the total slant vapor amount      *
!*                    VAP_SLANT_MDL and the column vapor amount CLMVAP.         *
!*                                         *
!*     SSH2O(I) (I = 1, ..., 60) - a pure scaling factor relative to the total  *
!*                    slant vapor amount of VAP_SLANT_MDL, and                  *
!*           SSH2O(I) = VAP_SLANT(I) / VAP_SLANT_MDL               *
!*                                         *
!*     SH2O = one value of SSH2O(I). SH2O is used during generation of the      *
!*           look-up table.                            *
!*                                         *
!*     VAPTT  = VAP_SLANT_MDL*SH2O, is the absolute total vapor amount on the   *
!*                    L-shaped path corresponding to a spectrum stored in the   *
!*                    look-up table.                                       *
!*                                         *
!*     CLMWVP = 0.5*(VAPTTA+VAPTTB)/G_VAP_EQUIV, is the retrieved column water  *
!*                    vapor amount from imaging spectrometer data.         *
!********************************************************************************
      SUBROUTINE init_speccal
 
      include 'COMMONS_INC.f'
 
!C  Common variables
      dimension h(25), t(25), p(25), vmr(25)
      dimension ssh2o(60)
      dimension wavobs(1024),fwhm(1024)
      dimension dp(25), pm(25), tm(25), vmrm(25)
      dimension finst2(100)
 
!C  Local variables
 
      COMMON /getinput1/ ih2ovp,ico2,io3,in2o,ico,ich4,io2,ino2
      COMMON /getinput3/ h,t,p,vmr,nb,nl,model,iaer,v,taer55,vrto3,sno2
      COMMON /getinput4/ wavobs,fwhm
      COMMON /getinput5/ nobs,hsurf,dlt,dlt2
      COMMON /getinput6/ wndow1,wndow2,wp94c,wndow3,wndow4,w1p14c
      COMMON /getinput7/ nb1,nb2,nbp94,nb3,nb4,nb1p14
      COMMON /geometry2/ gco2,go3,gn2o,gco,gch4,go2,ssh2o,totlo3,ggeom
      COMMON /model_adj1/ clmvap,q
 
      COMMON /init_speccal3/ nh2o
      COMMON /init_speccal5/ dp,pm,tm,vmrm
      COMMON /init_speccal6/ ist1,ied1,ist2,ied2,istp94,iedp94
      COMMON /init_speccal7/ ist3,ied3,ist4,ied4,ist1p14,ied1p14
      COMMON /init_speccal8/ wt1,wt2,wt3,wt4,ja
 
      COMMON /init_speccal10/ ncv2,ncvhf2,ncvtt2,istrt2,iend2,finst2
      COMMON /init_speccal11/ natot,nbtot,nctot,ndtot
 
      dimension g_vap(25), g_other(25)
      COMMON /geometry3/ g_vap, g_other, g_vap_equiv, vap_slant_mdl
 
      dimension o3cf(5001)
      COMMON /o3cf_init1/ o3cf
 
      dimension tran_o3_std(5001)
      COMMON /init_speccal16/ tran_o3_std
 
      dimension rno2cf(5001)
      COMMON /no2cf_init1/ rno2cf
 
      dimension tran_no2_std(5001)
      COMMON /init_speccal17/ tran_no2_std
 
      COMMON /model_adj4/ k_surf
      INTEGER start(2)/1,1/
      INTEGER cnt(2)/NP_HI,19/
 
 
      nh2o = 60         !Number of column water vapor values
!C
!C Initialization for high resolution spectral calculations -
!C First initialize arrays to smooth high resolution spectra (0.05 cm-1) to
!C    medium resolution spectra (0.2 nm):
!C
!C*** Note: The array WAVNO_HI is not used in actual computing, the array
!C*          should be removed from COMMONS_INC and the following DO LOOP
!C*          of this program, the purpose of keeping WAVNO_HI now is to
!C*          check array indices and make sure the correctness of the
!C*          indices ****************************************************
      DO i = 1, np_hi
         wavno_hi(i)  = 3000. + float(i-1)*dwavno   ! Wavenumber (cm-1)of high 
                                                    !    resolution spectrum,
         tran_hi(i)   = 1.0                         !    3000 - 18000 cm-1
      END DO
!C---        print*,'WAVNO_HI ,1, 10000,NP_HI = ',WAVNO_HI(1),
!C---     &     WAVNO_HI(10000), WAVNO_HI(NP_HI)
!C
!C
      DO i = 1, np_med
         wavln_med(i) = vstart + float(i-1)*dwavln  ! Wavelength of medium 
                                                    ! resolution spectrum,
         tran_med(i)  = 1.0                         ! FWHM=.2 nm, .56-3.1 um.
      END DO
!C-----
!C---        print*,'WAVLN_MED ,1, 10000,NP_MED = ',WAVLN_MED(1),
!C---     &     WAVLN_MED(10000), WAVLN_MED(NP_MED)
!C-----
!C
!C  NOTE:  WAVLN_STD starts from 0.3 um, instead of 0.56 um
      DO i = 1, np_std
         wavln_std(i)  = 0.3   + float(i-1)*dwavln  ! Wavelength of medium 
                                                    ! resolution spectrum,
         tran_std(i)   = 1.0                        ! FWHM=.2 nm, .3-3.1 um.
      END DO
!C-----
!C---        print*,'WAVLN_STD ,1, 10000,NP_STD = ',WAVLN_STD(1),
!C---     &     WAVLN_STD(10000), WAVLN_STD(NP_STD)
!C-----
!C
!C Note: The grids of WAVNO_HI do not match the grids of 10000./WAVLN_MED.
!C       INDEX_MED is a specially designed index for finding close matches
!C       between the two kinds of grids.
!C
      DO i = 1, np_med  
         index_med(i) = ( (10000./wavln_med(i) - 3000.)/dwavno + 1.)
      END DO
!C-----
!C---        print*,'INDEX_MED ,1, 10000,NP_MED = ',INDEX_MED(1),
!C---     &     INDEX_MED(10000), INDEX_MED(NP_MED)
!C-----
!C
!C Note:     WAVLN_MED_INDEX(I) is very close to WAVLN_MED(I),
!C       and WAVLN_MED_INDEX(I) >= WAVLN_MED(I)
!C
      DO i = 1, np_med
         wavln_med_index(i) = 10000. /(float(index_med(i)-1)*dwavno &
                                     + 3000.)
      END DO
!C-----
!C---        print*,'WAVLN_MED_INDEX ,1, 10000,NP_MED = ',WAVLN_MED_INDEX(1),
!C---     &     WAVLN_MED_INDEX(10000), WAVLN_MED_INDEX(NP_MED)
!C-----
 
 
      DO i = 1, np_med
         fwhm_wavno(i) = 10000.*dlt_med &
                         /(wavln_med_index(i)*wavln_med_index(i))
      END DO
!C-----
!C---        print*,'FWHM_WAVNO ,1, 10000,NP_MED = ',FWHM_WAVNO(1),
!C---     &     FWHM_WAVNO(10000), FWHM_WAVNO(NP_MED)
!C-----
!C
 
      DO i = 1, np_med
         ncvhf_wavno(i) = ( facdlt * fwhm_wavno(i) / dwavno + 1.)
      END DO
!C-----
!C---        print*,'NCVHF_WAVNO ,1, 10000,NP_MED = ',NCVHF_WAVNO(1),
!C---     &     NCVHF_WAVNO(10000), NCVHF_WAVNO(NP_MED)
!C-----
 
!C Initialize arrays for smoothing medium resolution spectrum (DLT_MED = 0.2 nm,
!C         and point spacing DWAVLN = 0.0001 micron) to coarser spectral
!C         resolution data from imaging spectrometers.
 
      DO i = 1, nobs
         ncvhf(i) = ( facdlt * fwhm(i) / dwavln + 1.)
      END DO
!C
!C parameters and arrays to smooth output surface reflectance spectrum
      wavcv2=facdlt*dlt2
 
!C Find the largest value in the FWHM array, and use this value in calculation
!C    of indices for smoothing output reflectance spectra. This smoothing
!C    algorithm should work well with grating spectrometers having nearly
!C    constant spectral resolutions, but not so well for prism spectrometers
!C    having variable spectral resolution.
      dwvavr = fwhm(1)
 
      DO i = 2, nobs
         IF(dwvavr.LT.fwhm(i)) dwvavr = fwhm(i)
      END DO
 
      rncv2=wavcv2/dwvavr
      ncv2=rncv2
      ncvhf2=ncv2+1
      ncvtt2=2*ncv2+1
 
      cons2=dlt2*sqrt(3.1415926/const1)
 
      IF (dlt2 .NE. 0.0) THEN
        sumins=0.0
        DO 585 i=ncvhf2,ncvtt2
          finst2(i)=exp(-const1*(float(i-ncvhf2)*dwvavr/dlt2)**2.0)
          sumins=sumins+finst2(i)
  585   CONTINUE
 
        DO 590 i=1,ncvhf2-1
          finst2(i)=finst2(ncvtt2-i+1)
          sumins=sumins+finst2(i)
  590   CONTINUE
 
        sumins=sumins*dwvavr
 
        DO 595 i=1,ncvtt2
          finst2(i)=finst2(i)*dwvavr/sumins
  595   CONTINUE
      ENDIF
 
      istrt2=ncvhf2
      iend2=nobs-ncvhf2
 
!C  number of channels of the four AVIRIS spectrometers.  These are used
!C  in removing null AVIRIS radiance values in the overlap portions of two
!C  adjacent spectrometers.
      nchnla=32
      nchnlb=64
      nchnlc=64
      nchnld=64
 
      natot=nchnla
      nbtot=nchnla+nchnlb
      nctot=nchnla+nchnlb+nchnlc
      ndtot=nchnla+nchnlb+nchnlc+nchnld
 
!C Resetting window wavelength positions and calculating weights for
!C  window and absorption channels used in 3-channel ratioing.
      iwndw1=findmatch(wavobs,nobs,wndow1)
      iwndw2=findmatch(wavobs,nobs,wndow2)
 
      wndow1=wavobs(iwndw1)
      wndow2=wavobs(iwndw2)
 
      jj=mod(nb1,2)
      IF(jj.EQ.0) nb1=nb1+1
      kk=mod(nb2,2)
      IF(kk.EQ.0) nb2=nb2+1
      nb1haf=(nb1-1)/2
      nb2haf=(nb2-1)/2
 
      ist1=iwndw1-nb1haf
      ied1=iwndw1+nb1haf
      ist2=iwndw2-nb2haf
      ied2=iwndw2+nb2haf
 
      iwp94c=findmatch(wavobs,nobs,wp94c)
      wp94c=wavobs(iwp94c)
 
      ll=mod(nbp94,2)
      IF(ll.EQ.0) nbp94=nbp94+1
      nb3haf=(nbp94-1)/2
      istp94=iwp94c-nb3haf
      iedp94=iwp94c+nb3haf
 
      wt1=(wndow2-wp94c)/(wndow2-wndow1)
      wt2=(wp94c-wndow1)/(wndow2-wndow1)
 
      iwndw4=findmatch(wavobs,nobs,wndow3)
      iwndw5=findmatch(wavobs,nobs,wndow4)
 
      wndow3=wavobs(iwndw4)
      wndow4=wavobs(iwndw5)
 
      jj=mod(nb3,2)
      IF(jj.EQ.0) nb3=nb3+1
      kk=mod(nb4,2)
      IF(kk.EQ.0) nb4=nb4+1
 
      nb4haf=(nb3-1)/2
      nb5haf=(nb4-1)/2
 
      ist3=iwndw4-nb4haf
      ied3=iwndw4+nb4haf
      ist4=iwndw5-nb5haf
      ied4=iwndw5+nb5haf
      iw1p14c=findmatch(wavobs,nobs,w1p14c)
 
      w1p14c=wavobs(iw1p14c)
      ll=mod(nb1p14,2)
      IF(ll.EQ.0) nb1p14=nb1p14+1
      nb6haf=(nb1p14-1)/2
      ist1p14=iw1p14c-nb6haf
      ied1p14=iw1p14c+nb6haf
 
      wt3=(wndow4-w1p14c)/(wndow4-wndow3)
      wt4=(w1p14c-wndow3)/(wndow4-wndow3)
 
!C Initialization for searching water vapor table (TRNTBL)
      ja = 30
!C
!C
!C Calculate medium resolution O3 transmittances (0.2 nm, 2-way) with
!C     a point spacing of 0.1 nm between 0.3 and 0.8 micron.
!C
      IF(io3.EQ.1) THEN
         DO i=1,no3pt
            tran_o3_std(i) = exp(-totlo3*o3cf(i))
         END DO
      END IF
 
!C
!C If O3 is not intended to be included in total atmospheric gaseous
!C     transmittance calculations, assigning TRAN_O3_STD = 1.0:
!C
      IF(io3.NE.1) THEN
         DO i=1,no3pt
            tran_o3_std(i) = 1.0
         END DO
      END IF
!C
!C Calculate medium resolution NO2 transmittances (0.2 nm, 2-way) with
!C     a point spacing of 0.1 nm between 0.3 and 0.8 micron.
!C
         no2pt  = no3pt
         vrtno2 = 5.0e+15
         vrtno2 = sno2 * vrtno2
 
         gno2   = go3
         totno2 = gno2 * vrtno2
 
      IF(ino2.EQ.1) THEN
         DO i=1,no2pt
            tran_no2_std(i) = exp(-totno2*rno2cf(i))
         END DO
      END IF
 
!C---Temp Code:
!C---      OPEN(57,file='zzzzzz_tst_NO2_trn',status='unknown')
!C---         DO I = 1, NO2PT
!C---          write(57,*)WAVLN_STD(I), TRAN_NO2_STD(I)
!C---         END DO
!C---      CLOSE(57)
!C---End of Temp Code
!
!C---      print*,' TOTNO2 = ', TOTNO2
!C---      print*,'VRTNO2, SNO2, GNO2 = ', VRTNO2, SNO2, GNO2
 
!C
!C If NO2 is not intended to be included in total atmospheric gaseous
!C     transmittance calculations, assigning TRAN_NO2_STD = 1.0:
!C
      IF(ino2.NE.1) THEN
         DO i=1,no2pt
            tran_no2_std(i) = 1.0
         END DO
      END IF
 
 
!C
!C Initial arrays for mean layer pressure and temperatures
 
      DO 320 i=1,nl
        dp(i)=p(i)-p(i+1)
        pm(i)=(p(i)+p(i+1))/2.0
        tm(i)=(t(i)+t(i+1))/2.0
  320 CONTINUE
 
!C
!C Calculate high resolution transmittances (0.05 cm-1) of CO2, N2O, CO,
!C     CH4, and O2 in the 0.56 - 3.1 micron range, and save values for
!C     calculating total atmospheric transmittances later.
!C     Because water vapor amounts are allowed to vary,
!C     the high resolution water vapor transmittances are calculated
!C     in subroutines TRAN_TABLE and TRANCAL. TRAN_TABLE provides variable
!C     water vapor amounts, and calls TRANCAL for the calculation of
!C     corresponding vapor transmittance spectrum.
!C
!C*** Note: The array WAVNO_HI is not used in actual computing, the array
!C*          should be removed from COMMONS_INC and the following DO LOOP
!C*          of this program, the purpose of keeping WAVNO_HI now is to
!C*          check array indices and make sure the correctness of the
!C*          indices.
!C*
!C Initialize the TRAN_HI_OTHERS array for high resolution spectrum:
        DO i = 1, np_hi
           tran_hi_others(i) = 1.0
        END DO
 
       ncid = ncopn('abscf_gas.nc',ncnowrit,ircode)
 
!C---------------------------------------
!C For CO2 transmittance calculation -
        IF(ico2.EQ.1) THEN
!C----        SCLCO2=1.0
!C----  On 2/7/2013 B.-C. Gao made the modification - Increased SCLCO2 i
!C----     from 1.0 to 1.1 to reflect the fact that the CO2 VMR reached the
!C----     2012 level of 390 ppmv.
!C----
          sclco2=1.1
 
          DO 322 i=1,nl
            vmrm(i)=sclco2*355.0*1.0e-06
!C Scale the VMRM by the two-way path geometrical factors. The geometric
!C           factors, G_OTHER, varies with atmospheric layer number for
!C           aircraft observational geometries.
            vmrm(i)= vmrm(i)*g_other(i)
  322     CONTINUE
  
  
       
!        NRHID = NCVID (NCID, 'waveno_hi_co2', IRCODE)
!       CALL NCVGT (NCID, NRHID, START(1), CNT(1), WAVNO_HI, IRCODE)
!       if (Ircode .ne.0) then
!          write(*,*) 'Error reading abscf_gas.nc: wavno_hi: rcode=',ircode
!          stop 0
!       end if
       nrhid = ncvid(ncid, 'abscf_co2', ircode)
       CALL ncvgt (ncid, nrhid, start, cnt, abscf, ircode)
        if (ircode .ne.0) then
          write(*,*) 'Error reading abscf_gas.nc: abscf_co2: rcode=',ircode
          stop 0
       end if
 
!       OPEN(31,file='abscf_co2_PC',status='old', &
!         form='unformatted',access='direct',recl=4*300000)
 
!          READ(31,REC=1) (WAVNO_HI(I), I = 1, 300000)
!
       DO j = k_surf, 19
!          READ(31,REC=J+1) (A(I), I = 1, 300000)
!
             tran_hi_others(:) = tran_hi_others(:) * exp( - abscf(:,j)*q* &
               dp(j-k_surf+1)*vmrm(j-k_surf+1) * 28.966 / &
               6.0225e+23 / 1.0e-06)
       END DO
!
!       CLOSE(31)
      ENDIF
 
!C--------------------------------------------
!C       For N2O transmittance calculation.
 
        IF(in2o.EQ.1) THEN
          DO 324 i=1,nl
            vmrm(i)=0.3*1.0e-06
            vmrm(i)= vmrm(i)*g_other(i)
  324     CONTINUE
  
!        NRHID = NCVID (NCID, 'waveno_hi_n2o', IRCODE)
!       CALL NCVGT (NCID, NRHID, START(1), CNT(1), WAVNO_HI, IRCODE)
!       if (Ircode .ne.0) then
!          write(*,*) 'Error reading abscf_gas.nc: wavno_hi: rcode=',ircode
!          stop 0
!       end if
       nrhid = ncvid(ncid, 'abscf_n2o', ircode)
       CALL ncvgt (ncid, nrhid, start, cnt, abscf, ircode)
        if (ircode .ne.0) then
          write(*,*) 'Error reading abscf_gas.nc: abscf_n2o: rcode=',ircode
          stop 0
       end if
 
!       OPEN(31,file='abscf_n2o_PC',status='old', &
!          form='unformatted',access='direct',recl=4*300000)
 
!          READ(31,REC=1) (WAVNO_HI(I), I = 1, 300000)
        
       DO j = k_surf, 19
!          READ(31,REC=J+1) (A(I), I = 1, 300000)
 
             tran_hi_others(:) = tran_hi_others(:) * exp( - abscf(:,j)*q* &
                dp(j-k_surf+1)*vmrm(j-k_surf+1) * 28.966 / &
                6.0225e+23 / 1.0e-06)
       END DO
 
!       CLOSE(31)
      ENDIF
 
!C--------------------------------------------
!C       For CO transmittance calculation.
        IF(ico.EQ.1) THEN
          DO 325 i=1,nl
            vmrm(i)=0.1*1.0e-06
            vmrm(i)= vmrm(i)*g_other(i)
  325     CONTINUE
  
!        NRHID = NCVID (NCID, 'waveno_hi_co', IRCODE)
!       CALL NCVGT (NCID, NRHID, START(1), CNT(1), WAVNO_HI, IRCODE)
!       if (Ircode .ne.0) then
!          write(*,*) 'Error reading abscf_gas.nc: wavno_hi: rcode=',ircode
!          stop 0
!       end if
       nrhid = ncvid(ncid, 'abscf_co', ircode)
       CALL ncvgt (ncid, nrhid, start, cnt, abscf, ircode)
        if (ircode .ne.0) then
          write(*,*) 'Error reading abscf_gas.nc: abscf_co: rcode=',ircode
          stop 0
       end if
 
!       OPEN(31,file='abscf_co_PC',status='old', &
!         form='unformatted',access='direct',recl=4*300000)
 
 !         READ(31,REC=1) (WAVNO_HI(I), I = 1, 300000)
 
       DO j = k_surf, 19
 !         READ(31,REC=J+1) (A(I), I = 1, 300000)
 
             tran_hi_others(:) = tran_hi_others(:) * exp( - abscf(:,j)*q* &
               dp(j-k_surf+1)*vmrm(j-k_surf+1) * 28.966 / &
               6.0225e+23 / 1.0e-06)
       END DO
 
 !      CLOSE(31)
      ENDIF
 
!C--------------------------------------------
!C       For CH4 transmittance calculation.
!C For assigning CH4 VMRM
!C NOTE: The scaling factor of 0.8 for the CH4 VMRS was obtained by comparing
!C       transmittance spectra calculated using our program, which is based on
!C       the Malkmus narrow band spectral model, with a ratioed spectrum
!C       provided by G. C. Toon at Jet Propulsion Laboratory (JPL). The JPL
!C       ratio spectrum was obtained by ratioing a high resolution (0.005 cm-1)
!C       solar spectrum measured at ground level against a solar spectrum
!C       measured at an altitude of approximately 35 km with the same Fourier
!C       Transform spectrometer. The high resolution ratio spectrum was
!C       degraded to a resolution of 10 nm during our derivation of the
!C       scaling factor for the CH4 VMRS.
 
        IF(ich4.EQ.1) THEN
!C***          SCLCH4=0.8
          sclch4=1.0
          DO 326 i=1,nl
            vmrm(i)=sclch4*1.6*1.0e-06
            vmrm(i)= vmrm(i)*g_other(i)
  326     CONTINUE
 
!        NRHID = NCVID (NCID, 'waveno_hi_ch4', IRCODE)
!       CALL NCVGT (NCID, NRHID, START(1), CNT(1), WAVNO_HI, IRCODE)
!       if (Ircode .ne.0) then
!          write(*,*) 'Error reading abscf_gas.nc: wavno_hi: rcode=',ircode
!          stop 0
!       end if
       nrhid = ncvid(ncid, 'abscf_ch4', ircode)
       CALL ncvgt (ncid, nrhid, start, cnt, abscf, ircode)
        if (ircode .ne.0) then
          write(*,*) 'Error reading abscf_gas.nc: abscf_ch4: rcode=',ircode
          stop 0
       end if
 
!       OPEN(31,file='abscf_ch4_PC',status='old', &
!          form='unformatted',access='direct',recl=4*300000)
 
!          READ(31,REC=1) (WAVNO_HI(I), I = 1, 300000)
 
       DO j = k_surf, 19
!          READ(31,REC=J+1) (A(I), I = 1, 300000)
 
              tran_hi_others(:) = tran_hi_others(:) * exp( - abscf(:,j)*q* &
               dp(j-k_surf+1)*vmrm(j-k_surf+1) * 28.966 / &
               6.0225e+23 / 1.0e-06)
       END DO
 
!       CLOSE(31)
      ENDIF
 
!C--------------------------------------------
!C       For O2 transmittance calculation.
        IF(io2.EQ.1) THEN
 
!C***Modified by Bo-Cai Gao on 2/7/2013 to increase O2 absorption
!C---  coefficients by the factor SCL_O2 for wavelengths > 1.2 micron
!C---  in order to model properly the atmospheric O2 band centered
!C---  near 1.265 micron.
 
          scl_o2  = 2.60
 
          DO 327 i=1,nl
            vmrm(i)=0.21
            vmrm(i)= vmrm(i)*g_other(i)
  327     CONTINUE
 
!        NRHID = NCVID (NCID, 'waveno_hi_o2', IRCODE)
!       CALL NCVGT (NCID, NRHID, START(1), CNT(1), WAVNO_HI, IRCODE)
!       if (Ircode .ne.0) then
!          write(*,*) 'Error reading abscf_gas.nc: wavno_hi: rcode=',ircode
!          stop 0
!       end if
       nrhid = ncvid(ncid, 'abscf_o2', ircode)
       CALL ncvgt (ncid, nrhid, start, cnt, abscf, ircode)
        if (ircode .ne.0) then
          write(*,*) 'Error reading abscf_gas.nc: abscf_o2: rcode=',ircode
          stop 0
       end if
 
!       OPEN(31,file='abscf_o2_PC',status='old', &
!         form='unformatted',access='direct',recl=4*300000)
 
!          READ(31,REC=1) (WAVNO_HI(I), I = 1, 300000)
 
       DO j = k_surf, 19
!          READ(31,REC=J+1) (A(I), I = 1, 300000)
 
!C***Modified by Bo-Cai Gao on 2/7/2013 to increase O2 absorption
!C---  coefficients by the factor SCL_O2 for wavelengths > 1.2 micron
!C---  i& < 1.3333 micron in order to model properly the atmospheric
!C---  O2 band centered near 1.265 micron.
 
             tran_hi_others(1:9000) = tran_hi_others(1:9000) * exp( - abscf(1:9000,j)*q* &
               dp(j-k_surf+1)*vmrm(j-k_surf+1) * 28.966 / &
               6.0225e+23 / 1.0e-06)
 
             tran_hi_others(9001:106600) = tran_hi_others(9001:106600) * exp( - abscf(9001:106600,j)*q* &
             dp(j-k_surf+1)*vmrm(j-k_surf+1)*scl_o2 *28.966/ &
               6.0225e+23 / 1.0e-06)
 
             tran_hi_others(106601:np_hi) = tran_hi_others(106601:np_hi) * exp( - abscf(106601:np_hi,j)*q* &
               dp(j-k_surf+1)*vmrm(j-k_surf+1) * 28.966 / &
               6.0225e+23 / 1.0e-06)
       END DO
 
       CLOSE(31)
      ENDIF
 
       CALL ncclos(ncid, rcode)
 
!C--------------------------------------------
!C End of calculation of high resolution transmittances for CO2, N2O, CO,
!C     CH4, and O2.
!C--------------------------------------------
!C Initial water vapor VMRs for repeated use in other subroutines and
!C    adjust layered water vapor VMRM with geometrical factors.
      DO i = 1, nl
         vmrm(i) = (vmr(i)+vmr(i+1))/2.0
!C Scale the VMRM by the two-way path geometrical factors. The geometric
!C           factors, G_VAP, varies with atmospheric layer number for
!C           aircraft observational geometries.
         vmrm(i) = vmrm(i)*g_vap(i)
      END DO
!C--------------------------------------------
!C
!C--D     WRITE(*,*)'NPSHIF=',NPSHIF,' DWAVLN=',DWAVLN
!C--D     WRITE(*,*)'NO3PT=',NO3PT,' VMIN=',VMIN,' ISTART=',ISTART
!C--D     WRITE(*,*)'IH2OLQ=',IH2OLQ,' RLQAMT=',RLQAMT,' NGASTT=',NGASTT
!C--D     WRITE(*,*)'NH2O=',NH2O,' VSTART=',VSTART,' VEND=',VEND
!C--D     WRITE(*,*)'IEND=',IEND
!C--D     WRITE(*,*)'ISTCAL=',ISTCAL,' IEDCAL=',IEDCAL
!C--D     DO 545 I=1,NL
!C--D 545   WRITE(*,*)I,DP(I),PM(I),TM(I),VMRM(I)
!C--D     WRITE(*,*)'IST1=',IST1,' IED1=',IED1,' IST2=',IST2,' IED2=',IED2
!C--D     WRITE(*,*)'ISTP94=',ISTP94,' IEDP94=',IEDP94
!C--D     WRITE(*,*)'IST3=',IST3,' IED3=',IED3,' IST4=',IST4,' IED4=',IED4
!C--D     WRITE(*,*)'IST1P14=',IST1P14,' IED1P14=',IED1P14
!C--D     WRITE(*,*)'WT1=',WT1,' WT2=',WT2,' WT3=',WT3,' WT4=',WT4,' JA=',JA
!C--D     WRITE(*,*)'NCV2=',NCV2,' NCVHF2=',NCVHF2,' NCVTT2=',NCVTT2,
!C--D    &' ISTRT2=',ISTRT2,' IEND2=',IEND2
!C--D     DO 544 I=1,30
!C--D 544   WRITE(*,*)'I=',I,' FINST2(I)=',FINST2(I)
!C--D     WRITE(*,*)'NATOT=',NATOT,' NBTOT=',NBTOT
!C--D     WRITE(*,*)'NCTOT=',NCTOT,' NDTOT=',NDTOT
     
      RETURN
      END
 
!********************************************************************************
!*                                             *
!*  Name: TRAN_TABLE                                   *
!*  Purpose: This subroutine generates a table consisting of 60 atmospheric     *
!*           transmittance spectra at the solar and observational               *
!*           geometry and with 60 column water vapor values. The table also     *
!*           includes the total amounts of column water vapor used in the       *
!*           calculations, and the 3-channel ratios calculated from the window  *
!*           and absorption channels in and around the 0.94- and 1.14-um water  *
!*           vapor bands.                              *
!*  Parameters: none                                   *
!*  Algorithm: For each of the 60 water vapor amounts, calculate the           *
!*             atmospheric transmittance, and save                 *
!*  Globals Used: NH2O   - number of column water vapor values             *
!*       VAPTT  - geometrically adjusted water vapor total         *
!*       R094   - channel ratio for .94 um region              *
!*       R114   - channel ratio for 1.14 um region             *
!*       TRNCAL - atmospheric transmittance spectra            *
!*  Global Output: VAPTOT() - array containing geometrically adjusted water     *
!*               vapor values                      *
!*        ROP94()  - array containing channel ratios for .94 um region *
!*        R1P14()  - array containing channel ratios for 1.14 um region*
!*        TRNTBL() - 2 dimensional array containing one transmittance  *
!*               spectrum for each column water vapor amount       *
!*  Return Codes: none                                                          *
!*  Special Considerations: none                                                *
!*                                         *
!********************************************************************************
      SUBROUTINE tran_table
 
      include 'COMMONS_INC.f'
 
!C  Common variables
      dimension trncal(1024)
      dimension wavobs(1024),fwhm(1024)
      dimension dp(25), pm(25), tm(25), vmrm(25)
      dimension ssh2o(60)
      dimension vaptot(60), r0p94(60), r1p14(60), trntbl(1024,60)
      
      dimension o3cf(5001)
      COMMON /o3cf_init1/ o3cf
 
      COMMON /getinput1/ ih2ovp,ico2,io3,in2o,ico,ich4,io2,ino2
      COMMON /getinput4/ wavobs,fwhm
      COMMON /getinput5/ nobs,hsurf,dlt,dlt2
      COMMON /init_speccal3/ nh2o
      COMMON /init_speccal5/ dp,pm,tm,vmrm
 
      COMMON /tran_table1/ sh2o,vaptot,r0p94,r1p14,trntbl
      COMMON /trancal1/ trncal,vaptt
      COMMON /geometry2/ gco2,go3,gn2o,gco,gch4,go2,ssh2o,totlo3,ggeom
      COMMON /chnlratio1/ r094,r114
 
!C For each water vapor amount, calculate the geometrically adjusted water
!C     vapor amount, the channel ratios, and the transmittance spectrum.
      vaptt = 0.
 
      DO 530 i=1,nh2o
        sh2o=ssh2o(i)
 
        CALL trancal
 
        vaptot(i)= vaptt
        r0p94(i) = r094
        r1p14(i) = r114
        DO 540 j=1,nobs
          trntbl(j,i)=trncal(j)
 540    CONTINUE
 530  CONTINUE
 
!C Write above calculated values to a file.
      nj=nh2o
      nk=float(nj)/10.
      njdiff=nj-nk*10
 
!RJH      DO 380 KK=1,NK
!        WRITE(35,77) (VAPTOT(I+(KK-1)*10),I=1,10)
!   77   FORMAT(7X,10(1X,E11.4))
!        WRITE(35,77) ( R0P94(I+(KK-1)*10),I=1,10)
!        WRITE(35,77) ( R1P14(I+(KK-1)*10),I=1,10)
!        DO 385 II=1,NOBS
!  385     WRITE(35,78) WAVOBS(II),(TRNTBL(II,I+(KK-1)*10),I=1,10)
!   78     FORMAT(1X,F6.4,10(1X,E11.4))
!  380 CONTINUE
 
!      WRITE(35,77) (VAPTOT(I+NK*10),I=1,NJDIFF)
!      WRITE(35,77) ( R0P94(I+NK*10),I=1,NJDIFF)
!      WRITE(35,77) ( R1P14(I+NK*10),I=1,NJDIFF)
!
!      IF(NJDIFF.GE.1) THEN
!        DO 386 II=1,NOBS
!  386     WRITE(35,78) WAVOBS(II),(TRNTBL(II,I+NK*10),I=1,NJDIFF)
!      ENDIF
!      CLOSE(35,IOSTAT=ISTAT)
 
      RETURN
      END
 
!********************************************************************************
!*                                         *
!*  Name: TRANCAL                                  *
!*  Purpose: This program calculates combined transmittances of H2O, CO2, O3,   *
!*      N2O, CO, CH4, and O2.                                                   *
!*  Parameters: none.                                  *
!*  Algorithm: The calculations were based on the line-by-line absorption       *
!*      parameters supplied by William R. Ridgway of NASA/GSFC.            *
!*  Global output:VAPTT  - geometrically adjusted water vapor total.           *
!*       R094   - channel ratio for 0.94 um region.            *
!*       R114   - channel ratio for 1.14 um region.            *
!*       TRNCAL - total transmittances of all gases that match the     *
!*                         resolutions of imaging spectrometers.               *
!*  Return Codes: none.                                *
!*  Special Considerations: The high resolution (0.05 cm-1) line-by-line        *
!*      absorption parameters cover the 0.555 - 3.33 micron spectral range      *
!*      (3000 - 18000 cm-1). The medium resolution ozone absorption             *
!*      coefficients covers the 0.3-0.8 um spectral range. The line-by-line     *
!*      high resolution spectra were first smoothed to medium resolution        *
!*      spectra (resolution = 0.2 nm, wavelength spacing = 0.1 nm) covering     *
!*      the 0.56 - 3.1 micron spectral region. The medium resolution spectra    *
!*      of O3 and other gases are combined (in SUBROUTINE TRAN_SMOOTH) to form  *
!*      a single medium resolution spectrum from 0.3 to 3.1 micron. This        *
!*      combined spectrum (medium resolution) is then smoothed to lower         *
!*      resolutions to match the resolutions of imaging spectrometers. The      *
!*      smoothing is also done in SUBROUTINE TRAN_SMOOTH.                       *
!*                                         *
!********************************************************************************
 
      SUBROUTINE trancal
 
      include 'COMMONS_INC.f'
      include 'netcdf.inc'
!C  Common variables
      dimension trncal(1024)
      dimension wavobs(1024),fwhm(1024)
      dimension dp(25), pm(25), tm(25), vmrm(25)
      dimension vaptot(60), r0p94(60), r1p14(60), trntbl(1024,60)
 
!C  Local variables
      dimension trans(1024)
      dimension trncv(1024)
      dimension trnstd_sm(1050)
      INTEGER   INDEX
      
      COMMON /getinput4/ wavobs,fwhm
      COMMON /getinput5/ nobs,hsurf,dlt,dlt2
      COMMON /model_adj1/ clmvap,q
      COMMON /init_speccal5/ dp,pm,tm,vmrm
      COMMON /init_speccal6/ ist1,ied1,ist2,ied2,istp94,iedp94
      COMMON /init_speccal7/ ist3,ied3,ist4,ied4,ist1p14,ied1p14
      COMMON /tran_table1/ sh2o,vaptot,r0p94,r1p14,trntbl
      COMMON /trancal1/ trncal,vaptt
 
      dimension g_vap(25), g_other(25)
      COMMON /geometry3/ g_vap, g_other, g_vap_equiv, vap_slant_mdl
 
      COMMON /model_adj4/ k_surf
      INTEGER FIRST/1/
      INTEGER start(2)/1,1/
      INTEGER cnt(2)/NP_HI,19/
      
!      REAL ABSCF(NP_HI,19)
      SAVE first
 
!C
!C Calculate water vapor transmittances with different scaling factors:
!C
 
       DO i = 1, np_hi
          tran_hi(i) = 1.0
       END DO
 
      if (first.eq.1) then
       ncid = ncopn('abscf_gas.nc',ncnowrit,ircode)
!        NRHID = NCVID (NCID, 'waveno_hi_h2o', IRCODE)
!       CALL NCVGT (NCID, NRHID, START(1), CNT(1), WAVNO_HI, IRCODE)
!       if (Ircode .ne.0) then
!          write(*,*) 'Error reading abscf_gas.nc: wavno_hi: rcode=',ircode
!          stop 0
!       end if
       nrhid = ncvid(ncid, 'abscf_h2o', ircode)
       CALL ncvgt (ncid, nrhid, start, cnt, abscf, ircode)
        if (ircode .ne.0) then
          write(*,*) 'Error reading abscf_gas.nc: abscf_h2o: rcode=',ircode
          stop 0
       end if
       CALL ncclos(ncid, rcode)
        
!       OPEN(31,file='abscf_h2o_PC',status='old', &
!         form='unformatted',access='direct',recl=4*300000)
 
!          READ(31,REC=1) (WAVNO_HI(I), I = 1, 300000)
 
!       DO J = 1, 19
!          READ(31,REC=J+1) (ABSCF(I,J), I = 1, 300000)
 
!       END DO
 
!       CLOSE(31)
       first=0
      endif
 
       DO j = k_surf, 19
 
             tran_hi(:) = tran_hi(:) * exp( - abscf(:,j) * q * &
               dp(j-k_surf+1)*vmrm(j-k_surf+1) * sh2o * 28.966 / &
               6.0225e+23 / 1.0e-06)
 
       END DO
 
!C Multiplying the high resolution water vapor transmittance with combined
!C    transmittances of CO2, N2O, CO, CH4, and O2:
!C
       tran_hi(:) = tran_hi(:)*tran_hi_others(:)
!       DO I = 1, NP_HI
!          TRAN_HI(I) = TRAN_HI(I) * TRAN_HI_OTHERS(I)
!       END DO
 
!C
!C Smooth the high resolution spectra to resolutions of measured spectrum
!C
      CALL tran_smooth
!C
!C---      DO 470 I=1,NOBS
!C---        WRITE(*,*)'I=',I,' TRNCAL(I)=',TRNCAL(I)
!C---  470 CONTINUE
 
!C Calculate 3-channel ratio values, R094 and R114, from simulated spectra.
      CALL chnlratio
 
!C Total amount of water vapor (in unit of cm) corresponding to the spectrum.
      vaptt = vap_slant_mdl * sh2o
!C--       print*, 'VAPTT= ', VAPTT,'VAP_SLANT_MDL= ',VAP_SLANT_MDL
 
 9999 RETURN
      END
 
!********************************************************************************
!*                                                                              *
!*  Name: TRAN_SMOOTH                                                           *
!*  Purpose: This program is to smooth the line-by-line high resolution         *
!*           spectrum to lower resolution spectrum that matches the resolutions *
!*           of imaging spectrometer data.                                      *
!*  Parameters: none.                                                           *
!*  Algorithm: The smoothing is done in two stages. The 1st stage is to smooth  *
!*             the high resolution spectrum to medium resolution spectrum at a  *
!*             constant FWHM (0.2 nm) and a constant wavelength interval        *
!*             (0.1 nm). The 2nd stage smoothing is to smooth the medium        *
!*             resolution spectrum to resolutions of input imaging spectrometer *
!*             data.                                                            *
!*  Globals used:  The global variables used are contained in the file          *
!*                         "COMMONS_INC"                                        *
!*  Global output: TRNCAL - total transmittances of all gases that match the    *
!*                          resolutions of imaging spectrometers.               *
!*  Return Codes:  none.                                                        *
!*                                                                              *
!********************************************************************************
 
      SUBROUTINE tran_smooth
 
      include 'COMMONS_INC.f'
 
      dimension tran_o3_std(5001)
      COMMON /init_speccal16/ tran_o3_std
 
      dimension tran_no2_std(5001)
      COMMON /init_speccal17/ tran_no2_std
 
      dimension wavobs(1024),fwhm(1024)
      COMMON /getinput4/ wavobs,fwhm
 
      dimension trncal(1024)
      COMMON /trancal1/ trncal,vaptt
 
      COMMON /getinput5/ nobs,hsurf,dlt,dlt2
 
!C First stage of smoothing - smooth line-by-line high resolution spectrum with
!C     over 300,000 points (point spacing of 0.05 cm-1) to medium resolution
!C     spectrum (resolution of 0.2 nm and point spacing of 0.1 nm) with about
!C     25,000 points.
!C
!C     The smoothing of line-by-line spectrum is done in wavenumber domain. For
!C     a spectrum with a constant 0.2 nm resolution in wavelength domain, it has
!C     variable resolution in wavenumber domain. This effect is properly taken
!C     care of in the design of smoothing functions.
!C
!C     Because the high resolution spectrum is in wavenumber units (cm-1), while
!C     the medium resolution spectrum is in wavelength units, the two kinds of
!C     grids do not automatically match. In order to match the grids, arrays
!C     of INDEX_MED and TRAN_MED_INDEX are specially designed. The desired
!C     medium resolution spectrum, TRAN_MED, at constant 0.1 nm point spacing
!C     and 0.2 nm resolution is obtained through linear interpolation of
!C     TRAN_MED_INDEX array.
!C
      DO 466 j =1, np_med
 
             tran_med_index(j)  = 0.0
             ncvtot_wavno = 2 * ncvhf_wavno(j) - 1
    
             sumins = 0.0
 
          DO 560 i = ncvhf_wavno(j), ncvtot_wavno
             finstr_wavno(i) = &
                exp( -const1*(float(i-ncvhf_wavno(j))*dwavno &
                /fwhm_wavno(j))**2.0)
             sumins = sumins + finstr_wavno(i)
 560      CONTINUE
 
          DO 565 i = 1, ncvhf_wavno(j)-1
             finstr_wavno(i) = finstr_wavno(ncvtot_wavno-i+1)
             sumins = sumins + finstr_wavno(i)
 565      CONTINUE
 
          sumins = sumins * dwavno
 
          DO 570 i = 1, ncvtot_wavno
             finstr_wavno(i) = finstr_wavno(i)*dwavno/sumins
 570      CONTINUE
  
!C*** !!!High resolution transmittances of CO2, N2O, CO, CH4, and O2 should
!C         also be calculated somewhere else (wavelength start = 0.56 micron).
!C***  Here assuming TCO2*TN2O*TCO*TCH4*TO2 is already calculated previously,
!C     i.e., TRANS(I) = TRAN_CO2(I)*TRAN_N2O(I)*TRAN_CO(I)*TRAN_CH4(I)*TRAN_O2(I)
!C      and TRAN_HI(I) = TRAN_HI_H2O(I) * TRANS(I), and TRAN_HI_H2O for varying
!C      water vapor amounts is calculated in this subroutine.
 
          DO 491 k = index_med(j)-(ncvhf_wavno(j)-1), &
                                 index_med(j)+ncvhf_wavno(j)-1
             tran_med_index(j) = tran_med_index(j) + tran_hi(k)*  &
                          finstr_wavno(k-index_med(j)+ncvhf_wavno(j))
 491      CONTINUE
 
 466  CONTINUE
 
!C
!C Linear interpolation to get TRAN_MED from TRAN_MED_INDEX:
!C     (Note that WAVLN_MED_INDEX(J) >= WAVLN_MED(J)    )
!C
         tran_med(1)      = tran_med_index(1)  
         tran_med(np_med) = tran_med_index(np_med)  
 
      DO j = 2, np_med-1
         IF(wavln_med_index(j).LE.wavln_med(j)) THEN
           tran_med(j) = tran_med_index(j)
         ELSE
           dlt  =  wavln_med_index(j) - wavln_med_index(j-1)
           fjm1 = (wavln_med_index(j) - wavln_med(j))        /dlt
           fj   = (wavln_med(j)       - wavln_med_index(j-1))/dlt
           tran_med(j) = fjm1*tran_med_index(j-1) + fj*tran_med_index(j)
!C---
!C---           print*,j,fjm1,fj
!C---
         END IF
      END DO
!C
!C--- Here multiplying O3 and NO2 spectra and other spectrum at medium resolution:
!C
       DO i = 1, np_std
          tran_std(i) = 1.
       END DO
 
       DO i = 1, no3pt
          tran_std(i) = tran_o3_std(i) * tran_no2_std(i)
       END DO
 
       DO i = npshif+1, np_std
          tran_std(i) = tran_std(i)*tran_med(i-npshif)
       END DO
 
 
!C The 2nd stage of smoothing - smooth the medium resolution spectrum (resolution
!C     of 0.2 nm and point spacing of 0.1 nm) with about 25,000 points to match
!C     the coarser and variable resolution spectrum from imaging spectrometers.
!C
!C Initialize some index parameters:
      ia = 1000
!C
      DO 1466 j =1, nobs
 
             trncal(j) = 0.0
             tran_ia     = 0.0
             tran_iap1   = 0.0
 
             ncvtot = 2 * ncvhf(j) - 1
!C---
!C---        print*,'J= ',j, 'NCVHF =', NCVHF(J), 'NCVTOT=',NCVTOT
    
!C Calculate instrumental response functions...
             sumins = 0.0
 
          DO 1560 i = ncvhf(j), ncvtot
             finstr(i) = &
                exp( -const1*(float(i-ncvhf(j))*dwavln &
                /fwhm(j))**2.0)
             sumins = sumins + finstr(i)
1560      CONTINUE
 
          DO 1565 i = 1, ncvhf(j)-1
             finstr(i) = finstr(ncvtot-i+1)
             sumins = sumins + finstr(i)
1565      CONTINUE
 
          sumins = sumins * dwavln
 
          DO 1570 i = 1, ncvtot
             finstr(i) = finstr(i)*dwavln/sumins
1570      CONTINUE
  
!C---          IF(j.eq.87) then
!C---           print*, ' J = 87 =', J
!C---           do i = 1, NCVTOT
!C---            print*,i, FINSTR(i)
!C---           end do
!C---          end if
 
!C  Index searching...
!C
          CALL hunt(wavln_std, np_std, wavobs(j), ia)
!C
!C---          print*,'J =', J, ' IA =', IA
!
!C  Smoothing...
!C
          DO 1491 k = ia-(ncvhf(j)-1), ia+ncvhf(j)-1
             tran_ia = tran_ia + tran_std(k)* &
                          finstr(k-ia+ncvhf(j))
 
!C---             IF(J.eq.1) then
!C---          print*,'j =', j, 'K = ',K, ' IA= ',IA,
!C---     &       K-IA+NCVHF(J),
!C---     &       FINSTR(K-IA+NCVHF(J))
!C---             End IF
 
1491      CONTINUE
 
            ia_p1 = ia + 1
          DO 1492 k = ia_p1-(ncvhf(j)-1), ia_p1+ncvhf(j)-1
             tran_iap1 = tran_iap1 + tran_std(k)* &
                          finstr(k-ia_p1+ncvhf(j))
1492      CONTINUE
!C
!C Linear interpolation to get TRNCAL from TRAN_IA and TRAN_IAP1:
!C
           dlt_ia  =  wavln_std(ia_p1) - wavln_std(ia)
           fia     = (wavln_std(ia_p1) - wavobs(j)) /dlt_ia
!C          FIA_P1  = (WAVOBS(J)     - WAVLN_STD(IA))/DLT_IA
           fia_p1  = 1. - fia
           trncal(j) = fia*tran_ia + fia_p1*tran_iap1
!C---
!C---       print*,'j=',j,'IA =',IA,'FIA =',FIA,'FIA_P1=',FIA_P1,DLT_IA
!C---
!C
1466  CONTINUE
 
 
 
      RETURN
      END
!********************************************************************************
!*                                                 *
!*  Name: CHNLRATIO                                *
!*  Purpose: Calculate 3-channel ratios.                           *
!*  Parameters: none                                   *
!*  Algorithm: The 0.94-um water vapor absorption channel is ratioed against    *
!*             the linear combination of two window channels near 0.86 and      *
!*             1.03 um to obtain one channel ratio for the 0.94-um band.           *
!*             Similar calculation is done for the 1.14-um water vapor band.    *
!*  Globals used: NB1,NB2,NBP94,NB3,NB4,NB1P14 - number of points used in       *
!*                          channel ratios for both the .94- and 1.14-um regions*
!*                IST1,IED1,IST2,IED2,ISTP94,IEDP94 - 3-channel ratioing        *
!*                          parameters for the 0.94-um water vapor band        *
!*                IST3,IED3,IST4,IED4,IST1P14,IED1P14 - 3-channel ratioing      *
!*                          parameters for the 1.14-um water vapor band.           *
!*       WT1,WT2,WT3,WT4,JA - Relative weights for the four window     *
!*                          channels used in channel-ratioing calculations. JA  *
!*             is an output parameter from a table searching       *
!*             routine.                        *
!*       TRNCAL -  Atmospheric transmittance spectra.              *
!*  Global output:R094,R114 - 3-channel ratio values for the 0.94- and 1.14-um  *
!*                          water vapor bands.                     *
!*  Return Codes: none                                 *
!*  Special Considerations: none                               *
!*                                         *
!********************************************************************************
 
      SUBROUTINE chnlratio
 
!C  Common variables
      dimension trncal(1024)
 
      COMMON /getinput7/ nb1,nb2,nbp94,nb3,nb4,nb1p14
      COMMON /init_speccal6/ ist1,ied1,ist2,ied2,istp94,iedp94
      COMMON /init_speccal7/ ist3,ied3,ist4,ied4,ist1p14,ied1p14
      COMMON /init_speccal8/ wt1,wt2,wt3,wt4,ja
      COMMON /trancal1/ trncal,vaptt
      COMMON /chnlratio1/ r094,r114
 
!C Calculate average of spectra over window and water vapor absorption regions.
      const1=0.0
      DO 560 i=ist1,ied1
        const1=const1+trncal(i)
  560 CONTINUE
      const1=const1/float(nb1)
 
      const2=0.0
      DO 570 i=ist2,ied2
        const2=const2+trncal(i)
  570 CONTINUE
      const2=const2/float(nb2)
 
      const3=0.0
      DO 575 i=istp94,iedp94
        const3=const3+trncal(i)
  575 CONTINUE
      const3=const3/float(nbp94)
 
      r094=const3/((wt1*const1) + (wt2*const2))
 
      const4=0.0
      DO 580 i=ist3,ied3 
        const4=const4+trncal(i)
  580 CONTINUE
      const4=const4/float(nb3)
 
      const5=0.0
      DO 590 i=ist4,ied4
        const5=const5+trncal(i)
  590 CONTINUE
      const5=const5/float(nb4)
 
      const6=0.0
      DO 595 i=ist1p14,ied1p14
        const6=const6+trncal(i)
  595 CONTINUE
      const6=const6/float(nb1p14)
 
      r114=const6/((wt3*const4) + (wt4*const5))
 
      RETURN
      END
 
!********************************************************************************
!*                                         *
!*  Name: PROCESS_CUBE                                 *
!*  Purpose: Processes an input cube one spectral slice at a time to derive     *
!*      surface reflectance for each spectrum and calculate the column     *
!*      water vapor amount for each pixel.  The derived surface reflectance*
!*      values are written to an output image file with the same dimensions*
!*      as the input image, and the column water vapor amounts are written *
!*      to a separate file as a single channel image.              *
!*  Parameters: none                                   *
!*  Algorithm: C programs are used to perform all of the file I/O.  A spectral  *
!*        slice is read from the input file.  Each spectrum in the slice   *
!*        is divided by the solar irradiance curve, then the water vapor   *
!*        ratio is calculated.  The corresponding ratio and associated     *
!*        transmission spectrum is used to derive the surface reflectance  *
!*        values.                                  *
!*                                         *
!*        The output is an image file in the same storage order format as  *
!*        the input image file.  It has a 512-byte SIPS header.  The values*
!*        are surface reflectance values multiplied by an input scale      *
!*             factor to make them integer*2 data.                 *
!*                                         *
!*        The output water vapor file is in BSQ format. It has a 512-byte  *
!*        SIPS header.  It contains 1 channel with each pixel value        *
!*        representing the column water vapor amount in units of cm        *
!*        multiplied by a scale factor of 1000 to make the values          *
!*        integer*2.                               *
!*                                             *
!*  Globals used: NSAMPS,NLINES,NBANDS - Input image dimensions used in I/O and *
!*           cube processing.                                  *
!*           FPIN - File pointer for input cube passed to rdslice()        *
!*           FPOCUB - File pointer for output cube passed to wtslice()          *
!*           FPOH2O - File pointer for water vapor image passed to wtslice()    *
!*           SPATIAL,ST_SAMPLE,ST_ROW,SAMPLES_TODO,ROWS_TODO,ISLICE,DIMS - Cube *
!*             dimensions used in rdslice() and wtslice().             *
!*           BUFFER - Holds the slice of data read by rdslice() and written by  *
!*             wtslice().                              *
!*           H2OBUF - Holds the water vapor image.                 *
!*  Global output: None.                                   *
!*  Return Codes: None.                                *
!*  Special Considerations: None.                          *
!*                                         *
!********************************************************************************
      SUBROUTINE process_cube
 
      use cubeio
 
      INTEGER              LUN_IN, LUN_OUT, LUN_VAP, I_RET
      COMMON /inout_units/ lun_in, lun_out, lun_vap
 
!C  Common variables
      dimension yy(1024)       ! Observed radiances
 
!C  parameters for cube I/O
      CHARACTER (LEN = 1000) :: FINAV,FOCUB,FOH2O
      INTEGER*2 BUFFER(8388608)
      INTEGER*2 H2OBUF(8388608)
      INTEGER SORDER,HDREC
 
      INTEGER I_Sample, J_Line
 
!C Local variables
      REAL TBUF(1024)
      INTEGER OFFSET
 
      INTEGER J, ISAMP, IBAND
      real*4  scalef,clmwvp
 
      dimension yirr(1024)
      COMMON /solar_irr1/yirr !RJH
      DATA pi/3.1415926535897/ !RJH
!C Putting these large arrays in global memory forces memory allocation at
!C program initialization.  Otherwise, the program chokes on the first
!C call to PROCESS_CUBE.
!C---      COMMON /DUMMYGLOB/ BUFFER,H2OBUF
 
      COMMON /proccube1/ yy
      COMMON /getinput11/hdrec,nsamps,nlines,nbands,sorder
      COMMON /getinput12/scalef
 
!C Commons for use with the C programs for cube I/O
      COMMON /outcube/ focub
      COMMON /incube/ finav
      COMMON /outh2ovap/ foh2o
      COMMON /rjhdebug/ const1,const2,const3,const4,const5,const6, &
             r094co, r114co, jb
      COMMON /init_speccal8/ wt1,wt2,wt3,wt4,ja
      COMMON /init_speccal6/ ist1,ied1,ist2,ied2,istp94,iedp94
       open(90,file='atrem_yy_pix512line2000_before.txt',form='formatted')
       open(91,file='atrem_yy_pix512line2000_after.txt',form='formatted')
 
!C Read input cube one slice at a time into the array BUFFER.
!C BUFFER contains a spectral slice of data.  The first NBANDS elements are
!C the first spectrum, the second NBANDS elements are the second spectrum, etc
!C Process each spectrum, then write to output file one slice at a time.
      DO 11 j=1,nlines
 
!C Let us know the progress of program.
!        WRITE(*,*)'LINE=',J
 
!C-Note: BUFFER is now passed in through COMMON /DUMMYGLOB/ BUFFER,H2OBUF
!C    for reducing the required total computer memory
!C---        CALL RD_SPECTRAL_SLICE(J,HDREC,NSAMPS,NLINES,NBANDS,
!C---     &                         SORDER)
!        CALL RD_SLICE(LUN_IN,NSAMPS,NBANDS,SORDER,BUFFER)
 
!C For each spectrum in the slice, assign it to the YY array, and use it in
!C REFLDRV to create a surface reflectance spectrum.
        DO 46 isamp= 1,nsamps
          offset = (isamp-1) * nbands
          DO 45 iband=1,nbands           ! fill YY array with this spectrum
            yy(iband) = buffer(offset + iband)
   45     CONTINUE
    
!C Write out sample input spectrum.
        IF((isamp.EQ.512).AND.(j.EQ.2000)) THEN
!C            WRITE(90,*)'Islice=3, Isamp=2, Spectrum before REFLDRV'
           DO 49 i=1,nbands
  49         WRITE(90,*) i,isamp,j,pi*yy(i)/(50*yirr(i))  !HICO CORRECTION
 
         ENDIF
 
          CALL app_refl_with_gas_removal_drv(clmwvp)
!C         IF((ISAMP.EQ.512).AND.(J.EQ.2000)) THEN
!C           WRITE(91,*)'Islice=3, Isamp=2, Spectrum before REFLDRV'
           DO  i=1,nbands
            WRITE(91,'(3i5,9F9.4,4I5)') i,isamp,j,yy(i), &
           r094co,r114co,const1,const2,const3, &
           const4,const5,const6,ja,jb,ist2,ied2
           end do
!C         ENDIF
!C---          CALL REFLDRV(CLMWVP)
 
!C The YY array now contains true surface reflectance values. Multiply
!C by input scale factor to convert values to integer*2 format.
          DO 55 iband=1,nbands
            tbuf(iband)=scalef*yy(iband)
            buffer(offset + iband) = nint(tbuf(iband))
   55     CONTINUE
 
!C Save the water vapor amount (also scaled by 1000 to convert to integer*2).
!C Later, H2OBUF will be written to the water vapor output file.
 
          h2obuf(isamp) = nint(1000.*clmwvp)
 
!C Write out sample processed spectrum.
!C--D         IF((ISAMP.EQ.2).AND.(J.EQ.3)) THEN
!C--D           WRITE(*,*)'Islice=3, Isamp=2, Spectrum after App_Refl_DRV'
!C--D           DO 48 I=1,NBANDS
!C--D             WRITE(*,*)'I=',I,' YY(I)=',YY(I)
!C--D  48         WRITE(*,*)'I=',I,' TBUF(I)=',TBUF(I)
!C--D         ENDIF
 
   46   CONTINUE
 
!C Write slice BUFFER to output file.  The output file is the same storage order
!C as the input file, and it always has a header that is one record (512 bytes).
!C-Note: BUFFER is now passed through COMMON /DUMMYGLOB/ BUFFER,H2OBUF
!C---        CALL WT_SPECTRAL_SLICE(FPOCUB,J,512,NSAMPS,NLINES,NBANDS,
!C---     &                SORDER)
         
!        CALL WT_SLICE(LUN_OUT,NSAMPS,NBANDS,SORDER,BUFFER)
!
!        CALL WT_LINE(LUN_VAP,NSAMPS,H2OBUF)
 
   11 CONTINUE
 
!C Write the water vapor image to a file.  The file is one spatial image
!C with one header record (512 bytes).
!C-Note: H2OBUF is now passed through COMMON /DUMMYGLOB/ BUFFER,H2OBUF
!C---      CALL WT_SPATIAL_SLICE(FPOH2O,1,512,NSAMPS,NLINES,NBANDS,0)
 
!      CALL CLOSEINFILE(LUN_IN)
!      CALL CLOSEOUTFILE(LUN_OUT)
!      CALL CLOSEVAPFILE(LUN_VAP)
 
!        close(90)
!        close(91)
      RETURN
      END
 
!********************************************************************************
!*                                         *
!*  Name: App_Refl_With_Gas_Removal_DRV                            *
!*  Purpose: To derive column water vapor amount and apparent reflectance curve *
!*            with atmospheric gaseous absorption effect removed from an input  *
!*            spectrum using a 3-channel ratioing technique.                    *
!*  Parameters: none                                   *
!*  Algorithm: the algorithm includes following procedures:            *
!*     1. An input radiance spectrum is divided by a solar radiance curve above *
!*            the atmosphere to derive 'APPARENT REFLECTANCE SPECTRUM T(Lambda)'*
!*     2. Two 3-channel ratios:  T(0.94 um)/(WT1*T(0.86)+WT2*T(1.02))          *
!*            and T(1.14 um)/(WT3*T(1.05)+WT4*T(1.23)), are calculated from     *
!*            the observed spectrum T(Lambda).                     *
!*     3. Based on the two ratios using look-up table procedures plus linear    *
!*            interpolation techniques, the column water vapor values and the   *
!*            atmospheric transmittance spectrum are derived.              *
!*     4. The apparent reflectance spectrum is divided by the calculated        *
!*            transmittance spectrum to remove atmospheric gaseous absorption   *
!*            features.                                                         *
!*     5. The output apprent reflectance spectrum with atmospheric gaseous      *
!*            absorption effect is smoothed if desired.                         *
!*  Globals used:                                  *
!*  Global output:                                     *
!*  Return Codes: none                                 *
!*  Special Considerations: None                               *
!*                                         *
!********************************************************************************
 
      SUBROUTINE app_refl_with_gas_removal_drv(CLMWVP)
 
!C  Common variables
      dimension vaptot(60), r0p94(60), r1p14(60), trntbl(1024,60)
      dimension yirr(1024)
      dimension yy(1024) 
      dimension rr(1024) 
      dimension rotot(1050), ttot(1050), stot(1050)
      dimension finst2(100)
      dimension ssh2o(60)
      dimension trncal(1024)
 
!C Local variables
      dimension trncv(1024)
 
      COMMON /getinput5/ nobs,hsurf,dlt,dlt2
      COMMON /getinput7/ nb1,nb2,nbp94,nb3,nb4,nb1p14
      COMMON /getinput8/ imn,idy,iyr,ih,im,is
      COMMON /geometry2/ gco2,go3,gn2o,gco,gch4,go2,ssh2o,totlo3,ggeom
      COMMON /init_speccal3/ nh2o
      COMMON /init_speccal6/ ist1,ied1,ist2,ied2,istp94,iedp94
      COMMON /init_speccal7/ ist3,ied3,ist4,ied4,ist1p14,ied1p14
      COMMON /init_speccal8/ wt1,wt2,wt3,wt4,ja
      COMMON /init_speccal10/ ncv2,ncvhf2,ncvtt2,istrt2,iend2,finst2
      COMMON /init_speccal11/ natot,nbtot,nctot,ndtot
      COMMON /tran_table1/ sh2o,vaptot,r0p94,r1p14,trntbl
      COMMON /trancal1/ trncal,vaptt
      COMMON /solar_irr1/yirr
      COMMON /proccube1/ yy
      COMMON /sixs1/ rotot, ttot, stot
      COMMON /rjhdebug/ const1,const2,const3,const4,const5,const6, &
             r094co, r114co, jb
      dimension g_vap(25), g_other(25)
      COMMON /geometry3/ g_vap, g_other, g_vap_equiv, vap_slant_mdl
 
!C Parameters for plane observations:
      CHARACTER (LEN = 80) :: NAME_INSTRU, NAMES(10)
      COMMON /getinput13/ name_instru, names
 
!C Arrays related to look-up table:
!C     VAPTOT: TOTAL SUN-SURFACE-SENSOR PATH WATER VAPOR IN UNITS OF CM
!C     R0P94 : Tc(0.94 um)/(WT1*Tc(0.86)+WT2*Tc(1.02))
!C     R1P14 : Tc(1.14 um)/(WT3*Tc(1.05)+WT4*Tc(1.23))
 
      dimension speca(1024),specb(1024),specav(1024)
 
    integer*4 itmp1
    real*4 xfa
    real*4 xfb
    real*4 xfc
    real*4 xfd
 
      DATA pi,dtorad /3.1415926,0.0174533/
 
!C
!C Loops for processing image data cube ...
      IF((name_instru.EQ.names(1)).AND.(iyr.LE.2009)) THEN
 
!C  Ratioing observed spectrum against solar irradiance curve.
!C  The radiance data from each year is scaled by a different constant.
!C    XFA = factor for A spectrometer
!C    XFB = factor for B spectrometer
!C    XFC = factor for C spectrometer
!C    XFD = factor for D spectrometer
!C
!C rewritten to handle each spectrometer differently
!C                                         9/14/1995 Roger N. Clark, U.S.G.S.
 
      IF(iyr.LE.1989) THEN
    xfa = 10.0
    xfb = 10.0
    xfc = 10.0
    xfd = 10.0
      ENDIF
 
      IF((iyr.GE.1990).AND.(iyr.LE.1991)) THEN
    xfa = 20.0
    xfb = 20.0
    xfc = 20.0
    xfd = 20.0
      ENDIF
 
      IF((iyr.GE.1992).AND.(iyr.LE.1994)) THEN
    xfa = 50.0
    xfb = 50.0
    xfc = 50.0
    xfd = 50.0
      ENDIF
 
      IF((iyr.GE.1995).AND.(iyr.LE.2009)) THEN
    xfa = 50.0
    xfb = 50.0
    xfc = 50.0
    xfd = 100.0
      ENDIF
 
!C   AVIRIS A Spectrometer (ch = 1-32)
 
      IF(nobs.GE.32) THEN
    itmp1 = 32
      ELSE
    itmp1 = nobs
      ENDIF
 
      DO 551 i=1,itmp1
          yy(i)=pi*yy(i)/(xfa*yirr(i))
  551 CONTINUE
 
!C   AVIRIS B Spectrometer (ch = 33-96)
 
      IF(nobs.GE.96) THEN
    itmp1 = 96
      ELSE if ((nobs.GE.33).and.(nobs.LT.96)) THEN
    itmp1 = nobs
      ELSE
    go to 559
      ENDIF
 
      DO 552 i=33,itmp1
          yy(i)=pi*yy(i)/(xfb*yirr(i))
  552 CONTINUE
 
!C   AVIRIS C Spectrometer (ch = 97-160)
 
      IF(nobs.GE.160) THEN
    itmp1 = 160
      ELSE if ((nobs.GE.97).and.(nobs.LT.160)) THEN
    itmp1 = nobs
      ELSE
    go to 559
      ENDIF
 
      DO 553 i=97,itmp1
          yy(i)=pi*yy(i)/(xfc*yirr(i))
  553 CONTINUE
 
!C   AVIRIS D Spectrometer (ch = 161-224)
 
      IF(nobs.GE.161) THEN
    itmp1 = nobs
      ELSE
    go to 559
      ENDIF
 
      DO 554 i=161,itmp1
          yy(i)=pi*yy(i)/(xfd*yirr(i))
  554 CONTINUE
 
!C branch point in case not all spectrometers are being done.
  559   CONTINUE
 
!C End of R. Clark contribution
      END IF
!C
!C Loops for processing AVIRIS data starting from year 2010---
      IF((name_instru.EQ.names(1)).AND.(iyr.GE.2010)) THEN
 
          f_av_1 = 30.
          f_av_2 = 60.
          f_av_3 = 120.
 
        DO i = 1, 110
           yy(i) = pi*yy(i) / (f_av_1 * yirr(i))
        END DO
 
        DO i = 111, 160
           yy(i) = pi*yy(i) / (f_av_2 * yirr(i))
        END DO
 
        DO i = 161, 224
           yy(i) = pi*yy(i) / (f_av_3 * yirr(i))
        END DO
 
      END IF
 
!C Loops for processing HYDICE data ...
!C
      IF(name_instru.EQ.names(2)) THEN
 
          f_hydice = 75.
 
        DO i = 1, nobs
           yy(i) = pi*yy(i) / (f_hydice * yirr(i))
        END DO
 
      END IF 
 
 
!C Loops for processing HSI data ...
!C
      IF(name_instru.EQ.names(3)) THEN
 
          f_hsi = 100.
 
        DO i = 1, nobs
           yy(i) = pi*yy(i) / (f_hsi * yirr(i))
        END DO
 
      END IF 
 
 
!C Loops for processing TRWIS-III data ...
!C
      IF(name_instru.EQ.names(4)) THEN
 
          f_trwis = 100.
 
        DO i = 1, nobs
           yy(i) = pi*yy(i) / (f_trwis * yirr(i))
        END DO
 
      END IF 
 
!C Loops for processing EO-1 Hyperion data ...
!C
      IF(name_instru.EQ.names(6)) THEN
 
          f_hyperion_a = 40.
          f_hyperion_b = 80.
 
        DO i = 1, 70
           yy(i) = pi*yy(i) / (f_hyperion_a * yirr(i))
        END DO
 
        DO i = 71, nobs 
           yy(i) = pi*yy(i) / (f_hyperion_b * yirr(i))
        END DO
 
      END IF
 
!C Loops for processing HICO data ...
!C
      IF(name_instru.EQ.names(7)) THEN
 
!C---          F_HICO = 50./1.32
          f_hico = 50.
 
        DO i = 1, nobs
           yy(i) = pi*yy(i) / (f_hico * yirr(i))
        END DO
 
      END IF
 
!C Loops for processing NIS (NEON Imaging Spectrometer) data ...
!C
      IF(name_instru.EQ.names(8)) THEN
 
          f_nis = 100.
 
        DO i = 1, nobs
           yy(i) = pi*yy(i) / (f_nis * yirr(i))
        END DO
 
      END IF 
 
!C Loops for processing PRISM Imaging Spectrometer data ...
!C
      IF(name_instru.EQ.names(9)) THEN
 
          f_prism = 100.
 
        DO i = 1, nobs
           yy(i) = pi*yy(i) / (f_prism * yirr(i))
        END DO
 
      END IF
 
!C
!C Calculating 3-channel ratios from an observed spectrum, using a
!C  look-up table procedure to derive the column amount of water vapor
!C  and to find the associated simulated transmittance spectrum.
      const1=0.0
      DO 560 i=ist1,ied1
        const1=const1+yy(i)
  560 CONTINUE
      const1=const1/float(nb1)
 
      const2=0.0
      DO 570 i=ist2,ied2
        const2=const2+yy(i)
  570 CONTINUE
      const2=const2/float(nb2)
 
      const3=0.0
      DO 575 i=istp94,iedp94
        const3=const3+yy(i)
  575 CONTINUE
      const3=const3/float(nbp94)
 
      r094co=const3/((wt1*const1) + (wt2*const2))
      r094c =r094co
 
      IF(r094co.GT.1.0) THEN
        const1=0.0
        DO 561 i=ist1,ied1
          const1=const1+1.0/yy(i)
  561   CONTINUE
        const1=const1/float(nb1)
 
        const2=0.0
        DO 571 i=ist2,ied2
          const2=const2+1.0/yy(i)
  571   CONTINUE
        const2=const2/float(nb2)
 
        const3=0.0
        DO 576 i=istp94,iedp94
          const3=const3+1.0/yy(i)
  576   CONTINUE
        const3=const3/float(nbp94)
 
        r094c=const3/((wt1*const1) + (wt2*const2))
      ENDIF
 
      const4=0.0
      DO 580 i=ist3,ied3
        const4=const4+yy(i)
  580 CONTINUE
      const4=const4/float(nb3)
 
      const5=0.0
      DO 590 i=ist4,ied4
        const5=const5+yy(i)
  590 CONTINUE
      const5=const5/float(nb4)
 
      const6=0.0
      DO 595 i=ist1p14,ied1p14
        const6=const6+yy(i)
  595 CONTINUE
      const6=const6/float(nb1p14)
 
      r114co=const6/((wt3*const4) + (wt4*const5))
      r114c =r114co
 
      IF(r114co.GT.1.0) THEN
        const4=0.0
        DO 581 i=ist3,ied3
          const4=const4+1.0/yy(i)
  581   CONTINUE
        const4=const4/float(nb3)
 
        const5=0.0
        DO 591 i=ist4,ied4
          const5=const5+1.0/yy(i)
  591   CONTINUE
        const5=const5/float(nb4)
 
        const6=0.0
        DO 596 i=ist1p14,ied1p14
          const6=const6+1.0/yy(i)
  596   CONTINUE
        const6=const6/float(nb1p14)
 
        r114c=const6/((wt3*const4) + (wt4*const5))
      ENDIF
 
      CALL hunt(r0p94,nh2o,r094c,ja)
      IF (ja.GT.0.AND.ja.LT.60) THEN
        dlta  =r0p94(ja+1)-r0p94(ja)
        fja   =(r0p94(ja+1)-r094c)/dlta
        fjap1 =(r094c-r0p94(ja))/dlta
 
        vaptta=fja*vaptot(ja)+fjap1*vaptot(ja+1)
        IF(r094co.GT.1.) vaptta=-vaptta
      ELSE
        IF(ja.LE.0) vaptta = vaptot(ja+1)
        IF(ja.GE.60) vaptta = vaptot(ja)
      ENDIF
 
      IF(r094co.LE.1.)THEN
        DO 900 i=1,nobs
          IF (ja.GT.0.AND.ja.LT.60) THEN
            speca(i)=fja*trntbl(i,ja)+fjap1*trntbl(i,ja+1)
          ELSE
            IF(ja.LE.0) speca(i)=trntbl(i,ja+1)
            IF(ja.GE.60) speca(i)=trntbl(i,ja)
          ENDIF
  900   CONTINUE
      ENDIF
 
      IF(r094co.GT.1.)THEN
        DO 901 i=1,nobs
          IF (ja.GT.0.AND.ja.LT.60) THEN
            speca(i)=1.0/(fja*trntbl(i,ja)+fjap1*trntbl(i,ja+1))
          ELSE
            IF(ja.LE.0) speca(i)=1.0/trntbl(i,ja+1)
            IF(ja.GE.60) speca(i)=1.0/trntbl(i,ja)
          ENDIF
  901   CONTINUE
      ENDIF
 
      jb = ja
      CALL hunt(r1p14,nh2o,r114c,jb)
      IF (jb.GT.0.AND.jb.LT.60) THEN
        dltb  =r1p14(jb+1)-r1p14(jb)
        fjb   =(r1p14(jb+1)-r114c)/dltb
        fjbp1 =(r114c-r1p14(jb))/dltb
  
        vapttb=fjb*vaptot(jb)+fjbp1*vaptot(jb+1)
        IF(r114co.GT.1.) vapttb=-vapttb
      ELSE
        IF(jb.LE.0) vapttb = vaptot(jb+1)
        IF(jb.GE.60) vapttb = vaptot(jb)
      ENDIF
 
      IF(r114co.LE.1.)THEN
        DO 910 i=1,nobs
          IF (jb.GT.0.AND.jb.LT.60) THEN
            specb(i)=fjb*trntbl(i,jb)+fjbp1*trntbl(i,jb+1)
          ELSE
            IF(jb.LE.0) specb(i)=trntbl(i,jb+1)
            IF(jb.GE.60) specb(i)=trntbl(i,jb)
          ENDIF
  910   CONTINUE
      ENDIF
 
      IF(r114co.GT.1.)THEN
        DO 911 i=1,nobs
          IF (jb.GT.0.AND.jb.LT.60) THEN
            specb(i)=1.0/(fjb*trntbl(i,jb)+fjbp1*trntbl(i,jb+1))
          ELSE
            IF(jb.LE.0) specb(i)=1.0/trntbl(i,jb+1)
            IF(jb.GE.60) specb(i)=1.0/trntbl(i,jb)
          ENDIF
  911   CONTINUE
      ENDIF
 
!C---      CLMWVP=0.5*(VAPTTA+VAPTTB)/GGEOM
      clmwvp=0.5*(vaptta+vapttb)/g_vap_equiv
!C
!C  Derivation of surface reflectances
      DO 915 i=1,nobs
        specav(i)=0.5*(speca(i)+specb(i))
!C---        TRNTMP=(YY(I)/SPECAV(I))-ROTOT(I)
!C---        YY(I)=TRNTMP/(TTOT(I)+(STOT(I)*TRNTMP))
        yy(i) = yy(i) / specav(i)
  915 CONTINUE
!C
!C--       print*, 'i, rotot(i), ttot(i), stot(i)'
!C--      DO I = 1, NOBS
!C--       print*, i, rotot(i), ttot(i), stot(i)
!C--      END DO
 
!C
!C Smooth the derived surface reflectance spectra
 
      IF(dlt2.GT.dlt) THEN
!C
!C       First, replace radiances <= 0 near the spectral overlapping parts of the
!C       four AVIRIS spectrometers by radiances in the nearby AVIRIS' channels.
        DO 920 i=natot-2,natot+2
          IF(yy(i).LE.0.0) yy(i)=yy(i-1)
  920   CONTINUE
 
        DO 921 i=nbtot-2,nbtot+2
          IF(yy(i).LE.0.0) yy(i)=yy(i-1)
  921   CONTINUE
 
        DO 922 i=nctot-4,nctot+4
          IF(yy(i).LE.0.0) yy(i)=yy(i-1)
  922   CONTINUE
 
        DO 923 i=ndtot-7,nobs
          IF(yy(i).LE.0.0) yy(i)=yy(i-1)
  923   CONTINUE
 
        DO 3480 i=istrt2,iend2
          trncv(i)=0.0
          ii=i-ncv2
          DO 3490 j=i-ncv2,i+ncv2
            trncv(i)=trncv(i)+yy(j)*finst2(j-ii+1)
 3490     CONTINUE
 3480   CONTINUE
        DO 3495 i=istrt2,iend2
          yy(i)=trncv(i)
 3495   CONTINUE
      ENDIF
 
      RETURN
      END
 
!********************************************************************************
!*                                                 *
!*  Name: LOCATE                                       *
!*  Purpose: given an array XX of length N, and given a value X, returns a value*
!*           J such that X is between XX(J) and XX(J+1).  XX must be monotonic, *
!*  Parameters: XX - monotonic array of values                     *
!*              N  - number of elements in XX                      *
!*              X  - value that will be matched to the XX array                *
!*              J  - index into the XX array where XX(J) <= X <= XX(J+1)           *
!*  Algorithm:  bisectional table searching, copied from Numerical Recipes.     *
!*  Globals used: none                                 *
!*  Global output: none                                *
!*  Return Codes: J=0 or J=N is returned to indicate that X is out of range     *
!*  Special Considerations: none                               *
!*                                         *
!********************************************************************************
      SUBROUTINE locate(xx,n,x,j)
      INTEGER j,n
      REAL x,xx(n)
      INTEGER jl,jm,ju
      jl=0
      ju=n+1
10    if(ju-jl.gt.1)then
        jm=(ju+jl)/2
        if((xx(n).ge.xx(1)).eqv.(x.ge.xx(jm)))then
          jl=jm
        else
          ju=jm
        endif
      goto 10
      endif
      if(x.eq.xx(1))then
        j=1
      else if(x.eq.xx(n))then
        j=n-1
      else
        j=jl
      endif
      return
      END
 
!********************************************************************************
!*                                                 *
!*  Name: CUBSPLN                                  *
!*  Purpose: an interface for performing cubic spline interpolations.          *
!*  Parameters: N - number of elements in XORGN, YORGN                 *
!*              XORGN - original x values                      *
!*              YORGN - original y values                      *
!*              XINT  - interpolated x values                      *
!*              YINT  - interpolated y values                      *
!*  Algorithm: Straight forward calculations                       *
!*  Globals used: NOBS - number of spectral points.                *
!*  Global output: none                                    *
!*  Return Codes: none                                 *
!*  Special Considerations: none                               *
!*                                         *
!********************************************************************************
 
      SUBROUTINE cubspln(N,XORGN,YORGN,XINT,YINT)
 
      dimension xorgn(1050),yorgn(1050),y2(1050)
      dimension xint(1024),yint(1024)
      INTEGER N                               !number of elements in XORGN,YORGN
      
      COMMON /getinput5/ nobs,hsurf,dlt,dlt2
!C
!C  YP1 and YP2 are two parameters specifying the values of 2nd derivatives at
!C     XORGN(1) and XORGN(N) and were used in cubic spline interpolations. The
!C     setting of both YP1 and YP2 to 2.0E30 is equivalent to set the 2nd
!C     derivatives at the two boundaries as zero, according to Numeric Recipes.
      yp1=2.0e30
      ypn=2.0e30
 
      CALL spline(xorgn,yorgn,n,yp1,ypn,y2)
 
      DO 320 i=1,nobs
        x=xint(i)
        CALL splint(xorgn,yorgn,y2,n,x,y)
        IF(y.LT.0.0) y=0.0
        yint(i)=y
 
  320 CONTINUE
 
      RETURN
      END
!
!********************************************************************************
!*                                             *
!*  Name: SPLINE                                       *
!*  Purpose: program for cubic spline interpolation --- copyed from Numerical   *
!*           Recipes, p.88-89                              *
!*  Parameters: X - x-values                               *
!*              Y - y-values                               *
!*              N - length of X                            *
!*              YP1 - 2nd derivative at X(1)                       *
!*              YPN - 2nd derivative at X(N)                       *
!*              Y2  - 2nd derivatives at all X positions                   *
!*  Algorithm: Given arrays X and Y of length N containing a tabulated function,*
!*      i.e.,Yj=f(Xj), with X1 < X2...<XN, and given values YP1 and YPN for     *
!*      the first derivative of the interpolating function at points 1         *
!*      and N, respectively, this routine returns an array Y2 of length N       *
!*      which contains the second derivatives of the interpolating function     *
!*      at the tabulated points Xj. If YP1 and/or YPN are equal to 1.0E30 or    *
!*      larger, the routine is signalled to set the corresponding boundary      *
!*      condition for a natural spline, with zero second derivative on          *
!*      that boundary.                                 *
!*  Globals used: none                                 *
!*  Global output: none                                *
!*  Return Codes: none                                 *
!*  Special Considerations: none                               *
!*                                         *
!********************************************************************************
 
      subroutine spline(x,y,n,yp1,ypn,y2)
      parameter(nmax=1050)
      integer n,i,k
      real x(n),y(n),y2(n),u(nmax)
      real yp1,ypn,sig,p,qn,un
      if (yp1.gt..99e30) then
        y2(1)=0.
        u(1)=0.
      else
        y2(1)=-0.5
        u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
      endif
      do 11 i=2,n-1
        sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
        p=sig*y2(i-1)+2.
        y2(i)=(sig-1.)/p
        u(i)=(6.*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1)) &
           /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
11    continue
      if (ypn.gt..99e30) then
        qn=0.
        un=0.
      else
        qn=0.5
        un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
      endif
      y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
      do 12 k=n-1,1,-1
        y2(k)=y2(k)*y2(k+1)+u(k)
12    continue
      return
      end
 
!********************************************************************************
!*                                                 *
!*  Name: SPLINT                                       *
!*  Purpose: calculate a cubic-spline interpolated Y value, -- copied from      *
!*           Numerical Recipes.                            *
!*  Parameters: XA  - original X array                         *
!*              YA  - original Y array                         *
!*              Y2A - 2nd derivative array from SUBROUTINE SPLINE          *
!*              N   - number of X elements                     *
!*              X   - an X position at which interpolation should be made       *
!*              Y   - interpolated y value at position X               *
!*  Algorithm: Given the arrays XA and YA of length N, which tabulate a function*
!*      (with the XAj's in order), and given the array Y2A, which is the output *
!*      from SPLINE above, and given a value of X, this routine returns a       *
!*      cubic-spline interpolated value Y.                     *
!*  Globals used: none                                 *
!*  Global output: none                                *
!*  Return Codes: none                                 *
!*  Special Considerations: SPLINE is called only once to process an entire     *
!*      tabulated function in arrays Xi and Yi. Once this has been done, values *
!*      of the interpolated function for any value of X are obtained by calls   *
!*      (as many as desired) to a separate routine SPLINT (for cubic spline     *
!*      interpolation")                                *
!*                                         *
!********************************************************************************
 
      subroutine splint(xa,ya,y2a,n,x,y)
      integer n,klo,khi,k
      real xa(n),ya(n),y2a(n)
      real x,y,h,a,b
      klo=1
      khi=n
1     if (khi-klo.gt.1) then
        k=(khi+klo)/2
        if(xa(k).gt.x)then
          khi=k
        else
          klo=k
        endif
      goto 1
      endif
      h=xa(khi)-xa(klo)
      if (h.eq.0.) stop 'bad xa input.'
      a=(xa(khi)-x)/h
      b=(x-xa(klo))/h
      y=a*ya(klo)+b*ya(khi)+ &
           ((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**2)/6.
      return
      end
 
!********************************************************************************
!*                                                 *
!*  Name: HUNT                                     *
!*  Purpose: finds the element in array XX that is closest to value X.  Array AA*
!*      must be monotonic, either increasing or decreasing.            *
!*  Parameters: XX  - array to search                          *
!*              N - number of elements in the array                *
!*              X - element to search for closest match                *
!*         JLO - index of the closest matching element             *
!*  Algorithm: this subroutine was copied from Numerical Recipes               *
!*  Globals used: none                                 *
!*  Global output: none                                *
!*  Return Codes: none                                 *
!*  Special Considerations: none                               *
!*                                         *
!********************************************************************************
 
      SUBROUTINE hunt(xx,n,x,jlo)
      INTEGER jlo,n
      REAL x,xx(n)
      INTEGER inc,jhi,jm
      LOGICAL ascnd
      ascnd=xx(n).ge.xx(1)
      if(jlo.le.0.or.jlo.gt.n)then
        jlo=0
        jhi=n+1
        goto 3
      endif
      inc=1
      if(x.ge.xx(jlo).eqv.ascnd)then
1       jhi=jlo+inc
        if(jhi.gt.n)then
          jhi=n+1
        else if(x.ge.xx(jhi).eqv.ascnd)then
          jlo=jhi
          inc=inc+inc
          goto 1
        endif
      else
        jhi=jlo
2       jlo=jhi-inc
        if(jlo.lt.1)then
          jlo=0
        else if(x.lt.xx(jlo).eqv.ascnd)then
          jhi=jlo
          inc=inc+inc
          goto 2
        endif
      endif
3     if(jhi-jlo.eq.1)then
        if(x.eq.xx(n))jlo=n-1
        if(x.eq.xx(1))jlo=1
        return
      endif
      jm=(jhi+jlo)/2
      if(x.ge.xx(jm).eqv.ascnd)then
        jlo=jm
      else
        jhi=jm
      endif
      goto 3
      END
 
!********************************************************************************
!*                                                 *
!*  Name: SUNCOR                                           *
!*  Purpose: computes at some reference time TZ the declination of the sun and  *
!*           the hour angle-TZ.                            *
!*  Algorithm: At any time T, the local solar hour angle is HAZ+T(GMT)-XLONG,   *
!*           where XLONG is the longitude measured positive west of Greenwich.  *
!*           Both time and angles are in radians. To compute azimuth and        *
!*           elevation at latitude XLAT, longitude XLONG, and time T, call      *
!*           HAZEL(HAZ+T-XLONG,DEC,AZ,EL,XLAT). This routine was copied        *
!*           from W. Mankin at NCAR, Boulder, CO.                  *
!*  Globals used: none                                 *
!*  Global output: none                                    *
!*  Return Codes: none                                 *
!*  Special Considerations: none                               *
!*                                         *
!********************************************************************************
 
      SUBROUTINE suncor(IDAY,MONTH,IYR,TZ,DEC,HAZ)
 
      jd=julian(iday,month,iyr)
      fjd=0.5+tz/6.283185307
      CALL solcor(jd,fjd,ras,dec,gsdt,bzero,pzero,solong)
      haz=gsdt-ras-tz
 
      RETURN
      END
 
!********************************************************************************
!*  Name:  SOLCOR                                  *
!*  Purpose: A routine for solar geometry calculations --- copied from         *
!*                  W. Mankin at NCAR.                         *
!*  Parameters: none                                   *
!*  Algorithm:  This routine was supplied by W. Mankin at NCAR, Boulder, CO     *
!*  Globals used:                                  *
!*  Global output:                                     *
!*  Return Codes: none                                 *
!*  Special Considerations: none                               *
!*                                         *
!********************************************************************************
      SUBROUTINE solcor(JD,FJD,RAS,DECS,GSDT,BZRO,P,SOLONG)
 
      pi=3.141592654
      d=(jd-2415020)+fjd
      iyr=d/365.25
      g=-.026601523+.01720196977*d-1.95e-15*d*d -2.*pi*iyr
      xlms=4.881627938+.017202791266*d+3.95e-15*d*d-2.*pi*iyr
      obl=.409319747-6.2179e-9*d
      ecc=.01675104-1.1444e-9*d
      f=d-365.25*iyr
      gsdt=1.739935476+2.*pi*f/365.25+1.342027e-4*d/365.25
      gsdt=gsdt+6.2831853*(fjd-0.5)
      xlts=xlms+2.*ecc*sin(g)+1.25*ecc*ecc*sin(2.*g)
      sndc=sin(xlts)*sin(obl)
      decs=asin(sndc)
      csra=cos(xlts)/cos(decs)
      ras=acos(csra)
      IF(sin(xlts).LT.0.) ras=2.*pi-ras
      omega=1.297906+6.66992e-7*d
      thetac=xlts-omega
      bzro=asin(.126199*sin(thetac))
      p=-atan(cos(xlts)*tan(obl))-atan(.127216*cos(thetac))
      xlmm=atan2(.992005*sin(thetac),cos(thetac))
      jdr=jd-2398220
      irot=(jdr+fjd)/25.38
      frot=(jdr+fjd)/25.38-irot
      solong=xlmm-2.*pi*frot+pi-3.e-4
      IF(solong.LT.0.) solong=solong+2.*pi
      IF(solong.GT.(2.*pi)) solong=solong-2.*pi
 
      RETURN
      END
 
 
!********************************************************************************
!*                                             *
!*  Name: JULIAN                                       *
!*  Purpose: computes the julian day                           *
!*  Parameters: IDAY - day of the year                         *
!*              MONTH - day of the year                        *
!*  Algorithm: This routine was supplied by William Mankin at NCAR, Boulder, CO *
!*  Globals used: none                                 *
!*  Global output:  none                                   *
!*  Return Codes: none                                 *
!*  Special Considerations: none                               *
!*                                         *
!********************************************************************************
 
      FUNCTION julian(IDAY,MONTH,IYR)
 
!C Local variables
      dimension md(12)
 
      DATA md/0,31,59,90,120,151,181,212,243,273,304,334/
 
      IF(iyr.LT.100) iyr=iyr+1900
      jyr=iyr-1600
      i1=jyr/400
      i2=(jyr-400*i1)/100
      i3=(jyr-400*i1-100*i2)/4
      julian=2305447+365*jyr+97*i1+24*i2+i3
      leap=0
      IF(mod(jyr,4).EQ.0) leap=1
      IF(mod(jyr,100).EQ.0) leap=0
      IF(mod(jyr,400).EQ.0) leap=1
 
!C     LEAP=1 if iyr is a leap year
      julian=julian+md(month)+iday
      IF(month.LE.2) julian=julian-leap
 
      RETURN
      END
 
!********************************************************************************
!*                                                 *
!*  Name: HAZEL                                        *
!*  Purpose: Calculates azimuth and elevation                      *
!*  Parameters: none                                   *
!*  Algorithm: This routine was supplied by William Mankin at NCAR, Boulder, CO *
!*  Globals used: None.                                    *
!*  Global output:   None.                             *
!*  Return Codes: None                                 *
!*  Special Considerations: None.                          *
!*                                         *
!********************************************************************************
      SUBROUTINE hazel (H,D,A,E,XLAT)
 
!C     H = HOUR ANGLE
!C     D = DECLINATION
!C     A = AZIMUTH
!C     E = ELEVATION
!C     XLAT = LATITUDE
!C     ALL ANGLES IN RADIANS
 
      pi = 3.14159265
      sne = sin(d)*sin(xlat)+cos(d)*cos(xlat)*cos(h)
      e=asin(sne)
      sna = cos(d)*sin(h)
      csa=(sin(xlat)*cos(h)*cos(d)-sin(d)*cos(xlat))
      a=atan2(sna,csa)+pi
      END
 
!********************************************************************************
!*                                                 *
!*  Name: FINDMATCH                                *
!*  Purpose: finds the closest match for ELEM in LIST                  *
!*  Parameters:  LIST - array of values to match ELEM to.  Elements is array    *
!*            should increase in value.                    *
!*  Algorithm: linearly compare ELEM to each element.  When ELEM is smaller     *
!*             than the LIST(I), then it is assumed to be closest to LIST(I-1)  *
!*             or LIST(I).  The one that has the smallest absolute difference   *
!*             to ELEM is returned as the closest match.                   *
!*  Globals used: none                                 *
!*  Global output: none                                *
!*  Return Codes: the closest matching element index                   *
!*  Special Considerations: none                               *
!*                                         *
!********************************************************************************
 
      FUNCTION findmatch(LIST,NOBS,ELEM)
 
      dimension list(1024)
      INTEGER   nobs
      REAL      elem,list
 
      DO 460 i=1,nobs
        IF(list(i).GT.elem) GOTO 470
  460 CONTINUE
  470 CONTINUE
      diff1=abs(list(i-1)-elem)
      diff2=abs(list(i)-elem)
      IF (diff1.LT.diff2) THEN
        findmatch=i-1
      ELSE
        findmatch=i
      ENDIF
 
      RETURN
      END
!C--------1---------2---------3---------4---------5---------6---------7--