C ***** File = iso3.f ******************************************** C * * C * f77 -r8 -o iso3a iso3.f iso3a.f isofft.f * C * or f77 -r8 -o iso3b iso3.f iso3b.f isofft.f * C * -------------------------------------------------------------- * C * Decay of Isotropic Turbulence * C * - Staggered Grid * C * * Nonlinear - Convection Form * C * Interpolation method with iso3a.f * C * Direct method with iso3b.f * C * - 4 point Finite Difference Method * C * - Adams-Bashforth Method * C * * Viscous - Laplacian Form * C * - Crank-Nicolson Method * C * * Pressure - Fractional Step Method * C ****************************************************************** PARAMETER (N=32) COMMON /STEP/ISTEP,TSTEP /RE/RE /DT/DT COMMON /ERRP/ERRP /ERRU/ERRU /ERRV/ERRV /ERRW/ERRW COMMON /ITRP/ITRP /ITRU/ITRU /ITRV/ITRV /ITRW/ITRW RE=41.1 ISTOP=300 DT=0.01 CALL SBMESH OPEN(10,file='iso-0.d',status='old',form='unformatted') CALL SBREAD CALL SBCOUR(COUR) CALL SBENGY(EK) CALL SBDISS(DS) CALL SBEKII(EKII,EKIII) CALL SBDSII(DSII,DSIII,PRDS) CALL SBDSDS(DSDS) RET=RE*EK**2/DS BEK=EK BDS=DS ISTART=ISTEP TSTART=TSTEP ISTOP=ISTEP+ISTOP WRITE(6,1100) WRITE(6,1000) ISTEP,TSTEP,EK,DS,RET,0.,0. & ,1.E3*EKII,1.E5*EKIII,1.E3*DSII,1.E5*DSIII,PRDS,DSDS C & ,COUR,ERRD,ERRP,ERRU,ERRV,ERRW 100 ISTEP=ISTEP+1 IADV=ISTEP-ISTART TSTEP=TSTART+DT*IADV CALL SBVISC CALL SBNONL CALL SBSORU CALL SBRHSP CALL SBSORP CALL SBCORR CALL SBENGY(EK) CALL SBDISS(DS) RET=RE*EK**2/DS CALL SBDIVU(ERRD) IF(MOD(ISTEP,5).ge.0.or.ISTEP.le.10) THEN CALL SBCOUR(COUR) CALL SBVINT CALL SBEKII(EKII,EKIII) CALL SBDSII(DSII,DSIII,PRDS) CALL SBDSDS(DSDS) WRITE(6,1000) ISTEP,TSTEP,EK,DS,RET,(EK-BEK)/DT,(DS-BDS)/DT & ,1.E3*EKII,1.E5*EKIII,1.E3*DSII,1.E5*DSIII,PRDS,DSDS C & ,COUR,ERRD,ERRP,ERRU,ERRV,ERRW ENDIF IF(ISTEP.GE.ISTOP) GO TO 200 IF(MOD(ISTEP,100).eq.0) THEN CALL SBSPEC CALL SBVINT C SBVINT is neccesary for iso3b.f WRITE(6,1100) ENDIF 1000 FORMAT(I5,F7.3,2F11.6,F10.3,2F11.6,4F11.5,2F11.6 C &,F8.3,5E10.2 &) 1100 FORMAT(/1X,4HSTEP,6X,1HT,5X,6HENERGY,5X,6HDISSIP,7X,3HRET &,6X,5HDK/DT,6X,5HDE/DT,3X,8HKII(E-3),2X,9HKIII(E-5) &,3X,8HDII(E-3),2X,9HDIII(E-5),7X,4HPRDS,7X,4HDSDS C &,4X,4HCOUR,6X,4HERRD,6X,4HERRP,6X,4HERRU,6X,4HERRV,6X,4HERRW &/) BEK=EK BDS=DS GO TO 100 200 CONTINUE OPEN(12,file='iso-3b.d',status='new',form='unformatted') CALL SBSAVE OPEN(14,file='iso-pb.d',status='new',form='formatted') CALL SBPLOT CALL SBSOUK STOP END C ---------------------------------------------------------------------- C Mesh & FDM Parameters C ---------------------------------------------------------------------- SUBROUTINE SBMESH PARAMETER (N=32) COMMON /DX/DX /C/CM3,CM2,CM1,C00,CP1,CP2,CP3 /CWV/CWV(N) COMMON /LSTM/M1(N),M2(N),M3(N) /LSTP/L1(N),L2(N),L3(N) PI=ACOS(-1.) DX=2.*PI/REAL(N) DX24=24.*DX AM2= 1./DX24 AM1=-27./DX24 A00= 27./DX24 AP1= -1./DX24 BM1= 1./DX24 B00=-27./DX24 BP1= 27./DX24 BP2= -1./DX24 CM3=AM2*BM1 CM2=AM2*B00+AM1*BM1 CM1=AM2*BP1+AM1*B00+A00*BM1 C00=AM2*BP2+AM1*BP1+A00*B00+AP1*BM1 CP1= AM1*BP2+A00*BP1+AP1*B00 CP2= A00*BP2+AP1*BP1 CP3= AP1*BP2 P2N=2.*PI/REAL(N) DO 20 J=1,N JWAVE=J-1 IF(JWAVE.GE.N/2) JWAVE=JWAVE-N WAV1=REAL(MOD(JWAVE,N)) WAV2=REAL(MOD(JWAVE,N)*2) WAV3=REAL(MOD(JWAVE,N)*3) CWV(J)=C00+2.*(CM1*COS(P2N*WAV1) & +CM2*COS(P2N*WAV2)+CM3*COS(P2N*WAV3)) 20 CONTINUE DO 10 I=1,N M3(I)=I-3 M2(I)=I-2 M1(I)=I-1 IF(M3(I).LT.1) M3(I)=M3(I)+N IF(M2(I).LT.1) M2(I)=M2(I)+N IF(M1(I).LT.1) M1(I)=M1(I)+N L3(I)=I+3 L2(I)=I+2 L1(I)=I+1 IF(L3(I).GT.N) L3(I)=L3(I)-N IF(L2(I).GT.N) L2(I)=L2(I)-N IF(L1(I).GT.N) L1(I)=L1(I)-N 10 CONTINUE RETURN END C ---------------------------------------------------------------------- C Data Read C ---------------------------------------------------------------------- SUBROUTINE SBREAD PARAMETER (N=32) COMMON /STEP/ISTEP,TSTEP COMMON /U/U(N,N,N) /W1/UP(N,N,N) COMMON /V/V(N,N,N) /W2/VP(N,N,N) COMMON /W/W(N,N,N) /W3/WP(N,N,N) COMMON /LSTM/M1(N),M2(N),M3(N) /LSTP/L1(N),L2(N),L3(N) READ(10) ISTEP,TSTEP READ(10) UP READ(10) VP READ(10) WP CLOSE(10) write(6,1200) ISTEP,TSTEP 1200 format(/'Data READ: STEP=',I7,10X,'T=',F10.3/) SU=0. TU=0. SV=0. TV=0. SW=0. TW=0. FN=1./REAL(N**3) DO 10 K=1,N DO 10 J=1,N DO 10 I=1,N U(I,J,K)= & (-UP(M1(I),J,K)+9.*UP(I,J,K)+9.*UP(L1(I),J,K)-UP(L2(I),J,K))/16. V(I,J,K)= & (-VP(I,M1(J),K)+9.*VP(I,J,K)+9.*VP(I,L1(J),K)-VP(I,L2(J),K))/16. W(I,J,K)= & (-WP(I,J,M1(K))+9.*WP(I,J,K)+9.*WP(I,J,L1(K))-WP(I,J,L2(K)))/16. SU=SU+U(I,J,K)**2 TU=TU+UP(I,J,K)**2 SV=SV+V(I,J,K)**2 TV=TV+VP(I,J,K)**2 SW=SW+W(I,J,K)**2 TW=TW+WP(I,J,K)**2 10 CONTINUE SU=FN*SU TU=FN*TU SV=FN*SV TV=FN*TV SW=FN*SW TW=FN*TW WRITE(6,1000) SU,SV,SW,(SU+SV+SW)/2., TU,TV,TW,(TU+TV+TW)/2. 1000 FORMAT(/1H ,' Interpolated: Ku=',F12.8 &,' Kv=',F12.8,' Kw=',F12.8,' K =',F12.8 & /1H ,' Pure(initial): Ku=',F12.8 &,' Kv=',F12.8,' Kw=',F12.8,' K =',F12.8/) RETURN END C ---------------------------------------------------------------------- C Viscous Terms C ---------------------------------------------------------------------- SUBROUTINE SBVISC PARAMETER (N=32) COMMON /U/U(N,N,N) /UF/UF(N,N,N) COMMON /V/V(N,N,N) /VF/VF(N,N,N) COMMON /W/W(N,N,N) /WF/WF(N,N,N) COMMON /RE/RE /DX/DX /DT/DT COMMON /LSTM/M1(N),M2(N),M3(N) /LSTP/L1(N),L2(N),L3(N) VIS=1./(2.*RE *12.*DX**2) DO 10 K=1,N DO 10 J=1,N DO 10 I=1,N UF(I,J,K)=U(I,J,K)/DT+VIS*( & -U(I,J,M2(K))+16.*U(I,J,M1(K)) -U(I,M2(J),K)+16.*U(I,M1(J),K) & -U(M2(I),J,K)+16.*U(M1(I),J,K) -90.*U(I,J,K) & +16.*U(L1(I),J,K)-U(L2(I),J,K) +16.*U(I,L1(J),K)-U(I,L2(J),K) & +16.*U(I,J,L1(K))-U(I,J,L2(K)) ) VF(I,J,K)=V(I,J,K)/DT+VIS*( & -V(I,J,M2(K))+16.*V(I,J,M1(K)) -V(I,M2(J),K)+16.*V(I,M1(J),K) & -V(M2(I),J,K)+16.*V(M1(I),J,K) -90.*V(I,J,K) & +16.*V(L1(I),J,K)-V(L2(I),J,K) +16.*V(I,L1(J),K)-V(I,L2(J),K) & +16.*V(I,J,L1(K))-V(I,J,L2(K)) ) WF(I,J,K)=W(I,J,K)/DT+VIS*( & -W(I,J,M2(K))+16.*W(I,J,M1(K)) -W(I,M2(J),K)+16.*W(I,M1(J),K) & -W(M2(I),J,K)+16.*W(M1(I),J,K) -90.*W(I,J,K) & +16.*W(L1(I),J,K)-W(L2(I),J,K) +16.*W(I,L1(J),K)-W(I,L2(J),K) & +16.*W(I,J,L1(K))-W(I,J,L2(K)) ) 10 CONTINUE RETURN END C ---------------------------------------------------------------------- C Nonlinear Terms C ---------------------------------------------------------------------- SUBROUTINE SBNONL COMMON /STEP/ISTEP,TSTEP IF(ISTEP.GT.1) THEN AB= 3./2. BB=-1./2. ELSE AB=1. BB=0. ENDIF CALL SBNONU(AB,BB) CALL SBNONV(AB,BB) CALL SBNONW(AB,BB) RETURN END C ---------------------------------------------------------------------- C Poisson Equation for Velocity C ---------------------------------------------------------------------- SUBROUTINE SBSORU PARAMETER (N=32) COMMON /U/U(N,N,N) /UF/UF(N,N,N) /ERRU/ERRU /ITRU/ITRU COMMON /V/V(N,N,N) /VF/VF(N,N,N) /ERRV/ERRV /ITRV/ITRV COMMON /W/W(N,N,N) /WF/WF(N,N,N) /ERRW/ERRW /ITRW/ITRW COMMON /W1/A(N,N,N) /W2/B(N,N,N) /RE/RE RE2=2.*RE DO 10 K=1,N DO 10 J=1,N DO 10 I=1,N A(I,J,K)=U(I,J,K) B(I,J,K)=-RE2*UF(I,J,K) 10 CONTINUE CALL SBSORF(ERRU,ITRU) DO 15 K=1,N DO 15 J=1,N DO 15 I=1,N U(I,J,K)=A(I,J,K) 15 CONTINUE DO 20 K=1,N DO 20 J=1,N DO 20 I=1,N A(I,J,K)=V(I,J,K) B(I,J,K)=-RE2*VF(I,J,K) 20 CONTINUE CALL SBSORF(ERRV,ITRV) DO 25 K=1,N DO 25 J=1,N DO 25 I=1,N V(I,J,K)=A(I,J,K) 25 CONTINUE DO 30 K=1,N DO 30 J=1,N DO 30 I=1,N A(I,J,K)=W(I,J,K) B(I,J,K)=-RE2*WF(I,J,K) 30 CONTINUE CALL SBSORF(ERRW,ITRW) DO 35 K=1,N DO 35 J=1,N DO 35 I=1,N W(I,J,K)=A(I,J,K) 35 CONTINUE RETURN END C ---------------------------------------------------------------------- SUBROUTINE SBSORF(ERR,ITR) PARAMETER (N=32) COMMON /W1/A(N,N,N) /W2/B(N,N,N) COMMON /DX/DX /RE/RE /DT/DT COMMON /LSTM/M1(N),M2(N),M3(N) /LSTP/L1(N),L2(N),L3(N) DEL=2.*RE/DT TX12=1./(12.*DX**2) C0=90.*TX12+DEL SUM=0. DO 10 K=1,N DO 10 J=1,N DO 10 I=1,N 10 SUM=SUM+B(I,J,K)**2 C WRITE(6,1100) ITR=0 100 ITR=ITR+1 SUMR=0. DO 50 K=1,N DO 50 J=1,N DO 50 I=1,N SLAP=TX12*( & -A(I,J,M2(K))+16.*A(I,J,M1(K)) -A(I,M2(J),K)+16.*A(I,M1(J),K) & -A(M2(I),J,K)+16.*A(M1(I),J,K) -90.*A(I,J,K) & +16.*A(L1(I),J,K)-A(L2(I),J,K) +16.*A(I,L1(J),K)-A(I,L2(J),K) & +16.*A(I,J,L1(K))-A(I,J,L2(K)) ) RESI=SLAP-DEL*A(I,J,K)-B(I,J,K) A(I,J,K)=A(I,J,K)+1.5*RESI/C0 SUMR=SUMR+RESI**2 50 CONTINUE ERR=SQRT(SUMR/SUM) IF(ITR.LT.100.AND.ERR.GT.1.E-10) GO TO 100 C WRITE(6,1000) ITR,SUM,SUMR,ERR C 1000 FORMAT(I8,3E13.4) C 1100 FORMAT(/5X,3HITR,10X,3HSUM,9X,4HSUMR,10X,3HERR) RETURN END C ---------------------------------------------------------------------- C R.H.S. of Poisson Equation for Pressure C ---------------------------------------------------------------------- SUBROUTINE SBRHSP PARAMETER (N=32) COMMON /U/U(N,N,N) /V/V(N,N,N) /W/W(N,N,N) /W1/Q(N,N,N) COMMON /DX/DX /DT/DT COMMON /LSTM/M1(N),M2(N),M3(N) /LSTP/L1(N),L2(N),L3(N) DTX24=24.*DX*DT DO 20 K=1,N DO 20 J=1,N DO 20 I=1,N Q(I,J,K)=( & (U(M2(I),J,K)-27.*U(M1(I),J,K)+27.*U(I,J,K)-U(L1(I),J,K)) & +(V(I,M2(J),K)-27.*V(I,M1(J),K)+27.*V(I,J,K)-V(I,L1(J),K)) & +(W(I,J,M2(K))-27.*W(I,J,M1(K))+27.*W(I,J,K)-W(I,J,L1(K))) & )/DTX24 20 CONTINUE RETURN END C ---------------------------------------------------------------------- C Solving Pressure Poisson Equation C ---------------------------------------------------------------------- SUBROUTINE SBSORP PARAMETER (N=32,NH=N/2) COMMON /ERRP/ERRP /ITRP/ITRP COMMON /W1/Q(N,N,N) /W2/S(N,N,N) /W3/QFFT(N,N,N) /P/P(N,N,N) COMMON /DX/DX /C/CM3,CM2,CM1,C00,CP1,CP2,CP3 /CWV/CWV(N) COMMON /LSTM/M1(N),M2(N),M3(N) /LSTP/L1(N),L2(N),L3(N) SUMQ=0. DO 80 K=1,N DO 80 J=1,N DO 80 I=1,N SUMQ=SUMQ+Q(I,J,K)**2 80 CONTINUE ISWFFT=1 IF(ISWFFT.EQ.1) THEN CALL SBRFFT DO 90 K=1,N/2 DO 90 J=1,N DO 90 I=1,N WAV=CWV(I)+CWV(J)+CWV(K) IF(ABS(WAV).LE.1.E-10) WAV=1. S(I,J,K )=QFFT(I,J,K )/WAV S(I,J,K+NH)=QFFT(I,J,K+NH)/WAV 90 CONTINUE CALL SBRIFT ELSE DO 10 K=1,N DO 10 J=1,N DO 10 I=1,N S(I,J,K)=P(I,J,K) 10 CONTINUE ENDIF C WRITE(6,1100) ITRP=0 100 ITRP=ITRP+1 SUMR=0. C03=3.*C00 DO 30 K=1,N DO 30 J=1,N DO 30 I=1,N RESI= & CM3*S(I,J,M3(K))+CM2*S(I,J,M2(K))+CM1*S(I,J,M1(K)) & +CM3*S(I,M3(J),K)+CM2*S(I,M2(J),K)+CM1*S(I,M1(J),K) & +CM3*S(M3(I),J,K)+CM2*S(M2(I),J,K)+CM1*S(M1(I),J,K) & +C03*S(I,J,K) & +CP1*S(L1(I),J,K)+CP2*S(L2(I),J,K)+CP3*S(L3(I),J,K) & +CP1*S(I,L1(J),K)+CP2*S(I,L2(J),K)+CP3*S(I,L3(J),K) & +CP1*S(I,J,L1(K))+CP2*S(I,J,L2(K))+CP3*S(I,J,L3(K)) & -Q(I,J,K) S(I,J,K)=S(I,J,K)-1.5*RESI/C03 SUMR=SUMR+RESI**2 30 CONTINUE ERRP=SQRT(SUMR/SUMQ) C WRITE(6,1000) ITRP,SUMQ,SUMR,ERRP IF(ITRP.LT.100.AND.ERRP.GT.1.E-10) GO TO 100 C 1000 FORMAT(I8,3E13.4) C 1100 FORMAT(/5X,3HITR,9X,4HSUMQ,9X,4HSUMR,9X,4HERRP) RETURN END C ---------------------------------------------------------------------- C Correction Step C ---------------------------------------------------------------------- SUBROUTINE SBCORR PARAMETER (N=32) COMMON /U/U(N,N,N) /V/V(N,N,N) /W/W(N,N,N) /P/P(N,N,N) COMMON /W2/S(N,N,N) /DT/DT /DX/DX COMMON /LSTM/M1(N),M2(N),M3(N) /LSTP/L1(N),L2(N),L3(N) DX24=24.*DX DO 10 K=1,N DO 10 J=1,N DO 10 I=1,N P(I,J,K)=S(I,J,K) U(I,J,K)=U(I,J,K)-DT & *(S(M1(I),J,K)-27.*S(I,J,K)+27.*S(L1(I),J,K)-S(L2(I),J,K))/DX24 V(I,J,K)=V(I,J,K)-DT & *(S(I,M1(J),K)-27.*S(I,J,K)+27.*S(I,L1(J),K)-S(I,L2(J),K))/DX24 W(I,J,K)=W(I,J,K)-DT & *(S(I,J,M1(K))-27.*S(I,J,K)+27.*S(I,J,L1(K))-S(I,J,L2(K)))/DX24 10 CONTINUE RETURN END C ---------------------------------------------------------------------- C Courant Number C ---------------------------------------------------------------------- SUBROUTINE SBCOUR(COUR) PARAMETER (N=32) COMMON /U/U(N,N,N) /V/V(N,N,N) /W/W(N,N,N) COMMON /DX/DX /DT/DT COMMON /LSTM/M1(N),M2(N),M3(N) /LSTP/L1(N),L2(N),L3(N) UMAX=0. DO 10 K=1,N DO 10 J=1,N DO 10 I=1,N US=(-U(M2(I),J,K)+9.*U(M1(I),J,K)+9.*U(I,J,K)-U(L1(I),J,K))/16. VS=(-V(I,M2(J),K)+9.*V(I,M1(J),K)+9.*V(I,J,K)-V(I,L1(J),K))/16. WS=(-W(I,J,M2(K))+9.*W(I,J,M1(K))+9.*W(I,J,K)-W(I,J,L1(K)))/16. UABS=ABS(US)+ABS(VS)+ABS(WS) UMAX=AMAX1(UMAX,UABS) 10 CONTINUE COUR=DT*UABS/DX RETURN END C ---------------------------------------------------------------------- C Kinetic Energy C ---------------------------------------------------------------------- SUBROUTINE SBENGY(EK) PARAMETER (N=32) COMMON /EK/EK11,EK22,EK33 COMMON /U/U(N,N,N) /V/V(N,N,N) /W/W(N,N,N) EK11=0. EK22=0. EK33=0. DO 10 K=1,N DO 10 J=1,N DO 10 I=1,N EK11=EK11+U(I,J,K)**2 EK22=EK22+V(I,J,K)**2 EK33=EK33+W(I,J,K)**2 10 CONTINUE AN=1./REAL(N**3) EK11=AN*EK11 EK22=AN*EK22 EK33=AN*EK33 EK=(EK11+EK22+EK33)/2. RETURN END C ---------------------------------------------------------------------- SUBROUTINE SBEKII(EKII,EKIII) PARAMETER (N=32) COMMON /EK/EK11,EK22,EK33 COMMON /W1/UP(N,N,N) /W2/VP(N,N,N) /W3/WP(N,N,N) COMMON /LSTM/M1(N),M2(N),M3(N) /LSTP/L1(N),L2(N),L3(N) EK12=0. EK23=0. EK31=0. c BK11=0. c BK22=0. c BK33=0. DO 10 K=1,N DO 10 J=1,N DO 10 I=1,N EK12=EK12+UP(I,J,K)*VP(I,J,K) EK23=EK23+VP(I,J,K)*WP(I,J,K) EK31=EK31+WP(I,J,K)*UP(I,J,K) c BK11=BK11+UP(I,J,K)**2 c BK22=BK22+VP(I,J,K)**2 c BK33=BK33+WP(I,J,K)**2 10 CONTINUE BN=1./REAL(N**3) EK12=BN*EK12 EK23=BN*EK23 EK31=BN*EK31 c BK11=BN*BK11 c BK22=BN*BK22 c BK33=BN*BK33 c WRITE(6,1000) BK11,BK22,BK33,(BK11+BK22+BK33)/2. c & ,EK11,EK22,EK33,(EK11+EK22+EK33)/2. c 1000 FORMAT(/1H ,' Interpolated: Ku=',F12.8 c &,' Kv=',F12.8,' Kw=',F12.8,' K =',F12.8 c & /1H ,' Pure(Consistent): Ku=',F12.8 c &,' Kv=',F12.8,' Kw=',F12.8,' K =',F12.8/) CALL SBTNSR(EK11,EK22,EK33,EK12,EK23,EK31,EKII,EKIII) RETURN END C ---------------------------------------------------------------------- C Dissipation Rate C ---------------------------------------------------------------------- SUBROUTINE SBDISS(DS) PARAMETER (N=32) COMMON /DS/DS11,DS22,DS33 /RE/RE /DX/DX COMMON /LSTM/M1(N),M2(N),M3(N) /LSTP/L1(N),L2(N),L3(N) COMMON /U/U(N,N,N) /V/V(N,N,N) /W/W(N,N,N) DS11=0. DS22=0. DS33=0. DO 10 K=1,N DO 10 J=1,N DO 10 I=1,N DUX=U(M2(I),J,K)-27.*U(M1(I),J,K)+27.*U(I,J,K)-U(L1(I),J,K) DUY=U(I,M1(J),K)-27.*U(I,J,K)+27.*U(I,L1(J),K)-U(I,L2(J),K) DUZ=U(I,J,M1(K))-27.*U(I,J,K)+27.*U(I,J,L1(K))-U(I,J,L2(K)) DS11=DS11+DUX**2+DUY**2+DUZ**2 DVX=V(M1(I),J,K)-27.*V(I,J,K)+27.*V(L1(I),J,K)-V(L2(I),J,K) DVY=V(I,M2(J),K)-27.*V(I,M1(J),K)+27.*V(I,J,K)-V(I,L1(J),K) DVZ=V(I,J,M1(K))-27.*V(I,J,K)+27.*V(I,J,L1(K))-V(I,J,L2(K)) DS22=DS22+DVX**2+DVY**2+DVZ**2 DWX=W(M1(I),J,K)-27.*W(I,J,K)+27.*W(L1(I),J,K)-W(L2(I),J,K) DWY=W(I,M1(J),K)-27.*W(I,J,K)+27.*W(I,L1(J),K)-W(I,L2(J),K) DWZ=W(I,J,M2(K))-27.*W(I,J,M1(K))+27.*W(I,J,K)-W(I,J,L1(K)) DS33=DS33+DWX**2+DWY**2+DWZ**2 10 CONTINUE AN=2./REAL(N**3)/RE/(24.*DX)**2 DS11=AN*DS11 DS22=AN*DS22 DS33=AN*DS33 DS=(DS11+DS22+DS33)/2. RETURN END C ---------------------------------------------------------------------- SUBROUTINE SBDSII(DSII,DSIII,PRDS) PARAMETER (N=32) COMMON /DS/DS11,DS22,DS33 /RE/RE /DX/DX COMMON /U/U(N,N,N) /W1/UP(N,N,N) COMMON /V/V(N,N,N) /W2/VP(N,N,N) COMMON /W/W(N,N,N) /W3/WP(N,N,N) COMMON /LSTM/M1(N),M2(N),M3(N) /LSTP/L1(N),L2(N),L3(N) DS12=0. DS23=0. DS31=0. c BS11=0. c BS22=0. c BS33=0. PRDS=0. DO 10 K=1,N DO 10 J=1,N DO 10 I=1,N c BUX=(UP(M2(I),J,K)-8.*UP(M1(I),J,K)+8.*UP(L1(I),J,K) c & -UP(L2(I),J,K))/12. c BVY=(VP(I,M2(J),K)-8.*VP(I,M1(J),K)+8.*VP(I,L1(J),K) c & -VP(I,L2(J),K))/12. c BWZ=(WP(I,J,M2(K))-8.*WP(I,J,M1(K))+8.*WP(I,J,L1(K)) c & -WP(I,J,L2(K)))/12. DUX=(U(M2(I),J,K)-27.*U(M1(I),J,K)+27.*U(I,J,K)-U(L1(I),J,K))/24. DUY=(UP(I,M2(J),K)-8.*UP(I,M1(J),K)+8.*UP(I,L1(J),K) & -UP(I,L2(J),K))/12. DUZ=(UP(I,J,M2(K))-8.*UP(I,J,M1(K))+8.*UP(I,J,L1(K)) & -UP(I,J,L2(K)))/12. DVX=(VP(M2(I),J,K)-8.*VP(M1(I),J,K)+8.*VP(L1(I),J,K) & -VP(L2(I),J,K))/12. DVY=(V(I,M2(J),K)-27.*V(I,M1(J),K)+27.*V(I,J,K)-V(I,L1(J),K))/24. DVZ=(VP(I,J,M2(K))-8.*VP(I,J,M1(K))+8.*VP(I,J,L1(K)) & -VP(I,J,L2(K)))/12. DWX=(WP(M2(I),J,K)-8.*WP(M1(I),J,K)+8.*WP(L1(I),J,K) & -WP(L2(I),J,K))/12. DWY=(WP(I,M2(J),K)-8.*WP(I,M1(J),K)+8.*WP(I,L1(J),K) & -WP(I,L2(J),K))/12. DWZ=(W(I,J,M2(K))-27.*W(I,J,M1(K))+27.*W(I,J,K)-W(I,J,L1(K)))/24. D12KK=DUX*DVX+DUY*DVY+DUZ*DVZ D23KK=DVX*DWX+DVY*DWY+DVZ*DWZ D31KK=DWX*DUX+DWY*DUY+DWZ*DUZ DS12=DS12+D12KK DS23=DS23+D23KK DS31=DS31+D31KK c BS11=BS11+BUX**2+DUY**2+DUZ**2 c BS22=BS22+DVX**2+BVY**2+DVZ**2 c BS33=BS33+DWX**2+DWY**2+BWZ**2 PRDS=PRDS & +DUX*(DUX**2+DUY**2+DUZ**2) & +DVY*(DVX**2+DVY**2+DVZ**2) & +DWZ*(DWX**2+DWY**2+DWZ**2) & +(DUY+DVX)*D12KK +(DVZ+DWY)*D23KK +(DWX+DUZ)*D31KK 10 CONTINUE c AN=2./REAL(N**3)/RE/DX**2 c BS11=AN*BS11 c BS22=AN*BS22 c BS33=AN*BS33 BN=2./REAL(N**3)/RE/DX**2 DS12=BN*DS12 DS23=BN*DS23 DS31=BN*DS31 CN=-2./REAL(N**3)/RE/DX**3 PRDS=CN*PRDS c WRITE(6,1000) BS11,BS22,BS33,(BS11+BS22+BS33)/2. c & ,DS11,DS22,DS33,(DS11+DS22+DS33)/2. c 1000 FORMAT(/1H ,' Interpolated: Du=',F12.8 c &,' Dv=',F12.8,' Dw=',F12.8,' D =',F12.8 c & /1H ,' Pure(Consistent): Du=',F12.8 c &,' Dv=',F12.8,' Dw=',F12.8,' D =',F12.8/) CALL SBTNSR(DS11,DS22,DS33,DS12,DS23,DS31,DSII,DSIII) RETURN END C ---------------------------------------------------------------------- SUBROUTINE SBDSDS(DSDS) PARAMETER (N=32) COMMON /RE/RE /DX/DX COMMON /U/U(N,N,N) /UF/DU(N,N,N) COMMON /V/V(N,N,N) /VF/DV(N,N,N) COMMON /W/W(N,N,N) /WF/DW(N,N,N) COMMON /LSTM/M1(N),M2(N),M3(N) /LSTP/L1(N),L2(N),L3(N) DSDS=0. DO 20 K=1,N DO 20 J=1,N DO 20 I=1,N DU(I,J,K)=U(M2(I),J,K)-27.*U(M1(I),J,K)+27.*U(I,J,K)-U(L1(I),J,K) DV(I,J,K)=V(I,M2(J),K)-27.*V(I,M1(J),K)+27.*V(I,J,K)-V(I,L1(J),K) DW(I,J,K)=W(I,J,M2(K))-27.*W(I,J,M1(K))+27.*W(I,J,K)-W(I,J,L1(K)) 20 CONTINUE DO 30 K=1,N DO 30 J=1,N DO 30 I=1,N DUXX=DU(M1(I),J,K)-27.*DU(I,J,K)+27.*DU(L1(I),J,K)-DU(L2(I),J,K) DUXY=DU(I,M1(J),K)-27.*DU(I,J,K)+27.*DU(I,L1(J),K)-DU(I,L2(J),K) DUXZ=DU(I,J,M1(K))-27.*DU(I,J,K)+27.*DU(I,J,L1(K))-DU(I,J,L2(K)) DVYX=DV(M1(I),J,K)-27.*DV(I,J,K)+27.*DV(L1(I),J,K)-DV(L2(I),J,K) DVYY=DV(I,M1(J),K)-27.*DV(I,J,K)+27.*DV(I,L1(J),K)-DV(I,L2(J),K) DVYZ=DV(I,J,M1(K))-27.*DV(I,J,K)+27.*DV(I,J,L1(K))-DV(I,J,L2(K)) DWZX=DW(M1(I),J,K)-27.*DW(I,J,K)+27.*DW(L1(I),J,K)-DW(L2(I),J,K) DWZY=DW(I,M1(J),K)-27.*DW(I,J,K)+27.*DW(I,L1(J),K)-DW(I,L2(J),K) DWZZ=DW(I,J,M1(K))-27.*DW(I,J,K)+27.*DW(I,J,L1(K))-DW(I,J,L2(K)) DSDS=DSDS+ DUXX**2 +DVYY**2 +DWZZ**2 & +2.*( DUXY**2+DUXZ**2 +DVYX**2+DVYZ**2 +DWZX**2+DWZY**2) 30 CONTINUE DO 40 K=1,N DO 40 J=1,N DO 40 I=1,N DU(I,J,K)=U(I,M1(J),K)-27.*U(I,J,K)+27.*U(I,L1(J),K)-U(I,L2(J),K) DV(I,J,K)=V(I,J,M1(K))-27.*V(I,J,K)+27.*V(I,J,L1(K))-V(I,J,L2(K)) DW(I,J,K)=W(M1(I),J,K)-27.*W(I,J,K)+27.*W(L1(I),J,K)-W(L2(I),J,K) 40 CONTINUE DO 50 K=1,N DO 50 J=1,N DO 50 I=1,N DUYY=DU(I,M2(J),K)-27.*DU(I,M1(J),K)+27.*DU(I,J,K)-DU(I,L1(J),K) DUYZ=DU(I,J,M2(K))-27.*DU(I,J,M1(K))+27.*DU(I,J,K)-DU(I,J,L1(K)) DVZX=DV(M2(I),J,K)-27.*DV(M1(I),J,K)+27.*DV(I,J,K)-DV(L1(I),J,K) DVZZ=DV(I,J,M2(K))-27.*DV(I,J,M1(K))+27.*DV(I,J,K)-DV(I,J,L1(K)) DWXX=DW(M2(I),J,K)-27.*DW(M1(I),J,K)+27.*DW(I,J,K)-DW(L1(I),J,K) DWXY=DW(I,M2(J),K)-27.*DW(I,M1(J),K)+27.*DW(I,J,K)-DW(I,L1(J),K) DSDS=DSDS+ DUYY**2+DVZZ**2+DWXX**2 & +2.*( DUYZ**2+DVZX**2+DWXY**2) 50 CONTINUE DO 60 K=1,N DO 60 J=1,N DO 60 I=1,N DU(I,J,K)=U(I,J,M1(K))-27.*U(I,J,K)+27.*U(I,J,L1(K))-U(I,J,L2(K)) DV(I,J,K)=V(M1(I),J,K)-27.*V(I,J,K)+27.*V(L1(I),J,K)-V(L2(I),J,K) DW(I,J,K)=W(I,M1(J),K)-27.*W(I,J,K)+27.*W(I,L1(J),K)-W(I,L2(J),K) 60 CONTINUE DO 70 K=1,N DO 70 J=1,N DO 70 I=1,N DUZZ=DU(I,J,M2(K))-27.*DU(I,J,M1(K))+27.*DU(I,J,K)-DU(I,J,L1(K)) DVXX=DV(M2(I),J,K)-27.*DV(M1(I),J,K)+27.*DV(I,J,K)-DV(L1(I),J,K) DWYY=DW(I,M2(J),K)-27.*DW(I,M1(J),K)+27.*DW(I,J,K)-DW(I,L1(J),K) DSDS=DSDS+ DUZZ**2+DVXX**2+DWYY**2 70 CONTINUE DSDS=-2.*DSDS/REAL(N**3)/RE**2/(24.*DX)**4 RETURN END C ---------------------------------------------------------------------- C Tnsor invariants C ---------------------------------------------------------------------- SUBROUTINE SBTNSR(T11,T22,T33,T12,T23,T31,TII,TIII) DIMENSION T(3,3) TI=T11+T22+T33 T(1,1)=T11/TI-1./3. T(2,2)=T22/TI-1./3. T(3,3)=T33/TI-1./3. T(1,2)=T12/TI T(2,3)=T23/TI T(3,1)=T31/TI T(2,1)=T(1,2) T(3,2)=T(2,3) T(1,3)=T(3,1) TII=0. DO 10 J=1,3 DO 10 I=1,3 10 TII=TII+T(I,J)*T(J,I) TII=TII/2. TIII=0. DO 20 K=1,3 DO 20 J=1,3 DO 20 I=1,3 20 TIII=TIII+T(I,J)*T(J,K)*T(K,I) TIII=ABS(TIII)/3. RETURN END C ---------------------------------------------------------------------- C Divergence check C ---------------------------------------------------------------------- SUBROUTINE SBDIVU(ERRD) PARAMETER (N=32) COMMON /U/U(N,N,N) /V/V(N,N,N) /W/W(N,N,N) COMMON /DX/DX COMMON /LSTM/M1(N),M2(N),M3(N) /LSTP/L1(N),L2(N),L3(N) DX24=24.*DX ERRD=0. DO 10 K=1,N DO 10 J=1,N DO 10 I=1,N DIV=( & U(M2(I),J,K)-27.*U(M1(I),J,K)+27.*U(I,J,K)-U(L1(I),J,K) & +V(I,M2(J),K)-27.*V(I,M1(J),K)+27.*V(I,J,K)-V(I,L1(J),K) & +W(I,J,M2(K))-27.*W(I,J,M1(K))+27.*W(I,J,K)-W(I,J,L1(K)) & )/DX24 ERRD=ERRD+DIV**2 10 CONTINUE ERRD=SQRT(ERRD/N**3) RETURN END C ---------------------------------------------------------------------- C Data Save C ---------------------------------------------------------------------- SUBROUTINE SBSAVE PARAMETER (N=32) COMMON /STEP/ISTEP,TSTEP COMMON /U/U(N,N,N) /V/V(N,N,N) /W/W(N,N,N) COMMON /P/P(N,N,N) /UF/PR(N,N,N) COMMON /LSTM/M1(N),M2(N),M3(N) /LSTP/L1(N),L2(N),L3(N) COMMON /DT/DT /RE/RE /C/CM3,CM2,CM1,C00,CP1,CP2,CP3 PS2=DT/2./RE C03=3.*C00 DO 10 K=1,N DO 10 J=1,N DO 10 I=1,N SLAP= & CM3*P(I,J,M3(K))+CM2*P(I,J,M2(K))+CM1*P(I,J,M1(K)) & +CM3*P(I,M3(J),K)+CM2*P(I,M2(J),K)+CM1*P(I,M1(J),K) & +CM3*P(M3(I),J,K)+CM2*P(M2(I),J,K)+CM1*P(M1(I),J,K) & +C03*P(I,J,K) & +CP1*P(L1(I),J,K)+CP2*P(L2(I),J,K)+CP3*P(L3(I),J,K) & +CP1*P(I,L1(J),K)+CP2*P(I,L2(J),K)+CP3*P(I,L3(J),K) & +CP1*P(I,J,L1(K))+CP2*P(I,J,L2(K))+CP3*P(I,J,L3(K)) PR(I,J,K)=P(I,J,K)+PS2*SLAP 10 CONTINUE WRITE(12) ISTEP,TSTEP WRITE(12) U WRITE(12) V WRITE(12) W WRITE(12) PR CLOSE(12) write(6,1200) ISTEP,TSTEP 1200 format(/'Data SAVED: STEP=',I7,10X,'T=',F10.3/) RETURN END C ---------------------------------------------------------------------- SUBROUTINE SBPLOT PARAMETER (N=32,NHP=N/2+1) COMMON /W1/UP(N,N,N) /UF/PR(N,N,N) COMMON /W2/VP(N,N,N) COMMON /W3/WP(N,N,N) /WF/EN(N,N,N) COMMON /DX/DX /LSTM/M1(N),M2(N),M3(N) /LSTP/L1(N),L2(N),L3(N) DX12=12.*DX DO 10 K=1,N DO 10 J=1,N DO 10 I=1,N DUY=(UP(I,M2(J),K)-8.*UP(I,M1(J),K)+8.*UP(I,L1(J),K) & -UP(I,L2(J),K))/DX12 DUZ=(UP(I,J,M2(K))-8.*UP(I,J,M1(K))+8.*UP(I,J,L1(K)) & -UP(I,J,L2(K)))/DX12 DVX=(VP(M2(I),J,K)-8.*VP(M1(I),J,K)+8.*VP(L1(I),J,K) & -VP(L2(I),J,K))/DX12 DVZ=(VP(I,J,M2(K))-8.*VP(I,J,M1(K))+8.*VP(I,J,L1(K)) & -VP(I,J,L2(K)))/DX12 DWX=(WP(M2(I),J,K)-8.*WP(M1(I),J,K)+8.*WP(L1(I),J,K) & -WP(L2(I),J,K))/DX12 DWY=(WP(I,M2(J),K)-8.*WP(I,M1(J),K)+8.*WP(I,L1(J),K) & -WP(I,L2(J),K))/DX12 OM1=DWY-DVZ OM2=DUZ-DWX OM3=DVX-DUY EN(I,J,K)=OM1**1+OM2**2+OM3**3 10 CONTINUE WRITE(14,1100) DO 100 I=1,N DO 100 J=1,N WRITE(14,1000) I-1,J-1,UP(I,J,NHP),PR(I,J,NHP),EN(I,J,NHP) 100 CONTINUE CLOSE(14) 1000 FORMAT(2I5,3F20.10) 1100 FORMAT(/'t=3, z=pi' &//10X,12X,8Hvelocity,12X,8Hpressure,11X,9Henstrophy &/4X,1HI,4X,1HJ &,10X,10HU(I,J,N/2),10X,10HP(I,J,N/2),10X,10HE(I,J,N/2)/) write(6,1200) 1200 format(/'Plot Data SAVED') RETURN END C ---------------------------------------------------------------------- SUBROUTINE SBSOUK PARAMETER (N=32,NHP=N/2+1) COMMON /W1/UP(N,N,N) /UF/PR(N,N,N) OPEN(16,file='upz16t3.d',status='old',form='formatted') SU1=0. SU2=0. SU3=0. SP1=0. SP2=0. SP3=0. DO 100 I=1,N DO 100 J=1,N READ(16,1000) UPR,PRR SU1=SU1+UP(I,J,NHP)**2 SU2=SU2+UPR**2 SU3=SU3+UPR*UP(I,J,NHP) SP1=SP1+PR(I,J,NHP)**2 SP2=SP2+PRR**2 SP3=SP3+PRR*PR(I,J,NHP) 100 CONTINUE CLOSE(16) SU1=SQRT(SU1)/REAL(N) SU2=SQRT(SU2)/REAL(N) SU3=SU3/REAL(N**2) SP1=SQRT(SP1)/REAL(N) SP2=SQRT(SP2)/REAL(N) SP3=SP3/REAL(N**2) WRITE(6,1200) 'R_uu(z=pi)=',SU3/(SU1*SU2),SU3,SU1,SU2 WRITE(6,1200) 'R_pp(z=pi)=',SP3/(SP1*SP2),SP3,SP1,SP2 OPEN(18,file='upy16t3.d',status='old',form='formatted') SU1=0. SU2=0. SU3=0. SP1=0. SP2=0. SP3=0. DO 110 I=1,N DO 110 K=1,N READ(18,1000) UPR,PRR SU1=SU1+UP(I,NHP,K)**2 SU2=SU2+UPR**2 SU3=SU3+UPR*UP(I,NHP,K) SP1=SP1+PR(I,NHP,K)**2 SP2=SP2+PRR**2 SP3=SP3+PRR*PR(I,NHP,K) 110 CONTINUE CLOSE(18) SU1=SQRT(SU1)/REAL(N) SU2=SQRT(SU2)/REAL(N) SU3=SU3/REAL(N**2) SP1=SQRT(SP1)/REAL(N) SP2=SQRT(SP2)/REAL(N) SP3=SP3/REAL(N**2) WRITE(6,1200) 'R_uu(y=pi)=',SU3/(SU1*SU2),SU3,SU1,SU2 WRITE(6,1200) 'R_pp(y=pi)=',SP3/(SP1*SP2),SP3,SP1,SP2 1000 FORMAT(2E34.18) 1200 FORMAT(/4X,A11,F12.6,' =(',F12.6,2('/',F12.6),')') RETURN END C ---------------------------------------------------------------------- C Velocity Interpolation ( U ==> P ) C ---------------------------------------------------------------------- SUBROUTINE SBVINT PARAMETER (N=32) COMMON /U/U(N,N,N) /W1/UP(N,N,N) COMMON /V/V(N,N,N) /W2/VP(N,N,N) COMMON /W/W(N,N,N) /W3/WP(N,N,N) COMMON /LSTM/M1(N),M2(N),M3(N) /LSTP/L1(N),L2(N),L3(N) DO 10 K=1,N DO 10 J=1,N DO 10 I=1,N UP(I,J,K)= & (-U(M2(I),J,K)+9.*U(M1(I),J,K)+9.*U(I,J,K)-U(L1(I),J,K))/16. VP(I,J,K)= & (-V(I,M2(J),K)+9.*V(I,M1(J),K)+9.*V(I,J,K)-V(I,L1(J),K))/16. WP(I,J,K)= & (-W(I,J,M2(K))+9.*W(I,J,M1(K))+9.*W(I,J,K)-W(I,J,L1(K)))/16. 10 CONTINUE RETURN END C ---------------------------------------------------------------------- C Energy Spectrum C ---------------------------------------------------------------------- SUBROUTINE SBSPEC PARAMETER (N=32,NH=N/2) DIMENSION WSPEC(NH-1),PSPEC(NH-1),OSPEC(NH-1) COMMON /U/U(N,N,N) /V/V(N,N,N) /W/W(N,N,N) /P/P(N,N,N) COMMON /W1/A(N,N,N) /W3/F(N,N,N) COMMON /LSTM/M1(N),M2(N),M3(N) /LSTP/L1(N),L2(N),L3(N) COMMON /DT/DT /DX/DX /RE/RE COMMON /C/CM3,CM2,CM1,C00,CP1,CP2,CP3 DO 10 K=1,N DO 10 J=1,N DO 10 I=1,N A(I,J,K)=W(I,J,K) 10 CONTINUE CALL SBRFFT DO 20 K=2,NH SUM=0. DO 30 J=1,N DO 30 I=1,N SUM=SUM+F(I,J,K)**2+F(I,J,K+NH)**2 30 CONTINUE WSPEC(K-1)=SUM/REAL(N**6) 20 CONTINUE PS2=DT/2./RE C03=3.*C00 DO 12 K=1,N DO 12 J=1,N DO 12 I=1,N SLAP= & CM3*P(I,J,M3(K))+CM2*P(I,J,M2(K))+CM1*P(I,J,M1(K)) & +CM3*P(I,M3(J),K)+CM2*P(I,M2(J),K)+CM1*P(I,M1(J),K) & +CM3*P(M3(I),J,K)+CM2*P(M2(I),J,K)+CM1*P(M1(I),J,K) & +C03*P(I,J,K) & +CP1*P(L1(I),J,K)+CP2*P(L2(I),J,K)+CP3*P(L3(I),J,K) & +CP1*P(I,L1(J),K)+CP2*P(I,L2(J),K)+CP3*P(I,L3(J),K) & +CP1*P(I,J,L1(K))+CP2*P(I,J,L2(K))+CP3*P(I,J,L3(K)) A(I,J,K)=P(I,J,K)-PS2*SLAP 12 CONTINUE CALL SBRFFT DO 22 K=2,NH SUM=0. DO 32 J=1,N DO 32 I=1,N SUM=SUM+F(I,J,K)**2+F(I,J,K+NH)**2 32 CONTINUE PSPEC(K-1)=SUM/REAL(N**6) 22 CONTINUE DX24=24.*DX DO 14 K=1,N DO 14 J=1,N DO 14 I=1,N A(I,J,K)= & (V(M1(I),J,K)-27.*V(I,J,K)+27.*V(L1(I),J,K)-V(L2(I),J,K))/DX24 &-(U(I,M1(J),K)-27.*U(I,J,K)+27.*U(I,L1(J),K)-U(I,L2(J),K))/DX24 14 CONTINUE CALL SBRFFT DO 24 K=2,NH SUM=0. DO 34 J=1,N DO 34 I=1,N SUM=SUM+F(I,J,K)**2+F(I,J,K+NH)**2 34 CONTINUE OSPEC(K-1)=SUM/REAL(N**6) 24 CONTINUE WRITE(6,1100) DO 40 K=1,NH-1 40 WRITE(6,1000) K,WSPEC(K),PSPEC(K),OSPEC(K) 1000 FORMAT(I5,3E15.6) 1100 FORMAT(/2X,3HK_3,9X,6HW_SPEC,9X,6HP_SPEC,9X,6HO_SPEC/) RETURN END C ---------------------------------------------------------------------- C FFT & Inverse FFT C ---------------------------------------------------------------------- SUBROUTINE SBRFFT PARAMETER (N=32,NH=N/2) DIMENSION WK(N) COMPLEX A(N),B(N) COMMON /W1/Q(N,N,N) /W3/F(N,N,N) C ----- z direction ----- DO 14 J=1,N DO 14 I=1,N DO 24 K=1,N A(K)=CMPLX( Q(I,J,K), 0. ) 24 CONTINUE CALL CFFT( A, N, 1, B, WK ) DO 34 K=1,NH F(I,J,K )= REAL(B(K)) F(I,J,K+NH)=AIMAG(B(K)) 34 CONTINUE 14 CONTINUE C ----- y direction ----- DO 12 K=1,NH DO 12 I=1,N DO 22 J=1,N A(J)=CMPLX( F(I,J,K), F(I,J,K+NH) ) 22 CONTINUE CALL CFFT( A, N, 1, B, WK ) DO 32 J=1,N F(I,J,K )= REAL(B(J)) F(I,J,K+NH)=AIMAG(B(J)) 32 CONTINUE 12 CONTINUE C ----- X direction ----- DO 10 K=1,NH DO 10 J=1,N DO 20 I=1,N A(I)=CMPLX( F(I,J,K), F(I,J,K+NH) ) 20 CONTINUE CALL CFFT( A, N, 1, B, WK ) DO 30 I=1,N F(I,J,K )= REAL(B(I)) F(I,J,K+NH)=AIMAG(B(I)) 30 CONTINUE 10 CONTINUE RETURN END C ---------------------------------------------------------------------- SUBROUTINE SBRIFT PARAMETER (N=32,NH=N/2) DIMENSION WK(N) COMPLEX A(N),B(N) COMMON /W2/S(N,N,N) C ----- X direction ----- DO 10 K=1,NH DO 10 J=1,N DO 20 I=1,N A(I)=CMPLX( S(I,J,K), S(I,J,K+NH) ) 20 CONTINUE CALL CFFT( A, N, 2, B, WK ) DO 30 I=1,N S(I,J,K )= REAL(B(I)) S(I,J,K+NH)=AIMAG(B(I)) 30 CONTINUE 10 CONTINUE C ----- y direction ----- DO 12 K=1,NH DO 12 I=1,N DO 22 J=1,N A(J)=CMPLX( S(I,J,K), S(I,J,K+NH) ) 22 CONTINUE CALL CFFT( A, N, 2, B, WK ) DO 32 J=1,N S(I,J,K )= REAL(B(J)) S(I,J,K+NH)=AIMAG(B(J)) 32 CONTINUE 12 CONTINUE C ----- z direction ----- DO 14 J=1,N DO 14 I=1,N DO 24 K=1,NH A(K)=CMPLX( S(I,J,K), S(I,J,K+NH) ) 24 CONTINUE A(NH+1)=CMPLX(0.,0.) DO 26 K=NH+2,N A(K)=CONJG(A(N+2-K)) 26 CONTINUE CALL CFFT( A, N, 2, B, WK ) DO 34 K=1,N S(I,J,K)= REAL(B(K)) 34 CONTINUE 14 CONTINUE RETURN END