C -*- Mode: Fortran; -*-
c-----------------------------------------------------------------
c     Ravi Samtaney
c     Copyright 2004
c     Princeton Plasma Physics Laboratory
c     All Rights Reserved
c-----------------------------------------------------------------
c     $Log: mhddriver.F,v $
c     Revision 1.4  2005/10/11 22:05:14  samtaney
c     Added more flags.
c
c     Revision 1.3  2005/10/02 19:09:32  samtaney
c     Put in several options.
c
c     Revision 1.2  2005/09/23 19:01:51  samtaney
c     Also writing tempAblation to time.txt.
c
c     Revision 1.1.1.1  2005/09/12 17:28:10  samtaney
c     Original source.
c
c     Revision 1.8  2005/08/12 19:59:46  samtaney
c     Added implicit none. New subroutine InitIO to open diagnostic files etc.
c
c     Revision 1.7  2005/08/11 20:01:38  samtaney
c     Allocating phi from izlo to izhi. Other very minor changes.
c
c     Revision 1.6  2004/12/09 14:39:12  samtaney
c     Added ttot as arg to mhdsolve. This is passed to the pellet source routine.
c
c     Revision 1.5  2004/12/04 15:45:40  samtaney
c     Added capability to include eqsrc. Also added calls to write 2D AVS files for
c     2D runs.
c
c     Revision 1.4  2004/09/24 19:35:31  samtaney
c     Added call to AVS output after initial conditions are specified.
c
c     Revision 1.3  2004/09/15 15:56:39  samtaney
c     Initialized hasxlo etc to zero. Also added conditional compilation for
c     ETAPERIODIC.
c
c     Revision 1.2  2004/06/18 22:30:14  samtaney
c     Intermediate checkin.
c
c     Revision 1.1.1.1  2004/05/28 19:33:27  samtaney
c     Original source - copied from 3D Unsplit code.
c
c-----------------------------------------------------------------
c
       subroutine MHDMain
       use mesh
       use mesh_common
c
      implicit none
#ifdef DYNAMIC
       ixlo=1-nghost
       iylo=1-nghost
       ixhi=nxlocal+nghost
       iyhi=nylocal+nghost
       hasxlo=0
       hasxhi=0
       hasylo=0
       hasyhi=0
       if(iprocx.eq.1)then
          ixlo=1
          hasxlo=1
       endif
       if(iprocx.eq.XPROCS)then
          ixhi=nxlocal
          hasxhi=1
       endif
#ifndef ETAPERIODIC
       if(iprocy.eq.1) then
          iylo=1
          hasylo=1
       endif
       if(iprocy.eq.YPROCS)then
          iyhi=nylocal
          hasyhi=1
       endif
#endif
c$$$       ilo=1
c$$$       ihi=nxlocal+1
c$$$       jlo=1
c$$$       jhi=nylocal+1
c$$$       klo=1
c$$$       khi=nzlocal+1
c     Don't need to do the following in the phi direction where we have
c     periodic boundary conditions
c$$$       if(iprocx.eq.1) ilo=2
c$$$       if(iprocy.eq.1) jlo=2
c$$$       if(iprocx.eq.XPROCS) ihi=nxlocal
c$$$       if(iprocy.eq.YPROCS) jhi=nylocal
c     
       nxlsize=nxlocal
       nylsize=nylocal
#endif
c     Allocate mesh
c     Face metrics (Is metrics the correct word for these?)
       allocate(dRdXi(ixlo:ixhi,iylo:iyhi+1))
       allocate(dZdXi(ixlo:ixhi,iylo:iyhi+1))
       allocate(dRdEta(ixlo:ixhi+1,iylo:iyhi))
       allocate(dZdEta(ixlo:ixhi+1,iylo:iyhi))
c
       allocate(dlEta(ixlo:ixhi,iylo:iyhi+1))
       allocate(dlXi(ixlo:ixhi+1,iylo:iyhi))
c
       allocate(Rc(ixlo:ixhi,iylo:iyhi))
       allocate(RJc(ixlo:ixhi,iylo:iyhi))
c
       allocate(REta(ixlo:ixhi,iylo:iyhi+1))
       allocate(RXi(ixlo:ixhi+1,iylo:iyhi))
c
       allocate(Zc(ixlo:ixhi,iylo:iyhi))
c
c       allocate(Zn(ixlo:ixhi+1,iylo:iyhi+1))
c       allocate(Rn(ixlo:ixhi+1,iylo:iyhi+1))
c
       allocate(JacI(ixlo:ixhi,iylo:iyhi))
c
       allocate(dRdXiC(ixlo:ixhi,iylo:iyhi))
       allocate(dRdEtaC(ixlo:ixhi,iylo:iyhi))
       allocate(dZdXiC(ixlo:ixhi,iylo:iyhi))
       allocate(dZdEtaC(ixlo:ixhi,iylo:iyhi))
c
c
       call MHDDriver
       end
c     
c-----------------------------------------------------------------
      subroutine MHDDriver
c     
      use mesh
      use mesh_common
c     use time
      use iounits
#ifdef PARALLEL
      use mpistuff
#endif
      use OptionFlags
      implicit none
c     
c=======Declarations=========
c      double precision,allocatable:: ux(:,:,:,:)
        double precision:: Zn(ixlo:ixhi+1,iylo:iyhi+1)
        double precision:: Rn(ixlo:ixhi+1,iylo:iyhi+1)
      double precision:: dt,ttot
c     
      integer:: istart, iter, maxiter, lastiter
      integer:: output_flag, binary_flag
      integer:: ndump, ndiag, nprof
      integer:: dump,new
      integer:: ipar
      integer:: tvdFlag

      integer:: timestep
c     
      namelist/inparams/maxiter,new, output_flag, binary_flag
      namelist/dumps/ndump,ndiag,nprof
      namelist/options/idealMHDFlag, fluxType, nonConFlag,
     &     cTransFlag,emCTFlag,eqsrcFlag,qModelFlag,
     &     tempBCFlag,launchType,staticBFlag
c      namelist/domain/RMajor, RMinor
c     
      open(16,file='mhd.inp',form='formatted')
      read(16,inparams)
c      write(6,inparams)
      read(16,dumps)
c      write(6,dumps)
c      read(16,domain) 
c      write(6,domain)
      read(16,options)
      close(16)
c     
c     
      lastiter=0
      istart=1
      dump=1
c     
c     
      call SetupDomain(Rn,Zn)
c      write(6,*) 'After domain setup'
c     
c     
#ifdef PARALLEL
      call MPI_Barrier(comm3d, ierr)
#endif
c     

#ifndef PARALLEL
c      write(ihis,*) 'AVS File at',timeStep, ttot
c$$$#ifdef TWO_D
c$$$      Call WriteAVSFile2D(ux,0)
c$$$#else
c$$$      Call WriteAVSFile(ux,0)
c$$$#endif
#else
      if(iprocx.eq.1.and.iprocy.eq.1.and.iprocz.eq.1) then
         write(ihis,*) 'Binary output at',istart
      endif
c      call WriteBinaryFileParallel(ux,istart)
#endif
c      do iter=istart,istart+maxiter-1,1
         call NewtonSolve
c     
c         if(iprocx.eq.1.and.iprocy.eq.1.and.iprocz.eq.1) then
c            write(6,*) 'iter=',iter
c            write(ihis,*) 'iter=',iter
c         endif
c     
c      enddo
c     
c     
      return
      end
c     
c-----------------------------------------------------------------------