      subroutine diffuse_distr
c
c    abcd:zeus3d.distr_evolve <--- Diffuses the already created particles
c          that have already undergone expansion in the shells.  The particles
c          are able to diffuse into neighbouring shell.
c
c                                                        DATE 07/01/02
c
c    written by: Ken Rice
c    modified 1:
c
c    modified by: Gang Li 02/21/2003    
c
c    modified by: Gang Li 03/11/2003: Add heavy ions
c    modified by: Gang Li 01/23/2004: 
c     1) Move code to after add_acc_ptcls.f in zeus3d.f
c     2) Write out current J = num /(MeV cm^2 s) at the shock front
c        to file dist_shock.dat (file 78) in the format
c        e, j, j_hi, r, t
c     3) Write out the content of distr_array(no_shells, p_int) to
c        file (96).
c     4) write out esc_number_up(p_int), x1a(shock_begin)+4.0D0*temp,
c        and time to file (88).
c     5) write out esc_number_up_hi(p_int), x1a(shock_begin)+4.0D0*temp,
c        and time to file (89).
c        
c
c
c  PURPOSE:
c
c
c  LOCAL VARIABLES:
c
c  EXTERNALS:
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-----------------------------------------------------------------------
c-------------- START OF PROGRAM/SUBROUTINE ----------------------------
c-----------------------------------------------------------------------
c
c      
       integer   maxshells, p_number
       real*8    max_p, mp,vo, p0 
c
       parameter (maxshells = 200, max_p = 10.0D0, p_number = 250, 
     *      p0 = 2.0, vo = 52483.25D0, mp = 1.67D-27) 
c
       real*8   get_number
       real*8      R_shell_x(maxshells)
       real*8      del_p, p, fi, num_ptcls  
       real*8      shock_speed, Machno, comp_rat            
       real*8      p_inj, p_cut, distr_array(maxshells, p_number)   
       real*8      new_ptcl_num(in), temp, temp0
       real*8      diffuse_array(maxshells, p_number)
       real*8      esc_dens_up(p_number), esc_dens_down(p_number),
     *             prev_acc_time, prev_expand_time(maxshells),
     *             prev_acc_posn
       real*8      esc_shell(in, p_number), dens_ptcls
       real*8      total_num, total_num_shells, total_num_esc_shell 
       real*8      shell_posn, esc_length, shell_vol
       real*8      kappa(p_number,maxshells)
       real*8      kappa_hi(p_number, maxshells)
       real*8      kappa_norm, Energy, vel
         
       real*8      esc_dens_shell(p_number, maxshells)
       real*8      esc_dens_shell_hi(p_number, maxshells)
       integer     shock_begin, shock_end, shell_posn_int
       integer     i , j , k, l, R_shell_i(maxshells), no_shells
       integer     p_int
       
       real*8      escape_location(p_number)
       real*8      vpar           ! particle speed.
       
c============= 02/02/2007, GL
       real*8      temp_esc(p_number), temp_esc_hi(p_number)
       real *8 delta_out
c============= 02/02/2007, GL       

c---- for heavy ions.

       real*8   charge, Amass, p_hi_inj, p_hi_cut
       real*8   distr_array_hi(maxshells,p_number)
       real*8   esc_dens_up_hi(p_number), esc_dens_down_hi(p_number)

       real*8   esc_shell_hi(in,p_number)
       
       real*8  esc_number_up(p_number), esc_number_up_hi(p_number)
       real*8  esc_number_down(p_number), esc_number_down_hi(p_number)
       real*8  temp1, temp2
       
       real*8  escape_location_hi(p_number)

       real*8  p_cut_old, p_hi_cut_old
       real*8  p_cut_expand, p_hi_cut_expand 

       real*8  test_alpha, test_fac, test_p

       real*8  factor_num 
       integer index1,index2, index3

       real *8 dV, den_sh, fi_hi, den_sh_hi ! den_sh is the density at shock front.
       real *8 den_sh_j, den_sh_hi_j

       real *8 energy_low, energy_hi
       real *8 Vsh
       real*8 totnum,tempR
       real*8 total_num_p(maxshells),total_num_hi(maxshells)
       
       real*8 kappa_norm0

       real*8 temp_counter1, temp_counter2

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

       real*8 c 

       integer escape_shell_num(1:250)
       integer escape_shell_num_hi(1:250)

       real*8 ipmax(maxshells),ipmin(maxshells)
       real*8 ipmax_hi(maxshells),ipmin_hi(maxshells)
       
       integer ichoice, ii

       real*8 temp_num(1:p_number),temp_num_hi(1:p_number), temp_dV


       integer i1
       real*8 temp_factor

       real*8 distr_all(250)

       real *8 temp_esc_up, temp_esc_up_hi

       real *8 tttt, ttt

       real *8 temp_check(250), temp_check2(250)

       common /p_checker/ipmax,ipmin,ipmax_hi,ipmin_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 /diff_coeff/ kappa,kappa_hi

       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 /diff_add_expand_ptcl/esc_number_up, esc_number_up_hi

       common /diff_expand/ p_cut_old, p_hi_cut_old,  p_cut_expand, 
     *      p_hi_cut_expand 

       common /escape_counter/ esc_ptcl_counter, esc_ptcl_counter_hi

       common /diag/ distr_all
c       
c-----------------------------------------------------------------------
c                            

       del_p = (max_p - p0)/(p_number*1.)  ! step of log(p) (p is in measure of p/mp) 
       
       print*, "--------Entering diffuse_distr.f --------"
       
       c = 3.0D8 

       Do j = 1, p_number
          p = p0 + (j-1)* del_p   
          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.
          write(28,*) Energy, kappa(j,no_shells)
       ENDDO

c----( debug 
c---- total_num_p(i) and total_num_hi(i) are the total proton number 
c---- and heavy ion number in the ith shell (integrated over all momentum).

        Do i = 1, no_shells
           total_num_p(i) = 0.
           total_num_hi(i) = 0.
          Do j = 1, p_number
             p = p0 + (j-1) * del_p
             total_num_p(i) = total_num_p(i) +  
     *            distr_array(i,j) *4./3. * pi *
     *            (Exp(p+del_p)**3 - Exp(p)**3)
             total_num_hi(i) =  total_num_hi(i) +  
     *            distr_array_hi(i,j) *4./3. * pi *
     *            (Exp(p+del_p)**3 - Exp(p)**3)
             escape_location(j) = 0.
             escape_location_hi(j) = 0.
          End Do
c          print*, " ---- before diffuse -----"
c          print*, " -- i, total_num_p(i) --", i, total_num_p(i)
c          print*, " -- i,total_num_hi(i) --", i, total_num_hi(i)
       End Do
   

c---- ) debug

c---- ( debug

       write(97,*)"--- time and r_shell_x(no_shells+1) ---", 
     * time, r_shell_x(no_shells+1)   

       Do j = 1, p_number
          temp_num(j) = 0.
       ENDDO
       temp_dV = 0.
       
       DO i =1, no_shells 
          DO j =1, p_number
             p = p0 + del_p * (j-1)
             Energy = Sqrt((1.67E-27*3.0E8**2.0)**2.0 + 
     *            (1.67E-27*Exp(p)*52483.25*3.0E8)**2.0) -
     *            1.67E-27*3.0E8**2.0 ! the kinetic energy.
             Energy = Energy/1.6E-19 !  (Energy in eV)
             factor_num = 4.*pi/3.* (Exp(p)**3 - Exp(p-del_p)**3) 
             dV = 4.*pi/3. *(  R_shell_x(i+1)**3 
     *         - R_shell_x(i)**3 ) ! no need to convert R_shell_x.
c             write(97,*) Energy,distr_array(i,j)* factor_num/dV,
c     *            distr_array_hi(i,j)*factor_num/dV

             temp_num(j) = temp_num(j) + distr_array(i,j)
             temp_dV = temp_dV + dV
             write(97,*) Energy,distr_array(i,j),
     *            distr_array(i,j)* factor_num/dV
c             print*, "j,energy,distr_array(no_shells,j)",
c     *            j,energy,distr_array(no_shells,j)
          ENDDO
       ENDDO


c----- write out time-integrated accelerated particles. 
c----- should be identical to fort.94 from add-acc_ptcls.f.

       Do j = 1, p_number
          temp_check(j) = 0.0
          temp_check2(j) = 0.0
       Enddo

       Do j = 1, p_number
          Do i = 1, no_shells ! i = no_shells + 1 for ptls in front of the shock
             temp_check(j) = temp_check(j) + distr_array(i,j)
          End Do
       End Do                   
       
       Do j = 1, p_number
          temp_check2(j) = temp_check(j) + esc_ptcl_counter(j)
       ENDDO

       Do j = 1, p_number
          p = p0 + del_p * (j-1)
          Energy = Sqrt((1.67E-27*3.0E8**2.0)**2.0 + 
     *         (1.67E-27*Exp(p)*52483.25*3.0E8)**2.0) -
     *         1.67E-27*3.0E8**2.0 ! the kinetic energy.
          Energy = Energy/1.6E-19 !  (Energy
          write(83,*) Energy, temp_check(j),temp_check2(j),distr_all(j),
     *         time
       Enddo

c---- ) debug

c---- Now we need to do something about ipmax(i=1:no_shells).
c---- The idea is that ipmax(i) must >= ipmax(i+1), for every i = 1 : no_shells -1 .
c---- This is to say, if the shock is still accelerating particles to higher energy, then
c---- turbulence will develop in all the previous shells as well.

       print*, "no_shells,ipmax(no_shells)", no_shells, ipmax(no_shells)
c       pause
       
       Do i1 = 1, no_shells -1  ! this DO END block is a rare case
          i = no_shells - i1

          if(ipmax(i) .lt. ipmax(i+1) ) then
             temp_factor = ( (R_shell_x(i+1)+R_shell_x(i))/
     *            (R_shell_x(i-1)+ R_shell_x(i)) )**3
             Do j = ipmax(i) + 1, ipmax(i+1)
c                kappa(j, i) = kappa(j, i+1) * temp_factor
                kappa(j, i) = kappa(j, i+1)
             Enddo
             ipmax(i) = ipmax(i+1)
          endif
       Enddo


       Do i1 = 1, no_shells -1  ! this DO END block is a rare case
          i = no_shells - i1

          if(ipmax_hi(i) .lt. ipmax_hi(i+1) ) then
             temp_factor = ( (R_shell_x(i+1)+R_shell_x(i))/
     *            (R_shell_x(i-1)+ R_shell_x(i)) )**3
             Do j = ipmax_hi(i) + 1, ipmax_hi(i+1)
c                kappa_hi(i, j) = kappa_hi(i+1, j) * temp_factor
                kappa_hi(i, j) = kappa_hi(i+1, j)
             Enddo
             ipmax_hi(i) = ipmax_hi(i+1)
          endif
       Enddo

       Do i =1, no_shells
          print*, " i, ipmax(i)", i, ipmax(i)
       ENDDO
c       pause

c---- initialization of escape_shell_num(j).
c---- Find the escaping location for different energies.

       Do j = 1, 250
          escape_shell_num(j) =  no_shells
          escape_shell_num_hi(j) =  no_shells
       ENDDO

       Do j = 1, 250
          Do i = 1, no_shells
             if ( j .le. ipmax(i) ) then
                escape_shell_num(j) = i
             endif

             if ( j .le. ipmax_hi(i) ) then
                escape_shell_num_hi(j) = i
             endif
          ENDDO
       ENDDO

       Do j =1, 250
          escape_location(j) =  R_shell_x(escape_shell_num(j) + 1) 
          escape_location_hi(j) = R_shell_x(escape_shell_num_hi(j)+ 1)
       ENDDO


c-- Alternatively, use the following simplified estimate ---------
c       Do j = 1, 250
c          escape_shell_num(j) =  no_shells
c          escape_shell_num_hi(j) =  no_shells
c          escape_location(j) =  R_shell_x(no_shells+1)
c          escape_location_hi(j) = R_shell_x(no_shells+1)
c       ENDDO
c-----------------------------------------------------------------

c---- Initialize the escaped particle ( Not used )
       do ii = 1, maxshells
          esc_dens_shell(p_int, ii) = 0.
          esc_dens_shell_hi(p_int, ii) = 0.
       enddo

c---- ( debug
       write(84,*)"time is", time
       Do j = 1, 250
c          print*, "j, escape_shell_num(j), no_shells", j, 
c     *         escape_shell_num(j),no_shells

          write(84,*) "j, escape_shell_num(j), no_shells", j, 
     *         escape_shell_num(j),no_shells

       ENDDO
c       pause

c---- ) debug

c---  now calculate the corresponding index for p_cut_old, p_cut_expand,p_cut
       
c       if ( p_cut .ne. 0.) then
c          index1 =  1 +  INT(( DLog(p_cut_old)- p0)/del_p + 0.5)
c          index2 =  1 +  INT(( DLog(p_cut_expand)- p0)/del_p + 0.5)
c          index3 =  1 +  INT(( DLog(p_cut)- p0)/del_p + 0.5)
c       else
c          index1 = 0
c          index2 = 0
c          index3 = 0
c       endif
c       print*, " *** index1, index2,index3 *** ", index1, index2, index3

c       pause

c                          
       Do i = 1, no_shells + 1  ! i = no_shells + 1 for ptls in front of the shock
          Do j = 1, p_number
             diffuse_array(i,j) = 0.0D0 ! The new array for the distribution  
                                ! after the diffusion.
             esc_dens_up(j) = 0.0D0 ! escaped distribution at upstream.
             esc_dens_down(j) = 0.0D0 ! escaped distribution at downstream.
          End Do
       End Do                   
c     
c--- ( del
       Do i = 1, in             ! # of position.
          Do j = 1, p_number         ! # of momentum.
             esc_shell(i, j) = 0.0D0
          End Do           
       End Do
c--- ) del

c----- Initialization of particles in up and dn stream.
       Do i = 1, p_number
          esc_number_up(i) = 0.
          esc_number_up_hi(i) = 0.
          esc_number_down(i) = 0.
          esc_number_down_hi(i) = 0.
       ENDDO
       

c===== added by Gang Li  02/02/2007

       do j = 1, p_number
          temp_esc(j)  = 0.
          temp_esc_hi(j)  = 0.
       enddo


c       delta_out = 1.7183  ! delta_out controls how many particles leak out when the shell is first generated.
       delta_out = 0.25  


       print*, " --- no_shells --", no_shells
       write(99,*) " time is ", time





c---- Now the real stuff...

       Do i = 1, no_shells  ! For the current shell and all the shells generated earlier.
          
          shell_posn = (R_shell_x(i + 1) + R_shell_x(i))/2.0D0 ! the loc. of the center of the ith shell.
          
          shell_posn_int = (R_shell_i(i+1)+R_shell_i(i))/2.0D0 ! the index of the center of the ith shell.
c     
          shell_vol = 4.0*pi*(R_shell_x(i+1)**3.0D0-
     *         R_shell_x(i)**3.0D0)/3.0D0 ! volume of the shell.    
c            
          print*, " shell no is ", i                                          
          write(99,*) " shell no is ", i                                          

          Do j = 1, p_number
             p_int = j
             p = p0 + (j-1) * del_p
             
             Energy = Sqrt((1.67E-27*3.0E8**2.0)**2.0 + 
     *            (1.67E-27*Exp(p)*52483.25*3.0E8)**2.0) -
     *            1.67E-27*3.0E8**2.0 ! the kinetic energy.
             
             Energy = Energy/1.6E-19 !  (Energy in eV)
c     
             vel = Exp(p)*52483.25/Sqrt(1.0 + Exp(p)*Exp(p)*
     *            52483.25*52483.25/(3.0E8*3.0E8)) ! particle velocity.
                                ! note v/c = (p/mc) / sqrt( 1 + (p/mc)^2 ) is in m/s.

             vel = vel/52483.25D0 ! Normalized velocity (in unit of the sound speed at 1 AU).

c             print*, "i,j,escape_shell_num(j)", i,j,escape_shell_num(j)
c             print*, "kappa(j,escape_shell_num(j)),kappa(j,i)",
c     *            kappa(j,escape_shell_num(j)), kappa(j,i)
            
             fi = distr_array(i, p_int) ! number of particles before the diffusion.

C            /added by Olga 
c	       if (no_shells .gt. 22) then
c	           delta_out =0.
c	       endif
C            /added by Olga 


c===========  add by Gang Li for a better approximation for escape 


             if ( i .eq. no_shells) then 
                temp_esc(p_int) =  fi  * 1. / (1.+delta_out)
                fi = (1. - 1./(1.+delta_out) )  * fi
             endif

c===========  end add by Gang Li

 

c----   should change kappa(j,1) to kappa(j, no_shells-i)
             
             kappa_norm = kappa(j,i)*2858055.0D0/1.5D11/1.5D11 ! get the normalised kappa.
                                ! kappa_norm is in the unit of AU * v0 where v0 
                                ! is the normolization speed.
             
             temp0 = 4.0D0*kappa_norm*(time - prev_acc_time) ! the length scale 4 
                                                            !the Green's function.
c             temp = temp0 + (kappa_norm/(shock_speed - 
c     *            v1(shock_end+2,1,5)) )**2.0D0 ! take into account the initial extended 
                                !   of the assumed delta function ditribution. 

             temp = temp0 + ((R_shell_x(i + 1) - R_shell_x(i))/4.)**2
                                ! temp is the [distance^2] in the code.

c --- ( debug
c             if (i .eq. no_shells) then
c                print*, "i,j,temp0, temp", i,j,temp0,temp
c                write(96,*) "i,j,temp0, temp", i,j,temp0,temp
c             endif 
c --- ) debug

c     if (i .eq. no_shells) then
c             write(99,*) j,energy,temp, fi/shell_vol
c     endif
c     
c     If (temp**0.5 .gt. 0.5D0*vel*(time-prev_acc_time)) Then ! Ensure ptcl does not 
c                temp = (0.5D0*vel*(time-prev_acc_time))**2. ! diffuse further than it
c                print*, " **** ballistically propagation happens ****"
c             EndIf              ! could propagate ballistically. 
c                             !  SHOULD change 0.5  to 1.0. 
                
             
             If (fi .NE. 0.0D0) Then ! fi = distr_array(i, p_int)
                                ! There are particles in shell j for momentum p_int.
c---  ( debug
c                print*, " --- start diff. [1 escape_shell_num(j) ] ---"
                
c                if ( p_int .ge. 150) then
c                   write(96,*),"---start diff.[escape] ---"
c                   write(96,*) "i,j,temp0, temp", i,j,temp0,temp
c                   write(96,*)" fi,sqrt(temp),j", fi,temp**0.5,j
c                endif
c---  ) debug
                
                total_num = 0.0D0 ! total number of ptcl.
                total_num_shells = 0.0D0 ! total number of ptcl in the shells.
                total_num_esc_shell = 0.0D0 ! total number of escaped ptcl.
                
                Do l = 1, escape_shell_num(j) ! only Re-distribute particles to those shells 
                                !  that have already had particles of this energy.
                   
c     Do l = 1, no_shells  ! only Re-distribute particles to those shells 
                   
                   call diffuse_distr_integ(num_ptcls,fi,shell_posn, ! get the number of particles 
     *                  R_shell_x(l), R_shell_x(l+1), temp) !in the lth shell. 
                   
c     write(22,*)"---num_ptcls from subroutine", num_ptcls
                   num_ptcls = get_number(fi,shell_posn, ! get the number of particles 
     *                  R_shell_x(l), R_shell_x(l+1), temp)
c     write(22,*)"---num_ptcls from code", num_ptcls
                   
ccc   print*,"fi,num_ptcls,temp,j",fi,num_ptcls,temp,j
ccc   print*,"l, shell_posn,R_shell_x(l,l+1)",l,
ccc   *                  shell_posn,R_shell_x(l), R_shell_x(l+1)
                   
c---  (- debug
c     if ( j .ge. 150) then
c     write(96,*)" fi,num_ptcls,temp,j",
c     *                     fi,num_ptcls,temp,j
c     write(96,*)" l,shell_posn,R_shell_x(l,l+1)",l,
c     *                  shell_posn,R_shell_x(l), R_shell_x(l+1)
c     endif
c---  -) debug 
                   
c------num_ptcls is the output: ptcl number (per unit momentum space)
c------between R_shell_x(l) and R_shell_x(l+1).
c------fi: input ptcl number in the shell (per unit momentum space) designated by shell_posn.
c------shell_posn: the shell in consideration.
c------R_shell_x(l): beginning loc. of  lth shell. 
c------R_shell_x(l): ending loc. of  lth shell.                              
c------temp: the length scale associated with p.
                   
                   diffuse_array(l,p_int) = diffuse_array(l,p_int) +
     *                  num_ptcls
                   
                   total_num = total_num + num_ptcls
                   total_num_shells = total_num_shells + num_ptcls 
c     
                End Do          ! (l = 1, escape_shell_num(j) )

c       print*, " ----i, j,total_num,fi----", i,j,total_num,fi
c     print*, "fi, total_num ", fi, total_num
                
c---  Now calculate kappa_norm0 which is used to calculate esc_length for jth momentum.
                
                kappa_norm0 = kappa(j, escape_shell_num(j) )
     *               *2858055.0D0/1.5D11/1.5D11 

c---  
c---  Calculate esc_length for this momentum.
c---  The index of the outmost shell for this momentum is 
c---  given by escape_shell_num(j).
c---  Then get the value for escape_location(p_int)
                


                if ( j .le. ipmax(escape_shell_num(j) ) ) then
                   esc_length = 2.0D0 * kappa_norm0/
     *                  (shock_speed-v1(shock_end+2,1,5)) ! Note e^4 = 55
                else
                   esc_length = 0.
                endif

                esc_length = 0.
               
c--------- adding the following by GL -------
                if ( escape_shell_num(j)  .eq. no_shells ) then ! we are in the first shell.
                   esc_length = 2.0D0 * kappa_norm0/
     *                  (shock_speed-v1(shock_end+2,1,5)) ! Note e^4 = 55
                else if ( escape_shell_num(j)  .lt. no_shells ) then
                   esc_length = 0.
                     endif
                endif
c---------end adding the following by GL --------                 

                escape_location(p_int)=R_shell_x(escape_shell_num(j)+1)
     *               + esc_length


c------ After testing the above way for escape_location, I decided
c------ to use the following value for escape_location(p_int).
c------ Basically, this is the location of the outer boundary of the shell
c------ for particles with momentum p_int. I expect this to be a good approx.
c------ GL, 02/20/2004.
c                escape_location(p_int)=R_shell_x(escape_shell_num(j)+1) 

                
c
c                print*, "R_shell_x(escape_shell_num(j)+1) ",
c     *               R_shell_x(escape_shell_num(j)+1) 

c                print*, "i,j,escape_shell_num(j)+1, esc_length",i,j,
c     *               escape_shell_num(j)+1,esc_length
                
c                print*, "j,esc_length,escape_location(p_int)",
c     *               p_int,esc_length, escape_location(p_int)

                write(85,*)"i,j,escape_shell_num(j)+1, esc_length",i,j,
     *               escape_shell_num(j)+1,esc_length
                
                write(85,*)"j,esc_length,escape_location(p_int)",
     *               p_int,esc_length, escape_location(p_int)


c----- Now let us calculate escaped particles.
                
                temp1 = escape_location(p_int)
                temp2 = min(x1a(ie-10), 3.5 * temp1)
c     temp2 = 10. * temp1
                
                call diffuse_distr_integ(num_ptcls, fi, shell_posn,
     *               temp1, temp2, temp) ! temp is the (length scale)^2 for the diffusion.
                
ccc   print*, " num_ptcls using diffuse_distr", num_ptcls
                
                num_ptcls = get_number(fi, shell_posn,
     *               temp1, temp2, temp) ! temp is the (length scale)^2 for the diffusion.
                
ccc   print*, " num_ptcls using get_number", num_ptcls
                
c     pause
                
ccc   print*, "esc_num,temp1,temp2",num_ptcls,temp1,temp2
ccc   print*, "fi,shell_posn,temp",fi,shell_posn,temp
                
ccc   if (j .ge. 150) then
ccc   write(96,*)"esc_#,temp1,temp2",num_ptcls,temp1,temp2
ccc   endif

                esc_number_up(p_int) = esc_number_up(p_int) + 
     *               num_ptcls 
                
                esc_ptcl_counter(p_int) = esc_ptcl_counter(p_int) + 
     *               num_ptcls 
                
                total_num = total_num + num_ptcls
                total_num_shells = total_num_shells + num_ptcls 

c--- let us put particles between R_shell_x(escape_shell_num(p_int)+1) 
c--- and  escape_location(p_int) to diffuse_array(l,p_int).
                
                temp1 = R_shell_x(escape_shell_num(p_int)+1) 
                temp2 = escape_location(p_int)
                
                if (temp1 .ne. temp2) then 
                   l = escape_shell_num(j)
                   call diffuse_distr_integ(num_ptcls,fi,shell_posn, ! get the number of particles 
     *                  temp1,temp2 , temp) !in the lth shell. 
                   
                   num_ptcls= get_number(fi,shell_posn, ! get the number of particles 
     *                  temp1,temp2 , temp) !in the lth shell. 
                   
                   diffuse_array(l,p_int) = diffuse_array(l,p_int) +
     *                  num_ptcls

c---- ( debug
                   if ( j .ge. 150) then
                      write(96,*) "-- ptcls in (outlayer, escape bc) --"
                      write(96,*)"num_ptcls,temp1,temp2",
     *                     num_ptcls,temp1,temp2 
                   endif
                   
ccc                   print*, " ----- particles between the outlayer and the
ccc     *                  the escape boundary ---- "
ccc                   print*, "num_ptcls", num_ptcls
ccc                   print*, " temp1,temp2", temp1,temp2
             
c---- ) debug      
                   total_num = total_num + num_ptcls
                   total_num_shells = total_num_shells + num_ptcls 
                   
                endif
                
c---- The following calculate particles leaking from the shock rear end. 
                
                k = R_shell_i(1) ! index of the location of the inner bd. of the first shell.                                
                temp1 = x1a(is+10) 
                temp2 = R_shell_x(1)
                
                if ( (temp2 - temp1) .lt. 3.* temp**0.5 ) then
                   temp1 =  0.01
                endif

                call diffuse_distr_integ(num_ptcls, fi, shell_posn,
     *               temp1, temp2, temp) 
                
                num_ptcls=get_number(fi, shell_posn,
     *               temp1, temp2, temp) 
                
                if ( j .ge. 150) then
                   write(96,*)"shell_posn,x1a(is+10),R_shell_x(1)",
     *                  shell_posn,x1a(is+10),R_shell_x(1)
                   write(96,*)"fi, num_ptcls, temp1,temp2",
     *                  fi,num_ptcls,temp1,temp2
                endif
                
c---- Behind the shock, the CME ejecta can reflect energetic particles back.
                
                Do l = 1, escape_shell_num(p_int) 
                   diffuse_array(l,p_int) = diffuse_array(l,p_int) +
     *                  num_ptcls * (R_shell_x(l+1)**3- R_shell_x(l)**3)
     *                  /  (R_shell_x(escape_shell_num(p_int)+1)**3
     *                  - R_shell_x(1)**3)
                Enddo

c---- The following for back propagation. 

                esc_number_down(p_int) = esc_number_down(p_int) 
     *               + num_ptcls

                total_num = total_num + num_ptcls 
                total_num_esc_shell = total_num_esc_shell + num_ptcls

                write(96,*) "j, fi, total_num, esc, esc/tot,time", 
     *               j, fi, total_num, esc_number_up(p_int), 
     *               esc_number_up(p_int)/total_num, time
                write(96,*) " **** end  diff ****"

c                print*, " p_int, escape_shell_num(p_int)", 
c     *               p_int, escape_shell_num(p_int)
c                pause

c                print*,"total_num(1+2)-distr_array(i,j),total_num,fi", 
c     *               total_num-distr_array(i,j),total_num,fi

c---   GL April 22nd, 2003.
                

c---  Now calculate the number density at the escape boundary @ shock front.
                call diff_distr_fn_posn(dens_ptcls,fi,shell_posn, ! get the density from the     
     *               x1a(shock_begin)+esc_length, temp) ! particle number fi.
c---  try the following.
                
                call diff_distr_fn_posn(dens_ptcls,fi,shell_posn, ! get the density from the     
     *               x1a(shock_begin), temp) ! particle number fi.
                
                esc_dens_up(p_int) = esc_dens_up(p_int) + 
     *               dens_ptcls

c---  Now calculate the number density at the inner escape boundary at the 1st shell.
                call diff_distr_fn_posn(dens_ptcls,fi,shell_posn, ! get the density from the     
     *               R_shell_x(1), temp) ! particle number fi.
                
                esc_dens_down(p_int) = esc_dens_down(p_int) + 
     *               dens_ptcls
                
                do ii = 1, no_shells
                   call diff_distr_fn_posn(dens_ptcls,fi,shell_posn, ! get the density from the     
     *                  R_shell_x(ii+1), temp) ! particle number fi.
                   
                   esc_dens_shell(p_int, ii) = esc_dens_shell(p_int, ii) 
     *                  + dens_ptcls
                enddo
                
             End Do             ! if (fi .NE. 0.) 
             
c     Write(46,*) " ***", Energy, esc_dens_up(j), 
c     *          diffuse_array(i,j)/shell_vol
             
          End Do                ! do j = 1, p_number
       End Do                   ! Do i = 1, no_shells
      

 
       
c---  ( debug

c     do j =1 ,250
c          print*, "j,escape_location(j)",j,escape_location(j)
c     enddo
c     pause

c---  ) debug
       
       Do i = 1, no_shells +1
          Do j = 1, p_number    !  Reset ptcl distributions as they
                                !  have now diffused and need to be included 
             
c---  ( debug
c     print*, "---distr_array (i,j) before diffuse ---", i,j,
c     *            distr_array(i,j)
c---  ) debug
             
             distr_array(i,j) = diffuse_array(i,j) ! in the distr array.
             
c---  ( debug
c     print*, "---distr_array (i,j) after diffuse ---", i,j,
c     *         distr_array(i,j)
c---  ) debug
             
          End Do
       End Do  

c--------------
       
       Do j = 1, p_number
          temp_num(j) = 0.
       ENDDO
       temp_dV = 0.
       
       DO i =1, no_shells 
          DO j =1, p_number
             p = p0 + del_p * (j-1)

             factor_num = 4.*pi/3.* (Exp(p)**3 - Exp(p-del_p)**3) 
             dV = 4.*pi/3. *(  R_shell_x(i+1)**3 
     *         - R_shell_x(i)**3 ) ! no need to convert R_shell_x.
             
             temp_num(j) = temp_num(j) + distr_array(i,j)
             temp_dV = temp_dV + dV
          ENDDO
       ENDDO
       
       Do j = 1, p_number
          p = p0 + del_p * (j-1)
          Energy = Sqrt((1.67E-27*3.0E8**2.0)**2.0 + 
     *         (1.67E-27*Exp(p)*52483.25*3.0E8)**2.0) -
     *         1.67E-27*3.0E8**2.0 ! the kinetic energy.
          Energy = Energy/1.6E-19 !  (Energ
          ttt = temp_num(j) + esc_number_up(j)
          tttt = temp_num(j) + esc_ptcl_counter(j)

          write(82,*) Energy,temp_num(j), ttt, tttt,
     *         temp_num(j)/temp_dV
       Enddo

c-------------


c      pause
c-----------------Now for heavy ions-------------------------------------------

       del_p = (max_p - p0)/(p_number *1.) ! step of log(p) (p is in measure of p/mp) 
c     
       Do i = 1, no_shells + 1
          Do j = 1, p_number
             diffuse_array(i,j) = 0.0D0 
             esc_dens_up_hi(j) = 0.0D0 
             esc_dens_down_hi(j) = 0.0D0 
          End Do
       End Do                   
c     
c--- ( del
       Do i = 1, in             ! # of position.
          Do j = 1, p_number         ! # of momentum.
             esc_shell_hi(i, j) = 0.0D0 
          End Do           
       End Do    
c--- ) del   


       Do i = 1, no_shells   ! For all the shells generated earlier.
          
          shell_posn = (R_shell_x(i + 1) + R_shell_x(i))/2.0D0 ! the loc. of the center of the ith shell.
          shell_posn_int = (R_shell_i(i+1)+R_shell_i(i))/2.0D0 ! the index of the center of the ith shell.
c     
          shell_vol = 4.0*pi*(R_shell_x(i+1)**3.0D0-
     *         R_shell_x(i)**3.0D0)/3.0D0 ! volume of the shell.    
c     

          Do j = 1, p_number
             p_int = j
c             print*, " p_int and j are", p_int, j
             p = p0 + (j-1) * del_p

             Energy = Sqrt((1.67E-27*3.0E8**2.0)**2.0 + 
     *            (1.67E-27*Exp(p)*52483.25*3.0E8)**2.0) -
     *            1.67E-27*3.0E8**2.0 ! the kinetic energy.
             Energy = Energy/1.6E-19 !  (Energy in eV)
c     
             vel = Exp(p)*52483.25/Sqrt(1.0 + Exp(p)*Exp(p)*
     *            52483.25*52483.25/(3.0E8*3.0E8)) ! particle velocity.
                                ! note v/c = (p/mc) / sqrt( 1 + (p/mc)^2 ) is in m/s.
             
             vel = vel/52483.25D0 ! Normalized velocity (in unit of the sound speed at 1 AU).
             
             fi = distr_array_hi(i, p_int) ! number of particles before the diffusion.

c===========add by Gang Li for a better approximation for escape 

C            /added by Olga 
c	       if (no_shells .gt. 22) then
c                   delta_out =0.
c	       endif
C            /added by Olga  


             if ( i .eq. no_shells) then 
                temp_esc_hi(p_int) =  fi * 1. / (1.+delta_out)
                fi  = (1. - 1./(1.+delta_out) ) * fi
             endif


c===========  end add by Gang Li



c---     
             kappa_norm = kappa_hi(j,i)*2858055.0D0/1.5D11/1.5D11 ! get the normalised kappa.
            
             temp0 = 4.0D0*kappa_norm*(time - prev_acc_time) ! the length scale squre 4 the 
                                                            ! Green's function.

c             temp = temp0 + (kappa_norm/(shock_speed - 
c     *            v1(shock_end+2,1,5)) )**2.0D0 ! take into account the initial extended 
c                                !   of the assumed delta function ditribution. 
             temp = temp0

             temp = temp0 + ((R_shell_x(i + 1) - R_shell_x(i))/5.)**2


c     temp_new = vel*(time-prev_acc_time)  ! the distence particle can move ballistically.
c     If (temp .gt. temp_new) Then ! Ensure ptcl does not 
c     print*, "--- temp, temp2---", temp, temp_new
c     temp = temp_new ! diffuse further than it
c     ! could propagate ballistically. 
c     EndIf              
             
             If (fi .NE. 0.0D0) Then ! so this momentum bin contains particles.
                
                total_num = 0.0D0 ! total number of ptcl.
                total_num_shells = 0.0D0 ! total number of ptcl in the shells.
                total_num_esc_shell = 0.0D0 ! total number of escaped ptcl.
                
                Do l = 1, no_shells 
                   
                   call diffuse_distr_integ(num_ptcls,fi,shell_posn, ! get the number of particles 
     *                  R_shell_x(l), R_shell_x(l+1), temp) !in the lth shell. 
                   
                   num_ptcls = get_number(fi,shell_posn, ! get the number of particles 
     *                  R_shell_x(l), R_shell_x(l+1), temp) !in the lth shell. 

c------num_ptcls is the output: ptcl number between R_shell_x(l) and R_shell_x(l+1).
c------fi: input ptcl number in the shell designated by shell_posn.
c------shell_posn: the shell in consideration.
c------R_shell_x(l): beginning loc. of  lth shell. 
c------R_shell_x(l): ending loc. of  lth shell.                              
c------temp: the length scale associated with p.
                   
                   diffuse_array(l,p_int) = diffuse_array(l,p_int)+
     *                  num_ptcls
                   
                   total_num = total_num + num_ptcls
                   total_num_shells = total_num_shells + num_ptcls 
c     
                End Do          ! (l = 1, no_shells)

c                print*, "total_num  - distr_array_hi(i,j)", 
c     *               total_num - distr_array_hi(i,j), total_num
c     
c--- Now calculate kappa_norm0 which is used to calculate esc_length for jth momentum.
                kappa_norm0 = kappa_hi(j,escape_shell_num_hi(j))
     *               *2858055.0D0/1.5D11/1.5D11 ! used to calculate esc_length

                if ( j .lt. ipmax_hi(no_shells) ) then
                   esc_length = 1.0D0 * kappa_norm0/
     *                  (shock_speed-v1(shock_end+2,1,5)) ! Note e^4 = 55
                else
                   esc_length = 0.
                endif

                esc_length = 0.

c--------- adding the following by GL  06/14/2007 -------
                if ( escape_shell_num_hi(j)  .eq. no_shells ) then ! we are in the first shell.
                   esc_length = 2.0D0 * kappa_norm0/
     *                  (shock_speed-v1(shock_end+2,1,5)) ! Note e^4 = 55
                else if ( escape_shell_num_hi(j)  .lt. no_shells ) then
                   esc_length = 0.
                     endif
                endif
c---------end adding the following by GL --------                 

                escape_location_hi(p_int)=
     *               R_shell_x(escape_shell_num_hi(p_int)+1)+esc_length

c------ After testing the above way for escape_location, I decided
c------ to use the following value for escape_location(p_int).
c------ Basically, this is the location of the outer boundary of the shell
c------ for particles with momentum p_int. I expect this to be a good approx.
c------ GL, 02/20/2004.
c                escape_location_hi(p_int)=
c     *               R_shell_x(escape_shell_num_hi(j)+1) 
                
                
cc     If ((DExp(p) .GE. p_hi_inj) .AND. ! only particles within the acceleration range has 
cc     *               (DExp(p) .LE. p_hi_cut)) Then ! non-zero esc_length.
                
cc     esc_length = 4.0D0*kappa_norm0/
cc     *                  (shock_speed-v1(shock_end+2,1,5))
cc     Else
cc     esc_length = 0.0D0
cc     EndIf  
cc     escape_location_hi(p_int) = x1a(shock_begin)+ esc_length
                
c-----Now let us calculate escaped particles.
                temp1 = escape_location_hi(p_int)
                temp2 = min(x1a(ie-10), 3.5 * temp1)
                
                call diffuse_distr_integ(num_ptcls, fi, shell_posn,
     *               temp1, temp2, temp) ! temp is the (length scale)^2 for the diffusion.
                
                num_ptcls =get_number(fi, shell_posn,
     *               temp1, temp2, temp) ! temp is the (length scale)^2 for the diffusion.

                esc_number_up_hi(p_int) = esc_number_up_hi(p_int) 
     *               + num_ptcls
                
                esc_ptcl_counter_hi(p_int) = esc_ptcl_counter_hi(p_int) 
     *               + num_ptcls
                
c---  let us put particles between R_shell_x(escape_shell_num(p_int)+1) 
c---  and  escape_location(p_int) back to diffuse_array(l,p_int)
                
                temp1 = R_shell_x(escape_shell_num_hi(p_int)+1) 
                temp2 = escape_location_hi(p_int)
                
                if (temp1 .ne. temp2) then 
                   l = escape_shell_num(p_int)
                   call diffuse_distr_integ(num_ptcls,fi,shell_posn, ! get the number of particles 
     *                  temp1,temp2 , temp) !in the lth shell. 

                   num_ptcls=get_number(fi,shell_posn, ! get the number of particles 
     *                  temp1,temp2 , temp) !in the lth shell. 

                   diffuse_array(l,p_int) = diffuse_array(l,p_int) +
     *                  num_ptcls
                endif

                total_num = total_num + num_ptcls 
                total_num_esc_shell = total_num_esc_shell + num_ptcls
                
c     diffuse_array(no_shells+1, p_int) =  num_ptcls + 
c     *               diffuse_array(no_shells+1, p_int)
                
c---  GL April 22nd, 2003.
                
c---- Now calculate the number density at the escape boundary @ shock front.   
                call diff_distr_fn_posn(dens_ptcls,fi,shell_posn,
     *               x1a(shock_begin)+esc_length, temp)             
                
                esc_dens_up_hi(p_int) = esc_dens_up_hi(p_int) + 
     *               dens_ptcls
                
c---  Now calculate the number density at the inner escape boundary at the 1st shell.
                call diff_distr_fn_posn(dens_ptcls,fi,shell_posn, ! get the density from the     
     *               R_shell_x(1), temp) ! particle number fi.
                
                esc_dens_down_hi(p_int) = esc_dens_down_hi(p_int) + 
     *               dens_ptcls
  
                do ii = 1, no_shells
                   call diff_distr_fn_posn(dens_ptcls,fi,shell_posn, ! get the density from the     
     *                  R_shell_x(ii), temp) ! particle number fi.
                   
                   esc_dens_shell_hi(p_int, ii) = 
     *                  esc_dens_shell_hi(p_int, ii) + dens_ptcls
                enddo

c------The following calculate the number of ptcl in the region $ R(t) < r < R(t) + \lambda_{esc}$.
                
c                k = R_shell_i(no_shells + 1) ! index of the location of the shock.
c                
c                call diffuse_distr_integ(num_ptcls, fi, shell_posn,
c     *               R_shell_x(no_shells+1),
c     *               R_shell_x(no_shells+1)+esc_length,temp)

c                esc_shell_hi(k, p_int)=esc_shell_hi(k,p_int)+num_ptcls
                
c                total_num = total_num + num_ptcls 
c                total_num_esc_shell = total_num_esc_shell + num_ptcls
c---------------------------

                k = R_shell_i(1) ! index of the location of the inner bd. of the first shell.                
                temp1 = x1a(is+10) 
                temp2 = R_shell_x(1)

                if ( (temp2 - temp1) .lt. 3.* temp**0.5 ) then
                   temp1 =  0.01
                endif


                call diffuse_distr_integ(num_ptcls, fi, shell_posn,
     *               temp1, temp2, temp) 

                num_ptcls =get_number(fi, shell_posn,
     *               temp1, temp2, temp) 
                

c---- Behind the shock, the CME ejecta can reflect energetic particles back.

                Do l = 1, escape_shell_num_hi(p_int) 
                   diffuse_array(l,p_int) = diffuse_array(l,p_int) +
     *                  num_ptcls * (R_shell_x(l+1)**3- R_shell_x(l)**3)
     *                  /  (R_shell_x(escape_shell_num(p_int)+1)**3
     *                  - R_shell_x(1)**3)
                Enddo

c---- The following for back propagation. 
                esc_number_down_hi(p_int) = esc_number_down_hi(p_int) 
     *               + num_ptcls
c----
                total_num = total_num + num_ptcls 
                total_num_esc_shell = total_num_esc_shell + num_ptcls
                
                if (abs(total_num-distr_array_hi(i,j))
     *               /distr_array_hi(i,j)  .ge. 5.D-5) then
                   
c                   print*,"i,j,total_num-distr_array_hi(i,j),distr_arr",
c     *                  i,j,total_num-distr_array_hi(i,j), 
c     *                  distr_array_hi(i,j)                   
c                   pause 
                endif
             Enddo             ! if (fi .NE. 0.) 
                
c             Write(66,*) Energy, esc_dens_up_hi(j), 
c     *            diffuse_array(i,j)/shell_vol
             
          End Do                ! do j = 1, p_number
       End Do                   ! Do i = 1, no_shells - 1 
       

       Do i = 1, no_shells + 1
          Do j = 1, p_number         !  Reset ptcl distributions as they
                                !  have now diffused and need to be included 

c---( debug
c          print*, "---distr_array_hi (i,j) before diffuse ---", i,j,
c     *            distr_array_hi(i,j)
c---  ) debug

             distr_array_hi(i,j) = diffuse_array(i,j) ! in the distr array.

c---( debug
c          print*, "---distr_array_hi (i,j) after diffuse ---", i,j,
c     *            distr_array_hi(i,j)
c---  ) debug

          End Do
       End Do  
c---
       write(13,*)"--- time is ---", time, r_shell_x(no_shells+1)
       write(14,*)"--- time is ---", time, r_shell_x(no_shells+1)




       temp_dV = 0.

       Do j = 1, p_number
          temp_num(j) = 0.          
          temp_num_hi(j) = 0.          
          temp_dV = 0.
       Enddo


       DO i =1, no_shells 
c          write(13,*)"--- shell no. ---", i
c          write(14,*)"--- shell no. ---", i
          DO j =1, p_number
             p = p0 + del_p * (j-1)
             Energy = Sqrt((1.67E-27*3.0E8**2.0)**2.0 + 
     *            (1.67E-27*Exp(p)*52483.25*3.0E8)**2.0) -
     *            1.67E-27*3.0E8**2.0 ! the kinetic energy.
             Energy = Energy/1.6E-19 !  (Energy in eV)
             
             factor_num = 4.*pi/3.* (Exp(p)**3 - Exp(p-del_p)**3) 
             dV = 4.*pi/3. *(  R_shell_x(i+1)**3 
     *         - R_shell_x(i)**3 ) ! no need to convert R_shell_x.

             temp_num(j) = temp_num(j) + distr_array(i,j)
             temp_num_hi(j) = temp_num_hi(j) + distr_array_hi(i,j)
             temp_dV = temp_dV + dV

             write(13,*) Energy,distr_array(i,j)/dV,  !distr_array(i,j)/dV * factor_num,
     *            distr_array_hi(i,j)/dV, distr_array(i,j),
     *            distr_array_hi(i,j)

             write(14,*) Energy,esc_dens_shell(j,i),
     *            esc_dens_shell_hi(j,i)

c             write(14,*) Energy,distr_array(i,j)/dV,
c     *            distr_array_hi(i,j)/dV

                                !write distr_array(no_shells, j) to file (14).
          ENDDO
       ENDDO

       Do j = 1, p_number
          p = p0 + del_p * (j-1)
          Energy = Sqrt((1.67E-27*3.0E8**2.0)**2.0 + 
     *         (1.67E-27*Exp(p)*52483.25*3.0E8)**2.0) -
     *         1.67E-27*3.0E8**2.0 ! the kinetic energy.
          Energy = Energy/1.6E-19 !  (Energy
          write(13,*) Energy,temp_num(j)/temp_dV,
     *         temp_num_hi(j)/temp_dV, temp_num(j), temp_num_hi(j)
       Enddo


       
c------ what does the following do?

c       if ( (index1 .gt. index2) .and. (index2 .gt. index3) ) then

c          test_alpha = log(esc_number_up(index3)/esc_number_up(index2))
c     *         / log(p_cut_expand / p_cut) 

c          test_fac=(p_cut**(1-test_alpha)-p_cut_expand**(1-test_alpha))
c     *         /( p_cut **(1-test_alpha) - p_cut_old**(1-test_alpha) )

c          do j = index3, index2
c             esc_number_up(j) =  test_fac * esc_number_up(j)
c          enddo

c          do j = index2+1, index1
c             test_p = exp( ( (j-1)*del_p + p0) )
c             esc_number_up(j) =  test_fac * esc_number_up(index2)
c     *            * (test_p/p_cut_expand) **(-test_alpha)
c          enddo

c       endif

c       p_cut_old = p_cut
c       p_hi_cut_old = p_hi_cut

c       do j = 1, p_number
c          print*, "j,esc_number_up(j)",j,esc_number_up(j)
c       enddo

c---
c--- Now calculate the current at the shock front and
c--- write out to file (78). Also write out to file (96)
c--- the distr_array(no_shells, p_int) and distr_array_hi(no_shells, p_int)
c--- Note, distr_array (no_shells,p_int) has unit of number/ (p^3).
c---
 
       do j = 1, p_number 
          p_int = j
          p = p0 + (j - 1) *del_p

          fi = distr_array(no_shells, p_int)
          
          fi_hi = distr_array_hi(no_shells, p_int)

          dV = 4.*pi/3. *(  R_shell_x(no_shells+1)**3 
     *         - R_shell_x(no_shells)**3 ) ! no need to convert R_shell_x.
          
          den_sh = fi / dV * 4./3.*pi*(Exp(p)**3-Exp(p-del_p)**3) 
          den_sh_hi = fi_hi/dV * 4./3.*pi*(Exp(p)**3-Exp(p-del_p)**3) 
          
c---- Now calculate the current f * v_sh / dE where dE is the energy gap. 

          energy_low = Sqrt((1.67E-27*3.0E8**2.0)**2.0 + ! particle kinetic energy/nucleon.
     *         (1.67E-27*Exp(p-del_p)*52483.25*3.0E8)**2.0) -
     *         1.67E-27*3.0E8**2.0
          energy_low = energy_low/1.6D-19 ! Convert to eV.
          
          energy = Sqrt((1.67E-27*3.0E8**2.0)**2.0 + ! particle kinetic energy/nucleon.
     *         (1.67E-27*Exp(p)*52483.25*3.0E8)**2.0) -
     *         1.67E-27*3.0E8**2.0
          energy = energy / 1.6D-19
          energy_hi = energy

          Vpar = (Exp(p)* 52483.25)* 938. / (energy/1E6 +938 )  ! particle speed

          Vsh = shock_speed * vo  ! shock velocity in m/s.
          
c          print*, " ---- Vpar, Vsh -----", Exp(p)*52483.25, Vpar, Vsh

c          den_sh_j = den_sh * Vsh/(energy_hi - energy_low) 
          den_sh_j = den_sh * Vpar/(energy_hi - energy_low) 
          den_sh_j = den_sh_j *1.D6 * 1.D2  ! to 1/cm^3 (cm/s) * (1/MeV)

c          den_sh_hi_j = den_sh_hi*Vsh / (energy_hi - energy_low)
          den_sh_hi_j = den_sh_hi*Vpar / (energy_hi - energy_low)
          den_sh_hi_j = den_sh_hi_j *1.D6 * 1.D2  ! to 1/cm^3 (cm/s) * (1/MeV)
          
                                ! den_sh_j is the current
                                ! den_sh_hi_j is the current for hvy ion.

          write(78,*) Energy, den_sh_j, den_sh_hi_j,
     *         x1a(shock_begin), time

          write(79,*) Energy, distr_array(no_shells, p_int), 
     *         distr_array_hi(no_shells, p_int), time
       enddo

c       pause


c====== added by Gang Li, Feb. 2, 2007 

       do j = 1, p_number
          esc_number_up(j) = esc_number_up(j) + temp_esc(j) 
          esc_number_up_hi(j) = esc_number_up_hi(j) + temp_esc_hi(j) 
       enddo

c====== end added by Gang Li, Jan. 

c---  Now write out esc_number_up(p_int), x1a(shock_begin)+ esc_length,
c---  and time to file (88).

       do j = 1, p_number 
          p_int = j
          p = p0 + (j - 1) *del_p
          energy = Sqrt((1.67E-27*3.0E8**2.0)**2.0 + ! particle kinetic energy/nucleon.
     *         (1.67E-27*Exp(p)*52483.25*3.0E8)**2.0) -
     *         1.67E-27*3.0E8**2.0
          energy = energy / 1.6D-19

c ----- Need to get the real number,rather than number/d^3p.
          temp_esc_up = esc_number_up(p_int) *  ! real particle number that leaked out.
     *         4 * pi /3. * (Exp(p)**3 - Exp(p-del_p)**3)

          temp_esc_up_hi = esc_number_up_hi(p_int) *  ! same as temp_esc_up, for heavyIon.
     *         4. * pi /3. * (Exp(p)**3 - Exp(p-del_p)**3)

          temp_counter1  = esc_ptcl_counter(p_int) * ! accumulated particle number that leaked out.
     *         4. * pi /3. * (Exp(p)**3 - Exp(p-del_p)**3)

          temp_counter2  = esc_ptcl_counter_hi(p_int) * !same as temp_counter1, for HeavyIon.
     *         4. * pi /3. * (Exp(p)**3 - Exp(p-del_p)**3)

c ----- Should times (1.5*13)**3 for the real number. No need for that.
 
          If (Exp(p) .LE. p_cut ) Then

             write(88,*) Energy, 
     *            temp_esc_up, temp_counter1,
     *            escape_location(p_int),time

          ELSE
             write(88,*) Energy, 
     *            temp_esc_up, temp_counter1,
     *            escape_location(p_int),time
          ENDIF

          write(38,*) Energy, esc_dens_up(p_int), 
     *         escape_location(p_int), time


          
          If (Exp(p) .GE. p_hi_cut ) Then
             write(29,*) Energy, esc_number_up_hi(p_int), 
     *            esc_ptcl_counter_hi(p_int),
     *            escape_location_hi(p_int),time

             write(89,*) Energy, 
     *            temp_esc_up_hi, temp_counter2,
     *            escape_location(p_int),time

          ELSE

             write(29,*) Energy, esc_number_up_hi(p_int), 
     *            esc_ptcl_counter_hi(p_int),
     *            escape_location_hi(p_int),time

             write(89,*) Energy, 
     *            temp_esc_up_hi, temp_counter2,
     *            escape_location(p_int),time
          ENDIF
          
          write(58,*) Energy, esc_dens_up_hi(p_int), 
     *         escape_location_hi(p_int), time
       enddo

c----(  move data-writing to file (68) from expand_distr.f to here. 
c---- Now let us calculate the total number of each shell after the expansion.

       totnum = 0.0
c     
       Do i = 1, no_shells

          Do j = 1, 250
             p = 2.0 + (j-1) * del_p
             totnum = totnum + 4.* pi/3. * distr_array(i,j)*
     *            ( Exp(p+del_p)**3 - Exp(p)**3 )
          End Do
          total_num_p(i) = totnum 
          
          totnum = 0.0
          Do j = 1, 250
             p = 2.0 + (j-1) * del_p
             totnum = totnum + 4.* pi/3. * distr_array_hi(i,j)*
     *            ( Exp(p+del_p)**3 - Exp(p)**3 )
          End Do
          total_num_hi(i) = totnum 
       ENDDO
       tempR= R_shell_x(no_shells+1)

       write(68,*) "--- time and R_shock ---", time, tempR
       write(69,*) "--- time and R_shock ---", time, tempR

       do i = 1, no_shells 

          dV = 4.*pi/3. *(  R_shell_x(i+1)**3 
     *         - R_shell_x(i)**3 )

          write(68,*) i, R_shell_i(i), R_shell_x(i),total_num_p(i)/dV,
     *         total_num_hi(i)/dV
          write(69,*) i, R_shell_i(i), R_shell_x(i),total_num_p(i),
     *         total_num_hi(i)

       enddo

c---- ) 

c       pause
c---  ( debugging.

       print*, "---- after diffusion ---"
ccc       do i = 1, no_shells  
ccc          print*, " -- total_num_p(i) --", total_num_p(i)
ccc          print*, " -- total_num_hi(i) --", total_num_hi(i)
ccc       enddo

ccc       print*, "---- for debuging purpose -----"
ccc       print*, "---- print out from 100 to 110, the following ---"
ccc       print*, "---- pint, energy,esc_number_up[hi](p_int),x,time" 
ccc       do p_int = 100, 110
ccc          write(*,*)p_int,esc_number_up(p_int),esc_number_up_hi(p_int), 
ccc     *         escape_location(p_int), time
ccc       enddo

c---  ) debugging.
c       pause

       print*, "******* Leaving diffuse_distr.f *******"       
       print*, "                                       "
c       pause
       return
       End


      Function get_number(fo, ri, r1, r2,temp)
      real*8 get_number
      real*8 num_ptcls, fo, ri, r1, r2, temp, pi    
      real*8 rsml,rbig, alpha_temp, xx1,xx2,xx3,test2
      
      pi = 3.14159265
      alpha_temp = temp**0.5
      
      if (r1 .lt. r2) then
         rsml = r1
         rbig = r2
      else
         rsml = r2
         rbig = r1
      endif


      if ( (rsml .lt. ri) .and. (ri .lt. rbig) ) then
         xx1 = abs( (rsml - ri) / alpha_temp )
         xx2 = abs( (rbig - ri) / alpha_temp )
         test2 = alpha_temp/ri / sqrt(4 *pi) * fo
         test2 = test2 * ( Dexp(-xx1**2) - Dexp(-xx2**2)  +
     *        sqrt(pi) * ri/alpha_temp *(erf(xx2) + erf(xx1) ))
         num_ptcls = test2
      else
         
         xx1 =  abs(rsml - ri) / alpha_temp 
         xx2 =  abs(rbig - ri) / alpha_temp 
         
         if ((Derf(xx2) .eq. 1.) .and. (Derf(xx1) .eq. 1.)) then
            test2 = alpha_temp/ri / sqrt(4 *pi) * fo
            test2 = test2 * (Dexp(-xx1**2) - Dexp(-xx2**2) +
     *           sqrt(pi) * ri/alpha_temp *abs(Derfc(xx1)-Derfc(xx2)))
            num_ptcls = test2
            
         else
            test2 = alpha_temp/ri / sqrt(4 *pi) * fo
            test2 = test2 * (Dexp(-xx1**2) - Dexp(-xx2**2) +
     *           sqrt(pi) * ri/alpha_temp *abs(Derf(xx2) - Derf(xx1)))
            num_ptcls = test2
         endif

      endif

c     endif

c      if (num_ptcls .eq. 0.) then
c         print*, " --- xx1,xx2, alpha_temp ---",xx1,xx2, alpha_temp 
c         print*, " --- Derfc(xx1),Derfc(xx2) ---",Derfc(xx1),Derfc(xx2)
c      endif

c     pause

      get_number = num_ptcls
      return
      
      End


      
