      parameter(it=164,jt=84,kt=132)

c j=wall normal direction, k=spanwise
      dimension u_64(jt,kt),v_64(jt,kt),w_64(jt,kt),fblend(jt)


c- section 0 ------initialization-----------------------------------------      

c 
c fix blending function (see Eq. 4 in TSFP5 paper)
      yblend=0.1
      ywidth=0.5
      do j=2,njm1
         yp_s=0.5*(yc(1,j,1)+yc(1,j-1,1))-yc(1,1,1)
         yp_n=yc(1,njm1,1)-0.5*(yc(1,j,1)+yc(1,j-1,1))
         yp1=min(yp_s,yp_n)
         arg=1.*(yp1-yblend)/ywidth
         fblend(j)=0.5*(1-tanh(arg))
         fblend(j)=max(0.1,fblend(j))
      end do


c- section 1 ------ read inlet data. Do this once per time step -------

c turbulent time scale=T in Eq. 22.8 in Lecture notes, Chapt. 22
      ustar1=1/25
      delta1=0.5
      tturb=0.05*delta1/ustar1
      a=exp(-dt(itstep)/tturb) !dt=computational time step
      b=(1.-a**2)**0.5

c filename='synth_uvw_non_const_dy_N.dat' where N=time step. When
c timestep=5001, let N=4999 and so on
      open(unit=35,file=filename)
      do k= 2,nkm1
      do j= 2,njm1
          read(35,*)u_64(j,k),v_64(j,k),w_64(j,k)
      end do
      end do


      do k= 2,nkm1
      do j= 2,njm1
          uprim=u_64(j,k)*fblend(j) !v'_1 in Eq. 22.6 in Lecture notes, Chapt. 22
          vprim=v_64(j,k)*fblend(j) !v'_2
          wprim=w_64(j,k)*fblend(j) !v'_3

c ufluct_inlet=V'_1^m in Eq. 22.7 in Lecture notes, Chapt. 22
          ufluct_inlet(j,k)=a*ufluct_inlet(j,k)+b*uprim
          phi(1,j,k,u)=ufluct_inlet(j,k)+umean(j) !=v_1 in Eq. 22.9 in Lecture notes, Chapt. 22
          phi(1,j,k,v)=a*phi(1,j,k,v)+b*vprim !=v_2 
          phi(1,j,k,w)=a*phi(1,j,k,w)+b*wprim !=v_3
       end do
       end do

