subroutine fresnelwik(nopr,iout,nair,nwater,csthetasun,snthetasun,&
   csthetaview,snthetaview,csphidiff,snphidiff,&
   thetasun,thetaview,phidiff,&
   thetai,fresnelnormal,fresnelpar,fresneleff)

! ******************
! See Wikipedia

! 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: fresnelnormal  fresnel reflection coefficient normal to 
!                                  incident-reflection plane
!                   fresnelpar     fresnel reflection coefficient parallel to
!                                  incident-reflection plane
!                   fresneleff     0.5 * ( fresnelnormal + fresnelpar )

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

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

! ******************
! thetai is equation (5) of Breon and Heriot (2006) paper
    a1=(csthetasun*csthetaview)+(snthetasun*snthetaview*csphidiff)
    a2=0.5*acos(a1)
    thetai=a2/convr

! ******************
! As per the Wikipedia article
    ratio=nair/nwater

    cs=cos(convr*thetai)
    sn=sin(convr*thetai)

! ******************
    a1=ratio*sn
    a2=1.00-(a1*a1)
    a3=sqrt(a2)

    a4=nair*cs

    a4=a4-(nwater*a3)
    a5=a4+(nwater*a3)

    a6=a4/a5

    fresnelnormal=abs(a6*a6)

! ******
    b1=nwater*cs

    b4=(nair*a3)-b1
    b5=(nair*a3)+b1

    b6=b4/b5

    fresnelpar=abs(b6*b6)

! ******
! The effective fresnel reflection
    fresneleff=0.5*(fresnelnormal+fresnelpar)

! ******
     if (nopr .eq. 1) then
      write(iout,100) nair,nwater,thetasun,thetaview,phidiff,&
      csthetasun,csthetaview,snthetasun,snthetaview,csphidiff,&
      thetai,fresnelnormal,fresnelpar,fresneleff
      100 format(/,2x,"fresnelwik: nair,nwater ",/,&
      2x,2(1x,f12.5),/,&
      2x,"fresnelwik: thetasun,thetaview,phidiff",/,&
      2x,3(1x,f12.5),/,&
      2x,"fresnelwik: csthetasun,csthetaview,snthetasun,",&
         "snthetaview,csphidiff",/,&
      2x,5(1x,f12.5),/,&
      2x,"fresnelwik: thetai ",/,&
      2x,1(1x,f12.5),/,&
      2x,"fresnelwik: fresnelnormal,fresnelpar",/,&
      2x,2(1x,f12.5),/,&
      2x,"fresnelwik: fresneleff",/,&
      2x,1(1x,f12.5))
     endif

! ******
      return
      end