module compute
  implicit none
  save
  private
  public :: getnumbers, stats
  
  integer(4), parameter :: stat = 1005
  
contains
  !-------------------------------------------------------------------------------------------------------------------
  subroutine stats()
    use data
    use io
    
    real(8) uu(my), vv(my), ww(my), um(my), factor(0:mx/2-1)
    real(8) umax, nu, tau
    integer i,j,k, ire, iim
    real(4), allocatable, dimension(:,:) :: uxz, vxz, wxz
    allocate(uxz(mx,mz),vxz(mx,mz),wxz(mx,mz))
    factor=2.0d0
    factor(0)=1.0d0
    
    uu = 0.0d0
    vv = 0.0d0
    ww = 0.0d0
    um = 0.0d0

    do j = 1,my
       call read_plane_xz(u,uxz,j)
       call read_plane_xz(v,vxz,j)
       call read_plane_xz(w,wxz,j)
       um(j)=uxz(1,1)
       do i = 0,mx/2-1
          ire=2*i+1
          iim=2*i+2
          do k = 1,mz
             if (i.gt.0.or.k.gt.1) then
                uu(j)=uu(j)+factor(i)*(uxz(ire,k)**2+uxz(iim,k)**2)
                vv(j)=vv(j)+factor(i)*(vxz(ire,k)**2+vxz(iim,k)**2)
                ww(j)=ww(j)+factor(i)*(wxz(ire,k)**2+wxz(iim,k)**2)
             endif
          enddo
       enddo
    enddo

    open(stat,file="stats.dat",form='formatted')
    do j = 1,my
       write(stat,'(5(1p1e15.6))') ypos(j),um(j),uu(j),vv(j),ww(j)
    enddo
    close(stat)

    ! Calculating wall shear stress, and related quantities

    umax=maxval(um) 
    nu=umax/Re      ! Re is based on umax
    tau=nu*abs(um(2)-um(1))/(ypos(2)-ypos(1))
    utau=sqrt(tau)
    retau=utau/nu

    write(*,*) "utau = ",utau, " nu = ",nu," Retau = ",retau

    deallocate(uxz,vxz,wxz)

  end subroutine stats
  !-------------------------------------------------------------------------------------------------------------------
  subroutine getnumbers()
    use data, only: dig1, dig2, dig3

    character(1) u, t, h
    integer i, j, k

    do i = 1,10
       write(u,'(I1.1)') (i-1)
       dig1(i)=u
       do j = 1,10
          write(t,'(I1.1)') (j-1)
          dig2((j-1)*10+i)=t//u
          do k = 1,10
             write(h,'(I1.1)') (k-1)
             dig3(i+10*((j-1)+10*(k-1)))=h//t//u
          enddo
       enddo
    enddo
             
  end subroutine getnumbers
  !--------------------------------------------------------------------------------------------------------------
  function ypos(j)
    use data, only: my, pi
    implicit none
    
    integer(4) j
    real(8) ypos
    
    ypos = -cos(pi*real(j-1,8)/real(my-1,8))
    
  end function ypos

end module compute