OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
viirs_aerosol_luts_nc4.f95
Go to the documentation of this file.
2 
3 private
4 
7 
8 public :: aero_412, aero_470, aero_650
9 public :: aero_412_abs, aero_470_abs
12 public :: new_intep, polint
13 
14 integer, parameter :: SGL = selected_real_kind(p=3)
15 
16 real, dimension(30) :: lut_raa
17 real, dimension(46) :: lut_vza
18 
19 real, dimension(12,46,30,10,8,20) :: nvalx412
20 real, dimension(12,46,30,10,4,24) :: nvalx470
21 real, dimension(12,46,30,10,24) :: nvalx650
22 
23 type :: viirs_aerosol_lut
24 
25  integer :: nsza
26  integer :: nvza
27  integer :: nraa
28  integer :: naot
29  integer :: nssa
30  integer :: nsfc
31 
32  real, dimension(:), allocatable :: sza
33  real, dimension(:), allocatable :: vza
34  real, dimension(:), allocatable :: raa
35  real, dimension(:), allocatable :: aot
36  real, dimension(:), allocatable :: ssa
37  real, dimension(:), allocatable :: sfc
38 
39  real, dimension(:,:,:,:,:,:), allocatable :: nvalx
40 end type viirs_aerosol_lut
41 
42 type(viirs_aerosol_lut) :: default_lut412, default_lut488, default_lut672
43 type(viirs_aerosol_lut) :: dust_lut412, dust_lut488, dust_lut672
44 
45 common /aottbl/ nvalx412, nvalx(12,46,30,10), &
46 & theta0(12), theta(46), &
47 & phi(30), sfc_ref412(20), tau(10), w0(8), w0_470(4), &
48 & nvalx650, sfc_ref650(24), sfc_ref470(24), &
49 & sfcprs(360,720), nvalx470
50 
51 contains
52 
53 integer function load_viirs_aerosol_luts(lut_file) result(status)
54  implicit none
55 
56  character(len=255), intent(in) :: lut_file
57  character(len=255) :: lut_type
58 
59 ! print *, 'Reading land LUT...'
60  lut_type = 'LAND_AEROSOL_FINE'
61  status = read_aerosol_lut_file(lut_file, lut_type, default_lut412, default_lut488, default_lut672)
62  if (status /= 0) then
63  print *, 'ERROR: Failed to read in default land aerosol model LUT: ', status
64  print *, 'File: ', trim(lut_file)
65  return
66  end if
67 ! print *, 'done.'
68 
69 ! print *, 'Reading dust LUT...'
70  lut_type = 'LAND_AEROSOL_DUST'
71  status = read_aerosol_lut_file(lut_file, lut_type, dust_lut412, dust_lut488, dust_lut672)
72  if (status /= 0) then
73  print *, 'ERROR: Failed to read in dust aerosol model LUT: ', status
74  print *, 'File: ', trim(lut_file)
75  return
76  end if
77 ! print *, 'done.'
78 
79 ! -- translate new LUT node arrays to old variabled in common block
80 ! -- copy over the data to the deprecated common block
81 ! -- still used in the vegetated retrieval -- see find_v_vegset.f
82  nvalx412 = default_lut412%nvalx(:,:,:,:,:,:)
83  nvalx470 = default_lut488%nvalx(:,:,:,:,:,:)
84  nvalx650 = default_lut672%nvalx(:,:,:,:,1,:)
85 
86  tau = default_lut412%aot
87  theta0 = default_lut412%sza
88  phi = default_lut412%raa
89  theta = default_lut412%vza
90  w0 = default_lut412%ssa
91 
92  w0_470 = default_lut488%ssa
93  sfc_ref412 = default_lut412%sfc
94  sfc_ref470 = default_lut488%sfc
95  sfc_ref650 = default_lut672%sfc
96 
97  return
98 
99 end function load_viirs_aerosol_luts
100 
101 subroutine unload_viirs_aerosol_luts(status)
102  implicit none
103 
104  integer,intent(inout) :: status
105 
106  status = 0
107 ! deallocate(nvalx412, nvalx470, nvalx650, stat=status)
108 ! if (status /= 0) then
109 ! print *, "ERROR: Failed to deallocate aerosol LUT arrays: ", status
110 ! return
111 ! end if
112 
113 end subroutine unload_viirs_aerosol_luts
114 
115 integer function read_aerosol_lut_file(lut_file, type, lut412, lut488, lut672) result(status)
116 
117 ! include 'hdf.f90'
118 ! include 'dffunc.f90'
119  use netcdf
120 
121  implicit none
122 
123  character(len=255), intent(in) :: lut_file
124  character(len=255), intent(in) :: type
125  type(viirs_aerosol_lut), intent(inout) :: lut412
126  type(viirs_aerosol_lut), intent(inout) :: lut488
127  type(viirs_aerosol_lut), intent(inout) :: lut672
128 
129  integer :: sd_id, sds_index, sds_id, rank, n_attrs, data_type
130  integer, dimension(1) :: start1, stride1, edges1, dim_sizes1
131  integer, dimension(5) :: edges5, start5, stride5, dim_sizes5
132  integer, dimension(6) :: edges6, start6, stride6, dim_sizes6
133  integer, dimension(32):: dimids
134  integer :: i, dimlen
135  integer :: xtype, ndims
136  integer :: nc_id
137  integer :: dim_id
138  integer :: dset_id
139  integer :: grp_id
140 
141  character(len=255) :: dset_name
142  character(len=255) :: attr_name
143  character(len=255) :: group_name
144  character(len=255) :: err_msg
145 
146  status = nf90_open(lut_file, nf90_nowrite, nc_id)
147  if (status /= nf90_noerr) then
148  print *, "ERROR: Failed to open deepblue lut_nc4 file: ", status
149  return
150  end if
151 
152  group_name = trim(type)
153  status = nf90_inq_ncid(nc_id, group_name, grp_id)
154  if (status /= nf90_noerr) then
155  print *, "ERROR: Failed to get ID of group "//trim(group_name)//": ", status
156  return
157  end if
158 
159 ! -- 412nm
160 
161  dset_name = 'SZA412_Nodes'
162  status = nf90_inq_varid(grp_id, dset_name, dset_id)
163  if (status /= nf90_noerr) then
164  print *, "ERROR: Failed to get ID of dataset "//trim(dset_name)//": ", status
165  return
166  end if
167  status = nf90_inquire_variable(grp_id, dset_id, dimids=dimids)
168  status = nf90_inquire_dimension(grp_id, dimids(1), len = dimlen)
169  err_msg = nf90_strerror(status)
170  if (status /= nf90_noerr) then
171  print *, "ERROR: Unable to get info on SDS "//trim(dset_name)//": ", status
172  status = -1
173  return
174  end if
175 
176  lut412%nsza = dimlen
177  allocate(lut412%sza(lut412%nsza), stat=status)
178  if (status /= 0) then
179  print *, "ERROR: Unable to allocate SZA412 data array: ", status
180  return
181  end if
182 
183  start1 = (/1/)
184  stride1 = (/1/)
185  edges1 = (/lut412%nsza/)
186  status = nf90_get_var(grp_id, dset_id, lut412%sza, start=start1, &
187  stride=stride1, count=edges1)
188  err_msg = nf90_strerror(status)
189  if (status /= nf90_noerr) then
190  print *, "ERROR: Failed to read dataset "//trim(dset_name)//": ", status
191  return
192  end if
193 
194  dset_name = 'VZA412_Nodes'
195  status = nf90_inq_varid(grp_id, dset_name, dset_id)
196  if (status /= nf90_noerr) then
197  print *, "ERROR: Failed to get ID of dataset "//trim(dset_name)//": ", status
198  return
199  end if
200  status = nf90_inquire_variable(grp_id, dset_id, dimids=dimids)
201  status = nf90_inquire_dimension(grp_id, dimids(1), len = dimlen)
202  err_msg = nf90_strerror(status)
203  if (status /= nf90_noerr) then
204  print *, "ERROR: Unable to get info on SDS "//trim(dset_name)//": ", status
205  status = -1
206  return
207  end if
208 
209  lut412%nvza = dimlen
210  allocate(lut412%vza(lut412%nvza), stat=status)
211  if (status /= 0) then
212  print *, "ERROR: Unable to allocate VZA412 data array: ", status
213  return
214  end if
215 
216  start1 = (/1/)
217  stride1 = (/1/)
218  edges1 = (/lut412%nvza/)
219  status = nf90_get_var(grp_id, dset_id, lut412%vza, start=start1, &
220  stride=stride1, count=edges1)
221  err_msg = nf90_strerror(status)
222  if (status /= nf90_noerr) then
223  print *, "ERROR: Failed to read dataset "//trim(dset_name)//": ", status
224  return
225  end if
226 
227  dset_name = 'RAA412_Nodes'
228  status = nf90_inq_varid(grp_id, dset_name, dset_id)
229  if (status /= nf90_noerr) then
230  print *, "ERROR: Failed to get ID of dataset "//trim(dset_name)//": ", status
231  return
232  end if
233  status = nf90_inquire_variable(grp_id, dset_id, dimids=dimids)
234  status = nf90_inquire_dimension(grp_id, dimids(1), len = dimlen)
235  err_msg = nf90_strerror(status)
236  if (status /= nf90_noerr) then
237  print *, "ERROR: Unable to get info on SDS "//trim(dset_name)//": ", status
238  status = -1
239  return
240  end if
241 
242  lut412%nraa = dimlen
243  allocate(lut412%raa(lut412%nraa), stat=status)
244  if (status /= 0) then
245  print *, "ERROR: Unable to allocate RAA412 data array: ", status
246  return
247  end if
248 
249  start1 = (/1/)
250  stride1 = (/1/)
251  edges1 = (/lut412%nraa/)
252  status = nf90_get_var(grp_id, dset_id, lut412%raa, start=start1, &
253  stride=stride1, count=edges1)
254  if (status /= nf90_noerr) then
255  print *, "ERROR: Failed to read dataset "//trim(dset_name)//": ", status
256  return
257  end if
258 
259  dset_name = 'AOT412_Nodes'
260  status = nf90_inq_varid(grp_id, dset_name, dset_id)
261  if (status /= nf90_noerr) then
262  print *, "ERROR: Failed to get ID of dataset "//trim(dset_name)//": ", status
263  return
264  end if
265  status = nf90_inquire_variable(grp_id, dset_id, dimids=dimids)
266  status = nf90_inquire_dimension(grp_id, dimids(1), len = dimlen)
267  err_msg = nf90_strerror(status)
268  if (status /= nf90_noerr) then
269  print *, "ERROR: Unable to get info on SDS "//trim(dset_name)//": ", status
270  status = -1
271  return
272  end if
273 
274  lut412%naot = dimlen
275  allocate(lut412%aot(lut412%naot), stat=status)
276  if (status /= 0) then
277  print *, "ERROR: Unable to allocate AOT412 data array: ", status
278  return
279  end if
280 
281  start1 = (/1/)
282  stride1 = (/1/)
283  edges1 = (/lut412%naot/)
284  status = nf90_get_var(grp_id, dset_id, lut412%aot, start=start1, &
285  stride=stride1, count=edges1)
286  err_msg = nf90_strerror(status)
287  if (status /= nf90_noerr) then
288  print *, "ERROR: Failed to read dataset "//trim(dset_name)//": ", status
289  return
290  end if
291 
292  dset_name = 'SSA412_Nodes'
293  status = nf90_inq_varid(grp_id, dset_name, dset_id)
294  if (status /= nf90_noerr) then
295  print *, "ERROR: Failed to get ID of dataset "//trim(dset_name)//": ", status
296  return
297  end if
298  status = nf90_inquire_variable(grp_id, dset_id, dimids=dimids)
299  status = nf90_inquire_dimension(grp_id, dimids(1), len = dimlen)
300  err_msg = nf90_strerror(status)
301  if (status /= nf90_noerr) then
302  print *, "ERROR: Unable to get info on SDS "//trim(dset_name)//": ", status
303  status = -1
304  return
305  end if
306 
307  lut412%nssa = dimlen
308  allocate(lut412%ssa (lut412%nssa ), stat=status)
309  if (status /= 0) then
310  print *, "ERROR: Unable to allocate SSA412 data array: ", status
311  return
312  end if
313 
314  start1 = (/1/)
315  stride1 = (/1/)
316  edges1 = (/lut412%nssa /)
317  status = nf90_get_var(grp_id, dset_id, lut412%ssa , start=start1, &
318  stride=stride1, count=edges1)
319  err_msg = nf90_strerror(status)
320  if (status /= nf90_noerr) then
321  print *, "ERROR: Failed to read dataset "//trim(dset_name)//": ", status
322  return
323  end if
324 
325  dset_name = 'SR412_Nodes'
326  status = nf90_inq_varid(grp_id, dset_name, dset_id)
327  if (status /= nf90_noerr) then
328  print *, "ERROR: Failed to get ID of dataset "//trim(dset_name)//": ", status
329  return
330  end if
331  status = nf90_inquire_variable(grp_id, dset_id, dimids=dimids)
332  status = nf90_inquire_dimension(grp_id, dimids(1), len = dimlen)
333  err_msg = nf90_strerror(status)
334  if (status /= nf90_noerr) then
335  print *, "ERROR: Unable to get info on SDS "//trim(dset_name)//": ", status
336  status = -1
337  return
338  end if
339 
340  lut412%nsfc = dimlen
341  allocate(lut412%sfc (lut412%nsfc ), stat=status)
342  if (status /= 0) then
343  print *, "ERROR: Unable to allocate SR412 data array: ", status
344  return
345  end if
346 
347  start1 = (/1/)
348  stride1 = (/1/)
349  edges1 = (/lut412%nsfc /)
350  status = nf90_get_var(grp_id, dset_id, lut412%sfc , start=start1, &
351  stride=stride1, count=edges1)
352  err_msg = nf90_strerror(status)
353  if (status /= nf90_noerr) then
354  print *, "ERROR: Failed to read dataset "//trim(dset_name)//": ", status
355  return
356  end if
357 
358  allocate(lut412%nvalx(lut412%nsza,lut412%nvza,lut412%nraa,lut412%naot,lut412%nssa, &
359  & lut412%nsfc), stat=status)
360  if (status /= 0) then
361  print *, "ERROR: Unable to allocate nvalx412 data array: ", status
362  return
363  end if
364 
365  dset_name = 'NVALX412'
366  status = nf90_inq_varid(grp_id, dset_name, dset_id)
367  if (status /= nf90_noerr) then
368  print *, "ERROR: Failed to get ID of dataset "//trim(dset_name)//": ", status
369  return
370  end if
371  status = nf90_inquire_variable(grp_id, dset_id, dimids=dimids)
372  status = nf90_inquire_dimension(grp_id, dimids(1), len = dimlen)
373  err_msg = nf90_strerror(status)
374  if (status /= nf90_noerr) then
375  print *, "ERROR: Unable to get info on SDS "//trim(dset_name)//": ", status
376  status = -1
377  return
378  end if
379 
380  lut412%nvalx = dimlen
381  if (status /= 0) then
382  print *, "ERROR: Unable to allocate SR412 data array: ", status
383  return
384  end if
385 
386  start6 = (/1,1,1,1,1,1/)
387  stride6 = (/1,1,1,1,1,1/)
388  edges6 = shape(lut412%nvalx)
389  status = nf90_get_var(grp_id, dset_id, lut412%nvalx , start=start6, &
390  stride=stride6, count=edges6)
391  err_msg = nf90_strerror(status)
392  if (status /= nf90_noerr) then
393  print *, "ERROR: Failed to read dataset "//trim(dset_name)//": ", status
394  return
395  end if
396 
397 
398 ! -- 488nm
399 
400  dset_name = 'SZA488_Nodes'
401  status = nf90_inq_varid(grp_id, dset_name, dset_id)
402  if (status /= nf90_noerr) then
403  print *, "ERROR: Failed to get ID of dataset "//trim(dset_name)//": ", status
404  return
405  end if
406  status = nf90_inquire_variable(grp_id, dset_id, dimids=dimids)
407  status = nf90_inquire_dimension(grp_id, dimids(1), len = dimlen)
408  err_msg = nf90_strerror(status)
409  if (status /= nf90_noerr) then
410  print *, "ERROR: Unable to get info on SDS "//trim(dset_name)//": ", status
411  status = -1
412  return
413  end if
414 
415  lut488%nsza = dimlen
416  allocate(lut488%sza(lut488%nsza), stat=status)
417  if (status /= 0) then
418  print *, "ERROR: Unable to allocate SZA488 data array: ", status
419  return
420  end if
421 
422  start1 = (/1/)
423  stride1 = (/1/)
424  edges1 = (/lut488%nsza/)
425  status = nf90_get_var(grp_id, dset_id, lut488%sza, start=start1, &
426  stride=stride1, count=edges1)
427  err_msg = nf90_strerror(status)
428  if (status /= nf90_noerr) then
429  print *, "ERROR: Failed to read dataset "//trim(dset_name)//": ", status
430  return
431  end if
432 
433  dset_name = 'VZA488_Nodes'
434  status = nf90_inq_varid(grp_id, dset_name, dset_id)
435  if (status /= nf90_noerr) then
436  print *, "ERROR: Failed to get ID of dataset "//trim(dset_name)//": ", status
437  return
438  end if
439  status = nf90_inquire_variable(grp_id, dset_id, dimids=dimids)
440  status = nf90_inquire_dimension(grp_id, dimids(1), len = dimlen)
441  err_msg = nf90_strerror(status)
442  if (status /= nf90_noerr) then
443  print *, "ERROR: Unable to get info on SDS "//trim(dset_name)//": ", status
444  status = -1
445  return
446  end if
447 
448  lut488%nvza = dimlen
449  allocate(lut488%vza(lut488%nvza), stat=status)
450  if (status /= 0) then
451  print *, "ERROR: Unable to allocate VZA488 data array: ", status
452  return
453  end if
454 
455  start1 = (/1/)
456  stride1 = (/1/)
457  edges1 = (/lut488%nvza/)
458  status = nf90_get_var(grp_id, dset_id, lut488%vza, start=start1, &
459  stride=stride1, count=edges1)
460  err_msg = nf90_strerror(status)
461  if (status /= nf90_noerr) then
462  print *, "ERROR: Failed to read dataset "//trim(dset_name)//": ", status
463  return
464  end if
465 
466  dset_name = 'RAA488_Nodes'
467  status = nf90_inq_varid(grp_id, dset_name, dset_id)
468  if (status /= nf90_noerr) then
469  print *, "ERROR: Failed to get ID of dataset "//trim(dset_name)//": ", status
470  return
471  end if
472  status = nf90_inquire_variable(grp_id, dset_id, dimids=dimids)
473  status = nf90_inquire_dimension(grp_id, dimids(1), len = dimlen)
474  err_msg = nf90_strerror(status)
475  if (status /= nf90_noerr) then
476  print *, "ERROR: Unable to get info on SDS "//trim(dset_name)//": ", status
477  status = -1
478  return
479  end if
480 
481  lut488%nraa = dimlen
482  allocate(lut488%raa(lut488%nraa), stat=status)
483  if (status /= 0) then
484  print *, "ERROR: Unable to allocate RAA488 data array: ", status
485  return
486  end if
487 
488  start1 = (/1/)
489  stride1 = (/1/)
490  edges1 = (/lut488%nraa/)
491  status = nf90_get_var(grp_id, dset_id, lut488%raa, start=start1, &
492  stride=stride1, count=edges1)
493  err_msg = nf90_strerror(status)
494  if (status /= nf90_noerr) then
495  print *, "ERROR: Failed to read dataset "//trim(dset_name)//": ", status
496  return
497  end if
498 
499  dset_name = 'AOT488_Nodes'
500  status = nf90_inq_varid(grp_id, dset_name, dset_id)
501  if (status /= nf90_noerr) then
502  print *, "ERROR: Failed to get ID of dataset "//trim(dset_name)//": ", status
503  return
504  end if
505  status = nf90_inquire_variable(grp_id, dset_id, dimids=dimids)
506  status = nf90_inquire_dimension(grp_id, dimids(1), len = dimlen)
507  err_msg = nf90_strerror(status)
508  if (status /= nf90_noerr) then
509  print *, "ERROR: Unable to get info on SDS "//trim(dset_name)//": ", status
510  status = -1
511  return
512  end if
513 
514  lut488%naot = dimlen
515  allocate(lut488%aot(lut488%naot), stat=status)
516  if (status /= 0) then
517  print *, "ERROR: Unable to allocate AOT488 data array: ", status
518  return
519  end if
520 
521  start1 = (/1/)
522  stride1 = (/1/)
523  edges1 = (/lut488%naot/)
524  status = nf90_get_var(grp_id, dset_id, lut488%aot, start=start1, &
525  stride=stride1, count=edges1)
526  err_msg = nf90_strerror(status)
527  if (status /= nf90_noerr) then
528  print *, "ERROR: Failed to read dataset "//trim(dset_name)//": ", status
529  return
530  end if
531 
532  dset_name = 'SSA488_Nodes'
533  status = nf90_inq_varid(grp_id, dset_name, dset_id)
534  if (status /= nf90_noerr) then
535  print *, "ERROR: Failed to get ID of dataset "//trim(dset_name)//": ", status
536  return
537  end if
538  status = nf90_inquire_variable(grp_id, dset_id, dimids=dimids)
539  status = nf90_inquire_dimension(grp_id, dimids(1), len = dimlen)
540  err_msg = nf90_strerror(status)
541  if (status /= nf90_noerr) then
542  print *, "ERROR: Unable to get info on SDS "//trim(dset_name)//": ", status
543  status = -1
544  return
545  end if
546 
547  lut488%nssa = dimlen
548  allocate(lut488%ssa(lut488%nssa ), stat=status)
549  if (status /= 0) then
550  print *, "ERROR: Unable to allocate SSA488 data array: ", status
551  return
552  end if
553 
554  start1 = (/1/)
555  stride1 = (/1/)
556  edges1 = (/lut488%nssa /)
557  status = nf90_get_var(grp_id, dset_id, lut488%ssa , start=start1, &
558  stride=stride1, count=edges1)
559  err_msg = nf90_strerror(status)
560  if (status /= nf90_noerr) then
561  print *, "ERROR: Failed to read dataset "//trim(dset_name)//": ", status
562  return
563  end if
564 
565  dset_name = 'SR488_Nodes'
566  status = nf90_inq_varid(grp_id, dset_name, dset_id)
567  if (status /= nf90_noerr) then
568  print *, "ERROR: Failed to get ID of dataset "//trim(dset_name)//": ", status
569  return
570  end if
571  status = nf90_inquire_variable(grp_id, dset_id, dimids=dimids)
572  status = nf90_inquire_dimension(grp_id, dimids(1), len = dimlen)
573  err_msg = nf90_strerror(status)
574  if (status /= nf90_noerr) then
575  print *, "ERROR: Unable to get info on SDS "//trim(dset_name)//": ", status
576  status = -1
577  return
578  end if
579 
580  lut488%nsfc = dimlen
581  allocate(lut488%sfc (lut488%nsfc ), stat=status)
582  if (status /= 0) then
583  print *, "ERROR: Unable to allocate SR488 data array: ", status
584  return
585  end if
586 
587  start1 = (/1/)
588  stride1 = (/1/)
589  edges1 = (/lut488%nsfc /)
590  status = nf90_get_var(grp_id, dset_id, lut488%sfc , start=start1, &
591  stride=stride1, count=edges1)
592  if (status /= nf90_noerr) then
593  print *, "ERROR: Failed to read dataset "//trim(dset_name)//": ", status
594  return
595  end if
596 
597  allocate(lut488%nvalx(lut488%nsza,lut488%nvza,lut488%nraa,lut488%naot, &
598  lut488%nssa, lut488%nsfc), stat=status)
599  if (status /= 0) then
600  print *, "ERROR: Unable to allocate nvalx488 data array: ", status
601  return
602  end if
603 
604  dset_name = 'NVALX488'
605  status = nf90_inq_varid(grp_id, dset_name, dset_id)
606  if (status /= nf90_noerr) then
607  print *, "ERROR: Failed to get ID of dataset "//trim(dset_name)//": ", status
608  return
609  end if
610  status = nf90_inquire_variable(grp_id, dset_id, dimids=dimids)
611  status = nf90_inquire_dimension(grp_id, dimids(1), len = dimlen)
612  err_msg = nf90_strerror(status)
613  if (status /= nf90_noerr) then
614  print *, "ERROR: Unable to get info on SDS "//trim(dset_name)//": ", status
615  status = -1
616  return
617  end if
618 
619  lut488%nvalx = dimlen
620  if (status /= 0) then
621  print *, "ERROR: Unable to allocate SR488 data array: ", status
622  return
623  end if
624 
625  start6 = (/1,1,1,1,1,1/)
626  stride6 = (/1,1,1,1,1,1/)
627  edges6 = shape(lut488%nvalx)
628  status = nf90_get_var(grp_id, dset_id, lut488%nvalx , start=start6, &
629  stride=stride6, count=edges6)
630  err_msg = nf90_strerror(status)
631  if (status /= nf90_noerr) then
632  print *, "ERROR: Failed to read dataset "//trim(dset_name)//": ", status
633  return
634  end if
635 
636 
637 ! -- 672nm
638 
639  dset_name = 'SZA672_Nodes'
640  status = nf90_inq_varid(grp_id, dset_name, dset_id)
641  if (status /= nf90_noerr) then
642  print *, "ERROR: Failed to get ID of dataset "//trim(dset_name)//": ", status
643  return
644  end if
645  status = nf90_inquire_variable(grp_id, dset_id, dimids=dimids)
646  status = nf90_inquire_dimension(grp_id, dimids(1), len = dimlen)
647  err_msg = nf90_strerror(status)
648  if (status /= nf90_noerr) then
649  print *, "ERROR: Unable to get info on SDS "//trim(dset_name)//": ", status
650  status = -1
651  return
652  end if
653 
654  lut672%nsza = dimlen
655  allocate(lut672%sza(lut672%nsza), stat=status)
656  if (status /= 0) then
657  print *, "ERROR: Unable to allocate SZA672 data array: ", status
658  return
659  end if
660 
661  start1 = (/1/)
662  stride1 = (/1/)
663  edges1 = (/lut672%nsza/)
664  status = nf90_get_var(grp_id, dset_id, lut672%sza, start=start1, &
665  stride=stride1, count=edges1)
666  err_msg = nf90_strerror(status)
667  if (status /= nf90_noerr) then
668  print *, "ERROR: Failed to read dataset "//trim(dset_name)//": ", status
669  return
670  end if
671 
672  dset_name = 'VZA672_Nodes'
673  status = nf90_inq_varid(grp_id, dset_name, dset_id)
674  if (status /= nf90_noerr) then
675  print *, "ERROR: Failed to get ID of dataset "//trim(dset_name)//": ", status
676  return
677  end if
678  status = nf90_inquire_variable(grp_id, dset_id, dimids=dimids)
679  status = nf90_inquire_dimension(grp_id, dimids(1), len = dimlen)
680  err_msg = nf90_strerror(status)
681  if (status /= nf90_noerr) then
682  print *, "ERROR: Unable to get info on SDS "//trim(dset_name)//": ", status
683  status = -1
684  return
685  end if
686 
687  lut672%nvza = dimlen
688  allocate(lut672%vza(lut672%nvza), stat=status)
689  if (status /= 0) then
690  print *, "ERROR: Unable to allocate VZA672 data array: ", status
691  return
692  end if
693 
694  start1 = (/1/)
695  stride1 = (/1/)
696  edges1 = (/lut672%nvza/)
697  status = nf90_get_var(grp_id, dset_id, lut672%vza, start=start1, &
698  stride=stride1, count=edges1)
699  err_msg = nf90_strerror(status)
700  if (status /= nf90_noerr) then
701  print *, "ERROR: Failed to read dataset "//trim(dset_name)//": ", status
702  return
703  end if
704 
705  dset_name = 'RAA672_Nodes'
706  status = nf90_inq_varid(grp_id, dset_name, dset_id)
707  if (status /= nf90_noerr) then
708  print *, "ERROR: Failed to get ID of dataset "//trim(dset_name)//": ", status
709  return
710  end if
711  status = nf90_inquire_variable(grp_id, dset_id, dimids=dimids)
712  status = nf90_inquire_dimension(grp_id, dimids(1), len = dimlen)
713  err_msg = nf90_strerror(status)
714  if (status /= nf90_noerr) then
715  print *, "ERROR: Unable to get info on SDS "//trim(dset_name)//": ", status
716  status = -1
717  return
718  end if
719 
720  lut672%nraa = dimlen
721  allocate(lut672%raa(lut672%nraa), stat=status)
722  if (status /= 0) then
723  print *, "ERROR: Unable to allocate RAA672 data array: ", status
724  return
725  end if
726 
727  start1 = (/1/)
728  stride1 = (/1/)
729  edges1 = (/lut672%nraa/)
730  status = nf90_get_var(grp_id, dset_id, lut672%raa, start=start1, &
731  stride=stride1, count=edges1)
732  err_msg = nf90_strerror(status)
733  if (status /= nf90_noerr) then
734  print *, "ERROR: Failed to read dataset "//trim(dset_name)//": ", status
735  return
736  end if
737 
738  dset_name = 'AOT672_Nodes'
739  status = nf90_inq_varid(grp_id, dset_name, dset_id)
740  if (status /= nf90_noerr) then
741  print *, "ERROR: Failed to get ID of dataset "//trim(dset_name)//": ", status
742  return
743  end if
744  status = nf90_inquire_variable(grp_id, dset_id, dimids=dimids)
745  status = nf90_inquire_dimension(grp_id, dimids(1), len = dimlen)
746  err_msg = nf90_strerror(status)
747  if (status /= nf90_noerr) then
748  print *, "ERROR: Unable to get info on SDS "//trim(dset_name)//": ", status
749  status = -1
750  return
751  end if
752 
753  lut672%naot = dimlen
754  allocate(lut672%aot(lut672%naot), stat=status)
755  if (status /= 0) then
756  print *, "ERROR: Unable to allocate AOT672 data array: ", status
757  return
758  end if
759 
760  start1 = (/1/)
761  stride1 = (/1/)
762  edges1 = (/lut672%naot/)
763  status = nf90_get_var(grp_id, dset_id, lut672%aot, start=start1, &
764  stride=stride1, count=edges1)
765  err_msg = nf90_strerror(status)
766  if (status /= nf90_noerr) then
767  print *, "ERROR: Failed to read dataset "//trim(dset_name)//": ", status
768  return
769  end if
770 
771  dset_name = 'SSA672_Nodes'
772  status = nf90_inq_varid(grp_id, dset_name, dset_id)
773  if (status /= nf90_noerr) then
774  print *, "ERROR: Failed to get ID of dataset "//trim(dset_name)//": ", status
775  return
776  end if
777  status = nf90_inquire_variable(grp_id, dset_id, dimids=dimids)
778  status = nf90_inquire_dimension(grp_id, dimids(1), len = dimlen)
779  err_msg = nf90_strerror(status)
780  if (status /= nf90_noerr) then
781  print *, "ERROR: Unable to get info on SDS "//trim(dset_name)//": ", status
782  status = -1
783  return
784  end if
785 
786  lut672%nssa = dimlen
787  allocate(lut672%ssa(lut672%nssa ), stat=status)
788  if (status /= 0) then
789  print *, "ERROR: Unable to allocate SSA672 data array: ", status
790  return
791  end if
792 
793  start1 = (/1/)
794  stride1 = (/1/)
795  edges1 = (/lut672%nssa /)
796  status = nf90_get_var(grp_id, dset_id, lut672%ssa , start=start1, &
797  stride=stride1, count=edges1)
798  err_msg = nf90_strerror(status)
799  if (status /= nf90_noerr) then
800  print *, "ERROR: Failed to read dataset "//trim(dset_name)//": ", status
801  return
802  end if
803 
804  dset_name = 'SR672_Nodes'
805  status = nf90_inq_varid(grp_id, dset_name, dset_id)
806  if (status /= nf90_noerr) then
807  print *, "ERROR: Failed to get ID of dataset "//trim(dset_name)//": ", status
808  return
809  end if
810  status = nf90_inquire_variable(grp_id, dset_id, dimids=dimids)
811  status = nf90_inquire_dimension(grp_id, dimids(1), len = dimlen)
812  err_msg = nf90_strerror(status)
813  if (status /= nf90_noerr) then
814  print *, "ERROR: Unable to get info on SDS "//trim(dset_name)//": ", status
815  status = -1
816  return
817  end if
818 
819  lut672%nsfc = dimlen
820  allocate(lut672%sfc(lut672%nsfc ), stat=status)
821  if (status /= 0) then
822  print *, "ERROR: Unable to allocate SR672 data array: ", status
823  return
824  end if
825 
826  start1 = (/1/)
827  stride1 = (/1/)
828  edges1 = (/lut672%nsfc /)
829  status = nf90_get_var(grp_id, dset_id, lut672%sfc , start=start1, &
830  stride=stride1, count=edges1)
831  err_msg = nf90_strerror(status)
832  if (status /= nf90_noerr) then
833  print *, "ERROR: Failed to read dataset "//trim(dset_name)//": ", status
834  return
835  end if
836 
837  allocate(lut672%nvalx(lut672%nsza,lut672%nvza,lut672%nraa,lut672%naot,lut672%nssa, &
838  lut672%nsfc), stat=status)
839  if (status /= 0) then
840  print *, "ERROR: Unable to allocate nvalx672 data array: ", status
841  return
842  end if
843 
844  dset_name = 'NVALX672'
845  status = nf90_inq_varid(grp_id, dset_name, dset_id)
846  if (status /= nf90_noerr) then
847  print *, "ERROR: Failed to get ID of dataset "//trim(dset_name)//": ", status
848  return
849  end if
850  status = nf90_inquire_variable(grp_id, dset_id, dimids=dimids)
851  status = nf90_inquire_dimension(grp_id, dimids(1), len = dimlen)
852  err_msg = nf90_strerror(status)
853  if (status /= nf90_noerr) then
854  print *, "ERROR: Unable to get info on SDS "//trim(dset_name)//": ", status
855  status = -1
856  return
857  end if
858 
859  lut672%nvalx = dimlen
860  if (status /= 0) then
861  print *, "ERROR: Unable to allocate SR672 data array: ", status
862  return
863  end if
864 
865  start6 = (/1,1,1,1,1,1/)
866  stride6 = (/1,1,1,1,1,1/)
867  edges6 = shape(lut672%nvalx)
868  status = nf90_get_var(grp_id, dset_id, lut672%nvalx , start=start6, &
869  stride=stride6, count=edges6)
870  err_msg = nf90_strerror(status)
871  if (status /= nf90_noerr) then
872  print *, "ERROR: Failed to read dataset "//trim(dset_name)//": ", status
873  return
874  end if
875 
876  status = nf90_close(nc_id)
877  if (status /= nf90_noerr) then
878  print *, "ERROR: Failed to close lut_nc4 file: ", status
879  return
880  end if
881 
882  return
883 end function read_aerosol_lut_file
884 
885 
886 !--------------------------------------------------------
887 !-- NOTE: if model_frac = 0.0, the aerosol model = imod
888 !-- if model_frac > 0.0, the aerosol model will be interpolated between
889 !-- imod and imod + 1, using the value of model_frac
890 !
891 subroutine aero_470(dflag, refl, x1, x2, x3, mm, nn, ll, ma, imod, &
892 & r470, tau_x470, tau_x470_flag, trflg, model_frac, debug)
893 
894  implicit none
895 
896  logical, intent(inout) :: dflag
897  real, intent(in) :: refl
898  real, intent(in) :: x1
899  real, intent(in) :: x2
900  real, intent(in) :: x3
901  integer, intent(in) :: mm
902  integer, intent(in) :: nn
903  integer, intent(in) :: ll
904  integer, intent(in) :: ma
905  integer, intent(in) :: imod
906  real, intent(in) :: r470
907  real, intent(inout) :: tau_x470
908  integer, intent(inout) :: tau_x470_flag
909  real, intent(inout) :: trflg
910  real, intent(in) :: model_frac
911  logical, intent(in) :: debug
912 
913  real, dimension(:), allocatable :: yy
914  real, dimension(8) :: yy2
915 
916  integer :: ii
917  real :: frac
918 
919  integer :: status
920 
921  tau_x470 = -999.0
922  tau_x470_flag = -999
923 
924  status = default_lut488%naot
925  status = 0
926  if (allocated(yy)) deallocate(yy)
927  allocate(yy(default_lut488%naot), stat=status)
928  if (status /= 0) then
929  print *, "ERROR: Failed to allocate array for reduced AOT 488 table: ", status
930  return
931  end if
932 
933  status = create_reduced_lut_aot(default_lut488, refl, x1,x2,x3, imod, &
934  & r470, model_frac, yy, debug)
935  if (status /= 0) then
936  deallocate(yy, stat=status)
937  dflag = .true.
938  return
939  end if
940 
941  if (refl <= yy(1)) then
942  tau_x470 = 0.06
943  if (trflg > 0.0) tau_x470 = 0.02
944  tau_x470_flag = -10
945  if (debug) print *, 'aero_470, hit low bound: ', refl, yy(1)
946  return
947  end if
948 
949 ! Reflc off the charts! Set AOT to max and set flag.
950  if (refl >= yy(size(yy))) then
951  tau_x470 = extrap(refl, yy, default_lut488%aot(1:default_lut488%naot), status)
952  if (status /= 0) then
953  if (status == 1) then
954  tau_x470 = 5.0
955  else
956  print *, 'ERROR: Failed to extrapolate AOT: ', status
957  return
958  end if
959  end if
960  if (tau_x470 > 5.0) tau_x470 = 5.0
961  tau_x470_flag = 1
962  if (debug) print *, 'aero_470, hit hi bound: ', refl, yy(10)
963  deallocate(yy, stat=status)
964  return
965  endif
966 
967 !
968 ! Check if the reflectance increase with AOT
969 !
970 
971  if (yy(1) < yy(2)) go to 650
972 
973  if (refl < yy(4)) return
974 
975  yy2(:) = -999.0
976  yy2(1:7) = yy(4:10) ! last cell is empty, how it was originally written!
977 
978  if (yy2(2) < yy2(1)) return
979 
980  ii = search(refl, yy2, status, frac=frac)
981  if (status /= 0) then
982  dflag = .true.
983  return
984  end if
985 
986  tau_x470 = frac*default_lut488%aot(ii+1+3) + (1.-frac)*default_lut488%aot(ii+3)
987  tau_x470_flag = 0
988 ! if (debug) print *, 'aero_470, exit 2358, aot: ', tau_x470
989  return
990 
991 650 continue
992 
993 !
994 ! Pass the monotonic order check
995 !
996  ii = search(refl, yy, status, frac=frac)
997  if (status /= 0) then
998  dflag = .true.
999  return
1000  end if
1001 
1002  tau_x470 = frac*default_lut488%aot(ii+1) + (1.-frac)*default_lut488%aot(ii)
1003  tau_x470_flag = 0
1004 ! if (debug) print *, 'aero_470, exit 2371, aot: ', tau_x470
1005 ! print *,'tau_x470 =', tau_x470
1006 
1007  deallocate(yy, stat=status)
1008  if (status /= 0) then
1009  print *, "WARNING: Failed to deallocate reduced AOT table: ", status
1010  end if
1011 
1012  return
1013 
1014 end subroutine aero_470
1015 
1016 
1017 !--------------------------------------------------------
1018 ! TODO: check/use/add status flag
1019 subroutine aero_650(dflag,refl,x1,x2,x3,mm,nn,ll,ma,r650,tau_x650, &
1020 & tau_x650_flag,tau_x470_flag,tau_x412,tau_x470,tau_x412_flag_91,trflg)
1021  implicit none
1022 
1023  logical, intent(inout) :: dflag
1024  real, intent(in) :: refl
1025  real, intent(in) :: x1
1026  real, intent(in) :: x2
1027  real, intent(in) :: x3
1028  integer, intent(in) :: mm
1029  integer, intent(in) :: nn
1030  integer, intent(in) :: ll
1031  integer, intent(in) :: ma
1032  real, intent(in) :: r650
1033  real, intent(inout) :: tau_x650
1034  integer, intent(inout) :: tau_x650_flag
1035  integer, intent(in) :: tau_x470_flag
1036  real, intent(in) :: tau_x412
1037  real, intent(in) :: tau_x470
1038  integer, intent(in) :: tau_x412_flag_91
1039  real, intent(in) :: trflg
1040  !logical, intent(in), optional :: debug
1041 
1042  real, dimension(:), allocatable :: yy
1043  real, dimension(8) :: yy2
1044  real, dimension(4) :: yy3
1045  real, dimension(6) :: yy5
1046 
1047  real :: frac
1048  real :: w0_x
1049  integer :: ii
1050 
1051  integer :: status
1052  logical :: debug
1053 
1054  tau_x650_flag = -999
1055  tau_x650 = -999.0
1056  debug = .false.
1057 
1058  if (allocated(yy)) deallocate(yy)
1059  allocate(yy(default_lut672%naot), stat=status)
1060  if (status /= 0) then
1061  print *, "ERROR: Failed to allocate array for reduced AOT 672 table: ", status
1062  return
1063  end if
1064 
1065  status = create_reduced_lut_aot(default_lut672, refl, x1,x2,x3, 1, &
1066  & r650, 1.0, yy, debug)
1067  if (status /= 0) then
1068  deallocate(yy, stat=status)
1069  dflag = .true.
1070  return
1071  end if
1072 
1073  if (refl <= yy(1) .and. yy(1) < yy(2)) then
1074  tau_x650 = 0.002
1075  if (trflg > 0.0) tau_x650 = 0.02
1076  tau_x650_flag = -10
1077  return
1078  end if
1079 
1080 ! Reflc off the charts! Set AOT to max and set flag.
1081  if (refl >= yy(10)) then
1082  tau_x650 = extrap(refl, yy, default_lut672%aot(1:default_lut672%naot), status)
1083  if (status /= 0) then
1084  if (status == 1) then
1085  tau_x650 = 5.0
1086  else
1087  print *, 'ERROR: Failed to extrapolate AOT: ', status
1088  return
1089  end if
1090  end if
1091  if (tau_x650 > 5.0) tau_x650 = 5.0
1092  w0_x = -999.
1093  tau_x650_flag = 1
1094  deallocate(yy, stat=status)
1095  return
1096  end if
1097 
1098 ! -- check for absorbing dust over bright surface
1099  if (refl >= yy(7)) then
1100  yy3(:) = yy(7:10)
1101 
1102  if (yy3(2) < yy3(1)) return
1103 
1104  ii = search(refl, yy3, status, frac=frac)
1105  if (status /= 0) then
1106  dflag = .true.
1107  return
1108  end if
1109 
1110  tau_x650 = frac*default_lut672%aot(ii+1+6) + (1.-frac)*default_lut672%aot(ii+6)
1111  tau_x650_flag = 0
1112  deallocate(yy, stat=status)
1113  return
1114 
1115  end if
1116 
1117 ! -- check if 470 or 412 retrievals were off the charts (refl > yy(10))
1118 ! if (tau_x470_flag > 0) go to 670
1119 ! if (tau_x412_flag_91 > 0) go to 680
1120 
1121 !
1122 ! Check if the reflectance increase with AOT
1123 !
1124  if (yy(1) < yy(2)) go to 650
1125 
1126  if (refl < yy(4)) return
1127 
1128  yy2(1:7) = yy(4:10)
1129 
1130  if (yy2(2) < yy2(1)) return
1131 
1132  ii = search(refl, yy2, status, frac=frac)
1133  if (status /= 0) then
1134  dflag = .true.
1135  return
1136  end if
1137 
1138  tau_x650 = frac*default_lut672%aot(ii+1+3) + (1.-frac)*default_lut672%aot(ii+3)
1139  tau_x650_flag = 0
1140  return
1141 
1142 670 continue
1143 ! if (refl < yy(8)) return
1144 !
1145 ! do i = 1, 3
1146 ! yy3(i) = yy(i+7)
1147 ! end do
1148 !
1149 ! if (yy3(2) < yy3(1)) return
1150 ! call search2(dflag,refl,yy3,3,ii,frac)
1151 ! if (dflag) return
1152 ! tau_x650 = frac*tau(ii+1+7) + (1.-frac)*tau(ii+7)
1153 ! tau_x650_flag = 0
1154 ! return
1155 
1156 680 continue
1157  if (refl < yy(5)) return
1158 
1159  yy5(:) = yy(5:10)
1160 
1161  if (yy5(2) < yy5(1)) return
1162 
1163  ii = search(refl, yy5, status, frac=frac)
1164  if (status /= 0) then
1165  dflag = .true.
1166  return
1167  end if
1168 
1169 ! if (yy3(1) > yy3(2)) print *,'yy=',refl,(yy(i),I=1,10)
1170  tau_x650 = frac*default_lut672%aot(ii+1+4) + (1.-frac)*default_lut672%aot(ii+4)
1171  tau_x650_flag = 0
1172 
1173  return
1174 
1175 650 continue
1176 !
1177 ! Pass the monotonic order check
1178 !
1179  ii = search(refl, yy, status, frac=frac)
1180  if (status /= 0) then
1181  dflag = .true.
1182  return
1183  end if
1184 
1185 ! print *,'after 2nd search'
1186  tau_x650 = frac*default_lut672%aot(ii+1) + (1.-frac)*default_lut672%aot(ii)
1187  tau_x650_flag = 0
1188 
1189  deallocate(yy, stat=status)
1190  if (status /= 0) then
1191  print *, "WARNING: Failed to deallocate reduced AOT table: ", status
1192  end if
1193 
1194  return
1195 
1196 end subroutine aero_650
1197 
1198 ! -- NOTE: if model_frac = 0.0, the aerosol model = imod
1199 ! -- if model_frac > 0.0, the aerosol model will be interpolated between
1200 ! -- imod and imod + 1, using the value of model_frac
1201 !
1202 !--------------------------------------------------------
1203  subroutine aero_412(dflag,refl,x1,x2,x3,mm,nn,ll,ma,imod,r412, &
1204  & tau_x412,tau_x412_flag,trflg,model_frac,debug)
1206  implicit none
1207 
1208 
1209  logical, intent(inout) :: dflag
1210  real, intent(in) :: refl
1211  real, intent(in) :: x1
1212  real, intent(in) :: x2
1213  real, intent(in) :: x3
1214  integer, intent(in) :: mm
1215  integer, intent(in) :: nn
1216  integer, intent(in) :: ll
1217  integer, intent(in) :: ma
1218  integer, intent(in) :: imod
1219  real, intent(in) :: r412
1220  real, intent(inout) :: tau_x412
1221  integer, intent(inout) :: tau_x412_flag
1222  real, intent(in) :: trflg
1223  real, intent(in) :: model_frac
1224  logical, intent(in) :: debug
1225 
1226  real, dimension(:), allocatable :: yy
1227  real, dimension(8) :: yy2
1228 
1229  integer :: ii
1230  real :: frac
1231 
1232  integer :: status
1233 
1234  tau_x412 = -999.0
1235  tau_x412_flag = -999
1236 
1237  if (allocated(yy)) deallocate(yy)
1238  allocate(yy(default_lut412%naot), stat=status)
1239  if (status /= 0) then
1240  print *, "ERROR: Failed to allocate array for reduced AOT 412 table: ", status
1241  return
1242  end if
1243 
1244  status = create_reduced_lut_aot(default_lut412, refl, x1,x2,x3, imod, &
1245  & r412, model_frac, yy, debug)
1246  if (status /= 0) then
1247  deallocate(yy, stat=status)
1248  dflag = .true.
1249  return
1250  end if
1251 
1252  if (refl <= yy(1)) then
1253  tau_x412 = 0.06
1254  if (trflg > 0.0) tau_x412 = 0.02
1255  tau_x412_flag = -10
1256  if (debug) print *, 'aero_412, hit low bound: ', refl, yy(1)
1257  return
1258  end if
1259 
1260 ! Reflc off the charts! Set AOT to max and set flag.
1261  if (refl >= yy(10)) then
1262  tau_x412 = extrap(refl, yy, default_lut412%aot(1:default_lut412%naot), status)
1263  if (status /= 0) then
1264  if (status == 1) then
1265  tau_x412 = 5.0
1266  else
1267  print *, 'ERROR: Failed to extrapolate AOT: ', status
1268  return
1269  end if
1270  end if
1271  if (tau_x412 > 5.0) tau_x412 = 5.0
1272 
1273  tau_x412_flag = 1
1274  if (debug) print *, 'aero_412, hit hi bound: ', refl, yy(10)
1275  deallocate(yy, stat=status)
1276  return
1277  end if
1278 !
1279 ! Check if the reflectance increase with AOT
1280 !
1281  if (yy(1) < yy(2)) go to 650
1282 
1283  if (refl < yy(4)) return
1284 
1285  yy2(:) = -999.0
1286  yy2(1:7) = yy(4:10)
1287 
1288  if (yy2(2) < yy2(1)) return
1289 
1290  ii = search(refl, yy2, status, frac=frac)
1291  if (status /= 0) then
1292  dflag = .true.
1293  return
1294  end if
1295 
1296  tau_x412 = frac*default_lut412%aot(ii+1+3) + (1.0-frac)*default_lut412%aot(ii+3)
1297  tau_x412_flag = 0
1298 ! if (debug) print *, 'aero_412, exit 2355, aot: ', tau_x412
1299  return
1300 
1301 650 continue
1302 !
1303 ! Pass the monotonic order check
1304 !
1305  ii = search(refl, yy, status, frac=frac)
1306  if (status /= 0) then
1307  dflag = .true.
1308  return
1309  end if
1310 
1311  tau_x412 = frac*default_lut412%aot(ii+1) + (1.-frac)*default_lut412%aot(ii)
1312  tau_x412_flag = 0
1313 ! if (debug) print *, 'aero_412, exit 2367, aot: ', tau_x412, refl
1314 
1315  deallocate(yy, stat=status)
1316  if (status /= 0) then
1317  print *, "WARNING: Failed to deallocate reduced AOT table: ", status
1318  end if
1319 
1320  return
1321 
1322 end subroutine aero_412
1323 
1324 !--------------------------------------------------------
1325 subroutine aero_412_abs(dflag,refl,x1,x2,x3,mm,nn,ll,r412,tau_x,w0_x)
1326  implicit none
1327 
1328  logical, intent(inout) :: dflag
1329  real, intent(in) :: refl
1330  real, intent(in) :: x1
1331  real, intent(in) :: x2
1332  real, intent(in) :: x3
1333  integer, intent(in) :: mm
1334  integer, intent(in) :: nn
1335  integer, intent(in) :: ll
1336  real, intent(in) :: r412
1337  real, intent(in) :: tau_x
1338  real, intent(inout) :: w0_x
1339 
1340  integer :: index_ii, index_ia
1341  real :: frac, frac_ia
1342  integer :: status
1343 
1344  real, dimension(:,:,:,:), allocatable :: nnvalxw
1345  real, dimension(:), allocatable :: yyw
1346 
1347  dflag = .false.
1348 
1349  if (allocated(nnvalxw)) deallocate(nnvalxw, yyw)
1350  allocate(nnvalxw(4,4,2,default_lut412%nssa), yyw(default_lut412%nssa), stat=status)
1351  if (status /= 0) then
1352  print *, "ERROR: Failed to allocate arrays: ", status
1353  dflag = .false.
1354  return
1355  end if
1356 
1357  index_ia = search(tau_x, default_lut412%aot, status, frac=frac_ia)
1358  if (status /= 0) then
1359  dflag = .false.
1360  return
1361  end if
1362 
1363  status = create_reduced_lut_ssa(default_lut412, refl, x1, x2, x3, index_ia, &
1364  & r412, frac_ia, yyw)
1365  if (status /= 0) then
1366  dflag = .true.
1367  return
1368  end if
1369 
1370  if (refl.le.yyw(1)) then
1371  w0_x = 0.82
1372  return
1373  endif
1374 
1375  if (refl.ge.yyw(8)) then
1376  w0_x = 1.0
1377  return
1378  endif
1379 
1380  index_ii = search(refl, yyw, status, frac=frac)
1381  if (status /= 0) then
1382  dflag = .false.
1383  return
1384  end if
1385  w0_x = frac*default_lut412%ssa(index_ii+1) + (1.-frac)*default_lut412%ssa(index_ii)
1386 
1387  deallocate(nnvalxw, yyw, stat=status)
1388  if (status /= 0) then
1389  print *, "WARNING: Failed to deallocate arrays: ", status
1390  return
1391  end if
1392 
1393  return
1394 
1395 end subroutine aero_412_abs
1396 
1397 !--------------------------------------------------------
1398 subroutine aero_470_abs(dflag2,refl,x1,x2,x3,mm,nn,ll,r470,tau_x,w0_x)
1399  implicit none
1400 
1401  logical, intent(inout) :: dflag2
1402  real, intent(in) :: refl
1403  real, intent(in) :: x1
1404  real, intent(in) :: x2
1405  real, intent(in) :: x3
1406  integer, intent(in) :: mm
1407  integer, intent(in) :: nn
1408  integer, intent(in) :: ll
1409  real, intent(in) :: r470
1410  real, intent(in) :: tau_x
1411  real, intent(inout) :: w0_x
1412 
1413  integer :: index_ii, index_ia
1414  real :: frac, frac_ia
1415  integer :: status
1416 
1417  real, dimension(:,:,:,:), allocatable :: nnvalxw
1418  real, dimension(:), allocatable :: yyw
1419 
1420  dflag2 = .false.
1421 
1422  if (allocated(nnvalxw)) deallocate(nnvalxw, yyw)
1423  allocate(nnvalxw(4,4,2,default_lut488%nssa), yyw(default_lut488%nssa), stat=status)
1424  if (status /= 0) then
1425  print *, "ERROR: Failed to allocate arrays: ", status
1426  dflag2 = .false.
1427  return
1428  end if
1429 
1430  index_ia = search(tau_x, default_lut488%aot, status, frac=frac_ia)
1431  if (status /= 0) then
1432  dflag2 = .false.
1433  return
1434  end if
1435 
1436  status = create_reduced_lut_ssa(default_lut488, refl, x1, x2, x3, index_ia, &
1437  & r470, frac_ia, yyw)
1438  if (status /= 0) then
1439  dflag2 = .true.
1440  return
1441  end if
1442 
1443  if (refl.le.yyw(1)) then
1444  w0_x = -999.
1445  return
1446  endif
1447 
1448  if (refl.ge.yyw(4)) then
1449  w0_x = 1.0
1450  return
1451  endif
1452 
1453  index_ii = search(refl, yyw, status, frac=frac)
1454  if (status /= 0) then
1455  dflag2 = .false.
1456  return
1457  end if
1458  w0_x = frac*default_lut488%ssa(index_ii+1) + (1.-frac)*default_lut488%ssa(index_ii)
1459 
1460  deallocate(nnvalxw, yyw, stat=status)
1461  if (status /= 0) then
1462  print *, "WARNING: Failed to deallocate arrays: ", status
1463  return
1464  end if
1465 
1466  return
1467 
1468 end subroutine aero_470_abs
1469 
1470 !--------------------------------------------------------
1471 !-- NOTE: if model_frac = 0.0, the aerosol model = imod
1472 !-- if model_frac > 0.0, the aerosol model will be interpolated between
1473 !-- imod and imod + 1, using the value of model_frac
1474 !
1475 subroutine aero_470_dust(dflag, refl, x1, x2, x3, mm, nn, ll, ma, imod, &
1476 & r470, tau_x470, tau_x470_flag, trflg, model_frac, debug)
1478  implicit none
1479 
1480  logical, intent(inout) :: dflag
1481  real, intent(in) :: refl
1482  real, intent(in) :: x1
1483  real, intent(in) :: x2
1484  real, intent(in) :: x3
1485  integer, intent(in) :: mm
1486  integer, intent(in) :: nn
1487  integer, intent(in) :: ll
1488  integer, intent(in) :: ma
1489  integer, intent(in) :: imod
1490  real, intent(in) :: r470
1491  real, intent(inout) :: tau_x470
1492  integer, intent(inout) :: tau_x470_flag
1493  real, intent(inout) :: trflg
1494  real, intent(in) :: model_frac
1495  logical, intent(in) :: debug
1496 
1497  real, dimension(:), allocatable :: yy
1498  real, dimension(8) :: yy2
1499 
1500  integer :: ii
1501  real :: frac
1502 
1503  integer :: status
1504 
1505  tau_x470 = -999.0
1506  tau_x470_flag = -999
1507 
1508  if (allocated(yy)) deallocate(yy)
1509  allocate(yy(dust_lut488%naot), stat=status)
1510  if (status /= 0) then
1511  print *, "ERROR: Failed to allocate array for reduced AOT DUST 488 table: ", status
1512  return
1513  end if
1514 
1515  status = create_reduced_lut_aot(dust_lut488, refl, x1,x2,x3, imod, &
1516  & r470, model_frac, yy, debug)
1517  if (status /= 0) then
1518  deallocate(yy, stat=status)
1519  dflag = .true.
1520  return
1521  end if
1522 
1523  if (refl <= yy(1)) then
1524  tau_x470 = 0.06
1525  if (trflg > 0.0) tau_x470 = 0.02
1526  tau_x470_flag = -10
1527  if (debug) print *, 'aero_470_dust, hit low bound: ', refl, yy(1)
1528  return
1529  end if
1530 
1531 ! Reflc off the charts! Set AOT to max and set flag.
1532  if (refl >= yy(size(yy))) then
1533  tau_x470 = extrap(refl, yy, dust_lut488%aot(1:dust_lut488%naot), status)
1534  if (status /= 0) then
1535  if (status == 1) then
1536  tau_x470 = 5.0
1537  else
1538  print *, 'ERROR: Failed to extrapolate AOT: ', status
1539  return
1540  end if
1541  end if
1542  if (tau_x470 > 5.0) tau_x470 = 5.0
1543  tau_x470_flag = 1
1544  if (debug) print *, 'aero_470_dust, hit hi bound: ', refl, yy(10)
1545  deallocate(yy, stat=status)
1546  return
1547  endif
1548 
1549 !
1550 ! Check if the reflectance increase with AOT
1551 !
1552 
1553  if (yy(1) < yy(2)) go to 650
1554 
1555  if (refl < yy(4)) return
1556 
1557  yy2(:) = -999.0
1558  yy2(1:7) = yy(4:10) ! last cell is empty, how it was originally written!
1559 
1560  if (yy2(2) < yy2(1)) return
1561 
1562  ii = search(refl, yy2, status, frac=frac)
1563  if (status /= 0) then
1564  dflag = .true.
1565  return
1566  end if
1567 
1568  tau_x470 = frac*dust_lut488%aot(ii+1+3) + (1.-frac)*dust_lut488%aot(ii+3)
1569  tau_x470_flag = 0
1570  if (debug) print *, 'aero_470_dust, exit 2358, aot: ', tau_x470
1571  return
1572 
1573 650 continue
1574 
1575 !
1576 ! Pass the monotonic order check
1577 !
1578  ii = search(refl, yy, status, frac=frac)
1579  if (status /= 0) then
1580  dflag = .true.
1581  return
1582  end if
1583 
1584  tau_x470 = frac*dust_lut488%aot(ii+1) + (1.-frac)*dust_lut488%aot(ii)
1585  tau_x470_flag = 0
1586  if (debug) print *, 'aero_470_dust, exit 2371, aot: ', tau_x470
1587 ! print *,'tau_x470 =', tau_x470
1588 
1589  deallocate(yy, stat=status)
1590  if (status /= 0) then
1591  print *, "WARNING: Failed to deallocate reduced AOT table: ", status
1592  end if
1593 
1594  return
1595 
1596 end subroutine aero_470_dust
1597 
1598 !--------------------------------------------------------
1599 ! TODO: check/use/add status flag
1600 subroutine aero_650_dust(dflag,refl,x1,x2,x3,mm,nn,ll,ma,r650,tau_x650, &
1601 & tau_x650_flag,tau_x470_flag,tau_x412,tau_x470,tau_x412_flag_91,trflg)
1602  implicit none
1603 
1604  logical, intent(inout) :: dflag
1605  real, intent(in) :: refl
1606  real, intent(in) :: x1
1607  real, intent(in) :: x2
1608  real, intent(in) :: x3
1609  integer, intent(in) :: mm
1610  integer, intent(in) :: nn
1611  integer, intent(in) :: ll
1612  integer, intent(in) :: ma
1613  real, intent(in) :: r650
1614  real, intent(inout) :: tau_x650
1615  integer, intent(inout) :: tau_x650_flag
1616  integer, intent(in) :: tau_x470_flag
1617  real, intent(in) :: tau_x412
1618  real, intent(in) :: tau_x470
1619  integer, intent(in) :: tau_x412_flag_91
1620  real, intent(in) :: trflg
1621  !logical, intent(in), optional :: debug
1622 
1623  real, dimension(:), allocatable :: yy
1624  real, dimension(8) :: yy2
1625  real, dimension(4) :: yy3
1626  real, dimension(6) :: yy5
1627 
1628  real :: frac
1629  real :: w0_x
1630  integer :: ii
1631 
1632  integer :: status
1633  logical :: debug
1634 
1635  tau_x650_flag = -999
1636  tau_x650 = -999.0
1637  debug = .false.
1638 
1639  if (allocated(yy)) deallocate(yy)
1640  allocate(yy(dust_lut672%naot), stat=status)
1641  if (status /= 0) then
1642  print *, "ERROR: Failed to allocate array for reduced AOT DUST 672 table: ", status
1643  return
1644  end if
1645 
1646  status = create_reduced_lut_aot(dust_lut672, refl, x1,x2,x3, 1, &
1647  & r650, 1.0, yy, debug)
1648  if (status /= 0) then
1649  deallocate(yy, stat=status)
1650  dflag = .true.
1651  return
1652  end if
1653 
1654  if (refl <= yy(1) .and. yy(1) < yy(2)) then
1655  tau_x650 = 0.002
1656  if (trflg > 0.0) tau_x650 = 0.02
1657  tau_x650_flag = -10
1658  return
1659  end if
1660 
1661 ! Reflc off the charts! Set AOT to max and set flag.
1662  if (refl >= yy(10)) then
1663  tau_x650 = extrap(refl, yy, dust_lut672%aot(1:dust_lut672%naot), status)
1664  if (status /= 0) then
1665  if (status == 1) then
1666  tau_x650 = 5.0
1667  else
1668  print *, 'ERROR: Failed to extrapolate AOT: ', status
1669  return
1670  end if
1671  end if
1672  if (tau_x650 > 5.0) tau_x650 = 5.0
1673  w0_x = -999.
1674  tau_x650_flag = 1
1675  deallocate(yy, stat=status)
1676  return
1677  end if
1678 
1679 ! -- check for absorbing dust over bright surface
1680  if (refl >= yy(7)) then
1681  yy3(:) = yy(7:10)
1682 
1683  if (yy3(2) < yy3(1)) return
1684 
1685  ii = search(refl, yy3, status, frac=frac)
1686  if (status /= 0) then
1687  dflag = .true.
1688  return
1689  end if
1690 
1691  tau_x650 = frac*dust_lut672%aot(ii+1+6) + (1.-frac)*dust_lut672%aot(ii+6)
1692  tau_x650_flag = 0
1693 
1694  return
1695 
1696  end if
1697 
1698 ! -- check if 470 or 412 retrievals were off the charts (refl > yy(10))
1699 ! if (tau_x470_flag > 0) go to 670
1700 ! if (tau_x412_flag_91 > 0) go to 680
1701 
1702 !
1703 ! Check if the reflectance increase with AOT
1704 !
1705  if (yy(1) < yy(2)) go to 650
1706 
1707  if (refl < yy(4)) return
1708 
1709  yy2(1:7) = yy(4:10)
1710 
1711  if (yy2(2) < yy2(1)) return
1712 
1713  ii = search(refl, yy2, status, frac=frac)
1714  if (status /= 0) then
1715  dflag = .true.
1716  return
1717  end if
1718 
1719  tau_x650 = frac*dust_lut672%aot(ii+1+3) + (1.-frac)*dust_lut672%aot(ii+3)
1720  tau_x650_flag = 0
1721  return
1722 
1723 670 continue
1724 ! if (refl < yy(8)) return
1725 !
1726 ! do i = 1, 3
1727 ! yy3(i) = yy(i+7)
1728 ! end do
1729 !
1730 ! if (yy3(2) < yy3(1)) return
1731 ! call search2(dflag,refl,yy3,3,ii,frac)
1732 ! if (dflag) return
1733 ! tau_x650 = frac*tau(ii+1+7) + (1.-frac)*tau(ii+7)
1734 ! tau_x650_flag = 0
1735 ! return
1736 
1737 680 continue
1738  if (refl < yy(5)) return
1739 
1740  yy5(:) = yy(5:10)
1741 
1742  if (yy5(2) < yy5(1)) return
1743 
1744  ii = search(refl, yy5, status, frac=frac)
1745  if (status /= 0) then
1746  dflag = .true.
1747  return
1748  end if
1749 
1750 ! if (yy3(1) > yy3(2)) print *,'yy=',refl,(yy(i),I=1,10)
1751  tau_x650 = frac*dust_lut672%aot(ii+1+4) + (1.-frac)*dust_lut672%aot(ii+4)
1752  tau_x650_flag = 0
1753 
1754  return
1755 
1756 650 continue
1757 !
1758 ! Pass the monotonic order check
1759 !
1760  ii = search(refl, yy, status, frac=frac)
1761  if (status /= 0) then
1762  dflag = .true.
1763  return
1764  end if
1765 
1766 ! print *,'after 2nd search'
1767  tau_x650 = frac*dust_lut672%aot(ii+1) + (1.-frac)*dust_lut672%aot(ii)
1768  tau_x650_flag = 0
1769 
1770  deallocate(yy, stat=status)
1771  if (status /= 0) then
1772  print *, "WARNING: Failed to deallocate reduced AOT table: ", status
1773  end if
1774 
1775  return
1776 
1777 end subroutine aero_650_dust
1778 
1779 ! -- NOTE: if model_frac = 0.0, the aerosol model = imod
1780 ! -- if model_frac > 0.0, the aerosol model will be interpolated between
1781 ! -- imod and imod + 1, using the value of model_frac
1782 !
1783 !--------------------------------------------------------
1784  subroutine aero_412_dust(dflag,refl,x1,x2,x3,mm,nn,ll,ma,imod,r412, &
1785  & tau_x412,tau_x412_flag,trflg,model_frac,debug)
1787  implicit none
1788 
1789 
1790  logical, intent(inout) :: dflag
1791  real, intent(in) :: refl
1792  real, intent(in) :: x1
1793  real, intent(in) :: x2
1794  real, intent(in) :: x3
1795  integer, intent(in) :: mm
1796  integer, intent(in) :: nn
1797  integer, intent(in) :: ll
1798  integer, intent(in) :: ma
1799  integer, intent(in) :: imod
1800  real, intent(in) :: r412
1801  real, intent(inout) :: tau_x412
1802  integer, intent(inout) :: tau_x412_flag
1803  real, intent(in) :: trflg
1804  real, intent(in) :: model_frac
1805  logical, intent(in) :: debug
1806 
1807  real, dimension(:), allocatable :: yy
1808  real, dimension(8) :: yy2
1809 
1810  integer :: ii
1811  real :: frac
1812 
1813  integer :: status
1814 
1815  tau_x412 = -999.0
1816  tau_x412_flag = -999
1817 
1818  if (allocated(yy)) deallocate(yy)
1819  allocate(yy(dust_lut412%naot), stat=status)
1820  if (status /= 0) then
1821  print *, "ERROR: Failed to allocate array for reduced AOT DUST 412 table: ", status
1822  return
1823  end if
1824 
1825  status = create_reduced_lut_aot(dust_lut412, refl, x1,x2,x3, imod, &
1826  & r412, model_frac, yy, debug)
1827  if (status /= 0) then
1828  deallocate(yy, stat=status)
1829  dflag = .true.
1830  return
1831  end if
1832 
1833  if (refl <= yy(1)) then
1834  tau_x412 = 0.06
1835  if (trflg > 0.0) tau_x412 = 0.02
1836  tau_x412_flag = -10
1837  if (debug) print *, 'aero_412_dust, hit low bound: ', refl, yy(1)
1838  return
1839  end if
1840 
1841 ! Reflc off the charts! Set AOT to max and set flag.
1842  if (refl >= yy(10)) then
1843  tau_x412 = extrap(refl, yy, dust_lut412%aot(1:dust_lut412%naot), status)
1844  if (status /= 0) then
1845  if (status == 1) then
1846  tau_x412 = 5.0
1847  else
1848  print *, 'ERROR: Failed to extrapolate AOT: ', status
1849  return
1850  end if
1851  end if
1852  if (tau_x412 > 5.0) tau_x412 = 5.0
1853 
1854  tau_x412_flag = 1
1855  if (debug) print *, 'aero_412_dust, hit hi bound: ', refl, yy(10)
1856  deallocate(yy, stat=status)
1857  return
1858  end if
1859 !
1860 ! Check if the reflectance increase with AOT
1861 !
1862  if (yy(1) < yy(2)) go to 650
1863 
1864  if (refl < yy(4)) return
1865 
1866  yy2(:) = -999.0
1867  yy2(1:7) = yy(4:10)
1868 
1869  if (yy2(2) < yy2(1)) return
1870 
1871  ii = search(refl, yy2, status, frac=frac)
1872  if (status /= 0) then
1873  dflag = .true.
1874  return
1875  end if
1876 
1877  tau_x412 = frac*dust_lut412%aot(ii+1+3) + (1.0-frac)*dust_lut412%aot(ii+3)
1878  tau_x412_flag = 0
1879  if (debug) print *, 'aero_412_dust, exit 2355, aot: ', tau_x412
1880  return
1881 
1882 650 continue
1883 !
1884 ! Pass the monotonic order check
1885 !
1886  ii = search(refl, yy, status, frac=frac)
1887  if (status /= 0) then
1888  dflag = .true.
1889  return
1890  end if
1891 
1892  tau_x412 = frac*dust_lut412%aot(ii+1) + (1.-frac)*dust_lut412%aot(ii)
1893  tau_x412_flag = 0
1894  if (debug) print *, 'aero_412_dust, exit 2367, aot: ', tau_x412, refl
1895 
1896  deallocate(yy, stat=status)
1897  if (status /= 0) then
1898  print *, "WARNING: Failed to deallocate reduced AOT table: ", status
1899  end if
1900 
1901  return
1902 
1903 end subroutine aero_412_dust
1904 
1905 !--------------------------------------------------------
1906 subroutine aero_412_abs_dust(dflag,refl,x1,x2,x3,mm,nn,ll,r412,tau_x,w0_x)
1907  implicit none
1908 
1909  logical, intent(inout) :: dflag
1910  real, intent(in) :: refl
1911  real, intent(in) :: x1
1912  real, intent(in) :: x2
1913  real, intent(in) :: x3
1914  integer, intent(in) :: mm
1915  integer, intent(in) :: nn
1916  integer, intent(in) :: ll
1917  real, intent(in) :: r412
1918  real, intent(in) :: tau_x
1919  real, intent(inout) :: w0_x
1920 
1921  integer :: index_ii, index_ia
1922  real :: frac, frac_ia
1923  integer :: status
1924 
1925  real, dimension(:,:,:,:), allocatable :: nnvalxw
1926  real, dimension(:), allocatable :: yyw
1927 
1928  dflag = .false.
1929 
1930  if (allocated(nnvalxw)) deallocate(nnvalxw, yyw)
1931  allocate(nnvalxw(4,4,2,dust_lut412%nssa), yyw(dust_lut412%nssa), stat=status)
1932  if (status /= 0) then
1933  print *, "ERROR: Failed to allocate arrays: ", status
1934  dflag = .false.
1935  return
1936  end if
1937 
1938  index_ia = search(tau_x, dust_lut412%aot, status, frac=frac_ia)
1939  if (status /= 0) then
1940  dflag = .false.
1941  return
1942  end if
1943 
1944  status = create_reduced_lut_ssa(dust_lut412, refl, x1, x2, x3, index_ia, &
1945  & r412, frac_ia, yyw)
1946  if (status /= 0) then
1947  dflag = .true.
1948  return
1949  end if
1950 
1951  if (refl.le.yyw(1)) then
1952  w0_x = 0.82
1953  return
1954  endif
1955 
1956  if (refl.ge.yyw(8)) then
1957  w0_x = 1.0
1958  return
1959  endif
1960 
1961  index_ii = search(refl, yyw, status, frac=frac)
1962  if (status /= 0) then
1963  dflag = .false.
1964  return
1965  end if
1966  w0_x = frac*dust_lut412%ssa(index_ii+1) + (1.-frac)*dust_lut412%ssa(index_ii)
1967 
1968  deallocate(nnvalxw, yyw, stat=status)
1969  if (status /= 0) then
1970  print *, "WARNING: Failed to deallocate arrays: ", status
1971  return
1972  end if
1973 
1974  return
1975 
1976 end subroutine aero_412_abs_dust
1977 
1978 !--------------------------------------------------------
1979 subroutine aero_470_abs_dust(dflag2,refl,x1,x2,x3,mm,nn,ll,r470,tau_x,w0_x)
1980  implicit none
1981 
1982  logical, intent(inout) :: dflag2
1983  real, intent(in) :: refl
1984  real, intent(in) :: x1
1985  real, intent(in) :: x2
1986  real, intent(in) :: x3
1987  integer, intent(in) :: mm
1988  integer, intent(in) :: nn
1989  integer, intent(in) :: ll
1990  real, intent(in) :: r470
1991  real, intent(in) :: tau_x
1992  real, intent(inout) :: w0_x
1993 
1994  integer :: index_ii, index_ia
1995  real :: frac, frac_ia
1996  integer :: status
1997 
1998  real, dimension(:,:,:,:), allocatable :: nnvalxw
1999  real, dimension(:), allocatable :: yyw
2000 
2001  dflag2 = .false.
2002 
2003  if (allocated(nnvalxw)) deallocate(nnvalxw, yyw)
2004  allocate(nnvalxw(4,4,2,dust_lut488%nssa), yyw(dust_lut488%nssa), stat=status)
2005  if (status /= 0) then
2006  print *, "ERROR: Failed to allocate arrays: ", status
2007  dflag2 = .false.
2008  return
2009  end if
2010 
2011  index_ia = search(tau_x, dust_lut488%aot, status, frac=frac_ia)
2012  if (status /= 0) then
2013  dflag2 = .false.
2014  return
2015  end if
2016 
2017  status = create_reduced_lut_ssa(dust_lut488, refl, x1, x2, x3, index_ia, &
2018  & r470, frac_ia, yyw)
2019  if (status /= 0) then
2020  dflag2 = .true.
2021  return
2022  end if
2023 
2024  if (refl.le.yyw(1)) then
2025  w0_x = -999.
2026  return
2027  endif
2028 
2029  if (refl.ge.yyw(4)) then
2030  w0_x = 1.0
2031  return
2032  endif
2033 
2034  index_ii = search(refl, yyw, status, frac=frac)
2035  if (status /= 0) then
2036  dflag2 = .false.
2037  return
2038  end if
2039  w0_x = frac*dust_lut488%ssa(index_ii+1) + (1.-frac)*dust_lut488%ssa(index_ii)
2040 
2041  deallocate(nnvalxw, yyw, stat=status)
2042  if (status /= 0) then
2043  print *, "WARNING: Failed to deallocate arrays: ", status
2044  return
2045  end if
2046 
2047  return
2048 
2049 end subroutine aero_470_abs_dust
2050 
2051 !--------------------------------------------------------
2052 integer function create_reduced_lut_aot(lut, refl, sza, vza, raa, imod, &
2053 & rXXX, model_frac, yy, debug) result(status)
2054 
2055  implicit none
2056 
2057  type(viirs_aerosol_lut) :: lut
2058  real, intent(in) :: refl
2059  real, intent(in) :: sza
2060  real, intent(in) :: vza
2061  real, intent(in) :: raa
2062  integer, intent(in) :: imod
2063  real, intent(in) :: rxxx
2064  real, intent(in) :: model_frac
2065  real, intent(inout),dimension(:) :: yy
2066  logical, intent(in),optional :: debug
2067 
2068  real, dimension(:,:,:,:), allocatable :: nnvalx, nnvalx1, nnvalx2
2069 
2070  integer :: index_ii, ii, jj, mbeg, nbeg
2071  integer :: ia, j, i
2072  real :: frac, xfrac, y, dy
2073 
2074  real, parameter :: pi = 3.14159
2075 
2076 ! if (debug) print *, 'aero_XXX, in: ', refl, sza, vza, raa, rXXX, imod, model_frac
2077 
2078  status = -1
2079 
2080  if (rxxx < 0.0) return
2081 
2082  index_ii = search(rxxx, lut%sfc, status, frac=frac)
2083  if (status /= 0) then
2084  return
2085  endif
2086 ! if (debug) print *, 'aero_XXX, sfc indx: ', rXXX, index_ii
2087 
2088  ii = search(raa, lut%raa, status, frac=xfrac)
2089  if (status /= 0) then
2090  return
2091  end if
2092 ! if (debug) print *, 'aero_XXX, raa: ', raa, ii, lut%raa(ii), lut%raa(ii+1)
2093 
2094  mbeg = search(sza, lut%sza, status)
2095  if (status /= 0) then
2096  print *, 'ERROR: Specified SZA not within table: ', sza, lut%sza(1), lut%sza(lut%nsza)
2097  return
2098  end if
2099  mbeg = max(0, mbeg-2)
2100  if (mbeg > lut%nsza-4) mbeg = lut%nsza - 4
2101 
2102  nbeg = search(vza, lut%vza, status)
2103  if (status /= 0) then
2104  print *, 'ERROR: Specified VZA not within table: ', vza, lut%vza(1), lut%vza(lut%nvza)
2105  return
2106  end if
2107  nbeg = max(0, nbeg-2)
2108  if (nbeg > lut%nvza-4) nbeg = lut%nvza - 4
2109 ! if (debug) print *, 'mbeg, nbeg: ', mbeg, nbeg
2110 
2111 ! -- interpolate between models if requested.
2112 ! if (present(model_frac)) then
2113  if (allocated(nnvalx)) deallocate(nnvalx, nnvalx1, nnvalx2)
2114  allocate(nnvalx(4,4,2,lut%naot), nnvalx1(4,4,2,lut%naot), &
2115  & nnvalx2(4,4,2,lut%naot), stat=status)
2116  if (status /= 0) then
2117  print *, "ERROR: Failed to allocate reduced LUT arrays: ", status
2118  return
2119  end if
2120 
2121  if (imod < lut%nssa) then ! temp check to avoid out-of-bounds issue with imod=nssa
2122  nnvalx1(:,:,:,:) = -999.0
2123  nnvalx2(:,:,:,:) = -999.0
2124  do ia = 1, lut%naot
2125  do i = 1, 4
2126  do j = 1, 4
2127  nnvalx1(i,j,1,ia) = lut%nvalx(mbeg+i,nbeg+j,ii,ia,imod,index_ii)* &
2128  & (1.0-frac) + lut%nvalx(mbeg+i,nbeg+j,ii,ia,imod,index_ii+1)*frac
2129  nnvalx1(i,j,2,ia) = lut%nvalx(mbeg+i,nbeg+j,ii+1,ia,imod,index_ii)* &
2130  & (1.0-frac) + lut%nvalx(mbeg+i,nbeg+j,ii+1,ia,imod,index_ii+1)*frac
2131 
2132  nnvalx2(i,j,1,ia) = lut%nvalx(mbeg+i,nbeg+j,ii,ia,imod+1,index_ii)* &
2133  & (1.0-frac) + lut%nvalx(mbeg+i,nbeg+j,ii,ia,imod+1,index_ii+1)*frac
2134  nnvalx2(i,j,2,ia) = lut%nvalx(mbeg+i,nbeg+j,ii+1,ia,imod+1,index_ii)* &
2135  & (1.0-frac) + lut%nvalx(mbeg+i,nbeg+j,ii+1,ia,imod+1,index_ii+1)*frac
2136 
2137  nnvalx(i,j,1,ia) = (1.0-model_frac) * nnvalx1(i,j,1,ia) + &
2138  & model_frac * nnvalx2(i,j,1,ia)
2139  nnvalx(i,j,2,ia) = (1.0-model_frac) * nnvalx1(i,j,2,ia) + &
2140  & model_frac * nnvalx2(i,j,2,ia)
2141  end do
2142  end do
2143  end do
2144 
2145 ! -- just use imod, no interpolation
2146  else
2147  nnvalx(:,:,:,:) = -999.0
2148  do ia = 1, lut%naot
2149  do i = 1, 4
2150  do j = 1, 4
2151  nnvalx(i,j,1,ia) = lut%nvalx(mbeg+i,nbeg+j,ii,ia,imod,index_ii)* &
2152  & (1.-frac) + lut%nvalx(mbeg+i,nbeg+j,ii,ia,imod,index_ii+1)*frac
2153  nnvalx(i,j,2,ia) = lut%nvalx(mbeg+i,nbeg+j,ii+1,ia,imod,index_ii)* &
2154  & (1.-frac) + lut%nvalx(mbeg+i,nbeg+j,ii+1,ia,imod,index_ii+1)*frac
2155  end do
2156  end do
2157  end do
2158  end if
2159 
2160 !--- interpolating AOT tables
2161  yy(:) = -999.0
2162  do ia = 1, lut%naot
2163  call new_intep(lut%sza, lut%vza, lut%raa, nnvalx, lut%nsza, lut%nvza, lut%nraa, ia, &
2164  & sza,vza,raa,y,dy,mbeg,nbeg,xfrac)
2165 
2166  yy(ia) = y/pi
2167  !print *,'tau, i/f=', tau(ia), y/pi, dy
2168  end do
2169 
2170  deallocate(nnvalx, nnvalx1, nnvalx2, stat=status)
2171  if (status /= 0) then
2172  print *, "WARNING: Failed to deallocate reduced LUT arrays: ", status
2173  end if
2174 
2175  return
2176 
2177 end function create_reduced_lut_aot
2178 
2179 !--------------------------------------------------------
2180 integer function create_reduced_lut_ssa(lut, refl, sza, vza, raa, index_ia, &
2181 & rXXX, aot_frac, yy, debug) result(status)
2182  implicit none
2183 
2184  type(viirs_aerosol_lut) :: lut
2185  real, intent(in) :: refl
2186  real, intent(in) :: sza
2187  real, intent(in) :: vza
2188  real, intent(in) :: raa
2189  integer, intent(in) :: index_ia
2190  real, intent(in) :: rxxx
2191  real, intent(in) :: aot_frac
2192  real, intent(inout),dimension(:) :: yy
2193  logical, intent(in),optional:: debug
2194 
2195  real, dimension(:,:,:,:), allocatable :: nnvalx
2196 
2197  real, parameter :: pi = 3.14159
2198 
2199  integer :: index_ii, ii, mbeg, nbeg
2200  integer :: iw, j, i
2201  real :: frac, xfrac, y, dy
2202  real :: dd1, dd2
2203 
2204 ! -- calculate the interpolation indices for surf. reflectance, RAA, SZA, VZA.
2205  index_ii = search(rxxx, lut%sfc, status, frac=frac)
2206  if (status /= 0) then
2207  print *, "ERROR: Failed to get interpolation index for surface: ", rxxx, status
2208  return
2209  end if
2210 
2211  ii = search(raa, lut%raa, status, frac=xfrac)
2212  if (status /= 0) then
2213  print *, "ERROR: Failed to get interpolation index for RAA: ", raa, status
2214  return
2215  end if
2216 
2217  mbeg = search(sza, lut%sza, status)
2218  if (status /= 0) then
2219  print *, "ERROR: Failed to get interpolation index for SZA: ", sza, status
2220  return
2221  end if
2222  mbeg = max(0, mbeg-2)
2223  if (mbeg > lut%nsza-4) mbeg = lut%nsza - 4
2224 
2225  nbeg = search(vza, lut%vza, status)
2226  if (status /= 0) then
2227  print *, "ERROR: Failed to get interpolation index for VZA: ", vza, status
2228  return
2229  end if
2230  nbeg = max(0, nbeg-2)
2231  if (nbeg > lut%nvza-4) nbeg = lut%nvza - 4
2232 
2233 ! -- perform interpolation of LUT as function of SSA.
2234  if (allocated(nnvalx)) deallocate(nnvalx)
2235  allocate(nnvalx(4,4,2,lut%nssa), stat=status)
2236  if (status /= 0) then
2237  print *, "ERROR: Failed to allocate nnvalx or yy array: ", status
2238  return
2239  end if
2240 
2241  do iw = 1, lut%nssa
2242  do i = 1, 4
2243  do j = 1, 4
2244  dd1 = lut%nvalx(mbeg+i,nbeg+j,ii,index_ia,iw,index_ii)* &
2245  & (1.-frac) + &
2246  & lut%nvalx(mbeg+i,nbeg+j,ii,index_ia,iw,index_ii+1)*frac
2247  dd2= lut%nvalx(mbeg+i,nbeg+j,ii,index_ia+1,iw,index_ii)* &
2248  & (1.-frac)+ &
2249  & lut%nvalx(mbeg+i,nbeg+j,ii,index_ia+1,iw,index_ii+1) &
2250  & *frac
2251 
2252  nnvalx(i,j,1,iw) = dd1* (1.-aot_frac) + dd2*aot_frac
2253 
2254  dd1 = lut%nvalx(mbeg+i,nbeg+j,ii+1,index_ia,iw,index_ii)*&
2255  & (1.-frac) + &
2256  & lut%nvalx(mbeg+i,nbeg+j,ii+1,index_ia,iw,index_ii+1)*frac
2257  dd2= lut%nvalx(mbeg+i,nbeg+j,ii+1,index_ia+1,iw,index_ii)*&
2258  & (1.-frac)+ &
2259  & lut%nvalx(mbeg+i,nbeg+j,ii+1,index_ia+1,iw,index_ii+1) &
2260  & *frac
2261 
2262  nnvalx(i,j,2,iw) = dd1* (1.-aot_frac) + dd2*aot_frac
2263 
2264  end do
2265  end do
2266  end do
2267 
2268 !-- interpolating W0 tables
2269  do iw = 1, lut%nssa
2270  call new_intep(lut%sza, lut%vza, lut%raa, nnvalx, &
2271  & lut%nsza, lut%nvza, lut%nraa, iw, sza,vza,raa,y, &
2272  & dy,mbeg,nbeg,xfrac)
2273 
2274  yy(iw) = y/pi
2275  end do
2276 
2277 end function create_reduced_lut_ssa
2278 
2279 !--------------------------------------------------------
2280  integer function search(xbar,x,status,frac) result(i)
2281 !
2282 ! purpose
2283 ! locate position in table of point at which interpolation is
2284 ! required
2285 !
2286 ! usage
2287 ! i = search (xbar, x, status, frac=frac)
2288 !
2289 ! description of parameters
2290 ! xbar - point at which interpolation is required
2291 ! x - array containing independent variable
2292 ! status - integer indicates success (0) or failure (/=0)
2293 ! frac - optional parameter for fraction.
2294 !
2295  implicit none
2296 
2297  real, intent(in) :: xbar
2298  real, dimension(:), intent(in) :: x
2299  integer, intent(inout) :: status
2300  real, intent(inout), optional :: frac
2301 
2302  integer :: n
2303 
2304  integer :: m, k
2305  integer :: icnt
2306  real, parameter :: b = 0.69314718 ! (ln(2))
2307 
2308  status = -1
2309 
2310  n = size(x)
2311  icnt = 0
2312  if (n < 2) then
2313  print *, "Search n is less than 2."
2314  return
2315  end if
2316 
2317  if (x(1) > x(2)) then
2318  print *, "Search table is not in increasing order."
2319  status = -1
2320  return
2321  end if
2322 
2323  m = int((log(float(n)))/b)
2324  i = 2**m
2325  if (i >= n) i = n-1
2326  k = i
2327  do
2328  k=k/2
2329 
2330  if (k == 0) icnt = icnt + 1
2331  if (icnt >= 2) exit ! starting an infinite loop, bail out!
2332 
2333 ! -- check if we found our indices and exit the loop if so.
2334 ! -- otherwise, reset i and loop again.
2335  if (xbar >= x(i) .AND. xbar < x(i+1)) then
2336  status = 0
2337  exit
2338  else
2339  if (xbar > x(i)) then
2340  i = i + k
2341  if (i >= n) then
2342  i = n - 1
2343  end if
2344  else
2345  i = i - k
2346  end if
2347  end if
2348 
2349  end do
2350 
2351 ! -- detect failed search above via icnt. if >= 2, do sequential search,
2352 ! -- else above search was successful.
2353  if (icnt >= 2 .OR. status /= 0) then
2354  do i = 1, n-1
2355  if (xbar >= x(i) .AND. xbar <= x(i+1)) then
2356  status = 0
2357  exit
2358  end if
2359  end do
2360  end if
2361 
2362 ! -- calculate fraction if desired and search was successful.
2363  if (status == 0 .AND. present(frac)) then
2364  if (i < n) then
2365  frac = (xbar-x(i))/ (x(i+1)-x(i))
2366  else
2367  frac = 1.0
2368  end if
2369  end if
2370 
2371  return
2372 
2373  end function search
2374 
2375 
2376 ! -- https://en.wikipedia.org/wiki/Extrapolation
2377  real function extrap(x, xx, yy, status) result(res)
2378  implicit none
2379 
2380  real, intent(in) :: x
2381  real, dimension(:), intent(in) :: xx
2382  real, dimension(:), intent(in) :: yy
2383  integer, intent(out) :: status
2384 
2385  real, dimension(2) :: r
2386  integer :: n, i
2387 
2388  res = -999.0
2389  status = -1
2390 
2391  n = size(xx, 1)
2392  status = linfit(xx(n-1:n), yy(n-1:n), r)
2393  if (status /= 0) then
2394  print *, "ERROR: linfit failed, skipping: ", status
2395  return
2396  end if
2397 
2398  res = r(1) + (x)*r(2)
2399 
2400 ! -- if extrapolation produces a negative AOT, exclude the last node point
2401 ! -- and try again.
2402  if (res < 3.5) then
2403  i = n-1
2404  do while (res < 3.5 .AND. i >= 2)
2405  status = linfit(xx(i-1:i), yy(i-1:i), r)
2406  if (status /= 0) then
2407  print *, "ERROR: linfit failed, skipping: ", status
2408  return
2409  end if
2410 
2411  res = r(1) + (x)*r(2)
2412  i = i-1
2413  end do
2414  if (i < 2) then
2415  status = 1
2416  return
2417  end if
2418  end if
2419 
2420  status = 0
2421  return
2422 
2423  end function extrap
2424 
2425  integer function linfit(x, y, r) result(status)
2426  implicit none
2427 
2428  real, dimension(:), intent(in) :: x
2429  real, dimension(:), intent(in) :: y
2430  real, dimension(2), intent(inout) :: r
2431 
2432  real :: sx, sy
2433  real :: sxx, syy, sxy
2434 
2435  integer :: i, n
2436 
2437  n = size(x)
2438  sx = 0.0
2439  sy = 0.0
2440  sxy = 0.0
2441  sxx = 0.0
2442  do i = 1, n
2443  sx = sx + x(i)
2444  sy = sy + y(i)
2445  sxx = sxx + (x(i) * x(i))
2446  sxy = sxy + (x(i) * y(i))
2447  syy = syy + (y(i) * y(i))
2448  end do
2449 
2450  r(2) = ((n*sxy) - (sx*sy))/((n*sxx)-(sx*sx))
2451  r(1) = (sy/n)-(r(2)*sx/n)
2452 
2453  status = 0
2454  return
2455 
2456  end function linfit
2457 
2458 !-------------------------------------------------------------
2459  subroutine new_intep(x1a,x2a,x3a,ya,m,n,l,ia,x1,x2,x3,y,dy, &
2460  & mbeg,nbeg,frac)
2461  implicit none
2462 
2463  real, dimension(:), intent(in) :: x1a
2464  real, dimension(:), intent(in) :: x2a
2465  real, dimension(:), intent(in) :: x3a
2466  real, dimension(:,:,:,:), intent(in) :: ya
2467  integer, intent(in) :: m, n, l
2468  integer, intent(in) :: ia
2469  real, intent(in) :: x1, x2, x3
2470  real, intent(inout) :: y, dy
2471  integer, intent(in) :: mbeg, nbeg
2472  real, intent(in) :: frac
2473 
2474  real, dimension(4) :: xx2a, xx1a
2475  real, dimension(4) :: yntmp, ymtmp
2476  real, dimension(2) :: yltmp
2477 
2478  integer :: j, k
2479 
2480  !dimension x1a(m),x2a(n),x3a(l),ya(4,4,2,10)
2481  !dimension xx2a(4), xx1a(4)
2482  !dimension yntmp(4),ymtmp(4),yltmp(2)
2483 
2484  do j=1,4
2485  do k=1,4
2486  yltmp(1)=ya(j,k,1,ia)
2487  yltmp(2)=ya(j,k,2,ia)
2488  yntmp(k) = yltmp(1)*(1.-frac) + yltmp(2)*frac
2489  xx2a(k) = x2a(k+nbeg)
2490  end do
2491  call polint(xx2a,yntmp,4,x2,ymtmp(j),dy)
2492  xx1a(j) = x1a(j+mbeg)
2493 
2494  end do
2495 
2496  call polint(xx1a,ymtmp,4,x1,y,dy)
2497 
2498  return
2499  end subroutine new_intep
2500 
2501 
2502  subroutine polint(xa,ya,n,x,y,dy)
2503  implicit none
2504  real, dimension(:), intent(in) :: xa
2505  real, dimension(:), intent(in) :: ya
2506  integer, intent(in) :: n
2507  real, intent(in) :: x
2508  real, intent(inout) :: y
2509  real, intent(inout) :: dy
2510 
2511  integer, parameter :: nmax = 50
2512  real, dimension(nmax) :: c
2513  real, dimension(nmax) :: d
2514 
2515  integer :: i, m, ns
2516  real :: ho, hp, dif, dift
2517  real :: w, den
2518 
2519  ns=1
2520  dif=abs(x-xa(1))
2521  do 11 i=1,n
2522  dift=abs(x-xa(i))
2523  if (dift.lt.dif) then
2524  ns=i
2525  dif=dift
2526  endif
2527  c(i)=ya(i)
2528  d(i)=ya(i)
2529 11 continue
2530  y=ya(ns)
2531  ns=ns-1
2532  do 13 m=1,n-1
2533  do 12 i=1,n-m
2534  ho=xa(i)-x
2535  hp=xa(i+m)-x
2536  w=c(i+1)-d(i)
2537  den=ho-hp
2538 ! if(den.eq.0.) print *,'den=',den
2539 ! if(den.eq.0.)pause
2540  den=w/den
2541  d(i)=hp*den
2542  c(i)=ho*den
2543 12 continue
2544  if (2*ns.lt.n-m)then
2545  dy=c(ns+1)
2546  else
2547  dy=d(ns)
2548  ns=ns-1
2549  endif
2550  y=y+dy
2551 13 continue
2552 
2553  return
2554 
2555  end subroutine polint
2556 
2557 
2558 end module viirs_aerosol_luts
subroutine, public aero_412_abs(dflag, refl, x1, x2, x3, mm, nn, ll, r412, tau_x, w0_x)
subroutine, public aero_470_abs_dust(dflag2, refl, x1, x2, x3, mm, nn, ll, r470, tau_x, w0_x)
subroutine, public aero_650_dust(dflag, refl, x1, x2, x3, mm, nn, ll, ma, r650, tau_x650, tau_x650_flag, tau_x470_flag, tau_x412, tau_x470, tau_x412_flag_91, trflg)
subroutine, public aero_412_abs_dust(dflag, refl, x1, x2, x3, mm, nn, ll, r412, tau_x, w0_x)
subroutine, public aero_412(dflag, refl, x1, x2, x3, mm, nn, ll, ma, imod, r412, tau_x412, tau_x412_flag, trflg, model_frac, debug)
subroutine, public aero_412_dust(dflag, refl, x1, x2, x3, mm, nn, ll, ma, imod, r412, tau_x412, tau_x412_flag, trflg, model_frac, debug)
subroutine search(dflag, xbar, x, n, i)
subroutine, public aero_470_dust(dflag, refl, x1, x2, x3, mm, nn, ll, ma, imod, r470, tau_x470, tau_x470_flag, trflg, model_frac, debug)
#define max(A, B)
Definition: main_biosmap.c:61
subroutine, public aero_470(dflag, refl, x1, x2, x3, mm, nn, ll, ma, imod, r470, tau_x470, tau_x470_flag, trflg, model_frac, debug)
subroutine, public unload_viirs_aerosol_luts(status)
string & trim(string &s, const string &delimiters)
Definition: EnvsatUtil.cpp:29
#define pi
Definition: vincenty.c:23
subroutine, public polint(xa, ya, n, x, y, dy)
subroutine, public new_intep(x1a, x2a, x3a, ya, m, n, l, ia, x1, x2, x3, y, dy, mbeg, nbeg, frac)
subroutine, public aero_650(dflag, refl, x1, x2, x3, mm, nn, ll, ma, r650, tau_x650, tau_x650_flag, tau_x470_flag, tau_x412, tau_x470, tau_x412_flag_91, trflg)
integer function, public load_viirs_aerosol_luts(lut_file)
#define abs(a)
Definition: misc.h:90
subroutine, public aero_470_abs(dflag2, refl, x1, x2, x3, mm, nn, ll, r470, tau_x, w0_x)