C PROGRAM OGNCG6 C C REC1G08 FROM /MHD2A/REC1G06 C COMPUTER SIMULATION OF EQUILIBRIUM C TWO DIMENSIONAL MHD MODEL C CARTESIAN COORDINATE FINITE RESISTIVITY C 1985.05.24 BY TATSUKI OGINO PARAMETER (NR=400,NZ=100,N1=NR+2,N2=N1*(NZ+2)) PARAMETER (NB=6,NBB=7,N3=N2*NB,N4=N2*NBB,N5=N2*4) PARAMETER (LAST=4,LLL=6,IBIG=0) PARAMETER (IIQ0=4,IIP0=4,IIS0=4,THR=8.00) PARAMETER (RMU0=0.010,EUD0=0.00,RO01=5.0E-4,PR01=3.56E-8) PARAMETER (NRP=100,ARU=10.0,ARU1=10.0,PMU0=1.0,DFU0=1.0) PARAMETER (EAT0=0.025,RAT1=0.2,RATT=0.2) PARAMETER (NIN=1,ITAP=11) PARAMETER (K2=N1*2,K3=N1*3,K4=N1*4,K5=N1*5,K6=N1*6,K7=N1*7) PARAMETER (KK2=K3*2,KK3=K3*3,KK4=K3*4,KK5=K3*5,KK6=K3*6) PARAMETER (N11=N1+1,NR1=NR+1,NR2=NR+2,NZ1=NZ+1,NZ2=NZ+2) PARAMETER (M2=N2*2,M3=N2*3,M4=N2*4,M5=N2*5) PARAMETER (NRZ=NR*NZ,MN1=N1*NZ1) DOUBLE PRECISION IT1,IT2,ITD,IT0 DIMENSION F(N3),FF(KK6),U(KK6),P(K7),ZMAX(LLL),ZMIN(LLL), 1 PP(N5),CJ(10),CP(10),FBB(6),HHX(N1) C M 1-RO 2-VR 3-VZ 4-P 5-BR 6-BZ 7-JY 8-JYZ C C CJ R1 DR DZ B0 R11 R22 AID AR1 GAM GRA DATA CJ/10.0,0.3,0.5,0.0,40.0,80.0,0.7,0.8,0.9,1.0/ C C CP R0 Z0 RA VR0 P0 GRA 07 VSW IBR IBZ C DATA CP/1203.,303.,16.,.044,2.68,1.35,7.0,.044,9.0,-0.0/ DATA CP/1203.,303.,16.,.044,2.68,1.35,7.0,.044,0.0,-1.5/ C REWIND 12 PI=3.1415926 GAM=2.0 GM2=0.5*(GAM-1.0) GRA=CP(6)*1.0E-6 PO0=CP(5)*1.0E-7 BIS=CP(10)*1.0E-4 CP(10)=0.0 CJ(2)=RATT VSW=CP(8) AR1=CP(3) AR2=AR1*AR1 HR=CP(1)/FLOAT(NR1) HZ=CP(2)/FLOAT(NZ1) T=0.5*HR*THR T1=0.5*T CJ(8)=AR1 CJ(9)=GAM CJ(10)=GRA C DR1=0.5*T1/HR DZ1=0.5*T1/HZ DR2=0.5*T/HR DZ2=0.5*T/HZ DR3=T1/(HR*HR) DZ3=T1/(HZ*HZ) DR4=T/(HR*HR) DZ4=T/(HZ*HZ) C RMU=RMU0*RO01 PMU=RMU0*PMU0 DFU=RMU0*DFU0 EUD=EUD0 RO02=RAT1*RO01 PR02=RAT1*PR01 C WRITE (6,12) LAST,NR,NZ,N1,N2,N3,N4,RMU,EUD, 1 HR,HZ,T,T1,RO01,PR01,GRA,DR2,DZ2,DR4,DZ4,EAT0, 2 (CP(I),I=1,10),(CJ(J),J=1,10) 12 FORMAT(1H ,5X,7I10/(1H ,5X,1P8E15.5)) C DO 20 I=1,N3 20 F(I)=0.0 DO 221 I=1,N1 221 HHX(I)=HR*FLOAT(I) DO 222 I=1,NB 222 FBB(I)=1.0 FBB(3)=-1.0 FBB(5)=-1.0 C C INITIAL CONDITION CALL EQUIP4(NR,NZ,NRP,RO01,PR01,CJ,CP,F,PP) GO TO 882 WRITE(6,212) (F(I),I=1,N3) 212 FORMAT(1H ,2X,1P10E12.4) 882 CONTINUE C JXG=N2 JYG=NZ2 JXN=NR2 LANK=20 IP=0 LB=1 LL=6 I1=II C CALL SCALM(NR,NZ,NRP,LB,LL,LLL,HR,HZ,AR2,ZMIN,ZMAX,F) DO 130 J=1,2 LB=1+3*(J-1) LL=3 C CALL GRAPE(JXG,JYG,LB,LL,ZMIN,ZMAX,I1,LLL,JXN,LANK,IP,F) 130 CONTINUE C C INITIAL CONDITION C IF(ITAP.EQ.999) GO TO 400 DO 410 I=1,ITAP C READ(13) F C IF(I.EQ.1) WRITE(12) F 410 CONTINUE 400 CONTINUE C TWO DIMENSIONAL CARTESIAN MODEL TIME=0.0 IIQ=0 IIP=0 IIS=0 C C C IT1=0.0D0 CALL CLOCK(IT0) C DO 100 II=1,LAST CALL CLOCK(IT2) IT2=IT2-IT0 ITD=IT2-IT1 WRITE(6,662) II,LAST,IT1,IT2,ITD 662 FORMAT(1H ,2I10,3F12.5) IT1=IT2 C IIQ=IIQ+1 IIP=IIP+1 IIS=IIS+1 C BOUNDARY CONDITION AT NZ=1 AND NZ=NZ2 DO 30 M=1,NB I3=N1*2 L1=N2*(M-1) I1=2+L1 I2=I1+MN1 P(2)=(4.0*F(I2-N1)-F(I2-I3))/3.0 U(2)=F(I1+N1)*FBB(M) DO 31 I=3,NR1 I1=I+L1 I2=I1+MN1 P(I)=(4.0*F(I2-N1-1)-F(I2-I3-2))/3.0 31 U(I)=F(I1+N1)*FBB(M) DO 32 I=2,NR1 I2=I+L1+MN1 F(I2)=P(I) 32 F(I2-MN1)=U(I) C BOUNDARY CONDITION AT NR=1 C FREE BOUNDARY AT NR=NR2 I1=N2*(M-1) DO 33 J=1,NZ2 I1=I1+N1 33 P(J)=(4.0*F(I1-1)-F(I1-2))/3.0 I1=N2*(M-1) DO 34 J=1,NZ2 I1=I1+N1 34 F(I1)=P(J) 30 CONTINUE DO 35 I=1,N2 F(I)=AMAX1(F(I),RO02) 35 CONTINUE C IF(II.NE.1) GO TO 136 IF(NIN.EQ.0) WRITE(12) F JXG=N2 JYG=NZ2 JXN=NR2 LANK=20 IP=0 LB=1 LL=6 I1=II C CALL SCALM(NR,NZ,NRP,LB,LL,LLL,HR,HZ,AR2,ZMIN,ZMAX,F) DO 134 J=1,2 LB=1+3*(J-1) LL=3 C CALL GRAPE(JXG,JYG,LB,LL,ZMIN,ZMAX,I1,LLL,JXN,LANK,IP,F) 134 CONTINUE 136 CONTINUE C C STEP OF J=1 J=1 C CURRENT L1=N1*(J-1) DO 381 I=1,NR1 I1=I+L1 I5=I1+M4 I6=I1+M5 J7=I+K6 P(J7)=0.5*((F(I5+N11)-F(I5+1)+F(I5+N1)-F(I5))/HZ 1 -(F(I6+N11)+F(I6+1)-F(I6+N1)-F(I6))/HR) C 2 -PP(I1) 381 CONTINUE DO 382 M=1,NB L4=L1+N2*(M-1) L2=N1*(M-1) L3=K3*(M-1) DO 382 I=1,NR1 I1=I+L4 J1=I+L2 J2=I+L3 P(J1)=0.25*(F(I1)+F(I1+1)+F(I1+N1)+F(I1+N11)) U(J2)=P(J1) 382 CONTINUE C C ZEROTH STEP XX4=0.5*HR*FLOAT(2*NRP-NR1-1) Z=0.5*HZ*FLOAT(2*J-2) DO 501 I=1,NR1 R=HHX(I)+XX4 RO2=R*R+Z*Z RO2=AMAX1(RO2,AR2) I1=I+L1 I2=I1+N2 I3=I2+N2 I4=I3+N2 I5=I4+N2 I6=I5+N2 J1=I J2=J1+N1 J3=J2+N1 J4=J3+N1 J5=J4+N1 J6=J5+N1 J7=J6+N1 JJ1=I X2=RO2/AR2-1.0 X2=ARU*X2*X2 X2=X2/(X2+1.0) C U(JJ1+KK5)=U(JJ1+KK5) 1 +DR1*(F(I3+N11)*F(I5+N11)-F(I2+N11)*F(I6+N11) 2 +F(I3+1)*F(I5+1)-F(I2+1)*F(I6+1) 3 -F(I3+N1)*F(I5+N1)+F(I2+N1)*F(I6+N1) 4 -F(I3)*F(I5)+F(I2)*F(I6)) U(JJ1+KK4)=U(JJ1+KK4) 1 -DZ1*(F(I3+N11)*F(I5+N11)-F(I2+N11)*F(I6+N11) 2 -F(I3+1)*F(I5+1)+F(I2+1)*F(I6+1) 3 +F(I3+N1)*F(I5+N1)-F(I2+N1)*F(I6+N1) 4 -F(I3)*F(I5)+F(I2)*F(I6)) X3=DR1*(F(I2+N11)+F(I2+1)-F(I2+N1)-F(I2)) 1 +DZ1*(F(I3+N11)-F(I3+1)+F(I3+N1)-F(I3)) U(JJ1+KK3)=U(JJ1+KK3)-P(J4)*GAM*X3 1 -DR1*P(J2)*(F(I4+N11)+F(I4+1)-F(I4+N1)-F(I4)) 2 -DZ1*P(J3)*(F(I4+N11)-F(I4+1)+F(I4+N1)-F(I4)) U(JJ1+KK2)=U(JJ1+KK2)+T1*((-P(J7)*P(J5))/P(J1) 1 -EUD*P(J3)-GRA*Z/RO2) 2 -DR1*P(J2)*(F(I3+N11)+F(I3+1)-F(I3+N1)-F(I3)) 3 -DZ1*P(J3)*(F(I3+N11)-F(I3+1)+F(I3+N1)-F(I3)) 4 -DZ1*(F(I4+N11)-F(I4+1)+F(I4+N1)-F(I4))/P(J1) U(JJ1+K3)=U(JJ1+K3)+T1*((P(J7)*P(J6))/P(J1) 1 -EUD*P(J2)-GRA*R/RO2) 2 -DR1*P(J2)*(F(I2+N11)+F(I2+1)-F(I2+N1)-F(I2)) 3 -DZ1*P(J3)*(F(I2+N11)-F(I2+1)+F(I2+N1)-F(I2)) 4 -DR1*(F(I4+N11)+F(I4+1)-F(I4+N1)-F(I4))/P(J1) U(JJ1)=U(JJ1) 1 -DR1*(F(I2+N11)*F(I1+N11)+F(I2+1)*F(I1+1) 2 -F(I2+N1)*F(I1+N1)-F(I2)*F(I1)) 3 -DZ1*(F(I3+N11)*F(I1+N11)-F(I3+1)*F(I1+1) 4 +F(I3+N1)*F(I1+N1)-F(I3)*F(I1)) 501 CONTINUE C DO 502 I=1,NR1 R=HHX(I)+XX4 RO2=R*R+Z*Z RO2=AMAX1(RO2,AR2) J1=I J2=J1+N1 J3=J2+N1 J4=J3+N1 J5=J4+N1 J6=J5+N1 J7=J6+N1 JJ1=I X2=RO2/AR2-1.0 X2=ARU*X2*X2 X2=X2/(X2+1.0) C U(JJ1+KK5)=U(JJ1+KK5)*X2+P(J6)*(1.0-X2) U(JJ1+KK4)=U(JJ1+KK4)*X2+P(J5)*(1.0-X2) U(JJ1+KK3)=U(JJ1+KK3)*X2+P(J4)*(1.0-X2) U(JJ1+KK2)=U(JJ1+KK2)*X2+P(J3)*(1.0-X2) U(JJ1+K3)=U(JJ1+K3)*X2+P(J2)*(1.0-X2) U(JJ1)=U(JJ1)*X2+P(J1)*(1.0-X2) U(JJ1)=AMAX1(U(JJ1),RO02) 502 CONTINUE C DO 621 M=1,NB L4=L1+N2*(M-1) L3=K3*(M-1) DO 621 I=1,N1 I1=I+L4 J1=I+L3 FF(J1)=F(I1) 621 CONTINUE C C STEP OF J=J DO 90 J=2,NZ1 C C CURRENT L1=N1*(J-1) DO 38 I=1,NR1 I1=I+L1 I5=I1+M4 I6=I1+M5 J7=I+K6 P(J7)=0.5*((F(I5+N11)-F(I5+1)+F(I5+N1)-F(I5))/HZ 1 -(F(I6+N11)+F(I6+1)-F(I6+N1)-F(I6))/HR) C 2 -PP(I1) 38 CONTINUE DO 39 M=1,NB L4=L1+N2*(M-1) L2=N1*(M-1) L3=N1+K3*(M-1) DO 39 I=1,NR1 I1=I+L4 J1=I+L2 J2=I+L3 P(J1)=0.25*(F(I1)+F(I1+1)+F(I1+N1)+F(I1+N11)) U(J2)=P(J1) 39 CONTINUE C C FIRST STEP XX4=0.5*HR*FLOAT(2*NRP-NR1-1) Z=0.5*HZ*FLOAT(2*J-2) DO 50 I=1,NR1 R=HHX(I)+XX4 RO2=R*R+Z*Z RO2=AMAX1(RO2,AR2) I1=I+L1 I2=I1+N2 I3=I2+N2 I4=I3+N2 I5=I4+N2 I6=I5+N2 J1=I J2=J1+N1 J3=J2+N1 J4=J3+N1 J5=J4+N1 J6=J5+N1 J7=J6+N1 JJ1=I+N1 X2=RO2/AR2-1.0 X2=ARU*X2*X2 X2=X2/(X2+1.0) C U(JJ1+KK5)=U(JJ1+KK5) 1 +DR1*(F(I3+N11)*F(I5+N11)-F(I2+N11)*F(I6+N11) 2 +F(I3+1)*F(I5+1)-F(I2+1)*F(I6+1) 3 -F(I3+N1)*F(I5+N1)+F(I2+N1)*F(I6+N1) 4 -F(I3)*F(I5)+F(I2)*F(I6)) U(JJ1+KK4)=U(JJ1+KK4) 1 -DZ1*(F(I3+N11)*F(I5+N11)-F(I2+N11)*F(I6+N11) 2 -F(I3+1)*F(I5+1)+F(I2+1)*F(I6+1) 3 +F(I3+N1)*F(I5+N1)-F(I2+N1)*F(I6+N1) 4 -F(I3)*F(I5)+F(I2)*F(I6)) X3=DR1*(F(I2+N11)+F(I2+1)-F(I2+N1)-F(I2)) 1 +DZ1*(F(I3+N11)-F(I3+1)+F(I3+N1)-F(I3)) U(JJ1+KK3)=U(JJ1+KK3)-P(J4)*GAM*X3 1 -DR1*P(J2)*(F(I4+N11)+F(I4+1)-F(I4+N1)-F(I4)) 2 -DZ1*P(J3)*(F(I4+N11)-F(I4+1)+F(I4+N1)-F(I4)) U(JJ1+KK2)=U(JJ1+KK2)+T1*((-P(J7)*P(J5))/P(J1) 1 -EUD*P(J3)-GRA*Z/RO2) 2 -DR1*P(J2)*(F(I3+N11)+F(I3+1)-F(I3+N1)-F(I3)) 3 -DZ1*P(J3)*(F(I3+N11)-F(I3+1)+F(I3+N1)-F(I3)) 4 -DZ1*(F(I4+N11)-F(I4+1)+F(I4+N1)-F(I4))/P(J1) U(JJ1+K3)=U(JJ1+K3)+T1*((P(J7)*P(J6))/P(J1) 1 -EUD*P(J2)-GRA*R/RO2) 2 -DR1*P(J2)*(F(I2+N11)+F(I2+1)-F(I2+N1)-F(I2)) 3 -DZ1*P(J3)*(F(I2+N11)-F(I2+1)+F(I2+N1)-F(I2)) 4 -DR1*(F(I4+N11)+F(I4+1)-F(I4+N1)-F(I4))/P(J1) U(JJ1)=U(JJ1) 1 -DR1*(F(I2+N11)*F(I1+N11)+F(I2+1)*F(I1+1) 2 -F(I2+N1)*F(I1+N1)-F(I2)*F(I1)) 3 -DZ1*(F(I3+N11)*F(I1+N11)-F(I3+1)*F(I1+1) 4 +F(I3+N1)*F(I1+N1)-F(I3)*F(I1)) 50 CONTINUE C DO 51 I=1,NR1 R=HHX(I)+XX4 RO2=R*R+Z*Z RO2=AMAX1(RO2,AR2) J1=I J2=J1+N1 J3=J2+N1 J4=J3+N1 J5=J4+N1 J6=J5+N1 J7=J6+N1 JJ1=I+N1 X2=RO2/AR2-1.0 X2=ARU*X2*X2 X2=X2/(X2+1.0) C U(JJ1+KK5)=U(JJ1+KK5)*X2+P(J6)*(1.0-X2) U(JJ1+KK4)=U(JJ1+KK4)*X2+P(J5)*(1.0-X2) U(JJ1+KK3)=U(JJ1+KK3)*X2+P(J4)*(1.0-X2) U(JJ1+KK2)=U(JJ1+KK2)*X2+P(J3)*(1.0-X2) U(JJ1+K3)=U(JJ1+K3)*X2+P(J2)*(1.0-X2) U(JJ1)=U(JJ1)*X2+P(J1)*(1.0-X2) U(JJ1)=AMAX1(U(JJ1),RO02) 51 CONTINUE C C C SECOND STEP C DO 62 M=1,NB L4=L1+N2*(M-1) L3=N1+K3*(M-1) DO 62 I=1,N1 I1=I+L4 J1=I+L3 FF(J1+N1)=F(I1+N1) FF(J1)=F(I1) 62 CONTINUE C C CURRENT L1=N1*(J-1) DO 68 I=2,NR1 I1=I+L1 I5=I+N1+KK4 I6=I5+K3 J7=I+K6 P(J7)=0.5*((U(I5)-U(I5-N1)+U(I5-1)-U(I5-N11))/HZ 1 -(U(I6)+U(I6-N1)-U(I6-1)-U(I6-N11))/HR) C 2 -0.25*(PP(I1)+PP(I1-1)+PP(I1-N1)+PP(I1-N11)) 68 CONTINUE DO 69 M=1,NB L2=N1+K3*(M-1) L3=N1*(M-1) DO 69 I=2,NR1 I1=I+L2 J1=I+L3 P(J1)=0.25*(U(I1)+U(I1-1)+U(I1-N1)+U(I1-N11)) 69 CONTINUE C C SECOND STEP XX4=0.5*HR*FLOAT(2*NRP-NR1-2) Z=0.5*HZ*FLOAT(2*J-3) DO 80 I=2,NR1 R=HHX(I)+XX4 RO2=R*R+Z*Z RO2=AMAX1(RO2,AR2) I1=I+N1 I2=I1+K3 I3=I2+K3 I4=I3+K3 I5=I4+K3 I6=I5+K3 J1=I J2=J1+N1 J3=J2+N1 J4=J3+N1 J5=J4+N1 J6=J5+N1 J7=J6+N1 JJ1=I+L1 X1=FF(I4)/(FF(I1)*PO0) X2=RO2/AR2-1.0 X2=ARU*X2*X2 X2=X2/(X2+1.0) EAT=EAT0*X2/SQRT(X1*X1*X1) C F(JJ1+M5)=F(JJ1+M5) 1 +DR2*(U(I3)*U(I5)-U(I2)*U(I6) 2 +U(I3-N1)*U(I5-N1)-U(I2-N1)*U(I6-N1) 3 -U(I3-1)*U(I5-1)+U(I2-1)*U(I6-1) 4 -U(I3-N11)*U(I5-N11)+U(I2-N11)*U(I6-N11)) 5 +EAT*(DR4*(FF(I6+1)-2.0*FF(I6)+FF(I6-1)) 6 +DZ4*(FF(I6+N1)-2.0*FF(I6)+FF(I6-N1))) F(JJ1+M4)=F(JJ1+M4) 1 -DZ2*(U(I3)*U(I5)-U(I2)*U(I6) 2 -U(I3-N1)*U(I5-N1)+U(I2-N1)*U(I6-N1) 3 +U(I3-1)*U(I5-1)-U(I2-1)*U(I6-1) 4 -U(I3-N11)*U(I5-N11)+U(I2-N11)*U(I6-N11)) 5 +EAT*(DR4*(FF(I5+1)-2.0*FF(I5)+FF(I5-1)) 6 +DZ4*(FF(I5+N1)-2.0*FF(I5)+FF(I5-N1))) X3=DR2*(U(I2)+U(I2-N1)-U(I2-1)-U(I2-N11)) 2 +DZ2*(U(I3)-U(I3-N1)+U(I3-1)-U(I3-N11)) F(JJ1+M3)=F(JJ1+M3)-P(J4)*GAM*X3 1 -DR2*P(J2)*(U(I4)+U(I4-N1)-U(I4-1)-U(I4-N11)) 2 -DZ2*P(J3)*(U(I4)-U(I4-N1)+U(I4-1)-U(I4-N11)) 3 +PMU*(DR4*(FF(I4+1)-2.0*FF(I4)+FF(I4-1)) 4 +DZ4*(FF(I4+N1)-2.0*FF(I4)+FF(I4-N1))) F(JJ1+M2)=F(JJ1+M2)+T*((-P(J7)*P(J5))/P(J1) 1 -EUD*P(J3)-GRA*Z/RO2) 2 -DR2*P(J2)*(U(I3)+U(I3-N1)-U(I3-1)-U(I3-N11)) 3 -DZ2*P(J3)*(U(I3)-U(I3-N1)+U(I3-1)-U(I3-N11)) 4 -DZ2*(U(I4)-U(I4-N1)+U(I4-1)-U(I4-N11))/P(J1) 5 +RMU*(DR4*(FF(I3+1)-2.0*FF(I3)+FF(I3-1)) 6 +DZ4*(FF(I3+N1)-2.0*FF(I3)+FF(I3-N1)))/FF(I1) F(JJ1+N2)=F(JJ1+N2)+T*((P(J7)*P(J6))/P(J1) 1 -EUD*P(J2)-GRA*R/RO2) 2 -DR2*P(J2)*(U(I2)+U(I2-N1)-U(I2-1)-U(I2-N11)) 3 -DZ2*P(J3)*(U(I2)-U(I2-N1)+U(I2-1)-U(I2-N11)) 4 -DR2*(U(I4)+U(I4-N1)-U(I4-1)-U(I4-N11))/P(J1) 5 +RMU*(DR4*(FF(I2+1)-2.0*FF(I2)+FF(I2-1)) 6 +DZ4*(FF(I2+N1)-2.0*FF(I2)+FF(I2-N1)))/FF(I1) F(JJ1)=F(JJ1) 1 -DR2*(U(I2)*U(I1)+U(I2-N1)*U(I1-N1) 2 -U(I2-1)*U(I1-1)-U(I2-N11)*U(I1-N11)) 3 -DZ2*(U(I3)*U(I1)-U(I3-N1)*U(I1-N1) 4 +U(I3-1)*U(I1-1)-U(I3-N11)*U(I1-N11)) 5 +DFU*(DR4*(FF(I1+1)-2.0*FF(I1)+FF(I1-1)) 6 +DZ4*(FF(I1+N1)-2.0*FF(I1)+FF(I1-N1))) 80 CONTINUE C DO 81 I=2,NR1 R=HHX(I)+XX4 RO2=R*R+Z*Z RO2=AMAX1(RO2,AR2) I1=I+N1 I2=I1+K3 I3=I2+K3 I4=I3+K3 I5=I4+K3 I6=I5+K3 JJ1=I+L1 X2=RO2/AR2-1.0 X2=ARU*X2*X2 X2=X2/(X2+1.0) C F(JJ1+M5)=F(JJ1+M5)*X2+FF(I6)*(1.0-X2) F(JJ1+M4)=F(JJ1+M4)*X2+FF(I5)*(1.0-X2) F(JJ1+M3)=F(JJ1+M3)*X2+FF(I4)*(1.0-X2) F(JJ1+M2)=F(JJ1+M2)*X2+FF(I3)*(1.0-X2) F(JJ1+N2)=F(JJ1+N2)*X2+FF(I2)*(1.0-X2) F(JJ1)=F(JJ1)*X2+FF(I1)*(1.0-X2) 81 CONTINUE C DO 82 I=2,NR1 JJ1=I+L1 JJ2=JJ1+N2 JJ3=JJ2+N2 JJ4=JJ3+N2 P(I+K3)=F(JJ4) P(I+K2)=F(JJ3) P(I+N1)=F(JJ2) P(I)=F(JJ1) 82 CONTINUE C DO 83 I=2,NR1 J1=I J2=J1+N1 J3=J2+N1 J4=J3+N1 JJ1=I+L1 C X1=P(J1) X2=P(J4) C Y2=AMAX1(X2,PR02) C IF(X1.GE.RO02) GO TO 83 X=AMAX1(X1/RO02,0.0) X=AMIN1(X,1.0) X3=P(J2)*P(J2)+P(J3)*P(J3) C IF(X3.LT.1.0) GO TO 84 C WRITE(12) F C WRITE(6,663) II,LAST,NR,NZ,X1,X2,X3 C 663 FORMAT(1H ,5X,4I10,1P3E15.5) C GO TO 300 C 84 CONTINUE X4=X2+GM2*X1*X3*(1.0-X) F(JJ1+M3)=AMAX1(X4,PR02) F(JJ1+M2)=F(JJ1+M2)*X F(JJ1+N2)=F(JJ1+N2)*X F(JJ1)=AMAX1(X1,RO02) 83 CONTINUE C DO 90 M=1,NB L2=K3*(M-1) DO 90 I=1,N1 I1=I+L2 U(I1)=U(I1+N1) FF(I1)=FF(I1+N1) 90 CONTINUE C C TIME=TIME+T IF(IIQ.NE.IIQ0) GO TO 401 IIQ=0 C X1=2.0 NR11=NR1-2*NRP DO 102 J=1,NZ2 L1=N1*(J-1) DO 102 I=1,NR2 X=(2.0*HHX(I)/HR-FLOAT(NR11+2))/FLOAT(NR11) X=(X1+1.0)*X+X1 X=AMAX1(-X,0.0) X=X*X I1=I+L1 F(I1+M5)=(1.0-X)*F(I1+M5)+BIS*X C IF(II.LE.256) F(I6)=(1.0-X)*F(I6)+0.0*X F(I1+M4)=(1.0-X)*F(I1+M4) F(I1+M3)=(1.0-X)*F(I1+M3)+PR01*X F(I1+M2)=(1.0-X)*F(I1+M2) F(I1+N2)=(1.0-X)*F(I1+N2)+VSW*X F(I1)=(1.0-X)*F(I1)+RO01*X 102 CONTINUE C X1=ARU1 XX4=0.5*HR*FLOAT(2*NRP-NR1-2) DO 104 J=1,NZ2 L1=N1*(J-1) Z=0.5*HZ*FLOAT(2*J-3) DO 104 I=1,NR2 R=HHX(I)+XX4 RO2=R*R+Z*Z RO2=AMAX1(RO2,AR2) X2=RO2/AR2-1.0 X2=AMAX1(X2,0.0) X2=X1*X2*X2 X2=X2/(X2+1.0) I1=I+L1 F(I1+M5)=F(I1+M5)*X2 F(I1+M4)=F(I1+M4)*X2 F(I1+M3)=F(I1+M3)*X2+PP(I1+M3)*(1.0-X2) F(I1+M2)=F(I1+M2)*X2+PP(I1+M2)*(1.0-X2) F(I1+N2)=F(I1+N2)*X2+PP(I1+N2)*(1.0-X2) F(I1)=F(I1)*X2+PP(I1)*(1.0-X2) 104 CONTINUE C 401 IF(IIP.NE.IIP0) GO TO 100 IIP=0 C JXG=N2 JYG=NZ2 JXN=NR2 LANK=20 IP=0 LB=1 LL=6 I1=II C CALL SCALM(NR,NZ,NRP,LB,LL,LLL,HR,HZ,AR2,ZMIN,ZMAX,F) DO 120 J=1,2 LB=1+3*(J-1) LL=3 C CALL GRAPE(JXG,JYG,LB,LL,ZMIN,ZMAX,I1,LLL,JXN,LANK,IP,F) 120 CONTINUE C C 501 IF(IIS.NE.IIS0) GO TO 100 IF(IIS.NE.IIS0) GO TO 100 IIS=0 WRITE(12) F 100 CONTINUE 300 CONTINUE 9 CONTINUE REWIND 12 STOP END SUBROUTINE EQUIP4(NR,NZ,NRP,RO01,PR01,CJ,CP,F,P) DIMENSION P(1),F(1),CJ(1),CP(1) C CJ 1-R1 2-DR 3-DZ 4-B0 5-RM 6-A2 7-AID 8-AR1 9-GAM 10-GRA C CP 1-R0 2-Z0 3-RA 4-VO0 5-P0 6-GRA 7-07 8-VSW 9-09 10-10 C VR0=CP(4) R11=CJ(5) R22=CJ(6) AID=CJ(7) AR1=CJ(8) GAM=CJ(9) GRA=CJ(10) RATT=CJ(2) P0=(GAM-1.0)*GRA/GAM BIS=CP(10)*1.0E-4 X1=1.0/(GAM-1.0) X2=GAM*X1 B0=CJ(4) RO02=RATT*RO01 PR02=PR01 C IC0=4 ICD=6 PI=3.1415926 NR1=NR+1 NR2=NR+2 N1=NR2 N11=N1+1 NB=8 NZ1=NZ+1 NZ2=NZ+2 N2=N1*NZ2 HR=CP(1)/FLOAT(NR1) HZ=CP(2)/FLOAT(NZ1) VSW=CP(8) C DO 10 J=1,NZ2 DO 10 I=1,NR2 R=0.5*HR*FLOAT(2*I-NR2-1+2*NRP) Z=0.5*HZ*FLOAT(2*J-3) RO2=R*R+Z*Z C RO2=AMAX1(RO2,AR2) RO1=SQRT(RO2) RO3=RO1*RO2 COST=R/RO1 SINT=Z/RO1 I1=I+N1*(J-1) I2=I1+N2 I3=I2+N2 I4=I3+N2 I5=I4+N2 I6=I5+N2 X=1.0 IF(R.GE.0.0) GO TO 12 IF(RO1.LT.R22) X=(RO1-R11)/(R22-R11) IF(RO1.LT.R11) X=0.0 GO TO 14 12 CONTINUE IF(ABS(Z).LT.R22) X=(ABS(Z)-R11)/(R22-R11) IF(ABS(Z).LT.R11) X=0.0 14 CONTINUE X1=RO01*(R11/RO1)*(VSW/VR0) X2=PR01*(R11*R11/RO2) X3=VR0*COST X4=VR0*SINT F(I1)=(1.0-X)*X1+X*RO01 F(I2)=((1.0-X)*X1*X3+X*VSW*RO01)/F(I1) F(I3)=(1.0-X)*X1*X4/F(I1) F(I4)=(1.0-X)*X2+X*PR01 F(I5)=-2.0*B0*R*Z/(RO2*RO2) F(I6)=B0*(R*R-Z*Z)/(RO2*RO2) F(I4)=AMAX1(F(I4),PR02) X=F(I1) F(I1)=AMAX1(F(I1),RO02) F(I2)=F(I2)*X/F(I1) F(I3)=F(I3)*X/F(I1) F(I6)=F(I6)+BIS P(I1)=F(I1) P(I2)=F(I2) P(I3)=F(I3) P(I4)=F(I4) 10 CONTINUE RETURN END C /SUBA/SUB1AA SUBROUTINE SCALM(NX,NZ,NXP,LB,LL,LLL,HX,HZ,AR2,ZMIN,ZMAX,P) DIMENSION ZMIN(LLL),ZMAX(LLL),P(1) NX2=NX+2 NZ2=NZ+2 N2=NX2*NZ2 L1=LB+LL-1 DO 10 J=LB,L1 I1=1+N2*(J-1) ZMIN(J)=P(I1)*(1.0-1.0E-5)-1.0E-10 ZMAX(J)=P(I1)*(1.0+1.0E-5)+1.0E-10 DO 10 K=1,NZ2 DO 10 I=1,NX2 Z=0.5*HZ*FLOAT(2*K-NZ2-1) X=0.5*HX*FLOAT(2*I-NX-3+2*NXP) RO=X*X+Z*Z IF(RO.LE.AR2) GO TO 10 I1=I+NX2*(K-1)+N2*(J-1) Y=P(I1) ZMIN(J)=AMIN1(ZMIN(J),Y) ZMAX(J)=AMAX1(ZMAX(J),Y) 10 CONTINUE RETURN END SUBROUTINE GRAPD(JXG,JYG,LB,LL,ZMIN,ZMAX,II,LLL,JXN,LANK,IP, 1 P) DIMENSION IA(132),JJG(61),ZMIN(LLL),ZMAX(LLL),P(1) DATA JJG/'.','1','2','3','4','5','6','7','8','9', 1 'A','B','C','D','E','F','G','H','I','J', 2 'K','L','M','N','O','P','Q','R','S','T', 3 'U','V','W','X','Y','Z','0','1','2','3', 4 'A','B','C','D','E','F','G','H','I','J', 5 'K','L','M','N','O','P','Q','R','S','T',' '/ L1=LB+LL-1 WRITE(6,522) II,LB,LL,JXN,JYG,LANK,IP,L1,LLL,JXG, 1 (ZMIN(I),ZMAX(I),I=LB,L1) 522 FORMAT(1H0,3X,3HII=,10I5/(1H ,3X,8E15.5)) C WRITE(6,554) DO 500 JJ12=1,JYG JJ=JYG+1-JJ12 DO 524 I=1,132 524 IA(I)=JJG(61) DO 506 K=LB,L1 M=1+K-LB YMIN=ZMIN(K) YMAX=ZMAX(K) X1=FLOAT(LANK)/(YMAX-YMIN) DO 506 I=1,JXN I1=I+JXN*(JJ-1)+JXN*JYG*(K-1) Y=P(I1) IF(Y.LT.YMIN) Y=YMIN IF(Y.GT.YMAX) Y=YMAX J1=IFIX(X1*(Y-YMIN))+1 IF(IP.NE.0) J1=IFIX(X1*(YMAX-Y))+1 IF(J1.GT.LANK) J1=LANK JX=I+JXN*(M-1)+M-1 IF(JX.LE.130) IA(JX)=JJG(J1) 506 CONTINUE WRITE(6,550) (IA(I),I=1,130) 500 CONTINUE 550 FORMAT(1H ,1X,130A1) 554 FORMAT(1H ) RETURN END FUNCTION CVMGP(X,Y,Z) CVMGP=Y IF(Z.GE.0.0) CVMGP=X RETURN END FUNCTION CVMGM(X,Y,Z) CVMGM=Y IF(Z.LT.0.0) CVMGM=X RETURN END SUBROUTINE GRAPE(JXG,JYG,LB,LL,ZMIN,ZMAX,II,LLL,JXN,LANK,IP, 1 P) DIMENSION IA(132),JJG(61),ZMIN(LLL),ZMAX(LLL),P(1) DATA JJG/'J','I','H','G','F','E','D','C','B','A', 1 '.','1','2','3','4','5','6','7','8','9', 2 '0','L','M','N','O','P','Q','R','S','T', 3 'U','V','W','X','Y','Z','7','8','9','0', 4 'A','B','C','D','E','F','G','H','I','J', 5 'K','L','M','N','O','P','Q','R','S','T',' '/ L1=LB+LL-1 LAN1=21 WRITE(6,522) II,LB,LL,JXN,JYG,LAN1,IP, 1 (ZMIN(I),ZMAX(I),I=LB,L1) 522 FORMAT(1H0,3X,3HII=,7I5/(1H ,3X,8E15.5)) C WRITE(6,554) DO 500 JJ12=1,JYG JJ=JYG+1-JJ12 DO 524 I=1,132 524 IA(I)=JJG(61) DO 506 K=LB,L1 M=1+K-LB YMIN=ZMIN(K) YMAX=ZMAX(K) YMAX=AMAX1(YMAX,-YMIN) YMIN=-YMAX X1=FLOAT(LAN1)/(YMAX-YMIN) DO 506 I=1,JXN I1=I+JXN*(JJ-1)+JXN*JYG*(K-1) Y=P(I1) IF(Y.LT.YMIN) Y=YMIN IF(Y.GT.YMAX) Y=YMAX J1=IFIX(X1*(Y-YMIN))+1 IF(IP.NE.0) J1=IFIX(X1*(YMAX-Y))+1 IF(J1.GT.LAN1) J1=LAN1 JX=I+JXN*(M-1)+M-1 IF(JX.LE.131) IA(JX)=JJG(J1) 506 CONTINUE WRITE(6,550) (IA(I),I=1,131) 500 CONTINUE 550 FORMAT(1H ,1X,131A1) 554 FORMAT(1H ) RETURN END SUBROUTINE CLOCK(IT) DOUBLE PRECISION IT X=0.0 Y=SECNDS(X) IT=1.00000*Y RETURN END