C *************************************************************** C * DNS37AB2 Ver.1 1994.09.13 (4th order accuracy) * C * Ver.3B 1997.07.12 (2nd order; 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='dns2' 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 MDQ1VN CALL MDQ1VD CALL MDQ1P 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/IP(NX) /IM/IM(NX) /KP/KP(NZ) /KM/KM(NZ) COMMON /TXZA/TXA,TZA /BXZA/BXA,BZA COMMON /CMPX/CM1X,C00X,CP1X /CMPZ/CM1Z,C00Z,CP1Z DX=18.D0/RET DZ= 9.D0/RET HX=NX*DX HZ=NZ*DZ BXA=1.D0/DX BZA=1.D0/DZ TXA=BXA/2.D0 TZA=BZA/2.D0 C --- List vectors to identify the neighbouring point DO 50 I=1,NX IP(I)=I+1 IM(I)=I-1 IF(IP(I).GT.NX) IP(I)=IP(I)-NX IF(IM(I).LT. 1) IM(I)=IM(I)+NX 50 CONTINUE DO 60 K=1,NZ KP(K)=K+1 KM(K)=K-1 IF(KP(K).GT.NZ) KP(K)=KP(K)-NZ IF(KM(K).LT. 1) KM(K)=KM(K)+NZ 60 CONTINUE C --- Parameters for finite-difference of Poisson equation CM1X= 1.D0/DX**2 C00X=-2.D0/DX**2 CP1X= 1.D0/DX**2 CM1Z= 1.D0/DZ**2 C00Z=-2.D0/DZ**2 CP1Z= 1.D0/DZ**2 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) /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) DO 20 J=1,NY DY(J)=-YV(J-1)+YV(J) 20 CONTINUE DO 30 J=1,NY-1 DYC(J)=-YP(J)+YP(J+1) 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) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /Y/YV(0:NY),YP(0:NY1) /D1VN/D1VNM(0:NY),D1VNP(0:NY) DO 10 J=1,NY-1 DYV=YP(J+1)-YP(J) D1VNM(J)=-1.D0/DYV D1VNP(J)= 1.D0/DYV 10 CONTINUE C For the Neumann B.C. at the wall (DP/DY=0) for Pressure D1VNM(0)= 0.D0 D1VNP(0)= 0.D0 D1VNM(NY)= 0.D0 D1VNP(NY)= 0.D0 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) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /Y/YV(0:NY),YP(0:NY1) COMMON /D0V/D0VM(NY-1),D0VP(NY-1) /D1V/D1VM(0:NY),D1VP(0:NY) DO 10 J=1,NY-1 DYV=YP(J+1)-YP(J) D0VM(J)=(YP(J+1)-YV(J))/DYV D0VP(J)=(YV(J)-YP(J))/DYV D1VM(J)=-1.D0/DYV D1VP(J)= 1.D0/DYV 10 CONTINUE Caution! Stencils are shifted at the wall point DYP0=YP(2)-YP(1) D1VM(0)= YP(2)/YP(1)/DYP0 D1VP(0)=-YP(1)/YP(2)/DYP0 DYPN=YP(NY)-YP(NY-1) D1VM(NY)= (1.-YP(NY))/(1.-YP(NY-1))/DYPN D1VP(NY)=-(1.-YP(NY-1))/(1.-YP(NY))/DYPN 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) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /Y/YV(0:NY),YP(0:NY1) /DY/DY(NY) COMMON /D0P/D0PM(NY),D0PP(NY) /D1P/D1PM(NY),D1PP(NY) DO 10 J=1,NY D0PM(J)=(YV(J)-YP(J))/DY(J) D0PP(J)=(YP(J)-YV(J-1))/DY(J) D1PM(J)=-1.D0/DY(J) D1PP(J)= 1.D0/DY(J) 10 CONTINUE RETURN END C----------------------------------------------------------------------- C--- FDM operators for viscous terms at V points --------------------- SUBROUTINE MDQ2V PARAMETER (NY=64,NY1=NY+1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /D2V/D2VM(NY-1),D2V0(NY-1),D2VP(NY-1) COMMON /D1V/D1VM(0:NY),D1VP(0:NY) /D1P/D1PM(NY),D1PP(NY) DO 10 J=1,NY-1 D2VM(J)=D1VM(J)*D1PM(J) D2V0(J)=D1VM(J)*D1PP(J)+D1VP(J)*D1PM(J+1) D2VP(J)=D1VP(J)*D1PP(J+1) 10 CONTINUE RETURN END C----------------------------------------------------------------------- C--- FDM operators for viscous terms at P points --------------------- SUBROUTINE MDQ2PV PARAMETER (NY=64,NY1=NY+1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /D2P/D2PM(2:NY),D2P0(NY),D2PP(NY-1) COMMON /D1V/D1VM(0:NY),D1VP(0:NY) /D1P/D1PM(NY),D1PP(NY) J=1 D2P0(J)=D1PM(J)*D1VM(J-1)+D1PP(J)*D1VM(J) D2PP(J)=D1PM(J)*D1VP(J-1)+D1PP(J)*D1VP(J) DO 10 J=2,NY-1 D2PM(J)=D1PM(J)*D1VM(J-1) D2P0(J)=D1PM(J)*D1VP(J-1)+D1PP(J)*D1VM(J) D2PP(J)= D1PP(J)*D1VP(J) 10 CONTINUE J=NY D2PM(J)=D1PM(J)*D1VM(J-1)+D1PP(J)*D1VM(J) D2P0(J)=D1PM(J)*D1VP(J-1)+D1PP(J)*D1VP(J) RETURN END C----------------------------------------------------------------------- C--- FDM operators for pressure's Poisson equation ----------------- SUBROUTINE MDQ2PP PARAMETER (NY=64,NY1=NY+1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /DPP/DPPM(2:NY),DPP0(NY),DPPP(NY-1) COMMON /D1VN/D1VNM(0:NY),D1VNP(0:NY) /D1P/D1PM(NY),D1PP(NY) DO 10 J=1,NY IF(J.GT. 1) DPPM(J)=D1PM(J)*D1VNM(J-1) DPP0(J)=D1PM(J)*D1VNP(J-1)+D1PP(J)*D1VNM(J) IF(J.LT.NY) DPPP(J)=D1PP(J)*D1VNP(J) 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) 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 COMMON /D1V/D1VM(0:NY),D1VP(0:NY) /D0P/D0PM(NY),D0PP(NY) COMMON /IP/IP(NX) /IM/IM(NX) /KP/KP(NZ) /KM/KM(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,IM(I),J)+U(K,I,J))*(-U(K,IM(I),J)+U(K,I,J)) CZ(K,I,J)=-TZA*(W(K,I,J)+W(K,IP(I),J))*(-U(K,I,J)+U(KP(K),I,J)) 10 CONTINUE C ... CY:-V^XU_Y at UV DO 30 J=0,NY IF(J.GT.0.AND.J.LT.NY) THEN DO 31 I=1,NX DO 31 K=1,NZ CY(K,I,J)=-(D1VM(J)*U(K,I,J)+D1VP(J)*U(K,I,J+1)) & *(V(K,I,J)+V(K,IP(I),J))/2.D0 31 CONTINUE ELSE DO 32 I=1,NX DO 32 K=1,NZ CY(K,I,J)=0.D0 32 CONTINUE ENDIF 30 CONTINUE C ... UF:-(U^XU_X)^X-(V^XU_Y)^Y-(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)=(CX(K,I,J)+CX(K,IP(I),J)+CZ(KM(K),I,J)+CZ(K,I,J))/2.D0 & +D0PM(J)*CY(K,I,J-1)+D0PP(J)*CY(K,I,J) 50 CONTINUE RETURN END C ---------------------------------------------------------------------- SUBROUTINE SBRNLV PARAMETER (NX=64,NY=64,NZ=64) 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 COMMON /D0V/D0VM(NY-1),D0VP(NY-1) COMMON /D0P/D0PM(NY),D0PP(NY) /D1P/D1PM(NY),D1PP(NY) COMMON /IP/IP(NX) /IM/IM(NX) /KP/KP(NZ) /KM/KM(NZ) 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)=-(D0PM(J)*V(K,I,J-1)+D0PP(J)*V(K,I,J)) & *(D1PM(J)*V(K,I,J-1)+D1PP(J)*V(K,I,J)) 22 CONTINUE C ... CX:-(U^Y*V_X) at UV, CZ:-(W^Y*V_Z) at VW DO 30 J=1,NY-1 DO 30 I=1,NX DO 30 K=1,NZ CX(K,I,J)=-(D0VM(J)*U(K,I,J)+D0VP(J)*U(K,I,J+1)) & *BXA*(-V(K,I,J)+V(K,IP(I),J)) CZ(K,I,J)=-(D0VM(J)*W(K,I,J)+D0VP(J)*W(K,I,J+1)) & *BZA*(-V(K,I,J)+V(KP(K),I,J)) 30 CONTINUE C ... VF:-(U^Y*V_X)^X-(V^Y*V_Y)^Y-(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)=(CX(K,IM(I),J)+CX(K,I,J)+CZ(KM(K),I,J)+CZ(K,I,J))/2.D0 & +D0VM(J)*CY(K,I,J)+D0VP(J)*CY(K,I,J+1) 50 CONTINUE RETURN END C ---------------------------------------------------------------------- SUBROUTINE SBRNLW PARAMETER (NX=64,NY=64,NZ=64) 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 COMMON /D1V/D1VM(0:NY),D1VP(0:NY) /D0P/D0PM(NY),D0PP(NY) COMMON /IP/IP(NX) /IM/IM(NX) /KP/KP(NZ) /KM/KM(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(K,I,J)+U(KP(K),I,J))*(-W(K,I,J)+W(K,IP(I),J)) CZ(K,I,J)=-TZA*(W(KM(K),I,J)+W(K,I,J))*(-W(KM(K),I,J)+W(K,I,J)) 10 CONTINUE C ... CY:-V^ZW_Y at VW DO 30 J=0,NY IF(J.GT.0.AND.J.LT.NY) THEN DO 31 I=1,NX DO 31 K=1,NZ CY(K,I,J)=-(D1VM(J)*W(K,I,J)+D1VP(J)*W(K,I,J+1)) & *(V(K,I,J)+V(KP(K),I,J))/2.D0 31 CONTINUE ELSE DO 32 I=1,NX DO 32 K=1,NZ CY(K,I,J)=0.D0 32 CONTINUE ENDIF 30 CONTINUE C ... WF:-(U^ZW_X)^X-(V^ZW_Y)^Y-(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)=(CX(K,IM(I),J)+CX(K,I,J)+CZ(K,I,J)+CZ(KP(K),I,J))/2.D0 & +D0PM(J)*CY(K,I,J-1)+D0PP(J)*CY(K,I,J) 50 CONTINUE RETURN END C ********************************************************************* C * SBR. VIS : VISCOUS TERM * C ********************************************************************* SUBROUTINE SBRVIS PARAMETER (NX=64,NY=64,NZ=64,NY1=NY+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) /RET/RET COMMON /UF/UF(NZ,NX,NY) /VF/VF(NZ,NX,0:NY) /WF/WF(NZ,NX,NY) COMMON /D2V/D2VM(NY-1),D2V0(NY-1),D2VP(NY-1) COMMON /D2P/D2PM(2:NY),D2P0(NY),D2PP(NY-1) COMMON /IP/IP(NX) /IM/IM(NX) /CMPX/CM1X,C00X,CP1X COMMON /KP/KP(NZ) /KM/KM(NZ) /CMPZ/CM1Z,C00Z,CP1Z DO 10 J=1,NY C00=C00X+C00Z+D2P0(J) IF(J.EQ.1) THEN DO 21 I=1,NX DO 21 K=1,NZ UF(K,I,J)=UF(K,I,J) &+(CM1X*U(K,IM(I),J)+CM1Z*U(KM(K),I,J)+C00*U(K,I,J) & +CP1Z*U(KP(K),I,J)+CP1X*U(K,IP(I),J)+D2PP(J)*U(K,I,J+1))/RET WF(K,I,J)=WF(K,I,J) &+(CM1X*W(K,IM(I),J)+CM1Z*W(KM(K),I,J)+C00*W(K,I,J) & +CP1Z*W(KP(K),I,J)+CP1X*W(K,IP(I),J)+D2PP(J)*W(K,I,J+1))/RET 21 CONTINUE ELSE IF(J.LT.NY) THEN DO 20 I=1,NX DO 20 K=1,NZ UF(K,I,J)=UF(K,I,J) &+(D2PM(J)*U(K,I,J-1)+CM1X*U(K,IM(I),J)+CM1Z*U(KM(K),I,J) & +C00*U(K,I,J) & +CP1Z*U(KP(K),I,J)+CP1X*U(K,IP(I),J)+D2PP(J)*U(K,I,J+1))/RET WF(K,I,J)=WF(K,I,J) &+(D2PM(J)*W(K,I,J-1)+CM1X*W(K,IM(I),J)+CM1Z*W(KM(K),I,J) & +C00*W(K,I,J) & +CP1Z*W(KP(K),I,J)+CP1X*W(K,IP(I),J)+D2PP(J)*W(K,I,J+1))/RET 20 CONTINUE ELSE DO 22 I=1,NX DO 22 K=1,NZ UF(K,I,J)=UF(K,I,J) &+(D2PM(J)*U(K,I,J-1)+CM1X*U(K,IM(I),J)+CM1Z*U(KM(K),I,J) & +C00*U(K,I,J)+CP1Z*U(KP(K),I,J)+CP1X*U(K,IP(I),J))/RET WF(K,I,J)=WF(K,I,J) &+(D2PM(J)*W(K,I,J-1)+CM1X*W(K,IM(I),J)+CM1Z*W(KM(K),I,J) & +C00*W(K,I,J)+CP1Z*W(KP(K),I,J)+CP1X*W(K,IP(I),J))/RET 22 CONTINUE ENDIF ENDIF 10 CONTINUE DO 30 J=1,NY-1 C00=C00X+C00Z+D2V0(J) DO 30 I=1,NX DO 30 K=1,NZ VF(K,I,J)=VF(K,I,J) &+(D2VM(J)*V(K,I,J-1)+CM1X*V(K,IM(I),J)+CM1Z*V(KM(K),I,J) & +C00*V(K,I,J) & +CP1Z*V(KP(K),I,J)+CP1X*V(K,IP(I),J)+D2VP(J)*V(K,I,J+1) & )/RET 30 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 /DT/DT /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) 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) 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 /D/DX,DZ COMMON /D1P/D1PM(NY),D1PP(NY) /IM/IM(NX) /KM/KM(NZ) DO 10 J=1,NY DO 10 I=1,NX DO 10 K=1,NZ Q(K,I,J)=((-U(K,IM(I),J)+U(K,I,J))/DX & +D1PM(J)*V(K,I,J-1)+D1PP(J)*V(K,I,J) & +(-W(KM(K),I,J)+W(K,I,J))/DZ )/DT 10 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) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /DPP/DPPM(2:NY),DPP0(NY),DPPP(NY-1) COMMON /P/P(NZ,NX,NY) /W3/Q(NZ,NX,NY) /ERRP/POIERR /ITRP/ITRP COMMON /IP/IP(NX) /IM/IM(NX) /CMPX/CM1X,C00X,CP1X COMMON /KP/KP(NZ) /KM/KM(NZ) /CMPZ/CM1Z,C00Z,CP1Z 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 20 J=1,NY C00=C00X+C00Z+DPP0(J) DPC=-1.5D0/C00 IF(J.EQ.1) THEN DO 31 I=1,NX DO 31 K=1,NZ RESI=CM1X*P(K,IM(I),J)+CM1Z*P(KM(K),I,J)+C00*P(K,I,J) & +CP1Z*P(KP(K),I,J)+CP1X*P(K,IP(I),J)+DPPP(J)*P(K,I,J+1) & -Q(K,I,J) P(K,I,J)=P(K,I,J)+RESI*DPC SUMR=SUMR+RESI**2 31 CONTINUE ELSE IF(J.LT.NY) THEN DO 30 I=1,NX DO 30 K=1,NZ RESI=DPPM(J)*P(K,I,J-1)+CM1X*P(K,IM(I),J)+CM1Z*P(KM(K),I,J) & +C00*P(K,I,J) & +CP1Z*P(KP(K),I,J)+CP1X*P(K,IP(I),J)+DPPP(J)*P(K,I,J+1) & -Q(K,I,J) P(K,I,J)=P(K,I,J)+RESI*DPC SUMR=SUMR+RESI**2 30 CONTINUE ELSE DO 32 I=1,NX DO 32 K=1,NZ RESI=DPPM(J)*P(K,I,J-1)+CM1X*P(K,IM(I),J)+CM1Z*P(KM(K),I,J) & +C00*P(K,I,J)+CP1Z*P(KP(K),I,J)+CP1X*P(K,IP(I),J) & -Q(K,I,J) P(K,I,J)=P(K,I,J)+RESI*DPC SUMR=SUMR+RESI**2 32 CONTINUE ENDIF ENDIF 20 CONTINUE POIERR=DSQRT(SUMR/SUMS) C WRITE(6,*) 'ITRP=',ITRP,' ERR=',POIERR 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) 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) /DT/DT /D/DX,DZ /DY/DY(NY) COMMON /IP/IP(NX) /KP/KP(NZ) /D1VN/D1VNM(0:NY),D1VNP(0:NY) 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)-TX*(-P(K,I,J)+P(K,IP(I),J)) W(K,I,J)=W(K,I,J)-TZ*(-P(K,I,J)+P(KP(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)-DT*(D1VNM(J)*P(K,I,J)+D1VNP(J)*P(K,I,J+1)) 20 CONTINUE RETURN END C ********************************************************************* C * SBR. UMR : MEAN & RMS VALUES * C ********************************************************************* SUBROUTINE SBRUMR PARAMETER (NX=64,NY=64,NZ=64,NY1=NY+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) /S12/S12L(0:NY),S12V(0:NY) COMMON /IP/IP(NX) /IM/IM(NX) /D0V/D0VM(NY-1),D0VP(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) 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 SUMUV=SUMUV+(D0VM(J)*U(K,I,J)+D0VP(J)*U(K,I,J+1)) & *(V(K,I,J)+V(K,IP(I),J)) 40 CONTINUE VM(J)=ARXZ*SUMV1 BVV2(J)=ARXZ*SUMV2 S12L(J)=ARXZ*SUMUV/2.D0 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) 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) 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/D0PM(NY),D0PP(NY) /D1P/D1PM(NY),D1PP(NY) COMMON /UR/UR(NY),VR(0:NY),WR(NY) COMMON /IM/IM(NX) /KM/KM(NZ) 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=(-U(K,IM(I),J)+U(K,I,J))/DX & +D1PM(J)*V(K,I,J-1)+D1PP(J)*V(K,I,J) & +(-W(KM(K),I,J)+W(K,I,J))/DZ UCP=(U(K,IM(I),J)+U(K,I,J))/2.D0 WCP=(W(KM(K),I,J)+W(K,I,J))/2.D0 VCP=D0PM(J)*V(K,I,J-1)+D0PP(J)*V(K,I,J) 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