C -*- Mode: Fortran; -*-
cc-----------------------------------------------------------------
c     Ravi Samtaney
c     Copyright 2004
c     Princeton Plasma Physics Laboratory
c     All Rights Reserved
c-----------------------------------------------------------------
c     $Log$

c
c-----------------------------------------------------------------
      subroutine SetupDomain(Rn,Zn)
      use mesh
      use mesh_common
      double precision:: Zn(ixlo:ixhi+1,iylo:iyhi+1)
      double precision:: Rn(ixlo:ixhi+1,iylo:iyhi+1)
      call ShapedTorusESC(Rn, Zn)
      return
      end
c
c----------------------------------------------------------------------
c     Read ESC equilibrium file for initial grid
      subroutine ShapedTorusESC(Rn,Zn)
      use mesh
      use mesh_common
      use iounits
      use profiles
      implicit none
      double precision:: Zn(ixlo:ixhi+1,iylo:iyhi+1)
      double precision:: Rn(ixlo:ixhi+1,iylo:iyhi+1)
      integer:: i,j, k
      integer:: ig, jg
      double precision:: twopi
      double precision:: RMajor, RMinor
      double precision, allocatable:: Rnode(:,:), Znode(:,:)
      double precision, allocatable:: RTmp(:,:), ZTmp(:,:)
      integer:: npsi, ntheta
c      double precision:: rma, zma
c
c      open(16,file='esc.txt',form='formatted')
c      open(16,file='eq3.txt',form='formatted')
      open(16,file='eq3r.txt',form='formatted')
      read(16,*) ntheta
      read(16,*) npsi       
c      write(6,*) 'ntheta, npsi=',ntheta,npsi

c
        allocate(psi(npsi))
        allocate(psibig(npsi))
        allocate(pprof(npsi))
        allocate(qprof(npsi))
        allocate(gprof(npsi))
        allocate(theta(ntheta))
        allocate(RNode(-1:npsi+2,-1:ntheta+2))
        allocate(ZNode(-1:npsi+2,-1:ntheta+2))
        Znode=0.D0
        Rnode=0.D0
        allocate(RTmp(ntheta,npsi))
        allocate(ZTmp(ntheta,npsi))
c
        read(16,*) rma
        read(16,*) zma
c        write(6,*) 'Rma, Zma',rma, zma
        
        read(16,*) psi
        do i=1,npsi,1
c           read(16,*) psi(i)
c           write(6,*) 'psi=',i,psi(i)
        enddo
c
        read(16,*) theta
        do i=1,ntheta,1
c           write(6,*) 'theta=',i,theta(i)
        enddo
c
        read(16,*) psibig
        do i=1,npsi,1
c           write(6,*) 'psibig=',i,psibig(i)
        enddo
c
        read(16,*) pprof
        do i=1,npsi,1
c           write(6,*) 'pprof=',i,pprof(i)
        enddo
c
        read(16,*) gprof
        do i=1,npsi,1
c           write(6,*) 'gprof=',i,gprof(i)
        enddo
c
        read(16,*) qprof
        do i=1,npsi,1
c           write(6,*) 'qprof=',i,qprof(i)
        enddo
c
        read(16,*) RTmp
        do i=1,npsi,1
           do j=1,ntheta-1,1
              Rnode(i,j)=RTmp(j+1,i)
c              write(6,*) 'RNode',i,j,Rnode(i,j)
           enddo
        enddo
c
        read(16,*) ZTmp
        close(16)
c
        do i=1,npsi,1
           do j=1,ntheta-1,1
              Znode(i,j)=ZTmp(j+1,i)
c              write(6,*) 'ZNode',i,j,Znode(i,j)
c              write(81,*) i,j,Rnode(i,j), Znode(i,j)
           enddo
c           write(81,*)
        enddo
c
c     
c     Periodic BC in eta
        Znode(-1:npsi+2,ntheta)=Znode(-1:npsi+2,1)
        Znode(-1:npsi+2,0)=Znode(-1:npsi+2,ntheta-1)
        Znode(-1:npsi+2,-1)=Znode(-1:npsi+2,ntheta-2)
        Znode(-1:npsi+2,ntheta+1)=Znode(-1:npsi+2,2)
        Znode(-1:npsi+2,ntheta+2)=Znode(-1:npsi+2,3)
c
        Rnode(-1:npsi+2,ntheta)=Rnode(-1:npsi+2,1)
        Rnode(-1:npsi+2,0)=Rnode(-1:npsi+2,ntheta-1)
        Rnode(-1:npsi+2,-1)=Rnode(-1:npsi+2,ntheta-2)
        Rnode(-1:npsi+2,ntheta+1)=Rnode(-1:npsi+2,2)
        Rnode(-1:npsi+2,ntheta+2)=Rnode(-1:npsi+2,3)
c
        
        twopi=8.D0*datan(1.D0)
c
        dxi=1.D0
        deta=1.D0
c
c	Get the global indices
        call iglobal(iprocx,1,ig)
        call jglobal(iprocy,1,jg)
        ig=ig-1
        jg=jg-1
c     
      do j=iylo,iyhi+1,1
         do i=ixlo,ixhi+1,1
            Zn(i,j)=Znode(i+ig,j+jg)
            Rn(i,j)=Rnode(i+ig,j+jg)
         enddo
      enddo
c
      do j=iylo,iyhi,1
         do i=ixlo,ixhi,1
            Zc(i,j)=0.25D0*(Znode(i+ig,j+jg)+
     &                      Znode(i+ig+1,j+jg)+
     &                      Znode(i+ig,j+jg+1)+
     &                      Znode(i+ig+1,j+jg+1))
            Rc(i,j)=0.25D0*(Rnode(i+ig,j+jg)+
     &                      Rnode(i+ig+1,j+jg)+
     &                      Rnode(i+ig,j+jg+1)+
     &                      Rnode(i+ig+1,j+jg+1))
c            write(71,*) i,j,Rc(i,j), Zc(i,j)
         enddo
      enddo
c
c$$$      i=1
c$$$      do j=iylo,iyhi,1
c$$$            Zc(i,j)=(Zn(i,j)+
c$$$     &                      Zn(i+1,j)+
c$$$     &                      Zn(i+1,j+1))/3.D0
c$$$            Rc(i,j)=(Rn(i,j)+
c$$$     &                      Rn(i+1,j)+
c$$$     &                      Rn(i+1,j+1))/3.D0
c$$$      enddo
c
c
      do j=iylo,iyhi,1
         do i=ixlo,ixhi+1,1
            dRdEta(i,j)=(RNode(i+ig,j+jg+1)-RNode(i+ig,j+jg))/deta
            dZdEta(i,j)=(ZNode(i+ig,j+jg+1)-ZNode(i+ig,j+jg))/deta
         enddo
      enddo
c
c     Length of xi face
      do j=iylo,iyhi,1
         do i=ixlo,ixhi+1,1
            dlXi(i,j)=dsqrt(dRdEta(i,j)**2+dZdEta(i,j)**2)
            RXi(i,j)=0.5D0*(RNode(i+ig,j+jg)+RNode(i+ig,j+jg+1))
         enddo
      enddo
c
      do j=iylo,iyhi+1,1
         do i=ixlo,ixhi,1
            dRdXi(i,j)=(RNode(i+ig+1,j+jg)-RNode(i+ig,j+jg))/dxi
            dZdXi(i,j)=(ZNode(i+ig+1,j+jg)-ZNode(i+ig,j+jg))/dxi
c            write(91,*) i,j, dRdXi(i,j), dZdXi(i,j)
         enddo
c         write(91,*)
      enddo
c
c     Length of eta face
      do j=iylo,iyhi+1,1
         do i=ixlo,ixhi,1
            dlEta(i,j)=dsqrt(dRdXi(i,j)**2+dZdXi(i,j)**2)
            REta(i,j)=0.5D0*(RNode(i+ig,j+jg)+RNode(i+ig+1,j+jg))
c            write(6,*) i+ig,j+jg, dRdXi(i,j), dZdXi(i,j), dlEta(i,j)
         enddo
      enddo
c     Inverse Jacobian = R_xi Z_eta - R_eta Z_xi
      do j=iylo,iyhi,1
         do i=ixlo,ixhi,1
            JacI(i,j)=0.25D0*(
     &           (dRdXi(i,j)+dRdXi(i,j+1))*
     &           (dZdEta(i,j)+dZdEta(i+1,j))-
     &           (dRdEta(i,j)+dRdEta(i+1,j))*
     &           (dZdXi(i,j)+dZdXi(i,j+1)))
            dRdXiC(i,j)=0.5D0*(dRdXi(i,j)+dRdXi(i,j+1))
            dZdXiC(i,j)=0.5D0*(dZdXi(i,j)+dZdXi(i,j+1))
            dRdEtaC(i,j)=0.5D0*(dRdEta(i,j)+dRdEta(i+1,j))
            dZdEtaC(i,j)=0.5D0*(dZdEta(i,j)+dZdEta(i+1,j))
            RJC(i,j)=RC(i,j)*Jaci(i,j)
c            write(51,*) i,j,dRdXiC(i,j), dZdXiC(i,j),
c     &           dRdEtaC(i,j), dZdEtaC(i,j)
         enddo
c         write(93,*) 
      enddo
c
        if(iprocx.eq.1.and.iprocy.eq.1.and.iprocz.eq.1) then
           write(ihis,*) 'MESH SIZE=',nx,ny
           write(ihis,*) 'LOCAL MESH SIZE=',nxlocal,nylocal
           write(ihis,*) 'Bounds xi=',ixlo,ixhi
           write(ihis,*) 'Bounds eta=',iylo,iyhi
           write(ihis,*) 
c
           write(ihis,*) 'DXi=',dxi
           write(ihis,*) 'DEta=',deta
c
           write(ihis,*) 
           write(ihis,*) 'MAX R=', maxval(RNode(1:npsi,1:ntheta))
           write(ihis,*) 'MIN R=', minval(RNode(1:npsi,1:ntheta))
           write(ihis,*) 'MAX Z=', maxval(ZNode(1:npsi,1:ntheta))
           write(ihis,*) 'MIN Z=', minval(ZNode(1:npsi,1:ntheta))
           write(ihis,*) 
        endif
c
c     Deallocate tmp variables
      deallocate(RNode)
      deallocate(ZNode)
      deallocate(RTmp)
      deallocate(ZTmp)
      return
      end
c     

c     



c----------------------------------------------------------------------
c     Read ESC equilibrium file for initial grid
      subroutine ComputeMetrics(Rn,Zn)
      use mesh
      use mesh_common
      use iounits
      use profiles
      implicit none
      double precision:: Zn(ixlo:ixhi+1,iylo:iyhi+1)
      double precision:: Rn(ixlo:ixhi+1,iylo:iyhi+1)
      integer:: i,j
c
      do j=iylo,iyhi,1
         do i=ixlo,ixhi,1
            Zc(i,j)=0.25D0*(Zn(i,j)+
     &                      Zn(i+1,j)+
     &                      Zn(i,j+1)+
     &                      Zn(i+1,j+1))
            Rc(i,j)=0.25D0*(Rn(i,j)+
     &                      Rn(i+1,j)+
     &                      Rn(i,j+1)+
     &                      Rn(i+1,j+1))
c            write(71,*) i,j,Rc(i,j), Zc(i,j)
         enddo
      enddo
c
c$$$      i=1
c$$$      do j=iylo,iyhi,1
c$$$            Zc(i,j)=(Zn(i,j)+
c$$$     &                      Zn(i+1,j)+
c$$$     &                      Zn(i+1,j+1))/3.D0
c$$$            Rc(i,j)=(Rn(i,j)+
c$$$     &                      Rn(i+1,j)+
c$$$     &                      Rn(i+1,j+1))/3.D0
c$$$      enddo
c
c
      do j=iylo,iyhi,1
         do i=ixlo,ixhi+1,1
            dRdEta(i,j)=(RN(i,j+1)-RN(i,j))/deta
            dZdEta(i,j)=(ZN(i,j+1)-ZN(i,j))/deta
         enddo
      enddo
c
c     Length of xi face
      do j=iylo,iyhi,1
         do i=ixlo,ixhi+1,1
            dlXi(i,j)=dsqrt(dRdEta(i,j)**2+dZdEta(i,j)**2)
            RXi(i,j)=0.5D0*(RN(i,j)+RN(i,j+1))
         enddo
      enddo
c
      do j=iylo,iyhi+1,1
         do i=ixlo,ixhi,1
            dRdXi(i,j)=(RN(i+1,j)-RN(i,j))/dxi
            dZdXi(i,j)=(ZN(i+1,j)-ZN(i,j))/dxi
c            write(91,*) i,j, dRdXi(i,j), dZdXi(i,j)
         enddo
c         write(91,*)
      enddo
c
c     Length of eta face
      do j=iylo,iyhi+1,1
         do i=ixlo,ixhi,1
            dlEta(i,j)=dsqrt(dRdXi(i,j)**2+dZdXi(i,j)**2)
            REta(i,j)=0.5D0*(RN(i,j)+RN(i+1,j))
c            write(6,*) i,j, dRdXi(i,j), dZdXi(i,j), dlEta(i,j)
         enddo
      enddo
c     Inverse Jacobian = R_xi Z_eta - R_eta Z_xi
      do j=iylo,iyhi,1
         do i=ixlo,ixhi,1
            JacI(i,j)=0.25D0*(
     &           (dRdXi(i,j)+dRdXi(i,j+1))*
     &           (dZdEta(i,j)+dZdEta(i+1,j))-
     &           (dRdEta(i,j)+dRdEta(i+1,j))*
     &           (dZdXi(i,j)+dZdXi(i,j+1)))
            dRdXiC(i,j)=0.5D0*(dRdXi(i,j)+dRdXi(i,j+1))
            dZdXiC(i,j)=0.5D0*(dZdXi(i,j)+dZdXi(i,j+1))
            dRdEtaC(i,j)=0.5D0*(dRdEta(i,j)+dRdEta(i+1,j))
            dZdEtaC(i,j)=0.5D0*(dZdEta(i,j)+dZdEta(i+1,j))
            RJC(i,j)=RC(i,j)*Jaci(i,j)
         enddo
c         write(93,*) 
      enddo
c
c$$$        if(iprocx.eq.1.and.iprocy.eq.1.and.iprocz.eq.1) then
c$$$           write(ihis,*) 'MESH SIZE=',nx,ny
c$$$           write(ihis,*) 'LOCAL MESH SIZE=',nxlocal,nylocal
c$$$           write(ihis,*) 'Bounds xi=',ixlo,ixhi
c$$$           write(ihis,*) 'Bounds eta=',iylo,iyhi
c$$$           write(ihis,*) 
c$$$c
c$$$           write(ihis,*) 'DXi=',dxi
c$$$           write(ihis,*) 'DEta=',deta
c$$$c
c$$$           write(ihis,*) 
c$$$           write(ihis,*) 'MAX R=', maxval(RN(1:nxlocal+1,1:nylocal+1))
c$$$           write(ihis,*) 'MIN R=', minval(RN(1:nxlocal+1,1:nylocal+1))
c$$$           write(ihis,*) 'MAX Z=', maxval(ZN(1:nxlocal+1,1:nylocal+1))
c$$$           write(ihis,*) 'MIN Z=', minval(ZN(1:nxlocal+1,1:nylocal+1))
c$$$           write(ihis,*) 
c$$$        endif
c
      return
      end
c     

c