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) REAL*16 C(0:NC),R0V COMMON /CDROP/ C,R0V IF (AXI.LT.0.3332103D0) THEN C(0) = 0.0Q0 C(1) = 0.0Q0 C(2) = 0.0Q0 C(3) = 0.0Q0 C(4) = 0.0Q0 C(5) = 0.0Q0 C(6) = 0.0Q0 C(7) = 0.0Q0 C(8) = 0.0Q0 C(9) = 0.0Q0 C(10) = 0.0Q0 ELSE IF (AXI.LT.0.5D0) THEN C(0) = (-0.0028Q0 / 0.1667897Q0) * (AXI - 0.3332103D0) C(1) = (0.003Q0 / 0.1667897Q0) * (AXI - 0.3332103D0) C(2) = (-0.0083Q0 / 0.1667897Q0) * (AXI - 0.3332103D0) C(3) = (0.0022Q0 / 0.1667897Q0) * (AXI - 0.3332103D0) C(4) = (-0.0003Q0 / 0.1667897Q0) * (AXI - 0.3332103D0) C(5) = (-0.0002Q0 / 0.1667897Q0) * (AXI - 0.3332103D0) C(6) = (0.0001Q0 / 0.1667897Q0) * (AXI - 0.3332103D0) C(7) = 0.0Q0 C(8) = 0.0Q0 C(9) = 0.0Q0 C(10) = 0.0Q0 ELSE IF (AXI.LT.0.75D0) THEN C(0) = -0.0028Q0 - 0.0176Q0 * (AXI - 0.5D0) C(1) = 0.003Q0 + 0.016Q0 * (AXI - 0.5D0) C(2) = -0.0083Q0 - 0.0508Q0 * (AXI - 0.5D0) C(3) = 0.0022Q0 + 0.014Q0 * (AXI - 0.5D0) C(4) = -0.0003Q0 - 0.0012Q0 * (AXI - 0.5D0) C(5) = -0.0002Q0 - 0.002Q0 * (AXI - 0.5D0) C(6) = 0.0001Q0 + 0.0008Q0 * (AXI - 0.5D0) C(7) = 0.0Q0 C(8) = -0.0004Q0 * (AXI - 0.5D0) C(9) = 0.0Q0 C(10) = 0.0001Q0 * (AXI - 0.5D0) ELSE IF (AXI.LT.1.0D0) THEN C(0) = -0.0072Q0 - 0.0248Q0 * (AXI - 0.75D0) C(1) = 0.007Q0 + 0.0192Q0 * (AXI - 0.75D0) C(2) = -0.021Q0 - 0.07Q0 * (AXI - 0.75D0) C(3) = 0.0057Q0 + 0.0172Q0 * (AXI - 0.75D0) C(4) = -0.0006Q0 + 0.0004Q0 * (AXI - 0.75D0) C(5) = -0.0007Q0 - 0.004Q0 * (AXI - 0.75D0) C(6) = 0.0003Q0 + 0.0012Q0 * (AXI - 0.75D0) C(7) = 0.0004Q0 * (AXI - 0.75D0) C(8) = -0.0001Q0 - 0.0008Q0 * (AXI - 0.75D0) C(9) = 0.0004Q0 * (AXI - 0.75D0) C(10) = 0.0001Q0 ELSE IF (AXI.LT.1.25D0) THEN C(0) = -0.0134Q0 - 0.0308Q0 * (AXI - 1.0D0) C(1) = 0.0118Q0 + 0.0248Q0 * (AXI - 1.0D0) C(2) = -0.0385Q0 - 0.0828Q0 * (AXI - 1.0D0) C(3) = 0.01Q0 + 0.0188Q0 * (AXI - 1.0D0) C(4) = -0.0005Q0 + 0.0036Q0 * (AXI - 1.0D0) C(5) = -0.0017Q0 - 0.006Q0 * (AXI - 1.0D0) C(6) = 0.0006Q0 + 0.0016Q0 * (AXI - 1.0D0) C(7) = 0.0001Q0 + 0.0008Q0 * (AXI - 1.0D0) C(8) = -0.0003Q0 - 0.0008Q0 * (AXI - 1.0D0) C(9) = 0.0001Q0 C(10) = 0.0001Q0 + 0.0004Q0 * (AXI - 1.0D0) ELSE IF (AXI.LT.1.5D0) THEN C(0) = -0.0211Q0 - 0.0344Q0 * (AXI - 1.25D0) C(1) = 0.018Q0 + 0.0268Q0 * (AXI - 1.25D0) C(2) = -0.0592Q0 - 0.0896Q0 * (AXI - 1.25D0) C(3) = 0.0147Q0 + 0.0164Q0 * (AXI - 1.25D0) C(4) = 0.0004Q0 + 0.008Q0 * (AXI - 1.25D0) C(5) = -0.0032Q0 - 0.008Q0 * (AXI - 1.25D0) C(6) = 0.001Q0 + 0.0012Q0 * (AXI - 1.25D0) C(7) = 0.0003Q0 + 0.002Q0 * (AXI - 1.25D0) C(8) = -0.0005Q0 - 0.0012Q0 * (AXI - 1.25D0) C(9) = 0.0001Q0 C(10) = 0.0002Q0 + 0.0008Q0 * (AXI - 1.25D0) ELSE IF (AXI.LT.1.75D0) THEN C(0) = -0.0297Q0 - 0.0364Q0 * (AXI - 1.5D0) C(1) = 0.0247Q0 + 0.0248Q0 * (AXI - 1.5D0) C(2) = -0.0816Q0 - 0.0904Q0 * (AXI - 1.5D0) C(3) = 0.0188Q0 + 0.0132Q0 * (AXI - 1.5D0) C(4) = 0.0024Q0 + 0.0116Q0 * (AXI - 1.5D0) C(5) = -0.0052Q0 - 0.0092Q0 * (AXI - 1.5D0) C(6) = 0.0013Q0 + 0.0008Q0 * (AXI - 1.5D0) C(7) = 0.0008Q0 + 0.0028Q0 * (AXI - 1.5D0) C(8) = -0.0008Q0 - 0.0016Q0 * (AXI - 1.5D0) C(9) = 0.0001Q0 - 0.0004Q0 * (AXI - 1.5D0) C(10) = 0.0004Q0 + 0.0012Q0 * (AXI - 1.5D0) ELSE IF (AXI.LT.2.0D0) THEN C(0) = -0.0388Q0 - 0.0372Q0 * (AXI - 1.75D0) C(1) = 0.0309Q0 + 0.02Q0 * (AXI - 1.75D0) C(2) = -0.1042Q0 - 0.0884Q0 * (AXI - 1.75D0) C(3) = 0.0221Q0 + 0.0092Q0 * (AXI - 1.75D0) C(4) = 0.0053Q0 + 0.0152Q0 * (AXI - 1.75D0) C(5) = -0.0075Q0 - 0.0096Q0 * (AXI - 1.75D0) C(6) = 0.0015Q0 C(7) = 0.0015Q0 + 0.004Q0 * (AXI - 1.75D0) C(8) = -0.0012Q0 - 0.0016Q0 * (AXI - 1.75D0) C(9) = -0.0008Q0 * (AXI - 1.75D0) C(10) = 0.0007Q0 + 0.0012Q0 * (AXI - 1.75D0) ELSE IF (AXI.LT.2.25D0) THEN C(0) = -0.0481Q0 - 0.0368Q0 * (AXI - 2.0D0) C(1) = 0.0359Q0 + 0.0168Q0 * (AXI - 2.0D0) C(2) = -0.1263Q0 - 0.0844Q0 * (AXI - 2.0D0) C(3) = 0.0244Q0 + 0.0044Q0 * (AXI - 2.0D0) C(4) = 0.0091Q0 + 0.0184Q0 * (AXI - 2.0D0) C(5) = -0.0099Q0 - 0.0088Q0 * (AXI - 2.0D0) C(6) = 0.0015Q0 - 0.0016Q0 * (AXI - 2.0D0) C(7) = 0.0025Q0 + 0.0044Q0 * (AXI - 2.0D0) C(8) = -0.0016Q0 - 0.0012Q0 * (AXI - 2.0D0) C(9) = -0.0002Q0 - 0.0016Q0 * (AXI - 2.0D0) C(10) = 0.001Q0 + 0.0012Q0 * (AXI - 2.0D0) ELSE IF (AXI.LT.2.5D0) THEN C(0) = -0.0573Q0 - 0.0368Q0 * (AXI - 2.25D0) C(1) = 0.0401Q0 + 0.0136Q0 * (AXI - 2.25D0) C(2) = -0.1474Q0 - 0.08Q0 * (AXI - 2.25D0) C(3) = 0.0255Q0 + 0.0012Q0 * (AXI - 2.25D0) C(4) = 0.0137Q0 + 0.02Q0 * (AXI - 2.25D0) C(5) = -0.0121Q0 - 0.008Q0 * (AXI - 2.25D0) C(6) = 0.0011Q0 - 0.0028Q0 * (AXI - 2.25D0) C(7) = 0.0036Q0 + 0.0048Q0 * (AXI - 2.25D0) C(8) = -0.0019Q0 - 0.0008Q0 * (AXI - 2.25D0) C(9) = -0.0006Q0 - 0.002Q0 * (AXI - 2.25D0) C(10) = 0.0013Q0 + 0.0016Q0 * (AXI - 2.25D0) ELSE IF (AXI.LT.2.75D0) THEN C(0) = -0.0665Q0 - 0.036Q0 * (AXI - 2.5D0) C(1) = 0.0435Q0 + 0.012Q0 * (AXI - 2.5D0) C(2) = -0.1674Q0 - 0.0756Q0 * (AXI - 2.5D0) C(3) = 0.0258Q0 - 0.0028Q0 * (AXI - 2.5D0) C(4) = 0.0187Q0 + 0.022Q0 * (AXI - 2.5D0) C(5) = -0.0141Q0 - 0.0064Q0 * (AXI - 2.5D0) C(6) = 0.0004Q0 - 0.0044Q0 * (AXI - 2.5D0) C(7) = 0.0048Q0 + 0.0052Q0 * (AXI - 2.5D0) C(8) = -0.0021Q0 C(9) = -0.0011Q0 - 0.0024Q0 * (AXI - 2.5D0) C(10) = 0.0017Q0 + 0.0016Q0 * (AXI - 2.5D0) ELSE IF (AXI.LT.3.0D0) THEN C(0) = -0.0755Q0 - 0.0352Q0 * (AXI - 2.75D0) C(1) = 0.0465Q0 + 0.0028Q0 * (AXI - 2.75D0) C(2) = -0.1863Q0 - 0.0708Q0 * (AXI - 2.75D0) C(3) = 0.0251Q0 - 0.0044Q0 * (AXI - 2.75D0) C(4) = 0.0242Q0 + 0.0228Q0 * (AXI - 2.75D0) C(5) = -0.0157Q0 - 0.0044Q0 * (AXI - 2.75D0) C(6) = -0.0007Q0 - 0.0056Q0 * (AXI - 2.75D0) C(7) = 0.0061Q0 + 0.0048Q0 * (AXI - 2.75D0) C(8) = -0.0021Q0 + 0.0004Q0 * (AXI - 2.75D0) C(9) = -0.0017Q0 - 0.0032Q0 * (AXI - 2.75D0) C(10) = 0.0021Q0 + 0.0012Q0 * (AXI - 2.75D0) ELSE IF (AXI.LT.3.25D0) THEN C(0) = -0.0843Q0 - 0.0348Q0 * (AXI - 3.0D0) C(1) = 0.0472Q0 + 0.006Q0 * (AXI - 3.0D0) C(2) = -0.204Q0 - 0.0668Q0 * (AXI - 3.0D0) C(3) = 0.024Q0 - 0.0072Q0 * (AXI - 3.0D0) C(4) = 0.0299Q0 + 0.0236Q0 * (AXI - 3.0D0) C(5) = -0.0168Q0 - 0.0028Q0 * (AXI - 3.0D0) C(6) = -0.0021Q0 - 0.0064Q0 * (AXI - 3.0D0) C(7) = 0.0073Q0 + 0.0044Q0 * (AXI - 3.0D0) C(8) = -0.002Q0 + 0.0016Q0 * (AXI - 3.0D0) C(9) = -0.0025Q0 - 0.0036Q0 * (AXI - 3.0D0) C(10) = 0.0024Q0 + 0.0012Q0 * (AXI - 3.0D0) ELSE IF (AXI.LT.3.5D0) THEN C(0) = -0.093Q0 - 0.0336Q0 * (AXI - 3.25D0) C(1) = 0.0487Q0 + 0.002Q0 * (AXI - 3.25D0) C(2) = -0.2207Q0 - 0.0628Q0 * (AXI - 3.25D0) C(3) = 0.0222Q0 - 0.0092Q0 * (AXI - 3.25D0) C(4) = 0.0358Q0 + 0.0244Q0 * (AXI - 3.25D0) C(5) = -0.0175Q0 - 0.0012Q0 * (AXI - 3.25D0) C(6) = -0.0037Q0 - 0.0076Q0 * (AXI - 3.25D0) C(7) = 0.0084Q0 + 0.0036Q0 * (AXI - 3.25D0) C(8) = -0.0016Q0 + 0.0016Q0 * (AXI - 3.25D0) C(9) = -0.0034Q0 - 0.0036Q0 * (AXI - 3.25D0) C(10) = 0.0027Q0 + 0.0012Q0 * (AXI - 3.25D0) ELSE IF (AXI.LT.4.0D0) THEN C(0) = -0.1014Q0 - 0.0346Q0 * (AXI - 3.5D0) C(1) = 0.0492Q0 - 0.002Q0 * (AXI - 3.5D0) C(2) = -0.2364Q0 - 0.0572Q0 * (AXI - 3.5D0) C(3) = 0.0199Q0 - 0.0102Q0 * (AXI - 3.5D0) C(4) = 0.0419Q0 + 0.0248Q0 * (AXI - 3.5D0) C(5) = -0.0178Q0 + 0.0014Q0 * (AXI - 3.5D0) C(6) = -0.0056Q0 - 0.0088Q0 * (AXI - 3.5D0) C(7) = 0.0093Q0 + 0.0028Q0 * (AXI - 3.5D0) C(8) = -0.0012Q0 + 0.0028Q0 * (AXI - 3.5D0) C(9) = -0.0043Q0 - 0.0042Q0 * (AXI - 3.5D0) C(10) = 0.003Q0 + 0.0004Q0 * (AXI - 3.5D0) ELSE IF (AXI.LT.4.5D0) THEN C(0) = -0.1187Q0 - 0.0282Q0 * (AXI - 4.0D0) C(1) = 0.0482Q0 - 0.0158Q0 * (AXI - 4.0D0) C(2) = -0.265Q0 - 0.0498Q0 * (AXI - 4.0D0) C(3) = 0.0148Q0 - 0.0084Q0 * (AXI - 4.0D0) C(4) = 0.0543Q0 + 0.0238Q0 * (AXI - 4.0D0) C(5) = -0.0171Q0 + 0.0036Q0 * (AXI - 4.0D0) C(6) = -0.01Q0 - 0.0092Q0 * (AXI - 4.0D0) C(7) = 0.0107Q0 + 0.0008Q0 * (AXI - 4.0D0) C(8) = 0.0002Q0 + 0.0032Q0 * (AXI - 4.0D0) C(9) = -0.0064Q0 - 0.0034Q0 * (AXI - 4.0D0) C(10) = 0.0032Q0 - 0.0002Q0 * (AXI - 4.0D0) ELSE C(0) = -0.1328Q0 C(1) = 0.0403Q0 C(2) = -0.2899Q0 C(3) = 0.0106Q0 C(4) = 0.0662Q0 C(5) = -0.0153Q0 C(6) = -0.0146Q0 C(7) = 0.0111Q0 C(8) = 0.0018Q0 C(9) = -0.0081Q0 C(10) = 0.0031Q0 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