! from interp.f90 from Talha on 2/04/08 MODULE SHARED_DATA ! This module defines all the global varicables and dynamic allocatable arrays IMPLICIT NONE ! INTEGER :: nx,ny,nz,nxt,nyt,nzt,iconbc,iperio,interx,intery,interz, & isgsm,irdger,iwtger,jconbc,ichann,nxo,nyo,nzo,nxi,nyi,nzi,nxto,nyto,nzto, & ic_length,ic_height,ic_width,iconbc_old,ibwest,ibeast REAL :: cony_old,dpmean,length,width,height,myu,nu,re,retau, & length_old,width_old,height_old, & roh,dt,time,rmsfli,rmsflm,rmsflo,rmsflc,rmsflb,dmsfl,del,pi,dpmean_old, & re_old,cony_new CHARACTER(LEN=30) :: unit4,unit21,unit24 REAL, DIMENSION(:,:,:), ALLOCATABLE :: u,v,w,p,varray,farray,test REAL, DIMENSION(:), ALLOCATABLE :: x,y,z,dx,dy,dz,xc,yc,zc,dxs,sys,dzs REAL, DIMENSION(:), ALLOCATABLE :: xo,yo,zo,xco,yco,zco REAL, DIMENSION(:,:), ALLOCATABLE :: sgc,facx,facy,facz INTEGER, DIMENSION(:,:), ALLOCATABLE :: i2,j2,k2 REAL, DIMENSION(:,:,:,:,:), ALLOCATABLE :: xinlet,yinlet ! END MODULE SHARED_DATA ! !***************************************************************************** !***************************************************************************** PROGRAM INTERPOLATION !***************************************************************************** !****************************************************************************** USE SHARED_DATA IMPLICIT NONE ! CALL INITIA CALL SAVEVF ! STOP END PROGRAM INTERPOLATION ! !****************************************************************************** SUBROUTINE INITIA !****************************************************************************** USE SHARED_DATA IMPLICIT NONE ! CALL INPUT_PARAMETERS CALL ARRAY_ALLOCATION CALL open_FILES CALL INITIALIZATION CALL DOMAIN_SIZE_new CALL NEW_GRID CALL READ_FLOW CALL INITIA_BC_UPDATE CALL INITIAL_OUT ! RETURN END SUBROUTINE INITIA !****************************************************************************** SUBROUTINE READ_FLOW !****************************************************************************** ! This procedure reads the velocit flow field to restart the simulation from ! previously calculated flow field USE SHARED_DATA IMPLICIT NONE ! OPEN(UNIT=4,FILE=unit4,STATUS='OLD') CALL READ_VEL_FIELD CLOSE(UNIT=4,STATUS='KEEP') ! ! Renormalisation ! !CALL RMSFL !CALL CMPNON CALL CHECK_PARA ! RETURN END SUBROUTINE READ_FLOW ! !****************************************************************************** SUBROUTINE READ_VEL_FIELD !****************************************************************************** !-------This procedure reads velocity and pressure fields for future runs. USE SHARED_DATA IMPLICIT NONE INTEGER :: i,j,k,m,N,nsmplg,idir1=2,idir2=2 CHARACTER(LEN=30) :: file_tecplot ! irdger=0 WRITE(*,*)'READVF' IF (irdger == 1) THEN READ(4) nsmplg READ(4,*)((sgc(i,J),i=1,nx+1),j=1,ny+1) END IF ! WRITE(*,*)'READVF: dpmean=',dpmean READ(4,*)length_old,height_old,width_old,re_old,dt,time,nxi,nyi,nzi,dpmean_old write(21,*)length_old,height_old,width_old,re_old,dt,time,nxi,nyi,nzi,dpmean_old write(*,*)length_old,height_old,width_old,re_old,dt,time,nxi,nyi,nzi,dpmean_old CALL OLD_GRID ! ! u velocity READ(4,*)(((test(i,j,k),i=-1,nxo+2),j=-1,nyo+2),k=-1,nzo+2) file_tecplot='u-old' CALL TECPLOT_3D(test,-1,2,nxo,nyo,nzo,1,file_tecplot,3,idir1,idir2) CALL INTERP (u,p,1) file_tecplot='u-new' CALL TECPLOT_3D(u,-1,2,nx,ny,nz,0,file_tecplot,3,idir1,idir2) ! ! v velocity READ(4,*)(((test(i,j,k),i=-1,nxo+2),j=-1,nyo+2),k=-1,nzo+2) file_tecplot='v-old' CALL TECPLOT_3D(test,-1,2,nxo,nyo,nzo,1,file_tecplot,3,idir1,idir2) CALL INTERP (v,p,2) file_tecplot='v-new' CALL TECPLOT_3D(v,-1,2,nx,ny,nz,0,file_tecplot,3,idir1,idir2) ! ! w velocity READ(4,*)(((test(i,j,k),i=-1,nxo+2),j=-1,nyo+2),k=-1,nzo+2) file_tecplot='w-old' CALL TECPLOT_3D(test,-1,2,nxo,nyo,nzo,1,file_tecplot,3,idir1,idir2) CALL INTERP (w,p,3) file_tecplot='w-new' CALL TECPLOT_3D(w,-1,2,nx,ny,nz,0,file_tecplot,3,idir1,idir2) ! ! p velocity READ(4,*)(((test(i,j,k),i=0,nxto),j=0,nyto),k=0,nzto) file_tecplot='p-old0' CALL TECPLOT_3D(test,-1,2,nxo,nyo,nzo,1,file_tecplot,3,idir1,idir2) CALL OLD_RESTRAT file_tecplot='p-old' CALL TECPLOT_3D(test,-1,2,nxo,nyo,nzo,1,file_tecplot,3,idir1,idir2) CALL INTERP (varray,p,0) file_tecplot='p-new' CALL TECPLOT_3D(p,0,1,nx,ny,nz,0,file_tecplot,3,idir1,idir2) ! ! psi velocity - maybe it required. PRINT*,'M=4.5',p(1,ny/2+1,nz/2) READ(4,*)(((test(i,j,k),i=-1,nxo+2),j=-1,nyo+2),k=-1,nzo+2) CALL INTERP (varray,p,4) ! OPEN(15,FILE='u1.dat',STATUS='UNKNOWN') j=NY/2+1 k=NZ/2 DO i=1,nxt WRITE(15,162) x(i),u(i,j,k),v(i,j,k),w(i,j,k),p(i,j,k) END DO CLOSE(15) ! OPEN(15,FILE='u2.dat',STATUS='UNKNOWN') i=NX/2 k=NZ/2 DO j=1,nyt WRITE(15,162) y(j),u(i,j,k),v(i,j,k),w(i,j,k),p(i,j,k) END DO CLOSE(15) ! OPEN(15,FILE='u3.dat',STATUS='UNKNOWN') i=NX/2 j=NY/2 DO k=1,nzt WRITE(15,162) z(k),u(i,j,k),v(i,j,k),w(i,j,k),p(i,j,k) END DO CLOSE(15) ! 121 FORMAT('title="',A10,'"') 114 FORMAT('variables="x","y","u","v","w","p"') 122 FORMAT('zone t="z1",i='I3',j='I3',f=point') 162 FORMAT(9E14.6) ! OPEN(15,FILE='up1-xy.dat',STATUS='UNKNOWN') WRITE(15,121) "up1-xy.dat" WRITE(15,114) WRITE(15,122) nxt,nyt k=nzt/2 DO j=1,nyt DO i=1,nxt WRITE(15,162) x(i),y(j),u(i,j,k),v(i,j,k),w(i,j,k),p(i,j,k) END DO END DO CLOSE(15) ! OPEN(15,FILE='up1-xz.dat',STATUS='UNKNOWN') WRITE(15,121) "up1-xz.dat" WRITE(15,114) WRITE(15,122) nxt,nzt j=6 DO k=1,nzt DO i=1,nxt WRITE(15,162) x(i),z(k),u(i,j,k),v(i,j,k),w(i,j,k),p(i,j,k) END DO END DO CLOSE(15) ! OPEN(15,FILE='up1-yz.dat',STATUS='UNKNOWN') WRITE(15,121) "up1-yz.dat" WRITE(15,114) WRITE(15,122) nzt,nyt i=NX/2 DO j=1,nyt DO k=1,nzt WRITE(15,162) z(k),y(j),u(i,j,k),v(i,j,k),w(i,j,k),p(i,j,k) END DO END DO CLOSE(15) ! RETURN END SUBROUTINE READ_VEL_FIELD ! !****************************************************************************** SUBROUTINE INTERP (temp1,temp2,m) !****************************************************************************** ! This procedure intepolates the pressure and velocity to new flow field. USE SHARED_DATA IMPLICIT NONE REAL, DIMENSION(-1:nx+2,-1:ny+2,-1:nz+2):: temp1 REAL, DIMENSION(0:nx+1,0:ny+1,0:nz+1):: temp2 INTEGER, INTENT(IN) :: m INTEGER :: i,j,k,jbeg REAL :: depdpm,dpmean1 ! CALL FACTOR(x,xc,xo,xco,facx,i2,nx,nxo,ic_length,length,length_old,interx,m,1) CALL FACTOR(y,yc,yo,yco,facy,j2,ny,nyo,ic_height,height,height_old,intery,m,2) CALL FACTOR(z,zc,zo,zco,facz,k2,nz,nzo,ic_width,width,width_old,interz,m,3) ! jbeg=1 IF(m.ne.2) jbeg=0 temp1=0.0e0 temp2=0.0e0 ! DO k=1,nz DO j=jbeg,nyt DO i=1,nx temp1(i,j,k) & =facx(i,1)*facy(j,1)*facz(k,1)*test(i2(i,1),j2(j,1),k2(k,1)) & +facx(i,2)*facy(j,1)*facz(k,1)*test(i2(i,2),j2(j,1),k2(k,1)) & +facx(i,1)*facy(j,2)*facz(k,1)*test(i2(i,1),j2(j,2),k2(k,1)) & +facx(i,2)*facy(j,2)*facz(k,1)*test(i2(i,2),j2(j,2),k2(k,1)) & +facx(i,1)*facy(j,1)*facz(k,2)*test(i2(i,1),j2(j,1),k2(k,2)) & +facx(i,2)*facy(j,1)*facz(k,2)*test(i2(i,2),j2(j,1),k2(k,2)) & +facx(i,1)*facy(j,2)*facz(k,2)*test(i2(i,1),j2(j,2),k2(k,2)) & +facx(i,2)*facy(j,2)*facz(k,2)*test(i2(i,2),j2(j,2),k2(k,2)) END DO END DO END DO 162 FORMAT(9E14.6) PRINT*,'m,temp2(1,ny/2+1,nz/2)',m,temp2(1,ny/2+1,nz/2) ! IF(iconbc_old == 0) THEN ! periodic bc in x DO i=-1,0 DO k=1,nz DO j=-1,ny+2 temp1(i,j,k)=temp1(i+nx,j,k) END DO END DO END DO DO i=1,2 DO k=1,nz DO j=-1,ny+2 temp1(i+nx,j,k)=temp1(i,j,k) END DO END DO END DO END IF ! periodic bc in z DO k=-1,0 DO j=-1,ny+2 DO i=-1,nx+2 temp1(i,j,k)=temp1(i,j,k+nz) END DO END DO END DO DO k=1,2 DO j=-1,ny+2 DO i=-1,nx+2 temp1(i,j,k+nz)=temp1(i,j,k) END DO END DO END DO ! change the pressure gradient pressure IF(M == 0) THEN do k=0,nzt DO j=0,nyt DO i=0,nxt temp2(i,j,k)=temp1(i,j,k) end do end do end do END IF RETURN END SUBROUTINE INTERP ! !***************************************************************************** SUBROUTINE INPUT_PARAMETERS !***************************************************************************** ! This procedure reads the input parameters form 'input' USE SHARED_DATA IMPLICIT NONE CHARACTER :: dummy*70,char_dummy*6 ! OPEN(1,FILE='input.par',STATUS='OLD') OPEN(15,FILE='input-check.dat') READ (1,300) dummy WRITE(15,300) dummy READ (1,300) dummy WRITE(15,300) dummy READ (1,300) dummy WRITE(15,300) dummy READ (1,*) nx,char_dummy WRITE(15,*) nx,char_dummy READ (1,*) NY,char_dummy WRITE(15,*) NY,char_dummy READ (1,*) NZ,char_dummy WRITE(15,*) NZ,char_dummy READ (1,*) nxo,char_dummy WRITE(15,*) nxo,char_dummy READ (1,*) nyo,char_dummy WRITE(15,*) nyo,char_dummy READ (1,*) nzo,char_dummy WRITE(15,*) nzo,char_dummy READ (1,*) interx,char_dummy WRITE(15,*) interx,char_dummy READ (1,*) intery,char_dummy WRITE(15,*) intery,char_dummy READ (1,*) interz,char_dummy WRITE(15,*) interz,char_dummy READ (1,*) iconbc,char_dummy WRITE(15,*) iconbc,char_dummy READ (1,*) jconbc,char_dummy WRITE(15,*) jconbc,char_dummy READ (1,*) ichann,char_dummy WRITE(15,*) ichann,char_dummy READ (1,*) RE,char_dummy WRITE(15,*) RE,char_dummy READ (1,*) retau,char_dummy WRITE(15,*) retau,char_dummy READ (1,*) cony_old,char_dummy WRITE(15,*) cony_old,char_dummy READ (1,*) cony_new,char_dummy WRITE(15,*) cony_new,char_dummy READ (1,*) ic_length,char_dummy WRITE(15,*) ic_length,char_dummy READ (1,*) ic_height,char_dummy WRITE(15,*) ic_height,char_dummy READ (1,*) ic_width,char_dummy WRITE(15,*) ic_width,char_dummy READ (1,*) iconbc_old,char_dummy WRITE(15,*) iconbc_old,char_dummy READ (1,*) ibwest,char_dummy WRITE(15,*) ibwest,char_dummy READ (1,*) ibeast,char_dummy WRITE(15,*) ibeast,char_dummy write(*,*)'ic_length=',ic_length write(*,*)'ic_height=',ic_height write(*,*)'ic_width=',ic_width CLOSE(1) CLOSE(15) ! pi=2.0e0*ASIN(1.0) nxt=nx+1 nyt=ny+1 nzt=nz+1 nxto=nxo+1 nyto=nyo+1 nzto=nzo+1 roh=1.0e0 myu=1.0e0/RE nu=myu/roh dpmean=0.0e0 irdger=0 iwtger=0 iperio=1 ! 300 FORMAT(A70) RETURN END SUBROUTINE INPUT_PARAMETERS ! !****************************************************************************** SUBROUTINE ARRAY_ALLOCATION !****************************************************************************** ! This subroutine allocates the size of allocatable arrays USE SHARED_DATA IMPLICIT NONE ! ALLOCATE(u(-1:nx+2,-1:ny+2,-1:nz+2),v(-1:nx+2,-1:ny+2,-1:nz+2), & w(-1:nx+2,-1:ny+2,-1:nz+2),p(-1:nx+2,-1:ny+2,-1:nz+2),sgc(0:nx+1,ny+1),& xinlet(-1:ny+2,-1:nzt,3,2,2),yinlet(-1:nx+2,-1:nzt,3,2,2), & varray(-1:nx+2,-1:ny+2,-1:nz+2),test(-1:nxo+2,-1:nyo+2,-1:nzo+2)) ! ALLOCATE (x(nxt),y(nyt),z(nzt),dx(-1:nx+2),dy(-1:ny+2),dz(-1:nz+2), & xc(0:nxt),yc(0:nyt),zc(0:nzt),dxs(nxt),sys(nyt),dzs(nzt)) ! ALLOCATE(xo(nxo+1),yo(nyo+1),zo(nzo+1)) ALLOCATE(xco(0:nxo+1),yco(0:nyo+1),zco(0:nzo+1)) ALLOCATE(i2(0:nxt,2),j2(0:nyt,2),k2(0:nzt,2)) ALLOCATE(facx(0:nxt,2),facy(0:nyt,2),facz(0:nzt,2)) RETURN END SUBROUTINE ARRAY_ALLOCATION ! !****************************************************************************** SUBROUTINE open_FILES !****************************************************************************** ! This procedure open files USE SHARED_DATA IMPLICIT NONE ! ! unit4='restrt-64_129_64.dat' unit4='restrt241h6.dat' unit21='logfile1h6.dat' unit24='restrt-in-64_129_64.dat' ! OPEN(UNIT=4,FILE=unit4,STATUS='UNKNOWN') OPEN(UNIT=21,FILE=unit21,STATUS='UNKNOWN') OPEN(UNIT=24,FILE=unit24,STATUS='UNKNOWN') ! RETURN END SUBROUTINE open_FILES ! !****************************************************************************** SUBROUTINE DOMAIN_SIZE_new !****************************************************************************** ! This procedure determines the size of computational domain on which USE SHARED_DATA IMPLICIT NONE ! IF(ichann == 1) THEN height=2.0e0 length=2.0e0*pi width=2.0e0/3.0e0*pi ! length=4.0e0*pi ! width=4.0e0/3.0e0*pi ! ! Boundary Layer Flow ! ELSE IF(ichann == 0) THEN width=1.0e0 length=2.40e1 height=8.0e0 ELSE del=height*0.5e0 length=1.28e1*del height=2.0e0 width=4.0e0*del END IF ! RETURN END SUBROUTINE DOMAIN_SIZE_new ! !******************************************************************************* SUBROUTINE INITIALIZATION !******************************************************************************* ! This procedure performs the initializtion of arrays USE SHARED_DATA IMPLICIT NONE ! xinlet=0.0e0 yinlet=0.0e0 U=0.0e0 V=0.0e0 W=0.0e0 P=0.0e0 varray=0.0e0 ! RETURN END SUBROUTINE INITIALIZATION ! !****************************************************************************** SUBROUTINE OLD_GRID !****************************************************************************** ! This procedure defines the computational grid for the old domain USE SHARED_DATA IMPLICIT NONE INTEGER :: i,j,k ! DO i=1,nxto xo(i)=length_old/REAL(nxo)*(i-1) END DO DO j=1,nyto yo(j)=height_old/REAL(nyo)*(j-1) END DO DO k=1,nzto zo(k)=width_old/REAL(nzo)*(k-1) END DO ! !----- Nonuniform grid for Y ! IF (cony_old /= 0.0e0) THEN DO j=1,nyto yo(j)=1.0e0+TANH(cony_old*(2.0e0*(j-1.0e0)/nyo-1.0e0))/TANH(cony_old) END DO ENDIF ! yo(1)=0.0e0 yo(nyto)=height_old ! DO i=1,nxo xco(i)=0.5e0*(xo(i)+xo(i+1)) END DO xco(0)=-xco(1) xco(nxt)=xo(nxt)+xco(1) DO j=1,nyo yco(j)=0.5e0*(yo(j)+yo(j+1)) END DO yco(0)=yo(1) yco(nyto)=yo(nyto) DO k=1,nzo zco(k)=0.5e0*(zo(k)+zo(k+1)) END DO zco(0)=-zco(1) zco(nzto)=zo(nzto)+zco(1) ! OPEN(15,FILE='x-old.dat',STATUS='UNKNOWN') DO i=1,nxto WRITE(15,*) i,xo(i),xco(i) END DO CLOSE(15) OPEN(15,FILE='y-old.dat',STATUS='UNKNOWN') DO j=1,nyto WRITE(15,*) j,yo(j),yco(j) END DO CLOSE(15) OPEN(15,FILE='z-old.dat',STATUS='UNKNOWN') DO k=1,nzto WRITE(15,*) k,zo(k),zco(k) END DO CLOSE(15) RETURN END SUBROUTINE OLD_GRID ! !****************************************************************************** SUBROUTINE NEW_GRID !****************************************************************************** ! This procedure defines the computational grid for the domain on which ! the velocity flow field to be interpolated USE SHARED_DATA IMPLICIT NONE INTEGER :: i,j,k ! DO i=1,nxt x(i)=length/REAL(nx)*(I-1) END DO ! DO j=1,nyt y(j)=height/REAL(ny)*(J-1) END DO ! DO k=1,nzt z(k)=width/REAL(nz)*(K-1) END DO ! !----- Nonuniform grid for Y ! IF (cony_new.NE.0.0e0) THEN IF(ichann == 1) THEN DO j=1,nyt y(j)=1.0e0+TANH(cony_new*(2.0e0*(J-1.0e0)/NY-1.0e0))/TANH(cony_new) END DO ELSE IF(ichann == 0) THEN DO j=1,nyt y(j)=(1.0e0+TANH(cony_new*(1.0e0*(J-1.0e0)/NY-1.0e0))/TANH(cony_new))*height END DO END IF ENDIF ! y(1)=0.0e0 y(nyt)=height ! DO i=1,nx xc(i)=0.5e0*(x(i)+x(i+1)) END DO xc(0)=-xc(1) xc(nxt)=x(nxt)+xc(1) ! DO j=1,ny yc(j)=0.5e0*(y(j)+y(j+1)) END DO yc(0)=y(1) yc(nyt)=y(nyt) ! DO k=1,nz zc(k)=0.5e0*(z(k)+z(k+1)) END DO zc(0)=-zc(1) zc(nzt)=z(nzt)+zc(1) ! OPEN(15,FILE='x-new.dat',STATUS='UNKNOWN') DO i=1,nxt WRITE(15,*) i,x(i),xc(i) END DO CLOSE(15) ! OPEN(15,FILE='y-new.dat',STATUS='UNKNOWN') DO j=1,nyt WRITE(15,*) j,y(j),yc(j) END DO CLOSE(15) ! OPEN(15,FILE='z-new.dat',STATUS='UNKNOWN') DO k=1,nzt WRITE(15,*) k,z(k),zc(k) END DO CLOSE(15) ! DO i=1,nx dx(i)=x(i+1)-x(i) END DO ! dx(-1)=dx(1) dx(0)=dx(1) dx(nx+1)=dx(nx) dx(nx+2)=dx(nx) ! DO j=1,ny dy(j)=y(j+1)-y(j) END DO ! dy(-1)=dy(1) dy(0)=dy(1) ! dy(ny+1)=dy(ny) dy(ny+2)=dy(ny) ! DO k=1,nz dz(k)=z(k+1)-z(k) END DO ! dz(-1)=dz(1) dz(0)=dz(1) dz(nz+1)=dz(nz) dz(nz+2)=dz(nz) ! DO i=1,nx+1 dxs(i)=0.5e0*(dx(I-1)+dx(i)) END DO ! DO j=1,ny+1 sys(j)=0.5e0*(dy(J-1)+dy(j)) END DO ! DO k=1,nz+1 dzs(k)=0.5e0*(dz(K-1)+dz(k)) END DO ! RETURN END SUBROUTINE NEW_GRID ! !****************************************************************************** SUBROUTINE INITIA_BC_UPDATE !****************************************************************************** ! This procedure sets the boundary conditions for flow field in INITIA ! subroutine. USE SHARED_DATA INTEGER :: i,j,k,iuvw,IS ! ! Wall B.C. ! DO k=1,nzt DO i=-1,nx+2 yinlet(i,k,1,1,1)=u(i,0,k) yinlet(i,k,2,1,1)=v(i,1,k) yinlet(i,k,3,1,1)=w(i,0,k) yinlet(i,k,1,2,1)=u(i,nyt,k) yinlet(i,k,2,2,1)=v(i,nyt,k) yinlet(i,k,3,2,1)=w(i,nyt,k) END DO END DO ! DO IS=1,2 DO iuvw=1,3 DO k=1,nzt DO i=-1,nx+2 yinlet(i,k,iuvw,IS,2)=yinlet(i,k,iuvw,IS,1) END DO END DO END DO END DO ! ! Periodic B.C. ! ! periodic IF(iconbc == 0) THEN DO i=-1,0 DO k=1,nz DO j=-1,ny+2 u(i,j,k)=u(i+nx,j,k) v(i,j,k)=v(i+nx,j,k) w(i,j,k)=w(i+nx,j,k) END DO END DO END DO ! DO i=1,2 DO k=1,nz DO j=-1,ny+2 u(i+nx,j,k)=u(i,j,k) v(i+nx,j,k)=v(i,j,k) w(i+nx,j,k)=w(i,j,k) END DO END DO END DO ! ! Convective B.C. ! ELSE IF(iconbc == 1) THEN DO k=1,nzt DO j=-1,ny+2 IF(iconbc_old == 1.AND.ibwest == 1) THEN v(0,j,k)=0.5e0*(v(0,j,k)+v(1,j,k)) w(0,j,k)=0.5e0*(w(0,j,k)+w(1,j,k)) p(0,j,k)=0.5e0*(p(0,j,k)+p(1,j,k)) END IF IF(iconbc_old == 1.AND.ibeast == 1) THEN v(nxt,j,k)=0.5e0*(v(nx,j,k)+v(nxt,j,k)) w(nxt,j,k)=0.5e0*(w(nx,j,k)+w(nxt,j,k)) p(nxt,j,k)=0.5e0*(p(nx,j,k)+p(nxt,j,k)) END IF xinlet(j,k,1,1,1)=u(1,j,k) xinlet(j,k,2,1,1)=v(0,j,k) xinlet(j,k,3,1,1)=w(0,j,k) xinlet(j,k,1,2,1)=u(nxt,j,k) xinlet(j,k,2,2,1)=v(nxt,j,k) xinlet(j,k,3,2,1)=w(nxt,j,k) END DO END DO DO IS=1,2 DO iuvw=1,3 DO k=1,nzt DO j=-1,ny+2 xinlet(j,k,iuvw,IS,2)=xinlet(j,k,iuvw,IS,1) END DO END DO END DO END DO END IF ! ! Update Z boundary values for periodic B.C. ! DO k=-1,0 DO j=1,nyt DO i=1,nxt u(i,j,k)=u(i,j,k+nz) v(i,j,k)=v(i,j,k+nz) w(i,j,k)=w(i,j,k+nz) END DO END DO END DO ! DO k=1,2 DO j=1,nyt DO i=1,nxt u(i,j,k+nz)=u(i,j,k) v(i,j,k+nz)=v(i,j,k) w(i,j,k+nz)=w(i,j,k) END DO END DO END DO ! RETURN END SUBROUTINE INITIA_BC_UPDATE ! !****************************************************************************** SUBROUTINE INITIAL_OUT !***************************************************************************** ! This procedure writes the output for initial flow field USE SHARED_DATA IMPLICIT NONE INTEGER :: nxh,nyh,nzh,i,j,k REAL :: uconv ! WRITE(21,150) 150 FORMAT('********************* INPUT PARAMETERS ********************'//) ! WRITE(21,151)length,height,width,roh,re,nx,ny,nz,nu 151 FORMAT(' length=',E9.3,' height=',E9.3,' width=',E9.3, & ' Roh=',E9.3,' Re=',E9.3,//' nx=',I5,' ny=',I5,' nz=',I5,// & ' nu=',E10.3,//'x, y and z are...') ! WRITE(21,1511)(x(i),i=1,nxt) WRITE(21,1512)(y(j),j=1,nyt) WRITE(21,1511)(z(k),k=1,nzt) WRITE(21,1513) 1513 FORMAT(' dx, and dy are...') ! WRITE(21,1511)(dx(i),i=-1,nx+2) WRITE(21,1512)(dy(j),j=-1,ny+2) WRITE(21,1511)(dz(k),k=-1,nz+2) 1511 FORMAT(5F15.5) 1512 FORMAT(6F12.5) ! nxh=NX/2 nyh=NY/2 nzh=NZ/2 ! WRITE(21,18)time WRITE(21,198)y(nyh+1),z(nzh+2) 18 FORMAT(' time=',E15.5/) 198 FORMAT(' @Y=',E15.5,2X,' @Z=',E15.5/) WRITE(21,16) 16 FORMAT(4X,'X',12X,'U @CL',12X,'V @CL',12X,'W @CL',12X,'P @CL',12X, & ' DPdx @CL') ! DO i=1,nx WRITE(21,17) xc(i),u(i,nyh,nzh),v(i,nyh,nzh),w(i,nyh,nzh), & p(i,nyh,nzh),(p(i,nyh,nzh)-p(I-1,nyh,nzh))/dxs(i) 17 FORMAT(1X,F12.5,2X,5E15.5) END DO ! WRITE(21,199)x(nxh+1),z(nzh+2) 199 FORMAT(' @X=',E15.5,2X,' @Z=',E15.5/) ! DO j=1,ny WRITE(21,17) yc(j),u(nxh,j,nzh),v(nxh,j,nzh),w(nxh,j,nzh),p(nxh,j,nzh) END DO ! WRITE(21,200)x(nxh+1),y(nyh+1) 200 FORMAT(' @X=',E15.5,2X,' @Y=',E15.5/) ! DO k=1,nz WRITE(21,17) zc(k),u(nxh,nyh,k),v(nxh,nyh,k),w(nxh,nyh,k),p(nxh,nyh,k) END DO ! CALL RMSFL ! WRITE(21,1904) rmsfli,rmsflm,rmsflo,dmsfl 1904 FORMAT(' mass flow at_inlet=',E15.8,/' mass flow at_middle=', & E15.8,/' mass flow at_outlet=',E15.8,/' Difference in massflow=',E15.8) ! ! Convective B.C. ! uconv=rmsfli/height/width/roh 121 FORMAT('title="',A10,'"') 114 FORMAT('variables="x","y","u","v","w","p"') 122 FORMAT('zone t="z1",i='I3',j='I3',f=point') 162 FORMAT(6E14.6) ! RETURN END SUBROUTINE INITIAL_OUT ! !****************************************************************************** SUBROUTINE RMSFL !****************************************************************************** !-------This procedure calculates mass flux at inlet and outlet ! and their difference. USE SHARED_DATA IMPLICIT NONE INTEGER :: i,j,k ! rmsfli=0.0e0 rmsflm=0.0e0 rmsflo=0.0e0 rmsflc=0.0e0 rmsflb=0.0e0 ! DO k=1,nz DO j=1,ny rmsfli=rmsfli+u(1,j,k)*dy(j)*dz(k)*roh rmsflm=rmsflm+u(NX/2+1,j,k)*dy(j)*dz(k)*roh rmsflo=rmsflo+u(nxt,j,k)*dy(j)*dz(k)*roh rmsflc=rmsflc+xinlet(j,k,1,2,2)*dy(j)*dz(k)*roh END DO END DO ! DO k=1,nz DO i=1,nx rmsflb=rmsflb+yinlet(i,k,2,1,2)*dx(i)*dz(k)*roh rmsflb=rmsflb-yinlet(i,k,2,2,2)*dx(i)*dz(k)*roh END DO END DO ! dmsfl=rmsfli+rmsflb-rmsflo ! RETURN END SUBROUTINE RMSFL ! !****************************************************************************** SUBROUTINE CHECK_PARA !****************************************************************************** ! This procedure checks the number of grid points,domain and renold number ! from restart file and compare these values with those given by input file. USE SHARED_DATA IMPLICIT NONE ! IF ((nxo /= NX).OR.(nyo /= NY).OR.(nzo /= NZ)) THEN WRITE(21,155) 155 FORMAT(' Input parameters different. Abort...') WRITE(21,*)'nxo=',nxo,' NX=',nx WRITE(21,*)'nyo=',nyo,' NY=',ny WRITE(21,*)'nzo=',nzo,' NZ=',nz ELSE IF (length_old /= length) THEN WRITE(21,155) WRITE(21,*)'length_old=',length_old,' length=',length ELSE IF (height_old /= height) THEN WRITE(21,155) WRITE(21,*)'height_old=',height_old,' height=',height ELSE IF (width_old /= width) THEN WRITE(21,155) WRITE(21,*)'width_old=',width_old,' width=',width END IF ! RETURN END SUBROUTINE CHECK_PARA ! !****************************************************************************** SUBROUTINE CMPNON !****************************************************************************** ! This oprocedure performs the renormalization of velocity field to ensure Umean (Bulk)=1 USE SHARED_DATA IMPLICIT NONE INTEGER :: i,j,k,rkstep REAL :: factu,factp,retau1 ! retau1=re/(rmsflo/width/height) WRITE(21,*)'Re=',re,' retau=',retau1 factu=retau1/re factp=factu*factu WRITE(21,*)'factu=',factu,' factp=',factp ! u=u*factu v=v*factu w=w*factu w=w*factp ! dpmean=dpmean*factp ! RETURN END SUBROUTINE CMPNON ! !****************************************************************************** SUBROUTINE SAVEVF !****************************************************************************** ! This procedure saves the velocity flow field at the end of program USE SHARED_DATA IMPLICIT NONE INTEGER :: i,j,k ! IF (iwtger == 1) THEN REWIND 24 WRITE(24,*)((sgc(i,J),i=1,nx+1),j=1,ny+1) END IF WRITE(*,*)length,height,width,re,dt,time,nx,ny,nz,dpmean ! GOTO 1 WRITE(24,*)length,height,width,re,dt,time,nx,ny,nz,dpmean WRITE(24,*)(((u(i,j,k),i=-1,nx+2),j=-1,ny+2),k=-1,nz+2) WRITE(24,*)(((v(i,j,k),i=-1,nx+2),j=-1,ny+2),k=-1,nz+2) WRITE(24,*)(((w(i,j,k),i=-1,nx+2),j=-1,ny+2),k=-1,nz+2) WRITE(24,*)(((p(i,j,k),i=0,nxt),j=0,nyt),k=0,nzt) WRITE(24,*)(((varray(i,j,k),i=-1,nx+2),j=-1,ny+2),k=-1,nz+2) 1 CONTINUE 121 FORMAT('title="',A10,'"') 113 FORMAT('variables="x","y","z","p"') 114 FORMAT('variables="x","y","u","v","w","p"') 122 FORMAT('zone t="z1",i='I3',j='I3',f=point') 123 FORMAT('zone t="z1",i='I3',j='I3',k='I3',f=point') 162 FORMAT(6E14.6) OPEN(15,FILE='up2-xy.dat',STATUS='UNKNOWN') WRITE(15,121) "up2-xy.dat" WRITE(15,114) WRITE(15,122) nxt,nyt k=nzt/2 DO j=1,nyt DO i=1,nxt WRITE(15,162) x(i),y(j),u(i,j,k),v(i,j,k),w(i,j,k),p(i,j,k) END DO END DO CLOSE(15) OPEN(15,FILE='up2-xz.dat',STATUS='UNKNOWN') WRITE(15,121) "up2-xz.dat" WRITE(15,114) WRITE(15,122) nxt,nzt j=6 DO k=1,nzt DO i=1,nxt WRITE(15,162) x(i),z(k),u(i,j,k),v(i,j,k),w(i,j,k),p(i,j,k) END DO END DO CLOSE(15) OPEN(15,FILE='up2-yz.dat',STATUS='UNKNOWN') WRITE(15,121) "up2-yz.dat" WRITE(15,114) WRITE(15,122) nzt,nyt i=NX/2 DO j=1,nyt DO k=1,nzt WRITE(15,162) z(k),y(j),u(i,j,k),v(i,j,k),w(i,j,k),p(i,j,k) END DO END DO CLOSE(15) OPEN(15,FILE='up3d.dat',STATUS='UNKNOWN') WRITE(15,121) "up3d.dat" WRITE(15,113) WRITE(15,123) NX/2,ny/2,nz+1 DO k=1,nz+1 DO j=1,ny/2 DO i=1,nx/2 WRITE(15,162) xc(i),yc(j),zc(k),p(i,j,k)+x(i) END DO END DO END DO CLOSE(15) ! RETURN END SUBROUTINE SAVEVF !***************************************************************************** SUBROUTINE FACTOR(x_tec,xc_tec,xo,xco,facx,i2,nx,nxo,ic_length,length, & length_old,interx,m,idir) !***************************************************************************** ! This procedure calculates the factors used for interpolation in y-direction ! USE SHARED_DATA IMPLICIT NONE INTEGER, INTENT(IN) :: nx,nxo,interx,m,idir,ic_length REAL x(nx+1),xc(0:nx+1),xo(nxo+1),xco(0:nxo+1),facx(0:nx+1,2) INTEGER i2(0:nx+1,2) INTEGER :: nxt,i,io,istagger=0,n_length,n2_length REAL :: X_tec(nx+1),xc_tec(0:nx+1),length,length_old,small=1.0E-3 nxt=nx+1 ! ! staggered or non-staggered IF((m >= 1.AND.m <= 3).AND.m == idir) istagger=1 x=x_tec xc=xc_tec IF(ic_length == 1.AND.((idir == 1.AND.length > length_old+small).& OR.(idir == 3.AND.length > length_old+small))) THEN write(*,*)'M,idir,ic_length=',m,idir,ic_length,length,length_old do i=1,nxt IF(x(i) > length_old+small) THEN n_length=x(i)/length_old ! write(*,*) x(i),n_length x(i)=x(i)-length_old*n_length END IF end do n2_length=n_length do i=1,nxt IF(xc(i) > length_old+small) THEN n_length=xc(i)/length_old n_length=min(n_length,n2_length) xc(i)=xc(i)-length_old*n_length END IF end do END IF ! no change outer: IF(interx == 0) THEN DO i=0,nxt i2(i,1)=I i2(i,2)=I facx(i,1)=1.0e0 facx(i,2)=0.0e0 END DO ! change ELSE IF(interx == 1) THEN ! staggered inner: IF(istagger == 1) THEN DO i=1,nxt DO io=1,nxo IF(x(i) >= xo(io).AND.x(i) /= xo(io+1)) THEN i2(i,1)=io i2(i,2)=io+1 facx(i,1)=(xo(io+1)-x(i))/(xo(io+1)-xo(io)) facx(i,2)=1.0e0-facx(i,1) END IF END DO END DO ! non-staggered ELSE IF(istagger == 0) THEN DO i=0,nxt DO io=0,nxo IF(xc(i) >= xco(io).AND.xc(i) /= xco(io+1)) THEN i2(i,1)=io i2(i,2)=io+1 facx(i,1)=(xco(io+1)-xc(i))/(xco(io+1)-xco(io)) facx(i,2)=1.0e0-facx(i,1) END IF END DO END DO END IF inner END IF outer RETURN END SUBROUTINE FACTOR ! !***************************************************** SUBROUTINE TECPLOT_3D(f,n_1,n_2,nx_tec,ny_tec,nz_tec,inew,& file_tecplot,idim,idir1,idir2) ! n_1: begins at n_1 ! n_2: ends at nx+n_2 or nxo+n_2 ! iold=0: new; 1: old ! file_tecplot: file name ! idim: 1, 2 or 3 dimension ! idir: 1 in x; 2 in y; 3 in z USE SHARED_DATA IMPLICIT NONE INTEGER :: i,j,k,n_1,n_2,nx_tec,ny_tec,nz_tec,nxt_tec,nyt_tec,nzt_tec,& inew,idim,idir,idir1,idir2,n_sec INTEGER, DIMENSION(3) :: ic_dim REAL :: f(n_1:nx_tec+n_2,n_1:ny_tec+n_2,n_1:nz_tec+n_2) REAL :: x_tec(nx_tec+1),y_tec(ny_tec+1),z_tec(nz_tec+1) CHARACTER(LEN=30) :: file_tecplot,file_tecplot_x n_sec=5 ! new ! IF(inew.eq.1) THEN ! nx_tec=nx ! ny_tec=ny ! nz_tec=nz ! old ! ELSE IF(inew.eq.o) THEN ! nx_tec=nxo ! ny_tec=nyo ! nz_tec=nzo ! END IF nxt_tec=nx_tec+1 nyt_tec=ny_tec+1 nzt_tec=nz_tec+1 IF(inew==0) THEN x_tec=x y_tec=y z_tec=z ELSE IF(inew==1) THEN x_tec=xo y_tec=yo z_tec=zo END IF ic_dim=0 IF(idim.ge.1) ic_dim(1)=1 IF(idim.ge.3) ic_dim(2)=1 IF(idim.eq.7) ic_dim(3)=1 ! 1 dimensional IF(ic_dim(1)==1) THEN idir=idir1 file_tecplot_x="" file_tecplot_x = TRIM(file_tecplot)//"-1d" IF(idir==1) file_tecplot_x = TRIM(file_tecplot_x)//"-x.plt" IF(idir==2) file_tecplot_x = TRIM(file_tecplot_x)//"-y.plt" IF(idir==3) file_tecplot_x = TRIM(file_tecplot_x)//"-z.plt" OPEN(15,FILE=file_tecplot_x,STATUS='UNKNOWN') write(*,*) file_tecplot_x IF(idir==1) THEN k=nz_tec/2 do i=1,nx_tec+1 write(15,162) x_tec(i),(f(i,j,k),j=1,ny_tec,ny_tec/n_sec) end do ELSE IF(idir==2) THEN k=nz_tec/2 do j=1,ny_tec+1 write(15,162) y_tec(j),(f(i,j,k),i=1,nx_tec,nx_tec/n_sec) end do ELSE IF(idir==3) THEN i=nx_tec/2 do k=1,nz_tec+1 write(15,162) z_tec(k),(f(i,j,k),j=1,ny_tec,ny_tec/n_sec) end do END IF CLOSE(15) END IF ! 2 dimensional IF(ic_dim(2)==1) THEN file_tecplot_x="" file_tecplot_x = TRIM(file_tecplot)//"-2d" idir=idir2 IF(idir==1) file_tecplot_x = TRIM(file_tecplot_x)//"-x.plt" IF(idir==2) file_tecplot_x = TRIM(file_tecplot_x)//"-y.plt" IF(idir==3) file_tecplot_x = TRIM(file_tecplot_x)//"-z.plt" OPEN(15,FILE=file_tecplot_x,STATUS='UNKNOWN') write(*,*) file_tecplot_x 121 FORMAT('title="',A10,'"') 111 FORMAT('variables="x","y","u"') 114 FORMAT('variables="x","y","u","v","w","p"') 122 FORMAT('zone t="z1",i='I3',j='I3',f=point') 162 FORMAT(9E14.6) ! IF(idir==3) THEN WRITE(15,121) "u2d-xy.plt" WRITE(15,111) WRITE(15,122) nxt_tec,nyt_tec k=nzt_tec/2 DO j=1,nyt_tec DO i=1,nxt_tec WRITE(15,162) x_tec(i),y_tec(j),f(i,j,k) END DO END DO ! ELSE IF(idir==2) THEN WRITE(15,121) "u2d-xz.plt" WRITE(15,111) WRITE(15,122) nxt_tec,nzt_tec j=6 DO k=1,nzt_tec DO i=1,nxt_tec WRITE(15,162) x_tec(i),z_tec(k),f(i,j,k) END DO END DO ! ELSE IF(idir==1) THEN WRITE(15,121) "u2d-yz.plt" WRITE(15,111) WRITE(15,122) nzt_tec,nyt_tec i=nx_tec/2 DO j=1,nyt_tec DO k=1,nzt_tec WRITE(15,162) z_tec(k),y_tec(j),f(i,j,k) END DO END DO END IF CLOSE(15) END IF IF(ic_dim(3)==1) THEN file_tecplot_x="" file_tecplot_x = TRIM(file_tecplot)//"-3d" file_tecplot_x = TRIM(file_tecplot_x)//".plt" OPEN(15,FILE=file_tecplot_x,STATUS='UNKNOWN') write(*,*) file_tecplot_x CLOSE(15) END IF RETURN END SUBROUTINE TECPLOT_3D SUBROUTINE OLD_RESTRAT ! Substracting the mean pressure component from the mean pressure gradient USE SHARED_DATA IMPLICIT NONE INTEGER :: i,j,k REAL :: dpdx_mean ! dpdx_mean=dpmean_old/length_old write(*,*)'dpmean,dpdx_mean=',dpmean_old,dpdx_mean do k=0,nzo+1 do j=0,nyo+1 do i=0,nxo+1 test(i,j,k)=test(i,j,k)-dpdx_mean*xo(i) end do end do end do RETURN END SUBROUTINE OLD_RESTRAT