OB.DAAC Logo
NASA Logo
Ocean Color Science Software

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