[lehoanganh07]

hi ,
I am a beginner in the coding. I will use fortran for my course. I'm trying to make the grid for my model, but when i ran, it had error. I know this forum is for OpenFOAM, but I really need help!
Could you help me to fix it!!
Thank you

program nozzle_mesh
implicit none
real,allocatable,dimension(:,::x,y
integer::i,j,num_xi,num_eta
real::d_xi,d_eta,xi,eta,le,wi,L1,L2,eta0
!eta0 is the hight of base
!L1 is the lenght of base
!L2 is the lenght of combustion chamber
!eta,xi is coordinate in computer plane, eta vetical axial, xi holizon axial
!d_xi; d_eta is distance between 2 node
!i,j is the node
!num_xi;num_eta is number of node
!beta_x;beta_y is the compressible coefficience of cluster grid
!xi0 is the position where cluster compressed in x direction
!xii;etaj distance to origin
real,parameter::e=2.71828
real::beta_x=5.0,beta_y=1.2,anpha,A,xi0,xii,etaj,c ,eta0,x0
open(1,"Nozlle.dat")
read*,le,hg,L1
read*,num_xi,num_eta
xi=hg
eta=le
call number(xi,num_xi,L1)
call number(eta,num_eta,eta0)
d_xi=xi/(num_xi-1)
d_eta=eta/(num_eta-1)
xi0=L1
allocate(x(num_xi,num_eta),y(num_xi,num_eta))
do j=1,num_xi
do i=1,num_eta
xii=d_xi*(j-1)
x0=log((1.0+(exp(beta_x)-1.0)*xi0)/(1.0+(exp(-beta_x)-1.0)*xi0))&
& /log(e)/2.0/beta_x
A=sinh(beta_x*x0)
x(i,j)=xi0*(sinh((xii-x0)*beta_x)+A)/A
if(xii<=L1 .and. i<=eta0/d_eta)y(i,j)=eta0
else if (xii<=L1 .and. i>eta0/d_eta) then
etaj=d_eta*(i-1)
c=log((beta_y+1)/(beta_y-1))
if (xii<=L1) anpha=0
if (xii>L1) anpha=0.5
y(i,j)=((beta_y+1)-(beta_y-1)*exp(-c*(etaj-1-anpha)/(1-anpha))) &
& /((2*anpha+1)*(1+exp(-c*(etaj-1-anpha))/(1-anpha)))
write(1,*),x(i,j),y(i,j)
end if
end do
write(1,*),x(i,j),y(i,j)
print*,x(i,j),y(i,j)
end do
end

//************************************************** *****************//
subroutine number(L,number,L1)
real,intent(inout)::L,number,L1
real::delta
integer::num1,num2
delta=L/(number-1)
num1=L1/delta
delta=L1/num1
num2=L2/delta
L2=num2*delta
L=L1+L2
number=num1+num2
return (L),(number),(L1)
end