C *************************************************************** C * DNS37AB4 Ver.1 1994.09.13 (4th order accuracy) * C * Ver.3C 1997.07.12 (corrected) * C * by Takeo KAJISHIMA, Dept. Mech. Eng., Osaka University * C * 2-1 Yamada-oka, Suita, Osaka 565-0871, Japan * C * TEL: +81 6 879 7249, FAX: + 81 6 879 7250 * C * E-mail: kajisima@mech.eng.osaka-u.ac.jp * C *************************************************************** PROGRAM DNS37AB4 PARAMETER (NX=64,NY=64,NZ=64) IMPLICIT DOUBLE PRECISION (A-H,O-Z) CHARACTER*3 CT(0:100) CHARACTER*4 FILE1 COMMON /FILE1/FILE1 /RMS/URMSX,URMSC,VRMSC,WRMSC COMMON /RET/RET /DT/DT /STEP/ISTEP /TSTEP/TSTEP /START/ISTART COMMON /DIV/DIVX /COU/COMX /ENE/ENE /UME/UME,UMX COMMON /ERRP/POIERR /ITRP/ITRP DATA CT/'000', &'001','002','003','004','005','006','007','008','009','010', &'011','012','013','014','015','016','017','018','019','020', &'021','022','023','024','025','026','027','028','029','030', &'031','032','033','034','035','036','037','038','039','040', &'041','042','043','044','045','046','047','048','049','050', &'051','052','053','054','055','056','057','058','059','060', &'061','062','063','064','065','066','067','068','069','070', &'071','072','073','074','075','076','077','078','079','080', &'081','082','083','084','085','086','087','088','089','090', &'091','092','093','094','095','096','097','098','099','100'/ FILE1='dns4' RET=300.D0 DT= 1.0D-4 ISKIP=40 ISAVE=100 ISTOP=ISKIP*ISAVE ICOUNT=0 ICT=0 CALL SBRMSH CALL SBRDTR(CT(ICT)) WRITE( 6,1000) CALL SBRUMR CALL SBRST1 CALL SBRCHK WRITE( 6,1100) ISTART,TSTEP,ITRP,POIERR,DIVX,COMX & ,UME,UMX,INT(RET*UME),INT(RET*UMX),ENE,URMSX,URMSC,VRMSC,WRMSC 1000 FORMAT(3X,4HSTEP,7X,1HT,1X,4HITRP,4X,6HPOIerr,4X,6HDIVmax &,2X,6HCOUmax,3X,5HUmean,4X,4HUmax,4X,3HRem,4X,3HRex &,4X,6HEnergy,3X,5HU'max,5X,3HU'c,5X,3HV'c,5X,3HW'c) 1100 FORMAT(I7,F8.4,I5,2E10.2,3F8.3,2I7,F10.5,4F8.4) C================================================================= 10 CONTINUE CALL SBRCON CALL SBRVIS IF(ICOUNT.GE.1) THEN CALL SBRPRE(1.5D0,-0.5D0) ELSE CALL SBRPRE(1.0D0, 0.0D0) ENDIF CALL SBRRHP IPSTOP=20 CALL SBRSOR(IPSTOP) CALL SBRCOR ICOUNT=ICOUNT+1 ISTEP=ISTART+ICOUNT TSTEP=TSTEP+DT IF(MOD(ISTEP,10).eq.0.OR.ICOUNT.LE.10) THEN CALL SBRUMR CALL SBRST1 CALL SBRCHK WRITE( 6,1100) ISTEP,TSTEP,ITRP,POIERR,DIVX,COMX & ,UME,UMX,INT(RET*UME),INT(RET*UMX),ENE,URMSX,URMSC,VRMSC,WRMSC ENDIF IF(MOD(ISTEP,ISKIP).EQ.0) THEN ICT=ICT+1 CALL SBRDTS(CT(ICT)) ENDIF IF(ICOUNT.GE.ISTOP) GO TO 20 IF(ISTEP.GE.10.AND.DIVX.GE.1.D+2) GO TO 30 GO TO 10 C================================================================= 20 CONTINUE STOP 30 CONTINUE WRITE(6,*) 'Stopped because of the numerical instability' STOP END C ********************************************************************* C * SBR. MSH : MESH PARAMETERS * C ********************************************************************* SUBROUTINE SBRMSH PARAMETER (NX=64,NY=64,NZ=64,NY1=NY+1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /H/HX,HZ /Y/YV(0:NY),YP(0:NY1) COMMON /D/DX,DZ /DY/DY(NY) /RET/RET /DT/DT CALL MESHY CALL MESHXZ WRITE(6,1000) NX,HX,DX,RET*HX,RET*DX,NZ,HZ,DZ,RET*HZ,RET*DZ &,NY,1.,DY(1),DY(NY/2),RET*DY(1),RET*DY(NY/2),YP(1),RET*YP(1) &,INT(RET),DT 1000 FORMAT(/10(1H*)," IN SBR.MSH ",10(1H*)/ &/5X,"NX=",I4,5X,"HX=",F7.3,3X,"DX=",F7.5 &,5X,"HX+=",F7.1,3X,"DX+=",F7.3 &/5X,"NZ=",I4,5X,"HZ=",F7.3,3X,"DZ=",F7.5 &,5X,"HZ+=",F7.1,3X,"DZ+=",F7.3 &/5X,"NY=",I4,5X,"HY=",F7.3,3X,"DYmin=",F7.5,3X,"DYmax=",F7.5 &/30X,"DYmin+=",F6.3,3X,"DYmax+=",F6.2 &/5X,"YP(1)=",F7.5,3X,"YP(1)+=",F6.2 &/5X,"RET=",I6/5X,"DT=",F9.5) CALL MDQLST CALL MDQ1VN CALL MDQ1VD CALL MDQ1P CALL MDQ1PW CALL MDQ2V CALL MDQ2PV CALL MDQ2PP RETURN END C--------------------------------------------------------------------- C----- Mesh generation ..... for X, Z directions ------------------ SUBROUTINE MESHXZ PARAMETER (NX=64,NY=64,NZ=64) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /H/HX,HZ /D/DX,DZ /RET/RET COMMON /IP/IP1(NX),IP2(NX),IP3(NX) /IM/IM1(NX),IM2(NX),IM3(NX) COMMON /KP/KP1(NZ),KP2(NZ),KP3(NZ) /KM/KM1(NZ),KM2(NZ),KM3(NZ) COMMON /TXZA/TXA,TZA /BXZA/BXA,BZA /USI/E0 COMMON /CMPX/CM3X,CM2X,CM1X,C00X,CP1X,CP2X,CP3X COMMON /CMPZ/CM3Z,CM2Z,CM1Z,C00Z,CP1Z,CP2Z,CP3Z COMMON /AMX/AM2X,AM1X,A00X,AP1X /AMZ/AM2Z,AM1Z,A00Z,AP1Z COMMON /BMX/BM1X,B00X,BP1X,BP2X /BMZ/BM1Z,B00Z,BP1Z,BP2Z DX=18.D0/RET DZ= 9.D0/RET HX=NX*DX HZ=NZ*DZ BXA=1.D0/(24.D0*DX) BZA=1.D0/(24.D0*DZ) E0 =1.D0/ 16.D0 TXA=E0*BXA TZA=E0*BZA C --- List vectors to identify the neighbouring point DO 50 I=1,NX IP1(I)=I+1 IP2(I)=I+2 IP3(I)=I+3 IM1(I)=I-1 IM2(I)=I-2 IM3(I)=I-3 IF(IP1(I).GT.NX) IP1(I)=IP1(I)-NX IF(IP2(I).GT.NX) IP2(I)=IP2(I)-NX IF(IP3(I).GT.NX) IP3(I)=IP3(I)-NX IF(IM1(I).LT. 1) IM1(I)=IM1(I)+NX IF(IM2(I).LT. 1) IM2(I)=IM2(I)+NX IF(IM3(I).LT. 1) IM3(I)=IM3(I)+NX 50 CONTINUE DO 60 K=1,NZ KP1(K)=K+1 KP2(K)=K+2 KP3(K)=K+3 KM1(K)=K-1 KM2(K)=K-2 KM3(K)=K-3 IF(KP1(K).GT.NZ) KP1(K)=KP1(K)-NZ IF(KP2(K).GT.NZ) KP2(K)=KP2(K)-NZ IF(KP3(K).GT.NZ) KP3(K)=KP3(K)-NZ IF(KM1(K).LT. 1) KM1(K)=KM1(K)+NZ IF(KM2(K).LT. 1) KM2(K)=KM2(K)+NZ IF(KM3(K).LT. 1) KM3(K)=KM3(K)+NZ 60 CONTINUE C --- Parameters for finite-difference of Poisson equation DX24=24.D0*DX AM2X= 1.D0/DX24 AM1X=-27.D0/DX24 A00X= 27.D0/DX24 AP1X= -1.D0/DX24 BM1X= 1.D0/DX24 B00X=-27.D0/DX24 BP1X= 27.D0/DX24 BP2X= -1.D0/DX24 CM3X=AM2X*BM1X CM2X=AM2X*B00X+AM1X*BM1X CM1X=AM2X*BP1X+AM1X*B00X+A00X*BM1X C00X=AM2X*BP2X+AM1X*BP1X+A00X*B00X+AP1X*BM1X CP1X= AM1X*BP2X+A00X*BP1X+AP1X*B00X CP2X= A00X*BP2X+AP1X*BP1X CP3X= AP1X*BP2X DZ24=24.D0*DZ AM2Z= 1.D0/DZ24 AM1Z=-27.D0/DZ24 A00Z= 27.D0/DZ24 AP1Z= -1.D0/DZ24 BM1Z= 1.D0/DZ24 B00Z=-27.D0/DZ24 BP1Z= 27.D0/DZ24 BP2Z= -1.D0/DZ24 CM3Z=AM2Z*BM1Z CM2Z=AM2Z*B00Z+AM1Z*BM1Z CM1Z=AM2Z*BP1Z+AM1Z*B00Z+A00Z*BM1Z C00Z=AM2Z*BP2Z+AM1Z*BP1Z+A00Z*B00Z+AP1Z*BM1Z CP1Z= AM1Z*BP2Z+A00Z*BP1Z+AP1Z*B00Z CP2Z= A00Z*BP2Z+AP1Z*BP1Z CP3Z= AP1Z*BP2Z RETURN END C--------------------------------------------------------------------- C----- Mesh generation ..... for Y direction ---------------------- SUBROUTINE MESHY PARAMETER (NX=64,NY=64,NZ=64,NY1=NY+1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /Y/YV(0:NY),YP(0:NY1) /YR/YR(0:NY1) COMMON /DY/DY(NY) /DYC/DYC(NY-1) ALG=0.95D0 AT=DLOG((1.D0+ALG)/(1.D0-ALG))/2.D0 YV(0)=0. DO 10 J=1,NY-1 ETA=AT*(-1.D0+2.D0*DBLE(J)/DBLE(NY)) YV(J)=(DTANH(ETA)/ALG+1.D0)/2.D0 10 CONTINUE YV(NY)=1.D0 DO 15 J=1,NY ETA=AT*(-1.D0+2.D0*(DBLE(J)-0.5D0)/DBLE(NY)) YP(J)=(DTANH(ETA)/ALG+1.D0)/2.D0 15 CONTINUE C ... Outer points (half mesh) YP(0)=2.D0*YV(0)-YP(1) YP(NY1)=2.D0*YV(NY)-YP(NY) YR(0)=YV(0) DO 20 J=1,NY YR(J)=YP(J) DY(J)=-YV(J-1)+YV(J) 20 CONTINUE YR(NY1)=YV(NY) DO 30 J=1,NY-1 DYC(J)=-YP(J)+YP(J+1) 30 CONTINUE RETURN END C----------------------------------------------------------------------- SUBROUTINE MDQLST PARAMETER (NY=64,NY1=NY+1,NYH1=NY-1,MM=4,MH=MM/2,MI=MH-1) COMMON /NSEV/NSV(0:NY),NEV(0:NY) /MSEV/MSV(0:NY),MEV(0:NY) COMMON /NSEP/NSP(NY),NEP(NY) /MSEP/MSP(NY),MEP(NY) DO 10 J=0,NY NSV(J)=-MI NEV(J)= MH JS=J+NSV(J) JE=J+NEV(J) IF(JS.LT. 0) NSV(J)=NSV(J)+(0-JS) IF(JE.GT.NY1) NEV(J)=NEV(J)-(JE-NY1) MSV(J)=-MI MEV(J)=+MH JS=J+MSV(J) JE=J+MEV(J) IF(JS.LT. 1) MSV(J)=MSV(J)+(1-JS) IF(JE.GT.NY) MEV(J)=MEV(J)-(JE-NY) 10 CONTINUE DO 30 J=1,NY NSP(J)=-MH NEP(J)=+MI JS=J+NSP(J) JE=J+NEP(J) IF(JS.LT. 0) NSP(J)=NSP(J)+(0-JS) IF(JE.GT.NY) NEP(J)=NEP(J)-(JE-NY) MSP(J)=-MH MEP(J)=+MI JS=J+MSP(J) JE=J+MEP(J) IF(JS.LT. 1) MSP(J)=MSP(J)+(1-JS) IF(JE.GT.NYH1) MEP(J)=MEP(J)-(JE-NYH1) 30 CONTINUE RETURN END C----------------------------------------------------------------------- C---- FDM operators at V points using data on P points and INSIDE walls SUBROUTINE MDQ1VN PARAMETER (NY=64,NY1=NY+1,MM=4,MH=MM/2,MI=MH-1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION E0V(-MI:MH,0:NY),E1V(-MI:MH,0:NY) COMMON /Y/YV(0:NY),YP(0:NY1) /YK/YK(-MI:MH) COMMON /NSEV/NSV(0:NY),NEV(0:NY) /MSEV/MSV(0:NY),MEV(0:NY) COMMON /D0VN/D0VN(-MI:MH,0:NY) /G0W/G0WA(0:MH),G0WB(-MH:0) COMMON /D1VN/D1VN(-MI:MH,0:NY) /G1W/G1WA(0:MH),G1WB(-MH:0) DO 10 J=0,NY NS=NSV(J) NE=NEV(J) DO 20 M=-MI,MH YK(M)=0.D0 E0V(M,J)=0.D0 E1V(M,J)=0.D0 20 CONTINUE YC=YV(J) DO 40 M=NS,NE 40 YK(M)=YP(J+M) PI0=1.D0 DO 50 M=NS,NE 50 PI0=PI0*(YC-YK(M)) PI1=0.D0 DO 60 L1=NS,NE PIS=1.D0 DO 70 M=NS,NE 70 IF(M.NE.L1) PIS=PIS*(YC-YK(M)) 60 PI1=PI1+PIS SUM0=0.D0 SUM1=0.D0 NPT=NE IF(J.EQ.0) NPT=NS DO 90 K=NS,NE IF(K.NE.NPT) THEN YCK=YK(K) PIK=1.D0 DO 30 M=NS,NE 30 IF(M.NE.K) PIK=PIK*(YCK-YK(M)) XD=YC-YCK E0V(K,J)=PI0/XD/PIK E1V(K,J)=(XD*PI1-PI0)/XD**2/PIK SUM0=SUM0+E0V(K,J) SUM1=SUM1+E1V(K,J) ENDIF 90 CONTINUE E0V(NPT,J)=1.D0-SUM0 E1V(NPT,J)=-SUM1 10 CONTINUE C ... Neumann Boundary Condition at the wall (dP/dn=0) for Pressure J0E=NEV(0) J1S=NY+NSV(NY) DO 100 J=0,NY DO 100 M=-MI,MH D0VN(M,J)=E0V(M,J) D1VN(M,J)=E1V(M,J) 100 CONTINUE DO 150 M=0,MH G0WA(M)=0.D0 G1WA(M)=0.D0 G0WB(-M)=0.D0 G1WB(-M)=0.D0 150 CONTINUE DO 120 J=0,NY IF(J+NSV(J).EQ.0) THEN DO 130 M=NSV(J)+1,NEV(J) JM=J+M IF(JM.LE.J0E) THEN D0VN(M,J)=D0VN(M,J)-D0VN(NSV(J),J)*E1V(JM,0)/E1V(NSV(0),0) D1VN(M,J)=D1VN(M,J)-D1VN(NSV(J),J)*E1V(JM,0)/E1V(NSV(0),0) ENDIF 130 CONTINUE D0VN(NSV(J),J)=0.D0 D1VN(NSV(J),J)=0.D0 G0WA(J)=E0V(-J,J)/E1V(0,0) G1WA(J)=E1V(-J,J)/E1V(0,0) ENDIF IF(J+NEV(J).EQ.NY1) THEN DO 140 M=NSV(J),NEV(J)-1 JM=J+M IF(JM.GE.J1S) THEN D0VN(M,J)=D0VN(M,J)-D0VN(NEV(J),J)*E1V(JM-NY,NY)/E1V(NEV(NY),NY) D1VN(M,J)=D1VN(M,J)-D1VN(NEV(J),J)*E1V(JM-NY,NY)/E1V(NEV(NY),NY) ENDIF 140 CONTINUE D0VN(NEV(J),J)=0.D0 D1VN(NEV(J),J)=0.D0 G0WB(J-NY)=E0V(NY1-J,J)/E1V(1,NY) G1WB(J-NY)=E1V(NY1-J,J)/E1V(1,NY) ENDIF 120 CONTINUE RETURN END C----------------------------------------------------------------------- C---- FDM operators at V points using data on P points and AT walls ---- SUBROUTINE MDQ1VD PARAMETER (NY=64,NY1=NY+1,MM=4,MH=MM/2,MI=MH-1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /NSEV/NSV(0:NY),NEV(0:NY) /MSEV/MSV(0:NY),MEV(0:NY) COMMON /Y/YV(0:NY),YP(0:NY1) /YR/YR(0:NY1) /YK/YK(-MI:MH) COMMON /D0V/D0V(-MI:MH,0:NY) /D1V/D1V(-MI:MH,0:NY) DO 10 J=0,NY NS=NSV(J) NE=NEV(J) DO 20 M=-MI,MH YK(M)=0.D0 D0V(M,J)=0.D0 D1V(M,J)=0.D0 20 CONTINUE YC=YV(J) DO 40 M=NS,NE 40 YK(M)=YR(J+M) PI0=1.D0 DO 50 M=NS,NE 50 PI0=PI0*(YC-YK(M)) PI1=0.D0 DO 60 L1=NS,NE PIS=1.D0 DO 70 M=NS,NE 70 IF(M.NE.L1) PIS=PIS*(YC-YK(M)) 60 PI1=PI1+PIS SUM0=0.D0 SUM1=0.D0 NPT=NE IF(J.EQ.0) NPT=NS DO 90 K=NS,NE IF(K.NE.NPT) THEN YCK=YK(K) PIK=1.D0 DO 30 M=NS,NE 30 IF(M.NE.K) PIK=PIK*(YCK-YK(M)) XD=YC-YCK D0V(K,J)=PI0/XD/PIK D1V(K,J)=(XD*PI1-PI0)/XD**2/PIK SUM0=SUM0+D0V(K,J) SUM1=SUM1+D1V(K,J) ENDIF 90 CONTINUE D0V(NPT,J)=1.D0-SUM0 D1V(NPT,J)= -SUM1 10 CONTINUE C For the Dirichlet Boundary Condition at the wall (U,W=0) for Velocity, C use MSV & MEV instead of NSV & NEV RETURN END C----------------------------------------------------------------------- C--- FDM operators for velocity at P points using data on V points ---- SUBROUTINE MDQ1P PARAMETER (NY=64,NY1=NY+1,NYH1=NY-1,MM=4,MH=MM/2,MI=MH-1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /NSEP/NSP(NY),NEP(NY) /MSEP/MSP(NY),MEP(NY) COMMON /Y/YV(0:NY),YP(0:NY1) /YK/YK(-MH:MI) COMMON /D0P/D0P(-MH:MI,NY) /D1P/D1P(-MH:MI,NY) DO 30 J=1,NY MS=NSP(J) ME=NEP(J) DO 20 M=-MH,MI YK(M)=0.D0 D0P(M,J)=0.D0 D1P(M,J)=0.D0 20 CONTINUE YC=YP(J) DO 40 M=MS,ME 40 YK(M)=YV(J+M) PI0=1.D0 DO 50 M=MS,ME 50 PI0=PI0*(YC-YK(M)) PI1=0.D0 DO 60 L1=MS,ME PIS=1.D0 DO 65 M=MS,ME 65 IF(M.NE.L1) PIS=PIS*(YC-YK(M)) 60 PI1=PI1+PIS SUM0=0.D0 SUM1=0.D0 MPT=ME IF(J.EQ.0) MPT=MS DO 90 K=MS,ME IF(K.NE.MPT) THEN YCK=YK(K) PIK=1.D0 DO 85 M=MS,ME 85 IF(M.NE.K) PIK=PIK*(YCK-YK(M)) XD=YC-YCK D0P(K,J)=PI0/XD/PIK D1P(K,J)=(XD*PI1-PI0)/XD**2/PIK SUM0=SUM0+D0P(K,J) SUM1=SUM1+D1P(K,J) ENDIF 90 CONTINUE D0P(MPT,J)=1.D0-SUM0 D1P(MPT,J)= -SUM1 30 CONTINUE C For the Dirichlet Boundary Condition at the wall (V=0) for Velocity, C use MSP & MEP instead of NSP & NEP RETURN END C----------------------------------------------------------------------- C--- FDM operators for dV/dY at the wall ------------------------------ SUBROUTINE MDQ1PW PARAMETER (NY=64,NY1=NY+1,MM=4,MH=MM/2) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /Y/YV(0:NY),YP(0:NY1) COMMON /D1W/D1W0(0:MH),D1W1(-MH:0) DIMENSION YK0(0:MH),YK1(-MH:0) DO 20 M=0,MH YK0(M)=YV(M) D1W0(M)=0.D0 YK1(-M)=YV(NY-M) D1W1(-M)=0.D0 20 CONTINUE YC0=YV(0) YC1=YV(NY) PI0=0.D0 PI1=0.D0 DO 60 L1=0,MH PIS0=1.D0 PIS1=1.D0 DO 65 M=0,MH IF(M.NE.L1) THEN PIS0=PIS0*(YC0-YK0(M)) PIS1=PIS1*(YC1-YK1(-M)) ENDIF 65 CONTINUE PI0=PI0+PIS0 PI1=PI1+PIS1 60 CONTINUE SUM0=0.D0 SUM1=0.D0 DO 90 K=1,MH YCK0=YK0( K) YCK1=YK1(-K) PIK0=1.D0 PIK1=1.D0 DO 85 M=0,MH IF(M.NE.K) THEN PIK0=PIK0*(YCK0-YK0( M)) PIK1=PIK1*(YCK1-YK1(-M)) ENDIF 85 CONTINUE XD0=YC0-YCK0 XD1=YC1-YCK1 D1W0( K)=PI0/XD0/PIK0 D1W1(-K)=PI1/XD1/PIK1 SUM0=SUM0+D1W0( K) SUM1=SUM1+D1W1(-K) 90 CONTINUE D1W0(0)=-SUM0 D1W1(0)=-SUM1 C ... Dirichlet Boundary Condition at the wall (V=0) for Velocity D1W0(0)=0.D0 D1W1(0)=0.D0 RETURN END C----------------------------------------------------------------------- C--- FDM operators for viscous terms at V points --------------------- SUBROUTINE MDQ2V PARAMETER (NY=64,NY1=NY+1,MM=4,MH=MM/2,MI=MH-1,MP=MH+MI) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /D2V/D2V(-MP:MP,0:NY) /LSEV/LSV(0:NY),LEV(0:NY) COMMON /D1V/D1V(-MI:MH,0:NY) /NSEV/NSV(0:NY),NEV(0:NY) COMMON /D1P/D1P(-MH:MI,NY) /MSEP/MSP(NY),MEP(NY) COMMON /D1W/D1W0(0:MH),D1W1(-MH:0) DO 10 J=0,NY LSV(J)=999 LEV(J)=-999 DO 20 M=-MP,MP D2V(M,J)=0.D0 20 CONTINUE DO 40 M1=NSV(J),NEV(J) JM=J+M1 C ... dV/dY is given at inner P point (inner) IF(JM.GE.1.AND.JM.LE.NY) THEN DO 30 M2=MSP(JM),MEP(JM) M=M1+M2 IF(M.LT.LSV(J)) LSV(J)=M IF(M.GT.LEV(J)) LEV(J)=M D2V(M,J)=D2V(M,J)+D1P(M2,JM)*D1V(M1,J) 30 CONTINUE ELSE C ... dV/dY is given at the wall IF(JM.EQ.0) THEN DO 32 M2=1,MH M=M1+M2 IF(M.LT.LSV(J)) LSV(J)=M D2V(M,J)=D2V(M,J)+D1W0(M2)*D1V(M1,J) 32 CONTINUE ENDIF IF(JM.EQ.NY1) THEN DO 34 M2=-MH,-1 M=M1+M2 -1 IF(M.GT.LEV(J)) LEV(J)=M D2V(M,J)=D2V(M,J)+D1W1(M2)*D1V(M1,J) 34 CONTINUE ENDIF ENDIF 40 CONTINUE 10 CONTINUE RETURN END C----------------------------------------------------------------------- C--- FDM operators for viscous terms at P points --------------------- SUBROUTINE MDQ2PV PARAMETER (NY=64,NY1=NY+1,MM=4,MH=MM/2,MI=MH-1,MP=MH+MI) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /D2P/D2P(-MP:MP,NY) /LSEP/LSP(NY),LEP(NY) COMMON /D1V/D1V(-MI:MH,0:NY) /MSEV/MSV(0:NY),MEV(0:NY) COMMON /D1P/D1P(-MH:MI,NY) /NSEP/NSP(NY),NEP(NY) DO 10 J=1,NY LSP(J)=999 LEP(J)=-999 DO 20 M=-MP,MP D2P(M,J)=0.D0 20 CONTINUE DO 40 M1=NSP(J),NEP(J) JM=J+M1 IF(JM.GE.0.AND.JM.LE.NY) THEN DO 30 M2=MSV(JM),MEV(JM) M=M1+M2 IF(M.LT.LSP(J)) LSP(J)=M IF(M.GT.LEP(J)) LEP(J)=M D2P(M,J)=D2P(M,J)+D1V(M2,JM)*D1P(M1,J) 30 CONTINUE ELSE WRITE(6,*) 'WARNING IN SBR MDQ2P: J=',J,' M1=',M1,' M2=',M2 ENDIF 40 CONTINUE 10 CONTINUE RETURN END C----------------------------------------------------------------------- C--- FDM operators for pressure's Poisson equation ----------------- SUBROUTINE MDQ2PP PARAMETER (NY=64,NY1=NY+1,MM=4,MH=MM/2,MI=MH-1,MP=MH+MI) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /DPP/DPP(-MP:MP,NY),DPC(NY) /D1VN/D1VN(-MI:MH,0:NY) COMMON /MSEV/MSV(0:NY),MEV(0:NY) /LSEP/LSP(NY),LEP(NY) COMMON /D1P/D1P(-MH:MI,NY) /MSEP/MSP(NY),MEP(NY) COMMON /CMPX/CM3X,CM2X,CM1X,C00X,CP1X,CP2X,CP3X COMMON /CMPZ/CM3Z,CM2Z,CM1Z,C00Z,CP1Z,CP2Z,CP3Z DATA ALPHA/1.5D0/ C00=C00X+C00Z DO 10 J=1,NY LSP(J)=999 LEP(J)=-999 DO 20 M=-MP,MP DPP(M,J)=0.D0 20 CONTINUE DO 40 M1=MSP(J),MEP(J) JM=J+M1 IF(JM.GE.0.AND.JM.LE.NY) THEN DO 30 M2=MSV(JM),MEV(JM) M=M1+M2 IF(M.LT.LSP(J)) LSP(J)=M IF(M.GT.LEP(J)) LEP(J)=M DPP(M,J)=DPP(M,J)+D1VN(M2,JM)*D1P(M1,J) 30 CONTINUE ELSE WRITE(6,*) 'WARNING IN SBR MDQ2P: J=',J,' M1=',M1,' M2=',M2 ENDIF 40 CONTINUE DPC(J)=-ALPHA/(DPP(0,J)+C00) 10 CONTINUE RETURN END C ********************************************************************* C * SBR. DTR : DATA READ * C ********************************************************************* SUBROUTINE SBRDTR(CT) PARAMETER (NX=64,NY=64,NZ=64,NY1=NY+1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) CHARACTER*3 CT CHARACTER*4 FILE1 COMMON /FILE1/FILE1 /START/ISTART /TSTEP/TSTEP COMMON /U/U(NZ,NX,NY) /V/V(NZ,NX,0:NY) /W/W(NZ,NX,NY) COMMON /P/P(NZ,NX,NY) OPEN(10,file=FILE1//CT//'u.d',status='old',form='formatted') READ(10,1000) ISTART,TSTEP READ(10,2000) U CLOSE(10) OPEN(10,file=FILE1//CT//'v.d',status='old',form='formatted') READ(10,1000) ISTART,TSTEP READ(10,2000) V CLOSE(10) OPEN(10,file=FILE1//CT//'w.d',status='old',form='formatted') READ(10,1000) ISTART,TSTEP READ(10,2000) W CLOSE(10) OPEN(10,file=FILE1//CT//'p.d',status='old',form='formatted') READ(10,1000) ISTART,TSTEP READ(10,3000) P CLOSE(10) 1000 FORMAT(I10,F12.6) 2000 FORMAT(8F10.6) 3000 FORMAT(8F10.5) WRITE(6,1200) ISTART,TSTEP 1200 FORMAT(/'Data READ',10X,'STEP=',I7,10X,'T=',F10.3) RETURN END C ********************************************************************* C * SBR. CON : NONLINEAR TERM * C ********************************************************************* SUBROUTINE SBRCON CALL SBRNLU CALL SBRNLV CALL SBRNLW RETURN END C ---------------------------------------------------------------------- SUBROUTINE SBRNLU PARAMETER (NX=64,NY=64,NZ=64,MM=4,MH=MM/2,MI=MH-1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /UF/UF(NZ,NX,NY) COMMON /U/U(NZ,NX,NY) /V/V(NZ,NX,0:NY) /W/W(NZ,NX,NY) COMMON /W1/CX(NZ,NX,NY) /W2/CY(NZ,NX,0:NY) /W3/CZ(NZ,NX,NY) COMMON /TXZA/TXA,TZA /BXZA/BXA,BZA /USI/E0 COMMON /MSEV/MSV(0:NY),MEV(0:NY) /MSEP/MSP(NY),MEP(NY) COMMON /D1V/D1V(-MI:MH,0:NY) /D0P/D0P(-MH:MI,NY) COMMON /IP/IP1(NX),IP2(NX),IP3(NX) /IM/IM1(NX),IM2(NX),IM3(NX) COMMON /KP/KP1(NZ),KP2(NZ),KP3(NZ) /KM/KM1(NZ),KM2(NZ),KM3(NZ) C ... CX:-U^XU_X at P, CZ:-W^XU_Z at P DO 10 J=1,NY DO 10 I=1,NX DO 10 K=1,NZ CX(K,I,J)=-TXA* & (-U(K,IM2(I),J)+ 9.D0*U(K,IM1(I),J)+ 9.D0*U(K,I,J)-U(K,IP1(I),J)) &*( U(K,IM2(I),J)-27.D0*U(K,IM1(I),J)+27.D0*U(K,I,J)-U(K,IP1(I),J)) CZ(K,I,J)=-TZA* & (-W(K,IM1(I),J)+ 9.D0*W(K,I,J)+ 9.D0*W(K,IP1(I),J)-W(K,IP2(I),J)) &*( U(KM1(K),I,J)-27.D0*U(K,I,J)+27.D0*U(KP1(K),I,J)-U(KP2(K),I,J)) 10 CONTINUE C ... CY:U_Y at UV DO 31 J=0,NY DO 31 I=1,NX DO 31 K=1,NZ CY(K,I,J)=0.D0 31 CONTINUE DO 30 J=1,NY-1 DO 30 M=MSV(J),MEV(J) JM=J+M D1J=D1V(M,J) DO 30 I=1,NX DO 30 K=1,NZ CY(K,I,J)=CY(K,I,J)+D1J*U(K,I,JM) 30 CONTINUE C ... CY:-V^XU_Y at UV DO 20 J=1,NY-1 DO 20 I=1,NX DO 20 K=1,NZ CY(K,I,J)=-CY(K,I,J)*E0 &*(-V(K,IM1(I),J)+ 9.D0*V(K,I,J)+ 9.D0*V(K,IP1(I),J)-V(K,IP2(I),J)) 20 CONTINUE C ... UF:-(U^XU_X)^X-(W^XU_Z)^Z at U DO 50 J=1,NY DO 50 I=1,NX DO 50 K=1,NZ UF(K,I,J)=E0*( &-CX(K,IM1(I),J)+9.D0*CX(K,I,J)+9.D0*CX(K,IP1(I),J)-CX(K,IP2(I),J) &-CZ(KM2(K),I,J)+9.D0*CZ(KM1(K),I,J)+9.D0*CZ(K,I,J)-CZ(KP1(K),I,J)) 50 CONTINUE C ... UF:-(U^XU_X)^X-(V^XU_Y)^Y-(W^XU_Z)^Z at U DO 60 J=1,NY DO 60 M=MSP(J),MEP(J) JM=J+M D0J=D0P(M,J) DO 60 I=1,NX DO 60 K=1,NZ UF(K,I,J)=UF(K,I,J)+D0J*CY(K,I,JM) 60 CONTINUE RETURN END C ---------------------------------------------------------------------- SUBROUTINE SBRNLV PARAMETER (NX=64,NY=64,NZ=64,MM=4,MH=MM/2,MI=MH-1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /VF/VF(NZ,NX,0:NY) COMMON /U/U(NZ,NX,NY) /V/V(NZ,NX,0:NY) /W/W(NZ,NX,NY) COMMON /W1/CX(NZ,NX,NY) /W2/CY(NZ,NX,0:NY) /W3/CZ(NZ,NX,NY) COMMON /TXZA/TXA,TZA /BXZA/BXA,BZA /USI/E0 COMMON /MSEV/MSV(0:NY),MEV(0:NY) /MSEP/MSP(NY),MEP(NY) COMMON /D0V/D0V(-MI:MH,0:NY) COMMON /D0P/D0P(-MH:MI,NY) /D1P/D1P(-MH:MI,NY) COMMON /IP/IP1(NX),IP2(NX),IP3(NX) /IM/IM1(NX),IM2(NX),IM3(NX) COMMON /KP/KP1(NZ),KP2(NZ),KP3(NZ) /KM/KM1(NZ),KM2(NZ),KM3(NZ) C ... CX:V^Y, CZ:V_Y at P DO 21 J=1,NY DO 21 I=1,NX DO 21 K=1,NZ CX(K,I,J)=0.D0 CZ(K,I,J)=0.D0 21 CONTINUE DO 20 J=1,NY DO 20 M=MSP(J),MEP(J) JM=J+M D0J=D0P(M,J) D1J=D1P(M,J) DO 20 I=1,NX DO 20 K=1,NZ CX(K,I,J)=CX(K,I,J)+D0J*V(K,I,JM) CZ(K,I,J)=CZ(K,I,J)+D1J*V(K,I,JM) 20 CONTINUE C ... CY:-(V^Y*V_Y)^Y at P DO 22 J=1,NY DO 22 I=1,NX DO 22 K=1,NZ CY(K,I,J)=-CX(K,I,J)*CZ(K,I,J) 22 CONTINUE C ... CX:U^Y at UV, CZ:W^Y at VW DO 31 J=1,NY-1 DO 31 I=1,NX DO 31 K=1,NZ CX(K,I,J)=0.D0 CZ(K,I,J)=0.D0 31 CONTINUE DO 30 J=1,NY-1 DO 30 M=MSV(J),MEV(J) JM=J+M D0J=D0V(M,J) DO 30 I=1,NX DO 30 K=1,NZ CX(K,I,J)=CX(K,I,J)+D0J*U(K,I,JM) CZ(K,I,J)=CZ(K,I,J)+D0J*W(K,I,JM) 30 CONTINUE C ... CX:-(U^Y*V_X) at UV, CZ:-(W^Y*V_Z) at VW DO 10 J=1,NY-1 DO 10 I=1,NX DO 10 K=1,NZ CX(K,I,J)=-CX(K,I,J)*BXA* & (V(K,IM1(I),J)-27.D0*V(K,I,J)+27.D0*V(K,IP1(I),J)-V(K,IP2(I),J)) CZ(K,I,J)=-CZ(K,I,J)*BZA* & (V(KM1(K),I,J)-27.D0*V(K,I,J)+27.D0*V(KP1(K),I,J)-V(KP2(K),I,J)) 10 CONTINUE C ... VF:-(U^Y*V_X)^X-(W^Y*V_Z)^Z at V DO 50 J=1,NY-1 DO 50 I=1,NX DO 50 K=1,NZ VF(K,I,J)=E0*( &-CX(K,IM2(I),J)+9.D0*CX(K,IM1(I),J)+9.D0*CX(K,I,J)-CX(K,IP1(I),J) &-CZ(KM2(K),I,J)+9.D0*CZ(KM1(K),I,J)+9.D0*CZ(K,I,J)-CZ(KP1(K),I,J)) 50 CONTINUE C ... VF:-(U^Y*V_X)^X-(V^Y*V_Y)^Y-(W^Y*V_Z)^Z at V DO 60 J=1,NY-1 DO 60 M=MSV(J),MEV(J) JM=J+M D0J=D0V(M,J) DO 60 I=1,NX DO 60 K=1,NZ VF(K,I,J)=VF(K,I,J)+D0J*CY(K,I,JM) 60 CONTINUE RETURN END C ---------------------------------------------------------------------- SUBROUTINE SBRNLW PARAMETER (NX=64,NY=64,NZ=64,MM=4,MH=MM/2,MI=MH-1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /WF/WF(NZ,NX,NY) COMMON /U/U(NZ,NX,NY) /V/V(NZ,NX,0:NY) /W/W(NZ,NX,NY) COMMON /W1/CX(NZ,NX,NY) /W2/CY(NZ,NX,0:NY) /W3/CZ(NZ,NX,NY) COMMON /TXZA/TXA,TZA /BXZA/BXA,BZA /USI/E0 COMMON /MSEV/MSV(0:NY),MEV(0:NY) /MSEP/MSP(NY),MEP(NY) COMMON /D1V/D1V(-MI:MH,0:NY) /D0P/D0P(-MH:MI,NY) COMMON /IP/IP1(NX),IP2(NX),IP3(NX) /IM/IM1(NX),IM2(NX),IM3(NX) COMMON /KP/KP1(NZ),KP2(NZ),KP3(NZ) /KM/KM1(NZ),KM2(NZ),KM3(NZ) C ... CX:-U^ZW_X at UW, CZ:-W^ZW_Z at P DO 10 J=1,NY DO 10 I=1,NX DO 10 K=1,NZ CX(K,I,J)=-TXA* & (-U(KM1(K),I,J)+ 9.D0*U(K,I,J)+ 9.D0*U(KP1(K),I,J)-U(KP2(K),I,J)) &*( W(K,IM1(I),J)-27.D0*W(K,I,J)+27.D0*W(K,IP1(I),J)-W(K,IP2(I),J)) CZ(K,I,J)=-TZA* & (-W(KM2(K),I,J)+ 9.D0*W(KM1(K),I,J)+ 9.D0*W(K,I,J)-W(KP1(K),I,J)) &*( W(KM2(K),I,J)-27.D0*W(KM1(K),I,J)+27.D0*W(K,I,J)-W(KP1(K),I,J)) 10 CONTINUE C ... CY:W_Y at VW DO 31 J=0,NY DO 31 I=1,NX DO 31 K=1,NZ CY(K,I,J)=0.D0 31 CONTINUE DO 30 J=1,NY-1 DO 30 M=MSV(J),MEV(J) JM=J+M D1J=D1V(M,J) DO 30 I=1,NX DO 30 K=1,NZ CY(K,I,J)=CY(K,I,J)+D1J*W(K,I,JM) 30 CONTINUE C ... CY:-V^ZW_Y at VW DO 20 J=1,NY-1 DO 20 I=1,NX DO 20 K=1,NZ CY(K,I,J)=-CY(K,I,J)*E0 &*(-V(KM1(K),I,J)+ 9.D0*V(K,I,J)+ 9.D0*V(KP1(K),I,J)-V(KP2(K),I,J)) 20 CONTINUE C ... WF:-(U^ZW_X)^X-(W^ZW_Z)^Z at W DO 50 J=1,NY DO 50 I=1,NX DO 50 K=1,NZ WF(K,I,J)=E0*( &-CX(K,IM2(I),J)+9.D0*CX(K,IM1(I),J)+9.D0*CX(K,I,J)-CX(K,IP1(I),J) &-CZ(KM1(K),I,J)+9.D0*CZ(K,I,J)+9.D0*CZ(KP1(K),I,J)-CZ(KP2(K),I,J)) 50 CONTINUE C ... WF:-(U^ZW_X)^X-(V^ZW_Y)^Y-(W^ZW_Z)^Z at W DO 60 J=1,NY DO 60 M=MSP(J),MEP(J) JM=J+M D0J=D0P(M,J) DO 60 I=1,NX DO 60 K=1,NZ WF(K,I,J)=WF(K,I,J)+D0J*CY(K,I,JM) 60 CONTINUE RETURN END C ********************************************************************* C * SBR. VIS : VISCOUS TERM * C ********************************************************************* SUBROUTINE SBRVIS PARAMETER (NX=64,NY=64,NZ=64,NY1=NY+1) PARAMETER (MM=4,MH=MM/2,MI=MH-1,MP=MH+MI) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /U/U(NZ,NX,NY) /V/V(NZ,NX,0:NY) /W/W(NZ,NX,NY) COMMON /UF/UF(NZ,NX,NY) /VF/VF(NZ,NX,0:NY) /WF/WF(NZ,NX,NY) COMMON /RET/RET /D2V/D2V(-MP:MP,0:NY) /D2P/D2P(-MP:MP,NY) COMMON /LSEV/LSV(0:NY),LEV(0:NY) /LSEP/LSP(NY),LEP(NY) COMMON /IP/IP1(NX),IP2(NX),IP3(NX) /IM/IM1(NX),IM2(NX),IM3(NX) COMMON /KP/KP1(NZ),KP2(NZ),KP3(NZ) /KM/KM1(NZ),KM2(NZ),KM3(NZ) COMMON /CMPX/CM3X,CM2X,CM1X,C00X,CP1X,CP2X,CP3X COMMON /CMPZ/CM3Z,CM2Z,CM1Z,C00Z,CP1Z,CP2Z,CP3Z C00=C00X+C00Z C ... UF:U_XX+U_ZZ, WF:W_XX+W_ZZ DO 10 J=1,NY DO 10 I=1,NX DO 10 K=1,NZ UF(K,I,J)=UF(K,I,J) &+(CM3X*U(K,IM3(I),J)+CM2X*U(K,IM2(I),J)+CM1X*U(K,IM1(I),J) & +CM3Z*U(KM3(K),I,J)+CM2Z*U(KM2(K),I,J)+CM1Z*U(KM1(K),I,J) & +C00*U(K,I,J) & +CP1Z*U(KP1(K),I,J)+CP2Z*U(KP2(K),I,J)+CP3Z*U(KP3(K),I,J) & +CP1X*U(K,IP1(I),J)+CP2X*U(K,IP2(I),J)+CP3X*U(K,IP3(I),J))/RET WF(K,I,J)=WF(K,I,J) &+(CM3X*W(K,IM3(I),J)+CM2X*W(K,IM2(I),J)+CM1X*W(K,IM1(I),J) & +CM3Z*W(KM3(K),I,J)+CM2Z*W(KM2(K),I,J)+CM1Z*W(KM1(K),I,J) & +C00*W(K,I,J) & +CP1Z*W(KP1(K),I,J)+CP2Z*W(KP2(K),I,J)+CP3Z*W(KP3(K),I,J) & +CP1X*W(K,IP1(I),J)+CP2X*W(K,IP2(I),J)+CP3X*W(K,IP3(I),J))/RET 10 CONTINUE C .. UF:U_XX+U_YY+U_ZZ, WF:W_XX+W_YY+W_ZZ DO 20 J=1,NY DO 20 L=LSP(J),LEP(J) JL=J+L BNU=D2P(L,J)/RET DO 20 I=1,NX DO 20 K=1,NZ UF(K,I,J)=UF(K,I,J)+BNU*U(K,I,JL) WF(K,I,J)=WF(K,I,J)+BNU*W(K,I,JL) 20 CONTINUE C ... VF:V_XX+V_ZZ DO 30 J=1,NY-1 DO 30 I=1,NX DO 30 K=1,NZ VF(K,I,J)=VF(K,I,J) &+(CM3X*V(K,IM3(I),J)+CM2X*V(K,IM2(I),J)+CM1X*V(K,IM1(I),J) & +CM3Z*V(KM3(K),I,J)+CM2Z*V(KM2(K),I,J)+CM1Z*V(KM1(K),I,J) & +C00*V(K,I,J) & +CP1Z*V(KP1(K),I,J)+CP2Z*V(KP2(K),I,J)+CP3Z*V(KP3(K),I,J) & +CP1X*V(K,IP1(I),J)+CP2X*V(K,IP2(I),J)+CP3X*V(K,IP3(I),J))/RET 30 CONTINUE C ... VF:V_XX+V_YY+V_ZZ DO 40 J=1,NY-1 DO 40 L=LSV(J),LEV(J) JL=J+L BNU=D2V(L,J)/RET DO 40 I=1,NX DO 40 K=1,NZ VF(K,I,J)=VF(K,I,J)+BNU*V(K,I,JL) 40 CONTINUE RETURN END C ********************************************************************* C * SBR. PRE : PREDICTION STEP * C ********************************************************************* SUBROUTINE SBRPRE(AB,BB) PARAMETER (NX=64,NY=64,NZ=64) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /U/U(NZ,NX,NY) /V/V(NZ,NX,0:NY) /W/W(NZ,NX,NY) COMMON /UF/UF(NZ,NX,NY) /VF/VF(NZ,NX,0:NY) /WF/WF(NZ,NX,NY) COMMON /UB/UB(NZ,NX,NY) /VB/VB(NZ,NX,NY-1) /WB/WB(NZ,NX,NY) COMMON /DT/DT DAB=AB*DT DBB=BB*DT DO 10 J=1,NY DO 10 I=1,NX DO 10 K=1,NZ UF(K,I,J)=UF(K,I,J)+2.D0 U(K,I,J)=U(K,I,J)+DAB*UF(K,I,J)+DBB*UB(K,I,J) W(K,I,J)=W(K,I,J)+DAB*WF(K,I,J)+DBB*WB(K,I,J) UB(K,I,J)=UF(K,I,J) WB(K,I,J)=WF(K,I,J) 10 CONTINUE DO 20 J=1,NY-1 DO 20 I=1,NX DO 20 K=1,NZ V(K,I,J)=V(K,I,J)+DAB*VF(K,I,J)+DBB*VB(K,I,J) VB(K,I,J)=VF(K,I,J) 20 CONTINUE RETURN END C ********************************************************************* C * SBR. RHP : R.H.S. OF POISSON EQ. * C ********************************************************************* C --- R.H.S. FOR SCALER POTANTIAL ------------------------------------- SUBROUTINE SBRRHP PARAMETER (NX=64,NY=64,NZ=64,MM=4,MH=MM/2,MI=MH-1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /U/U(NZ,NX,NY) /V/V(NZ,NX,0:NY) /W/W(NZ,NX,NY) COMMON /W3/Q(NZ,NX,NY) /DT/DT COMMON /D1P/D1P(-MH:MI,NY) /MSEP/MSP(NY),MEP(NY) COMMON /IP/IP1(NX),IP2(NX),IP3(NX) /IM/IM1(NX),IM2(NX),IM3(NX) COMMON /KP/KP1(NZ),KP2(NZ),KP3(NZ) /KM/KM1(NZ),KM2(NZ),KM3(NZ) COMMON /AMX/AM2X,AM1X,A00X,AP1X /AMZ/AM2Z,AM1Z,A00Z,AP1Z DO 10 J=1,NY DO 10 I=1,NX DO 10 K=1,NZ Q(K,I,J)= & AM2X*U(K,IM2(I),J)+AM1X*U(K,IM1(I),J) & +A00X*U(K,I,J)+AP1X*U(K,IP1(I),J) & +AM2Z*W(KM2(K),I,J)+AM1Z*W(KM1(K),I,J) & +A00Z*W(K,I,J)+AP1Z*W(KP1(K),I,J) 10 CONTINUE DO 20 J=1,NY DO 20 M=MSP(J),MEP(J) JM=J+M D1J=D1P(M,J) DO 20 I=1,NX DO 20 K=1,NZ Q(K,I,J)=Q(K,I,J)+D1J*V(K,I,JM) 20 CONTINUE DO 30 J=1,NY DO 30 I=1,NX DO 30 K=1,NZ Q(K,I,J)=Q(K,I,J)/DT 30 CONTINUE RETURN END C --- R.H.S. FOR PRESSURE --------------------------------------------- SUBROUTINE SBRRHW PARAMETER (NX=64,NY=64,NZ=64,MM=4,MH=MM/2,MI=MH-1,MP=MH+MI) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /VISWALL/S0W(NZ,NX),S1W(NZ,NX) COMMON /UF/UF(NZ,NX,NY) /VF/VF(NZ,NX,0:NY) /WF/WF(NZ,NX,NY) COMMON /V/V(NZ,NX,0:NY) COMMON /W3/Q(NZ,NX,NY) COMMON /D1P/D1P(-MH:MI,NY) COMMON /NSEP/NSP(NY),NEP(NY) /MSEP/MSP(NY),MEP(NY) COMMON /IP/IP1(NX),IP2(NX),IP3(NX) /IM/IM1(NX),IM2(NX),IM3(NX) COMMON /KP/KP1(NZ),KP2(NZ),KP3(NZ) /KM/KM1(NZ),KM2(NZ),KM3(NZ) COMMON /AMX/AM2X,AM1X,A00X,AP1X /AMZ/AM2Z,AM1Z,A00Z,AP1Z COMMON /D2V/D2V(-MP:MP,0:NY) /LSEV/LSV(0:NY),LEV(0:NY) COMMON /G1W/G1WA(0:MH),G1WB(NY-MH:NY) /RET/RET DO 10 J=1,NY DO 10 I=1,NX DO 10 K=1,NZ Q(K,I,J)= & AM2X*UF(K,IM2(I),J)+AM1X*UF(K,IM1(I),J) & +A00X*UF(K,I,J)+AP1X*UF(K,IP1(I),J) & +AM2Z*WF(KM2(K),I,J)+AM1Z*WF(KM1(K),I,J) & +A00Z*WF(K,I,J)+AP1Z*WF(KP1(K),I,J) 10 CONTINUE DO 20 J=1,NY DO 20 M=MSP(J),MEP(J) JM=J+M D1J=D1P(M,J) DO 20 I=1,NX DO 20 K=1,NZ Q(K,I,J)=Q(K,I,J)+D1J*VF(K,I,JM) 20 CONTINUE C --- For Wall B.C. (Viscous Terms) ------------- DO 30 I=1,NX DO 30 K=1,NZ S0W(K,I)=0.D0 DO 31 L=LSV(0),LEV(0) S0W(K,I)=S0W(K,I)+D2V(L,0)*V(K,I,L) 31 CONTINUE S0W(K,I)=S0W(K,I)/RET S1W(K,I)=0.D0 DO 32 L=LSV(NY),LEV(NY) S1W(K,I)=S1W(K,I)+D2V(L,NY)*V(K,I,NY+L) 32 CONTINUE S1W(K,I)=S1W(K,I)/RET 30 CONTINUE DO 40 J=1,NY DO 40 M=NSP(J),NEP(J) JM=J+M IF(JM.LE.MH) THEN D1J=D1P(M,J)*G1WA(JM) DO 41 I=1,NX DO 41 K=1,NZ Q(K,I,J)=Q(K,I,J)-D1J*S0W(K,I) 41 CONTINUE ENDIF IF(JM.GE.NY-MH) THEN D1J=D1P(M,J)*G1WB(JM-NY) DO 42 I=1,NX DO 42 K=1,NZ Q(K,I,J)=Q(K,I,J)-D1J*S1W(K,I) 42 CONTINUE ENDIF 40 CONTINUE RETURN END C ********************************************************************* C * SBR. SOR : S.O.R. SCHEME FOR POISSON EQ. * C ********************************************************************* SUBROUTINE SBRSOR(ISOR) PARAMETER (NX=64,NY=64,NZ=64,NY1=NY+1) PARAMETER (MM=4,MH=MM/2,MI=MH-1,MP=MH+MI) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /DPP/DPP(-MP:MP,NY),DPC(NY) /LSEP/LSP(NY),LEP(NY) COMMON /P/P(NZ,NX,NY) /W3/Q(NZ,NX,NY) /ERRP/POIERR /ITRP/ITRP COMMON /IP/IP1(NX),IP2(NX),IP3(NX) /IM/IM1(NX),IM2(NX),IM3(NX) COMMON /KP/KP1(NZ),KP2(NZ),KP3(NZ) /KM/KM1(NZ),KM2(NZ),KM3(NZ) COMMON /CMPX/CM3X,CM2X,CM1X,C00X,CP1X,CP2X,CP3X COMMON /CMPZ/CM3Z,CM2Z,CM1Z,C00Z,CP1Z,CP2Z,CP3Z C00=C00X+C00Z ITRP=0 SUMS=0.D0 DO 10 J=1,NY DO 10 I=1,NX DO 10 K=1,NZ SUMS=SUMS+Q(K,I,J)**2 10 CONTINUE 100 ITRP=ITRP+1 SUMR=0. DO 30 J=1,NY DO 30 I=1,NX DO 30 K=1,NZ RESI=CM3X*P(K,IM3(I),J)+CM2X*P(K,IM2(I),J)+CM1X*P(K,IM1(I),J) & +CM3Z*P(KM3(K),I,J)+CM2Z*P(KM2(K),I,J)+CM1Z*P(KM1(K),I,J) & +C00*P(K,I,J) & +CP1Z*P(KP1(K),I,J)+CP2Z*P(KP2(K),I,J)+CP3Z*P(KP3(K),I,J) & +CP1X*P(K,IP1(I),J)+CP2X*P(K,IP2(I),J)+CP3X*P(K,IP3(I),J) & -Q(K,I,J) DO 40 M=LSP(J),LEP(J) RESI=RESI+DPP(M,J)*P(K,I,J+M) 40 CONTINUE P(K,I,J)=P(K,I,J)+RESI*DPC(J) SUMR=SUMR+RESI**2 30 CONTINUE POIERR=DSQRT(SUMR/SUMS) IF(ITRP.GE.ISOR.OR.POIERR.LT.1.D-5) GO TO 200 GO TO 100 200 CONTINUE RETURN END C ********************************************************************* C * SBR. COR : CORRECTION STEP * C ********************************************************************* SUBROUTINE SBRCOR PARAMETER (NX=64,NY=64,NZ=64,NY1=NY+1,MM=4,MH=MM/2,MI=MH-1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /U/U(NZ,NX,NY) /V/V(NZ,NX,0:NY) /W/W(NZ,NX,NY) COMMON /P/P(NZ,NX,NY) COMMON /DT/DT /D/DX,DZ /DY/DY(NY) COMMON /D1VN/D1VN(-MI:MH,0:NY) /MSEV/MSV(0:NY),MEV(0:NY) COMMON /IP/IP1(NX),IP2(NX),IP3(NX) /IM/IM1(NX),IM2(NX),IM3(NX) COMMON /KP/KP1(NZ),KP2(NZ),KP3(NZ) /KM/KM1(NZ),KM2(NZ),KM3(NZ) COMMON /BMX/BM1X,B00X,BP1X,BP2X /BMZ/BM1Z,B00Z,BP1Z,BP2Z TX=DT/DX TZ=DT/DZ DO 10 J=1,NY DO 10 I=1,NX DO 10 K=1,NZ U(K,I,J)=U(K,I,J)-DT*(BM1X*P(K,IM1(I),J) & +B00X*P(K,I,J)+BP1X*P(K,IP1(I),J)+BP2X*P(K,IP2(I),J)) W(K,I,J)=W(K,I,J)-DT*(BM1Z*P(KM1(K),I,J) & +B00Z*P(K,I,J)+BP1Z*P(KP1(K),I,J)+BP2Z*P(KP2(K),I,J)) 10 CONTINUE DO 20 J=1,NY-1 DO 20 M=MSV(J),MEV(J) JM=J+M TY=DT*D1VN(M,J) DO 20 I=1,NX DO 20 K=1,NZ V(K,I,J)=V(K,I,J)-TY*P(K,I,JM) 20 CONTINUE RETURN END C ********************************************************************* C * SBR. UMR : MEAN & RMS VALUES * C ********************************************************************* SUBROUTINE SBRUMR PARAMETER (NX=64,NY=64,NZ=64,NY1=NY+1,MM=4,MH=MM/2,MI=MH-1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /U/U(NZ,NX,NY) /V/V(NZ,NX,0:NY) /W/W(NZ,NX,NY) COMMON /UM/UM(NY),VM(0:NY),WM(NY) COMMON /S12/S12L(0:NY),S12V(0:NY) COMMON /IP/IP1(NX),IP2(NX),IP3(NX) /IM/IM1(NX),IM2(NX),IM3(NX) COMMON /MSEV/MSV(0:NY),MEV(0:NY) /D0V/D0V(-MI:MH,0:NY) COMMON /BVU/BVU2(NY),BVU3(NY),BVU4(NY) COMMON /BVV/BVV2(0:NY),BVV3(0:NY),BVV4(0:NY) COMMON /BVW/BVW2(NY),BVW3(NY),BVW4(NY) ARXZ=1./DBLE(NX*NZ) DO 20 J=1,NY SUMU1=0.D0 SUMU2=0.D0 SUMW1=0.D0 SUMW2=0.D0 DO 10 I=1,NX DO 10 K=1,NZ SUMU1=SUMU1+U(K,I,J) SUMU2=SUMU2+U(K,I,J)**2 SUMW1=SUMW1+W(K,I,J) SUMW2=SUMW2+W(K,I,J)**2 10 CONTINUE UM(J)=ARXZ*SUMU1 BVU2(J)=ARXZ*SUMU2 WM(J)=ARXZ*SUMW1 BVW2(J)=ARXZ*SUMW2 20 CONTINUE DO 50 J=1,NY-1 SUMV1=0.D0 SUMV2=0.D0 SUMUV=0.D0 DO 40 I=1,NX DO 40 K=1,NZ SUMV1=SUMV1+V(K,I,J) SUMV2=SUMV2+V(K,I,J)**2 U12=0.D0 DO 30 M=MSV(J),MEV(J) U12=U12+D0V(M,J)*U(K,I,J+M) 30 CONTINUE V12=(-V(K,IM1(I),J)+9.D0*(V(K,I,J)+V(K,IP1(I),J))-V(K,IP2(I),J)) & /16.D0 SUMUV=SUMUV+U12*V12 40 CONTINUE VM(J)=ARXZ*SUMV1 BVV2(J)=ARXZ*SUMV2 S12L(J)=ARXZ*SUMUV 50 CONTINUE RETURN END C ********************************************************************* C * SBR. STT : TURBULENCE STATISTICS * C ********************************************************************* C --- Mean Velocity and Turbulence Intensity---------------------------- SUBROUTINE SBRST1 PARAMETER (NX=64,NY=64,NZ=64,NY1=NY+1) PARAMETER (MM=4,MH=MM/2,MI=MH-1,MP=MI+MH) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /UM/UM(NY),VM(0:NY),WM(NY) /PM/PM(NY) COMMON /UR/UR(NY),VR(0:NY),WR(NY) /PR/PR(NY) COMMON /ENE/ENE /UME/UME,UMX /RMS/URMSX,URMSC,VRMSC,WRMSC COMMON /DY/DY(NY) /DYC/DYC(NY-1) COMMON /BVU/BVU2(NY),BVU3(NY),BVU4(NY) COMMON /BVV/BVV2(0:NY),BVV3(0:NY),BVV4(0:NY) COMMON /BVW/BVW2(NY),BVW3(NY),BVW4(NY) COMMON /BVP/BVP2(NY) ENE=0.D0 UME=0.D0 UMX=0.D0 URMSX=0.D0 DO 10 J=1,NY UR(J)=DSQRT(BVU2(J)-UM(J)**2) WR(J)=DSQRT(BVW2(J)-WM(J)**2) PR(J)=DSQRT(BVP2(J)-PM(J)**2) ENE=ENE+DY(J)*(UR(J)**2+WR(J)**2) UME=UME+DY(J)*UM(J) UMX=DMAX1(UMX,UM(J)) URMSX=DMAX1(URMSX,UR(J)) 10 CONTINUE DO 20 J=1,NY-1 VR(J)=DSQRT(BVV2(J)-VM(J)**2) ENE=ENE+DYC(J)*VR(J)**2 20 CONTINUE URMSC=(UR(NY/2)+UR(NY/2+1))/2.D0 VRMSC=VR(NY/2) WRMSC=(WR(NY/2)+WR(NY/2+1))/2.D0 ENE=ENE/2.D0 RETURN END C ********************************************************************* C * SBR. CHK : CHECK of DIVERGENCE and COURANT-NUMBER * C ********************************************************************* SUBROUTINE SBRCHK PARAMETER (NX=64,NY=64,NZ=64,MM=4,MH=MM/2,MI=MH-1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /U/U(NZ,NX,NY) /V/V(NZ,NX,0:NY) /W/W(NZ,NX,NY) COMMON /D/DX,DZ /DY/DY(NY) /DT/DT /DIV/DIVX /COU/COMX COMMON /D0P/D0P(-MH:MI,NY) /D1P/D1P(-MH:MI,NY) COMMON /MSEP/MSP(NY),MEP(NY) COMMON /UR/UR(NY),VR(0:NY),WR(NY) COMMON /IP/IP1(NX),IP2(NX),IP3(NX) /IM/IM1(NX),IM2(NX),IM3(NX) COMMON /KP/KP1(NZ),KP2(NZ),KP3(NZ) /KM/KM1(NZ),KM2(NZ),KM3(NZ) COMMON /AMX/AM2X,AM1X,A00X,AP1X /AMZ/AM2Z,AM1Z,A00Z,AP1Z DIVX=0.D0 COMX=0.D0 DO 40 J=1,NY DNOMAL=(DX*DY(J)*DZ)**(1.D0/3.D0)/UR(J) DO 40 I=1,NX DO 40 K=1,NZ DIV=AM2X*U(K,IM2(I),J)+AM1X*U(K,IM1(I),J) & +A00X*U(K,I,J)+AP1X*U(K,IP1(I),J) & +AM2Z*W(KM2(K),I,J)+AM1Z*W(KM1(K),I,J) & +A00Z*W(K,I,J)+AP1Z*W(KP1(K),I,J) UCP=(-U(K,IM2(I),J)+9.D0*U(K,IM1(I),J) & +9.D0*U(K,I,J)-U(K,IP1(I),J))/16.D0 WCP=(-W(KM2(K),I,J)+9.D0*W(KM1(K),I,J) & +9.D0*W(K,I,J)-W(KP1(K),I,J))/16.D0 VCP=0.D0 DO 20 M=MSP(J),MEP(J) VCP=VCP+D0P(M,J)*V(K,I,J+M) DIV=DIV+D1P(M,J)*V(K,I,J+M) 20 CONTINUE DIVX=DMAX1(DIVX,DNOMAL*DIV) COU=DT*(DABS(UCP)/DX+DABS(VCP)/DY(J)+DABS(WCP)/DZ) COMX=DMAX1(COMX,COU) 40 CONTINUE RETURN END C ********************************************************************* C * SBR. DTS : DATA SAVE * C ********************************************************************* SUBROUTINE SBRDTS(CT) PARAMETER (NX=64,NY=64,NZ=64,NY1=NY+1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) CHARACTER*3 CT CHARACTER*4 FILE1 COMMON /FILE1/FILE1 /STEP/ISTEP /TSTEP/TSTEP COMMON /U/U(NZ,NX,NY) /V/V(NZ,NX,0:NY) /W/W(NZ,NX,NY) COMMON /P/P(NZ,NX,NY) OPEN(10,file=FILE1//CT//'u.d',status='new',form='formatted') WRITE(10,1000) ISTEP,TSTEP WRITE(10,2000) U CLOSE(10) OPEN(10,file=FILE1//CT//'v.d',status='new',form='formatted') WRITE(10,1000) ISTEP,TSTEP WRITE(10,2000) V CLOSE(10) OPEN(10,file=FILE1//CT//'w.d',status='new',form='formatted') WRITE(10,1000) ISTEP,TSTEP WRITE(10,2000) W CLOSE(10) OPEN(10,file=FILE1//CT//'p.d',status='new',form='formatted') WRITE(10,1000) ISTEP,TSTEP WRITE(10,3000) P CLOSE(10) 1000 FORMAT(I10,F12.6) 2000 FORMAT(8F10.6) 3000 FORMAT(8F10.5) WRITE(6,1200) ISTEP,TSTEP 1200 FORMAT(/'Data SAVE',10X,'STEP=',I7,10X,'T=',F10.3/) RETURN END