OB.DAAC Logo
NASA Logo
Ocean Color Science Software

ocssw V2022
ocn_misr.f
Go to the documentation of this file.
1 c***********************************************************************
2 c..................program ocn.....................................
3 c
4 c
5 c
6  implicit real*8(a-h,o-z)
7  real*8 ez(4),tha(26),pha(47),thp(25), phcnd(46),thcnd(25)
8  real*8 ftrx(16,25,46),pti(8,25,46),xz(3),athcnd(25),aphcnd(46)
9  real*8 wav(20),sig_hg(10),solz(25)
10  integer*4 izx(2),iwav(20)
11 c
12  data ez/0.5d0,0.5d0,0.0d0,0.0d0/
13 c
14  character*60 dir1
15  character*60 dir2
16  character*80 temp
17  character*2 vc
18  character*2 xw
19  character*2 solzc
20  character*4 xwc
21 c
22 c
23 c print*,'entering the main program'
24 c
25  open(5,file='ocn.dir',status='old')
26  read(5,'(a)')dir1
27  read(5,'(a)')dir2
28  len1=index(dir1,' ')-1
29  len2=index(dir2,' ')-1
30 c
31  pi = dacos(-1.0d0 )
32  con=pi/180.0d0
33  iprob=1
34  opt=0.1546d0
35  xsal=0.0d0
36  open(unit=3,file=dir1(1:len1)//'input.dat',status='old')
37  read(3,23)temp
38  read(3,24)iopt,nww,ibgn,iend,isig,isolz
39  read(3,23)temp
40  read(3,25)dlth,dlph
41  read(3,23)temp
42  read(3,26) xr,xi,v,thcel,phcel
43  read(3,23)temp
44  read(3,15)(wav(i),i=1,nww)
45  read(3,23)temp
46  read(3,16)ntha,(tha(i),i=1,ntha)
47  read(3,23)temp
48  read(3,35)nsig,(sig_hg(i),i=1,nsig)
49 c
50  do i=1,nww
51  iwav(i)=int((wav(i)+0.0001)*1000.0 )
52  enddo
53 c
54  if(isolz .eq.1)then
55  read(3,23)temp
56  read(3,16)nsolz,(solz(i),i=1,nsolz)
57  endif
58 
59 c print*,'read complete'
60 c
61  nloop=1
62  if(isolz .eq. 1)then
63  nloop=2
64  endif
65 c
66 c print*,'nloop',nloop
67 c print*,solz
68 
69  do iloop=1,nloop
70 
71  phip=0.0d0
72  nphi=180.0d0/dlph+2.0d0+0.01d0
73  nthcnd=ntha-1
74  nphcnd=nphi-1
75  nthp=nthcnd
76  call rtang(tha,pha,thcnd,phcnd,thp,dlth,dlph,solz,ntha,nphi,
77  1 nthcnd,nphcnd,nthp,nloop)
78 
79 c print*,'thp',thp
80 
81 c
82  do ii=ibgn,iend
83  do mm=1,nsig
84  if(isig.eq.1)then
85  v=sig_hg(mm)**2/5.34d-3
86  endif
87  if(iopt.eq.1)then
88  call riwat(wav(ii),xsal,xr,xi)
89  endif
90  rewind(3)
91  jsg=sig_hg(mm)*100.0
92  iv=(v+0.005d0)*100.0
93  call convtc(jsg,2,vc)
94  call convtc(ii,2,xw)
95  call convtc(iwav(ii),4,xwc)
96  call convtc(isolz,2,solzc)
97 c
98  open(6,file=
99  1 dir2(1:len2)//'owl'//xwc//'x'//vc//'solz'//solzc//'.dat',
100  2 status='unknown')
101 c
102 c print*,dir2(1:len2)//'ol'//xw//'x'//vc//'solz'//solzc//'.dat'
103 
104  open(7,
105  1 file=
106  2 dir2(1:len2)//'sumrywl'//xwc//'x'//vc//'solz'//solzc//'.dat',
107  3 status='unknown')
108 
109  print*,file
110 
111  open(20,file=
112  1 dir2(1:len2)//'rufrtwl'//xwc//'x'//vc//'solz'//solzc//'.dat',
113  2 access='sequential',form='unformatted',
114  3 status='unknown')
115 c
116 c print*,file
117 
118 c stop
119 c
120  izx(1) = nthcnd
121  izx(2) = nphcnd
122  xz(1) = xr
123  xz(2) = xi
124  xz(3) = v
125  do i=1,nthcnd
126  athcnd(i) = thcnd(i)
127  enddo
128  do i=1,nphcnd
129  aphcnd(i) = phcnd(i)
130  enddo
131 c
132  write (20) xz,izx,athcnd,aphcnd
133 c
134  write(7,510)v
135 c
136  do k=1,nthp,1
137  thetap = thp(k)
138  amuo=dcos(thp(k)*con)
139  call rufocn(ez,thp(k),phip,xr,xi,v,pha,tha,phcnd,
140  1 thcnd,thcel,phcel,opt,ftrx,pti,pi,con,
141  2 ntha,nphi,iprob,nthcnd,nphcnd)
142 c
143  write (20) thetap,ftrx,pti
144 c
145  enddo
146  close(6)
147  close(7)
148  close(20)
149  enddo
150  enddo
151  enddo
152 c***********************************************************************
153 c......format statements................................................
154 c
155 15 format(20e10.3)
156 16 format(i5,26f5.1)
157 23 format(a)
158 24 format(6i8)
159 25 format(2f10.3)
160 26 format(2e10.3,4f10.4)
161 27 format(2f10.1)
162 28 format(12f6.1,8x)
163 29 format(10f11.4)
164 30 format(1x,2i3,1p8d15.4/1x,6x,1p8d15.4/)
165 35 format(i5,10f7.3)
166 500 format(1x,f6.2,1x,f6.2,2x,1p4d11.3,1x,1p3d11.3,1x,1p2d11.3)
167 505 format(1x, ' ')
168 510 format('Cox-Monk Rough Ocean (Wind Direction Independent)'/
169  1 'wind speed (m/sec)'/f10.2)
170 c
171  stop
172  end
173 c***********************************************************************
174 
175  subroutine convtc (num,nchar, loc)
176 c***********************************************************************
177 c library subroutine
178 c
179 c module name: convtc
180 c
181 c version : 1.0
182 c
183 c programmer : lucy liu and michael peng, stx, august 1990
184 c
185 c purpose: to convert first nchar digits of an integer into a character
186 c or a string. (note: nchar doesn't include the sign if num
187 c is negative)
188 c
189 c calling sequence: convtc (num,nchar, loc)
190 c
191 c subroutines called: none
192 c
193 c intrinsic functions used: char, ichar, iabs
194 c
195 c common blocks used: none
196 c
197 c arguments and local variables:
198 c
199 c name type i/o descriptions
200 c --------- ---- --- ----------------------------------------------
201 c num i*4 i integer number for conversion
202 c nchar i*4 i number of chars to be converted
203 c loc c*x o char variable or array (string)
204 c istart i*4 start position of integral char string:
205 c 1 for positive integer; 2 for negative integer
206 c irem i*4 absolute value or remainder by 10
207 c itemp i*4 temporary buffer
208 c***********************************************************************
209 c
210 c --- arguments
211 c
212  integer num, nchar
213  character loc(1)
214 c
215 c --- local variables
216 c
217  integer istart, irem, itemp
218 c
219  irem = iabs(num )
220  istart = 1
221 c
222 c --- keep the sign if the integer is negative
223 c
224  if (num.lt.0) then
225  nchar = nchar + 1
226  loc(1) = '-'
227  istart = istart + 1
228  endif
229 c
230 c --- convert digit by digit starting with the least significant digit
231 c
232  do 100 i = nchar, istart, -1
233  itemp = irem / 10
234  loc(i) = char(irem - itemp * 10 + ichar('0'))
235  irem = itemp
236  100 continue
237 c
238 c --- put an asterisk at the first position if the actual number of
239 c --- digits in the integer is found greater than nchar
240 c
241  if (irem.ne.0) loc(1) = '*'
242 c
243  return
244  end
245 c*****************************************************************************
246 c
247  subroutine probl(zx,zy,v,prob,pi,sig,iprob)
248  implicit real *8(a-h), real*8(o-z)
249 c*****************************************************************************
250 c
251  csum=0.0d0
252  aaa=1.0d0
253  go to(100,200),iprob
254 100 continue
255  sigsqx = 0.003d0 + 5.12d-3 *v
256 c sigsq_hg=5.34d-3*v !use HG relationship
257 c sigsqx=sigsq_hg
258  sigsqy = sigsqx
259  sig=dsqrt(sigsqx)
260  go to 250
261 200 continue
262  aaa=2.0d0
263  sigsqx = 0.003d0 + 1.92d-3*v
264  sigsqy = 3.16d-3*v
265 c
266 250 continue
267  sigx = dsqrt(sigsqx)
268  sigy = dsqrt(sigsqy)
269  if(dabs(zx ) .lt. 1.0d-06) zx =0.0d0
270  if(dabs(zy ) .lt. 1.0d-06) zy =0.0d0
271  zeta = zx/sigx
272  eata = zy/sigy
273 c
274 c determine the limiting value of slope
275 c that would yield prob1=1.0d-6
276  pqr= -dlog( aaa*pi*sigx*sigy*1.0d-6)
277  eee=(zeta**2 + eata**2 )/aaa
278  if(eee .gt. pqr) go to 81
279 c
280 71 continue
281  if(iprob .eq.1) go to 300
282 c
283 c following paramters are cox-munk's wind dependent coeff.
284  c21 = 0.01d0 - 0.0086d0 * v
285  c03 = 0.04d0 - 0.033d0 * v
286  c40 = 0.40d0
287  c22 = 0.12d0
288  c04 = 0.23d0
289 c
290  c1 = c21*(zeta**2-1.0d0)*eata *0.5d0
291  c2 = (c03*(eata**3 -3.0d0*eata))/6.0d0
292  c3 = (c40*( zeta**4 - 6.0d0*zeta**2 + 3.0d0))/24.0d0
293  c4 = ( c22*(zeta**2-1.0d0)*(eata**2-1.0d0))/4.0d0
294  c5 = (c04 * (eata**4 - 6.0d0*eata**2 + 3.0d0) )/24.0d0
295 c
296 72 continue
297  csum = -c1-c2+c3+c4+c5
298 300 continue
299  prob1 = dexp(-(zeta**2 + eata**2)/aaa)/( aaa*pi*sigx*sigy)
300  prob = prob1*(1.0d0+csum)
301 c
302 81 continue
303 c
304  return
305  end
306 c*********************************************************************
307 c
308  subroutine rufocn(ei,thpd,phpd,xr,xi,v,pha,tha,phcnd,thcnd,
309  1 thcel,phcel,opt, ftrx,pti,pi,con,ntha,nphi,iprob,
310  2 nthcnd,nphcnd)
311 c*********************************************************************
312  implicit real*8(a-h), real*8(o-z)
313 c
314  real*8 ftrx(16,25,46),pti(8,25,46),a(6)
315  real*8 ei(4),zint(4,25,46),zmio(4),tha(26),pha(47),thx(15),
316  1 phx(35),rint(4,25,46),phsl(3),thsl(3),sang(25,46),snside(5),
317  2 phnl(3,3),thnl(3,3),coside(5),side(5),cosang(6),angtri(6),
318  3 tmio(4),xint(4),tzm(4,4),tint(4,4),phcnd(46),apti(4),thcnd(25),
319  4 frac(10),axz(16),amul(2),amulsq(2),xi_hgcnt(25,46)
320 c*********************************************************************
321  data frac/0.0d0,0.1d0,0.2d0,0.3d0,0.4d0,0.5d0,0.6d0,
322  10.7d0,0.8d0,0.9d0/
323 c*********************************************************************
324 c
325 c write(6,*)'entering rufocn sub program'
326  csum=0.0d0
327  nmtha=ntha-1
328  nmphi=nphi-1
329  ithpd = thpd
330  thp =thpd*con
331  php =phpd*con
332  costhp=dcos(thp)
333  sinthp=dsqrt(1.0d0-costhp**2)
334  domgan=0.0d0
335 c
336 c computation of rough ocean begins
337  nn = 1 !no lambertian fraction
338  pfluxs=0.0d0
339  pfluxz = 0.0d0
340 c
341  do 50 ii=1,nmtha,1 !begin loop over theta (big cell)
342  thcntd=thcnd(ii)
343  thcnt = thcntd*con
344  cthcnt = dcos(thcnt)
345  sthcnt=dsqrt(1.0d0-cthcnt**2)
346 c
347 c divide the big cell into small cells of delta theta (thcel deg.)
348 c and delta phi (phcel deg.)
349 c star constructing small cells in theta space
350  i1=1
351  thx(i1)=tha(ii)
352 c
353 102 continue
354  i1 = i1 + 1
355  thx(i1)=thx(i1-1)+thcel
356  dif1 = dabs(tha(ii+1) - thx(i1) )
357  if(dif1 .lt. 1.0d-6 .or. thx(i1) .gt. tha(ii+1) ) go to 101
358  go to 102
359 c
360 101 continue
361  thx(i1) = tha(ii+1)
362  i1m1=i1-1
363 c
364  do 60 jj=1,nmphi !begin loop over phi (big cell)
365  phcntd= phcnd(jj)
366  iphcnt=phcntd
367  phcnt = phcntd*con
368  cphcnt=dcos(phcnt)
369  sphcnt=dsqrt(1.0d0-cphcnt**2)
370 c
371 c solid angle of the big cell
372  sang3=(dcos(tha(ii)*con)-dcos(tha(ii+1)*con))*
373  1 (pha(jj+1)-pha(jj))*con
374 c
375 c initialize the buffers apti and tint
376  do ij=1,4
377  apti(ij) = 0.0d0
378  do ik=1,4
379  tint(ij,ik)=0.0d0
380  enddo
381  enddo
382 c
383 c start constructing small cells in phi space
384  i2=1
385  phx(i2)=pha(jj)
386 c
387 103 continue
388  i2 = i2 + 1
389  phx(i2)=phx(i2-1)+phcel
390  dif2 = dabs(pha(jj+1) - phx(i2) )
391  if( dif2 .lt. 1.0d-10 .or. phx(i2) .gt. pha(jj+1) )goto 104
392  go to 103
393 c
394 104 continue
395  phx(i2) = pha(jj+1)
396  i2m1=i2-1
397 c
398 c construct a small cell around the center of the big cell
399  phx(i2+1)=phcntd-phcel/2.0d0
400  phx(i2+2)=phcntd+phcel/2.0d0
401  thx(i1+1)=thcntd-thcel/2.0d0
402  thx(i1+2)=thcntd+thcel/2.0d0
403 c
404 c now loop over two times. first to compute the intensity
405 c values at the center of the big cell to get the point value
406 c for the center of the big cell and then compute the intensity
407 c values at the center of the small cells to get the average
408 c value for the big cell
409  do 250 mmnn=1,2
410  if(mmnn .eq.1) goto 260
411  limth1=1
412  limth2=i1m1
413  limfi1=1
414  limfi2=i2m1
415  goto 270
416 c
417 260 continue
418  limth1=i1+1
419  limth2=(i1+2)-1
420  limfi1=i2+1
421  limfi2=(i2+2)-1
422 270 continue
423 c
424 c find the mid point of the small cell and express the extremities
425 c as well as the mid point in radians
426 c begin the loop over theta of the small cells
427  do 10 i=limth1,limth2
428  thsd = (thx(i)+thx(i+1))/2.0d0
429  if(thx(i) .lt. 0.0d0) thx(i)=0.0d0
430  thsl(1) = thx(i)*con
431  thsl(2) = thx(i+1)*con
432  thsl(3) = thsd*con
433  ths = thsl(3)
434  cosths =dcos(ths)
435  sinths=dsqrt(1.0d0-cosths**2)
436  amul(1) = dcos(thsl(1))
437  amul(2) = dcos(thsl(2))
438  amulsq(1) = amul(1) ** 2
439  amulsq(2) = amul(2) ** 2
440 c
441 c begin the loop over phi of the small cells
442  do 11 j=limfi1,limfi2
443 c first initalialize the buffers xint and rint
444  do im=1,4
445  xint(im)=0.0d0
446  rint(im,i,j)=0.0d0
447  enddo
448  phsd = (phx(j)+phx(j+1))/2.0d0
449  if(phx(j+1) .gt. 180.0d0) phx(j+1)=180.0d0
450  if(phx(j) .lt. 0.0d0)phx(j)=0.0d0
451  phsl(1) = phx(j)*con
452  phsl(2) = phx(j+1)*con
453  phsl(3) = phsd*con
454  phs = phsl(3)
455  delphi = phsl(2)-phsl(1)
456 c initialize the buffer thnl and phnl for the 'normal' vector
457  do j1=1,3
458  do j2=1,3
459  thnl(j1,j2)=0.0d0
460  phnl(j1,j2)=0.0d0
461  enddo
462  enddo
463  if(thpd .gt. 1.0d-3) go to 74
464 c
465 73 continue
466  thnl1=thsl(1)/2.0d0
467  thnl2=thsl(2)/2.0d0
468  domgan=(dcos(thnl1)-dcos(thnl2) )*(phsl(2)-phsl(1))
469  sang(i,j)=(dcos(thsl(1))-dcos(thsl(2)))*
470  1 (phsl(2)-phsl(1) )
471  thnl(3,3)=thsl(3)/2.0d0
472  phnl(3,3)=phsl(3)
473  go to 70
474 74 continue
475  if(thsd .gt.1.0d-3) go to 75
476  call solidz(thp,php,thsl,phsl,domgan,sang(i,j),con)
477  thnl(3,3)=thp/2.0d0
478  phnl(3,3)=php
479  go to 70
480 c
481 75 continue
482 c
483  iphsd = phsd+1.0d-5
484 c
485 c calculate the solid angles of the small cell
486  sang(i,j)=(dcos(thx(i)*con)-dcos(thx(i+1)*con))*
487  1 (phx(j+1)-phx(j))*con
488 c
489 c find the normals that would reflect the ray to the
490 c four corners as well as the mid point of the small cell
491  do 31 ix=1,3
492  do 32 iy=1,3
493  thnl(ix,iy)=0.0d0
494  phnl(ix,iy)=0.0d0
495  if(iy .eq. 3 .and. ix .ne. 3) go to 32
496  if(ix .eq. 3 .and. iy .ne. 3) go to 32
497  cthsl=dcos(thsl(ix))
498  sthsl=dsqrt(1.0d0-cthsl**2)
499 c compute the scattering angle assuming the direction
500 c of the photon as phi=0
501 c costwo=costhp*cthsl+sinthp*sthsl*dcos(php-phsl(iy))
502  cos_scat=-costhp*cthsl+sinthp*sthsl*dcos(php-phsl(iy))
503  costwo=-cos_scat
504  if(costwo.ge.1.0d0)costwo=1.0d0
505  if(costwo.le. -1.0d0)costwo=-1.0d0
506  sintwo=dsqrt(1.0d0-costwo**2)
507  scath = 0.5d0*dacos(costwo)
508  cscath=dcos(scath)
509  if(cscath.ge.1.0d0)cscath=1.0d0
510  if(cscath.le. -1.0d0)cscath=-1.0d0
511  sscath=dsqrt(1.0d0-cscath**2)
512  cosai1 = 1.0d0
513  diff =dabs(php-phsl(iy))
514  cons = dabs(dsin(diff) )
515  if(sintwo.lt.1.0d-10)goto 37
516  cosai1=(cthsl-costhp*costwo)/(sintwo*sinthp)
517 37 continue
518  if(cosai1 .gt. 1.0d0 ) cosai1=1.0d0
519  if(cosai1 .lt. -1.00001d0) write(6,143) ix,iy,cosai1
520  if(cosai1 .lt. -1.0d0 ) cosai1=-1.0d0
521  ai1 =dacos(cosai1)
522  costhn=costhp*cscath+sinthp*sscath*cosai1
523  if(costhn.ge.1.0d0)costhn=1.0d0
524  if(costhn.le. -1.0d0)costhn=-1.0d0
525  sinthn=dsqrt(1.0d0-costhn**2)
526  thnl(ix,iy) = dacos(costhn)
527 c
528  if(thnl(ix,iy) .lt.1.0d-10 .and. scath .le. thp)
529  1 cosxx=1.0d0
530  if(thnl(ix,iy) .lt.1.0d-10.and. scath .gt.thp)
531  1 cosxx=-1.0d0
532  if(thnl(ix,iy).lt. 1.0d-10 ) go to 150
533  cosxx=1.0d0
534 c
535  if(sinthn.le.0.0d0)goto 150
536  cosxx=(cscath-costhp*costhn)/( sinthp*sinthn )
537 150 continue
538  if(cosxx .gt. 1.0d0) cosxx=1.0d0
539  if(cosxx .lt.-1.0d0) cosxx = -1.0d0
540  xx=dacos(cosxx)
541  phnl(ix,iy) = xx-php
542 32 continue
543 31 continue
544 c compute the solid angle of the normal cone
545  call znorml(thnl,phnl,domgan )
546  if(iphsd .eq. 0 .or. iphsd .eq. 180) domgan=2.0d0*domgan
547  if(iphsd .eq. 0 .or. iphsd .eq.180)sang(i,j)=
548  1 2.0d0*sang(i,j)
549 70 continue
550 c compute the x and y component of the slope and the dzxdzy
551  thn = thnl(3,3)
552  phn = phnl(3,3)
553  cosphn=dcos(phn)
554  sinphn=dsqrt(1.0d0-cosphn**2)
555  thnd=thn/con
556  phnd=phn/con
557  cosnnn =dcos(thn)
558  zx =dtan(thn)*sinphn
559  zy =dtan(thn)*cosphn
560  dzxdzy=domgan/cosnnn**3
561 c compute the slope probabilty and the shadowing factors
562  prob=0.0d0
563  probx=0.0d0
564  proby=0.0d0
565  xi_hg_corr=0.0d0
566  call probl(zx,zy,v,probx,pi,sig,iprob)
567  call ufunc(pi,cosths,sig,uuff_ths)
568  call ufunc(pi,costhp,sig,uuff_thp)
569  prob=probx*uuff_ths*uuff_thp
570  if(mmnn.eq.1)then
571  call probl_hg(pi,costhp,php,cthcnt,phcnt,sig,
572  1 alfa_hg,proby,xi_hg)
573  call ufunc(pi,cthcnt,sig,uuff_thcnt)
574  call ufunc(pi,costhp,sig,uuff_thp)
575  xi_hg_corr=xi_hg*uuff_thcnt*uuff_thp
576  xi_hgcnt(ii,jj)=xi_hg_corr
577  endif
578 c
579 c initialize the buff tzm
580  do ij=1,4
581  do ik=1,4
582  tzm(ij,ik)=0.0d0
583  enddo
584  enddo
585 c
586 c compute the elements of the reflection matrix
587  call zangl(thp,php,ths,phs,thn,phn,xr,xi,zmio,
588  1 ei,cost,thcnt,phcnt,tzm)
589 c
590 c determine the lambertian and non-lambertian fractions
591  r1=frac(nn)
592  r2=(1.0d0-r1)
593 c
594 c calculate the intensity at the center of the small cell
595  do 15 im=1,4
596  rint(im,i,j) =(cost/(cosnnn*cosths))*
597  1 prob*dzxdzy*zmio(im)*r2
598  if(im.eq.1 .or. im.eq.2)rint(im,i,j)=rint(im,i,j)+
599  1 r1*costhp*sang(i,j)/(2.0d0*pi)
600  xint(im)=rint(im,i,j)/sang(i,j)
601  if (ithpd.eq.0) go to 15
602  if(iphsd.eq.0 .or. iphsd.eq.180)
603  1 rint(im,i,j)=0.5d0*rint(im,i,j)
604 15 continue
605  xin=xint(1)+xint(2)
606  if(mmnn .eq. 1) go to 26
607  pfluxs=pfluxs+(rint(1,i,j)+rint(2,i,j) )*cosths*2.0d0
608  pfluxz=pfluxz+(amulsq(1)-amulsq(2))*delphi*xin
609 c
610  cccc = ((cost*prob*dzxdzy)/(cosnnn*cosths*sang3))
611  do ij=1,4
612  do ik=1,4
613  tint(ij,ik)=tint(ij,ik)+tzm(ij,ik)*cccc
614  enddo
615  enddo
616  go to 27
617 26 continue
618 c save the intensity value in apti buffer
619  do ij=1,4
620  apti(ij)=xint(ij)
621  enddo
622 27 continue
623 c
624 11 continue
625 c end of summation over phi interval of the big cell
626 10 continue
627 c end of summation over theta interval of the big cell
628 250 continue
629 c end of the loop for point value and average value calculations
630 c
631 c calculate the average intensity over the big cell
632 c intialize the buffer zint
633  do im=1,4
634  zint(im,ii,jj)=0.0d0
635  enddo
636  do im=1,4
637  do i=1,i1m1
638  do j=1,i2m1
639  zint(im,ii,jj)=zint(im,ii,jj)+rint(im,i,j)/sang3
640  enddo
641  enddo
642  enddo
643 c
644 c THE PROGRAM HAS BEEN MODIFIED SO THAT THE SUN VECTOR IS
645 C ALONG THE DIRECTION OF THE PHOTON. NO NEED TO SWAP ZERO DEGREE PHI WITH 180 DEG. PHI
646 C ll =nmphi-jj+1
647  ll=jj
648  do ij = 1,4
649  pti(ij,ii,ll)=apti(ij) !point value at the center of the big cell
650  pti(ij+4,ii,ll) = zint(ij,ii,jj) ! average value
651  do ik = 1,4
652  kk = (ij-1)*4 + ik
653  ftrx(kk,ii,ll) = tint(ij,ik)*r2
654  if(kk.eq.1 .or.kk.eq.6)ftrx(kk,ii,ll)=ftrx(kk,ii,ll)+
655  1 (costhp*r1)/pi
656  enddo
657  enddo
658 
659  a(1)=0.5d0*(ftrx(1,ii,ll)+ftrx(2,ii,ll))
660  a(2)=0.5d0*(ftrx(5,ii,ll)+ftrx(6,ii,ll))
661  a(3)=0.5d0*(ftrx(9,ii,ll)+ftrx(10,ii,ll))
662  a(4)=0.5d0*(ftrx(13,ii,ll)+ftrx(14,ii,ll))
663  a(5)=a(1)+a(2)
664  a(6)=a(5)*pi
665  pti(5,ii,ll)=a(1)
666  pti(6,ii,ll)=a(2)
667  pti(7,ii,ll)=a(3)
668  pti(8,ii,ll)=a(4)
669  refl=pfluxs/dcos(thp)
670  reflz = pfluxz / dcos(thp)
671  azn=zint(1,ii,jj)+zint(2,ii,jj)
672  bzn=azn*pi
673 c write(6,701)ii,jj,nmtha,(zint(im,ii,jj),im=1,4),azn,bzn,sang3
674 701 format('ii,jj,nmtha,zint,azn,bzn,sang3',3i4/1p7e11.3)
675 60 continue
676 c write(6,146)
677  if(thcntd .gt. 1.0d-1) go to 50
678  do 62 ib=1,16
679  axz(ib)= (ftrx(ib,ii,1)+ftrx(ib,ii,nmphi))/2.0
680 62 continue
681  do 64 ia=1,nmphi
682  do 63 ib=1,16
683  ftrx(ib,ii,ia)=axz(ib)
684 63 continue
685  pti(5,ii,ia)=( ftrx(1,ii,ia)+ftrx(2,ii,ia) )/2.0
686  pti(6,ii,ia)=( ftrx(5,ii,ia)+ftrx(6,ii,ia) )/2.0
687  pti(7,ii,ia)=( ftrx(9,ii,ia)+ftrx(10,ii,ia))/2.0
688  pti(8,ii,ia)=( ftrx(13,ii,ia)+ftrx(14,ii,ia))/2.0
689 64 continue
690 50 continue
691  rr1=r1*100.0d0
692  rr2=r2*100.0d0
693 c write(6,149) thpd,refl
694  write(7,515)
695  write(7,151)thpd,refl
696 c
697 c call rufint(pi,con,amuo,thcnd,phcnd,opt,pti,
698 c 1 xi_hgcnt,nthcnd,nphcnd)
699 c
700  call rufint_new(pi,con,costhp,thcnd,phcnd,opt,pti,
701  1 xi_hgcnt,nthcnd,nphcnd)
702 c55 continue
703 c***********************************************************************
704 c......format statements................................................
705 c
706 100 format(1x,8i5)
707 136 format( /1x,1p11d10.2/ )
708 139 format(1x,' the avrg.value of stokes parameter for big cell no.',
709  1 5x,i2,',',i2, 'are')
710 143 format(1x,' the cosine of',2i2, ' angle is large',1pd15.5)
711 145 format(t2,'ang',t11,'mu',t17,'phi',t25,'i sub l',t38,'i sub r',
712  1 t52,' u',t64,' v',t76,'i',t88,'p',t97,' pi*i ',t110,'prob' /)
713 146 format(1x ' ' /)
714 147 format(1x,20x,1p4d12.4)
715 148 format(1x,f5.2,1x,f6.2,1x,1p4d11.3,1p5d12.4 )
716 149 format(1x,' sea surface reflectivity for solar'
717  1 ,' zenith angle = ',0pf8.2,' is =',1p2d12.3/)
718 151 format(f10.1,2f10.3)
719 155 format(1x,'fraction of incident flux reflected according to
720  1lambert law (in percent)',0pf7.1/
721  2 1x,'fraction of incident flux reflected according to fresnel
722  3 law (inpercent)',0pf7.1//)
723 180 format(1x,1p8e14.4/1x,1p8e14.4//)
724 200 format(1x,1p10e12.4)
725 500 format(1x,0pf7.2,0pf8.4,0pf7.1,0p2f7.2,1p9e11.4)
726 515 format(t7,'sza',t15,'h. ref1')
727 600 format(1x,2i5,1p4e15.4)
728 c************************************************************************
729  return
730  end
731 c************************************************************************
732 c
733  subroutine rtang(tha,pha,thcnd,phcnd,thp,dlth,dlph,solz,ntha,nphi,
734  1 nthcnd,nphcnd,nthp,nloop)
735 c
736  implicit real*8 (a-h,o-z)
737  real*8 tha(26),pha(47),thp(25),phcnd(46),thcnd(25),solz(25)
738 c*****************************************************************************
739 c
740  do i=1,nthcnd
741  thcnd(i)=(tha(i+1)+tha(i))/2.0d0
742  thp(i)=thcnd(i)
743  if(nloop .eq.2)then
744  thp(i)=solz(i)
745  endif
746  enddo
747 c
748  pha(1)=0.0d0
749  pha(2)=pha(1)+dlph/2.0d0
750  phcnd(1)=0.0d0
751  do i=2,nphcnd
752  pha(i+1)=pha(i)+dlph
753  phcnd(i)=(pha(i+1)+pha(i))/2.0d0
754  if(pha(i+1).gt.180.0d0)pha(i+1)=180.0d0
755  if(phcnd(i).ge.180.0d0)phcnd(i)=180.0d0
756  enddo
757 c
758  return
759  end
760 c************************************************************************
761 c
762  subroutine rufint(pi,con,amuo,thcnd,phcnd,opt,pti,
763  1 xi_hgcnt,nthcnd,nphcnd)
764 c************************************************************************
765 c
766  implicit real * 8 (a-h,o-z)
767  real*8 thcnd(25),phcnd(46),xi_hgcnt(25,46)
768  real * 8 pti(8,25,46)
769 c************************************************************************
770 c
771  write(7,*)'welcome to rufint'
772  do 30 kk=1,1
773  if (kk .eq. 1) n=4
774  if (kk .eq. 2) n=8
775  l = (kk-1) * 4 + 1
776  do 20 i=1,nthcnd
777  do 10 j=1,nphcnd
778  jr=nphcnd-j+1
779  amu = dcos(thcnd(i) * con)
780  pol = 0.0d0
781  xin = pti(l,i,j) + pti(l+1,i,j)
782  xpt = opt
783  xinpi = xin * pi
784  xinsea = xinpi * dexp(-opt/amuo)
785  xintop = xinsea * dexp(-xpt/amu)
786  if (xin .gt. 1.0d-10) poln=sqrt((pti(l,i,j) - pti(l+1,i,j))
787  1 **2 + pti(l+2,i,j)**2 + pti(l+3,i,j)**2 )
788  if (xin .gt. 1.0d-10) pol = poln/xin
789  if (pti(l,i,j) .gt. pti(l+1,i,j) ) pol = -pol
790  write (7,500) thcnd(i),phcnd(j),(pti(m,i,j),m=l,n),xin,pol,
791  1 xinpi,xi_hgcnt(i,jr)*pi,xintop
792  10 continue
793  write (7,505)
794  20 continue
795  30 continue
796 c***********************************************************************
797 c..................format statments.....................................
798 c
799  500 format(1x,f6.2,1x,f6.2,2x,1p4d11.3,1x,1p3d11.3,1x,1p3d11.3)
800  505 format(1x,' ')
801 c.......................................................................
802  return
803  end
804 c******************************************************************************
805 c
806  subroutine solidz(thp,php,thsl,phsl,domgan,domega,con)
807  implicit real * 8(a-h,o-z)
808  real*8 thsl(3),phsl(3)
809 c******************************************************************************
810 c
811  costhp=dcos(thp)
812  sinthp=dsqrt(1.0d0-costhp**2)
813  costh2=dcos(thsl(2) )
814  sn=thp/2.0d0
815  sn1=(thp+thsl(2) )/2.0d0
816  costci=costhp*costh2
817  sintci=dsqrt(1.0d0-costci**2)
818  chi=0.50d0*dacos(costci)
819  chid=chi/con
820 c write(6,126)costhp,sinthp,costh2,sn,sn1,costci,
821 c 1 sintci,chi,chid
822  coschi=dcos(chi)
823  sinchi=dsqrt(1.0d0-coschi**2)
824  cosai=(costh2-costhp*costci)/(sinthp*sintci)
825  ai=dacos(cosai)
826  cossn=dcos(sn)
827  sinsn=dsqrt(1.0d0-cossn**2)
828  cossn1=dcos(sn1)
829  sinsn1=dsqrt(1.0d0-cossn1**2)
830  coss1=cossn*coschi+sinsn*sinchi*cosai
831  coss2=cossn1* coschi+sinsn1*sinchi*cosai
832  s1=dacos(coss1)
833  s2=dacos(coss2)
834  s3=sn1-sn
835  ss=(s1+s2+s3)/2.0d0
836  s1d=s1/con
837  s2d=s2/con
838  s3d=s3/con
839  ssd=ss/con
840 c write(6,126)ai,coss1,coss2,s1,s2,s3,ss,s1d,s2d,
841 c 1 s3d
842  st=dtan(ss/2.0d0)*dtan((ss-s1)/2.0d0)*
843  1 dtan((ss-s2)/2.0d0)*dtan((ss-s3)/2.0d0)
844  se=dsqrt(1.0d0/(1.0d0+st) )
845  domgan=4.0d0*dacos(se)
846 c write(6,126)st,se,domgan
847  r1=thsl(2)-thsl(3)
848  r3=dacos(costh2*costh2)
849  rr=r1+r3/2.0d0
850  sr=dtan(rr/2.0d0)*dtan((rr-r1)/2.0d0)**2 *
851  1 dtan((rr-r3)/2.0d0)
852  re=dsqrt(1.0d0/(1.0d0+sr) )
853  domega=4.0d0*dacos(re)
854  r1d=r1/con
855  r3d=r3/con
856 c write(6,126)r1,r3,rr,sr,re,domega,r1d,r3d
857 126 format(1x,1p10d12.4/)
858 
859  return
860  end
861 c******************************************************************************
862 c
863  subroutine zangl(thp,php,ths,phs,thn,phn,xr,xi,tmio,ei,cost,
864  1 thcnt,phcnt,tzm )
865 c
866  implicit real*8(a-h), real*8(o-z)
867  real*8 zmio(4),ei(4),tmio(4),z(4,4),t(4,4),tzm(4,4)
868 c******************************************************************************
869 c
870 c calculate the scattering angle (in the scattering plane)
871 c write(6,*)'entering sub zangl'
872  xzx=-1.0d-10
873  cost = dcos(thn)*dcos(thp) + dsin(thn)*dsin(thp)*dcos(php-phn)
874  thetac= dacos(cost)
875 c calculate the elements of fresnel matrix
876 c write(6,*)'begin fresnel matrix'
877  sinsq = (1.0d0-cost**2)
878  sinsq2 = sinsq**2
879  a = (xr**2-xi**2 -1.0d0 + cost**2)
880  b = (-2.0d0 * xr * xi)
881  r = dsqrt(a**2 + b**2)
882  tmr = cost**2*r
883  xp = dsqrt(2.0d0*r + 2.0d0*a)
884  xmm = (2.0d0*r - 2.0d0*a)
885  if(xmm .lt. 0.0d0 .and. xmm .gt. xzx) xmm=-xmm
886  xm=dsqrt(xmm)
887  dnmr1 = (sinsq2 + tmr + cost*sinsq*xp)
888  qrmu = (sinsq2 - tmr)/dnmr1
889  qimu = (cost*sinsq*xm)/dnmr1
890  dnmr2 = (cost**2 + r + cost*xp)
891  rrr = (cost**2 - r)/dnmr2
892  rri = (cost*xm)/dnmr2
893  rer = qrmu*rrr - qimu*rri
894  rei = qimu*rrr + qrmu*rri
895  r11 = rer**2 + rei**2
896  r22 = rrr**2 + rri**2
897  r33 = rer*rrr + rei*rri
898  r34 = rri*rer - rei*rrr
899 c
900 c
901 c write(6,*)'begin calculation of i1 and i2'
902 c calculate the angle between the meridian and scattering plane
903  costot = dcos(2.0d0*thetac)
904  sintot = dsin(2.0d0*thetac)
905  cosai1 = 1.0d0
906  cosai2 = 1.0d0
907  if(thp .le.1.0d-6 .or. ths .le. 1.0d-6) go to 65
908  cons = dsin(php-phs)
909  cons = dabs(cons)
910  if(cons .lt.1.0d-10) go to 65
911  cosai1 = (dcos(ths) - dcos(thp)*costot)/(dsin(thp)*sintot)
912  cosai2 = (dcos(thp) - dcos(ths)*costot)/(dsin(ths)*sintot)
913 65 continue
914  ai1 = dacos(cosai1)
915  ai2 = dacos(cosai2)
916 c calculate the elements of matrix zm (= l(pi-ai2)r(thetac)l(-ai1))
917 c first the matrix element of rl=r(thetac)*(-ai1)
918 c write(6,*)'rotation through -ai1'
919  ssqai1 = dsin(ai1)**2
920  csqai1 = dcos(ai1)**2
921  stoai1 = dsin(2.0d0*ai1)
922  ctoai1 = dcos(2.0d0*ai1)
923  ssqai2 = dsin(ai2)**2
924  csqai2 = dcos(ai2)**2
925  stoai2 = dsin(2.0d0*ai2)
926  ctoai2 = dcos(2.0d0*ai2)
927  rl11 = r11 * csqai1
928  rl21 = r22 * ssqai1
929  rl31 = r33 * stoai1
930  rl41 =-r34 * stoai1
931  rl12 = r11 * ssqai1
932  rl22 = r22 * csqai1
933  rl32 =-r33 * stoai1
934  rl42 = r34 * stoai1
935  rl13 =-0.5d0*r11*stoai1
936  rl23 = 0.5d0*r22*stoai1
937  rl33 = r33 * ctoai1
938  rl43 =-r34 * ctoai1
939  rl14 = 0.d0
940  rl24 = 0.d0
941  rl34 = r34
942  rl44 = r33
943 c write(6,*)'rotation through ai2'
944 c matrix elements of zm = l(pi-ai2)*rl(4*4)
945  z(1,1)= csqai2*rl11 + ssqai2*rl21 - 0.5d0*stoai2*rl31
946  z(2,1)= ssqai2*rl11 + csqai2*rl21 + 0.5d0*stoai2*rl31
947  z(3,1) = stoai2*rl11- stoai2*rl21+ ctoai2*rl31
948  z(4,1) = rl41
949  z(1,2)= csqai2*rl12 + ssqai2*rl22 - 0.5d0*stoai2*rl32
950  z(2,2)= ssqai2*rl12 + csqai2*rl22 + 0.5d0*stoai2*rl32
951  z(3,2) = stoai2*rl12- stoai2*rl22+ ctoai2*rl32
952  z(4,2) = rl42
953  z(1,3)= csqai2*rl13 + ssqai2*rl23 - 0.5d0*stoai2*rl33
954  z(2,3)= ssqai2*rl13 + csqai2*rl23 + 0.5d0*stoai2*rl33
955  z(3,3) = stoai2*rl13- stoai2*rl23+ ctoai2*rl33
956  z(4,3) = rl43
957  z(1,4)= csqai2*rl14 + ssqai2*rl24 - 0.5d0*stoai2*rl34
958  z(2,4)= ssqai2*rl14 + csqai2*rl24 + 0.5d0*stoai2*rl34
959  z(3,4) = stoai2*rl14- stoai2*rl24+ ctoai2*rl34
960  z(4,4) = rl44
961 c transform the stokes parameter to the center of the big cell
962 c write(6,*)'transform to the center of the big cell'
963 c
964  costi1 = 1.0d0
965  costi2 = 1.0d0
966  if(ths .le.1.0d-6 .or. thcnt .le. 1.0d-6) go to 75
967  cons = dsin(phcnt-phs)
968  cons = dabs(cons)
969  if(cons .lt.1.0d-10) go to 75
970  thtx = dcos(thcnt)*dcos(ths) + dsin(thcnt)*dsin(ths)*
971  1 dcos(phcnt-phs)
972  tht = dacos(thtx)
973  thty = dsin(tht)
974  costi1 = (dcos(thcnt)-dcos(ths)*thtx)/(dsin(ths)*thty)
975  costi2 = (dcos(ths) - dcos(thcnt)*thtx) / (dsin(thcnt)*thty)
976 75 continue
977 c write(6,*)'computation of tht done'
978  ti1 = dacos(costi1)
979  ti2 = dacos(costi2)
980  ssqti1= dsin(ti1)**2
981  csqti1= dcos(ti1)**2
982  stoti1= dsin(2.0d0*ti1)
983  ctoti1= dcos(2.0d0*ti1)
984  ssqti2= dsin(ti2)**2
985  csqti2= dcos(ti2)**2
986  stoti2= dsin(2.0d0*ti2)
987  ctoti2= dcos(2.0d0*ti2)
988 c write(6,*)'compute the element of t matrix'
989  t(1,1)= (csqti2*csqti1 +ssqti2*ssqti1 - 0.5d0*stoti2*stoti1)
990  t(2,1)= (ssqti2*csqti1 +csqti2*ssqti1 + 0.5d0*stoti2*stoti1)
991  t(3,1)=stoti2*csqti1 - stoti2*ssqti1 + ctoti2*stoti1
992  t(4,1)= 0.0d0
993  t(1,2)= (csqti2*ssqti1 +ssqti2*csqti1 + 0.5d0*stoti2*stoti1)
994  t(2,2)= (ssqti2*ssqti1 +csqti2*csqti1 - 0.5d0*stoti2*stoti1)
995  t(3,2)=stoti2*ssqti1 - stoti2*csqti1 - ctoti2*stoti1
996  t(4,2)=0.0d0
997  t(1,3)= (-0.5d0*csqti2*stoti1+0.5d0*ssqti2*stoti1-
998  1 0.5d0*stoti2*ctoti1)
999  t(2,3)= (-0.5d0*ssqti2*stoti1+0.5d0*csqti2*stoti1+
1000  1 0.5d0*stoti2*ctoti1)
1001  t(3,3)= (-0.5d0*stoti2*stoti1-0.5d0*stoti2*stoti1+
1002  1 ctoti2*ctoti1)
1003  t(4,3)= 0.0d0
1004  t(1,4)= 0.0d0
1005  t(2,4)= 0.0d0
1006  t(3,4)= 0.0d0
1007  t(4,4)= 0.0d0
1008  do 10 i=1,4
1009  do 10 j=1,4
1010  do 10 k=1,4
1011  tzm(i,j) = tzm(i,j) + t(i,k)*z(k,j)
1012  10 continue
1013  do 20 i=1,4
1014  tmio(i) = 0.0d0
1015  do 20 j=1,4
1016  tmio(i) = tmio(i) + tzm(i,j)*ei(j)
1017 20 continue
1018 c write(6,*)'leaving sub..zangle'
1019 c
1020  return
1021  end
1022 c******************************************************************************
1023 c
1024  subroutine znorml(thnl,phnl,domgan)
1025 c
1026  implicit real*8(a-h), real*8(o-z)
1027  real*8 thnl(3,3),phnl(3,3), coside(5),snside(5), cosang(6),
1028  1 angtri(6),side(5)
1029 c
1030 c******************************************************************************
1031 c
1032  pi = dacos(-1.0d0)
1033 c construct two spherical triangles using four normals(corresponding
1034 c to four reflected beams along the four corners of the small cell
1035  coside(1) = dcos(thnl(1,1))*dcos(thnl(2,1))+dsin(thnl(1,1))*
1036  1 dsin(thnl(2,1))*dcos(phnl(1,1) -phnl(2,1))
1037  coside(2) = dcos(thnl(2,1))*dcos(thnl(2,2))+dsin(thnl(2,1))*
1038  1 dsin(thnl(2,2))*dcos(phnl(2,1) -phnl(2,2))
1039  coside(3) = dcos(thnl(2,2))*dcos(thnl(1,2))+dsin(thnl(2,2))*
1040  1 dsin(thnl(1,2))*dcos(phnl(2,2) -phnl(1,2))
1041  coside(4) = dcos(thnl(1,2))*dcos(thnl(1,1))+dsin(thnl(1,2))*
1042  1 dsin(thnl(1,1))*dcos(phnl(1,2) -phnl(1,1))
1043  coside(5) = dcos(thnl(1,1))*dcos(thnl(2,2))+dsin(thnl(1,1))*
1044  1 dsin(thnl(2,2))*dcos(phnl(1,1) -phnl(2,2))
1045  do 33 iz=1,5
1046  if(coside(iz) .gt. 1.000000d0)then
1047  coside(iz)=1.000000d0
1048  endif
1049  if(coside(iz) .lt. -1.000000d0)then
1050  coside(iz)=-1.000000d0
1051  endif
1052  side(iz) = dacos(coside(iz))
1053  snside(iz)=dsin(side(iz))
1054 33 continue
1055  do 36 k=1,6
1056 36 cosang(k) = 1.0d0
1057  if(side(1).lt.1.0d-10 .or. side(2).lt.1.0d-10) cosang(5)=-1.0d0
1058  if(side(1).lt.1.0d-10 .or. side(2).lt.1.0d-10) go to 45
1059  cosang(1) = (coside(1) -coside(2)*coside(5))/
1060  1 (snside(2)*snside(5))
1061  cosang(2) = (coside(2) -coside(1)*coside(5))/
1062  1 (snside(1)*snside(5))
1063  cosang(5) = (coside(5) -coside(1)*coside(2))/
1064  1 (snside(1)*snside(2))
1065 45 continue
1066  if(side(3).lt.1.0d-10 .or. side(4).lt.1.0d-10) cosang(6)=-1.0d0
1067  if(side(3).lt.1.0d-10 .or. side(4).lt.1.0d-10) go to 55
1068  cosang(3) = (coside(3) -coside(4)*coside(5))/
1069  1 (snside(4)*snside(5))
1070  cosang(4) = (coside(4) -coside(3)*coside(5))/
1071  1 (snside(3)*snside(5))
1072  cosang(6) = (coside(5) -coside(3)*coside(4))/
1073  1 (snside(3)*snside(4))
1074 55 continue
1075  aaaaa=-1.0d0
1076  do 34 iz=1,6
1077  difabs=dabs(dabs(cosang(iz) ) -1.0d0 )
1078  if(difabs .lt. 1.0d-10 .and. cosang(iz) .gt. 1.0d0)
1079  1cosang(iz)=1.0d0
1080  if(difabs .lt. 1.0d-10 .and. cosang(iz) .lt.aaaaa )
1081  1cosang(iz)=-1.0d0
1082 34 angtri(iz)=dacos(cosang(iz))
1083 c find the area of the triangles
1084  area1 =(angtri(1)+angtri(2)+angtri(5) - pi)
1085  area2 =(angtri(3)+angtri(4)+angtri(6) - pi)
1086  domgan = area1+area2
1087  if(phnl(1,1).lt.1.0d-6)then
1088 c write(6,*)'coside',coside
1089 c write(6,*)'side',side
1090 c write(6,*)'cosang',cosang
1091 c write(6,*)'angtri',angtri
1092 c write(6,*)'area1,area2',area1,area2
1093  endif
1094  return
1095  end
1096 c**********************************************************************
1097 c
1098  subroutine riwat(wl,xsal,nr,ni)
1099 C
1100  implicit real*8 (a-h,o-z)
1101  real*8 twl(62),tnr(62),tni(62)
1102  real*8 nr,ni,nrc,nic
1103 c******************************************************************************
1104 c
1105 C input parameters: wl=wavelength (in micrometers)
1106 C xsal=salinity (in ppt), if xsal<0 then 34.3ppt
1107 C by default
1108 C output parameters: nr=index of refraction of sea water
1109 C ni=extinction coefficient of sea water
1110 C
1111 C
1112 C Indices of refraction for pure water from Hale and Querry,
1113 C Applied Optique, March 1973, Vol. 12, No. 3, pp. 555-563
1114  data twl/
1115  s 0.250,0.275,0.300,0.325,0.345,0.375,0.400,0.425,0.445,0.475,
1116  s 0.500,0.525,0.550,0.575,0.600,0.625,0.650,0.675,0.700,0.725,
1117  s 0.750,0.775,0.800,0.825,0.850,0.875,0.900,0.925,0.950,0.975,
1118  s 1.000,1.200,1.400,1.600,1.800,2.000,2.200,2.400,2.600,2.650,
1119  s 2.700,2.750,2.800,2.850,2.900,2.950,3.000,3.050,3.100,3.150,
1120  s 3.200,3.250,3.300,3.350,3.400,3.450,3.500,3.600,3.700,3.800,
1121  s 3.900,4.000/
1122  data tnr/
1123  s 1.362,1.354,1.349,1.346,1.343,1.341,1.339,1.338,1.337,1.336,
1124  s 1.335,1.334,1.333,1.333,1.332,1.332,1.331,1.331,1.331,1.330,
1125  s 1.330,1.330,1.329,1.329,1.329,1.328,1.328,1.328,1.327,1.327,
1126  s 1.327,1.324,1.321,1.317,1.312,1.306,1.296,1.279,1.242,1.219,
1127  s 1.188,1.157,1.142,1.149,1.201,1.292,1.371,1.426,1.467,1.483,
1128  s 1.478,1.467,1.450,1.432,1.420,1.410,1.400,1.385,1.374,1.364,
1129  s 1.357,1.351/
1130  data tni/
1131  s 3.35e-08,2.35e-08,1.60e-08,1.08e-08,6.50e-09,
1132  s 3.50e-09,1.86e-09,1.30e-09,1.02e-09,9.35e-10,
1133  s 1.00e-09,1.32e-09,1.96e-09,3.60e-09,1.09e-08,
1134  s 1.39e-08,1.64e-08,2.23e-08,3.35e-08,9.15e-08,
1135  s 1.56e-07,1.48e-07,1.25e-07,1.82e-07,2.93e-07,
1136  s 3.91e-07,4.86e-07,1.06e-06,2.93e-06,3.48e-06,
1137  s 2.89e-06,9.89e-06,1.38e-04,8.55e-05,1.15e-04,
1138  s 1.10e-03,2.89e-04,9.56e-04,3.17e-03,6.70e-03,
1139  s 1.90e-02,5.90e-02,1.15e-01,1.85e-01,2.68e-01,
1140  s 2.98e-01,2.72e-01,2.40e-01,1.92e-01,1.35e-01,
1141  s 9.24e-02,6.10e-02,3.68e-02,2.61e-02,1.95e-02,
1142  s 1.32e-02,9.40e-03,5.15e-03,3.60e-03,3.40e-03,
1143  s 3.80e-03,4.60e-03/
1144  i=2
1145  10 if (wl.lt.twl(i)) goto 20
1146  if (i.lt.62) then
1147  i=i+1
1148  goto 10
1149  endif
1150  20 continue
1151  xwl=twl(i)-twl(i-1)
1152  yr=tnr(i)-tnr(i-1)
1153  yi=tni(i)-tni(i-1)
1154  nr=tnr(i-1)+(wl-twl(i-1))*yr/xwl
1155  ni=tni(i-1)+(wl-twl(i-1))*yi/xwl
1156 c
1157 c Correction to be applied to the index of refraction and to the extinction
1158 c coefficients of the pure water to obtain the ocean water one (see for
1159 c example Friedman). By default, a typical sea water is assumed
1160 c (Salinity=34.3ppt, Chlorinity=19ppt) as reported by Sverdrup.
1161 c In that case there is no correction for the extinction coefficient between
1162 c 0.25 and 4 microns. For the index of refraction, a correction of +0.006
1163 c has to be applied (McLellan). For a chlorinity of 19.0ppt the correction
1164 c is a linear function of the salt concentration. Then, in 6S users are able
1165 c to enter the salt concentration (in ppt).
1166 c REFERENCES:
1167 c Friedman D., Applied Optics, 1969, Vol.8, No.10, pp.2073-2078.
1168 c McLellan H.J., Elements of physical Oceanography, Pergamon Press, Inc.,
1169 c New-York, 1965, p 129.
1170 c Sverdrup H.V. et al., The Oceans (Prentice-Hall, Inc., Englewood Cliffs,
1171 c N.J., 1942, p 173.
1172 
1173  nrc=0.006
1174  nic=0.000
1175  nr=nr+nrc*(xsal/34.3)
1176  ni=ni+nic*(xsal/34.3)
1177  return
1178  end
1179 c***********************************************************************
1180 c
1181  subroutine ufunc(pi,xmu,sig,uu)
1182  implicit real*8 (a-h,o-z)
1183 c******************************************************************************
1184 c
1185  gnu=xmu/(sig*dsqrt(1.0d0-xmu**2))
1186  gnusq=gnu**2
1187  ff1=dexp(-gnusq)/(dsqrt(pi)*gnu)
1188  ff2=(1.0d0-erfz(gnu))
1189  ffnu=0.5d0*(ff1-ff2)
1190  uu=1.0d0/(1.0d0+ffnu)
1191 c
1192 c write(*,*)'gnu,ff1,ff2,ffnu,uu',
1193 c 1 gnu,ff1,ff2,ffnu,uu
1194  return
1195  end
1196 c**********************************************************************
1197 
1198  function erfz(x)
1199  implicit real*8 (a-h,o-z)
1200 c******************************************************************************
1201  if(x.lt.0.d0)then
1202  erfz=-gammp(.5d0,x**2)
1203  else
1204  erfz=gammp(.5d0,x**2)
1205  endif
1206 c
1207  return
1208  end
1209 c***********************************************************************
1210 
1211  function gammp(a,x)
1212  implicit real*8 (a-h,o-z)
1213 c******************************************************************************
1214  if(x.lt.0.d0.or.a.le.0.d0)pause
1215  if(x.lt.a+1.d0)then
1216  call gser(gamser,a,x,gln)
1217  gammp=gamser
1218  else
1219  call gcf(gammcf,a,x,gln)
1220  gammp=1.d0-gammcf
1221  endif
1222 c
1223  return
1224  end
1225 c**********************************************************************
1226 
1227  subroutine gser(gamser,a,x,gln)
1228  implicit real*8 (a-h,o-z)
1229  parameter(itmax=100,eps=3.0d-7)
1230 c******************************************************************************
1231 c
1232  gln=gammln(a)
1233  if(x.le.0.d0)then
1234  if(x.lt.0.d0)pause
1235  gamser=0.d0
1236  return
1237  endif
1238  ap=a
1239  sum=1.d0/a
1240  del=sum
1241  do 11 n=1,itmax
1242  ap=ap+1.d0
1243  del=del*x/ap
1244  sum=sum+del
1245  if(dabs(del).lt.dabs(sum)*eps)go to 1
1246 11 continue
1247  pause 'a too large, itmax too small'
1248 1 gamser=sum*dexp(-x+a*dlog(x)-gln)
1249  return
1250  end
1251 c**********************************************************************
1252 c
1253  subroutine gcf(gammcf,a,x,gln)
1254  implicit real*8 (a-h,o-z)
1255  parameter(itmax=100,eps=3.0d-7)
1256 c******************************************************************************
1257 c
1258  gln=gammln(a)
1259  gold=0.d0
1260  a0=1.d0
1261  a1=x
1262  b0=0.d0
1263  b1=1.d0
1264  fac=1.d0
1265  do 11 n=1,itmax
1266  an=dfloat(n)
1267  ana=an-a
1268  a0=(a1+a0*ana)*fac
1269  b0=(b1+b0*ana)*fac
1270  anf=an*fac
1271  a1=x*a0+anf*a1
1272  b1=x*b0+anf*b1
1273  if(a1.ne.0.d0)then
1274  fac=1.d0/a1
1275  g=b1*fac
1276  if(dabs((g-gold)/g).lt.eps)go to 1
1277  gold=g
1278  endif
1279 11 continue
1280  pause 'a too large, itmax too small'
1281 1 continue
1282  gammcf=dexp(-x+a*dlog(x)-gln)*g
1283  return
1284  end
1285 c**********************************************************************
1286 c
1287  function gammln(xx)
1288 c
1289  implicit real*8 (a-h,o-z)
1290  real*8 cof(6),stp,half,one,fpf,x,tmp,ser
1291 c******************************************************************************
1292  data cof,stp/76.18009173d0,-86.50532033d0,24.01409822d0,
1293  * -1.231739516d0,.120858003d-2,-.536382d-5,2.50662827465d0/
1294  data half,one,fpf/0.5d0,1.0d0,5.5d0/
1295 c******************************************************************************
1296 c
1297  x=xx-one
1298  tmp=x+fpf
1299  tmp=(x+half)*log(tmp)-tmp
1300  ser=one
1301  do 11 j=1,6
1302  x=x+one
1303  ser=ser+cof(j)/x
1304 11 continue
1305  gammln=tmp+dlog(stp*ser)
1306  return
1307  end
1308 c***********************************************************************
1309 c
1310  subroutine probl_hg(pi,xmup,phirpx,xmux,phirx,sig,alfa,prob,
1311  1 xi_hg)
1312 c
1313  implicit real*8 (a-h,o-z)
1314 c******************************************************************************
1315 c compute the slope probability dist using hg function
1316 c
1317 c in rough ocean program phi is defined wrt to the vector
1318 c pointing toward the sun. The phi angle is corrected in the
1319 c end.....za
1320  xn1=1.334d0
1321  angl=dacos(xmux)
1322  xmu=dcos(pi-angl)
1323  phir=pi-phirx
1324 c
1325  cos2chi=-xmup*xmu-dsqrt((1.0d0-xmup**2)*(1.0d0-xmu**2))*
1326  1 dcos(phirp-phir)
1327 c
1328  alfa=(1.0d0+cos2chi)/(xmup-xmu)**2
1329  arg_exp=(1.0d0-2.0d0*alfa)/sig**2
1330  exp_term=dexp(arg_exp)
1331  prob=exp_term/(pi*sig**2)
1332 c
1333  chii=dacos(cos2chi)/2.0d0
1334  sinchii=dsin(chii)
1335  sinchit=sinchii/xn1
1336  chit=dasin(sinchit)
1337  thitp=chii+chit
1338  thitm=chii-chit
1339  sinthip=dsin(thitp)
1340  sinthim=dsin(thitm)
1341  sinthipsq=sinthip**2
1342  sinthimsq=sinthim**2
1343  costhipsq=1.0d0-sinthipsq
1344  costhimsq=1.0d0-sinthimsq
1345  tanthipsq=sinthipsq/costhipsq
1346  tanthimsq=sinthimsq/costhimsq
1347  if(dabs(chii).lt.1.0d-6)then
1348  rho_pls=((1.0d0-xn1)/(1.0d0+xn1))**2
1349  else
1350  rho_pls=0.5d0*(tanthimsq/tanthipsq +
1351  1 sinthimsq/sinthipsq)
1352  endif
1353  xi_hg=(rho_pls*alfa**2*prob)/dabs(xmu)
1354 c
1355 c
1356 c write(6,99)xmup,phirp,xmu,phir,sig
1357 99 format('xmup,phirp,xmu,phir,sig',5f10.4)
1358 c write(6,100)cos2chi,alfa,sig,arg_exp,exp_term,prob
1359 100 format('cos2chi,alfa,sig,arg_exp,exp_term,prob'/1p6e12.4)
1360 c write(6,101)sinthipsq,sinthimsq,costhipsq,costhimsq
1361 101 format('sinthipsq,sinthimsq,costhipsq,costhimsq'/1p4e12.4)
1362 c write(6,102)tanthipsq,tanthimsq,rho_pls,xi_hg
1363 102 format('tanthipsq,tanthimsq,rho_pls,xi_hg'/1p4e12.4)
1364 
1365 c
1366  return
1367  end
1368 c***********************************************************************
1369 c
1370  subroutine rufint_new (pi,con,amuo,thcnd,phcnd,opt,pti,
1371  1 xi_hgcnt,nthcnd,nphcnd)
1372 c
1373  implicit real * 8 (a-h,o-z)
1374  real*8 thcnd(25),phcnd(46),xi_hgcnt(25,46)
1375  real * 8 pti(8,25,46)
1376 c************************************************************************
1377 c
1378  write(7,100)
1379  do 20 i=1,nthcnd
1380  do 10 j=1,nphcnd
1381  jr=nphcnd-j+1
1382  amu = dcos(thcnd(i)*con)
1383  pol = 0.0d0
1384  xin1=pti(1,i,j)+pti(2,i,j)
1385  xin2=pti(5,i,j)+pti(6,i,j)
1386  xpt=opt
1387  xi_hgpi=xi_hgcnt(i,jr)*pi
1388  xin1pi=xin1*pi
1389  xin1sea=xin1pi*dexp(-opt/amuo)
1390  xin1top =xin1sea*dexp(-xpt/amu)
1391  xin2pi=xin2*pi
1392  xin2sea=xin2pi*dexp(-opt/amuo)
1393  xin2top =xin2sea*dexp(-xpt/amu)
1394  write (7,500) thcnd(i),phcnd(j),xin1pi,xi_hgpi,
1395  1 xin2pi,xin2top
1396  10 continue
1397  write (7,505)
1398  20 continue
1399 c***********************************************************************
1400 c..................format statments.....................................
1401 c
1402  100 format(' the phi I_ZA I_HG',
1403  1 ' Avg_I_ZA Avg_I_ZA_Top')
1404  500 format(1x,f6.2,1x,f6.2,2x,1p4d13.4)
1405  505 format(1x,' ')
1406 c
1407  return
1408  end
1409 c************************************************************************
gammln
function gammln(xx)
Definition: ocn.f:1270
index
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
erfz
function erfz(x)
Definition: ocn.f:1181
gser
subroutine gser(gamser, a, x, gln)
Definition: ocn.f:1210
rufint
subroutine rufint(pi, con, amuo, thcnd, phcnd, opt, pti, xi_hgcnt, nthcnd, nphcnd)
Definition: ocn.f:741
rufint_new
subroutine rufint_new(pi, con, amuo, thcnd, phcnd, opt, pti, xi_hgcnt, nthcnd, nphcnd)
Definition: ocn.f:1354
angl
subroutine angl
Definition: angl.f:2
gammp
function gammp(a, x)
Definition: ocn.f:1194
zangl
subroutine zangl(thp, php, ths, phs, thn, phn, xr, xi, tmio, ei, cost, thcnt, phcnt, tzm)
Definition: ocn.f:842
rtang
subroutine rtang(tha, pha, thcnd, phcnd, thp, dlth, dlph, ntha, nphi, nthcnd, nphcnd, nthp)
Definition: ocn.f:715
znorml
subroutine znorml(thnl, phnl, domgan)
Definition: ocn.f:1002
real
#define real
Definition: DbAlgOcean.cpp:26
probl_hg
subroutine probl_hg(pi, xmup, phirpx, xmux, phirx, sig, alfa, prob, xi_hg)
Definition: ocn.f:1294
riwat
subroutine riwat(wl, xsal, nr, ni)
Definition: ocn.f:1076
gcf
subroutine gcf(gammcf, a, x, gln)
Definition: ocn.f:1236
fac
#define fac
Definition: atmocor1_land.c:111
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
solidz
subroutine solidz(thp, php, thsl, phsl, domgan, domega, con)
Definition: ocn.f:784
re
#define re
Definition: l1_czcs_hdf.c:701
diff
subroutine diff(x, conec, n, dconecno, dn, dconecmk, units, u, inno, i, outno, o, input, deriv)
Definition: ffnet.f:205
rufocn
subroutine rufocn(ei, thpd, phpd, xr, xi, v, pha, tha, phcnd, thcnd, thcel, phcel, opt, ftrx, pti, pi, con, ntha, nphi, iprob, nthcnd, nphcnd)
Definition: ocn.f:291
ufunc
subroutine ufunc(pi, xmu, sig, uu)
Definition: ocn.f:1159
pi
#define pi
Definition: vincenty.c:23
ei
float ei(float x)
probl
subroutine probl(zx, zy, v, prob, pi, sig, iprob)
Definition: ocn.f:228
rint
double rint(double)
convtc
subroutine convtc(num, nchar, loc)
Definition: ocn.f:156


ABOUT

Our Team

Contact Us

DATA

Get Started

Find Data

Use Data

TOOLS

SeaBASS

SeaDAS

Giovanni Data Visualizer

Worldview Mapping Interface

Mission Quality Monitor

HELP

Level 1 & 2 Browser Help

Level 3 & 4 Browser Help

File Search Help

Search & Download Help

Earthdata Forum

Join Mailing List

How to Cite Data

CONNECT WITH US
COMPLIANCE

Web Privacy Policy

Data & Information Policy

Communications Policy

Accessibility Policy

Freedom of Information Act

USA.gov

Responsible NASA Official: Sean Bailey
Curator: OceanColor Webmaster