OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
atrem_app_refl_plus_gas_removal_2013.f
Go to the documentation of this file.
1 !------------------------------------------------------------------------------|
2 ! The September 2012 version of atrem code, 'atrem_for_PRISM.f' (atrem for JPL |
3 ! PRISM Imaging Spectrometer, was obtained through modifying a 2011 |
4 ! version of 'atrem_f90_cubeio.f'. |
5 ! The March 2013 version of atrem code, 'atrem_for_PRISM_2013.f', was updated |
6 ! from the September 2012 version of atrem code. The upgrades include: |
7 ! a) replacing the MODTRAN 3.5 solar irradiance curve with a new solar curve|
8 ! constructed from the solar irradiance curve of Thuillier et al. 2004 |
9 ! ATLAS 3 (<=644.7 nm) and that of Kurucz 2005 (> 644.7 nm) built in |
10 ! MODTRAN 5.2. This newly contructed solar curve is more suited for |
11 ! modeling imaging imaging spectrometer data below 450 nm and with a |
12 ! spectral resolution finer than 5 nm; |
13 ! b) replacing high spectral resolution gas absorption tables with new |
14 ! tables calculated with the HITRAN2008+ line database and a fast |
15 ! line-by-line code originally developed by Bill Ridgway. After such |
16 ! upgrades, the retrieved surface reflectances near 1.45 and 1.95 micron |
17 ! are improved significantly; |
18 ! c) modifying volume mixing ratios of atmospheric carbon dioxide; |
19 ! d) applying a scaling factor for improved modeling of atmospheric |
20 ! oxygen band absorption centered near 1.265 micron. |
21 !------------------------------------------------------------------------------|
22 ! The 2001 version of ATREM code uses the new 'cubeio.f90' for I/O operations. |
23 ! The output data cubes have no headers, unlike the previous generations |
24 ! of ATREM codes containing 512 bytes of headers in output data cubes. |
25 !------------------------------------------------------------------------------|
26 ********************************************************************************
27 * *
28 * Name: atrem_app_refl_plus_gas_removal_2013 April, 2013 *
29 * (ATmosphere REMoval Program) *
30 * *
31 * Author: Bo-Cai Gao, Remote Sensing Division, Code 7230, *
32 * Naval Research Laboratory, Washington, DC 20375 USA *
33 * *
34 * Purpose: To derive surface reflectances from spectral imaging data collected *
35 * by the Airborne Visible Infrared Imaging Spectrometer (AVIRIS), *
36 * HYDICE, HSI, PHILLS, HICO, NIS, PRISM, and possibly other *
37 * airborne and spaceborne imaging spectrometers, and to derive a *
38 * column water vapor image for each scene. *
39 * *
40 * General Principles: In order to derive surface reflectances from image data *
41 * cube, a thorough compensation for the atmospheric absorption and *
42 * scattering is required. The spatial and temporal variations of *
43 * atmospheric water vapor amounts pose difficulties in removing *
44 * water vapor absorption features in data cube using standard *
45 * atmospheric models. In this algorithm, the amount of water vapor *
46 * on a pixel-by-pixel basis is derived from data cube themselves *
47 * using the 0.94- and the 1.14-um water vapor bands and a *
48 * three-channel ratioing technique. The derived water vapor *
49 * values are then used for modeling water vapor absorption effects *
50 * in the entire 0.4-2.5 um region. The absorption effects of *
51 * well mixed atmospheric gases, such as CO2, N2O, CH4, and O2, *
52 * are modeled using standard atmospheric models. A line-by-line *
53 * code developed by W. Ridgway at NASA Goddard Space Flight Center *
54 * is used to generate a database of gaseous absorption *
55 * coefficients at high spectral resolution. This database and *
56 * ozone cross sections and NO2 cross sections are used in *
57 * calculating transmittances of eight atmospheric gases (H2O, CO2 *
58 * O3, N2O, CO, CH4, O2, and NO2). The scattering effect is modeled *
59 * using a modified version of 6S code (Vermote et al. 1996). *
60 * *
61 * Algorithm: 1. An input parameter file is read in and global data are *
62 * initialized. *
63 * 2. The solar zenith angle is derived based on the flight time *
64 * and latitude and longitude of the scene. *
65 * 3. A table consisting of transmittance spectra at the solar *
66 * and observational geometry but with varying amounts of water *
67 * vapor is created. *
68 * The transmittance spectra for water vapor (H2O), carbon *
69 * dioxide (CO2), ozone (O3), nitrous oxide (N2O), carbon *
70 * monoxide (CO), methane (CH4), oxygen (O2), and nitrogen *
71 * dioxide (NO2) in the 0.4-2.5 micron region are calculated. *
72 * 4. Aerosol and molecular scattering terms are calculated using *
73 * a modified version of the 6S code. *
74 * 5. The image data cube (2D spatial and 1D spectral) are accessed *
75 * one spectral slice at a time from an image file in any storage*
76 * order, Band SeQuential (BSQ), Band Interleaved by Pixel (BIP),*
77 * and Band Interleaved by Line (BIL). Each measured radiance *
78 * spectrum is divided by the solar irradiance curve above the *
79 * atmosphere to get apparent reflectances. *
80 * 6. The integrated water vapor amount for each apparent *
81 * reflectance spectrum is derived from the 0.94- and the 1.14-um*
82 * water vapor absorption bands using a 3-channel ratioing *
83 * technique and a look-up table procedure. *
84 * 7. The surface reflectances are derived using another look-up *
85 * table procedure. Routines written in C language are used to *
86 * perform all the input/output operations. *
87 * *
88 * Limitations: *
89 * 1. This algorithm works best for areas where the surface elevation *
90 * does not vary by more than 1 km. *
91 * 2. The algorithm works well only for measurements taken on clear *
92 * days. For hazy days with large aerosol optical depths, the *
93 * algorithm does not work nearly as well. *
94 * 3. The algorithm does not work well over dark areas such as rivers *
95 * and lakes. It does not perform corrections for "atmospheric *
96 * adjacency" effects. *
97 * 4. The algorithm assumes horizontal surfaces having Lambertian *
98 * reflectances. *
99 * 5. Additional work is needed to port this algorithm to different *
100 * computing systems. This algorithm was developed on an SGI *
101 * workstation and used FORTRAN statements to read several binary *
102 * files. The binary file format on different computers can be *
103 * different. The record lengths defined on an SGI machine are *
104 * different from most other computers. Some of the Dec Alpha *
105 * workstations use 8 bytes (64 bits) to store a C-pointer, instead *
106 * of 4 bytes (32 bits). All these factors need to be taken into *
107 * considerations when porting this algorithm to other computers. *
108 * *
109 * User Input: *
110 * 1. Name - Imaging spectrometer name, e.g., AVIRIS, HYDICE, HSI, PHYLLS. *
111 * 2. Plane altitude in km, above the sea level. *
112 * 3. Date - The month, day, and year that the data were collected. *
113 * 4. Time - The hour, minute, second (GMT) that the data were *
114 * collected. *
115 * 5. Latitude - The mean latitude of the measurement area in degrees, *
116 * minutes, seconds, and hemisphere. *
117 * 6. Longitude - The mean longitude of the measurement area in *
118 * degrees, minutes, seconds, and hemisphere. *
119 * 7. View zenith angle - the angle in degrees, minutes, and seconds *
120 * between the viewing direction and nadir *
121 * 8. View azimuth angle in degree, minutes and seconds. The view azimuth *
122 * angle is defined as the angle between local North *
123 * the view direction projected on the ground *
124 * (clockwise count). This angle can vary from 0 *
125 * to 360 deg. Ths definition here is consistent with*
126 * navigator's convention, but not consistent with *
127 * astronomer's convention. Here is an example: for *
128 * a downward looking instrument pointed off-nadir *
129 * to the West, the view azimuthal angle is 90 deg. *
130 * Conversely, if the instrument is pointed toward *
131 * East, the view azimuthal angle is 270 deg. *
132 * *
133 * 9. Resolution of input spectra in nm (~ 10 nm for AVIRIS) or "0" *
134 * if the FWHM values (in micron) are present in the wavelength file *
135 * 10. Filename of spectrometer wavelength table - 2 or 3 column ASCII data. *
136 * The first column is the channel number, the second column is *
137 * the wavelength and the optional third column is the FWHM value for *
138 * each channel. *
139 * 11. Channel ratio parameters - Center positions and widths of window and *
140 * water vapor channels used in channel ratios for deriving *
141 * the amount of water vapor from imaging data. Each of the window *
142 * and water vapor channels typically consists of several narrow *
143 * imaging spectrometer channels. *
144 * 12. Atmospheric model - Temperature, pressure, water vapor volume mixing *
145 * ratio profiles. The model can be either a predefined model or a *
146 * user defined model. *
147 * 13. Gases - An indication of which of the 8 gases (H2O, CO2, O3, *
148 * N2O, CO, CH4, O2, and NO2) should be included in the spectral *
149 * calculations. *
150 * 14. Ozone - Column amount of O3 - changes with season and *
151 * latitude. Typical value is .34 atm-cm. *
152 * 15. Aerosol model number and visibility when measurements were taken. *
153 * 16. Aerosol optical depth at 550 nm (Optional, only if visibility V *
154 * is set to 0) *
155 * 17. Surface elevation - The average surface elevation of an imaging scene *
156 * in km. *
157 * 18. Filename of input image file *
158 * 19. Dimensions of input image. *
159 * 20. Filename of output image file in same storage order as input file *
160 * 21. Resolution of output surface reflectance spectra in micron *
161 * 22. Scale factor for output reflectance values *
162 * 23. Filename of output water vapor image. *
163 * 24. Filename of output atmospheric transmittance look-up table *
164 * *
165 * Output: *
166 * a. Output to user specified files: *
167 * 1. Surface reflectance cube data - retrieved from the image data cube *
168 * 2. Water vapor image - an image of the derived column water vapor *
169 * amounts at each pixel. *
170 * 3. Transmittance Look-up table - consisting of channel ratio values *
171 * for the .94- and 1.14-um channels corresponding to 60 water vapor *
172 * amounts, and the associated atmospheric transmittance spectra. *
173 * b. Output to standard output: *
174 * 1. Debugging information - if the source is compiled with the *
175 * debug flag set (-d_lines for the ULTRIX compiler), *
176 * then debug information is written out. *
177 * 2. Error messages - if any of the user input is invalid, or there is *
178 * an I/O error, a message is written out, and the program halts. *
179 * 3. Progress indicator - a message, which shows the number of the *
180 * spectral slices that have been processed, to give the user an *
181 * indication of the progress of the program execution. *
182 * *
183 * Special Considerations: *
184 * a. Make sure that the imaging spectrometer's channel positions specified *
185 * in the input wavelength table are correct. The incorrect channel *
186 * positions will introduce sharp features in derived surface *
187 * reflectance spectra. *
188 * b. The output reflectance cube file has the same size as the input *
189 * data cube file. Make sure there is enough space in the file *
190 * system for the output cube before running this program. *
191 * *
192 * Change History: *
193 * ATREM 1.2 - uses full-width half-max values to smooth solar irradiance *
194 * curve and a new scaling factor for methane *
195 * ATREM 1.3 - used new solar irradiance curve *
196 * - supports 1992 AVIRIS data *
197 * ATREM 1.3.1 - allows 0 and above for output resolution *
198 * ATREM 2.0 - allows variable scale factors for each spectrometer *
199 * (code submitted by Roger Clark, USGS) *
200 * - user can input output scale factor for reflectance values *
201 * - user can input radiance file header size in bytes *
202 * - new solar irradiance curve (Neckel and Labs plus ATMOS) *
203 * used (Green and Gao, 1993). *
204 * ATREM 3.0 - replaced the 5S code by the 6S code for modeling atmospheric *
205 * scattering effects and for modeling measurements from *
206 * low-altitude aircrafts. ATREM 3.0, like previous versions *
207 * of ATREM, only models nadir viewing geometry. *
208 * - changed algorithms for calculating atmospheric gaseous *
209 * transmittances for the upward surface-sensor path *
210 * - users need to specify aircraft altitude (km) in input file *
211 * ATREM 4.0 - replaced the band model by a line-by-line-based algorithm *
212 * for calculating atmospheric gaseous transmittances. This *
213 * allows the modeling of imaging spectrometer data at *
214 * spectral resolutions better than 10 nm. *
215 * - increased the buffer size to 1024x1024 in order to handle *
216 * images larger than the AVIRIS' size (614x512). *
217 * - modified the algorithm to allow off-nadir pointing geometry. *
218 * - users need to specify view zenith angle and view azimuth *
219 * angle. The view azimuth angle is defined as the angle *
220 * between local North and the view direction projected on the *
221 * ground (clockwise count). This angle can vary from 0 to *
222 * 360 deg. *
223 * - replaced the solar irradiance curve (Neckel and Labs plus *
224 * ATMOS) by the high spectral resolution solar irradiance *
225 * curve from MODTRAN 3.5 released in December of 1996. *
226 * ATREM 4.1 - included atmospheric NO2 in transmittance calculations. *
227 * *
228 * Acknowledgments: *
229 * This work was partially supported over several years by grants *
230 * from NASA Jet Propulsion Laboratory, California Institute of *
231 * Technology, from NASA Headquarters, and from the Office of Naval *
232 * Research. Special thanks goes to Kathy Heidebrecht at the Center *
233 * for the Study of Earth from Space, University of Colorado at *
234 * Boulder, Colorado, for supporting the development of the earlier *
235 * versions of ATREM codes. Special thanks also goes to William L. *
236 * Ridgway for providing the line-by-line gaseous absorption database *
237 * used in the present version of ATREM. *
238 * *
239 * References: *
240 * Gao, B.-C., K. Heidebrecht, and A. F. H. Goetz, Derivation of scaled *
241 * surface reflectances from AVIRIS data, Remote Sens. Env., 44, *
242 * 165-178, 1993. *
243 * Gao, B.-C., and C. O. Davis, Development of an operational algorithm *
244 * for removing atmospheric effects from HYDICE and HSI data, *
245 * in SPIE'96 Conference Proceedings, Vol. 2819, 45-55, 1996. *
246 * Gao, B.-C., and A. F. H. Goetz, Column atmospheric water vapor and *
247 * vegetation liquid water retrievals from airborne imaging *
248 * spectrometer data, J. Geophys. Res., 95, 3549-3564, 1990. *
249 * Goetz, A. F. H., and M. Herring, The high resolution imaging *
250 * spectrometer (HIRIS) for Eos, IEEE Trans. Geosci. Remote Sens., 27,*
251 * 136-144, 1989. *
252 * Goetz, A. F. H., G. Vane, J. Solomon, and B. N. Rock, Imaging *
253 * spectrometry for Earth remote sensing, Science, 228, 1147-1153,1985*
254 * Green, R. O., and B.-C. Gao, A Proposed Update to the Solar Irradiance*
255 * Spectrum used in LOWTRAN and MODTRAN, in Summaries of the Fourth *
256 * Annual JPL Airborne Geoscience Workshop, October 25-29, (Editor, *
257 * R. O. Green), JPL Publ. 93-26, Vol. 1, pp. 81-84, Jet Propul. Lab, *
258 * Pasadena, Calif., 1993. *
259 * Kneizys, F. X., E. P. Shettle, L. W. Abreu, J. H. Chetwynd, G. P. *
260 * Anderson, W. O. Gallery, J. E. A. Selby, and S. A. Clough, Users *
261 * guide to LOWTRAN7, AFGL-TR-8-0177, Air Force Geophys. Lab., *
262 * Bedford, Mass., 1988. *
263 * Iqbal, M., An Introduction To Solar Radiation, Academic, San Diego, *
264 * Calif., 1983. *
265 * Malkmus, W., Random Lorentz band model with exponential-tailed S line *
266 * intensity distribution function, J. Opt. Soc. Am., 57, 323-329,1967*
267 * Press, W. H., B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, *
268 * Numerical Recipes-The ART of Scientific Computing, Cambridge *
269 * University Press, 1986. *
270 * Rothman, L. S., et al., The HITRAN 2008 molecular spectroscopic *
271 * database, JQSRT, 110, 533-572, 2009. *
272 * Solomon, S., R. W. Portmann, R. W. Sanders, J. S. Daniel, W. Madsen, *
273 * B. Bartram, and E. G. Dutton, On the role of nitrogen dioxide in *
274 * the absorption of solar radiation, J. Geophys. Res., 104, *
275 * 12047-12058, 1999. *
276 * Tanre, D., C. Deroo, P. Duhaut, M. Herman, J. J. Morcrette, J. Perbos,*
277 * and P. Y. Deschamps, Description of a computer code to simulate *
278 * the satellite signal in the solar spectrum: the 5S code, Int. *
279 * J. Remote Sens., 11, 659-668, 1990. *
280 * Tanre, D., C. Deroo, P. Duhaut, M. Herman, J. J. Morcrette, J. Perbos,*
281 * and P. Y. Deschamps, Simulation of the satellite signal in the *
282 * solar spectrum (5S), Users' Guide (U. S. T. De Lille, 59655 *
283 * Villeneu d'Ascq, France: Laboratoire d'Optique Atmospherique), *
284 * 1986. *
285 * Thuillier, G., et al., Solar irradiance reference spectra for two *
286 * solar active levels, Adv. Space Res., 34, 256-261, 2004. *
287 * Vane, G., R. O. Green, T. G. Chrien, H. T. Enmark, E. G. Hansen, and *
288 * W. M. Porter, The Airborne Visible/Infrared Imaging Spectrometer, *
289 * Remote Sens. Env., 44, 127-143, 1993. *
290 * Vane, G. (Ed), Airborne visible/infrared imaging spectrometer *
291 * (AVIRIS), JPL Publ. 87-38, Jet Propul. Lab, Pasadena, Calif., 1987.*
292 * Vermote, E., D. Tanre, J. L. Deuze, M. Herman, and J. J. Morcrette, *
293 * Second simulation of the satellite signal in the solar spectrum *
294 * (6S), 6S User's Guide Version 1, NASA-GSFC, Greenbelt, Maryland, *
295 * 134 pages, 1994. *
296 * *
297 ********************************************************************************
298 C--------1---------2---------3---------4---------5---------6---------7--
299  character (len=8) :: start_date, stop_date
300  character (len=10), dimension(12) :: timing
301 
302  call date_and_time(date=start_date, time=timing(1))
303 
304  WRITE(*,*)' atrem_app_refl_plus_gas_removal_2013'
305  WRITE(*,*)' (atrem for apparent reflectance + gas absorption '
306  WRITE(*,*)' effect removal from imaging spectrometer Data )'
307  WRITE(*,*)' April 2013'
308  WRITE(*,*)' '
309 
310 C get user input
311  CALL get_input
312  call date_and_time(time=timing(2))
313 
314 C adjust bottom layer boundary of the model atmosphere if the surface elevation
315 C of the imaging scene > 0
316  CALL model_adj
317 
318 C calculate solar and observational geometric factors
319  CALL geometry
320  call date_and_time(time=timing(3))
321 
322 C simulating scattering effects using the 6S code (Vermote & Tanre et al., 1996)
323 C--- CALL SSSSSS
324 C--- call date_and_time(time=timing(4))
325 
326 C initialize parameters for transmittance calculations
327  CALL init_speccal
328  call date_and_time(time=timing(5))
329 
330 C obtain solar irradiances above the atmosphere
331  CALL solar_irr_pc
332  call date_and_time(time=timing(6))
333 
334 C calculate transmittance spectra with different water vapor values and save
335 C them into a table
336  CALL tran_table
337  call date_and_time(time=timing(7))
338 
339 C process image data cube to derive surface reflectances and column water
340 C vapor values
341  CALL process_cube
342  call date_and_time(date=stop_date, time=timing(8))
343 
344  WRITE(*,*)'*** ATREM: Finished processing image. ***'
345 
346  WRITE(*,*) ' '
347  write(*,*)' ***TIMING SUMMARY**'
348  WRITE(*,*) ' '
349  write(*,*)'Started on : ', start_date
350  WRITE(*,*) ' '
351  write(*,*)'Started time : ', timing(1)
352  write(*,*)'Finished GET_INPUT : ', timing(2)
353  write(*,*)'Finished GEOMETRY : ', timing(3)
354 C--- write(*,*)'Finished 6S : ', timing(4)
355  write(*,*)'Finished INIT_SPECCAL : ', timing(5)
356  write(*,*)'Finished SOLAR_IRR_PC : ', timing(6)
357  write(*,*)'Finished TRAN_TABLE : ', timing(7)
358  write(*,*)'Finished PROCESS_CUBE at: ', timing(8)
359  WRITE(*,*) ' '
360  write(*,*)'Stopped on : ', stop_date
361 
362  stop
363  END
364 
365 ********************************************************************************
366 * *
367 * Name: GET_INPUT *
368 * Purpose: Reads and verifies all user provided input. *
369 * Parameters: none *
370 * Algorithm: Data is read in from standard input. The user is not prompted *
371 * interactively, so data should be redirected from an input file. *
372 * The data is validated after it is read in. If the data is *
373 * invalid, a message is printed and the program stops. *
374 * Globals used: TPVMR - two dimensional array containing 6 predefined *
375 * atmospheric models. The models contain values for *
376 * the number of atmospheric layer boundaries, altitude,*
377 * temperature, pressure, and water vapor volume mixing *
378 * ratio. *
379 * Global output: IH2OVP, ICO2, IO3, IN2O, ICO, ICH4, IO2 - set to 0 to *
380 * indicate that the gas should NOT be used in the *
381 * calculations, and set to 1 if the gas should be used *
382 * H(), T(), P(), VMR(), NB, NL, MODEL - altitude (km), *
383 * temperature (K), pressure (atm), water vapor volume *
384 * mixing ratio (ppm), number of atmospheric layer *
385 * boundaries, number of atmospheric layers (NB-1), *
386 * and model number for the selected atmospheric model. *
387 * VRTO3 - column amount of O3. *
388 * WAVOBS() - wavelengths for all channels. *
389 * FWHM() - resolutions for each channel. *
390 * NOBS - number of AVIRIS wavelengths. *
391 * HSURF - the mean surface elevation of the imaging scene *
392 * DLT, DLT2 - resolution, in units of nm, of input *
393 * spectra and resolution of output surface reflectance *
394 * spectra. If DLT2>DLT, output spectra are smoothed *
395 * using a gaussian function. *
396 * WNDOW1, WNDOW2, WP94C - center wavelength positions of two *
397 * broad window channels and one broad .94-um water *
398 * vapor channel. *
399 * WNDOW3, WNDOW4, W1P14C - center wavelength positions of two *
400 * broad window channels and one broad 1.14-um water *
401 * vapor channel. *
402 * NB1, NB2, NBP94, NB3, NB4, NB1P14 - number of individual *
403 * narrow channels that form the corresponding broad *
404 * window and water vapor absorption channels. *
405 * IMN, IDY, IYR, IH, IM, IS - month, day, year, hour, minute, *
406 * and second of measurement. *
407 * XLATD, XLATM, XLATS, LATHEM - degrees, minutes, seconds, and*
408 * hemisphere of latitude of measured area. *
409 * XLONGD, XLONGM, XLONGS, LNGHEM - degrees, minutes, seconds, *
410 * and hemisphere of latitude of measured area. *
411 * NAME_INSTRU: Imaging spectrometer name, e.g., AVIRIS, HYDICE*
412 * XPSS: = HSURF, an interface for using 6S. *
413 * XPPP: plane height (km, above the sea level). *
414 * *
415 * Return Codes: none *
416 * Special Considerations: None *
417 * *
418 ********************************************************************************
419 
420  SUBROUTINE get_input
421 
422  use cubeio
423 
424  include 'COMMONS_INC'
425 
426  INTEGER LUN_IN, LUN_OUT, LUN_VAP, I_RET
427  COMMON /inout_units/ lun_in, lun_out, lun_vap
428 
429 
430 C Common variables
431  dimension h(25), t(25), p(25), vmr(25)
432  dimension wavobs(1024),fwhm(1024)
433  dimension tpvmr(7,81)
434  CHARACTER*1 LATHEM, LNGHEM
435 
436  CHARACTER (LEN = 1000) :: FINAV,FOCUB,FOH2O
437  INTEGER SORDER,HDREC
438 
439 C Local variables
440  INTEGER ANS
441  CHARACTER (LEN = 1000) :: FINPWV,FTPVMR
442  LOGICAL GOOD_DATA
443  CHARACTER (LEN = 1000) :: FOUT1
444  INTEGER DIMS(4)
445  INTEGER ST_ORDER
446 
447  COMMON /getinput1/ ih2ovp,ico2,io3,in2o,ico,ich4,io2,ino2
448  COMMON /getinput3/ h,t,p,vmr,nb,nl,model,iaer,v,taer55,vrto3,sno2
449  COMMON /getinput4/ wavobs,fwhm
450  COMMON /getinput5/ nobs,hsurf,dlt,dlt2
451  COMMON /getinput6/ wndow1,wndow2,wp94c,wndow3,wndow4,w1p14c
452  COMMON /getinput7/ nb1,nb2,nbp94,nb3,nb4,nb1p14
453  COMMON /getinput8/ imn,idy,iyr,ih,im,is
454  COMMON /getinput9/ xlatd,xlatm,xlats,lathem
455  COMMON /getinput10/xlongd,xlongm,xlongs,lnghem
456  COMMON /getinput11/hdrec,nsamps,nlines,nbands,sorder
457  COMMON /getinput12/scalef
458  COMMON /tpvmr_init1/ tpvmr
459 
460 C Commons for use with the C programs for cube I/O
461  COMMON /outcube/ focub
462  COMMON /incube/ finav
463  COMMON /outh2ovap/ foh2o
464 
465 C Parameters for names of imaging spectrometers
466  CHARACTER (LEN = 80) :: NAME_INSTRU, NAMES(10)
467  COMMON /getinput13/ name_instru, names
468 
469  REAL XPSS, XPPP
470  COMMON /getinput14/ xpss, xppp
471 
472  REAL XVIEWD, XVIEWM, XVIEWS
473  REAL XAZMUD, XAZMUM, XAZMUS
474  COMMON /getinput15/ xviewd,xviewm,xviews, xazmud,xazmum,xazmus
475 
476  good_data = .true. ! initialize flag for good data
477 
478 C Get the name of imaging spectrometer, e.g., AVIRIS, HYDICE
479  READ(5,627) name_instru
480  627 FORMAT(a10)
481 C
482 C
483 C***Temp code for assigning names of imaging spectrometers. Based on these
484 C names, different scale factors should be used for different instrument
485 C when converting measured radiances to standard radiance units.
486 C*** The coding here may need to be moved to the file "COMMONS_INC"
487  names(1) = 'AVIRIS'
488  names(2) = 'HYDICE'
489  names(3) = 'HSI'
490  names(4) = 'TRWIS-III'
491  names(5) = 'PHYLLS'
492  names(6) = 'Hyperion'
493  names(7) = 'HICO'
494  names(8) = 'NIS'
495  names(9) = 'PRISM'
496  names(10)= 'OTHERS?'
497 C***End of temp coding ------------
498 
499  print*, 'Instrument Name: ', name_instru
500 
501 C Get the plane altitude (in km, above sea level). Must > HSURF.
502  READ(*,*) xppp
503 
504 C-- print*, 'XPPP = ', XPPP
505 
506 C Get date and time of measurements. Format: MM DD YYYY HH MM SS
507  READ (5,*) imn,idy,iyr,ih,im,is
508  IF((imn.LT.1).OR.(imn.GT.12).OR.(idy.LT.1).OR.(idy.GT.31).OR.
509  & (iyr.LT.1987).OR.(iyr.GT.2020)) THEN
510  WRITE(*,*)'Invalid date:',imn,idy,iyr
511  WRITE(*,*)'Format is MM DD YYYY where 0<MM<12, 0<DD<32,
512  &YYYY>1986.'
513  good_data = .false.
514  ENDIF
515  IF((ih.LT.0).OR.(ih.GT.24).OR.(im.LT.0).OR.(im.GT.60).OR.
516  & (is.LT.0).OR.(is.GT.60)) THEN
517  WRITE(*,*)'Invalid time:',ih,im,is
518  WRITE(*,*)'Format is HH MM SS where 0<HH<24, 0<MM<60, SS<60.'
519  good_data = .false.
520  ENDIF
521 
522 C Get latitude of measurements. Format: degrees minutes seconds
523  READ (5,*) xlatd,xlatm,xlats
524  IF((xlatd.GT.90).OR.(xlatm.GT.60).OR.(xlats.GT.60)) THEN
525  WRITE(*,*)'Invalid latitude:',xlatd,xlatm,xlats
526  WRITE(*,*)'Format: degrees minutes seconds'
527  WRITE(*,*)'Valid values are degrees < 90, minutes < 60,
528  &seconds<60.'
529  good_data = .false.
530  ENDIF
531 
532 C Get hemisphere of latitude. Format: "N" or "S"
533  READ (5,127) lathem
534  127 FORMAT(a1)
535  IF((lathem.NE.'N').AND.(lathem.NE.'S')) THEN
536  WRITE(*,*)'Invalid hemisphere value:',lathem
537  WRITE(*,*)'Valid values are "N" or "S".'
538  good_data = .false.
539  ENDIF
540 
541 C Get longitude of measurements. Format: degrees minutes seconds
542  READ (5,*) xlongd,xlongm,xlongs
543  IF((xlongd.GT.180).OR.(xlongm.GT.60).OR.(xlongs.GT.60)) THEN
544  WRITE(*,*)'Invalid longitude:',xlongd,xlongm,xlongs
545  WRITE(*,*)'Format: degrees minutes seconds'
546  WRITE(*,*)'Valid values are degrees < 180 minutes < 60,
547  &seconds<60.'
548  good_data = .false.
549  ENDIF
550 
551 C Get hemisphere of longitude. Format: "E" or "W"
552  READ (5,127) lnghem
553  IF((lnghem.NE.'E').AND.(lnghem.NE.'W')) THEN
554  WRITE(*,*)'Invalid hemisphere:',lnghem
555  WRITE(*,*)'Valid values are "E" or "W".'
556  good_data = .false.
557  ENDIF
558 
559 C Get view zenith angle. Format: degrees minutes seconds
560  READ (5,*) xviewd,xviewm,xviews
561  IF((xviewd.GT.90).OR.(xviewm.GT.60).OR.(xviews.GT.60)) THEN
562  WRITE(*,*)'Invalid latitude:',xviewd,xviewm,xviews
563  WRITE(*,*)'Format: degrees minutes seconds'
564  WRITE(*,*)'Valid values are degrees < 90, minutes < 60,
565  &seconds<60.'
566  good_data = .false.
567  ENDIF
568 
569 C Get relative azimuth angle between the solar ray and view direction with
570 C both projected on the ground. Format: degrees minutes seconds
571  READ (5,*) xazmud,xazmum,xazmus
572  IF((xazmud.GT.360.).OR.(xazmum.GT.60).OR.(xazmus.GT.60)) THEN
573  WRITE(*,*)'Invalid view azimuth ang.:',xazmud,xazmum,xazmus
574  WRITE(*,*)'Format: degrees minutes seconds'
575  WRITE(*,*)'Valid values are degrees < 360, minutes < 60,
576  &seconds<60.'
577  good_data = .false.
578  ENDIF
579 
580 C Get spectral resolutions (FWHM) in micron. If the value is 0, that means
581 C read the FWHM data from the wavelength file. Otherwise, a non-zero value
582 C will be used as a the constant FWHM value.
583  READ(*,*)dlt
584  IF (dlt.GT.0.0) THEN
585  IF((dlt.LT.0.0005).OR.(dlt.GT.0.1)) THEN
586  WRITE(*,*)'Invalid spectral resolution value:',dlt
587  WRITE(*,*)'Valid values are 0.0005-0.1 micron.'
588  good_data = .false.
589  ELSE
590 C Initialize the FWHM array
591 C*** DLTNEW=DLT/1000.
592  DO 521 i=1,1024
593  fwhm(i)=dlt
594  521 CONTINUE
595  ENDIF
596  ENDIF
597 
598 C Get imaging spectrometer wavelength file name
599  READ(*,85) finpwv
600  85 FORMAT(a1000)
601  OPEN(11,file=finpwv,status='OLD',iostat=istat)
602  IF(istat.NE.0) THEN
603  WRITE(*,*)'wavelength file did not open
604  &successfully:',finpwv
605  good_data = .false.
606  ELSE
607 
608 C Read wavelength data. For 1991 and before, the file contains the channel
609 C number and the wavelength for each channel. For 1992, the file also contains
610 C a third column with the resolution (FWHM) in micron for each channel.
611 
612  nobs = 1024 !initially set to max allowable # of observations
613  DO 510 i=1,nobs
614  IF (dlt.GT.0.) THEN
615  READ(11,*,END=520) X,WAVOBS(I)
616  ELSE
617  READ(11,*,END=520) X,WAVOBS(I), FWHM(I)
618  ENDIF
619  510 CONTINUE
620  520 nobs=i-1 !set to actual # of observations in input file
621  CLOSE(11)
622  ENDIF
623 
624 C Determine if channel ratio parameters should be read. Format: 0=no 1=yes
625  READ(*,*) ans
626  IF((ans.NE.0).AND.(ans.NE.1)) THEN
627  WRITE(*,*)'Invalid value to indicate whether the channel ratio
628  & parameters should be read:',ans
629  WRITE(*,*)'Valid values: 0=no 1=yes.'
630  good_data = .false.
631  ELSE
632  IF (ans .EQ. 1) THEN
633 C Get ranges on curve for the two atmospheric windows surrounding the .94-um
634 C water vapor absorption feature and the center point of the .94-um water
635 C vapor absorption feature. Enter:
636 C 1. the midpoint of first window (0.6-2.5)
637 C 2. number of points to average for first window (1-10)
638 C 3. the midpoint of second window (0.6-2.5)
639 C 4. number of points to average for second window (1-10)
640 C 5. the midpoint of water vapor absorption feature (0.6-2.5)
641 C 6. the number of points to average absorption feature (1-30)
642 
643  114 READ(*,*)wndow1, nb1, wndow2, nb2, wp94c, nbp94
644 
645  IF((wndow1.LT.0.6).OR.(wndow1.GT.2.5)) THEN
646  WRITE(*,*)'Invalid wavelength position for first atmospheric
647  & window in the .94-um region:',wndow1
648  WRITE(*,*)'Valid values are 0.6-2.5.'
649  good_data = .false.
650  ENDIF
651 
652  IF((nb1.LT.1).OR.(nb1.GT.50)) THEN
653  WRITE(*,*)'Invalid number of channels for first wavelength
654  &position in the .94-um region:',nb1
655  WRITE(*,*)'Valid values are 1-50.'
656  good_data = .false.
657  ENDIF
658 
659  IF((wndow2.LT.0.6).OR.(wndow2.GT.2.5)) THEN
660  WRITE(*,*)'Invalid wavelength position for second atmospheric
661  & window in the .94-um region:',wndow2
662  WRITE(*,*)'Valid values are 0.6-2.5.'
663  good_data = .false.
664  ENDIF
665 
666  IF((nb2.LT.1).OR.(nb2.GT.50)) THEN
667  WRITE(*,*)'Invalid number of channels for second wavelength
668  &position in the .94-um region:',nb2
669  WRITE(*,*)'Valid values are 1-50.'
670  good_data = .false.
671  ENDIF
672 
673  IF((wp94c.LE.wndow1).OR.(wp94c.GE.wndow2)) THEN
674  WRITE(*,*)'Invalid water vapor wavelength position for
675  & the .94-um region:',wp94c
676  WRITE(*,*)'Valid range is:',wndow1,' < value < ',
677  &wndow2,'.'
678  good_data = .false.
679  ENDIF
680 
681  IF((nbp94.LT.1).OR.(nbp94.GT.90)) THEN
682  WRITE(*,*)'Invalid number of channels for water vapor
683  &wavelength position in the .94-um region:',nbp94
684  WRITE(*,*)'Valid values are 1-90.'
685  good_data = .false.
686  ENDIF
687 
688 C Get ranges on curve for the two atmospheric windows surrounding the 1.14-um
689 C water vapor absorption feature and the center point of the 1.14-um water
690 C vapor absorption feature. Enter:
691 C 1. the midpoint of third window (0.6-2.5)
692 C 2. number of points to average for third window (1-10)
693 C 3. the midpoint of fourth window (0.6-2.5)
694 C 4. number of points to average for fourth window (1-10)
695 C 5. the midpoint of 1.14-um absorption feature (0.6-2.5)
696 C 6. the number of points to average for the absorption feature (1-30)
697 
698 
699  READ(*,*)wndow3, nb3, wndow4, nb4, w1p14c, nb1p14
700 
701  IF((wndow3.LT.0.6).OR.(wndow3.GT.2.5)) THEN
702  WRITE(*,*)'Invalid wavelength position for first atmospheric
703  & window in the 1.14-um region:',wndow3
704  WRITE(*,*)'Valid values are 0.6-2.5'
705  good_data = .false.
706  ENDIF
707 
708  IF((nb3.LT.1).OR.(nb3.GT.50)) THEN
709  WRITE(*,*)'Invalid number of channels for first wavelength
710  &position in the 1.14-um region:',nb3
711  WRITE(*,*)'Valid values are 1-50.'
712  good_data = .false.
713  ENDIF
714 
715  IF((wndow4.LT.0.6).OR.(wndow4.GT.2.5)) THEN
716  WRITE(*,*)'Invalid wavelength position for second atmospheric
717  & window in the 1.14-um region:',wndow4
718  WRITE(*,*)'Valid values are 0.6-2.5'
719  good_data = .false.
720  ENDIF
721 
722  IF((nb4.LT.1).OR.(nb4.GT.50)) THEN
723  WRITE(*,*)'Invalid number of channels for second wavelength
724  &position in the 1.14-um region:',nb4
725  WRITE(*,*)'Valid values are 1-50.'
726  good_data = .false.
727  ENDIF
728 
729  IF((w1p14c.LE.wndow3).OR.(w1p14c.GE.wndow4)) THEN
730  WRITE(*,*)'Invalid water vapor wavelength position for
731  &1.14-um:',w1p14c
732  WRITE(*,*)'Valid range is:',wndow3,' < value <',
733  &wndow4,'.'
734  good_data = .false.
735  ENDIF
736 
737  IF((nb1p14.LT.1).OR.(nb1p14.GT.110)) THEN
738  WRITE(*,*)'Invalid number of channels for water vapor
739  &wavelength position in the 1.14-um region:',nb1p14
740  WRITE(*,*)'Valid values are 1-110.'
741  good_data = .false.
742  ENDIF
743  ELSE
744 
745 C Initialize default atmospheric window and water vapor absorption regions
746 C for 3-channel ratioing calculations.
747  wndow1 = 0.865
748  nb1 = 3
749  wndow2 = 1.030
750  nb2 = 3
751  wp94c = 0.940
752  nbp94 = 5
753  wndow3 = 1.050
754  nb3 = 3
755  wndow4 = 1.235
756  nb4 = 3
757  w1p14c = 1.1375
758  nb1p14 = 7
759  ENDIF
760  ENDIF
761 
762 C Determine which atmospheric model to use
763 C 1 = tropical
764 C 2 = mid latitude summer
765 C 3 = mid latitude winter
766 C 4 = subarctic summer
767 C 5 = subarctic winter
768 C 6 = US standard 1962
769 C 7 = user defined model
770  READ(*,*)model
771  IF((model.LT.1).OR.(model.GT.7)) THEN
772  WRITE(*,*)'Invalid atmospheric model value:',model
773  WRITE(*,*)'Valid values are 1-7.'
774  good_data = .false.
775  model = 1 ! init model to a valid value for use below
776  ENDIF
777 
778  IF (model.EQ.7) THEN
779 C Get file name of atmospheric model information
780  READ(*,101) ftpvmr
781  101 FORMAT(a1000)
782  OPEN(20,file=ftpvmr,status='OLD',iostat=istat)
783  IF(istat.NE.0) THEN
784  WRITE(*,*)'Atmospheric model file did not open
785  & successfully:',ftpvmr
786  good_data = .false.
787  ELSE
788 C Read in number of boundaries, temperature (K), pressure (mb), and water
789 C vapor volume mixing ratio (VMR, in units of parts per million (ppm)) profiles
790  READ(20,*) nb
791  IF((nb.LE.0).OR.(nb.GT.25)) THEN
792  WRITE(*,*)'Invalid number of boundaries:',nb
793  WRITE(*,*)'Valid values are 1-25.'
794  good_data = .false.
795  ELSE
796  DO 300 i=1,nb
797  READ(20,*) h(i),p(i),t(i),vmr(i)
798  300 CONTINUE
799  CLOSE(20,iostat=istat)
800  ENDIF
801  ENDIF
802  ELSE
803 C Initialize NB, H, P, T, VMR from predefined atmospheric model array
804  nb = tpvmr(model,1)
805  DO 120 i=1,nb
806  h(i) = tpvmr(model,2+(4*(i-1)))
807  p(i) = tpvmr(model,3+(4*(i-1)))
808  t(i) = tpvmr(model,4+(4*(i-1)))
809  vmr(i) = tpvmr(model,5+(4*(i-1)))
810  120 CONTINUE
811  ENDIF
812  nl=nb-1
813 
814  DO i = nb+1, 25
815  h(i) = 1000.
816  p(i) = 0.0
817  t(i) = 300.
818  vmr(i) = 0.0
819  END DO
820 
821 C Determine if various gases should be included in the calculations. Format:
822 C 1=yes or 0=no. The order should be:
823 C 1. water vapor
824 C 2. carbon dioxide
825 C 3. ozone
826 C 4. nitrous oxide
827 C 5. carbon monoxide
828 C 6. methane
829 C 7. oxygen
830 C 8. nitrogen dioxide
831 
832  READ(*,*)ih2ovp, ico2, io3, in2o, ico, ich4, io2, ino2
833  IF((ih2ovp.NE.0).AND.(ih2ovp.NE.1)) THEN
834  WRITE(*,*)'Invalid selection for H2O vapor:',ih2ovp
835  WRITE(*,*)'Valid values: 0=no, 1=yes.'
836  good_data = .false.
837  ENDIF
838 
839  IF((ico2.NE.0).AND.(ico2.NE.1)) THEN
840  WRITE(*,*)'Invalid selection for CO2:',ico2
841  WRITE(*,*)'Valid values: 0=no, 1=yes.'
842  good_data = .false.
843  ENDIF
844 
845  IF((io3.NE.0).AND.(io3.NE.1)) THEN
846  WRITE(*,*)'Invalid selection for O3:',io3
847  WRITE(*,*)'Valid values: 0=no, 1=yes.'
848  good_data = .false.
849  ENDIF
850 
851  IF((in2o.NE.0).AND.(in2o.NE.1)) THEN
852  WRITE(*,*)'Invalid selection for N2O:',in2o
853  WRITE(*,*)'Valid values: 0=no, 1=yes.'
854  good_data = .false.
855  ENDIF
856 
857  IF((ico.NE.0).AND.(ico.NE.1)) THEN
858  WRITE(*,*)'Invalid selection for CO:',ico
859  WRITE(*,*)'Valid values: 0=no, 1=yes.'
860  good_data = .false.
861  ENDIF
862 
863  IF((ich4.NE.0).AND.(ich4.NE.1)) THEN
864  WRITE(*,*)'Invalid selection for CH4:',ich4
865  WRITE(*,*)'Valid values: 0=no, 1=yes.'
866  good_data = .false.
867  ENDIF
868 
869  IF((io2.NE.0).AND.(io2.NE.1)) THEN
870  WRITE(*,*)'Invalid selection for O2:',io2
871  WRITE(*,*)'Valid values: 0=no, 1=yes.'
872  good_data = .false.
873  ENDIF
874 
875  IF((ino2.NE.0).AND.(ino2.NE.1)) THEN
876  WRITE(*,*)'Invalid selection for NO2:',ino2
877  WRITE(*,*)'Valid values: 0=no, 1=yes.'
878  good_data = .false.
879  ENDIF
880 
881 C Read in the column ozone amount and scaling factor for NO2 column amount.
882 C The typical ozone amount is: 0.28-0.55 (atm-cm). The built-in NO2
883 C column amount is 5.0E+15 molecules/cm^2.
884  READ(*,*)vrto3, sno2
885  IF((vrto3.LT.0.1).OR.(vrto3.GT.0.6)) THEN
886  WRITE(*,*)'Invalid vertical column amount of ozone:',vrto3
887  WRITE(*,*)'Valid values are 0.1-0.6.'
888  WRITE(*,*)'Reasonable values are 0.28-0.55.'
889  good_data = .false.
890  ENDIF
891 
892 C Get aerosol model, visibility, or aerosol optical depth at 0.55 um (TAER55)
893 C for SSSSS(). If V > 0, no TAER55 is needed in the input file. If V = 0,
894 C TAER55 is needed.
895 C IAER: 0=no aerosol, set V = -1
896 C IAER: 1=continental aerosol, 2=maritime aerosol, 3=urban aerosol, set V>5 km.
897 
898  READ(*,*) iaer,v
899  IF((iaer.LT.0).OR.(iaer.GT.6)) THEN
900  WRITE(*,*)'Invalid aerosol model value:',iaer
901  WRITE(*,*)'Valid values are 0-3.'
902  good_data = .false.
903  ENDIF
904 
905  IF(iaer.EQ.0) THEN
906  v = -1
907  taer55=0.0
908  ELSE
909  IF((v.LT.0).OR.(v.GT.300)) THEN
910  WRITE(*,*)'Invalid visibility value:',v
911  WRITE(*,*)'Value must be greater than 0 and less than 300.'
912  good_data = .false.
913  ENDIF
914  ENDIF
915 
916  IF(v.EQ.0) THEN
917  READ(*,*) taer55
918  IF((taer55.GT.10.).OR.(taer55.LE.0)) THEN
919  WRITE(*,*)'Invalid aerosol optical depth at 550 nm: ',taer55
920  WRITE(*,*)'Valid values are between 0 and 10.'
921  good_data = .false.
922  ENDIF
923  ENDIF
924 
925 C Get the mean surface elevation. Range: 0.0 - H(NB-1) km
926  READ(*,*)hsurf
927 C--
928  xpss = hsurf
929 C--
930  IF((hsurf.LT.0).OR.(hsurf.GT.h(nb-1))) THEN
931  WRITE(*,*)'Invalid surface elevation value:',hsurf
932  WRITE(*,*)'Value must be less than the maximum elevation in
933  & the atmospheric model.'
934  good_data = .false.
935  ENDIF
936 
937 C Test to see if plane altitude is greater than bottom surface elevation
938  IF(xppp .LT. hsurf) THEN
939  WRITE(*,*)'Invalid plane altitude:',xppp
940  WRITE(*,*)'Plane altitude must be greater than the bottom
941  & surface elevation.'
942  good_data = .false.
943  ENDIF
944 
945 C Read in file name of data cube and get data dimensions in DIMS, ST_ORDER
946 C If no SIPS header on data, then default dimensions are header size = 0,
947 C number of samples = 614, number of lines = 512, number of bands = 224,
948 C storage order = BIL
949  READ(*,85) finav
950 C------ CALL OPENINFILE(IRET,DIMS,ST_ORDER)
951  lun_in = 47
952  CALL openinfile(lun_in, i_ret)
953 
954  IF(i_ret.NE.0) THEN
955  WRITE(*,*)'Input image file did not open
956  & successfully:',finav
957  good_data = .false.
958  ENDIF
959 
960 C Determine if cube dimensions should be read in. If not, the dimensions
961 C in the image file header will be used.
962  READ(*,*) ans
963  IF((ans.NE.0).AND.(ans.NE.1)) THEN
964  WRITE(*,*)'Invalid value to indicate whether the cube
965  & dimensions should be read:',ans
966  WRITE(*,*)'Valid values: 0=no 1=yes.'
967  good_data = .false.
968  ELSE
969  IF (ans .EQ. 1) THEN
970 C Read the header size in bytes, samples, lines, bands, and
971 C storage order (0=BSQ, 1=BIP, 2=BIL)'
972  READ(*,*) hdrec, nsamps, nlines, nbands, sorder
973  IF ((hdrec.LT.0).OR.(nsamps.LE.0).OR.
974  & (nlines.LE.0).OR.(nbands.LE.0).OR.
975  & (sorder.LE.0).OR.(sorder.GT.2)) THEN
976  WRITE(*,*)'Invalid cube parameters:',hdrec,nsamps,nlines,
977  &nbands,sorder
978  WRITE(*,*)'Values must be greater than or equal to 1 and
979  &the storage order must be less than 3.'
980  WRITE(*,*)'Storage Order (0=BSQ, 1=BIP, 2=BIL)'
981  WRITE(*,*)' 0 = BSQ is not supported in this version of ATREM'
982  good_data = .false.
983  ENDIF
984  ELSE
985  hdrec = dims(1)
986  nsamps = dims(2)
987  nlines = dims(3)
988  nbands = dims(4)
989  sorder = st_order
990  ENDIF
991  ENDIF
992 
993 C Get the output reflectance cube file name. The output cube always has a 512
994 C byte header record.
995  READ(*,85) focub
996 C--- CALL OPENOUTFILE(IRET,512,NSAMPS,NLINES,NBANDS,SORDER)
997  lun_out = 48
998  CALL openoutfile(lun_out, i_ret)
999 
1000  IF(i_ret.NE.0) THEN
1001  WRITE(*,*)'Output image file did not open
1002  & successfully:',focub
1003  good_data = .false.
1004  ENDIF
1005 
1006 C Get the resolution of output spectra. Format: 0-100 nm.
1007  READ(*,*)dlt2
1008  IF((dlt2.LT.0.).OR.(dlt2.GT.100.)) THEN
1009  WRITE(*,*)'Invalid output resolution value:',dlt2
1010  WRITE(*,*)'Valid values are 0-100 nm.'
1011  good_data = .false.
1012  ENDIF
1013 
1014 C Get the scaling factor for output reflectance values. Range: 1. to 32000.
1015  READ(*,*)scalef
1016  IF((scalef.LT.1.).OR.(scalef.GT.32000.)) THEN
1017  WRITE(*,*)'Invalid output resolution value:',scalef
1018  WRITE(*,*)'Valid values are 1-32000 nm.'
1019  good_data = .false.
1020  ENDIF
1021 
1022 C Get the output water vapor file name. The water vapor output has a 512-byte
1023 C header, is always one band, and is in BSQ format (BSQ=0).
1024  READ(*,85) foh2o
1025 C--- CALL OPENVAPFILE(IRET,512,NSAMPS,NLINES,1,0)
1026  lun_vap = 49
1027  CALL openvapfile(lun_vap, i_ret)
1028 
1029  IF(i_ret.NE.0) THEN
1030  WRITE(*,*)'Output water vapor file did not open
1031  & successfully:',foh2o
1032  good_data = .false.
1033  ENDIF
1034 
1035 C Get the file name for output table consisting of atmospheric
1036 C transmission spectra at the solar and observational geometry but
1037 C with varying amounts of water vapor.
1038 
1039  READ(*,85) fout1
1040  OPEN(35,file=fout1,status='UNKNOWN',iostat=istat)
1041  IF(istat.NE.0) THEN
1042  WRITE(*,*)'Output file did not open successfully:',fout1
1043  good_data = .false.
1044  ENDIF
1045 
1046  IF(good_data .NEQV. .true.) THEN
1047  WRITE(*,*)'**** ERROR: Invalid input detected. ****'
1048  stop
1049  ENDIF
1050 
1051 C--D WRITE(*,*)'IH2OVP=',IH2OVP,' ICO2=',ICO2,' IO3=',IO3,' IN2O=',IN2O
1052 C--D WRITE(*,*)'ICO=',ICO,' ICH4=',ICH4,' IO2=',IO2
1053 C--D WRITE(*,*)'H T P VMR'
1054 C--D DO 998 I=1,25
1055 C--D 998 WRITE(*,*)H(I),' ',T(I),' ',P(I),' ',VMR(I)
1056 C--D WRITE(*,*)'NB=',NB,' NL=',NL,' MODEL=',MODEL
1057 C--D WRITE(*,*)' IAER=',IAER,' V=',V,' VRTO3=',VRTO3
1058 C--D DO 997 I=1,NOBS
1059 C--D 997 WRITE(*,*)'I=',I,' WAVOBS(I)=',WAVOBS(I),' FWHM(I)=',FWHM(I)
1060 C--D WRITE(*,*)'NOBS=',NOBS,' HSURF=',HSURF,' DLT=',DLT,' DLT2=',DLT2
1061 C--D WRITE(*,*)'SCALEF=',SCALEF
1062 C--D WRITE(*,*)'WNDOW1=',WNDOW1,' WNDOW2=',WNDOW2,' WP94C=',WP94C
1063 C--D WRITE(*,*)'WNDOW3=',WNDOW3,' WNDOW4=',WNDOW4,' W1P14C=',W1P14C
1064 C--D WRITE(*,*)'NB1=',NB1,' NB2=',NB2,' NBP94=',NBP94
1065 C--D WRITE(*,*)'NB3=',NB3,' NB4=',NB4,' NB1P14=',NB1P14
1066 C--D WRITE(*,*)'IMN=',IMN,' IDY=',IDY,' IYR=',IYR
1067 C--D WRITE(*,*)'IH=',IH,' IM=',IM,' IS=',IS
1068 C--D WRITE(*,*)'XLATD=',XLATD,' XLATM=',XLATM,' XLATS=',XLATS
1069 C--D WRITE(*,147)LATHEM
1070 C--D 147 FORMAT('LATHEM=',A1)
1071 C--D WRITE(*,*)'XLONGD=',XLONGD,' XLONGM=',XLONGM,' XLONGS=',XLONGS
1072 C--D WRITE(*,148)LNGHEM
1073 C--D 148 FORMAT('LNGHEM=',A1)
1074 C--D WRITE(*,*)'HDREC=',HDREC,' NSAMPS=',NSAMPS,' NLINES=',NLINES
1075 C--D WRITE(*,*)'NBANDS=',NBANDS,' SORDER=',SORDER
1076 
1077  RETURN
1078  END
1079 
1080 ********************************************************************************
1081 * *
1082 * Name: MODEL_ADJ *
1083 * Purpose: resets the bottom boundary of the input model if the surface *
1084 * elevation is greater than 0, and calculate the column water vapor *
1085 * amount in the selected model. *
1086 * Parameters: none *
1087 * Algorithm: If the surface elevation > 0, the bottom layer temperature and *
1088 * water vapor volume mixing ratio are obtained through linear *
1089 * interpolation, while the bottom layer pressure is obtained *
1090 * through exponential interpolation. *
1091 * Globals used: H, T, P, VMR, NB, NL, HSURF - these values are adjusted if *
1092 * HSURF > 0 *
1093 * Global output: CLMVAP - Column water vapor amount in unit of cm. *
1094 * Q - Number of molecules above the surface at one *
1095 * atmosphere in units of molecules/cm**2 *
1096 * Return Codes: none *
1097 * Special Considerations: none *
1098 * *
1099 ********************************************************************************
1100  SUBROUTINE model_adj
1101 
1102  include 'COMMONS_INC'
1103 
1104 C Common variables
1105  dimension h(25), t(25), p(25), vmr(25)
1106 
1107  COMMON /getinput3/ h,t,p,vmr,nb,nl,model,iaer,v,taer55,vrto3,sno2
1108  COMMON /getinput5/ nobs,hsurf,dlt,dlt2
1109  COMMON /model_adj1/ clmvap,q
1110 
1111  dimension hp(25), tp(25), pp(25), vmrp(25)
1112  COMMON /model_adj2/ hp, tp, pp, vmrp
1113  COMMON /model_adj3/ k_plane, dvap_plane, dvap_layer,
1114  & dp_plane, dp_layer, clmvapp
1115  COMMON /model_adj4/ k_surf
1116 
1117  REAL XPSS, XPPP
1118  COMMON /getinput14/ xpss, xppp
1119 
1120 C H = Layer boundary height, SL=Slant path length, DSLODH=DSL/DH
1121 
1122  re=6380.0
1123 C Q=# of molecules above the surface at one atmosphere (molecules/cm**2)
1124  q=2.152e25
1125 
1126 C-- print*, 'original model atmosphere ...'
1127 C-- DO I = 1, 25
1128 C-- print*, 'H,T,P,VMR = ', H(I), T(I), P(I), VMR(I)
1129 C-- END DO
1130 
1131 C--- Convert the atmospheric pressure from "mb" to "atm.":
1132  DO i = 1, nb
1133  p(i) = p(i) / 1013.
1134  END DO
1135 C--- End of conversion
1136 
1137 C Convert the VMR from the ppm unit in the model to absolute unit.
1138  DO 310 i=1,nb
1139  310 vmr(i)=vmr(i)*1.0e-06
1140 
1141 C Do special processing if surface altitude is greater than 0
1142 C--- IF(HSURF.NE.0.0) THEN
1143 C Reset H(I) to a smaller value if HSURF.EQ.H(I) to avoid possible
1144 C problems in table searching using LOCATE
1145  DO 7455 i=1,nb
1146  7455 IF(hsurf.EQ.h(i)) hsurf=h(i)+0.0001
1147 
1148 C Determine index of H() such that H(index) < HSURF < H(index+1)
1149  CALL locate(h,nb,hsurf,k)
1150 C
1151 C K_SURF is an index relative to the original model atmosphere (0 - 100 km)
1152  k_surf = k
1153 C--- print*, 'K_SURF = ', K_SURF
1154 C
1155  IF(k.EQ.0) THEN
1156  WRITE(6,5237)
1157  5237 FORMAT(2x,'***WARNING: Surface elevation smaller then lowest'/
1158  & 2x,'boundary of the model atmosphere.')
1159 
1160  k_surf = 1
1161  GOTO 5255
1162  ENDIF
1163 C
1164  IF(k.GT.0) THEN
1165  dhk=h(k+1)-h(k)
1166  dhs=hsurf -h(k)
1167 
1168 C linear interpolation for surface temperature (TSURF) and VMR ( VMRS)
1169  tsurf=t(k)+(dhs/dhk)*(t(k+1)-t(k))
1170  vmrs =vmr(k)+(dhs/dhk)*(vmr(k+1)-vmr(k))
1171 
1172 C exponential interpolation for surface pressure (PSURF)
1173  psurf=p(k)*exp(-alog(p(k)/p(k+1))*dhs/dhk)
1174  h(1)=hsurf
1175  p(1)=psurf
1176  t(1)=tsurf
1177  vmr(1)=vmrs
1178 
1179  nb=nb-k+1
1180  nl=nb-1
1181 C--- print*,'NL = ', NL, ' NB = ', NB, 'in Model_Adj'
1182 C
1183  DO 5240 i=2,nb
1184  h(i)=h(i+k-1)
1185  p(i)=p(i+k-1)
1186  t(i)=t(i+k-1)
1187  vmr(i)=vmr(i+k-1)
1188  5240 CONTINUE
1189 C Zero out pressures and VMRS of top atmospheric layers.
1190  DO 5245 i=nb+1,25
1191  h(i)=1000.
1192  p(i)=0.0
1193  t(i)=300.
1194  vmr(i)=0.0
1195  5245 CONTINUE
1196  ENDIF
1197 C--- ENDIF
1198 
1199  5255 CONTINUE
1200 
1201  amtvrt=0.0
1202 
1203  DO 350 i=1,nl
1204  damtvt=q*(p(i)-p(i+1))*(vmr(i)+vmr(i+1))/2.0
1205  amtvrt=amtvrt+damtvt
1206  350 CONTINUE
1207 
1208  clmvap=amtvrt/3.34e+22
1209 
1210  WRITE(*,*)'Column vapor amount in model atmosphere from ground'
1211  WRITE(*,*)' to space = ', clmvap, ' cm'
1212 
1213 C--- WRITE(*,*)'In MODEL_ADJ, NB = ',NB,' NL = ',NL
1214 C--- print*, 'After adjusting for elevated surface ...'
1215 C--- DO I = 1, 25
1216 C--- print*, 'H,T,P,VMR = ', H(I), T(I), P(I), VMR(I)
1217 C--- END DO
1218 
1219 C
1220 C
1221 C Setting the upward atmospheric path's T, P, and VMR profiles:
1222 C
1223 C 1st duplicate the entire atmospheric profiles from the downward path
1224 C to the upward path
1225 C
1226  DO i = 1, 25
1227  hp(i) = h(i)
1228  tp(i) = t(i)
1229  pp(i) = p(i)
1230  vmrp(i) = vmr(i)
1231  END DO
1232 C
1233 C
1234  hplane = xppp
1235 C Set the highest plane altitude to the upper bound of model atmosphere
1236 C--- IF(HPLANE.GT.100.0) HPLANE = 100. - 0.0001
1237  IF(hplane.GE.100.0) hplane = 100. - 0.0001
1238 C
1239 C Do special processing if the plane height (HPLANE) is greater than HP(1)
1240  IF(hplane.GT.hp(1)) THEN
1241 C Reset Plane altitude HPLANE (= XPPP) to a larger value if
1242 C HPLANE.EQ.HP(I) to avoid possible problems in table
1243 C searching using LOCATE
1244  DO 7456 i=1,25
1245  7456 IF(hplane.EQ.hp(i)) hplane=hp(i)-0.0001
1246 
1247 C Determine index of HP() such that HP(index) < HPLANE < H(index+1)
1248  CALL locate(hp,nb,hplane,kk)
1249 
1250  IF(kk.EQ.0) THEN
1251  WRITE(6,5239)
1252  5239 FORMAT(2x,'***WARNING: Plane altitude less then lowest'/
1253  & 2x,'boundary of the model atmosphere.')
1254  GOTO 5256
1255  ENDIF
1256 C
1257  IF(kk.GT.0) THEN
1258  dhkk = hp(kk+1) - hp(kk)
1259  dhss = hplane - hp(kk)
1260 
1261 C linear interpolation for plane temperature (TPLANE) and VMR ( VMRSP)
1262  tplane = tp(kk) + (dhss/dhkk)*(tp(kk+1)-tp(kk))
1263  vmrsp = vmrp(kk) + (dhss/dhkk)*(vmrp(kk+1)-vmrp(kk))
1264 
1265 C exponential interpolation for plane pressure (PPLANE)
1266  pplane = pp(kk)*exp(-alog(pp(kk)/pp(kk+1))*dhss/dhkk)
1267  hp(kk+1) = hplane
1268  pp(kk+1) = pplane
1269  tp(kk+1) = tplane
1270  vmrp(kk+1) = vmrsp
1271 C
1272 C Zero out pressures and VMRP of top atmospheric layers.
1273  IF(kk.LT.24) THEN
1274  DO i=kk+2,25
1275  hp(i)=1000.
1276  pp(i)=0.0
1277  tp(i)=300.
1278  vmrp(i)=0.0
1279  END DO
1280  END IF
1281 
1282  ENDIF
1283  ENDIF
1284 
1285  5256 CONTINUE
1286 
1287  amtvrtp=0.0
1288 
1289  DO 357 i=1,kk
1290  damtvtp=q*(pp(i)-pp(i+1))*(vmrp(i)+vmrp(i+1))/2.0
1291  amtvrtp=amtvrtp+damtvtp
1292  357 CONTINUE
1293 
1294  clmvapp=amtvrtp/3.34e+22
1295 
1296  WRITE(*,*)'Column vapor below plane (CLMVAPP) = ',
1297  & clmvapp, ' cm'
1298 
1299 C--- WRITE(*,*)'In MODEL_ADJ, NB = ',NB,' KK = ', KK
1300 C--- print*, 'After further adjusting for plane height ...'
1301 C--- DO I = 1, 25
1302 C--- print*, 'HP,TP,PP,VMRP = ', HP(I), TP(I), PP(I), VMRP(I)
1303 C--- END DO
1304 
1305 C--- Indices and parameters for the plane layer
1306  k_plane = kk
1307 
1308  dvap_plane = q*(pp(k_plane) - pp(k_plane+1))*
1309  & (vmrp(k_plane) + vmrp(k_plane+1))/2.0 / 3.34e+22
1310 
1311  dvap_layer = q*(p(k_plane) - p(k_plane+1))*
1312  & (vmr(k_plane) + vmr(k_plane+1))/2.0 / 3.34e+22
1313 
1314  dp_plane = pp(k_plane) - pp(k_plane+1)
1315  dp_layer = p(k_plane) - p(k_plane+1)
1316 
1317 C--- print*, 'K_PLANE, DVAP_PLANE, DVAP_LAYER = ',
1318 C--- & K_PLANE, DVAP_PLANE, DVAP_LAYER
1319 C--- print*, 'DP_PLANE, DP_LAYER = ', DP_PLANE, DP_LAYER
1320 
1321 
1322  RETURN
1323  END
1324 ********************************************************************************
1325 * *
1326 * Name: GEOMETRY *
1327 * Purpose: Calculates the solar and the observational geometric factors. *
1328 * Parameters: none *
1329 * Algorithm: The solar geometry was obtained based on the latitude, longitude,*
1330 * GMT time using programs written by W. Mankin at National Center *
1331 * for Atmospheric Research in Boulder, Colorado. The *
1332 * geometric factors for CO2, O3, N2O, CO, CH4, and O2 are based *
1333 * only on the solar and observational angles. Sixty artificial *
1334 * geometric factors for H2O are set up to produce a transmittance *
1335 * table for different atmospheric water vapor amounts. *
1336 * Globals used: VRTO3 - Column O3 amount in units of atm-cm *
1337 * IMN,IDY,IYR,IH,IM,IS - time and date of data measurements *
1338 * XLATD,XLATM,XLATS,LATHEM - Latitude of measured area *
1339 * XLONGD,XLONGM,XLONGS,LNGHEM - Longitude of measured area *
1340 * CLMVAP - Column water vapor in unit of cm in the model atmosphere *
1341 * Global output: *
1342 * SOLZNI,SOLAZ,OBSZNI,OBSPHI,IDAY - Solar zenith angle, solar azimuth *
1343 * angle, observational zenith angle, observational azimuth angle, *
1344 * and the day number in the year *
1345 * GCO2,GO3,GN2O,GCO,GCH4,GO2,SSH2O,GGEOM,TOTLO3 - Geometric factors for *
1346 * different gases, and total O3 amount in the sun-surface ray path. *
1347 * The geometric factor is defined as: if the vertical column amount *
1348 * of the gas is equal 1, then GGAS is equal to the total amount of *
1349 * the gas in the combined Sun-surface-sensor ray path. *
1350 * Return Codes: none *
1351 * Special Considerations: none *
1352 * *
1353 ********************************************************************************
1354 
1355  SUBROUTINE geometry
1356 
1357  include 'COMMONS_INC'
1358 
1359  dimension vapvrt(60), vap_slant(60)
1360  dimension md(12)
1361  dimension ssh2o(60)
1362  dimension h(25), t(25), p(25), vmr(25)
1363  CHARACTER*1 LATHEM,LNGHEM
1364 
1365 C The VAPVRT array contains 60 column water vapor values used for generating
1366 C a table of transmittance spectra with different amount of water vapor.
1367 C The values in VAPVRT is designed so that there is approximately 2% change
1368 C in the .94-um H2O channel ratio for column water vapor in the range .4-6 cm.
1369 
1370  DATA vapvrt/.00, .02, .06, .11, .16, .21, .26, .31, .36, .40,
1371  & .43, .46, .50, .54, .58, .62, .66, .70, .75, .80,
1372  & .86, .92, .98, 1.06,1.14, 1.22, 1.3, 1.4, 1.5, 1.6,
1373  & 1.7, 1.8, 1.9, 2.05, 2.2, 2.35, 2.55, 2.75, 2.95, 3.2,
1374  & 3.5, 3.8, 4.1, 4.4, 4.7, 5.0, 5.3, 5.6, 6.0, 6.4,
1375  & 7.0, 7.7, 8.5, 9.4,10.4, 11.6, 13.0, 15.0, 25.0, 50./
1376 
1377  DATA md/0,31,59,90,120,151,181,212,243,273,304,334/
1378 
1379  COMMON /getinput3/ h,t,p,vmr,nb,nl,model,iaer,v,taer55,vrto3,sno2
1380 
1381  COMMON /getinput8/ imn,idy,iyr,ih,im,is
1382  COMMON /getinput9/ xlatd,xlatm,xlats,lathem
1383  COMMON /getinput10/xlongd,xlongm,xlongs,lnghem
1384  COMMON /model_adj1/ clmvap,q
1385  COMMON /geometry1/ solzni,solaz,obszni,obsphi,iday
1386  COMMON /geometry2/ gco2,go3,gn2o,gco,gch4,go2,ssh2o,totlo3,ggeom
1387 
1388  COMMON /model_adj3/ k_plane, dvap_plane, dvap_layer,
1389  & dp_plane, dp_layer, clmvapp
1390 
1391  dimension g_vap(25), g_other(25)
1392  COMMON /geometry3/ g_vap, g_other, g_vap_equiv, vap_slant_mdl
1393 
1394  REAL XPSS, XPPP
1395  COMMON /getinput14/ xpss, xppp
1396 
1397  REAL XVIEWD, XVIEWM, XVIEWS
1398  REAL XAZMUD, XAZMUM, XAZMUS
1399  COMMON /getinput15/ xviewd,xviewm,xviews, xazmud,xazmum,xazmus
1400 C
1401  hplane = xppp
1402 C
1403 C VAP_SLANT is a new array for containing two-way total vapor amounts
1404 C Here the number "2" can be changed to other numbers, e.g., 2.5,
1405 C without any major effects in derived water vapor values from
1406 C imaging spectrometer data.
1407 C
1408  DO i = 1, 60
1409  vap_slant(i) = vapvrt(i) * 2.0
1410  END DO
1411 C
1412  xh=ih
1413  xm=im
1414  xs=is
1415  tt=6.28318*((xh)/24+xm/1440+xs/86400)
1416 
1417  xlat = xlatd + xlatm/60.0 + xlats/3600.0
1418  xlong = xlongd + xlongm/60.0 + xlongs/3600.0
1419 
1420 C northern hemisphere is positive, and western hemisphere is positive
1421 C (NOTE: this is not standard, but has been internally coded in this program)
1422  IF(lathem.NE.'N') xlat = -xlat
1423  IF(lnghem.NE.'W') xlong = -xlong
1424 
1425 325 xlatr = xlat/57.2958
1426  xlongr = xlong/57.2958
1427 
1428  CALL suncor(idy,imn,iyr,tt,dec,haz)
1429  CALL hazel(haz+tt-xlongr,dec,solaz,el,xlatr)
1430 C---Note: DEC, SOLAZ,and EL DEC are in radiant units
1431 
1432  IF (el .LE. 0.) THEN
1433  WRITE(*,*)'ERROR: Sun is below the horizon!!!'
1434  WRITE(*,*)'Check input date, time, latitude and longitude.'
1435  stop
1436  ENDIF
1437 
1438 C converting observational zenith angle and relative azimuth angle into degrees
1439  obszni = xviewd + xviewm/60.0 + xviews/3600.0
1440  obsphi = xazmud + xazmum/60.0 + xazmus/3600.0
1441 
1442  print*,'OBSZNI = ',obszni, ' OBSPHI = ',obsphi, ' degrees'
1443 
1444  obszni = obszni / 57.2958
1445  obsphi = obsphi / 57.2958
1446 
1447  solzni = 90.0/57.2958 - el
1448 
1449  ggeom = 1./cos(solzni) + 1./cos(obszni)
1450 
1451  gco2 = ggeom
1452 
1453  go3 = ggeom
1454  IF(hplane.LT.27.) go3 = ggeom - 1./cos(obszni)
1455 
1456  gn2o = ggeom
1457  gco = ggeom
1458  gch4 = ggeom
1459  go2 = ggeom
1460 
1461 C--- print*, 'GGEOM, GCO2, GO3, GN2O, GCO, GCH4, GO2 = '
1462 C--- print*, GGEOM, GCO2, GO3, GN2O, GCO, GCH4, GO2
1463 
1464  totlo3 = go3 * vrto3
1465 
1466 C--- print*, 'TOTLO3 = ', TOTLO3
1467 
1468 C Initialize newly created geometrical factors for each atmospheric
1469 C layers (here G_VAP and G_OTHER are true geometrical factors that
1470 C account for the actual Sun-surface-plane path lengths in the
1471 C model atmosphere):
1472 C
1473 C---For layers below the plane layer---
1474  DO i = 1, k_plane - 1
1475  g_vap(i) = ggeom
1476  g_other(i) = ggeom
1477  END DO
1478 C
1479 C---For layers above the plane layer
1480  DO i = k_plane + 1, 25
1481  g_vap(i) = ggeom - 1./cos(obszni)
1482  g_other(i) = ggeom - 1./cos(obszni)
1483  END DO
1484 C
1485 C---Special treatment for the plane layer to take account the
1486 C "shorter" upward path length
1487  g_vap(k_plane) = ggeom - 1./cos(obszni)
1488  & + dvap_plane/dvap_layer/cos(obszni)
1489  g_other(k_plane) = ggeom - 1./cos(obszni)
1490  & + dp_plane/dp_layer/cos(obszni)
1491 
1492 C--- print*, ' G_VAP, G_OTHER, I ='
1493 C--- DO I = 1, 25
1494 C--- print*, G_VAP(I), G_OTHER(I), I
1495 C--- END DO
1496 
1497 C Calculate the water vapor SCALING factor relative to the total amount
1498 C of water vapor in the model atmosphere in the L-shaped
1499 C Sun-surface-plane ray path.
1500 C
1501  vap_slant_mdl = clmvap/cos(solzni) + clmvapp/cos(obszni)
1502  print*, 'VAP_SLANT_MDL =', vap_slant_mdl, ' cm'
1503 C
1504 C The "equivalent" geometrical factor corresponding to the total
1505 C slant vapor amount of VAP_SLANT_MDL':
1506 C
1507  g_vap_equiv = vap_slant_mdl / clmvap
1508 C--- print*, 'G_VAP_EQUIV = ', G_VAP_EQUIV
1509 
1510  DO 310 i=1,60
1511  ssh2o(i) = vap_slant(i) / vap_slant_mdl
1512 C--- print*, 'SSH2O(I), I = ', SSH2O(I), I
1513  310 CONTINUE
1514 
1515 C Calculate the number of days that have passed in this year. Take leap year
1516 C into account.
1517  iday = md(imn) + idy
1518  lpyr = iyr - (4 * (iyr/4))
1519  IF((lpyr.EQ.0).AND.(iday.GT.59).AND.(imn.NE.2)) iday = iday + 1
1520 
1521 C--D WRITE(*,*)'SOLZNI=',SOLZNI,' SOLAZ=',SOLAZ,' OBSZNI=',OBSZNI
1522 C--D WRITE(*,*)'OBSPHI=',OBSPHI,' IDAY = ',IDAY
1523 C--D WRITE(*,*)'GCO2=',GCO2,' GO3=',GO3,' GN2O=',GN2O,' GCO=',GCO
1524 C--D WRITE(*,*)'GCH4=',GCH4,' GO2=',GO2,' TOTLO3=',TOTLO3
1525 C--D WRITE(*,*)'GGEOM=',GGEOM
1526 C--D DO 311 I=1,60
1527 C--D 311 WRITE(*,*)'I=',I,' SSH2O(I)',SSH2O(I)
1528 
1529  RETURN
1530  END
1531 
1532 ********************************************************************************
1533 * *
1534 * Name: INIT_SPECCAL *
1535 * Purpose: initialize global data for spectrum calculations. *
1536 * Parameters: none *
1537 * Algorithm: initialize data. *
1538 * Globals used: AH2O, APH2O, BH2O, BPH2O, SODLT, SOGAM, O3CF - Band model *
1539 * parameters for spectral calculations. *
1540 * WNDOW1, WNDOW2, WP94C, WNDOW3, WNDOW4, W1P14C - Center *
1541 * positions of window and water vapor absorption *
1542 * channels used in 3-channel ratio calculations. *
1543 * NB1, NB2, NBP94, NB3, NB4, NB1P14 - Number of narrow channels *
1544 * to form broader window and absorption channels. *
1545 * Global output: *
1546 * IH2OLQ,RLQAMT,NGASTT,NH2O,VSTART,VEND - Flag for including liquid *
1547 * water, liquid water amount (cm), total number of gases *
1548 * (typically 8), number of water vapor values, starting *
1549 * and ending wavelengths in internal calculations. *
1550 * NO3PT,NCV,NCVHAF,NCVTOT,VMIN,ISTART,IEND - Number of O3 abs. coef. *
1551 * points, parameters for gaussian function and spectral *
1552 * calculations. *
1553 * ISTCAL,IEDCAL,DP,PM,TM,VMRM - Parameters for spectral calculations *
1554 * IST1,IED1,IST2,IED2,ISTP94,IEDP94 - 3-channel ratioing parameters for *
1555 * the 0.94-um water vapor band. *
1556 * IST3,IED3,IST4,IED4,IST1P14,IED1P14 - 3-channel ratioing parameters for *
1557 * the 1.14-um water vapor band. *
1558 * WT1,WT2,WT3,WT4,JA - Relative weights for the four window channels *
1559 * used in channel-ratioing calculations. JA is a *
1560 * output parameter from a table searching routine. *
1561 * NCV2,NCVHF2,NCVTT2,ISTRT2,IEND2,FINST2 - Parameters for smoothing *
1562 * output reflectance spectra. *
1563 * NATOT,NBTOT,NCTOT,NDTOT - Number of channels for the four AVIRIS' *
1564 * grating spectrometers (A, B, C, and D). *
1565 * Return Codes: None. *
1566 * Special Considerations: Some parameters may need to be fine-tuned. *
1567 * *
1568 ********************************************************************************
1569 * Notes about water vapor VMRS and related quantities: *
1570 * *
1571 * VAPVRT(60) - a table containing 60 column vapor values (in unit of cm) *
1572 * *
1573 * VAP_SLANT(I) = VAPVRT(I) * 2.0, VAP_SLANT is a new table for containing *
1574 * two-way total vapor amounts. Here the number "2" can be *
1575 * changed to other numbers, e.g., 2.5, without major *
1576 * effects on retrieved water vapor values. *
1577 * *
1578 * G_VAP(I = 1,..., NL) = true vapor geometric factor for each layer in *
1579 * the model atmosphere (after adjusting for the elevated *
1580 * surface. *
1581 * *
1582 * VMRM(I) = VMRM(I)*G_VAP(I). The VMRS are multiplied by the geometrical *
1583 * factor. We can calculate the vapor transmittance on the *
1584 * Sun-surface-sensor path by assuming a vertical path in *
1585 * the model atmosphere with geometric-factor-adjusted VMRS. *
1586 * *
1587 * CLMVAP = vertical column amount from ground to space in model atmosphere*
1588 * CLMVAPP = vertical column amount from ground to aircraft or satellite *
1589 * sensor in model atmosphere *
1590 * Q = 2.152E25 = # of molecules above the surface at one atmosphere *
1591 * (in unit of molecules/cm**2) *
1592 * *
1593 * VAP_SLANT_MDL= CLMVAP/COS(SOLZNI) + CLMVAPP/COS(OBSZNI) = total amount *
1594 * of water vapor in the model atmosphere in the L-shaped *
1595 * Sun-surface-plane ray path. *
1596 * *
1597 * G_VAP_EQUIV = VAP_SLANT_MDL / CLMVAP = the "equivalent" geometrical *
1598 * factor corresponding to the total slant vapor amount *
1599 * VAP_SLANT_MDL and the column vapor amount CLMVAP. *
1600 * *
1601 * SSH2O(I) (I = 1, ..., 60) - a pure scaling factor relative to the total *
1602 * slant vapor amount of VAP_SLANT_MDL, and *
1603 * SSH2O(I) = VAP_SLANT(I) / VAP_SLANT_MDL *
1604 * *
1605 * SH2O = one value of SSH2O(I). SH2O is used during generation of the *
1606 * look-up table. *
1607 * *
1608 * VAPTT = VAP_SLANT_MDL*SH2O, is the absolute total vapor amount on the *
1609 * L-shaped path corresponding to a spectrum stored in the *
1610 * look-up table. *
1611 * *
1612 * CLMWVP = 0.5*(VAPTTA+VAPTTB)/G_VAP_EQUIV, is the retrieved column water *
1613 * vapor amount from imaging spectrometer data. *
1614 ********************************************************************************
1615  SUBROUTINE init_speccal
1616 
1617  include 'COMMONS_INC'
1618 
1619 C Common variables
1620  dimension h(25), t(25), p(25), vmr(25)
1621  dimension ssh2o(60)
1622  dimension wavobs(1024),fwhm(1024)
1623  dimension dp(25), pm(25), tm(25), vmrm(25)
1624  dimension finst2(100)
1625 
1626 C Local variables
1627 
1628  COMMON /getinput1/ ih2ovp,ico2,io3,in2o,ico,ich4,io2,ino2
1629  COMMON /getinput3/ h,t,p,vmr,nb,nl,model,iaer,v,taer55,vrto3,sno2
1630  COMMON /getinput4/ wavobs,fwhm
1631  COMMON /getinput5/ nobs,hsurf,dlt,dlt2
1632  COMMON /getinput6/ wndow1,wndow2,wp94c,wndow3,wndow4,w1p14c
1633  COMMON /getinput7/ nb1,nb2,nbp94,nb3,nb4,nb1p14
1634  COMMON /geometry2/ gco2,go3,gn2o,gco,gch4,go2,ssh2o,totlo3,ggeom
1635  COMMON /model_adj1/ clmvap,q
1636 
1637  COMMON /init_speccal3/ nh2o
1638  COMMON /init_speccal5/ dp,pm,tm,vmrm
1639  COMMON /init_speccal6/ ist1,ied1,ist2,ied2,istp94,iedp94
1640  COMMON /init_speccal7/ ist3,ied3,ist4,ied4,ist1p14,ied1p14
1641  COMMON /init_speccal8/ wt1,wt2,wt3,wt4,ja
1642 
1643  COMMON /init_speccal10/ ncv2,ncvhf2,ncvtt2,istrt2,iend2,finst2
1644  COMMON /init_speccal11/ natot,nbtot,nctot,ndtot
1645 
1646  dimension g_vap(25), g_other(25)
1647  COMMON /geometry3/ g_vap, g_other, g_vap_equiv, vap_slant_mdl
1648 
1649  dimension o3cf(5001)
1650  COMMON /o3cf_init1/ o3cf
1651 
1652  dimension tran_o3_std(5001)
1653  COMMON /init_speccal16/ tran_o3_std
1654 
1655  dimension rno2cf(5001)
1656  COMMON /no2cf_init1/ rno2cf
1657 
1658  dimension tran_no2_std(5001)
1659  COMMON /init_speccal17/ tran_no2_std
1660 
1661  COMMON /model_adj4/ k_surf
1662 
1663 C
1664  nh2o = 60 !Number of column water vapor values
1665 C
1666 C Initialization for high resolution spectral calculations -
1667 C First initialize arrays to smooth high resolution spectra (0.05 cm-1) to
1668 C medium resolution spectra (0.2 nm):
1669 C
1670 C*** Note: The array WAVNO_HI is not used in actual computing, the array
1671 C* should be removed from COMMONS_INC and the following DO LOOP
1672 C* of this program, the purpose of keeping WAVNO_HI now is to
1673 C* check array indices and make sure the correctness of the
1674 C* indices ****************************************************
1675  DO i = 1, np_hi
1676  wavno_hi(i) = 3000. + float(i-1)*dwavno ! Wavenumber (cm-1)of high
1677  ! resolution spectrum,
1678  tran_hi(i) = 1.0 ! 3000 - 18000 cm-1
1679  END DO
1680 C--- print*,'WAVNO_HI ,1, 10000,NP_HI = ',WAVNO_HI(1),
1681 C--- & WAVNO_HI(10000), WAVNO_HI(NP_HI)
1682 C
1683 C
1684  DO i = 1, np_med
1685  wavln_med(i) = vstart + float(i-1)*dwavln ! Wavelength of medium
1686  ! resolution spectrum,
1687  tran_med(i) = 1.0 ! FWHM=.2 nm, .56-3.1 um.
1688  END DO
1689 C-----
1690 C--- print*,'WAVLN_MED ,1, 10000,NP_MED = ',WAVLN_MED(1),
1691 C--- & WAVLN_MED(10000), WAVLN_MED(NP_MED)
1692 C-----
1693 C
1694 C NOTE: WAVLN_STD starts from 0.3 um, instead of 0.56 um
1695  DO i = 1, np_std
1696  wavln_std(i) = 0.3 + float(i-1)*dwavln ! Wavelength of medium
1697  ! resolution spectrum,
1698  tran_std(i) = 1.0 ! FWHM=.2 nm, .3-3.1 um.
1699  END DO
1700 C-----
1701 C--- print*,'WAVLN_STD ,1, 10000,NP_STD = ',WAVLN_STD(1),
1702 C--- & WAVLN_STD(10000), WAVLN_STD(NP_STD)
1703 C-----
1704 C
1705 C Note: The grids of WAVNO_HI do not match the grids of 10000./WAVLN_MED.
1706 C INDEX_MED is a specially designed index for finding close matches
1707 C between the two kinds of grids.
1708 C
1709  DO i = 1, np_med
1710  index_med(i) = ( (10000./wavln_med(i) - 3000.)/dwavno + 1.)
1711  END DO
1712 C-----
1713 C--- print*,'INDEX_MED ,1, 10000,NP_MED = ',INDEX_MED(1),
1714 C--- & INDEX_MED(10000), INDEX_MED(NP_MED)
1715 C-----
1716 C
1717 C Note: WAVLN_MED_INDEX(I) is very close to WAVLN_MED(I),
1718 C and WAVLN_MED_INDEX(I) >= WAVLN_MED(I)
1719 C
1720  DO i = 1, np_med
1721  wavln_med_index(i) = 10000. /(float(index_med(i)-1)*dwavno
1722  & + 3000.)
1723  END DO
1724 C-----
1725 C--- print*,'WAVLN_MED_INDEX ,1, 10000,NP_MED = ',WAVLN_MED_INDEX(1),
1726 C--- & WAVLN_MED_INDEX(10000), WAVLN_MED_INDEX(NP_MED)
1727 C-----
1728 C
1729 
1730  DO i = 1, np_med
1731  fwhm_wavno(i) = 10000.*dlt_med
1732  & /(wavln_med_index(i)*wavln_med_index(i))
1733  END DO
1734 C-----
1735 C--- print*,'FWHM_WAVNO ,1, 10000,NP_MED = ',FWHM_WAVNO(1),
1736 C--- & FWHM_WAVNO(10000), FWHM_WAVNO(NP_MED)
1737 C-----
1738 C
1739 
1740  DO i = 1, np_med
1741  ncvhf_wavno(i) = ( facdlt * fwhm_wavno(i) / dwavno + 1.)
1742  END DO
1743 C-----
1744 C--- print*,'NCVHF_WAVNO ,1, 10000,NP_MED = ',NCVHF_WAVNO(1),
1745 C--- & NCVHF_WAVNO(10000), NCVHF_WAVNO(NP_MED)
1746 C-----
1747 
1748 C Initialize arrays for smoothing medium resolution spectrum (DLT_MED = 0.2 nm,
1749 C and point spacing DWAVLN = 0.0001 micron) to coarser spectral
1750 C resolution data from imaging spectrometers.
1751 
1752  DO i = 1, nobs
1753  ncvhf(i) = ( facdlt * fwhm(i) / dwavln + 1.)
1754  END DO
1755 C
1756 C parameters and arrays to smooth output surface reflectance spectrum
1757  wavcv2=facdlt*dlt2
1758 
1759 C Find the largest value in the FWHM array, and use this value in calculation
1760 C of indices for smoothing output reflectance spectra. This smoothing
1761 C algorithm should work well with grating spectrometers having nearly
1762 C constant spectral resolutions, but not so well for prism spectrometers
1763 C having variable spectral resolution.
1764  dwvavr = fwhm(1)
1765 
1766  DO i = 2, nobs
1767  IF(dwvavr.LT.fwhm(i)) dwvavr = fwhm(i)
1768  END DO
1769 
1770  rncv2=wavcv2/dwvavr
1771  ncv2=rncv2
1772  ncvhf2=ncv2+1
1773  ncvtt2=2*ncv2+1
1774 
1775  cons2=dlt2*sqrt(3.1415926/const1)
1776 
1777  IF (dlt2 .NE. 0.0) THEN
1778  sumins=0.0
1779  DO 585 i=ncvhf2,ncvtt2
1780  finst2(i)=exp(-const1*(float(i-ncvhf2)*dwvavr/dlt2)**2.0)
1781  sumins=sumins+finst2(i)
1782  585 CONTINUE
1783 
1784  DO 590 i=1,ncvhf2-1
1785  finst2(i)=finst2(ncvtt2-i+1)
1786  sumins=sumins+finst2(i)
1787  590 CONTINUE
1788 
1789  sumins=sumins*dwvavr
1790 
1791  DO 595 i=1,ncvtt2
1792  finst2(i)=finst2(i)*dwvavr/sumins
1793  595 CONTINUE
1794  ENDIF
1795 
1796  istrt2=ncvhf2
1797  iend2=nobs-ncvhf2
1798 
1799 C number of channels of the four AVIRIS spectrometers. These are used
1800 C in removing null AVIRIS radiance values in the overlap portions of two
1801 C adjacent spectrometers.
1802  nchnla=32
1803  nchnlb=64
1804  nchnlc=64
1805  nchnld=64
1806 
1807  natot=nchnla
1808  nbtot=nchnla+nchnlb
1809  nctot=nchnla+nchnlb+nchnlc
1810  ndtot=nchnla+nchnlb+nchnlc+nchnld
1811 
1812 C Resetting window wavelength positions and calculating weights for
1813 C window and absorption channels used in 3-channel ratioing.
1814  iwndw1=findmatch(wavobs,nobs,wndow1)
1815  iwndw2=findmatch(wavobs,nobs,wndow2)
1816 
1817  wndow1=wavobs(iwndw1)
1818  wndow2=wavobs(iwndw2)
1819 
1820  jj=mod(nb1,2)
1821  IF(jj.EQ.0) nb1=nb1+1
1822  kk=mod(nb2,2)
1823  IF(kk.EQ.0) nb2=nb2+1
1824  nb1haf=(nb1-1)/2
1825  nb2haf=(nb2-1)/2
1826 
1827  ist1=iwndw1-nb1haf
1828  ied1=iwndw1+nb1haf
1829  ist2=iwndw2-nb2haf
1830  ied2=iwndw2+nb2haf
1831 
1832  iwp94c=findmatch(wavobs,nobs,wp94c)
1833  wp94c=wavobs(iwp94c)
1834 
1835  ll=mod(nbp94,2)
1836  IF(ll.EQ.0) nbp94=nbp94+1
1837  nb3haf=(nbp94-1)/2
1838  istp94=iwp94c-nb3haf
1839  iedp94=iwp94c+nb3haf
1840 
1841  wt1=(wndow2-wp94c)/(wndow2-wndow1)
1842  wt2=(wp94c-wndow1)/(wndow2-wndow1)
1843 
1844  iwndw4=findmatch(wavobs,nobs,wndow3)
1845  iwndw5=findmatch(wavobs,nobs,wndow4)
1846 
1847  wndow3=wavobs(iwndw4)
1848  wndow4=wavobs(iwndw5)
1849 
1850  jj=mod(nb3,2)
1851  IF(jj.EQ.0) nb3=nb3+1
1852  kk=mod(nb4,2)
1853  IF(kk.EQ.0) nb4=nb4+1
1854 
1855  nb4haf=(nb3-1)/2
1856  nb5haf=(nb4-1)/2
1857 
1858  ist3=iwndw4-nb4haf
1859  ied3=iwndw4+nb4haf
1860  ist4=iwndw5-nb5haf
1861  ied4=iwndw5+nb5haf
1862  iw1p14c=findmatch(wavobs,nobs,w1p14c)
1863 
1864  w1p14c=wavobs(iw1p14c)
1865  ll=mod(nb1p14,2)
1866  IF(ll.EQ.0) nb1p14=nb1p14+1
1867  nb6haf=(nb1p14-1)/2
1868  ist1p14=iw1p14c-nb6haf
1869  ied1p14=iw1p14c+nb6haf
1870 
1871  wt3=(wndow4-w1p14c)/(wndow4-wndow3)
1872  wt4=(w1p14c-wndow3)/(wndow4-wndow3)
1873 
1874 C Initialization for searching water vapor table (TRNTBL)
1875  ja = 30
1876 C
1877 C
1878 C Calculate medium resolution O3 transmittances (0.2 nm, 2-way) with
1879 C a point spacing of 0.1 nm between 0.3 and 0.8 micron.
1880 C
1881  IF(io3.EQ.1) THEN
1882  DO i=1,no3pt
1883  tran_o3_std(i) = exp(-totlo3*o3cf(i))
1884  END DO
1885  END IF
1886 
1887 C
1888 C If O3 is not intended to be included in total atmospheric gaseous
1889 C transmittance calculations, assigning TRAN_O3_STD = 1.0:
1890 C
1891  IF(io3.NE.1) THEN
1892  DO i=1,no3pt
1893  tran_o3_std(i) = 1.0
1894  END DO
1895  END IF
1896 C
1897 C Calculate medium resolution NO2 transmittances (0.2 nm, 2-way) with
1898 C a point spacing of 0.1 nm between 0.3 and 0.8 micron.
1899 C
1900  no2pt = no3pt
1901  vrtno2 = 5.0e+15
1902  vrtno2 = sno2 * vrtno2
1903 
1904  gno2 = go3
1905  totno2 = gno2 * vrtno2
1906 
1907  IF(ino2.EQ.1) THEN
1908  DO i=1,no2pt
1909  tran_no2_std(i) = exp(-totno2*rno2cf(i))
1910  END DO
1911  END IF
1912 
1913 C---Temp Code:
1914 C--- OPEN(57,file='zzzzzz_tst_NO2_trn',status='unknown')
1915 C--- DO I = 1, NO2PT
1916 C--- write(57,*)WAVLN_STD(I), TRAN_NO2_STD(I)
1917 C--- END DO
1918 C--- CLOSE(57)
1919 C---End of Temp Code
1920 
1921 C--- print*,' TOTNO2 = ', TOTNO2
1922 C--- print*,'VRTNO2, SNO2, GNO2 = ', VRTNO2, SNO2, GNO2
1923 
1924 C
1925 C If NO2 is not intended to be included in total atmospheric gaseous
1926 C transmittance calculations, assigning TRAN_NO2_STD = 1.0:
1927 C
1928  IF(ino2.NE.1) THEN
1929  DO i=1,no2pt
1930  tran_no2_std(i) = 1.0
1931  END DO
1932  END IF
1933 
1934 
1935 C
1936 C Initial arrays for mean layer pressure and temperatures
1937 
1938  DO 320 i=1,nl
1939  dp(i)=p(i)-p(i+1)
1940  pm(i)=(p(i)+p(i+1))/2.0
1941  tm(i)=(t(i)+t(i+1))/2.0
1942  320 CONTINUE
1943 
1944 C
1945 C Calculate high resolution transmittances (0.05 cm-1) of CO2, N2O, CO,
1946 C CH4, and O2 in the 0.56 - 3.1 micron range, and save values for
1947 C calculating total atmospheric transmittances later.
1948 C Because water vapor amounts are allowed to vary,
1949 C the high resolution water vapor transmittances are calculated
1950 C in subroutines TRAN_TABLE and TRANCAL. TRAN_TABLE provides variable
1951 C water vapor amounts, and calls TRANCAL for the calculation of
1952 C corresponding vapor transmittance spectrum.
1953 C
1954 C*** Note: The array WAVNO_HI is not used in actual computing, the array
1955 C* should be removed from COMMONS_INC and the following DO LOOP
1956 C* of this program, the purpose of keeping WAVNO_HI now is to
1957 C* check array indices and make sure the correctness of the
1958 C* indices.
1959 C*
1960 C Initialize the TRAN_HI_OTHERS array for high resolution spectrum:
1961  DO i = 1, np_hi
1962  tran_hi_others(i) = 1.0
1963  END DO
1964 
1965 C---------------------------------------
1966 C For CO2 transmittance calculation -
1967  IF(ico2.EQ.1) THEN
1968 C---- SCLCO2=1.0
1969 C---- On 2/7/2013 B.-C. Gao made the modification - Increased SCLCO2 i
1970 C---- from 1.0 to 1.1 to reflect the fact that the CO2 VMR reached the
1971 C---- 2012 level of 390 ppmv.
1972 C----
1973  sclco2=1.1
1974 
1975  DO 322 i=1,nl
1976  vmrm(i)=sclco2*355.0*1.0e-06
1977 C Scale the VMRM by the two-way path geometrical factors. The geometric
1978 C factors, G_OTHER, varies with atmospheric layer number for
1979 C aircraft observational geometries.
1980  vmrm(i)= vmrm(i)*g_other(i)
1981  322 CONTINUE
1982 
1983  OPEN(31,file='abscf_co2_PC',status='old',
1984  & form='unformatted',access='direct',recl=4*300000)
1985 
1986  READ(31,rec=1) (wavno_hi(i), i = 1, 300000)
1987 
1988  DO j = k_surf, 19
1989  READ(31,rec=j+1) (a(i), i = 1, 300000)
1990 
1991  DO i = 1, np_hi
1992  tran_hi_others(i) = tran_hi_others(i) * exp( - a(i)*q*
1993  & dp(j-k_surf+1)*vmrm(j-k_surf+1) * 28.966 /
1994  & 6.0225e+23 / 1.0e-06)
1995  END DO
1996  END DO
1997 
1998  CLOSE(31)
1999  ENDIF
2000 
2001 C--------------------------------------------
2002 C For N2O transmittance calculation.
2003 
2004  IF(in2o.EQ.1) THEN
2005  DO 324 i=1,nl
2006  vmrm(i)=0.3*1.0e-06
2007  vmrm(i)= vmrm(i)*g_other(i)
2008  324 CONTINUE
2009 
2010  OPEN(31,file='abscf_n2o_PC',status='old',
2011  & form='unformatted',access='direct',recl=4*300000)
2012 
2013  READ(31,rec=1) (wavno_hi(i), i = 1, 300000)
2014 
2015  DO j = k_surf, 19
2016  READ(31,rec=j+1) (a(i), i = 1, 300000)
2017 
2018  DO i = 1, np_hi
2019  tran_hi_others(i) = tran_hi_others(i) * exp( - a(i)*q*
2020  & dp(j-k_surf+1)*vmrm(j-k_surf+1) * 28.966 /
2021  & 6.0225e+23 / 1.0e-06)
2022  END DO
2023  END DO
2024 
2025  CLOSE(31)
2026  ENDIF
2027 
2028 C--------------------------------------------
2029 C For CO transmittance calculation.
2030  IF(ico.EQ.1) THEN
2031  DO 325 i=1,nl
2032  vmrm(i)=0.1*1.0e-06
2033  vmrm(i)= vmrm(i)*g_other(i)
2034  325 CONTINUE
2035 
2036  OPEN(31,file='abscf_co_PC',status='old',
2037  & form='unformatted',access='direct',recl=4*300000)
2038 
2039  READ(31,rec=1) (wavno_hi(i), i = 1, 300000)
2040 
2041  DO j = k_surf, 19
2042  READ(31,rec=j+1) (a(i), i = 1, 300000)
2043 
2044  DO i = 1, np_hi
2045  tran_hi_others(i) = tran_hi_others(i) * exp( - a(i)*q*
2046  & dp(j-k_surf+1)*vmrm(j-k_surf+1) * 28.966 /
2047  & 6.0225e+23 / 1.0e-06)
2048  END DO
2049  END DO
2050 
2051  CLOSE(31)
2052  ENDIF
2053 
2054 C--------------------------------------------
2055 C For CH4 transmittance calculation.
2056 C For assigning CH4 VMRM
2057 C NOTE: The scaling factor of 0.8 for the CH4 VMRS was obtained by comparing
2058 C transmittance spectra calculated using our program, which is based on
2059 C the Malkmus narrow band spectral model, with a ratioed spectrum
2060 C provided by G. C. Toon at Jet Propulsion Laboratory (JPL). The JPL
2061 C ratio spectrum was obtained by ratioing a high resolution (0.005 cm-1)
2062 C solar spectrum measured at ground level against a solar spectrum
2063 C measured at an altitude of approximately 35 km with the same Fourier
2064 C Transform spectrometer. The high resolution ratio spectrum was
2065 C degraded to a resolution of 10 nm during our derivation of the
2066 C scaling factor for the CH4 VMRS.
2067 
2068  IF(ich4.EQ.1) THEN
2069 C*** SCLCH4=0.8
2070  sclch4=1.0
2071  DO 326 i=1,nl
2072  vmrm(i)=sclch4*1.6*1.0e-06
2073  vmrm(i)= vmrm(i)*g_other(i)
2074  326 CONTINUE
2075 
2076  OPEN(31,file='abscf_ch4_PC',status='old',
2077  & form='unformatted',access='direct',recl=4*300000)
2078 
2079  READ(31,rec=1) (wavno_hi(i), i = 1, 300000)
2080 
2081  DO j = k_surf, 19
2082  READ(31,rec=j+1) (a(i), i = 1, 300000)
2083 
2084  DO i = 1, np_hi
2085  tran_hi_others(i) = tran_hi_others(i) * exp( - a(i)*q*
2086  & dp(j-k_surf+1)*vmrm(j-k_surf+1) * 28.966 /
2087  & 6.0225e+23 / 1.0e-06)
2088  END DO
2089  END DO
2090 
2091  CLOSE(31)
2092  ENDIF
2093 
2094 C--------------------------------------------
2095 C For O2 transmittance calculation.
2096  IF(io2.EQ.1) THEN
2097 
2098 C***Modified by Bo-Cai Gao on 2/7/2013 to increase O2 absorption
2099 C--- coefficients by the factor SCL_O2 for wavelengths > 1.2 micron
2100 C--- in order to model properly the atmospheric O2 band centered
2101 C--- near 1.265 micron.
2102 
2103  scl_o2 = 2.60
2104 
2105  DO 327 i=1,nl
2106  vmrm(i)=0.21
2107  vmrm(i)= vmrm(i)*g_other(i)
2108  327 CONTINUE
2109 
2110  OPEN(31,file='abscf_o2_PC',status='old',
2111  & form='unformatted',access='direct',recl=4*300000)
2112 
2113  READ(31,rec=1) (wavno_hi(i), i = 1, 300000)
2114 
2115  DO j = k_surf, 19
2116  READ(31,rec=j+1) (a(i), i = 1, 300000)
2117 
2118 C***Modified by Bo-Cai Gao on 2/7/2013 to increase O2 absorption
2119 C--- coefficients by the factor SCL_O2 for wavelengths > 1.2 micron
2120 C--- i& < 1.3333 micron in order to model properly the atmospheric
2121 C--- O2 band centered near 1.265 micron.
2122 
2123  DO i = 1, 9000
2124  tran_hi_others(i) = tran_hi_others(i) * exp( - a(i)*q*
2125  & dp(j-k_surf+1)*vmrm(j-k_surf+1) * 28.966 /
2126  & 6.0225e+23 / 1.0e-06)
2127  END DO
2128 
2129  DO i = 9001, 106600
2130  tran_hi_others(i) = tran_hi_others(i) * exp( - a(i)*q*
2131  & dp(j-k_surf+1)*vmrm(j-k_surf+1)*scl_o2 *28.966/
2132  & 6.0225e+23 / 1.0e-06)
2133  END DO
2134 
2135  DO i = 106601, np_hi
2136  tran_hi_others(i) = tran_hi_others(i) * exp( - a(i)*q*
2137  & dp(j-k_surf+1)*vmrm(j-k_surf+1) * 28.966 /
2138  & 6.0225e+23 / 1.0e-06)
2139  END DO
2140  END DO
2141 
2142  CLOSE(31)
2143  ENDIF
2144 
2145 C--------------------------------------------
2146 C End of calculation of high resolution transmittances for CO2, N2O, CO,
2147 C CH4, and O2.
2148 C--------------------------------------------
2149 C Initial water vapor VMRs for repeated use in other subroutines and
2150 C adjust layered water vapor VMRM with geometrical factors.
2151  DO i = 1, nl
2152  vmrm(i) = (vmr(i)+vmr(i+1))/2.0
2153 C Scale the VMRM by the two-way path geometrical factors. The geometric
2154 C factors, G_VAP, varies with atmospheric layer number for
2155 C aircraft observational geometries.
2156  vmrm(i) = vmrm(i)*g_vap(i)
2157  END DO
2158 C--------------------------------------------
2159 C
2160 C--D WRITE(*,*)'NPSHIF=',NPSHIF,' DWAVLN=',DWAVLN
2161 C--D WRITE(*,*)'NO3PT=',NO3PT,' VMIN=',VMIN,' ISTART=',ISTART
2162 C--D WRITE(*,*)'IH2OLQ=',IH2OLQ,' RLQAMT=',RLQAMT,' NGASTT=',NGASTT
2163 C--D WRITE(*,*)'NH2O=',NH2O,' VSTART=',VSTART,' VEND=',VEND
2164 C--D WRITE(*,*)'IEND=',IEND
2165 C--D WRITE(*,*)'ISTCAL=',ISTCAL,' IEDCAL=',IEDCAL
2166 C--D DO 545 I=1,NL
2167 C--D 545 WRITE(*,*)I,DP(I),PM(I),TM(I),VMRM(I)
2168 C--D WRITE(*,*)'IST1=',IST1,' IED1=',IED1,' IST2=',IST2,' IED2=',IED2
2169 C--D WRITE(*,*)'ISTP94=',ISTP94,' IEDP94=',IEDP94
2170 C--D WRITE(*,*)'IST3=',IST3,' IED3=',IED3,' IST4=',IST4,' IED4=',IED4
2171 C--D WRITE(*,*)'IST1P14=',IST1P14,' IED1P14=',IED1P14
2172 C--D WRITE(*,*)'WT1=',WT1,' WT2=',WT2,' WT3=',WT3,' WT4=',WT4,' JA=',JA
2173 C--D WRITE(*,*)'NCV2=',NCV2,' NCVHF2=',NCVHF2,' NCVTT2=',NCVTT2,
2174 C--D &' ISTRT2=',ISTRT2,' IEND2=',IEND2
2175 C--D DO 544 I=1,30
2176 C--D 544 WRITE(*,*)'I=',I,' FINST2(I)=',FINST2(I)
2177 C--D WRITE(*,*)'NATOT=',NATOT,' NBTOT=',NBTOT
2178 C--D WRITE(*,*)'NCTOT=',NCTOT,' NDTOT=',NDTOT
2179 
2180  RETURN
2181  END
2182 
2183 ********************************************************************************
2184 * *
2185 * Name: TRAN_TABLE *
2186 * Purpose: This subroutine generates a table consisting of 60 atmospheric *
2187 * transmittance spectra at the solar and observational *
2188 * geometry and with 60 column water vapor values. The table also *
2189 * includes the total amounts of column water vapor used in the *
2190 * calculations, and the 3-channel ratios calculated from the window *
2191 * and absorption channels in and around the 0.94- and 1.14-um water *
2192 * vapor bands. *
2193 * Parameters: none *
2194 * Algorithm: For each of the 60 water vapor amounts, calculate the *
2195 * atmospheric transmittance, and save *
2196 * Globals Used: NH2O - number of column water vapor values *
2197 * VAPTT - geometrically adjusted water vapor total *
2198 * R094 - channel ratio for .94 um region *
2199 * R114 - channel ratio for 1.14 um region *
2200 * TRNCAL - atmospheric transmittance spectra *
2201 * Global Output: VAPTOT() - array containing geometrically adjusted water *
2202 * vapor values *
2203 * ROP94() - array containing channel ratios for .94 um region *
2204 * R1P14() - array containing channel ratios for 1.14 um region*
2205 * TRNTBL() - 2 dimensional array containing one transmittance *
2206 * spectrum for each column water vapor amount *
2207 * Return Codes: none *
2208 * Special Considerations: none *
2209 * *
2210 ********************************************************************************
2211  SUBROUTINE tran_table
2212 
2213  include 'COMMONS_INC'
2214 
2215 C Common variables
2216  dimension trncal(1024)
2217  dimension wavobs(1024),fwhm(1024)
2218  dimension dp(25), pm(25), tm(25), vmrm(25)
2219  dimension ssh2o(60)
2220  dimension vaptot(60), r0p94(60), r1p14(60), trntbl(1024,60)
2221 
2222  dimension o3cf(5001)
2223  COMMON /o3cf_init1/ o3cf
2224 
2225  COMMON /getinput1/ ih2ovp,ico2,io3,in2o,ico,ich4,io2,ino2
2226  COMMON /getinput4/ wavobs,fwhm
2227  COMMON /getinput5/ nobs,hsurf,dlt,dlt2
2228  COMMON /init_speccal3/ nh2o
2229  COMMON /init_speccal5/ dp,pm,tm,vmrm
2230 
2231  COMMON /tran_table1/ sh2o,vaptot,r0p94,r1p14,trntbl
2232  COMMON /trancal1/ trncal,vaptt
2233  COMMON /geometry2/ gco2,go3,gn2o,gco,gch4,go2,ssh2o,totlo3,ggeom
2234  COMMON /chnlratio1/ r094,r114
2235 
2236 C For each water vapor amount, calculate the geometrically adjusted water
2237 C vapor amount, the channel ratios, and the transmittance spectrum.
2238  vaptt = 0.
2239 
2240  DO 530 i=1,nh2o
2241  sh2o=ssh2o(i)
2242 
2243  CALL trancal
2244 
2245  vaptot(i)= vaptt
2246  r0p94(i) = r094
2247  r1p14(i) = r114
2248  DO 540 j=1,nobs
2249  trntbl(j,i)=trncal(j)
2250  540 CONTINUE
2251  530 CONTINUE
2252 
2253 C Write above calculated values to a file.
2254  nj=nh2o
2255  nk=float(nj)/10.
2256  njdiff=nj-nk*10
2257 
2258  DO 380 kk=1,nk
2259  WRITE(35,77) (vaptot(i+(kk-1)*10),i=1,10)
2260  77 FORMAT(7x,10(1x,e11.4))
2261  WRITE(35,77) ( r0p94(i+(kk-1)*10),i=1,10)
2262  WRITE(35,77) ( r1p14(i+(kk-1)*10),i=1,10)
2263  DO 385 ii=1,nobs
2264  385 WRITE(35,78) wavobs(ii),(trntbl(ii,i+(kk-1)*10),i=1,10)
2265  78 FORMAT(1x,f6.4,10(1x,e11.4))
2266  380 CONTINUE
2267 
2268  WRITE(35,77) (vaptot(i+nk*10),i=1,njdiff)
2269  WRITE(35,77) ( r0p94(i+nk*10),i=1,njdiff)
2270  WRITE(35,77) ( r1p14(i+nk*10),i=1,njdiff)
2271 
2272  IF(njdiff.GE.1) THEN
2273  DO 386 ii=1,nobs
2274  386 WRITE(35,78) wavobs(ii),(trntbl(ii,i+nk*10),i=1,njdiff)
2275  ENDIF
2276  CLOSE(35,iostat=istat)
2277 
2278  RETURN
2279  END
2280 
2281 ********************************************************************************
2282 * *
2283 * Name: TRANCAL *
2284 * Purpose: This program calculates combined transmittances of H2O, CO2, O3, *
2285 * N2O, CO, CH4, and O2. *
2286 * Parameters: none. *
2287 * Algorithm: The calculations were based on the line-by-line absorption *
2288 * parameters supplied by William R. Ridgway of NASA/GSFC. *
2289 * Global output:VAPTT - geometrically adjusted water vapor total. *
2290 * R094 - channel ratio for 0.94 um region. *
2291 * R114 - channel ratio for 1.14 um region. *
2292 * TRNCAL - total transmittances of all gases that match the *
2293 * resolutions of imaging spectrometers. *
2294 * Return Codes: none. *
2295 * Special Considerations: The high resolution (0.05 cm-1) line-by-line *
2296 * absorption parameters cover the 0.555 - 3.33 micron spectral range *
2297 * (3000 - 18000 cm-1). The medium resolution ozone absorption *
2298 * coefficients covers the 0.3-0.8 um spectral range. The line-by-line *
2299 * high resolution spectra were first smoothed to medium resolution *
2300 * spectra (resolution = 0.2 nm, wavelength spacing = 0.1 nm) covering *
2301 * the 0.56 - 3.1 micron spectral region. The medium resolution spectra *
2302 * of O3 and other gases are combined (in SUBROUTINE TRAN_SMOOTH) to form *
2303 * a single medium resolution spectrum from 0.3 to 3.1 micron. This *
2304 * combined spectrum (medium resolution) is then smoothed to lower *
2305 * resolutions to match the resolutions of imaging spectrometers. The *
2306 * smoothing is also done in SUBROUTINE TRAN_SMOOTH. *
2307 * *
2308 ********************************************************************************
2309 
2310  SUBROUTINE trancal
2311 
2312  include 'COMMONS_INC'
2313 
2314 C Common variables
2315  dimension trncal(1024)
2316  dimension wavobs(1024),fwhm(1024)
2317  dimension dp(25), pm(25), tm(25), vmrm(25)
2318  dimension vaptot(60), r0p94(60), r1p14(60), trntbl(1024,60)
2319 
2320 C Local variables
2321  dimension trans(1024)
2322  dimension trncv(1024)
2323  dimension trnstd_sm(1050)
2324  INTEGER INDEX
2325 
2326  COMMON /getinput4/ wavobs,fwhm
2327  COMMON /getinput5/ nobs,hsurf,dlt,dlt2
2328  COMMON /model_adj1/ clmvap,q
2329  COMMON /init_speccal5/ dp,pm,tm,vmrm
2330  COMMON /init_speccal6/ ist1,ied1,ist2,ied2,istp94,iedp94
2331  COMMON /init_speccal7/ ist3,ied3,ist4,ied4,ist1p14,ied1p14
2332  COMMON /tran_table1/ sh2o,vaptot,r0p94,r1p14,trntbl
2333  COMMON /trancal1/ trncal,vaptt
2334 
2335  dimension g_vap(25), g_other(25)
2336  COMMON /geometry3/ g_vap, g_other, g_vap_equiv, vap_slant_mdl
2337 
2338  COMMON /model_adj4/ k_surf
2339 
2340 C
2341 C Calculate water vapor transmittances with different scaling factors:
2342 C
2343  DO i = 1, np_hi
2344  tran_hi(i) = 1.0
2345  END DO
2346 
2347  OPEN(31,file='abscf_h2o_PC',status='old',
2348  & form='unformatted',access='direct',recl=4*300000)
2349 
2350  READ(31,rec=1) (wavno_hi(i), i = 1, 300000)
2351 
2352  DO j = k_surf, 19
2353  READ(31,rec=j+1) (a(i), i = 1, 300000)
2354 
2355  DO i = 1, np_hi
2356  tran_hi(i) = tran_hi(i) * exp( - a(i) * q *
2357  & dp(j-k_surf+1)*vmrm(j-k_surf+1) * sh2o * 28.966 /
2358  & 6.0225e+23 / 1.0e-06)
2359  END DO
2360  END DO
2361 
2362  CLOSE(31)
2363 
2364 
2365 C Multiplying the high resolution water vapor transmittance with combined
2366 C transmittances of CO2, N2O, CO, CH4, and O2:
2367 C
2368  DO i = 1, np_hi
2369  tran_hi(i) = tran_hi(i) * tran_hi_others(i)
2370  END DO
2371 
2372 C
2373 C Smooth the high resolution spectra to resolutions of measured spectrum
2374 C
2375  CALL tran_smooth
2376 C
2377 C--- DO 470 I=1,NOBS
2378 C--- WRITE(*,*)'I=',I,' TRNCAL(I)=',TRNCAL(I)
2379 C--- 470 CONTINUE
2380 
2381 C Calculate 3-channel ratio values, R094 and R114, from simulated spectra.
2382  CALL chnlratio
2383 
2384 C Total amount of water vapor (in unit of cm) corresponding to the spectrum.
2385  vaptt = vap_slant_mdl * sh2o
2386 C-- print*, 'VAPTT= ', VAPTT,'VAP_SLANT_MDL= ',VAP_SLANT_MDL
2387 
2388  9999 RETURN
2389  END
2390 
2391 ********************************************************************************
2392 * *
2393 * Name: TRAN_SMOOTH *
2394 * Purpose: This program is to smooth the line-by-line high resolution *
2395 * spectrum to lower resolution spectrum that matches the resolutions *
2396 * of imaging spectrometer data. *
2397 * Parameters: none. *
2398 * Algorithm: The smoothing is done in two stages. The 1st stage is to smooth *
2399 * the high resolution spectrum to medium resolution spectrum at a *
2400 * constant FWHM (0.2 nm) and a constant wavelength interval *
2401 * (0.1 nm). The 2nd stage smoothing is to smooth the medium *
2402 * resolution spectrum to resolutions of input imaging spectrometer *
2403 * data. *
2404 * Globals used: The global variables used are contained in the file *
2405 * "COMMONS_INC" *
2406 * Global output: TRNCAL - total transmittances of all gases that match the *
2407 * resolutions of imaging spectrometers. *
2408 * Return Codes: none. *
2409 * *
2410 ********************************************************************************
2411 
2412  SUBROUTINE tran_smooth
2413 C
2414  include 'COMMONS_INC'
2415 
2416  dimension tran_o3_std(5001)
2417  COMMON /init_speccal16/ tran_o3_std
2418 
2419  dimension tran_no2_std(5001)
2420  COMMON /init_speccal17/ tran_no2_std
2421 
2422  dimension wavobs(1024),fwhm(1024)
2423  COMMON /getinput4/ wavobs,fwhm
2424 
2425  dimension trncal(1024)
2426  COMMON /trancal1/ trncal,vaptt
2427 
2428  COMMON /getinput5/ nobs,hsurf,dlt,dlt2
2429 
2430 C First stage of smoothing - smooth line-by-line high resolution spectrum with
2431 C over 300,000 points (point spacing of 0.05 cm-1) to medium resolution
2432 C spectrum (resolution of 0.2 nm and point spacing of 0.1 nm) with about
2433 C 25,000 points.
2434 C
2435 C The smoothing of line-by-line spectrum is done in wavenumber domain. For
2436 C a spectrum with a constant 0.2 nm resolution in wavelength domain, it has
2437 C variable resolution in wavenumber domain. This effect is properly taken
2438 C care of in the design of smoothing functions.
2439 C
2440 C Because the high resolution spectrum is in wavenumber units (cm-1), while
2441 C the medium resolution spectrum is in wavelength units, the two kinds of
2442 C grids do not automatically match. In order to match the grids, arrays
2443 C of INDEX_MED and TRAN_MED_INDEX are specially designed. The desired
2444 C medium resolution spectrum, TRAN_MED, at constant 0.1 nm point spacing
2445 C and 0.2 nm resolution is obtained through linear interpolation of
2446 C TRAN_MED_INDEX array.
2447 C
2448  DO 466 j =1, np_med
2449 
2450  tran_med_index(j) = 0.0
2451  ncvtot_wavno = 2 * ncvhf_wavno(j) - 1
2452 
2453  sumins = 0.0
2454 
2455  DO 560 i = ncvhf_wavno(j), ncvtot_wavno
2456  finstr_wavno(i) =
2457  & exp( -const1*(float(i-ncvhf_wavno(j))*dwavno
2458  & /fwhm_wavno(j))**2.0)
2459  sumins = sumins + finstr_wavno(i)
2460  560 CONTINUE
2461 
2462  DO 565 i = 1, ncvhf_wavno(j)-1
2463  finstr_wavno(i) = finstr_wavno(ncvtot_wavno-i+1)
2464  sumins = sumins + finstr_wavno(i)
2465  565 CONTINUE
2466 
2467  sumins = sumins * dwavno
2468 
2469  DO 570 i = 1, ncvtot_wavno
2470  finstr_wavno(i) = finstr_wavno(i)*dwavno/sumins
2471  570 CONTINUE
2472 
2473 C*** !!!High resolution transmittances of CO2, N2O, CO, CH4, and O2 should
2474 C also be calculated somewhere else (wavelength start = 0.56 micron).
2475 C*** Here assuming TCO2*TN2O*TCO*TCH4*TO2 is already calculated previously,
2476 C i.e., TRANS(I) = TRAN_CO2(I)*TRAN_N2O(I)*TRAN_CO(I)*TRAN_CH4(I)*TRAN_O2(I)
2477 C and TRAN_HI(I) = TRAN_HI_H2O(I) * TRANS(I), and TRAN_HI_H2O for varying
2478 C water vapor amounts is calculated in this subroutine.
2479 
2480  DO 491 k = index_med(j)-(ncvhf_wavno(j)-1),
2481  & index_med(j)+ncvhf_wavno(j)-1
2482  tran_med_index(j) = tran_med_index(j) + tran_hi(k)*
2483  & finstr_wavno(k-index_med(j)+ncvhf_wavno(j))
2484  491 CONTINUE
2485 
2486  466 CONTINUE
2487 
2488 C
2489 C Linear interpolation to get TRAN_MED from TRAN_MED_INDEX:
2490 C (Note that WAVLN_MED_INDEX(J) >= WAVLN_MED(J) )
2491 C
2492  tran_med(1) = tran_med_index(1)
2493  tran_med(np_med) = tran_med_index(np_med)
2494 
2495  DO j = 2, np_med-1
2496  IF(wavln_med_index(j).LE.wavln_med(j)) THEN
2497  tran_med(j) = tran_med_index(j)
2498  ELSE
2499  dlt = wavln_med_index(j) - wavln_med_index(j-1)
2500  fjm1 = (wavln_med_index(j) - wavln_med(j)) /dlt
2501  fj = (wavln_med(j) - wavln_med_index(j-1))/dlt
2502  tran_med(j) = fjm1*tran_med_index(j-1) + fj*tran_med_index(j)
2503 C---
2504 C--- print*,j,fjm1,fj
2505 C---
2506  END IF
2507  END DO
2508 C
2509 C--- Here multiplying O3 and NO2 spectra and other spectrum at medium resolution:
2510 C
2511  DO i = 1, np_std
2512  tran_std(i) = 1.
2513  END DO
2514 
2515  DO i = 1, no3pt
2516  tran_std(i) = tran_o3_std(i) * tran_no2_std(i)
2517  END DO
2518 
2519  DO i = npshif+1, np_std
2520  tran_std(i) = tran_std(i)*tran_med(i-npshif)
2521  END DO
2522 
2523 
2524 C The 2nd stage of smoothing - smooth the medium resolution spectrum (resolution
2525 C of 0.2 nm and point spacing of 0.1 nm) with about 25,000 points to match
2526 C the coarser and variable resolution spectrum from imaging spectrometers.
2527 C
2528 C Initialize some index parameters:
2529  ia = 1000
2530 C
2531  DO 1466 j =1, nobs
2532 
2533  trncal(j) = 0.0
2534  tran_ia = 0.0
2535  tran_iap1 = 0.0
2536 
2537  ncvtot = 2 * ncvhf(j) - 1
2538 C---
2539 C--- print*,'J= ',j, 'NCVHF =', NCVHF(J), 'NCVTOT=',NCVTOT
2540 
2541 C Calculate instrumental response functions...
2542  sumins = 0.0
2543 
2544  DO 1560 i = ncvhf(j), ncvtot
2545  finstr(i) =
2546  & exp( -const1*(float(i-ncvhf(j))*dwavln
2547  & /fwhm(j))**2.0)
2548  sumins = sumins + finstr(i)
2549 1560 CONTINUE
2550 
2551  DO 1565 i = 1, ncvhf(j)-1
2552  finstr(i) = finstr(ncvtot-i+1)
2553  sumins = sumins + finstr(i)
2554 1565 CONTINUE
2555 
2556  sumins = sumins * dwavln
2557 
2558  DO 1570 i = 1, ncvtot
2559  finstr(i) = finstr(i)*dwavln/sumins
2560 1570 CONTINUE
2561 
2562 C--- IF(j.eq.87) then
2563 C--- print*, ' J = 87 =', J
2564 C--- do i = 1, NCVTOT
2565 C--- print*,i, FINSTR(i)
2566 C--- end do
2567 C--- end if
2568 
2569 C Index searching...
2570 C
2571  CALL hunt(wavln_std, np_std, wavobs(j), ia)
2572 C
2573 C--- print*,'J =', J, ' IA =', IA
2574 
2575 C Smoothing...
2576 C
2577  DO 1491 k = ia-(ncvhf(j)-1), ia+ncvhf(j)-1
2578  tran_ia = tran_ia + tran_std(k)*
2579  & finstr(k-ia+ncvhf(j))
2580 
2581 C--- IF(J.eq.1) then
2582 C--- print*,'j =', j, 'K = ',K, ' IA= ',IA,
2583 C--- & K-IA+NCVHF(J),
2584 C--- & FINSTR(K-IA+NCVHF(J))
2585 C--- End IF
2586 
2587 1491 CONTINUE
2588 
2589  ia_p1 = ia + 1
2590  DO 1492 k = ia_p1-(ncvhf(j)-1), ia_p1+ncvhf(j)-1
2591  tran_iap1 = tran_iap1 + tran_std(k)*
2592  & finstr(k-ia_p1+ncvhf(j))
2593 1492 CONTINUE
2594 C
2595 C Linear interpolation to get TRNCAL from TRAN_IA and TRAN_IAP1:
2596 C
2597  dlt_ia = wavln_std(ia_p1) - wavln_std(ia)
2598  fia = (wavln_std(ia_p1) - wavobs(j)) /dlt_ia
2599 C FIA_P1 = (WAVOBS(J) - WAVLN_STD(IA))/DLT_IA
2600  fia_p1 = 1. - fia
2601  trncal(j) = fia*tran_ia + fia_p1*tran_iap1
2602 C---
2603 C--- print*,'j=',j,'IA =',IA,'FIA =',FIA,'FIA_P1=',FIA_P1,DLT_IA
2604 C---
2605 C
2606 1466 CONTINUE
2607 
2608 C
2609 C
2610  RETURN
2611  END
2612 ********************************************************************************
2613 * *
2614 * Name: CHNLRATIO *
2615 * Purpose: Calculate 3-channel ratios. *
2616 * Parameters: none *
2617 * Algorithm: The 0.94-um water vapor absorption channel is ratioed against *
2618 * the linear combination of two window channels near 0.86 and *
2619 * 1.03 um to obtain one channel ratio for the 0.94-um band. *
2620 * Similar calculation is done for the 1.14-um water vapor band. *
2621 * Globals used: NB1,NB2,NBP94,NB3,NB4,NB1P14 - number of points used in *
2622 * channel ratios for both the .94- and 1.14-um regions*
2623 * IST1,IED1,IST2,IED2,ISTP94,IEDP94 - 3-channel ratioing *
2624 * parameters for the 0.94-um water vapor band *
2625 * IST3,IED3,IST4,IED4,IST1P14,IED1P14 - 3-channel ratioing *
2626 * parameters for the 1.14-um water vapor band. *
2627 * WT1,WT2,WT3,WT4,JA - Relative weights for the four window *
2628 * channels used in channel-ratioing calculations. JA *
2629 * is an output parameter from a table searching *
2630 * routine. *
2631 * TRNCAL - Atmospheric transmittance spectra. *
2632 * Global output:R094,R114 - 3-channel ratio values for the 0.94- and 1.14-um *
2633 * water vapor bands. *
2634 * Return Codes: none *
2635 * Special Considerations: none *
2636 * *
2637 ********************************************************************************
2638 
2639  SUBROUTINE chnlratio
2640 
2641 C Common variables
2642  dimension trncal(1024)
2643 
2644  COMMON /getinput7/ nb1,nb2,nbp94,nb3,nb4,nb1p14
2645  COMMON /init_speccal6/ ist1,ied1,ist2,ied2,istp94,iedp94
2646  COMMON /init_speccal7/ ist3,ied3,ist4,ied4,ist1p14,ied1p14
2647  COMMON /init_speccal8/ wt1,wt2,wt3,wt4,ja
2648  COMMON /trancal1/ trncal,vaptt
2649  COMMON /chnlratio1/ r094,r114
2650 
2651 C Calculate average of spectra over window and water vapor absorption regions.
2652  const1=0.0
2653  DO 560 i=ist1,ied1
2654  const1=const1+trncal(i)
2655  560 CONTINUE
2656  const1=const1/float(nb1)
2657 
2658  const2=0.0
2659  DO 570 i=ist2,ied2
2660  const2=const2+trncal(i)
2661  570 CONTINUE
2662  const2=const2/float(nb2)
2663 
2664  const3=0.0
2665  DO 575 i=istp94,iedp94
2666  const3=const3+trncal(i)
2667  575 CONTINUE
2668  const3=const3/float(nbp94)
2669 
2670  r094=const3/((wt1*const1) + (wt2*const2))
2671 
2672  const4=0.0
2673  DO 580 i=ist3,ied3
2674  const4=const4+trncal(i)
2675  580 CONTINUE
2676  const4=const4/float(nb3)
2677 
2678  const5=0.0
2679  DO 590 i=ist4,ied4
2680  const5=const5+trncal(i)
2681  590 CONTINUE
2682  const5=const5/float(nb4)
2683 
2684  const6=0.0
2685  DO 595 i=ist1p14,ied1p14
2686  const6=const6+trncal(i)
2687  595 CONTINUE
2688  const6=const6/float(nb1p14)
2689 
2690  r114=const6/((wt3*const4) + (wt4*const5))
2691 
2692  RETURN
2693  END
2694 
2695 ********************************************************************************
2696 * *
2697 * Name: PROCESS_CUBE *
2698 * Purpose: Processes an input cube one spectral slice at a time to derive *
2699 * surface reflectance for each spectrum and calculate the column *
2700 * water vapor amount for each pixel. The derived surface reflectance*
2701 * values are written to an output image file with the same dimensions*
2702 * as the input image, and the column water vapor amounts are written *
2703 * to a separate file as a single channel image. *
2704 * Parameters: none *
2705 * Algorithm: C programs are used to perform all of the file I/O. A spectral *
2706 * slice is read from the input file. Each spectrum in the slice *
2707 * is divided by the solar irradiance curve, then the water vapor *
2708 * ratio is calculated. The corresponding ratio and associated *
2709 * transmission spectrum is used to derive the surface reflectance *
2710 * values. *
2711 * *
2712 * The output is an image file in the same storage order format as *
2713 * the input image file. It has a 512-byte SIPS header. The values*
2714 * are surface reflectance values multiplied by an input scale *
2715 * factor to make them integer*2 data. *
2716 * *
2717 * The output water vapor file is in BSQ format. It has a 512-byte *
2718 * SIPS header. It contains 1 channel with each pixel value *
2719 * representing the column water vapor amount in units of cm *
2720 * multiplied by a scale factor of 1000 to make the values *
2721 * integer*2. *
2722 * *
2723 * Globals used: NSAMPS,NLINES,NBANDS - Input image dimensions used in I/O and *
2724 * cube processing. *
2725 * FPIN - File pointer for input cube passed to rdslice() *
2726 * FPOCUB - File pointer for output cube passed to wtslice() *
2727 * FPOH2O - File pointer for water vapor image passed to wtslice() *
2728 * SPATIAL,ST_SAMPLE,ST_ROW,SAMPLES_TODO,ROWS_TODO,ISLICE,DIMS - Cube *
2729 * dimensions used in rdslice() and wtslice(). *
2730 * BUFFER - Holds the slice of data read by rdslice() and written by *
2731 * wtslice(). *
2732 * H2OBUF - Holds the water vapor image. *
2733 * Global output: None. *
2734 * Return Codes: None. *
2735 * Special Considerations: None. *
2736 * *
2737 ********************************************************************************
2738  SUBROUTINE process_cube
2739 
2740  use cubeio
2741 
2742  INTEGER LUN_IN, LUN_OUT, LUN_VAP, I_RET
2743  COMMON /inout_units/ lun_in, lun_out, lun_vap
2744 
2745 C Common variables
2746  dimension yy(1024) ! Observed radiances
2747 
2748 C parameters for cube I/O
2749  CHARACTER (LEN = 1000) :: FINAV,FOCUB,FOH2O
2750  INTEGER*2 BUFFER(8388608)
2751  INTEGER*2 H2OBUF(8388608)
2752  INTEGER SORDER,HDREC
2753 
2754  INTEGER I_Sample, J_Line
2755 
2756 C Local variables
2757  REAL TBUF(1024)
2758  INTEGER OFFSET
2759 
2760  INTEGER J, ISAMP, IBAND
2761  real*4 scalef,clmwvp
2762 
2763 C Putting these large arrays in global memory forces memory allocation at
2764 C program initialization. Otherwise, the program chokes on the first
2765 C call to PROCESS_CUBE.
2766 C--- COMMON /DUMMYGLOB/ BUFFER,H2OBUF
2767 
2768  COMMON /proccube1/ yy
2769  COMMON /getinput11/hdrec,nsamps,nlines,nbands,sorder
2770  COMMON /getinput12/scalef
2771 
2772 C Commons for use with the C programs for cube I/O
2773  COMMON /outcube/ focub
2774  COMMON /incube/ finav
2775  COMMON /outh2ovap/ foh2o
2776 
2777 C Read input cube one slice at a time into the array BUFFER.
2778 C BUFFER contains a spectral slice of data. The first NBANDS elements are
2779 C the first spectrum, the second NBANDS elements are the second spectrum, etc
2780 C Process each spectrum, then write to output file one slice at a time.
2781  DO 11 j=1,nlines
2782 
2783 C Let us know the progress of program.
2784  WRITE(*,*)'LINE=',j
2785 
2786 C-Note: BUFFER is now passed in through COMMON /DUMMYGLOB/ BUFFER,H2OBUF
2787 C for reducing the required total computer memory
2788 C--- CALL RD_SPECTRAL_SLICE(J,HDREC,NSAMPS,NLINES,NBANDS,
2789 C--- & SORDER)
2790  CALL rd_slice(lun_in,nsamps,nbands,sorder,buffer)
2791 
2792 C For each spectrum in the slice, assign it to the YY array, and use it in
2793 C REFLDRV to create a surface reflectance spectrum.
2794  DO 46 isamp= 1,nsamps
2795  offset = (isamp-1) * nbands
2796  DO 45 iband=1,nbands ! fill YY array with this spectrum
2797  yy(iband) = buffer(offset + iband)
2798  45 CONTINUE
2799 
2800 C Write out sample input spectrum.
2801 C--D IF((ISAMP.EQ.2).AND.(J.EQ.3)) THEN
2802 C--D WRITE(*,*)'Islice=3, Isamp=2, Spectrum before REFLDRV'
2803 C--D DO 49 I=1,NBANDS
2804 C--D 49 WRITE(*,*)'I=',I,' YY(I)=',YY(I)
2805 C--D ENDIF
2806 
2807  CALL app_refl_with_gas_removal_drv(clmwvp)
2808 C--- CALL REFLDRV(CLMWVP)
2809 
2810 C The YY array now contains true surface reflectance values. Multiply
2811 C by input scale factor to convert values to integer*2 format.
2812  DO 55 iband=1,nbands
2813  tbuf(iband)=scalef*yy(iband)
2814  buffer(offset + iband) = nint(tbuf(iband))
2815  55 CONTINUE
2816 
2817 C Save the water vapor amount (also scaled by 1000 to convert to integer*2).
2818 C Later, H2OBUF will be written to the water vapor output file.
2819 
2820  h2obuf(isamp) = nint(1000.*clmwvp)
2821 
2822 C Write out sample processed spectrum.
2823 C--D IF((ISAMP.EQ.2).AND.(J.EQ.3)) THEN
2824 C--D WRITE(*,*)'Islice=3, Isamp=2, Spectrum after App_Refl_DRV'
2825 C--D DO 48 I=1,NBANDS
2826 C--D WRITE(*,*)'I=',I,' YY(I)=',YY(I)
2827 C--D 48 WRITE(*,*)'I=',I,' TBUF(I)=',TBUF(I)
2828 C--D ENDIF
2829 
2830  46 CONTINUE
2831 
2832 C Write slice BUFFER to output file. The output file is the same storage order
2833 C as the input file, and it always has a header that is one record (512 bytes).
2834 C-Note: BUFFER is now passed through COMMON /DUMMYGLOB/ BUFFER,H2OBUF
2835 C--- CALL WT_SPECTRAL_SLICE(FPOCUB,J,512,NSAMPS,NLINES,NBANDS,
2836 C--- & SORDER)
2837 
2838  CALL wt_slice(lun_out,nsamps,nbands,sorder,buffer)
2839 
2840  CALL wt_line(lun_vap,nsamps,h2obuf)
2841 
2842  11 CONTINUE
2843 
2844 C Write the water vapor image to a file. The file is one spatial image
2845 C with one header record (512 bytes).
2846 C-Note: H2OBUF is now passed through COMMON /DUMMYGLOB/ BUFFER,H2OBUF
2847 C--- CALL WT_SPATIAL_SLICE(FPOH2O,1,512,NSAMPS,NLINES,NBANDS,0)
2848 
2849  CALL closeinfile(lun_in)
2850  CALL closeoutfile(lun_out)
2851  CALL closevapfile(lun_vap)
2852 
2853  RETURN
2854  END
2855 
2856 ********************************************************************************
2857 * *
2858 * Name: App_Refl_With_Gas_Removal_DRV *
2859 * Purpose: To derive column water vapor amount and apparent reflectance curve *
2860 * with atmospheric gaseous absorption effect removed from an input *
2861 * spectrum using a 3-channel ratioing technique. *
2862 * Parameters: none *
2863 * Algorithm: the algorithm includes following procedures: *
2864 * 1. An input radiance spectrum is divided by a solar radiance curve above *
2865 * the atmosphere to derive 'APPARENT REFLECTANCE SPECTRUM T(Lambda)'*
2866 * 2. Two 3-channel ratios: T(0.94 um)/(WT1*T(0.86)+WT2*T(1.02)) *
2867 * and T(1.14 um)/(WT3*T(1.05)+WT4*T(1.23)), are calculated from *
2868 * the observed spectrum T(Lambda). *
2869 * 3. Based on the two ratios using look-up table procedures plus linear *
2870 * interpolation techniques, the column water vapor values and the *
2871 * atmospheric transmittance spectrum are derived. *
2872 * 4. The apparent reflectance spectrum is divided by the calculated *
2873 * transmittance spectrum to remove atmospheric gaseous absorption *
2874 * features. *
2875 * 5. The output apprent reflectance spectrum with atmospheric gaseous *
2876 * absorption effect is smoothed if desired. *
2877 * Globals used: *
2878 * Global output: *
2879 * Return Codes: none *
2880 * Special Considerations: None *
2881 * *
2882 ********************************************************************************
2883 
2884  SUBROUTINE app_refl_with_gas_removal_drv(CLMWVP)
2885 
2886 C Common variables
2887  dimension vaptot(60), r0p94(60), r1p14(60), trntbl(1024,60)
2888  dimension yirr(1024)
2889  dimension yy(1024)
2890  dimension rr(1024)
2891  dimension rotot(1050), ttot(1050), stot(1050)
2892  dimension finst2(100)
2893  dimension ssh2o(60)
2894  dimension trncal(1024)
2895 
2896 C Local variables
2897  dimension trncv(1024)
2898 
2899  COMMON /getinput5/ nobs,hsurf,dlt,dlt2
2900  COMMON /getinput7/ nb1,nb2,nbp94,nb3,nb4,nb1p14
2901  COMMON /getinput8/ imn,idy,iyr,ih,im,is
2902  COMMON /geometry2/ gco2,go3,gn2o,gco,gch4,go2,ssh2o,totlo3,ggeom
2903  COMMON /init_speccal3/ nh2o
2904  COMMON /init_speccal6/ ist1,ied1,ist2,ied2,istp94,iedp94
2905  COMMON /init_speccal7/ ist3,ied3,ist4,ied4,ist1p14,ied1p14
2906  COMMON /init_speccal8/ wt1,wt2,wt3,wt4,ja
2907  COMMON /init_speccal10/ ncv2,ncvhf2,ncvtt2,istrt2,iend2,finst2
2908  COMMON /init_speccal11/ natot,nbtot,nctot,ndtot
2909  COMMON /tran_table1/ sh2o,vaptot,r0p94,r1p14,trntbl
2910  COMMON /trancal1/ trncal,vaptt
2911  COMMON /solar_irr1/yirr
2912  COMMON /proccube1/ yy
2913  COMMON /sixs1/ rotot, ttot, stot
2914 
2915  dimension g_vap(25), g_other(25)
2916  COMMON /geometry3/ g_vap, g_other, g_vap_equiv, vap_slant_mdl
2917 
2918 C Parameters for plane observations:
2919  CHARACTER (LEN = 80) :: NAME_INSTRU, NAMES(10)
2920  COMMON /getinput13/ name_instru, names
2921 
2922 C Arrays related to look-up table:
2923 C VAPTOT: TOTAL SUN-SURFACE-SENSOR PATH WATER VAPOR IN UNITS OF CM
2924 C R0P94 : Tc(0.94 um)/(WT1*Tc(0.86)+WT2*Tc(1.02))
2925 C R1P14 : Tc(1.14 um)/(WT3*Tc(1.05)+WT4*Tc(1.23))
2926 
2927  dimension speca(1024),specb(1024),specav(1024)
2928 
2929  integer*4 itmp1
2930  real*4 xfa
2931  real*4 xfb
2932  real*4 xfc
2933  real*4 xfd
2934 
2935  DATA pi,dtorad /3.1415926,0.0174533/
2936 
2937 C
2938 C Loops for processing image data cube ...
2939  IF((name_instru.EQ.names(1)).AND.(iyr.LE.2009)) THEN
2940 
2941 C Ratioing observed spectrum against solar irradiance curve.
2942 C The radiance data from each year is scaled by a different constant.
2943 C XFA = factor for A spectrometer
2944 C XFB = factor for B spectrometer
2945 C XFC = factor for C spectrometer
2946 C XFD = factor for D spectrometer
2947 C
2948 C rewritten to handle each spectrometer differently
2949 C 9/14/1995 Roger N. Clark, U.S.G.S.
2950 
2951  IF(iyr.LE.1989) THEN
2952  xfa = 10.0
2953  xfb = 10.0
2954  xfc = 10.0
2955  xfd = 10.0
2956  ENDIF
2957 
2958  IF((iyr.GE.1990).AND.(iyr.LE.1991)) THEN
2959  xfa = 20.0
2960  xfb = 20.0
2961  xfc = 20.0
2962  xfd = 20.0
2963  ENDIF
2964 
2965  IF((iyr.GE.1992).AND.(iyr.LE.1994)) THEN
2966  xfa = 50.0
2967  xfb = 50.0
2968  xfc = 50.0
2969  xfd = 50.0
2970  ENDIF
2971 
2972  IF((iyr.GE.1995).AND.(iyr.LE.2009)) THEN
2973  xfa = 50.0
2974  xfb = 50.0
2975  xfc = 50.0
2976  xfd = 100.0
2977  ENDIF
2978 
2979 C AVIRIS A Spectrometer (ch = 1-32)
2980 
2981  IF(nobs.GE.32) THEN
2982  itmp1 = 32
2983  ELSE
2984  itmp1 = nobs
2985  ENDIF
2986 
2987  DO 551 i=1,itmp1
2988  yy(i)=pi*yy(i)/(xfa*yirr(i))
2989  551 CONTINUE
2990 
2991 C AVIRIS B Spectrometer (ch = 33-96)
2992 
2993  IF(nobs.GE.96) THEN
2994  itmp1 = 96
2995  ELSE if ((nobs.GE.33).and.(nobs.LT.96)) THEN
2996  itmp1 = nobs
2997  ELSE
2998  go to 559
2999  ENDIF
3000 
3001  DO 552 i=33,itmp1
3002  yy(i)=pi*yy(i)/(xfb*yirr(i))
3003  552 CONTINUE
3004 
3005 C AVIRIS C Spectrometer (ch = 97-160)
3006 
3007  IF(nobs.GE.160) THEN
3008  itmp1 = 160
3009  ELSE if ((nobs.GE.97).and.(nobs.LT.160)) THEN
3010  itmp1 = nobs
3011  ELSE
3012  go to 559
3013  ENDIF
3014 
3015  DO 553 i=97,itmp1
3016  yy(i)=pi*yy(i)/(xfc*yirr(i))
3017  553 CONTINUE
3018 
3019 C AVIRIS D Spectrometer (ch = 161-224)
3020 
3021  IF(nobs.GE.161) THEN
3022  itmp1 = nobs
3023  ELSE
3024  go to 559
3025  ENDIF
3026 
3027  DO 554 i=161,itmp1
3028  yy(i)=pi*yy(i)/(xfd*yirr(i))
3029  554 CONTINUE
3030 
3031 C branch point in case not all spectrometers are being done.
3032  559 CONTINUE
3033 
3034 C End of R. Clark contribution
3035  END IF
3036 C
3037 C Loops for processing AVIRIS data starting from year 2010---
3038  IF((name_instru.EQ.names(1)).AND.(iyr.GE.2010)) THEN
3039 
3040  f_av_1 = 30.
3041  f_av_2 = 60.
3042  f_av_3 = 120.
3043 
3044  DO i = 1, 110
3045  yy(i) = pi*yy(i) / (f_av_1 * yirr(i))
3046  END DO
3047 
3048  DO i = 111, 160
3049  yy(i) = pi*yy(i) / (f_av_2 * yirr(i))
3050  END DO
3051 
3052  DO i = 161, 224
3053  yy(i) = pi*yy(i) / (f_av_3 * yirr(i))
3054  END DO
3055 
3056  END IF
3057 
3058 C Loops for processing HYDICE data ...
3059 C
3060  IF(name_instru.EQ.names(2)) THEN
3061 
3062  f_hydice = 75.
3063 
3064  DO i = 1, nobs
3065  yy(i) = pi*yy(i) / (f_hydice * yirr(i))
3066  END DO
3067 
3068  END IF
3069 
3070 
3071 C Loops for processing HSI data ...
3072 C
3073  IF(name_instru.EQ.names(3)) THEN
3074 
3075  f_hsi = 100.
3076 
3077  DO i = 1, nobs
3078  yy(i) = pi*yy(i) / (f_hsi * yirr(i))
3079  END DO
3080 
3081  END IF
3082 
3083 
3084 C Loops for processing TRWIS-III data ...
3085 C
3086  IF(name_instru.EQ.names(4)) THEN
3087 
3088  f_trwis = 100.
3089 
3090  DO i = 1, nobs
3091  yy(i) = pi*yy(i) / (f_trwis * yirr(i))
3092  END DO
3093 
3094  END IF
3095 
3096 C Loops for processing EO-1 Hyperion data ...
3097 C
3098  IF(name_instru.EQ.names(6)) THEN
3099 
3100  f_hyperion_a = 40.
3101  f_hyperion_b = 80.
3102 
3103  DO i = 1, 70
3104  yy(i) = pi*yy(i) / (f_hyperion_a * yirr(i))
3105  END DO
3106 
3107  DO i = 71, nobs
3108  yy(i) = pi*yy(i) / (f_hyperion_b * yirr(i))
3109  END DO
3110 
3111  END IF
3112 
3113 C Loops for processing HICO data ...
3114 C
3115  IF(name_instru.EQ.names(7)) THEN
3116 
3117 C--- F_HICO = 50./1.32
3118  f_hico = 50.
3119 
3120  DO i = 1, nobs
3121  yy(i) = pi*yy(i) / (f_hico * yirr(i))
3122  END DO
3123 
3124  END IF
3125 
3126 C Loops for processing NIS (NEON Imaging Spectrometer) data ...
3127 C
3128  IF(name_instru.EQ.names(8)) THEN
3129 
3130  f_nis = 100.
3131 
3132  DO i = 1, nobs
3133  yy(i) = pi*yy(i) / (f_nis * yirr(i))
3134  END DO
3135 
3136  END IF
3137 
3138 C Loops for processing PRISM Imaging Spectrometer data ...
3139 C
3140  IF(name_instru.EQ.names(9)) THEN
3141 
3142  f_prism = 100.
3143 
3144  DO i = 1, nobs
3145  yy(i) = pi*yy(i) / (f_prism * yirr(i))
3146  END DO
3147 
3148  END IF
3149 
3150 C
3151 C Calculating 3-channel ratios from an observed spectrum, using a
3152 C look-up table procedure to derive the column amount of water vapor
3153 C and to find the associated simulated transmittance spectrum.
3154  const1=0.0
3155  DO 560 i=ist1,ied1
3156  const1=const1+yy(i)
3157  560 CONTINUE
3158  const1=const1/float(nb1)
3159 
3160  const2=0.0
3161  DO 570 i=ist2,ied2
3162  const2=const2+yy(i)
3163  570 CONTINUE
3164  const2=const2/float(nb2)
3165 
3166  const3=0.0
3167  DO 575 i=istp94,iedp94
3168  const3=const3+yy(i)
3169  575 CONTINUE
3170  const3=const3/float(nbp94)
3171 
3172  r094co=const3/((wt1*const1) + (wt2*const2))
3173  r094c =r094co
3174 
3175  IF(r094co.GT.1.0) THEN
3176  const1=0.0
3177  DO 561 i=ist1,ied1
3178  const1=const1+1.0/yy(i)
3179  561 CONTINUE
3180  const1=const1/float(nb1)
3181 
3182  const2=0.0
3183  DO 571 i=ist2,ied2
3184  const2=const2+1.0/yy(i)
3185  571 CONTINUE
3186  const2=const2/float(nb2)
3187 
3188  const3=0.0
3189  DO 576 i=istp94,iedp94
3190  const3=const3+1.0/yy(i)
3191  576 CONTINUE
3192  const3=const3/float(nbp94)
3193 
3194  r094c=const3/((wt1*const1) + (wt2*const2))
3195  ENDIF
3196 
3197  const4=0.0
3198  DO 580 i=ist3,ied3
3199  const4=const4+yy(i)
3200  580 CONTINUE
3201  const4=const4/float(nb3)
3202 
3203  const5=0.0
3204  DO 590 i=ist4,ied4
3205  const5=const5+yy(i)
3206  590 CONTINUE
3207  const5=const5/float(nb4)
3208 
3209  const6=0.0
3210  DO 595 i=ist1p14,ied1p14
3211  const6=const6+yy(i)
3212  595 CONTINUE
3213  const6=const6/float(nb1p14)
3214 
3215  r114co=const6/((wt3*const4) + (wt4*const5))
3216  r114c =r114co
3217 
3218  IF(r114co.GT.1.0) THEN
3219  const4=0.0
3220  DO 581 i=ist3,ied3
3221  const4=const4+1.0/yy(i)
3222  581 CONTINUE
3223  const4=const4/float(nb3)
3224 
3225  const5=0.0
3226  DO 591 i=ist4,ied4
3227  const5=const5+1.0/yy(i)
3228  591 CONTINUE
3229  const5=const5/float(nb4)
3230 
3231  const6=0.0
3232  DO 596 i=ist1p14,ied1p14
3233  const6=const6+1.0/yy(i)
3234  596 CONTINUE
3235  const6=const6/float(nb1p14)
3236 
3237  r114c=const6/((wt3*const4) + (wt4*const5))
3238  ENDIF
3239 
3240  CALL hunt(r0p94,nh2o,r094c,ja)
3241  IF (ja.GT.0.AND.ja.LT.60) THEN
3242  dlta =r0p94(ja+1)-r0p94(ja)
3243  fja =(r0p94(ja+1)-r094c)/dlta
3244  fjap1 =(r094c-r0p94(ja))/dlta
3245 
3246  vaptta=fja*vaptot(ja)+fjap1*vaptot(ja+1)
3247  IF(r094co.GT.1.) vaptta=-vaptta
3248  ELSE
3249  IF(ja.LE.0) vaptta = vaptot(ja+1)
3250  IF(ja.GE.60) vaptta = vaptot(ja)
3251  ENDIF
3252 
3253  IF(r094co.LE.1.)THEN
3254  DO 900 i=1,nobs
3255  IF (ja.GT.0.AND.ja.LT.60) THEN
3256  speca(i)=fja*trntbl(i,ja)+fjap1*trntbl(i,ja+1)
3257  ELSE
3258  IF(ja.LE.0) speca(i)=trntbl(i,ja+1)
3259  IF(ja.GE.60) speca(i)=trntbl(i,ja)
3260  ENDIF
3261  900 CONTINUE
3262  ENDIF
3263 
3264  IF(r094co.GT.1.)THEN
3265  DO 901 i=1,nobs
3266  IF (ja.GT.0.AND.ja.LT.60) THEN
3267  speca(i)=1.0/(fja*trntbl(i,ja)+fjap1*trntbl(i,ja+1))
3268  ELSE
3269  IF(ja.LE.0) speca(i)=1.0/trntbl(i,ja+1)
3270  IF(ja.GE.60) speca(i)=1.0/trntbl(i,ja)
3271  ENDIF
3272  901 CONTINUE
3273  ENDIF
3274 
3275  jb = ja
3276  CALL hunt(r1p14,nh2o,r114c,jb)
3277  IF (jb.GT.0.AND.jb.LT.60) THEN
3278  dltb =r1p14(jb+1)-r1p14(jb)
3279  fjb =(r1p14(jb+1)-r114c)/dltb
3280  fjbp1 =(r114c-r1p14(jb))/dltb
3281 
3282  vapttb=fjb*vaptot(jb)+fjbp1*vaptot(jb+1)
3283  IF(r114co.GT.1.) vapttb=-vapttb
3284  ELSE
3285  IF(jb.LE.0) vapttb = vaptot(jb+1)
3286  IF(jb.GE.60) vapttb = vaptot(jb)
3287  ENDIF
3288 C
3289  IF(r114co.LE.1.)THEN
3290  DO 910 i=1,nobs
3291  IF (jb.GT.0.AND.jb.LT.60) THEN
3292  specb(i)=fjb*trntbl(i,jb)+fjbp1*trntbl(i,jb+1)
3293  ELSE
3294  IF(jb.LE.0) specb(i)=trntbl(i,jb+1)
3295  IF(jb.GE.60) specb(i)=trntbl(i,jb)
3296  ENDIF
3297  910 CONTINUE
3298  ENDIF
3299 C
3300  IF(r114co.GT.1.)THEN
3301  DO 911 i=1,nobs
3302  IF (jb.GT.0.AND.jb.LT.60) THEN
3303  specb(i)=1.0/(fjb*trntbl(i,jb)+fjbp1*trntbl(i,jb+1))
3304  ELSE
3305  IF(jb.LE.0) specb(i)=1.0/trntbl(i,jb+1)
3306  IF(jb.GE.60) specb(i)=1.0/trntbl(i,jb)
3307  ENDIF
3308  911 CONTINUE
3309  ENDIF
3310 
3311 C--- CLMWVP=0.5*(VAPTTA+VAPTTB)/GGEOM
3312  clmwvp=0.5*(vaptta+vapttb)/g_vap_equiv
3313 C
3314 C Derivation of surface reflectances
3315  DO 915 i=1,nobs
3316  specav(i)=0.5*(speca(i)+specb(i))
3317 C--- TRNTMP=(YY(I)/SPECAV(I))-ROTOT(I)
3318 C--- YY(I)=TRNTMP/(TTOT(I)+(STOT(I)*TRNTMP))
3319  yy(i) = yy(i) / specav(i)
3320  915 CONTINUE
3321 C
3322 C-- print*, 'i, rotot(i), ttot(i), stot(i)'
3323 C-- DO I = 1, NOBS
3324 C-- print*, i, rotot(i), ttot(i), stot(i)
3325 C-- END DO
3326 
3327 C
3328 C Smooth the derived surface reflectance spectra
3329 
3330  IF(dlt2.GT.dlt) THEN
3331 C
3332 C First, replace radiances <= 0 near the spectral overlapping parts of the
3333 C four AVIRIS spectrometers by radiances in the nearby AVIRIS' channels.
3334  DO 920 i=natot-2,natot+2
3335  IF(yy(i).LE.0.0) yy(i)=yy(i-1)
3336  920 CONTINUE
3337 
3338  DO 921 i=nbtot-2,nbtot+2
3339  IF(yy(i).LE.0.0) yy(i)=yy(i-1)
3340  921 CONTINUE
3341 
3342  DO 922 i=nctot-4,nctot+4
3343  IF(yy(i).LE.0.0) yy(i)=yy(i-1)
3344  922 CONTINUE
3345 
3346  DO 923 i=ndtot-7,nobs
3347  IF(yy(i).LE.0.0) yy(i)=yy(i-1)
3348  923 CONTINUE
3349 
3350  DO 3480 i=istrt2,iend2
3351  trncv(i)=0.0
3352  ii=i-ncv2
3353  DO 3490 j=i-ncv2,i+ncv2
3354  trncv(i)=trncv(i)+yy(j)*finst2(j-ii+1)
3355  3490 CONTINUE
3356  3480 CONTINUE
3357  DO 3495 i=istrt2,iend2
3358  yy(i)=trncv(i)
3359  3495 CONTINUE
3360  ENDIF
3361 
3362  RETURN
3363  END
3364 
3365 ********************************************************************************
3366 * *
3367 * Name: LOCATE *
3368 * Purpose: given an array XX of length N, and given a value X, returns a value*
3369 * J such that X is between XX(J) and XX(J+1). XX must be monotonic, *
3370 * Parameters: XX - monotonic array of values *
3371 * N - number of elements in XX *
3372 * X - value that will be matched to the XX array *
3373 * J - index into the XX array where XX(J) <= X <= XX(J+1) *
3374 * Algorithm: bisectional table searching, copied from Numerical Recipes. *
3375 * Globals used: none *
3376 * Global output: none *
3377 * Return Codes: J=0 or J=N is returned to indicate that X is out of range *
3378 * Special Considerations: none *
3379 * *
3380 ********************************************************************************
3381  SUBROUTINE locate(xx,n,x,j)
3382  INTEGER j,n
3383  REAL x,xx(n)
3384  INTEGER jl,jm,ju
3385  jl=0
3386  ju=n+1
3387 10 if(ju-jl.gt.1)then
3388  jm=(ju+jl)/2
3389  if((xx(n).ge.xx(1)).eqv.(x.ge.xx(jm)))then
3390  jl=jm
3391  else
3392  ju=jm
3393  endif
3394  goto 10
3395  endif
3396  if(x.eq.xx(1))then
3397  j=1
3398  else if(x.eq.xx(n))then
3399  j=n-1
3400  else
3401  j=jl
3402  endif
3403  return
3404  END
3405 
3406 ********************************************************************************
3407 * *
3408 * Name: CUBSPLN *
3409 * Purpose: an interface for performing cubic spline interpolations. *
3410 * Parameters: N - number of elements in XORGN, YORGN *
3411 * XORGN - original x values *
3412 * YORGN - original y values *
3413 * XINT - interpolated x values *
3414 * YINT - interpolated y values *
3415 * Algorithm: Straight forward calculations *
3416 * Globals used: NOBS - number of spectral points. *
3417 * Global output: none *
3418 * Return Codes: none *
3419 * Special Considerations: none *
3420 * *
3421 ********************************************************************************
3422 
3423  SUBROUTINE cubspln(N,XORGN,YORGN,XINT,YINT)
3424 
3425  dimension xorgn(1050),yorgn(1050),y2(1050)
3426  dimension xint(1024),yint(1024)
3427  INTEGER N !number of elements in XORGN,YORGN
3428 
3429  COMMON /getinput5/ nobs,hsurf,dlt,dlt2
3430 C
3431 C YP1 and YP2 are two parameters specifying the values of 2nd derivatives at
3432 C XORGN(1) and XORGN(N) and were used in cubic spline interpolations. The
3433 C setting of both YP1 and YP2 to 2.0E30 is equivalent to set the 2nd
3434 C derivatives at the two boundaries as zero, according to Numeric Recipes.
3435  yp1=2.0e30
3436  ypn=2.0e30
3437 
3438  CALL spline(xorgn,yorgn,n,yp1,ypn,y2)
3439 
3440  DO 320 i=1,nobs
3441  x=xint(i)
3442  CALL splint(xorgn,yorgn,y2,n,x,y)
3443  IF(y.LT.0.0) y=0.0
3444  yint(i)=y
3445 C
3446  320 CONTINUE
3447 
3448  RETURN
3449  END
3450 
3451 ********************************************************************************
3452 * *
3453 * Name: SPLINE *
3454 * Purpose: program for cubic spline interpolation --- copyed from Numerical *
3455 * Recipes, p.88-89 *
3456 * Parameters: X - x-values *
3457 * Y - y-values *
3458 * N - length of X *
3459 * YP1 - 2nd derivative at X(1) *
3460 * YPN - 2nd derivative at X(N) *
3461 * Y2 - 2nd derivatives at all X positions *
3462 * Algorithm: Given arrays X and Y of length N containing a tabulated function,*
3463 * i.e.,Yj=f(Xj), with X1 < X2...<XN, and given values YP1 and YPN for *
3464 * the first derivative of the interpolating function at points 1 *
3465 * and N, respectively, this routine returns an array Y2 of length N *
3466 * which contains the second derivatives of the interpolating function *
3467 * at the tabulated points Xj. If YP1 and/or YPN are equal to 1.0E30 or *
3468 * larger, the routine is signalled to set the corresponding boundary *
3469 * condition for a natural spline, with zero second derivative on *
3470 * that boundary. *
3471 * Globals used: none *
3472 * Global output: none *
3473 * Return Codes: none *
3474 * Special Considerations: none *
3475 * *
3476 ********************************************************************************
3477 
3478  subroutine spline(x,y,n,yp1,ypn,y2)
3479  parameter(nmax=1050)
3480  integer n,i,k
3481  real x(n),y(n),y2(n),u(nmax)
3482  real yp1,ypn,sig,p,qn,un
3483  if (yp1.gt..99e30) then
3484  y2(1)=0.
3485  u(1)=0.
3486  else
3487  y2(1)=-0.5
3488  u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
3489  endif
3490  do 11 i=2,n-1
3491  sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
3492  p=sig*y2(i-1)+2.
3493  y2(i)=(sig-1.)/p
3494  u(i)=(6.*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1))
3495  * /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
3496 11 continue
3497  if (ypn.gt..99e30) then
3498  qn=0.
3499  un=0.
3500  else
3501  qn=0.5
3502  un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
3503  endif
3504  y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
3505  do 12 k=n-1,1,-1
3506  y2(k)=y2(k)*y2(k+1)+u(k)
3507 12 continue
3508  return
3509  end
3510 
3511 ********************************************************************************
3512 * *
3513 * Name: SPLINT *
3514 * Purpose: calculate a cubic-spline interpolated Y value, -- copied from *
3515 * Numerical Recipes. *
3516 * Parameters: XA - original X array *
3517 * YA - original Y array *
3518 * Y2A - 2nd derivative array from SUBROUTINE SPLINE *
3519 * N - number of X elements *
3520 * X - an X position at which interpolation should be made *
3521 * Y - interpolated y value at position X *
3522 * Algorithm: Given the arrays XA and YA of length N, which tabulate a function*
3523 * (with the XAj's in order), and given the array Y2A, which is the output *
3524 * from SPLINE above, and given a value of X, this routine returns a *
3525 * cubic-spline interpolated value Y. *
3526 * Globals used: none *
3527 * Global output: none *
3528 * Return Codes: none *
3529 * Special Considerations: SPLINE is called only once to process an entire *
3530 * tabulated function in arrays Xi and Yi. Once this has been done, values *
3531 * of the interpolated function for any value of X are obtained by calls *
3532 * (as many as desired) to a separate routine SPLINT (for cubic spline *
3533 * interpolation") *
3534 * *
3535 ********************************************************************************
3536 
3537  subroutine splint(xa,ya,y2a,n,x,y)
3538  integer n,klo,khi,k
3539  real xa(n),ya(n),y2a(n)
3540  real x,y,h,a,b
3541  klo=1
3542  khi=n
3543 1 if (khi-klo.gt.1) then
3544  k=(khi+klo)/2
3545  if(xa(k).gt.x)then
3546  khi=k
3547  else
3548  klo=k
3549  endif
3550  goto 1
3551  endif
3552  h=xa(khi)-xa(klo)
3553  if (h.eq.0.) pause 'bad xa input.'
3554  a=(xa(khi)-x)/h
3555  b=(x-xa(klo))/h
3556  y=a*ya(klo)+b*ya(khi)+
3557  * ((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**2)/6.
3558  return
3559  end
3560 
3561 ********************************************************************************
3562 * *
3563 * Name: HUNT *
3564 * Purpose: finds the element in array XX that is closest to value X. Array AA*
3565 * must be monotonic, either increasing or decreasing. *
3566 * Parameters: XX - array to search *
3567 * N - number of elements in the array *
3568 * X - element to search for closest match *
3569 * JLO - index of the closest matching element *
3570 * Algorithm: this subroutine was copied from Numerical Recipes *
3571 * Globals used: none *
3572 * Global output: none *
3573 * Return Codes: none *
3574 * Special Considerations: none *
3575 * *
3576 ********************************************************************************
3577 
3578  SUBROUTINE hunt(xx,n,x,jlo)
3579  INTEGER jlo,n
3580  REAL x,xx(n)
3581  INTEGER inc,jhi,jm
3582  LOGICAL ascnd
3583  ascnd=xx(n).ge.xx(1)
3584  if(jlo.le.0.or.jlo.gt.n)then
3585  jlo=0
3586  jhi=n+1
3587  goto 3
3588  endif
3589  inc=1
3590  if(x.ge.xx(jlo).eqv.ascnd)then
3591 1 jhi=jlo+inc
3592  if(jhi.gt.n)then
3593  jhi=n+1
3594  else if(x.ge.xx(jhi).eqv.ascnd)then
3595  jlo=jhi
3596  inc=inc+inc
3597  goto 1
3598  endif
3599  else
3600  jhi=jlo
3601 2 jlo=jhi-inc
3602  if(jlo.lt.1)then
3603  jlo=0
3604  else if(x.lt.xx(jlo).eqv.ascnd)then
3605  jhi=jlo
3606  inc=inc+inc
3607  goto 2
3608  endif
3609  endif
3610 3 if(jhi-jlo.eq.1)then
3611  if(x.eq.xx(n))jlo=n-1
3612  if(x.eq.xx(1))jlo=1
3613  return
3614  endif
3615  jm=(jhi+jlo)/2
3616  if(x.ge.xx(jm).eqv.ascnd)then
3617  jlo=jm
3618  else
3619  jhi=jm
3620  endif
3621  goto 3
3622  END
3623 
3624 ********************************************************************************
3625 * *
3626 * Name: SUNCOR *
3627 * Purpose: computes at some reference time TZ the declination of the sun and *
3628 * the hour angle-TZ. *
3629 * Algorithm: At any time T, the local solar hour angle is HAZ+T(GMT)-XLONG, *
3630 * where XLONG is the longitude measured positive west of Greenwich. *
3631 * Both time and angles are in radians. To compute azimuth and *
3632 * elevation at latitude XLAT, longitude XLONG, and time T, call *
3633 * HAZEL(HAZ+T-XLONG,DEC,AZ,EL,XLAT). This routine was copied *
3634 * from W. Mankin at NCAR, Boulder, CO. *
3635 * Globals used: none *
3636 * Global output: none *
3637 * Return Codes: none *
3638 * Special Considerations: none *
3639 * *
3640 ********************************************************************************
3641 
3642  SUBROUTINE suncor(IDAY,MONTH,IYR,TZ,DEC,HAZ)
3643 
3644  jd=julian(iday,month,iyr)
3645  fjd=0.5+tz/6.283185307
3646  CALL solcor(jd,fjd,ras,dec,gsdt,bzero,pzero,solong)
3647  haz=gsdt-ras-tz
3648 
3649  RETURN
3650  END
3651 
3652 ********************************************************************************
3653 * Name: SOLCOR *
3654 * Purpose: A routine for solar geometry calculations --- copied from *
3655 * W. Mankin at NCAR. *
3656 * Parameters: none *
3657 * Algorithm: This routine was supplied by W. Mankin at NCAR, Boulder, CO *
3658 * Globals used: *
3659 * Global output: *
3660 * Return Codes: none *
3661 * Special Considerations: none *
3662 * *
3663 ********************************************************************************
3664  SUBROUTINE solcor(JD,FJD,RAS,DECS,GSDT,BZRO,P,SOLONG)
3665 
3666  pi=3.141592654
3667  d=(jd-2415020)+fjd
3668  iyr=d/365.25
3669  g=-.026601523+.01720196977*d-1.95e-15*d*d -2.*pi*iyr
3670  xlms=4.881627938+.017202791266*d+3.95e-15*d*d-2.*pi*iyr
3671  obl=.409319747-6.2179e-9*d
3672  ecc=.01675104-1.1444e-9*d
3673  f=d-365.25*iyr
3674  gsdt=1.739935476+2.*pi*f/365.25+1.342027e-4*d/365.25
3675  gsdt=gsdt+6.2831853*(fjd-0.5)
3676  xlts=xlms+2.*ecc*sin(g)+1.25*ecc*ecc*sin(2.*g)
3677  sndc=sin(xlts)*sin(obl)
3678  decs=asin(sndc)
3679  csra=cos(xlts)/cos(decs)
3680  ras=acos(csra)
3681  IF(sin(xlts).LT.0.) ras=2.*pi-ras
3682  omega=1.297906+6.66992e-7*d
3683  thetac=xlts-omega
3684  bzro=asin(.126199*sin(thetac))
3685  p=-atan(cos(xlts)*tan(obl))-atan(.127216*cos(thetac))
3686  xlmm=atan2(.992005*sin(thetac),cos(thetac))
3687  jdr=jd-2398220
3688  irot=(jdr+fjd)/25.38
3689  frot=(jdr+fjd)/25.38-irot
3690  solong=xlmm-2.*pi*frot+pi-3.e-4
3691  IF(solong.LT.0.) solong=solong+2.*pi
3692  IF(solong.GT.(2.*pi)) solong=solong-2.*pi
3693 
3694  RETURN
3695  END
3696 
3697 
3698 ********************************************************************************
3699 * *
3700 * Name: JULIAN *
3701 * Purpose: computes the julian day *
3702 * Parameters: IDAY - day of the year *
3703 * MONTH - day of the year *
3704 * Algorithm: This routine was supplied by William Mankin at NCAR, Boulder, CO *
3705 * Globals used: none *
3706 * Global output: none *
3707 * Return Codes: none *
3708 * Special Considerations: none *
3709 * *
3710 ********************************************************************************
3711 
3712  FUNCTION julian(IDAY,MONTH,IYR)
3713 
3714 C Local variables
3715  dimension md(12)
3716 
3717  DATA md/0,31,59,90,120,151,181,212,243,273,304,334/
3718 C
3719  IF(iyr.LT.100) iyr=iyr+1900
3720  jyr=iyr-1600
3721  i1=jyr/400
3722  i2=(jyr-400*i1)/100
3723  i3=(jyr-400*i1-100*i2)/4
3724  julian=2305447+365*jyr+97*i1+24*i2+i3
3725  leap=0
3726  IF(mod(jyr,4).EQ.0) leap=1
3727  IF(mod(jyr,100).EQ.0) leap=0
3728  IF(mod(jyr,400).EQ.0) leap=1
3729 
3730 C LEAP=1 if iyr is a leap year
3731  julian=julian+md(month)+iday
3732  IF(month.LE.2) julian=julian-leap
3733 
3734  RETURN
3735  END
3736 
3737 ********************************************************************************
3738 * *
3739 * Name: HAZEL *
3740 * Purpose: Calculates azimuth and elevation *
3741 * Parameters: none *
3742 * Algorithm: This routine was supplied by William Mankin at NCAR, Boulder, CO *
3743 * Globals used: None. *
3744 * Global output: None. *
3745 * Return Codes: None *
3746 * Special Considerations: None. *
3747 * *
3748 ********************************************************************************
3749  SUBROUTINE hazel (H,D,A,E,XLAT)
3750 
3751 C H = HOUR ANGLE
3752 C D = DECLINATION
3753 C A = AZIMUTH
3754 C E = ELEVATION
3755 C XLAT = LATITUDE
3756 C ALL ANGLES IN RADIANS
3757 
3758  pi = 3.14159265
3759  sne = sin(d)*sin(xlat)+cos(d)*cos(xlat)*cos(h)
3760  e=asin(sne)
3761  sna = cos(d)*sin(h)
3762  csa=(sin(xlat)*cos(h)*cos(d)-sin(d)*cos(xlat))
3763  a=atan2(sna,csa)+pi
3764  END
3765 
3766 ********************************************************************************
3767 * *
3768 * Name: FINDMATCH *
3769 * Purpose: finds the closest match for ELEM in LIST *
3770 * Parameters: LIST - array of values to match ELEM to. Elements is array *
3771 * should increase in value. *
3772 * Algorithm: linearly compare ELEM to each element. When ELEM is smaller *
3773 * than the LIST(I), then it is assumed to be closest to LIST(I-1) *
3774 * or LIST(I). The one that has the smallest absolute difference *
3775 * to ELEM is returned as the closest match. *
3776 * Globals used: none *
3777 * Global output: none *
3778 * Return Codes: the closest matching element index *
3779 * Special Considerations: none *
3780 * *
3781 ********************************************************************************
3782 
3783  FUNCTION findmatch(LIST,NOBS,ELEM)
3784 
3785  dimension list(1500)
3786  INTEGER nobs
3787  REAL elem,list
3788 
3789  DO 460 i=1,nobs
3790  IF(list(i).GT.elem) GOTO 470
3791  460 CONTINUE
3792  470 CONTINUE
3793  diff1=abs(list(i-1)-elem)
3794  diff2=abs(list(i)-elem)
3795  IF (diff1.LT.diff2) THEN
3796  findmatch=i-1
3797  ELSE
3798  findmatch=i
3799  ENDIF
3800 
3801  RETURN
3802  END
3803 C--------1---------2---------3---------4---------5---------6---------7--
app_refl_with_gas_removal_drv
subroutine app_refl_with_gas_removal_drv(CLMWVP)
Definition: atrem_app_refl_plus_gas_removal_2013.f:2885
list
list(APPEND LIBS ${PGSTK_LIBRARIES}) add_executable(atteph_info_modis atteph_info_modis.c) target_link_libraries(atteph_info_modis $
Definition: CMakeLists.txt:7
cubeio::closeoutfile
subroutine closeoutfile(LUN_OUT)
Definition: cubeio.f90:122
suncor
subroutine suncor(IDAY, MONTH, IYR, TZ, DEC, HAZ)
Definition: atrem_app_refl_f90_cubeio.f:3375
cubeio::wt_line
subroutine wt_line(LUN_VAP, NSAMPS, H2OBUF)
Definition: cubeio.f90:100
tran_table
subroutine tran_table
Definition: atrem_app_refl_f90_cubeio.f:2155
model_adj
subroutine model_adj
Definition: atrem_app_refl_f90_cubeio.f:1074
hunt
int32_t hunt(double *xx, int32_t n, double x, int32_t jlo)
spline
subroutine spline(s, x, y, n, in, t, il, iu, vl, vu, e, u)
Definition: phs.f:1348
get_input
subroutine get_input
Definition: atrem_app_refl_f90_cubeio.f:395
cubspln
subroutine cubspln(N, XORGN, YORGN, XINT, YINT)
Definition: atrem_app_refl_f90_cubeio.f:3156
process_cube
subroutine process_cube
Definition: atrem_app_refl_f90_cubeio.f:2682
real
#define real
Definition: DbAlgOcean.cpp:26
bzero
void bzero()
cubeio::rd_slice
subroutine rd_slice(LUN_IN, NSAMPS, NBANDS, SORDER, BUFFER)
Definition: cubeio.f90:56
trancal
subroutine trancal
Definition: atrem_app_refl_f90_cubeio.f:2254
cubeio::wt_slice
subroutine wt_slice(LUN_OUT, NSAMPS, NBANDS, SORDER, BUFFER)
Definition: cubeio.f90:77
parameter
README for MOD_PR03(V6.1.0) 2. POINTS OF CONTACT it can be either SDP Toolkit or MODIS Packet for Terra input files The orbit validation configuration parameter(LUN 600281) must be either "TRUE" or "FALSE". It needs to be "FALSE" when running in Near Real Time mode
init_speccal
subroutine init_speccal
Definition: atrem_app_refl_f90_cubeio.f:1589
cubeio::closevapfile
subroutine closevapfile(LUN_VAP)
Definition: cubeio.f90:129
cubeio::closeinfile
subroutine closeinfile(LUN_IN)
Definition: cubeio.f90:114
re
#define re
Definition: l1_czcs_hdf.c:701
findmatch
function findmatch(LIST, NOBS, ELEM)
Definition: atrem_app_refl_f90_cubeio.f:3516
chnlratio
subroutine chnlratio
Definition: atrem_app_refl_f90_cubeio.f:2583
jd
Definition: jd.py:1
tt
subroutine tt(NMAX, NCHECK)
Definition: ampld.lp.f:1852
tran_smooth
subroutine tran_smooth
Definition: atrem_app_refl_f90_cubeio.f:2356
locate
subroutine locate(xx, n, x, j)
Definition: atrem_app_refl_f90_cubeio.f:3114
cubeio
Definition: cubeio.f90:1
solar_irr_pc
void solar_irr_pc()
Definition: solar_irr_PC.f:36
names
Definition: names.f90:1
geometry
subroutine geometry
Definition: atrem_app_refl_f90_cubeio.f:1329
splint
subroutine splint(xa, ya, y2a, n, x, y)
Definition: atrem_app_refl_f90_cubeio.f:3270
hazel
subroutine hazel(H, D, A, E, XLAT)
Definition: atrem_app_refl_f90_cubeio.f:3482
f
#define f
Definition: l1_czcs_hdf.c:702
leap
logical function leap(YEAR)
Definition: leap.f:10
julian
integer function julian(DAY, MONTH, YEAR)
Definition: julian.f:2
abs
#define abs(a)
Definition: misc.h:90
solcor
subroutine solcor(JD, FJD, RAS, DECS, GSDT, BZRO, P, SOLONG)
Definition: atrem_app_refl_f90_cubeio.f:3397


ABOUT

Our Team

Contact Us

DATA

Get Started

Find Data

Use Data

TOOLS

SeaBASS

SeaDAS

Giovanni Data Visualizer

Worldview Mapping Interface

Mission Quality Monitor

HELP

Level 1 & 2 Browser Help

Level 3 & 4 Browser Help

File Search Help

Search & Download Help

Earthdata Forum

Join Mailing List

How to Cite Data

CONNECT WITH US
COMPLIANCE

Web Privacy Policy

Data & Information Policy

Communications Policy

Accessibility Policy

Freedom of Information Act

USA.gov

Responsible NASA Official: Sean Bailey
Curator: OceanColor Webmaster