c======================================================================|
c     solar wind model (6-ref.long. calculation)
c     last update 2015.7.23
c======================================================================|
      implicit none
      integer :: i,j,iang,nang
      integer :: idp_in,idp_prop,idp_out
      real*8 :: angref0,touts,dtr
      real*8 :: angref(30)
      integer :: instop,instopfin,ifnpo,idprop
      character*100 :: fnin,fnout,fnin1,fnout1,fntmp(12),fntmp2(12)
      character*100 :: fdirtmp
      real*8,allocatable :: angref6(:)
      real*8 :: xmin1,xmax1,xear
      integer :: ix
!!!!! elena : namelist should be before executable statements
      namelist /INPARA1/instop,fnin,fnout,ifnpo,dtr,touts,idprop
     1     ,fdirtmp,angref,idp_in,idp_prop,idp_out,xmin1,xmax1
     
c---setting
      angref(:)=-1.e10
      
      open(unit=50,file='namelist',status='old',form='formatted')
      read(50,nml=INPARA1)
      close(50)

c----set angle
      do i=1,30
       if(angref(i).eq.-1.e10) goto 319
      enddo
319   continue
      nang=i-1
      
      allocate(angref6(nang))
!!!! elena :  dimensions mismatch 6/30      
      angref6(1:nang)=angref(1:nang)

c----set filename
      do iang=1,nang
       angref0=angref6(iang)
       write(fntmp(iang),'(a4,i3.3,a7)')
     1       '/tmp',int(angref0),'_in.txt'
       write(fntmp2(iang),'(a4,i3.3,a8)')
     1       '/tmp',int(angref0),'_out.txt'
       fntmp(iang)=trim(fdirtmp)//fntmp(iang)
       fntmp2(iang)=trim(fdirtmp)//fntmp2(iang)
      enddo

c---prepare dataset at angref0
      if(idp_in.eq.1)
     1 call sub_indat(fnin,fntmp(1:nang),angref6,nang,dtr,instopfin
     1   ,ifnpo,idprop)
      if(instop.eq.0) call count_lines(fntmp(1),instop)

c---solar wind propagation along angref0
       do iang=1,nang         
       angref0=angref6(iang)
       fnin1=fntmp(iang)
       fnout1=fntmp2(iang)
       if(idp_prop.eq.1)
     1  call swmodel(angref0,instop,fnin1,fnout1,ifnpo,dtr,touts,idprop
     1   ,xmax1,xmin1)
      enddo
c---read & select dataset
      if(idp_out.eq.1)
     1 call sub_slct(fnin,fnout,fntmp2(1:nang),angref6,nang
     1     ,touts,dtr,idprop,ifnpo)

      stop
      end
c======================================================================|
      subroutine swmodel(angref0,instop,fnin,fnout,ifnpo,dtr,touts0
     1   ,idprop ,xmax1,xmin1)  !--REV0318
c======================================================================|
      implicit double precision (a-h,o-z)
      real*8,allocatable,dimension(:) :: x,xm,dx,dxm,ro,pr,vx,vy,by,bx
     1  ,bxm,sc,scm,dsc,dscm,dv,gx,gxm,rr,rrm,drr,drrm,vz,bz,gy,gz
      integer :: idprop  !--REV0318
      integer :: xin,xout,ifnpo,fnum,fnum2,fnpo,merr
      real*8 :: touts0,touts
      real*8 :: p8bd(8),parain(13),paranm(8)
      real*8,allocatable :: jsinref(:),parainref(:,:)
      real*8 :: bx0at1au,angref0,ns,gm
      character*100 :: fnin,fnout
!!!!! elena
      real*8 jjs

      namelist /INPARA2/npout,bx0at1au,gm 
      open(unit=50,file='namelist',status='old',form='formatted')
      read(50,nml=INPARA2)
      close(50)
       
      touts=touts0*0.5
      ix=int((xmax1-xmin1)*400.+1.e-10)
      xear=(1.-xmin1)*400.

      allocate(x(ix),xm(ix),dx(ix),dxm(ix)
     1   ,ro(ix),pr(ix),vx(ix),vy(ix),by(ix),bx(ix),bxm(ix)
     1   ,sc(ix),scm(ix),dsc(ix),dscm(ix),dv(ix),gx(ix),gxm(ix)
     1   ,rr(ix),rrm(ix),drr(ix),drrm(ix)
     1   ,vz(ix),bz(ix),gy(ix),gz(ix))

c----------------------------------------------------------------------|
c     prologue
c----------------------------------------------------------------------|
c----------------------------------------------------------------------|
c   time control parametersn
c   instop : number of time steps taken from input data
c   (nstop : number of total time steps for the run, = instop)
c   touts : data output's rate (ex.12 with dtr=300 -> 12x300=3600 sec.)
c   npout : prosess print rate (<nstop!!)
c   dtr : real time step (ex. 300 sec.)
c      instop=131600 !***
c      touts=12.
c      npout=100
c      dtr=300.0
      bx0=bx0at1au/16.34d0
c----------------------------------------------------------------------|
c  file open
c  input file(32):solar wind at input, in/out positions (real value)
      fnum=32
      open(fnum,file=trim(fnin),status='old')

c  output file(30):solar wind at specified position (real value)
      fnum2=30
      open(fnum2,file=trim(fnout)
     &     ,status='unknown',form='formatted')

c  monitor file (ex. fnpo=61 -> fort.61)
      fnpo=ifnpo
c----------------------------------------------------------------------|
c   parameters
c   margin: 1 in MLW, 2 in Roe, 4 in CIP
      margin=2
c----------------------------------------------------------------------|
      nstop=instop
      allocate(jsinref(instop),parainref(6,instop))
      paranm=(/10.,10000.,400.,400.,400.,1.,1.,1./)
c----------------------------------------------------------------------|
c  initialize counters
      time  = 0.0
      timep = 0.0
      ns    = 0.
      merr  = 0
c----------------------------------------------------------------------|
c   setup numerical model (grid, initial conditions, etc.)
      call model(idf,ro,pr,vx,vy,vz,by,bz,vstar,bx,bxm,gm,gx,gxm,gy,gz
     &   ,rr,rrm,sc,scm,dv,margin,x,ix,bx0,xear,xmin1,xmax1)
      if(idprop.eq.1) 
     &  p8bd(:)=(/ro(1),pr(1),vx(1),vy(1),vz(1),bx(1),by(1),bz(1)/)
      if(idprop.eq.-1)
     &  p8bd(:)=(/ro(ix),pr(ix),-vx(ix),vy(ix),vz(ix)
     &  ,bx(ix),by(ix),bz(ix)/)
      vx=vx*float(idprop) !--REV0318
      call bnd3(margin,ro,pr,vx,vy,vz,by,bz,gm,vstar,ix,p8bd,idprop) !--REV0318 
c----------------------------------------------------------------------|
c     ready
      call grdrdy(dx,xm,dxm,x,ix)
      call scrdy(dsc,dscm,sc,scm,dx,dxm,ix)
      call scrdy(drr,drrm,rr,rrm,dx,dxm,ix)
      call bndsc(margin,sc,dsc,scm,dscm,ix)
c----------------------------------------------------------------------|
      l=0.
      js=0.d0
      k=10
c----------------------------------------------------------------------|
c     time integration
c----------------------------------------------------------------------|
1000  continue  
         ns = ns+1.
c----------------------------------------------------------------------|
c     obtain time spacing
      call cfl_m3(dt,merr,gm,bx,ro,pr,vx,vy,vz,by,bz,dx,ix)

      dt0=dtr/375000.d0
      if(dt0.lt.dt) then
        dt=dt0
      else
        write(fnpo,*) 'Time Spacing is NOT Correct! : dt0 > dt'
        write(fnpo,*) 'dt=',dt
        write(fnpo,*) 'dt0=',dt0       
      endif
      if (merr.ne.0) goto 9999


c----data read & input----------
      read(fnum,*) jjs,(parain(i),i=1,13)
!!!!! elena : js is defined as REAL in the file
      js=int(jjs)
      parain(3)=parain(3)*float(idprop)  !--REV0318
      xin=minval(minloc(abs(parain(9)-x)))
      if(parain(9).lt.x(xin).and.idprop.eq.1) xin=xin-1 !--REV0318
      if(parain(9).gt.x(xin).and.idprop.eq.-1) xin=xin+1 !--REV0318

      jsinref(ns)=float(js)
      parainref(1:2,ns)=parain(9:10)
      parainref(3,ns)=x(xin)
      parainref(4,ns)=parain(3) !km/s
      
      call datin_varied3(ro,pr,vx,vy,vz,bx,by,bz
     &        ,vstar,ns,ix,parain,xin,paranm ,idprop,x)  !--REV0318

c----initial setting
      if(ns.eq.1.) js0=js
      if(ns.eq.1.) write(fnpo,*) 'js=',js
      if(ns.eq.1.) then
        parainref(5,1)=parain(10)
      else
        dang0=parain(10)-parainref(2,ns-1)
        if(dang0.le.-180.) dang0=dang0+360.
        parainref(5,ns)=parainref(5,ns-1)+dang0
      endif
        parainref(6,ns)=parain(13)

c----------------------------------------------------------------------|
c     solve hydrodynamic equations
      call roe_m_bg(ro,pr,vx,vy,by,bx,bxm,dt,gm
     &            ,gx,dsc,scm,dv,rr,rrm,drr,dx,ix)
      call bnd3_2(margin,ro,pr,vx,vy,vz,by,bz,gm,vstar,ix,idprop,x) !--REV0318
c----------------------------------------------------------------------|
c     data output 
       
c---monitoring
       if(minval(pr).lt.0.) write(6,*) ns
       if(minval(pr).lt.0.) write(6,*) minval(pr)
       if(minval(pr).lt.0.) stop 'Pr<0' 
       if(idprop.eq.1.and.minval(vx).lt.0.) stop 'Vx opposite'
       if(idprop.eq.-1.and.maxval(vx).gt.0.) then
         write(6,1115) vx
         stop 'Vx opposite'
       endif
1115   format(15f10.7)
c----first dt estimation
      if (mod(int(ns),int(touts)).eq.0) then
       js=js0+ns*dtr *float(idprop) !--REV0318
       xout=minval(minloc(abs(parain(11)-x)))
       dtp1=(parain(11)-parain(9))*1.5d11/400./1.d3
       
!       write(6,*) 'first dt estimation', dtp1
       
       if((js-dtp1.le.js0.and.idprop.eq.1) !--REV0318
     1    .or.(js+dtp1.ge.js0.and.idprop.eq.-1)) then !--REV0318
         write(fnum2,704) js,1.,1.,100.,1.,1.,1.d-5,1.d-5,1.d-5
     &    ,0.d0,float(js0)
       else
         itin=minval(minloc(abs((js-dtp1*float(idprop))-jsinref))) !--REV0318
         
         !!!!! elena
         if(parainref(4,itin).eq.0.0)then        
           write(6,*) '!!!second dt',itin, ns, parain(12), parain(10)
           itin = ns
         endif
         
c----second dt estimation
         dtp1=(parain(11)-parain(9))*1.5d11/parainref(4,itin)/1.d3
                      
        if((js-dtp1.le.js0.and.idprop.eq.1)
     1    .or.(js+dtp1.ge.js0.and.idprop.eq.-1)) then
         write(fnum2,704) js,1.,1.,100.,1.,1.,1.d-5,1.d-5,1.d-5
     &    ,0.d0,float(js0)
         goto 1001
        endif
        
        itin_tmp=itin
        itin=minval(minloc(abs((js-dtp1*float(idprop))-jsinref)))
        
        !!!!! elena
        if (parainref(4,itin).eq.0.0) itin=itin_tmp
c----output
         dang0=-(parain(12)-parainref(2,itin))
         dang0=mod(dang0,360.)
         if(dang0.ge.180.) dang0=dang0-360.
         
         if(angref0.ge.0.)
     &     dtp2=(parain(11)-x(xout))*1.5d11/vx(xout)/400.d3
     &      -(parainref(1,itin)-parainref(3,itin))*1.5d11 !x_datain-x_in
     &        /parainref(4,itin)/1.d3  
          !!!!  elena   
          if (parainref(4,itin).eq.0.0) then ! NaN 
            write(6,*) '!!!!!! NaN in dtp2 ', angref0
            write(6,*) parain(11), x(xout),vx(xout)           
            write(6,*) ns, itin, xout 
            write(6,*) parain(12), parain(10)
            stop           
          endif
          
         write(fnum2,704) js,ro(xout)*paranm(1)
     &      ,pr(xout)*1937./ro(xout)*paranm(2)
     &      ,vx(xout)*paranm(3)*float(idprop)  !--REV0318
     &      ,vy(xout)*paranm(4),vz(xout)*paranm(5)
     &	    ,bx(xout)*paranm(6)*16.34,by(xout)*paranm(7)*16.34 
     &      ,bz(xout)*paranm(8)*16.34
     &	    ,dtp2*float(idprop),jsinref(itin)-parainref(6,itin)
       endif
        
704    format(i15.1, 9e18.8,f18.1)

1001   continue
      endif

      if (mod(int(ns),int(nstop/npout)).eq.0) then
      	l=l+1.
      	write(fnpo,'(i6,a1,i6,a25,f7.1,a1)') 
     1    l,'/',npout,' :processed ! (angref=',angref0,')'
      endif
c----------------------------------------------------------------------|
      if (ns .lt. nstop) goto 1000
c----------------------------------------------------------------------|
c     epilogue
c----------------------------------------------------------------------|
9999  continue

      write(fnpo,*) 'js=',js
      CLOSE (fnum, STATUS = 'KEEP')
      CLOSE (fnum2, STATUS = 'KEEP')
      deallocate(jsinref,parainref)

      deallocate(x,xm,dx,dxm,ro,pr,vx,vy,by,bx,bxm,sc,scm,dsc,dscm
     1   ,dv,gx,gxm,rr,rrm,drr,drrm,vz,bz,gy,gz)

      return
      end
c======================================================================|
      subroutine model(idf,ro,pr,vx,vy,vz,by,bz
     &  ,vstar,bx,bxm,gm,gx,gxm,gy,gz,rr,rrm,sc,scm,dv,margin,x,ix,bx0
     &  ,xear,xmin1,xmax1)
c======================================================================|
      implicit double precision (a-h,o-z)
c----------------------------------------------------------------------|
      dimension x(ix),dxm(ix)
      dimension ro(ix),pr(ix),vx(ix),vy(ix),by(ix)
      dimension sc(ix),scm(ix),dv(ix)
      dimension bx(ix),bxm(ix)
      dimension gx(ix),gxm(ix)
      dimension rr(ix),rrm(ix)
      dimension vz(ix),bz(ix)
      dimension gy(ix),gz(ix)
c----------------------------------------------------------------------|
c   parameters
c----------------------------------------------------------------------|
      gmst=0.0055248d0
      vstar=0.d0

      rstar=xmin1
      rmax=xmax1
c-----------------------------------------------------------------------
c     grid
c-----------------------------------------------------------------------
      dx0=(rmax-rstar)/real(ix-margin*2)
c-----------------------------------------------------------------------
c      dxm,x
      do i=1,ix
         dxm(i)=dx0
      enddo
      izero=margin
      x(izero)=rstar-dxm(izero)/2.d0
      do i=izero+1,ix
         x(i) = x(i-1)+dxm(i-1)
      enddo
      do i=izero-1,1,-1
         x(i) = x(i+1)-dxm(i)
      enddo
c----------------------------------------------------------------------|
c       gravity
c----------------------------------------------------------------------|
      do i=1,ix
        gx(i)=-1.d0*gmst/x(i)**2d0
        gxm(i)=-1.d0*gmst/(x(i)+0.5d0*dxm(i))**2
        gy(i)=0.d0
        gz(i)=0.d0
      enddo
c----------------------------------------------------------------------|
c   rotation
c----------------------------------------------------------------------|
      do i=1,ix
        rr(i)=x(i)
        rrm(i)=x(i)+dxm(i)/2.d0
      enddo
c-----------------------------------------------------------------------|
c   initial temperature, density, pressure distributions
c-----------------------------------------------------------------------|
      ro=1./x**2
      do i=1,ix
        te=1.d0/x(i)**0.79
        pr(i)=ro(i)*te/1937.d0
      enddo
      vx(:)=1.; vy(:)=0.; vz(:)=0.
c----------------------------------------------------------------------|
c       magnetic field
c----------------------------------------------------------------------|
      pi = acos(-1.0d0)
      do i=1,ix
        by(i)=0.d0
	bz(i)=0.d0
        bx(i)=bx0*x(xear)**2/x(i)**2
        bxm(i)=bx0*x(xear)**2/(x(i)+dxm(i)/2.d0)**2
      enddo
c----------------------------------------------------------------------|
c       cross section of flux tube
c----------------------------------------------------------------------|
      do i=1,ix
        sc(i)=x(i)**2/x(xear)**2
        scm(i)=(x(i)+dxm(i)/2.)**2/x(xear)**2
        dv(i)=sc(i)
      enddo
      return
      end
c======================================================================|
      subroutine bdcnsx(mbnd,margin,qq,q0,ix)
c======================================================================|
c
c NAME  bdcnsx
c
c PURPOSE
c    apply constant-value boundary condition
c
c INPUTS & OUTPUTS
c    qq(ix): [double] variable
c
c OUTPUTS
c    None
c
c INPUTS
c    ix: [integer] dimension size
c    margin: [integer] margin, i.e. # of grid points outside the boundary
c    mbnd: [integer] If mbnd=0, smaller 'i' side. 
c                    If mbnd=1, larger  'i' side.
c    q0: [double] constant boundary value to be taken 
c
c HISTORY
c    written 2002-3-1 T. Yokoyama
c
c----------------------------------------------------------------------|
      implicit double precision (a-h,o-z)
      dimension qq(ix)
c----------------------------------------------------------------------|
      if (mbnd.eq.0) then 
        ibnd=1+margin
        do i=1,margin
          qq(ibnd-i) = q0
        enddo
      else
        ibnd=ix-margin
        do i=1,margin
          qq(ibnd+i) = q0
        enddo
      endif      
      return
      end
c======================================================================|
      subroutine bdfrdx(mbnd,margin,qq,dxm,ix)
c======================================================================|
c
c NAME  bdfrdx
c
c PURPOSE
c    apply free boundary condition
c    values are extended to have constant gradient
c
c INPUTS & OUTPUTS
c    qq(ix): [double] variable
c
c OUTPUTS
c    None
c
c INPUTS
c    ix: [integer] dimension size
c    margin: [integer] margin, i.e. # of grid points outside the boundary
c    mbnd: [integer] If mbnd=0, smaller 'i' side. 
c                    If mbnd=1, larger  'i' side.
c
c HISTORY
c    written 2002-3-1 T. Yokoyama
c
c----------------------------------------------------------------------|
      implicit double precision (a-h,o-z)
      dimension qq(ix),dxm(ix)
c----------------------------------------------------------------------|
      if (mbnd.eq.0) then 
        ibnd=1+margin
        dqq=(qq(ibnd+1)-qq(ibnd))/dxm(ibnd)
        do i=1,margin 
          qq(ibnd-i) = qq(ibnd-i+1)-dqq*dxm(ibnd-i)
        enddo
      else
        ibnd=ix-margin
        dqq=(qq(ibnd)-qq(ibnd-1))/dxm(ibnd-1)
        do i=1,margin
          qq(ibnd+i) = qq(ibnd+i-1)+dqq*dxm(ibnd+i-1)
        enddo
      endif
      
      return
      end
c======================================================================|
      subroutine bdfrex(mbnd,margin,qq,ix)
c======================================================================|
c
c NAME  bdfrex
c
c PURPOSE
c    apply free boundary condition
c
c INPUTS & OUTPUTS
c    qq(ix): [double] variable
c
c OUTPUTS
c    None
c
c INPUTS
c    ix: [integer] dimension size
c    margin: [integer] margin, i.e. # of grid points outside the boundary
c    mbnd: [integer] If mbnd=0, smaller 'i' side. 
c                    If mbnd=1, larger  'i' side.
c
c HISTORY
c    written 2002-3-1 T. Yokoyama
c
c----------------------------------------------------------------------|
      implicit double precision (a-h,o-z)
      dimension qq(ix)
c----------------------------------------------------------------------|
      if (mbnd.eq.0) then 
        ibnd=1+margin
        do i=1,margin 
          qq(ibnd-i) = qq(ibnd)
        enddo
      else
        ibnd=ix-margin
        do i=1,margin
          qq(ibnd+i) = qq(ibnd)
        enddo
      endif
      return
      end
c======================================================================|
      subroutine bnd3(margin,ro,pr,vx,vy,vz,by,bz,gm,vstar,ix 
     1  ,p8bd,idprop) !--REV0318 
c======================================================================|
c     apply boundary condition 
c----------------------------------------------------------------------|      
      implicit double precision (a-h,o-z)
      dimension ro(ix),pr(ix),vx(ix),vy(ix),by(ix)
      dimension vz(ix),bz(ix)
      real*8 p8bd(8)
      integer :: idprop  !--REV0318 
c----------------------------------------------------------------------|      
      ro0=p8bd(1)
      pr0=p8bd(2)
      vstar=p8bd(4)
      if(idprop.eq.1)then !--REV0318
       call bdcnsx(0,margin,ro,ro0,ix)
       call bdcnsx(0,margin,pr,pr0,ix)
       call bdfrex(0,margin,vx,ix)
       call bdcnsx(0,margin,vy,vstar,ix)
       call bdfrex(0,margin,by,ix)
       call bdfrex(0,margin,vz,ix)
       call bdfrex(0,margin,bz,ix)
       call bdfrex(1,margin,ro,ix)
       call bdfrex(1,margin,pr,ix)
       call bdfrex(1,margin,vx,ix)
       call bdfrex(1,margin,vy,ix)
       call bdfrex(1,margin,by,ix)
       call bdfrex(1,margin,vz,ix)
       call bdfrex(1,margin,bz,ix)
      else
c       call bdcnsx(1,margin,ro,ro0,ix)
c       call bdcnsx(1,margin,pr,pr0,ix)
       call bdfrex(1,margin,ro,ix)
       call bdfrex(1,margin,pr,ix)
       call bdfrex(1,margin,vx,ix)
c       call bdcnsx(1,margin,vy,vstar,ix)
       call bdfrex(1,margin,vy,ix)
       call bdfrex(1,margin,by,ix)
       call bdfrex(1,margin,vz,ix)
       call bdfrex(1,margin,bz,ix)
       call bdfrex(0,margin,ro,ix)
       call bdfrex(0,margin,pr,ix)
       call bdfrex(0,margin,vx,ix)
       call bdfrex(0,margin,vy,ix)
       call bdfrex(0,margin,by,ix)
       call bdfrex(0,margin,vz,ix)
       call bdfrex(0,margin,bz,ix)
      endif
      return
      end
c======================================================================|
      subroutine bnd3_2(margin,ro,pr,vx,vy,vz,by,bz,gm,vstar,ix
     1   ,idprop,x)
c======================================================================|
c     apply boundary condition 
c----------------------------------------------------------------------|      
      implicit double precision (a-h,o-z)
      dimension ro(ix),pr(ix),vx(ix),vy(ix),by(ix)
      dimension vz(ix),bz(ix)
      dimension x(ix)
      integer :: idprop !--REV0318
c----------------------------------------------------------------------|      
      if(idprop.eq.1) then
       call bdfrex(1,margin,ro,ix)
       call bdfrex(1,margin,pr,ix)
       call bdfrex(1,margin,vx,ix)
       call bdfrex(1,margin,vy,ix)
       call bdfrex(1,margin,by,ix)
       call bdfrex(1,margin,vz,ix)
       call bdfrex(1,margin,bz,ix)
      else
       call bdfrex(0,margin,ro,ix)
c       ro(1)=ro(3)*x(3)**2/x(1)**2
c       ro(2)=ro(3)*x(3)**2/x(2)**2
       call bdfrex(0,margin,pr,ix)
       call bdfrex(0,margin,vx,ix)
       call bdfrex(0,margin,vy,ix)
       call bdfrex(0,margin,by,ix)
       call bdfrex(0,margin,vz,ix)
       call bdfrex(0,margin,bz,ix)
c       vx(1:2)=-1.
      endif
      return
      end
c======================================================================|
      subroutine bndsc(margin,sc,dsc,scm,dscm,ix)
c======================================================================|
c     apply boundary condition 
c----------------------------------------------------------------------|      
      implicit double precision (a-h,o-z)
      dimension sc(ix),dsc(ix)
      dimension scm(ix),dscm(ix)
c----------------------------------------------------------------------|      
      call bdfrex(0,margin,dsc,ix)
      call bdfrex(0,margin-1,dscm,ix)
      call bdfrex(1,margin,dsc,ix)
      call bdfrex(1,margin,dscm,ix)
      return
      end
c======================================================================|
      subroutine cfl_m3(dt,merr,gm,bx,ro,pr,vx,vy,vz,by,bz,dx,ix)
c======================================================================|
c 
c NAME  cfl_m
c
c PURPOSE
c    determine time step such that it satisfies CFL condition.
c        * MHD equations
c
c OUTPUTS
c    dt: [double] delta time
c    merr: [integer] error code, merr=0 is nominal.
c
c INPUTS
c    ix: [integer] dimension size
c    ro(ix): [double] density
c    pr(ix): [double] pressure
c    vx(ix): [double] velocity
c    vy(ix): [double] velocity
c    bx(ix): [double] magnetic field
c    by(ix): [double] magnetic field
c    dx(ix): [double] grid spacing
c    gm: [double] polytropic index gamma
c
c HISTORY
c    written 2002-3-1 T. Yokoyama
c
c----------------------------------------------------------------------|
c     determine time step such that it satisfies cfl condition.
c----------------------------------------------------------------------|      
      implicit double precision (a-h,o-z)
      dimension dx(ix)
      dimension ro(ix)
      dimension pr(ix)
      dimension vx(ix)
      dimension vy(ix)
      dimension bx(ix)
      dimension by(ix)
      dimension dtq(ix)
      dimension vz(ix),bz(ix)
c----------------------------------------------------------------------|
      pi = acos(-1.0d0)
      pi4i=0.25/pi
      dtmin=2.0e-10
!      safety=0.8 !cf. 0.4 --original 
      safety=0.4
c----------------------------------------------------------------------|
      dt=1.e20
      imin   = 0
      do i=2,ix-1
         onero = 1.0/ro(i)
         v2 = vx(i)*vx(i)+vy(i)*vy(i)+vz(i)*vz(i)
         ca2 = (bx(i)*bx(i)+by(i)*by(i)+bz(i)*bz(i))*pi4i*onero
         cs2 = gm*pr(i)*onero
         dtcfl = dx(i)/sqrt(v2+cs2+ca2)
         dtq(i)=safety*dtcfl
      enddo
      do i=2,ix-1
        if(dtq(i).lt.dt) then
          imin=i
          dt=dtq(i)
        endif
      enddo
!      write(6,*) minval(pr)
c----------------------------------------------------------------------|
c     write the point where dt is smaller than critical value    
c----------------------------------------------------------------------|
      merr=0
      if (dt.lt.dtmin) then
         merr=9001
         write(6,*) '  ### stop due to small dt, less than dtmin ###'
         write(6,620) dt,dtmin,imin
 620     format('   dt = ',1pe10.3,'  < ',1pe10.3,' @ i =',i5) 
      endif
      return
      end
c======================================================================|
      subroutine datin_varied3(ro,pr,vx,vy,vz,bx,by,bz,
     &		vstar,ns,ix,parain,xear,paranm ,idprop,x) !--REV0318
c======================================================================|
      implicit double precision (a-h,o-z)
      dimension ro(ix),pr(ix),vx(ix),vy(ix),vz(ix)
      dimension bx(ix),by(ix),bz(ix)
      dimension x(ix)
      real*8 parain(12),paranm(8)
      integer :: xear,idprop
      if(idprop.eq.1)then
cc       ro(1:xear)=parain(1)/paranm(1)
cc       pr(1:xear)=parain(2)/paranm(2)*parain(1)/paranm(1)/1937.d0
       ro(1:xear)=parain(1)/paranm(1)/x(1:xear)**2*x(xear)**2
       pr(1:xear)=parain(2)/paranm(2)*ro(1:xear)/1937.d0
       vx(1:xear)=parain(3)/paranm(3)
       vy(1:xear)=parain(4)/paranm(4)
cc       vz(1:xear)=parain(5)/paranm(5)
       by(1:xear)=parain(7)/16.34d0
cc       bz(1:xear)=parain(8)/16.34d0
      else
       ro(xear:ix)=parain(1)/paranm(1)
       pr(xear:ix)=parain(2)/paranm(2)*parain(1)/paranm(1)/1937.d0
cc       ro(xear:ix)=parain(1)/paranm(1)/x(xear:ix)**2*x(xear)**2
cc       pr(xear:ix)=parain(2)/paranm(2)/ro(xear:ix)/1937.d0
       vx(xear:ix)=parain(3)/paranm(3)
       vy(xear:ix)=parain(4)/paranm(4)
cc       vz(xear:ix)=parain(5)/paranm(5)
       by(xear:ix)=parain(7)/16.34d0
cc       bz(xear:ix)=parain(8)/16.34d0
      endif
      return
      end
c======================================================================|
      subroutine datin_const3(ro,pr,vx,vy,vz,bx,by,bz,vstar,ns,ix)
c======================================================================|
      implicit double precision (a-h,o-z)
      dimension ro(ix),pr(ix),vx(ix),vy(ix),vz(ix)
      dimension bx(ix),by(ix),bz(ix)
      ro(1:xear)=1.d0
      pr(1:xear)=1.d0/1937.d0
      vx(1:xear)=1.d0
      vy(1:xear)=vstar
      vz(1:xear)=0.d0
      bx(1:xear)=0.0d0
      by(1:xear)=0.0d0
      bz(1:xear)=0.0d0
      ro(xear+1)=ro(xear)
      pr(xear+1)=pr(xear)
      vx(xear+1)=vx(xear)
      vy(xear+1)=vy(xear)
      vz(xear+1)=vz(xear)
      bx(xear+1)=bx(xear)
      by(xear+1)=by(xear)
      bz(xear+1)=bz(xear)
      return
      end
c======================================================================|
      subroutine grdrdy(dx,xm,dxm,x,ix)
c======================================================================|
c
c NAME  grdrdy
c
c PURPOSE
c    calculate coordinate of mid-grid points and
c    grid spacing on the grid points
c
c OUTPUTS
c    dx(ix),dxm(ix): [double] grid spacing
c    xm(ix): [double] coordinate
c
c INPUTS
c    x(ix): [double] coordinate
c    ix: [integer] dimension size
c
c HISTORY
c    written 2002-3-1 T. Yokoyama
c
c----------------------------------------------------------------------|
      implicit double precision (a-h,o-z)
      dimension dx(ix),dxm(ix)
      dimension x(ix),xm(ix)
c----------------------------------------------------------------------|
      do i=1,ix-1
         dxm(i)=x(i+1)-x(i)
      enddo
      dxm(ix)=dxm(ix-1)

      do i=2,ix-1
         dx(i)  = 0.5*(dxm(i-1)+dxm(i))
      enddo
      dx(1)=dx(2)
      dx(ix)=dx(ix-1)

      do i=1,ix-1
         xm(i)=0.5*(x(i)+x(i+1))
      enddo
      xm(ix)=xm(ix-1)+dx(ix-1)
      return
      end
c======================================================================|
      subroutine roe_m_bg(ro,pr,vx,vy,by,bx,bxm,dt,gm
     &            ,gx,dsc,scm,dv,rr,rrm,drr,dx,ix)
c======================================================================|
c
c NAME  roe_m_bg
c
c PURPOSE
c    solve eqs. by modified Roe + MUSCL-TVD  method with effects of
c        * MHD
c        * axial symmetry
c        * non-uniform poloidal magnetic field
c        * gravity
c
c INPUTS & OUTPUTS
c    ro(ix): [double] density
c    pr(ix): [double] pressure
c    vx(ix): [double] velocity 
c    vy(ix): [double] velocity 
c    by(ix): [double] magnetic field
c
c OUTPUTS
c    None
c
c INPUTS
c    NOTE: ??m(ix) is the variable array defined at grid bounds
c
c    bx(ix), bxm(ix) : [double] magnetic field
c    gx(ix), gxm(ix) : [double] gravity
c    gy(ix), gym(ix) : [double] gravity
c    gz(ix), gzm(ix) : [double] gravity
c    scm(ix) : [double] cross section
c    dsc(ix), dscm(ix) : [double] cross section gradient
c    rr(ix), rrm(ix) : [double] distance from rotation axis
c    drr(ix), drrm(ix) : [double] distance gradient from rotation axis
c    dx(ix) : [double] grid spacing
c    gm: [double] polytropic index gamma
c    dt: [double] delta time
c    ix: [integer] dimension size
c
c HISTORY
c    written 2002-3-1 T. Yokoyama based on N. Fukuda's code
c
c----------------------------------------------------------------------|
      implicit double precision (a-h,o-z)
      dimension dx(ix)
      dimension ro(ix),pr(ix),vx(ix),vy(ix),by(ix)
      dimension bx(ix),bxm(ix)
      dimension dsc(ix),scm(ix),dv(ix)
      dimension rosc(ix),eesc(ix),rxsc(ix),rysc(ix),bysc(ix)
      dimension rosch(ix),eesch(ix),rxsch(ix),rysch(ix),bysch(ix)
      dimension fro(ix),fee(ix),frx(ix),fry(ix),fby(ix)
      dimension roh(ix),prh(ix),vxh(ix),vyh(ix),byh(ix)
      dimension row(ix,2),prw(ix,2),vxw(ix,2),vyw(ix,2)
     &         ,bxw(ix,2),byw(ix,2)
      dimension gx(ix)
      dimension rr(ix),drr(ix),rrm(ix)
c----------------------------------------------------------------------|      
c     numerical parameters
      pi = acos(-1.0d0)
      pi4=4.0d0*pi
      pi8=8.0d0*pi
      pi4i=1.0d0/pi4
      pi8i=5.0d-1*pi4i
c----------------------------------------------------------------------|
c     computation of conservative variables w(i,l)
      do i=1,ix
         rosc(i)=dv(i)*ro(i)
         rxsc(i)=dv(i)*ro(i)*vx(i)
         rysc(i)=dv(i)*ro(i)*vy(i)*rr(i)
         bysc(i)=dv(i)*by(i)/rr(i)
         v2=vx(i)**2+vy(i)**2
         b2=bx(i)**2+by(i)**2
         eesc(i)=dv(i)*(pr(i)/(gm-1.0d0) +0.5d0*ro(i)*v2 + pi8i*b2)
      enddo
c----------------------------------------------------------------------|
c     proceed half step
c     computation of 1st order flux f(i,l)
c----------------------------------------------------------------------|
      do i=1,ix-1
         row(i,1)=ro(i)
         prw(i,1)=pr(i)
         vxw(i,1)=vx(i)
         vyw(i,1)=vy(i)
         bxw(i,1)=bxm(i)
         byw(i,1)=by(i)
         row(i,2)=ro(i+1)
         prw(i,2)=pr(i+1)
         vxw(i,2)=vx(i+1)
         vyw(i,2)=vy(i+1)
         bxw(i,2)=bxm(i)
         byw(i,2)=by(i+1)
      enddo

      call roeflux_m(fro,fee,frx,fry,fby,gm,row,prw,vxw,vyw,bxw,byw,ix)
      do i=1,ix-1
        fro(i)=scm(i)*fro(i)
        fee(i)=scm(i)*fee(i)
        frx(i)=scm(i)*frx(i)
        fry(i)=scm(i)*fry(i)*rrm(i)
        fby(i)=scm(i)*fby(i)/rrm(i)
      enddo

      do i=2,ix-1
         rosch(i)=rosc(i)+0.5d0*dt*( (fro(i-1)-fro(i))/dx(i) )
         eesch(i)=eesc(i)+0.5d0*dt*( (fee(i-1)-fee(i))/dx(i) )
         rxsch(i)=rxsc(i)+0.5d0*dt*( (frx(i-1)-frx(i))/dx(i) )
         rysch(i)=rysc(i)+0.5d0*dt*( (fry(i-1)-fry(i))/dx(i) )
         bysch(i)=bysc(i)+0.5d0*dt*( (fby(i-1)-fby(i))/dx(i) )
      enddo

      do i=2,ix-1
         see=dv(i)*ro(i)*vx(i)*gx(i)
         eesch(i)=eesch(i)+0.5d0*dt*see
         b2=bx(i)**2+by(i)**2
         srx=dv(i)
     &   *(ro(i)*gx(i)+(ro(i)*vy(i)**2-by(i)**2*pi4i)/rr(i)*drr(i))
     &       +(pr(i)+b2*pi8i)*dsc(i)
         rxsch(i)=rxsch(i)+0.5d0*dt*srx
      enddo

c     computation of basic variables on half step

      do i=2,ix-1
         roh(i)=rosch(i)/dv(i)
         vxh(i)=rxsch(i)/rosch(i)
         vyh(i)=rysch(i)/rosch(i)/rr(i)
         byh(i)=bysch(i)/dv(i)*rr(i)
         v2=vxh(i)**2+vyh(i)**2
         b2= bx(i)**2+byh(i)**2
         prh(i)=(gm-1.0d0)* (eesch(i)/dv(i)-0.5d0*roh(i)*v2 -b2*pi8i)
      enddo

c----------------------------------------------------------------------|
c     proceed full step
c     computation of 2nd order flux f(i,l)
c----------------------------------------------------------------------|
      call tvdminmod(roh,row,ix)
      call tvdminmod(prh,prw,ix)
      call tvdminmod(vxh,vxw,ix)
      call tvdminmod(vyh,vyw,ix)
      call tvdminmod(byh,byw,ix)

      call roeflux_m(fro,fee,frx,fry,fby,gm,row,prw,vxw,vyw,bxw,byw,ix)
      do i=2,ix-2
        fro(i)=scm(i)*fro(i)
        fee(i)=scm(i)*fee(i)
        frx(i)=scm(i)*frx(i)
        fry(i)=scm(i)*fry(i)*rrm(i)
        fby(i)=scm(i)*fby(i)/rrm(i)
      enddo

      do i=2,ix-2
         rosc(i)=rosc(i)+dt*( (fro(i-1)-fro(i))/dx(i) )
         eesc(i)=eesc(i)+dt*( (fee(i-1)-fee(i))/dx(i) )
         rxsc(i)=rxsc(i)+dt*( (frx(i-1)-frx(i))/dx(i) )
         rysc(i)=rysc(i)+dt*( (fry(i-1)-fry(i))/dx(i) )
         bysc(i)=bysc(i)+dt*( (fby(i-1)-fby(i))/dx(i) )
      enddo

      do i=3,ix-2
         see=dv(i)*roh(i)*vxh(i)*gx(i)
         eesc(i)=eesc(i)+dt*see
         b2=bx(i)**2+byh(i)**2
         srx=dv(i)
     &   *(roh(i)*gx(i)+(roh(i)*vyh(i)**2-byh(i)**2*pi4i)/rrm(i)*drr(i))
     &       +(prh(i)+b2*pi8i)*dsc(i)
         rxsc(i)=rxsc(i)+dt*srx
      enddo

c----------------------------------------------------------------------|
c     computation of basic variables on full step
      do i=3,ix-2
         ro(i)=rosc(i)/dv(i)
         vx(i)=rxsc(i)/rosc(i)
         vy(i)=rysc(i)/rosc(i)/rr(i)
         by(i)=bysc(i)/dv(i)*rr(i)
         v2=vx(i)**2+vy(i)**2
         b2=bx(i)**2+by(i)**2
         pr(i)=(gm-1.0d0)* (eesc(i)/dv(i)-0.5d0*ro(i)*v2-pi8i*b2)
      enddo
c----------------------------------------------------------------------|
      return
      end
c======================================================================|
      subroutine roeflux_m(fro,fee,frx,fry,fby
     &                         ,gm,row,prw,vxw,vyw,bxw,byw,ix)
c======================================================================|
c
c NAME  roeflux_m
c
c PURPOSE
c    derive numerical flux by solving the linearized Riemann problem
c        * MHD
c
c INPUTS & OUTPUTS
c    None
c
c OUTPUTS
c    fro(ix): [double] density flux
c    fee(ix): [double] total-energy flux
c    frx(ix): [double] momentum flux
c    fry(ix): [double] momentum flux
c    fby(ix): [double] magnetic field flux
c
c INPUTS
c    row(ix,2): [double] density at cell boundary
c    prw(ix,2): [double] pressure at cell boundary
c    vxw(ix,2): [double] velocity at cell boundary
c    vyw(ix,2): [double] velocity at cell boundary
c    byw(ix,2): [double] magnetic field at cell boundary
c    gm: [double] polytropic index gamma
c    ix: [integer] dimension size
c
c HISTORY
c    written 2002-3-1 T. Yokoyama based on N. Fukuda's code
c
c----------------------------------------------------------------------|
      implicit double precision (a-h,o-z)
      dimension row(ix,2),prw(ix,2),vxw(ix,2)
      dimension vyw(ix,2),byw(ix,2),bxw(ix,2)
      dimension fro(ix),fee(ix),frx(ix),fry(ix),fby(ix)
c----------------------------------------------------------------------|
      pi = acos(-1.0d0)
      pi4=4.0d0*pi
      pi4i=1.0d0/pi4
      pi8i=5.0d-1*pi4i

      do i=1,ix-1
         rhol=row(i,1)
         vxl=vxw(i,1)
         vyl=vyw(i,1)
         bxl=bxw(i,1)
         byl=byw(i,1)
         prl=prw(i,1)
         rhor=row(i,2)
         vxr=vxw(i,2)
         vyr=vyw(i,2)
         bxr=bxw(i,2)
         byr=byw(i,2)
         prr=prw(i,2)
c-----roe's variable
      sr0=sqrt(rhol)
      sr1=sqrt(rhor)
      sri=1.0d0/(sr0+sr1)
      rhobar=sr0*sr1
      vxbar=(sr0*vxl+sr1*vxr)*sri
      vybar=(sr0*vyl+sr1*vyr)*sri
      bxbar=(sr0*bxr+sr1*bxl)*sri
      bybar=(sr0*byr+sr1*byl)*sri
      hl=0.5d0*(vxl**2+vyl**2)+gm*prl/((gm-1.0d0)*rhol)
     1  +(bxbar**2+byl**2)/(pi4*rhol)
      hr=0.5d0*(vxr**2+vyr**2)+gm*prr/((gm-1.0d0)*rhor)
     1  +(bxbar**2+byr**2)/(pi4*rhor)
      hbar=(sr0*hl+sr1*hr)*sri
      byave=(byl+byr)/2.0d0
c-----characteristic speed
      delb2=(gm-2.0d0)/(gm-1.0d0)
     1     *((byr-byl)**2)*sri**2*pi8i
      cs2=(gm-1.0d0)*(hbar-0.5d0*(vxbar**2+vybar**2)
     1   -delb2-(bxbar**2+bybar**2)*pi4i/rhobar)
      astar2=(gm-1.0d0)
     1      *(hbar-0.5d0*(vxbar**2+vybar**2)-delb2)
     2      -(gm-2.0d0)*(bxbar**2+bybar**2)*pi4i/rhobar
      ca2=bxbar**2/(pi4*rhobar)
c      cfast2=0.5d0*(astar2+sqrt(astar2**2-4.0d0*cs2*ca2))
      cbr2=(bybar**2)*pi4i/rhobar
      cfast2=0.5d0*(astar2+sqrt(cbr2*(astar2+cs2+ca2)+(cs2-ca2)**2))
      cslow2=cs2*ca2/cfast2
      cfast=sqrt(cfast2)
      cslow=sqrt(cslow2)
      ca=sqrt(ca2)
      cs=sqrt(cs2)
c----- for singular points
      epsi=1.0d-12
      sgr=bybar**2-epsi
      sp=0.5d0+sign(0.5d0,sgr)
      betay=sp*bybar*sqrt(1.0d0/(bybar**2+1.0d0-sp))
     1     +sqrt(0.5d0)*(1.0d0-sp)
      betaz=sqrt(0.5d0)*(1.0d0-sp)
      eps2=1.0d-12
      sgr2=(bybar**2)/(pi4*rhobar)+abs(ca2-cs2)-eps2
      sp2=0.5d0+sign(0.5d0,sgr2)
      cfca=max(0.0d0,cfast2-ca2)
      cfcs=max(0.0d0,cfast2-cslow2)
      cfa=max(0.0d0,cfast2-cs2)
      alphf=sp2*sqrt(cfca/(cfcs+1.0d0-sp2))+1.0d0-sp2
      alphs=sp2*sqrt(cfa/(cfcs+1.0d0-sp2))
      sgnbx=sign(1.0d0,bxbar)
c----- eigen value & entropy condition
      eeps=(vxr-vxl+abs(vxr-vxl))*2.5d-1 
      elpf=-max(abs(vxbar+cfast),eeps)
      elmf=-max(abs(vxbar-cfast),eeps)
      elps=-max(abs(vxbar+cslow),eeps)
      elms=-max(abs(vxbar-cslow),eeps)
      elpa=-max(abs(vxbar+ca),eeps)
      elma=-max(abs(vxbar-ca),eeps)
      elze=-max(abs(vxbar),eeps)
c     elmax=max(abs(elpf),abs(elmf))
c----- amplitude;w's
      drho=rhor-rhol
      du21=rhobar*(vxr-vxl)
      du31=rhobar*(vyr-vyl)
      du41=0
      du6=byr-byl
      du5=0.
      t1=betay*du6+betaz*du5
      t2=(prr-prl+(byave*du6)*pi4i
     1  +(gm-2.0d0)*(bybar*du6)*pi4i)/(gm-1.0d0)
      t3=betaz*du6-betay*du5
      s1=du21
      s2=betay*du31+betaz*du41
      s3=betaz*du31-betay*du41
      p11=alphs*cfast*sqrt(pi4/rhobar)
      p12=-alphf*cs2/cfast*sqrt(pi4/rhobar)
      p21=alphf*(cfast2-cs2*(gm-2.0d0)/(gm-1.0d0))
      p22=alphs*(cslow2-cs2*(gm-2.0d0)/(gm-1.0d0))
      q11=alphf*cfast
      q12=alphs*cslow
      q21=-alphs*ca*sgnbx
      q22=alphf*cs*sgnbx
      detp=p11*p22-p12*p21
      detq=q11*q22-q12*q21
c     chkdp=cs2*cfast/(gm-1.0d0)*sqrt(pi4/rhobar)
c     chkdq=cfast*cs*sgnbx
      wpf=0.5d0*((p22*t1-p12*t2)/detp+(q22*s1-q12*s2)/detq)
      wmf=0.5d0*((p22*t1-p12*t2)/detp-(q22*s1-q12*s2)/detq)
      wps=0.5d0*((-p21*t1+p11*t2)/detp+(-q21*s1+q11*s2)/detq)
      wms=0.5d0*((-p21*t1+p11*t2)/detp-(-q21*s1+q11*s2)/detq)
      wpa=0.5d0*(sqrt(rhobar*pi4i)*t3-sgnbx*s3)
      wma=0.5d0*(sqrt(rhobar*pi4i)*t3+sgnbx*s3)
      wze=drho-alphf*(wpf+wmf)-alphs*(wps+wms)
c----- flux
      fluxlro=rhol*vxl
      fluxlrx=rhol*vxl*vxl+prl+(-bxbar**2+byl**2)*pi8i
      fluxlry=rhol*vxl*vyl-bxbar*byl*pi4i
      fluxlby=vxl*byl-vyl*bxbar
      fluxlee=rhol*vxl*hl-bxbar*(bxbar*vxl+byl*vyl)*pi4i

      fluxrro=rhor*vxr
      fluxrrx=rhor*vxr*vxr+prr+(-bxbar**2+byr**2)*pi8i
      fluxrry=rhor*vxr*vyr-bxbar*byr*pi4i
      fluxrby=vxr*byr-vyr*bxbar
      fluxree=rhor*vxr*hr-bxbar*(bxbar*vxr+byr*vyr)*pi4i
c----- components of the eigen vectors
      rpfro=alphf
      rpfrx=alphf*(vxbar+cfast)
      rpfry=alphf*vybar-alphs*betay*ca*sgnbx
      rpfby=alphs*betay*cfast*sqrt(pi4/rhobar)
      rpfee=alphf*(0.5d0*(vxbar**2+vybar**2)
     1           +delb2+cfast*vxbar+cfast2/(gm-1.0d0)
     2           +(cfast2-cs2)*(gm-2.0d0)/(gm-1.0d0))
     3    -alphs*ca*(betay*vybar)*sgnbx
      rmfro=alphf
      rmfrx=alphf*(vxbar-cfast)
      rmfry=alphf*vybar+alphs*betay*ca*sgnbx
      rmfby=rpfby
      rmfee=alphf*(0.5d0*(vxbar**2+vybar**2)
     1           +delb2-cfast*vxbar +cfast2/(gm-1.0d0)
     2           +(cfast2-cs2)*(gm-2.0d0)/(gm-1.0d0))
     3    +alphs*ca*(betay*vybar)*sgnbx
      rpsro=alphs
      rpsrx=alphs*(vxbar+cslow)
      rpsry=alphs*vybar+cs*sgnbx*alphf*betay
      rpsby=-sqrt(pi4/rhobar)*cs2*alphf*betay/cfast
      rpsee=alphs*(0.5d0*(vxbar**2+vybar**2)
     1           +delb2+cslow*vxbar+cslow2/(gm-1.0d0)
     2           +(cslow2-cs2)*(gm-2.0d0)/(gm-1.0d0))
     3    +alphf*cs*(betay*vybar)*sgnbx
      rmsro=alphs
      rmsrx=alphs*(vxbar-cslow)
      rmsry=alphs*vybar-cs*sgnbx*alphf*betay
      rmsby=rpsby
      rmsee=alphs*(0.5d0*(vxbar**2+vybar**2)
     1           +delb2-cslow*vxbar+cslow2/(gm-1.0d0)
     2           +(cslow2-cs2)*(gm-2.0d0)/(gm-1.0d0))
     3    -alphf*cs*(betay*vybar)*sgnbx
      rparo=0.0d0
      rparx=0.0d0
      rpary=-sgnbx*betaz
      rpaby=sqrt(pi4/rhobar)*betaz
      rpaee=-(betaz*vybar)*sgnbx
      rmaro=0.0d0
      rmarx=0.0d0
      rmary=-rpary
      rmaby=rpaby
      rmaee=-rpaee
      rzero=1.0d0
      rzerx=vxbar
      rzery=vybar
      rzeby=0.0d0
      rzeee=0.5d0*(vxbar**2+vybar**2)+delb2

c-----computation of f(i+1/2,j)
      fro(i)=0.5d0*(fluxlro+fluxrro
     1       +elpf*wpf*rpfro +elmf*wmf*rmfro 
     &       +elps*wps*rpsro +elms*wms*rmsro
     2       +elpa*wpa*rparo +elma*wma*rmaro +elze*wze*rzero)

      fee(i)=0.5d0*(fluxlee+fluxree
     1       +elpf*wpf*rpfee +elmf*wmf*rmfee 
     &       +elps*wps*rpsee +elms*wms*rmsee
     2       +elpa*wpa*rpaee +elma*wma*rmaee +elze*wze*rzeee)

      frx(i)=0.5d0*(fluxlrx+fluxrrx
     1       +elpf*wpf*rpfrx +elmf*wmf*rmfrx 
     &       +elps*wps*rpsrx +elms*wms*rmsrx
     2       +elpa*wpa*rparx +elma*wma*rmarx +elze*wze*rzerx)

      fry(i)=0.5d0*(fluxlry+fluxrry
     1       +elpf*wpf*rpfry +elmf*wmf*rmfry 
     &       +elps*wps*rpsry +elms*wms*rmsry
     2       +elpa*wpa*rpary +elma*wma*rmary +elze*wze*rzery)

      fby(i)=0.5d0*(fluxlby+fluxrby
     1       +elpf*wpf*rpfby +elmf*wmf*rmfby 
     &       +elps*wps*rpsby +elms*wms*rmsby
     2       +elpa*wpa*rpaby +elma*wma*rmaby +elze*wze*rzeby)

      enddo

      return
      end
c======================================================================|
      subroutine scrdy(dsc,dscm,sc,scm,dx,dxm,ix)
c======================================================================|
c
c NAME  scrrdy
c
c PURPOSE
c    calculate cross section derivatives
c
c OUTPUTS
c    dsc(ix), dscm(ix) : [double] cross section derivative
c
c INPUTS
c    sc(ix), scm(ix) : [double] cross section
c    dx(ix),dxm(ix): [double] grid spacing
c    ix: [integer] dimension size
c
c HISTORY
c    written 2002-3-1 T. Yokoyama
c
c----------------------------------------------------------------------|
      implicit double precision (a-h,o-z)
      dimension dsc(ix),dscm(ix)
      dimension sc(ix),scm(ix)
      dimension dx(ix),dxm(ix)
c----------------------------------------------------------------------|

      do i=2,ix-1
          dsc(i) = (scm(i) - scm(i-1))/dx(i)
      enddo

      do i=2,ix-2
        dscm(i)= (sc(i+1) - sc(i))/dxm(i)
      enddo

      return
      end
c======================================================================|
      subroutine tvdminmod(da,daw,ix)
c======================================================================|
c
c NAME  tvdminmod
c
c PURPOSE
c    Interporate the physical variables based on MUSCL
c    using 'min-mod' function as a limitter
c
c INPUTS & OUTPUTS
c    None
c
c OUTPUTS
c    daw(ix,2): [double] variable at cell boundary
c
c INPUTS
c    da(ix): [double] physical variable
c    ix: [integer] dimension size
c
c HISTORY
c    written 2002-3-1 T. Yokoyama based on N. Fukuda's code
c
c----------------------------------------------------------------------|
      implicit double precision (a-h,o-z)

      dimension da(ix)
      dimension daw(ix,2)
c----------------------------------------------------------------------|
c     define limiter functions
      flmt(a,b)=max(0.0d0,min(b*sign(1.0d0,a),abs(a)))*sign(1.0d0,a)
c----------------------------------------------------------------------|
      do i=2,ix-2
         daw(i,1)=da(i)+0.5*flmt(da(i+1)-da(i),da(i)-da(i-1))
         daw(i,2)=da(i+1)-0.5*flmt(da(i+1)-da(i),da(i+2)-da(i+1))
      enddo
      return
      end