OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
find_v_vegset_nc4.f
Go to the documentation of this file.
1  subroutine find_v_veg(month,season,realbuf,tmpvg,
2  1 r412sv,r470sv,gzflg,outbufvg,tau_x470_flag) !modified by CH 11/18/2016
3 c. subroutine myd_aero_veg(imod,season,itmp,jtmp,realbuf,tmpvg,
4 c. 1 outbufvg)
5 
6 c... Main subroutines for aerosol retrieval over vegetated surfaces.
7 c... Written by Myeong-Jae Jeong (MJ)
8 c... Last modified Aug 08, 2011
9 c...
10 c... Inputs:
11 c... season: 1:DJF; 2:MAM; 3:JJA; 4:SON
12 c... imod: aerosol model index for 470nm --> no longer an input
13 c... nvalx470: TOA refl LUT for 470nm
14 c... nvalx650: TOA refl LUT for 650nm
15 c... xlcvr_2: MODIS land cover (IGBP)
16 c... regid_2: Region Index (to choose a set of aerosol models
17 c... depending on regions and land cover)
18 c... realbuf: see below for definitions
19 c... tmpvg:
20 c...
21 c... Output:
22 c... outbufvg: see the end of this subroutine for definitions
23 c...
24 
25  use viirs_aerosol_luts, only: aero_470, aero_650
26 
27 c include 'aottbl.inc'
28  include 'newaottbl.inc'
29 
30  parameter(nx2=3600, ny2=1800) ! landcover data dimension
31 
32  logical dflag, debug, do_sv, do_nir
33 
34  integer bflag ! jlee added
35  integer month, season, gzflg
36 
37  real realbuf(26), tmpvg(7), outbufvg(21), xnvalm6(6)
38  real r412sv, r470sv
39 
40 c real nvalx470(10,46,30,10,4,24),nvalx650(10,46,30,10,24)
41  real nval(10,46,30), yy(10), yyw(8) !tau(10),
42 c. real xnvalm6(6), realbuf(13), outbuf(20)
43 c real*4 xlcvr_2(nx2,ny2)
44  real tau_x470_1, tau_x470_2, tau_x470_3
45  real tau_x470sv_1, tau_x470sv_2, tau_x470sv_3
46  integer tau_x470_flag, tau_x650_flag, tau_ini_flag, tau_x470_flag_ini
47  integer tau_x470_flag1, tau_x470_flag2, tau_x470_flag3
48  integer tau_x470sv_flag1, tau_x470sv_flag2, tau_x470sv_flag3
49 
50 c real theta0(10), theta(46), phi(30)
51 c real sfc_ref412(20), sfc_ref470(24), sfc_ref650(24)
52  real*4 r412db,r470db,r650db,cl_flag
53  real xtau(3),ssa(3),qa_flag(4),aot_mod(6) !,w0_470(4)
54  character*4 w0_name470(4)
55  character*12 aer_tab(10)
56  real*4 ctharr(3) ! 08/05/2011
57  integer imodarr(3) ! 08/05/2011
58  real modfrac(3) ! 16 January 2018 JLee
59 
60 c common /angle_node/ theta0, theta, phi
61 c common /sfcref_node/ sfc_ref412, sfc_ref470, sfc_ref650
62  common /fname_node/ aer_tab, w0_name470
63  data pi /3.14159/
64 C data theta0 /0.0,8.0,16.0,24.0,32.0,40.0,48.0,56.0,64.0,72.0/
65 C data tau /0.0, 0.1, 0.3, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5/
66 C data w0_470 /0.91, 0.94, 0.96, 0.99/
67 C data sfc_ref412 /1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,
68 C 1 11.,12.,13.,14.,15.,16.,17.,18.,19.,20./
69 C data sfc_ref470 /1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,
70 C 1 11.,12.,13.,14.,15.,16.,17.,18.,19.,20.,
71 C 2 21.,22.,23.,24./
72 C data sfc_ref650 /1.,3.,5.,7.,9.,11.,13.,15.,17.,19.,21.,
73 C 1 23.,25.,27.,29.,31.,33.,35.,37.,39.,41.,
74 C 2 43.,45.,47./
75 
76 c. data aer_tab /'table_0.0aot', 'table_0.1aot', 'table_0.3aot',
77 c. 1 'table_0.5aot', 'table_1.0aot', 'table_1.5aot',
78 c. 2 'table_2.0aot', 'table_2.5aot', 'table_3.0aot',
79 c. 3 'table_3.5aot'/
80  data w0_name470 /'0.91','0.94','0.96','0.99'/
81 
82  bflag = 0 ! jlee added
83  debug = .false.
84 
85 c-------------------------------------------
86 c... define nodes for angles
87 c do i = 1, 46
88 c theta(i) = 2.*float(i-1)
89 c enddo
90 
91 c do i = 1, 30
92 c phi(i) = 6. + 6.*float(i-1)
93 c enddo
94 
95  mm = 12 ! solar zenith
96  nn = 46 ! satellite zenith
97  ll = 30 ! rel. azimuth
98  ma = 10 ! tau
99  mw = 8 ! ssa
100 c-------------------------------------------
101 
102 c-------------------------------------------
103 c... IGBP Landcover data
104  xlatbeg=-90.0
105  xlonbeg=-180.0
106  xyintv2=0.10
107 c-------------------------------------------
108 
109 c... ************************************************
110 c... ************************************************
111 c... come back later to set these threshold values, if necessary...(MJ) Jun 15, 2010
112  xmnsfc1=1.0 ! allowed min. sfc refl. for B1
113  xmxsfc1=47.0 ! allowed max. sfc refl. for B1
114  xmnsfc3=1.0 ! allowed min. sfc refl. for B3
115  xmxsfc3=24.0 ! allowed max. sfc refl. for B3
116 c... ************************************************
117 c... ************************************************
118 
119 c.? xlat5 = realbuf(1)
120 c.? xlong5 = realbuf(2)
121  xlat = realbuf(1)
122  xlong = realbuf(2)
123  sza = realbuf(3)
124  xthet = realbuf(4)
125  xphi = realbuf(5)
126  xnvalm6(1)=tmpvg(1) ! (B1)
127  xnvalm6(2)=tmpvg(2) ! (B3)
128  xnvalm6(3)=tmpvg(3) ! (B8)
129  xnvalm6(4)=tmpvg(4) ! 2.1um (B7)
130  xnvalm6(5)=tmpvg(5) ! 850nm (B2)
131  xnvalm6(6)=tmpvg(6) ! 1.2um (B5)
132 c band26 =realbuf(10) ! 1.38um(B26)
133 c band02 =xnvalm6(5)
134 c band05 =xnvalm6(6)
135 c band07 =tmpvg(7) ! jlee commented out
136 c band08 =xnvalm6(3) ! (B8)
137 c band03 =xnvalm6(2) ! (B3)
138 c band01 =xnvalm6(1) ! (B1)
139 c cl_flag =realbuf(13)
140  r412db =realbuf(24)*100.
141  r470db =realbuf(25)*100.
142  r650db =realbuf(26)*100.
143 
144 c print *, itmp, jtmp, realbuf, lsf ! test
145 
146 c... ============================================================
147 c -- sun glint mask
148 c
149  cc = 3.14159/180.
150  psi = acos(cos(sza*cc)*cos(xthet*cc) +
151  1 sin(sza*cc)*sin(xthet*cc)*cos(xphi*cc))
152  glint_ang = psi/cc
153 
154 c if (abs(psi/cc).lt.35.0) go to 10
155 
156 c -- scattering angle (scat_ang)
157 c
158  psi = acos(cos(sza*cc)*cos(xthet*cc) -
159  1 sin(sza*cc)*sin(xthet*cc)*cos(xphi*cc))
160  scat_ang = 180. - psi/cc
161 
162 c if (scat_ang .gt. 158.) go to 10
163 c if (scat_ang .gt. 175.) go to 10
164 
165 c... ######################################
166  lcvr=1 ! initialization
167  ioprg=5 ! ititialization
168 c... get lcvr (IGBP landcover type; integer) data here
169  idx=int((xlong-xlonbeg)/xyintv2) + 1
170  idy=int((xlat-xlatbeg)/xyintv2) + 1
171 
172  if(idx.ge.1.and.idx.le.nx2.and.idy.ge.1.and.idy.le.ny2) then
173  sfc_typ = xlcvr_2(idx,idy)
174  xreg_id = regid_2(idx,idy) ! 08/05/2011
175  else
176 c sfc_typ = 1.0*lcvr ! @@@@@ just for test
177  print *, 'lcvr data out of bound: idx,nx2,idy,ny2: ',idx,nx2,idy,ny2,xlat,xlong
178 c noob=noob+1
179  stop ! @@@@@
180  endif
181  lcvr=int(sfc_typ)
182  ioprg=int(xreg_id) ! 08/05/2011
183 
184 c. if(itmp.ge.100.and.itmp.lt.110.and.jtmp.ge.200.and.
185 c. 1 jtmp.lt.210) print *, sfc_typ, lcvr, xreg_id, ioprg
186 c... ######################################
187 
188  !xmu = cos(sza * 3.14159/180.)
189 
190 c... refl unit conversion (pi*L/F/xmu --> L/F)
191  !do i = 1, 6
192  ! xnvalm6(i) = xnvalm6(i) * xmu/3.14159
193  !enddo
194 
195 c if (xphi.gt.179.99) go to 10
196 c if (xphi.gt.179.99) return ! choose this or below (MJ)
197  if (xphi.gt.179.99) xphi=179.99 ! choose this or above (MJ)
198  if (xphi.lt.6.0) xphi = 6.
199 
200  x1 = sza
201  x2 = xthet
202  x3 = xphi
203 
204 c... ============================================================
205 
206  refl1= xnvalm6(3) ! 412 nm
207  refl6= xnvalm6(1) ! 650 nm
208  refl3= xnvalm6(2) ! 470 nm
209  !refl21=band07 ! 2.1um ! jlee commneted out
210  refl21 = xnvalm6(4) !2.1um ! jlee added
211  refl865 = xnvalm6(5) !865 nm !jlee added
212  refln21 = xnvalm6(4)*3.14159/cos(sza*cc) ! 2.1 um, normalized ! jlee added
213  refln865 = xnvalm6(5)*3.14159/cos(sza*cc)! 865 nm, normalized ! jlee added
214 
215  refn672_rc = -999.0
216  refn865_rc = -999.0
217 
218  call calc_rc_ref672(sza,xthet,xphi,refl6,refn672_rc)
219  call calc_rc_ref865(sza,xthet,xphi,refl865,refn865_rc)
220 
221  sirndvi=(xnvalm6(6)-xnvalm6(4))/(xnvalm6(6)+xnvalm6(4)) !jlee added, ndvi_swir
222  rc_ndvi =(refn865_rc-refn672_rc)/(refn865_rc+refn672_rc) !jlee added, Rayleigh-corrected ndvi_vis
223  toa_ndvi = (refl865-refl6)/(refl865+refl6)
224  airmass = 1.0/cos(sza*cc)*1.0/cos(vza*cc) !jlee added
225  !print *, toa_ndvi, rc_ndvi
226 c... ============================================================
227 
228 !===================================================================
229 c---------------------------------------------------------------
230 c Initialization
231 c---------------------------------------------------------------
232 c --intermediate parameters
233  w0_x = -999.
234  w0_int = -999.
235  tau_x412 = -999.
236  tau_x412_91 = -999.
237  tau_x470 = -999.
238  tau_x650 = -999.
239  xxrat = -999.
240  xxrat2 = -999.
241  aot = -999.
242 c -- output parameters
243  tau550 = -999.
244  alpha = -999.
245 
246  do i = 1, 3
247  xtau(i) = -999.
248  ssa(i) = -999.
249  enddo
250 
251  do i = 1, 4
252  qa_flag(i) = 0
253  enddo
254 
255  do i = 1, 6
256  aot_mod(i) = -999.
257  enddo
258 
259  do i = 1,21
260  outbufvg(i) = -999.
261  enddo
262 c-----------------------------------------------------------------------
263 c... estimate surface reflectance at 470nm & 650nm
264 c-----------------------------------------------------------------------
265  xeas_b1 = -999.0 ! initialization
266  xeas_b3 = -999.0
267 !==============================================================
268 !jlee added
269  ! do_sv = false and do_nir = false: use swir method
270  ! do_sv = false and do_nir = true: use nir method
271  ! do_sv = true and do_nir = false: use 2.2 um database, alpha from swir method
272  ! do_sv = true and do_nir = true: use 2.2 um database, alpha from nir method
273  do_sv = .false.
274  do_nir = .false.
275  if (gzflg.eq.15.or.gzflg.eq.19) then !swir method only over india
276  do_sv = .false.
277  do_nir = .false.
278  elseif (lcvr.eq.12.and.rc_ndvi.lt.0.4.and.(ioprg.eq.1.or.ioprg.eq.3)) then !use nir method over bright croplands in North and South America
279  do_sv = .false.
280  do_nir = .true.
281  else !default, swir method only
282  do_sv = .false.
283  do_nir = .false.
284  endif
285 
286  if (do_nir) then !this condition is also used other parts of the code below
287  call get_sfcrfl_nir_vis(season,lcvr,refn865_rc,rc_ndvi,
288  1 xeas_b1,xeas_b3)
289  else
290  call get_sfcrfl_swir_vis(season,lcvr,refln21,rc_ndvi,
291  1 xeas_b1,xeas_b3)
292  endif
293 
294  if (xeas_b1.lt.xmnsfc1.and.xeas_b1.gt.-900.0) xeas_b1=xmnsfc1
295  if (xeas_b3.lt.xmnsfc3.and.xeas_b3.gt.-900.0) xeas_b3=xmnsfc3
296 
297  refl = refl6
298  r650 = xeas_b1
299  dflag = .false.
300  trflg = 1
301  call aero_650(dflag,refl,x1,x2,x3,mm,nn,ll,ma,
302  1 r650,tau_ini,tau_ini_flag,0,0.0,0.0,0,trflg)!tau_x470_flag,tau_x412,tau_x470,tau_x412_flag_91
303 
304 ! call aero_650veg(refl,x1,x2,x3,mm,nn,ll,ma,
305 ! 1 r650,tau_ini,tau_ini_flag,
306 ! 2 tau_x470_flag_ini,tau_x470_ini)
307 
308 ! imod=imodarr(1)
309 ! refl = refl3
310 ! r470 = xeAs_B3
311 !
312 ! call aero_470veg(refl,x1,x2,x3,mm,nn,ll,ma,imod,r470,
313 ! 1 tau_ini,tau_ini_flag)
314 
315  if (tau_ini.gt.0.3) then
316  xeas_b1 = -999.0
317  xeas_b3 = -999.0
318 
319  rc_ndvi_c = rc_ndvi+0.257*tau_ini*rc_ndvi
320  1 +0.172*(airmass-1.0)*tau_ini*rc_ndvi
321 
322  if (scat_ang.lt.120) then
323  rc_ndvi_c = rc_ndvi+0.257*tau_ini*rc_ndvi
324  1 +0.172*(airmass-1.0)*tau_ini*rc_ndvi
325  2 +0.025*(30.0-(scat_ang-90.0))*tau_ini*rc_ndvi
326  endif
327 
328  if (do_nir) then
329  call get_sfcrfl_nir_vis(season,lcvr,refn865_rc,rc_ndvi_c,
330  1 xeas_b1,xeas_b3)
331 
332 !! swir sfc for high aod over cropland
333 ! refl = refl6
334 ! r650 = xeAs_B1
335 ! call aero_650veg(refl,x1,x2,x3,mm,nn,ll,ma,
336 ! 1 r650,tau_ini,tau_ini_flag,
337 ! 2 tau_x470_flag_ini,tau_x470_ini)
338 !
339 ! if (tau_ini.gt.0.4) then
340 ! call get_sfcrfl_swir_vis(season,lcvr,refln21,rc_ndvi_c,
341 ! 1 xeAs_B1,xeAs_B3)
342 ! endif
343  else
344  call get_sfcrfl_swir_vis(season,lcvr,refln21,rc_ndvi_c,
345  1 xeas_b1,xeas_b3)
346  endif
347 
348  if (xeas_b1.lt.xmnsfc1.and.xeas_b1.gt.-900.0) xeas_b1=xmnsfc1
349  if (xeas_b3.lt.xmnsfc3.and.xeas_b3.gt.-900.0) xeas_b3=xmnsfc3
350  endif! if (tau_ini.gt.0.3) then
351 
352  if (xeas_b3.gt.24.0) return
353 
354 !end jlee added
355 !===================================================================
356 
357 c sfc_typ = -999.
358 
359 c print *, 'xeAs_B1,xmnsfc1,xeAs_B3,xmnsfc3',xeAs_B1,xmnsfc1,xeAs_B3,xmnsfc3
360 
361 c if(xeAs_B1.lt.xmnsfc1.or.xeAs_B1.gt.xmxsfc1.or.
362 c 1 xeAs_B3.lt.xmnsfc3.or.xeAs_B3.gt.xmxsfc3) return ! return to main w/o retrieval
363 c print *, 'xeAs_B1,xeAs_B3',xeAs_B1,xeAs_B3,refl21,lcvr
364 
365 c-----------------------------------------------------------------------
366 
367 c if(itmp.ge.100.and.itmp.lt.200.and.
368 c 1 jtmp.ge.200.and.jtmp.lt.300)
369 c if(xeAs_B1.gt.0.and.xeAs_B3.gt.0)
370 c 2 print *, 'sfc refl= ',xeAs_B1,xeAs_B3,r650db,r470db
371 c... write outputs on fort.23 for tests !MJ@@@@@
372 c. if(xeAs_B1.gt.0.and.xeAs_B1.lt.47.and.xeAs_B3.gt.0.and.
373 c. 1 xeAs_B3.lt.24)
374 c. 2write(23,2323) itmp,jtmp,xlat,xlong,sza,xthet,xphi,lcvr,
375 c. 3 (tmpvg(kk),kk=1,6),xeAs_B1,xeAs_B3,r650db,r470db,r412db,
376 c. 4 sirndvi,toa_ndvi,ioprg,lcvr
377 c2323 format(2(i4,2x),2(f9.4,2x),3(f8.3,2x),i3,6(2x,e14.6),
378 c 1 5(2x,f8.3),2(2x,f9.4),2(2x,i3))
379 
380 c... Begin aerosol retrieval .....
381 c... getting thresholds for selecting a set of aerosol models
382 c... for retrieval depending on regions (08/05/2011)
383 !... below 37 lines --> @@new@@
384  ctharr(:) = -999.0
385  imodarr(:) = -999 !1:0.91, 2:0.94, 3:0.96, 4:0.995
386  modfrac(:) = 0.0 ! aerosol model fraction of imodarr()+1
387  if(ioprg.eq.1) then ! N. America
388  if (season.eq.1) then
389  ctharr(1)=0.2
390  ctharr(2)=0.5
391  ctharr(3)=1.0
392  imodarr(1)=3
393  imodarr(2)=2
394  imodarr(3)=2
395  modfrac(:)=0.5
396  else
397  ctharr(1)=0.2
398  ctharr(2)=0.5
399  ctharr(3)=0.8
400  imodarr(1)=4
401  imodarr(2)=4
402  imodarr(3)=2
403  modfrac(:)=0.0
404  modfrac(3)=0.5
405  endif
406  elseif(ioprg.eq.2) then ! China
407  if (season.eq.3) then
408  ctharr(1)=0.2
409  ctharr(2)=0.5
410  ctharr(3)=1.0
411  imodarr(1)=3
412  imodarr(2)=3
413  imodarr(3)=2
414  modfrac(:)=0.0
415  else
416  ctharr(1)=0.2
417  ctharr(2)=0.5
418  ctharr(3)=1.0
419  imodarr(1)=3
420  imodarr(2)=2
421  imodarr(3)=1
422  modfrac(1)=0.0
423  modfrac(2:3)=0.5
424  endif
425  elseif(ioprg.eq.9) then ! Korea and Japan
426  if (season.eq.3) then
427  ctharr(1)=0.2
428  ctharr(2)=0.5
429  ctharr(3)=1.0
430  imodarr(1)=3
431  imodarr(2)=3
432  imodarr(3)=2
433  modfrac(:)=0.0
434  else
435  ctharr(1)=0.2
436  ctharr(2)=0.5
437  ctharr(3)=1.5
438  imodarr(1)=3
439  imodarr(2)=2
440  imodarr(3)=1
441  modfrac(1)=0.0
442  modfrac(2:3)=0.5
443  endif
444  elseif(ioprg.eq.3) then ! S. America
445  if (season.eq.1.) then
446  ctharr(1)=0.2
447  ctharr(2)=0.5
448  ctharr(3)=1.0
449  imodarr(1)=3
450  imodarr(2)=3
451  imodarr(3)=3
452  modfrac(:)=0.0
453  else
454  ctharr(1)=0.2
455  ctharr(2)=0.5
456  ctharr(3)=1.0
457  imodarr(1)=2
458  imodarr(2)=2
459  imodarr(3)=2
460  modfrac(:)=0.5
461  endif
462  elseif(ioprg.eq.4) then !India
463  if (season.eq.3) then ! Summer
464  ctharr(1)=0.2
465  ctharr(2)=0.5
466  ctharr(3)=1.0
467  imodarr(1)=3
468  imodarr(2)=2
469  imodarr(3)=2
470  modfrac(:)=0.5
471  else
472  ctharr(1)=0.2
473  ctharr(2)=0.5
474  ctharr(3)=1.0
475  imodarr(1)=3
476  imodarr(2)=2
477  imodarr(3)=1
478  modfrac(:)=0.0
479  modfrac(3)=0.5
480  end if
481  elseif(ioprg.eq.5) then ! N. Africa
482  ctharr(1)=0.2
483  ctharr(2)=0.5
484  ctharr(3)=1.0
485  imodarr(1)=2
486  imodarr(2)=1
487  imodarr(3)=1
488  modfrac(1)=0.0
489  modfrac(2:3)=0.5
490  elseif(ioprg.eq.6) then ! S. Africa
491  if (season.eq.3.or.season.eq.4) then
492  ctharr(1)=0.2
493  ctharr(2)=0.5
494  ctharr(3)=1.0
495  imodarr(1)=1
496  imodarr(2)=1
497  imodarr(3)=1
498  modfrac(:)=0.5
499  elseif (season.eq.1) then
500  ctharr(1)=0.2
501  ctharr(2)=0.5
502  ctharr(3)=1.0
503  imodarr(1)=2
504  imodarr(2)=2
505  imodarr(3)=1
506  modfrac(:)=0.5
507  modfrac(3)=0.2
508  else
509  ctharr(1)=0.2
510  ctharr(2)=0.5
511  ctharr(3)=1.0
512  imodarr(1)=2
513  imodarr(2)=2
514  imodarr(3)=1
515  modfrac(:)=0.5
516  end if
517  elseif(ioprg.eq.7) then ! SE Asia
518  if (season.eq.1) then ! Winter
519  ctharr(1)=0.2
520  ctharr(2)=0.5
521  ctharr(3)=1.0
522  imodarr(1)=1
523  imodarr(2)=1
524  imodarr(3)=1
525  modfrac(:)=0.5
526  else
527  ctharr(1)=0.2
528  ctharr(2)=0.5
529  ctharr(3)=1.5
530  imodarr(1)=2
531  imodarr(2)=1
532  imodarr(3)=1
533  modfrac(1:2)=0.5
534  modfrac(3)=0.0
535  end if
536  elseif(ioprg.eq.8.or.ioprg.eq.10) then ! Europe, Central Asia
537  if (season.eq.3) then
538  ctharr(1)=0.2
539  ctharr(2)=0.5
540  ctharr(3)=1.0
541  imodarr(1)=3
542  imodarr(2)=3
543  imodarr(3)=3
544  modfrac(:)=0.0
545  else
546  ctharr(1)=0.2
547  ctharr(2)=0.5
548  ctharr(3)=1.0
549  imodarr(1)=3
550  imodarr(2)=3
551  imodarr(3)=2
552  modfrac(:)=0.0
553  modfrac(3)=0.5
554  endif
555  elseif(ioprg.eq.13) then ! C. America
556  ctharr(1)=0.2
557  ctharr(2)=0.5
558  ctharr(3)=1.0
559  imodarr(1)=2
560  imodarr(2)=1
561  imodarr(3)=1
562  modfrac(1)=0.0
563  modfrac(2:3)=0.5
564  else !default
565  ctharr(1)=0.2
566  ctharr(2)=0.5
567  ctharr(3)=1.0
568  imodarr(1)=2
569  imodarr(2)=2
570  imodarr(3)=1
571  modfrac(:)=0.0
572  modfrac(3)=0.5
573  endif
574 
575 c---------------------------------------------------------------
576 c Screen for pixels outside reasonable ranges of reflectance
577 c---------------------------------------------------------------
578 c... will need to re-open the following lines .... (MJ)
579  if (refl1.gt.0.0.and.refl1.lt.0.09.and.
580  1 refl6.gt.0.0.and.refl6.lt.0.14) go to 11
581  if (toa_ndvi.gt.0.1.and.
582  1 refln21.gt.0.01.and.refln21.le.0.25) go to 11 ! remove water/bright sfc contamination
583  return ! return to main w/o retrieval
584 
585 11 continue
586 
587 c--------------------------------------------------------
588 c Input surface reflectance at 470 nm
589 c--------------------------------------------------------
590 c Input surface reflectance at 470 nm
591  r470 = xeas_b3
592 c
593 c Input toa reflectance (L/F) at 470 nm
594  refl = refl3
595 c x3 = xphi ! ori. position (moved up)
596 
597 c Retrieving 470 nm AOT
598 c imod = 1 ! w0 = 0.91
599 c imod = 2 ! w0 = 0.94
600 c imod = 3 ! w0 = 0.96
601 c imod = 4 ! w0 = 0.995
602 !c... below 40 lines --> AOT-dependent aerosol model selection (new@@@)
603  as_21=-999.0
604  tau_x470 = -999.
605  tau_x470_1 = -999.
606  tau_x470_2 = -999.
607  tau_x470_3 = -999.
608  tau_x470_flag=0
609  tau_x650_flag=0
610  tau_x470_flag1=0
611  tau_x470_flag2=0
612  tau_x470_flag3=0
613  cth0=ctharr(1)
614  cth1=ctharr(2)
615  cth2=ctharr(3)
616 c print *, ctharr
617 c print *, cth0, cth1, cth2
618  imod=imodarr(1)
619 
620  call aero_470(dflag,refl,x1,x2,x3,mm,nn,ll,ma,
621  1 imod,r470,tau_x470_1,tau_x470_flag1,trflg,0.0,debug)!model_frac
622 
623  go to 333 ! skip 2.1 um SFC re-calculation !jlee added
624 
625 c... **********************************************
626 c... **********************************************
627 c.... new addition: calculate 2.1um sfc refl. using
628 c... first guess of 470nm AOT
629 c... Aug. 12, 2011
630 
631  tau_x4701_00=tau_x470_1
632  xeas_b3_00=xeas_b3
633  tref21in=xnvalm6(4)
634  alpest=1.5
635  aot21est=tau_x470_1*(466.0/2110.)**alpest
636  if(aot21est.lt.0.0.or.aot21est.gt.0.5) return ! no retr. for too-heavy aerosol
637  if(aot21est.ge.0.0.and.aot21est.lt.0.06) go to 333
638 
639 c.
640  call calc_sfc21(x1,x2,x3,tref21in,aot21est,as_21)
641 c.
642 c... refl21r=3.14159*As_21/100.0/xmu
643  refl21r=as_21/100.0
644 c print *, 'refl21r: ', refl21r, As_21
645  if(refl21r.lt.0) return
646 
647  call get_sfcrfl_veg(season,lcvr,refl21r,rc_ndvi,xeas_b1,xeas_b3)
648  if(xeas_b1.lt.xmnsfc1.or.xeas_b1.gt.xmxsfc1.or.
649  1 xeas_b3.lt.xmnsfc3.or.xeas_b3.gt.xmxsfc3) return ! return to main w/o retrieval
650  r470=xeas_b3
651 c... **********************************************
652 c... **********************************************
653  call aero_470(dflag,refl,x1,x2,x3,mm,nn,ll,ma,
654  1 imod,r470,tau_x470_1,tau_x470_flag1,trflg,0.0,debug)
655 
656 333 continue
657 c... ********
658 c if(itmp.ge.800.and.itmp.lt.820.and.
659 c 1 jtmp.ge.900.and.jtmp.lt.920)
660 c 2 print *,100.*refl21,As_21,aot21est,tau_x4701_00,tau_x4701,
661 c 3 xeAs_B3_00, xeAs_B3
662 c... ********
663 
664  if(xeas_b1.lt.xmnsfc1.or.xeas_b1.gt.xmxsfc1.or.
665  1 xeas_b3.lt.xmnsfc3.or.xeas_b3.gt.xmxsfc3) return ! return to main w/o retrieval
666 
667 ! 16 January 2018 JLee
668  dflag = .false.
669  trflg = 1 ! default AOD = 0.02 if TOA reflectance < min(lut)
670  call aero_470(dflag,refl,x1,x2,x3,mm,nn,ll,ma,
671  1 imodarr(1),r470,tau_x470_1,tau_x470_flag1,trflg,modfrac(1),debug)
672  call aero_470(dflag,refl,x1,x2,x3,mm,nn,ll,ma,
673  1 imodarr(2),r470,tau_x470_2,tau_x470_flag2,trflg,modfrac(2),debug)
674  call aero_470(dflag,refl,x1,x2,x3,mm,nn,ll,ma,
675  1 imodarr(3),r470,tau_x470_3,tau_x470_flag3,trflg,modfrac(3),debug)
676 
677  tau_x470 = tau_x470_1
678  tau_x470_flag = tau_x470_flag1
679 
680  if(tau_x470_1.ge.cth0.and.tau_x470_1.lt.cth1) then
681  f_1=(tau_x470_1-cth0)/(cth1-cth0)
682  tau_x470=(1.0-f_1)*tau_x470_1+f_1*tau_x470_2
683  tau_x470_flag=tau_x470_flag2
684  endif
685 
686  if(tau_x470_1.ge.cth1.and.tau_x470_1.lt.cth2) then
687  f_1=(tau_x470_2-cth1)/(cth2-cth1)
688  tau_x470=(1.0-f_1)*tau_x470_2+f_1*tau_x470_3
689  tau_x470_flag=tau_x470_flag3
690  endif
691 
692  if(tau_x470_1.ge.cth2) then
693  tau_x470=tau_x470_3
694  tau_x470_flag=tau_x470_flag3
695  endif
696 
697  ! aod using 2.2 um database
698  if (do_sv.and.r470sv.gt.0.0.and.r470sv.lt.24.0) then
699  if (r470sv.lt.1.0) r470sv = 1.0
700  call aero_470(dflag,refl,x1,x2,x3,mm,nn,ll,ma,
701  1 imodarr(1),r470sv,tau_x470sv_1,tau_x470sv_flag1,trflg,modfrac(1),debug)
702  call aero_470(dflag,refl,x1,x2,x3,mm,nn,ll,ma,
703  1 imodarr(2),r470sv,tau_x470sv_2,tau_x470sv_flag2,trflg,modfrac(2),debug)
704  call aero_470(dflag,refl,x1,x2,x3,mm,nn,ll,ma,
705  1 imodarr(3),r470sv,tau_x470sv_3,tau_x470sv_flag3,trflg,modfrac(3),debug)
706 
707  tau_x470sv = tau_x470sv_1
708  tau_x470sv_flag = tau_x470sv_flag1
709 
710  if(tau_x470sv_1.ge.cth0.and.tau_x470sv_1.lt.cth1) then
711  f_1=(tau_x470sv_1-cth0)/(cth1-cth0)
712  tau_x470sv=(1.0-f_1)*tau_x470sv_1+f_1*tau_x470sv_2
713  tau_x470sv_flag=tau_x470sv_flag2
714  endif
715 
716  if(tau_x470sv_1.ge.cth1.and.tau_x470sv_1.lt.cth2) then
717  f_1=(tau_x470sv_2-cth1)/(cth2-cth1)
718  tau_x470sv=(1.0-f_1)*tau_x470sv_2+f_1*tau_x470sv_3
719  tau_x470sv_flag=tau_x470sv_flag3
720  endif
721 
722  if(tau_x470sv_1.ge.cth2) then
723  tau_x470sv=tau_x470sv_3
724  tau_x470sv_flag=tau_x470sv_flag3
725  endif
726  endif
727 ! end JLee added
728 
729 ! if(tau_x4701.ge.0.0.and.tau_x4701.lt.cth0) then
730 ! tau_x470=tau_x4701
731 ! tau_x470_flag=tau_x470_flag1
732 ! elseif(tau_x4701.ge.cth0.and.tau_x4701.lt.cth1) then
733 ! imod=imodarr(2)
734 ! call aero_470veg(refl,x1,x2,x3,mm,nn,ll,ma,
735 ! 1 imod,r470,tau_x4702,tau_x470_flag2)
736 ! f_1=(tau_x4701-cth0)/(cth1-cth0)
737 ! tau_x470=f_1*tau_x4702+(1.0-f_1)*tau_x4701
738 ! tau_x470_flag=tau_x470_flag2
739 ! elseif(tau_x4701.ge.cth1.and.tau_x4701.lt.cth2) then
740 ! imod=imodarr(2)
741 ! call aero_470veg(refl,x1,x2,x3,mm,nn,ll,ma,
742 ! 1 imod,r470,tau_x4702,tau_x470_flag2)
743 ! imod=imodarr(3)
744 ! call aero_470veg(refl,x1,x2,x3,mm,nn,ll,ma,
745 ! 1 imod,r470,tau_x4703,tau_x470_flag3)
746 ! f_2=(tau_x4701-cth1)/(cth2-cth1)
747 ! tau_x470=f_2*tau_x4703+(1.0-f_2)*tau_x4702
748 ! tau_x470_flag=tau_x470_flag3
749 ! elseif(tau_x4701.ge.cth2) then
750 ! imod=imodarr(3)
751 ! call aero_470veg(refl,x1,x2,x3,mm,nn,ll,ma,
752 ! 1 imod,r470,tau_x4703,tau_x470_flag3)
753 ! tau_x470=tau_x4703
754 ! tau_x470_flag=tau_x470_flag3
755 ! else
756 ! tau_x470=tau_x4701
757 ! tau_x470_flag=tau_x470_flag1
758 !c tau_x470=-999.0
759 ! endif
760 
761 c;--------------------------------------------<
762 c. original call for 470nm AOT retrieval (before 08/05/2011)
763 c. imod=3
764 c. call aero_470veg(refl,x1,x2,x3,mm,nn,ll,ma, ! commented out
765 c. 1 imod,r470,tau_x470,tau_x470_flag) ! 08/05/2011
766 
767 !c. @@@test@@@--> 5lines
768  if(tau_x470.gt.5.or.tau_x470.lt.0) then
769  tau_x470=-999.
770  tau_x650=-999.
771  return
772  endif
773 c print *, 'hereC'
774 c print *,'tau_x470 =', tau_x470
775 
776 
777 c--------------------------------------------------------
778 c Retrieving 650 nm AOT --> added Feb 26, 2010
779 c--------------------------------------------------------
780 c Input surface reflectance at 650 nm
781  r650 = xeas_b1
782 
783 !jlee added
784 ! if (test_type.eq.1) then
785 ! r650 = xeAs_B1
786 ! elseif (test_type.eq.2) then
787 ! r650 = xeAs_B1 - 4.67*tau_x470
788 ! elseif (test_type.eq.3) then
789 ! r650 = xeAs_B1
790 ! if (rc_ndvi.lt.0.4) then
791 ! r650 = xeAs_B1 - 4.67*tau_x470
792 ! endif
793 ! elseif (test_type.eq.4) then
794 ! r650 = xeAs_B1
795 ! if (rc_ndvi.lt.0.4.and.lcvr.eq.12) then
796 ! r650 = xeAs_B1 - 4.67*tau_x470
797 ! endif
798 ! else
799 ! endif
800 !end jlee added
801 c
802 c Input toa reflectance (L/F) at 650 nm
803  refl = refl6
804 
805 c Retrieving 650 nm AOT
806  call aero_650(dflag,refl,x1,x2,x3,mm,nn,ll,ma,
807  1 r650,tau_x650,tau_x650_flag,0,0.0,0.0,0,trflg)!tau_x470_flag,tau_x412,tau_x470,tau_x412_flag_91
808 
809 ! call aero_650veg(refl,x1,x2,x3,mm,nn,ll,ma,
810 ! 1 r650,tau_x650,tau_x650_flag,
811 ! 2 tau_x470_flag,tau_x470)
812 
813 c print *, 'hereD'
814 c print *,'tau_x650 =', tau_x650, r650
815 
816 !c. @@@test@@@--> 5lines
817  if(tau_x650.gt.5.or.tau_x650.lt.0) then
818  tau_x470=-999.
819  tau_x650=-999.
820  return
821  endif
822 
823 ! -- only set bright surface flag over N.America (ioprg==1) and
824 ! -- croplands (lcvr == 12)
825  if (ioprg .eq. 1 .and. lcvr .eq. 12) then
826  if (xeas_b1.gt.17.or.xeas_b3.gt.17) then
827  bflag = 1
828  endif
829  endif
830 
831 c--------------------------------------------------------
832 
833 c -- tweak AOT over S.Africa, Mongu site.
834  if (ioprg .eq. 6) then
835  if (season .eq .3) then
836  tau_x470 = tau_x470 * 1.25
837  tau_x650 = tau_x650 * 1.25
838  end if
839 ! if (month .ge. 9 .AND. month .le. 11) then
840 ! tau_x470 = tau_x470 * 1.2
841 ! tau_x650 = tau_x650 * 1.2
842 ! end if
843  end if
844 
845 c -- tweak over Mukdahan, SE Asia site.
846  if (ioprg .eq. 7) then
847 ! if (month .ge. 3 .AND. month .le. 5) then
848 ! if (tau_x470 .gt. 0.7) then
849 ! tau_x470 = tau_x470 * 1.1
850 ! tau_x650 = tau_x650 * 1.1
851 ! end if
852 ! end if
853 
854 ! if (season .eq. 1) then !DJF
855 ! tau_x470 = tau_x470 * 1.1
856 ! tau_x650 = tau_x650 * 1.1
857 ! end if
858 
859 ! if (month .ge. 9 .AND. month .le. 11) then
860 ! tau_x470 = tau_x470 * 1.1
861 ! tau_x650 = tau_x650 * 1.1
862 ! end if
863  end if
864 
865 ! jlee modified 650->672, 466->488
866 c... original set up (temporary)
867 c...==============================================
868 c... MJ re-defines output
869 
870  alpha=alog10(tau_x470/tau_x650)/alog10(672./488.)
871 
872 c if(tau_x470.lt.0.0) tau_x470=-999.0
873 c if(tau_x650.lt.0.0) tau_x650=-999.0
874 c if(tau_x470.lt.0.0.and.tau_x650.lt.0.0) alpha=-999.0
875 c if(alpha.ge.-0.5.and.alpha.le.3.5.and.tau_x470.gt.0.and.
876 c & tau_x650.gt.0.) then
877 
878  if (tau_x470.lt.1.0) then
879  if(alpha.ge.-0.5.and.alpha.le.3.5) then
880  tau550=tau_x470*(488.0/500.0)**alpha
881 c xtau(1)=tau_x470*(466.0/412.0)**alpha ! AOT 412nm (extrapol.)
882  else
883  if((tau_x470_flag.eq.1.or.tau_x470_flag.eq.-10).and.
884  1 (tau_x650_flag.eq.1.or.tau_x650_flag.eq.-10)) then
885  alpha=1.0
886  tau550=(tau_x470+tau_x650)/2.0
887  tau_x470=tau550*(500./488.)**alpha
888  tau_x650=tau550*(500./672.)**alpha
889  elseif((tau_x470_flag.eq.1.or.tau_x470_flag.eq.-10).and.
890  1 tau_x650_flag.eq.0) then
891  alpha=1.0
892  tau550=tau_x650*(672./500.)**alpha
893  tau_x470=tau_x650*(672./488.)**alpha
894  elseif(tau_x470_flag.eq.0.and.
895  1 (tau_x650_flag.eq.1.or.tau_x650_flag.eq.-10)) then
896  alpha=1.0
897  tau550=tau_x470*(488./500.)**alpha
898  tau_x650=tau_x470*(488./672.)**alpha
899  else
900  alpha=1.0
901  tau550=(tau_x470+tau_x650)/2.0
902  tau_x470=tau550*(500./488.)**alpha
903  tau_x650=tau550*(500./672.)**alpha
904  endif
905  endif
906  else ! if (tau_x470.lt.1.0) !start jlee added
907  if(alpha.ge.-0.5.and.alpha.le.3.5.and.
908  1 tau_x470_flag.eq.0.and.tau_x650_flag.eq.0) then
909  tau550=tau_x470*(488.0/500.0)**alpha
910 c xtau(1)=tau_x470*(466.0/412.0)**alpha ! AOT 412nm (extrapol.)
911  else
912  if((tau_x470_flag.eq.1.or.tau_x470_flag.eq.-10).and.
913  1 (tau_x650_flag.eq.1.or.tau_x650_flag.eq.-10)) then
914  alpha=1.0
915  tau550=(tau_x470+tau_x650)/2.0
916  tau_x470=tau550*(500./488.)**alpha
917  tau_x650=tau550*(500./672.)**alpha
918  elseif(tau_x470_flag.eq.-10.and.tau_x650_flag.eq.0) then
919  alpha=1.0
920  tau550=tau_x650*(672./500.)**alpha
921  tau_x470=tau_x650*(672./488.)**alpha
922  elseif((tau_x470_flag.eq.1.or.tau_x470_flag.eq.0).and.
923  1 tau_x650_flag.eq.0) then
924  tau550=tau_x470
925  alpha = 1.8 !thick smoke
926  tau_x470=tau550*(500./488.)**alpha
927  tau_x650=tau550*(500./672.)**alpha
928  elseif(tau_x470_flag.eq.0.and.
929  1 (tau_x650_flag.eq.1.or.tau_x650_flag.eq.-10)) then
930  alpha=1.0
931  tau550=tau_x470*(488./500.)**alpha
932  tau_x650=tau_x470*(488./672.)**alpha
933  else
934  alpha=1.0
935  tau550=(tau_x470+tau_x650)/2.0
936  tau_x470=tau550*(500./488.)**alpha
937  tau_x650=tau550*(500./672.)**alpha
938  endif
939  endif
940  endif !if (tau_x470.lt.1.5) ! end jlee added
941 
942  if (do_sv) then
943  tau550=tau_x470sv*(488./500.)**alpha
944  tau_x470 = tau_x470sv
945  tau_x470_flag = tau_x470sv_flag
946  endif
947 c...==============================================
948 !end jlee modified
949 
950 c...==========================================================
951 c... MJ re-defines output
952 c if(tau_x470.lt.0.0) tau_x470=-999.0
953 c if(tau_x650.lt.0.0) tau_x650=-999.0
954 c if(tau_x470.lt.0.0.or.tau_x650.lt.0.0) return
955 c-------------------------------------------------------------
956 c Set output buffer --> based on find_v.f
957 c-------------------------------------------------------------
958 ! 22 December 2017 decided not to report Veg SSA
959 ! To do: report SSA as a diagnostic variable in the extended output
960 ! read(w0_name470(imod),'(f4.2)') ssatmp
961 ! ssa(2)=ssatmp
962 ! ssa(3)=0.976
963  xtau(1)=-999.0 ! fillvalue at 412nm for over-veg retr.
964  xtau(2)=tau_x470
965  xtau(3)=tau_x650
966  do i=1,3
967  outbufvg(i) = xtau(i)
968  outbufvg(i+3) = ssa(i)
969  enddo
970  outbufvg(7) = tau550
971  outbufvg(8) = alpha
972  outbufvg(9) = 1.0*tau_x470_flag ! used to be "r412"
973  outbufvg(10) = 1.0*tau_x650_flag
974  outbufvg(11) = r470 ! over-veg. sfc. refl.
975  outbufvg(12) = r650 ! over-veg. sfc. refl.
976  outbufvg(13) = xthet
977  outbufvg(14) = scat_ang
978  outbufvg(15) = sfc_typ
979 
980  outbufvg(16) = xlat
981  outbufvg(17) = xlong
982  outbufvg(18) = sirndvi
983  outbufvg(19) = rc_ndvi
984  outbufvg(20) = xreg_id
985  outbufvg(21) = bflag
986 
987 c print *,'outbufvg(1,2,3) ',outbufvg(1),outbufvg(2),outbufvg(3)
988  return
989  end
990 c... --------------------------------------
991 c... The end of subroutine myd_aero_veg
992 c... (main subroutine for aerosol retr. over vegetation)
993 c... ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
994 
995 c... ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
996  subroutine get_sfcrfl_veg(season,lcvr,refln21,rc_ndvi,
997  1 xeAs_B1,xeAs_B3)
998 c...
999 c... subroutine to estimate surface reflectance at 470nm & 650nm over
1000 c... vegetated surfaces using 2.1um TOA reflectance
1001 c... program written by Myeong-Jae Jeong (MJ)
1002 c... last modified Jun 15, 2010
1003 c...
1004 c... Inputs
1005 c... season: season ID (1: MAM; 2: JJA; 3: SON; 4: DJF)
1006 c... refln21: 2.1um TOA reflectance (normalized, Ipi/fcos(sza))
1007 c... xnvalm6: TOA reflectances at other bands
1008 c... lcvr: MODIS land cover (IGBP)
1009 c...
1010 c... Output
1011 c... xeAs_B1: estimated surface refl. for band 1 (LER in %)
1012 c... xeAs_B3: estimated surface refl. for band 3 (LER in %)
1013 c...
1014 c... 02/11/2014, VIIRS surface parameterization for general vegetation areas, modified by J. Lee
1015 
1016  integer season
1017 
1018  sb7=100.0*refln21
1019 
1020  if(season.eq.2) then ! MAM 2012
1021 
1022  if((lcvr.gt.0.and.lcvr.le.12).or.lcvr.eq.14) then
1023  xeas_b1 = 0.283 + 0.509*sb7 + 0.0038*sb7**2.0
1024  sb1 = xeas_b1
1025  xeas_b3 = 0.810 + 0.535*sb1
1026  elseif(lcvr.eq.13) then
1027  xeas_b1 = 0.762*sb7
1028  sb1 = xeas_b1
1029  xeas_b3 = 0.718*sb1
1030  else
1031  return
1032  endif
1033 
1034  elseif(season.eq.3) then ! JJA 2012
1035 
1036  if((lcvr.gt.0.and.lcvr.le.12).or.lcvr.eq.14) then
1037  xeas_b1 = 0.295 + 0.405*sb7 + 0.0095*sb7**2.0
1038  sb1 = xeas_b1
1039  xeas_b3 = 0.605 + 0.507*sb1
1040  elseif(lcvr.eq.13) then
1041  xeas_b1 = 0.765*sb7
1042  sb1 = xeas_b1
1043  xeas_b3 = 0.664*sb1
1044  else
1045  return
1046  endif
1047 
1048  elseif(season.eq.4) then ! SON 2012
1049 
1050  if((lcvr.gt.0.and.lcvr.le.12).or.lcvr.eq.14) then
1051  xeas_b1 = 0.298 + 0.430*sb7 + 0.0084*sb7**2.0
1052  sb1 = xeas_b1
1053  xeas_b3 = 0.494 + 0.524*sb1
1054  elseif(lcvr.eq.13) then
1055  xeas_b1 = 0.784*sb7
1056  sb1 = xeas_b1
1057  xeas_b3 = 0.665*sb1
1058  else
1059  return
1060  endif
1061 
1062  elseif(season.eq.1) then ! DJF 2013
1063 
1064  if((lcvr.gt.0.and.lcvr.le.12).or.lcvr.eq.14) then
1065  xeas_b1 = 0.283 + 0.509*sb7 + 0.0038*sb7**2.0
1066  sb1 = xeas_b1
1067  xeas_b3 = 0.810 + 0.535*sb1
1068  elseif(lcvr.eq.13) then
1069  xeas_b1 = 0.762*sb7
1070  sb1 = xeas_b1
1071  xeas_b3 = 0.718*sb1
1072  else
1073  return
1074  endif
1075 
1076  else
1077  return
1078  endif
1079 
1080  return
1081  end
1082 c... --------------------------------------
1083 c... The end of subroutine get_sfcrfl_veg
1084 c... ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1085  subroutine get_sfcrfl_veg_high(season,lcvr,refln21,rc_ndvi,
1086  1 xeAs_B1,xeAs_B3)
1087 c...
1088 c... subroutine to estimate surface reflectance at 470nm & 650nm over
1089 c... vegetated surfaces using 2.1um TOA reflectance
1090 c... program written by Myeong-Jae Jeong (MJ)
1091 c... last modified Jun 15, 2010
1092 c...
1093 c... Inputs
1094 c... season: season ID (1: DJF; 2: MAM; 3: JJA; 4: SON)
1095 c... refln21: 2.1um TOA reflectance
1096 c... xnvalm6: TOA reflectances at other bands
1097 c... lcvr: MODIS land cover (IGBP)
1098 c...
1099 c... Output
1100 c... xeAs_B1: estimated surface refl. for band 1 (LER in %)
1101 c... xeAs_B3: estimated surface refl. for band 3 (LER in %)
1102 c...
1103 c... 02/11/2014, VIIRS surface parameterization for general vegetation areas, modified by J. Lee
1104 
1105  integer season
1106 
1107  sb7=100.0*refln21
1108 
1109  if(season.eq.2) then ! MAM 2012
1110 
1111  if((lcvr.gt.0.and.lcvr.le.12).or.lcvr.eq.14) then
1112  xeas_b1 = -0.542 +0.575*sb7 +0.0022*sb7**2
1113  sb1 = xeas_b1
1114  xeas_b3 = +0.342 +0.598*sb1
1115  elseif(lcvr.eq.13) then
1116  xeas_b1 = -1.201 +0.785*sb7
1117  sb1 = xeas_b1
1118  xeas_b3 = +0.401 +0.679*sb1
1119  else
1120  return
1121  endif
1122 
1123  elseif(season.eq.3) then ! JJA 2012
1124 
1125  if((lcvr.gt.0.and.lcvr.le.12).or.lcvr.eq.14) then
1126  xeas_b1 = -0.039 +0.401*sb7 +0.0108*sb7**2
1127  sb1 = xeas_b1
1128  xeas_b3 = +0.568 +0.572*sb1
1129  elseif(lcvr.eq.13) then
1130  xeas_b1 = -2.147 +0.976*sb7
1131  sb1 = xeas_b1
1132  xeas_b3 = +0.713 +0.620*sb1
1133  else
1134  return
1135  endif
1136 
1137  elseif(season.eq.4) then ! SON 2012
1138 
1139  if((lcvr.gt.0.and.lcvr.le.12).or.lcvr.eq.14) then
1140  xeas_b1 = -0.125 +0.466*sb7 +0.0076*sb7**2
1141  sb1 = xeas_b1
1142  xeas_b3 = +0.438 +0.518*sb1
1143  elseif(lcvr.eq.13) then
1144  xeas_b1 = -1.428 +0.839*sb7
1145  sb1 = xeas_b1
1146  xeas_b3 = +0.233 +0.648*sb1
1147  else
1148  return
1149  endif
1150 
1151  elseif(season.eq.1) then ! DJF 2013
1152 
1153  if((lcvr.gt.0.and.lcvr.le.12).or.lcvr.eq.14) then
1154  xeas_b1 = -0.542 +0.575*sb7 +0.0022*sb7**2
1155  sb1 = xeas_b1
1156  xeas_b3 = +0.342 +0.598*sb1
1157  elseif(lcvr.eq.13) then
1158  xeas_b1 = -1.201 +0.785*sb7
1159  sb1 = xeas_b1
1160  xeas_b3 = +0.401 +0.679*sb1
1161  else
1162  return
1163  endif
1164 
1165  else
1166  return
1167  endif
1168 
1169  return
1170  end
1171 c... --------------------------------------
1172 c... The end of subroutine get_sfcrfl_veg_high
1173 c... ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1174  subroutine get_sfcrfl_swir_vis(season,lcvr,refln21,rc_ndvi,
1175  1 xeAs_B1,xeAs_B3)
1176 c...
1177 c... subroutine to estimate surface reflectance at 470nm & 650nm over
1178 c... vegetated surfaces using 2.25 um TOA reflectance and Rayleigh-corrected NDVI
1179 c... program written by Jaehwa Lee
1180 c... last modified March 24, 2014
1181 c...
1182 c... Inputs
1183 c... season: season ID (1: DJF; 2: MAM; 3: JJA; 4: SON)
1184 c... lcvr: MODIS land cover (IGBP)
1185 c... refln21: 2.1um TOA reflectance
1186 c... rc_ndvi: Rayleigh-corrected NDVI
1187 c...
1188 c... Output
1189 c... xeAs_B1: estimated surface refl. for band 1 (LER in %)
1190 c... xeAs_B3: estimated surface refl. for band 3 (LER in %)
1191 c...
1192 c... 03/24/2014, Coeffs. w.r.t. Rayleigh-corrected NDVI are applied.
1193  integer season
1194 
1195  real, dimension(20) :: coeff0_red
1196  real, dimension(20) :: coeff1_red
1197  real, dimension(20) :: coeff0_blue
1198  real, dimension(20) :: coeff1_blue
1199 
1200  rc_ndvi = max(rc_ndvi, 0.1)
1201  rc_ndvi = min(rc_ndvi, 1.0)
1202  sb7=100.0*refln21
1203 
1204  if(season.eq.1.or.season.eq.2) then ! DJF, MAM
1205 
1206  if(lcvr.gt.0.and.lcvr.ne.7.and.lcvr.ne.13) then
1207  coeff0_red(1:20) = (/ 0.000, 0.000, 0.000, 0.000, 0.000,
1208  1 0.000, 0.000, 0.000, 0.000, 0.000,
1209  1 0.000, 0.000, 0.000, 0.000, 0.000,
1210  1 0.000, 0.000, 0.000, 0.000, 0.000/)
1211  coeff1_red(1:20) = (/ 0.688, 0.688, 0.688, 0.679, 0.632,
1212  1 0.604, 0.592, 0.578, 0.561, 0.538,
1213  1 0.519, 0.508, 0.502, 0.491, 0.468,
1214  1 0.426, 0.382, 0.337, 0.000, 0.000/)
1215  coeff0_blue(1:20) =(/ 0.000, 0.000, 0.000, 0.000, 0.000,
1216  1 0.000, 0.000, 0.000, 0.000, 0.000,
1217  1 0.000, 0.000, 0.000, 0.000, 0.000,
1218  1 0.000, 0.000, 0.000, 0.000, 0.000/)
1219  coeff1_blue(1:20) =(/ 0.522, 0.522, 0.522, 0.534, 0.530,
1220  1 0.551, 0.567, 0.596, 0.618, 0.633,
1221  1 0.651, 0.676, 0.703, 0.734, 0.771,
1222  1 0.818, 0.885, 0.906, 0.000, 0.000/)
1223 
1224  ndvi_idx = min(floor(rc_ndvi*20)+1, 18)
1225 
1226  xeas_b1 = coeff0_red(ndvi_idx) + coeff1_red(ndvi_idx)*sb7
1227  sb1 = xeas_b1
1228  xeas_b3 = coeff0_blue(ndvi_idx) + coeff1_blue(ndvi_idx)*sb1
1229 
1230 ! ndvi_idx = min(floor((rc_ndvi-0.025)*20+1), 18)
1231 ! coeff0_red_fin = coeff0_red(ndvi_idx)
1232 ! 1 +(coeff0_red(ndvi_idx+1)-coeff0_red(ndvi_idx))/0.05
1233 ! 1 *(rc_ndvi-(0.025+(ndvi_idx-1)*0.05))
1234 ! coeff1_red_fin = coeff1_red(ndvi_idx)
1235 ! 1 +(coeff1_red(ndvi_idx+1)-coeff1_red(ndvi_idx))/0.05
1236 ! 1 *(rc_ndvi-(0.025+(ndvi_idx-1)*0.05))
1237 ! coeff0_blue_fin = coeff0_blue(ndvi_idx)
1238 ! 1 +(coeff0_blue(ndvi_idx+1)-coeff0_blue(ndvi_idx))/0.05
1239 ! 1 *(rc_ndvi-(0.025+(ndvi_idx-1)*0.05))
1240 ! coeff1_blue_fin = coeff1_blue(ndvi_idx)
1241 ! 1 +(coeff1_blue(ndvi_idx+1)-coeff1_blue(ndvi_idx))/0.05
1242 ! 1 *(rc_ndvi-(0.025+(ndvi_idx-1)*0.05))
1243 !
1244 ! xeAs_B1 = coeff0_red_fin + coeff1_red_fin*sb7
1245 ! sb1 = xeAs_B1
1246 ! xeAs_B3 = coeff0_blue_fin + coeff1_blue_fin*sb1
1247 
1248  elseif(lcvr.eq.7) then !open shrublands
1249  coeff0_red(1:20) = (/ 0.000, 0.000, 0.000, 0.000, 0.000,
1250  1 0.000, 0.000, 0.000, 0.000, 0.000,
1251  1 0.000, 0.000, 0.000, 0.000, 0.000,
1252  1 0.000, 0.000, 0.000, 0.000, 0.000/)
1253  coeff1_red(1:20) = (/ 0.735, 0.735, 0.735, 0.724, 0.716,
1254  1 0.709, 0.692, 0.647, 0.601, 0.565,
1255  1 0.555, 0.535, 0.537, 0.513, 0.498,
1256  1 0.000, 0.000, 0.000, 0.000, 0.000/)
1257  coeff0_blue(1:20) =(/ 0.000, 0.000, 0.000, 0.000, 0.000,
1258  1 0.000, 0.000, 0.000, 0.000, 0.000,
1259  1 0.000, 0.000, 0.000, 0.000, 0.000,
1260  1 0.000, 0.000, 0.000, 0.000, 0.000/)
1261  coeff1_blue(1:20) =(/ 0.527, 0.527, 0.527, 0.508, 0.526,
1262  1 0.540, 0.555, 0.555, 0.606, 0.631,
1263  1 0.744, 0.800, 0.809, 0.810, 0.813,
1264  1 0.000, 0.000, 0.000, 0.000, 0.000/)
1265 
1266  ndvi_idx = min(floor(rc_ndvi*20)+1, 15)
1267 
1268  xeas_b1 = coeff0_red(ndvi_idx) + coeff1_red(ndvi_idx)*sb7
1269  sb1 = xeas_b1
1270  xeas_b3 = coeff0_blue(ndvi_idx) + coeff1_blue(ndvi_idx)*sb1
1271 
1272  elseif(lcvr.eq.13) then !urban and built-up
1273  coeff0_red(1:20) = (/ 0.000, 0.000, 0.000, 0.000, 0.000,
1274  1 0.000, 0.000, 0.000, 0.000, 0.000,
1275  1 0.000, 0.000, 0.000, 0.000, 0.000,
1276  1 0.000, 0.000, 0.000, 0.000, 0.000/)
1277  coeff1_red(1:20) = (/ 0.903, 0.903, 0.903, 0.860, 0.818,
1278  1 0.751, 0.721, 0.677, 0.642, 0.627,
1279  1 0.604, 0.581, 0.561, 0.530, 0.483,
1280  1 0.436, 0.404, 0.000, 0.000, 0.000/)
1281  coeff0_blue(1:20) =(/ 0.000, 0.000, 0.000, 0.000, 0.000,
1282  1 0.000, 0.000, 0.000, 0.000, 0.000,
1283  1 0.000, 0.000, 0.000, 0.000, 0.000,
1284  1 0.000, 0.000, 0.000, 0.000, 0.000/)
1285  coeff1_blue(1:20) =(/ 0.693, 0.693, 0.693, 0.729, 0.735,
1286  1 0.725, 0.728, 0.725, 0.729, 0.738,
1287  1 0.749, 0.768, 0.793, 0.805, 0.818,
1288  1 0.838, 0.862, 0.000, 0.000, 0.000/)
1289 
1290  ndvi_idx = min(floor(rc_ndvi*20)+1, 17)
1291 
1292  xeas_b1 = coeff0_red(ndvi_idx) + coeff1_red(ndvi_idx)*sb7
1293  sb1 = xeas_b1
1294  xeas_b3 = coeff0_blue(ndvi_idx) + coeff1_blue(ndvi_idx)*sb1
1295 
1296  else
1297  return
1298  endif
1299 
1300  elseif(season.eq.3) then ! JJA
1301 
1302  if(lcvr.gt.0.and.lcvr.ne.7.and.lcvr.ne.13) then
1303  coeff0_red(1:20) = (/ 0.000, 0.000, 0.000, 0.000, 0.000,
1304  1 0.000, 0.000, 0.000, 0.000, 0.000,
1305  1 0.000, 0.000, 0.000, 0.000, 0.000,
1306  1 0.000, 0.000, 0.000, 0.000, 0.000/)
1307  coeff1_red(1:20) = (/ 0.692, 0.692, 0.692, 0.681, 0.684,
1308  1 0.653, 0.626, 0.602, 0.593, 0.580,
1309  1 0.564, 0.540, 0.519, 0.496, 0.471,
1310  1 0.441, 0.395, 0.332, 0.260, 0.000/)
1311  coeff0_blue(1:20) =(/ 0.000, 0.000, 0.000, 0.000, 0.000,
1312  1 0.000, 0.000, 0.000, 0.000, 0.000,
1313  1 0.000, 0.000, 0.000, 0.000, 0.000,
1314  1 0.000, 0.000, 0.000, 0.000, 0.000/)
1315  coeff1_blue(1:20) =(/ 0.522, 0.522, 0.522, 0.516, 0.490,
1316  1 0.487, 0.499, 0.516, 0.530, 0.541,
1317  1 0.556, 0.583, 0.620, 0.676, 0.767,
1318  1 0.827, 0.899, 0.926, 0.665, 0.000/)
1319 
1320  ndvi_idx = min(floor(rc_ndvi*20)+1, 18)
1321 
1322  xeas_b1 = coeff0_red(ndvi_idx) + coeff1_red(ndvi_idx)*sb7
1323  sb1 = xeas_b1
1324  xeas_b3 = coeff0_blue(ndvi_idx) + coeff1_blue(ndvi_idx)*sb1
1325 
1326  elseif(lcvr.eq.7) then
1327  coeff0_red(1:20) = (/ 0.000, 0.000, 0.000, 0.000, 0.000,
1328  1 0.000, 0.000, 0.000, 0.000, 0.000,
1329  1 0.000, 0.000, 0.000, 0.000, 0.000,
1330  1 0.000, 0.000, 0.000, 0.000, 0.000/)
1331  coeff1_red(1:20) = (/ 0.735, 0.735, 0.735, 0.724, 0.716,
1332  1 0.709, 0.692, 0.647, 0.601, 0.565,
1333  1 0.555, 0.535, 0.537, 0.513, 0.498,
1334  1 0.000, 0.000, 0.000, 0.000, 0.000/)
1335  coeff0_blue(1:20) =(/ 0.000, 0.000, 0.000, 0.000, 0.000,
1336  1 0.000, 0.000, 0.000, 0.000, 0.000,
1337  1 0.000, 0.000, 0.000, 0.000, 0.000,
1338  1 0.000, 0.000, 0.000, 0.000, 0.000/)
1339  coeff1_blue(1:20) =(/ 0.527, 0.527, 0.527, 0.508, 0.526,
1340  1 0.540, 0.555, 0.555, 0.606, 0.631,
1341  1 0.744, 0.800, 0.809, 0.810, 0.813,
1342  1 0.000, 0.000, 0.000, 0.000, 0.000/)
1343 
1344  ndvi_idx = min(floor(rc_ndvi*20)+1, 15)
1345 
1346  xeas_b1 = coeff0_red(ndvi_idx) + coeff1_red(ndvi_idx)*sb7
1347  sb1 = xeas_b1
1348  xeas_b3 = coeff0_blue(ndvi_idx) + coeff1_blue(ndvi_idx)*sb1
1349 
1350  elseif(lcvr.eq.13) then !urban and built-up
1351  coeff0_red(1:20) = (/ 0.000, 0.000, 0.000, 0.000, 0.000,
1352  1 0.000, 0.000, 0.000, 0.000, 0.000,
1353  1 0.000, 0.000, 0.000, 0.000, 0.000,
1354  1 0.000, 0.000, 0.000, 0.000, 0.000/)
1355  coeff1_red(1:20) = (/ 0.901, 0.901, 0.901, 0.874, 0.844,
1356  1 0.811, 0.802, 0.778, 0.742, 0.708,
1357  1 0.677, 0.646, 0.615, 0.574, 0.526,
1358  1 0.481, 0.427, 0.362, 0.000, 0.000/)
1359  coeff0_blue(1:20) =(/ 0.000, 0.000, 0.000, 0.000, 0.000,
1360  1 0.000, 0.000, 0.000, 0.000, 0.000,
1361  1 0.000, 0.000, 0.000, 0.000, 0.000,
1362  1 0.000, 0.000, 0.000, 0.000, 0.000/)
1363  coeff1_blue(1:20) =(/ 0.653, 0.653, 0.653, 0.688, 0.698,
1364  1 0.692, 0.699, 0.700, 0.701, 0.695,
1365  1 0.710, 0.770, 0.819, 0.847, 0.865,
1366  1 0.895, 0.951, 1.045, 0.000, 0.000/)
1367 
1368  ndvi_idx = min(floor(rc_ndvi*20)+1, 18)
1369 
1370  xeas_b1 = coeff0_red(ndvi_idx) + coeff1_red(ndvi_idx)*sb7
1371  sb1 = xeas_b1
1372  xeas_b3 = coeff0_blue(ndvi_idx) + coeff1_blue(ndvi_idx)*sb1
1373 
1374  else
1375  return
1376  endif
1377 
1378  elseif(season.eq.4) then ! SON
1379 
1380  if(lcvr.gt.0.and.lcvr.ne.7.and.lcvr.ne.13) then
1381  coeff0_red(1:20) = (/ 0.000, 0.000, 0.000, 0.000, 0.000,
1382  1 0.000, 0.000, 0.000, 0.000, 0.000,
1383  1 0.000, 0.000, 0.000, 0.000, 0.000,
1384  1 0.000, 0.000, 0.000, 0.000, 0.000/)
1385  coeff1_red(1:20) = (/ 0.686, 0.686, 0.686, 0.677, 0.610,
1386  1 0.572, 0.564, 0.559, 0.549, 0.534,
1387  1 0.526, 0.517, 0.511, 0.500, 0.475,
1388  1 0.439, 0.397, 0.337, 0.238, 0.000/)
1389  coeff0_blue(1:20) =(/ 0.000, 0.000, 0.000, 0.000, 0.000,
1390  1 0.000, 0.000, 0.000, 0.000, 0.000,
1391  1 0.000, 0.000, 0.000, 0.000, 0.000,
1392  1 0.000, 0.000, 0.000, 0.000, 0.000/)
1393  coeff1_blue(1:20) =(/ 0.506, 0.506, 0.506, 0.489, 0.485,
1394  1 0.490, 0.497, 0.519, 0.560, 0.587,
1395  1 0.605, 0.634, 0.668, 0.703, 0.736,
1396  1 0.782, 0.843, 0.859, 0.622, 0.000/)
1397 
1398  ndvi_idx = min(floor(rc_ndvi*20)+1, 18)
1399 
1400  xeas_b1 = coeff0_red(ndvi_idx) + coeff1_red(ndvi_idx)*sb7
1401  sb1 = xeas_b1
1402  xeas_b3 = coeff0_blue(ndvi_idx) + coeff1_blue(ndvi_idx)*sb1
1403 
1404  elseif(lcvr.eq.7) then
1405  coeff0_red(1:20) = (/ 0.000, 0.000, 0.000, 0.000, 0.000,
1406  1 0.000, 0.000, 0.000, 0.000, 0.000,
1407  1 0.000, 0.000, 0.000, 0.000, 0.000,
1408  1 0.000, 0.000, 0.000, 0.000, 0.000/)
1409  coeff1_red(1:20) = (/ 0.735, 0.735, 0.735, 0.724, 0.716,
1410  1 0.709, 0.692, 0.647, 0.601, 0.565,
1411  1 0.555, 0.535, 0.537, 0.513, 0.498,
1412  1 0.000, 0.000, 0.000, 0.000, 0.000/)
1413  coeff0_blue(1:20) =(/ 0.000, 0.000, 0.000, 0.000, 0.000,
1414  1 0.000, 0.000, 0.000, 0.000, 0.000,
1415  1 0.000, 0.000, 0.000, 0.000, 0.000,
1416  1 0.000, 0.000, 0.000, 0.000, 0.000/)
1417  coeff1_blue(1:20) =(/ 0.527, 0.527, 0.527, 0.508, 0.526,
1418  1 0.540, 0.555, 0.555, 0.606, 0.631,
1419  1 0.744, 0.800, 0.809, 0.810, 0.813,
1420  1 0.000, 0.000, 0.000, 0.000, 0.000/)
1421 
1422  ndvi_idx = min(floor(rc_ndvi*20)+1, 15)
1423 
1424  xeas_b1 = coeff0_red(ndvi_idx) + coeff1_red(ndvi_idx)*sb7
1425  sb1 = xeas_b1
1426  xeas_b3 = coeff0_blue(ndvi_idx) + coeff1_blue(ndvi_idx)*sb1
1427 
1428  elseif(lcvr.eq.13) then !urban and built-up
1429  coeff0_red(1:20) = (/ 0.000, 0.000, 0.000, 0.000, 0.000,
1430  1 0.000, 0.000, 0.000, 0.000, 0.000,
1431  1 0.000, 0.000, 0.000, 0.000, 0.000,
1432  1 0.000, 0.000, 0.000, 0.000, 0.000/)
1433  coeff1_red(1:20) = (/ 0.854, 0.854, 0.854, 0.846, 0.824,
1434  1 0.780, 0.746, 0.699, 0.664, 0.647,
1435  1 0.632, 0.606, 0.576, 0.540, 0.510,
1436  1 0.481, 0.425, 0.000, 0.000, 0.000/)
1437  coeff0_blue(1:20) =(/ 0.000, 0.000, 0.000, 0.000, 0.000,
1438  1 0.000, 0.000, 0.000, 0.000, 0.000,
1439  1 0.000, 0.000, 0.000, 0.000, 0.000,
1440  1 0.000, 0.000, 0.000, 0.000, 0.000/)
1441  coeff1_blue(1:20) =(/ 0.700, 0.700, 0.700, 0.729, 0.722,
1442  1 0.709, 0.709, 0.711, 0.711, 0.732,
1443  1 0.751, 0.769, 0.789, 0.808, 0.829,
1444  1 0.858, 0.917, 0.000, 0.000, 0.000/)
1445 
1446  ndvi_idx = min(floor(rc_ndvi*20)+1, 17)
1447 
1448  xeas_b1 = coeff0_red(ndvi_idx) + coeff1_red(ndvi_idx)*sb7
1449  sb1 = xeas_b1
1450  xeas_b3 = coeff0_blue(ndvi_idx) + coeff1_blue(ndvi_idx)*sb1
1451 
1452  else
1453  return
1454  endif
1455 
1456  else
1457  return
1458  endif
1459 
1460  return
1461  end
1462 
1463 c... --------------------------------------
1464 c... The end of subroutine get_sfcrfl_swir_vis
1465 c... ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1466 
1467  subroutine get_sfcrfl_nir_vis(season,lcvr,refn865_rc,rc_ndvi,
1468  1 xeAs_B1,xeAs_B3)
1469 c...
1470 c... subroutine to estimate surface reflectance at 470nm & 650nm over
1471 c... vegetated surfaces using 0.86 um TOA reflectance
1472 c... program written by Jaehwa Lee
1473 c... last modified March 24, 2014
1474 c...
1475 c... Inputs
1476 c... season: season ID (1: DJF; 2: MAM; 3: JJA; 4: SON)
1477 c... lcvr: MODIS land cover (IGBP)
1478 c... refn865_rc: Rayleigh-corrected 865 nm reflectance
1479 c... rc_ndvi: Rayleigh-corrected NDVI
1480 c...
1481 c... Output
1482 c... xeAs_B1: estimated surface refl. for band 1 (LER in %)
1483 c... xeAs_B3: estimated surface refl. for band 3 (LER in %)
1484 c...
1485 c... 03/24/2014, Coeffs. w.r.t. Rayleigh-corrected NDVI are applied.
1486  integer season
1487 
1488  real, dimension(20) :: coeff0_red
1489  real, dimension(20) :: coeff1_red
1490  real, dimension(20) :: coeff0_blue
1491  real, dimension(20) :: coeff1_blue
1492 
1493  rc_ndvi = max(rc_ndvi, 0.1)
1494  rc_ndvi = min(rc_ndvi, 1.0)
1495  sb2=100.0*refn865_rc
1496 
1497  if(season.eq.1.or.season.eq.2) then ! MAM, DJF
1498 
1499  if(lcvr.gt.0.and.lcvr.ne.13) then
1500  coeff0_red(1:20) = (/ 0.000, 0.000, 0.000, 0.000, 0.000,
1501  1 0.000, 0.000, 0.000, 0.000, 0.000,
1502  1 0.000, 0.000, 0.000, 0.000, 0.000,
1503  1 0.000, 0.000, 0.000, 0.000, 0.000/)
1504  coeff1_red(1:20) = (/ 0.770, 0.770, 0.770, 0.681, 0.617,
1505  1 0.548, 0.487, 0.430, 0.378, 0.330,
1506  1 0.286, 0.244, 0.207, 0.171, 0.135,
1507  1 0.102, 0.075, 0.056, 0.000, 0.000/)
1508  coeff0_blue(1:20) =(/ 0.000, 0.000, 0.000, 0.000, 0.000,
1509  1 0.000, 0.000, 0.000, 0.000, 0.000,
1510  1 0.000, 0.000, 0.000, 0.000, 0.000,
1511  1 0.000, 0.000, 0.000, 0.000, 0.000/)
1512  coeff1_blue(1:20) =(/ 0.406, 0.406, 0.406, 0.379, 0.344,
1513  1 0.308, 0.277, 0.251, 0.226, 0.204,
1514  1 0.183, 0.163, 0.144, 0.124, 0.102,
1515  1 0.081, 0.064, 0.050, 0.000, 0.000/)
1516 
1517  ndvi_idx = min(floor((rc_ndvi-0.025)*20+1), 18)
1518  coeff0_red_fin = coeff0_red(ndvi_idx)
1519  1 +(coeff0_red(ndvi_idx+1)-coeff0_red(ndvi_idx))/0.05
1520  1 *(rc_ndvi-(ndvi_idx*0.05-0.025))
1521  coeff1_red_fin = coeff1_red(ndvi_idx)
1522  1 +(coeff1_red(ndvi_idx+1)-coeff1_red(ndvi_idx))/0.05
1523  1 *(rc_ndvi-(ndvi_idx*0.05-0.025))
1524  coeff0_blue_fin = coeff0_blue(ndvi_idx)
1525  1 +(coeff0_blue(ndvi_idx+1)-coeff0_blue(ndvi_idx))/0.05
1526  1 *(rc_ndvi-(ndvi_idx*0.05-0.025))
1527  coeff1_blue_fin = coeff1_blue(ndvi_idx)
1528  1 +(coeff1_blue(ndvi_idx+1)-coeff1_blue(ndvi_idx))/0.05
1529  1 *(rc_ndvi-(ndvi_idx*0.05-0.025))
1530 
1531  xeas_b1 = coeff0_red_fin + coeff1_red_fin*sb2
1532  xeas_b3 = coeff0_blue_fin + coeff1_blue_fin*sb2
1533  elseif(lcvr.eq.13) then
1534  coeff0_red(1:20) = (/ 0.000, 0.000, 0.000, 0.000, 0.000,
1535  1 0.000, 0.000, 0.000, 0.000, 0.000,
1536  1 0.000, 0.000, 0.000, 0.000, 0.000,
1537  1 0.000, 0.000, 0.000, 0.000, 0.000/)
1538  coeff1_red(1:20) = (/ 0.739, 0.739, 0.739, 0.673, 0.604,
1539  1 0.537, 0.478, 0.422, 0.372, 0.326,
1540  1 0.282, 0.236, 0.197, 0.162, 0.129,
1541  1 0.094, 0.000, 0.000, 0.000, 0.000/)
1542  coeff0_blue(1:20) =(/ 0.000, 0.000, 0.000, 0.000, 0.000,
1543  1 0.000, 0.000, 0.000, 0.000, 0.000,
1544  1 0.000, 0.000, 0.000, 0.000, 0.000,
1545  1 0.000, 0.000, 0.000, 0.000, 0.000/)
1546  coeff1_blue(1:20) =(/ 0.516, 0.516, 0.516, 0.482, 0.438,
1547  1 0.386, 0.344, 0.304, 0.269, 0.240,
1548  1 0.212, 0.184, 0.157, 0.130, 0.103,
1549  1 0.080, 0.000, 0.000, 0.000, 0.000/)
1550 
1551  ndvi_idx = min(floor((rc_ndvi-0.025)*20+1), 16)
1552  coeff0_red_fin = coeff0_red(ndvi_idx)
1553  1 +(coeff0_red(ndvi_idx+1)-coeff0_red(ndvi_idx))/0.05
1554  1 *(rc_ndvi-(ndvi_idx*0.05-0.025))
1555  coeff1_red_fin = coeff1_red(ndvi_idx)
1556  1 +(coeff1_red(ndvi_idx+1)-coeff1_red(ndvi_idx))/0.05
1557  1 *(rc_ndvi-(ndvi_idx*0.05-0.025))
1558  coeff0_blue_fin = coeff0_blue(ndvi_idx)
1559  1 +(coeff0_blue(ndvi_idx+1)-coeff0_blue(ndvi_idx))/0.05
1560  1 *(rc_ndvi-(ndvi_idx*0.05-0.025))
1561  coeff1_blue_fin = coeff1_blue(ndvi_idx)
1562  1 +(coeff1_blue(ndvi_idx+1)-coeff1_blue(ndvi_idx))/0.05
1563  1 *(rc_ndvi-(ndvi_idx*0.05-0.025))
1564 
1565  xeas_b1 = coeff0_red_fin + coeff1_red_fin*sb2
1566  xeas_b3 = coeff0_blue_fin + coeff1_blue_fin*sb2
1567  else
1568  return
1569  endif
1570 
1571  elseif(season.eq.3) then ! JJA
1572 
1573  if(lcvr.gt.0.and.lcvr.ne.13) then
1574  coeff0_red(1:20) = (/ 0.000, 0.000, 0.000, 0.000, 0.000,
1575  1 0.000, 0.000, 0.000, 0.000, 0.000,
1576  1 0.000, 0.000, 0.000, 0.000, 0.000,
1577  1 0.000, 0.000, 0.000, 0.000, 0.000/)
1578  coeff1_red(1:20) = (/ 0.764, 0.764, 0.764, 0.703, 0.619,
1579  1 0.552, 0.489, 0.431, 0.380, 0.331,
1580  1 0.288, 0.247, 0.207, 0.171, 0.137,
1581  1 0.107, 0.080, 0.058, 0.041, 0.000/)
1582  coeff0_blue(1:20) =(/ 0.000, 0.000, 0.000, 0.000, 0.000,
1583  1 0.000, 0.000, 0.000, 0.000, 0.000,
1584  1 0.000, 0.000, 0.000, 0.000, 0.000,
1585  1 0.000, 0.000, 0.000, 0.000, 0.000/)
1586  coeff1_blue(1:20) =(/ 0.398, 0.398, 0.398, 0.358, 0.300,
1587  1 0.269, 0.245, 0.223, 0.202, 0.180,
1588  1 0.161, 0.144, 0.128, 0.113, 0.099,
1589  1 0.084, 0.069, 0.052, 0.037, 0.000/)
1590 
1591  ndvi_idx = min(floor((rc_ndvi-0.025)*20+1), 19)
1592  coeff0_red_fin = coeff0_red(ndvi_idx)
1593  1 +(coeff0_red(ndvi_idx+1)-coeff0_red(ndvi_idx))/0.05
1594  1 *(rc_ndvi-(ndvi_idx*0.05-0.025))
1595  coeff1_red_fin = coeff1_red(ndvi_idx)
1596  1 +(coeff1_red(ndvi_idx+1)-coeff1_red(ndvi_idx))/0.05
1597  1 *(rc_ndvi-(ndvi_idx*0.05-0.025))
1598  coeff0_blue_fin = coeff0_blue(ndvi_idx)
1599  1 +(coeff0_blue(ndvi_idx+1)-coeff0_blue(ndvi_idx))/0.05
1600  1 *(rc_ndvi-(ndvi_idx*0.05-0.025))
1601  coeff1_blue_fin = coeff1_blue(ndvi_idx)
1602  1 +(coeff1_blue(ndvi_idx+1)-coeff1_blue(ndvi_idx))/0.05
1603  1 *(rc_ndvi-(ndvi_idx*0.05-0.025))
1604 
1605  xeas_b1 = coeff0_red_fin + coeff1_red_fin*sb2
1606  xeas_b3 = coeff0_blue_fin + coeff1_blue_fin*sb2
1607  elseif(lcvr.eq.13) then
1608  coeff0_red(1:20) = (/ 0.000, 0.000, 0.000, 0.000, 0.000,
1609  1 0.000, 0.000, 0.000, 0.000, 0.000,
1610  1 0.000, 0.000, 0.000, 0.000, 0.000,
1611  1 0.000, 0.000, 0.000, 0.000, 0.000/)
1612  coeff1_red(1:20) = (/ 0.745, 0.745, 0.745, 0.681, 0.610,
1613  1 0.546, 0.485, 0.429, 0.373, 0.325,
1614  1 0.281, 0.241, 0.204, 0.170, 0.139,
1615  1 0.111, 0.087, 0.000, 0.000, 0.000/)
1616  coeff0_blue(1:20) =(/ 0.000, 0.000, 0.000, 0.000, 0.000,
1617  1 0.000, 0.000, 0.000, 0.000, 0.000,
1618  1 0.000, 0.000, 0.000, 0.000, 0.000,
1619  1 0.000, 0.000, 0.000, 0.000, 0.000/)
1620  coeff1_blue(1:20) =(/ 0.492, 0.492, 0.492, 0.471, 0.422,
1621  1 0.376, 0.338, 0.300, 0.260, 0.226,
1622  1 0.197, 0.183, 0.169, 0.144, 0.121,
1623  1 0.100, 0.083, 0.000, 0.000, 0.000/)
1624 
1625  ndvi_idx = min(floor((rc_ndvi-0.025)*20+1), 17)
1626  coeff0_red_fin = coeff0_red(ndvi_idx)
1627  1 +(coeff0_red(ndvi_idx+1)-coeff0_red(ndvi_idx))/0.05
1628  1 *(rc_ndvi-(ndvi_idx*0.05-0.025))
1629  coeff1_red_fin = coeff1_red(ndvi_idx)
1630  1 +(coeff1_red(ndvi_idx+1)-coeff1_red(ndvi_idx))/0.05
1631  1 *(rc_ndvi-(ndvi_idx*0.05-0.025))
1632  coeff0_blue_fin = coeff0_blue(ndvi_idx)
1633  1 +(coeff0_blue(ndvi_idx+1)-coeff0_blue(ndvi_idx))/0.05
1634  1 *(rc_ndvi-(ndvi_idx*0.05-0.025))
1635  coeff1_blue_fin = coeff1_blue(ndvi_idx)
1636  1 +(coeff1_blue(ndvi_idx+1)-coeff1_blue(ndvi_idx))/0.05
1637  1 *(rc_ndvi-(ndvi_idx*0.05-0.025))
1638 
1639  xeas_b1 = coeff0_red_fin + coeff1_red_fin*sb2
1640  xeas_b3 = coeff0_blue_fin + coeff1_blue_fin*sb2
1641  else
1642  return
1643  endif
1644 
1645  elseif(season.eq.4) then ! SON
1646 
1647  if(lcvr.gt.0.and.lcvr.ne.13) then
1648  coeff0_red(1:20) = (/ 0.000, 0.000, 0.000, 0.000, 0.000,
1649  1 0.000, 0.000, 0.000, 0.000, 0.000,
1650  1 0.000, 0.000, 0.000, 0.000, 0.000,
1651  1 0.000, 0.000, 0.000, 0.000, 0.000/)
1652  coeff1_red(1:20) = (/ 0.762, 0.762, 0.762, 0.691, 0.617,
1653  1 0.549, 0.488, 0.432, 0.380, 0.332,
1654  1 0.288, 0.248, 0.210, 0.175, 0.141,
1655  1 0.110, 0.083, 0.063, 0.043, 0.000/)
1656  coeff0_blue(1:20) =(/ 0.000, 0.000, 0.000, 0.000, 0.000,
1657  1 0.000, 0.000, 0.000, 0.000, 0.000,
1658  1 0.000, 0.000, 0.000, 0.000, 0.000,
1659  1 0.000, 0.000, 0.000, 0.000, 0.000/)
1660  coeff1_blue(1:20) =(/ 0.395, 0.395, 0.395, 0.342, 0.309,
1661  1 0.276, 0.247, 0.223, 0.204, 0.185,
1662  1 0.166, 0.150, 0.136, 0.119, 0.101,
1663  1 0.085, 0.070, 0.054, 0.028, 0.000/)
1664 
1665  ndvi_idx = min(floor((rc_ndvi-0.025)*20+1), 19)
1666  coeff0_red_fin = coeff0_red(ndvi_idx)
1667  1 +(coeff0_red(ndvi_idx+1)-coeff0_red(ndvi_idx))/0.05
1668  1 *(rc_ndvi-(ndvi_idx*0.05-0.025))
1669  coeff1_red_fin = coeff1_red(ndvi_idx)
1670  1 +(coeff1_red(ndvi_idx+1)-coeff1_red(ndvi_idx))/0.05
1671  1 *(rc_ndvi-(ndvi_idx*0.05-0.025))
1672  coeff0_blue_fin = coeff0_blue(ndvi_idx)
1673  1 +(coeff0_blue(ndvi_idx+1)-coeff0_blue(ndvi_idx))/0.05
1674  1 *(rc_ndvi-(ndvi_idx*0.05-0.025))
1675  coeff1_blue_fin = coeff1_blue(ndvi_idx)
1676  1 +(coeff1_blue(ndvi_idx+1)-coeff1_blue(ndvi_idx))/0.05
1677  1 *(rc_ndvi-(ndvi_idx*0.05-0.025))
1678 
1679  xeas_b1 = coeff0_red_fin + coeff1_red_fin*sb2
1680  xeas_b3 = coeff0_blue_fin + coeff1_blue_fin*sb2
1681  elseif(lcvr.eq.13) then
1682  coeff0_red(1:20) = (/ 0.000, 0.000, 0.000, 0.000, 0.000,
1683  1 0.000, 0.000, 0.000, 0.000, 0.000,
1684  1 0.000, 0.000, 0.000, 0.000, 0.000,
1685  1 0.000, 0.000, 0.000, 0.000, 0.000/)
1686  coeff1_red(1:20) = (/ 0.735, 0.735, 0.735, 0.672, 0.601,
1687  1 0.532, 0.473, 0.419, 0.369, 0.324,
1688  1 0.283, 0.243, 0.204, 0.169, 0.141,
1689  1 0.119, 0.093, 0.000, 0.000, 0.000/)
1690  coeff0_blue(1:20) =(/ 0.000, 0.000, 0.000, 0.000, 0.000,
1691  1 0.000, 0.000, 0.000, 0.000, 0.000,
1692  1 0.000, 0.000, 0.000, 0.000, 0.000,
1693  1 0.000, 0.000, 0.000, 0.000, 0.000/)
1694  coeff1_blue(1:20) =(/ 0.492, 0.492, 0.492, 0.482, 0.429,
1695  1 0.371, 0.330, 0.293, 0.258, 0.232,
1696  1 0.209, 0.185, 0.159, 0.136, 0.119,
1697  1 0.103, 0.083, 0.000, 0.000, 0.000/)
1698 
1699  ndvi_idx = min(floor((rc_ndvi-0.025)*20+1), 17)
1700  coeff0_red_fin = coeff0_red(ndvi_idx)
1701  1 +(coeff0_red(ndvi_idx+1)-coeff0_red(ndvi_idx))/0.05
1702  1 *(rc_ndvi-(ndvi_idx*0.05-0.025))
1703  coeff1_red_fin = coeff1_red(ndvi_idx)
1704  1 +(coeff1_red(ndvi_idx+1)-coeff1_red(ndvi_idx))/0.05
1705  1 *(rc_ndvi-(ndvi_idx*0.05-0.025))
1706  coeff0_blue_fin = coeff0_blue(ndvi_idx)
1707  1 +(coeff0_blue(ndvi_idx+1)-coeff0_blue(ndvi_idx))/0.05
1708  1 *(rc_ndvi-(ndvi_idx*0.05-0.025))
1709  coeff1_blue_fin = coeff1_blue(ndvi_idx)
1710  1 +(coeff1_blue(ndvi_idx+1)-coeff1_blue(ndvi_idx))/0.05
1711  1 *(rc_ndvi-(ndvi_idx*0.05-0.025))
1712 
1713  xeas_b1 = coeff0_red_fin + coeff1_red_fin*sb2
1714  xeas_b3 = coeff0_blue_fin + coeff1_blue_fin*sb2
1715  else
1716  return
1717  endif
1718 
1719  else
1720  return
1721  endif
1722 
1723  return
1724  end
1725 
1726 c... --------------------------------------
1727 c... The end of subroutine get_sfcrfl_nir_vis
1728 c... ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1729 
1730  subroutine aero_470veg(refl,x1,x2,x3,mm,nn,ll,ma,
1731  1 imod,r470,tau_x470,tau_x470_flag)
1732 
1733  use viirs_aerosol_luts, only: new_intep
1734 
1735  include 'aottbl.inc'
1736  include 'newaottbl.inc'
1737 
1738 c common /angle_node/ theta0, theta, phi
1739 c common /sfcref_node/ sfc_ref412, sfc_ref470, sfc_ref650
1740 c real theta0(10), theta(46), phi(30), tau(10)
1741 c real sfc_ref412(20), sfc_ref470(24), sfc_ref650(24)
1742  real nnvalx(4,4,2,10), yy(10), yy2(8)
1743 c real nvalx470(10,46,30,10,4,24)
1744  integer tau_x470_flag, imod
1745  data pi /3.14159/
1746 c data tau /0.0, 0.1, 0.3, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5/
1747 
1748 c print *, 'aero_470veg in: ', refl,x1,x2,x3,mm,nn,ll,ma,imod,r470,tau_x470,tau_x470_flag
1749  index_ii = r470
1750  if (index_ii.lt.0) return ! orig.
1751 
1752  frac = (r470-sfc_ref470(index_ii))/
1753  1 (sfc_ref470(index_ii+1)-sfc_ref470(index_ii))
1754 
1755  if (index_ii.lt.1.or.index_ii.gt.24)
1756  1 print *,'index_iir470 = ', index_ii,xlat,xlong
1757  if (frac.lt.0.0.or.frac.gt.1.0)
1758  1 print *,'frac on sfc470=', frac
1759  if (index_ii.lt.1.or.index_ii.gt.24) then ! MJ added
1760  tau_x470_flag = 9
1761  tau_x470=-999.0
1762  return
1763  endif
1764 
1765 
1766  call search(dflag2,x3,phi,ll,ii)
1767  xfrac = (x3-phi(ii))/(phi(ii+1)-phi(ii))
1768 
1769  nsm=1
1770  dif=x1-theta0(1)
1771  do i=1,mm
1772  dift=x1-theta0(i)
1773  if (dift.gt.0. .and. dift.lt.dif) then
1774  nsm=i
1775  dif=dift
1776  endif
1777  enddo
1778  mbeg = nsm - 2
1779  if (mbeg.le.0) then
1780  mbeg = 0
1781  else if (mbeg.gt.mm-4) then
1782  mbeg = mm-4
1783  endif
1784 
1785  nsn=1
1786  dif=x2-theta(1)
1787  do i=1,nn
1788  dift=x2-theta(i)
1789  if (dift.gt.0. .and. dift.lt.dif) then
1790  nsn=i
1791  dif=dift
1792  endif
1793  enddo
1794  nbeg = nsn - 2
1795  if (nbeg.le.0) then
1796  nbeg = 0
1797  else if (nbeg.gt.nn-4) then
1798  nbeg = nn-4
1799  endif
1800 
1801 c write(6,*) "frac =",frac
1802  do ia = 1, 10
1803  do i = 1, 4
1804  do j = 1, 4
1805  nnvalx(i,j,1,ia) = nvalx470(mbeg+i,nbeg+j,ii,ia,imod,index_ii)*
1806  1 (1.-frac) + nvalx470(mbeg+i,nbeg+j,ii,ia,imod,index_ii+1)*frac
1807  nnvalx(i,j,2,ia) = nvalx470(mbeg+i,nbeg+j,ii+1,ia,imod,index_ii)*
1808  1 (1.-frac) + nvalx470(mbeg+i,nbeg+j,ii+1,ia,imod,index_ii+1)*frac
1809  enddo
1810  enddo
1811  enddo
1812 
1813 c--- interpolating AOT tables
1814 
1815  do 105 ia = 1, 10
1816 
1817  call new_intep(theta0, theta, phi, nnvalx, mm, nn, ll, ia,
1818  1 x1,x2,x3,y,dy,mbeg,nbeg,xfrac)
1819 
1820  yy(ia) = y/pi
1821 c print *,'tau, i/f=', tau(ia), y/pi, dy
1822  105 continue
1823  if (refl.le.yy(1)) then
1824  tau_x470 = 0.02
1825  tau_x470_flag = 2 ! MJ added
1826  return
1827  endif
1828 
1829  if (refl.ge.yy(10)) then
1830  xxrat = 0.8
1831  tau_x470 = 3.5
1832  tau_x470_flag = 1
1833  return
1834  endif
1835 c
1836 c
1837 c Check if the reflectance increase with AOT
1838 c
1839 
1840  if (yy(1).lt.yy(2)) go to 650
1841 
1842  if (refl.lt.yy(4)) return
1843 
1844  do i = 1, 7
1845  yy2(i) = yy(i+3)
1846  enddo
1847 
1848  if (yy2(2).lt.yy2(1)) return
1849  call search2(dflag2,refl,yy2,7,index_ii,frac)
1850 c print *,'after 1st search'
1851 c if (yy2(1).gt.yy2(2)) print *,'yy=',refl,(yy(i),I=1,10)
1852  tau_x = frac*tau(index_ii+1+3) + (1.-frac)*tau(index_ii+3)
1853  tau_x470 = tau_x
1854  return
1855 
1856 650 continue
1857 c
1858 c Pass the monotonic order check
1859 c
1860  call search2(dflag2,refl,yy,10,index_ii,frac)
1861 c print *,'after 2nd search'
1862  tau_x = frac*tau(index_ii+1) + (1.-frac)*tau(index_ii)
1863 
1864  tau_x470 = tau_x
1865 
1866 c print *,'tau_x470 =', tau_x470
1867 
1868  return
1869  end
1870 c... --------------------------------------
1871 c... The end of subroutine aero_470veg
1872 c... ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1873 
1874 
1875 
1876 c... ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1877 c--------------------------------------------------------
1878  subroutine aero_650veg(refl,x1,x2,x3,mm,nn,ll,ma,
1879  1 r650,tau_x650,tau_x650_flag,
1880  2 tau_x470_flag,tau_x470)
1881 
1882  use viirs_aerosol_luts, only: new_intep
1883 
1884  include 'aottbl.inc'
1885  include 'newaottbl.inc'
1886 
1887 c common /angle_node/ theta0, theta, phi
1888 c common /sfcref_node/ sfc_ref412, sfc_ref470, sfc_ref650
1889 
1890 c real theta0(10), theta(46), phi(30), tau(10)
1891 c real sfc_ref412(20), sfc_ref470(24), sfc_ref650(24)
1892  real nnvalx(4,4,2,10), yy(10), yy2(8), yy3(3), yy5(6)
1893  real tau_x650, tau_x412, tau_x470
1894  real refl,x1,x2,x3,r650
1895 c real nvalx650(10,46,30,10,24)
1896  integer tau_x470_flag, tau_x650_flag, tau_x412_flag_91
1897  logical dflag2
1898  data pi /3.14159/
1899 c data tau /0.0, 0.1, 0.3, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5/
1900 
1901  index_ii = (r650+1.)/2.
1902 
1903  frac = (r650-sfc_ref650(index_ii))/
1904  1 (sfc_ref650(index_ii+1)-sfc_ref650(index_ii))
1905 
1906 c if (index_ii.lt.1.or.index_ii.gt.24)
1907 c 1 print *,'index_iir650 = ', index_ii,xlat,xlong
1908 c if (frac.lt.0.0.or.frac.gt.1.0)
1909 c 1 print *,'frac on sfc=', frac
1910 
1911 c... --------------------------------------------------
1912 c... MJ added 12Feb2010 @@@@@
1913 c... --------------------------------------------------
1914 c. if(index_ii.lt.2.or.index_ii.gt.12) then ! ok ver.
1915  if(index_ii.lt.1.or.index_ii.gt.23) then
1916  tau_x650_flag = 9
1917  tau_x650 = -999.0
1918  return
1919  endif ! if(index_ii.lt.1.or.index_ii.gt.23) then
1920 c... --------------------------------------------------
1921 
1922  if (index_ii.lt.1) then
1923  index_ii = 1
1924  frac = 0.
1925  endif
1926 
1927  call search(dflag2,x3,phi,ll,ii)
1928  xfrac = (x3-phi(ii))/(phi(ii+1)-phi(ii))
1929 
1930  nsm=1
1931  dif=x1-theta0(1)
1932  do i=1,mm
1933  dift=x1-theta0(i)
1934  if (dift.gt.0. .and. dift.lt.dif) then
1935  nsm=i
1936  dif=dift
1937  endif
1938  enddo
1939  mbeg = nsm - 2
1940  if (mbeg.le.0) then
1941  mbeg = 0
1942  else if (mbeg.gt.mm-4) then
1943  mbeg = mm-4
1944  endif
1945 
1946  nsn=1
1947  dif=x2-theta(1)
1948  do i=1,nn
1949  dift=x2-theta(i)
1950  if (dift.gt.0. .and. dift.lt.dif) then
1951  nsn=i
1952  dif=dift
1953  endif
1954  enddo
1955  nbeg = nsn - 2
1956  if (nbeg.le.0) then
1957  nbeg = 0
1958  else if (nbeg.gt.nn-4) then
1959  nbeg = nn-4
1960  endif
1961 
1962 c write(6,*) "frac =",frac
1963  do ia = 1, 10
1964  do i = 1, 4
1965  do j = 1, 4
1966  nnvalx(i,j,1,ia) = nvalx650(mbeg+i,nbeg+j,ii,ia,index_ii)*
1967  1 (1.-frac) + nvalx650(mbeg+i,nbeg+j,ii,ia,index_ii+1)*frac
1968  nnvalx(i,j,2,ia) = nvalx650(mbeg+i,nbeg+j,ii+1,ia,index_ii)*
1969  1 (1.-frac) + nvalx650(mbeg+i,nbeg+j,ii+1,ia,index_ii+1)*frac
1970  enddo
1971  enddo
1972  enddo
1973 
1974 c--- interpolating AOT tables
1975 
1976  do 600 ia = 1, 10
1977 
1978  call new_intep(theta0, theta, phi, nnvalx, mm, nn, ll, ia,
1979  1 x1,x2,x3,y,dy,mbeg,nbeg,xfrac)
1980 
1981  yy(ia) = y/pi
1982 c print *,'tau, i/f=', tau(ia), y/pi, dy
1983  600 continue
1984 
1985 
1986  if (refl.le.yy(1).and.yy(1).lt.yy(2)) then
1987  tau_x650_flag = 2 ! MJ added
1988  tau_x650 = 0.02
1989  return
1990  endif
1991 
1992  if (refl.ge.yy(10)) then
1993  tau_x650 = 3.5
1994  w0_x = -999.
1995  tau_x650_flag = 1
1996  return
1997  endif
1998 
1999 c!!! if (tau_x470_flag.gt.0) go to 670 ! ori open 20100216 @@@@@
2000 c if (tau_x412_flag_91.gt.0) go to 680
2001 c
2002 c
2003 c Check if the reflectance increase with AOT
2004 c
2005  if (yy(1).lt.yy(2)) go to 650
2006 
2007 c... *****************************
2008 c... *****************************
2009 c... *****************************
2010  tau_x650_flag = 4 ! mj added 08/05/2011
2011  return ! mj added 20100216 @@@@@
2012  ! no retrieval for non-monotonic
2013  ! changes in LUT reflectance
2014 c... needs investigation to comment this out!
2015 c... *****************************
2016 c... *****************************
2017 c... *****************************
2018 
2019  if (refl.lt.yy(4)) return
2020 
2021  do i = 1, 7
2022  yy2(i) = yy(i+3)
2023  enddo
2024 
2025  if (yy2(2).lt.yy2(1)) return
2026  call search2(dflag2,refl,yy2,7,index_ii,frac)
2027 c print *,'after 1st search'
2028 c if (yy2(1).gt.yy2(2)) print *,'yy=',refl,(yy(i),I=1,10)
2029  tau_x = frac*tau(index_ii+1+3) + (1.-frac)*tau(index_ii+3)
2030  tau_x650 = tau_x
2031  return
2032 
2033 670 continue
2034  if (refl.lt.yy(8)) return
2035 
2036  do i = 1, 3
2037  yy3(i) = yy(i+7)
2038  enddo
2039 
2040  if (yy3(2).lt.yy3(1)) return
2041  call search2(dflag2,refl,yy3,3,index_ii,frac)
2042 c print *,'after 1st search'
2043 c if (yy3(1).gt.yy3(2)) print *,'yy=',refl,(yy(i),I=1,10)
2044  tau_x = frac*tau(index_ii+1+7) + (1.-frac)*tau(index_ii+7)
2045  tau_x650 = tau_x
2046  return
2047 
2048 680 continue
2049  if (refl.lt.yy(5)) return
2050 
2051  do i = 1, 6
2052  yy5(i) = yy(i+4)
2053  enddo
2054 
2055  if (yy5(2).lt.yy5(1)) return
2056  call search2(dflag2,refl,yy5,6,index_ii,frac)
2057 c print *,'after 1st search'
2058 c if (yy3(1).gt.yy3(2)) print *,'yy=',refl,(yy(i),I=1,10)
2059  tau_x = frac*tau(index_ii+1+4) + (1.-frac)*tau(index_ii+4)
2060  tau_x650 = tau_x
2061  return
2062 
2063 650 continue
2064 c
2065 c Pass the monotonic order check
2066 c
2067  call search2(dflag2,refl,yy,10,index_ii,frac)
2068 c print *,'after 2nd search'
2069  tau_x = frac*tau(index_ii+1) + (1.-frac)*tau(index_ii)
2070 
2071  tau_x650 = tau_x
2072 
2073  return
2074  end
2075 c... --------------------------------------
2076 c... The end of subroutine aero_650veg
2077 c... ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2078 
2079 c **********************************************************************************************************
2080 
2081 
2082 c **********************************************************************************************************
2083  subroutine cmb_stddb_veg(outbuf, outbufvg, tmpvg, idalg)
2084 c...
2085 c... A subroutine to combine standard Deep Blue and over-vegetation
2086 c... retrievals, depending on ndvi, ndvi_swir, refl at 2.1um, etc.
2087 c... Written by Myeong-Jae Jeong (MJ)
2088 c... Ver. 0.5 Aug 8, 2011
2089 c...
2090 c.... Inputs:
2091 c... outbuf(20) - standard DeepBlue output
2092 c... outbufvg(20) - over-vegetation retr. output
2093 c... tmpvg(6) - TOA reflectances
2094 c...
2095 c... Output:
2096 c... outfub(20) - combined output
2097 c... idalg - index for algorithm in effect
2098 c...
2099 c... idalg=0 : standard Deep Blue
2100 c... idalg=1 : over-vegetation retrieval
2101 c... idalg=2 : mixture of standard DB and over-veg. retrievals
2102 
2103  real outbuf(20), outbufvg(20), tmpvg(6)
2104  real outbuftmp(20), xfrdb
2105  real c_ndsir1,c_ndvi1,c_rf21_1,c_rf21_2,c_rf21_3
2106  integer idalg
2107 
2108  c_ndsir1=0.1 ! NDVI_swir minimum threshold
2109  c_ndvi1=0.1 ! NDVI minimum threshold
2110  c_rf21_1=0.01 ! 2.1um refl. threshold 1
2111  c_rf21_2=0.25 ! 2.1um refl. threshold 2
2112  c_rf21_3=0.35 ! 2.1um refl. threshold 3
2113 
2114 c... Note. Also try to employ NDVI_swir thresholds later to combine
2115 c... standard DeepBlue and over-vegeation retrievals... (MJ)
2116 
2117 c... Initialization ... set to standard DeepBlue retrieval results
2118  do i=1,20
2119  outbuftmp(i)=outbuf(i)
2120  enddo
2121  idalg=0 ! standard DeepBlue
2122 
2123 
2124 c if(outbufvg(18).ge.c_ndsir1.and.outbufvg(19).ge.c_ndvi1) then
2125 c if(tmpvg(4).ge.c_rf21_1.and.tmpvg(4).lt.c_rf21_2) then ! veg.
2126 c do i=1,8
2127 c outbuftmp(i)=outbufvg(i)
2128 c enddo
2129 c do i=11,14
2130 c outbuftmp(i)=outbufvg(i)
2131 c enddo
2132 c outbuftmp(15) = 100.0
2133 c idalg=1
2134 c elseif(tmpvg(4).ge.c_rf21_2.and.tmpvg(4).lt.c_rf21_3) then
2135 cc xfrdb=(tmpvg(4)-c_rf21_2)/(c_rf21_3-c_rf21_2)
2136 cc if(outbufvg(7).ge.0.and.outbuf(7).ge.0) then ! veg+DB
2137 cc do i=2,3
2138 cc outbuftmp(i)=(1.0-xfrdb)*outbufvg(i)+xfrdb*outbuf(i)
2139 cc outbuftmp(i+3)=((1.0-xfrdb)*outbufvg(i+3)*outbufvg(i)
2140 cc 1 + xfrdb*outbuf(i+3)*outbuf(i))
2141 cc 2 /(outbuftmp(i))
2142 cc enddo
2143 cc outbuftmp(7)=(1.0-xfrdb)*outbufvg(7)+xfrdb*outbuf(7)
2144 cc outbuftmp(1)=-999.0 ! set aot412nm to -999.0
2145 cc outbuftmp(4)=-999.0 ! set ssa412nm to -999.0
2146 cc outbuftmp(8)=alog(outbuftmp(2)/outbuftmp(3))/
2147 cc 1 alog(650./466.)
2148 cc if(outbufvg(11).gt.0.and.outbuf(11).gt.0) then
2149 cc outbuftmp(11)=((1.0-xfrdb)*outbufvg(11)*outbufvg(2)
2150 cc 1 +xfrdb*outbuf(11)*outbuf(2))
2151 cc 2 /(outbuftmp(2))
2152 cc outbuftmp(12)=((1.0-xfrdb)*outbufvg(12)*outbufvg(3)
2153 cc 1 +xfrdb*outbuf(12)*outbuf(3))
2154 cc 2 /(outbuftmp(3))
2155 cc else
2156 cc outbuftmp(11)=-999.0
2157 cc outbuftmp(12)=-999.0
2158 cc endif
2159 cc idalg=2
2160 c if(outbufvg(7).ge.0.and.outbuf(7).lt.0) then
2161 c do i=1,8
2162 c outbuftmp(i)=outbufvg(i)
2163 c enddo
2164 c do i=11,14
2165 c outbuftmp(i)=outbufvg(i)
2166 c enddo
2167 c outbuftmp(15) = 100.0
2168 c idalg=1
2169 c else
2170 c do i=1,8
2171 c outbuftmp(i)=outbuf(i)
2172 c enddo
2173 c do i=11,14
2174 c outbuftmp(i)=outbuf(i)
2175 c enddo
2176 c idalg=0
2177 c endif
2178 c endif
2179 c endif
2180 
2181 c... Finalization ... push the results to outbuf
2182  do i=1,20
2183  outbuf(i)=outbuftmp(i)
2184  enddo
2185 
2186  return
2187  end
2188 c... --------------------------------------
2189 c... The end of subroutine cmb_stddb_veg
2190 c... ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2191 c **********************************************************************************************************
2192 
2193 
2194 c **********************************************************************************************************
2195  subroutine calc_sfc21(sza,xthet,xphi,refl21,aot21,As_21)
2196 c...
2197 c... ----------------------------------------------------------
2198 c... written by Myeong-Jae Jeong (MJ) at GWNU, Korea
2199 c... ver. 1.0 Aug 12, 2011
2200 c...
2201 c... subroutine to calculate surface reflectance at 2.1um.
2202 c... code inherited from "sfc470_5wav.f" by N.C. Hsu
2203 c...
2204 c... Inputs
2205 c... sza : solar zenith angle
2206 c... xthet : viewing zenith angle
2207 c... xphi : relative azimuth angle
2208 c... refl21: 2.1um TOA reflectance
2209 c... aot21 : 2.1um AOT
2210 c...
2211 c... Output
2212 c... As_21 : 2.1 surface reflectance (%)
2213 c...
2214 c... ----------------------------------------------------------
2215 c... refl21 = refl21 * xmu/3.14159
2216 
2217  include 'aottbl.inc'
2218  include 'sfc21tbl.inc'
2219 
2220 c real tau2(4), theta0(10), theta(46), phi(30)
2221  real tau2(4)
2222  real nnvalx(4,4,2), rr0x(4,4,2)
2223  real ttx(4,4,2), ssx(4,4,2)
2224  logical dflag3
2225  integer index_ii
2226 
2227  data tau2 /0.0, 0.1, 0.3, 0.5/
2228  data pi /3.14159/
2229 c data theta0 /0.0,8.0,16.0,24.0,32.0,40.0,48.0,56.0,64.0,72.0/
2230 
2231  x1 = sza
2232  x2 = xthet
2233  x3 = xphi
2234  refl7 = refl21 ! 2.1 micron
2235 
2236 c do i = 1, 46
2237 c theta(i) = 2.*float(i-1)
2238 c enddo
2239 
2240 c do i = 1, 30
2241 c phi(i) = 6. + 6.*float(i-1)
2242 c enddo
2243 
2244  mm = 10 ! solar zenith
2245  nn = 46 ! satellite zenith
2246  ll = 30 ! rel. azimuth
2247  ma2= 4 ! tau2 (aot at 2.1um)
2248 
2249 c -- Derive SFC
2250 
2251 c-- interpolating for 2.1 micron channel
2252 
2253  nsm=1
2254  dif=x1-theta0(1)
2255  do i=1,mm
2256  dift=x1-theta0(i)
2257  if (dift.gt.0. .and. dift.lt.dif) then
2258  nsm=i
2259  dif=dift
2260  endif
2261  enddo
2262  mbeg = nsm - 2
2263  if (mbeg.le.0) then
2264  mbeg = 0
2265  else if (mbeg.gt.mm-4) then
2266  mbeg = mm-4
2267  endif
2268 
2269  nsn=1
2270  dif=x2-theta(1)
2271  do i=1,nn
2272  dift=x2-theta(i)
2273  if (dift.gt.0. .and. dift.lt.dif) then
2274  nsn=i
2275  dif=dift
2276  endif
2277  enddo
2278  nbeg = nsn - 2
2279  if (nbeg.le.0) then
2280  nbeg = 0
2281  else if (nbeg.gt.nn-4) then
2282  nbeg = nn-4
2283  endif
2284 c...
2285  aot2=aot21
2286  call search(dflag3,aot2,tau2,ma2,index_ii)
2287  frac = (aot2-tau2(index_ii))
2288  1 /(tau2(index_ii+1)-tau2(index_ii))
2289  if(frac.lt.0.or.frac.gt.1)
2290  1 print *,'aot2, frac, index_ii=',aot2, frac, index_ii
2291 
2292  call search(dflag3,x3,phi,ll,ii)
2293  xfrac = (x3-phi(ii))/(phi(ii+1)-phi(ii))
2294 
2295 c print *, 'nsm, nsn, ii, aot2, frac, index_ii'
2296 c print *, nsm, nsn, ii, aot2, frac, index_ii
2297  do 123 i = 1, 4
2298  do 124 j = 1, 4
2299  nnvalx(i,j,1) = nvalx21(mbeg+i,nbeg+j,ii,index_ii)*
2300  1 (1.-frac) + nvalx21(mbeg+i,nbeg+j,ii,index_ii+1)*frac
2301  nnvalx(i,j,2) = nvalx21(mbeg+i,nbeg+j,ii+1,index_ii)*
2302  1 (1.-frac) + nvalx21(mbeg+i,nbeg+j,ii+1,index_ii+1)*frac
2303 
2304  rr0x(i,j,1) = r0x_21(mbeg+i,nbeg+j,ii,index_ii)*
2305  1 (1.-frac) + r0x_21(mbeg+i,nbeg+j,ii,index_ii+1)*frac
2306  rr0x(i,j,2) = r0x_21(mbeg+i,nbeg+j,ii+1,index_ii)*
2307  1 (1.-frac) + r0x_21(mbeg+i,nbeg+j,ii+1,index_ii+1)*frac
2308 
2309  ttx(i,j,1) = tx_21(mbeg+i,nbeg+j,ii,index_ii)*
2310  1 (1.-frac) + tx_21(mbeg+i,nbeg+j,ii,index_ii+1)*frac
2311  ttx(i,j,2) = tx_21(mbeg+i,nbeg+j,ii+1,index_ii)*
2312  1 (1.-frac) + tx_21(mbeg+i,nbeg+j,ii+1,index_ii+1)*frac
2313 
2314  ssx(i,j,1) = sx_21(mbeg+i,nbeg+j,ii,index_ii)*
2315  1 (1.-frac) + sx_21(mbeg+i,nbeg+j,ii,index_ii+1)*frac
2316  ssx(i,j,2) = sx_21(mbeg+i,nbeg+j,ii+1,index_ii)*
2317  1 (1.-frac) + sx_21(mbeg+i,nbeg+j,ii+1,index_ii+1)*frac
2318 
2319 124 continue
2320 123 continue
2321 c enddo
2322 c enddo
2323 
2324 c---
2325 
2326  call new_intepsf(theta0, theta, phi, nnvalx, mm, nn, ll,
2327  1 x1,x2,x3,y,dy,mbeg,nbeg,xfrac)
2328 
2329  rr = y/pi
2330 
2331  call new_intepsf(theta0, theta, phi, rr0x, mm, nn, ll,
2332  1 x1,x2,x3,y,dy,mbeg,nbeg,xfrac)
2333 
2334  rr0 = y/pi
2335 
2336  call new_intepsf(theta0, theta, phi, ttx, mm, nn, ll,
2337  1 x1,x2,x3,y,dy,mbeg,nbeg,xfrac)
2338 
2339  tt = y/pi
2340 
2341  call new_intepsf(theta0, theta, phi, ssx, mm, nn, ll,
2342  1 x1,x2,x3,y,dy,mbeg,nbeg,xfrac)
2343 
2344 c! SS = y/pi
2345  ss = y ! chg20090818
2346 
2347  dff = refl7 - rr0
2348  as_21 = 100.* dff / (tt + ss*dff)
2349 
2350 c print *,nday,sza,xthet,xphi,refl3,refl7,As,As_21
2351 c print *,'refl7,RR0,dff,TT,SS,As_21= ',
2352 c 1 refl7,RR0,dff,TT,SS,As_21
2353 
2354  return
2355  end
2356 c... --------------------------------------
2357 c... The end of subroutine calc_sfc21
2358 c... ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2359 c **********************************************************************************************************
2360 
2361 c **********************************************************************************************************
2362  subroutine calc_rc_ref672(sza,xthet,xphi,ref_toa,ref_rc)
2363 c...
2364 c... ----------------------------------------------------------
2365 c... written by Jaehwa Lee
2366 c... ver. 1.0 Mar. 21, 2014
2367 c...
2368 c... subroutine to calculate Rayleigh-corrected reflectance
2369 c... code inherited from "calc_sfc21.f" by MJ
2370 c...
2371 c... Inputs
2372 c... sza : solar zenith angle
2373 c... xthet : viewing zenith angle
2374 c... xphi : relative azimuth angle
2375 c... ref_toa: TOA reflectance (I/F)
2376 c...
2377 c... Output
2378 c... ref_rc : Rayleigh-corrected reflectance (normalized)
2379 c...
2380 c... ----------------------------------------------------------
2381 
2382  include 'aottbl.inc'
2383  include 'sfc21tbl.inc'
2384 
2385 c real tau2(4), theta0(10), theta(46), phi(30)
2386  real tau2(4)
2387  real nnvalx(4,4,2), rr0x(4,4,2)
2388  real ttx(4,4,2), ssx(4,4,2)
2389  logical dflag3
2390 
2391  data pi /3.14159/
2392 
2393  x1 = sza
2394  x2 = xthet
2395  x3 = xphi
2396 
2397  mm = 10 ! solar zenith
2398  nn = 46 ! satellite zenith
2399  ll = 30 ! rel. azimuth
2400  ma2= 4 ! tau2 (aot at 2.1um)
2401 
2402 c-- interpolating for 2.1 micron channel
2403 
2404  nsm=1
2405  dif=x1-theta0(1)
2406  do i=1,mm
2407  dift=x1-theta0(i)
2408  if (dift.gt.0. .and. dift.lt.dif) then
2409  nsm=i
2410  dif=dift
2411  endif
2412  enddo
2413  mbeg = nsm - 2
2414  if (mbeg.le.0) then
2415  mbeg = 0
2416  else if (mbeg.gt.mm-4) then
2417  mbeg = mm-4
2418  endif
2419 
2420  nsn=1
2421  dif=x2-theta(1)
2422  do i=1,nn
2423  dift=x2-theta(i)
2424  if (dift.gt.0. .and. dift.lt.dif) then
2425  nsn=i
2426  dif=dift
2427  endif
2428  enddo
2429  nbeg = nsn - 2
2430  if (nbeg.le.0) then
2431  nbeg = 0
2432  else if (nbeg.gt.nn-4) then
2433  nbeg = nn-4
2434  endif
2435 
2436  call search(dflag3,x3,phi,ll,ii)
2437  xfrac = (x3-phi(ii))/(phi(ii+1)-phi(ii))
2438 
2439  do 123 i = 1, 4
2440  do 124 j = 1, 4
2441  nnvalx(i,j,1) = nvalx672(mbeg+i,nbeg+j,ii,1)
2442  nnvalx(i,j,2) = nvalx672(mbeg+i,nbeg+j,ii+1,1)
2443 
2444  rr0x(i,j,1) = r0x_672(mbeg+i,nbeg+j,ii,1)
2445  rr0x(i,j,2) = r0x_672(mbeg+i,nbeg+j,ii+1,1)
2446 
2447  ttx(i,j,1) = tx_672(mbeg+i,nbeg+j,ii,1)
2448  ttx(i,j,2) = tx_672(mbeg+i,nbeg+j,ii+1,1)
2449 
2450  ssx(i,j,1) = sx_672(mbeg+i,nbeg+j,ii,1)
2451  ssx(i,j,2) = sx_672(mbeg+i,nbeg+j,ii+1,1)
2452 
2453 124 continue
2454 123 continue
2455 
2456  call new_intepsf(theta0, theta, phi, nnvalx, mm, nn, ll,
2457  1 x1,x2,x3,y,dy,mbeg,nbeg,xfrac)
2458 
2459  rr = y/pi
2460 
2461  call new_intepsf(theta0, theta, phi, rr0x, mm, nn, ll,
2462  1 x1,x2,x3,y,dy,mbeg,nbeg,xfrac)
2463 
2464  rr0 = y/pi
2465 
2466  call new_intepsf(theta0, theta, phi, ttx, mm, nn, ll,
2467  1 x1,x2,x3,y,dy,mbeg,nbeg,xfrac)
2468 
2469  tt = y/pi
2470 
2471  call new_intepsf(theta0, theta, phi, ssx, mm, nn, ll,
2472  1 x1,x2,x3,y,dy,mbeg,nbeg,xfrac)
2473 
2474  ss = y
2475 
2476  dff = ref_toa - rr0
2477  ref_rc = dff / (tt + ss*dff)
2478 
2479  return
2480  end
2481 c... --------------------------------------
2482 c... The end of subroutine calc_rc_ref672
2483 c... ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2484 c **********************************************************************************************************
2485 
2486 c **********************************************************************************************************
2487  subroutine calc_rc_ref865(sza,xthet,xphi,ref_toa,ref_rc)
2488 c...
2489 c... ----------------------------------------------------------
2490 c... written by Jaehwa Lee
2491 c... ver. 1.0 Mar. 21, 2014
2492 c...
2493 c... subroutine to calculate Rayleigh-corrected reflectance
2494 c... code inherited from "calc_sfc21.f" by MJ
2495 c...
2496 c... Inputs
2497 c... sza : solar zenith angle
2498 c... xthet : viewing zenith angle
2499 c... xphi : relative azimuth angle
2500 c... ref_toa: TOA reflectance (I/F)
2501 c...
2502 c... Output
2503 c... ref_rc : Rayleigh-corrected reflectance (normalized)
2504 c...
2505 c... ----------------------------------------------------------
2506 
2507  include 'aottbl.inc'
2508  include 'sfc21tbl.inc'
2509 
2510 c real tau2(4), theta0(10), theta(46), phi(30)
2511  real tau2(4)
2512  real nnvalx(4,4,2), rr0x(4,4,2)
2513  real ttx(4,4,2), ssx(4,4,2)
2514  logical dflag3
2515 
2516  data pi /3.14159/
2517 
2518  x1 = sza
2519  x2 = xthet
2520  x3 = xphi
2521 
2522  mm = 10 ! solar zenith
2523  nn = 46 ! satellite zenith
2524  ll = 30 ! rel. azimuth
2525  ma2= 4 ! tau2 (aot at 2.1um)
2526 
2527 c-- interpolating for 2.1 micron channel
2528 
2529  nsm=1
2530  dif=x1-theta0(1)
2531  do i=1,mm
2532  dift=x1-theta0(i)
2533  if (dift.gt.0. .and. dift.lt.dif) then
2534  nsm=i
2535  dif=dift
2536  endif
2537  enddo
2538  mbeg = nsm - 2
2539  if (mbeg.le.0) then
2540  mbeg = 0
2541  else if (mbeg.gt.mm-4) then
2542  mbeg = mm-4
2543  endif
2544 
2545  nsn=1
2546  dif=x2-theta(1)
2547  do i=1,nn
2548  dift=x2-theta(i)
2549  if (dift.gt.0. .and. dift.lt.dif) then
2550  nsn=i
2551  dif=dift
2552  endif
2553  enddo
2554  nbeg = nsn - 2
2555  if (nbeg.le.0) then
2556  nbeg = 0
2557  else if (nbeg.gt.nn-4) then
2558  nbeg = nn-4
2559  endif
2560 
2561  call search(dflag3,x3,phi,ll,ii)
2562  xfrac = (x3-phi(ii))/(phi(ii+1)-phi(ii))
2563 
2564  do 123 i = 1, 4
2565  do 124 j = 1, 4
2566  nnvalx(i,j,1) = nvalx865(mbeg+i,nbeg+j,ii,1)
2567  nnvalx(i,j,2) = nvalx865(mbeg+i,nbeg+j,ii+1,1)
2568 
2569  rr0x(i,j,1) = r0x_865(mbeg+i,nbeg+j,ii,1)
2570  rr0x(i,j,2) = r0x_865(mbeg+i,nbeg+j,ii+1,1)
2571 
2572  ttx(i,j,1) = tx_865(mbeg+i,nbeg+j,ii,1)
2573  ttx(i,j,2) = tx_865(mbeg+i,nbeg+j,ii+1,1)
2574 
2575  ssx(i,j,1) = sx_865(mbeg+i,nbeg+j,ii,1)
2576  ssx(i,j,2) = sx_865(mbeg+i,nbeg+j,ii+1,1)
2577 
2578 124 continue
2579 123 continue
2580 
2581  call new_intepsf(theta0, theta, phi, nnvalx, mm, nn, ll,
2582  1 x1,x2,x3,y,dy,mbeg,nbeg,xfrac)
2583 
2584  rr = y/pi
2585 
2586  call new_intepsf(theta0, theta, phi, rr0x, mm, nn, ll,
2587  1 x1,x2,x3,y,dy,mbeg,nbeg,xfrac)
2588 
2589  rr0 = y/pi
2590 
2591  call new_intepsf(theta0, theta, phi, ttx, mm, nn, ll,
2592  1 x1,x2,x3,y,dy,mbeg,nbeg,xfrac)
2593 
2594  tt = y/pi
2595 
2596  call new_intepsf(theta0, theta, phi, ssx, mm, nn, ll,
2597  1 x1,x2,x3,y,dy,mbeg,nbeg,xfrac)
2598 
2599  ss = y
2600 
2601  dff = ref_toa - rr0
2602  ref_rc = dff / (tt + ss*dff)
2603 
2604  return
2605  end
2606 c... --------------------------------------
2607 c... The end of subroutine calc_rc_ref865
2608 c... ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2609 c **********************************************************************************************************
2610 
2611  subroutine new_intepsf(x1a,x2a,x3a,ya,m,n,l,x1,x2,x3,y,dy,
2612  1 mbeg,nbeg,frac)
2613 
2614  use viirs_aerosol_luts, only: polint
2615 
2616  parameter(nmax=46,mmax=10,lmax=30)
2617  dimension x1a(m),x2a(n),x3a(l),ya(4,4,2)
2618  dimension xx2a(4), xx1a(4)
2619  dimension yntmp(4),ymtmp(4),yltmp(2)
2620 
2621  do 12 j=1,4
2622  do 11 k=1,4
2623  yltmp(1)=ya(j,k,1)
2624  yltmp(2)=ya(j,k,2)
2625  yntmp(k) = yltmp(1)*(1.-frac) + yltmp(2)*frac
2626  xx2a(k) = x2a(k+nbeg)
2627 11 continue
2628  call polint(xx2a,yntmp,4,x2,ymtmp(j),dy)
2629  xx1a(j) = x1a(j+mbeg)
2630 12 continue
2631  call polint(xx1a,ymtmp,4,x1,y,dy)
2632  return
2633  end
2634 
calc_rc_ref865
subroutine calc_rc_ref865(sza, xthet, xphi, ref_toa, ref_rc)
Definition: find_v_vegset_nc4.f:2488
get_sfcrfl_veg_high
subroutine get_sfcrfl_veg_high(season, lcvr, refln21, rc_ndvi, xeAs_B1, xeAs_B3)
Definition: find_v_vegset_nc4.f:1087
calc_sfc21
subroutine calc_sfc21(sza, xthet, xphi, refl21, aot21, As_21)
Definition: find_v_vegset_nc4.f:2196
real
#define real
Definition: DbAlgOcean.cpp:26
find_v_veg
subroutine find_v_veg(month, season, realbuf, tmpvg, r412sv, r470sv, gzflg, outbufvg, tau_x470_flag)
Definition: find_v_vegset_nc4.f:3
new_intepsf
subroutine new_intepsf(x1a, x2a, x3a, ya, m, n, l, x1, x2, x3, y, dy, mbeg, nbeg, frac)
Definition: find_v_vegset_nc4.f:2613
aero_470veg
subroutine aero_470veg(refl, x1, x2, x3, mm, nn, ll, ma, imod, r470, tau_x470, tau_x470_flag)
Definition: find_v_vegset_nc4.f:1732
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
viirs_aerosol_luts
Definition: viirs_aerosol_luts_nc4.f95:1
search
subroutine search(dflag, xbar, x, n, i)
Definition: find_v_viirs_nc4.f:3401
search2
subroutine search2(dflag, xbar, x, n, i, fx)
Definition: find_v_viirs_nc4.f:3452
max
#define max(A, B)
Definition: main_biosmap.c:61
tt
subroutine tt(NMAX, NCHECK)
Definition: ampld.lp.f:1852
endif
endif() set(LIBS $
Definition: CMakeLists.txt:6
get_sfcrfl_nir_vis
subroutine get_sfcrfl_nir_vis(season, lcvr, refn865_rc, rc_ndvi, xeAs_B1, xeAs_B3)
Definition: find_v_vegset_nc4.f:1469
calc_rc_ref672
subroutine calc_rc_ref672(sza, xthet, xphi, ref_toa, ref_rc)
Definition: find_v_vegset_nc4.f:2363
viirs_aerosol_luts::aero_470
subroutine, public aero_470(dflag, refl, x1, x2, x3, mm, nn, ll, ma, imod, r470, tau_x470, tau_x470_flag, trflg, model_frac, debug)
Definition: viirs_aerosol_luts_nc4.f95:893
min
#define min(A, B)
Definition: main_biosmap.c:62
pi
#define pi
Definition: vincenty.c:23
viirs_aerosol_luts::polint
subroutine, public polint(xa, ya, n, x, y, dy)
Definition: viirs_aerosol_luts_nc4.f95:2503
eq
#define eq(a, b)
Definition: rice.h:167
viirs_aerosol_luts::new_intep
subroutine, public new_intep(x1a, x2a, x3a, ya, m, n, l, ia, x1, x2, x3, y, dy, mbeg, nbeg, frac)
Definition: viirs_aerosol_luts_nc4.f95:2461
aero_650veg
subroutine aero_650veg(refl, x1, x2, x3, mm, nn, ll, ma, r650, tau_x650, tau_x650_flag, tau_x470_flag, tau_x470)
Definition: find_v_vegset_nc4.f:1881
get_sfcrfl_veg
subroutine get_sfcrfl_veg(season, lcvr, refln21, rc_ndvi, xeAs_B1, xeAs_B3)
Definition: find_v_vegset_nc4.f:998
get_sfcrfl_swir_vis
subroutine get_sfcrfl_swir_vis(season, lcvr, refln21, rc_ndvi, xeAs_B1, xeAs_B3)
Definition: find_v_vegset_nc4.f:1176
viirs_aerosol_luts::aero_650
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)
Definition: viirs_aerosol_luts_nc4.f95:1021
cmb_stddb_veg
subroutine cmb_stddb_veg(outbuf, outbufvg, tmpvg, idalg)
Definition: find_v_vegset_nc4.f:2084


ABOUT

Our Team

Contact Us

DATA

Get Started

Find Data

Use Data

TOOLS

SeaBASS

SeaDAS

Giovanni Data Visualizer

Worldview Mapping Interface

Level 3 Time Series Plotter

Request a SeaHawk Image

Mission Quality Monitor

HELP

Level 1 & 2 Browser Help

Level 3 & 4 Browser Help

File Search Help

Search & Download Help

Earthdata Forum

Join Mailing List

How to Cite Data

CONNECT WITH US
COMPLIANCE

Web Privacy Policy

Data & Information Policy

Communications Policy

Accessibility Policy

Freedom of Information Act

USA.gov

Responsible NASA Official: Sean Bailey
Curator: OceanColor Webmaster