PROGRAM EMP
        CHARACTER*10 FIL1,KWORD
        DIMENSION A(10,500),X(10,3)  !定義兩個數(shù)組
        REAL X
        !TNT的參數(shù)
        VD=6930 !m/s
          DEN=1.630 !g/cm3
!          PH=21*10E9 !pa
      !  write(*,*)' AN or other EXPLOSIVE? (AN=1)'
   !     r=(vd*vd/1608/1608+1.)**.5 !銨油多方指數(shù)的計算
        

        !DEN,TIE分別為炸藥的密度和厚度
          DEN=1.63
          TIE=20
      
        !DEP,TIP分別為平板的密度和厚度
          DEP=1.7
          TIP=8

        !Q表示質(zhì)量比
        Q=(DEN*TIE)/DEP/TIP

          !R即為等熵指數(shù)K 單位分別為kg/m3,m/s,pa
          !或者為g/cm3,km/s,Gpa
          R=2.72
        !R=DEN0*VD**2/PH-1.0

        !角度
          CTA=30
        
          !插值點數(shù)
          N=20

          EPS=0.001
        WRITE(*,3)Q,R,TIE !Q質(zhì)量比R多方指數(shù)及GUMMA
        WRITE(*,7)VD,DEN   !炸藥的爆速和密度
7       FORMAT(6H   Vd=,F7.2,5H(m/s),10H DENSITYe=,F5.3,9H(g/cm**3))
3       FORMAT(4H  R=,F7.4,9H   GAMMA=,F5.3,22H   Explosive thickness,
     & F8.4,4H(mm))
        WRITE(*,5)
5       FORMAT(' Angle(DEG)    Vp(m/s)      H(mm)     P(GPa)     X(mm)')
        DEM=DEN/2.
          SS=1.0172
  !      CALL MMM(DEM,VD,R,EPS,CTA,SS,CTAB)        !CTA(指的是角度)子例行程序
        CALL PDENFRI(N,2.2292,4.2945,R,A,0.,0.,1.)  !務(wù)必注意這兩步的作用
        DO 10 I=N,1,-1
        DO 10 J=1,10
10      A(J,I+1)=A(J,I)
        DO 20 J=1,10
20      A(J,1)=0.
        A(7,1)=3.1415926/2.
        DD=-180./3.1415926*A(1,1)                        !角度轉(zhuǎn)化
        DV=-VD*2.*SIN(A(1,1)/2.)                         !泰勒公式計算的拋擲速度
        DHH=-TIE*A(9,1)                                  !炸藥的厚度乘以角度
        PH=DEN*(VD/1000.)**2/(R+1.)                      !爆轟波頭的壓力
        WRITE(*,4)DD,DV,DHH,PH,A(8,1)*Tie
4       FORMAT(10(1X,F9.4))
200     DO 30 J=1,10
        X(J,1)=A(J,1)
30      X(J,2)=A(J,2)
        CALL RIM(X,R,-1,Q,EPS,P) !調(diào)用邊界插值函數(shù)
        DO 40 I=1,10
40      A(I,1)=X(I,3)
        DD=-180./3.1415926*A(1,1)
        DV=-VD*2.*SIN(A(1,1)/2.)
        DHH=-TIE*A(9,1)
!          WRITE(*,4)P P的數(shù)值
        WRITE(*,4)DD,DV,DHH,P*PH,A(8,1)*Tie
        N=N-1
        DO 50 I=1,N
        DO 60 J=1,10
        X(J,1)=A(J,I)
60      X(J,2)=A(J,I+2)
        CALL INSIDE(X,EPS)  !調(diào)用內(nèi)部插值函數(shù)
        DO 70 J=1,10
70      A(J,I+1)=X(J,3)
50      CONTINUE
        IF(N.EQ.1) GOTO 100
        GOTO 200
100     WRITE(*,*)'  Anther Calculation(Y/N)?'
!       READ(*,'(A)')KWORD
!        IF(KWORD.EQ.'Y') GOTO 909
!        WRITE(*,77)FIL1
        STOP

       
        END

C      plane-fluid subrountine
      !*******************************
      !等熵的壓力關(guān)系
        !*******************************

        FUNCTION PS(X)!密度我變量
        V=1.63/X                                                                                                                                                   
        A=373.8 !A,B,C GPA
        B=2.747
        CC=0.734
        R1=4.15
        R2=0.90
        W=0.3
        PS=(A*EXP(-R1*V)+B*EXP(-R2*V)+CC/(V**(W+1)))*1.0
        RETURN
        END

        !*******************************
      !等熵的比內(nèi)能函數(shù)
        !*******************************

        FUNCTION ES(X) !X指的是密度
        V=1.63/X
        A=373.8 !A,B,C GPA
        B=2.747
        CC=0.734
        R1=4.15
        R2=0.90
        W=0.3
        ES=A*EXP(-R1*V)/R1+B*EXP(-R2*V)/R2+CC/(V**W)/W
        RETURN
        END

        !*******************************
      !等熵的聲速函數(shù)
        !*******************************

      !X指的是密度
        FUNCTION CS(X)! 聲速的函數(shù)     
        V=1.63/X
        A=373.8 !A,B,C GPA
        B=2.747
        CC=0.734
        R1=4.15
        R2=0.90
        W=0.3
        DPV=-A*R1*EXP(-R1*V)-B*R2*EXP(-R2*V)-(W+1)*CC/V**(W+1)
        CS2=-V**2/X*DPV
        CS=SQRT(CS2)
        RETURN
        END
  
      !*******************************
        !普朗特-邁耶函數(shù)
      !*******************************

      !X指的是密度
         REAL FUNCTION F(X) !普朗特函數(shù)需要修改CS時含有變量的函數(shù)
        REAL X
!         XI=ES(X)+PS(X)/X
  !     XI0=ES(DENH)+PS(DENH)/DENH+0.5*qH**2
         XX=LOG(SQRT(30-XI(X)))!XI0需要帶入
         F=-SQRT((2.0*EXP(2.0*XX)-CS(X)**2)/CS(X)**2)
         END


       !*******************************
         !定義了焓的函數(shù)
       !*******************************

         FUNCTION XI(X)
       XI=ES(X)+PS(X)/X
         END

       !*******************************
       !辛普生積分函數(shù)
       !*******************************

       FUNCTION SIMPSON(F,A,B,N)
         H=(B-A)/(2.0*N)
         S=F(A)-F(B)
         DO I=1,N
          S1=F(A+(2.0*I-1)*H)
          S2=F(A+2.0*I*H)
          S=S+4.0*S1+2.0*S2
         END DO
         SIMPSON=H/3.0*S
        END

      !*********************************
        !利用辛普生積分求得馬赫角
        !*********************************

      FUNCTION V(X)
        REAL X,V
        EXTERNAL F,XI
        V=SIMPSON(F(X),30.0,XI(X),20)
        RETURN
        END

      !*******************************
      !普邁函數(shù)的反函數(shù)
        !*******************************

        FUNCTION VN(Q,EPS) !反函數(shù)!Q角度R等熵指數(shù)需要修改
        Y0=0.
        Y1=2.
        IF(V(Y1).LE.Q)THEN
        Y0=Y1
        Y1=5.
        END IF
        IF(V(Y1).LE.Q)THEN
        Y0=Y1
        Y1=10.
        END IF
        IF(V(Y1).LE.Q)THEN
        Y0=Y1
        Y1=15.
        END IF
        IF(V(Y1).LE.Q)THEN
        Y0=Y1
        Y1=20.
        END IF
        IF(V(Y1).LE.Q)THEN
        Y0=Y1
        Y1=30.
        END IF
1        YY=(Y0+Y1)/2.
        IF(V(YY).LE.Q)THEN
        Y0=YY   
        ELSE
        Y1=YY       !VN得出來的是a=log(i0-i)
        END IF
        IF(ABS(Y1-Y0).LT.(.01*EPS)) GOTO 2
        GOTO 1
2        VN=YY
        RETURN
        END

      !*******************************
      !正弦函數(shù)的反函數(shù)
        !*******************************

        FUNCTION ASINN(X)
           ASINN=ATAN(1./X)
        RETURN
        END

      !*******************************
      !焓的反函數(shù)
        !*******************************
       FUNCTION VI(Q,EPS)
          Y0=0.
        Y1=2.
        IF(XI(Y1).LE.Q)THEN
        Y0=Y1
        Y1=5.
        END IF
        IF(XI(Y1).LE.Q)THEN
        Y0=Y1
        Y1=10.
        END IF
        IF(XI(Y1).LE.Q)THEN
        Y0=Y1
        Y1=15.
        END IF
        IF(XI(Y1).LE.Q)THEN
        Y0=Y1
        Y1=20.
        END IF
        IF(XI(Y1).LE.Q)THEN
        Y0=Y1
        Y1=30.
        END IF
3        YY=(Y0+Y1)/2.
        IF(XI(YY).LE.Q)THEN
        Y0=YY   
        ELSE
        Y1=YY       !VN得出來的是a=log(i0-i)
        END IF
        IF(ABS(Y1-Y0).LT.(.01*EPS)) GOTO 4
        GOTO 3
4        VI=YY
        RETURN
        END

       !*******************************
       !內(nèi)部插值函數(shù) 已修改 還需要修改焓0，VN
       !*******************************

        SUBROUTINE INSIDE(X,EPS)
        DIMENSION X(10,3)
        X(1,3)=(X(10,2)-X(10,1)+X(1,1)+X(1,2))/2.
        X(10,3)=X(1,3)-X(1,1)+X(10,1)
        X(5,3)=VN(X(10,3),EPS)
          X(3,3)=30.0-EXP(2.0*X(5,3))!X(3,0)指的是I0
          X(2,3)=VI(X(3,3),EPS)
          X(6,3)=SQRT(2.0*EXP(2.0*X(5,3))/CS(X(2,3))**2)
        X(7,3)=ASINN(X(6,3))
        A=1./TAN(X(1,1)+X(7,1))
        B=1./TAN(X(1,2)-X(7,2))
        X(9,3)=(X(8,2)-X(8,1)+A*X(4,1)-B*X(9,2))/(A-B)
        X(8,3)=X(8,1)+A*(X(9,3)-X(9,1))
        RETURN
        END

        !*******************************
        !邊界插值函數(shù)
        !*******************************

        SUBROUTINE RIM(X,R,K,Q,EPS,PRESS)
        DIMENSION X(10,3)
        A=1./TAN(X(1,2)+K*X(7,2))
        F=X(1,2)+K*X(7,2)
        PK1=PS(X(2,1))*Q/(R+1.)
        PK=PK1
        DO 10 I=1,2
        SINF=SIN(X(1,1)-F)
       X(1,3)=X(1,1)-K*(X(3,2)-X(3,1)-A*(X(4,2)-X(4,1)))*SIN(F)/SINF
     & *PK
        X(10,3)=FLOAT(K)*(X(1,3)-X(1,1))+X(10,2)
        X(5,3)=VN(X(10,3),EPS)
          X(3,3)=30.0-EXP(2.0*X(5,3)) !30.0=XI0
          X(2,3)=VI(X(3,3),EPS)            !VI是焓的反函數(shù)
          X(6,3)=SQRT(2.0*EXP(2.0*X(5,3))/(CS(X(2,3))**2))
        X(7,3)=ASINN(X(6,3))

        PK3=PK1+(X(3,3)-X(3,1))/0.5/(X(2,1)+X(2,3))*Q/(R+1.)
        PK=.5*PK1+.5*PK3
10        CONTINUE
        X(8,3)=(SIN(X(1,3))-SIN(X(1,1)))*K/PK+X(8,1)
        X(9,3)=(COS(X(1,1))-COS(X(1,3)))*K/PK+X(9,1)
        PRESS=PK/(Q/(R+1.))!計算P的表達式
        RETURN
        END

          !*******************************
        !密度流初值選取 已修改
          !*******************************

        SUBROUTINE PDENFRI(N,DEN0,DENM,R,X,Q0,X0,Y0)
        DIMENSION X(10,N)
        DO 1 I=1,N
        X(2,I)=DEN0+(DENM-DEN0)/(N-1)*(I-1)
        IF(X(2,I).LE..1E-6)THEN
        X(2,I)=DEN0+0.1E-6
        END IF
          X(3,I)=ES(X(2,I))+PS(X(2,I))/X(2,I)
        X(4,I)=VI(X(3,I),EPS)
          X(5,I)=LOG(SQRT(30.0-X(3,I)))
          X(6,I)=SQRT(2.0*EXP(2.0*X(5,I)))
          X(7,I)=ASIN(X(5,I))
          X(8,I)=X0
          X(9,I)=Y0
!          X(10,I)=V(X(5,1))
          X(1,I)=V(X(5,I))+Q0
1        CONTINUE
        RETURN
        END

Originally posted by 20033657 at 2010-10-11 11:18:17:
In the call to SIMPSON, actual argument #1 does not match the type and kind of the corresponding dummy argument.什么原因
程序如下
        PROGRAM EMP
        CHARACTER*10 FIL1,KWORD
        DIM ...