C PROGRAM OGNCB7 C C OGNCB72 M3RTM83 FROM /MHD3/RECTM80 C COMPUTER SIMULATION OF MAGNETOSPHERIC EQUILIBRIUM C THREE DIMENSIONAL MHD MODEL C CARTESIAN COORDINATE FINITE RESISTIVITY 45 DEGREE BOUDARY C 1983.08.18 BY TATSUKI OGINO C PARAMETER (NX=60,NY=30,NZ=30,N1=NX+2,N2=N1*(NY+2),N3=N2*(NZ+2)) PARAMETER (NB=8,NBB=11,N4=N3*NB,N5=N3*NBB,N6=N3*18) PARAMETER (LAST=4,NIL=1,LLL=18,IBIG=0) PARAMETER (IIQ0=4,IIP0=4,IIS0=4,THX=8.00) PARAMETER (RMUU=0.005,EUD0=0.00,RO01=5.0E-4,PR01=1.78E-8) PARAMETER (NXP=0,ARUU=100.0,PMU0=1.0,DFU0=1.0) PARAMETER (EATT=0.010,RRAT=0.2,RRA1=0.2) PARAMETER (N11=N1+1,N12=N2+1,N21=N1+N2,N22=N21+1) PARAMETER (NX1=NX+1,NX2=NX+2,NY1=NY+1,NY2=NY+2) PARAMETER (NZ1=NZ+1,NZ2=NZ+2,NXZ=NX*NZ) PARAMETER (NIN=0,ITAP=1,MN1=N1*NY1,MN2=N2*NZ1) PARAMETER (M2=N3*2,M3=N3*3,M4=N3*4,M5=N3*5,M6=N3*6,M7=N3*7) PARAMETER (M8=N3*8,M9=N3*9,M10=N3*10,M11=N3*11,M12=N3*12) PARAMETER (K2=N2*2,K3=N2*3,K4=N2*4,K5=N2*5,K6=N2*6) PARAMETER (K7=N2*7,K8=N2*8,K9=N2*9,K10=N2*10,K11=N2*11) PARAMETER (KK2=K3*2,KK3=K3*3,KK4=K3*4,KK5=K3*5) PARAMETER (KK6=K3*6,KK7=K3*7,KK8=K3*8) PARAMETER (NP1=NX2/2,NP2=NP1*NY1,NP3=NP2*NZ1,NP4=NP3*3) PARAMETER (NP11=NP1+1,NP12=NP2+1,NP21=NP1+NP2,NP22=NP21+1) C DOUBLE PRECISION IT1,IT2,ITD,IT0 DIMENSION F(N4),FF(KK8),U(KK8),P(K11),HHX(N1),FBB(16), 1 PP(NP4),ZMAX(18),ZMIN(18),CJ(10),CP(11) C M 1RO 2VR 3VO 4VZ 5PR 6BR 7BO 8BZ 9JR 10JO 11JZ C C CJ R1 DR DZ B0 RM A2 AID AR1 GAM GRA DATA CJ/10.0,0.3,0.5,1.0,1.9,-0.01,0.7,0.8,0.9,1.0/ C C CP XL YL ZL RA VO0 P0 GRA 07 VSW IBR IBZ C DATA CP/61.0,31.0,31.0,4.0,6.81,2.68,1.35,7.0,.044,9.0,-0.0/ DATA CP/61.0,31.0,31.0,3.5,0.00,2.68,1.35,7.0,.044,0.0,1.50/ C REWIND 12 C DO 300 III=1,1 EAT0=EATT RMU0=RMUU RMU0=RMUU ARU=ARUU PI=3.1415926 GAM=5.0/3.0 GM1=2.0/(GAM-1) CP(6)=10.0*(GAM-1.0)*CP(7)/GAM GRA=CP(7)*1.0E-6 PO0=CP(6)*1.0E-7 BIS=CP(11)*1.0E-4 VSW=CP(9) AR1=CP(4) AR2=AR1*AR1 HX=CP(1)/FLOAT(NX1) HY=CP(2)/FLOAT(NY1) HZ=CP(3)/FLOAT(NZ1) T=0.5*HX*THX T1=0.5*T CJ(8)=AR1 CJ(9)=GAM CJ(10)=GRA C DX1=0.25*T1/HX DY1=0.25*T1/HY DZ1=0.25*T1/HZ DX2=0.25*T/HX DY2=0.25*T/HY DZ2=0.25*T/HZ DX3=T1/(HX*HX) DY3=T1/(HY*HY) DZ3=T1/(HZ*HZ) DX4=T/(HX*HX) DY4=T/(HY*HY) DZ4=T/(HZ*HZ) C RMU=RMU0*RO01 PMU=RMU0*PMU0 DFU=RMU0*DFU0 EUD=EUD0 RO02=RRAT*RO01 PR02=RRAT*PR01 C WRITE (6,12) III,LAST,NX,NY,NZ,N1,N2,N3,N4,N5,N6,EAT0,RMU0,ARU, 1 EUD,RRAT,HX,HY,HZ,T,T1,RO01,PR01,GRA,DX2,DY2,DZ2,DX4,DY4,DZ4, 2 (CP(I),I=1,11),(CJ(J),J=1,10) 12 FORMAT(1H ,5X,11I10/(1H ,5X,1P8E15.5)) C DO 20 I=1,N4 20 F(I)=0.0 DO 221 I=1,N1 221 HHX(I)=HX*FLOAT(I) DO 222 I=1,NB FBB(I+NB)=1.0 222 FBB(I)=1.0 FBB(4)=-1.0 FBB(6)=-1.0 FBB(7)=-1.0 FBB(3+NB)=-1.0 FBB(7+NB)=-1.0 C C INITIAL CONDITION CALL EQUIB7(NX,NY,NZ,NXP,RO01,PR01,RRA1,CJ,CP,F,PP) C GO TO 882 C WRITE(6,212) (F(I),I=1,N4) 212 FORMAT(1H ,2X,1P10E12.4) C 882 CONTINUE C C C INITIAL CONDITION C IF(NIN.EQ.0) GO TO 400 DO 410 I=1,ITAP C READ(12) F 410 CONTINUE 400 CONTINUE C C THREE DIMENSIONAL CARTESIAN MODEL TIME=0.0 IIQ=0 IIP=0 IIS=0 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 C BOUNDARY CONDITION AT NY=1 AND NY=NY2 DO 30 M=1,NB I3=N2*2 DO 31 J=2,NY1 L1=N1*(J-1)+N3*(M-1) I1=2+L1 I2=I1+MN2 P(2)=F(I2-N2-1) U(2)=F(I1+N2)*FBB(M) DO 311 I=3,NX1 I1=I+L1 I2=I1+MN2 P(I)=(4.0*F(I2-N2-1)-F(I2-I3-2))/3.0 311 U(I)=F(I1+N2)*FBB(M) DO 31 I=2,NX1 I2=I+L1+MN2 F(I2)=P(I) 31 F(I2-MN2)=U(I) I3=N1*2 DO 32 K=1,NZ2 L1=N2*(K-1)+N3*(M-1) I1=2+L1 I2=I1+MN1 P(2)=F(I2-N1-1) U(2)=F(I1+N1)*FBB(M+NB) DO 321 I=3,NX1 I1=I+L1 I2=I1+MN1 P(I)=(4.0*F(I2-N1-1)-F(I2-I3-2))/3.0 321 U(I)=F(I1+N1)*FBB(M+NB) DO 32 I=2,NX1 I2=I+L1+MN1 F(I2)=P(I) 32 F(I2-MN1)=U(I) C BOUNDARY CONDITION AT NX=NX2 DO 33 K=1,NZ2 I1=N1+N2*(K-1)+N3*(M-1) DO 331 J=1,NY2 P(J)=(4.0*F(I1-1)-F(I1-2))/3.0 331 I1=I1+N1 I1=N1+N2*(K-1)+N3*(M-1) DO 33 J=1,NY2 F(I1)=P(J) 33 I1=I1+N1 30 CONTINUE C IF(NIN.EQ.0.AND.II.EQ.1) WRITE(12) F C C STEP OF K=1 K=1 C C CURRENT DO 3810 J=1,NY1 L1=N1*(J-1) L2=N1*(J-1) DO 3810 I=1,NX1 I6=I+L1+M5 I7=I6+N3 I8=I7+N3 J9=I+L2+K8 P(J9+K2)=0.25*((F(I7+N22)+F(I7+N12)+F(I7+N11)+F(I7+1) 1 -F(I7+N21)-F(I7+N2)-F(I7+N1)-F(I7))/HX 2 -(F(I6+N22)-F(I6+N12)+F(I6+N11)-F(I6+1) 3 +F(I6+N21)-F(I6+N2)+F(I6+N1)-F(I6))/HY) P(J9+N2)=0.25*((F(I6+N22)+F(I6+N12)-F(I6+N11)-F(I6+1) 1 +F(I6+N21)+F(I6+N2)-F(I6+N1)-F(I6))/HZ 2 -(F(I8+N22)+F(I8+N12)+F(I8+N11)+F(I8+1) 3 -F(I8+N21)-F(I8+N2)-F(I8+N1)-F(I8))/HX) P(J9)=0.25*((F(I8+N22)-F(I8+N12)+F(I8+N11)-F(I8+1) 1 +F(I8+N21)-F(I8+N2)+F(I8+N1)-F(I8))/HY 2 -(F(I7+N22)+F(I7+N12)-F(I7+N11)-F(I7+1) 3 +F(I7+N21)+F(I7+N2)-F(I7+N1)-F(I7))/HZ) 3810 CONTINUE DO 3812 J=1,NY1 L1=NP1*(J-1)+NP2*(K-1) L2=N1*(J-1)+K8 DO 3814 I=1,NP1 J9=I+L2 I1=I+L1 P(J9+K2)=P(J9+K2)-PP(I1+2*NP3) P(J9+N2)=P(J9+N2)-PP(I1+NP3) P(J9)=P(J9)-PP(I1) 3814 CONTINUE DO 3812 I=2,NP1 J9=I+NP1-1+L2 I1=NP1-I+1+L1 P(J9+K2)=P(J9+K2)+PP(I1+2*NP3) P(J9+N2)=P(J9+N2)+PP(I1+NP3) P(J9)=P(J9)-PP(I1) 3812 CONTINUE C DO 3811 M=1,NB DO 3811 J=1,NY1 L1=N1*(J-1)+N3*(M-1) L2=N1*(J-1)+N2*(M-1) L3=N1*(J-1)+K3*(M-1) DO 3811 I=1,NX1 I1=I+L1 J1=I+L2 J2=I+L3 P(J1)=0.125*(F(I1)+F(I1+1)+F(I1+N1)+F(I1+N11) 1 +F(I1+N2)+F(I1+N12)+F(I1+N21)+F(I1+N22)) U(J2)=P(J1) 3811 CONTINUE C C ZEROTH STEP XX4=0.5*HX*FLOAT(2*NXP-NX1-1) Z=0.5*HZ*FLOAT(2*K-2) DO 5010 J=1,NY1 L1=N1*(J-1) L2=N1*(J-1) L3=L2 Y=0.5*HY*FLOAT(2*J-2) C DO 5011 I=1,NX1 X=HHX(I)+XX4 RO2=X*X+Y*Y+Z*Z RO2=MAX(RO2,AR2) RO3=RO2*SQRT(RO2) I1=I+L1 I2=I1+N3 I3=I2+N3 I4=I3+N3 I5=I4+N3 I6=I5+N3 I7=I6+N3 I8=I7+N3 J1=I+L2 J2=J1+N2 J3=J2+N2 J4=J3+N2 J5=J4+N2 J6=J5+N2 J7=J6+N2 J8=J7+N2 J9=J8+N2 J10=J9+N2 J11=J10+N2 JJ1=I+L3 C U(JJ1+KK7)=P(J8) 1 +DX1*(F(I4+N22)*F(I6+N22)-F(I2+N22)*F(I8+N22) 2 +F(I4+N12)*F(I6+N12)-F(I2+N12)*F(I8+N12) 3 +F(I4+N11)*F(I6+N11)-F(I2+N11)*F(I8+N11) 4 +F(I4+1)*F(I6+1)-F(I2+1)*F(I8+1) 5 -F(I4+N21)*F(I6+N21)+F(I2+N21)*F(I8+N21) 6 -F(I4+N2)*F(I6+N2)+F(I2+N2)*F(I8+N2) 7 -F(I4+N1)*F(I6+N1)+F(I2+N1)*F(I8+N1) 8 -F(I4)*F(I6)+F(I2)*F(I8)) 1 -DY1*(F(I3+N22)*F(I8+N22)-F(I4+N22)*F(I7+N22) 2 -F(I3+N12)*F(I8+N12)+F(I4+N12)*F(I7+N12) 3 +F(I3+N11)*F(I8+N11)-F(I4+N11)*F(I7+N11) 4 -F(I3+1)*F(I8+1)+F(I4+1)*F(I7+1) 5 +F(I3+N21)*F(I8+N21)-F(I4+N21)*F(I7+N21) 6 -F(I3+N2)*F(I8+N2)+F(I4+N2)*F(I7+N2) 7 +F(I3+N1)*F(I8+N1)-F(I4+N1)*F(I7+N1) 8 -F(I3)*F(I8)+F(I4)*F(I7)) U(JJ1+KK6)=P(J7) 1 +DZ1*(F(I3+N22)*F(I8+N22)-F(I4+N22)*F(I7+N22) 2 +F(I3+N12)*F(I8+N12)-F(I4+N12)*F(I7+N12) 3 -F(I3+N11)*F(I8+N11)+F(I4+N11)*F(I7+N11) 4 -F(I3+1)*F(I8+1)+F(I4+1)*F(I7+1) 5 +F(I3+N21)*F(I8+N21)-F(I4+N21)*F(I7+N21) 6 +F(I3+N2)*F(I8+N2)-F(I4+N2)*F(I7+N2) 7 -F(I3+N1)*F(I8+N1)+F(I4+N1)*F(I7+N1) 8 -F(I3)*F(I8)+F(I4)*F(I7)) 1 -DX1*(F(I2+N22)*F(I7+N22)-F(I3+N22)*F(I6+N22) 2 +F(I2+N12)*F(I7+N12)-F(I3+N12)*F(I6+N12) 3 +F(I2+N11)*F(I7+N11)-F(I3+N11)*F(I6+N11) 4 +F(I2+1)*F(I7+1)-F(I3+1)*F(I6+1) 5 -F(I2+N21)*F(I7+N21)+F(I3+N21)*F(I6+N21) 6 -F(I2+N2)*F(I7+N2)+F(I3+N2)*F(I6+N2) 7 -F(I2+N1)*F(I7+N1)+F(I3+N1)*F(I6+N1) 8 -F(I2)*F(I7)+F(I3)*F(I6)) U(JJ1+KK5)=P(J6) 1 +DY1*(F(I2+N22)*F(I7+N22)-F(I3+N22)*F(I6+N22) 2 -F(I2+N12)*F(I7+N12)+F(I3+N12)*F(I6+N12) 3 +F(I2+N11)*F(I7+N11)-F(I3+N11)*F(I6+N11) 4 -F(I2+1)*F(I7+1)+F(I3+1)*F(I6+1) 5 +F(I2+N21)*F(I7+N21)-F(I3+N21)*F(I6+N21) 6 -F(I2+N2)*F(I7+N2)+F(I3+N2)*F(I6+N2) 7 +F(I2+N1)*F(I7+N1)-F(I3+N1)*F(I6+N1) 8 -F(I2)*F(I7)+F(I3)*F(I6)) 1 -DZ1*(F(I4+N22)*F(I6+N22)-F(I2+N22)*F(I8+N22) 2 +F(I4+N12)*F(I6+N12)-F(I2+N12)*F(I8+N12) 3 -F(I4+N11)*F(I6+N11)+F(I2+N11)*F(I8+N11) 4 -F(I4+1)*F(I6+1)+F(I2+1)*F(I8+1) 5 +F(I4+N21)*F(I6+N21)-F(I2+N21)*F(I8+N21) 6 +F(I4+N2)*F(I6+N2)-F(I2+N2)*F(I8+N2) 7 -F(I4+N1)*F(I6+N1)+F(I2+N1)*F(I8+N1) 8 -F(I4)*F(I6)+F(I2)*F(I8)) C X3=DX1*(F(I2+N22)+F(I2+N12)+F(I2+N11)+F(I2+1) 1 -F(I2+N21)-F(I2+N2)-F(I2+N1)-F(I2)) 2 +DY1*(F(I3+N22)-F(I3+N12)+F(I3+N11)-F(I3+1) 3 +F(I3+N21)-F(I3+N2)+F(I3+N1)-F(I3)) 4 +DZ1*(F(I4+N22)+F(I4+N12)-F(I4+N11)-F(I4+1) 5 +F(I4+N21)+F(I4+N2)-F(I4+N1)-F(I4)) U(JJ1+KK4)=P(J5)-P(J5)*GAM*X3 1 -DX1*P(J2)*(F(I5+N22)+F(I5+N12)+F(I5+N11)+F(I5+1) 2 -F(I5+N21)-F(I5+N2)-F(I5+N1)-F(I5)) 3 -DY1*P(J3)*(F(I5+N22)-F(I5+N12)+F(I5+N11)-F(I5+1) 4 +F(I5+N21)-F(I5+N2)+F(I5+N1)-F(I5)) 5 -DZ1*P(J4)*(F(I5+N22)+F(I5+N12)-F(I5+N11)-F(I5+1) 6 +F(I5+N21)+F(I5+N2)-F(I5+N1)-F(I5)) U(JJ1+KK3)=P(J4)+T1*((P(J9)*P(J7)-P(J10)*P(J6))/P(J1) 1 -EUD*P(J4)-GRA*Z/RO3) 2 -DX1*P(J2)*(F(I4+N22)+F(I4+N12)+F(I4+N11)+F(I4+1) 3 -F(I4+N21)-F(I4+N2)-F(I4+N1)-F(I4)) 4 -DY1*P(J3)*(F(I4+N22)-F(I4+N12)+F(I4+N11)-F(I4+1) 5 +F(I4+N21)-F(I4+N2)+F(I4+N1)-F(I4)) 6 -DZ1*P(J4)*(F(I4+N22)+F(I4+N12)-F(I4+N11)-F(I4+1) 7 +F(I4+N21)+F(I4+N2)-F(I4+N1)-F(I4)) 8 -DZ1*(F(I5+N22)+F(I5+N12)-F(I5+N11)-F(I5+1) 9 +F(I5+N21)+F(I5+N2)-F(I5+N1)-F(I5))/P(J1) U(JJ1+KK2)=P(J3)+T1*((P(J11)*P(J6)-P(J9)*P(J8))/P(J1) 1 -EUD*P(J3)-GRA*Y/RO3) 2 -DX1*P(J2)*(F(I3+N22)+F(I3+N12)+F(I3+N11)+F(I3+1) 3 -F(I3+N21)-F(I3+N2)-F(I3+N1)-F(I3)) 4 -DY1*P(J3)*(F(I3+N22)-F(I3+N12)+F(I3+N11)-F(I3+1) 5 +F(I3+N21)-F(I3+N2)+F(I3+N1)-F(I3)) 6 -DZ1*P(J4)*(F(I3+N22)+F(I3+N12)-F(I3+N11)-F(I3+1) 7 +F(I3+N21)+F(I3+N2)-F(I3+N1)-F(I3)) 8 -DY1*(F(I5+N22)-F(I5+N12)+F(I5+N11)-F(I5+1) 9 +F(I5+N21)-F(I5+N2)+F(I5+N1)-F(I5))/P(J1) U(JJ1+K3)=P(J2)+T1*((P(J10)*P(J8)-P(J11)*P(J7))/P(J1) 1 -EUD*P(J2)-GRA*X/RO3) 2 -DX1*P(J2)*(F(I2+N22)+F(I2+N12)+F(I2+N11)+F(I2+1) 3 -F(I2+N21)-F(I2+N2)-F(I2+N1)-F(I2)) 4 -DY1*P(J3)*(F(I2+N22)-F(I2+N12)+F(I2+N11)-F(I2+1) 5 +F(I2+N21)-F(I2+N2)+F(I2+N1)-F(I2)) 6 -DZ1*P(J4)*(F(I2+N22)+F(I2+N12)-F(I2+N11)-F(I2+1) 7 +F(I2+N21)+F(I2+N2)-F(I2+N1)-F(I2)) 8 -DX1*(F(I5+N22)+F(I5+N12)+F(I5+N11)+F(I5+1) 9 -F(I5+N21)-F(I5+N2)-F(I5+N1)-F(I5))/P(J1) U(JJ1)=P(J1) 1 -DX1*(F(I2+N22)*F(I1+N22)+F(I2+N12)*F(I1+N12) 2 +F(I2+N11)*F(I1+N11)+F(I2+1)*F(I1+1) 3 -F(I2+N21)*F(I1+N21)-F(I2+N2)*F(I1+N2) 4 -F(I2+N1)*F(I1+N1)-F(I2)*F(I1)) 5 -DY1*(F(I3+N22)*F(I1+N22)-F(I3+N12)*F(I1+N12) 6 +F(I3+N11)*F(I1+N11)-F(I3+1)*F(I1+1) 7 +F(I3+N21)*F(I1+N21)-F(I3+N2)*F(I1+N2) 8 +F(I3+N1)*F(I1+N1)-F(I3)*F(I1)) 1 -DZ1*(F(I4+N22)*F(I1+N22)+F(I4+N12)*F(I1+N12) 2 -F(I4+N11)*F(I1+N11)-F(I4+1)*F(I1+1) 3 +F(I4+N21)*F(I1+N21)+F(I4+N2)*F(I1+N2) 4 -F(I4+N1)*F(I1+N1)-F(I4)*F(I1)) 5011 CONTINUE C DO 5201 I=1,NX1 X=HHX(I)+XX4 RO2=X*X+Y*Y+Z*Z I1=I+L2 I2=I1+N2 I3=I2+N2 I4=I3+N2 I5=I4+N2 I6=I5+N2 I7=I6+N2 I8=I7+N2 J1=I+L3 X2=MAX(RO2/AR2-1.0,0.0) X2=ARU*X2*X2 X2=X2/(X2+1.0) U(J1+KK7)=U(J1+KK7)*X2+P(I8)*(1.0-X2) U(J1+KK6)=U(J1+KK6)*X2+P(I7)*(1.0-X2) U(J1+KK5)=U(J1+KK5)*X2+P(I6)*(1.0-X2) U(J1+KK4)=U(J1+KK4)*X2+P(I5)*(1.0-X2) U(J1+KK3)=U(J1+KK3)*X2+P(I4)*(1.0-X2) U(J1+KK2)=U(J1+KK2)*X2+P(I3)*(1.0-X2) U(J1+K3)=U(J1+K3)*X2+P(I2)*(1.0-X2) U(J1)=U(J1)*X2+P(I1)*(1.0-X2) 5201 CONTINUE 5010 CONTINUE C DO 620 M=1,NB L1=N3*(M-1) L2=K3*(M-1) DO 620 I=1,N2 I1=I+L1 J1=I+L2 FF(J1)=F(I1) 620 CONTINUE C C STEP OF K=K DO 90 K=2,NZ1 C C CURRENT DO 38 J=1,NY1 L1=N1*(J-1)+N2*(K-1) L2=N1*(J-1) DO 38 I=1,NX1 I6=I+L1+M5 I7=I6+N3 I8=I7+N3 J9=I+L2+K8 P(J9+K2)=0.25*((F(I7+N22)+F(I7+N12)+F(I7+N11)+F(I7+1) 1 -F(I7+N21)-F(I7+N2)-F(I7+N1)-F(I7))/HX 2 -(F(I6+N22)-F(I6+N12)+F(I6+N11)-F(I6+1) 3 +F(I6+N21)-F(I6+N2)+F(I6+N1)-F(I6))/HY) P(J9+N2)=0.25*((F(I6+N22)+F(I6+N12)-F(I6+N11)-F(I6+1) 1 +F(I6+N21)+F(I6+N2)-F(I6+N1)-F(I6))/HZ 2 -(F(I8+N22)+F(I8+N12)+F(I8+N11)+F(I8+1) 3 -F(I8+N21)-F(I8+N2)-F(I8+N1)-F(I8))/HX) P(J9)=0.25*((F(I8+N22)-F(I8+N12)+F(I8+N11)-F(I8+1) 1 +F(I8+N21)-F(I8+N2)+F(I8+N1)-F(I8))/HY 2 -(F(I7+N22)+F(I7+N12)-F(I7+N11)-F(I7+1) 3 +F(I7+N21)+F(I7+N2)-F(I7+N1)-F(I7))/HZ) 38 CONTINUE DO 382 J=1,NY1 L1=NP1*(J-1)+NP2*(K-1) L2=N1*(J-1)+K8 DO 384 I=1,NP1 J9=I+L2 I1=I+L1 P(J9+K2)=P(J9+K2)-PP(I1+2*NP3) P(J9+N2)=P(J9+N2)-PP(I1+NP3) P(J9)=P(J9)-PP(I1) 384 CONTINUE DO 382 I=2,NP1 J9=I+NP1-1+L2 I1=NP1-I+1+L1 P(J9+K2)=P(J9+K2)+PP(I1+2*NP3) P(J9+N2)=P(J9+N2)+PP(I1+NP3) P(J9)=P(J9)-PP(I1) 382 CONTINUE C DO 381 M=1,NB DO 381 J=1,NY1 L1=N1*(J-1)+N2*(K-1)+N3*(M-1) L2=N1*(J-1)+N2*(M-1) L3=N1*(J-1)+N2+K3*(M-1) DO 381 I=1,NX1 I1=I+L1 J1=I+L2 J2=I+L3 P(J1)=0.125*(F(I1)+F(I1+1)+F(I1+N1)+F(I1+N11) 1 +F(I1+N2)+F(I1+N12)+F(I1+N21)+F(I1+N22)) U(J2)=P(J1) 381 CONTINUE C C FIRST STEP XX4=0.5*HX*FLOAT(2*NXP-NX1-1) Z=0.5*HZ*FLOAT(2*K-2) DO 50 J=1,NY1 L1=N1*(J-1)+N2*(K-1) L2=N1*(J-1) L3=L2+N2 Y=0.5*HY*FLOAT(2*J-2) C DO 501 I=1,NX1 X=HHX(I)+XX4 RO2=X*X+Y*Y+Z*Z RO2=MAX(RO2,AR2) RO3=RO2*SQRT(RO2) I1=I+L1 I2=I1+N3 I3=I2+N3 I4=I3+N3 I5=I4+N3 I6=I5+N3 I7=I6+N3 I8=I7+N3 J1=I+L2 J2=J1+N2 J3=J2+N2 J4=J3+N2 J5=J4+N2 J6=J5+N2 J7=J6+N2 J8=J7+N2 J9=J8+N2 J10=J9+N2 J11=J10+N2 JJ1=I+L3 C U(JJ1+KK7)=P(J8) 1 +DX1*(F(I4+N22)*F(I6+N22)-F(I2+N22)*F(I8+N22) 2 +F(I4+N12)*F(I6+N12)-F(I2+N12)*F(I8+N12) 3 +F(I4+N11)*F(I6+N11)-F(I2+N11)*F(I8+N11) 4 +F(I4+1)*F(I6+1)-F(I2+1)*F(I8+1) 5 -F(I4+N21)*F(I6+N21)+F(I2+N21)*F(I8+N21) 6 -F(I4+N2)*F(I6+N2)+F(I2+N2)*F(I8+N2) 7 -F(I4+N1)*F(I6+N1)+F(I2+N1)*F(I8+N1) 8 -F(I4)*F(I6)+F(I2)*F(I8)) 1 -DY1*(F(I3+N22)*F(I8+N22)-F(I4+N22)*F(I7+N22) 2 -F(I3+N12)*F(I8+N12)+F(I4+N12)*F(I7+N12) 3 +F(I3+N11)*F(I8+N11)-F(I4+N11)*F(I7+N11) 4 -F(I3+1)*F(I8+1)+F(I4+1)*F(I7+1) 5 +F(I3+N21)*F(I8+N21)-F(I4+N21)*F(I7+N21) 6 -F(I3+N2)*F(I8+N2)+F(I4+N2)*F(I7+N2) 7 +F(I3+N1)*F(I8+N1)-F(I4+N1)*F(I7+N1) 8 -F(I3)*F(I8)+F(I4)*F(I7)) U(JJ1+KK6)=P(J7) 1 +DZ1*(F(I3+N22)*F(I8+N22)-F(I4+N22)*F(I7+N22) 2 +F(I3+N12)*F(I8+N12)-F(I4+N12)*F(I7+N12) 3 -F(I3+N11)*F(I8+N11)+F(I4+N11)*F(I7+N11) 4 -F(I3+1)*F(I8+1)+F(I4+1)*F(I7+1) 5 +F(I3+N21)*F(I8+N21)-F(I4+N21)*F(I7+N21) 6 +F(I3+N2)*F(I8+N2)-F(I4+N2)*F(I7+N2) 7 -F(I3+N1)*F(I8+N1)+F(I4+N1)*F(I7+N1) 8 -F(I3)*F(I8)+F(I4)*F(I7)) 1 -DX1*(F(I2+N22)*F(I7+N22)-F(I3+N22)*F(I6+N22) 2 +F(I2+N12)*F(I7+N12)-F(I3+N12)*F(I6+N12) 3 +F(I2+N11)*F(I7+N11)-F(I3+N11)*F(I6+N11) 4 +F(I2+1)*F(I7+1)-F(I3+1)*F(I6+1) 5 -F(I2+N21)*F(I7+N21)+F(I3+N21)*F(I6+N21) 6 -F(I2+N2)*F(I7+N2)+F(I3+N2)*F(I6+N2) 7 -F(I2+N1)*F(I7+N1)+F(I3+N1)*F(I6+N1) 8 -F(I2)*F(I7)+F(I3)*F(I6)) U(JJ1+KK5)=P(J6) 1 +DY1*(F(I2+N22)*F(I7+N22)-F(I3+N22)*F(I6+N22) 2 -F(I2+N12)*F(I7+N12)+F(I3+N12)*F(I6+N12) 3 +F(I2+N11)*F(I7+N11)-F(I3+N11)*F(I6+N11) 4 -F(I2+1)*F(I7+1)+F(I3+1)*F(I6+1) 5 +F(I2+N21)*F(I7+N21)-F(I3+N21)*F(I6+N21) 6 -F(I2+N2)*F(I7+N2)+F(I3+N2)*F(I6+N2) 7 +F(I2+N1)*F(I7+N1)-F(I3+N1)*F(I6+N1) 8 -F(I2)*F(I7)+F(I3)*F(I6)) 1 -DZ1*(F(I4+N22)*F(I6+N22)-F(I2+N22)*F(I8+N22) 2 +F(I4+N12)*F(I6+N12)-F(I2+N12)*F(I8+N12) 3 -F(I4+N11)*F(I6+N11)+F(I2+N11)*F(I8+N11) 4 -F(I4+1)*F(I6+1)+F(I2+1)*F(I8+1) 5 +F(I4+N21)*F(I6+N21)-F(I2+N21)*F(I8+N21) 6 +F(I4+N2)*F(I6+N2)-F(I2+N2)*F(I8+N2) 7 -F(I4+N1)*F(I6+N1)+F(I2+N1)*F(I8+N1) 8 -F(I4)*F(I6)+F(I2)*F(I8)) C X3=DX1*(F(I2+N22)+F(I2+N12)+F(I2+N11)+F(I2+1) 1 -F(I2+N21)-F(I2+N2)-F(I2+N1)-F(I2)) 2 +DY1*(F(I3+N22)-F(I3+N12)+F(I3+N11)-F(I3+1) 3 +F(I3+N21)-F(I3+N2)+F(I3+N1)-F(I3)) 4 +DZ1*(F(I4+N22)+F(I4+N12)-F(I4+N11)-F(I4+1) 5 +F(I4+N21)+F(I4+N2)-F(I4+N1)-F(I4)) U(JJ1+KK4)=P(J5)-P(J5)*GAM*X3 1 -DX1*P(J2)*(F(I5+N22)+F(I5+N12)+F(I5+N11)+F(I5+1) 2 -F(I5+N21)-F(I5+N2)-F(I5+N1)-F(I5)) 3 -DY1*P(J3)*(F(I5+N22)-F(I5+N12)+F(I5+N11)-F(I5+1) 4 +F(I5+N21)-F(I5+N2)+F(I5+N1)-F(I5)) 5 -DZ1*P(J4)*(F(I5+N22)+F(I5+N12)-F(I5+N11)-F(I5+1) 6 +F(I5+N21)+F(I5+N2)-F(I5+N1)-F(I5)) U(JJ1+KK3)=P(J4)+T1*((P(J9)*P(J7)-P(J10)*P(J6))/P(J1) 1 -EUD*P(J4)-GRA*Z/RO3) 2 -DX1*P(J2)*(F(I4+N22)+F(I4+N12)+F(I4+N11)+F(I4+1) 3 -F(I4+N21)-F(I4+N2)-F(I4+N1)-F(I4)) 4 -DY1*P(J3)*(F(I4+N22)-F(I4+N12)+F(I4+N11)-F(I4+1) 5 +F(I4+N21)-F(I4+N2)+F(I4+N1)-F(I4)) 6 -DZ1*P(J4)*(F(I4+N22)+F(I4+N12)-F(I4+N11)-F(I4+1) 7 +F(I4+N21)+F(I4+N2)-F(I4+N1)-F(I4)) 8 -DZ1*(F(I5+N22)+F(I5+N12)-F(I5+N11)-F(I5+1) 9 +F(I5+N21)+F(I5+N2)-F(I5+N1)-F(I5))/P(J1) U(JJ1+KK2)=P(J3)+T1*((P(J11)*P(J6)-P(J9)*P(J8))/P(J1) 1 -EUD*P(J3)-GRA*Y/RO3) 2 -DX1*P(J2)*(F(I3+N22)+F(I3+N12)+F(I3+N11)+F(I3+1) 3 -F(I3+N21)-F(I3+N2)-F(I3+N1)-F(I3)) 4 -DY1*P(J3)*(F(I3+N22)-F(I3+N12)+F(I3+N11)-F(I3+1) 5 +F(I3+N21)-F(I3+N2)+F(I3+N1)-F(I3)) 6 -DZ1*P(J4)*(F(I3+N22)+F(I3+N12)-F(I3+N11)-F(I3+1) 7 +F(I3+N21)+F(I3+N2)-F(I3+N1)-F(I3)) 8 -DY1*(F(I5+N22)-F(I5+N12)+F(I5+N11)-F(I5+1) 9 +F(I5+N21)-F(I5+N2)+F(I5+N1)-F(I5))/P(J1) U(JJ1+K3)=P(J2)+T1*((P(J10)*P(J8)-P(J11)*P(J7))/P(J1) 1 -EUD*P(J2)-GRA*X/RO3) 2 -DX1*P(J2)*(F(I2+N22)+F(I2+N12)+F(I2+N11)+F(I2+1) 3 -F(I2+N21)-F(I2+N2)-F(I2+N1)-F(I2)) 4 -DY1*P(J3)*(F(I2+N22)-F(I2+N12)+F(I2+N11)-F(I2+1) 5 +F(I2+N21)-F(I2+N2)+F(I2+N1)-F(I2)) 6 -DZ1*P(J4)*(F(I2+N22)+F(I2+N12)-F(I2+N11)-F(I2+1) 7 +F(I2+N21)+F(I2+N2)-F(I2+N1)-F(I2)) 8 -DX1*(F(I5+N22)+F(I5+N12)+F(I5+N11)+F(I5+1) 9 -F(I5+N21)-F(I5+N2)-F(I5+N1)-F(I5))/P(J1) U(JJ1)=P(J1) 1 -DX1*(F(I2+N22)*F(I1+N22)+F(I2+N12)*F(I1+N12) 2 +F(I2+N11)*F(I1+N11)+F(I2+1)*F(I1+1) 3 -F(I2+N21)*F(I1+N21)-F(I2+N2)*F(I1+N2) 4 -F(I2+N1)*F(I1+N1)-F(I2)*F(I1)) 5 -DY1*(F(I3+N22)*F(I1+N22)-F(I3+N12)*F(I1+N12) 6 +F(I3+N11)*F(I1+N11)-F(I3+1)*F(I1+1) 7 +F(I3+N21)*F(I1+N21)-F(I3+N2)*F(I1+N2) 8 +F(I3+N1)*F(I1+N1)-F(I3)*F(I1)) 1 -DZ1*(F(I4+N22)*F(I1+N22)+F(I4+N12)*F(I1+N12) 2 -F(I4+N11)*F(I1+N11)-F(I4+1)*F(I1+1) 3 +F(I4+N21)*F(I1+N21)+F(I4+N2)*F(I1+N2) 4 -F(I4+N1)*F(I1+N1)-F(I4)*F(I1)) 501 CONTINUE C DO 520 I=1,NX1 X=HHX(I)+XX4 RO2=X*X+Y*Y+Z*Z I1=I+L2 I2=I1+N2 I3=I2+N2 I4=I3+N2 I5=I4+N2 I6=I5+N2 I7=I6+N2 I8=I7+N2 J1=I+L3 X2=MAX(RO2/AR2-1.0,0.0) X2=ARU*X2*X2 X2=X2/(X2+1.0) U(J1+KK7)=U(J1+KK7)*X2+P(I8)*(1.0-X2) U(J1+KK6)=U(J1+KK6)*X2+P(I7)*(1.0-X2) U(J1+KK5)=U(J1+KK5)*X2+P(I6)*(1.0-X2) U(J1+KK4)=U(J1+KK4)*X2+P(I5)*(1.0-X2) U(J1+KK3)=U(J1+KK3)*X2+P(I4)*(1.0-X2) U(J1+KK2)=U(J1+KK2)*X2+P(I3)*(1.0-X2) U(J1+K3)=U(J1+K3)*X2+P(I2)*(1.0-X2) U(J1)=U(J1)*X2+P(I1)*(1.0-X2) 520 CONTINUE 50 CONTINUE C C SECOND STEP C DO 62 M=1,NB L1=N2*(K-1)+N3*(M-1) L2=N2+K3*(M-1) DO 62 I=1,N2 I1=I+L1 J1=I+L2 FF(J1+N2)=F(I1+N2) FF(J1)=F(I1) 62 CONTINUE C C CURRENT DO 68 J=2,NY1 L1=N1*(J-1)+N2*(K-1) L2=N1*(J-1) L3=L2+N2 DO 68 I=2,NX1 I6=I+L3+KK5 I7=I6+K3 I8=I7+K3 J9=I+L2+K8 P(J9+K2)=0.25*((U(I7)+U(I7-N1)+U(I7-N2)+U(I7-N21) 1 -U(I7-1)-U(I7-N11)-U(I7-N12)-U(I7-N22))/HX 2 -(U(I6)-U(I6-N1)+U(I6-N2)-U(I6-N21) 3 +U(I6-1)-U(I6-N11)+U(I6-N12)-U(I6-N22))/HY) P(J9+N2)=0.25*((U(I6)+U(I6-N1)-U(I6-N2)-U(I6-N21) 1 +U(I6-1)+U(I6-N11)-U(I6-N12)-U(I6-N22))/HZ 2 -(U(I8)+U(I8-N1)+U(I8-N2)+U(I8-N21) 3 -U(I8-1)-U(I8-N11)-U(I8-N12)-U(I8-N22))/HX) P(J9)=0.25*((U(I8)-U(I8-N1)+U(I8-N2)-U(I8-N21) 1 +U(I8-1)-U(I8-N11)+U(I8-N12)-U(I8-N22))/HY 2 -(U(I7)+U(I7-N1)-U(I7-N2)-U(I7-N21) 3 +U(I7-1)+U(I7-N11)-U(I7-N12)-U(I7-N22))/HZ) 68 CONTINUE DO 682 J=2,NY1 L1=NP1*(J-1)+NP2*(K-1) L2=N1*(J-1)+K8 DO 684 I=2,NP1 J9=I+L2 I1=I+L1 I2=I1+NP3 I3=I2+NP3 P(J9+K2)=P(J9+K2)-0.125*(PP(I3)+PP(I3-1)+PP(I3-NP1)+PP(I3-NP2) 1 +PP(I3-NP11)+PP(I3-NP12)+PP(I3-NP21)+PP(I3-NP22)) P(J9+N2)=P(J9+N2)-0.125*(PP(I2)+PP(I2-1)+PP(I2-NP1)+PP(I2-NP2) 1 +PP(I2-NP11)+PP(I2-NP12)+PP(I2-NP21)+PP(I2-NP22)) P(J9)=P(J9)-0.125*(PP(I1)+PP(I1-1)+PP(I1-NP1)+PP(I1-NP2) 1 +PP(I1-NP11)+PP(I1-NP12)+PP(I1-NP21)+PP(I1-NP22)) 684 CONTINUE DO 682 I=2,NP1 J9=I+NP1-1+L2 I1=NP1-I+2+L1 I2=I1+NP3 I3=I2+NP3 P(J9+K2)=P(J9+K2)+0.125*(PP(I3)+PP(I3-1)+PP(I3-NP1)+PP(I3-NP2) 1 +PP(I3-NP11)+PP(I3-NP12)+PP(I3-NP21)+PP(I3-NP22)) P(J9+N2)=P(J9+N2)+0.125*(PP(I2)+PP(I2-1)+PP(I2-NP1)+PP(I2-NP2) 1 +PP(I2-NP11)+PP(I2-NP12)+PP(I2-NP21)+PP(I2-NP22)) P(J9)=P(J9)-0.125*(PP(I1)+PP(I1-1)+PP(I1-NP1)+PP(I1-NP2) 1 +PP(I1-NP11)+PP(I1-NP12)+PP(I1-NP21)+PP(I1-NP22)) 682 CONTINUE C DO 681 M=1,NB DO 681 J=2,NY1 L2=N1*(J-1)+N2*(M-1) L3=N1*(J-1)+N2+K3*(M-1) DO 681 I=2,NX1 I1=I+L3 P(I+L2)=0.125*(U(I1)+U(I1-1)+U(I1-N1)+U(I1-N11) 1 +U(I1-N2)+U(I1-N12)+U(I1-N21)+U(I1-N22)) 681 CONTINUE C C SECOND STEP XX4=0.5*HX*FLOAT(2*NXP-NX1-2) Z=0.5*HZ*FLOAT(2*K-3) DO 80 J=2,NY1 L1=N1*(J-1)+N2*(K-1) L2=N1*(J-1) L3=L2+N2 Y=0.5*HY*FLOAT(2*J-3) C DO 801 I=2,NX1 X=HHX(I)+XX4 RO2=X*X+Y*Y+Z*Z RO2=MAX(RO2,AR2) RO3=RO2*SQRT(RO2) I1=I+L3 I2=I1+K3 I3=I2+K3 I4=I3+K3 I5=I4+K3 I6=I5+K3 I7=I6+K3 I8=I7+K3 J1=I+L2 J2=J1+N2 J3=J2+N2 J4=J3+N2 J5=J4+N2 J6=J5+N2 J7=J6+N2 J8=J7+N2 J9=J8+N2 J10=J9+N2 J11=J10+N2 JJ1=I+L1 C X1=FF(I5)/(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 X1=EAT*(DX4*(FF(I8+1)-2.0*FF(I8)+FF(I8-1)) 1 +DY4*(FF(I8+N1)-2.0*FF(I8)+FF(I8-N1)) 2 +DZ4*(FF(I8+N2)-2.0*FF(I8)+FF(I8-N2))) F(JJ1+M7)=FF(I8)+X1 1 +DX2*(U(I4)*U(I6)-U(I2)*U(I8) 2 +U(I4-N1)*U(I6-N1)-U(I2-N1)*U(I8-N1) 3 +U(I4-N2)*U(I6-N2)-U(I2-N2)*U(I8-N2) 4 +U(I4-N21)*U(I6-N21)-U(I2-N21)*U(I8-N21) 5 -U(I4-1)*U(I6-1)+U(I2-1)*U(I8-1) 6 -U(I4-N11)*U(I6-N11)+U(I2-N11)*U(I8-N11) 7 -U(I4-N12)*U(I6-N12)+U(I2-N12)*U(I8-N12) 8 -U(I4-N22)*U(I6-N22)+U(I2-N22)*U(I8-N22)) 1 -DY2*(U(I3)*U(I8)-U(I4)*U(I7) 2 -U(I3-N1)*U(I8-N1)+U(I4-N1)*U(I7-N1) 3 +U(I3-N2)*U(I8-N2)-U(I4-N2)*U(I7-N2) 4 -U(I3-N21)*U(I8-N21)+U(I4-N21)*U(I7-N21) 5 +U(I3-1)*U(I8-1)-U(I4-1)*U(I7-1) 6 -U(I3-N11)*U(I8-N11)+U(I4-N11)*U(I7-N11) 7 +U(I3-N12)*U(I8-N12)-U(I4-N12)*U(I7-N12) 8 -U(I3-N22)*U(I8-N22)+U(I4-N22)*U(I7-N22)) X1=EAT*(DX4*(FF(I7+1)-2.0*FF(I7)+FF(I7-1)) 1 +DY4*(FF(I7+N1)-2.0*FF(I7)+FF(I7-N1)) 2 +DZ4*(FF(I7+N2)-2.0*FF(I7)+FF(I7-N2))) F(JJ1+M6)=FF(I7)+X1 1 +DZ2*(U(I3)*U(I8)-U(I4)*U(I7) 2 +U(I3-N1)*U(I8-N1)-U(I4-N1)*U(I7-N1) 3 -U(I3-N2)*U(I8-N2)+U(I4-N2)*U(I7-N2) 4 -U(I3-N21)*U(I8-N21)+U(I4-N21)*U(I7-N21) 5 +U(I3-1)*U(I8-1)-U(I4-1)*U(I7-1) 6 +U(I3-N11)*U(I8-N11)-U(I4-N11)*U(I7-N11) 7 -U(I3-N12)*U(I8-N12)+U(I4-N12)*U(I7-N12) 8 -U(I3-N22)*U(I8-N22)+U(I4-N22)*U(I7-N22)) 1 -DX2*(U(I2)*U(I7)-U(I3)*U(I6) 2 +U(I2-N1)*U(I7-N1)-U(I3-N1)*U(I6-N1) 3 +U(I2-N2)*U(I7-N2)-U(I3-N2)*U(I6-N2) 4 +U(I2-N21)*U(I7-N21)-U(I3-N21)*U(I6-N21) 5 -U(I2-1)*U(I7-1)+U(I3-1)*U(I6-1) 6 -U(I2-N11)*U(I7-N11)+U(I3-N11)*U(I6-N11) 7 -U(I2-N12)*U(I7-N12)+U(I3-N12)*U(I6-N12) 8 -U(I2-N22)*U(I7-N22)+U(I3-N22)*U(I6-N22)) X1=EAT*(DX4*(FF(I6+1)-2.0*FF(I6)+FF(I6-1)) 1 +DY4*(FF(I6+N1)-2.0*FF(I6)+FF(I6-N1)) 2 +DZ4*(FF(I6+N2)-2.0*FF(I6)+FF(I6-N2))) F(JJ1+M5)=FF(I6)+X1 1 +DY2*(U(I2)*U(I7)-U(I3)*U(I6) 2 -U(I2-N1)*U(I7-N1)+U(I3-N1)*U(I6-N1) 3 +U(I2-N2)*U(I7-N2)-U(I3-N2)*U(I6-N2) 4 -U(I2-N21)*U(I7-N21)+U(I3-N21)*U(I6-N21) 5 +U(I2-1)*U(I7-1)-U(I3-1)*U(I6-1) 6 -U(I2-N11)*U(I7-N11)+U(I3-N11)*U(I6-N11) 7 +U(I2-N12)*U(I7-N12)-U(I3-N12)*U(I6-N12) 8 -U(I2-N22)*U(I7-N22)+U(I3-N22)*U(I6-N22)) 1 -DZ2*(U(I4)*U(I6)-U(I2)*U(I8) 2 +U(I4-N1)*U(I6-N1)-U(I2-N1)*U(I8-N1) 3 -U(I4-N2)*U(I6-N2)+U(I2-N2)*U(I8-N2) 4 -U(I4-N21)*U(I6-N21)+U(I2-N21)*U(I8-N21) 5 +U(I4-1)*U(I6-1)-U(I2-1)*U(I8-1) 6 +U(I4-N11)*U(I6-N11)-U(I2-N11)*U(I8-N11) 7 -U(I4-N12)*U(I6-N12)+U(I2-N12)*U(I8-N12) 8 -U(I4-N22)*U(I6-N22)+U(I2-N22)*U(I8-N22)) C X3=DX2*(U(I2)+U(I2-N1)+U(I2-N2)+U(I2-N21) 1 -U(I2-1)-U(I2-N11)-U(I2-N12)-U(I2-N22)) 2 +DY2*(U(I3)-U(I3-N1)+U(I3-N2)-U(I3-N21) 3 +U(I3-1)-U(I3-N11)+U(I3-N12)-U(I3-N22)) 4 +DZ2*(U(I4)+U(I4-N1)-U(I4-N2)-U(I4-N21) 5 +U(I4-1)+U(I4-N11)-U(I4-N12)-U(I4-N22)) F(JJ1+M4)=FF(I5)-P(J5)*GAM*X3 1 -DX2*P(J2)*(U(I5)+U(I5-N1)+U(I5-N2)+U(I5-N21) 2 -U(I5-1)-U(I5-N11)-U(I5-N12)-U(I5-N22)) 3 -DY2*P(J3)*(U(I5)-U(I5-N1)+U(I5-N2)-U(I5-N21) 4 +U(I5-1)-U(I5-N11)+U(I5-N12)-U(I5-N22)) 5 -DZ2*P(J4)*(U(I5)+U(I5-N1)-U(I5-N2)-U(I5-N21) 6 +U(I5-1)+U(I5-N11)-U(I5-N12)-U(I5-N22)) 7 +PMU*(DX4*(FF(I5+1)-2.0*FF(I5)+FF(I5-1)) 8 +DY4*(FF(I5+N1)-2.0*FF(I5)+FF(I5-N1)) 9 +DZ4*(FF(I5+N2)-2.0*FF(I5)+FF(I5-N2))) F(JJ1+M3)=FF(I4)+T*((P(J9)*P(J7)-P(J10)*P(J6))/P(J1) 1 -EUD*P(J4)-GRA*Z/RO3) 2 -DX2*P(J2)*(U(I4)+U(I4-N1)+U(I4-N2)+U(I4-N21) 3 -U(I4-1)-U(I4-N11)-U(I4-N12)-U(I4-N22)) 4 -DY2*P(J3)*(U(I4)-U(I4-N1)+U(I4-N2)-U(I4-N21) 5 +U(I4-1)-U(I4-N11)+U(I4-N12)-U(I4-N22)) 6 -DZ2*P(J4)*(U(I4)+U(I4-N1)-U(I4-N2)-U(I4-N21) 7 +U(I4-1)+U(I4-N11)-U(I4-N12)-U(I4-N22)) 8 -DZ2*(U(I5)+U(I5-N1)-U(I5-N2)-U(I5-N21) 9 +U(I5-1)+U(I5-N11)-U(I5-N12)-U(I5-N22))/P(J1) 1 +RMU*(DX4*(FF(I4+1)-2.0*FF(I4)+FF(I4-1)) 2 +DY4*(FF(I4+N1)-2.0*FF(I4)+FF(I4-N1)) 3 +DZ4*(FF(I4+N2)-2.0*FF(I4)+FF(I4-N2)))/FF(I1) F(JJ1+M2)=FF(I3)+T*((P(J11)*P(J6)-P(J9)*P(J8))/P(J1) 1 -EUD*P(J3)-GRA*Y/RO3) 2 -DX2*P(J2)*(U(I3)+U(I3-N1)+U(I3-N2)+U(I3-N21) 3 -U(I3-1)-U(I3-N11)-U(I3-N12)-U(I3-N22)) 4 -DY2*P(J3)*(U(I3)-U(I3-N1)+U(I3-N2)-U(I3-N21) 5 +U(I3-1)-U(I3-N11)+U(I3-N12)-U(I3-N22)) 6 -DZ2*P(J4)*(U(I3)+U(I3-N1)-U(I3-N2)-U(I3-N21) 7 +U(I3-1)+U(I3-N11)-U(I3-N12)-U(I3-N22)) 8 -DY2*(U(I5)-U(I5-N1)+U(I5-N2)-U(I5-N21) 9 +U(I5-1)-U(I5-N11)+U(I5-N12)-U(I5-N22))/P(J1) 1 +RMU*(DX4*(FF(I3+1)-2.0*FF(I3)+FF(I3-1)) 2 +DY4*(FF(I3+N1)-2.0*FF(I3)+FF(I3-N1)) 3 +DZ4*(FF(I3+N2)-2.0*FF(I3)+FF(I3-N2)))/FF(I1) F(JJ1+N3)=FF(I2)+T*((P(J10)*P(J8)-P(J11)*P(J7))/P(J1) 1 -EUD*P(J2)-GRA*X/RO3) 2 -DX2*P(J2)*(U(I2)+U(I2-N1)+U(I2-N2)+U(I2-N21) 3 -U(I2-1)-U(I2-N11)-U(I2-N12)-U(I2-N22)) 4 -DY2*P(J3)*(U(I2)-U(I2-N1)+U(I2-N2)-U(I2-N21) 5 +U(I2-1)-U(I2-N11)+U(I2-N12)-U(I2-N22)) 6 -DZ2*P(J4)*(U(I2)+U(I2-N1)-U(I2-N2)-U(I2-N21) 7 +U(I2-1)+U(I2-N11)-U(I2-N12)-U(I2-N22)) 8 -DX2*(U(I5)+U(I5-N1)+U(I5-N2)+U(I5-N21) 9 -U(I5-1)-U(I5-N11)-U(I5-N12)-U(I5-N22))/P(J1) 1 +RMU*(DX4*(FF(I2+1)-2.0*FF(I2)+FF(I2-1)) 2 +DY4*(FF(I2+N1)-2.0*FF(I2)+FF(I2-N1)) 3 +DZ4*(FF(I2+N2)-2.0*FF(I2)+FF(I2-N2)))/FF(I1) F(JJ1)=FF(I1) 1 -DX2*(U(I2)*U(I1)+U(I2-N1)*U(I1-N1) 2 +U(I2-N2)*U(I1-N2)+U(I2-N21)*U(I1-N21) 3 -U(I2-1)*U(I1-1)-U(I2-N11)*U(I1-N11) 4 -U(I2-N12)*U(I1-N12)-U(I2-N22)*U(I1-N22)) 5 -DY2*(U(I3)*U(I1)-U(I3-N1)*U(I1-N1) 6 +U(I3-N2)*U(I1-N2)-U(I3-N21)*U(I1-N21) 7 +U(I3-1)*U(I1-1)-U(I3-N11)*U(I1-N11) 8 +U(I3-N12)*U(I1-N12)-U(I3-N22)*U(I1-N22)) 1 -DZ2*(U(I4)*U(I1)+U(I4-N1)*U(I1-N1) 2 -U(I4-N2)*U(I1-N2)-U(I4-N21)*U(I1-N21) 3 +U(I4-1)*U(I1-1)+U(I4-N11)*U(I1-N11) 4 -U(I4-N12)*U(I1-N12)-U(I4-N22)*U(I1-N22)) 5 +DFU*(DX4*(FF(I1+1)-2.0*FF(I1)+FF(I1-1)) 6 +DY4*(FF(I1+N1)-2.0*FF(I1)+FF(I1-N1)) 7 +DZ4*(FF(I1+N2)-2.0*FF(I1)+FF(I1-N2))) 801 CONTINUE C DO 820 I=2,NX1 X=HHX(I)+XX4 RO2=X*X+Y*Y+Z*Z I1=I+L3 I2=I1+K3 I3=I2+K3 I4=I3+K3 I5=I4+K3 I6=I5+K3 I7=I6+K3 I8=I7+K3 J1=I+L1 X2=MAX(RO2/AR2-1.0,0.0) X2=ARU*X2*X2 X2=X2/(X2+1.0) F(J1+M7)=F(J1+M7)*X2+FF(I8)*(1.0-X2) F(J1+M6)=F(J1+M6)*X2+FF(I7)*(1.0-X2) F(J1+M5)=F(J1+M5)*X2+FF(I6)*(1.0-X2) F(J1+M4)=F(J1+M4)*X2+FF(I5)*(1.0-X2) F(J1+M3)=F(J1+M3)*X2+FF(I4)*(1.0-X2) F(J1+M2)=F(J1+M2)*X2+FF(I3)*(1.0-X2) F(J1+N3)=F(J1+N3)*X2+FF(I2)*(1.0-X2) F(J1)=F(J1)*X2+FF(I1)*(1.0-X2) 820 CONTINUE C DO 822 I=2,NX1 I1=I+L1 I2=I1+N3 I3=I2+N3 I4=I3+N3 I5=I4+N3 P(I+K4)=F(I5) P(I+K3)=F(I4) P(I+K2)=F(I3) P(I+N2)=F(I2) P(I)=F(I1) 822 CONTINUE C DO 824 I=2,NX1 I1=I+L1 X1=P(I) X5=P(I+K4) Y1=AMAX1(X1,RO02) Y5=AMAX1(X5,PR02) X2=P(I+N2)*P(I+N2)+P(I+K2)*P(I+K2)+P(I+K3)*P(I+K3) X2=AMAX1(X2,1.0E-14) X=(X1*X2+GM1*(X5-Y5))/(Y1*X2) X=AMAX1(X,0.0) X=SQRT(X) F(I1+M4)=Y5 F(I1+M3)=P(I+K3)*X F(I1+M2)=P(I+K2)*X F(I1+N3)=P(I+N2)*X F(I1)=Y1 824 CONTINUE 80 CONTINUE C DO 830 M=1,NB L1=K3*(M-1) DO 830 I=1,N2 I1=I+L1 U(I1)=U(I1+N2) FF(I1)=FF(I1+N2) 830 CONTINUE 90 CONTINUE C C TIME=TIME+T C IF(IIQ.NE.IIQ0) GO TO 401 IIQ=0 C X1=2.0 DO 102 K=1,NZ2 DO 102 J=1,NY2 L1=N1*(J-1)+N2*(K-1) DO 102 I=1,NX2 X=(2.0*HHX(I)/HX-FLOAT(1+NX2))/FLOAT(NX1) X=(X1+1.0)*X+X1 C X=CVMGP(X*X,0.0,-X) IF(X.GT.0.0) THEN X=0.0 ELSE X=X*X ENDIF I1=I+L1 F(I1+M7)=(1.0-X)*F(I1+M7)+BIS*X F(I1+M6)=(1.0-X)*F(I1+M6) F(I1+M5)=(1.0-X)*F(I1+M5) F(I1+M4)=(1.0-X)*F(I1+M4)+PR01*X F(I1+M3)=(1.0-X)*F(I1+M3) F(I1+M2)=(1.0-X)*F(I1+M2) F(I1+N3)=(1.0-X)*F(I1+N3)+VSW*X F(I1)=(1.0-X)*F(I1)+RO01*X 102 CONTINUE C 401 IF(IIP.NE.IIP0) GO TO 100 IIP=0 C WRITE(12) F WRITE(6,661) II,LAST,NIN,ITAP 661 FORMAT(1H ,5X,4I10) C JYG=NY1+NZ1 JYG2=JYG-2 JXG=N1*JYG JXN=NX2 LANK=36 IP=0 LB=1 LL=8 X=2.0*AR2 DO 250 M=1,LL I3=M IF(M.GT.NB) I3=1 DO 250 I=1,NX2 DO 252 K=2,NZ2 I1=I+N2*(K-1)+N3*(I3-1) I2=I+N1*(NY1+K-2)+JXG*(M-1) FF(I2)=0.5*(F(I1)+F(I1+N1)) 252 CONTINUE DO 250 J=2,NY2 I1=I+N1*(J-1)+N3*(I3-1) I2=I+N1*(NY1+1-J)+JXG*(M-1) FF(I2)=0.5*(F(I1)+F(I1+N2)) 250 CONTINUE I1=II C CALL SCALM(NX,JYG2,NXP,LB,LL,LLL,HX,HZ,X,ZMIN,ZMAX,FF) LL1=2 LL2=LL/LL1 DO 150 J=1,LL2 LB=1+LL1*(J-1) C CALL GRAPD(JXG,JYG,LB,LL1,ZMIN,ZMAX,I1,LLL,JXN,LANK,IP,FF) 150 CONTINUE C JYG=NZ2 JYG2=JYG-2 JXG=NY2*JYG JXN=NY2 LANK=36 ARR1=0.1 IP=0 DO 260 IL=1,NIL LB=1 LL=8 I=NX2/2+1+(NX2*IL)/(2*(NIL+1)) DO 262 M=1,LL I3=M IF(M.GT.NB) I3=1 DO 262 K=1,NZ2 DO 262 J=1,NY2 I1=I+N1*(J-1)+N2*(K-1)+N3*(I3-1) I2=J+NY2*(K-1)+JXG*(M-1) FF(I2)=F(I1) 262 CONTINUE I1=I C CALL SCALM(NY,JYG2,NXP,LB,LL,LLL,HX,HZ,ARR1,ZMIN,ZMAX,FF) LL1=4 LL2=LL/LL1 DO 260 J=1,LL2 LB=1+LL1*(J-1) C CALL GRAPD(JXG,JYG,LB,LL1,ZMIN,ZMAX,I1,LLL,JXN,LANK,IP,FF) 260 CONTINUE C 100 CONTINUE 300 CONTINUE 9 CONTINUE REWIND 12 STOP END SUBROUTINE EQUIB7(NX,NY,NZ,NXP,RO01,PR01,RRAT,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: AID=CJ(7) AR1=CJ(8) GAM=CJ(9) GRA=CJ(10) P0=(GAM-1.0)*GRA/GAM BIS=CP(11)*1.0E-4 X1=1.0/(GAM-1.0) X2=GAM*X1 B0=CJ(4) RO02=RRAT*RO01 PR02=PR01 C: IC0=4 ICD=6 PI=3.1415926 NX1=NX+1 NX2=NX+2 NY1=NY+1 NY2=NY+2 NZ1=NZ+1 NZ2=NZ+2 N1=NX2 N2=N1*NY2 N3=N2*NZ2 N11=N1+1 N12=N2+1 N21=N1+N2 N22=N21+1 NB=8 HX=CP(1)/FLOAT(NX1) HY=CP(2)/FLOAT(NY1) HZ=CP(3)/FLOAT(NZ1) VSW=CP(9) C: DO 10 K=1,NZ2 DO 10 J=1,NY2 VR1=VSW DO 10 I=1,NX2 X=0.5*HX*FLOAT(2*I-NX2-1+2*NXP) Y=0.5*HY*FLOAT(2*J-3) Z=0.5*HZ*FLOAT(2*K-3) RO2=X*X+Y*Y+Z*Z RO1=SQRT(RO2) RO3=RO1*RO2 I1=I+N1*(J-1)+N2*(K-1) I2=I1+N3 I3=I2+N3 I4=I3+N3 I5=I4+N3 I6=I5+N3 I7=I6+N3 I8=I7+N3 F(I1)=1.0/RO3 C F(I1)=RO1**(-X1) IF(F(I1).LT.RO02) F(I1)=RO02 F(I5)=1.0E-7*CP(6)/RO2 C F(I5)=P0*RO1**(-X2) IF(F(I5).LT.PR02) F(I5)=PR02 F(I6)=-2.0*B0*X*Z/(RO2*RO2) F(I7)=-2.0*B0*Y*Z/(RO2*RO2) F(I8)=B0*(X*X+Y*Y-Z*Z)/(RO2*RO2) F(I6)=-3.0*B0*X*Z/(RO2*RO3) F(I7)=-3.0*B0*Y*Z/(RO2*RO3) F(I8)=B0*(X*X+Y*Y-2.0*Z*Z)/(RO2*RO3) PX1=SQRT(F(I6)*F(I6)+F(I7)*F(I7)+F(I8)*F(I8)) IF(ABS(VSW).LT.1.0E-8) GO TO 10 Y1=(RO01-CJ(1)*(PX1/VSW)**2)/RO01 IF(Y1.LE.0.0) Y1=0.0 Y1=VSW*SQRT(Y1) VR1=AMIN1(VR1,Y1) F(I2)=VR1 XX=VR1/VSW F(I1)=F(I1)+(RO01-RO02)*XX F(I5)=F(I5)+(PR01-PR02)*XX F(I8)=F(I8)+BIS 10 CONTINUE C NP1=NX2/2 NP2=NP1*NY1 NP3=NP2*NZ1 C CURRENT DO 38 K=1,NZ1 DO 38 J=1,NY1 L1=N1*(J-1)+N2*(K-1)+N3*5 L2=NP1*(J-1)+NP2*(K-1) DO 38 I=1,NP1 I6=I+L1 I7=I6+N3 I8=I7+N3 J9=I+L2 P(J9+2*NP3)=0.25*((F(I7+N22)+F(I7+N12)+F(I7+N11)+F(I7+1) 1 -F(I7+N21)-F(I7+N2)-F(I7+N1)-F(I7))/HX 2 -(F(I6+N22)-F(I6+N12)+F(I6+N11)-F(I6+1) 3 +F(I6+N21)-F(I6+N2)+F(I6+N1)-F(I6))/HY) P(J9+NP3)=0.25*((F(I6+N22)+F(I6+N12)-F(I6+N11)-F(I6+1) 1 +F(I6+N21)+F(I6+N2)-F(I6+N1)-F(I6))/HZ 2 -(F(I8+N22)+F(I8+N12)+F(I8+N11)+F(I8+1) 3 -F(I8+N21)-F(I8+N2)-F(I8+N1)-F(I8))/HX) P(J9)=0.25*((F(I8+N22)-F(I8+N12)+F(I8+N11)-F(I8+1) 1 +F(I8+N21)-F(I8+N2)+F(I8+N1)-F(I8))/HY 2 -(F(I7+N22)+F(I7+N12)-F(I7+N11)-F(I7+1) 3 +F(I7+N21)+F(I7+N2)-F(I7+N1)-F(I7))/HZ) 38 CONTINUE RETURN END SUBROUTINE CLOCK(IT) DOUBLE PRECISION IT X=0.0 Y=SECNDS(X) IT=1.00000*Y RETURN END