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 mm. * C * 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,AXI) PARAMETER (NC=10, NG=60) IMPLICIT REAL*8 (A-H,O-Z) REAL*8 X(NG),W(NG),C(0:NC) COMMON /CDROP/ C,R0V 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