       subroutine gridx1
c
c    dac:zeus3d.gridx1 <----------------------- initialises the x1a-grid
c    from jms:zeus3d.gridi                                     may, 1992
c
c    written by: David Clarke
c    modified 1:
c
c  PURPOSE:  Initialises the x1a-grid as specified by the parameters in
c  the namelist "ggen1".
c
c  LOCAL VARIABLES:
c
c  EXTERNALS:
c    ROUND
c    MODULO
c
c-----------------------------------------------------------------------
c
c
c-----------------------------------------------------------------------
c------------------------ IMPLICIT STATEMENTS --------------------------
c-----------------------------------------------------------------------
c
       implicit      none
c
c-----------------------------------------------------------------------
c----------------------------- PARAMETERS ------------------------------
c-----------------------------------------------------------------------
c
c      Primary parameters:
c
       integer       in      , jn      , kn      , ijkn    , nxpx
     1             , nypx    , nxrd    , nyrd    , niov    , nios
     2             , ntsl    , mn      , nmd     , isig
       real*8        pi  , tiny  , huge  , one  , mone  , haf  , zro
c
       parameter     ( in = 5005, jn = 1, kn = 10, ijkn = 5005)
       parameter     ( nxpx = 1, nypx = 1, nxrd =   1, nyrd =   1 )
       parameter     ( niov =  20, nios =  20, ntsl =1000 )
       parameter     ( mn   =   1, nmd  =   1 )
       parameter     ( pi   = 3.1415926535898 )
       parameter     ( isig = 9, tiny = 1.0d-99, huge = 1.0d+99 )
       parameter     ( one  = 1.0d+00, mone = -1.0d+00
     1                               , haf  =  0.5d+00, zro = 0.0d+00)
c
c      Secondary parameters:
c
       integer        idim    , jdim
     1              , nt01    , nt02    , nt03    , nt04    , nt05
     2              , nt06    , nt07    , nt08    , nt09    , nt10
     3              , nt11    , nt12    , nt13    , nt14    , nt15
     4              , nt16    , nt17    , nt18    , nt19    , nt20
       parameter     ( idim =   kn, jdim =   in )
       parameter     ( nt01 = 600, nt02 = 600, nt03 =   7, nt04 =   5 )
       parameter     ( nt05 =  18, nt06 =   1, nt07 =  10, nt08 =  12 )
       parameter     ( nt09 =   5, nt10 =  54, nt11 =   8 )
       parameter     ( nt12 =        16, nt17 =         4 )
       parameter     ( nt13 = nt12 +  4, nt18 = nt17 +  1 )
       parameter     ( nt14 = nt13     , nt19 = nt18      )
       parameter     ( nt15 = nt14 + 30, nt20 = nt19 +  3 )
       parameter     ( nt16 = nt14 + 42 )
c
c-----------------------------------------------------------------------
c--------------------------- ROOT VARIABLES ----------------------------
c-----------------------------------------------------------------------
c
c============================ DECLARATIONS =============================
c
c------------------------- boundary variables --------------------------
c
c      Inner and outer i-boundaries
c
       integer       niib    (  jn,  kn), noib    (  jn,  kn)
     1             , niib2   (  jn,  kn), noib2   (  jn,  kn)
     2             , niib3   (  jn,  kn), noib3   (  jn,  kn)
     3             , niib23  (  jn,  kn), noib23  (  jn,  kn)
c
       real*8        diib1   (  jn,  kn), doib1   (  jn,  kn)
     1             , diib2   (  jn,  kn), doib2   (  jn,  kn)
     2             , v1iib1  (  jn,  kn), v1oib1  (  jn,  kn)
     3             , v1iib2  (  jn,  kn), v1oib2  (  jn,  kn)
     4             , v2iib1  (  jn,  kn), v2oib1  (  jn,  kn)
     5             , v2iib2  (  jn,  kn), v2oib2  (  jn,  kn)
     6             , v3iib1  (  jn,  kn), v3oib1  (  jn,  kn)
     7             , v3iib2  (  jn,  kn), v3oib2  (  jn,  kn)
       real*8        eiib1   (  jn,  kn), eoib1   (  jn,  kn)
     1             , eiib2   (  jn,  kn), eoib2   (  jn,  kn)
       real*8        b1iib1  (  jn,  kn), b1oib1  (  jn,  kn)
     1             , b1iib2  (  jn,  kn), b1oib2  (  jn,  kn)
     2             , b2iib1  (  jn,  kn), b2oib1  (  jn,  kn)
     3             , b2iib2  (  jn,  kn), b2oib2  (  jn,  kn)
     4             , b3iib1  (  jn,  kn), b3oib1  (  jn,  kn)
     5             , b3iib2  (  jn,  kn), b3oib2  (  jn,  kn)
     6             , emf1iib1(  jn,  kn), emf1oib1(  jn,  kn)
     7             , emf1iib2(  jn,  kn), emf1oib2(  jn,  kn)
     7             , emf1iib3(  jn,  kn), emf1oib3(  jn,  kn)
     8             , emf2iib1(  jn,  kn), emf2oib1(  jn,  kn)
       real*8        emf2iib2(  jn,  kn), emf2oib2(  jn,  kn)
     1             , emf2iib3(  jn,  kn), emf2oib3(  jn,  kn)
     1             , emf3iib1(  jn,  kn), emf3oib1(  jn,  kn)
     2             , emf3iib2(  jn,  kn), emf3oib2(  jn,  kn)
     3             , emf3iib3(  jn,  kn), emf3oib3(  jn,  kn)
     4             , a1iib1  (  jn,  kn)
     5             , a1iib2  (  jn,  kn)
     6             , a1iib3  (  jn,  kn)
     6             , a1iib4  (  jn,  kn)
     7             , a2iib1  (  jn,  kn)
       real*8        a2iib2  (  jn,  kn)
     1             , a2iib3  (  jn,  kn)
     1             , a2iib4  (  jn,  kn)
     1             , a3iib1  (  jn,  kn)
     2             , a3iib2  (  jn,  kn)
     3             , a3iib3  (  jn,  kn)
     4             , a3iib4  (  jn,  kn)
c
c      Inner and outer k-boundaries
c
       integer       nikb    (  in,  jn), nokb    (  in,  jn)
     1             , nikb1   (  in,  jn), nokb1   (  in,  jn)
     2             , nikb2   (  in,  jn), nokb2   (  in,  jn)
     3             , nikb12  (  in,  jn), nokb12  (  in,  jn)
c
       real*8        dikb1   (  in,  jn), dokb1   (  in,  jn)
     1             , dikb2   (  in,  jn), dokb2   (  in,  jn)
     2             , v1ikb1  (  in,  jn), v1okb1  (  in,  jn)
     3             , v1ikb2  (  in,  jn), v1okb2  (  in,  jn)
     4             , v2ikb1  (  in,  jn), v2okb1  (  in,  jn)
     5             , v2ikb2  (  in,  jn), v2okb2  (  in,  jn)
     6             , v3ikb1  (  in,  jn), v3okb1  (  in,  jn)
     7             , v3ikb2  (  in,  jn), v3okb2  (  in,  jn)
       real*8        eikb1   (  in,  jn), eokb1   (  in,  jn)
     1             , eikb2   (  in,  jn), eokb2   (  in,  jn)
       real*8        b1ikb1  (  in,  jn), b1okb1  (  in,  jn)
     1             , b1ikb2  (  in,  jn), b1okb2  (  in,  jn)
     2             , b2ikb1  (  in,  jn), b2okb1  (  in,  jn)
     3             , b2ikb2  (  in,  jn), b2okb2  (  in,  jn)
     4             , b3ikb1  (  in,  jn), b3okb1  (  in,  jn)
     5             , b3ikb2  (  in,  jn), b3okb2  (  in,  jn)
     6             , emf1ikb1(  in,  jn), emf1okb1(  in,  jn)
     7             , emf1ikb2(  in,  jn), emf1okb2(  in,  jn)
     8             , emf1ikb3(  in,  jn), emf1okb3(  in,  jn)
     9             , emf2ikb1(  in,  jn), emf2okb1(  in,  jn)
       real*8        emf2ikb2(  in,  jn), emf2okb2(  in,  jn)
     1             , emf2ikb3(  in,  jn), emf2okb3(  in,  jn)
     2             , emf3ikb1(  in,  jn), emf3okb1(  in,  jn)
     3             , emf3ikb2(  in,  jn), emf3okb2(  in,  jn)
     3             , emf3ikb3(  in,  jn), emf3okb3(  in,  jn)
c
c------------------------- display variables ---------------------------
c
       integer       iswdis  , idisdmp , ndis    , mdis
       real*8        tdis    , dtdis
c
       character*2   disvar  (niov)
       integer       idismn  (nios), jdismn  (nios), kdismn  (nios)
     1             , idismx  (nios), jdismx  (nios), kdismx  (nios)
     2             , idisdir (nios), ldis    (nios)
c
c-------------------------- field variables ----------------------------
c
       real*8        d      (  in,  jn,  kn), v1     (  in,  jn,  kn)
     1             , v2     (  in,  jn,  kn), v3     (  in,  jn,  kn)
     2             , e      (  in,  jn,  kn)
     4             , b1     (  in,  jn,  kn), b2     (  in,  jn,  kn)
     5             , b3     (  in,  jn,  kn)
c
c---------------------- essential grid variables -----------------------
c
       integer       ismn    , is      , iemx    , ie      , ione
     1             , jsmn    , js      , jemx    , je      , jone
     2             , ksmn    , ks      , kemx    , ke      , kone
     3             , ia      , ja      , ka      , igcon
       real*8        x1scale , x2scale , x3scale , x1fac   , x2fac
     1             , x3fac
c
       integer       istretch(   4), jstretch(   4), kstretch(   4)
       real*8        x1a     (  in), x2a     (  jn), x3a     (  kn)
c
c---------------------------- HDF variables ----------------------------
c
       integer       iswhdf  , ihdfdmp , nhdf
       real*8        thdf    , dthdf
c
       character*2   hdfvar  (niov)
c
c--------------------------- pixel variables ---------------------------
c
       integer       iswpix  , ipixdmp , npix    , mpix    , ipixmm
     1             , ncpix   , rorhpix , iref    , jref
       real*8        tpix    , dtpix
c
       character*2   pixvar  (niov)
       integer       nlpix   (niov)
       integer       ipixdir (nios), lpix    (nios), npi     (nios)
     1             , npj     (nios)
       real*8        pixmin  (niov), pixmax  (niov)
       real*8        x1pxmn  (nios), x2pxmn  (nios), x3pxmn  (nios)
     1             , x1pxmx  (nios), x2pxmx  (nios), x3pxmx  (nios)
c
c------------------------- 1-D NCAR variables --------------------------
c
       integer       iswplt1 , iplt1dmp, nplt1   , mplt1   , iplt1mm
     1             , corl
       real*8        tplt1   , dtplt1
c
       character*2   plt1var (niov)
       integer       nlplt1  (niov)
       integer       iplt1   (nios), jplt1   (nios), kplt1   (nios)
     1             , iplt1dir(nios)
       real*8        plt1min (niov), plt1max (niov)
c
c------------------------- 2-D NCAR variables --------------------------
c
       integer       iswplt2 , iplt2dmp, nplt2   , mplt2   , iplt2mm
     1             , numcl   , incvx   , incvy
       real*8        tplt2   , dtplt2  , vscale
c
       character*2   plt2var (niov)
       integer       nlplt2  (niov)
       integer       iplt2dir(nios), lplt2   (nios)
       real*8        plt2min (niov), plt2max (niov)
       real*8        x1p2mn  (nios), x1p2mx  (nios), x2p2mn  (nios)
     1             , x2p2mx  (nios), x3p2mn  (nios), x3p2mx  (nios)
c
c--------------------------- RADIO variables ---------------------------
c
       integer       iswrad  , iraddmp , nrad    , iradmm  , ncrad
     1             , rorhrad , incpx   , incpy   , itype   , ntheta
     2             , nphi    , icnvlv  , eorb    , porf    , bworb
     3             , nisc    , njsc
       real*8        trad    , dtrad   , x1rdmn  , x2rdmn  , x3rdmn
     1             , x1rdmx  , x2rdmx  , x3rdmx  , theta   , thetamin
     2             , thetamax, dtheta  , phi     , phimin  , phimax
     3             , dphi    , alpha   , freq    , dsw     , btsw
     4             , dlower  , elower  , selower , pblower , shlower
     5             , vvlower , dupper  , eupper  , seupper , pbupper
     6             , shupper , vvupper , bmajor  , bminor  , bpa
     7             , cpb     , icut    , pcut    , vlmin   , pscale
c
       character*2   radvar  (niov)
       integer       nlrad   (niov)
       real*8        radmin  (niov), radmax  (niov)
c
c-------------------------- sundry variables ---------------------------
c
       integer       iswint  , iswdmp  , iresdmp , ilogdmp , incdmp
     1             , iord    , iordb1  , iordb2  , iordb3  , iordd
     2             , iorde   , iords1  , iords2  , iords3  , ippa
     3             , istp    , istpb1  , istpb2  , istpb3  , istpd
     4             , istpe   , istps1  , istps2  , istps3  , ifsen
     5             , nhy     , nlim    , nwarn   , ix1x2x3 , icool
     6             , itote   , imxd    , nmat    , iwiggle , ibpol
     7             , ijetmn  , ijetmx  , jjetmn  , jjetmx  , kjetmn
     8             , kjetmx  , igpiib
     9             , igpoib  , igpijb  , igpojb  , igpikb  , igpokb
       real*8        tdmp    , b1floor , b2floor , b3floor , dfloor
     1             , efloor  , v1floor , v2floor , v3floor , time
     2             , tlim    , trem    , tsave   , ttotal  , tused
     3             , dt      , dtrat   , dtmin   , dtnew   , dtcs
     4             , dtal    , dtqq    , dtqqi2  , dtv1    , dtv2
     5             , dtv3    , courno  , qcon    , qlin    , omega
     6             , thetai  , thetaj  , thetak  , rmode   , mamb
     7             , zsigma  , wnz     , ojet    , sig1    , sig2
     8             , sig3    , xibgen  , rjet    , g       , ptmass
     9             , x1ptm   , x2ptm   , x3ptm
c
       integer       niso    (  mn)
       real*8        ciso    (  mn), rmetal  (  mn), gam     (  mn)
     1             , gamma   (  mn), rgas    (  mn)
c
c------------------------ time slice variables -------------------------
c
       integer       iswtsl  , iswtslp , itsldmp , itslpdmp
       real*8        ttsl    , ttslp   , dttsl   , dttslp  , tslpmn
     1             , tslpmx
c
       real*8        timetsl (ntsl), massp   (ntsl), x1cmp   (ntsl)
     1             , x2cmp   (ntsl), x3cmp   (ntsl), eintp   (ntsl)
     2             , ek1p    (ntsl), ek2p    (ntsl), ek3p    (ntsl)
     3             , ekp     (ntsl), epotp   (ntsl), tep     (ntsl)
     4             , l1p     (ntsl), l2p     (ntsl), l3p     (ntsl)
     5             , dminp   (ntsl), dmaxp   (ntsl), seminp  (ntsl)
     6             , semaxp  (ntsl), pminp   (ntsl), pmaxp   (ntsl)
     7             , v1minp  (ntsl), v1maxp  (ntsl), v2minp  (ntsl)
     8             , v2maxp  (ntsl), v3minp  (ntsl), v3maxp  (ntsl)
     9             , vminp   (ntsl), vmaxp   (ntsl), dtcsp   (ntsl)
       real*8        dtv1p   (ntsl), dtv2p   (ntsl), dtv3p   (ntsl)
     1             , dtqqp   (ntsl), dtnewp  (ntsl), emag1p  (ntsl)
     2             , emag2p  (ntsl), emag3p  (ntsl), emagp   (ntsl)
     3             , b1flxp  (ntsl), b2flxp  (ntsl), b3flxp  (ntsl)
     4             , dvbp    (ntsl), b1minp  (ntsl), b1maxp  (ntsl)
     5             , b2minp  (ntsl), b2maxp  (ntsl), b3minp  (ntsl)
     6             , b3maxp  (ntsl), bminp   (ntsl), bmaxp   (ntsl)
     7             , dvbminp (ntsl), dvbmaxp (ntsl), dtalp   (ntsl)
c
c------------------------- user dump variables -------------------------
c
       integer       iswusr  , iusrdmp
       real*8        tusr    , dtusr
c
c--------------------------- voxel variables ---------------------------
c
       integer       iswvox  , ivoxdmp , nvox    , ivoxmm  , ncvox
     1             , rorhvox , nvi     , nvj     , nvk
       real*8        tvox    , dtvox   , x1vxmn  , x2vxmn  , x3vxmn
     1             , x1vxmx  , x2vxmx  , x3vxmx
c
       character*2   voxvar  (niov)
       integer       nlvox   (niov)
       real*8        voxmin  (niov), voxmax  (niov)
c
c======================== EQUIVALENCE STATEMENTS =======================
c
c------------------------------- scalars -------------------------------
c
       integer       insc    (nt01)
       real*8        rlsc    (nt02)
c
       equivalence
     1   (insc(  1),iswdis  ),(insc(  2),idisdmp ),(insc(  3),ndis    )
     2  ,(insc(  4),mdis    )
       equivalence
     1   (insc( 31),ismn    ),(insc( 32),is      ),(insc( 33),iemx    )
     2  ,(insc( 34),ie      ),(insc( 35),ione    ),(insc( 36),jsmn    )
     3  ,(insc( 37),js      ),(insc( 38),jemx    ),(insc( 39),je      )
     4  ,(insc( 40),jone    ),(insc( 41),ksmn    ),(insc( 42),ks      )
     5  ,(insc( 43),kemx    ),(insc( 44),ke      ),(insc( 45),kone    )
     6  ,(insc( 46),ia      ),(insc( 47),ja      ),(insc( 48),ka      )
     7  ,(insc( 49),igcon   )
       equivalence
     1   (insc( 91),istretch),(insc( 95),jstretch),(insc( 99),kstretch)
       equivalence
     1   (insc(121),iswhdf  ),(insc(122),ihdfdmp ),(insc(123),nhdf    )
       equivalence
     1   (insc(151),iswpix  ),(insc(152),ipixdmp ),(insc(153),npix    )
     2  ,(insc(154),mpix    ),(insc(155),ipixmm  ),(insc(156),ncpix   )
     3  ,(insc(157),rorhpix ),(insc(158),iref    ),(insc(159),jref    )
       equivalence
     1   (insc(181),iswplt1 ),(insc(182),iplt1dmp),(insc(183),nplt1   )
     2  ,(insc(184),mplt1   ),(insc(185),iplt1mm ),(insc(186),corl    )
       equivalence
     1   (insc(211),iswplt2 ),(insc(212),iplt2dmp),(insc(213),nplt2   )
     2  ,(insc(214),mplt2   ),(insc(215),iplt2mm ),(insc(216),numcl   )
     3  ,(insc(217),incvx   ),(insc(218),incvy   )
       equivalence
     1   (insc(241),iswrad  ),(insc(242),iraddmp ),(insc(243),nrad    )
     2  ,(insc(244),iradmm  ),(insc(245),ncrad   ),(insc(246),rorhrad )
     3  ,(insc(247),incpx   ),(insc(248),incpy   ),(insc(249),itype   )
     4  ,(insc(250),ntheta  ),(insc(251),nphi    ),(insc(252),icnvlv  )
     5  ,(insc(253),eorb    ),(insc(254),porf    ),(insc(255),bworb   )
     6  ,(insc(256),nisc    ),(insc(257),njsc    )
       equivalence
     1   (insc(331),iswint  ),(insc(332),iswdmp  ),(insc(333),iresdmp )
     2  ,(insc(334),ilogdmp ),(insc(335),incdmp  ),(insc(336),iord    )
     3  ,(insc(337),iordb1  ),(insc(338),iordb2  ),(insc(339),iordb3  )
     4  ,(insc(340),iordd   ),(insc(341),iorde   ),(insc(342),iords1  )
     5  ,(insc(343),iords2  ),(insc(344),iords3  ),(insc(345),ippa    )
     6  ,(insc(346),istp    ),(insc(347),istpb1  ),(insc(348),istpb2  )
     7  ,(insc(349),istpb3  ),(insc(350),istpd   ),(insc(351),istpe   )
     8  ,(insc(352),istps1  ),(insc(353),istps2  ),(insc(354),istps3  )
     9  ,(insc(355),ifsen   ),(insc(356),nhy     ),(insc(357),nlim    )
       equivalence
     1   (insc(358),nwarn   ),(insc(359),ix1x2x3 ),(insc(360),icool   )
     2  ,(insc(361),itote   ),(insc(362),imxd    ),(insc(363),nmat    )
     3  ,(insc(364),iwiggle ),(insc(365),ibpol   ),(insc(366),ijetmn  )
     4  ,(insc(367),ijetmx  ),(insc(368),jjetmn  ),(insc(369),jjetmx  )
     5  ,(insc(370),kjetmn  ),(insc(371),kjetmx  )
       equivalence
     1   (insc(375),igpiib  )
     2  ,(insc(376),igpoib  ),(insc(377),igpijb  ),(insc(378),igpojb  )
     3  ,(insc(379),igpikb  ),(insc(380),igpokb  )
       equivalence
     1   (insc(451),iswtsl  ),(insc(452),iswtslp ),(insc(453),itsldmp )
     2  ,(insc(454),itslpdmp)
       equivalence
     1   (insc(481),iswusr  ),(insc(482),iusrdmp )
       equivalence
     1   (insc(511),iswvox  ),(insc(512),ivoxdmp ),(insc(513),nvox    )
     2  ,(insc(514),ivoxmm  ),(insc(515),ncvox   ),(insc(516),rorhvox )
     3  ,(insc(517),nvi     ),(insc(518),nvj     ),(insc(519),nvk     )
c
       equivalence
     1   (rlsc(  1),tdis    ),(rlsc(  2),dtdis   )
       equivalence
     1   (rlsc( 31),x1scale ),(rlsc( 32),x2scale ),(rlsc( 33),x3scale )
     2  ,(rlsc( 34),x1fac   ),(rlsc( 35),x2fac   ),(rlsc( 36),x3fac   )
       equivalence
     1   (rlsc(121),thdf    ),(rlsc(122),dthdf   )
       equivalence
     1   (rlsc(151),tpix    ),(rlsc(152),dtpix   )
       equivalence
     1   (rlsc(181),tplt1   ),(rlsc(182),dtplt1  )
       equivalence
     1   (rlsc(211),tplt2   ),(rlsc(212),dtplt2  ),(rlsc(213),vscale  )
       equivalence
     1   (rlsc(241),trad    ),(rlsc(242),dtrad   ),(rlsc(243),x1rdmn  )
     2  ,(rlsc(244),x2rdmn  ),(rlsc(245),x3rdmn  ),(rlsc(246),x1rdmx  )
     3  ,(rlsc(247),x2rdmx  ),(rlsc(248),x3rdmx  ),(rlsc(249),theta   )
     4  ,(rlsc(250),thetamin),(rlsc(251),thetamax),(rlsc(252),dtheta  )
     5  ,(rlsc(253),phi     ),(rlsc(254),phimin  ),(rlsc(255),phimax  )
     6  ,(rlsc(256),dphi    ),(rlsc(257),alpha   ),(rlsc(258),freq    )
     7  ,(rlsc(259),dsw     ),(rlsc(260),btsw    ),(rlsc(261),dlower  )
     8  ,(rlsc(262),elower  ),(rlsc(263),selower ),(rlsc(264),pblower )
     9  ,(rlsc(265),shlower ),(rlsc(266),vvlower ),(rlsc(267),dupper  )
       equivalence
     1   (rlsc(268),eupper  ),(rlsc(269),seupper ),(rlsc(270),pbupper )
     2  ,(rlsc(271),shupper ),(rlsc(272),vvupper ),(rlsc(273),bmajor  )
     3  ,(rlsc(274),bminor  ),(rlsc(275),bpa     ),(rlsc(276),cpb     )
     4  ,(rlsc(277),icut    ),(rlsc(278),pcut    ),(rlsc(279),vlmin   )
     5  ,(rlsc(280),pscale  )
       equivalence
     1   (rlsc(331),tdmp    ),(rlsc(332),b1floor ),(rlsc(333),b2floor )
     2  ,(rlsc(334),b3floor ),(rlsc(335),dfloor  ),(rlsc(336),efloor  )
     3  ,(rlsc(337),v1floor ),(rlsc(338),v2floor ),(rlsc(339),v3floor )
     4  ,(rlsc(340),time    ),(rlsc(341),tlim    ),(rlsc(342),trem    )
     5  ,(rlsc(343),tsave   ),(rlsc(344),ttotal  ),(rlsc(345),tused   )
     6  ,(rlsc(346),dt      ),(rlsc(347),dtrat   ),(rlsc(348),dtmin   )
     7  ,(rlsc(349),dtnew   ),(rlsc(350),dtcs    ),(rlsc(351),dtal    )
     8  ,(rlsc(352),dtv1    ),(rlsc(353),dtv2    ),(rlsc(354),dtv3    )
     1  ,(rlsc(355),dtqq    ),(rlsc(356),dtqqi2  ),(rlsc(357),courno  )
       equivalence
     1   (rlsc(358),qcon    ),(rlsc(359),qlin    ),(rlsc(360),omega   )
     2  ,(rlsc(361),thetai  ),(rlsc(362),thetaj  ),(rlsc(363),thetak  )
     3  ,(rlsc(364),rmode   ),(rlsc(365),mamb    ),(rlsc(366),zsigma  )
     4  ,(rlsc(367),wnz     ),(rlsc(368),ojet    ),(rlsc(369),sig1    )
     5  ,(rlsc(370),sig2    ),(rlsc(371),sig3    ),(rlsc(372),xibgen  )
     6  ,(rlsc(373),rjet    ),(rlsc(374),g       ),(rlsc(375),ptmass  )
     6  ,(rlsc(376),x1ptm   ),(rlsc(377),x2ptm   ),(rlsc(378),x3ptm   )
       equivalence
     1   (rlsc(451),ttsl    ),(rlsc(452),ttslp   ),(rlsc(453),dttsl   )
     2  ,(rlsc(454),dttslp  ),(rlsc(455),tslpmn  ),(rlsc(456),tslpmx  )
       equivalence
     1   (rlsc(481),tusr    ),(rlsc(482),dtusr   )
       equivalence
     1   (rlsc(511),tvox    ),(rlsc(512),dtvox   ),(rlsc(513),x1vxmn  )
     2  ,(rlsc(514),x1vxmx  ),(rlsc(515),x2vxmn  ),(rlsc(516),x2vxmx  )
     3  ,(rlsc(517),x3vxmn  ),(rlsc(518),x3vxmx  )
c
c----------------------------- 1-D arrays ------------------------------
c
       character*2   c2niov  (niov,nt03)
       integer       inniov  (niov,nt04), innios  (nios,nt05)
     1             , inmn    (  mn,nt06)
       real*8        rlniov  (niov,nt07), rlnios  (nios,nt08)
     1             , rlmn    (  mn,nt09)
     2             , rlntsl  (ntsl,nt10)
c
       equivalence
     1    (c2niov  (1,  1), disvar  ), (c2niov  (1,  2), hdfvar  )
     2  , (c2niov  (1,  3), pixvar  ), (c2niov  (1,  4), plt1var )
     3  , (c2niov  (1,  5), plt2var ), (c2niov  (1,  6), radvar  )
     4  , (c2niov  (1,  7), voxvar  )
c
       equivalence
     1    (inniov  (1,  1), nlpix   ), (inniov  (1,  2), nlplt1  )
     2  , (inniov  (1,  3), nlplt2  ), (inniov  (1,  4), nlrad   )
     3  , (inniov  (1,  5), nlvox   )
       equivalence
     1    (innios  (1,  1), idisdir ), (innios  (1,  2), ldis    )
     2  , (innios  (1,  3), idismn  ), (innios  (1,  4), jdismn  )
     3  , (innios  (1,  5), kdismn  ), (innios  (1,  6), idismx  )
     4  , (innios  (1,  7), jdismx  ), (innios  (1,  8), kdismx  )
     5  , (innios  (1,  9), ipixdir ), (innios  (1, 10), lpix    )
     6  , (innios  (1, 11), npi     ), (innios  (1, 12), npj     )
     7  , (innios  (1, 13), iplt1   ), (innios  (1, 14), jplt1   )
     8  , (innios  (1, 15), kplt1   ), (innios  (1, 16), iplt1dir)
     9  , (innios  (1, 17), iplt2dir), (innios  (1, 18), lplt2   )
       equivalence
     1    (inmn    (1,  1), niso    )
c
       equivalence
     1    (rlniov  (1,  1), pixmin  ), (rlniov  (1,  2), pixmax  )
     2  , (rlniov  (1,  3), plt1min ), (rlniov  (1,  4), plt1max )
     3  , (rlniov  (1,  5), plt2min ), (rlniov  (1,  6), plt2max )
     4  , (rlniov  (1,  7), radmin  ), (rlniov  (1,  8), radmax  )
     5  , (rlniov  (1,  9), voxmin  ), (rlniov  (1, 10), voxmax  )
       equivalence
     1    (rlnios  (1,  1), x1pxmn  ), (rlnios  (1,  2), x2pxmn  )
     2  , (rlnios  (1,  3), x3pxmn  ), (rlnios  (1,  4), x1pxmx  )
     3  , (rlnios  (1,  5), x2pxmx  ), (rlnios  (1,  6), x3pxmx  )
     4  , (rlnios  (1,  7), x1p2mn  ), (rlnios  (1,  8), x2p2mn  )
     5  , (rlnios  (1,  9), x3p2mn  ), (rlnios  (1, 10), x1p2mx  )
     6  , (rlnios  (1, 11), x2p2mx  ), (rlnios  (1, 12), x3p2mx  )
       equivalence
     1    (rlmn    (1,  1), ciso    ), (rlmn    (1,  2), rmetal  )
     2  , (rlmn    (1,  3), gam     ), (rlmn    (1,  4), gamma   )
     3  , (rlmn    (1,  5), rgas    )
       equivalence
     1    (rlntsl  (1,  1), timetsl ), (rlntsl  (1,  2), massp   )
     2  , (rlntsl  (1,  3), x1cmp   ), (rlntsl  (1,  4), x2cmp   )
     3  , (rlntsl  (1,  5), x3cmp   ), (rlntsl  (1,  6), eintp   )
     4  , (rlntsl  (1,  7), ek1p    ), (rlntsl  (1,  8), ek2p    )
     5  , (rlntsl  (1,  9), ek3p    ), (rlntsl  (1, 10), ekp     )
     6  , (rlntsl  (1, 11), epotp   ), (rlntsl  (1, 12), tep     )
     7  , (rlntsl  (1, 13), l1p     ), (rlntsl  (1, 14), l2p     )
     8  , (rlntsl  (1, 15), l3p     ), (rlntsl  (1, 16), dminp   )
     9  , (rlntsl  (1, 17), dmaxp   ), (rlntsl  (1, 18), seminp  )
       equivalence
     1    (rlntsl  (1, 19), semaxp  ), (rlntsl  (1, 20), pminp   )
     2  , (rlntsl  (1, 21), pmaxp   ), (rlntsl  (1, 22), v1minp  )
     3  , (rlntsl  (1, 23), v1maxp  ), (rlntsl  (1, 24), v2minp  )
     4  , (rlntsl  (1, 25), v2maxp  ), (rlntsl  (1, 26), v3minp  )
     5  , (rlntsl  (1, 27), v3maxp  ), (rlntsl  (1, 28), vminp   )
     6  , (rlntsl  (1, 29), vmaxp   ), (rlntsl  (1, 30), dtcsp   )
     7  , (rlntsl  (1, 31), dtv1p   ), (rlntsl  (1, 32), dtv2p   )
     8  , (rlntsl  (1, 33), dtv3p   ), (rlntsl  (1, 34), dtqqp   )
     9  , (rlntsl  (1, 35), dtnewp  ), (rlntsl  (1, 36), emag1p  )
       equivalence
     1    (rlntsl  (1, 37), emag2p  ), (rlntsl  (1, 38), emag3p  )
     2  , (rlntsl  (1, 39), emagp   ), (rlntsl  (1, 40), b1flxp  )
     3  , (rlntsl  (1, 41), b2flxp  ), (rlntsl  (1, 42), b3flxp  )
     4  , (rlntsl  (1, 43), dvbp    ), (rlntsl  (1, 44), b1minp  )
     5  , (rlntsl  (1, 45), b1maxp  ), (rlntsl  (1, 46), b2minp  )
     6  , (rlntsl  (1, 47), b2maxp  ), (rlntsl  (1, 48), b3minp  )
     7  , (rlntsl  (1, 49), b3maxp  ), (rlntsl  (1, 50), bminp   )
     8  , (rlntsl  (1, 51), bmaxp   ), (rlntsl  (1, 52), dvbminp )
     9  , (rlntsl  (1, 53), dvbmaxp ), (rlntsl  (1, 54), dtalp   )
c
c----------------------------- 2-D arrays ------------------------------
c
       integer       injnkn  (  jn,  kn,nt11)
       real*8        rljnkn  (  jn,  kn,nt16)
c
       equivalence
     1   (injnkn(1,1,      1),niib    ),(injnkn(1,1,      2),noib    )
     2  ,(injnkn(1,1,      3),niib2   ),(injnkn(1,1,      4),noib2   )
     3  ,(injnkn(1,1,      5),niib3   ),(injnkn(1,1,      6),noib3   )
     4  ,(injnkn(1,1,      7),niib23  ),(injnkn(1,1,      8),noib23  )
       equivalence
     1   (rljnkn(1,1,      1),diib1   ),(rljnkn(1,1,      2),doib1   )
     2  ,(rljnkn(1,1,      3),diib2   ),(rljnkn(1,1,      4),doib2   )
     3  ,(rljnkn(1,1,      5),v1iib1  ),(rljnkn(1,1,      6),v1oib1  )
     4  ,(rljnkn(1,1,      7),v1iib2  ),(rljnkn(1,1,      8),v1oib2  )
     5  ,(rljnkn(1,1,      9),v2iib1  ),(rljnkn(1,1,     10),v2oib1  )
     6  ,(rljnkn(1,1,     11),v2iib2  ),(rljnkn(1,1,     12),v2oib2  )
     7  ,(rljnkn(1,1,     13),v3iib1  ),(rljnkn(1,1,     14),v3oib1  )
     8  ,(rljnkn(1,1,     15),v3iib2  ),(rljnkn(1,1,     16),v3oib2  )
       equivalence
     1   (rljnkn(1,1,nt12+ 1),eiib1   ),(rljnkn(1,1,nt12+ 2),eoib1   )
     2  ,(rljnkn(1,1,nt12+ 3),eiib2   ),(rljnkn(1,1,nt12+ 4),eoib2   )
       equivalence
     1   (rljnkn(1,1,nt14+ 1),b1iib1  ),(rljnkn(1,1,nt14+ 2),b1oib1  )
     2  ,(rljnkn(1,1,nt14+ 3),b1iib2  ),(rljnkn(1,1,nt14+ 4),b1oib2  )
     3  ,(rljnkn(1,1,nt14+ 5),b2iib1  ),(rljnkn(1,1,nt14+ 6),b2oib1  )
     4  ,(rljnkn(1,1,nt14+ 7),b2iib2  ),(rljnkn(1,1,nt14+ 8),b2oib2  )
     5  ,(rljnkn(1,1,nt14+ 9),b3iib1  ),(rljnkn(1,1,nt14+10),b3oib1  )
     6  ,(rljnkn(1,1,nt14+11),b3iib2  ),(rljnkn(1,1,nt14+12),b3oib2  )
     7  ,(rljnkn(1,1,nt14+13),emf1iib1),(rljnkn(1,1,nt14+14),emf1oib1)
     8  ,(rljnkn(1,1,nt14+15),emf1iib2),(rljnkn(1,1,nt14+16),emf1oib2)
     8  ,(rljnkn(1,1,nt14+17),emf1iib3),(rljnkn(1,1,nt14+18),emf1oib3)
       equivalence
     9   (rljnkn(1,1,nt14+19),emf2iib1),(rljnkn(1,1,nt14+20),emf2oib1)
     1  ,(rljnkn(1,1,nt14+21),emf2iib2),(rljnkn(1,1,nt14+22),emf2oib2)
     2  ,(rljnkn(1,1,nt14+23),emf2iib3),(rljnkn(1,1,nt14+24),emf2oib3)
     3  ,(rljnkn(1,1,nt14+25),emf3iib1),(rljnkn(1,1,nt14+26),emf3oib1)
     4  ,(rljnkn(1,1,nt14+27),emf3iib2),(rljnkn(1,1,nt14+28),emf3oib2)
     5  ,(rljnkn(1,1,nt14+29),emf3iib3),(rljnkn(1,1,nt14+30),emf3oib3)
     6  ,(rljnkn(1,1,nt14+31),a1iib1  ),(rljnkn(1,1,nt14+32),a1iib2  )
     7  ,(rljnkn(1,1,nt14+33),a1iib3  ),(rljnkn(1,1,nt14+34),a1iib4  )
     1  ,(rljnkn(1,1,nt14+35),a2iib1  ),(rljnkn(1,1,nt14+36),a2iib2  )
       equivalence
     8   (rljnkn(1,1,nt14+37),a2iib3  ),(rljnkn(1,1,nt14+38),a2iib4  )
     9  ,(rljnkn(1,1,nt14+39),a3iib1  ),(rljnkn(1,1,nt14+40),a3iib2  )
     1  ,(rljnkn(1,1,nt14+41),a3iib3  ),(rljnkn(1,1,nt14+42),a3iib4  )
c
       integer       inknin
       real*8        rlknin
c
       integer       ininjn  (  in,  jn,nt11)
       real*8        rlinjn  (  in,  jn,nt15)
c
       equivalence
     1   (ininjn(1,1,      1),nikb    ),(ininjn(1,1,      2),nokb    )
     2  ,(ininjn(1,1,      3),nikb1   ),(ininjn(1,1,      4),nokb1   )
     3  ,(ininjn(1,1,      5),nikb2   ),(ininjn(1,1,      6),nokb2   )
     4  ,(ininjn(1,1,      7),nikb12  ),(ininjn(1,1,      8),nokb12  )
       equivalence
     1   (rlinjn(1,1,      1),dikb1   ),(rlinjn(1,1,      2),dokb1   )
     2  ,(rlinjn(1,1,      3),dikb2   ),(rlinjn(1,1,      4),dokb2   )
     3  ,(rlinjn(1,1,      5),v1ikb1  ),(rlinjn(1,1,      6),v1okb1  )
     4  ,(rlinjn(1,1,      7),v1ikb2  ),(rlinjn(1,1,      8),v1okb2  )
     5  ,(rlinjn(1,1,      9),v2ikb1  ),(rlinjn(1,1,     10),v2okb1  )
     6  ,(rlinjn(1,1,     11),v2ikb2  ),(rlinjn(1,1,     12),v2okb2  )
     7  ,(rlinjn(1,1,     13),v3ikb1  ),(rlinjn(1,1,     14),v3okb1  )
     8  ,(rlinjn(1,1,     15),v3ikb2  ),(rlinjn(1,1,     16),v3okb2  )
       equivalence
     1   (rlinjn(1,1,nt12+ 1),eikb1   ),(rlinjn(1,1,nt12+ 2),eokb1   )
     2  ,(rlinjn(1,1,nt12+ 3),eikb2   ),(rlinjn(1,1,nt12+ 4),eokb2   )
       equivalence
     1   (rlinjn(1,1,nt14+ 1),b1ikb1  ),(rlinjn(1,1,nt14+ 2),b1okb1  )
     2  ,(rlinjn(1,1,nt14+ 3),b1ikb2  ),(rlinjn(1,1,nt14+ 4),b1okb2  )
     3  ,(rlinjn(1,1,nt14+ 5),b2ikb1  ),(rlinjn(1,1,nt14+ 6),b2okb1  )
     4  ,(rlinjn(1,1,nt14+ 7),b2ikb2  ),(rlinjn(1,1,nt14+ 8),b2okb2  )
     5  ,(rlinjn(1,1,nt14+ 9),b3ikb1  ),(rlinjn(1,1,nt14+10),b3okb1  )
     6  ,(rlinjn(1,1,nt14+11),b3ikb2  ),(rlinjn(1,1,nt14+12),b3okb2  )
     7  ,(rlinjn(1,1,nt14+13),emf1ikb1),(rlinjn(1,1,nt14+14),emf1okb1)
     8  ,(rlinjn(1,1,nt14+15),emf1ikb2),(rlinjn(1,1,nt14+16),emf1okb2)
     9  ,(rlinjn(1,1,nt14+17),emf1ikb3),(rlinjn(1,1,nt14+18),emf1okb3)
       equivalence
     1   (rlinjn(1,1,nt14+19),emf2ikb1),(rlinjn(1,1,nt14+20),emf2okb1)
     2  ,(rlinjn(1,1,nt14+21),emf2ikb2),(rlinjn(1,1,nt14+22),emf2okb2)
     3  ,(rlinjn(1,1,nt14+23),emf2ikb3),(rlinjn(1,1,nt14+24),emf2okb3)
     4  ,(rlinjn(1,1,nt14+25),emf3ikb1),(rlinjn(1,1,nt14+26),emf3okb1)
     5  ,(rlinjn(1,1,nt14+27),emf3ikb2),(rlinjn(1,1,nt14+28),emf3okb2)
     5  ,(rlinjn(1,1,nt14+29),emf3ikb3),(rlinjn(1,1,nt14+30),emf3okb3)
c
c----------------------------- 3-D arrays ------------------------------
c
       real*8        rlinjnkn(  in,  jn,  kn,nt20)
c
       equivalence
     1    ( rlinjnkn(1,1,1,      1), d        )
     2  , ( rlinjnkn(1,1,1,      2), v1       )
     3  , ( rlinjnkn(1,1,1,      3), v2       )
     4  , ( rlinjnkn(1,1,1,      4), v3       )
     5  , ( rlinjnkn(1,1,1,nt17+ 1), e        )
     7  , ( rlinjnkn(1,1,1,nt19+ 1), b1       )
     8  , ( rlinjnkn(1,1,1,nt19+ 2), b2       )
     9  , ( rlinjnkn(1,1,1,nt19+ 3), b3       )
c
c============================= COMMON BLOCKS ===========================
c
       common / rootc    /
     1               c2niov
       common / rooti    /
     1               insc    , inniov  , innios  , inmn    , injnkn
     2             , inknin  , ininjn
       common / rootr    /
     1               rlsc    , rlniov  , rlnios  , rlmn    , rlntsl
     2             , x1a     , x2a     , x3a     , rljnkn  , rlknin
     3             , rlinjn  , rlinjnkn
c
c-----------------------------------------------------------------------
c-------------------------- COOLING VARIABLES --------------------------
c-----------------------------------------------------------------------
c
       real*8        tfac    , tmin    , trat
       real*8        cfunc   ( nmd), dcfnc   ( nmd)
       real*8        xheat   (  in)
c
c----------------------------- Common blocks ---------------------------
c
       common / coolcomr /
     1               tfac    , tmin    , trat    , cfunc   , dcfnc
     2             , xheat
c
c-----------------------------------------------------------------------
c---------------------- ADDITIONAL GRID VARIABLES ----------------------
c-----------------------------------------------------------------------
c
       integer       ismnm2  , ismnm1  , ismnp1  , ismnp2  , ismnp3
     1             , ism2    , ism1    , isp1    , isp2    , isp3
     2             , iemxm2  , iemxm1  , iemxp1  , iemxp2  , iemxp3
     3             , iem2    , iem1    , iep1    , iep2    , iep3
     4             , nx1zmx  , nx1z
       integer       jsmnm2  , jsmnm1  , jsmnp1  , jsmnp2  , jsmnp3
     1             , jsm2    , jsm1    , jsp1    , jsp2    , jsp3
     2             , jemxm2  , jemxm1  , jemxp1  , jemxp2  , jemxp3
     3             , jem2    , jem1    , jep1    , jep2    , jep3
     4             , nx2zmx  , nx2z
       integer       ksmnm2  , ksmnm1  , ksmnp1  , ksmnp2  , ksmnp3
     1             , ksm2    , ksm1    , ksp1    , ksp2    , ksp3
     2             , kemxm2  , kemxm1  , kemxp1  , kemxp2  , kemxp3
     3             , kem2    , kem1    , kep1    , kep2    , kep3
     4             , nx3zmx  , nx3z
       real*8        onei    , onej    , onek
c
       real*8                        x1ah    (  in), x1an    (  in)
     1             , dx1a    (  in), dx1ah   (  in), dx1an   (  in)
     2             , dvl1a   (  in), dvl1ah  (  in), dvl1an  (  in)
     3             , g2a     (  in), g2ah    (  in), g2an    (  in)
     4             , g31a    (  in), g31ah   (  in), g31an   (  in)
     5             , dg2ad1  (  in)
     6             , dg31ad1 (  in)
       real*8        x1ai    (  in), x1ahi   (  in), x1ani   (  in)
     1             , dx1ai   (  in), dx1ahi  (  in), dx1ani  (  in)
     2             , dvl1ai  (  in), dvl1ahi (  in), dvl1ani (  in)
     3             , g2ai    (  in), g2ahi   (  in), g2ani   (  in)
     4             , g31ai   (  in), g31ahi  (  in), g31ani  (  in)
       real*8        x1b     (  in), x1bh    (  in), x1bn    (  in)
     1             , dx1b    (  in), dx1bh   (  in), dx1bn   (  in)
     2             , dvl1b   (  in), dvl1bh  (  in), dvl1bn  (  in)
     3             , g2b     (  in), g2bh    (  in), g2bn    (  in)
     4             , g31b    (  in), g31bh   (  in), g31bn   (  in)
     5             , dg2bd1  (  in)
     6             , dg31bd1 (  in)
       real*8        x1bi    (  in), x1bhi   (  in), x1bni   (  in)
     1             , dx1bi   (  in), dx1bhi  (  in), dx1bni  (  in)
     2             , dvl1bi  (  in), dvl1bhi (  in), dvl1bni (  in)
     3             , g2bi    (  in), g2bhi   (  in), g2bni   (  in)
     4             , g31bi   (  in), g31bhi  (  in), g31bni  (  in)
       real*8        vg1     (  in)
       real*8        f1dx1a  (  in), f2dx1a  (  in), f3dx1a  (  in)
     1             , f4dx1a  (  in), f5dx1a  (  in), f6dx1a  (  in)
     2             , f7dx1a  (  in), f8dx1a  (  in)
       real*8        f1dx1b  (  in), f2dx1b  (  in), f3dx1b  (  in)
     1             , f4dx1b  (  in), f5dx1b  (  in), f6dx1b  (  in)
     2             , f7dx1b  (  in), f8dx1b  (  in), f9dx1b  (  in)
c
       real*8                        x2ah    (  jn), x2an    (  jn)
     1             , dx2a    (  jn), dx2ah   (  jn), dx2an   (  jn)
     2             , dvl2a   (  jn), dvl2ah  (  jn), dvl2an  (  jn)
     3             , g32a    (  jn), g32ah   (  jn), g32an   (  jn)
     4             , dg32ad2 (  jn)
       real*8        x2ai    (  jn), x2ahi   (  jn), x2ani   (  jn)
     1             , dx2ai   (  jn), dx2ahi  (  jn), dx2ani  (  jn)
     2             , dvl2ai  (  jn), dvl2ahi (  jn), dvl2ani (  jn)
     3             , g32ai   (  jn), g32ahi  (  jn), g32ani  (  jn)
       real*8        x2b     (  jn), x2bh    (  jn), x2bn    (  jn)
     1             , dx2b    (  jn), dx2bh   (  jn), dx2bn   (  jn)
     2             , dvl2b   (  jn), dvl2bh  (  jn), dvl2bn  (  jn)
     3             , g32b    (  jn), g32bh   (  jn), g32bn   (  jn)
     4             , dg32bd2 (  jn)
       real*8        x2bi    (  jn), x2bhi   (  jn), x2bni   (  jn)
     1             , dx2bi   (  jn), dx2bhi  (  jn), dx2bni  (  jn)
     2             , dvl2bi  (  jn), dvl2bhi (  jn), dvl2bni (  jn)
     3             , g32bi   (  jn), g32bhi  (  jn), g32bni  (  jn)
       real*8        vg2     (  jn)
       real*8        f1dx2a  (  jn), f2dx2a  (  jn), f3dx2a  (  jn)
     1             , f4dx2a  (  jn), f5dx2a  (  jn), f6dx2a  (  jn)
     2             , f7dx2a  (  jn), f8dx2a  (  jn)
       real*8        f1dx2b  (  jn), f2dx2b  (  jn), f3dx2b  (  jn)
     1             , f4dx2b  (  jn), f5dx2b  (  jn), f6dx2b  (  jn)
     2             , f7dx2b  (  jn), f8dx2b  (  jn), f9dx2b  (  jn)
c
       real*8                        x3ah    (  kn), x3an    (  kn)
     1             , dx3a    (  kn), dx3ah   (  kn), dx3an   (  kn)
     2             , dvl3a   (  kn), dvl3ah  (  kn), dvl3an  (  kn)
       real*8        x3ai    (  kn), x3ahi   (  kn), x3ani   (  kn)
     1             , dx3ai   (  kn), dx3ahi  (  kn), dx3ani  (  kn)
     2             , dvl3ai  (  kn), dvl3ahi (  kn), dvl3ani (  kn)
       real*8        x3b     (  kn), x3bh    (  kn), x3bn    (  kn)
     1             , dx3b    (  kn), dx3bh   (  kn), dx3bn   (  kn)
     2             , dvl3b   (  kn), dvl3bh  (  kn), dvl3bn  (  kn)
       real*8        x3bi    (  kn), x3bhi   (  kn), x3bni   (  kn)
     1             , dx3bi   (  kn), dx3bhi  (  kn), dx3bni  (  kn)
     2             , dvl3bi  (  kn), dvl3bhi (  kn), dvl3bni (  kn)
       real*8        vg3     (  kn)
       real*8        f1dx3a  (  kn), f2dx3a  (  kn), f3dx3a  (  kn)
     1             , f4dx3a  (  kn), f5dx3a  (  kn), f6dx3a  (  kn)
     2             , f7dx3a  (  kn), f8dx3a  (  kn)
       real*8        f1dx3b  (  kn), f2dx3b  (  kn), f3dx3b  (  kn)
     1             , f4dx3b  (  kn), f5dx3b  (  kn), f6dx3b  (  kn)
     2             , f7dx3b  (  kn), f8dx3b  (  kn), f9dx3b  (  kn)
c
c----------------------------- Common blocks ---------------------------
c
       common / gridcomi /
     1               ismnm2  , ismnm1  , ismnp1  , ismnp2  , ismnp3
     2             , ism2    , ism1    , isp1    , isp2    , isp3
     3             , iemxm2  , iemxm1  , iemxp1  , iemxp2  , iemxp3
     4             , iem2    , iem1    , iep1    , iep2    , iep3
     5             , nx1zmx  , nx1z
       common / gridcomi /
     1               jsmnm2  , jsmnm1  , jsmnp1  , jsmnp2  , jsmnp3
     2             , jsm2    , jsm1    , jsp1    , jsp2    , jsp3
     3             , jemxm2  , jemxm1  , jemxp1  , jemxp2  , jemxp3
     4             , jem2    , jem1    , jep1    , jep2    , jep3
     5             , nx2zmx  , nx2z
       common / gridcomi /
     1               ksmnm2  , ksmnm1  , ksmnp1  , ksmnp2  , ksmnp3
     2             , ksm2    , ksm1    , ksp1    , ksp2    , ksp3
     3             , kemxm2  , kemxm1  , kemxp1  , kemxp2  , kemxp3
     4             , kem2    , kem1    , kep1    , kep2    , kep3
     5             , nx3zmx  , nx3z
c
       common / gridcomr /
     1               onei    , onej    , onek
c
       common / gridcomr /
     1               x1ah    , x1an    , dx1a    , dx1ah   , dx1an
     2             , dvl1a   , dvl1ah  , dvl1an  , g2a     , g2ah
     3             , g2an    , g31a    , g31ah   , g31an   , dg2ad1
     4             , dg31ad1 , x1ai    , x1ahi   , x1ani   , dx1ai
     5             , dx1ahi  , dx1ani  , dvl1ai  , dvl1ahi , dvl1ani
     6             , g2ai    , g2ahi   , g2ani   , g31ai   , g31ahi
     7             , g31ani  , x1b     , x1bh    , x1bn    , dx1b
     8             , dx1bh   , dx1bn   , dvl1b   , dvl1bh  , dvl1bn
     9             , g2b     , g2bh    , g2bn    , g31b    , g31bh
       common / gridcomr /
     1               g31bn   , dg2bd1  , dg31bd1 , x1bi    , x1bhi
     2             , x1bni   , dx1bi   , dx1bhi  , dx1bni  , dvl1bi
     3             , dvl1bhi , dvl1bni , g2bi    , g2bhi   , g2bni
     4             , g31bi   , g31bhi  , g31bni  , vg1     , f1dx1a
     5             , f2dx1a  , f3dx1a  , f4dx1a  , f5dx1a  , f6dx1a
     6             , f7dx1a  , f8dx1a  , f1dx1b  , f2dx1b  , f3dx1b
     7             , f4dx1b  , f5dx1b  , f6dx1b  , f7dx1b  , f8dx1b
     8             , f9dx1b
c
       common / gridcomr /
     1               x2ah    , x2an    , dx2a    , dx2ah   , dx2an
     2             , dvl2a   , dvl2ah  , dvl2an  , g32a    , g32ah
     3             , g32an   , dg32ad2 , x2ai    , x2ahi   , x2ani
     4             , dx2ai   , dx2ahi  , dx2ani  , dvl2ai  , dvl2ahi
     5             , dvl2ani , g32ai   , g32ahi  , g32ani  , x2b
     6             , x2bh    , x2bn    , dx2b    , dx2bh   , dx2bn
     7             , dvl2b   , dvl2bh  , dvl2bn  , g32b    , g32bh
     8             , g32bn   , dg32bd2 , x2bi    , x2bhi   , x2bni
     9             , dx2bi   , dx2bhi  , dx2bni  , dvl2bi  , dvl2bhi
       common / gridcomr /
     1               dvl2bni , g32bi   , g32bhi  , g32bni  , vg2
     2             , f1dx2a  , f2dx2a  , f3dx2a  , f4dx2a  , f5dx2a
     3             , f6dx2a  , f7dx2a  , f8dx2a  , f1dx2b  , f2dx2b
     4             , f3dx2b  , f4dx2b  , f5dx2b  , f6dx2b  , f7dx2b
     5             , f8dx2b  , f9dx2b
c
       common / gridcomr /
     1               x3ah    , x3an    , dx3a    , dx3ah   , dx3an
     2             , dvl3a   , dvl3ah  , dvl3an  , x3ai    , x3ahi
     3             , x3ani   , dx3ai   , dx3ahi  , dx3ani  , dvl3ai
     4             , dvl3ahi , dvl3ani , x3b     , x3bh    , x3bn
     5             , dx3b    , dx3bh   , dx3bn   , dvl3b   , dvl3bh
     6             , dvl3bn  , x3bi    , x3bhi   , x3bni   , dx3bi
     7             , dx3bhi  , dx3bni  , dvl3bi  , dvl3bhi , dvl3bni
     8             , vg3     , f1dx3a  , f2dx3a  , f3dx3a  , f4dx3a
     9             , f5dx3a  , f6dx3a  , f7dx3a  , f8dx3a  , f1dx3b
       common / gridcomr /
     1               f2dx3b  , f3dx3b  , f4dx3b  , f5dx3b  , f6dx3b
     2             , f7dx3b  , f8dx3b  , f9dx3b
c
c-----------------------------------------------------------------------
c--------------------- ADDITIONAL RADIO VARIABLES ----------------------
c-----------------------------------------------------------------------
c
c
c-----------------------------------------------------------------------
c--------------------- ADDITIONAL SUNDRY VARIABLES ---------------------
c-----------------------------------------------------------------------
c
       character*2   id
       character*8   logfile , resfile , plt1file, plt2file, tslfile
     1             , tslpfile, disfile , usrfile
       integer       ioin    , iotty   , ioplt1  , ioplt2  , iolog
     1             , iodmp   , iopix   , iousr   , iotsl   , iotslp
     2             , iovox   , iodis   , iorad   , igetmn  , igetmx
     3             , jgetmn  , jgetmx  , kgetmn  , kgetmx  , iputmn
     4             , iputmx  , jputmn  , jputmx  , kputmn  , kputmx
     5             , iaddz   , jaddz   , kaddz   , irestart
       real*8        dtdmp
c
c----------------------------- Common blocks ---------------------------
c
       common / sundcomc /
     1               id      , logfile , resfile , plt1file, plt2file
     2             , tslfile , tslpfile, disfile , usrfile
       common / sundcomi /
     1               ioin    , iotty   , ioplt1  , ioplt2  , iolog
     2             , iodmp   , iopix   , iousr   , iotsl   , iotslp
     3             , iovox   , iodis   , iorad   , igetmn  , igetmx
     4             , jgetmn  , jgetmx  , kgetmn  , kgetmx  , iputmn
     5             , iputmx  , jputmn  , jputmx  , kputmn  , kputmx
     6             , iaddz   , jaddz   , kaddz   , irestart
       common / sundcomr /
     1               dtdmp
c
c-----------------------------------------------------------------------
c-------------------- SCRATCH SPACE (WORKER ARRAYS) --------------------
c-----------------------------------------------------------------------
c
c      1-D worker arrays
c
       real*8        wa1d    (ijkn), wb1d    (ijkn), wc1d    (ijkn)
     1             , wd1d    (ijkn), we1d    (ijkn), wf1d    (ijkn)
     2             , wg1d    (ijkn), wh1d    (ijkn), wi1d    (ijkn)
     3             , wj1d    (ijkn), wk1d    (ijkn), wl1d    (ijkn)
     4             , wm1d    (ijkn), wn1d    (ijkn), wo1d    (ijkn)
     5             , wp1d    (ijkn), wq1d    (ijkn), wr1d    (ijkn)
     6             , ws1d    (ijkn), wt1d    (ijkn), wu1d    (ijkn)
     7             , wv1d    (ijkn), ww1d    (ijkn), wx1d    (ijkn)
     8             , wy1d    (ijkn), wz1d    (ijkn)
c
c      2-D worker arrays
c
       real*8        wa2d    (idim,jdim), wa2dt   (jdim,idim)
     1             , wb2d    (idim,jdim), wb2dt   (jdim,idim)
     1             , wc2d    (idim,jdim), wc2dt   (jdim,idim)
     2             , wd2d    (idim,jdim), wd2dt   (jdim,idim)
     1             , we2d    (idim,jdim), we2dt   (jdim,idim)
     2             , wf2d    (idim,jdim), wf2dt   (jdim,idim)
c
c      3-D worker arrays
c
       real*8        wa3d    (  in,  jn,  kn), wb3d    (  in,  jn,  kn)
     1             , wc3d    (  in,  jn,  kn), wd3d    (  in,  jn,  kn)
     2             , we3d    (  in,  jn,  kn), wf3d    (  in,  jn,  kn)
     3             , wg3d    (  in,  jn,  kn)
c
       equivalence   ( wa2d    , wa2dt    ), ( wb2d    , wb2dt    )
     1             , ( wc2d    , wc2dt    ), ( wd2d    , wd2dt    )
     2             , ( we2d    , we2dt    ), ( wf2d    , wf2dt    )
       equivalence   ( wa3d    , wa2d     ), ( wb3d    , wb2d     )
     1             , ( wc3d    , wc2d     ), ( wd3d    , wd2d     )
     2             , ( we3d    , we2d     ), ( wf3d    , wf2d     )
c
c---------------------------- Common Blocks ----------------------------
c
       common / scratchr /
     1               wa1d    , wb1d    , wc1d    , wd1d    , we1d
     2             , wf1d    , wg1d    , wh1d    , wi1d    , wj1d
     3             , wk1d    , wl1d    , wm1d    , wn1d    , wo1d
     4             , wp1d    , wq1d    , wr1d    , ws1d    , wt1d
     5             , wu1d    , wv1d    , ww1d    , wx1d    , wy1d
     6             , wz1d
       common / scratchr /
     1               wa3d    , wb3d    , wc3d    , wd3d    , we3d
     2             , wf3d
     3             , wg3d
c
       integer       mgit
       real*8        gerr
       parameter     ( mgit = 20, gerr = 1.0e-6 )
c
       logical       lgrid
       integer       nbl     , igrid   , imin    , imax    , ibeg
     1             , iend    , iter    , i
       real*8        x1min   , x1max   , x1rat   , dx1min  , dfndx1r
     1             , x1r     , deltx1r , gerrx1r
       real*8        fn      , small
       character*32  namelist (   2,1000)
       external      nlsdac01, nlsdac02, nlsdac03, nlsdac04, nlsdac05
     1             , nlsdac06, nlsdac07, nlsdac08, nlsdac09, nlsdac10
     2             , nlsdac11, nlsdac12, nlsdac13, nlsdac14, nlsdac15
     3             , nlsdac16, nlsdac17, nlsdac18, nlsdac19, nlsdac20
     4             , nlsdac21, nlsdac22, nlsdac23, nlsdac24
c
c      Data statements
c
       data          small    / 0.001 /
c
c      External statements
c
       real*8        modulo
       external      round   , modulo
c
c-----------------------------------------------------------------------
c
       if (in .le. 6) then
c
c      Set up a token x1-grid with one zone for problems with 1-symmetry
c
         ismn   = 1
         iemx   = 1
         ione   = 0
         x1a(1) = 0.0
         return
       endif
c
c      Set up a full x1-grid for problems where in > 6.  Read in blocks
c  of x1 grid zones.  Loop over read statement until all blocks are read
c  (signalled by reading in "lgrid" = .true.).  Each block can be zoned
c  independently of the others.  The only constraint is that the
c  starting position of the current block ("x1min") must correspond to
c  the ending position ("x1max") of the previous block.  Over the
c  current block, the zone widths are either ratioed or scaled, as
c  explained below.
c
c      Case 1:  RATIOED ZONES
c
c      Ratioed zones means that the ratio dx(i+1)/dx(i) is constant over
c  a specified region of the grid.  The zones can then be computed using
c  the zoning equation:
c
c      x1(I) = x1min + SUM OVER i UP TO I [dx1min*x1rat**i]
c
c  Not all parameters are independent.  Since x1max is specified, then
c
c      x1max = x1(N) = x1min + SUM OVER i UP TO N [dx1min*x1rat**i]
c
c  will allow "dx1min" to be computed from "x1rat" ("igrid" = +1, -1),
c  or "x1rat" to be computed from "dx1min" ("igrid" = +2, -2).
c
c      Case 2:  SCALED ZONES
c
c      Scaled zones are zones where dx is proportional to x. Thus, zones
c  with uniform shapes in spherical coordinates can be achieved with
c  this option.  In this case, the ratio x(i+1)/x(i) is constant over
c  the specified region of the grid.  The zones are computed using the
c  zoning equation:
c
c      x1(I) = x1min * x1rat**(I)
c
c  Since x1max = x1(N) is specified, N and "x1rat" are not independent.
c
c      x1max = x1min * x1rat**(N)
c
c  Thus, in this case ("igrid" = 3) the grid can be generated from N,
c  "x1min", and "x1max" and the user specified values of "x1rat" and
c  "dx1min" are ignored.
c
c   nbl        number of active zones in block being generated
c   x1min      x1a(imin), bottom position of block
c   x1max      x1a(imax), top    position of block
c   x1scale    arbitrary scaling factor for x1min and x1max
c   igrid      selects zoning type.
c              = 0  => block has already been set (irestart=1)
c              =+1  => (ratioed) use input "x1rat" to compute "dx1min",
c                      where x1a(imin+1) = x1a(imin) + dx1min.
c              =-1  => (ratioed) use input "x1rat" to compute "dx1min",
c                      where x1a(imax) = x1a(imax-1) + dx1min.
c              =+2  => (ratioed) use input "dx1min" to compute "x1rat",
c                      where x1a(imin+1) = x1a(imin) + dx1min.
c              =-2  => (ratioed) use input "dx1min" to compute "x1rat",
c                      where x1a(imax) = x1a(imax-1) + dx1min.
c              = 3  => (scaled) compute "x1rat" and "dx1min" from "nbl".
c   lgrid      =.false. => read another block (namelist card).
c              =.true.  => all blocks are read in.  Do not look for
c                          another namelist card.
c   imax,imin  indices of top and bottom of block
c
       ismn   = 3
       imax   = ismn - 1
       call nlsdac01 ( namelist, 1000 )
       namelist (1,   1) = 'nbl                          i s'
       namelist (1,   2) = 'x1min                        r s'
       namelist (1,   3) = 'x1max                        r s'
       namelist (1,   4) = 'x1scale                      r s'
       namelist (1,   5) = 'igrid                        i s'
       namelist (1,   6) = 'x1rat                        r s'
       namelist (1,   7) = 'dx1min                       r s'
       namelist (1,   8) = 'lgrid                        l s'
       nbl     = 2
       x1min   = 0.0
       x1max   = 0.0
       x1scale = 1.0
       igrid   = 1
       x1rat   = 1.0
       dx1min  = 0.0
       lgrid   = .false.
c
10     continue
       call nlsdac02 ( ioin    , 'ggen1   ', namelist, 1000,    8 )
       call nlsdac04 ( nbl     , 'nbl     ', namelist, 1000 )
       call nlsdac05 ( x1min   , 'x1min   ', namelist, 1000 )
       call nlsdac05 ( x1max   , 'x1max   ', namelist, 1000 )
       call nlsdac05 ( x1scale , 'x1scale ', namelist, 1000 )
       call nlsdac04 ( igrid   , 'igrid   ', namelist, 1000 )
       call nlsdac05 ( x1rat   , 'x1rat   ', namelist, 1000 )
       call nlsdac05 ( dx1min  , 'dx1min  ', namelist, 1000 )
       call nlsdac06 ( lgrid   , 'lgrid   ', namelist, 1000 )
       if (iolog .gt. 0) then
       call nlsdac24 ( iolog   , 'ggen1   ', namelist, 1000 )
       endif
c
       imin   = imax + 1
       imax   = imax + nbl
       if (imax .gt. in-3) then
         if (iolog .gt. 0) write (iolog, 2010) imax, in
         if (iotty .gt. 0) write (iotty, 2010) imax, in
         stop
       endif
c
c      Set limits of do-loops so that boundary zones will be computed
c  with the active zones.
c
       ibeg = imin
       if (imin .eq. ismn) ibeg = imin - 2
       iend = imax + 1
       if (lgrid)          iend = imax + 3
c
c  IGRID = 0:  This portion of the grid has already been set (this is
c              a restarted job: irestart = 1).  Write grid to log file
c              (zlnnid) and go on to next block.
c
       if (igrid .eq. 0) then
         if (iolog .gt. 0) then
           write (iolog, 2050) ibeg, iend
           write (iolog, 2060) resfile
           write (iolog, 2100) (x1a(i),i=ibeg,iend)
         endif
         if (.not. lgrid) go to 10
         go to 80
       endif
c
c  IGRID = 1:  For ratioed zones, compute "dx1min" from given value of
c              "x1rat".
c
       if ( abs(igrid) .eq. 1)  then
         if (igrid .eq. -1) x1rat = 1.0 / x1rat
         if (x1rat .eq. 1.0)  then
           dx1min = ( x1max - x1min ) /  real(nbl)
         else
           dx1min = ( x1max - x1min ) * ( x1rat      - 1.0 )
     1                                / ( x1rat**nbl - 1.0 )
         endif
       endif
c
c  IGRID = 2:  For ratioed zones, compute "x1rat" from given value of
c              "dx1min".
c
c  A Newton Raphson iteration is used to find the root ("x1rat") of the
c  function:
c
c    fn(x1r) = (x1max-x1min) - dx1min * [(x1r)**nbl - 1] / [x1r - 1] = 0
c
c  dfndx1r = derivative of fn wrt x1r
c  deltx1r = correction term to x1r
c  gerrx1r = relative error of correction
c
       if ( abs(igrid) .eq. 2) then
         x1r = 1.01
         do 20 iter=1,mgit
           fn      = x1max - x1min
     1             - dx1min * ( x1r**nbl - 1.0 ) / ( x1r - 1.0 )
           dfndx1r = dx1min * ( x1r**nbl - 1.0 ) / ( x1r - 1.0 )**2
     1             - nbl * dx1min * x1r**(nbl-1) / ( x1r - 1.0 )
           deltx1r = - fn / dfndx1r
           gerrx1r = abs ( deltx1r / x1r)
           x1r     = x1r + deltx1r
           if (gerrx1r .lt. gerr) go to 30
20       continue
         if (iolog .gt. 0) write (iolog, 2020) imin, x1r, deltx1r, fn
         if (iotty .gt. 0) write (iotty, 2020) imin, x1r, deltx1r, fn
         stop
c
30       continue
         x1rat = x1r
       endif
       if (igrid .eq. -2) then
         dx1min = dx1min * x1rat**(nbl-1)
         x1rat  = 1.0 / x1rat
       endif
c
c  IGRID = 3:  For scaled zones, compute "x1rat", the ratio of adjacent
c              zone-centre coordinates.
c
       if (igrid .eq. 3) then
         if (x1min .le. 0.0) then
           if (iolog .gt. 0) write (iolog, 2030) x1min
           if (iotty .gt. 0) write (iotty, 2030) x1min
           stop
         endif
         x1rat = ( x1max / x1min )**( 1.0 /  real(nbl) )
       endif
c
c      Default = error in input of "igrid".  Issue error message and
c  terminate execution (only legal values of "igrid" are -2, -1, 0, 1,
c  2, and 3).
c
       if ( (igrid .lt. -2) .or. (igrid .gt. 3) ) then
         if (iolog .gt. 0) write (iolog, 2040) imin, igrid
         if (iotty .gt. 0) write (iotty, 2040) imin, igrid
         stop
       endif
c
c      Set "x1a" grid from i = "ibeg" to "iend", using known values of
c  "x1min" and "x1rat".
c
       if (igrid .le. 2) then
         if (x1rat .eq. 1.0) then
           do 40 i=ibeg,iend
             x1a(i) = x1min + dx1min *  real(i-imin)
40         continue
         else
           do 50 i=ibeg,iend
             x1a(i) = x1min + dx1min * ( x1rat**(i-imin) - 1.0 )
     1                               / ( x1rat           - 1.0 )
50         continue
         endif
       endif
       if (igrid .eq. 3) then
         do 60 i=ibeg,iend
           x1a(i) = x1min * x1rat**(i-imin)
60       continue
       endif
c
c      Write grid to log file (zlnnid).
c
       if (iolog .gt. 0) then
         write (iolog, 2050) ibeg, iend
         if (x1rat .ge. 1.0) write (iolog, 2070) dx1min, x1rat
         if (x1rat .lt. 1.0) write (iolog, 2080) dx1min, x1rat
         write (iolog, 2090) x1scale
         write (iolog, 2100) (x1a(i),i=ibeg,iend)
       endif
c
c      Scale the grid.
c
       do 70 i=ibeg,iend
         x1a (i) = x1scale * x1a(i)
70     continue
c
c      Go back and read another block of x1 grid zones if "lgrid" is
c  still .false..
c
       if (.not. lgrid) go to 10
80     continue
       iemx = imax
       ione = 1
c
c      Grid in 1-direction is complete.  Round values for "x1a" grid to
c  "isig" significant figures.
c
       call round ( x1a, in, isig )
c
c      Grid consistency check.  Make sure origin always occurs at
c  x1a(ismn).
c
       if ( (x1a(ismn-1) .le. 0.0) .and. (x1a(ismn) .ne. 0.0) ) then
         if (iolog .gt. 0) write (iolog, 2110) x1a(ismn)
         if (iotty .gt. 0) write (iotty, 2110) x1a(ismn)
         stop
       endif
c
c-----------------------------------------------------------------------
c----------------------- Write format statements -----------------------
c-----------------------------------------------------------------------
c
2010   format('GRIDX1  : ***** ERROR ***** Number of zones in '
     1       ,' 1-direction exceeds array bounds',/
     2       ,'GRIDX1  : imax = ',i4,'  in = ',i4,/
     3       ,'GRIDX1  : ABORTS!')
2020   format('GRIDX1  : ***** ERROR ***** Newton-Raphson did not'
     1       ,' converge for x1rat.',/
     2       ,'GRIDX1  : imin = ',i3,' x1r = ',1pe12.5,' deltx1r = '
     3       ,1e12.5,' fn = ',1e12.5,/
     4       ,'GRIDX1  : ABORTS!')
2030   format('GRIDX1  : ***** ERROR ***** Hey chucklehead, you '
     1       ,'can''t have x1min = ',1pg12.5,/
     2       ,'GRIDX1  : with the scaled zone option! (igrid = 3)',/
     3       ,'GRIDX1  : ABORTS!')
2040   format('GRIDX1  : ***** ERROR ***** Illegal value of igrid '
     1       ,'during x1 grid generation.',/
     2       ,'GRIDX1  : imin = ',i3,'  igrid = ',i2,/
     3       ,'GRIDX1  : ABORTS!')
2050   format('GRIDX1  :',/
     1       ,'GRIDX1  : --------------------> x1a (face-centred) '
     2       ,'grid, zones ',i4,' through ',i4)
2060   format('GRIDX1  : These zones read from restart dump ',a8)
2070   format('GRIDX1  : grid parameters:  dx1min  = ',1pg12.5
     1       ,', x1rat = ',g12.5)
2080   format('GRIDX1  : grid parameters:  dx1max  = ',1pg12.5
     1       ,', x1rat = ',g12.5)
2090   format('GRIDX1  :                   x1scale = ',1pg12.5)
2100   format('GRIDX1  : ',1p5g14.6)
2110   format('GRIDX1  : ***** ERROR ***** Ill-defined grid, zero of '
     1       ,'x1 not at x1a(ismn).',/
     2       ,'GRIDX1  : x1a(ismn) = ',1pe12.5,/
     3       ,'GRIDX1  : ABORTS!')
c
       return
       end
