atrem_app_refl_plus_gas_removal_l2.f90

Go to the documentation of this file.

!------------------------------------------------------------------------------|
! 2015 - Modifications for l2gen including speed enhancements by R. Healy      |
!        (richard.healy@nasa.gov)                                              |
! 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                               *
!*                                         *
!********************************************************************************
 
!INTERFACE
!  SUBROUTINE ecdf(xcdf, ycdf, nbin, xs, nsamp) BIND(C)
!    USE, INTRINSIC :: ISO_C_BINDING, ONLY: C_INT, C_FLOAT
!    IMPLICIT NONE
!    FLOAT(C_FLOAT) :: xcdf,ycdf,xs
!    INTEGER(C_INT) :: nsamp,nbin
!  END SUBROUTINE addnums
!END INTERFACE
      
      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(nobs_max),fwhm(nobs_max)
      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,ifullcalc,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
 
 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 ------------
 
  85  FORMAT(a1000)
 
!C        ELSE
 
!C Initialize default atmospheric window and water vapor absorption regions
!C     for 3-channel ratioing calculations.
!C          WNDOW1 = 0.865
!C          NB1    = 3
!C          WNDOW2 = 1.030
!C          NB2    = 3
!C          WP94C  = 0.940
!C          NBP94  = 5
!C          WNDOW3 = 1.050
!C          NB3    = 3
!C          WNDOW4 = 1.235
!C          NB4    = 3
!C          W1P14C = 1.1375
!C          NB1P14 = 7
!C        ENDIF
!C      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
      model = 1
!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)))
          !Convert the atmospheric pressure from "mb" to "atm.":
          p(i) = tpvmr(model,3+(4*(i-1))) / 1013.
          t(i) = tpvmr(model,4+(4*(i-1)))
          !Convert the VMR from the ppm unit in the model to absolute unit
          vmr(i) = tpvmr(model,5+(4*(i-1)))*1.0e-06
  120   CONTINUE
      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
 
!C      READ(*,*)IH2OVP, ICO2, IO3, IN2O, ICO, ICH4, IO2, INO2
 
    ih2ovp = 1
    ico2   = 0
    io3    = 0
    in2o   = 0
    ico    = 0
    ich4   = 0
    io2    = 0
    ino2   = 0
 
      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,ifullcalc,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.":
!      Moved this conversion to GET_INPUT and init_tpvmr - rjh 1/8/2016
!      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.
!      Moved this conversion to GET_INPUT and init_tpvmr - rjh 1/8/2016
!      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(91,*)'Column vapor amount in model atmosphere from ground'
      WRITE(91,*)'       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(91,*)'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(nh2o_max), vap_slant(nh2o_max)
      dimension md(12)
      dimension ssh2o(nh2o_max)
      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
      COMMON /geometry4/vap_slant_mdl
      REAL   MU,MU0, SSH2O_S(NH2O_MAX,2)
      COMMON /geometry5/mu,mu0,ssh2o_s
      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, nh2o_max
        vap_slant(i) = vapvrt(i) * 2.0
      END DO
 
      obszni = senzn_l2
      obsphi = senaz_l2
      solzni = solzn_l2
 
      obszni = obszni / radeg
      obsphi = obsphi / radeg
 
      solzni = solzni/radeg
 
!      write(91,*) 'KPLANE=',k_plane
 
      mu0 = 1./cos(solzni)
      mu  = 1./cos(obszni)
      ggeom  = mu0 + mu
      write(91,*) 'GGEOM =',ggeom,' OBSZNI = ',obszni, ' OBSPHI = ',obsphi, &
       'solzni=',solzni,' degrees :: MU0, MU = ',mu0, mu
      write(91,*) 'f1a,f2a,f1b,f2b=',f1a,f2a,f1b,f2b
      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, vrto3
 
!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)
      vap_sol = clmvapp*mu0
      vap_sen = clmvapp*mu
!      write(91,*) 'VAP_SLANT_MDL =', VAP_SLANT_MDL, ' cm ', SLANT_MDL
!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, vap_slant_mdl, clmvap
       write(91,*) 'VAP_SOL,VAP_SEN = ', vap_sol, vap_sen, wtrvpr
 
      DO 310 i=1,nh2o_max
        ssh2o(i) = vap_slant(i) / vap_slant_mdl
        write(91,*) 'SSH2O(I), I = ', ssh2o(i), i, vap_slant_mdl, wtrvpr, isplitp
        if (wtrvpr.gt.0) THEN
            ssh2o_s(i,1) = vapvrt(i) / vap_sol
            ssh2o_s(i,2) = vapvrt(i) / vap_sen
            write(91,*) 'SSH2O_S(I,1), I = ', ssh2o_s(i,1), i, solzni*radeg
            write(91,*) 'SSH2O_S(I,2), I = ', ssh2o_s(i,2), i, obszni*radeg
        else
            if (isplitp.ne.0) then
                write(6,*) 'ATREM: Split paths is not working because WaterVapor is 0'
                stop 0
            end if
        end if
  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(nh2o_max)
      dimension wavobs(nobs_max),fwhm(nobs_max)
      dimension dp(25), pm(25), tm(25), vmrm(25)
      dimension finst2(100)
      dimension sumcf(np_hi)
 
!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,ifullcalc,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
 
      REAL, ALLOCATABLE :: ABSCF_CO2(:,:),ABSCF_O2(:,:),ABSCF_N2O(:,:),ABSCF_CH4(:,:),ABSCF_CO(:,:)
      dimension g_vap(25), g_other(25)
      COMMON /geometry3/ g_vap, g_other, g_vap_equiv
      COMMON /geometry4/vap_slant_mdl
 
      dimension o3cf(no3pt)
      COMMON /o3cf_init1/ o3cf
 
      dimension tran_o3_std(no3pt)
      COMMON /init_speccal16/ tran_o3_std
 
      dimension rno2cf(no3pt)
      COMMON /no2cf_init1/ rno2cf
 
      dimension tran_no2_std(no3pt)
      COMMON /init_speccal17/ tran_no2_std
 
      COMMON /model_adj4/ k_surf
      INTEGER start(2)/1,1/
      INTEGER cnt(2)/NP_HI,19/
      CHARACTER(len=4096) :: filename
      INTEGER FINDMATCH
      external findMatch
      INTEGER FIRST/1/
      SAVE first,abscf_co2,abscf_o2,abscf_n2o,abscf_ch4,abscf_co
 
      IF (first.eq.1) THEN
 
!      CALL get_environment_variable("OCDATAROOT", datadir)
!      WRITE (*,*) 'OCDATAROOT='//TRIM(homedir)
 
      nh2o = nh2o_max         !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
 
!      TRAN_HI(:,:)  = 1.0
 
!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,
!C         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,
!C         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)
          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
      END IF
 
      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)
      iwndw1c = findmatch(wavobs,nobs,wndow1)
      write(91,*) 'IWNDW1=',iwndw1,iwndw1c,' list[0]=',wavobs(1)
      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)
 
      ALLOCATE ( abscf_co2(np_hi,19) )
      ALLOCATE ( abscf_n2o(np_hi,19) )
      ALLOCATE ( abscf_co(np_hi,19) )
      ALLOCATE ( abscf_ch4(np_hi,19) )
      ALLOCATE ( abscf_o2(np_hi,19) )
 
       write(filename(1:dln),'(a)') datpath(1:dln)
       write(filename(dln+1:),'(a)') 'abscf_gas.nc'
       ncid = ncopn(filename,ncnowrit,ircode)
 
       nrhid = ncvid(ncid, 'abscf_co2', ircode)
       CALL ncvgt (ncid, nrhid, start, cnt, abscf_co2, ircode)
        if (ircode .ne.0) then
          write(*,*) 'Error reading abscf_gas.nc: abscf_co2: rcode=',ircode
          stop 0
       end if
 
       nrhid = ncvid(ncid, 'abscf_n2o', ircode)
       CALL ncvgt (ncid, nrhid, start, cnt, abscf_n2o, ircode)
        if (ircode .ne.0) then
          write(*,*) 'Error reading abscf_gas.nc: abscf_n2o: rcode=',ircode
          stop 0
       end if
 
       nrhid = ncvid(ncid, 'abscf_co', ircode)
       CALL ncvgt (ncid, nrhid, start, cnt, abscf_co, ircode)
        if (ircode .ne.0) then
          write(*,*) 'Error reading abscf_gas.nc: abscf_co: rcode=',ircode
          stop 0
       end if
 
       nrhid = ncvid(ncid, 'abscf_ch4', ircode)
       CALL ncvgt (ncid, nrhid, start, cnt, abscf_ch4, ircode)
        if (ircode .ne.0) then
          write(*,*) 'Error reading abscf_gas.nc: abscf_ch4: rcode=',ircode
          stop 0
       end if
 
       nrhid = ncvid(ncid, 'abscf_o2', ircode)
       CALL ncvgt (ncid, nrhid, start, cnt, abscf_o2, ircode)
        if (ircode .ne.0) then
          write(*,*) 'Error reading abscf_gas.nc: abscf_o2: rcode=',ircode
          stop 0
       end if
 
    CALL ncclos(ncid, rcode)
 
    END IF ! IFIRST = 1
!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:
 
       tran_hi_others(:) = 1.0
 
 
!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
  
       sumcf(:) = 0.0
       DO j = k_surf, 19
          sumcf(:) = sumcf(:) - abscf_co2(:,j)*dp(j-k_surf+1)*vmrm(j-k_surf+1)
       END DO
 
        tran_hi_others(:)  = tran_hi_others(:)*exp(sumcf(:)*q* 28.966 / &
                              6.0225e+23 / 1.0e-06)
      END IF
 
!C--------------------------------------------
!C       For N2O transmittance calculation.
 
      IF(in2o.EQ.1) THEN
!          if (first.eq.1) print*,"Doing N2O..."
          DO 324 i=1,nl
            vmrm(i)=0.3*1.0e-06
            vmrm(i)= vmrm(i)*g_other(i)
  324     CONTINUE
  
 
      sumcf(:) = 0
       DO j = k_surf, 19
           sumcf(:) = sumcf(:) - abscf_n2o(:,j)*dp(j-k_surf+1)*vmrm(j-k_surf+1)
       END DO
 
        tran_hi_others(:)  = tran_hi_others(:)*exp(sumcf(:)*q* 28.966 / &
                              6.0225e+23 / 1.0e-06)
!        if (first.eq.1) then
!        do i=1,NP_HI
!        print*,i," Tran_hi_others=",TRAN_HI_OTHERS(i)
!        end do
!        endif
      END IF
 
!C--------------------------------------------
!C       For CO transmittance calculation.
      IF(ico.EQ.1) THEN
!          if (first.eq.1) print*,"Doing CO..."
          DO 325 i=1,nl
            vmrm(i)=0.1*1.0e-06
            vmrm(i)= vmrm(i)*g_other(i)
  325     CONTINUE
  
       sumcf(:) = 0.0
       DO j = k_surf, 19
           sumcf(:) = sumcf(:) - abscf_co(:,j)*dp(j-k_surf+1)*vmrm(j-k_surf+1)
       END DO
        tran_hi_others(:)  = tran_hi_others(:)*exp(sumcf(:)*q* 28.966 / &
                              6.0225e+23 / 1.0e-06)
!        if (first.eq.1) then
!        do i=1,NP_HI
!        print*,i," Tran_hi_others=",TRAN_HI_OTHERS(i)
!        end do
!        endif
      END IF
 
!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
 
 
       sumcf(:) = 0.0
       DO j = k_surf, 19
          sumcf(:) = sumcf(:) - abscf_ch4(:,j)*dp(j-k_surf+1)*vmrm(j-k_surf+1)
       END DO
 
        tran_hi_others(:)  = tran_hi_others(:)*exp(sumcf(:)*q* 28.966 / &
                              6.0225e+23 / 1.0e-06)
 
      END IF
 
!C--------------------------------------------
!C       For O2 transmittance calculation.
      IF(io2.EQ.1) THEN
            sumcf(:) = 0
!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
 
 
!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.
       DO j = k_surf, 19
          sumcf(:) = sumcf(:) - abscf_o2(:,j)*dp(j-k_surf+1)*vmrm(j-k_surf+1)
 
       END DO
       tran_hi_others(1:9000) = tran_hi_others(1:9000) * exp(sumcf(1:9000)*q* 28.966 / &
                              6.0225e+23 / 1.0e-06)
       tran_hi_others(9001:106600) = tran_hi_others(9001:106600) * exp(sumcf(9001:106600)*q*scl_o2* 28.966 / &
                              6.0225e+23 / 1.0e-06)
       tran_hi_others(106601:np_hi) = tran_hi_others(106601:np_hi) * exp(sumcf(106601:np_hi)*q* 28.966 / &
                              6.0225e+23 / 1.0e-06)
 
      END IF
 
 
!C--------------------------------------------
!C End of calculation of high resolution transmittances for CO2, N2O, CO,
!C     CH4, and O2.
!C--------------------------------------------
!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
     
      first = 0
      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                                                *
!*                                         *
!********************************************************************************
!********************************************************************************
!*  TRANCAL combined with TRAN_TABLE by R. Healy 4/28/2015                      *
!*  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 tran_table
 
      include 'COMMONS_INC.f'
 
!C  Common variables
      dimension dp(25), pm(25), tm(25), vmrm(25), g_vap(25),g_other(25)
      dimension h(25), t(25), p(25), vmr(25)
      dimension wavobs(nobs_max),fwhm(nobs_max)
      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,ifullcalc,hsurf,dlt,dlt2
      COMMON /init_speccal3/ nh2o
      COMMON /init_speccal5/ dp,pm,tm,vmrm
      COMMON /model_adj1/ clmvap,q
      COMMON /model_adj4/ k_surf
      COMMON /model_adj3/ k_plane, dvap_plane, dvap_layer, &
                          dp_plane, dp_layer, clmvapp
 
      dimension ssh2o(nh2o_max)
      dimension vaptot(nh2o_max), r0p94(nh2o_max), r1p14(nh2o_max), trntbl(nobs_max,nh2o_max), &
                 tran_kd(nobs_max,nh2o_max), diff_tran(nobs_max,nh2o_max),trh2(nobs_max,nh2o_max), &
                 trntblo(nobs_max)
      COMMON /tran_table1/ sh2o,vaptot,r0p94,r1p14,trntbl,tran_kd, diff_tran,trntblo
      ! Split path transmittance
      ! Just need the values for the four table entries calculated on call to get_atrem
      ! rjh - 1/11/2016
      REAL   TRAN_HI_SA(NP_HI,2),TRAN_HI_SAP1(NP_HI,2),TRAN_HI_SB(NP_HI,2),TRAN_HI_SBP1(NP_HI,2), &
              trntbl_s(nobs_max,2),trntbl_so(nobs_max,2)
      COMMON /tran_tables/tran_hi_sa,tran_hi_sap1,tran_hi_sb,tran_hi_sbp1
      COMMON /tran_table_l2gen/trntbl_s,trntbl_so
      !COMMON /TRANCAL1/ VAPTT
 
      REAL   MU,MU0,SSH2O_S(NH2O_MAX,2),VMRM_S(25,2)
      COMMON /geometry5/mu,mu0,ssh2o_s
 
      COMMON /geometry2/ gco2,go3,gn2o,gco,gch4,go2,ssh2o,totlo3,ggeom
      COMMON /geometry3/ g_vap, g_other, g_vap_equiv
      COMMON /geometry4/vap_slant_mdl
 
      INTEGER start(3)
      INTEGER cnt(3)
      CHARACTER(len=4096) :: filename
      dimension sumcf(np_hi), sumcf_s(np_hi,2)
      REAL ABSCF_H2O(NP_HI,19)
 
      REAL, ALLOCATABLE :: TG(:), TKCDF(:,:,:),DIFFT(:,:),SUM_KD(:),F1(:),F2(:),TKCDFC(:,:,:)
      REAL, ALLOCATABLE :: WAVNO_HI(:)
      CHARACTER*31 GNAM, BNAM,LNAM
 
      REAL, ALLOCATABLE :: TDP(:),TVMRM(:), TRAN_HI(:,:)
      SAVE tdp, tvmrm
 
      INTEGER FIRST/1/,IDOSMOOTH/0/
      REAL WV/4.1294627/
      LOGICAL DOINLINEKD/.TRUE./
      SAVE first, abscf_h2o, tg, tkcdf, difft,n_g, wv, idosmooth, doinlinekd
      data nl/19/,ng/7/
 
       ALLOCATE( tran_hi(np_hi,nh2o_max) )
       ALLOCATE( wavno_hi(np_hi) )
 
       if (first.eq.1) then
 
           if (ifullcalc.eq.1) then
              print*,'ATREM: **WARNING** Full_calc is on.  Processing will be extremely slow...'
           else
              print*,'ATREM: Full_calc is off.  Processing using k-distribution method...'
           endif
 
           idosmooth = ico2 + in2o + ico + ich4 + io2 + io3 + ino2
 
           write(filename(1:dln),'(a)') datpath(1:dln)
           write(filename(dln+1:),'(a)') 'abscf_gas.nc'
 
           start(1) = 1
           start(2) = 1
           start(3) = 1
           cnt(1)   = np_hi
           cnt(2)   = 19
           cnt(3)   = 0
           ncid = ncopn(filename,ncnowrit,ircode)
           nrhid = ncvid(ncid, 'abscf_h2o', ircode)
           CALL ncvgt (ncid, nrhid, start, cnt, abscf_h2o, ircode)
            if (ircode .ne.0) then
              write(*,*) 'Error reading abscf_gas.nc: abscf_h2o: rcode=',ircode
              stop 0
           end if
           cnt(2)   = 1
           nrhid = ncvid(ncid, 'waveno', ircode)
           CALL ncvgt (ncid, nrhid, start, cnt, wavno_hi, ircode)
            if (ircode .ne.0) then
              write(*,*) 'Error reading abscf_gas.nc: waveno: rcode=',ircode
              stop 0
           end if
 
           CALL ncclos(ncid, rcode)
 
           if (.not.doinlinekd) then
               write(filename(1:dln),'(a)') datpath(1:dln)
               write(filename(dln+1:),'(a)') 'hico_atrem_h2o_coef7.nc'
 
               ncid = ncopn(filename,ncnowrit,ircode)
 
               CALL ncinq(ncid,ndims,nvars,natts,irecdim,ircode)
 
               if (ndims.NE.3) then
                  write(*,*) 'Error Wrong number of dims in file'//trim(filename)// &
                  'Expected 4, but got ',ndims
                  stop 0
               end if
 
               CALL ncdinq(ncid, 1, gnam, n_g, ircode)
 
               CALL ncdinq(ncid, 2, bnam, n_b, ircode)
 
               CALL ncdinq(ncid, 3, lnam, n_l, ircode)
 
               !CALL NCDINQ(NCID, 4, LNAM, N_H, IRCODE)
 
               if (n_b.NE.nobs) then
                  write(*,*) ' Bands expected (',nobs,') do not match number of bands in file'// &
                 '(',n_b,') for file='//filename
                  stop 0
               end if
 
    !           if (N_H.NE.NH2O) then
    !              write(*,*) ' NH2O expected (',NH2O,') do not match NH2O in file'// &
    !             '(',N_H,') for file='//filename
    !              stop 0
    !           end if
            else
                n_g = 7
                n_l = 19
                n_b = nobs
            endif ! doinlinekd
 
           ALLOCATE( tg(n_g) )
           ALLOCATE( tkcdf(n_l,n_b,n_g) )
           !ALLOCATE( DIFFT(N_B,N_H) )
           ALLOCATE( f1(n_b) )
           ALLOCATE( f2(n_b) )
           ALLOCATE( sum_kd(n_b))
 
           ALLOCATE( tvmrm(n_l) )
           ALLOCATE( tdp(n_l) )
 
           if (doinlinekd) then
               tg(1:7) = (/ 0.0, 0.379, 0.5106, 0.81, 0.9548, 0.9933, 1.0/)
           else
 
               start(1) = 1
               start(2) = 1
               start(3) = 1
               cnt(1)   = n_g
               cnt(2)   = 0
               cnt(3)   = 0
 
               nrhid = ncvid(ncid, 'g', ircode)
               CALL ncvgt (ncid, nrhid, start, cnt, tg, ircode)
                if (ircode .ne.0) then
                  write(*,*) 'Error reading '//filename//': rcode=',ircode
                  stop 0
               end if
 
               start(1) = 1
               start(2) = 1
               start(3) = 1
               cnt(1)   = n_l
               cnt(2)   = n_b
               cnt(3)   = n_g
 
               nrhid = ncvid(ncid, 'k_h2o', ircode)
               CALL ncvgt (ncid, nrhid, start, cnt, tkcdf, ircode)
                if (ircode .ne.0) then
                  write(*,*) 'Error reading '//filename//': rcode=',ircode
                  stop 0
               end if
 
    !           start(1) = 1
    !           start(2) = 1
    !           start(3) = 1
    !           cnt(1)   = N_B
    !           cnt(2)   = N_H
    !           cnt(3)   = 0
 
    !           NRHID = NCVID (NCID, 'diffT', IRCODE)
    !           CALL NCVGT (NCID, NRHID, START, CNT, DIFFT, IRCODE)
    !            if (ircode .ne.0) then
    !              write(*,*) 'Error reading '//filename//': rcode=',ircode
    !              stop 0
    !           end if
 
            endif
 
 
      endif
 
! Moved from INIT_SPECAL because it seems more logical here. - rjh 01/04/2016
!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.
         if (isplitp.ne.0) then      ! perform the split path calculations
             vmrm_s(i,1) = vmrm(i)*mu0
             if (i.lt.k_plane) THEN
                vmrm_s(i,2) = vmrm(i)*mu
             else
                if (i.eq.k_plane) THEN
                    vmrm_s(i,2) = vmrm(i)*mu*dvap_plane/dvap_layer
                else
                    vmrm_s(i,2) = 0
                end if
             END IF
             WRITE(91,*) i,' VMRM SOL, OBS=',vmrm_s(i,1),vmrm_s(i,2)
         end if
         vmrm(i) = vmrm(i)*g_vap(i)
         WRITE(91,*) 'VMRM=',vmrm(i)
      END DO
 
      if (ifullcalc.ne.0.or.first.eq.1) then
 
        !C For each water vapor amount, calculate the geometrically adjusted water
        !C     vapor amount, the channel ratios, and the transmittance spectrum.
              sumcf(:)    = 0
              IF (isplitp.NE.0) THEN
                sumcf_s(:,1)    = 0
                sumcf_s(:,2)    = 0
              END IF
              DO j = k_surf, 19
                 sumcf(:) = sumcf(:) - abscf_h2o(:,j)*dp(j-k_surf+1)*vmrm(j-k_surf+1)
                 IF (isplitp.NE.0) THEN
                    sumcf_s(:,1) = sumcf_s(:,1) - abscf_h2o(:,j)*dp(j-k_surf+1)*vmrm_s(j-k_surf+1,1)
                    sumcf_s(:,2) = sumcf_s(:,2) - abscf_h2o(:,j)*dp(j-k_surf+1)*vmrm_s(j-k_surf+1,2)
                 END IF
              END DO
 
              DO  i=1,nh2o
        !C Calculate water vapor transmittances with different scaling factors:
        !C
        !C Multiplying the high resolution water vapor transmittance with combined
        !C    transmittances of CO2, N2O, CO, CH4, and O2:
        !C
 
               tran_hi(:,i) = exp(sumcf(:)*ssh2o(i) * q * 28.966 / &
                            6.0225e+23 / 1.0e-06)            ! *TRAN_HI_OTHERS(:) - Done in TRAN_SMOOTH_OTHERS
        !C Total amount of water vapor (in unit of cm) corresponding to the spectrum.
               vaptot(i)= vap_slant_mdl * ssh2o(i)
!         write(91,*) 'VAPTT=',VAPTOT(I),SSH2O(I),VAP_SLANT_MDL
 
           END DO
 
           if (isplitp.ne.0) then      ! perform the split path calculations
                DO j=1,2
                    tran_hi_sa(:,j) = exp(sumcf_s(:,j)  *ssh2o_s(ja_l2,j) * q * 28.966 / &
                            6.0225e+23 / 1.0e-06)
                    tran_hi_sap1(:,j) = exp(sumcf_s(:,j)*ssh2o_s(ja_l2+1,j) * q * 28.966 / &
                            6.0225e+23 / 1.0e-06)
 
                    tran_hi_sb(:,j) = exp(sumcf_s(:,j)  *ssh2o_s(jb_l2,j) * q * 28.966 / &
                            6.0225e+23 / 1.0e-06)
                    tran_hi_sbp1(:,j) = exp(sumcf_s(:,j)*ssh2o_s(jb_l2+1,j) * q * 28.966 / &
                            6.0225e+23 / 1.0e-06)
                END DO
!                    write(91,*) 'SPLIT: JA,JB=',JA_L2,JB_L2,'TRAN_HI=',TRAN_HI(1,JA_L2), &
!                    ' TRAN_HI_SA:',TRAN_HI_SA(1,1),TRAN_HI_SA(1,2), &
!                    TRAN_HI_SAP1(1,1),TRAN_HI_SAP1(1,2),TRAN_HI_SB(1,1),TRAN_HI_SB(1,2), &
!                    TRAN_HI_SBP1(1,1),TRAN_HI_SBP1(1,2)
!                    write(91,*)  'SSH2O JAJB:',SSH2O_S(JA_L2,1),SSH2O_S(JA_L2,2),SSH2O_S(JB_L2,1),SSH2O_S(JB_L2,2), ' SOL=',    &
!                    SUMCF_S(1,1), SSH2O_S(JA_L2,1), SUMCF_S(1,1)  *SSH2O_S(JA_L2,1) * Q * 28.966 / 6.0225E+23 / 1.0E-06,' SEN=', &
!                    SUMCF_S(1,2), SSH2O_S(JA_L2,2), SUMCF_S(1,2)  *SSH2O_S(JA_L2,2) * Q * 28.966 / 6.0225E+23 / 1.0E-06,' SOL=', &
!                    exp(SUMCF_S(1,1)  *SSH2O_S(JA_L2,1) * Q * 28.966 / 6.0225E+23 / 1.0E-06),' SEN=', &
!                    exp(SUMCF_S(1,2)  *SSH2O_S(JA_L2,2) * Q * 28.966 / 6.0225E+23 / 1.0E-06)
 
           endif
 
     endif
     if (ifullcalc.eq.0) then
! The following code calculates the k-distrubition coefficients
! as opposed to reading them in from a netcdf file.
! -rjh 9/23/2015
        if (first.eq.1.AND.doinlinekd) then
           ALLOCATE( tkcdfc(n_l,n_b,n_g) )
           call kdist_gas_abs(tkcdfc,abscf_h2o,np_hi,wavno_hi,wavobs,n_b)
           tkcdf = tkcdfc
        end if
        !  Perform the k-distribution calculation Trapezoidal integral for transmittance over the NOBS(NBANDS) bands
           DO  i=1,nh2o
              tran_kd(:,i) = 1.0
 
              DO j = k_surf, 19
 
                 sum_kd(:) = 0.0
 
                 f1(:) = exp(-tkcdf(j,:,1)  *dp(j-k_surf+1)*vmrm(j-k_surf+1)*ssh2o(i))
 
                 DO k = 1, n_g-1
                    f2(:) = exp(-tkcdf(j,:,k+1)*dp(j-k_surf+1)*vmrm(j-k_surf+1)*ssh2o(i))
                    sum_kd(:) = sum_kd(:)+(f1(:) + f2(:))*(tg(k+1) - tg(k))/2.0
                    f1(:) = f2(:)
                 END DO
 
                 tran_kd(:,i) = tran_kd(:,i)*sum_kd(:)
 
              END DO
 
        !C Total amount of water vapor (in unit of cm) corresponding to the spectrum.
               vaptot(i)= vap_slant_mdl * ssh2o(i)
           END DO
 
     endif
 
      if (ifullcalc.ne.0.or.first.eq.1) then
 
!C
!C Smooth the high resolution spectra to resolutions of measured spectrum
!C
 
            CALL tran_smooth(tran_hi)
 
            if (first.eq.1) then
                 DO  i=1,nh2o
                     trh2(:,i) = trntbl(:,i)
 !                DO J=1,NOBS
 !                   write(6,*) 'RJH: TRH2: ',TRH2(J,I)
 !                END DO
                 END DO
            endif
 
       endif
 
       trntblo(:) = 1.0
       trntbl_so(:,:) = 1.0
       ! Only call the TRAN_SMOOTH_OTHERS if the other trace gases were selected
       IF (idosmooth.GT.0) THEN
         CALL tran_smooth_others
       ENDIF
 
       if (ifullcalc.eq.0) then
 
           if (first.eq.1) then
!               do k=1,NOBS
!                    print*,K," FAST)TRNTBL=",TRAN_KD(k,60)*DIFF_TRAN(k,1),TRNTBLO(k)
!               end do
           end if
           DO  i=1,nh2o
               trntbl(:,i) = tran_kd(:,i)*diff_tran(:,i)*trntblo(:)
 
!                if (first.eq.1) then
!                 DO J=1,NOBS
!                    write(6,*) 'RJH: DIFF_TRAN: ',I,J,TRNTBL(J,I),TRH2(J,I),TRAN_KD(J,I),DIFF_TRAN(J,I)
!                 END DO
!                endif
           END DO
 
       ELSE
           if (first.eq.1) then
               do k=1,nobs
                    print*,k," SLOW)TRNTBL=",trntbl(k,60),trntblo(k)
               end do
           end if
 
           DO  i=1,nh2o
               trntbl(:,i) = trntbl(:,i)*trntblo(:)
           END DO
 
           trntbl_s(:,:) = trntbl_s(:,:)*trntbl_so(:,:)
 
       endif
 
!C Calculate 3-channel ratio values, R094 and R114, from simulated spectra.
       CALL chnlratio
 
       first=0
       DEALLOCATE(tran_hi,wavno_hi)
      RETURN
      END
 
      subroutine linterp(x,y,n,xi,yi)
      real*4 x(n),y(n),xi,yi
      integer n
 
      klo=1
      khi=n
      k=0
      do while (khi-klo.gt.1)
        k=(khi+klo)/2
        if(x(k).gt.xi)then
          khi=k
        else
          klo=k
        endif
      end do
 
      if (khi.eq.klo) then  !extrapolate
         klo = klo - 1
         if (klo.le.0) then
            !print*,"Ooops..."
            klo=1
            !stop 0
         endif
         yi = y(khi) + (xi-x(khi))/(x(khi)-x(klo))*(y(khi)-y(klo))
      else
         yi = y(klo) + (xi-x(klo))/(x(khi)-x(klo))*(y(khi)-y(klo))
      endif
 
      end
 
      SUBROUTINE kdist_gas_abs(kcdf,abscf,np_hi,waveno,wavobs,nwave)
 
      ! R. Healy 9/2015
      ! Reconstruct Transmittances from the k coefficients (from Amir Ibrahim's Matlab Code)
 
      real*4 :: abscf(np_hi,19),waveno(np_hi),wavobs(nwave),kcdf(19,nwave,7)
      real*4 alayers(np_hi,19)
      REAL, ALLOCATABLE :: diflam(:),UV_lam(:),IR_lam(:),uv_dlam(:),ir_dlam(:),lam(:), &
                             dwave(:),dwn(:),wn(:),g_i(:),k_i(:),g(:),kg(:),logk(:),kint(:,:), &
                             k7(:,:,:),kk(:)
      INTEGER, ALLOCATABLE :: wavel_window(:,:)
      INTEGER binnum(19)
      real*4 g7(7)/ 0,    0.379,    0.6,    0.81,    0.9548,    0.9933,    1/
      character*1 x /''/
      integer*1 inul
      equivalence(inul,x)
 
      data q/2.15199993e+25/  ! Q=# of molecules above the surface at one atmosphere (molecules/cm**2)
 
      do i=1,19
 
        alayers(:,i) = abscf(:,i)*q*28.966/ &
                            6.0225e+23 / 1.0e-06
      end do
! wavobs - array of center wavelengths of instrument
! DLAM_cent - delta center wavelengths
 
      ALLOCATE( diflam(nwave) )
 
      dlam = 9999
      wmin = 9999
      wmax = -9999
      do i=1,nwave-1
        diflam(i) = wavobs(i+1) - wavobs(i)
        if (dlam.gt.diflam(i)) dlam = diflam(i)
        if (wmin.gt.wavobs(i)) wmin = wavobs(i)
        if (wmax.lt.wavobs(i)) wmax = wavobs(i)
      end do
 
      diflam(nwave) = diflam(nwave-1)
 
      if (wmin.gt.wavobs(nwave)) wmin = wavobs(nwave)
      if (wmax.lt.wavobs(nwave)) wmax = wavobs(nwave)
 
      ndxwuv = (wmin-0.0001 - 0.3)/dlam + 1
      ndxwir = (3.1- (wmax+0.0001))/diflam(nwave) + 1
      ndxtot = nwave+ndxwuv+ndxwir
 
      ALLOCATE( uv_lam(ndxwuv) )
      ALLOCATE( ir_lam(ndxwir) )
      ALLOCATE( uv_dlam(ndxwuv) )
      ALLOCATE( ir_dlam(ndxwir) )
      ALLOCATE( lam(ndxtot) )
      ALLOCATE( dwave(ndxtot) )
      ALLOCATE( dwn(ndxtot) )
      ALLOCATE( wn(ndxtot) )
      ALLOCATE( wavel_window(ndxtot,2) )
 
      k = ndxtot
      swav = 0.3
 
      do i=1,ndxwuv
        uv_lam(i) = swav
        uv_dlam(i) = diflam(1)
        swav = swav + dlam
        lam(k) = uv_lam(i)
        dwn(k) = 10000.*dlam/(lam(k)**2)
         wn(k) = 10000./lam(k)
        k = k - 1
      end do
      swav = wmin
      ks=k
      do i=1,nwave
        lam(k) = swav
        dwn(k) = 10000.*diflam(i)/(lam(k)**2)
        wn(k) = 10000./lam(k)
        swav = swav + diflam(i)
        k = k - 1
      end do
 
      swav = wmax+0.0001
 
      do i=1,ndxwir
        ir_lam(i) = swav
        ir_dlam(i) = diflam(nwave)
        swav = swav + diflam(nwave)
        lam(k) = ir_lam(i)
        dwn(k) = 10000.*diflam(nwave)/(lam(k)**2)
         wn(k) = 10000./lam(k)
        k = k - 1
      end do
 
      ALLOCATE( kint(ndxtot,7) )
      ALLOCATE( k7(19,ndxtot,7) )
 
      iw = nwave ! index for wavelngths for output array kcdf
 
      do i=ndxwir+1,ndxtot-ndxwuv
 
         kwavdn = -999
         kwavup = -999
 
         do k=1,np_hi
           if (kwavdn.lt.0.and.waveno(k).gt.(wn(i)-dwn(i)/2.0)) kwavdn = k
           if (kwavup.lt.0.and.waveno(k).gt.(wn(i)+dwn(i)/2.0)) kwavup = k
         end do
 
         wavel_window(i,1) = kwavdn
         wavel_window(i,2) = kwavup
 
         if (kwavup.lt.1.or.kwavdn.lt.1) then
            nsamp = 1
         else
            nsamp = kwavup - kwavdn + 1
         end if
 
         nbins = nsamp
         ALLOCATE( g_i(nbins+1) )
         ALLOCATE( k_i(nbins+1) )
         ALLOCATE( kk(nsamp+1) )
 
         do k=1,19
            if (kwavdn.lt.1.or.kwavup.lt.1) then
               kk(:) = 0
            else
               kk(1:nsamp) = alayers(kwavdn:kwavup,k)
            end if
 
            binnum(k) = nbins
            call ecdf(k_i,g_i,binnum(k),kk,nsamp)
 
            ALLOCATE( logk(nbins+1) )
            logk(:) = log10(k_i(:))
 
            do l=1,7
 
 
              call linterp(g_i,logk,binnum(k),g7(l),kint(i,l))
              k7(k,i,l) = 10**kint(i,l)
 
              if (i.gt.ndxwir.and.i.le.ndxtot-ndxwuv) then
                  if (isnan(k7(k,i,l))) then
                     k7(k,i,l) = 0
                  endif
                  kcdf(k,iw,l) = k7(k,i,l)
              end if
            end do
 
            DEALLOCATE(logk)
 
         end do
         DEALLOCATE (g_i, stat=ialloerr)
         DEALLOCATE (k_i)
         DEALLOCATE (kk)
 
         if (i.gt.ndxwir.and.i.le.ndxtot-ndxwuv) iw = iw - 1
      end do
 
      DEALLOCATE (k7)
      DEALLOCATE (kint)
      DEALLOCATE( uv_lam )
      DEALLOCATE( ir_lam )
      DEALLOCATE( uv_dlam )
      DEALLOCATE( ir_dlam )
      DEALLOCATE( lam )
      DEALLOCATE( dwave )
      DEALLOCATE( dwn )
      DEALLOCATE( wn )
      DEALLOCATE( wavel_window )
 
      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(TRAN_HI)
 
      include 'COMMONS_INC.f'
 
      REAL :: TRAN_HI(NP_HI,NH2O_MAX)
      dimension wavobs(nobs_max),fwhm(nobs_max)
      COMMON /getinput4/ wavobs,fwhm
 
      dimension vaptot(nh2o_max), r0p94(nh2o_max), r1p14(nh2o_max), trntbl(nobs_max,nh2o_max), &
                 tran_kd(nobs_max,nh2o_max), diff_tran(nobs_max,nh2o_max),trntblo(nobs_max)
      COMMON /tran_table1/ sh2o,vaptot,r0p94,r1p14,trntbl,tran_kd, diff_tran,trntblo
 
      dimension tran_ia(nh2o_max),tran_iap1(nh2o_max)
 
      COMMON /getinput5/ nobs,ifullcalc,hsurf,dlt,dlt2
      REAL, ALLOCATABLE :: TRAN_STD(:,:),TRAN_MED(:,:)
 
      REAL   TRAN_HI_SA(NP_HI,2),TRAN_HI_SAP1(NP_HI,2),TRAN_HI_SB(NP_HI,2),TRAN_HI_SBP1(NP_HI,2), &
              trntbl_s(nobs_max,2),trntbl_so(nobs_max,2)
      COMMON /tran_tables/tran_hi_sa,tran_hi_sap1,tran_hi_sb,tran_hi_sbp1
      COMMON /tran_table_l2gen/trntbl_s,trntbl_so
      REAL   TRAN_MED_INDEX_SA(NP_MED,2),TRAN_MED_INDEX_SB(NP_MED,2), &
              tran_med_index_sap1(np_med,2),tran_med_index_sbp1(np_med,2), &
              tran_med_sa(np_med,2),tran_med_sap1(np_med,2),tran_med_sb(np_med,2),tran_med_sbp1(np_med,2), &
              tran_std_sa(np_std,2),tran_std_sap1(np_std,2),tran_std_sb(np_std,2),tran_std_sbp1(np_std,2),             &
              tran_ia_sa(2),tran_ia_sap1(2),tran_ia_sb(2),tran_ia_sbp1(2), &
              tran_iap1_sa(2),tran_iap1_sap1(2),tran_iap1_sb(2),tran_iap1_sbp1(2), &
              trntblsa(2),trntblsap1(2),trntblsb(2),trntblsbp1(2)
      COMMON /tran_tables1/tran_med_index_sa,tran_med_index_sap1,tran_med_index_sb,tran_med_index_sbp1, &
              tran_med_sa,tran_med_sap1,tran_med_sb,tran_med_sbp1, &
              tran_std_sa,tran_std_sap1,tran_std_sb,tran_std_sbp1
 
      INTEGER FIRST/1/,IA(NOBS_MAX)
      SAVE first, ia
 
!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
    IF (first.eq.1) THEN !RJH
      DO 466 j =1, np_med
 
             ncvtot_wavno = 2 * ncvhf_wavno(j) - 1
    
             sumins = 0.0
 
          DO 560 i = ncvhf_wavno(j), ncvtot_wavno
             finstr_wavno(i,j) = &
                exp( -const1*(float(i-ncvhf_wavno(j))*dwavno &
                /fwhm_wavno(j))**2)
             sumins = sumins + finstr_wavno(i,j)
 560      CONTINUE
 
          DO 565 i = 1, ncvhf_wavno(j)-1
             finstr_wavno(i,j) = finstr_wavno(ncvtot_wavno-i+1,j)
             sumins = sumins + finstr_wavno(i,j)
 565      CONTINUE
 
          sumins = sumins * dwavno
 
          DO 570 i = 1, ncvtot_wavno
             finstr_wavno(i,j) = finstr_wavno(i,j)*dwavno/sumins
 570      CONTINUE
  
 
 466  CONTINUE
 
 
!C  Index searching...
!C
       DO j=1,nobs
         !C Calculate instrumental response functions...
 
          sumins = 0.0
 
          ncvtot = 2 * ncvhf(j) - 1
 
          DO 1560 i = ncvhf(j), ncvtot
             finstr(i,j) = &
                exp( -const1*(float(i-ncvhf(j))*dwavln &
                /fwhm(j))**2)
             sumins = sumins + finstr(i,j)
1560      CONTINUE
 
          DO 1565 i = 1, ncvhf(j)-1
             finstr(i,j) = finstr(ncvtot-i+1,j)
             sumins = sumins + finstr(i,j)
1565      CONTINUE
 
          sumins = sumins * dwavln
 
          DO 1570 i = 1, ncvtot
             finstr(i,j) = finstr(i,j)*dwavln/sumins
1570      CONTINUE
 
 
          CALL hunt(wavln_std, np_std, wavobs(j), ia(j))
       END DO
 
      first = 0
      END IF
 
       ALLOCATE(tran_std(np_std,nh2o_max))
       ALLOCATE(tran_med(np_med,nh2o_max))
!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.
 
     tran_med_index(:,:)  = 0.0
     IF(isplitp.NE.0) THEN
         tran_med_index_sa(:,:)    = 0.0
         tran_med_index_sap1(:,:)  = 0.0
         tran_med_index_sb(:,:)    = 0.0
         tran_med_index_sbp1(:,:)  = 0.0
     END IF
 
 
     DO  j =1, np_med
 
          ndx1 = index_med(j)-(ncvhf_wavno(j)-1)
          ndx2 = index_med(j)+ ncvhf_wavno(j)-1
          DO 491 k = ndx1,ndx2
             tran_med_index(j,:) = tran_med_index(j,:) + tran_hi(k,:)*  &
                          finstr_wavno(k-index_med(j)+ncvhf_wavno(j),j)
              IF(isplitp.NE.0) THEN
                 tran_med_index_sa(j,:) = tran_med_index_sa(j,:) + tran_hi_sa(k,:)*  &
                          finstr_wavno(k-index_med(j)+ncvhf_wavno(j),j)
                 tran_med_index_sap1(j,:) = tran_med_index_sap1(j,:) + tran_hi_sap1(k,:)*  &
                          finstr_wavno(k-index_med(j)+ncvhf_wavno(j),j)
                 tran_med_index_sb(j,:) = tran_med_index_sb(j,:) + tran_hi_sb(k,:)*  &
                          finstr_wavno(k-index_med(j)+ncvhf_wavno(j),j)
                 tran_med_index_sbp1(j,:) = tran_med_index_sbp1(j,:) + tran_hi_sbp1(k,:)*  &
                          finstr_wavno(k-index_med(j)+ncvhf_wavno(j),j)
              END IF
 491      CONTINUE
     END DO
 
!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,:)
         IF(isplitp.NE.0) THEN
            tran_med_sa(1,:)        = tran_med_index_sa(1,:)
            tran_med_sa(np_med,:)   = tran_med_index_sa(np_med,:)
            tran_med_sap1(1,:)      = tran_med_index_sap1(1,:)
            tran_med_sap1(np_med,:) = tran_med_index_sap1(np_med,:)
            tran_med_sb(1,:)        = tran_med_index_sb(1,:)
            tran_med_sb(np_med,:)   = tran_med_index_sb(np_med,:)
            tran_med_sbp1(1,:)      = tran_med_index_sbp1(1,:)
            tran_med_sbp1(np_med,:) = tran_med_index_sbp1(np_med,:)
         ENDIF
 
      DO j = 2, np_med-1
         IF(wavln_med_index(j).LE.wavln_med(j)) THEN
           tran_med(j,:) = tran_med_index(j,:)
           IF(isplitp.NE.0) THEN
                tran_med_sa(j,:)   = tran_med_index_sa(j,:)
                tran_med_sap1(j,:) = tran_med_index_sap1(j,:)
                tran_med_sb(j,:)   = tran_med_index_sb(j,:)
                tran_med_sbp1(j,:) = tran_med_index_sbp1(j,:)
           ENDIF
         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,:)
           IF(isplitp.NE.0) THEN
                tran_med_sa(j,:)   = fjm1*tran_med_index_sa(j-1,:)   + fj*tran_med_index_sa(j,:)
                tran_med_sap1(j,:) = fjm1*tran_med_index_sap1(j-1,:) + fj*tran_med_index_sap1(j,:)
                tran_med_sb(j,:)   = fjm1*tran_med_index_sb(j-1,:)   + fj*tran_med_index_sb(j,:)
                tran_med_sbp1(j,:) = fjm1*tran_med_index_sbp1(j-1,:) + fj*tran_med_index_sbp1(j,:)
!                write(91,*) J,' TRAN_MED_SAP1=',TRAN_MED_INDEX(J,JA_L2),FJM1,FJ,TRAN_MED_INDEX(J-1,JA_L2), &
!                TRAN_MED_INDEX_SAP1(J,JA_L2)
!                write(91,*) J,' TRAN_MED_SAP1_SOL=',TRAN_MED_SAP1(J,1),FJM1,FJ,TRAN_MED_INDEX_SAP1(J-1,1), &
!                TRAN_MED_INDEX_SAP1(J,1)
           ENDIF
!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 = npshif+1, np_std
          tran_std(i,:) = tran_std(i,:)*tran_med(i-npshif,:)
       END DO
 
       IF(isplitp.NE.0) THEN
          DO i = 1, np_std
           tran_std_sa(i,:)   = 1.
           tran_std_sap1(i,:) = 1.
           tran_std_sb(i,:)   = 1.
           tran_std_sbp1(i,:) = 1.
          END DO
          DO i = npshif+1, np_std
            tran_std_sa(i,:)   = tran_std_sa(i,:)  *tran_med_sa(i-npshif,:)
            tran_std_sap1(i,:) = tran_std_sap1(i,:)*tran_med_sap1(i-npshif,:)
            tran_std_sb(i,:)   = tran_std_sb(i,:)  *tran_med_sb(i-npshif,:)
            tran_std_sbp1(i,:) = tran_std_sbp1(i,:)*tran_med_sbp1(i-npshif,:)
          END DO
!          DO I = 1,NP_STD
!           write(91,*) I,' TRAN_STD_SA_SOL=',TRAN_STD_SA(I,1)
!            write(91,*) I,' TRAN_STD_SAP1_SOL=',TRAN_STD_SAP1(I,1)
!          END DO
       ENDIF
 
 
 
!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:
!C
      DO 1466 j =1, nobs
 
             trntbl(j,:)       = 0.0
             tran_ia(:)        = 0.0
             tran_iap1(:)      = 0.0
 
       IF(isplitp.NE.0) THEN
             trntbl_s(j,:)     = 0.0
             tran_ia_sa(:)     = 0.0
             tran_iap1_sa(:)   = 0.0
             tran_ia_sap1(:)   = 0.0
             tran_iap1_sap1(:) = 0.0
             tran_ia_sb(:)     = 0.0
             tran_iap1_sb(:)   = 0.0
             tran_ia_sbp1(:)   = 0.0
             tran_iap1_sbp1(:) = 0.0
       ENDIF
!C---
!C---        print*,'J= ',j, 'NCVHF =', NCVHF(J), 'NCVTOT=',NCVTOT
    
!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
!C---          print*,'J =', J, ' IA =', IA
 
!
!C  Smoothing...
!C
          DO 1491 k = ia(j)-(ncvhf(j)-1), ia(j)+ncvhf(j)-1
             tran_ia(:) = tran_ia(:) + tran_std(k,:)* &
                          finstr(k-ia(j)+ncvhf(j),j)
             IF(isplitp.NE.0) THEN
                 tran_ia_sa(:) = tran_ia_sa(:)     + tran_std_sa(k,:)*   &
                          finstr(k-ia(j)+ncvhf(j),j)
                 tran_ia_sap1(:) = tran_ia_sap1(:) + tran_std_sap1(k,:)* &
                          finstr(k-ia(j)+ncvhf(j),j)
                 tran_ia_sb(:) = tran_ia_sb(:)     + tran_std_sb(k,:)*   &
                          finstr(k-ia(j)+ncvhf(j),j)
                 tran_ia_sbp1(:) = tran_ia_sbp1(:) + tran_std_sbp1(k,:)* &
                          finstr(k-ia(j)+ncvhf(j),j)
!                 write(91,*) K,' TRAN_IA=',TRAN_IA(JA_L2),TRAN_STD(K,JA_L2),FINSTR(K-IA(J)+NCVHF(J),J)
!                 write(91,*) K,' TRAN_IA_SA_SOL=',TRAN_IA_SA(1),TRAN_STD_SA(K,1),FINSTR(K-IA(J)+NCVHF(J),J)
!                 write(91,*) K,' TRAN_IA_SAP1_SOL=',TRAN_IA_SAP1(1),TRAN_STD_SAP1(K,1),FINSTR(K-IA(J)+NCVHF(J),J)
             ENDIF
 
!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(j) + 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),j)
             IF(isplitp.NE.0) THEN
                tran_iap1_sa(:) = tran_iap1_sa(:)     + tran_std_sa(k,:)*   &
                          finstr(k-ia_p1+ncvhf(j),j)
                tran_iap1_sap1(:) = tran_iap1_sap1(:) + tran_std_sap1(k,:)* &
                          finstr(k-ia_p1+ncvhf(j),j)
                tran_iap1_sb(:) = tran_iap1_sb(:)     + tran_std_sb(k,:)*   &
                          finstr(k-ia_p1+ncvhf(j),j)
                tran_iap1_sbp1(:) = tran_iap1_sbp1(:) + tran_std_sbp1(k,:)* &
                          finstr(k-ia_p1+ncvhf(j),j)
             ENDIF
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(j))
           fia     = (wavln_std(ia_p1) - wavobs(j)) /dlt_ia
!C          FIA_P1  = (WAVOBS(J)     - WAVLN_STD(IA))/DLT_IA
           fia_p1  = 1. - fia
           trntbl(j,:) = fia*tran_ia(:) + fia_p1*tran_iap1(:)
           IF(isplitp.NE.0) THEN
            trntblsa(:)   = fia*tran_ia_sa(:)   + fia_p1*tran_iap1_sa(:)
            trntblsap1(:) = fia*tran_ia_sap1(:) + fia_p1*tran_iap1_sap1(:)
            trntblsb(:)   = fia*tran_ia_sb(:)   + fia_p1*tran_iap1_sb(:)
            trntblsbp1(:) = fia*tran_ia_sbp1(:) + fia_p1*tran_iap1_sbp1(:)
            trntbl_s(j,:) = 0.5*(f1a*trntblsa(:) + f2a*trntblsap1(:) + &
                                 f1b*trntblsb(:) + f2b*trntblsbp1(:) )
!            write(91,*) J,' TRNTBL_TOT=',TRNTBL(J,JA_L2),FIA,FIA_P1,TRAN_IA(JA_L2),TRAN_IAP1(JA_L2)
!            write(91,*) J,' TRNTBL_SOL_SA=',TRNTBLSA(1),TRNTBLSAP1(1),F1A,F2A,F1B,F2B,TRAN_IA_SA(1), &
!                             TRAN_IAP1_SA(1),TRAN_IA_SAP1(1),TRAN_IAP1_SAP1(1)
!            write(91,*) J,JA_L2,JB_L2,' TRNTBL_SOL=',TRNTBL_S(J,1),TRNTBL(J,1),F1A,F2A,F1B,F2B,TRNTBLSA(1),TRNTBLSAP1(1), &
!                             TRNTBLSB(1),TRNTBLSBP1(1)
!            write(91,*) J,JA_L2,JB_L2,' TRNTBL_SEN=',TRNTBL_S(J,2),TRNTBL(J,2),F1A,F2A,F1B,F2B,TRNTBLSA(2),TRNTBLSAP1(2), &
!                             TRNTBLSB(2),TRNTBLSBP1(2)
           ENDIF
           IF (ifullcalc.eq.0) THEN
           do k=1,nh2o_max
             !! TRAN_KD(J,K) = TRAN_KD(J,K)*DIFF_TRAN(J,K)
             diff_tran(j,k) =  (trntbl(j,k)-tran_kd(j,k))/tran_kd(j,k) + 1
             !!DTT =  (TRAN_KD(J,K)-TRNTBL(J,K))/TRNTBL(J,K)
             !write(*,*) 'RJH: TRAN2: ',J, K, TRAN_KD(J,K),TRNTBL(J,K),DIFF_TRAN(J,K) !!,100*DTT
           end do
 
           end if
!C---
!C---       print*,'j=',j,'IA =',IA,'FIA =',FIA,'FIA_P1=',FIA_P1,DLT_IA
!C---
!C
1466  CONTINUE
 
 
      DEALLOCATE(tran_std,tran_med)
      RETURN
      END
 
!********************************************************************************
!*                                                                              *
!*  Name: TRAN_SMOOTH_OTHERS                                                    *
!*  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.                                                        *
!*                                                                              *
!* Added by R. Healy (richard.healy@nasa.gov) 9/29/2015                         *
!*                                                                              *
!* This version smooths the other gases, not H2O                                *
!********************************************************************************
 
      SUBROUTINE tran_smooth_others
 
      include 'COMMONS_INC.f'
 
      dimension tran_o3_std(no3pt)
      COMMON /init_speccal16/ tran_o3_std
 
      dimension tran_no2_std(no3pt)
      COMMON /init_speccal17/ tran_no2_std
 
      dimension wavobs(nobs_max),fwhm(nobs_max)
      COMMON /getinput4/ wavobs,fwhm
 
      dimension vaptot(nh2o_max), r0p94(nh2o_max), r1p14(nh2o_max), trntbl(nobs_max,nh2o_max), &
                 tran_kd(nobs_max,nh2o_max), diff_tran(nobs_max,nh2o_max),trntblo(nobs_max)
      COMMON /tran_table1/ sh2o,vaptot,r0p94,r1p14,trntbl,tran_kd, diff_tran,trntblo
      REAL   TRAN_HI_SA(NP_HI,2),TRAN_HI_SAP1(NP_HI,2),TRAN_HI_SB(NP_HI,2),TRAN_HI_SBP1(NP_HI,2), &
              trntbl_s(nobs_max,2),trntbl_so(nobs_max,2)
      COMMON /tran_tables/tran_hi_sa,tran_hi_sap1,tran_hi_sb,tran_hi_sbp1
      COMMON /tran_table_l2gen/trntbl_s,trntbl_so
      REAL   MU,MU0, SSH2O_S(NH2O_MAX,2)
      COMMON /geometry5/mu,mu0,ssh2o_s
 
      real*4 tran_ia,tran_iap1, tran_std_o(np_std),tran_std_os(np_std,2),am(np_std),tran_ias(2),tran_iap1s(2)
 
      COMMON /getinput5/ nobs,ifullcalc,hsurf,dlt,dlt2
      REAL TRAN_MED_O(NP_MED),TRAN_MED_INDEX_O(NP_MED)        !Transmittance of medium resolution data
      REAL TRAN_MED_OS(NP_MED,2),TRAN_MED_INDEX_OS(NP_MED,2)
      INTEGER FIRST/1/,IA(NOBS_MAX)
      SAVE first, ia
 
!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
    IF (first.eq.1) THEN !RJH
      DO 466 j =1, np_med
 
             ncvtot_wavno = 2 * ncvhf_wavno(j) - 1
 
             sumins = 0.0
 
          DO 560 i = ncvhf_wavno(j), ncvtot_wavno
             finstr_wavno(i,j) = &
                exp( -const1*(float(i-ncvhf_wavno(j))*dwavno &
                /fwhm_wavno(j))**2)
             sumins = sumins + finstr_wavno(i,j)
 560      CONTINUE
 
          DO 565 i = 1, ncvhf_wavno(j)-1
             finstr_wavno(i,j) = finstr_wavno(ncvtot_wavno-i+1,j)
             sumins = sumins + finstr_wavno(i,j)
 565      CONTINUE
 
          sumins = sumins * dwavno
 
          DO 570 i = 1, ncvtot_wavno
             finstr_wavno(i,j) = finstr_wavno(i,j)*dwavno/sumins
 570      CONTINUE
 
 
 466  CONTINUE
 
 
         tran_med_o(:)  = 1.0                         ! FWHM=.2 nm, .56-3.1 um.
 
 
!C  Index searching...
!C
       DO j=1,nobs
         !C Calculate instrumental response functions...
 
          sumins = 0.0
 
          ncvtot = 2 * ncvhf(j) - 1
 
          DO 1560 i = ncvhf(j), ncvtot
             finstr(i,j) = &
                exp( -const1*(float(i-ncvhf(j))*dwavln &
                /fwhm(j))**2)
             sumins = sumins + finstr(i,j)
1560      CONTINUE
 
          DO 1565 i = 1, ncvhf(j)-1
             finstr(i,j) = finstr(ncvtot-i+1,j)
             sumins = sumins + finstr(i,j)
1565      CONTINUE
 
          sumins = sumins * dwavln
 
          DO 1570 i = 1, ncvtot
             finstr(i,j) = finstr(i,j)*dwavln/sumins
1570      CONTINUE
 
 
          CALL hunt(wavln_std, np_std, wavobs(j), ia(j))
       END DO
 
      END IF ! IFIRST = 0
 
     airmass = 1.0/mu0 + 1.0/mu;
 
!C*** !!!High resolution transmittances for CO2, N2O, CO, CH4, and O2 should
!C         also be calculated for wavelength start = 0.56 micron.
 
     tran_med_index_o(:)  = 0.0
     IF(isplitp.NE.0) tran_med_index_os(:,:)  = 0.0
 
     DO  j =1, np_med
 
          ndx1 = index_med(j)-(ncvhf_wavno(j)-1)
          ndx2 = index_med(j)+ ncvhf_wavno(j)-1
          DO 491 k = ndx1,ndx2
             tran_med_index_o(j) = tran_med_index_o(j) + tran_hi_others(k)*  &
                          finstr_wavno(k-index_med(j)+ncvhf_wavno(j),j)
 
              IF(isplitp.NE.0) THEN
                airm = -log(tran_hi_others(k))/airmass
                tran_hi_others_sol = exp(-airm/mu0)
                tran_hi_others_sen = exp(-airm/mu)
 
                tran_med_index_os(j,1) = tran_med_index_os(j,1) + tran_hi_others_sol*  &
                          finstr_wavno(k-index_med(j)+ncvhf_wavno(j),j)
                tran_med_index_os(j,2) = tran_med_index_os(j,2) + tran_hi_others_sen*  &
                          finstr_wavno(k-index_med(j)+ncvhf_wavno(j),j)
 
              ENDIF
!C            if (FIRST.eq.1) print*,j," tran_med_index:", TRAN_MED_INDEX_O(J)
 491      CONTINUE
     END DO
 
!C
!C Linear interpolation to get TRAN_MED from TRAN_MED_INDEX_O:
!C     (Note that WAVLN_MED_INDEX(J) >= WAVLN_MED(J)    )
!C
         tran_med_o(1)      = tran_med_index_o(1)
         tran_med_o(np_med) = tran_med_index_o(np_med)
         IF(isplitp.NE.0) THEN
            tran_med_os(1,:)      = tran_med_index_os(1,:)
            tran_med_os(np_med,:) = tran_med_index_os(np_med,:)
         END IF
 
      DO j = 2, np_med-1
         IF(wavln_med_index(j).LE.wavln_med(j)) THEN
           tran_med_o(j) = tran_med_index_o(j)
           IF(isplitp.NE.0) tran_med_os(j,:) = tran_med_index_os(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_o(j) = fjm1*tran_med_index_o(j-1) + fj*tran_med_index_o(j)
           IF(isplitp.NE.0) &
             tran_med_os(j,:) = fjm1*tran_med_index_os(j-1,:) + fj*tran_med_index_os(j,:)
 
         END IF
      END DO
 
!C
!C--- Here multiplying O3 and NO2 spectra and other spectrum at medium resolution:
!C
       tran_std_o(:) = 1.
       IF(isplitp.NE.0) tran_std_os(:,:) = 1.
 
       tran_std_o(1:no3pt) = tran_o3_std(1:no3pt) * tran_no2_std(1:no3pt)
       IF(isplitp.NE.0) THEN
           am(1:no3pt) = -log(tran_std_o(1:no3pt))/airmass
           tran_std_os(1:no3pt,1) = exp(-am(1:no3pt)/mu0)
           tran_std_os(1:no3pt,2) = exp(-am(1:no3pt)/mu)
       END IF
 
       tran_std_o(npshif+1:np_std) = tran_std_o(npshif+1:np_std)*tran_med_o(1:np_std-npshif)
       IF(isplitp.NE.0) THEN
           tran_std_os(npshif+1:np_std,:) = tran_std_os(npshif+1:np_std,:)*tran_med_os(1:np_std-npshif,:)
!       do j=NPSHIF+1,NP_STD
!        write(91,*) j, ' TRAN_STD_OS_SOL=',TRAN_STD_OS(j,1),TRAN_MED_OS(j-NPSHIF,1),TRAN_STD_O(j),TRAN_MED_O(J)
!       write(91,*) j, ' TRAN_STD_OS_SEN=',TRAN_STD_OS(j,2),TRAN_MED_OS(j-NPSHIF,2),TRAN_STD_O(j),TRAN_MED_O(J)
!       end do
       ENDIF
!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:
!C
      DO 1466 j =1, nobs
 
             tran_ia         = 0.0
             tran_iap1       = 0.0
             tran_ias(:)     = 0.0
             tran_iap1s(:)   = 0.0
 
!
!C  Smoothing...
!C
          DO 1491 k = ia(j)-(ncvhf(j)-1), ia(j)+ncvhf(j)-1
             tran_ia = tran_ia + tran_std_o(k)* &
                          finstr(k-ia(j)+ncvhf(j),j)
            IF(isplitp.NE.0) &
                tran_ias(:) = tran_ias(:) + tran_std_os(k,:)* &
                          finstr(k-ia(j)+ncvhf(j),j)
 
 
1491      CONTINUE
 
            ia_p1 = ia(j) + 1
          DO 1492 k = ia_p1-(ncvhf(j)-1), ia_p1+ncvhf(j)-1
             tran_iap1 = tran_iap1 + tran_std_o(k)* &
                          finstr(k-ia_p1+ncvhf(j),j)
           IF(isplitp.NE.0) &
             tran_iap1s(:) = tran_iap1s(:) + tran_std_os(k,:)* &
                          finstr(k-ia_p1+ncvhf(j),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(j))
           fia     = (wavln_std(ia_p1) - wavobs(j)) /dlt_ia
!C          FIA_P1  = (WAVOBS(J)     - WAVLN_STD(IA))/DLT_IA
           fia_p1  = 1. - fia
           trntblo(j) = fia*tran_ia + fia_p1*tran_iap1
           IF(isplitp.NE.0) THEN
            trntbl_so(j,:) = fia*tran_ias(:) + fia_p1*tran_iap1s(:)
!            write(91,*) 'TRNTBL_SO:',TRNTBL_SO(J,1),TRNTBL_SO(J,2),TRAN_IAS(1), &
!             TRAN_IAS(2),TRAN_IAP1S(1),TRAN_IAP1S(2)
           ENDIF
 
!C---
!C---       print*,'j=',j,'IA =',IA,'FIA =',FIA,'FIA_P1=',FIA_P1,DLT_IA
!C---
!C
1466  CONTINUE
 
 
      first = 0
 
      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
      parameter(nobs_max=1024,nh2o_max=60)
 
!C  Common variables
      dimension const1(nh2o_max),const2(nh2o_max),const3(nh2o_max)
      dimension const4(nh2o_max),const5(nh2o_max),const6(nh2o_max)
 
      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
      dimension vaptot(nh2o_max), r0p94(nh2o_max), r1p14(nh2o_max),trntbl(nobs_max,nh2o_max), &
                 tran_kd(nobs_max,nh2o_max), diff_tran(nobs_max,nh2o_max),trntblo(nobs_max)
      COMMON /tran_table1/ sh2o,vaptot,r0p94,r1p14,trntbl, tran_kd, diff_tran,trntblo
 
!C Calculate average of spectra over window and water vapor absorption regions.
      const1(:)=0.0
      DO 560 i=ist1,ied1
        const1(:)=const1(:)+trntbl(i,:)
  560 CONTINUE
      const1(:)=const1(:)/float(nb1)
 
      const2=0.0
      DO 570 i=ist2,ied2
        const2(:)=const2(:)+trntbl(i,:)
  570 CONTINUE
      const2(:)=const2(:)/float(nb2)
 
      const3(:)=0.0
      DO 575 i=istp94,iedp94
        const3(:)=const3(:)+trntbl(i,:)
  575 CONTINUE
      const3(:)=const3(:)/float(nbp94)
 
      r0p94(:)=const3(:)/((wt1*const1(:)) + (wt2*const2(:)))
 
      const4(:)=0.0
      DO 580 i=ist3,ied3 
        const4(:)=const4(:)+trntbl(i,:)
  580 CONTINUE
      const4(:)=const4(:)/float(nb3)
 
      const5(:)=0.0
      DO 590 i=ist4,ied4
        const5(:)=const5(:)+trntbl(i,:)
  590 CONTINUE
      const5(:)=const5(:)/float(nb4)
 
      const6(:)=0.0
      DO 595 i=ist1p14,ied1p14
        const6(:)=const6(:)+trntbl(i,:)
  595 CONTINUE
      const6(:)=const6(:)/float(nb1p14)
 
      r1p14(:)=const6(:)/((wt3*const4(:)) + (wt4*const5(:)))
 
      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)
      parameter(nobs_max=1024)
 
      dimension xorgn(1050),yorgn(1050),y2(1050)
      dimension xint(nobs_max),yint(nobs_max)
      INTEGER N                               !number of elements in XORGN,YORGN
 
      COMMON /getinput5/ nobs,ifullcalc,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: 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                               *
!*                                         *
!********************************************************************************
 
      INTEGER FUNCTION findmatch(LIST,NOBS,ELEM)
      parameter(nobs_max=1024)
      dimension list(nobs_max)
      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--