module shared_data
!
  Implicit None
! Read Parameters
  Integer :: nx,ny,nz
  Real :: leng,hei,wid,re,dt,time,dpdx_mean
  Real, Dimension(:,:,:), Allocatable :: u,v,w,p
!
! Grid Parameters
  Real, Dimension(:), Allocatable :: x,y,z,xc,yc,zc
  Real :: cony
!
end module shared_data
!******************************************************************************
Program check_restrt
  Implicit None
!
!
  Call read_restrt
  Call grid_channel
  Call write_field
!
  STOP
End Program check_restrt
!******************************************************************************
Subroutine read_restrt
  Use shared_data
  Implicit None
!
!  This subroutine reads the variables from the binary restart file
!
  Open(Unit=22, File='restrt-nocontrol.rst', Status='OLD', Form='UNFormatTED')
  Read(22)leng,hei,wid,re,dt,time,nx,ny,nz,dpdx_mean
!
  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))
  u = 0.0e0;v = 0.0e0; w = 0.0e0; p = 0.0e0
!
  Read(22) u
  Read(22) v
  Read(22) w
  Read(22) p
!
  Close(22)
!
Return
End Subroutine read_restrt
!******************************************************************************
Subroutine grid_channel
!
  Use shared_data
  Implicit None
  Integer :: i,j,k
!
!  This subroutine defines the computational grid for the old domain
!
  Allocate(x(nx+1),y(ny+1),z(nz+1),xc(0:nx+1),yc(0:ny+1),zc(0:nz+1))
!
! cony sets spacing in y direction, warning: this is not always the same!
  cony = 2.04e0
!
  Do i = 1,nx+1
    x(i) = 1.0e0/Real(nx)*(i-1)
  End Do
  Do j = 1,ny+1
    y(j) = 1.0e0/Real(ny)*(j-1)
  End Do
  Do k = 1,nz+1
    z(k) = 1.0e0/Real(nz)*(k-1)
  End Do
!
! Sets up non-uniform grid for y
!
  If (cony > 1.0e-6) Then
    Do j = 1,ny+1
      y(j) = TANH(cony*(2.0e0*(j-1.0e0)/ny-1.0e0))/TANH(cony)
    End Do
    y=(y+1.0e0)*0.5e0
  End If
!
  y(1) = 0.0e0
  y(ny+1) = 1.0e0
!
  Do i = 1,nx
    xc(i) = 0.5e0*(x(i)+x(i+1))
  End Do
  xc(0) = -xc(1)
  xc(nx+1) = x(nx+1)+xc(1)
!
  Do j = 1,ny
    yc(j) = 0.5e0*(y(j)+y(j+1))
  End Do
  yc(0) = y(1)
  yc(ny+1) = y(ny+1)

  Do k = 1,nz
    zc(k) = 0.5e0*(z(k)+z(k+1))
  End Do
  zc(0) = -zc(1)
  zc(nz+1) = z(nz+1)+zc(1)
!
  x=x*leng; y=y*hei; z=z*wid
  xc=xc*leng; yc=yc*hei; zc=zc*wid
!
Return
End Subroutine grid_channel
!******************************************************************************
Subroutine write_field
!
  Use shared_data
  Implicit None
  Integer :: i,j,k
!
! This subroutine writes certain data, to check read
!
! Prints Data in xy plane
  Open(Unit=15, File='restrt-check-xy.plt', Status='Unknown')
  Write(15,121) "restrt-check-xy.plt"
  Write(15,114)
  Write(15,122) nx+1,ny+1
!
  k=nz+1/2
  Do j = 1,ny+1; Do i = 1,nx+1
    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)
!
! Prints Data in xz plane
  Open(15, File='restrt-check-xz.plt', Status='Unknown')
  Write(15,121) "restrt-check-xz.plt"
  Write(15,115)
  Write(15,122) nx+1,nz+1
!
  j=1
  Do k = 1,nz+1; Do i = 1,nx+1
    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)
!
! Prints Data in yz plane
  Open(15, File='restrt-check-yz.plt', Status='Unknown')
  Write(15,121) "restrt-check-yz.plt"
  Write(15,116)
  Write(15,122) nz+1,ny+1
!
  i=nx/2
  Do j = 1,ny+1; Do k = 1,nz+1
    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)
!
 121    Format('title="',a10,'"')
 114    Format('variables="x","y","u","v","w","p"')
 115    Format('variables="x","z","u","v","w","p"')
 116    Format('variables="z","y","u","v","w","p"')
 122    Format('zone t="z1",i='I3',j='I3',f=point')
 162    Format(6E14.6)
!
Return
End Subroutine write_field
!******************************************************************************