Pess
25-06-2011, 15:44
Ho un problema nell'esecuzione di un programma Fortran che serve a "decomprimere" ed ordinare Database sulla turbolenza (lo potete trovare qua (http://mp0806.cineca.it/icfd.php?uuid=3f9d4d7eac6746a18e832b99de119d6f) insieme ai database citati); tali file decompressi mi serviranno per una esercitazione per l'università, come da regola di sezione dunque questa è una "fase preliminare" per ottenere i dati su cui lavorare e non il lavoro da svolgere...
Il programmino in questione è stato scritto in Fortran dal ricercatore che ha raccolto i database ed in teoria basterebbe eseguirlo (mettendo al rigo 72 il nome del file da decomprimere) ma è da 2 giorni interi che ci provo senza successo.
Premetto che ho usato vari compilatori su varie macchine e quello con cui mi trovo meglio è FTN95 (http://www.silverfrost.com/11/ftn95/ftn95_fortran_95_for_windows.aspx) (è free per uso privato) che mi ha permesso di fare un debug ed individuare l'errore:
"Runtime error from program:c:\dmt\postbl-uvwp3d-af.exe
Run-time Error
*** Error 158, Unformatted record is too short for input list
main - in file postbl-uvwp3d-af.f at line 74 [+00d4]"
Come si può correggere? Sono assolutamente a digiuno di Fortran e non so dove mettere le mani! Ho provato con molto altri compilatori ed IDE (Force 2.0, gfortran...) ma il risultato è sempre un errore di Runtime... Ho fatto prove su due diversi PC di qualche annetto fa con Windows XP e sul notebook in firma...
Ringrazio in anticipo chiunque vorrà aiutarmi...
EDIT: Riporto il codice per intero con la riga che dà l'errore evidenziata in rosso; l'errore ovviamente vale anche per gli altri punti in cui compare questa fantomatica T0 come parametro di lista per i comandi READ e WRITE...
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C !
C ******************************************************* !
C * POSTPROCESSOR: * !
C * SPATIALLY DEVELOPING TURBULENT BOUNDARY LAYER * !
C ******************************************************* !
C !
C AUTHOR: A.FERRANTE !
C DATE: 21 AUGUST, 2006 !
C INPUT: uvwp3d.ibm.# (Unformatted IBM-P4 file) !
C OUTPUT: field2dy.plt field2dx.pld (Formatted Tecplot files) !
C !
C CODE DESCRIPTION: !
C READS INSTANTANEOUS 3D FLUID VELOCITY COMPONENTS U,V,W !
C AND PRESSURE P FROM uvwp3d.ibm.# !
C WRITES X,Z,U,V,W,P ON field2dy.plt !
C WRITES Y,Z,U,V,W,P ON field2dx.plt !
C !
C COMPILING LINE ON IBM-P4+: !
C xlf90_r -qfixed -qrealsize=8 -qsuppress=cmpmsg !
C !
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
MODULE CPARA
INTEGER, PARAMETER :: IMP=514,JMP=258,KMP=98
INTEGER, PARAMETER :: ISLICE=484, JSLICE=129
REAL, PARAMETER :: LX=20.0, LY=5.0, LZ=3.6
INTEGER :: NCUBE,NXZ,NT,IM,JM,KM
REAL :: DX,DY,DZ
END MODULE CPARA
!
MODULE CBOLA
REAL, PARAMETER :: REDEL=8000., ATL=0.66
REAL, DIMENSION(:), ALLOCATABLE :: ZETA,ZI,FZZ
REAL, PARAMETER :: VKC=0.41,BLL=5.2,PIWL=0.5
END MODULE CBOLA
!
MODULE CFLD
REAL, DIMENSION(:), ALLOCATABLE ::
. UFLD,VFLD,WFLD,DUMMY,PRESS
REAL, DIMENSION(:,:,:), ALLOCATABLE :: U, V, W
END MODULE CFLD
!
!----------------------------------------------------------------------
PROGRAM POSTPRO
C
C I= 1,IM+2 ARRAYS DIMENSION -- (IMP=IM+2) SIMILAR FOR J AND K
C I= 2,IM+1 PHYSICAL GRID POINTS IM
C DX=LX/IM, DY=LY/JM, DZ=LZ/KM
C
C U(I,J,K) u component of velocity in (i+.5,j ,k)
C V v component of velocity in (i ,j+.5,k)
C W w component of velocity in (i ,j ,k+.5)
C
C PRESS(I,J,K) pressure in (i,j,k)
C
C UFLD(:) u component of velocity in (i,j,k)
C VFLD v //
C WFLD w //
C
C X,Y,Z (I ,J ,K ) [-DX/2:LX+DX/2] {1,IMP}
C XS,YS,ZS (I+.5,J+.5,K+.5) [ 0:LX+DX ] {1,IMP}
C
USE cpara
USE cfld
USE cbola
C
CALL OP_WR
C
DO 1 IXYZ=1,2 ! X and Y planes
C
NGPH1=25
OPEN(NGPH1,FILE='uvwp3d.ibm.06160',FORM='UNFORMATTED')
REWIND(NGPH1)
READ(NGPH1) T0
BACKSPACE(NGPH1)
C
CALL ALLOC_VAR(IXYZ)
C
IM = IMP-2
JM = JMP-2
KM = KMP-2
DX = LX/FLOAT(IM)
DY = LY/FLOAT(JM)
DZ = LZ/FLOAT(KM)
C
CALL BLUREF
C
NT=0
C
10 CONTINUE ! LOOP ON TIME NT
NT=NT+1
C
READ(NGPH1,END=100) T0
BACKSPACE(NGPH1)
C
CALL READU(IXYZ)
print *,'nt=',NT
C
C ***AVERAGE VELOCITY FIELD IN GRID POINTS I,J,K***
C
CALL COPYFIELD(IXYZ)
CALL AVGRPO(IXYZ)
C
C ***WRITE IN OUTPUT VELOCITY AND PRESSURE FIELD***
C
CALL WRSLICE(IXYZ)
C
GOTO 10
C
100 CONTINUE
C
CALL DEALLOC_VAR
CLOSE(NGPH1)
C
1 CONTINUE ! IPLANE=1,2 FOR X AND Y
C
WRITE(*,*)'I-SLICE X=',DX*(ISLICE-1.5)
WRITE(*,*)'J-SLICE Y=',DY*(JSLICE-1.5)
STOP
END
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C S U B R O U T I N E O P _ W R
C
C PURPOSE:
C --------
C THIS SUBROUTINE OPENS INPUT FILES AND OUTPUT FILES
C AND WRITES HEADERS OF OUTPUT FILES.
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE OP_WR
C
USE cpara
USE cfld
USE cbola
c
OPEN( 9 ,FILE='field2dx.plt')
OPEN(11 ,FILE='field2dy.plt')
c OPEN(13 ,FILE='field2dz.plt')
WRITE( 9,'(35HVARIABLES = "Y" "Z" "U" "V" "W" "P")')
WRITE(11,'(35HVARIABLES = "X" "Z" "U" "V" "W" "P")')
c WRITE(13,'(35HVARIABLES = "X" "Y" "U" "V" "W" "P")')
RETURN
END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C S U B R O U T I N E A L L O C _ V A R
C
C PURPOSE:
C --------
C THIS SUBROUTINE ALLOCATES ALL THE ARRAYS
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE ALLOC_VAR(IXYZ)
!
USE cpara
USE cfld
USE cbola
!
IF(IXYZ==1) THEN
NX = 4 ! AF-5MAR02
NY = JMP
NZ = KMP
NCUBE = NX*NY*NZ
ALLOCATE (DUMMY(IMP*KMP))
ALLOCATE (ZETA(KMP),ZI(KMP),FZZ(KMP))
ALLOCATE (UFLD(NCUBE),VFLD(NCUBE),WFLD(NCUBE)
. ,PRESS(NCUBE))
ALLOCATE (U(1:4,0:NY,0:NZ)
. ,V(1:4,0:NY,0:NZ)
. ,W(1:4,0:NY,0:NZ))
ELSEIF(IXYZ==2) THEN
NX = IMP
NY = 4
NZ = KMP
NCUBE = NX*NY*NZ
ALLOCATE (DUMMY(IMP*KMP))
ALLOCATE (ZETA(KMP),ZI(KMP),FZZ(KMP))
ALLOCATE (UFLD(NCUBE),VFLD(NCUBE),WFLD(NCUBE)
. ,PRESS(NCUBE))
ALLOCATE (U(0:NX,1:4,0:NZ)
. ,V(0:NX,1:4,0:NZ)
. ,W(0:NX,1:4,0:NZ))
ENDIF
RETURN
END
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C S U B R O U T I N E D E A L L O C _ V A R
C
C PURPOSE:
C --------
C THIS SUBROUTINE ALLOCATES ALL THE ARRAYS
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE DEALLOC_VAR
!
USE cpara
USE cfld
USE cbola
!
DEALLOCATE (DUMMY)
DEALLOCATE (ZETA,ZI,FZZ)
DEALLOCATE (UFLD,VFLD,WFLD,PRESS)
DEALLOCATE (U,V,W)
RETURN
END
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE BLUREF
C
C COMPUTES THE STRETCHED MESH FUNCTION: ZETA(K)
C REF. FERRANTE & ELGHOBASHI, JCP 2004
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
USE CPARA
USE CFLD
USE CBOLA
C
C COMPUTE REFERENCE VELOCITY PROFILE u+=f_w(z+)
C
TGLZ =TANH(ATL*LZ)
QTGLZ=1./TGLZ
QTLZ =1./(ATL*LZ)
DO K=1,KM+1
ZI(K)=(K-1.5)*DZ
ZETA(K)=LZ*(1.-TANH(ATL*(LZ-ZI(K)))*QTGLZ)
FZZ(K)=TGLZ*QTLZ*(COSH(ATL*(LZ-ZI(K))))**2
ENDDO
RETURN
END
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C S U B R O U T I N E R E A D I N
C
C PURPOSE:
C --------
C THIS SUBROUTINE READS IN VARIOUS VARIABLES
C FOR THE GRAPHICS PROGRAM TO USE.
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE READU(IXYZ)
!
USE cpara
USE cfld
!
IK(I,K) = I + (K-1)*IMP
IJK(I,J,K)= I-ISLICE+2 + (K-1)*4 + (J-1)*4*KMP
C
print*,'in readu'
C
NGPH1=25
IEOF=0
IF (IXYZ==1) THEN ! X PLANE (YZ)
READ(NGPH1,END=100) T0
WRITE(*,*) 'T=',T0
NXZ = IMP*KMP
C
DO 20 J=1,JMP
READ(NGPH1,END=100) (DUMMY(NPCNT),NPCNT=1,NXZ)
DO 21 K=1,KMP
DO 21 I=ISLICE-1,ISLICE+2
IJKL=IJK(I,J,K)
IKL=IK(I,K)
UFLD(IJKL)=DUMMY(IKL)
21 CONTINUE
20 IEOF = IEOF+1
WRITE(*,*) 'READ U'
c
DO 30 J=1,JMP
READ(NGPH1,END=100) (DUMMY(NPCNT),NPCNT=1,NXZ)
DO 31 K=1,KMP
DO 31 I=ISLICE-1,ISLICE+2
IJKL=IJK(I,J,K)
IKL=IK(I,K)
VFLD(IJKL)=DUMMY(IKL)
31 CONTINUE
30 IEOF = IEOF+1
WRITE(*,*) 'READ V'
c
DO 40 J=1,JMP
READ(NGPH1,END=100) (DUMMY(NPCNT),NPCNT=1,NXZ)
DO 41 K=1,KMP
DO 41 I=ISLICE-1,ISLICE+2
IJKL=IJK(I,J,K)
IKL=IK(I,K)
WFLD(IJKL)=DUMMY(IKL)
41 CONTINUE
40 IEOF = IEOF+1
WRITE(*,*) 'READ W'
c
DO 42 J=1,JMP
READ(NGPH1,END=100) (DUMMY(NPCNT),NPCNT=1,NXZ)
DO 43 K=1,KMP
DO 43 I=ISLICE-1,ISLICE+2
IJKL=IJK(I,J,K)
IKL=IK(I,K)
PRESS(IJKL)=DUMMY(IKL)
43 CONTINUE
42 IEOF = IEOF+1
WRITE(*,*) 'READ P'
c
ELSEIF(IXYZ==2) THEN ! Y PLANE (XZ)
C
READ(NGPH1,END=100) T0
WRITE(*,*) 'T=',T0
NXZ = IMP*KMP
C
DO 50 J=1,JMP
IF((J==JSLICE-1).OR.(J==JSLICE).OR.(J==JSLICE+1).OR.(J==JSLICE+2))
. THEN
READ(NGPH1,END=100) (UFLD(NPCNT+NXZ*(J-JSLICE+1)) ,NPCNT=1,NXZ)
ELSE
READ(NGPH1,END=100) (DUMMY(NPCNT),NPCNT=1,NXZ)
ENDIF
50 IEOF = IEOF+1
WRITE(*,*) 'READ U'
c
DO 60 J=1,JMP
IF((J==JSLICE-1).OR.(J==JSLICE).OR.(J==JSLICE+1).OR.(J==JSLICE+2))
. THEN
READ(NGPH1,END=100) (VFLD(NPCNT+NXZ*(J-JSLICE+1)) ,NPCNT=1,NXZ)
ELSE
READ(NGPH1,END=100) (DUMMY(NPCNT),NPCNT=1,NXZ)
ENDIF
60 IEOF = IEOF+1
WRITE(*,*) 'READ V'
c
DO 70 J=1,JMP
IF((J==JSLICE-1).OR.(J==JSLICE).OR.(J==JSLICE+1).OR.(J==JSLICE+2))
. THEN
READ(NGPH1,END=100) (WFLD(NPCNT+NXZ*(J-JSLICE+1)) ,NPCNT=1,NXZ)
ELSE
READ(NGPH1,END=100) (DUMMY(NPCNT),NPCNT=1,NXZ)
ENDIF
70 IEOF = IEOF+1
WRITE(*,*) 'READ W'
C
DO 75 J=1,JMP
IF((J==JSLICE-1).OR.(J==JSLICE).OR.(J==JSLICE+1).OR.(J==JSLICE+2))
. THEN
READ(NGPH1,END=100) (PRESS(NPCNT+NXZ*(J-JSLICE+1)) ,NPCNT=1,NXZ)
ELSE
READ(NGPH1,END=100) (DUMMY(NPCNT),NPCNT=1,NXZ)
ENDIF
75 IEOF = IEOF+1
WRITE(*,*) 'READ P'
C
ENDIF !IF(IXYZ==1,2,3)
C
100 CONTINUE
C
print*,'IEOF=',IEOF
print*,'out readu'
C
RETURN
END
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C S U B R O U T I N E C O P Y F I E L D
C
C PURPOSE:
C --------
C THIS SUBROUTINE COPIES THE VELOCITY UFLD TO U,V,W(I,J,K)
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE COPYFIELD(IXYZ)
C
USE cpara
USE cfld
C ***SPACE LOOPS TO COPY THE FIELD***
IF(IXYZ==1) THEN
C
DO 100 K = 1,KMP
DO 100 J = 1,JMP
DO 100 I = 1,4
NPCNT = I + 4*KMP*(J-1) + 4*(K-1)
U(I,J,K) = UFLD(NPCNT)
V(I,J,K) = VFLD(NPCNT)
W(I,J,K) = WFLD(NPCNT)
100 CONTINUE
C
ELSEIF (IXYZ==2) THEN
C ***SPACE LOOPS TO COPY THE FIELD***
DO 200 K = 1,KMP
DO 200 J = 1,4
DO 200 I = 1,IMP
NPCNT = I + NXZ*(J-1) + IMP*(K-1)
U(I,J,K) = UFLD(NPCNT)
V(I,J,K) = VFLD(NPCNT)
W(I,J,K) = WFLD(NPCNT)
200 CONTINUE
C
ENDIF
C
RETURN
END
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C S U B R O U T I N E A V G R P O
C
C PURPOSE:
C --------
C THIS SUBROUTINE AVERAGES THE VELOCITY COMPONENTS
C IN THE PRESSURE GRID POINTS FROM THE STAGGERED VEL.
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE AVGRPO(IXYZ)
C
USE cpara
USE cfld
c
IF(IXYZ==1) THEN
C ***AVERAGE VELOCITY COMPONENTS***
DO 100 K = 2,KMP
DO 100 J = 2,JMP
DO 100 I = 2,4
NPCNT = I + 4*KMP*(J-1) + 4*(K-1)
UFLD(NPCNT) = .5 * (U(I-1,J,K) + U(I,J,K))
VFLD(NPCNT) = .5 * (V(I,J-1,K) + V(I,J,K))
WFLD(NPCNT) = .5 * (W(I,J,K-1) + W(I,J,K))
c
IF (J.EQ.JM+1) THEN ! PERIODIC BOUNDARY CONDITION IN J
NP1 = I + 4*(K-1)
UFLD(NP1) = UFLD(NPCNT)
VFLD(NP1) = VFLD(NPCNT)
WFLD(NP1) = WFLD(NPCNT)
ENDIF
IF (K.EQ.2) THEN ! NO SLIP BC AT THE WALL U(1)=0 FOR K=1
NP1 = I + NXZ*(J-1)
UFLD(NP1) = 0.
VFLD(NP1) = 0.
WFLD(NP1) = 0.
ENDIF
100 CONTINUE
C
ELSEIF(IXYZ==2) THEN
C ***AVERAGE VELOCITY COMPONENTS***
DO 200 K = 2,KMP
DO 200 J = 2,4
DO 200 I = 2,IMP
NPCNT = I + NXZ*(J-1) + IMP*(K-1)
UFLD(NPCNT) = .5 * (U(I-1,J,K) + U(I,J,K))
VFLD(NPCNT) = .5 * (V(I,J-1,K) + V(I,J,K))
WFLD(NPCNT) = .5 * (W(I,J,K-1) + W(I,J,K))
IF (K.EQ.2) THEN ! NO SLIP BC AT THE WALL U(1)=0 FOR K=1
NP1 = I + NXZ*(J-1)
UFLD(NP1) = 0.
VFLD(NP1) = 0.
WFLD(NP1) = 0.
ENDIF
C
200 CONTINUE
C
ENDIF
RETURN
END
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C S U B R O U T I N E W R S L I C E
C
C PURPOSE:
C --------
C THIS SUBROUTINE WRITES THE VELOCITY
C AND PRESSURE FIELD AT THE GRID POINTS ON X OR Y PLANE.
C THE VELOCITY FIELD IS AVERAGED AT THE GRID POINTS.
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE WRSLICE(IXYZ)
C
USE cpara
USE cfld
USE cbola
C ***WRITING OUTPUT***
print*,'writing velocity slice...'
C ***TIME LOOP***
IF(IXYZ==1) THEN
I=2
C
WRITE(9,200) T0,JM,KM
C ***SPACE LOOPS***
DO 90 K = 1,KM ! AF-28SEP02 AFTER REMOVING THE AVFIELD AND LEAVING AVGRPO
Z = ZETA(K)
IF(K==1) Z=0.0
DO 90 J = 2,JM+1
Y = DY*(J-1.5)
NPCNT = I + 4*KMP*(J-1) + 4*(K-1)
WRITE(9,300) Y,Z
. ,UFLD(NPCNT),VFLD(NPCNT),WFLD(NPCNT)
. ,PRESS(NPCNT)
90 CONTINUE
C
ELSEIF(IXYZ==2) THEN
C
J=2
WRITE(11,200) T0,IM,KM
C ***SPACE LOOPS***
DO 100 K = 1,KM
Z=ZETA(K)
IF(K==1) Z=0.0
DO 100 I = 2,IM+1
X = DX*(I-1.5)
NPCNT = I + NXZ*(J-1) + IMP*(K-1)
WRITE(11,300) X,Z
. ,UFLD(NPCNT),VFLD(NPCNT),WFLD(NPCNT)
. ,PRESS(NPCNT)
100 CONTINUE
C
ENDIF
C
200 FORMAT('ZONE T="T =',G14.5,'"'
. ,' F=POINT, I=',I4,' J=',I4)
300 FORMAT( 1X,6G14.5)
C
RETURN
END
C================================================================================
Il programmino in questione è stato scritto in Fortran dal ricercatore che ha raccolto i database ed in teoria basterebbe eseguirlo (mettendo al rigo 72 il nome del file da decomprimere) ma è da 2 giorni interi che ci provo senza successo.
Premetto che ho usato vari compilatori su varie macchine e quello con cui mi trovo meglio è FTN95 (http://www.silverfrost.com/11/ftn95/ftn95_fortran_95_for_windows.aspx) (è free per uso privato) che mi ha permesso di fare un debug ed individuare l'errore:
"Runtime error from program:c:\dmt\postbl-uvwp3d-af.exe
Run-time Error
*** Error 158, Unformatted record is too short for input list
main - in file postbl-uvwp3d-af.f at line 74 [+00d4]"
Come si può correggere? Sono assolutamente a digiuno di Fortran e non so dove mettere le mani! Ho provato con molto altri compilatori ed IDE (Force 2.0, gfortran...) ma il risultato è sempre un errore di Runtime... Ho fatto prove su due diversi PC di qualche annetto fa con Windows XP e sul notebook in firma...
Ringrazio in anticipo chiunque vorrà aiutarmi...
EDIT: Riporto il codice per intero con la riga che dà l'errore evidenziata in rosso; l'errore ovviamente vale anche per gli altri punti in cui compare questa fantomatica T0 come parametro di lista per i comandi READ e WRITE...
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C !
C ******************************************************* !
C * POSTPROCESSOR: * !
C * SPATIALLY DEVELOPING TURBULENT BOUNDARY LAYER * !
C ******************************************************* !
C !
C AUTHOR: A.FERRANTE !
C DATE: 21 AUGUST, 2006 !
C INPUT: uvwp3d.ibm.# (Unformatted IBM-P4 file) !
C OUTPUT: field2dy.plt field2dx.pld (Formatted Tecplot files) !
C !
C CODE DESCRIPTION: !
C READS INSTANTANEOUS 3D FLUID VELOCITY COMPONENTS U,V,W !
C AND PRESSURE P FROM uvwp3d.ibm.# !
C WRITES X,Z,U,V,W,P ON field2dy.plt !
C WRITES Y,Z,U,V,W,P ON field2dx.plt !
C !
C COMPILING LINE ON IBM-P4+: !
C xlf90_r -qfixed -qrealsize=8 -qsuppress=cmpmsg !
C !
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
MODULE CPARA
INTEGER, PARAMETER :: IMP=514,JMP=258,KMP=98
INTEGER, PARAMETER :: ISLICE=484, JSLICE=129
REAL, PARAMETER :: LX=20.0, LY=5.0, LZ=3.6
INTEGER :: NCUBE,NXZ,NT,IM,JM,KM
REAL :: DX,DY,DZ
END MODULE CPARA
!
MODULE CBOLA
REAL, PARAMETER :: REDEL=8000., ATL=0.66
REAL, DIMENSION(:), ALLOCATABLE :: ZETA,ZI,FZZ
REAL, PARAMETER :: VKC=0.41,BLL=5.2,PIWL=0.5
END MODULE CBOLA
!
MODULE CFLD
REAL, DIMENSION(:), ALLOCATABLE ::
. UFLD,VFLD,WFLD,DUMMY,PRESS
REAL, DIMENSION(:,:,:), ALLOCATABLE :: U, V, W
END MODULE CFLD
!
!----------------------------------------------------------------------
PROGRAM POSTPRO
C
C I= 1,IM+2 ARRAYS DIMENSION -- (IMP=IM+2) SIMILAR FOR J AND K
C I= 2,IM+1 PHYSICAL GRID POINTS IM
C DX=LX/IM, DY=LY/JM, DZ=LZ/KM
C
C U(I,J,K) u component of velocity in (i+.5,j ,k)
C V v component of velocity in (i ,j+.5,k)
C W w component of velocity in (i ,j ,k+.5)
C
C PRESS(I,J,K) pressure in (i,j,k)
C
C UFLD(:) u component of velocity in (i,j,k)
C VFLD v //
C WFLD w //
C
C X,Y,Z (I ,J ,K ) [-DX/2:LX+DX/2] {1,IMP}
C XS,YS,ZS (I+.5,J+.5,K+.5) [ 0:LX+DX ] {1,IMP}
C
USE cpara
USE cfld
USE cbola
C
CALL OP_WR
C
DO 1 IXYZ=1,2 ! X and Y planes
C
NGPH1=25
OPEN(NGPH1,FILE='uvwp3d.ibm.06160',FORM='UNFORMATTED')
REWIND(NGPH1)
READ(NGPH1) T0
BACKSPACE(NGPH1)
C
CALL ALLOC_VAR(IXYZ)
C
IM = IMP-2
JM = JMP-2
KM = KMP-2
DX = LX/FLOAT(IM)
DY = LY/FLOAT(JM)
DZ = LZ/FLOAT(KM)
C
CALL BLUREF
C
NT=0
C
10 CONTINUE ! LOOP ON TIME NT
NT=NT+1
C
READ(NGPH1,END=100) T0
BACKSPACE(NGPH1)
C
CALL READU(IXYZ)
print *,'nt=',NT
C
C ***AVERAGE VELOCITY FIELD IN GRID POINTS I,J,K***
C
CALL COPYFIELD(IXYZ)
CALL AVGRPO(IXYZ)
C
C ***WRITE IN OUTPUT VELOCITY AND PRESSURE FIELD***
C
CALL WRSLICE(IXYZ)
C
GOTO 10
C
100 CONTINUE
C
CALL DEALLOC_VAR
CLOSE(NGPH1)
C
1 CONTINUE ! IPLANE=1,2 FOR X AND Y
C
WRITE(*,*)'I-SLICE X=',DX*(ISLICE-1.5)
WRITE(*,*)'J-SLICE Y=',DY*(JSLICE-1.5)
STOP
END
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C S U B R O U T I N E O P _ W R
C
C PURPOSE:
C --------
C THIS SUBROUTINE OPENS INPUT FILES AND OUTPUT FILES
C AND WRITES HEADERS OF OUTPUT FILES.
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE OP_WR
C
USE cpara
USE cfld
USE cbola
c
OPEN( 9 ,FILE='field2dx.plt')
OPEN(11 ,FILE='field2dy.plt')
c OPEN(13 ,FILE='field2dz.plt')
WRITE( 9,'(35HVARIABLES = "Y" "Z" "U" "V" "W" "P")')
WRITE(11,'(35HVARIABLES = "X" "Z" "U" "V" "W" "P")')
c WRITE(13,'(35HVARIABLES = "X" "Y" "U" "V" "W" "P")')
RETURN
END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C S U B R O U T I N E A L L O C _ V A R
C
C PURPOSE:
C --------
C THIS SUBROUTINE ALLOCATES ALL THE ARRAYS
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE ALLOC_VAR(IXYZ)
!
USE cpara
USE cfld
USE cbola
!
IF(IXYZ==1) THEN
NX = 4 ! AF-5MAR02
NY = JMP
NZ = KMP
NCUBE = NX*NY*NZ
ALLOCATE (DUMMY(IMP*KMP))
ALLOCATE (ZETA(KMP),ZI(KMP),FZZ(KMP))
ALLOCATE (UFLD(NCUBE),VFLD(NCUBE),WFLD(NCUBE)
. ,PRESS(NCUBE))
ALLOCATE (U(1:4,0:NY,0:NZ)
. ,V(1:4,0:NY,0:NZ)
. ,W(1:4,0:NY,0:NZ))
ELSEIF(IXYZ==2) THEN
NX = IMP
NY = 4
NZ = KMP
NCUBE = NX*NY*NZ
ALLOCATE (DUMMY(IMP*KMP))
ALLOCATE (ZETA(KMP),ZI(KMP),FZZ(KMP))
ALLOCATE (UFLD(NCUBE),VFLD(NCUBE),WFLD(NCUBE)
. ,PRESS(NCUBE))
ALLOCATE (U(0:NX,1:4,0:NZ)
. ,V(0:NX,1:4,0:NZ)
. ,W(0:NX,1:4,0:NZ))
ENDIF
RETURN
END
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C S U B R O U T I N E D E A L L O C _ V A R
C
C PURPOSE:
C --------
C THIS SUBROUTINE ALLOCATES ALL THE ARRAYS
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE DEALLOC_VAR
!
USE cpara
USE cfld
USE cbola
!
DEALLOCATE (DUMMY)
DEALLOCATE (ZETA,ZI,FZZ)
DEALLOCATE (UFLD,VFLD,WFLD,PRESS)
DEALLOCATE (U,V,W)
RETURN
END
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE BLUREF
C
C COMPUTES THE STRETCHED MESH FUNCTION: ZETA(K)
C REF. FERRANTE & ELGHOBASHI, JCP 2004
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
USE CPARA
USE CFLD
USE CBOLA
C
C COMPUTE REFERENCE VELOCITY PROFILE u+=f_w(z+)
C
TGLZ =TANH(ATL*LZ)
QTGLZ=1./TGLZ
QTLZ =1./(ATL*LZ)
DO K=1,KM+1
ZI(K)=(K-1.5)*DZ
ZETA(K)=LZ*(1.-TANH(ATL*(LZ-ZI(K)))*QTGLZ)
FZZ(K)=TGLZ*QTLZ*(COSH(ATL*(LZ-ZI(K))))**2
ENDDO
RETURN
END
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C S U B R O U T I N E R E A D I N
C
C PURPOSE:
C --------
C THIS SUBROUTINE READS IN VARIOUS VARIABLES
C FOR THE GRAPHICS PROGRAM TO USE.
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE READU(IXYZ)
!
USE cpara
USE cfld
!
IK(I,K) = I + (K-1)*IMP
IJK(I,J,K)= I-ISLICE+2 + (K-1)*4 + (J-1)*4*KMP
C
print*,'in readu'
C
NGPH1=25
IEOF=0
IF (IXYZ==1) THEN ! X PLANE (YZ)
READ(NGPH1,END=100) T0
WRITE(*,*) 'T=',T0
NXZ = IMP*KMP
C
DO 20 J=1,JMP
READ(NGPH1,END=100) (DUMMY(NPCNT),NPCNT=1,NXZ)
DO 21 K=1,KMP
DO 21 I=ISLICE-1,ISLICE+2
IJKL=IJK(I,J,K)
IKL=IK(I,K)
UFLD(IJKL)=DUMMY(IKL)
21 CONTINUE
20 IEOF = IEOF+1
WRITE(*,*) 'READ U'
c
DO 30 J=1,JMP
READ(NGPH1,END=100) (DUMMY(NPCNT),NPCNT=1,NXZ)
DO 31 K=1,KMP
DO 31 I=ISLICE-1,ISLICE+2
IJKL=IJK(I,J,K)
IKL=IK(I,K)
VFLD(IJKL)=DUMMY(IKL)
31 CONTINUE
30 IEOF = IEOF+1
WRITE(*,*) 'READ V'
c
DO 40 J=1,JMP
READ(NGPH1,END=100) (DUMMY(NPCNT),NPCNT=1,NXZ)
DO 41 K=1,KMP
DO 41 I=ISLICE-1,ISLICE+2
IJKL=IJK(I,J,K)
IKL=IK(I,K)
WFLD(IJKL)=DUMMY(IKL)
41 CONTINUE
40 IEOF = IEOF+1
WRITE(*,*) 'READ W'
c
DO 42 J=1,JMP
READ(NGPH1,END=100) (DUMMY(NPCNT),NPCNT=1,NXZ)
DO 43 K=1,KMP
DO 43 I=ISLICE-1,ISLICE+2
IJKL=IJK(I,J,K)
IKL=IK(I,K)
PRESS(IJKL)=DUMMY(IKL)
43 CONTINUE
42 IEOF = IEOF+1
WRITE(*,*) 'READ P'
c
ELSEIF(IXYZ==2) THEN ! Y PLANE (XZ)
C
READ(NGPH1,END=100) T0
WRITE(*,*) 'T=',T0
NXZ = IMP*KMP
C
DO 50 J=1,JMP
IF((J==JSLICE-1).OR.(J==JSLICE).OR.(J==JSLICE+1).OR.(J==JSLICE+2))
. THEN
READ(NGPH1,END=100) (UFLD(NPCNT+NXZ*(J-JSLICE+1)) ,NPCNT=1,NXZ)
ELSE
READ(NGPH1,END=100) (DUMMY(NPCNT),NPCNT=1,NXZ)
ENDIF
50 IEOF = IEOF+1
WRITE(*,*) 'READ U'
c
DO 60 J=1,JMP
IF((J==JSLICE-1).OR.(J==JSLICE).OR.(J==JSLICE+1).OR.(J==JSLICE+2))
. THEN
READ(NGPH1,END=100) (VFLD(NPCNT+NXZ*(J-JSLICE+1)) ,NPCNT=1,NXZ)
ELSE
READ(NGPH1,END=100) (DUMMY(NPCNT),NPCNT=1,NXZ)
ENDIF
60 IEOF = IEOF+1
WRITE(*,*) 'READ V'
c
DO 70 J=1,JMP
IF((J==JSLICE-1).OR.(J==JSLICE).OR.(J==JSLICE+1).OR.(J==JSLICE+2))
. THEN
READ(NGPH1,END=100) (WFLD(NPCNT+NXZ*(J-JSLICE+1)) ,NPCNT=1,NXZ)
ELSE
READ(NGPH1,END=100) (DUMMY(NPCNT),NPCNT=1,NXZ)
ENDIF
70 IEOF = IEOF+1
WRITE(*,*) 'READ W'
C
DO 75 J=1,JMP
IF((J==JSLICE-1).OR.(J==JSLICE).OR.(J==JSLICE+1).OR.(J==JSLICE+2))
. THEN
READ(NGPH1,END=100) (PRESS(NPCNT+NXZ*(J-JSLICE+1)) ,NPCNT=1,NXZ)
ELSE
READ(NGPH1,END=100) (DUMMY(NPCNT),NPCNT=1,NXZ)
ENDIF
75 IEOF = IEOF+1
WRITE(*,*) 'READ P'
C
ENDIF !IF(IXYZ==1,2,3)
C
100 CONTINUE
C
print*,'IEOF=',IEOF
print*,'out readu'
C
RETURN
END
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C S U B R O U T I N E C O P Y F I E L D
C
C PURPOSE:
C --------
C THIS SUBROUTINE COPIES THE VELOCITY UFLD TO U,V,W(I,J,K)
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE COPYFIELD(IXYZ)
C
USE cpara
USE cfld
C ***SPACE LOOPS TO COPY THE FIELD***
IF(IXYZ==1) THEN
C
DO 100 K = 1,KMP
DO 100 J = 1,JMP
DO 100 I = 1,4
NPCNT = I + 4*KMP*(J-1) + 4*(K-1)
U(I,J,K) = UFLD(NPCNT)
V(I,J,K) = VFLD(NPCNT)
W(I,J,K) = WFLD(NPCNT)
100 CONTINUE
C
ELSEIF (IXYZ==2) THEN
C ***SPACE LOOPS TO COPY THE FIELD***
DO 200 K = 1,KMP
DO 200 J = 1,4
DO 200 I = 1,IMP
NPCNT = I + NXZ*(J-1) + IMP*(K-1)
U(I,J,K) = UFLD(NPCNT)
V(I,J,K) = VFLD(NPCNT)
W(I,J,K) = WFLD(NPCNT)
200 CONTINUE
C
ENDIF
C
RETURN
END
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C S U B R O U T I N E A V G R P O
C
C PURPOSE:
C --------
C THIS SUBROUTINE AVERAGES THE VELOCITY COMPONENTS
C IN THE PRESSURE GRID POINTS FROM THE STAGGERED VEL.
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE AVGRPO(IXYZ)
C
USE cpara
USE cfld
c
IF(IXYZ==1) THEN
C ***AVERAGE VELOCITY COMPONENTS***
DO 100 K = 2,KMP
DO 100 J = 2,JMP
DO 100 I = 2,4
NPCNT = I + 4*KMP*(J-1) + 4*(K-1)
UFLD(NPCNT) = .5 * (U(I-1,J,K) + U(I,J,K))
VFLD(NPCNT) = .5 * (V(I,J-1,K) + V(I,J,K))
WFLD(NPCNT) = .5 * (W(I,J,K-1) + W(I,J,K))
c
IF (J.EQ.JM+1) THEN ! PERIODIC BOUNDARY CONDITION IN J
NP1 = I + 4*(K-1)
UFLD(NP1) = UFLD(NPCNT)
VFLD(NP1) = VFLD(NPCNT)
WFLD(NP1) = WFLD(NPCNT)
ENDIF
IF (K.EQ.2) THEN ! NO SLIP BC AT THE WALL U(1)=0 FOR K=1
NP1 = I + NXZ*(J-1)
UFLD(NP1) = 0.
VFLD(NP1) = 0.
WFLD(NP1) = 0.
ENDIF
100 CONTINUE
C
ELSEIF(IXYZ==2) THEN
C ***AVERAGE VELOCITY COMPONENTS***
DO 200 K = 2,KMP
DO 200 J = 2,4
DO 200 I = 2,IMP
NPCNT = I + NXZ*(J-1) + IMP*(K-1)
UFLD(NPCNT) = .5 * (U(I-1,J,K) + U(I,J,K))
VFLD(NPCNT) = .5 * (V(I,J-1,K) + V(I,J,K))
WFLD(NPCNT) = .5 * (W(I,J,K-1) + W(I,J,K))
IF (K.EQ.2) THEN ! NO SLIP BC AT THE WALL U(1)=0 FOR K=1
NP1 = I + NXZ*(J-1)
UFLD(NP1) = 0.
VFLD(NP1) = 0.
WFLD(NP1) = 0.
ENDIF
C
200 CONTINUE
C
ENDIF
RETURN
END
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C S U B R O U T I N E W R S L I C E
C
C PURPOSE:
C --------
C THIS SUBROUTINE WRITES THE VELOCITY
C AND PRESSURE FIELD AT THE GRID POINTS ON X OR Y PLANE.
C THE VELOCITY FIELD IS AVERAGED AT THE GRID POINTS.
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE WRSLICE(IXYZ)
C
USE cpara
USE cfld
USE cbola
C ***WRITING OUTPUT***
print*,'writing velocity slice...'
C ***TIME LOOP***
IF(IXYZ==1) THEN
I=2
C
WRITE(9,200) T0,JM,KM
C ***SPACE LOOPS***
DO 90 K = 1,KM ! AF-28SEP02 AFTER REMOVING THE AVFIELD AND LEAVING AVGRPO
Z = ZETA(K)
IF(K==1) Z=0.0
DO 90 J = 2,JM+1
Y = DY*(J-1.5)
NPCNT = I + 4*KMP*(J-1) + 4*(K-1)
WRITE(9,300) Y,Z
. ,UFLD(NPCNT),VFLD(NPCNT),WFLD(NPCNT)
. ,PRESS(NPCNT)
90 CONTINUE
C
ELSEIF(IXYZ==2) THEN
C
J=2
WRITE(11,200) T0,IM,KM
C ***SPACE LOOPS***
DO 100 K = 1,KM
Z=ZETA(K)
IF(K==1) Z=0.0
DO 100 I = 2,IM+1
X = DX*(I-1.5)
NPCNT = I + NXZ*(J-1) + IMP*(K-1)
WRITE(11,300) X,Z
. ,UFLD(NPCNT),VFLD(NPCNT),WFLD(NPCNT)
. ,PRESS(NPCNT)
100 CONTINUE
C
ENDIF
C
200 FORMAT('ZONE T="T =',G14.5,'"'
. ,' F=POINT, I=',I4,' J=',I4)
300 FORMAT( 1X,6G14.5)
C
RETURN
END
C================================================================================