subroutine breon(nopr,iout,&
   thetasun,thetaview,phisun,phiview,windaz,windspeed,nreal,nim,fresnel,&
   phidiff,windazdiff,rsp,rspdiv,shadow,&
   m11h,m12h,m33h,m44h,&
   reflm11,reflm12,reflm33,reflm44,&
   reflm11b,reflm12b,reflm33b,reflm44b,&
   compareb,gc,exvalsig,cs4,betadeg,ratiom12m11)

! See Breon and Henriot, JGR, 2006
! Spaceborene observations of ocean glint reflectance and modeling of 
! wave slope distributions
! vol 111, C06005, doi:10.1029/2005JC003343

! ******************************
! Input
! thetasun     solar zenith angle (degrees)
! thetaview    sensor zenith angle (degrees)
! phisun       solar azimuth angle (degrees)
! phiview      ssensor azimuth angle (degrees)
! windaz       wind azimuth angle
! windspeed    m/sec
! nreal,nim    real and imaginary indices
! phidiff      relative azimuth angle between solar and sensor azimuths, phiview-phisun
!              for specular reflection phidiff is near 180 degrees
! windazdiff   difference in wind and solar azimuths, as per Breon instructions
!                   windazdiff=windaz+180.0-phisun

! Output
! fresnel      fresnel reflection coefficient
! rsp          the brdf (equation (4) of Breon)

! ******************************
! Example

! Have the sun in y=0 xz plane in a x-y-z grid, in the negative x side of the grid
!    phisun=180.0
! Have the sensor in y=0 xz plane in the x-y-z grid, in the positive x side of the grid
!    phiview=0.0
! Have the wind speed be 8.0 m/sec
!    windspeed=8.0
! Have the wind going from west to east (windaz=0.0 in z=0  x-y cartesian grid)
!    windaz=0.0
! For a solar zenith angle (thetasun=17.1 degrees) greater than 
!    thetasun=17.1
! the sensor zenith angle (thetaview=13.52 degrees), the
!    thetaview=13.52
! ocean reflective surface "faces the sun" 
! See Page 2 of 10 of Breon's paper for the geometry discussion
! (If thetasun=thetaview, then the reflective surface is flat in the x-y palne)
! Have windazdiff=0 (windazdiff=windaz+180.0-phisun=0.0+180.0-180.0=0.0)
! For this case, The angle between the x-y plane upwind wind unit vector (-1,0,0) and the solar
! direction unit vector (-1,0,0) is 0.0.
!    windazdiff=0.0

! This particular case yields Zx=-0.3127, Zy=0.0; Zup=Zx=-0.3127, Zcr=0.0
! tilt angle beta=1.79  (0.5 times (17.1 - 13.52))
! zeta and eta are -0.1924 and 0.0000
! equation(7) pval is is 8.04
! (The Gram Charlier part of equation(7) is gc=1.11)
! Final equation(4) Rsp is 0.136

! ******************************
    real :: nreal,nim
    real :: nair,nwater
    real :: omega,omegaprime
    real :: fresnel,fresnelsayerval

! Will use hafermanmatrix.f90 for the reflection matrix elements
    real :: m11h,m12h,m33h,m44h
    real :: reflm11,reflm12,reflm33,reflm44
    real :: ratiom12m11
    real :: exval,exvalsig
    real :: shadow

! ******************************
! Helpful for conversions
    radius=6370.0
!   pi=3.14159265
    pi=atan(1.)*4.
    twopi=2.0*pi
    circum=2.00*pi*radius
! radians per degree
    convr=pi/180.0
! km per degree (110.0)
    conv=circum/360.0

! ******************************
! Define cos and sin terms

     s1=convr*thetasun
    snthetasun=sin(s1)
    csthetasun=cos(s1)

     s1=convr*thetaview
    snthetaview=sin(s1)
    csthetaview=cos(s1)

     s1=convr*phisun
    snphisun=sin(s1)
    csphisun=cos(s1)

     s1=convr*phiview
    snphiview=sin(s1)
    csphiview=cos(s1)

! Values not used
!    s1=convr*windaz
!   snwindaz=sin(s1)
!   cswindaz=cos(s1)

! Careful
! For the specular direction phidiff is near 180.0 degrees
     s1=convr*phidiff
    snphidiff=sin(s1)
    csphidiff=cos(s1)

! Careful
! You calculate windazdiff outside of breon.f90
     s1=convr*windazdiff
    snwindazdiff=sin(s1)
    cswindazdiff=cos(s1)

! ******************************
! Calculate the Fresnel coefficient (should be near 0.02)

! oops, this does not fresnel varying from 0.02
! Gives correct value for small solar zenith angles
!       phir=thetasun-thetaview
!       s1=convr*phir
!     cos2chi=(csthetasun*csthetaview)+(snthetaview*snthetasun*cos(s1))
!     if (cos2chi .gt. 1.0) then 
!      cos2chi=0.99999999999
!     endif
!     if (cos2chi .lt. -1.0) then 
!      cos2chi=-0.99999999999
!     endif
!     coschi=sqrt(0.5*(1+cos2chi))
!     sinchi=sqrt(0.5*(1-cos2chi))

! Calculate the fresnel coefficient (on the order of 0.02 for O2 A-band)
       noprf=0
      if (nopr .eq. 1) then
       noprf=1
      endif

! hmmm, fresnel is near 0.02 (for small solar zenith angles) and
! also for larger solar zenith angles (booh)
!     call fresnelval(noprf,iout,nreal,nim,coschi,sinchi,fresnel)

! ******************************
! Calculate fresnel coefficients noraml and parallel to 
! the plane of incidence-reflection
    nair=1.00029
    nwater=nreal

! From wikipedia
!  call fresnelwik(noprf,iout,nair,nwater,csthetasun,snthetasun,&
!  csthetaview,snthetaview,csphidiff,snphidiff,&
!  thetasun,thetaview,phidiff,&
!  thetai,fresnelnormal,fresnelpar,fresneleff)

! By Sayer Oxford AOPP doc
   call fresnelsayer(nopr,iout,nair,nwater,csthetasun,snthetasun,&
   csthetaview,snthetaview,csphidiff,snphidiff,&
   thetasun,thetaview,phidiff,&
   omega,omegaprime,fresnelsayerval)

   fresnel=fresnelsayerval

! ******************************
! Equation (1) of Breon
    a1=-(snthetasun+(snthetaview*csphidiff))
    a2=csthetasun+csthetaview
   zx=a1/a2

    b1=-(snthetaview*snphidiff)
   zy=b1/a2

! ******************************
! Equation (2) of Breon
! Tilt angle
    sum= (zx*zx) + (zy*zy)
   tnbeta=sqrt(sum)
   betarad=atan(tnbeta)
   betadeg=betarad/convr

   csbeta=cos(betarad)
   cs4=csbeta**4.00

! ******************************
! Equation (3) of Breon
! Notice use of windazdiff values
! zup and zcr are upwind and crosswind slopes
   zup=(zx*cswindazdiff)+(zy*snwindazdiff)
   zcr=(zy*cswindazdiff)-(zx*snwindazdiff)

! ******************************
! Equation (9) of Breon
   sigma2up=1.0e-3+(3.16e-3*windspeed)
   sigma2cr=3.0e-3+(1.85e-3*windspeed)

   sigmaup=sqrt(sigma2up)
   sigmacr=sqrt(sigma2cr)

! can compare to gisscoxmunk
    sigma2breon=sigmaup*sigmacr

! ******************************
! Equation (7) of Breon
   zeta=zup/sigmaup
   eta=zcr/sigmacr
   zetaeta=-((zeta*zeta)+(eta*eta))/2.0
   exval=exp(zetaeta)

   p0=(1.0/(twopi*sigmaup*sigmacr))*exval

! for diagnostics
   exvalsig=exp(zetaeta)/sigma2breon

   c21=-9.0e-4*(windspeed*windspeed)
   c03=-0.45/(1.0+(exp(7.0-windspeed)))
   c40=0.3
   c04=0.4
   c22=0.12

   zeta2=zeta*zeta
   zeta4=zeta2*zeta2
   eta2=eta*eta
   eta3=eta2*eta
   eta4=eta2*eta2

   term1=-(c21/2.0)*(zeta2-1.0)*eta
   term2=-(c03/6.0)*(eta3-(3.0*eta))
   term3= (c40/24.0)*(zeta4-(6.0*zeta2)+3.0)
   term4= (c04/24.0)*(eta4-(6.0*eta2)+3.0)
   term5= (c22/4.0)*((zeta2-1.0)*(eta2-1.0))
   termsum=term1+term2+term3+term4+term5

! Gram-Charlier series
   gc=1.00+termsum

! The full pdf
   pval=p0*gc

! ******************************
! Equation (7) of Breon
   top=pi*fresnel*pval
   denom=4.0*csthetasun*csthetaview*cs4

! The reflectance generated by specular reflection
! at the ocean surface
   rsp=top/denom

! Can compare to coxmunk_new.f90 comparenew 
   compareb=rsp*(cs4/(gc*fresnel))

! ******************************
! Wrote several routines to calculate the reflection matrix terms
! based on several papers

! hmmm, m12 is of opposite sign as m12 of hafermanmatris.f90
! and otarmatrix.f90
!  call collett(nopr,iout,&
!  thetasun,thetaview,phisun,phiview,nreal,nim,betadeg,&
!  m11,m12,m33,m44,ratiom12m11,ratiom33m11)

! otarmatrix and hafermanmatrix.f90 have same
! m11,m12,m33,m44 values
!  call otarmatrix(nopr,iout,&
!  thetasun,thetaview,phisun,phiview,nreal,nim,betadeg,&
!  m11,m12,m33,m44,ratiom12m11,ratiom33m11)

! Gave result similar to gisscoxmunk when you calculated ratiom12m11
! and used it below
! Note that haferman reflection matrix elements are given by the 
! m11h,m12h,m33h,m44h values
   call hafermanmatrix(nopr,iout,&
   thetasun,thetaview,phisun,phiview,nreal,nim,betadeg,&
   m11h,m12h,m33h,m44h,ratiom12m11,ratiom33m11)

! See equation 87 and 88 of Mischchenko and Travis
! for the shadow term (seems to be 1.00)
   call mtshadow(nopr,iout,&
   thetasun,thetaview,phisun,phiview,nreal,nim,betadeg,&
   sigma2breon,shadow)

! Include the shadow term
   rspuse=rsp*shadow

! The final sun glint reflection matrix
! Caution: This may be an approximation
   reflm11=rspuse
   reflm12=rspuse*ratiom12m11
   reflm33=rspuse*ratiom33m11
   reflm44=rspuse*ratiom33m11

! rsp divided by fresnel
   rspdiv=rspuse/fresnel

! These values should be equal to reflm11,reflm12 etc
   reflm11b=rspdiv*m11h
   reflm12b=rspdiv*m12h
   reflm33b=rspdiv*m33h
   reflm44b=rspdiv*m44h

! ******************************
   if (nopr .eq. 1) then 
    write(iout,100) thetasun,thetaview,phisun,phiview
    100 format(/,2x,'breon: thetasun,thetaview,phisun,phiview',/,&
    2x,4(1x,f12.5))
    write(iout,105) phidiff
    105 format(2x,'breon: phidiff, near 180 for specular direction ',/,&
    2x,1(1x,f12.5))
    write(iout,110) windaz,windspeed
    110 format(2x,'breon: windaz,windspeed',/,&
    2x,2(1x,f12.5))
    write(iout,115) windazdiff
    115 format(2x,'breon: windazdiff=windaz+180-phisun',/,&
    2x,1(1x,f12.5))
    write(iout,120) nreal,nim,fresnel
    120 format(2x,'breon: nreal,nim,fresnel',/,&
    2x,3(1x,f12.5))
    write(iout,125) zx,zy
    125 format(2x,'breon: zx,zy ',/,&
    2x,2(1x,f12.5))
    write(iout,130) sum,tnbeta,betarad,betadeg,cs4
    130 format(2x,'breon:  sum,tnbeta,betarad,betadeg,cs4',/,&
    2x,5(1x,f12.5))
    write(iout,135) zup,zcr
    135 format(2x,'breon: zup,zcr',/,&
    2x,2(1x,f12.5))
    write(iout,140) sigmaup,sigmacr,sigma2breon
    140 format(2x,'breon: sigmaup,sigmacr,sigma2breon',/,&
    2x,3(1x,f12.5))
    write(iout,145) zeta,eta
    145 format(2x,'breon: zeta,eta',/,&
    2x,2(1x,f12.5))
    write(iout,150) zetaeta
    150 format(2x,'breon: zetaeta',/,&
    2x,1(1x,f12.5))
    write(iout,155) exval,exvalsig
    155 format(2x,'breon: exval,exvalsig',/,&
    2x,2(1x,f12.5))
    write(iout,160) p0
    160 format(2x,'breon: p0',/,&
    2x,1(1x,f12.5))
    write(iout,162) c21,c03,c40,c04,c22
    162 format(2x,'breon: c21,c03,c40,c04,c22 ',/,&
    2x,5(1x,f12.5))
    write(iout,165) term1,term2,term3,term4,term5
    165 format(2x,'breon: term1,term2,term3,term4,term5',/,&
    2x,5(1x,f12.5))
    write(iout,170) termsum
    170 format(2x,'breon: termsum ',/,&
    2x,1(1x,f12.5))
    write(iout,175) gc
    175 format(2x,'breon: gc ',/,&
    2x,1(1x,f12.5))
    write(iout,180) pval
    180 format(2x,'breon: pval=p0*gc',/,&
    2x,1(1x,f12.5))
    write(iout,185) top,denom 
    185 format(2x,'breon: top,denom',/,&
    2x,1(1x,f12.5))
    write(iout,190) rsp
    190 format(2x,'breon: rsp',/,&
    2x,1(1x,f12.5))
    write(iout,195) shadow
    195 format(2x,'breon: shadow',/,&
    2x,1(1x,f12.5))
    write(iout,200) rspuse
    200 format(2x,'breon: rspuse',/,&
    2x,1(1x,f12.5))
    write(iout,205) m11h,m12h,m33h,m44h
    205 format(2x,'breon: m11h,m12h,m33h,m44h',/,&
    2x,4(1x,f12.5))
    write(iout,210) reflm11,reflm12,reflm33,reflm44
    210 format(2x,'breon: reflm11,reflm12,reflm33,reflm44',/,&
    2x,4(1x,f12.5))
    write(iout,215) rspdiv
    215 format(2x,'breon: rspdiv',/,&
    2x,1(1x,f12.5))
    write(iout,220) reflm11b,reflm12b,reflm33b,reflm44b
    220 format(2x,'breon: reflm11b,reflm12b,reflm33b,reflm44b',/,&
    2x,4(1x,f12.5))
    write(iout,225) compareb
    225 format(2x,'breon: compareb and coxmunk_new.f90 comparenew',/,&
    2x,'breon: shoud be similar',/,&
    2x,1(1x,f12.5))
   endif

! ******************************
  return
  end