program gnomonic

      implicit none
      include 'param_o.h'

      integer,parameter::kstrd=4,lstrd=kstrd,km=(im-1)*kstrd+1
      real,dimension(1-kstrd:km+kstrd,1-lstrd:km+lstrd,nm)::xa,ya,za
      real::xc,zc,a,dlm,dph,lmt,ph,lm0,ph0,d
      integer::k,l

      a=1./(sqrt(3.))
      dlm=3.1415926*0.5/(km-1)
      dph=dlm
      lm0=3.1415926*0.25
      ph0=-lm0

      do k=1,km
      do l=1,km
        lmt=lm0+(k-1)*dlm
        ph=ph0+(l-1)*dph
        xc=a/tan(lmt)
        zc=a*tan(ph)
        d=sqrt(xc**2+zc**2+a**2)

        xa(k,l,2)=xc/d
        ya(k,l,2)=a/d
        za(k,l,2)=zc/d
	if(za(k,l,2).lt.-0.9)then
	  print *,k,l,zc,d,xc,a
	  stop
        endif
!	print *,k,l,xa(k,l,2),ya(k,l,2),za(k,l,2)
      enddo
      enddo


      do l=1,(km+1)/2
        xa((km+1)/2,l,2)=0.
	za((km+1)/2,l,2)=-sqrt(1.-ya((km+1)/2,l,2)**2)
      enddo

      do l=1,(km+1)/2
        xa(l,l,2)=-za(l,l,2)
	ya(l,l,2)=sqrt(1.-2.*xa(l,l,2)**2)
      enddo

      do l=1,(km+1)/2
        ya(l,1,2)=-za(l,1,2)
	xa(l,1,2)=sqrt(1.-2.*ya(l,1,2)**2)
      enddo

      xa(1,1,2)=1./sqrt(3.)
      ya(1,1,2)=xa(1,1,2)
      za(1,1,2)=-xa(1,1,2)

      do l=2,(km+1)/2
      do k=1,l-1
        ya(k,l,2)=ya(l,k,2)
	xa(k,l,2)=-za(l,k,2)
	za(k,l,2)=-xa(l,k,2)
      enddo
      enddo
!
      do k=1,(km+1)/2
      do l=1,(km+1)/2
        ya(km-l+1,k,2)=ya(k,l,2)
	xa(km-l+1,k,2)=za(k,l,2)
	za(km-l+1,k,2)=-xa(k,l,2)
!
	ya(km-k+1,km-l+1,2)=ya(k,l,2)
	xa(km-k+1,km-l+1,2)=-xa(k,l,2)
	za(km-k+1,km-l+1,2)=-za(k,l,2)
!!
        ya(l,km-k+1,2)=ya(k,l,2)
        xa(l,km-k+1,2)=-za(k,l,2)
        za(l,km-k+1,2)=xa(k,l,2)

      enddo
      enddo

      do k=1,km
      do l=1,km
!	print *,k,l,xa(k,l,2),ya(k,l,2),za(k,l,2)
      enddo
      enddo

      xa(:,:,1)= ya(:,:,2)
      ya(:,:,1)=-xa(:,:,2)
      za(:,:,1)= za(:,:,2)

      xa(:,:,3)=-xa(:,:,1)
      ya(:,:,3)=-ya(:,:,1)
      za(:,:,3)= za(:,:,1)

      xa(:,:,4)= ya(:,:,1)
      ya(:,:,4)=-xa(:,:,1)
      za(:,:,4)= za(:,:,1)

      xa(:,:,5)=-za(:,:,1)
      ya(:,:,5)= ya(:,:,1)
      za(:,:,5)= xa(:,:,1)

      xa(:,:,6)= za(:,:,1)
      ya(:,:,6)= ya(:,:,1)
      za(:,:,6)= -xa(:,:,1)

      open(11,file='cbdt_gn.dat',form='unformatted')
      write(11)xa(1:km:2,1:km:2,1) 
      write(11)ya(1:km:2,1:km:2,1) 
      write(11)za(1:km:2,1:km:2,1) 
      close(11)

       open(12,file='gnom.dat',form='unformatted')
       write(12)xa,ya,za
       close(12)

      end program gnomonic