       subroutine wave_intens
c
c    dac:zeus3d.wave_intens <---------------- Calculates the wave intensity
c          (see Gordon et al. 1999) upstream of a shock for a given momentum.
c
c    Written by : Ken Rice
c
c                           Modified : 04/19/2000
c
c    Modified by:  Gang Li (Extensive changes)
c  EXTERNALS: [NONE]
c
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
c-----------------------------------------------------------------------
c-------------- START OF PROGRAM/SUBROUTINE ----------------------------
c
c  This subroutine is all written in cgs units and is taken from 
c  Gordon et al., 1999.
c
c-----------------------------------------------------------------------
c
       integer maxshells
       real*8  max_p, p0
c
       parameter (maxshells = 200, max_p = 10.0D0, p0 = 2.0)
c
       real*8 p, psf, del_p, vel, gamma_rel, energy
       real*8 B, muo
       real*8 vo, VA, Vsw, Vsh, V, Vpr, Vd
       real*8 n, mp, rg
       real*8 Lo, k, ko        
       real*8 injectrate, intensnorm
       real*8 alpha_wave, Io
       real*8 pinj, acc_time, acc_time_hi
       real*8 c, ech, gyrof
       real*8 psi, beta
       real*8 shock_speed, Machno, comp_rat
       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
c
       real*8 Iplus(250,maxshells), Ipluso(250), Ibohm(250)
       real*8 fo(250,maxshells)
       real*8 kappa(250,maxshells),kappa_hi(250,maxshells)
       real*8 kappa_perp(250,maxshells), kappa_perp_hi(250,maxshells)
       real*8 s_ratio, s(250)
       real*8 wave_comp_rat
       real*8 x, dx, xintegral(250,1000), uintegral, u, Uo
c
       integer p_int
       integer shock_begin, shock_end
       integer i, j, l, index, ik
       real *8 Amass, charge,p_hi_inj, p_hi_cut, p_hi_cut2
       real *8 p_max, p_max_hi,time_all, psf_hi
       real *8 energy_inj
       real *8 Iplustemp, Iplusotemp

       real*8   esc_shell(in, 250), R_shell_x(maxshells)
       integer  R_shell_i(maxshells), no_shells

       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.
       real*8 temp_p, temp_k, temp_I ! for heavy ions.
       real*8 temp, gammag
       
       real*8 alpha_bohm
       real*8 alphaK

       real*8 shock_start_time

       integer j_p_cut_chk, j_p_cut
       integer j_p_cut_chk_hi, j_p_cut_hi

       real*8 to
       real*8 time_initslap
       integer i_perp

c----- the following is added for the function get_injection_speed(....)       
       real*8 get_injection_speed
       real*8 theta_BN, deltaB,B0,uspeed,comp,alpha_inj,s2d, 
     X      lambda_c, vinj, tempVd
       real*8 lambda_perp, lambda_par
c----- the following is added for the function get_rl(....)       
       real*8 get_rl
       real*8 tempP, temp_rl, temp_kappa_perp, temp_kappa_par, tempt

       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 /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 /diff_coeff/ kappa,kappa_hi
       common /heavyIon/ charge, Amass, p_hi_inj, p_hi_cut
       common /p_checker/ipmax,ipmin,ipmax_hi,ipmin_hi

       common /wave_action/ Iplus
       common /time_init/ time_initslap
       common /kappa_perp/ i_perp
       common /zeus3d_wave_intens/ alphaK, theta_BN
       real *8 shock_init_loc
c
c---------------------------------------------------------------------- 
c
       print*, "-------Entering wave_intens.f -------"
       print*, " -- time is --", time 
       print*, "*** no_shells is  ***", no_shells
c
       print*, "amass, charge", amass, charge
       muo = 4.0D0*pi*1.0D-7       
c        
       del_p = (max_p - p0)/250.0D0
c
       vo = 52483.25   ! Velocity normalization in MKS units.
       to = 2858055.0  ! time normalization in seconds.
c                        
       alpha_bohm = 2.

       gammag = 5./3.


       shock_init_loc = 0.025 ! the "real" shock initialization point.

c       alphaK = 0.001

c----- set j_p_cut_chk = j_p_cut
       j_p_cut_chk = j_p_cut
       j_p_cut_chk_hi = j_p_cut_hi

       call parker_field(B, x1a(shock_end+10), v1(shock_end+10,1,1))
c       print*, "B_parker", B
       B = Sqrt(B1(shock_end+2,1,5)**2.0+B3(shock_end+2,1,5)**2.0)
c       print*, "B", B
c
       B = B*5.38E-9                ! Unnormalized B in Tesla.
c       B0 = 5.38D-9 is the scale of magnetic field through B^2/2 ~ pressure.

       n = d(shock_end+5,1,5)*5.0D6   ! Density in m^-3.
c       n0 = 5.0D6 is the normalized number density at 0.1 AU.
       mp = 1.67D-27                ! Proton mass in kg.                   
c                        
       VA = Sqrt(B*B/(muo*mp*n)) ! Alfven speed in m/s.
c
       Lo = 0.025               ! Correlation length in AU.
       ko = 2.0*pi/(Lo*1.5D11)  ! ko in m^-1.
c     
       c = 3.0D8                ! Speed of light in m/s.
       ech = 1.6D-19            ! Electron charge in coulombs
c     
       gyrof = ech*B/mp         ! Non-relativistic Gyrofrequency in MKS.
c       print*, gyrof,'    gyrof'
c
       psi = 0.0*pi/180.0       ! Shock obliquity
c     
       Vsw = V1(shock_end+10,1,1)*vo ! Solar wind speed in m s^-1.
c     Vsw = 400000.0D0
       
       Vsh = shock_speed*vo     ! Shock speed in MKS.
c     Vsh = 540000.0D0 

       V = Vsh - Vsw            ! upstream speed in shock frame.
       If (VA*cos(psi) .LT. V) Then
          Vpr = V - VA*cos(psi) ! Scattering center speed.
       Else
          Vpr = 0.0
       End If  
       Vpr = V     
c     
       Vd = (shock_speed-v1(shock_begin-5,1,1))*vo
c     
c       print*, "Vsw, Vsh, V, VA, Vpr, Vd, x1a(shock_begin)"
c       print*, Vsw, Vsh, V, VA, Vpr, Vd, x1a(shock_begin)
c   
       wave_comp_rat = Vpr/Vd
c     
       beta = 3.0*wave_comp_rat/(wave_comp_rat - 1.0) 
c       print*, "beta, comp_rat, wave_comp_rat"

       beta = 3.0*comp_rat/(comp_rat - 1.0)
       print*, "beta, comp_rat, wave_comp_rat",
     *      beta, comp_rat, wave_comp_rat

       print*, " shock_begin , shock_end",shock_begin,shock_end
c     
       injectrate = 0.01D0*n*V  !simple approach on injection.
c                        
       alpha_wave = 1.0D0                 ! Fraction of magnetic energy density in waves.
       Io = 8. * pi * alpha_wave/(5.0*ko)*B**2.0/(2.0*muo) ! Background intensity.

       temp = 2.0/((gammag - 1.0)*Machno*Machno)  
       temp = 1.0 - 1.0/(comp_rat*comp_rat) + temp

       p_inj = 3.0*0.5*(gammag - 1.0)
       p_inj = p_inj*(shock_speed - v1(shock_end+2,1,5))**3
       p_inj = p_inj/(2.0*gammag
     *           *(shock_speed - v1(shock_end+2, 1, 5)))
c
       p_inj = Sqrt(p_inj*temp)    !  Lower boundary of injection momentum.
c
c       p_inj = 3. * (shock_speed-v1(shock_end+2,1,5))/(comp_rat -1.)
c------ for parallel case, take zeta = 3u/ (v*(s-1))  = 0.5 to decide u
       
       print*, "*** p_inj, temp, comp_rat ***", p_inj, temp, comp_rat
       print*, " ---- revised p_inj for para case ----", p_inj, comp_rat

       print*, "--- shock_speed,v1 --",shock_speed, v1(shock_end+2,1,5)
       print*, "--- v0 ---", vo
c      pause  ! GL

       p_inj = shock_speed + p_inj   

       if (i_perp .eq. 1) then  ! redecide p_inj
          
c ---------use function get_injection_speed(theta_BN,deltaB,B0,uspeed,comp,alpha, beta)
c     -------- theta_BN == the angle between shock normal and B0 in degree.
c     -------- deltaB == the turbulent field delta_B.
c     -------- B0 == ambient B.
c     -------- uspeed == flow speed.
c     -------- alpha == lambda_2D / lambda_slab 
c     -------- s2d == (b_2D/b)^2
c ---------------------------------------------------------------------------

c       theta_BN = 80.  ! the angle in degree.

          B0 =  B * 1.D4        ! in Gauss
          print*, "x1a(shock_end+10)", x1a(shock_end+10)
          deltaB = 13.6**(0.5) * 1.D-5 * (x1a(shock_end+10)/1.)**(-3/2.) ! in gauss
          uspeed = (shock_speed - v1(shock_end+10, 1,1)) *vo/1000. ! upstream flow speed in shock frame in km/s
          tempVd = (shock_speed-v1(shock_begin-5,1,1))*vo/1000. ! dnstream flow speed in shock frame in km/s.
          comp = uspeed/tempVd  ! compression ratio.
          alpha_inj = 0.1       ! lambda_2D / lambda_slab 
          s2d = 0.8
          lambda_c = Lo * 1.5 D8 ! correlation length in km (assuming not change with distance).
          
          vinj = get_injection_speed(theta_BN,deltaB,B0,uspeed,
     X         comp,alpha_inj, s2d, lambda_c, lambda_perp, lambda_par) ! in km/s
          
          p_inj = vinj * 1000. / vo ! normalized p_inj (vo is in m/s)
          
          print*, "********** p_inj for perp case **********", p_inj
       endif

c     pause
       
       pinj = p_inj*mp*vo       ! un-normalized pinj
c     
       intensnorm = 8.0*gyrof*Vpr**2.0*mp/pinj
       intensnorm = intensnorm/(VA*pinj*beta*injectrate*cos(psi))
       intensnorm = intensnorm*mp/pinj
c
       If (intensnorm .LE. 0.0) Then
          intensnorm = 1.0E30   ! Make intensnorm very large - no wave growth.
       End If  
c     
c       write(43,*) "time is", time
c       write(44,*) "time is", time, x1a(shock_begin)

c       print*, time
       
c---- decide kappa perp.
       Do i = 1, 250
          p = p0 + del_p * (i-1)
          psf = Exp(p)*mp*vo    ! momentum in MKS.
          vel = psf/mp*(1/Sqrt(1.0 + (psf/(mp*c))**2.0)) 
          
          kappa_perp(i,no_shells) = (R_shell_x(no_shells)/1.0)**(2/3.) 
     *         * 1/3. * alphaK * ( psf/(mp*c)*0.938)**(1/3.) 
     *         * 0.4 * 1.5D11 * vel
c     print*, "R_shell_x(no_shells),vel",R_shell_x(no_shells), vel
c     print*, "i, kappa_perp",i, kappa_perp(i, no_shells)
          
c     print*,"---lambda_perp, lambda_par---", lambda_perp, lambda_par
          
          
          IF ( i_perp .eq. 1) then 
             
             tempP = psf/(mp*c) * 938. ! momentum in MeV/c.
             temp_rl = get_rl(B0,tempP, lambda_c,alpha_inj,s2d,deltaB,
     X            temp_kappa_perp, temp_kappa_par, tempt)
             
             print*,"--temp_rl,temp_kappa_perp,temp_kappa_par--"

             print*, temp_rl, temp_kappa_perp, temp_kappa_par

             print*, "---- i, kappa_perp(i,no_shells), t",
     X            i, kappa_perp(i,no_shells), tempt
             
             kappa_perp(i,no_shells) = temp_kappa_perp  * 
     X            sin(theta_BN/90. * 3.14/2.)**2/ 5.  !  devide by 5 to agree with the 
c lambda = 0.4 AU * P^(1/3)*r^(2/3) assumption made in Li 2003 and 2005.
          ENDIF 
      enddo
      
c     pause
      
      Do i = 1, 250
         p = p0 + del_p * (i-1)
         psf = Exp(p)*mp*vo     ! momentum in MKS.
c     
         vel = psf/mp*(1/Sqrt(1.0 + (psf/(mp*c))**2.0)) ! in real unit.
         gamma_rel = 1.0/Sqrt(1.0 - vel**2.0/c**2.0)  
c     
         k = mp*gyrof/(psf)     ! the wave vector whose resonance p is psf.
                                ! k is the same for rel. or non-rel. case.
         
c         rg = gamma_rel*vel/gyrof ! Need gamma_rel since gyrof is Non-relativistic.	     
         rg =  1./ k
c     
         energy = Sqrt((mp*c**2.0)**2.0+(mp*Exp(p)*vo*c)**2.0)- !kinetic energy
     *        mp*c**2.0
         energy = energy/1.6E-19 ! in eV.
         
c     print*, "gamma_rel", gamma_rel, psf * 3.*1.D8/1.6D-13, i
                                !     psf is in (MeV).
c------- The weakest Iplus can NOT be smaller than  B^2/(2 \mu_0 k)
         If (k .LT. ko) Then
            Ipluso(i) = Io      ! Energy containing range.
         Else
            Ipluso(i) = Io*(k/ko)**(-5.0/3.0)
         End If          
         
c----------The strongest Iplus can NOT exceed 24 pi/32 * B^2/(2 \mu_0 k)
         Ibohm(i) = alpha_bohm * 24.* pi/32.*(B**2.0/(muo*2.0*k)) ! The Bohm approximation
         
         Ibohm(i) = Ibohm(i) /2. ! &&&&&&& ----- special treatment
         
         s(i) = k/(mp*gyrof/pinj) ! ratio of k to the K>, K> corresponds to pinj
         s_ratio = k/(mp*gyrof/pinj)
c     print*, "--- s(i) ---", s(i)
         
         If (s_ratio .GE. 1.0) Then
            Iplus(i,no_shells) = pi*beta*gyrof*injectrate*VA*
     *           cos(psi)*pinj/(Vpr**2 * k**2) *
     *           (1.0/(beta-4.0)-(psf/pinj)**2.0/(beta-2.0))
     *           + Ipluso(i)
c     print*, " --------- YOU SHOULD NOT SEE THIS ------------"
c     print*, "--- i,k > mp*Omega/pinj--",  i,k, mp*gyrof/pinj
         ELSE                   !  So, s(i) .LT. 1.0
            Iplus(i,no_shells)= 2*pi*beta*gyrof*injectrate*VA*cos(psi)*
     *           psf*(psf/pinj)**(3.0-beta)/( (beta-4.)*(beta-2.) )
     *           /(Vpr**2 * k**2) +  Ipluso(i)
            fo(i,no_shells) = (psf/pinj)**(-1.0*beta)
         ENDIF
         Iplus(i,no_shells) = Iplus(i,no_shells) / 2.0 ! &&&&&&& ----- special treatment
         
c          print*," i,Iplus(i,no_shells),Ibohm", 
c     *         i,Iplus(i,no_shells),Ibohm(i)

c          print*," k, ko", k, ko

          If (Iplus(i,no_shells)/(Ibohm(i)+Ipluso(i)).GT. 1.0) Then
c             print*," i,Iplus(i,no_shells),Ibohm", 
c     *            i,Iplus(i,no_shells),Ibohm(i)
             Iplus(i,no_shells) = Ibohm(i)+Ipluso(i)
          End If
          
          kappa(i,no_shells) = psf**2.0*vel/
     *         (ech**2.0*2.0*muo*Iplus(i,no_shells))
       End Do

c-------Need to decide the maximum momentum carefully.

       time_all = x1a(shock_end+5) *1.5D11 / (shock_speed*vo)       
       print*, " ---time_all is--", time_all

       print*, "no of shells", no_shells

       if (no_shells .eq. 1) then
          time_initslap =  ( x1a(is) - shock_init_loc)  ! 0.025 ==> 4 solar radii 
     *      *1.5D11 / (shock_speed*vo)
       endif

       time_all = (x1a(shock_end+5) - x1a(is) ) 
     *      *1.5D11 /(shock_speed*vo)  + time_initslap

       print*, " time_all is", time_all

       time_all = (time - shock_start_time) * to + time_initslap

       print*, " time_all is", time_all

       acc_time = 0.
       
       print*, "-------time all is ********", time_all

c       pause

       p_max = 0.
       
c       psf = Exp(p)*mp*vo       ! momentum.

       if (i_perp .ne. 1) then
          goto 16
       endif

       DO i = 1, 250
          p = p0 + (i-1)*del_p   
          
          if ( Exp(p) .ge. p_inj) then
             acc_time = acc_time+2.* beta/Vpr**2.0*
     *            kappa_perp(i,no_shells)* del_p ! See Drury's paper for the factor 2.
c     *            (exp(p+del_p/2.0D0)-exp(p-del_p/2.0D0))/exp(p)
c             print*, "i,p, acc_time,time_all", i,p, acc_time, time_all
c             print*, "beta,kappa(i,1),Vpr", beta,kappa(i,1), Vpr

             if ( (acc_time .gt. time_all) .and. (p_max .eq. 0.)) then
                p_max = Exp(p)
c                print*, "p_max, p_cut",p_max,p_cut
                p_cut = p_max
                j_p_cut = i
                goto 16
             endif
c             print*, "acc_time", acc_time
          endif
          
       ENDDO


 16    CONTINUE
        print*, "p_cut",p_cut
c        pause


       IF ( i_perp .eq. 1) then 
          goto 17
       ENDIF

       DO i = 1, 250
          p = p0 + (i-1)*del_p   
          
          if ( Exp(p) .ge. p_inj) then
             acc_time = acc_time+2.* beta/Vpr**2.0*kappa(i,no_shells)*  ! See Drury's paper for the factor 2.
     *            del_p
c     *            (exp(p+del_p/2.0D0)-exp(p-del_p/2.0D0))/exp(p)
             print*, "i,p, acc_time,time_all", i,p, acc_time, time_all
c             print*, "beta,kappa(i,1),Vpr", beta,kappa(i,1), Vpr

             if ( (acc_time .gt. time_all) .and. (p_max .eq. 0.)) then
                p_max = Exp(p)
c                print*, "p_max, p_cut",p_max,p_cut
                p_cut = p_max
                j_p_cut = i
                goto 15
             endif
c             print*, "acc_time", acc_time
          endif
          
       ENDDO
       
 15    continue
c       pause
c-------Since no spectra transfer involves, Reset Iplus(i,1) to Iplus0(i)
c-------when p<p_inj or p>p_cut.
c-------Reset the kappa(i,1) after reset Iplus(i,1)

       DO i = 1, 250
          p = p0 + (i-1)*del_p   
          psf = Exp(p)*mp*vo    ! momentum.
          vel = psf/mp*(1/Sqrt(1.0 + (psf/(mp*c))**2.0)) ! in real unit.          
          
          If ((Exp(p) .LT. p_inj) .OR. (Exp(p) .GT. p_cut)) Then
             Iplus(i,no_shells) = Ipluso(i)
          ENDIF
          
          kappa(i,no_shells) = psf**2.0*vel/
     *         (ech**2.0*2.0*muo*Iplus(i,no_shells))
       ENDDO

 17    continue

c-------Since Iplus is known for all $k$ now, we next 
c-------find out kappa_hi(i,no_shells). Note, the resonant index j for 
c-------kappa_hi(i,no_shells) is given by (Amass/charge) * p_proton. 
c-------Also note there is a (A/Q)**2 dependence of kappa for heavy Ions.
       
       energy_inj =  Sqrt((mp*c**2.0)**2.0+(mp*vo*p_inj*c)**2.0)- !kinetic energy
     *      mp*c**2.0
       energy_inj = energy_inj/1.6E-19 ! in eV.
       
       print*, "The injection energy in eV is", energy_inj
c       pause !GL 

       DO i = 1, 250
          p = p0 + (i-1) *del_p
          psf = Exp(p)*mp*vo    ! momentum per nucleon.
          vel = psf/mp*(1/Sqrt(1.0 + (psf/(mp*c))**2.0)) ! in real unit.          
          
          energy = Sqrt((mp*c**2.0)**2.0+(psf*c)**2.0)- !kinetic energy
     *         mp*c**2.0
          energy = energy/1.6E-19 ! in eV.
          
          psf_hi = Exp(p) * (Amass/charge) ! psf_hi is the resonant proton momentum for heavy ions.          
          k = mp*gyrof /psf_hi

c     print*, "psf_hi,  p", psf_hi, Exp(p)
c     print*, "(Amass/charge), Log",(Amass/charge),
c     *         DLog(Amass/charge)

          index = i + INT(( DLog(psf_hi)-p)/del_p + 0.5)
          
          if (index .ge. 250) then  ! the corresponding turbulence must be due to solar wind.
             temp_p = Exp(p0 + index*del_p) * mp * vo 
             temp_k = (charge/Amass)*(ech * B)/ temp_p 
             temp_k = mp*gyrof /temp_p

             if ( temp_k .lt. ko) then
                temp_I = Io
             else 
                temp_I = Io*(temp_k/ko)**(-5.0/3.0)
             endif
             kappa_hi(i,no_shells) = psf**2.0*vel/
     *            (ech**2.0*2.0*muo*temp_I) * (Amass/charge)**2.
c---- ( debug
c             print*, "i,index,kappa_hi(i,no_shells)",
c     *            i,index,kappa_hi(i,no_shells)
c---- ) debug

          else
             kappa_hi(i,no_shells) = psf**2.0*vel/
     *            (ech**2.0*2.0*muo*Iplus(index,no_shells))          
     *            * (Amass/charge)**2.
          endif
       ENDDO

c-----Now decide the maximum heavy ion momentum.
       
       p_hi_cut2 =  p_cut * (charge/Amass)
       p_max_hi = 0.
       acc_time_hi = 0.
       
       p_hi_inj = p_inj


       DO i = 1, 250
          p = p0 + (i-1)* del_p   
          
c     print*,"kappa(i,1),kappa_hi(i,1),i",kappa(i,1),kappa_hi(i,1),i
          if ( Exp(p) .ge. p_inj) then 
             acc_time_hi = acc_time_hi + beta/Vpr**2.0*kappa_hi(i,
     *            no_shells)*
     *            (exp(p+del_p/2.0D0)-exp(p-del_p/2.0D0))/exp(p)
             
             if ( (acc_time_hi .gt. time_all) 
     *            .and. (p_max_hi .eq. 0.)) then
                p_max_hi = Exp(p-del_p)
                p_hi_cut = p_max_hi
                goto 20
             endif
          endif
       ENDDO
       
 20    CONTINUE

c----( debug
c       print*, "i,acc_time_hi,time_all", i,acc_time_hi, time_all
c       print*, "p_hi_cut, p_hi_cut2", p_hi_cut, p_hi_cut2
c----) debug


       print*, "The proton injection and maximum momentum (in mp*c)", 
     *      p_inj * vo/c, p_cut*vo/c
       print*, "The hvyIon injection and maximum momentum (in mp*c)", 
     *      p_hi_inj*vo/c, p_hi_cut*vo/c
       print*, "*** p_hi_cut, p_hi_cut2 ***", p_hi_cut, p_hi_cut2

       write(90,*) "***********************************"
       write(90,*) "*** no_shells is ***", no_shells
       write(90,*) "*** p_hi_cut, p_hi_cut2 ***", p_hi_cut, p_hi_cut2
       write(90,*) "*** Normal case: p_hi_cut2 <p_hi_cut ***"
       write(90,*) "***********************************"
       if ( p_hi_cut2 .lt. p_hi_cut) then
          p_hi_cut = p_hi_cut2
       endif

       j_p_cut_hi = INT (0.5+(Dlog(p_hi_cut) - p0) / del_p) + 1

c       print*, " j_p_cut_hi", j_p_cut_hi
c       print*, " j_p_cut", j_p_cut
c       print*, "(Dlog(p_cut) - p0) / del_p ",(Dlog(p_cut) - p0) / del_p 
       j_p_cut = INT ( 0.5+(Dlog(p_cut) - p0) / del_p) + 1
c       print*, " j_p_cut", j_p_cut

c       pause

       write(90,*) "*** p_hi_cut2 > p_hi_cut1 ***"

       energy = Sqrt((mp*c**2.0)**2.0+(mp*p_cut*vo*c)**2.0)- !kinetic energy
     *      mp*c**2.0
       energy = energy/1.6E-19  ! in eV.
       
       print*, "p_inj,p_cut, energy_cut",p_inj,p_cut, energy
       
c       pause

c---- write out data at the shock location.

       DO i = 1, 250
          p = p0 + (i-1)* del_p   
          psf = Exp(p)*mp*vo    ! momentum.
          k = mp*gyrof/psf
c          
          vel = psf/mp*(1/Sqrt(1.0 + (psf/(mp*c))**2.0)) ! in real unit.
          gamma_rel = 1.0/Sqrt(1.0 - vel**2.0/c**2.0)  
          rg = gamma_rel*vel/gyrof 

          energy = Sqrt((mp*c**2.0)**2.0+(mp*Exp(p)*vo*c)**2.0)- !kinetic energy
     *         mp*c**2.0
          energy = energy/1.6E-19 ! in eV.
          If (s(i) .NE. 0.0) Then
             Write(43,111) k*Vsw/2.0/pi, energy,
     *            2.0E18*Ipluso(i)*2.0*pi/Vsw, 
     *            2.0e18*Iplus(i,no_shells)*2.0*pi/Vsw, 
     *            2.0E18*Ibohm(i)*2.0*pi/Vsw
             Write(44,112) k*Vsw/2.0/pi,energy, kappa_perp(i,no_shells), 
c     *            kappa_hi(i,no_shells),4.*pi/3.*rg*vel/alpha_bohm,
     *            kappa(i,no_shells),
     *            4.*pi/3.*rg*vel/alpha_bohm,x1a(shock_begin)
          Else
             Write(43,111) 0.0, 0.0, 0.0, 0.0, 0.0
             Write(44,112) 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
          End If
       ENDDO

       ipmin(no_shells) = INT(( DLog(p_inj)- 2.0)/del_p ) + 2
       ipmax(no_shells) = INT(( DLog(p_cut)- 2.0)/del_p ) + 1
       ipmin_hi(no_shells) = INT(( DLog(p_hi_inj)- 2.0)/del_p ) +2
       ipmax_hi(no_shells) = INT(( DLog(p_hi_cut)- 2.0)/del_p ) +1
       
        print*, " ipmin, ipmax ", ipmin(no_shells), ipmax(no_shells)
       print*, " ipmin_hi, ipmax_hi ", 
     *      ipmin_hi(no_shells), ipmax_hi(no_shells)

       Do j = ipmax(no_shells) - 10, ipmax(no_shells) + 10
          
          print*, " j, kappa(ipmax(j,no_shells )", 
     *         j,kappa(j, no_shells )
       Enddo

c       pause

c----( debug       
c       do j = ipmin(no_shells) - 1, ipmax(no_shells) + 1
c          print*, " j,kappa(j,no_shells) ", j,kappa(j,no_shells)
c       enddo
c       do  j = ipmin_hi(no_shells) -1 , ipmax_hi(no_shells) +2
c          print*, "j, kappa_hi(j,no_shells) ", j,kappa_hi(j,no_shells)
c       enddo
c----) debug       

c--------- Reset kappa(j,no_shells) to happa_perp(i,no_shells) for later diffusion.
c--------- Gang Li, Sep. 19th, 2005.
c--------- Need to do the same for heavy ions. (need Q/A dependence).

       if (i_perp .eq. 1) then 
          Do i = 1, 250
             kappa(i,no_shells) = kappa_perp(i,no_shells)
          Enddo
       endif

       print*, " -- time is --", time 

c---- Increase the 

       print*, "---***--- Inc. the turb. power acr. the shock ---**---"
       Do j = 1, 250
          Iplus(j,no_shells) =   Iplus(j,no_shells) * (comp_rat**2)
          kappa(j, no_shells) =  kappa(j, no_shells)/ (comp_rat**2)
          kappa_hi(j,no_shells)= kappa_hi(j,no_shells)/(comp_rat**2)
       ENDDO


       print*, "******* Leaving wave_intens.f *******"
       print*, "---comp_rat ----", comp_rat
c       pause
       print*, "                                       "
       return
       
 111   format (5 (2x,E9.4) )
 112   format (6 (2x,E9.4) )

       end

      Function get_injection_speed(theta,deltaB,B0,u,comp,
     X     alpha, s,lambda_c, lambda_perp, lambda_par)      

      real*8 get_injection_speed
c ---------------------------------------------------------------------------
c----        by Gang Li
c----        Sep., 13th, 2005.
c----        This code returns the injection speed for perp. shock case.
c----        See Zank et al. 2004.
c -------------------   Input    -----------------------------------------
c     -------- theta_BN == the angle between shock normal and B0 in degree.
c     -------- deltaB == the turbulent field delta_B.
c     -------- B0 == ambient B.
c     -------- u == flow speed.
c     -------- alpha == lambda_2D / lambda_slab 
c     -------- s == (b_2D/b)^2
c -------------------   Output  ------------------------------------------
c     -------- injection momentum p_inj  (in MeV/c ?)
c----   Code description.
c----
c----   y = lambda_{xx} / r_l
c----   x = rl/lambda_c
c----   rl/lambda_c. 

      integer j
      
      real*8 theta,deltaB, B0,comp, u,alpha,s,lambda_c
      real*8 v_inj              ! the return value of momentum in SGI unit.
      
      real*8  beta, phi
      real*8  oneAU
      
      real*8  lambda_par, lambda_perp
      real*8  kesi, kesi1, kesi2, kesi_pre
      real*8  t,rl,x, sg, q, A, D, z,y,v
      real*8  ratio, factor, solution

c ----- need the following input

      oneAU = 1.5D8             ! 1AU =1.5E8 km      
      phi =  theta/90. * (3.14/2.) ! theta in radian.
      beta = (deltaB/B0)**2 ! beta is the ratio of deltaB/B0

      kesi_pre = 10.            ! initialize solution for kesi


c      alpha = 0.05               ! = lambda_2d/lambda_SLAB
c      beta = 1.                 ! = b^2/B^2_0
c      s = 0.8                   ! = b^2_2d/b^2

c      oneAU = 1.5E8             ! 1AU =1.5E8 km
c      u  = 400.  ! flow speed 400 km/s
c      comp_ratio = 4. ! compression ratio 

c      print*, " alpha, beta, s, u, comp"
c      print*, alpha, beta, s, u, comp
c      print*, "B0,deltaB, lambda_c", B0,deltaB,lambda_c

c ----- need the following input

c      B0 = 4.12E-5  ! magnetic field at 1 AU
c      lambda_c = 0.03 * 1.5 E8 ! correlation length at 1 AU in km.

      do j = 0, 1300  
         
         t =  1D-4 * exp(j* 0.01) ! t is the kinetic energy in MeV
         rl = 31. /B0 * (2.13E-3 * t * (1+t/(2.*938.) ) )**0.5 ! Larmor radius in km.
         
         x = rl /lambda_c       ! normalized rl
         
         sg = 0.7468 * x
         q = 5. /3. * sg**2 / (1+ sg**2 -(1+sg**2)**(1/6.))
         A = (1+sg**2)**(5/6.) - 1.
         D = 7./9. * A / ( (q + 1./3.) * (q + 7./3.)  )
         
         z = 2.343 /( (1-s) * beta ) * x**(1/3.) * (1+D)
         
         if( z/1.732 .lt. 1.) then
            factor = 4.33 * z/1.732
         else
            factor = 3.09
         endif
         
         y = 0.36 * (s * alpha * beta)**(2/3.)*z**(1/3.)*(1.339/x)
     X        *(1+ 0.1102*(1-s) *beta**(1/3.) / (alpha *s)**(2/3.)
     X        /(z**(1/3.))* factor ) 
         
         lambda_perp = y        ! lambda_perp / rl
         lambda_par =  z * (1.339/x) ! lambda_par / rl
         
         ratio = lambda_perp / lambda_par
         
c     --now calculate kesi.
         
         v = ((t/938.) * ( 2+ t/938.))**0.5/(1+t/938.)
     X        * 3.E5            ! particle speed in km/s
         
c         write (30,*) t, rl/OneAU,lambda_perp*rl/OneAU,
c     X        lambda_par*rl/OneAU, D            
         
         kesi  = 3 * u /v  *( 
     X        1 / (comp_ratio-1)**2 + 
     X        (1+ ((lambda_par - lambda_perp)*cos(phi))**2)
     X        *sin(phi)**2 
     X        / ( lambda_perp*sin(phi)**2 
     X        + lambda_par*cos(phi)**2 )**2 )**0.5 
         
         if ( (kesi_pre .gt. 1.) .and. (kesi .le. 1.)) then
            solution = t
            kesi1 = kesi_pre
            kesi2 = kesi
            goto 10
         endif
         kesi_pre = kesi
      enddo
         
 10   print*, "***solution***",solution, kesi1, kesi2

      v_inj = ((solution / 938.) * (2 + solution / 938.))**0.5/
     X     (1 + solution / 938.) ! injection speed in unit of c.
      
      get_injection_speed = v_inj * 3.D5 ! injection speed in unit of km/s.
      lambda_perp = lambda_perp * rl ! lambda_perp in km.
      lambda_par =  lambda_par * rl ! lambda_par in km.

      return 
      
      end
      

      Function get_rl(B0,tempP, lambda_c, alpha, s, deltaB,
     X     temp_kappa_perp, temp_kappa_par, tempt) 

c---- B0  == magnatic field in gauss. 
c---- tempP ==  momentum in MeV/c.
c---- lambda_c == correlation length in km.
c---- alpha == lambda_2d/lambda_SLAB
c---- s == b^2_2d/b^2
c---- deltaB == turbulent magnetic field in gauss.
      
      real*8 tempP, tempt
      real*8 get_rl, rl, temp_kappa_perp, temp_kappa_par
      real*8 deltaB, B0,alpha,s,lambda_c, beta
      real*8 t, x, sg, q, A, D, z, y, v
      

      beta = (deltaB/B0)**2 ! beta is the ratio of deltaB/B0
      t = (tempP**2 + 938**2)**0.5 - 938.  ! particle energy in MeV.
      v = ((t/ 938.) * (2 + t / 938.))**0.5/
     X     (1 + t / 938.) * 3.D5 ! particle speed in km/s.
      
      rl = 31. /B0 * (2.13E-3 * t * (1+t/(2.*938.) ) )**0.5 ! Larmor radius in km.
      
      x = rl /lambda_c          ! normalized rl
      
      sg = 0.7468 * x
      q = 5. /3. * sg**2 / (1+ sg**2 -(1+sg**2)**(1/6.))
      A = (1+sg**2)**(5/6.) - 1.
      D = 7./9. * A / ( (q + 1./3.) * (q + 7./3.)  )
      
      z = 2.343 /( (1-s) * beta ) * x**(1/3.) * (1+D)
      
      if( z/1.732 .lt. 1.) then
         factor = 4.33 * z/1.732
      else
         factor = 3.09
      endif
      
      y = 0.36 * (s * alpha * beta)**(2/3.)*z**(1/3.)*(1.339/x)
     X     *(1+ 0.1102*(1-s) *beta**(1/3.) / (alpha *s)**(2/3.)
     X     /(z**(1/3.))* factor ) 

c      print*, " ----y,z, v-----", y,z, v
      

      get_rl = rl
      temp_kappa_perp = y * rl * v / 3. * 1.D6 ! kappa_perp in m^2/s
      temp_kappa_par =  z * (1.339/x)  * rl * v /3. * 1.D6 ! kapa_par in m^2/s
      temp_t = t

      return 
      end
