Quiero calcular el valor z como la coordenada en el rango de x: -50 ~ 50 e y: -50 ~ 50 como el código siguiente.

program test
implicit none
! --- [local entities]
      real*8            :: rrr,th,U0,amp,alp,Ndiv
      real*8            :: pi,alpR,NR,Rmin,Rmax,z
      integer           :: ir, i, j

do i=0, 50
  do j=0, 50
      th=datan2(i,j)
      pi=datan(1.d0)*4.d0
!
      Ndiv= 24.d0            !! Number of circumferential division
      alp = 90.d0/180.d0*pi  !! phase [rad]
      U0  = 11.4d0           !! average velocity
      amp = 0.5d0            !! amplitude of velocity
      Rmin = 10              !! [m] 
      Rmax = 50              !! [m]
      NR = 6.d0              !! Number of radial division
!
      rrr=dsqrt(i**2+j**2)
      ir=int((rrr-Rmin)/(Rmax-Rmin)*NR)
      alpR=2.d0*pi/dble(Ndiv)*dble(mod(ir,2))
      z=U0*(1.d0+amp*dsin(0.5d0*Ndiv*th+alp+alpR))

  write(*,*) 'i, j, z'
  write(*,*) i, j, z 
 end do
end do

stop
end program test

Pero no pude hacer que funcione como el error a continuación. Creo que porque i, j están en datan (i, j). ¿Cómo debo cambiar este código?

test.f90:10.16:

      th=datan2(i,j)
                1
Error: 'y' argument of 'datan2' intrinsic at (1) must be REAL
test.f90:21.16:

      rrr=dsqrt(i**2+j**2)
                1
Error: 'x' argument of 'dsqrt' intrinsic at (1) must be REAL