SUBROUTINE DA01A(Y,E,A,B,C,D,N0,X0,DX,MOP) C RUNGE KUTTA MERSON (K ENTRIES TO DYBDX)WITH EST.TRUNC.ERROR E(1) C CALLS S/R DYBDX(Y,F,N,X,) TO SET UP F(4)=FUNC(Y(I),X),I,J=1,N IMPLICIT REAL*8(A-H,O-Z) DIMENSION Y(10),E(10),A(10),B(10),C(10),D(10) N=N0 Z=X0 H=DX/3. C WRITE(6,50)X0 50 FORMAT(1H ,'T4= ',D13.4) 51 FORMAT(1H ,'T5= ',D13.4) C DO 1 I=1,N D(I)=Y(I) 1 CONTINUE C X=Z DO 41 J=1,5 CALL DYBDX(Y,E,N,MOP,X) C DO 21 I=1,N GO TO (11,12,13,14,15),J 11 A(I)=H*E(I) Y(I)=D(I)+A(I) GO TO 21 12 B(I)=H*E(I) Y(I)=D(I)+(A(I)+B(I))*0.5 GO TO 21 13 B(I)=H*E(I) Y(I)=D(I)+(A(I)+B(I)*3.)*0.375 GO TO 21 14 C(I)=H*E(I) Y(I)=D(I)+(A(I)-B(I)*3.+C(I)*4.)*1.5 GO TO 21 15 Z =D(I) D(I)=H*E(I) Y(I)=Z+(A(I)+C(I)*4.+D(I))*0.5 E(I)=(A(I)*2.-B(I)*9.+C(I)*8.-D(I))/10. 21 CONTINUE C GO TO (31,41,33,34,41),J 31 X=Z+H GO TO 41 33 X=Z+DX*0.5 GO TO 41 34 X=Z+DX 41 CONTINUE X0=X C WRITE(6,51)X0 RETURN END C DYBDX FOR RUNGE KUTTA ROUTINES SUBROUTINE DYBDX(Y,F,N,M,T) REAL*8 Y,F,T,AK,X COMMON /b_2/AMPSSP(10),AMP(10),OHM(10),RT(10), 1 FK(10,10),RK(10,10) COMMON /b_3/AMPS(10),VELZ(10),P1(10),P2(10),P3(10),P4(10), 1 Z(10),VOL(10,800),IND(10),DT COMMON /b_5/TOTVOL(10),ARCON(10),PROVEL(6,5),ALPHA(10) DIMENSION Y(10),F(10),AK(10,10),X(10) MA=M*N I=1 DO 1 J=1,N,1 DO 1 K=1,M,1 AK(K,J)=DBLE(RK(K,J)) 1 I=I+1 DO 15 I=1,N,1 15 X(I)=-(Y(I)*DBLE(RT(I))+(Y(I)-DBLE(AMPSSP(I)))*DBLE(OHM(I))) DO 11 I=1,N,1 F(I)=0.0D0 DO 11 J=1,N,1 11 F(I)=F(I)+AK(I,J)*X(J) RETURN END C CCC C FL=KEMAINDR.QUENCH.DEKM(NO) C CC C MODIFIED CULHAM QUENCH INCLUDING RUNGE KUTTA C SUBROUTINE DEKM(YY,ES,TIME,DT,NC,A1) REAL*8 Y,E,AA,BB,CC,DD,T,DTO,A2,T2 COMMON/ERR/EB1,EB2 DIMENSION YY(10),ES(10),DY(10),A1(10) DIMENSION Y(10),E(10),AA(10),BB(10),CC(10),DD(10),A2(10) DO 1 I=1,NC,1 A2(I)=DBLE(A1(I)) 1 Y(I)=DBLE(YY(I)) T=DBLE(TIME) DTO=DBLE(DT) N=NC M=NC KC=1 20 AM=0.0 IO=0 DO 30 KL=1,N,1 AM=AM+Y(KL) 30 A2(KL)=Y(KL) AM=AM/FLOAT(NC) C WRITE(6,33)TIME 33 FORMAT(1H ,'T1= ',F10.5) 31 FORMAT(1H ,'T2= ',F10.5) 32 FORMAT(1H ,'T3= ',F10.5) CALL DA01A(Y,E,AA,BB,CC,DD,N,T,DTO,M) C WRITE(6,31)TIME DO 10 IND=1,NC,1 1001 FORMAT(3D13.6,1E13.6,2D13.6,3I10) C WRITE(6,35)Y(IND) 35 FORMAT(1H ,'Y(IND)= ',D13.4) IF(DABS(Y(IND)).LT.0.1) GO TO 10 E(IND)=DABS(E(IND)/Y(IND)) C WRITE(6,34)E(IND) 34 FORMAT(1H ,'E(IND)= ',D13.4) IF(E(IND).GT.ABS(EB1)) GO TO 200 IF(E(IND).LT.ABS(EB2)) IO=IO+1 702 FORMAT(3X,3F10.4,2I10) 10 CONTINUE IF(IO.EQ.NC) DTO=DTO*2.0D0 C WRITE(6,36)DTO 36 FORMAT(1H ,'DTO= ',D13.4) GO TO 400 200 T=T-DTO DTO=DTO/2.0D0 KC=KC+1 C WRITE(6,37)T,DTO 37 FORMAT(1H ,'T,DTO= ',2D13.4) IF(KC.LE.3) GO TO 17 KC=1 WRITE(6,656)(E(IND),IND=1,NC) WRITE(6,23)DTO DTO=2.0D0*DTO T=T+DTO GO TO 400 17 DO 14 KL=1,NC 14 Y(KL)=A2(KL) GO TO 20 C 400 T2=TIME+DT *** WARNING NDG ON VAX *** 400 T2=DBLE(TIME)+DBLE(DT) 701 FORMAT(3F10.4,1I10) C DTEST=1.-SNGL(T/DABS(T2)) C IF(ABS(DTEST).LE..00001)DT=SNGL(T)-TIME IF(T.GE.DABS(T2)) DT=SNGL(T)-TIME KC=1 C WRITE(6,39)T2,T,DT 39 FORMAT(1H ,'T2,T,DT= ',2D13.4,E13.4) IF(T.GE.DABS(T2)) GO TO 500 C IF(ABS(DTEST).LE..00001)GO TO 500 DO 42 I=1,N,1 42 A2(I)=Y(I) GO TO 20 500 DO 40 KC=1,NC,1 40 YY(KC)=SNGL(Y(KC)) TIME=TIME+DT C WRITE(6,32)TIME RETURN 600 CONTINUE WRITE(6,601) 601 FORMAT(2X,' DT IN DEKM CANNOT BE MATCHED TO TIME INTVAL') GO TO 500 650 RETURN 656 FORMAT(' ERRORS='10G11.3) 23 FORMAT(' INTERVAL TOO SMALL:'G15.6) END