OB.DAAC Logo
NASA Logo
Ocean Color Science Software

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

Level 3 Time Series Plotter

Request a SeaHawk Image

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