!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
      Program metrcs
!     ******************************************************************
!     *                                                                *
!     *  Convert Map Coordinates into Cartesian                        *
!     *                                                                *
!     ******************************************************************

       implicit none
!----------------------------------------------------------------------
       include 'param_o.h'  
       integer,Parameter::i0=im+im-1                          
       integer,Parameter::i1=i0+i0-1
       integer,Parameter::img=i1+i0-1                    &
      , j0=i0,j1=i1,jmg=img                              &
      , iml=i0,jml=iml,lmg=img                           &
      , n1grid=img/2,nagrid=n1grid/3,ntay=30,nn=128                          
       real,parameter::flon0=270.  &
      , flat0=-90.,scen=1
       real::JUMBLE(nn),WORK(nn),RA(0:nn),QA(0:nn)
       integer,dimension(jmg)::imin,imax
       real,dimension(lmg)::s
       real,dimension(3)::xe
       real,dimension(0:img+1,0:img+1,2)::xg,yg,zg
       real,dimension(0:iml+1,0:iml+1,nm)::xa,ya,za
       real,dimension(0:im+1,0:im+1,nm)::xh,yh,zh,xv,yv,zv
       integer,dimension(0:im+1,0:jm+1,nm)::imod
       logical,dimension(0:im+1,0:jm+1,nm)::lv
       real,Dimension(0:im+1,0:jm+1,nm)::DXVDX,DXVDY,DYVDX,DYVDY,  &
                DZVDX,DZVDY 
       real::qs11,qs12,qs22,sqhc,dxhdx,dyhdx,dzhdx,dxhdy,dyhdy,dzhdy
       real,dimension(0:im+1,0:jm+1,nm)::sqv,sqh,rsqv,rsqh,q11,q12,q22
       real,dimension(0:im+1,0:jm+1,nm,4)::qd11,qd12,qd21,qd22
       real,dimension(0:im+1,0:jm+1,nm)::qbv11,qbv12,qbv22
       real,dimension(0:im+1,0:jm+1,nm)::vcosp,vsinp,vlam,vphi,  &
                                     hcosp,hsinp,hlam,hphi

!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
      real,Dimension(igm,jgm):: GLON , GLAT, HG, &
                XMP0, YMP0
      integer,dimension(igm,jgm)::KMP0
!-----------------------------------------------------------------------
      real,Dimension(0:im+1, 0:jm+1,nm):: u2us,v2us , &
                u2vs, v2vs,qh11,qh12,qh22  
!-----------------------------------------------------------------------
      Logical NCOR (0:im + 1, 0:jm + 1)  
       real,dimension(0:im+1,0:jm+1,nm,4,2,2)::qvh,qhv,qvh2,qhv2
!-----------------------------------------------------------------------
      Character * 01 cres  
      Character * 37 cfile1  
      Character * 37 cfile2  
      real:: mfac,dmap,xmap,ymap,atang1,csgm,sngm
      real::pihlf,pi,rad,dlat,dlon,xlon,ylat,csal,snal,csbt  &
          ,snbt,x0,y0,z0,x1,y1,z1,x2,y2,z2,vlamda,vcosl,vsinl,ax,ay,bx &
	  ,by,det,rdet
      integer::j,i,l,il,jl,n,imd,k,it,jt,nres,list1,list2,kmap
      CHARACTER (len=39) :: tfile=& 
          '../data_in/grid/SOCT0004.DAT'
      real,dimension(2)::xm
      real,dimension(0:im+1,0:jm+1,nm)::rarea
      real,dimension(0:im+1,0:jm+1,nm,4)::ds
      real,parameter::a=6.371e6

!-----------------------------------------------------------------------

      mfac=0.5
!-----------------------------------------------------------------------
!
! Initialize grid generation
!
      CALL insoct2(flon0, flat0, scen, tfile)
!         Call Inoct(n1grid,nagrid,ntay,flon0,flat0,scen,RA,QA,nn  &
!                      ,JUMBLE,WORK)
!
! Define control indices for global arrays
!
      Do j=1,j0-1
        IMIN(j)=i0+1-j
        IMAX(j)=i1+j-1
      End Do
      Do j=j0,j1
        IMIN(j)=1
        IMAX(j)=img
      End Do
      Do j=j1+1,jmg
        IMIN(j)=j-j1+1
        IMAX(j)=img-j+j1
      End Do
      xm(1)=-0.13
      xm(2)=-0.87
      call Otoc1(xm,xe)
      call Otoc2(xm,xe)
      call ctoo(xe,xm,kmap)

!
! Read_in Global Absolute Coordinates (XG,YG,ZG)
!
      dmap=2./(img-1)

      S(1)=-1
      Do l=2,lmg
        S(l)=S(l-1)+dmap
      End Do

      Do j=1,jmg
      Do i=IMIN(j),IMAX(j)
        xm(1)=S(i)
        xm(2)=S(j)
          Call Otoc1(xm,XE)
        XG(i,j,1)=XE(1)
        YG(i,j,1)=XE(2)
        ZG(i,j,1)=XE(3)
        
          Call Otoc2(xm,XE)
        XG(i,j,2)=XE(1)
        YG(i,j,2)=XE(2)
        ZG(i,j,2)=XE(3)

      End Do
      End Do
!
! Define Local Arrays in the Corners of the Northern Hemisphere
! 
! Domain 6
!
      Do j=1,j0 
        Do i=IMIN(j),i0
          Xa(i,j,6)=XG(i,j,1)
          Ya(i,j,6)=YG(i,j,1)
          Za(i,j,6)=ZG(i,j,1)
        End Do 
      End Do 

      Do j=1,j0-1 
        Do i=1,IMIN(j)-1
          Xa(i,j,6)=XG(i1+j-1,j0+1-i,2)
          Ya(i,j,6)=YG(i1+j-1,j0+1-i,2)
          Za(i,j,6)=ZG(i1+j-1,j0+1-i,2)
        End Do 
      End Do 
! 
! Domain 8
!
      Do j=1,j0 
        Do i=i1,IMAX(j)
            il=i-i1+1
          Xa(il,j,8)=XG(i,j,1)
          Ya(il,j,8)=YG(i,j,1)
          Za(il,j,8)=ZG(i,j,1)
        End Do 
      End Do 

      Do j=1,j0-1
        Do i=IMAX(j)+1,img
            il=i-i1+1
          Xa(il,j,8)=XG(i0+1-j,i-i1+1,2)
          Ya(il,j,8)=YG(i0+1-j,i-i1+1,2)
          Za(il,j,8)=ZG(i0+1-j,i-i1+1,2)
        End Do 
      End Do 
! 
! Domain 12
!
      Do j=j1,jmg
        Do i=IMIN(j),i0
            jl=j-j1+1
          Xa(i,jl,12)=XG(i,j,1)
          Ya(i,jl,12)=YG(i,j,1)
          Za(i,jl,12)=ZG(i,j,1)
        End Do 
      End Do 

      Do j=j1+1,jmg
        Do i=1,IMIN(j)-1
            jl=j-j1+1
          Xa(i,jl,12)=XG(i1+img-j,j1+i-1,2)
          Ya(i,jl,12)=YG(i1+img-j,j1+i-1,2)
          Za(i,jl,12)=ZG(i1+img-j,j1+i-1,2)
        End Do 
      End Do 
!
! Domain 14
!
      Do j=j1,jmg
        Do i=i1,IMAX(j)
           il=i-i1+1
           jl=j-j1+1
          Xa(il,jl,14)=XG(i,j,1)
          Ya(il,jl,14)=YG(i,j,1)
          Za(il,jl,14)=ZG(i,j,1)
        End Do 
      End Do 

      Do j=j1+1,jmg
        Do i=IMAX(j)+1,img
            il=i-i1+1
            jl=j-j1+1
          Xa(il,jl,14)=XG(jl,jmg-il+1,2)
          Ya(il,jl,14)=YG(jl,jmg-il+1,2)
          Za(il,jl,14)=ZG(jl,jmg-il+1,2)
        End Do 
      End Do 
!
!  Define All Other Local Arrays 
!  
      Xa(1:iml,1:jml,1)=XG(i0:i1,1:j0,2)
      Ya(1:iml,1:jml,1)=YG(i0:i1,1:j0,2)
      Za(1:iml,1:jml,1)=ZG(i0:i1,1:j0,2)

      Xa(1:iml,1:jml,2)=XG(1:i0,j0:j1,2)
      Ya(1:iml,1:jml,2)=YG(1:i0,j0:j1,2)
      Za(1:iml,1:jml,2)=ZG(1:i0,j0:j1,2)

      Xa(1:iml,1:jml,3)=XG(i0:i1,j0:j1,2)
      Ya(1:iml,1:jml,3)=YG(i0:i1,j0:j1,2)
      Za(1:iml,1:jml,3)=ZG(i0:i1,j0:j1,2)

      Xa(1:iml,1:jml,4)=XG(i1:img,j0:j1,2)
      Ya(1:iml,1:jml,4)=YG(i1:img,j0:j1,2)
      Za(1:iml,1:jml,4)=ZG(i1:img,j0:j1,2)

      Xa(1:iml,1:jml,5)=XG(i0:i1,j1:jmg,2)
      Ya(1:iml,1:jml,5)=YG(i0:i1,j1:jmg,2)
      Za(1:iml,1:jml,5)=ZG(i0:i1,j1:jmg,2)

      Xa(1:iml,1:jml,7)=XG(i0:i1,1:j0,1)
      Ya(1:iml,1:jml,7)=YG(i0:i1,1:j0,1)
      Za(1:iml,1:jml,7)=ZG(i0:i1,1:j0,1)

      Xa(1:iml,1:jml,9)=XG(1:i0,j0:j1,1)
      Ya(1:iml,1:jml,9)=YG(1:i0,j0:j1,1)
      Za(1:iml,1:jml,9)=ZG(1:i0,j0:j1,1)

      Xa(1:iml,1:jml,10)=XG(i0:i1,j0:j1,1)
      Ya(1:iml,1:jml,10)=YG(i0:i1,j0:j1,1)
      Za(1:iml,1:jml,10)=ZG(i0:i1,j0:j1,1)

      Xa(1:iml,1:jml,11)=XG(i1:img,j0:j1,1)
      Ya(1:iml,1:jml,11)=YG(i1:img,j0:j1,1)
      Za(1:iml,1:jml,11)=ZG(i1:img,j0:j1,1)

      Xa(1:iml,1:jml,13)=XG(i0:i1,j1:jmg,1)
      Ya(1:iml,1:jml,13)=YG(i0:i1,j1:jmg,1)
      Za(1:iml,1:jml,13)=ZG(i0:i1,j1:jmg,1)

      Call Bocoh(Xa,iml,jml)
      Call Bocoh(Ya,iml,jml)
      Call Bocoh(Za,iml,jml)

      csal=Cos(alfa)
      snal=Sin(alfa)
      csbt=Cos(beta)
      snbt=Sin(beta)
      csgm=cos(gamm)
      sngm=sin(gamm)

      do n=1,nm
      do i=0,iml+1
      do j=0,iml+1
          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 

!
!  Define XH,YH,ZH
!
       Do n=1,14
         XH(1:im,1:jm,n)=Xa(1:iml:2,1:jml:2,n)
         YH(1:im,1:jm,n)=Ya(1:iml:2,1:jml:2,n)
         ZH(1:im,1:jm,n)=Za(1:iml:2,1:jml:2,n)
       End Do
       Call Bocoh(XH,im,jm)
       Call Bocoh(YH,im,jm)
       Call Bocoh(ZH,im,jm)

       do n=1,nm
       do i=0,im+1
       do j=0,jm+1
        HCOSP (i, j, n) = Sqrt (Xh (i,j, n) **2 + Yh (i,j, n) **2)  
        HSINP (i, j, n) = Zh (i,j, n)  
        HLAM (i, j, n) = Atang1 (Yh (i,j, n), Xh (i,j, n) )  
        HPHI (i, j, n) = ASIN (HSINP (i, j, n) )  
       enddo
       enddo
       enddo
!       print *,'center of tile 3,',hlam((im-1)/2,(jm-1)/2,3)*180./3.14,   &
!              hphi((im-1)/2,(jm-1)/2,3)*180./3.14,gamm
!       stop

!       print *,hlam(9,13,9),yh(9,13,9),xh(9,13,9)
!       stop

!
!  Define XV,YV,ZV
!
       Do n=1,nm
         XV(0:im,0:jm,n)=Xa(0:iml+1:2,0:jml+1:2,n)
         YV(0:im,0:jm,n)=Ya(0:iml+1:2,0:jml+1:2,n)
         ZV(0:im,0:jm,n)=Za(0:iml+1:2,0:jml+1:2,n)
       End Do

       do n=1,nm
       do i=0,im
       do j=0,jm
        VCOSP (i, j, n) = Sqrt (Xv (i,j, n) **2 + Yv (i,j, n) **2)  
        VSINP (i, j, n) = Zv (i,j, n)  
        VLAM (i, j, n) = Atang1 (Yv (i,j, n), Xv (i,j, n) )  
        VPHI (i, j, n) = ASIN (VSINP (i, j, n) )  
       enddo
       enddo
       enddo

!
!  Correction of coefficients in the corners
!
      IMOD(:,:,:)=1

      IMOD(1,1, 1)=0
      IMOD(1,1, 2)=0
      IMOD(1,1, 7)=0
      IMOD(1,1, 8)=0
      IMOD(1,1, 9)=0
      IMOD(1,1,12)=0
      IMOD(im,1,6 )=0
      IMOD(im,1,7 )=0
      IMOD(im,1,11)=0
      IMOD(im,1,14)=0
      IMOD(im,1, 1)=0
      IMOD(1,jm, 2)=0
      IMOD(im,1, 4)=0
      IMOD(1,jm, 5)=0
      IMOD(1,jm, 6)=0
      IMOD(1,jm, 9)=0
      IMOD(1,jm,13)=0
      IMOD(1,jm,14)=0
      IMOD(im,jm, 4)=0
      IMOD(im,jm, 5)=0
      IMOD(im,jm, 8)=0
      IMOD(im,jm,11)=0
      IMOD(im,jm,12)=0
      IMOD(im,jm,13)=0

      LV(:,:,:)=.true.

      LV(0 ,0 ,1)=.false.
      LV(im,0 ,1)=.false.
      LV(0 ,0 ,2)=.false.
      LV(0 ,jm,2)=.false.
      LV(0 ,jm,5)=.false.
      LV(im,jm,5)=.false.
      LV(im,jm,4)=.false.
      LV(im,0 ,4)=.false.

      LV(0 ,0 ,8)=.false.
      LV(im,jm,8)=.false.
      LV(im,0 ,14)=.false.
      LV(0 ,jm,14)=.false.
      LV(im,jm,12)=.false.
      LV(0 ,0 ,12)=.false.
      LV(im,0 ,6 )=.false.
      LV(0 ,jm,6 )=.false.

      LV(0 ,0 ,7)=.false.
      LV(im,0 ,7)=.false.
      LV(0 ,0 ,9)=.false.
      LV(0 ,jm,9)=.false.
      LV(0 ,jm,13)=.false.
      LV(im,jm,13)=.false.
      LV(im,jm,11)=.false.
      LV(im,0 ,11)=.false.

      Do n = 1, nm  
      Do j = 0, jm  
      Do i = 0, im  
        DXVDX (i,j,n) =(XH(i+1,j,n)-XH(i,j+1,n))*mfac 
        DYVDX (i,j,n) =(YH(i+1,j,n)-YH(i,j+1,n))*mfac 
        DZVDX (i,j,n) =(ZH(i+1,j,n)-ZH(i,j+1,n))*mfac 
        DXVDY (i,j,n) =(XH(i+1,j+1,n)-XH(i,j,n))*mfac 
        DYVDY (i,j,n) =(YH(i+1,j+1,n)-YH(i,j,n))*mfac 
        DZVDY (i,j,n) =(ZH(i+1,j+1,n)-ZH(i,j,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)

        If (LV(i, j,n) ) Then  
            SQV (i, j,n) = Sqrt (qs11 * qs22 - qs12 * qs12)  
            RSQV (i, j,n) = 1. / SQV (i, j,n)  
            Q11 (i, j, n) = qs22 * RSQV (i, j,n) **2  
            Q22 (i, j, n) = qs11 * RSQV (i, j,n) **2  
            Q12 (i, j, n) = - qs12 * RSQV (i, j,n) **2  
        else
            SQV (i, j,n)  =  0. 
            RSQV (i, j,n) = 0.
            Q11 (i, j, n) = 0.
            Q22 (i, j, n) = 0.
            Q12 (i, j, n) = 0.
       EndIf  
!       print *,i,j,n,sqv(i,j,n),qs11,qs22,qs12
      EndDo  
      EndDo  
      EndDo  

      call bocov1(sqv,im,jm)
      call bocov2(q11,q22,im,jm)
      call bocov3(q12,im,jm)
      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)
!DRAGAN      call qintc_2(xh,yh,zh,xv,yv,zv,im,jm,nm,mfac)

!
! Grid components: H-points
!
       sqhc=SQV(1,1,7)*2./3.

      do n=1,nm
      Do j = 1, jm  
      Do i = 1, im  
        imd=IMOD(i,j,n)
        dxhdx=(XV(i,j-1,n)-XV(i-1,j,n))*mfac 
        dyhdx=(YV(i,j-1,n)-YV(i-1,j,n))*mfac 
        dzhdx=(ZV(i,j-1,n)-ZV(i-1,j,n))*mfac 
        dxhdy=(XV(i,j,n)-XV(i-1,j-1,n))*mfac 
        dyhdy=(YV(i,j,n)-YV(i-1,j-1,n))*mfac 
        dzhdy=(ZV(i,j,n)-ZV(i-1,j-1,n))*mfac 
        qs11 = dxhdx**2 + dyhdx**2 + dzhdx**2  
        qs22 = dxhdy**2 + dyhdy**2 + dzhdy**2  
        qs12 = dxhdx * dxhdy + dyhdx * dyhdy + dzhdx * dzhdy  
!        SQH(i,j,n)=Sqrt(qs11*qs22-qs12*qs12)*imd  &
!                   +    sqhc*(1-imd)
       sqh(i,j,n)=0.25*(sqv(i-1,j-1,n)+sqv(i,j,n)+sqv(i,j-1,n)  &
	     +sqv(i-1,j,n))
!        print *,i,j,n,sqh(i,j,n)
       rarea(i,j,n)=1./(2.*sqh(i,j,n)*a*a)
       ds(i,j,n,1)=0.5*sqrt(qs11)*a
       ds(i,j,n,3)=0.5*sqrt(qs22)*a
      EndDo  
      EndDo  
      EndDo  

      call bocoh(sqh,im,jm)
      call bocoh(rarea,im,jm)

      do n=1,nm
      Do j = 1, jm  
      Do i = 1, im  
        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  
        ds(i,j,n,2)=0.25*sqrt(qs11)*a
        ds(i,j,n,4)=0.25*sqrt(qs22)*a
      EndDo  
      EndDo  
      EndDo  

      do n=1,nm
      do i=1,im1
      do j=1,jm1
      do k=1,4
        if(k.eq.1)then
          it=i-1
          jt=j
        else if(k.eq.2)then
          it=i+1
          jt=j
        else if(k.eq.3)then
          it=i
          jt=j-1
        else
          it=i
          jt=j+1
        endif

        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
         enddo

      do n=1,nm
      do i=0,im+1
      do j=0,jm+1
       if(sqh(i,j,n).ne.0)then
        rsqh(i,j,n)=1./sqh(i,j,n)
       else
        rsqh(i,j,n)=0.
       endif
      enddo
      enddo
      enddo

      do n=1,nm
      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))
!	print *,i,j,n,qbv11(i,j,n),qbv12(i,j,n),qbv22(i,j,n)
      enddo
      enddo
      enddo

!
! Outfile for Absolute Coordintes
!

       print *,"hsinp(1,1,1),",hsinp(1,1,1)
         list1 = 11  
	 nres=1
       write (cres, 900) nres  
  900 format(i1.1)  

      cfile1 = "../data_in/grid/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
      Close (list1)  

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

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

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


!test
       print * , 'Wrote ', cfile1, ' file'
       print *,'center of tile 3,',hlam((im-1)/2,(jm-1)/2,3),   &
              hphi((im-1)/2,(jm-1)/2,3)
!test
!-----------------------------------------------------------------------
!
! Prepare Graphical Output
!
!
! Define GLON and GLAT
!

      pihlf=2.*Atan(1.)
      pi=2.*pihlf
      rad=pi/180.

      dlat=   pi/(jgm-1)
      dlon=2 *pi/igm

      Do i=1,igm
        xlon=(i-1)*dlon
      Do j=1,jgm
        GLON(i,j)=xlon
      End Do
      End Do

      Do j=1,jgm
       ylat = - pihlf + (j - 1) * dlat  
      Do i=1,igm
        GLAT(i,j)=ylat
      End Do
      End Do

      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
        XE(1)=x0
        XE(2)=y0
        XE(3)=z0
        Call Ctoo(XE,xm,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
!	print *,i,j,xmap,ymap,kmap
!	stop
      End Do
      End Do

!
! 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  
!print *,i,j,n,u2us(i,j,n),v2us(i,j,n),u2vs(i,j,n),v2vs(i,j,n)
       enddo  
       enddo  
       enddo  
!-----------------------------------------------------------------------
!
! Outfile for graphics
!

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

       cfile2 = "../data_in/grid/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'  
!       open(11,file="ocdt.dat",form="unformatted")
!       write(11)xa
!       write(11)ya
!       write(11)za
!       close(11)
!test
!-----------------------------------------------------------------------
       Stop  
       END program metrcs
