subroutine fresnelsayer(nopr,iout,nair,nwater,csthetasun,snthetasun,&
   csthetaview,snthetaview,csphidiff,snphidiff,&
   thetasun,thetaview,phidiff,&
   omega,omegaprime,fresnelsayerval)

! ******************
! See Sayer 2007 AOPP Oxford doc "A SEA SURFACE REFLECTANCE MODEL
! SUITABLE FOR USE WITH AATSR AEROSOL RETRIEVAL"

! Input from breon.f90
! cs and sn are cosines and sines
! thetasun is incident solar zenith angle 
! thetaview is incident sensor zenith angle 
! phidiff is near 180 degrees (near specular reflection)

! output parameter: fresnelsayerval

! ******************
      real :: nair,nwater
      real :: ratio,cs,sn
      real :: a1,a2,a3,a4,a5,a6
      real :: b1,b2,b3,b4,b5,b6
      real :: fresnelsayerval

! ******************************
! Helpful for conversions
!   pi=3.14159265
    pi=atan(1.)*4.
    twopi=2.0*pi
! radians per degree
    convr=pi/180.0

! ******************
! omega is thetai in equation (5) of Breon and Heriot (2006) paper
! Same as equation 12 of Sayer
    a1=(csthetasun*csthetaview)+(snthetasun*snthetaview*csphidiff)
    a2=0.5*acos(a1)

! In degrees
    omega=a2/convr

    csomega=cos(convr*omega)
    snomega=sin(convr*omega)

! ******************
    ratio=nair/nwater

! Calculate omegaprime using Snell's law
    a1=ratio*snomega
    a2=asin(a1)

! In degrees
    omegaprime=a2/convr

    a3=convr*(omega-omegaprime)
    a4=convr*(omega+omegaprime)

    b1=sin(a3)
    b2=sin(a4)

    b3=tan(a3)
    b4=tan(a4)

    a5=(b1/b2)**2.0
    a6=(b3/b4)**2.0

    fresnelsayerval=0.5*(a5+a6)

! ******
     if (nopr .eq. 1) then
      write(iout,100) nair,nwater,thetasun,thetaview,phidiff,&
      csthetasun,csthetaview,snthetasun,snthetaview,csphidiff,&
      omega,omegaprime,fresnelsayerval
      100 format(/,2x,"fresnelsayer: nair,nwater ",/,&
      2x,2(1x,f12.5),/,&
      2x,"fresnelsayer: thetasun,thetaview,phidiff",/,&
      2x,3(1x,f12.5),/,&
      2x,"fresnelsayer: csthetasun,csthetaview,snthetasun,",&
         "snthetaview,csphidiff",/,&
      2x,5(1x,f12.5),/,&
      2x,"fresnelsayer: omega,omegaprime ",/,&
      2x,2(1x,f12.5),/,&
      2x,"fresnelsayer: fresnelsayerval",/,&
      2x,1(1x,f12.5))
     endif

! ******
      return
      end