C********************************************************************
C                                                                   *
C   Computation of the coefficients for the generalized Chebyshev   *
C   shape by interpolating the values given in Table 1 of Chuang    *
C   and Beard (1990). It is assumed that the particle radius AXI    *
C   is given in m.                                                  *
C                                                                   *
C   Changed by RMAY to take AXI in m instead of mm                  *
C   Reference:                                                      *
C   Chuang, C.C. and K.V. Beard, 1990: A numerical model for the    *
C     equilibrium shape of electrified raindrops. J. Atmos. Sci.,   *
C     47(11), 1374-1389.                                            *
C                                                                   *
C********************************************************************

      SUBROUTINE DROP (RAT,AXI0)
      PARAMETER (NC=10, NG=60)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 X(NG),W(NG),C(0:NC),AXI
      COMMON /CDROP/ C,R0V
      AXI = AXI0 * 1000.
      IF (AXI.LT.0.3332103D0) THEN
         C(0) = 0.0D0
         C(1) = 0.0D0
         C(2) = 0.0D0
         C(3) = 0.0D0
         C(4) = 0.0D0
         C(5) = 0.0D0
         C(6) = 0.0D0
         C(7) = 0.0D0
         C(8) = 0.0D0
         C(9) = 0.0D0
         C(10) =0.0D0
      ELSE IF (AXI.LT.0.5D0) THEN
         C(0) = (-0.0028D0 / 0.1667897D0) * (AXI - 0.3332103D0)
         C(1) = (0.003D0 / 0.1667897D0) * (AXI - 0.3332103D0)
         C(2) = (-0.0083D0 / 0.1667897D0) * (AXI - 0.3332103D0)
         C(3) = (0.0022D0 / 0.1667897D0) * (AXI - 0.3332103D0)
         C(4) = (-0.0003D0 / 0.1667897D0) * (AXI - 0.3332103D0)
         C(5) = (-0.0002D0 / 0.1667897D0) * (AXI - 0.3332103D0)
         C(6) = (0.0001D0 / 0.1667897D0) * (AXI - 0.3332103D0)
         C(7) = 0.0D0
         C(8) = 0.0D0
         C(9) = 0.0D0
         C(10) = 0.0D0
      ELSE IF (AXI.LT.0.75D0) THEN
         C(0) = -0.0028D0 - 0.0176D0 * (AXI - 0.5D0)
         C(1) = 0.003D0 + 0.016D0 * (AXI - 0.5D0)
         C(2) = -0.0083D0 - 0.0508D0 * (AXI - 0.5D0)
         C(3) = 0.0022D0 + 0.014D0 * (AXI - 0.5D0)
         C(4) = -0.0003D0 - 0.0012D0 * (AXI - 0.5D0)
         C(5) = -0.0002D0 - 0.002D0 * (AXI - 0.5D0)
         C(6) = 0.0001D0 + 0.0008D0 * (AXI - 0.5D0)
         C(7) = 0.0D0
         C(8) = -0.0004D0 * (AXI - 0.5D0)
         C(9) = 0.0D0
         C(10) = 0.0001D0 * (AXI - 0.5D0)
      ELSE IF (AXI.LT.1.0D0) THEN
         C(0) = -0.0072D0 - 0.0248D0 * (AXI - 0.75D0)
         C(1) = 0.007D0 + 0.0192D0 * (AXI - 0.75D0)
         C(2) = -0.021D0 - 0.07D0 * (AXI - 0.75D0)
         C(3) = 0.0057D0 + 0.0172D0 * (AXI - 0.75D0)
         C(4) = -0.0006D0 + 0.0004D0 * (AXI - 0.75D0)
         C(5) = -0.0007D0 - 0.004D0 * (AXI - 0.75D0)
         C(6) = 0.0003D0 + 0.0012D0 * (AXI - 0.75D0)
         C(7) = 0.0004D0 * (AXI - 0.75D0)
         C(8) = -0.0001D0 - 0.0008D0 * (AXI - 0.75D0)
         C(9) = 0.0004D0 * (AXI - 0.75D0)
         C(10) = 0.0001D0
      ELSE IF (AXI.LT.1.25D0) THEN
         C(0) = -0.0134D0 - 0.0308D0 * (AXI - 1.0D0)
         C(1) = 0.0118D0 + 0.0248D0 * (AXI - 1.0D0)
         C(2) = -0.0385D0 - 0.0828D0 * (AXI - 1.0D0)
         C(3) = 0.01D0 + 0.0188D0 * (AXI - 1.0D0)
         C(4) = -0.0005D0 + 0.0036D0 * (AXI - 1.0D0)
         C(5) = -0.0017D0 - 0.006D0 * (AXI - 1.0D0)
         C(6) = 0.0006D0 + 0.0016D0 * (AXI - 1.0D0)
         C(7) = 0.0001D0 + 0.0008D0 * (AXI - 1.0D0)
         C(8) = -0.0003D0 - 0.0008D0 * (AXI - 1.0D0)
         C(9) = 0.0001D0
         C(10) = 0.0001D0 + 0.0004D0 * (AXI - 1.0D0)
      ELSE IF (AXI.LT.1.5D0) THEN
         C(0) = -0.0211D0 - 0.0344D0 * (AXI - 1.25D0)
         C(1) = 0.018D0 + 0.0268D0 * (AXI - 1.25D0)
         C(2) = -0.0592D0 - 0.0896D0 * (AXI - 1.25D0)
         C(3) = 0.0147D0 + 0.0164D0 * (AXI - 1.25D0)
         C(4) = 0.0004D0 + 0.008D0 * (AXI - 1.25D0)
         C(5) = -0.0032D0 - 0.008D0 * (AXI - 1.25D0)
         C(6) = 0.001D0 + 0.0012D0 * (AXI - 1.25D0)
         C(7) = 0.0003D0 + 0.002D0 * (AXI - 1.25D0)
         C(8) = -0.0005D0 - 0.0012D0 * (AXI - 1.25D0)
         C(9) = 0.0001D0
         C(10) = 0.0002D0 + 0.0008D0 * (AXI - 1.25D0)
      ELSE IF (AXI.LT.1.75D0) THEN
         C(0) = -0.0297D0 - 0.0364D0 * (AXI - 1.5D0)
         C(1) = 0.0247D0 + 0.0248D0 * (AXI - 1.5D0)
         C(2) = -0.0816D0 - 0.0904D0 * (AXI - 1.5D0)
         C(3) = 0.0188D0 + 0.0132D0 * (AXI - 1.5D0)
         C(4) = 0.0024D0 + 0.0116D0 * (AXI - 1.5D0)
         C(5) = -0.0052D0 - 0.0092D0 * (AXI - 1.5D0)
         C(6) = 0.0013D0 + 0.0008D0 * (AXI - 1.5D0)
         C(7) = 0.0008D0 + 0.0028D0 * (AXI - 1.5D0)
         C(8) = -0.0008D0 - 0.0016D0 * (AXI - 1.5D0)
         C(9) = 0.0001D0 - 0.0004D0 * (AXI - 1.5D0)
         C(10) = 0.0004D0 + 0.0012D0 * (AXI - 1.5D0)
      ELSE IF (AXI.LT.2.0D0) THEN
         C(0) = -0.0388D0 - 0.0372D0 * (AXI - 1.75D0)
         C(1) = 0.0309D0 + 0.02D0 * (AXI - 1.75D0)
         C(2) = -0.1042D0 - 0.0884D0 * (AXI - 1.75D0)
         C(3) = 0.0221D0 + 0.0092D0 * (AXI - 1.75D0)
         C(4) = 0.0053D0 + 0.0152D0 * (AXI - 1.75D0)
         C(5) = -0.0075D0 - 0.0096D0 * (AXI - 1.75D0)
         C(6) = 0.0015D0
         C(7) = 0.0015D0 + 0.004D0 * (AXI - 1.75D0)
         C(8) = -0.0012D0 - 0.0016D0 * (AXI - 1.75D0)
         C(9) = -0.0008D0 * (AXI - 1.75D0)
         C(10) = 0.0007D0 + 0.0012D0 * (AXI - 1.75D0)
      ELSE IF (AXI.LT.2.25D0) THEN
         C(0) = -0.0481D0 - 0.0368D0 * (AXI - 2.0D0)
         C(1) = 0.0359D0 + 0.0168D0 * (AXI - 2.0D0)
         C(2) = -0.1263D0 - 0.0844D0 * (AXI - 2.0D0)
         C(3) = 0.0244D0 + 0.0044D0 * (AXI - 2.0D0)
         C(4) = 0.0091D0 + 0.0184D0 * (AXI - 2.0D0)
         C(5) = -0.0099D0 - 0.0088D0 * (AXI - 2.0D0)
         C(6) = 0.0015D0 - 0.0016D0 * (AXI - 2.0D0)
         C(7) = 0.0025D0 + 0.0044D0 * (AXI - 2.0D0)
         C(8) = -0.0016D0 - 0.0012D0 * (AXI - 2.0D0)
         C(9) = -0.0002D0 - 0.0016D0 * (AXI - 2.0D0)
         C(10) = 0.001D0 + 0.0012D0 * (AXI - 2.0D0)
      ELSE IF (AXI.LT.2.5D0) THEN
         C(0) = -0.0573D0 - 0.0368D0 * (AXI - 2.25D0)
         C(1) = 0.0401D0 + 0.0136D0 * (AXI - 2.25D0)
         C(2) = -0.1474D0 - 0.08D0 * (AXI - 2.25D0)
         C(3) = 0.0255D0 + 0.0012D0 * (AXI - 2.25D0)
         C(4) = 0.0137D0 + 0.02D0 * (AXI - 2.25D0)
         C(5) = -0.0121D0 - 0.008D0 * (AXI - 2.25D0)
         C(6) = 0.0011D0 - 0.0028D0 * (AXI - 2.25D0)
         C(7) = 0.0036D0 + 0.0048D0 * (AXI - 2.25D0)
         C(8) = -0.0019D0 - 0.0008D0 * (AXI - 2.25D0)
         C(9) = -0.0006D0 - 0.002D0 * (AXI - 2.25D0)
         C(10) = 0.0013D0 + 0.0016D0 * (AXI - 2.25D0)
      ELSE IF (AXI.LT.2.75D0) THEN
         C(0) = -0.0665D0 - 0.036D0 * (AXI - 2.5D0)
         C(1) = 0.0435D0 + 0.012D0 * (AXI - 2.5D0)
         C(2) = -0.1674D0 - 0.0756D0 * (AXI - 2.5D0)
         C(3) = 0.0258D0 - 0.0028D0 * (AXI - 2.5D0)
         C(4) = 0.0187D0 + 0.022D0 * (AXI - 2.5D0)
         C(5) = -0.0141D0 - 0.0064D0 * (AXI - 2.5D0)
         C(6) = 0.0004D0 - 0.0044D0 * (AXI - 2.5D0)
         C(7) = 0.0048D0 + 0.0052D0 * (AXI - 2.5D0)
         C(8) = -0.0021D0
         C(9) = -0.0011D0 - 0.0024D0 * (AXI - 2.5D0)
         C(10) = 0.0017D0 + 0.0016D0 * (AXI - 2.5D0)
      ELSE IF (AXI.LT.3.0D0) THEN
         C(0) = -0.0755D0 - 0.0352D0 * (AXI - 2.75D0)
         C(1) = 0.0465D0 + 0.0028D0 * (AXI - 2.75D0)
         C(2) = -0.1863D0 - 0.0708D0 * (AXI - 2.75D0)
         C(3) = 0.0251D0 - 0.0044D0 * (AXI - 2.75D0)
         C(4) = 0.0242D0 + 0.0228D0 * (AXI - 2.75D0)
         C(5) = -0.0157D0 - 0.0044D0 * (AXI - 2.75D0)
         C(6) = -0.0007D0 - 0.0056D0 * (AXI - 2.75D0)
         C(7) = 0.0061D0 + 0.0048D0 * (AXI - 2.75D0)
         C(8) = -0.0021D0 + 0.0004D0 * (AXI - 2.75D0)
         C(9) = -0.0017D0 - 0.0032D0 * (AXI - 2.75D0)
         C(10) = 0.0021D0 + 0.0012D0 * (AXI - 2.75D0)
      ELSE IF (AXI.LT.3.25D0) THEN
         C(0) = -0.0843D0 - 0.0348D0 * (AXI - 3.0D0)
         C(1) = 0.0472D0 + 0.006D0 * (AXI - 3.0D0)
         C(2) = -0.204D0 - 0.0668D0 * (AXI - 3.0D0)
         C(3) = 0.024D0 - 0.0072D0 * (AXI - 3.0D0)
         C(4) = 0.0299D0 + 0.0236D0 * (AXI - 3.0D0)
         C(5) = -0.0168D0 - 0.0028D0 * (AXI - 3.0D0)
         C(6) = -0.0021D0 - 0.0064D0 * (AXI - 3.0D0)
         C(7) = 0.0073D0 + 0.0044D0 * (AXI - 3.0D0)
         C(8) = -0.002D0 + 0.0016D0 * (AXI - 3.0D0)
         C(9) = -0.0025D0 - 0.0036D0 * (AXI - 3.0D0)
         C(10) = 0.0024D0 + 0.0012D0 * (AXI - 3.0D0)
      ELSE IF (AXI.LT.3.5D0) THEN
         C(0) = -0.093D0 - 0.0336D0 * (AXI - 3.25D0)
         C(1) = 0.0487D0 + 0.002D0 * (AXI - 3.25D0)
         C(2) = -0.2207D0 - 0.0628D0 * (AXI - 3.25D0)
         C(3) = 0.0222D0 - 0.0092D0 * (AXI - 3.25D0)
         C(4) = 0.0358D0 + 0.0244D0 * (AXI - 3.25D0)
         C(5) = -0.0175D0 - 0.0012D0 * (AXI - 3.25D0)
         C(6) = -0.0037D0 - 0.0076D0 * (AXI - 3.25D0)
         C(7) = 0.0084D0 + 0.0036D0 * (AXI - 3.25D0)
         C(8) = -0.0016D0 + 0.0016D0 * (AXI - 3.25D0)
         C(9) = -0.0034D0 - 0.0036D0 * (AXI - 3.25D0)
         C(10) = 0.0027D0 + 0.0012D0 * (AXI - 3.25D0)
      ELSE IF (AXI.LT.4.0D0) THEN
         C(0) = -0.1014D0 - 0.0346D0 * (AXI - 3.5D0)
         C(1) = 0.0492D0 - 0.002D0 * (AXI - 3.5D0)
         C(2) = -0.2364D0 - 0.0572D0 * (AXI - 3.5D0)
         C(3) = 0.0199D0 - 0.0102D0 * (AXI - 3.5D0)
         C(4) = 0.0419D0 + 0.0248D0 * (AXI - 3.5D0)
         C(5) = -0.0178D0 + 0.0014D0 * (AXI - 3.5D0)
         C(6) = -0.0056D0 - 0.0088D0 * (AXI - 3.5D0)
         C(7) = 0.0093D0 + 0.0028D0 * (AXI - 3.5D0)
         C(8) = -0.0012D0 + 0.0028D0 * (AXI - 3.5D0)
         C(9) = -0.0043D0 - 0.0042D0 * (AXI - 3.5D0)
         C(10) = 0.003D0 + 0.0004D0 * (AXI - 3.5D0)
      ELSE IF (AXI.LT.4.5D0) THEN
         C(0) = -0.1187D0 - 0.0282D0 * (AXI - 4.0D0)
         C(1) = 0.0482D0 - 0.0158D0 * (AXI - 4.0D0)
         C(2) = -0.265D0 - 0.0498D0 * (AXI - 4.0D0)
         C(3) = 0.0148D0 - 0.0084D0 * (AXI - 4.0D0)
         C(4) = 0.0543D0 + 0.0238D0 * (AXI - 4.0D0)
         C(5) = -0.0171D0 + 0.0036D0 * (AXI - 4.0D0)
         C(6) = -0.01D0 - 0.0092D0 * (AXI - 4.0D0)
         C(7) = 0.0107D0 + 0.0008D0 * (AXI - 4.0D0)
         C(8) = 0.0002D0 + 0.0032D0 * (AXI - 4.0D0)
         C(9) = -0.0064D0 - 0.0034D0 * (AXI - 4.0D0)
         C(10) = 0.0032D0 - 0.0002D0 * (AXI - 4.0D0)
      ELSE
         C(0) = -0.1328D0
         C(1) = 0.0403D0
         C(2) = -0.2899D0
         C(3) = 0.0106D0
         C(4) = 0.0662D0
         C(5) = -0.0153D0
         C(6) = -0.0146D0
         C(7) = 0.0111D0
         C(8) = 0.0018D0
         C(9) = -0.0081D0
         C(10) = 0.0031D0
      ENDIF
      CALL GAUSS (NG,0,0,X,W)
      S=0D0
      V=0D0
      DO I=1,NG
         XI=DACOS(X(I))
         WI=W(I)
         RI=1D0+C(0)
         DRI=0D0
         DO N=1,NC
            XIN=XI*N
            RI=RI+C(N)*DCOS(XIN)
            DRI=DRI-C(N)*N*DSIN(XIN)
         ENDDO
         SI=DSIN(XI)
         CI=X(I)
         RISI=RI*SI
         S=S+WI*RI*DSQRT(RI*RI+DRI*DRI)
         V=V+WI*RI*RISI*(RISI-DRI*CI)
      ENDDO
      RS=DSQRT(S*0.5D0)
      RV=(V*3D0*0.25D0)**(1D0/3D0)
      IF (DABS(RAT-1D0).GT.1D-8) RAT=RV/RS
      R0V=1D0/RV
!      WRITE (6,1000) R0V
!      DO N=0,NC
!         WRITE (6,1001) N,C(N)
!      ENDDO
! 1000 FORMAT ('r_0/r_ev=',F7.4)
! 1001 FORMAT ('c_',I2,'=',F7.4)
      RETURN
      END
