subroutine otarmatrix(nopr,iout,&
   thetasun,thetaview,phisun,phiview,nreal,nim,betadeg,&
   m11,m12,m33,m44,ratiom12m11)

! Caution: seems your implementation is good only for the 
!          scattering plane (with sun vector and observation vector in same plane)

! ******************************
! See Ota "Matrix formulations of radiative transfer including
! the polarization effect in a coupled atmosphere-ocean system"
! vol 111, pgs 878-894, 2010

! The full 4x4 Mueller-Stokes matrix (MSmatrix) (A.1) is 
!    m11  m12  0   0 
!    m12  m11  0   0
!    0    0    m33 0
!    0    0    0   m33

! For input Stokes vector V0=(s0,s1,s2,s3)T 
!  T, do transpose to get column vector
! You get output vector V1=(soR,s1R,s2R,s3R)T
! With    V1 = MSmatrix V0

! ******************************
! Input
! thetasun     solar zenith angle (degrees)
! thetaview    sensor zenith angle (degrees)
! phisun       solar azimuth angle (degrees)
! phiview      ssensor azimuth angle (degrees)
! betadeg      surface tilt angle (degrees)
! nreal,nim    real and imaginary indices

! Output
! m11,m12      the 1,1 and 1,2 matrix elements of the reflection matrix
! m33,m44      the 1,1 and 1,2 matrix elements of the reflection matrix

! ******************************
    real :: nreal,nim
    real :: nreal2
    real :: m11,m12,m33,m44
    real :: ratiom12m11

! ******************************
! 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

! ******************************
! Asuume you can use previous info to calculate the
! reflection matrix elements

    thetanew=thetasun-betadeg
     s1=convr*thetanew
    csbeta=cos(s1)
    snbeta=sin(s1)
    sn2=snbeta*snbeta

    nreal2=nreal*nreal

! equation A.2
    a1=nreal2*csbeta
    a2=sqrt(nreal2-sn2)
    rp=(a1-a2)/(a1+a2)
    rp2=rp*rp

! equation A.3
    b1=csbeta
    rperp=(b1-a2)/(b1+a2)
    rperp2=rperp*rperp

    m11=0.5*(rp2+rperp2)
    m12=0.5*(rp2-rperp2)

    m33=rp*rperp
    m44=m33

   ratiom12m11=m12/m11
   ratiom33m11=m33/m11

! ******************************
   if (nopr .eq. 1) then 
    write(iout,100) thetasun,betadeg,nreal,thetanew
    100 format(/,2x,'otarmatrix: thetasun,betadeg,nreal,thetanew',/,&
    2x,5(1x,f10.4))
    write(iout,105) rp,rperp
    105 format(2x,'otarmatrix: rp,rperp',/,&
    2x,2(1x,f10.4))
    write(iout,107) rp2,rperp2
    107 format(2x,'otarmatrix: rp2,rperp2',/,&
    2x,2(1x,f10.4))
    write(iout,115) m11,m12,m33,m44
    115 format(2x,'otarmatrix: m11,m12,m33,m44',/,&
    2x,4(1x,f10.4))
    write(iout,120) ratiom12m11,ratiom33m11
    120 format(2x,'otarmatrix: ratiom12m11,ratiom33m11',/,&
    2x,2(1x,f10.4))
   endif

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