1 2 PROGRAM CFDF 3 4 PARAMETER ( N = 50, niters = 10 ) 5 6 REAL VX(N,N,N), VY(N,N,N), VZ(N,N,N) 7 REAL VX2(N,N,N), VY2(N,N,N), VZ2(N,N,N) 8 REAL P(N,N,N) 9 REAL FX(N,N,N), FY(N,N,N), FZ(N,N,N) 10 REAL AX(N,N,N), AY(N,N,N), AZ(N,N,N) 11 12 INTEGER iter 13 14C Initialize arrays 15 16 CALL initialize(N, VX, VY, VZ, VX2, VY2, VZ2, P, FX, FY, FZ, 17 . AX, AY, AZ) 18 19C Apply the stencil a few times 20 21 DO iter=1,niters 22 CALL cfdStencil(N, VX, VY, VZ, VX2, VY2, VZ2, P, 23 . FX, FY, FZ, AX, AY, AZ) 24 END DO 25 26 STOP 27 END 28 29 SUBROUTINE cfdStencil(N, VX, VY, VZ, VX2, VY2, VZ2, P, 30 . FX, FY, FZ, AX, AY, AZ) 31 32 INTEGER N 33 REAL VX(N,N,N), VY(N,N,N), VZ(N,N,N) 34 REAL VX2(N,N,N), VY2(N,N,N), VZ2(N,N,N) 35 REAL P(N,N,N) 36 REAL FX(N,N,N), FY(N,N,N), FZ(N,N,N) 37 REAL AX(N,N,N), AY(N,N,N), AZ(N,N,N) 38 39 PARAMETER ( delta_t = 0.001, recip_rho = 1.0e-3, 40 . eta = 1.0e-6, c1 = 0.1, c2 = 0.1 ) 41 42 DO i=3,N-2 43 DO j=3,N-2 44 DO k=3,N-2 45 VX2(i,j,k) = VX(i,j,k)+delta_t*(recip_rho*(eta* 46 . c1 * (-90*VX(i,j,k)-VX(i-2,j,k)+16*VX(i-1,j,k) 47 . +16*VX(i+1,j,k)-VX(i+2,j,k)-VX(i,j-2,k)+16*VX(i,j-1,k) 48 . +16*VX(i,j+1,k)-VX(i,j+2,k)-VX(i,j,k-2)+16*VX(i,j,k-1) 49 . +16*VX(i,j,k+1)-VX(i,j,k+2))+c2*(P(i-2,j,k) 50 . -8*P(i-1,j,k)+8*P(i+1,j,k)+P(i+2,j,k))+FX(i,j,k)) 51 . -AX(i,j,k)) 52 VY2(i,j,k) = VY(i,j,k)+delta_t*(recip_rho*(eta* 53 . c1 * (-90*VY(i,j,k)-VY(i-2,j,k)+16*VY(i-1,j,k) 54 . +16*VY(i+1,j,k)-VY(i+2,j,k)-VY(i,j-2,k)+16*VY(i,j-1,k) 55 . +16*VY(i,j+1,k)-VY(i,j+2,k)-VY(i,j,k-2)+16*VY(i,j,k-1) 56 . +16*VY(i,j,k+1)-VY(i,j,k+2))+c2*(P(i,j-2,k) 57 . -8*P(i,j-1,k)+8*P(i,j+1,k)+P(i,j+2,k))+FY(i,j,k)) 58 . -AY(i,j,k)) 59 VZ2(i,j,k) = VZ(i,j,k)+delta_t*(recip_rho*(eta* 60 . c1 * (-90*VZ(i,j,k)-VZ(i-2,j,k)+16*VZ(i-1,j,k) 61 . +16*VZ(i+1,j,k)-VZ(i+2,j,k)-VZ(i,j-2,k)+16*VZ(i,j-1,k) 62 . +16*VZ(i,j+1,k)-VZ(i,j+2,k)-VZ(i,j,k-2)+16*VZ(i,j,k-1) 63 . +16*VZ(i,j,k+1)-VZ(i,j,k+2))+c2*(P(i,j,k-2) 64 . -8*P(i,j,k-1)+8*P(i,j,k+1)+P(i,j,k+2))+FZ(i,j,k)) 65 . -AZ(i,j,k)) 66 END DO 67 END DO 68 END DO 69 70 RETURN 71 END 72 73 74 75 SUBROUTINE initialize(N, VX, VY, VZ, VX2, VY2, VZ2, P, 76 . FX, FY, FZ, AX, AY, AZ) 77 78 INTEGER N 79 REAL VX(N,N,N), VY(N,N,N), VZ(N,N,N) 80 REAL VX2(N,N,N), VY2(N,N,N), VZ2(N,N,N) 81 REAL P(N,N,N) 82 REAL FX(N,N,N), FY(N,N,N), FZ(N,N,N) 83 REAL AX(N,N,N), AY(N,N,N), AZ(N,N,N) 84 85 DO i=3,N-2 86 DO j=3,N-2 87 DO k=3,N-2 88 VX(i,j,k) = 0 89 VY(i,j,k) = 0 90 VZ(i,j,k) = 0 91 P(i,j,k) = 0 92 FX(i,j,k) = 0 93 FY(i,j,k) = 0 94 FZ(i,j,k) = 0 95 AX(i,j,k) = 0 96 AY(i,j,k) = 0 97 AZ(i,j,k) = 0 98 END DO 99 END DO 100 END DO 101 102 END 103 104