      program mflu    
c     (1) Simple, time-explicit, staggered mesh, 1-D hydro code for ideal,
c         compressible, isothermal multicomponent flows using van Leer (1977) advection.
c     (2) The code is as short as possible, at the price of generality.
c     (3) To avoid introducing many control parameters, path control is
c         done via "if (1.eq.1)" and "if (1.eq.2)". Check!
c     (4) Two "viruses" vi1,vi2 infect all subroutines: they serve as
c         global parameters, to avoid common blocks. 
c     (5) au1,au2 are auxiliary variables that frequently change meaning.
c     (6) All arrays run from 1 to ig+4 (with 1| never used). 
c     (7) Basic fields are density (r) and momentum density (ru), where r
c         is zone-centered (symbol :) and ru is face-centered (symbol |).
c     (8) The mesh indices are chosen in accordance with Zeus 2-D. "I" are
c         interior points, "B" are true boundaries, "G" are ghost cells:
c
c     1  1  2  2  3  3  4... i  i i+1 ...ig+2 ig+2 ig+3 ig+3 ig+4 ig+4
c     |  :  |  :  |  :  |    |  :  |       |    :    |    :    |    :
c        G  G  B  B  I  I    I  I  I       I    I    B    B    G    G
c
c
c     (p) Viktor Votruba 2003-2005
c     2003 Major change of one-component code flu.f to two-component fluid
c     2004 Fixed problems with courant time step, and corrected problems
c          with procedure for radiative acceleration.  
c     2005 Developed algorithm for calculation of friction term
c     2006 I found strong decoupling instability problem

      implicit NONE
      integer i,ig,coo,pro,iou,ble,blep,bri,brip,iau1,iau2,grb
      parameter (ig=1600)
      double precision a,ap,a2,ap2,xmi,xma,dt,ti,tio,tie,ddt,
     & cfl,rc,au1,au2,vi1,vi2,vth,radc,abu
      double precision tiop,tiep
      double precision r(ig+4),ru(ig+4),ua(ig+4),ub(ig+4),
     & xa(ig+4),xb(ig+4),dfa(ig+4),dfb(ig+4),dva(ig+4),
     & dvb(ig+4),g(ig+4),coss(ig+4)
      double precision rgl(ig+4),kp(ig+4),rpi(ig+4)
      double precision rp(ig+4),rup(ig+4),uap(ig+4),ubp(ig+4),
     & gp(ig+4),cosse(ig+4)
      double precision rn,mn,sn,lum
      open (1,file='inp.dat',status='unknown')
      open (2,file='out.dat',status='unknown')
      open (4,file='flua.dat',status='unknown')
      open (5,file='flub.dat',status='unknown')    
      open (9,file='scratch.dat',status='unknown')
      open (10,file='header',status='unknown')
      pro=8 ! setup
      call norm  (rn,mn,sn,vth,radc)
      call set1  (ig,coo,pro,xmi,xma,ua,uap,ub,ubp,
     & 	 xa,xb,g,a,a2,ap,ap2,cfl,iou,ble,blep,bri,
     &   brip,ddt,ti,tio,tie,vi1,vi2,rn,mn,sn,radc,grb)
      call grid  (ig,coo,pro,xmi,xma,xa,xb,dfa,dfb,
     &            dva,dvb,vi1,vi2,grb)
      call set2  (ig,coo,pro,xa,g,coss,vi1,vi2)
      call init  (ig,coo,pro,r,rp,ru,rup,xa,xb,xmi,xma,
     &            a,ap,vi1,vi2,rn,mn,sn)
C      call stat1 (ig,r,ru,rp,rup,xb,1,rn,mn,sn) 
      call out   (ig,pro,r,ru,rp,rup,xa,xb,iou,ddt,ti,
     &             tio,tie,vi1,vi2,rgl,rpi,rn,mn,sn)
      do ! time loop
      call speed (ig,r,ru,ua,ub,vi1,vi2,rn,mn,sn)
      call speed (ig,rp,rup,uap,ubp,vi1,vi2,rn,mn,sn)
      call clock (ig,ub,ubp,r,rp,xa,xb,a,ap,cfl,dt,
     &            ti,vi1,vi2,rn,mn,sn)
      call friction (ig,r,rp,ru,rup,ua,uap,ub,ubp,xa,xb,
     &      a,ap,dt,kp,rn,mn,sn,rgl,rpi,vi1,vi2,vth,abu)
      call leer1 (ig,ua,r,xa,xb,dfa,dvb,dt,vi1,vi2)
      call leer1 (ig,uap,rp,xa,xb,dfa,dvb,dt,vi1,vi2)
      call leer2 (ig,ub,ru,xa,xb,dfb,dva,dt,vi1,vi2)
      call leer2 (ig,ubp,rup,xa,xb,dfb,dva,dt,vi1,vi2)
      call radia (ig,coo,pro,r,ru,rp,ua,ub,xa,xb,dt,ti,coss,
     &            vi1,vi2,rn,mn,sn,rgl,radc,abu)
      call press (ig,r,ru,xb,a2,dt,vi1,vi2)
      call press (ig,rp,rup,xb,ap2,dt,vi1,vi2)
      call grav  (ig,r,ru,g,dt,vi1,vi2,rn,mn,sn,lum)
      call grav  (ig,rp,rup,g,dt,vi1,vi2,rn,mn,sn,lum)

c      call abbo  (ig,pro,r,ru,dt,ti,vi1,vi2)
      call visco (ig,r,ru,dt,xa,xb,ua,vi1,vi2)    
      call visco (ig,rp,rup,dt,xa,xb,ua,vi1,vi2)
      call bound (ig,r,ru,rp,rup,ble,blep,bri,brip,a,ap,
     &            kp,rgl,vi1,vi2,rn,mn,sn)
      call out   (ig,pro,r,ru,rp,rup,xa,xb,iou,ddt,ti,tio,tie,
     &            vi1,vi2,rgl,rpi,rn,mn,sn)
      enddo
      stop
      end
      
c     ----------------------------------------------------        
c     Modification: normalization routine using relation
c     GM=1, Msdot/4*Pi=1
c     parameter M_m,R_m,M_dot from Owocki (1984)
c     (12.5. 2004)

      subroutine norm (rn,mn,sn,vth,radc)
      double precision rn,mn,sn
      double precision rs,ms,r_m,m_m,ms_dot,m_dot,ls,lum,T,vth
      double precision pi,ge,c,se,sig,nu,kb,radc,mp,kappa
      double precision alpha,k,delta
      double precision buf1,buf2,D
      parameter (rs=6.96d10,ms=1.99d33,ls=3.9d33,sig=5.67051d-5,
     & pi=3.141592653589,kb=1.38d-16,mp=1.67d-24,ai=12)
      parameter (T=28500,r_m=37.d0,m_m=90.d0)
      parameter (alpha=0.5,k=0.270,delta=0.03)
      parameter (kappa=6.6726d-8)
      rn=r_m*rs
      se=0.346d0
      c=2.99d10
      lum=4*pi*rn**2*sig*T**4
C     Eddington parameter Ge
      Ge=se*lum/(4*pi*c*kappa*m_m*ms)
      sn=sqrt(rn**3/((1.d0-Ge)*m_m*ms*kappa))
      nu=2/(1.+3*0.73+0.5*0.16)
      a=sqrt(kb*T/(nu*mp))*sn/rn
      vth=sqrt(2*kb*T/(ai*mp))*sn/rn
      buf1=((1.-Ge)*m_m*ms*kappa)**(1.-alpha)*(1.36d7*(r_m*rs)**2)
     & **delta*(1.-alpha)**(alpha-1.)*alpha**(-alpha)*(alpha+1)*
     & 0.9*(1.-0.1**(alpha+1.)) 
      D=1.2/(1.4*mp*1d11)
      ms_dot=(buf1*(2.*pi/D)**delta*(se*vth*rn/sn/(4.*pi))**alpha/(k*Ge*
     & m_m*ms*kappa*2**(-delta/2.)))**(1./(delta-alpha))
c      ms_dot=4*pi/(se*vth*rn/sn)*(se/(4.*pi))**(1./alpha)*
c     & ((1.-alpha)/alpha)**((1.-alpha)/alpha)*(k*alpha)**(1./alpha)*
c     & (lum/c)**(1./alpha)*(kappa*m_m*ms*(1.-Ge))**((alpha-1.)/alpha)     
c      ms_dot=m_dot*ms/(365.25d0*24.d0*3600.d0)
      mn=ms_dot/(4.d0*pi)*sn
      se=se*mn/rn**2
      radc=lum/rn**2*sn**3/mn*se/(4.*pi*c*sn/rn)*k/
     & (vth**alpha*se**alpha) 
      write(10,*)'Parametry hvezdy: '
      write(10,*)'Polomer [R*]:',r_m,'Hmotnost [M*]: ',
     & m_m,' Teplota:',T
      write(10,*)'Svitivost [L*]:',lum/ls,' Ztrata hmoty [M*]:',
     & ms_dot*365.25d0*24.d0*3600/ms,' Ge:',Ge
      write(10,*) 
      write(10,*)'Normovaci konstanty',rn,mn,sn
      write(10,*)'Rn: ',rn,' Mn:',mn,' Sn: ',sn
      write(10,*)'a: ',a*rn/sn,' vth',vth*rn/sn
      write(10,*)'alpha: ',alpha,' delta: ',delta,' k: ',k
      write(10,*)'Radiative scalling constant',radc
      close(10);
      return
      end
c     -----------------------------------------------------
      subroutine leer1 (ig,ua,r,xa,xb,dfa,dvb,dt,vi1,vi2)
c     van leer (1977) advection for density r (b grid)
      implicit NONE
      integer i,ii,ig
      integer iau1,iau2
      double precision au1,au2,vi1,vi2
      double precision dd,dm,dp,vl,dt
      double precision ua(ig+4),r(ig+4),xa(ig+4),xb(ig+4),
     & dfa(ig+4),dvb(ig+4),ri(ig+4)
c     center index of upwinded 3-point molecule 
      do i=3,ig+3
      if (ua(i).gt.0.d0) then
         ii=i-1
      else
         ii=i
      endif
c     van Leer derivative
      dm=(r(ii)-r(ii-1))/(xb(ii)-xb(ii-1))
      dp=(r(ii+1)-r(ii))/(xb(ii+1)-xb(ii))
      dd=dm*dp
      vl=0.d0
      if (dd.gt.0.d0) vl=2.d0*dd/(dm+dp)
c     interface interpolants (absolute coordinates required)
      ri(i)=r(ii)+vl*(xa(i)-xb(ii)-ua(i)*dt/2.d0)
      enddo
c     zone loop for advection
      do i=3,ig+2
      r(i)=r(i)-dt/dvb(i)*
     &        (ri(i+1)*ua(i+1)*dfa(i+1)-ri(i)*ua(i)*dfa(i))
      enddo
      return
      end

c     -----------------------------------------------------
      subroutine leer2 (ig,ub,ru,xa,xb,dfb,dva,dt,vi1,vi2)
c     van Leer (1977) advection for momentum density ru (a grid)
      implicit NONE
      integer i,ii,ig
      integer iau1,iau2
      double precision au1,au2,vi1,vi2
      double precision dd,dm,dp,vl,dt
      double precision ub(ig+4),ru(ig+4),xa(ig+4),xb(ig+4),
     & dfb(ig+4),dva(ig+4),rui(ig+4)
c     center index of upwinded 3-point molecule 
      do i=3,ig+2
      if (ub(i).gt.0.d0) then
         ii=i
      else
         ii=i+1
      endif
c     van Leer derivative
      dm=(ru(ii)-ru(ii-1))/(xa(ii)-xa(ii-1))
      dp=(ru(ii+1)-ru(ii))/(xa(ii+1)-xa(ii))
      dd=dm*dp
      vl=0.d0
      if (dd.gt.0.d0) vl=2.d0*dd/(dm+dp)
c     interface interpolants
      rui(i)=ru(ii)+vl*(xb(i)-xa(ii)-ub(i)*dt/2.d0)
      enddo
c     zone loop for advection
      do i=4,ig+2
      ru(i)=ru(i)-dt/dva(i)*
     &        (rui(i)*ub(i)*dfb(i)-rui(i-1)*ub(i-1)*dfb(i-1))
      enddo
      return
      end

c     -----------------------------------------------------
      subroutine press (ig,r,ru,xb,a2,dt,vi1,vi2)
c     pressure force
      implicit NONE
      integer i,ig
      integer iau1,iau2
      double precision au1,au2,vi1,vi2
      double precision dt,a2
      double precision r(ig+4),ru(ig+4),xb(ig+4)
      do i=4,ig+2
      ru(i)=ru(i)-a2*dt*(r(i)-r(i-1))/(xb(i)-xb(i-1))
      enddo
      return
      end

c     -----------------------------------------------------
      subroutine grav (ig,r,ru,g,dt,vi1,vi2,rn,mn,sn,lum)
c     gravitational acceleration
      implicit NONE
      integer i,ig
      integer iau1,iau2
      double precision au1,au2,vi1,vi2
      double precision dt
      double precision r(ig+4),ru(ig+4),g(ig+4)
      double precision rn,mn,sn,lum
      do i=4,ig+2
      ru(i)=ru(i)-dt*(r(i-1)+r(i))/2.d0*(g(i))
      enddo
      return
      end

c     ------------------------------------------------------
c     Modification of routine for computing radiative force
c     CAK (1975). We use dimensionless form with scaleable factor 
C     for interactions with absorbing ions only
      subroutine radia
     &   (ig,coo,pro,r,ru,rp,ua,ub,xa,xb,dt,ti,coss,vi1,vi2,
     &   rn,mn,sn,rgl,radc,abu)
      implicit NONE
      integer i,ig,coo,pro,ic,iau1,iau2
      double precision au1,au2,vi1,vi2,dt,ti,rc,pi,radc
      double precision r(ig+4),ru(ig+4),ua(ig+4),ub(ig+4),
     & xa(ig+4),xb(ig+4),rp(ig+4),coss(ig+4),rgl(ig+4)
      double precision rn,mn,sn
      double precision alp,del,k
      double precision f,w,mu2,pom,buff
      double precision ne,abu,vth,mp    
      parameter (alp=0.50,del=0.03,k=0.270)
      parameter (mp=1.672d-24)
      parameter (pi=3.14159265358979)
c     Sobolev force in flow direction ("point star approx")
c     For this time i assume alpha = 0.5
      buff=3.45228215
      do i=4,ig+2      
c     For simulation of pulsating shell (Porter, Skouza 1999)
C      
c      if (xa(i).ge.1.3) then
c      buff=0.0
c      else
      buff=radc
c      endif       
      rgl(i)=buff*sqrt(abu)/abu*(2/(r(i-1)+r(i))* !for alpha=0.5
     &       abs((ub(i)-ub(i-1))/(xb(i)-xb(i-1))))**alp
      enddo
c     flux drop-off for spherical symmetry
      do i=4,ig+2
      rgl(i)=rgl(i)*0.5*(r(i-1)+r(i))/xa(i)**2
      enddo
      do i=4,ig+2
c     we add correction factor [au2] Abbot (1980) 
      ne=(rp(i)*mn/rn**3)/(mp)*1.0d-11 
      mu2=coss(i)
      w=0.5d0*(1.-sqrt(mu2))
      pom=(ne/w)**del
c      write(*,*)((rp(i)*mn/rn**3)/(mp))/w
c     finite cone correction factor after CAK eq. [au1] (50)
      au1=abs((xa(i)/ua(i)*(ub(i)-ub(i-1))/(xb(i)-xb(i-1))))-1.d0 ! sigma
      f=((1.d0+au1)**(alp+1.d0)-(1.d0+au1*mu2)**(alp+1))/
     & ((1.d0-mu2)*(alp+1.d0)*au1*(1.d0+au1)**alp)
      rgl(i)=rgl(i)*f
      enddo
c     finally, acceleration
      do i=4,ig+2      
      ru(i)=ru(i)+dt*rgl(i)
      enddo
      return
      end
      
      
c     New routine for computing the friction force due to
c     Coulombic interactions, for technical purposes
c     we compute dimensionless form at this routine      
c     March 2005: add some changes in normalization
C     March 3: Chandrasekhar function was numericaly wrong
c     March 4: add exponential+linear grid, default turn off
c     April 20: Complete new routine, with semi implicit method
c     for calculation of friction
c     !Some danger with exponential function in part (1) of curve
      subroutine friction (ig,r,rp,ru,rup,ua,uap,ub,ubp,
     & xa,xb,a,ap,dt,kp,rn,mn,sn,rgl,rpi,vi1,vi2,vth,abu)
      implicit NONE
      integer i,j,ig,iau
      double precision r(ig+4),rp(ig+4),ru(ig+4),rup(ig+4),
     & ua(ig+4),uap(ig+4)
      double precision ub(ig+4),ubp(ig+4),xa(ig+4),xb(ig+4)
      double precision rgl(ig+4)
      double precision a,ap,vi1,vi2,dt
      double precision rn,mn,sn
      double precision rip(ig+4),rpi(ig+4),kp(ig+4)
      double precision vth,T,L,qp,qi,k,ts
      double precision ai,mp,e,rd
      double precision au1,au2,au3,au4,au5,au6,au7,au8
      double precision xpio(ig+4),xpin(ig+4)
      double precision g,buff,pi,abu
      double precision a1,b1,c1
      integer sng
      logical result
      parameter (pi=3.14159265358979,T=28500,k=1.3806e-16 )
      parameter (qp=1, qi=3, ai=12, mp=1.672e-24, e=4.8032e-10 )
      parameter (a1=-3.127,b1=6.5653,c1=-3.8439)
      au1=sqrt(-(4*a1*c1-b1*b1))
      au7=vth*sqrt(1+ai)
      sng=1
      do i=4,ig+2
        au2= (ua(i)-uap(i))
        if ((ua(i)-uap(i)).le.0) then
      	au2=(uap(i)-ua(i))
	sng=-1
	endif
        xpio(i)=sqrt(1/(ai+1))*(au2)/vth
        au4=2.*a1*xpio(i)+b1-au1
        au5=2.*a1*xpio(i)+b1+au1
        rip(i)=(0.5*(rp(i-1)+rp(i))+1/ai*0.5*(r(i-1)+r(i)))
        L=log(24.*pi*sqrt(mp/(rip(i)*mn/rn**3))*((k*T)/
     & (4.*pi*e**2))**1.5)
        kp(i)=((4.*pi*(L)/(k*T)*(au2)/abs(au2)*((qp*qi)**2)*(e**2))*
     &   e**2)
        au3=-(0.5*(rp(i-1)+rp(i)+r(i)+r(i-1))*mn/rn**3)/mp*
     &  sn**2/rn/(mp*ai)*kp(i)/vth       
      if (abs(xpio(i)).lt.0.15) then
          au3=au3*2./(3.*sqrt(pi))
          if ((au3*dt).le.(-70.)) then
	  xpin(i)=1.d-35
	  else
      	  xpin(i)=xpio(i)*exp(au3*dt)
          endif	  
 
      elseif (xpio(i).gt.1.85) then
          au3=au3*0.5
      	  xpin(i)=((xpio(i)**3)+3.*au3*dt)**(0.333)
      else
	  au6=exp(au3*dt*au1)
	  xpin(i)=0.5/a1*((au5+au4*au6)/(au5-au4*au6)*au1-b1)
      endif
      enddo
      do i=4,ig+2
      buff=au7*0.5*(r(i-1)+r(i))*(rp(i-1)+rp(i))/
     & (r(i-1)+rp(i-1)+r(i)+rp(i))
      rpi(i)=buff*(xpio(i)-xpin(i))
      abu=0.0127
      ru(i)=ru(i)-sng*rpi(i)
      rup(i)=rup(i)+sng*rpi(i)
      enddo
      return      
      end
      
c     -----------------------------------------------------
      subroutine speed (ig,r,ru,ua,ub,vi1,vi2,rn,mn,sn)
c     linear interpolation of advection speeds onto a,b grids
      implicit NONE
      integer i,ig
      integer iau1,iau2
      double precision rn,mn,sn
      double precision au1,au2,vi1,vi2
      double precision r(ig+4),ru(ig+4),ua(ig+4),ub(ig+4)
c     For r-advection. Boundaries must be included here, (i) for
c     advection itself, (ii) to calculate ub below
      do i=3,ig+3
      ua(i)=2.d0*ru(i)/(r(i-1)+r(i))
      enddo
c     For ru-advection, 5 point interpolation is better than
c     3-point interpolation 0.5*(ru(i)+ru(i+1))/r(i)
      do i=3,ig+2
      ub(i)=(ua(i)+ua(i+1))/2.d0
      enddo
      return
      end

c     -----------------------------------------------------
c     Modification for two components fluid, we set time step
c     as smallest number from predicted time step for first
c     component dt a second component dttp from courant stability
c     condition (14.06.2004)
c     add density of the ions as input parameter due to abbot waves
      subroutine clock (ig,ub,ubp,r,rp,xa,xb,a,ap,cfl,dt,ti,
     & vi1,vi2,rn,mn,sn)
c     Courant time step. Negative Courant number == fixed time step
      implicit NONE
      integer i,ig
      integer iau1,iau2
      double precision au1,au2,vi1,vi2,vth
      double precision a,ap,cfl,dt,ti,dtt,dttp
      double precision ub(ig+4),ubp(ig+4),xa(ig+4),xb(ig+4),r(ig+4)
      double precision rp(ig+4),rip(ig+4),kp(ig+4)
      double precision abu,pi,k,mp,ai,e,qp,qi,T,L
      double precision rn,mn,sn
      if (cfl.lt.0.d0) then
         dt=-cfl
      else
         dt=1.d30	 
         do i=3,ig+2
	 dtt=cfl*(xa(i+1)-xa(i))/(abs(ub(i))+a) !without Abbot waves
	 dttp=cfl*(xa(i+1)-xa(i))/(abs(ubp(i))+ap) !change for passive plasma ubp(i) don't account
	 if (dttp.lt.dtt) dtt=dttp
         if (dtt.lt.dt) dt=dtt
         enddo
      endif
      ti=ti+dt
c      write(*,*)'dt:',dt,' ti:',ti
      return
      end

c     -----------------------------------------------------
      subroutine set1 (ig,coo,pro,xmi,xma,ua,uap,ub,ubp,
     & 		      xa,xb,g,a,a2,ap,ap2, cfl,iou,ble,blep,
     &       bri,brip,ddt,ti,tio,tie,vi1,vi2,rn,mn,sn,radc,grb)
c     fixes problem parameters
      implicit NONE
      integer i,ig,coo,pro,iou,ble,blep,bri,brip,grb
      integer iau1,iau2
      double precision au1,au2,vi1,vi2
      double precision a,a2,ap,ap2,cfl,xmi,xma,ddt,ti,tio,tie
      double precision ua(ig+4),uap(ig+4),ub(ig+4),ubp(ig+4),
     & 		       xa(ig+4),xb(ig+4),g(ig+4)
      double precision rn,mn,sn
      double precision mp,T,k,ai,radc
c     Spherical stellar wind. See pro.eq.7
      parameter (T=28500,mp=1.67d-24,k=1.38d-16,ai=12)
c     vi1 is force proportionality constant      
         vi1=2.
c     vi2 si mass loos in units of critical one
         vi2=1.0
         a=0.1
         ap=0.09999
         write(*,*)ap,a
         ble=1
	 blep=2
         bri=2
	 brip=2
         cfl=0.01
         coo=2
         tie=25.d0
         xmi=1.d0
         xma=5.d0
         iou=100
         grb=2
c     General settings
      do i=1,ig+4
      ua(i)=0.d0
      uap(i)=0.d0
      ub(i)=0.d0
      ubp(i)=0.d0
      enddo
      a2=a*a
      ap2=ap*ap
      tio=0.d0
      ti=0.d0
      ddt=tie/float(iou)
      return
      end

c     -----------------------------------------------------
      subroutine set2 (ig,coo,pro,xa,g,coss,vi1,vi2)
c     fixes problem parameters that depend on the spatial grid
      implicit NONE
      integer i,ig,coo,pro
      integer iau1,iau2
      double precision au1,au2,vi1,vi2
      double precision xa(ig+4),g(ig+4),coss(ig+4)
c     CAK: stellar wind
         do i=1,ig+4
         au1=1.d0/xa(i)**2
         g(i)=au1
         coss(i)=1.d0-au1 ! mu_core^2 for FD correction factor
         enddo
      return
      end

c     -----------------------------------------------------
      subroutine init (ig,coo,pro,r,rp,ru,rup,xa,xb,
     & 		       xmi,xma,a,ap,vi1,vi2,rn,mn,sn)
c     fixes initial conditions for r and ru
      implicit NONE
      integer i,ig,coo,pro
      integer iau1,iau2
      double precision au1,au2,au3,au4,au5,au6,vi1,vi2
      double precision xmi,xma,a,ap,abu
      double precision r(ig+4),ru(ig+4),rp(ig+4),rup(ig+4),
     & xa(ig+4),xb(ig+4),u(ig+4),up(ig+4)
      double precision rn,mn,sn
      parameter(abu=0.0127)
c     Spherical CAK stellar wind
          au3=0.
	  au4=0.
	  do i=1,ig+4
	  if (au3.lt.ap) then
	    au5=xb(i)
	    au3=ap/10.*exp((au5-xmi)/ap**2)
          else
	    au3=ap+2.*(xb(i)-au5)
	  endif
	  au4=au3+0.1*sn/rn
	  r(i)=0.01/au4/xb(i)**2
	  rp(i)=0.99/au3/xb(i)**2
	  ru(i)=r(i)*au4
	  rup(i)=rp(i)*au3
	  enddo
	  	  
          r(1)=r(2)
	  rp(1)=rp(2)
	  r(ig+4)=r(ig+3)
	  rp(ig+4)=rp(ig+3)
      	  ru(ig+4)=ru(ig+3)	 
	  rup(ig+4)=rup(ig+3)
          ru(2)=ru(3)
	  rup(2)=rup(3)	 
      return
      end
      
c     -----------------------------------------------------
      subroutine stat1 (ig,r,ru,rp,rup,xb,io,rn,mn,sn)
c     Reads or writes input data
c     March 2005: Fixes some bugs with mesh grid
      integer ig,io
      double precision r(ig+4),ru(ig+4)
      double precision rp(ig+4),rup(ig+4)
      double precision xb(ig+4),u(ig+4),up(ig+4)
      double precision abu,a,ap
      double precision rn,mn,sn
      parameter (abu=0.0127)
      if (io.eq.1) then
         do i=1,ig+4
      	 read(1,*)r(i),ru(i),rp(i),rup(i)
	 enddo
      else
         do i=1,ig+4
         write (2,*) r(i),ru(i),rp(i),rup(i)
         enddo
      endif
      return
      end

c     -----------------------------------------------------
      subroutine bound (ig,r,ru,rp,rup,ble,blep,bri,brip,a,ap,
     &                 kp,rgl,vi1,vi2,rn,mn,sn)
c     Simple boundary conditions for r and ru. If a characteristic
c     leaves the mesh, a quantity (r or ru) is extrapolated, if a
c     characteristic enters the mesh, a quantity (r or ru) is fixed
c     via boundary conditions. This is the Courant stability theorem.
c     NO Riemann invariants are calculated.
c     NO automatic detection of transonic changes.
c     ble=1, bri=1 for subsonic inflow, outflow; 
c     ble=2, bri=2 for supersonic inflow, outflow.
      implicit NONE
      integer i,ig,ble,blep,bri,brip
      integer iau1,iau2
      double precision a,ap,abu,mp,pi
      double precision au1,au2,vi1,vi2
      double precision r(ig+4),ru(ig+4),rp(ig+4),rup(ig+4)
      double precision kp(ig+4),rgl(ig+4)
      double precision rn,mn,sn
      parameter (abu=0.0127,mp=1.67d-24,pi=3.14159265358979)
c     inner or lower boundary conditions
      if (ble.eq.1) then
         r(2) =r(2)        ! bc 
         ru(3)=ru(4)       ! ex
      else if (ble.eq.2) then
         r(2) =r(2)        ! bc
         ru(3)=ru(3)       ! bc
      endif

      if (blep.eq.1) then
        rp(2)=rp(2)       !bc
	rup(3)=rup(4)     !ex
      else if (blep.eq.2) then
	rp(2)=rp(2)       !bc
	rup(3)=rup(3)     !bc
      endif
c     outer or upper boundary condition
      if (bri.eq.1) then
         r(ig+3) =r(ig+2)  ! ex
         ru(ig+3)=ru(ig+3) ! bc
      else if (bri.eq.2) then
         r(ig+3) =r(ig+2)  ! ex
         ru(ig+3)=ru(ig+2) ! ex  
      endif
      if (brip.eq.1) then
         rp(ig+3)=rp(ig+2)   ! ex
      	 rup(ig+3)=rup(ig+3) ! bc
      else if (brip.eq.2) then
         rp(ig+3)=rp(ig+2)   ! ex
         rup(ig+3)=rup(ig+2) ! ex
      endif	 	
        
c     ghost cell reduction
      r(1)=r(2)
      rp(1)=rp(2)
      rup(2)=rup(3)
      ru(2)=ru(3)
      rup(1)=rup(2)
      ru(1)=ru(2)
      r(ig+4)=r(ig+3)
      rp(ig+4)=rp(ig+3)
      ru(ig+4)=ru(ig+3)
      rup(ig+4)=rup(ig+3)
      return
      end

c     -----------------------------------------------------
      subroutine grid
     &  (ig,coo,pro,xmi,xma,xa,xb,dfa,dfb,dva,dvb,vi1,vi2,grb)
c     Spatial mesh. Scalars (r) reside at cell center (termed b),
c     vectors (ru) at cell interfaces (termed a)
      implicit NONE
      integer i,ig,coo,pro,grb
      integer iau1,iau2,buff1,buff2,buff3
      double precision au1,au2,vi1,vi2
      double precision xmi,xma,rd1,rd2
      double precision xa(ig+4),xb(ig+4),dfa(ig+4),dfb(ig+4),
     & dva(ig+4),dvb(ig+4)
      parameter(rd1=1.4,rd2=2.1)
      parameter(buff1=300,buff2=300,buff3=300)
      if (grb.eq.1) then
c      Linear grid
      do i=1,ig+4
      xa(i)=xmi+(xma-xmi)*(i-3.d0)/float(ig)
      enddo
       elseif (grb.eq.2) then
c      Exponencial grid and linear
       iau1=int(ig*9/10)
       do i=3,ig+1
       xa(ig+4-i)=xma-(xma-1)*(log(float(i-2))/log(float(ig)))
       enddo
       xa(1)=xma-(xma-1)*(log(float(ig+1))/log(float(ig)))
       xa(2)=xma-(xma-1)*(log(float(ig+2))/log(float(ig)))
       do i=1,(ig+4-iau1)
       xa(iau1+i)=xa(iau1)+(xma-xa(iau1))*i/float(ig+2-iau1)
       enddo
       else
c     Manually adaptive linear grid :-) for part of decoupling
      do i=1,buff1+3
      xa(i)=xmi+(rd1-xmi)*(i-3)/float(buff1)
      enddo
      do i=(buff1+4),(buff1+buff2+3)
      xa(i)=rd1+(rd2-rd1)*(i-buff1-3)/float(buff2)
      enddo
      do i=(buff2+buff1+4),(buff1+buff2+buff3+4)
      xa(i)=rd2+(xma-rd2)*(i-(buff1+buff2)-3)/float(buff3)
      enddo
      endif

      do i=1,ig+3
      xb(i)=(xa(i)+xa(i+1))/2.d0
      enddo
      xb(ig+4)=2.d0*xb(ig+3)-xb(ig+2)
c     control volumes and interface areas. Cartesian coordinates
      if (coo.eq.0) then
         do i=1,ig+4
         dfa(i)=1.d0
         dfb(i)=1.d0
         enddo
         if (pro.eq.4) then ! Laval nozzle area function
            do i=1,ig+4
            if (xa(i).le.0.5d0) then
               dfa(i)=1.d0+2.d0*(0.5d0-xa(i))**2
               dfb(i)=1.d0+2.d0*(0.5d0-xb(i))**2
            else
               dfa(i)=1.d0+6.d0*(0.5d0-xa(i))**2
               dfb(i)=1.d0+6.d0*(0.5d0-xb(i))**2
            endif
            enddo
         endif ! Laval
         do i=2,ig+4
         dva(i)=(xb(i)-xb(i-1))*dfa(i)
         enddo
         dva(1)=dva(2)
         do i=1,ig+3
         dvb(i)=(xa(i+1)-xa(i))*dfb(i)
         enddo
         dvb(ig+4)=dvb(ig+3)
      else if (coo.eq.2) then
c     spherical coordinates. A factor 4*pi is dropped, as it would cancel!
         do i=1,ig+4
         dfa(i)=xa(i)**2
         dfb(i)=xb(i)**2
         enddo
         do i=2,ig+4
         dva(i)=(xb(i)**3-xb(i-1)**3)/3.d0
         enddo
         dva(1)=dva(2)
         do i=1,ig+3
         dvb(i)=(xa(i+1)**3-xa(i)**3)/3.d0
         enddo
         dvb(ig+4)=dvb(ig+3)
      endif ! cartesian vs spherical
      return
      end

c     -----------------------------------------------------
      subroutine out (ig,pro,r,ru,rp,rup,xa,xb,iou,ddt,ti,
     & tio,tie,vi1,vi2,rgl,rpi,rn,mn,sn)
c     writes output to file
      implicit NONE
      integer i,ig,iou,pro,io
      integer iau1,iau2
      double precision au1,au2,vi1,vi2
      double precision dtodx,a,ddt,ti,tio,tie
      double precision r(ig+4),ru(ig+4),xa(ig+4),xb(ig+4),u(ig+4)
      double precision rp(ig+4),rup(ig+4),rgl(ig+4),rpi(ig+4)
      double precision rn,mn,sn
c     output control
c      write(*,*)tio,'  ',ddt,ti
      if (ti.lt.(tio)) return
      tio=tio+ddt
c     data writing to file
      do i=2,ig+3
      write (4,*)xb(i),0.5*(r(i-1)+r(i)),
     & 2.*ru(i)/(r(i-1)+r(i))*rn/sn/100000
      write (5,*)xb(i),0.5*(rp(i-1)+rp(i)),
     & 2.*rup(i)/(rp(i-1)+rp(i))*rn/sn/100000
      enddo
c     program termination
      if (ti.ge.tie) then
      call stat1 (ig,r,ru,rp,rup,xb,2,rn,mn,sn) ! output for input
c     close (1)
      close (2)
      close (5)
c      close (7)
c      close (8)
      close (9)
      stop
      endif
c     print to screen
      print*,tio
      return
      end
      
c     -----------------------------------------------------
      subroutine abbo (ig,pro,r,ru,dt,ti,vi1,vi2)
c     Abbott wave perturbations in line driven winds
      implicit NONE
      integer i,ig,pro
      integer iau1,iau2
      double precision au1,au2,vi1,vi2
      double precision r(ig+4),ru(ig+4),dt,ti,pi
      if (pro.ne.9) return
      pi=3.1415927
      iau1=450  ! location
      au1=1.0  ! amplitude
      au2=0.002  ! period
      ru(iau1)=ru(iau1)+r(iau1)*dt*au1*2.*pi/au2*cos(2.*pi*ti/au2)
      return
      end

c     -----------------------------------------------------
      subroutine visco (ig,r,ru,dt,xa,xb,ua,vi1,vi2)
c     Planar artificial pressure after von Neumann and Richtmyer.
c     See Roache, comput. fluid dynamics, p. 232; paper by Noh; book
c     by Richtmyer & Morton; paper by Schulz; paper by Norman & Winkler
      implicit NONE
      integer i,ig,iau1,iau2
      double precision au1,au2,vi1,vi2,vis,dt
      double precision r(ig+4),ru(ig+4),xa(ig+4),xb(ig+4),ua(ig+4),
     & q(ig+4)
c     artificial viscosity, centered on scalar mesh. Only for
c     compression (Rosenbluth).
      vis=0.2
      do i=3,ig+2
      q(i)=r(i)*vis*min(0.,ua(i+1)-ua(i))**2
      enddo
c     momentum equation
      do i=4,ig+2
      ru(i)=ru(i)-dt*(q(i)-q(i-1))/(xb(i)-xb(i-1))
      enddo
      return
      end
      
      

