OB.DAAC Logo
NASA Logo
Ocean Color Science Software

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


ABOUT

Our Team

Contact Us

DATA

Get Started

Find Data

Use Data

TOOLS

SeaBASS

SeaDAS

Giovanni Data Visualizer

Worldview Mapping Interface

Mission Quality Monitor

HELP

Level 1 & 2 Browser Help

Level 3 & 4 Browser Help

File Search Help

Search & Download Help

Earthdata Forum

Join Mailing List

How to Cite Data

CONNECT WITH US
COMPLIANCE

Web Privacy Policy

Data & Information Policy

Communications Policy

Accessibility Policy

Freedom of Information Act

USA.gov

Responsible NASA Official: Sean Bailey
Curator: OceanColor Webmaster