OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
readin.f
Go to the documentation of this file.
1  subroutine readin(ilm,irh,isd,itau)
2 c
3 c subroutine readin reads the opens files and reads input from
4 c unit 21
5 c***********************************************************************
6 c.....include common blocks
7  implicit real*8 (a-h,o-z)
8  include 'common_all.cmn'
9 c***********************************************************************
10  integer*4 izx(2)
11  real*8 xzx(3),athcnd(25),aphcnd(46),ftrx(16,25,46),
12  1 ptix(8,25,46)
13 c
14 c character*1 ciwind,ciconc
15  character*1 ciconc
16  character*2 ciwind
17  character*2 citau,crh,cset
18  character*2 cisd,cilmd,cwind
19  character*4 cilm
20  character*90 xname1,xname2,xname3,xname4,xname5,file_name
21  character*90 xname4a,xname5a,xname6
22  character*80 txt,temp
23 c
24  len1=index(dir1,' ')-1
25  len2=index(dir2,' ')-1
26  len3=index(dir3,' ')-1
27  len4=index(dir4,' ')-1
28  len5=index(dir5,' ')-1
29 c
30 c if(wwl(ilm).lt.1.0)then
31  iwl=wwl(ilm)*1.0d3+0.01d0
32 c else
33 c iwl=wwl(ilm)*1.0d4+0.1d0
34 c endif
35  call convtc(idnint(wwl(ilm)*1.0d3+0.01),4,cilm)
36  write(*,*) wwl(ilm),cilm
37  call convtc(isd,2,cisd)
38  call convtc(itau,2,citau)
39  call convtc(iwind-1,2,ciwind)
40 c call convtc(iconc,1,ciconc)
41  call convtc(iww(iwind),2,cwind)
42  call convtc(krhum(irh),2,crh)
43  call convtc(iset,2,cset) !iset is passed through the common block zeroa
44 c
45  if(iref.eq.1)ciwind='00'
46  if(iref.eq.0)ciwind='99'
47  if(iset.eq.0)crh='00' !iset is read in sub. inusr and passed through common zeroa
48 c
49  xname1='rt2_wl'//cilm//'x'//ciwind
50  xname2='sd'//cisd//'rh'//crh//'ta'//citau//'_set'//cset
51 c xname3=cilm//ciwind//ciconc//cisd//citau
52  xname4a='rt1_'//'wl'//cilm//'sd'//cisd//'rh'//crh//'ta'//citau
53  nm4a=index(xname4a,' ')-1
54  xname4=xname4a(1:nm4a)//'_set'//cset
55  xname5a='rt1_sp'//'wl'//cilm//'sd'//cisd//'rh'//crh//'ta'//citau
56  nm5a=index(xname5a,' ')-1
57  xname5=xname5a(1:nm5a)//'_set'//cset
58 c
59 c
60  nm1=index(xname1,' ')-1
61  nm2=index(xname2,' ')-1
62 c nm3=index(xname3,' ')-1
63  nm4=index(xname4,' ')-1
64  nm5=index(xname5,' ')-1
65 c
66 c
67  num_arg = iargc()
68  if (num_arg .lt. 1) then
69  file_name = 'dumpza'
70  else
71  call getarg(1, file_name)
72  end if
73 c
74  open(3,file=
75  1 dir3(1:len3)//'/'//xname1(1:nm1)//xname2(1:nm2)//'dn.dat',
76  2 status='unknown')
77  open(4,file=
78  1 dir3(1:len3)//'/'//xname1(1:nm1)//xname2(1:nm2)//'up.dat',
79  2 status='unknown')
80  open(6,file=file_name,status='unknown')
81  open(21,file=dir1(1:len1)//'/'//xname4(1:nm4)//'.dat',
82  1 status='old')
83  open(24,file=
84  1 dir3(1:len3)//'/'//xname1(1:nm1)//xname2(1:nm2)//'.dat',
85  2 access='direct',recl=1915*4,
86  3 form='unformatted',status='unknown')
87 c
88  open(53,access='direct',
89  1 form='unformatted',status='scratch',recl=18400*4)
90  open(54,access='direct',
91  1 form='unformatted',status='scratch',recl=18400*4)
92  open(55,access='direct',
93  1 form='unformatted',status='scratch',recl=18400*4)
94  open(64,access='direct',
95  1 form='unformatted',status='scratch',recl=18400*4)
96  open(65,access='direct',
97  1 form='unformatted',status='scratch',recl=18400*4)
98  open(71,access='direct',
99  1 form='unformatted',status='scratch',recl=18400*4)
100  open(72,access='direct',
101  1 form='unformatted',status='scratch',recl=18400*4)
102 c
103  if(icrft.eq.1)then
104  open(13,file=
105  1 dir3(1:len3)//'/'//xname1(1:nm1)//xname2(1:nm2)//'crftdn.dat',
106  2 status='unknown')
107  open(14,file=
108  1 dir3(1:len3)//'/'//xname1(1:nm1)//xname2(1:nm2)//'crftup.dat',
109  2 status='unknown')
110  open(15,file=
111  1 dir3(1:len3)//'/'//xname1(1:nm1)//xname2(1:nm2)//'crft.dat',
112  2 access='direct',recl=1915*4,
113  3 form='unformatted',status='unknown')
114  open(73,access='direct',
115  1 form='unformatted',status='scratch',recl=18400*4)
116  open(74,access='direct',
117  1 form='unformatted',status='scratch',recl=18400*4)
118  endif
119 c
120  if(iactflx.eq.1)then
121  open(30,file=
122  1 dir3(1:len3)//'/'//xname1(1:nm1)//xname2(1:nm2)//'flx.dat',
123  2 status='unknown')
124  endif
125 c
126  if(isurf.eq.1)then
127  open(16,file=
128  1 dir3(1:len3)//'/'//xname1(1:nm1)//xname2(1:nm2)//'surfup.dat',
129  2 status='unknown')
130  if(iref.eq.2)then
131  open(41,file=
132  1 dir3(1:len3)//'/'//xname1(1:nm1)//xname2(1:nm2)//'dirup.dat',
133  2 status='unknown')
134  open(42,file=
135  1 dir3(1:len3)//'/'//xname1(1:nm1)//xname2(1:nm2)//'ocnup.dat',
136  2 status='unknown')
137  open(43,file=
138  1 dir3(1:len3)//'/'//xname1(1:nm1)//xname2(1:nm2)//'skyup.dat',
139  2 status='unknown')
140  endif
141  endif
142 c
143 c rewind(7)
144 c read(7,*)xday,xmonth,xyear
145  xday=-777.0
146  xmonth=-777.0
147  xyear=-777.0
148 c
149  read(21,125)temp
150  read(21,126)wvlth,psrfc,rho,xozn,tautot,deltau,nolyr
151  read(21,125)temp
152  read(21,126)tr,tm,ta,tcar,twat,tozn,ifc
153  read(21,125)temp
154  ht(1)=84.0d0
155  pl(1)=5.3040d-6
156  do i=1,nolyr
157  read(21,128)ht(i+1),pl(i+1),dtrr(i),dtmm(i),dtaa(i),dtot(i)
158 c write(6,128)ht(i+1),pl(i+1),dtrr(i),dtmm(i),dtaa(i),dtot(i)
159  enddo
160  press=pl(nolyr+1)
161 c
162  if(ipsudo.eq.1)then
163  open(2,file=
164  1 dir1(1:len1)//'/'//xname5(1:nm5)//'.dat',
165  2 status='old')
166  read(2,125)temp
167  read(2,145)nmodl
168  read(2,125)temp
169  do i=1,nmodl
170  read(2,150)htp(i),ppo(i),taur(i),taum(i),tauabs(i)
171  enddo
172 c
173  endif
174 c
175  qs=1.0d-15
176  qt=1.0d-15
177  do 70 i=1,1801
178  do 60 j=1,4
179  t(i,j) = 0.0d0
180  60 continue
181  70 continue
182  if(ifc.eq.1)then
183 c
184  xname6='phs_twl'//cilm//'sd'//cisd//'rh'//crh//'_set'//cset
185  nm6=index(xname6,' ')-1
186 
187  open(5,file=
188  1 dir2(1:len2)//'/'//xname6(1:nm6)//'.dat',
189  1 status='old')
190 c
191  read (5,1090) txt
192  read (5,1090) txt
193  read (5,1125) ifunc,mfunc
194  if(ifunc.eq.1)then
195  read (5,1090) txt
196  read (5,1100) wvlth2,refr1,refi1
197  read (5,1090) txt
198  read (5,1100) rmin1,rmax1
199  read (5,1090) txt
200  read (5,1100) delr1,delx1
201  read (5,1090) txt
202  read (5,1100) rc,cc,xnu
203  elseif(ifunc.eq.2)then
204  read (5,1090) txt
205  read (5,1100) wvlth2,refr1,refi1
206  read (5,1090) txt
207  read (5,1100) rmin1,rmax1
208  read (5,1090) txt
209  read (5,1100) delr1,delx1
210  read (5,1090) txt
211  read (5,1100) xalf,xgam,xa,xb
212  elseif(ifunc.eq.3)then
213  read (5,1090) txt
214  read (5,1100) wvlth2,refr1,refi1
215  read (5,1090) txt
216  read (5,1100) rmin1,rmax1
217  read (5,1090) txt
218  read (5,1100) delr1,delx1
219  read (5,1090) txt
220  read (5,1100) rg1,sig1,xnum1
221  elseif(ifunc.eq.4)then
222  read (5,1090) txt
223  read (5,1100) wvlth2,refr1,refi1,refr2,refi2
224  read (5,1090) txt
225  read (5,1100) rmin1,rmax1,rmin2,rmax2
226  read (5,1090) txt
227  read (5,1100) delr1,delx1,delr2,delx2
228  read (5,1090) txt
229  read (5,1100) rg1,sig1,xnum1,rg2,sig2,xnum2
230  elseif(ifunc.eq.5)then
231  read (5,1090) txt
232  read (5,1100) wvlth2,refr1,refi1,refr2,refi2,refr3,refi3
233  read (5,1090) txt
234  read (5,1100) rmin1,rmax1,rmin2,rmax2,rmin3,rmax3
235  read (5,1090) txt
236  read (5,1100) delr1,delx1,delr2,delx2,delr3,delx3
237  read (5,1090) txt
238  read (5,1130) rg1,sig1,xnum1,rg2,sig2,xnum2,rg3,sig3,xnum3
239  endif
240  read (5,1090) txt
241  read (5,1100) r11,r22,r33,r44,reff,veff
242  read (5,1090) txt
243  read (5,1095) ccn,bsr,salb,asf,qs,qt
244 c convert the wavelength in cm
245  wvlth2=wvlth2*1.0d-4
246 c
247  dw=abs(wvlth2-wvlth)
248 cbaf if(dw .ge.1.0e-12)then
249  if(dw .ge.2.0e-8)then
250 c write(6,1105)wvlth,wvlth2
251  write(*,*) "Wavelength mismatch ",dw,wvlth,wvlth2
252  stop
253  endif
254 c
255 c....read the prevoiusly computed t matrix
256  do i=1,1801
257  read(5,1120)(t(i,j),j=1,4),thd(i)
258  enddo
259 c
260  endif
261  if(iref.eq.2)then
262  open(18,
263  1 access='sequential',form='unformatted',status='scratch')
264  open(19,
265  1 access='sequential',form='unformatted',status='scratch')
266 c
267 c open(20,file=
268 c 1 dir4(1:len4)//'/'//'ocn_refl_wl'//cilm//'x'//cwind//'.dat',
269 c 2 access='sequential',form='unformatted',status='old')
270 c
271  open(20,file=
272  1 dir4(1:len4)//'/'//'rufrtwl'//cilm//'x'//cwind//'solz00.dat',
273  2 access='sequential',form='unformatted',status='old')
274 c
275 c read the rough ocean inputs
276  rewind(20)
277  read (20) xzx,izx,athcnd,aphcnd
278  write(*,*)
279  1 dir4(1:len4)//'/'//'ocn_refl_wl'//cilm//'x'//cwind//'.dat'
280  write(*,901)xzx
281  write(*,902)izx
282  write(*,903)athcnd
283  write(*,904)aphcnd
284 901 format('xzx',3f8.3)
285 902 format('izx',2i5)
286 903 format('athcnd',10f7.1)
287 904 format('aphcnd',10f7.1)
288  mtha = izx(1)
289  mphi = izx(2)
290  write (19) (xzx(i),i=1,3),mtha,mphi
291 c
292  do l=1,mtha
293  read (20) thetap,ftrx,ptix
294 c write (18) thetap,(((ftrx(i,j,k),i=1,16),j=1,mtha),
295 c 1 k=1,mphi)
296  do k=1,mphi
297  do j=1,mtha
298  do i=1,16
299  txx(i,k,l,j)=ftrx(i,j,k)
300  enddo
301  enddo
302  enddo
303 c
304 c write (19) (((ptix(i,j,k),i=1,8),j=1,mtha),k=1,mphi)
305  enddo
306  rewind(18)
307  rewind(19)
308  rewind(20)
309  endif
310 c
311 c if iref=3 then read the brdf of the ground surface
312 c
313  if(iref.eq.3)then
314  open(17,file=dir0(1:len0)//'/'//'brdf.dat',
315  1 status='old')
316 c
317  read(17,'(a)')temp
318  read(17,*)nsolbrdf,nthbrdf,nphbrdf,nwbrdf
319  do is=1,nsolbrdf
320  read(17,160)solbrdf(is)
321 160 format(t6,f4.1)
322  read(17,165)(wbrdf(i),i=1,nwbrdf)
323 165 format(t22,f5.1,t35,f5.1,t48,f5.1,t61,f5.1,
324  1 t74,f5.1,t87,f5.1)
325 c check the input wavelength against the wavelengths
326 c in wbrdf
327  mwbrdf=0
328  do i=1,nwbrdf
329  diffwbrdf=dabs(wvlth*1.0d7-wbrdf(i))
330  if(diffwbrdf.lt.0.1d-3)mwbrdf=i
331  enddo
332  if(mwbrdf.eq.0)then
333 c write(6,*)'input wavelength does not matches',
334 c 1 ' with the wavelengths in the buffer',
335 c 2 ' wbrdf'
336 c write(6,*)'input wavelength',wvlth*1.0d7,' nm'
337 c write(6,*)'wavelengths in the buffer wbrdf:',wbrdf
338  call closeunits
339  stop
340  endif
341  do j=1,nphbrdf
342  m=nphbrdf-j+1
343  do k=1,nthbrdf
344  read(17,170)thbrdf(k),phbrdf(j),
345  1 (brdfy(l),l=1,nwbrdf)
346 170 format(f6.2,2x,f6.2,10d13.4)
347  brdfx(is,k,m)=brdfy(mwbrdf)
348  enddo
349  enddo
350  enddo
351 c
352  endif
353 c
354 c initialize the buffer & fill the header block
355 c
356  if(ipol.eq.0)rho=0.0d0
357 c
358  do i=1,4625
359  bfr1(i)=-777.0
360  bfr2(i)=-777.0
361  bfr3(i)=-777.0
362  enddo
363 c
364  bfr1(1)=xyear
365  bfr1(2)=xmonth
366  bfr1(3)=xday
367  bfr1(4)=-777.0
368  bfr1(5)=wvlth*1.0d4
369  bfr1(6)=refr1
370  bfr1(7)=refi1
371  bfr1(8)=rg1
372  bfr1(9)=sig1
373  bfr1(10)=rmin1
374  bfr1(11)=rmax1
375  bfr1(12)=delx1
376  bfr1(13)=reff
377  bfr1(14)=r11
378  bfr1(15)=r22
379  bfr1(16)=r33
380  bfr1(17)=r44
381  bfr1(18)=salb
382  bfr1(19)=asf
383  bfr1(20)=qs
384  bfr1(21)=qt
385  bfr1(22)=xzx(3)
386  bfr1(23)=xzx(1)
387  bfr1(24)=xzx(2)
388  bfr1(25)=tr
389  bfr1(26)=tm
390  bfr1(27)=ta
391  bfr1(28)=twat
392  bfr1(29)=tautot
393  bfr1(30)=rho
394  bfr1(31)=ccn
395  bfr1(32)=bsr
396  bfr1(33)=float(ipol)
397  bfr1(34)=float(ipsudo)
398  bfr1(41)=dfloat(ifoam)
399  bfr1(42)=dfloat(iwatr)
400 c bfr1(43)=chl(iconc)
401  bfr1(44)=albwat(ilm)
402  bfr1(45)=dfloat(iref)
403  bfr1(46)=psrfc
404  bfr1(47)=dfloat(ifc)
405  bfr1(48)=dfloat(ifunc)
406  bfr1(49)=dfloat(mfunc)
407  goto(210,220,230,230,230),ifunc
408 210 continue
409  bfr1(51)=rc
410  bfr1(52)=cc
411  bfr1(53)=xnu
412  goto 260
413 220 continue
414  bfr1(54)=xalf
415  bfr1(55)=xgam
416  bfr1(56)=xa
417  bfr1(57)=xb
418  goto 260
419 230 continue
420  bfr1(58)=xnum1
421  if(ifunc.eq.3)goto 260
422  bfr1(59)=rg2
423  bfr1(60)=sig2
424  bfr1(61)=xnum2
425  bfr1(62)=delx2
426  bfr1(63)=refr2
427  bfr1(64)=refi2
428  bfr1(65)=rmin2
429  bfr1(66)=rmax2
430  if(ifunc.eq.4)goto 260
431  bfr1(67)=rg3
432  bfr1(68)=sig3
433  bfr1(69)=xnum3
434  bfr1(70)=delx3
435  bfr1(71)=refr3
436  bfr1(72)=refi3
437  bfr1(73)=rmin3
438  bfr1(74)=rmax3
439 260 continue
440  bfr1(76)=itrans !flag for computing diffused transmiss
441 c
442 c*****format statements******************************************************
443  100 format(2i4,1x,1p6e11.3)
444  125 format(a)
445  126 format(6d12.4,i4)
446  128 format(4x,6d12.4)
447  145 format(i5)
448  150 format(5x,1p6e12.4)
449  195 format (1x,'xzx',1p10d11.2)
450  196 format (1x,'izx',2i5)
451  197 format (1x,'athcnd',1p10d11.2)
452  198 format (1x,'aphcnd',1p10d11.2)
453  199 format (1x,'thetap',1p10d11.2)
454  200 format (1x,'ftrx',1p10d11.2)
455 c 205 format('xzx1-3,mtha,mphi'/1p3e11.3,0p2i4)
456  201 format (1x,'pti',1p10d11.2)
457  505 format (1x,'end of file encountered , check the program')
458  1090 format(a)
459  1095 format(1p4e11.3,1p2e12.4)
460  1100 format(1p7e11.3)
461  1105 format(1x,'input wavelength and phase matrix wavelength'
462  1 ' do not match.'/1x,'input wavelength=',1pe12.4,
463  2 'phase matrix wavelength=',1pe12.4)
464  1120 format(1p2d18.12,1p2d19.12,0pf6.2)
465  1125 format(i5,11x,i5)
466  1130 format(9f8.5)
467  1160 format(a)
468 c****************************************************************************
469  return
470  end
471 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
subroutine closeunits
Definition: closeunits.f:2
subroutine ccn(ifunc, irgm)
Definition: phs.f:940
subroutine readin(ilm, irh, isd, itau)
Definition: readin.f:2
#define real
Definition: DbAlgOcean.cpp:26
float rc(float x, float y)
#define abs(a)
Definition: misc.h:90
subroutine convtc(num, nchar, loc)
Definition: ocn.f:156