      SUBROUTINE mkgrid

      IMPLICIT NONE

      INCLUDE 'param.txt'
      INCLUDE 'param2.txt'
      INCLUDE 'common.txt'
      
      REAL pi
      
      LOGICAL done

      real x,dx

      real gam1,ksi1,gam2,ksi2,kk1,kk2,h
      integer ns
      real eta,zuni(0:kmax),runi(0:imax)
      real xacc
      parameter (xacc=1.0E-6)
      real lambda,theta
      external funcd
      real sech,rtnewt

      REAL yplus
c--------------------------------------------------------------------

      pi=4.0*ATAN(1.0)


c     theta-grid...


      dphi=2.0*pi/jmax


c     r-grid...

      IF(1.EQ.1)THEN  ! non-uniform grid

         DO i=0,imax
            runi(i)=i*(lr/2.0)/imax
         END DO

         ns=imax/1.5
         h=.55*(lr/2.0)
         eta=runi(ns)
         gam1=0.5
         ksi1=0.5*(lr/2.0)
         kk1=h/(ksi1*(1.0-tanh(gam1*(ksi1-eta))/tanh(gam1*ksi1)))

         DO i=0,ns
            ru(i)=kk1*ksi1*(1.0-tanh(gam1*(ksi1-runi(i)))/tanh(gam1*ksi1
     $           ))
         END DO

         kk2  = ((lr/2.0)-h)/((lr/2.0)-eta)
         ksi2 = (lr/2.0)

         lambda=kk1*gam1*ksi1*(sech(gam1*(ksi1-eta)))**2/tanh(gam1*ksi1)
      
         theta=rtnewt(funcd,0.5,5.0,xacc,kk2,lambda)
         gam2=theta/(2.0*(ksi2-eta))
      
         DO i=ns+1,imax
            ru(i)=h+kk2*(ksi2-eta)*(1.0-tanh(gam2*(ksi2-runi(i)))
     $           /tanh(gam2*(ksi2-eta)))
         END DO

         ru(0)=0.0
      
      ELSE   ! uniform grid

         ru(0)=0
         DO i=1,imax
            x= 1.*i/imax
            dx= 1.0
            ru(i)=ru(i-1)+dx
         END DO
         DO i=0,imax
            ru(i)=ru(i)*0.5*lr/ru(imax)
         END DO
      
      END IF

      DO i=1,imax
         rp(i)=(ru(i)+ru(i-1))/2.0
         dr(i)=(ru(i)-ru(i-1))
      END DO

      ru(i1)=ru(imax) + (ru(imax)-ru(imax-1))

      dr(i1) = dr(imax)
      dr(0)  = dr(1)

      rp(i1)=2.0*ru(imax) - rp(imax)

      rp(0)=ru(0) - (ru(1)-ru(0))/2.0


      done=.FALSE.
      ijet=i1
      DO i=1,imax
         IF(ru(i) .GT. 0.5 .AND. .NOT. done)THEN 
            ijet=i-1
            done=.TRUE.
         END IF
      END DO
            
c     z-grid...

      IF(1.eq.0)THEN    ! non-uniform grid



         DO k=0,kmax
            zuni(k)=k*lz/kmax
         END DO

c      ns=kmax/3.25
c      h=.65*lz
c      eta=zuni(ns)
c      gam1=7.

         ns=kmax/2.2
         h=.5*lz
         eta=zuni(ns)
         gam1=3.0


         ksi1=0.01
         kk1=h/(ksi1*(1.0-tanh(gam1*(ksi1-eta))/tanh(gam1*ksi1)))

         DO k=0,ns
            zw(k)=kk1*ksi1*(1.0-tanh(gam1*(ksi1-zuni(k)))/tanh(gam1*ksi1
     $           ))
         END DO

         kk2=(lz-h)/(lz-eta)
         ksi2=lz

         lambda=kk1*gam1*ksi1*(sech(gam1*(ksi1-eta)))**2/tanh(gam1*ksi1)

         theta=rtnewt(funcd,0.25,8.0,xacc,kk2,lambda)
         gam2=theta/(2.0*(ksi2-eta))
      
         DO k=ns+1,kmax
            zw(k)=h+kk2*(ksi2-eta)*(1.0-tanh(gam2*(ksi2-zuni(k)))
     $           /tanh(gam2*(ksi2-eta)))
         END DO

      ELSE   ! uniform grid

         zw(0)=0.  
         DO  k=1,kmax
            x = 1.*k/kmax 
            dx = 1.0
            zw(k)=zw(k-1)+dx 
         END DO
         DO k=0,kmax
            zw(k)=zw(k)*lz/zw(kmax)
         END DO

      END IF


      DO k=1,kmax
         zp(k)=(zw(k)+zw(k-1))/2.0
         dz(k)=(zw(k)-zw(k-1))
      END DO

      zw(k1)=zw(kmax) + (zw(kmax)-zw(kmax-1))

      dz(k1) = dz(kmax)
      dz(0)  = dz(1)

      zp(k1)=2.0*zw(kmax) - zp(kmax)

      zp(0)=zw(0) - (zw(1)-zw(0))/2.0


c *** Calculate mixing length of the subgrid scale motions:

      DO i=1,imax
         DO k=1,kmax

            delta(i,k) = ( rp(i) * dphi )**2 +
     &                   (     dr(i)    )**2 +
     &                   (     dz(k)    )**2 

            delbr(i,k) = ( rp(i) * dphi )**2 +
     &                   (     dr(i)    )**2 +
     &                   (     dz(k)    )**2

            delht(i,k) = ( 2*rp(i)*dphi )**2 +
     &                   (     dr(i)    )**2 +
     &                   (    2*dz(k)   )**2

         END DO
      END DO


      DO i=1,imax
         DO k=1,kmax

            lmix(i,k) = Cs * delta(i,k)
            
         END DO
      END DO

c...output r-phi grid to file to use in the program interpolating
c...turbulence in a box onto a cylinderical domain

      OPEN(15,FILE='grid.out',STATUS='UNKNOWN')

      WRITE(15,10) ijet,dphi
 10   FORMAT(1(I5),F17.12)

      DO i=0,i1
         WRITE(15,20) ru(i),rp(i)
 20      FORMAT(2(F17.12))
      END DO

      CLOSE(15)


      RETURN
      END

c---------------------------------------------------------------------
c---------------------------------------------------------------------

      
      function sech(x)
      real x
      
      sech=1.0/cosh(x)

      return
      end



      function rtnewt(funcd,x1,x2,xacc,k,l)

      integer jmax
      real rtnewt,x1,x2,xacc
      external funcd
      parameter (jmax=200)
      integer j
      real df,dx,f,k,l

      rtnewt=.5*(x1+x2)
      do j=1,jmax
         call funcd(rtnewt,k,l,f,df)
         dx=f/df
         rtnewt=rtnewt-dx
         if((x1-rtnewt)*(rtnewt-x2) .lt. 0.0) pause 'out brackets'
         if(abs(dx) .lt. xacc) return
      end do
      pause 'rtnewt exceeded maximum iterations'
      end 



      subroutine funcd(x,k,l,f,df)
      real k,l,x,f,df

      f=k*x - l*sinh(x)
      df=k-l*cosh(x)

      return
      end
