OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
afrt_rt2_read.f
Go to the documentation of this file.
1  subroutine af_rt2_read(adir1,adir2,adir3,adir4,adir5,ailm,aisd,airh,
2  1 aiwind,aitau,aiset,airef,aipol,aiww,akrhum,awwl,aalbwat,
3  2 ifunc,mfunc,refr,refi,rmin,rmax,delr,delx,rg,sig,npar,
4  3 r11,r22,r33,r44,reff,veff,ccn,bsr,salb,asf,qscat,qext,
5  4 t,thd,xr,xi,v,thcel,phcel,txx,pti,wvlth,psrfc,rho,
6  5 xozn,tautot,deltau,trp,tmp,tap,tcar,twat,tozn,
7  6 htlvl,pplvl,dtrr,dtmm,dtaa,dtot,htdv,pdv,
8  7 taur,taum,taua,ifc,nmodl,nolyr,bfr)
9 c
10 c subroutine readin reads the opens files and reads input from
11 c unit 21
12 c***********************************************************************
13 c.....include common blocks
14  implicit real*8 (a-h,o-z)
15 c include 'toa_rt2.cmn'
16 c parameter(nw=100,nth=26,nsz=25,nph=46,nsg=10,nstk=4,ntrx=16,npti=8,
17 c 1 nlyr=200,ntf=1801,nwnd=8)
18  parameter(nw=250,nth=91,nsz=90,nph=181,nsg=10,nstk=4,ntrx=16,npti=8,
19  1 nlyr=200,ntf=1801,nwnd=8,nmd=3)
20 c***********************************************************************
21  integer*4 izx(2),aiww(10),akrhum(10),ifc,nolyr,nmodl,ifunc,mfunc
22  integer*4 ailm,airh,aisd,aitau,aiwind,aiset,airef,aipol,ref
23  real*8 refr(nmd),refi(nmd),rmin(nmd),rmax(nmd),delr(nmd),delx(nmd),
24  1 rg(nmd),sig(nmd),npar(nmd),r11,r22,r33,r44,reff,veff,ccn,bsr,
25  2 salb,asf,qscat,qext,t(ntf,nstk),thd(ntf),txx(ntrx,nph,nth,nsz),
26  3 pti(npti,nph,nth,nsz),wvlth,psrfc,rho,xozn,
27  4 tautot,deltau,trp,tmp,tap,tcar,twat,tozn,
28  5 htlvl(nlyr),pplvl(nlyr),dtrr(nlyr),dtmm(nlyr),
29  6 dtaa(nlyr),dtot(nlyr),htdv(nlyr),pdv(nlyr),
30  7 taur(nlyr),taum(nlyr),taua(nlyr)
31  real*4 bfr(nlyr),aalbwat(nwnd)
32  real*4 bfr1(nlyr),bfr2(nlyr),bfr3(nlyr),gfr2(nlyr),gfr3(nlyr)
33  real*8 xzx(3),athcnd(nth),aphcnd(nph),ftrx(ntrx,nth,nph),
34  1 ptix(npti,nth,nph),awwl(nlyr)
35  real*8 solbrdf(nsz),thbrdf(nsz),phbrdf(nph),wbrdf(10)
36  real*8 brdfy(10),brdfx(nsz,nth,nph),fiosurf(4,2*nsz,nph)
37  real*8 surfzu(nsz,nth,nph)
38  integer*4 nsolbrdf,nthbrdf,nphbrdf,nwbrdf
39  character*255 adir0,adir1,adir2,adir3,adir4,adir5
40 
41 c
42 c character*1 ciwind,ciconc
43  character*1 ciconc
44  character*2 ciwind
45  character*2 citau,crh,cset
46  character*2 cisd,cilmd,cwind
47  character*4 cilm
48  character*90 xname1,xname2,xname3,xname4,xname5,file_name
49  character*90 xname4a,xname5a,xname6
50  character*80 txt,temp
51 c
52  len1=index(adir1,' ')-1
53  len2=index(adir2,' ')-1
54  len3=index(adir3,' ')-1
55  len4=index(adir4,' ')-1
56  len5=index(adir5,' ')-1
57 c
58 c if(awwl(ilm).lt.1.0)then
59  iwl=awwl(ilm)*1.0d3+0.01d0
60 c else
61 c iwl=awwl(ilm)*1.0d4+0.1d0
62 c endif
63 
64  iwind = aiwind
65  iset = aiset
66  call convtc(idnint(awwl(ailm)*1.0d3+0.01),4,cilm)
67 c write(*,*) awwl(ailm),cilm
68  call convtc(aisd,2,cisd)
69  call convtc(aitau,2,citau)
70  call convtc(aiwind,2,ciwind)
71 c call convtc(iconc,1,ciconc)
72  call convtc(aiww(iwind),2,cwind)
73  call convtc(akrhum(airh),2,crh)
74  call convtc(aiset,2,cset) !iset is passed through the common block zeroa
75 c
76  if(airef.eq.1)ciwind='00'
77  if(airef.eq.0)ciwind='99'
78  if(aiset.eq.0)crh='00' !iset is read in sub. inusr and passed through common zeroa
79 c
80  xname1='rt2_wl'//cilm//'x'//ciwind
81  xname2='sd'//cisd//'rh'//crh//'ta'//citau//'_set'//cset
82 c xname3=cilm//ciwind//ciconc//cisd//citau
83  xname4a='rt1_'//'wl'//cilm//'sd'//cisd//'rh'//crh//'ta'//citau
84  nm4a=index(xname4a,' ')-1
85  xname4=xname4a(1:nm4a)//'_set'//cset
86  xname5a='rt1_sp'//'wl'//cilm//'sd'//cisd//'rh'//crh//'ta'//citau
87  nm5a=index(xname5a,' ')-1
88  xname5=xname5a(1:nm5a)//'_set'//cset
89 c
90 c
91  nm1=index(xname1,' ')-1
92  nm2=index(xname2,' ')-1
93 c nm3=index(xname3,' ')-1
94  nm4=index(xname4,' ')-1
95  nm5=index(xname5,' ')-1
96 c
97 c
98  num_arg = iargc()
99  if (num_arg .lt. 1) then
100  file_name = 'dumpza'
101  else
102  call getarg(1, file_name)
103  end if
104 c
105 c open(6,file=file_name,status='unknown')
106  open(21,file=adir1(1:len1)//'/'//xname4(1:nm4)//'.dat',
107  1 status='old')
108 c
109 c rewind(7)
110 c read(7,*)xday,xmonth,xyear
111  xday=-777.0
112  xmonth=-777.0
113  xyear=-777.0
114 c
115  read(21,125)temp
116  read(21,126)wvlth,psrfc,rho,c,tautot,deltau,nolyr
117  read(21,125)temp
118  read(21,126)trp,tmp,tap,tcar,twat,tozn,ifc
119  read(21,125)temp
120  htlvl(1)=84.0d0
121  pplvl(1)=5.3040d-6
122  do i=1,nolyr
123  read(21,128)htlvl(i+1),pplvl(i+1),dtrr(i),dtmm(i),dtaa(i),dtot(i)
124 c write(6,128)ht(i+1),pl(i+1),dtrr(i),dtmm(i),dtaa(i),dtot(i)
125  enddo
126  press=pplvl(nolyr+1)
127 c
128  if(ipsudo.eq.1)then
129  open(2,file=
130  1 adir1(1:len1)//'/'//xname5(1:nm5)//'.dat',
131  2 status='old')
132  read(2,125)temp
133  read(2,145)nmodl
134  read(2,125)temp
135  do i=1,nmodl
136  read(2,150)htdv(i),pdv(i),taur(i),taum(i),taua(i)
137  enddo
138 c
139  endif
140 c
141  qscat=1.0d-15
142  qext=1.0d-15
143  do 70 i=1,1801
144  do 60 j=1,4
145  t(i,j) = 0.0d0
146  60 continue
147  70 continue
148  if(ifc.eq.1)then
149 c
150  xname6='phs_twl'//cilm//'sd'//cisd//'rh'//crh//'_set'//cset
151  nm6=index(xname6,' ')-1
152 
153  open(5,file=
154  1 adir2(1:len2)//'/'//xname6(1:nm6)//'.dat',
155  1 status='old')
156 c
157  read (5,1090) txt
158  read (5,1090) txt
159  read (5,1125) ifunc,mfunc
160  if(ifunc.eq.1)then
161  read (5,1090) txt
162  read (5,1100) wvlth2,refr(1),refi(1)
163  read (5,1090) txt
164  read (5,1100) rmin(1),rmax(1)
165  read (5,1090) txt
166  read (5,1100) delr(1),delx(1)
167  read (5,1090) txt
168  read (5,1100) rc,cc,xnu
169  elseif(ifunc.eq.2)then
170  read (5,1090) txt
171  read (5,1100) wvlth2,refr(1),refi(1)
172  read (5,1090) txt
173  read (5,1100) rmin(1),rmax(1)
174  read (5,1090) txt
175  read (5,1100) delr(1),delx(1)
176  read (5,1090) txt
177  read (5,1100) xalf,xgam,xa,xb
178  elseif(ifunc.eq.3)then
179  read (5,1090) txt
180  read (5,1100) wvlth2,refr(1),refi(1)
181  read (5,1090) txt
182  read (5,1100) rmin(1),rmax(1)
183  read (5,1090) txt
184  read (5,1100) delr(1),delx(1)
185  read (5,1090) txt
186  read (5,1100) rg(1),sig(1),npar(1)
187  elseif(ifunc.eq.4)then
188  read (5,1090) txt
189  read (5,1100) wvlth2,refr(1),refi(1),refr(2),refi(2)
190  read (5,1090) txt
191  read (5,1100) rmin(1),rmax(1),rmin(2),rmax(2)
192  read (5,1090) txt
193  read (5,1100) delr(1),delx(1),delr(2),delx(2)
194  read (5,1090) txt
195  read (5,1100) rg(1),sig(1),npar(1),rg(2),sig(2),npar(2)
196  elseif(ifunc.eq.5)then
197  read (5,1090) txt
198  read (5,1100) wvlth2,refr(1),refi(1),refr(2),refi(2),refr(3),refi(3)
199  read (5,1090) txt
200  read (5,1100) rmin(1),rmax(1),rmin(2),rmax(2),rmin(3),rmax(3)
201  read (5,1090) txt
202  read (5,1100) delr(1),delx(1),delr(2),delx(2),delr(3),delx(3)
203  read (5,1090) txt
204  read (5,1130) rg(1),sig(1),npar(1),rg(2),sig(2),npar(2),rg(3),sig(3),npar(3)
205  endif
206  read (5,1090) txt
207  read (5,1100) r11,r22,r33,r44,reff,veff
208  read (5,1090) txt
209  read (5,1095) ccn,bsr,salb,asf,qscat,qext
210 c convert the wavelength in cm
211  wvlth2=wvlth2*1.0d-4
212 c
213  dw=abs(wvlth2-wvlth)
214 cbaf if(dw .ge.1.0e-12)then
215  if(dw .ge.2.0e-8)then
216 c write(6,1105)wvlth,wvlth2
217  write(*,*) "Wavelength mismatch ",dw,wvlth,wvlth2
218  stop
219  endif
220 c
221 c....read the prevoiusly computed t matrix
222  do i=1,1801
223  read(5,1120)(t(i,j),j=1,4),thd(i)
224  enddo
225 c
226  endif
227  if(airef.eq.2)then
228 c open(18,
229 c 1 access='sequential',form='unformatted',status='scratch')
230 c open(19,
231 c 1 access='sequential',form='unformatted',status='scratch')
232  open(20,file=
233  1 adir4(1:len4)//'/'//'ocn_refl_wl'//cilm//'x'//cwind//'.dat',
234  2 access='sequential',form='unformatted',status='old')
235 c read the rough ocean inputs
236  rewind(20)
237  read (20) xzx,izx,athcnd,aphcnd
238  xr = xzx(2)
239  xi = xzx(3)
240  v = xzx(1)
241  thcel = izx(1)
242  phcel = izx(2)
243 c write(*,*) adir4(1:len4)//'/'//'ocn_refl_wl'//cilm//'x'//cwind//'.dat'
244 c write(*,901)xzx
245 c write(*,902)izx
246 c write(*,903)athcnd
247 c write(*,904)aphcnd
248 901 format('xzx',3f8.3)
249 902 format('izx',2i5)
250 903 format('athcnd',10f7.1)
251 904 format('aphcnd',10f7.1)
252  mtha = izx(1)
253  mphi = izx(2)
254 c write (19) (xzx(i),i=1,3),mtha,mphi
255 c
256  do l=1,mtha
257  read (20) thetap,ftrx,ptix
258 c write (18) thetap,(((ftrx(i,j,k),i=1,16),j=1,mtha),
259 c 1 k=1,mphi)
260  do k=1,mphi
261  do j=1,mtha
262  do i=1,16
263  txx(i,k,j,l)=ftrx(i,j,k)
264  enddo
265  do i=1,8
266  pti(i,k,j,l)=ptix(i,j,k)
267  enddo
268  enddo
269  enddo
270 c
271 c write (19) (((ptix(i,j,k),i=1,8),j=1,mtha),k=1,mphi)
272  enddo
273 c rewind(18)
274 c rewind(19)
275  rewind(20)
276  endif
277 c
278 c if airef=3 then read the brdf of the ground surface
279 c
280  if(airef.eq.3)then
281  open(17,file=adir0(1:len0)//'/'//'brdf.dat',
282  1 status='old')
283 c
284  read(17,'(a)')temp
285  read(17,*)nsolbrdf,nthbrdf,nphbrdf,nwbrdf
286  do is=1,nsolbrdf
287  read(17,160)solbrdf(is)
288 160 format(t6,f4.1)
289  read(17,165)(wbrdf(i),i=1,nwbrdf)
290 165 format(t22,f5.1,t35,f5.1,t48,f5.1,t61,f5.1,
291  1 t74,f5.1,t87,f5.1)
292 c check the input wavelength against the wavelengths
293 c in wbrdf
294  mwbrdf=0
295  do i=1,nwbrdf
296  diffwbrdf=dabs(wvlth*1.0d7-wbrdf(i))
297  if(diffwbrdf.lt.0.1d-3)mwbrdf=i
298  enddo
299  if(mwbrdf.eq.0)then
300 c write(6,*)'input wavelength does not matches',
301 c 1 ' with the wavelengths in the buffer',
302 c 2 ' wbrdf'
303 c write(6,*)'input wavelength',wvlth*1.0d7,' nm'
304 c write(6,*)'wavelengths in the buffer wbrdf:',wbrdf
305 c call closeunits
306  stop
307  endif
308  do j=1,nphbrdf
309  m=nphbrdf-j+1
310  do k=1,nthbrdf
311  read(17,170)thbrdf(k),phbrdf(j),
312  1 (brdfy(l),l=1,nwbrdf)
313 170 format(f6.2,2x,f6.2,10d13.4)
314  brdfx(is,k,m)=brdfy(mwbrdf)
315  enddo
316  enddo
317  enddo
318 c
319  endif
320 c
321 c initialize the buffer & fill the header block
322 c
323  if(aipol.eq.0)rho=0.0d0
324 c
325  do i=1,nlyr
326  bfr1(i)=-777.0
327  bfr2(i)=-777.0
328  bfr3(i)=-777.0
329  enddo
330 c
331  bfr1(1)=xyear
332  bfr1(2)=xmonth
333  bfr1(3)=xday
334  bfr1(4)=-777.0
335  bfr1(5)=wvlth*1.0d4
336  bfr1(6)=refr(1)
337  bfr1(7)=refi(1)
338  bfr1(8)=rg(1)
339  bfr1(9)=sig(1)
340  bfr1(10)=rmin(1)
341  bfr1(11)=rmax(1)
342  bfr1(12)=delx(1)
343  bfr1(13)=reff
344  bfr1(14)=r11
345  bfr1(15)=r22
346  bfr1(16)=r33
347  bfr1(17)=r44
348  bfr1(18)=salb
349  bfr1(19)=asf
350  bfr1(20)=qscat
351  bfr1(21)=qext
352  bfr1(22)=xzx(3)
353  bfr1(23)=xzx(1)
354  bfr1(24)=xzx(2)
355  bfr1(25)=trp
356  bfr1(26)=tmp
357  bfr1(27)=tap
358  bfr1(28)=twat
359  bfr1(29)=tautot
360  bfr1(30)=rho
361  bfr1(31)=ccn
362  bfr1(32)=bsr
363  bfr1(33)=float(ipol)
364  bfr1(34)=float(ipsudo)
365  bfr1(41)=dfloat(ifoam)
366  bfr1(42)=dfloat(iwatr)
367 c bfr1(43)=chl(iconc)
368  bfr1(44)=aalbwat(ilm)
369  bfr1(45)=dfloat(airef)
370  bfr1(46)=psrfc
371  bfr1(47)=dfloat(ifc)
372  bfr1(48)=dfloat(ifunc)
373  bfr1(49)=dfloat(mfunc)
374  goto(210,220,230,230,230),ifunc
375 210 continue
376  bfr1(51)=rc
377  bfr1(52)=cc
378  bfr1(53)=xnu
379  goto 260
380 220 continue
381  bfr1(54)=xalf
382  bfr1(55)=xgam
383  bfr1(56)=xa
384  bfr1(57)=xb
385  goto 260
386 230 continue
387  bfr1(58)=npar(1)
388  if(ifunc.eq.3)goto 260
389  bfr1(59)=rg(2)
390  bfr1(60)=sig(2)
391  bfr1(61)=npar(2)
392  bfr1(62)=delx(2)
393  bfr1(63)=refr(2)
394  bfr1(64)=refi(2)
395  bfr1(65)=rmin(2)
396  bfr1(66)=rmax(2)
397  if(ifunc.eq.4)goto 260
398  bfr1(67)=rg(3)
399  bfr1(68)=sig(3)
400  bfr1(69)=npar(3)
401  bfr1(70)=delx(3)
402  bfr1(71)=refr(3)
403  bfr1(72)=refi(3)
404  bfr1(73)=rmin(3)
405  bfr1(74)=rmax(3)
406 260 continue
407  bfr1(76)=itrans !flag for computing diffused transmiss
408 
409  do i=1,100
410  bfr(i)=bfr1(i)
411  enddo
412 
413  close (2)
414  close (5)
415  close (17)
416  close (20)
417  close (21)
418 
419 c
420 c*****format statements******************************************************
421  100 format(2i4,1x,1p6e11.3)
422  125 format(a)
423  126 format(6d12.4,i4)
424  128 format(4x,6d12.4)
425  145 format(i5)
426  150 format(5x,1p6e12.4)
427  195 format (1x,'xzx',1p10d11.2)
428  196 format (1x,'izx',2i5)
429  197 format (1x,'athcnd',1p10d11.2)
430  198 format (1x,'aphcnd',1p10d11.2)
431  199 format (1x,'thetap',1p10d11.2)
432  200 format (1x,'ftrx',1p10d11.2)
433 c 205 format('xzx1-3,mtha,mphi'/1p3e11.3,0p2i4)
434  201 format (1x,'pti',1p10d11.2)
435  505 format (1x,'end of file encountered , check the program')
436  1090 format(a)
437  1095 format(1p4e11.3,1p2e12.4)
438  1100 format(1p7e11.3)
439  1105 format(1x,'input wavelength and phase matrix wavelength'
440  1 ' do not match.'/1x,'input wavelength=',1pe12.4,
441  2 'phase matrix wavelength=',1pe12.4)
442  1120 format(1p2d18.12,1p2d19.12,0pf6.2)
443  1125 format(i5,11x,i5)
444  1130 format(9f8.5)
445  1160 format(a)
446 c****************************************************************************
447  return
448  end
449 c***********************************************************************
an array had not been initialized Several spelling and grammar corrections were which is read from the appropriate MCF the above metadata values were hard coded A problem calculating the average background DN for SWIR bands when the moon is in the space view port was corrected The new algorithm used to calculate the average background DN for all reflective bands when the moon is in the space view port is now the same as the algorithm employed by the thermal bands For non SWIR changes in the averages are typically less than Also for non SWIR the black body DNs remain a backup in case the SV DNs are not available For SWIR the changes in computed averages were larger because the old which used the black body suffered from contamination by the micron leak As a consequence of the if SV DNs are not available for the SWIR the EV pixels will not be the granule time is used to identify the appropriate tables within the set given for one LUT the first two or last two tables respectively will be used for the interpolation If there is only one LUT in the set of it will be treated as a constant LUT The manner in which Earth View data is checked for saturation was changed Previously the raw Earth View DNs and Space View DNs were checked against the lookup table values contained in the table dn_sat The change made is to check the raw Earth and Space View DNs to be sure they are less than the maximum saturation value and to check the Space View subtracted Earth View dns against a set of values contained in the new lookup table dn_sat_ev The metadata configuration and ASSOCIATEDINSTRUMENTSHORTNAME from the MOD02HKM product The same metatdata with extensions and were removed from the MOD021KM and MOD02OBC products ASSOCIATEDSENSORSHORTNAME was set to MODIS in all products These changes are reflected in new File Specification which users may consult for exact the pow functions were eliminated in Emissive_Cal and Emissive bands replaced by more efficient code Other calculations throughout the code were also made more efficient Aside from a few round off there was no difference to the product The CPU time decreased by about for a day case and for a night case A minor bug in calculating the uncertainty index for emissive bands was corrected The frame index(0-based) was previously being used the frame number(1-based) should have been used. There were only a few minor changes to the uncertainty index(maximum of 1 digit). 3. Some inefficient arrays(Sigma_RVS_norm_sq) were eliminated and some code lines in Preprocess_L1A_Data were moved into Process_OBCEng_Emiss. There were no changes to the product. Required RAM was reduced by 20 MB. Now
#define real
Definition: DbAlgOcean.cpp:26
subroutine af_rt2_read(adir1, adir2, adir3, adir4, adir5, ailm, aisd, air aiwind, aitau, aiset, airef, aipol, aiww, akrhum, awwl, aalbwat, ifunc, mfunc, refr, refi, rmin, rmax, delr, delx, rg, sig, npar, r11, r22, r33, r44, reff, veff, ccn, bsr, salb, asf, qscat, qext, t, thd, xr, xi, v, thcel, phcel, txx, pti, wvlth, psrfc, rho, xozn, tautot, deltau, trp, tmp, tap, tcar, twat, tozn, htlvl, pplvl, dtrr, dtmm, dtaa, dtot, htdv, pdv, taur, taum, taua, ifc, nmodl, nolyr, bfr)
Definition: afrt_rt2_read.f:9
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
float rc(float x, float y)
#define abs(a)
Definition: misc.h:90
subroutine convtc(num, nchar, loc)
Definition: ocn.f:156