      program zeus3d
c
c    dac:zeus3d <------------------------- multi-physics, 3-D fluid code
c    from jms:zeus2d (NCSA),                               january, 1990
c         mln,dac:zeus04 (NCSA),
c         mln,dac:b2 (LANL),
c         mln:a2 (LANL)
c
c    Written by: David Clarke
c    Modified 1: Robert Fiedler, 2/95, so that the calls to the various
c                timing routines work properly on Convex, Sun, SGI,
c                IBM, and HP machines even if compiler flags are used
c                to promote the default precision to 64-bits.
c
c    Modified by: Gang Li 08/2006 
c                 Delete ishock input since it is in gen_bndy.f
c
c  PURPOSE:  The main program orchestrates the calls to the major
c  subroutines in the hydrodynamical loop, and the subroutines which
c  handle the output.
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-boundares
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     
c-----------------------------------------------------------------------
c------------------ Start of Main program ------------------------------
c-----------------------------------------------------------------------
c
       integer       nzones  , nzceff  , nzcact  , cpuhours, cpumins
     1             , cpusecs , wchours , wcmins  , wcsecs
       real*8        zcspk   , zcseff  , zcsact  , cputot  , wctot
     1             , cpu1    , cpu2    , ptime
       real*4        tarray  (   2)
       integer*4     iarray  (   3)
c
       integer i,j,maxshells
       real*8    max_p
c
       parameter (maxshells = 200, max_p = 10.0D0)
c
       real*8   R_shell_x(maxshells)
       integer  R_shell_i(maxshells)
       real*8   esc_shell(in,250)
       real*8   shock_speed, Machno, comp_rat
       real*8   prev_shock_posn, prev_shock_time
       real*8   p_cut_bohm, p_cut_gorlee,p_cut_gorlee_hi
       real*8   p_inj, p_cut, distr_array(maxshells,250)
       real*8   esc_dens_up(250),esc_dens_down(250)
       real*8   prev_acc_time, prev_expand_time(maxshells), 
     *       prev_acc_posn
       real*8   arrive_distr(250,10) 
       integer  no_shells, shell_dt_count
       integer  shock_begin, shock_end     
       real*8   charge, Amass, p_hi_inj, p_hi_cut
       real*8   distr_array_hi(maxshells,250),esc_dens_up_hi(250),
     *      esc_dens_down_hi(250),esc_shell_hi(in,250)
c
       real *8 to, vo
       real *8 p_sml, energy_sml, p_big, energy_big, mp, c

       integer  N_shell_dt      ! controls the time steps for 
                                ! generating a new shell.

       integer i_shock_strength ! controls the shock strength.
       integer i_heavy          ! controls CNO or IRON.

       integer i_perp           ! controls if this is a perp. shock.
       integer i_restart        ! controls if this is a restart.

       integer R_shell_i_prev(maxshells)! used to calculate the shell position

       real*8 esc_ptcl_counter(250) ! time-integrated escaped particles.
       real*8 esc_ptcl_counter_hi(250) ! time-integrated escaped particles.


       real*8 ipmax(maxshells),ipmin(maxshells)
                                ! ipmax(j) is the maximum proton momentum index
                                ! for the jth shell.
       real*8 ipmax_hi(maxshells),ipmin_hi(maxshells)
                                ! ipmax(j) is the maximum HvyIon momentum index
                                ! for the jth shell.

       
       integer j_p_cut_chk, j_p_cut
       integer j_p_cut_chk_hi, j_p_cut_hi
       
       real*8   erfcc
       real*8   distr_array_all(250)
       real*8   distr_all(250)

       real*8  shock_start_time

       real*8  alphaK, theta_BN
       
       common /wave_add/j_p_cut_chk,j_p_cut,j_p_cut_chk_hi,j_p_cut_hi

       common /shell/  esc_shell, R_shell_x, R_shell_i, no_shells 
       common /shock1/ prev_shock_posn, prev_shock_time
       common /shock2/ shock_begin, shock_end, shock_speed,
     *                 Machno, comp_rat  
       common /ptcls/  p_inj, p_cut, distr_array, esc_dens_up,
     *                 esc_dens_down, prev_acc_time,
     *      prev_expand_time, prev_acc_posn
       common /shocktime/ shock_start_time

       common /arrive/ arrive_distr  
       common /heavyIon/ charge, Amass, p_hi_inj, p_hi_cut
       common /heavyIon2/distr_array_hi,esc_dens_up_hi,esc_dens_down_hi,
     *      esc_shell_hi

       common /kappa_perp/ i_perp
c       common /gl_save/ shell_dt_count

       common /R_shell_int/R_shell_i_prev

       common /p_checker/ipmax,ipmin,ipmax_hi,ipmin_hi
       common /escape_counter/ esc_ptcl_counter, esc_ptcl_counter_hi
       common /total_ptcl /distr_array_all       
       common /diag/distr_all
       common /zeus3d_wave_intens/ alphaK, theta_BN

       integer ctoi
       external ctoi 

       real*8 res_time
       common /dump_time/res_time

       common /shell_control/ shell_dt_count, N_shell_dt  

c     
c      External statements.
c
       real*4        etime
       external      etime   , itime
       external      empty   , renam
c
c      The following external statements are defined only if the EDITOR
c  alias is aliased to something other than "empty".  This prevents
c  "empty" from being declared in more than one external statement.
c
       external      gen_bndy  
       external      mstart
       external      srcstep
       external      transprt
       external      nudt
       external      dataio

       integer isolar
       real *8 time_0
       real *8 tmtreT


c
c-----------------------------------------------------------------------
c
c      Initialise cpu and wall clocks.  "cputot" will be the total cpu
c  time (in seconds) consumed by this job.  "wctot" will be the total
c  elapsed wall clock time (in seconds) since the job began.
c

       cputot = etime ( tarray )
       call itime ( iarray )
       wctot  = iarray(1) * 3600. + iarray(2) * 60. + iarray(3)
c
c      Initialise execution and make dump of initial model.
c
c --- Jan. 11th, 2006 GL: important parameters in strange unit. 
       to = 2858055.0   ! time unit  (in seconds)
       vo = 52483.25D0  ! speed unit (in m/sec) 
       mp = 1.67D-27    ! mass of proton (in gram)
       c = 3.0D8        ! speed of light (in m/sec)



       write(*,*) "input iheavy (1 for CNO, 2 for IRON)"
       write(*,*) "Is this a CONTINUED run (1:yes,0:no)"
       read(*,*) i_heavy
       read(*,*) i_restart

       write(*,*) "input i_perp (0: Bohm, 1:perp) "
       read(*,*) i_perp         ! i_perp = 0

       write(*,*) "input alphaK (0.001)"
       read(*,*) alphaK !        alphaK = 0.001

       write(*,*) "input theta_BN (80.)"
       read(*,*) theta_BN !      theta_BN = 80 degree

       call mstart

       res_time = 0.
       iresdmp = ctoi ( resfile, 3, 4 )
       print*, "iresdmp is", iresdmp
       print*, "time is", time
       res_time = time
c       pause

       if (iresdmp .eq. 0) then ! usrdmp file zr01IH 
          res_time = 0.
          isolar = 0
          print*, "---going to 2000-----"
          go to 2000
       else 
          isolar = 1
          print*,"res_time is",res_time 
          open (27,file="start_time.dat")
          write(27,*) res_time, x1a(is)
          close(27)
       endif

c---- initial values for our own code. --GL
       j_p_cut = 0                
       j_p_cut_hi = 0 
       j_p_cut_chk = 0
       j_p_cut_chk_hi = 0

       N_shell_dt = 150 ! the initial N_shell_dt

c     
c ---  Initialise ptcl acc variables and arrays.        
c
       Do i = 1, maxshells
         R_shell_x(i) = 0.1
         R_shell_i(i) = is + 3    
         Do j = 1, 250 
	    esc_shell(i,j) = 0.0
         End Do
       End Do 

       Do i = 1, maxshells
         R_shell_i_prev(i) = 0    
       End Do 
c
       Do i = 1, 10
         Do j = 1, 250
	   arrive_distr(j,i) = 0.0D0
	 EndDo
       EndDo	  

c ---initialize distr_array
       Do i =1, maxshells
          Do j =1, 250
             distr_array(i,j) = 0.  ! for proton
             distr_array_hi(i,j) = 0. ! heavy ion
          ENDDO
       ENDDO

c ---initialize distr_array_all
          Do j =1, 250
             distr_array_all(j) = 0.
          ENDDO

c ---initialize esc_ptcl_counter (p_number)
       Do j =1, 250
          esc_ptcl_counter(j) = 0.
          esc_ptcl_counter_hi(j) = 0.
          distr_all(j) = 0.
       ENDDO

c
c       shock_speed = 10.9       ! Initial shock speed.       
       prev_shock_posn = x1a(is) ! Shock starting position.
       prev_shock_time = time   ! Shock starting time.
       prev_acc_time = time     ! Acceleration starting time.
       prev_acc_posn = x1a(is)  ! Acceleration starting position.
       R_shell_x(1) = x1a(is)
       R_shell_i(1) = is
       R_shell_x(2) = x1a(is)
       R_shell_i(2) = is

       shock_begin = 0

       shock_start_time = time
c
c --------------------------------------------------------------------  
c  Open output files
c --------------------------------------------------------------------
c

       open(13,file='ptcls_in_shell.dat',form='formatted',
     *         status='unknown')
c
       open(14,file='distr_in_shell.dat',form='formatted',
     *         status='unknown')
c
       open(35,file='obs_ptcls.dat',form='formatted',
     *         status='unknown')
c
       open(36,file='posn_distr.dat',form='formatted',
     *         status='unknown')
c
       open(37,file='time_distr.dat',form='formatted',
     *         status='unknown')
c
       open(38,file='esc_distr.dat',form='formatted',
     *         status='unknown')                   
c
       open(39,file='distr_1AU_time.dat',form='formatted',
     *         status='unknown')   
c
       open(43,file='waveintens.dat',form='formatted',
     *         status='unknown')
c    
       open(44,file='kappa.dat',form='formatted',
     *         status='unknown') 
c
       open(45,file='shock_momenta.dat',form='formatted',
     *         status='unknown') 
c

c       open(46,file='f_acc.dat',form='formatted', 
c     *         status='unknown')   

c
       open(47,file='distr_fn_out.dat',form='formatted',
     *         status='unknown')  
c    
       open(48,file='shock_posn_comp.dat',form='formatted',
     *      status='unknown') 

       open(55,file='momenta-hi.dat',form='formatted',
     *         status='unknown') 

       open(58,file='esc_distr-hi.dat',form='formatted',
     *         status='unknown')                   

       open(65,file='momenta-bohm.dat',form='formatted',
     *         status='unknown') 

       open(66,file='f_acc-hi.dat',form='formatted', 
     *         status='unknown')   

       open(67,file='kappa-wave-I.dat',form='formatted',
     *         status='unknown') 

       open(68,file='shell-location.dat',form='formatted',
     *         status='unknown') 

       open(69,file='shell-location2.dat',form='formatted',
     *         status='unknown') 

       open(78,file='current_at_shock.dat',form='formatted', 
     *         status='unknown')   

       open(79,file='dist_at_shock.dat',form='formatted', 
     *         status='unknown')   

       open(88,file='esc_num_up.dat',form='formatted', 
     *         status='unknown')   

       open(89,file='esc_num_up_hi.dat',form='formatted', 
     *         status='unknown')   
       
       open(90, file ='zeus.log', status='unknown')   

       open(91, file ="gl.input", status='unknown')  ! no use now. -- GL

       open(92, file ="gl.output", status='unknown') ! no use now. -- GL

       no_shells = 1
       shell_dt_count = 0

c   --- ZEUS part.

 2000  continue
       if (iolog .gt. 0) then
         if (nwarn .eq. 0) write (iolog, 2010)
         if (nwarn .eq. 1) write (iolog, 2020)
         if (nwarn .gt. 1) write (iolog, 2030) nwarn
       endif
       if (iotty .gt. 0) then
         if (nwarn .eq. 0) write (iotty, 2010)
         if (nwarn .eq. 1) write (iotty, 2020)
         if (nwarn .gt. 1) write (iotty, 2030) nwarn
       endif
       nzceff  = 0
       nzcact  = 0
       zcspk   = 0.0
       tused   = 1.0
       cpu1    = etime ( tarray )

       time_0 = dtdmp
       print*, "--- dtdmp ---", dtdmp
c       pause

       if (isolar .eq. 0) then
          print*, "----going to 10 ----"
          go to 10
       endif

c ---initialize distr_array, esc_num_d

       if (i_restart .eq. 1) then
          print*, "---- Read continution from gl.input -----"
          
          read (91,*) no_shells,shell_dt_count,shock_speed
          read (91,*) prev_shock_posn,prev_shock_time
          read (91,*) prev_acc_time, prev_acc_posn

          Do i = 1, maxshells
             read(91,*) R_shell_x(i), R_shell_i(i), 
     *            prev_expand_time(i)
          End Do 
          
          Do i =1, maxshells
             Do j =1, 250
                read(91,*) distr_array(i,j),distr_array_hi(i,j)
             ENDDO
          ENDDO

       print*, "----successful read in gl.input -----"
       endif
       
c---------------
       
       if (ifsen .eq. 1) go to 20
c     
c-----------------------------------------------------------------------
c--------------------------  START OF MAIN LOOP  -----------------------
c-----------------------------------------------------------------------
c
c      Execution ends when INTCHK returns a value of 1 for "ifsen".
c
       print*, " ------ In zeus3d.f -------"

c---- only the ratio of amass/charge matters.

       if (i_heavy .eq. 0) then ! For checking code.
          amass = 1. 
          charge = 1.
       elseif (i_heavy .eq. 1) then ! select CNO.
c---  for CNO
          amass = 14. 
          charge = 6.57
       elseif ( i_heavy .eq. 2) then ! select IRON.
c---  for IRON
          amass = 56.
          charge = 14.
       endif

       write(90,*) " --- This file is generated at the begining "
       write(90,*) " --- of running xzeus -----"
       write(90,*) " i_heavy is: ", i_heavy
       write(90,*) " amass and charge are ", amass, charge 
       close (90)

10     continue
c----------- [FLAG]: most of the following are ZEUS part.

       call gen_bndy  ! boundary condition subroutine need to be changed for different shocks.
       call empty
       call empty
       call srcstep
       call empty
       call transprt
       call empty
       call nudt
       call empty
       call dataio

       nhy    = nhy   + 1
       time   = time  + dt
c
c       print*, 'nhy=',nhy,' & dt=',dt,' & time=',time
c     
       if (isolar .eq. 0) then
          if ( time .gt. time_0) then
             do i = 1, 450      ! boundary at 5000  write solar wind speed to data file.
c                print*, "----gamma is ---", gamma(1)
c                pause

                tmtreT = e(i*10, 1, 5)*(gamma(1)-1.)/(6.*d(i*10, 1, 5))
                write(22,*)x1a(i*10), d(i*10, 1, 5),
     X               v1(i*10,1,5)*vo/1000., tmtreT, time 
                write(23,*)x1a(i*10), b1(i*10, 1, 5),
     X               b3(i*10,1,5), time
             enddo
             time_0 = time_0 + dtdmp
             print*, "----writing solar wind speed to a dump file --"
             print*, "iresdmp is", iresdmp
          endif
c          print*, "---going to 2001-----"
          go to 2001 
       endif

c-----------END [FLAG].

c---   The following deals with shells. GL

       shell_dt_count = shell_dt_count + 1
c       
c       print*, "shell_dt_count, N_shell_dt",shell_dt_count, N_shell_dt
 
       If ((shell_dt_count.gt.0) .and. 
     *      (shell_dt_count.lt. N_shell_dt)) Then
          call shell_track
	  If (mod(shell_dt_count,5*N_shell_dt) .eq. 0) then
            call posn_distr(0.8D0)
	  EndIf  	    
       EndIf
c     
c
       If ((shell_dt_count .eq. N_shell_dt).or. 
     *      (shell_dt_count.eq.0)) Then
          call shell_track	 
          
          print*, "prev_acc_time ",prev_acc_time 
          print*, "prev_acc_posn", prev_acc_posn
          print*, 'time = ', time
          print*, "---time,shock_begin, x1a(shock_begin)---", 
     *         time, shock_begin, x1a(shock_begin)
          if( (prev_acc_posn .gt. x1a(shock_begin)) .and.
     *         (shell_dt_count.ne.0) )then
             print*, "N_shell_dt", N_shell_dt
             print*, "prev_acc_posn, x1a(shock_begin)", 
     *            prev_acc_posn, x1a(shock_begin)
c             pause
          endif

          If (shell_dt_count .eq. N_shell_dt) Then ! create a new shell.
             
             call wave_intens   ! get the wave intensity
             call expand_distr  ! Calculate the simultaneous expansion with the solar wind.
             call add_acc_ptcls ! At the shock front, add the newly generated particles.
             call diffuse_distr ! Calculate the change of distribution due to diffusion.
c------------
             
             print*, " --- After subroutine calls ----" 
             
c---- Check how far are we from the sun and update N_shell_dt if necessary.
             
             if (x1a(shock_begin) .lt. 0.35 ) then
                N_shell_dt = 120 ! time steps for generating a new shell.
             else
                N_shell_dt = 150 ! time steps for generating a new shell.
             endif
             
             prev_acc_time = time
             prev_acc_posn = x1a(shock_begin)
             
             print*, "shock_speed,shock_posn",shock_speed, prev_acc_posn
             print*, "no_shells", no_shells
c     
c---- write location of the shock, maximum and minimum proton energy
c---- to file, shock_momenta.dat(45)

             p_sml = (mp*vo*c*p_inj)/1.6D-19 ! Inject  mom. in eV
             p_big = (mp*vo*c*p_cut)/1.6D-19 ! Maximum mom. in eV
             
             energy_sml = (Sqrt((mp*c**2.0)**2.0+(mp*vo*c*p_inj)**2.0)- !Inject kinetic 
     *            mp*c**2.0)/1.6D-19 !energy in eV
             
             energy_big = (Sqrt((mp*c**2.0)**2.0+(mp*vo*c*p_cut)**2.0)- !kinetic energy
     *            mp*c**2.0)/1.6D-19 !energy in eV
             
             write(45,2110)x1a(shock_begin),(time-res_time)*to/3600.,
     *            p_sml,p_big,energy_sml, energy_big
             
             print*, " *** ---------***"
             
             print*, x1a(shock_begin),(time-res_time)*to/3600.,
     *            p_sml,p_big,energy_sml, energy_big
             
             print*, " *** ---------***"
             
c---- write location of the shock, maximum and minimum proton energy
c---- to file, shock_momenta.dat(45)

            write(65,2110)x1a(shock_begin),(time-1.00003167)*to/3600.,
     *           p_sml,p_cut_bohm

c---- write location of the shock, maximum and minimum heavy ion energy
c---- to file, shock_momenta.dat(55)

             p_sml = (mp*vo*c*p_hi_inj)/1.6D-19 ! Inject  mom. in eV
             p_big = (mp*vo*c*p_hi_cut)/1.6D-19 ! Maximum mom. in eV
           
             energy_sml = (Sqrt((mp*c**2.0)**2+(mp*vo*c*p_hi_inj)**2.0)- !Inject kinetic 
     *            mp*c**2.0)/1.6D-19 !energy in eV
             
             energy_big = (Sqrt((mp*c**2.0)**2+(mp*vo*c*p_hi_cut)**2.0)- !kinetic energy
     *            mp*c**2.0)/1.6D-19 !energy in eV
             
             write(55,2110)x1a(shock_begin),(time-res_time)*to/3600.,
     *            p_sml,p_big,energy_sml, energy_big
            
c---- write shock location, comp_rat, speed, time
c---- to file, shock_posn_comp.dat (48)

             print*,"output in shock_posn*.dat",time,res_time 
             write(48,2100)x1a(shock_begin),comp_rat,
     *            shock_speed*vo/1000.,v1(shock_end+2,1,5)*vo/1000.,	    
     *            (time - res_time)*to/3600.
             
             write(*,*) "x1a(shock_begin),shock_speed,comp_rat",
     *           x1a(shock_begin),shock_speed*vo/1000.,comp_rat
             
             no_shells = no_shells + 1
             shell_dt_count = 0  
c--- ( debug
c            pause
c--- ) debug

          EndIf
       EndIf

c---- ZEUS part of the code.       	 	              
 2001 continue

      nzones = nx1z * nx2z * nx3z
      cpu2   = cpu1
      cpu1   = etime ( tarray )
      cpu1   = tarray(1)
      cpu2   = cpu1 - cpu2
      zcspk  =   max ( zcspk,  real(nzones) / cpu2 )
      nzceff = nzceff + nx1zmx * nx2zmx * nx3zmx
      nzcact = nzcact + nzones
      if (ifsen .eq. 1) go to 20
      go to 10
 20   continue
c     
      close(47)
      close(46)
      close(45)
      close(48)
      close(44)
      close(43)
      close(39)
      close(38)
      close(37)
      close(36)
      close(35)
      close(54)
      close(55)
      
c     
c-----------------------------------------------------------------------
c---------------------------  END OF MAIN LOOP  ------------------------
c---------------------------  ZEUS part of the code  ------------------------
c-----------------------------------------------------------------------
c
c      Perform any user-specified clean-up tasks.
c
       call empty
c
c      Determine number of zone-cycles per second.
c
       zcseff   =  real(nzceff) / tused
       zcsact   =  real(nzcact) / tused
c
c      Final wall clock time, expressed in hours, minutes, and seconds.
c
       call itime (iarray)
       wctot  = iarray(1) * 3600. + iarray(2) * 60. + iarray(3) - wctot
       wchours  =  int ( wctot / 3600.0 )
       wcmins   =  int ( wctot / 60.0   ) - 60   * wchours
       wcsecs   =  int ( wctot + 0.5    ) - 3600 * wchours
     1                                    - 60   * wcmins
c
c      cpu usage, expressed in hours, minutes, and seconds.
c
       cpusecs  = etime ( tarray )
       cputot   = tarray(1) - cputot
       cpuhours =  int ( cputot / 3600.0 )
       cpumins  =  int ( cputot / 60.0   ) - 60   * cpuhours
       cpusecs  =  int ( cputot + 0.5    ) - 3600 * cpuhours
     1                                     - 60   * cpumins
c
c      Percent of real time that job was serviced.
c
       if (wctot .le. 0.0) wctot =   max ( cputot, tiny )
       ptime    = 100.0 * cputot / wctot
c
c      Final messages to crt and log.
c
       if (iolog .gt. 0) then
         if (nwarn .eq. 0) write (iolog, 2040)
         if (nwarn .eq. 1) write (iolog, 2050)
         if (nwarn .gt. 1) write (iolog, 2060) nwarn
         write (iolog, 2070) zcspk, zcseff, zcsact
c         write (iolog, 2080) wchours, wcmins, wcsecs, cpuhours
c     1                     , cpumins, cpusecs, ptime
         write (iolog, 2080) wctot, cputot, ptime
       endif
       if (iotty .gt. 0) then
         if (nwarn .eq. 0) write (iotty, 2040)
         if (nwarn .eq. 1) write (iotty, 2050)
         if (nwarn .gt. 1) write (iotty, 2060) nwarn
         write (iotty, 2070) zcspk, zcseff, zcsact
c         write (iotty, 2080) wchours, wcmins, wcsecs, cpuhours
c     1                     , cpumins, cpusecs, ptime
         write (iotty, 2080) wctot, cputot, ptime
       endif
c
c      Close ascii files, and then AMF!
c
       if (iolog .gt. 0) then
         close ( unit=iolog )
         ilogdmp = ilogdmp + 1
         write (logfile, 2090) 'zl', ilogdmp, id
         call renam ( 'zltemp', logfile(1:6) )
       endif
       if (dttsl .gt. 0.0) close ( unit=iotsl )
c
c-----------------------------------------------------------------------
c----------------------- Write format statements -----------------------
c-----------------------------------------------------------------------
c
2010   format('ZEUS3D  :',/
     1       ,'ZEUS3D  : Setup complete with no warnings issued, '
     2       ,'entering main loop...',/,'ZEUS3D  :')
2020   format('ZEUS3D  :',/
     1       ,'ZEUS3D  : Setup complete with 1 warning issued, '
     2       ,'entering main loop...',/,'ZEUS3D  :')
2030   format('ZEUS3D  :',/
     1       ,'ZEUS3D  : Setup complete with ',i2,' warnings issued, '
     2       ,'entering main loop...',/,'ZEUS3D  :')
2040   format('ZEUS3D  :',/,'ZEUS3D  : Execution terminating with no '
     1       ,'warnings issued.')
2050   format('ZEUS3D  :',/,'ZEUS3D  : Execution terminating with 1 '
     1       ,'warning issued.')
2060   format('ZEUS3D  :',/,'ZEUS3D  : Execution terminating with ',i4
     1       ,' warnings issued.')
2070   format('ZEUS3D  : Peak      zone-cycles per second for this '
     1       ,'job = ',1pe12.5,/
     2       ,'ZEUS3D  : Effective zone-cycles per second for this '
     3       ,'job = ',1pe12.5,/
     4       ,'ZEUS3D  : Average   zone-cycles per second for this '
     5       ,'job = ',1pe12.5)
c2080   format('ZEUS3D  : Real time =',i3,':',i2.2,':',i2.2
c     1       ,', CPU time =',i3,':',i2.2,':',i2.2,' (',f6.2,'%)')
2080   format('ZEUS3D  : Wall clock time = ', f12.2
     1       ,', CPU time = ', f12.2, ' (',f6.2,'%)')
2090   format(a2,i2.2,a2)
c
       stop

 2100  format(2x, f10.5, 2x, f8.4, 2(2x,f10.5), 2x, f12.5)
 2110      format( 2(f10.4,2x), 4(e12.6,2x) )


       print*, " ****** Leave zeus3d.f *******"
       end

      Function erfcc(x)
      Real *8 erfcc, x
      real *8 t, z
      z = abs(x)
      t = 1. /(1. + 0.5 * z)
      erfcc = t * Dexp(-z*z - 1.26551223 + t * ( 1.00002368 + 
     *     t* (0.39409196 + t * (0.09678418 + t* (-0.18628806 + 
     *     t* (-1.13520398 + t * (1.48851587 + t * (-0.82215223 +
     *     t* 0.17087277))))))))
      if ( x .lt. 0.) erfcc = 2. - erfcc
      return
      end
