C -*- Mode: Fortran; -*-
c-----------------------------------------------------------------
c       Ravi Samtaney
c       Copyright 2004
c       Princeton Plasma Physics Laboratory
c       All Rights Reserved
c-----------------------------------------------------------------
c       $Log$
c
c-----------------------------------------------------------------
c Subroutine to set up the initial conditions
c Read output from Jsolver which gives an equilibrium initial
c condition
c---------------------------------------------------------------------------
c
c
      subroutine SetParameters
      use mesh
      use mesh_common
      use properties
      use BoundaryPressTemp
      use iounits
      use OptionFlags
      use Profiles
      use MagAxisToroidalField
      implicit none
c
c
	character*24 fdate
        double precision:: alpha, pi
c        
        integer:: i,j,ig,jg
        double precision:: jactmp
c        double precision:: Rp, Zp
c        double precision:: xp0,yp0,zp0,xp1,yp1,zp1, chk
        double precision:: beta, bxi
c

c
c     Scale factor
        alpha=0.375D0
c        alpha=1.D0
c
        pi=4.D0*datan(1.D0)
c
        if(iprocx.eq.1.and.iprocy.eq.1.and.iprocz.eq.1) then
           write(ihis,*) 'Scale factor for B Field', alpha
           write(ihis,*) 'Scale factor for pressure', alpha*alpha
           write(ihis,*) 
        endif
c     Magnetic Axis Toroidal Field
        Btaxis=gprof(1)*alpha
c        Btaxis=0.D0
c        write(6,*) 'Magnetic Axis Toroidal Field=', Btaxis
        if(iprocx.eq.1.and.iprocy.eq.1.and.iprocz.eq.1) then
           write(ihis,*) 'Magnetic Axis Toroidal Field x R=', Btaxis
c           write(ihis,*) 'Magnetic Axis Toroidal Field=', Btaxis/rma
        endif
c
c
        psifactor=(psibig(nx+1)-psibig(1))/(nx)
        psifactor=psifactor*0.03698/(psibig(nx)-psibig(1))
c
c     Reference quantities
        wmol=2.D0
        RUniversal=8314.5D0
c     rgas=RUniversal/wmol
c     This is a non-dim gas constant - Eos is p=rho T
        rgas=1.D0
c     Permeability of free space in Henry/m or kg-m/Coulomb^2
        mu0=4.D0*pi * 1.D-07
c     Reference Magnetic field in Tesla (kg/s/coulomb)
c        Bref=0.25D0
        Bref=1.D0
c     CDXU
c        Bref=0.375D0
c     DIII-D
c        Bref=2.0D0
c     ITER
c        Bref=5.0D0
c      Bref=BToroidalMax
c
c     Alternative pref is Bref^2/mu0 (in N/m^2)
      pref=Bref**2/mu0
c     Reference no. of molecules /m^3
      n0ref=1.5D19
c     Mass of each proton in Kg
      Mi=1.6605 D-27
c     Reference density rhoRef=n0ref*Mi*wmol
c      rhoRef=1.6605*wmol*1.D-08
      rhoRef=Mi*wmol*n0Ref
c     Reference temperature in K
      TRef=pref/(rhoRef*RUniversal/wmol)
c     Convert K to eV
      K2eV=1.1604D04
c     Reference temperature in eV
      TRefeV=Tref/K2eV
c
c     Reference velocity is the alfven velocity defined 
c        as Bref/sqrt(rhoRef*mu0)
      velRef=Bref/dsqrt(rhoRef*mu0)
c     
c     Reference length scale is 1m
      lref=1.D0
c
c     Reference time scale tauRef
      tauRef=lref/velRef
c
      pambient=0.03D0/mu0*alpha*alpha
c
      if(iprocx.eq.1.and.iprocy.eq.1.and.iprocz.eq.1) then
         write(ihis,*) 'Reference Quantities'
         write(ihis,*) 'Density (kg/m^3)', rhoRef
         write(ihis,*) 'Number density (/m^3)',n0Ref
         write(ihis,*) 'Pressure (Pa)', pRef
         write(ihis,*) 'Temperature (K)', TRef
         write(ihis,*) 'Temperature (eV)', TRefeV
         write(ihis,*) 'Magnetic field (T)', Bref
         write(ihis,*) 'Velocity (m/s)', velRef
         write(ihis,*) 'Energy  (J/m^3)', rhoRef*velRef**2
         write(ihis,*) 'Length (m)', lref
         write(ihis,*) 'time (s)', tauRef
         write(ihis,*) 
         write(ihis,*) 'Ambient pressure (Pa)', pambient
         write(ihis,*) 'PSIFACTOR=', psifactor
         write(ihis,*) 
         write(ihis,*) 'Plasma gamma', gamma
         write(ihis,*) 'Universal gas constant (J/kmol/K)', RUniversal
         write(ihis,*) 'Molecular weight', wmol
         write(ihis,*) 'Gas constant (nondim)', rgas
         write(ihis,*) 'EoS: p=rho T'
         write(ihis,*) 
      endif
c
        pbdry=(pprof(nx+1)/mu0/pRef)
     &                *alpha*alpha
     &                +pambient/pRef
c     Taking density as unity
        TempBdry=pbdry
c
        return
        end
c-----------------------------------------------------------------------
      subroutine ComputeConserved(ux)
      use mesh
      use mesh_common
      use properties
      use BoundaryPressTemp
      use iounits
      use OptionFlags
      use MagAxisToroidalField
      use Profiles
      implicit none
      double precision:: ux(IXLO:IXHI,IYLO:IYHI,nvar)
c      double precision:: Zn(ixlo:ixhi+1,iylo:iyhi+1)
c      double precision:: Rn(ixlo:ixhi+1,iylo:iyhi+1)
c
c
c=======Declarations=========
        double precision:: rho,press, Bi, Bj, Bk, u,v,w
c
	character*24 fdate
        double precision:: alpha, pi
c        
        integer:: i,j,ig,jg
        double precision::  jactmp
c        double precision:: Rp, Zp
c        double precision:: xp0,yp0,zp0,xp1,yp1,zp1, chk
        double precision:: beta, bxi
c

c
c     Scale factor
        alpha=0.375D0
c        alpha=1.D0
c
        pi=4.D0*datan(1.D0)
c
c$$$        if(iprocx.eq.1.and.iprocy.eq.1.and.iprocz.eq.1) then
c$$$           write(ihis,*) 'Scale factor for B Field', alpha
c$$$           write(ihis,*) 'Scale factor for pressure', alpha*alpha
c$$$           write(ihis,*) 
c$$$        endif
c     Magnetic Axis Toroidal Field
        Btaxis=gprof(1)*alpha
c        Btaxis=0.D0
c        write(6,*) 'Magnetic Axis Toroidal Field=', Btaxis
c        if(iprocx.eq.1.and.iprocy.eq.1.and.iprocz.eq.1) then
c           write(ihis,*) 'Magnetic Axis Toroidal Field x R=', Btaxis
c           write(ihis,*) 'Magnetic Axis Toroidal Field=', Btaxis/rma
c        endif
c
c
        psifactor=(psibig(nx+1)-psibig(1))/(nx)
        psifactor=psifactor*0.03698/(psibig(nx+1)-psibig(1))
c
c     Reference quantities
        wmol=2.D0
        RUniversal=8314.5D0
c     rgas=RUniversal/wmol
c     This is a non-dim gas constant - Eos is p=rho T
        rgas=1.D0
c     Permeability of free space in Henry/m or kg-m/Coulomb^2
        mu0=4.D0*pi * 1.D-07
c     Reference Magnetic field in Tesla (kg/s/coulomb)
c        Bref=0.25D0
        Bref=1.D0
c     CDXU
c        Bref=0.375D0
c     DIII-D
c        Bref=2.0D0
c     ITER
c        Bref=5.0D0
c      Bref=BToroidalMax
c
c     Alternative pref is Bref^2/mu0 (in N/m^2)
      pref=Bref**2/mu0
c     Reference no. of molecules /m^3
      n0ref=1.5D19
c     Mass of each proton in Kg
      Mi=1.6605 D-27
c     Reference density rhoRef=n0ref*Mi*wmol
c      rhoRef=1.6605*wmol*1.D-08
      rhoRef=Mi*wmol*n0Ref
c     Reference temperature in K
      TRef=pref/(rhoRef*RUniversal/wmol)
c     Convert K to eV
      K2eV=1.1604D04
c     Reference temperature in eV
      TRefeV=Tref/K2eV
c
c     Reference velocity is the alfven velocity defined 
c        as Bref/sqrt(rhoRef*mu0)
      velRef=Bref/dsqrt(rhoRef*mu0)
c     
c     Reference length scale is 1m
      lref=1.D0
c
c     Reference time scale tauRef
      tauRef=lref/velRef
c
      pambient=0.03D0/mu0*alpha*alpha
c
        call iglobal(iprocx,1,ig)
        call jglobal(iprocy,1,jg)
c
           do j=iylo,iyhi,1
              do i=ixlo,ixhi,1
                 bk=0.5D0*(gprof(ig+i-1)+gprof(ig+i))*alpha
     &                /Rc(i,j)-Btaxis/Rc(i,j)
                 bi=psifactor*(dRdEtaC(i,j))
     &                /(Rc(i,j)*Jaci(i,j))*alpha
                 bj=psifactor*dZdEtaC(i,j)
     &                /Rc(i,j)
     &                /JacI(i,j)*alpha
c                 bk=0.D0
                 beta=-bi*dZdXiC(i,j)/Jaci(i,j)+bj*dRdXiC(i,j)/Jaci(i,j)
                 bxi=bi*dZdEtaC(i,j)-bj*dRdEtaC(i,j)
c                 write(6,*) 'Beta',i,j, beta, beta/(bk+Btaxis/Rc(i,j))
c$$$                 write(6,*) 'R x Beta', i,j,beta*Rc(i,j)*Jaci(i,j)
                 if(bxi.gt.1.D-13) then
                    write(6,*) 'Bxi=',i,j,bxi
                 endif
c$$$                 write(6,*) 'Bxi2=',i,j,bj*(-dZdXiC(i,j))
c$$$     &                -bi*(dRdXiC(i,j))
                 press=0.5D0*(pprof(ig+i-1)+pprof(ig+i))/psifactor
     &                *alpha*alpha
     &                +pambient/pRef
c
                 ux(i,j,1)=Bi
                 ux(i,j,2)=Bj
                 ux(i,j,3)=Bk
                 ux(i,j,4)=press+0.5D0*(bi*bi+bj*bj+bk*bk)
              enddo
           enddo
c
        return
        end