!+ Subroutine to calculate single scattering after Takano and Liou.
!
! Method:
!	The cross-sections of extinction and scattering and the 
!	asymmetry of hexagonal ice crystals are calculated using
!	a fit to Takano and Liou's calculations based on geometrical
!	optics. Note that this is based entirely on the fit in their
!	paper and therefore may not be very well sampled. The fit
!	to the asymmetry does not work well near 3 microns.
!
! Current Owner of Code: J. M. Edwards
!
! History:
!	Version		Date			Comment
!	1.0		12-04-95		First Version under RCS
!						(J. M. Edwards)
!
! Description of Code:
!   FORTRAN 77  with extensions listed in documentation.
!
!- ---------------------------------------------------------------------
      SUBROUTINE TAKANO_LIOU_SCATTER(WAVELENGTH, RADIUS, ASPECT
     &   , N_RELATIVE
     &   , XS_SCATTER, XS_EXT, ASYMMETRY
     &   )
!
!
      IMPLICIT NONE
!
!
!     INCLUDE HEADER FILES.
      INCLUDE 'incl/math_cnst.ccf'
!
!     DUMMY ARGUMENTS
      REAL*8	!, INTENT(IN)
     &     WAVELENGTH
!             WAVELENGTH
     &   , RADIUS
!             RADIUS OF CIRCUMSCRIPTION
     &   , ASPECT
!             ASPECT RATIO OF CRYSTAL
      COMPLEX*16	!, INTENT(IN)
     &     N_RELATIVE
!             RELATIVE REFRACTIVE INDEX
      REAL*8	!, INTENT(OUT)
     &     XS_SCATTER
!             SCATTERING CROSS_SECTION
     &   , XS_EXT
!             EXTINCTIVE CROSS_SECTION
     &   , ASYMMETRY
!             ASYMMTRY OF SCATTERING
!
!     LOCAL VARIABLES
      REAL*8
     &     AREA
!             SURFACE AREA OF CRYSTAL
     &   , RADIUS_EQV
!             EQUIVALENT RADIUS OF CRYSTAL
     &   , MI_SIZE
!             IM(REF. INDEX)* SIZE PARAMETER
     &   , OMEGA
!             SINGLE SCATTERING ALBEDO
     &   , Z
!             TAKANO AND LIOU'S Z
     &   , INV_REFRACT
!             DIFFERENCE OF INVERSE
!                REFRACTIVE INDICES
      REAL*8
     &     SQRT_3_BY_2
!             CONSTANT
     &   , E0
!             FITTING PARAMETER FOR ASYMMETRY
     &   , E1
!             FITTING PARAMETER FOR ASYMMETRY
     &   , F0
!             FITTING PARAMETER FOR ASYMMETRY
     &   , F1
!             FITTING PARAMETER FOR ASYMMETRY
     &   , E
!             FITTING PARAMETER FOR ASYMMETRY
     &   , F
!             FITTING PARAMETER
!
!     THESE VALUSE WERE DERIVED TAKING THE EFFECTIVE RADIUS IN MICRONS.
      DATA
     &     E0/1.4750D-02/
     &   , E1/8.6611D-01/
     &   , F0/-2.4652D-01/
     &   , F1/3.7026D-01/
!
!
      SQRT_3_BY_2=SQRT(3.0D+00)/2.0D+00
!
      AREA=6.0D+00*RADIUS*RADIUS*(SQRT_3_BY_2+2.0D+00*ASPECT)
      RADIUS_EQV=SQRT(AREA/(4.0D+00*PI))
      MI_SIZE=2.0D+00*PI*RADIUS*DIMAG(N_RELATIVE)/WAVELENGTH
!
!     WITH GEOMETRICAL OPTICS THE EXTINCTIVE CROSS-SECTION IS TWICE
!     THE GEOMETRICAL CROSS-SECTION, BUT WHEN THE RANDOM ORIENTATION OF
!     THE CRYSTALS IS ACCOUNTED FOR WE GET HALF THE SURFACE AREA.
      XS_EXT=5.0D-01*AREA
!
!     THE SINGLE SCATTERING ALBEDO IS FITTED BY TAKANO AND LIOU AS
!     A QUARTIC OVER THE RANGE .4 TO 2.5 MICRONS. WE EXTEND THIS BY 
!     TAKING IT TO BE 1 BELOW .4 MICRONS.
      IF (WAVELENGTH.LT.4.0D-07) THEN
         OMEGA=1.0D+00
      ELSE
         Z=2.0D+00*MI_SIZE
     &      *(SQRT(2.7D+01)*ASPECT)/(SQRT(3.0D+00)+4.0D+00*ASPECT)
         OMEGA=1.0D+00+Z*(-1.1128D+00+Z*(2.5576D+00+Z*(-5.6257D+00
     &      +Z*5.9498D+00)))
!        CHECK THAT OMEGA LIES WITHIN ITS RANGE OF VALIDITY AND
!        USE THE FORM BASED ON ADT IF NOT. WE TEST FOR OMEGA EXCEEDING
!        1.01 TO ENABLE US TO SEPARATE ABJECT FAILURE FROM POSSIBLE
!        SMALL ERRORS IN THE FIT WHER OMEGA SHOULD BE CLOSE TO 1.0.
         IF ( (OMEGA.GT.1.01D+00).OR.(OMEGA.LT.0.75D+00) )THEN
!           OUTSIDE THE RANGE OF VALIDITY OF THE QUARTIC, MODIFIED ADT
!           IS USED. THE SIZE PARAMETER HERE IS DEFINED IN TERMS OF
!           THE EFFECTIVE RADIUS, SO MI_SIZE IS RESCALED.
            MI_SIZE=MI_SIZE*RADIUS_EQV/RADIUS
            OMEGA=5.0D-01+(MI_SIZE+3.0D+00)
     &         /(6.0D+00*(MI_SIZE+1.0D+00)**3)
         ENDIF
!        NOW DEAL WITH ANY SMALL EXCESS ABOVE 1.
         OMEGA=MIN(1.0D+00, OMEGA)
      ENDIF
      XS_SCATTER=OMEGA*XS_EXT
!
!     THE ASYMMETRY IS DERIVED FROM LINEAR REGRESSION IN TERMS OF
!     LOG(RADIUS_EQV) OVER THE DATA FOR HEXAGONS AT .55 AND 2.2 MICRONS.
!     THE DEPENDENCE ON THE REFRACTIVE INDEX SEEMS TO FIT THE AVERAGED
!     DATA GIVEN FOR THE CLOUD DISTRIBUTIONS EXCEPT AT 3 MICRONS
!     WHERE THE ABSORPTION BECOMES SIGNIFICANT.
      INV_REFRACT=1.0D+00/DBLE(N_RELATIVE)
      E=E0+E1*INV_REFRACT
      F=F0+F1*INV_REFRACT
      ASYMMETRY=E+F*LOG(1.0D+06*RADIUS_EQV)
      ASYMMETRY=MAX(0.0D+00, MIN(ASYMMETRY, 9.6D-01))
!
!
!
      RETURN
      END
