      subroutine bram
C...Translated by Pacific-Sierra Research VAST-90 2.06G2  10:55:53   6/14/01
C...Switches: -p4 -yb

c ************ ************ ************ ************ ************
c
c   routine to set initial values of many variables
c   BRAM  for 'M. Brambilla' who wrote the original version
c   for lower hybrid waves.
c
c ************ ************ ************ ************ ************
c
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
!      use Vparam
      use Viogrl
      use Viorlh
      use Vcraylh
!      use Viogrl, ONLY: calfa, px, sumpx, bk, refl, rtot, trans,
!     1   aksupw, ddk, nnk, joutrf, ireout, jumpw, jtest, ispec, ngra,
!     2   freqcy, width, dwall, heigth, dnedge, dndx, dxpl, phmin, phmax,
!     3   itrans, radius, nwguid, nmod, nphase, isimm, ndim0, ndim1
!      use Viorlh, ONLY: nprofn, psimhd, dpsimhd, rmaj, rrho, triang,
!     1   elong, nmhd, denc, tempec, atm, azi, conc, tempic, nspec,
!     2   mainsp, pavtne, pavtte, pavtti, pavpp, beta, betap, ali3,
!     3   ppsi, pjavg, pne, pte, nprofs, nprof, iedge, ipropt, pni, pti,
!     4   psinm, tani, nspect, omega, anzin, anzsup, anzinf, anplin,
!     5   anpsup, anpinf, pwcpl, pwtra, specfac, power, powers, thgril,
!     6   delthe, phase, rgrill, nnkpar, nnkpol, poldk, nthin, maxref,
!     7   igrill, islofa, irayio, ioread, nraypic, deykpr, psistep,
!     8   pkexpnt, nequip, iprint, indvar, iploxy, igraph, iplopr, jwfr,
!     9   pwe, pwe0, pwi, pwii0, pwev, pwiv, wkarea, pwal, ffast, tite,
!     .   alpene1, alpte1, frcte, frcne, frcti, gamte1, gamti1, alpti1,
!     1   epserr, epser1, epspfw, epspsw, epsray, rzero, rplasm,
!     2   bzero, btsign, tcurr, fbsc, rshift, dtrian, qsaf0, qsafa,
!     3   gamene1, presfr, iout, iout1, iout2, nharm, nequil, isafa,
!     4   polasp
      use kind_spec

      implicit none
C-----------------------------------------------
C   G l o b a l   P a r a m e t e r s
C-----------------------------------------------
C-----------------------------------------------
C   L o c a l   P a r a m e t e r s
C-----------------------------------------------
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      integer :: i, jout6, npart, nx, ny, j, ifrang
      real(kind=dp) :: anzacc, edgabs, bzedge, aksups, zpsi, zdpsi,
     1    zdivne,
     1   zdivte, zne, zti, zdecay, zdivni, zdivti
C-----------------------------------------------
C   E x t e r n a l   F u n c t i o n s
C-----------------------------------------------
      real(kind=dp),EXTERNAL::furmaj,furrho,futria,fuelon,prone,
     1   prote,proti
C-----------------------------------------------
c
c  parameters
c  Nspect-->nspectm, Nion-->nions,  Nua-->nuam,  Nta-->ntam, Lrz-->lrza,
c  Nthp11-->nthp11a,  Npnt-->npnta,   Mpnt-->mpnta
c
c------------ --------------
c   common blocks of Comiogrl
c---------------------------
c
c----------------------------
c       common block Comiorlh(1-8)
c
c  ipropt         integer   #option for n-T profile calculation: 1 if
c                        # calculated from prontp; else =0.


c--------SCC  5/23/95
c*******  TKM  06/97, 2/00
c*************************
c********    TKM    11/98
c*******************************
c  #   SCC 3/27/92   SCC 4/13/95
c---------------------------
c   common block for Comraylh
c-------------------------
c
c
c$$$$$$$$     TKM    11/98
c*********    TKM   2/97, 6/99, 9/99
c*********    TKM   2/97
c*******   TK Mau   8/98
c*************************
c  SCC  3/7/95

c
c ************ ************ ************ ************ ************
c
c      call dropfile(0)
c
c  title & date
c
      write (iout, 5000)
 5000 format('***',2x,"FW/SW raytracing code with current drive",/,5x,
     1   "Derived from RAYLH                 ",/,5x,
     2   "By Brambilla:22/11/85, BH:10/4/90, TKM/SCC/BH:11/21/91",/)
c
c ===============================================
c  plasma parameters
c ===============================================
c  toroidal magnetic field on axis (tesla)
      bzero = 2.17
      btsign = 1.
c  main and plasma radius (m)
      rzero = 1.68
      rplasm = 0.39
c  toroidal current (a)
      tcurr = 3.20E5
c  central electron density (cm-3)
c  from parameterization if <= 0
c      ignored if ifrang not equal to 0
      denc = 2.E13
c ------------------------------ (main only)
c  specify central (ifrang=0) or line averaged (ifrang=1) density
c      ifrang = -1 for exponentially parametrized profiles
c      if ifrang>0 specify the density per frange zdnfrg (cm-3)
C        ifrang = 0
C        zdnfrg = 1.4e18
c  simulation of the scrape-off layer (iedge=1)
c        iedge = 0
c ------------------------------ (end)
c  central electron temperature (kev)
c  from parameterization if <= 0
c      ignored if ifrang not equal to 0
      tempec = 1.200
c  n.of ion species  (up to nions species allowed)
      nspec = 1
c  reference species
      mainsp = 1
c  mass and charges (atomic units)
      atm(1) = 2.
      azi(1) = 1.
c        atm(2) = 1.
c        azi(2) = 1.
c  concentrations conc(j,nsp)=n(j,nsp)/conc(j,nsp)=n(e); charge neutrality will be

c          implemented by re-evaluating n(j,nsp)
c  central ion temperatures (kev)
c  all ion species will have same temperature;
c  if they are different, one should input from curray_go or eqdskex
      tempic(1) = 0.9
c        tempic(2) = 20.
c  nprof = n. of points in the radial profile
c      accepted if 11 <= nprof <= nprofn
c**bh    nprof = 41
c  specify here if a standard case has to be run:
c        nequil=1 - data entered by the user
c                   (equilibrium code or diagnostic)
c        nequil=0 - circular or elliptic cross-section
c                   (specify the plasma radius (m.), the
c                   poloidal aspect ratio and the relative
c                   shift of the magnetic axis)
      nequil = 0
c  poloidal cross-section dimension and shape
c        rshift - shift of the magnetic axis (m)
c        polasp - poloidal aspect ratio (b/a)
c                 (1): limit on the magnetic axis
c                 (2): of the external magn. surface
c        dtrian - triangularity of the external magn. surf.
c               - ignored if nequil = 1
      rshift = 0.0
      polasp(1) = 1.
      polasp(2) = 1.
      dtrian = 0.
c****************9/15/91   SCC
c    epserr=tolerence of error in dispersion relation in raytracing
c********************
      epserr = 0.01
c
c ===============================================
c  lh heating parameters
c ===============================================
c   frequency (hz)
      freqcy = 2.45E9
      omega = 2.*3.14159265359*freqcy
c  total power coupled (mw)
      power = 1.
      powers(1) = 1.
      powers(2) = 0.
c  thgril - central position of the antenna (deg.)
      thgril(:nspect) = 0.
c  igrill = 1 - grill code called
c  igrill = 0 - grill code not called
c  igrill = -1 - fit to P(n_parallel)
      igrill = 0
c  wave emitted by the antenna
c   islofa = 1 : slow wave
c   islofa = -1 : fast wave (igrill = 0)
      islofa = 1
c
c ===============================================
c  grill parameters  (ignored if igrill = 0)
c ===============================================
c number of waveguides
      nwguid = 8
c waveguide distance, cm.
      width = 2.1
c waveguide separation, cm.
      dwall = 0.80
c heigth of the guides
      heigth = 16.5
c total number of modes
      nmod = 3
c density gradient at the plasma edge, cm-4
      dndx = 1.E11
c density at the plasma edge, cm-3
      dnedge = 1.E11
c effective absorption for nonaccessible waves
c   between npar=1 and npar=anzacc
      anzacc = 0.
      edgabs = 0.1
c magnetic field near the antenna
      bzedge = 1.8
c distance of the plasma from the wall, cm.
      dxpl = 0.00
c------------------------------------------------------------------
c phases for output
c        nphase > 0 - uniform excitation
c        nphase < 0 - excitations specified through subroutine
c                     config(calfa,ncall)
c                     in this case specify wether the config.
c                     is simmetric (isimm=1) or not (isimm=2)
      nphase = 1
      phmin = 180.
      phmax = 180.
c------------------------------------------------------------------
c output control
c       = 1  -  only reflection coefficients
c       = 2  -  also complex coefficients
      jout6 = 1
c  repetition of the waveguides parameters for each case
c     ireout = n. of cases - yes
c     ireout = 0 - no
      ireout = 0
c  max n-parallel for output
c remember: aksups should not be smaller than the largest of anzsup
      aksupw = 7.
      aksups = 7.
c  n. of points and mesh steps
      nnk(1) = 300
      ddk(1) = (aksups - 1.)/float(nnk(1))
      nnk(2) = 200
      ddk(2) = 0.1
c  jumpw = jump for print
      jumpw = 5
c  npart = n. of step in the unaccessible region
      npart = 5
c
c ===============================================
c  discretisation of the power spectre
c      for ray tracing purposes
c ===============================================
c  nnkpar(i) - n of values of nparallel in interval i.
c  anzsup - largest absolute value of k parallel
c  anzinf - smallest absolute value of k parallel
c    if nnkpar=1  anzinf(1) will be taken
c ****** actual
      nnkpar(1) = 40
      anzinf(1) = 1.05
      anzsup(1) = 3.00
c ******
      nnkpar(2) = 0
      anzinf(2) = -2.
      anzsup(2) = -4.
c ******
      nnkpar(3:nspect) = 0
      anzinf(3:nspect) = 0.0
      anzsup(3:nspect) = 0.0
c  nthin - n. of rays in the poloidal cross-section
c    (should not exceed ngr)
      nthin = 1
c  nequip - max. n. of equiphase surfaces with output
c    ignored if igraph not zero
c  deykpr - step of eikon for printing
      nequip = 200
      deykpr = 2.
c ************ ************ ************ ************
c independent variable
c  indvar = 0 - poloidal phase
c  indvar = 1 - toroidal phase
      indvar = 1
c ************ ************ ************ ************
c  maxref = n. of boundary reflections allowed
      maxref = 1
c ************ ************ ************ ************
c details of print
c  iprint = 0 - output suppressed
c  iprint = 1 - normal output
c         = 2 - more detailed output
c         = 3 - det. output at each call of output1(x)
c               for tests (no graphical output)
      iprint = 0
c  iplopr = 1  -  plots of the profiles
      iplopr = 0
c ************ ************ ************ ************
c  igraph = 0  no graphical output
      igraph = 0
      if (igraph /= 0) then
c ************ ************ ************ ************
c  wavefronts to be plotted every jwfr point
c      ignored if nthin=1
         jwfr = 1
c ************ ************ ************ ************
c  list of plots to be made
         iploxy(:,:nions) = 0
c ===================================================================
c  npolxy(nxvarb,nyvarb) different from zero produces a plot
c                        according to the following table:
c ====================================================================
c  nxvarb=1  indep. var. seikon  (phase in the poloidal plane)
c         2              psi (magnetic surfaces label)
c ====================================================================
c  nyvarb=1  dep. var.   skp
c         2              skpar
c         3              skr
c         4              skth
c         5              skx
c         6              skz
c         7              spower
c         8              sphi
c ===================================================================
         iploxy(2,1) = 1
         iploxy(2,2) = 1
c        iploxy(2,8) = 1
      endif
c   parameters for spectre
c        iabs(ngra) = n of graphs of spectre per page
c        ngra.gt.0    integral of px also plotted
c        ngra.lt.0    only px plotted
c  ignored if igrill = 0
      ngra = 1
c ===============================================
c  mhd configuration parameters
c ===============================================
c  nmhd = n. of points in the radial profiles
c      of rmaj(psi), rrho(psi), triang(psi), ellipt(psi)
c      accepted if 11 <= nmhd <= nprofn
c
      nmhd = 11
c
      dpsimhd = 1./(nmhd - 1.)
c
      do i = 1, nmhd
         psimhd(i) = (i - 1.)*dpsimhd
         zpsi = psimhd(i)
         rmaj(i) = furmaj(zpsi)
         rrho(i) = furrho(zpsi)
         triang(i) = futria(zpsi)
         elong(i) = fuelon(zpsi)
      end do
c
c ************ ************ ************ ************ ************
c
c  standard case with analytic profiles for density & temperatures
c
c ************ ************ ************ ************ ************
c
      conc(:nprofs,:nspec) = 0.
      pni(:nprofs,:nspec) = 0.
      pti(:nprofs,:nspec) = 0.
      conc(:nprofs,mainsp) = 1.
c
      zdpsi = 1./float(nprofs - 1)
c     print *,'nprofs = ',nprofs,' zdpsi = ',zdpsi,' ifrang = ',ifrang
c ======================================================================
c ====== only with the exponential version of the profiles =============
      zdivne = 1.
      zdivte = 1.
      if (ifrang == (-1)) then
         zdivne = 1.E14
         zdivte = 1.E-3
         zne = prone(0.)
         denc = 1.E14*zne
         tempec = 1.E-3*prote(0.)
         zti = min(1.0_dp,0.60_dp*zne + 0.40_dp)
         tempic(1) = zti*tempec
      endif
c ======================================================================
      do j = 1, nprofs
c        print *,'j = ',j
         psi = float(j - 1)*zdpsi
c        print *,'psi = ',psi
         ppsi(j) = psi
         pte(j) = prote(psi)*zdivte
         pne(j) = prone(psi)*zdivne
         pjavg(j) = pte(j)**1.5
         do i = 1, nspec
            pni(j,i) = pne(j)*conc(j,i)
            pti(j,i) = proti(psi)
         end do
      end do
      tempic(:nspec) = pti(1,:nspec)
c
c ************ ************ ************ ************
c  profile averaged density
c
C     if(ifrang.eq.1)  then
C        zsum = pne(1)
C     do 510  n=2,nprofs,2
C 510    zsum = zsum + 4.*pne(n) + pne(n+1)
C        zfact = zdpsi*zsum/(3.*pne(1))
C        zdivne = 0.5*zdnfrg*float(ifrang)/zfact
C     do 520  n=1,nprofs
C 520    pne(n) = zdivne*pne(n)
C        denc = pne(1)
C     end if
c
c ************ ************ ************ ************
c  edge simulation
c
      if (iedge /= 0) then
         do j = nprofs + 1, nprof
            ppsi(j) = float(j - 1)*zdpsi
         end do
         zdivne = pne(nprofs)
         zdivte = pte(nprofs)
         do j = nprofs + 1, nprof
            zdecay = exp((-float(j - nprofs)/2.))
            pne(j) = zdivne*zdecay
            pte(j) = zdivte*zdecay
            pjavg(j) = 0.0
            pti(j,:nspec) = pti(nprofs,:nspec)*zdecay
            pni(j,:nspec) = pni(nprofs,:nspec)*zdecay
         end do
c  extend MHD
         nmhd = 11
         dpsimhd = ppsi(nprof)/(nmhd - 1.)
         do i = 1, nmhd
            psimhd(i) = (i - 1.)*dpsimhd
            zpsi = psimhd(i)
            rmaj(i) = furmaj(zpsi)
            rrho(i) = furrho(zpsi)
            triang(i) = futria(zpsi)
            elong(i) = fuelon(zpsi)
         end do
      endif
c
c ************ ************ ************ ************

      return
      end subroutine bram
