subroutine cornermod2(xa,ya,za,k,im,jm)

      implicit none
      integer,parameter::nm=6
      integer::k,im,jm
      real,dimension(-1:2*im+1,-1:2*im+1,nm)::xa,ya,za
      integer::i,j,n,c1,m2,n2,km,iax,jax
      real,dimension(3)::p1,p2,p3,p4,pv
      real::cc1,cc2,b1,b2,a11,a12,a13,a21,a22,a23,xi,yi,zi
      real,dimension(k,k,nm,4,4,3)::pp
      real,dimension(k,k,nm,4,4,2)::ppll
      integer,dimension(k,k,4,4,3)::cind
      real,dimension(k,k,4,4,4,4)::cintp
      real,dimension(0:im+1,0:jm+1,nm)::hlat,hlon
      real,dimension(k,k,4,2)::mos
      real::pihlf,atang1,atang,tmpsum

      pihlf=2.*Atan(1.)
      km=2*im-1

       Do n=1,nm
       do j=0,jm+1
         jax=2*j-1
       do i=0,im+1
         iax=2*i-1
         xi=Xa(iax,jax,n)
         yi=Ya(iax,jax,n)
         zi=Za(iax,jax,n)
         HLON(i,j,n)=Atang1(yi,xi)
         HLAT(i,j,n)=Atang(zi,Sqrt(xi**2+yi**2),pihlf)
       End Do
       End Do
       End Do

      do n=1,nm
      do i=1,k
      do j=1,k
      do c1=1,4
        if(c1.eq.1)then
	  m2=2*(i-1)+1
	  n2=2*(j-1)+1
        else if(c1.eq.2)then
	  m2=km+2*(i-1-k)
	  n2=2*(j-1)+1
        else if(c1.eq.3)then
	  m2=2*(i-1)+1
	  n2=km+2*(j-1-k)
        else
	  m2=km+2*(i-1-k)
	  n2=km+2*(j-1-k)
        endif

	p1(1)=xa(m2,n2,n)
	p1(2)=ya(m2,n2,n)
	p1(3)=za(m2,n2,n)
	p2(1)=xa(m2+2,n2,n)
	p2(2)=ya(m2+2,n2,n)
	p2(3)=za(m2+2,n2,n)
	
	p3(1)=xa(m2,n2+2,n)
	p3(2)=ya(m2,n2+2,n)
	p3(3)=za(m2,n2+2,n)

	p4(1)=xa(m2+2,n2+2,n)
	p4(2)=ya(m2+2,n2+2,n)
	p4(3)=za(m2+2,n2+2,n)

	a11=p1(1)-p4(1)
	a12=p1(2)-p4(2)
	a13=p1(3)-p4(3)

	a21=p2(1)-p3(1)
	a22=p2(2)-p3(2)
	a23=p2(3)-p3(3)

        if(a11.ne.0)then
           cc1=a12*a21-a22*a11
	   cc2=a13*a21-a23*a11

	   b1=-cc2/cc1
	   b2=(cc2/cc1*a12-a13)/a11
        else
	   b1=-a13/a12
	   b2=(a13*a22-a23*a12)/(a12*a21)
        endif

	   if(p1(3).gt.0)then
	     pv(3)=1./(sqrt(b1**2+b2**2+1))
           else
	     pv(3)=-1./(sqrt(b1**2+b2**2+1))
           endif

	   pv(1)=b2*pv(3)
    	   pv(2)=b1*pv(3)

	xa(m2+1,n2+1,n)=pv(1)
	ya(m2+1,n2+1,n)=pv(2)
	za(m2+1,n2+1,n)=pv(3)

	call extend1(p4,pv,pp(i,j,n,c1,1,:))
	call extend1(p3,pv,pp(i,j,n,c1,2,:))
	call extend1(p2,pv,pp(i,j,n,c1,3,:))
	call extend1(p1,pv,pp(i,j,n,c1,4,:))

       if(n.eq.1)then
	call getmos(p3,p2,pp(i,j,n,c1,3,:),pp(i,j,n,c1,2,:),mos(i,j,c1,1))
	call getmos(p1,p4,pp(i,j,n,c1,1,:),pp(i,j,n,c1,4,:),mos(i,j,c1,2))
       endif

	enddo
	enddo
	enddo
	enddo

	call pptoll(pp,ppll,k)
	call cnind(pp,cind,k,im,jm)
	call cnintpco(ppll,cind,cintp,hlat,hlon,k,im,jm)
	
	print *,'Write to file cnmod.dat'
	open(11,file='../data_in/grid/cnmod.dat',form='unformatted')
	write(11)cind,cintp,mos
	close(11)

	end subroutine cornermod2