OB.DAAC Logo
NASA Logo
Ocean Color Science Software

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