c  Mahc_Evolve: Numerical Evolution of the General Relativistic Hydro Equations
c  Copyright (C) 2001 Mark Miller 

/*@@
  @file      MahcParam.F
  @date      December 1999
  @author     Mark Miller
  @desc
     Sets all the control parameters in MAHC.
  @enddesc
@@*/

#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Arguments.h"
#include "cctk_DefineThorn.h"

/*@@
  @routine    MahcParam
  @date       December 1999
  @author     Mark Miller
  @desc
     Sets all the control parameters in MAHC.
  @enddesc
  @calls
  @calledby
  @history
  @endhistory
@@*/

      subroutine MahcParam(CCTK_FARGUMENTS)
      implicit none
      DECLARE_CCTK_FARGUMENTS
      DECLARE_CCTK_PARAMETERS
      integer CCTK_Equals
#ifdef CACTUSEOS_EOS_BASE
#include "EOS_Base.inc"
#endif


      mahc_eos_k = eos_k
      mahc_eos_gamma = eos_gamma
      mahc_densmin = densmin
      mahc_perc_ptol = perc_ptol
      mahc_del_ptol = del_ptol

      if(CCTK_EQUALS(hydro_order,'fluxlim')) then
         if ( (mahc_stencil_size .lt. 2) .or.
     .        (cctk_nghostzones(1) .lt. 2) .or.
     .        (cctk_nghostzones(2) .lt. 2) .or.
     .        (cctk_nghostzones(3) .lt. 2) ) then
            write(*,*) 'MahcParam: if you want to use fluxlim, you need'
            write(*,*) '   to make stencil and ghost bigger than 2!'
            call CCTK_WARN(0,"MahcParam")
         endif
      endif

      if(CCTK_EQUALS(hydro_order,'upwind')) then
         if ( (mahc_stencil_size .gt. 1) .or.
     .        (cctk_nghostzones(1) .gt. 1) .or.
     .        (cctk_nghostzones(2) .gt. 1) .or.
     .        (cctk_nghostzones(3) .gt. 1) ) then
            write(*,*) 'MahcParam: if you want to use upwind, you need'
            write(*,*) '   to make stencil and ghost 1!'
            call CCTK_WARN(0,"MahcParam")
         endif
      endif


c     SET MAHC_EPS_MIN
c     this number is set as a multiple of the eps of the atmosphere, which
c     depends on the following parameters in the .par file:
c     eos_gamma,eos_k,centden,ns_atmos_fact

      if (mahc_eos_gamma.ne.1.0d0) then
         mahc_eps_min = epsmin_fact * mahc_eos_k *
     .        ((centden/ns_atmos_fact) ** (mahc_eos_gamma - 1.0)) /
     .        ((mahc_eos_gamma - 1.0) + 1.0d-29)
      else
         write(*,*)'using special form for mahc_eps_min'
         mahc_eps_min = 1.0d-29
      endif

c     set the mahc_primvar_solver variable
      if(CCTK_EQUALS(mahc_primsolve_type,"Pressure")) then
         mahc_primvar_solver = 1
      elseif(CCTK_EQUALS(mahc_primsolve_type,"Sigma")) then
         mahc_primvar_solver = 2
      endif

c     clear out mahc_prim_energyadd, which keeps track of energy 
c     additions in primative variable solver (sigma type)
      mahc_prim_energyadd = 0.0d0

c     interface with EOS_Base
c$$$#ifdef CACTUSEOS_EOS_BASE
c$$$      if(CCTK_EQUALS(mahc_eos,"tables")) then
c$$$         call EOS_Handle(mahc_eos_handle,mahc_eos_table_type)
c$$$         if(mahc_eos_handle .lt. 0) then
c$$$            write(*,*) 'Error:MahcParam: could not find EOS handle:',
c$$$     .                 mahc_eos_table_type
c$$$            call CCTK_WARN(0,"MahcParam")
c$$$         endif
c$$$      endif
c$$$#else
c$$$      mahc_eos_handle = 0
c$$$#endif

      return
      end