!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
      program metrics  
!     ******************************************************************
!     *                                                                *
!     *  convert map coordinates into cartesian                        *
!     *                                                                *
!     ******************************************************************

      implicit none

      integer, parameter :: m=64, np=m+1

      real, dimension(3,-1:np,-1:np) :: xtab
      common /qpan/ xtab
!----------------------------------------------------------------------
      include 'param_o.h'  
!----------------------------------------------------------------------
      integer, parameter:: kstrd=4
      integer, parameter:: lstrd=kstrd, km=(im-1)*kstrd+1
!-----------------------------------------------------------------------
      real,dimension(:,:,:),allocatable:: xa,ya,za
      real,dimension(:,:),allocatable:: glon,glat,hg,xmp0,ymp0
      integer,dimension(:,:),allocatable:: kmp0
      real,dimension(:,:),allocatable:: sqh2,rsqh2
      real,dimension(:,:,:),allocatable:: dxvdx,dxvdy,dyvdx,dyvdy,dzvdx,dzvdy
      real,dimension(:,:,:),allocatable:: q11,q12,q22,sqv,sqh,rsqv,rsqh
      real,dimension(:,:,:),allocatable:: vcosp,vsinp,vlam,vphi
      real,dimension(:,:,:),allocatable:: hcosp,hsinp,hlam,hphi
      real,dimension(:,:,:),allocatable:: u2us,v2us,u2vs,v2vs
      real,dimension(:,:,:),allocatable:: qbv11,qbv12,qbv22

      real,dimension(:,:,:,:),allocatable:: qd11,qd12,qd21,qd22,ds
      real,dimension(:,:,:,:,:,:),allocatable::qvh,qhv,qvh2,qhv2

      real,dimension(:,:,:),allocatable:: rarea,qh11,qh12,qh22
      real,dimension(:,:),allocatable:: rarea2
      real,dimension(:,:,:),allocatable:: ds2

      real,dimension(:,:,:),allocatable::xh,yh,zh,xv,yv,zv

      real,dimension(:,:,:),allocatable:: sqv_norm
!-----------------------------------------------------------------------
!      real,dimension(:,:),allocatable:: sqvl2
!-----------------------------------------------------------------------
      real,parameter:: a=6.371e6
!-----------------------------------------------------------------------
      logical ncor (0:im+1,0:jm+1)  
!-----------------------------------------------------------------------
      real xs (km), xmp (km, km), ymp (km, km)  
!-----------------------------------------------------------------------
      real xm (2),xc(3),dxcdxm(3,2),dxmdxc(2,3)  
!-----------------------------------------------------------------------
      character * 37 cfile1  
      character * 37 cfile2  
      real:: mfac
      real, parameter:: pi =3.1415926 
      real, parameter:: pihlf =pi*0.5 
      real, parameter:: rad =pi/180.  

      real:: dz, x
      integer:: k,kaux,i,j,l,n,it,jt,list1,list2,kmap
      real:: snal,csal,csbt,snbt,csgm,sngm
      real:: x0,y0,z0,x1,y1,z1,x2,y2,z2
      real:: qs11,qs22,qs12,dxhdx,dyhdx,dzhdx,dxhdy,dyhdy,dzhdy
      real:: dlat,dlon,xlon,ylat,xmap,ymap
      real:: vlamda,vcosl,vsinl,ax,ay,bx,by,det,rdet
      real:: sqv2,rsqv2
      real:: sqv_max

      real, external :: atang1
!-----------------------------------------------------------------------

      allocate(xa(1-kstrd:km+kstrd,1-lstrd:km+lstrd,nm))
      allocate(ya(1-kstrd:km+kstrd,1-lstrd:km+lstrd,nm))
      allocate(za(1-kstrd:km+kstrd,1-lstrd:km+lstrd,nm))

      allocate(glon (igm, jgm))
      allocate(glat (igm, jgm))
      allocate(hg   (igm, jgm))
      allocate(xmp0 (igm, jgm))
      allocate(ymp0 (igm, jgm))
      allocate(kmp0 (igm, jgm))

      allocate(sqh2 (0:im+1, 0:jm+1))
      allocate(rsqh2(0:im+1, 0:jm+1))

      allocate(dxvdx(0:im+1,0:jm+1,nm))
      allocate(dxvdy(0:im+1,0:jm+1,nm))
      allocate(dyvdx(0:im+1,0:jm+1,nm))
      allocate(dyvdy(0:im+1,0:jm+1,nm))
      allocate(dzvdx(0:im+1,0:jm+1,nm))
      allocate(dzvdy(0:im+1,0:jm+1,nm))

      allocate(q11 (0:im+1,0:jm+1,nm))
      allocate(q12 (0:im+1,0:jm+1,nm))
      allocate(q22 (0:im+1,0:jm+1,nm))
      allocate(sqv (0:im+1,0:jm+1,nm))
      allocate(sqh (0:im+1,0:jm+1,nm))
      allocate(rsqv(0:im+1,0:jm+1,nm))
      allocate(rsqh(0:im+1,0:jm+1,nm))
      allocate(sqv_norm(0:im+1,0:jm+1,nm))

!      allocate(sqvl2(-1:im+2,-1:jm+2))

      allocate(vcosp (0:im+1,0:jm+1,nm))
      allocate(vsinp (0:im+1,0:jm+1,nm))
      allocate(vlam  (0:im+1,0:jm+1,nm))
      allocate(vphi  (0:im+1,0:jm+1,nm))
      allocate(hcosp (0:im+1,0:jm+1,nm))
      allocate(hsinp (0:im+1,0:jm+1,nm))
      allocate(hlam  (0:im+1,0:jm+1,nm))
      allocate(hphi  (0:im+1,0:jm+1,nm))

      allocate(u2us (0:im+1, 0:jm+1,nm))
      allocate(v2us (0:im+1, 0:jm+1,nm))
      allocate(u2vs (0:im+1, 0:jm+1,nm))
      allocate(v2vs (0:im+1, 0:jm+1,nm))

      allocate(qbv11 (0:im+1,0:jm+1,nm))
      allocate(qbv12 (0:im+1,0:jm+1,nm))
      allocate(qbv22 (0:im+1,0:jm+1,nm))

      allocate(qvh(0:im+1,0:jm+1,nm,4,2,2))
      allocate(qhv(0:im+1,0:jm+1,nm,4,2,2))
      allocate(qvh2(0:im+1,0:jm+1,nm,4,2,2))
      allocate(qhv2(0:im+1,0:jm+1,nm,4,2,2))

      allocate(qd11(0:im+1,0:jm+1,nm,4))
      allocate(qd12(0:im+1,0:jm+1,nm,4))
      allocate(qd21(0:im+1,0:jm+1,nm,4))
      allocate(qd22(0:im+1,0:jm+1,nm,4))
      allocate(ds  (0:im+1,0:jm+1,nm,4))

      allocate(rarea(0:im+1,0:jm+1,nm))
      allocate(qh11 (0:im+1,0:jm+1,nm))
      allocate(qh12 (0:im+1,0:jm+1,nm))
      allocate(qh22 (0:im+1,0:jm+1,nm))

      allocate(rarea2 (0:im+1,0:jm+1))

      allocate(ds2  (0:im+1,0:jm+1,4))

      allocate(xh(0:im+1,0:im+1,nm))
      allocate(yh(0:im+1,0:im+1,nm))
      allocate(zh(0:im+1,0:im+1,nm))
      allocate(xv(0:im+1,0:im+1,nm))
      allocate(yv(0:im+1,0:im+1,nm))
      allocate(zv(0:im+1,0:im+1,nm))
 
      mfac=0.5

!-----------------------------------------------------------------------
!
! initialize smooth cube
!
      call infin3  
!
! read table
!
      open (unit = 10, file = 'round.dat', &
       form = 'unformatted')
      read (10) xtab  
      close (10)  
      print * , 'read round.dat file'  
!-----------------------------------------------------------------------
!
!  define map space
!

      dz=2./(km-1)  
      x=-1.  

      do k=1,km/2  
        xs(k)=x  
        x=x+dz  
      enddo  

      xs(km/2+1)=0.  
      kaux=km/2  

      do k=km/2+2,km  
        xs (k)=-xs(kaux)  
        kaux=kaux-1  
      enddo  

      do k=1,km  
      do l=1,km  
        xmp(k,l)=xs(k)  
      enddo  
      enddo  

      do k=1,km  
      do l=1,km  
        ymp(k,l)=xs(l)  
      enddo  
      enddo  
!-----------------------------------------------------------------------
!
! derive absolute cartesian coordinates
!
      do n=1,nm  
      do k=1,km  
      do l=1,km  
        xm(1)= xmp(k,l)  
        xm(2)= ymp(k,l)  
          call xmtoxc (xm, xc, dxcdxm, n)  
        xa(k,l,n) = xc (1)  
        ya(k,l,n) = xc (2)  
        za(k,l,n) = xc (3)  
      enddo  
      enddo  
      enddo  

!
! use gnomonic cube, comment out to use conformal or smoothed cube
!

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

      call bococ (xa)  
      call bococ (ya)  
      call bococ (za)  

! rotation
!       alfa = 45. * rad  
!       beta = 54. * rad  
!       alfa = 2. * rad  
!       beta = 0. * rad  
    
      csal=cos(alfa)
      snal=sin(alfa)
      csbt=cos(beta)
      snbt=sin(beta)
      csgm=cos(gamm)
      sngm=sin(gamm)

      do n=1,nm
      do i=1-kstrd,km+kstrd
      do j=1-kstrd,km+kstrd
          x0=xa(i,j,n)
          y0=ya(i,j,n)
          z0=za(i,j,n)
          x1= x0*csal+y0*snal
          y1=-x0*snal+y0*csal
          z1=z0
          x2=x1
          y2= y1*csbt+z1*snbt
          z2=-y1*snbt+z1*csbt
          x1= x2*csgm+y2*sngm
          y1=-x2*sngm+y2*csgm
          z1=z2
          xa(i,j,n)=x1
          ya(i,j,n)=y1
          za(i,j,n)=z1
      enddo
      enddo
      enddo 

!-----------------------------------------------------------------------
!
! grid components: v-points
!
      do n=1,nm  
        l=1+lstrd/2
      do j=1,jm-1  
        k=1+kstrd/2
      do i=1,im-1  
        dxvdx(i,j,n)=( xa(k+2,l-2,n)-xa(k-2,l+2,n) )*mfac 
        dyvdx(i,j,n)=( ya(k+2,l-2,n)-ya(k-2,l+2,n) )*mfac 
        dzvdx(i,j,n)=( za(k+2,l-2,n)-za(k-2,l+2,n) )*mfac 
        dxvdy(i,j,n)=( xa(k+2,l+2,n)-xa(k-2,l-2,n) )*mfac 
        dyvdy(i,j,n)=( ya(k+2,l+2,n)-ya(k-2,l-2,n) )*mfac 
        dzvdy(i,j,n)=( za(k+2,l+2,n)-za(k-2,l-2,n) )*mfac 

        qs11 = dxvdx(i,j,n)**2+dyvdx(i,j,n)**2+dzvdx(i,j,n)**2
        qs22 = dxvdy(i,j,n)**2+dyvdy(i,j,n)**2+dzvdy(i,j,n)**2 
     
        qs12 = dxvdx(i,j,n)*dxvdy(i,j,n) &
             + dyvdx(i,j,n)*dyvdy(i,j,n) &
             + dzvdx(i,j,n)*dzvdy(i,j,n)

          sqv2=qs11*qs22-qs12*qs12
!###########################################################
          if(sqv2==0.) then
            print *,'METRCS: sqv2,i,j=',sqv2,i,j
          end if
!###########################################################
          rsqv2=1./sqv2
            q11(i,j,n)= qs22*rsqv2
            q22(i,j,n)= qs11*rsqv2
            q12(i,j,n)=-qs12*rsqv2
          sqv(i,j,n)=sqrt (sqv2)
          rsqv(i,j,n)=1./sqv(i,j,n)

        k=k+kstrd  
      enddo  
        l=l+lstrd  
      enddo  
      enddo  

            call bocov1_cb(sqv,im,jm)
            call bocov1_cb(rsqv,im,jm)
            call bocov2_cb(q11,q22,im,jm)
            call bocov3_cb(q12,im,jm)

         sqv_max=maxval(sqv)

         sqv_norm=sqv/sqv_max

!################################################
!        sqvl2=0.
!      do i=0,im
!      do j=0,jm
!        sqvl2(i,j)=sqv(i,j,1)
!      end do
!      end do
! 
!      do i=1,im-1
!        sqvl2(i  ,-1)=sqvl2(i,2   )
!        sqvl2(i,jm+1)=sqvl2(i,jm-2)
!      end do
!
!      do j=1,jm-1
!        sqvl2(-1  ,j)=sqvl2(2   ,j)
!        sqvl2(im+1,j)=sqvl2(im-2,j)
!      end do
!
!################################################

      do n = 1, nm  
        l = 1 - lstrd / 2  
      do j = 0, jm  
        k = 1 - kstrd / 2  
      do i = 0, im  
        xv(i,j,n)=xa(k,l,n)
        yv(i,j,n)=ya(k,l,n)
        zv(i,j,n)=za(k,l,n)
        vcosp (i, j, n) = sqrt (xa (k, l, n) **2 + ya (k, l, n) **2)  
        vsinp (i, j, n) = za (k, l, n)  
        vlam (i, j, n) = atang1 (ya (k, l, n), xa (k, l, n) )  
        vphi (i, j, n) = asin (vsinp (i, j, n) )  
        k = k + kstrd  
      enddo  
        l = l + lstrd  
      enddo  
      enddo  

!
! grid components: h-points
!
      do j=1,jm  
      do i=1,im  
      if (i.eq.1.and.j.eq.1.or.i.eq.1.and.j.eq.jm.or.i.eq.im.and.j.eq.1 &
      .or.i.eq.im.and.j.eq.jm) then
         ncor(i,j)=.false.  
      else 
         ncor(i,j)=.true.  
      endif  
      enddo  
      enddo  

      n = 1  
        l = 1  
      do j=1,jm  
        k = 1  
      do i=1,im  
        dxhdx=( xa(k+2,l-2,n)-xa(k-2,l+2,n) )*mfac 
        dyhdx=( ya(k+2,l-2,n)-ya(k-2,l+2,n) )*mfac 
        dzhdx=( za(k+2,l-2,n)-za(k-2,l+2,n) )*mfac 
        dxhdy=( xa(k+2,l+2,n)-xa(k-2,l-2,n) )*mfac 
        dyhdy=( ya(k+2,l+2,n)-ya(k-2,l-2,n) )*mfac 
        dzhdy=( za(k+2,l+2,n)-za(k-2,l-2,n) )*mfac 
        qs11 = dxhdx**2 + dyhdx**2 + dzhdx**2  
        qs22 = dxhdy**2 + dyhdy**2 + dzhdy**2  
        qs12 = dxhdx*dxhdy + dyhdx*dyhdy + dzhdx*dzhdy  
      if (ncor (i,j) ) then  
         sqh2(i,j) =sqrt(qs11*qs22 - qs12*qs12)  
      else  
         sqh2(i,j) =sqv(1,1,1)*3./4.  
      endif  
        rsqh2 (i, j) = 1./sqh2 (i,j)  
        rarea2(i,j)=1./(2.*sqh2(i,j)*a*a)
        ds2(i,j,1)=0.5*sqrt(qs11)*a
        ds2(i,j,3)=0.5*sqrt(qs22)*a
        k=k+kstrd  
      enddo  
        l=l+lstrd  
      enddo  

      do n=1,nm  
        l=1-lstrd 
      do j=0,jm+1  
        k=1-kstrd  
      do i=0,im+1  
        xh(i,j,n)=xa(k,l,n)
        yh(i,j,n)=ya(k,l,n)
        zh(i,j,n)=za(k,l,n)
        hcosp(i,j,n)=sqrt(xa(k,l,n)**2+ya(k,l,n)**2)  
        hsinp(i,j,n)=za  (k,l,n)  
        hlam (i,j,n)=atang1 (ya(k,l,n),xa(k,l,n) )  
        hphi (i,j,n)=asin (hsinp (i,j,n) )  
        k=k+kstrd  
      enddo  
        l=l+lstrd  
      enddo  
      enddo  

      call qvhqhv(xh,yh,zh,xv,yv,zv,q11,q12,q22,qvh,qhv,im,jm,nm,mfac)
      call qvhqhv2(xh,yh,zh,xv,yv,zv,q11,q12,q22,qh11,qh12,qh22,qvh2,  &
                  qhv2,im,jm,nm,mfac)
!!!!!!!      call qintc(xh,yh,zh,xv,yv,zv,q11,q12,q22,im,jm,nm,mfac)

         
         call qintc_2(xh,yh,zh,xv,yv,zv,im,jm,nm,mfac)

      n=1
      do i = 1, im  
      do j = 1, jm  
        dxhdx=(xv(i,j-1,n)+xv(i,j,n)-xv(i-1,j,n)-xv(i-1,j-1,n))*0.5 
        dyhdx=(yv(i,j-1,n)+yv(i,j,n)-yv(i-1,j,n)-yv(i-1,j-1,n))*0.5 
        dzhdx=(zv(i,j-1,n)+zv(i,j,n)-zv(i-1,j,n)-zv(i-1,j-1,n))*0.5 
        dxhdy=(xv(i,j,n)+xv(i-1,j,n)-xv(i,j-1,n)-xv(i-1,j-1,n))*0.5 
        dyhdy=(yv(i,j,n)+yv(i-1,j,n)-yv(i,j-1,n)-yv(i-1,j-1,n))*0.5
        dzhdy=(zv(i,j,n)+zv(i-1,j,n)-zv(i,j-1,n)-zv(i-1,j-1,n))*0.5
        qs11 = dxhdx**2 + dyhdx**2 + dzhdx**2  
        qs22 = dxhdy**2 + dyhdy**2 + dzhdy**2  
        ds2(i,j,2)=0.25*sqrt(qs11)*a
        ds2(i,j,4)=0.25*sqrt(qs22)*a
      enddo  
      enddo  

!-----------------------------------------------------------------------
         do n=1,nm  
           sqh (:,:,n) = sqh2 (:,:)  
           rsqh(:,:,n) = rsqh2(:,:)  
           rarea(:,:,n)=rarea2
           ds(:,:,n,:)=ds2
         enddo  
!-----------------------------------------------------------------------

         n=1
      do i=1,im1
      do j=1,jm1
      do k=1,4
        select case(k)
          case(1)
            it=i-1
            jt=j
          case(2) 
            it=i+1
            jt=j
          case(3) 
            it=i
            jt=j-1
          case default  
            it=i
            jt=j+1
        end select

          qd11(i,j,n,k)=q11(it,jt,n)*(dxvdx(it,jt,n)*dxvdx(i,j,n)  &
             +dyvdx(it,jt,n)*dyvdx(i,j,n)+dzvdx(it,jt,n)*dzvdx(i,j,n)) &
             +q12(it,jt,n)*(dxvdy(it,jt,n)*dxvdx(i,j,n) &
             +dyvdy(it,jt,n)*dyvdx(i,j,n)+dzvdy(it,jt,n)*dzvdx(i,j,n))
          qd12(i,j,n,k)=q12(it,jt,n)*(dxvdx(it,jt,n)*dxvdx(i,j,n)  &
             +dyvdx(it,jt,n)*dyvdx(i,j,n)+dzvdx(it,jt,n)*dzvdx(i,j,n)) &
             +q22(it,jt,n)*(dxvdy(it,jt,n)*dxvdx(i,j,n) &
             +dyvdy(it,jt,n)*dyvdx(i,j,n)+dzvdy(it,jt,n)*dzvdx(i,j,n))
          qd21(i,j,n,k)=q11(it,jt,n)*(dxvdx(it,jt,n)*dxvdy(i,j,n)  &
             +dyvdx(it,jt,n)*dyvdy(i,j,n)+dzvdx(it,jt,n)*dzvdy(i,j,n)) &
             +q12(it,jt,n)*(dxvdy(it,jt,n)*dxvdy(i,j,n) &
             +dyvdy(it,jt,n)*dyvdy(i,j,n)+dzvdy(it,jt,n)*dzvdy(i,j,n))
          qd22(i,j,n,k)=q12(it,jt,n)*(dxvdx(it,jt,n)*dxvdy(i,j,n)  &
             +dyvdx(it,jt,n)*dyvdy(i,j,n)+dzvdx(it,jt,n)*dzvdy(i,j,n)) &
             +q22(it,jt,n)*(dxvdy(it,jt,n)*dxvdy(i,j,n) &
             +dyvdy(it,jt,n)*dyvdy(i,j,n)+dzvdy(it,jt,n)*dzvdy(i,j,n))

         enddo
         enddo
         enddo


      do i=0,im
      do j=0,jm
        qbv11(i,j,n)=0.5*(q11(i,j,n)+q22(i,j,n)+2.*q12(i,j,n))
        qbv12(i,j,n)=0.5*(q22(i,j,n)-q11(i,j,n))
        qbv22(i,j,n)=0.5*(q11(i,j,n)+q22(i,j,n)-2.*q12(i,j,n))
      enddo
      enddo

!
! For cube
!
      do n=2,nm
        qd11(:,:,n,:)=qd11(:,:,1,:)
        qd12(:,:,n,:)=qd12(:,:,1,:)
        qd21(:,:,n,:)=qd21(:,:,1,:)
        qd22(:,:,n,:)=qd22(:,:,1,:)
        qbv11(:,:,n)= qbv11(:,:,1) 
        qbv12(:,:,n)= qbv12(:,:,1) 
        qbv22(:,:,n)= qbv22(:,:,1) 
      end do
     

!
! outfile for absolute coordintes
!

      list1 = 11  

      call bocoh_cb(sqh,im,jm)
      call bocoh_cb(rsqh,im,jm)
      call bocoh_cb(hcosp,im,jm)
      call bocoh_cb(hsinp,im,jm)
      call bocoh_cb(hlam,im,jm)
      call bocoh_cb(hphi,im,jm)


      print *,"writing file cubeinit ..."
      cfile1 = "gridinit.dat"  
      open (unit = list1, file = cfile1, form = "unformatted")  
      write (list1) sqv,sqh,q11,q12,q22,dxvdx,dxvdy,dyvdx, &
         dyvdy,dzvdx,dzvdy,rsqv,vcosp,vsinp,vlam,vphi,rsqh,hcosp, &
         hsinp,hlam,hphi,qd11,qd12,qd21,qd22,qvh,qhv,qbv11,qbv12,qbv22, &
         qvh2,qhv2,qh11,qh12,qh22,sqv_norm
      close (list1)  
      print * , 'wrote ', cfile1, ' file'

      open(list1,file="hvphi.dat", form="unformatted")
      write(list1)hphi,vphi
      close(list1)

      open(list1,file="hpos.dat", form="unformatted")
      write(list1)hlam,hphi
      close(list1)

      open(list1,file="fzeff.dat", form="unformatted")
      write(list1)rarea,ds
      close(list1)

!-----------------------------------------------------------------------
!
! prepare graphical output
!
!
! define glon and glat
!
       dlat = pi / (jgm - 1)  

       dlon = 2 * pi / igm  
       print *,dlon*180/pi,dlat*180/pi
       do i = 1, igm  
       xlon = (i - 1) * dlon  
       do j = 1, jgm  
       glon (i, j) = xlon  
       enddo  

       enddo  
       do j = 1, jgm  
       ylat = - pihlf + (j - 1) * dlat  
       do i = 1, igm  
       glat (i, j) = ylat  
       enddo  

       enddo  
       csal = cos (alfa)  
       snal = sin (alfa)  
       csbt = cos (beta)  
       snbt = sin (beta)  
!
! find xmp0, ymp0, kmp0
!
       do j = 1, jgm  
       do i = 1, igm  
       xlon = glon (i, j)  
       ylat = glat (i, j)  
       x2 = cos (ylat) * cos (xlon)  
       y2 = cos (ylat) * sin (xlon)  
       z2 = sin (ylat)  
       x0 = x2 * csgm - y2 * sngm
       y0 = x2 * sngm + y2 * csgm
       z0 = z2  
       x1 = x0  
       y1 = y0 * csbt - z0 * snbt  
       z1 = y0 * snbt + z0 * csbt  
       x0 = x1 * csal - y1 * snal  
       y0 = x1 * snal + y1 * csal  
       z0 = z1  
       xc (1) = x0  
       xc (2) = y0  
       xc (3) = z0  
       call xctoxm (xc, xm, dxmdxc, kmap)  
       xmap = xm (1)  
       ymap = xm (2)  
       xmap = (1. + xmap) / 2.  
       ymap = (1. + ymap) / 2.  
       xmp0 (i, j) = xmap  
       ymp0 (i, j) = ymap  
       kmp0 (i, j) = kmap  
       enddo  
       enddo  
!
! prepare conversion of covariant back to spherical winds
!
       do n = 1, nm  
       do j = 1, jm1  
       do i = 1, im1  
       vlamda = vlam (i, j, n)  
       vcosl = cos (vlamda)  
       vsinl = sin (vlamda)  
       ax = - vsinl * dxvdx (i, j, n) + vcosl * dyvdx (i, j, n)  
       ay = - vsinl * dxvdy (i, j, n) + vcosl * dyvdy (i, j, n)  
       bx = - vcosl * vsinp (i, j, n) * dxvdx (i, j, n) - vsinl * vsinp ( &
        i, j, n) * dyvdx (i, j, n) + vcosp (i, j, n) * dzvdx (i, j, n)
       by = - vcosl * vsinp (i, j, n) * dxvdy (i, j, n) - vsinl * vsinp ( &
        i, j, n) * dyvdy (i, j, n) + vcosp (i, j, n) * dzvdy (i, j, n)
       det = ax * by - ay * bx  
       rdet = 1. / det  
       u2us (i, j, n) = by * rdet  
       v2us (i, j, n) = - bx * rdet  
       u2vs (i, j, n) = - ay * rdet  
       v2vs (i, j, n) = ax * rdet  
       enddo  
       enddo  
       enddo  
!-----------------------------------------------------------------------
!
! outfile for graphics
!

       print *,igm,jgm,im,jm,nm
       list2 = 12  

       cfile2 = "grph.dat"  
       open (unit = list2, file = cfile2, form = "unformatted")  
       write (list2) glon,glat,xmp0,ymp0,kmp0,u2us,v2us,u2vs,v2vs
       close (list2)  
!test

       print * , 'wrote ', cfile2, ' file'  
       k=(im-1)/20
       print *,'k,km,',k,im
       open(11,file="cbdt.dat",form="unformatted")
       write(11)sqh
       close(11)

       print *,'corner locations:'

       print *,hlam(1,1,1),hphi(1,1,1)
       print *,hlam(im,1,1),hphi(im,1,1)
       print *,hlam(1,jm,1),hphi(1,jm,1)
       print *,hlam(im,jm,1),hphi(im,jm,1)
       print *,hlam(1,1,3),hphi(1,1,3)
       print *,hlam(im,1,3),hphi(im,1,3)
       print *,hlam(1,jm,3),hphi(1,jm,3)
       print *,hlam(im,jm,3),hphi(im,jm,3)
       print *,xh(4,3,1),yh(4,3,1),zh(4,3,1)
!      call qintc_2(xh,yh,zh,xv,yv,zv,im,jm,nm,mfac)
       
!test
!-----------------------------------------------------------------------
       stop  
       end  program metrics