subroutine mtshadow(nopr,iout,&
   thetasun,thetaview,phisun,phiview,nreal,nim,betadeg,&
   sigma2,shadow)

! *********************
! MT refers to Mishchenko and Travis, "Satellite retrieval
! of aerosol properties over the ocean using polarization as
! well as intensity sunlight"
! JGR, vol 102, No. D14, pages 16989-17013, July 27, 1997

! ******************************
! 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
! sigma2       mean square slope

! Output
! shadow

! ******************************
    real :: nreal,nim
    real :: mu,mu2
    real :: nu,nu2
    real :: shadow

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

! ******************************
! note that you use betadeg (the tilt angle)
!hmmm, this could be in error

! nu is  scattering normal angle
! nu2 is incident normal angle
    mu=thetaview-betadeg
    mu2=thetasun-betadeg

    sigma=sqrt(sigma2)

    cmu=cos(convr*mu)
    cmu2=cos(convr*mu2)

    factor=sqrt(2.00)
    factor2=sqrt(2.0/pi)

     b1=sqrt(1.00-(cmu*cmu))
    nu=cmu/(sigma*factor*b1)

     b2=sqrt(1.00-(cmu2*cmu2))
    nu2=cmu2/(sigma*factor*b2)

     c1=(sigma/cmu)*factor2*b1
     a1=c1*exp(-nu*nu)
     ernu=erfc(nu)
    fnu=0.5*(a1-ernu)

     c2=(sigma/cmu2)*factor2*b2
     a2=c2*exp(-nu2*nu2)
     ernu2=erfc(nu2)
    fnu2=0.5*(a2-ernu2)

    shadow=1.00/(1.00+fnu+fnu2)

! ******************************
   if (nopr .eq. 1) then 
    write(iout,100) thetasun,thetaview,betadeg
    100 format(/,2x,'mtshadow: thetasun,thetaview,betadeg',/,&
    2x,3(1x,f10.4))
    write(iout,102) mu,mu2
    102 format(2x,'mtshadow: mu,mu2',/,&
    2x,2(1x,f10.4))
    write(iout,103) cmu,cmu2
    103 format(2x,'mtshadow: cmu,cmu2',/,&
    2x,2(1x,f10.4))
    write(iout,105) sigma2,sigma
    105 format(2x,'mtshadow: sigma2,sigma',/,&
    2x,2(1x,f10.4))
    write(iout,107) nu,nu2
    107 format(2x,'mtshadow: nu,nu2',/,&
    2x,2(1x,f10.4))
    write(iout,108) b1,b2
    108 format(2x,'mtshadow: b1,b2',/,&
    2x,2(1x,f10.4))
    write(iout,109) ernu,ernu2
    109 format(2x,'mtshadow: ernu,ernu2',/,&
    2x,2(1x,f10.4))
    write(iout,110) fnu,fnu2
    110 format(2x,'mtshadow: fnu,fnu2',/,&
    2x,2(1x,f10.4))
    write(iout,115) shadow
    115 format(2x,'mtshadow: shadow',/,&
    2x,4(1x,f10.4))
   endif

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