program lager
c     Lagrange solver after von Neumann & Richtmyer (1950) or Noh (1987)
c     The grid is as in hydro2.f
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     Again, 1| is never used. 
c     Parcel speed and location is u(|) and  x(|)
c     Their density and mass  is   r(:) and dm(:)
c     History:
c        23 Aug 01   my first Lagrange solver ever
c        24 Aug 01   supply code added to solver; Riemann tube
c        25 Aug 01   double shock, CAK wind
      implicit NONE
      integer i,ig,id,ima,iou,pro,ble,bri
      parameter (ig=200)
      integer iau1,iau2
      real au1,au2,vi1,vi2
      real a,a2,dt,cfl,ti,tie,tio,xmi,xma
      real u(ig+4),x(ig+4),r(ig+4),dm(ig+4)
      open (5,file='data.hyd3',status='unknown')
c     -----------------------------------------------------
c     MAIN CODE
      pro=7
      call setup    (id,pro,xmi,xma,ima,
     $     a,a2,cfl,ti,tie,tio,iou,ble,bri,vi1,vi2)
      call initial1 (ig,ima,pro,a,u,r,vi1,vi2)
      call initial2 (ig,id,pro,ima,xmi,xma,u,x,r,dm,vi1,vi2)
      call output   (ig,iou,ti,tie,tio,u,x,r,dm,vi1,vi2)
c*    time loop
      do i=1,1e8
      call clock    (ig,u,x,a,cfl,dt,ti,vi1,vi2)
      call lagrange (ig,id,pro,dt,a2,u,x,r,dm,vi1,vi2)
      call boundary (ig,a,r,u,ble,bri,vi1,vi2)
      call output   (ig,iou,ti,tie,tio,u,x,r,dm,vi1,vi2)
      enddo
      stop
      end

c     -----------------------------------------------------
      subroutine clock (ig,u,x,a,cfl,dt,ti,vi1,vi2)
c     Courant time step
      implicit NONE
      integer i,ig
      integer iau1,iau2
      real au1,au2,vi1,vi2
      real a,cfl,dt,ti
      real u(ig+4),x(ig+4)
      if (cfl.lt.0.) then
         dt=-cfl
      else
         dt=1.e20
         do i=3,ig+2
         dt=min(dt,cfl*(x(i+1)-x(i))/(a+max(abs(u(i+1)),abs(u(i)))))
         enddo
      endif
      ti=ti+dt
      return
      end

c     -----------------------------------------------------
      subroutine setup (id,pro,xmi,xma,ima,
     $     a,a2,cfl,ti,tie,tio,iou,ble,bri,vi1,vi2)
c     fixes problem parameters
      implicit NONE
      integer id,ima,pro,iou,ble,bri
      integer iau1,iau2
      real au1,au2,vi1,vi2
      real a,a2,cfl,xmi,xma,ti,tie,tio
      if (pro.eq.1) then ! Riemann shock tube
         a=1.
         ble=1
         bri=2
         cfl=0.5
         id=1
         tie=0.2
         ima=1
         xmi=0.
         xma=1.
         iou=50
      else if (pro.eq.2) then ! Riemann double shock
         a=0.3
         ble=2
         bri=2
         cfl=0.1
         id=1
         tie=0.45
         ima=1
         xmi=0.
         xma=0.5
         iou=50
      else if (pro.eq.7) then ! stationary CAK wind
c        Model wind from accretion disk. Normalization is (ru)_cri=1
c        vi1=2*sqrt(g_c*r_c*u_c) is force proportionality constant
c        vi2 is mass loss rate in units of critical mass loss
         vi1=2.*sqrt(2./3./sqrt(3.)) ! never change. Note r_c*u_c=1
         vi2=.8 ! this picks the stationary solution wanted
         a=0.2
         ble=1
         bri=2
         cfl=0.1
         id=1
         tie=30.
         ima=2
         xmi=-3.
         xma=3.
         iou=50
      endif ! different problems
c     Problem independent settings
      ti=0.
      tio=0.
      a2=a*a
      return
      end

c     ---------------------------------------
      subroutine initial1 (ig,ima,pro,a,u,r,vi1,vi2)
c     Fix initial fields u and r. x and dm are done in initial2
      implicit NONE
      integer i,ig,ima,pro
      integer iau1,iau2
      real au1,au2,vi1,vi2
      real a,u(ig+4),r(ig+4)
c     Riemann shock tube
      if (pro.eq.1) then
      do i=1,ig+4
      u(i)=0.
      r(i)=1.
      if (i.ge.0.4*ig+2) r(i)=0.1
      enddo
c     Double shock
      else if (pro.eq.2) then
      do i=1,3*ig/10+2
      r(i)=1.
      u(i)=1.
      enddo
      do i=3*ig/10+3,ig+4
      r(i)=1.
      u(i)=0.
      enddo
c     CAK accretion disk wind
      else if (pro.eq.7) then
      do i=1,ig+4
      u(i)=max(a/10.,2.*(i-0.5*ig-2)/float(ig))
      r(i)=vi2/u(i)
      enddo
      endif
      return
      end

c     ---------------------------------------
      subroutine initial2 (ig,id,pro,ima,xmi,xma,u,x,r,dm,vi1,vi2)
c     Fixes intial dm and r fields, eq (3.1) in Noh (1987). Supply density
      implicit NONE
      integer i,ig,id,ima,pro
      integer iau1,iau2
      real au1,au2,vi1,vi2
      real xmi,xma,dmm
      real u(ig+4),x(ig+4),r(ig+4),dm(ig+4)
      if (ima.eq.1) then 
c     Assuming constant parcel mass, parcel location is derived. First, the
c     total mass in [xmi,xma] is calculated assuming that initially,
c     r(i)=r(x(i)), with linear x(i). Later, r(i)=r(m(i)) is re-defined.
      dmm=0.
      au2=xmi
      do i=2,ig+2 ! any large number does for mass integral
      au1=au2
      au2=au1+(xma-xmi)/float(ig)
      dmm=dmm+r(i)*(au2**id-au1**id)
      enddo
c     Initializing parcels with this mass
      dmm=dmm/float(ig) ! now ig counts
      do i=2,ig+3
      dm(i)=dmm
      enddo
c     Calculating parcel boundaries
      au1=1./float(id)
      x(3)=xmi
      x(2)=(x(3)**id-dm(2)/r(2))**au1
      do i=4,ig+3
      x(i)=(x(i-1)**id+dm(i-1)/r(i-1))**au1
c     mapping r(x(i)) -> r(m(i)): don't separate this off into an extra loop!
c     Note that x=x(r(x))
      if (pro.eq.1) then
      r(i)=1.
      if (x(i).ge.0.6*xmi+0.4*xma) r(i)=0.1
      else if (pro.eq.2) then
      u(i)=1.
      if (x(i).ge.0.7*xmi+0.3*xma) u(i)=0.
      else if (pro.eq.7) then
      endif
      enddo
c**   2nd method: parcel locations (equidist, etc) given, masses derived
      else if (ima.eq.2) then
      do i=2,ig+4
      x(i)=xmi+(xma-xmi)*(i-3.)/float(ig)
      enddo
      do i=2,ig+3
      dm(i)=r(i)*(x(i+1)**id-x(i)**id)
      enddo
      endif ! ima
c     now gravity...
      return
      end

c     -----------------------------------------------------
      subroutine output (ig,iou,ti,tie,tio,u,x,r,dm,vi1,vi2)
c     writes output to file
      implicit NONE
      integer i,ig,it,iou
      integer iau1,iau2
      real au1,au2,vi1,vi2
      real ti,tie,tio
      real u(ig+4),x(ig+4),r(ig+4),dm(ig+4)
      if (ti.lt.tio) return
      tio=tio+tie/iou
      do i=3,ig+3
      write (5,*) x(i),u(i),r(i)
      enddo
      if (ti.ge.tie) then
         close (5)
         stop
      endif
      return
      end

c     -----------------------------------------------------
      subroutine boundary (ig,a,r,u,ble,bri,vi1,vi2)
c     boundary conditions for u and r.
c     Simple upwind boundary conditions: ble=1/bri=1 for subsonic
c     inflow/outflow; ble=2/bri=2 for supersonic inflow/outflow
      implicit NONE
      integer i,ig,ble,bri
      integer iau1,iau2
      real au1,au2,vi1,vi2
      real a,r(ig+4),u(ig+4)
c     left boundary conditions
      if (ble.eq.1) then      ! subsonic inflow: fix 1, extrapolate 1
         r(2)=r(2)
         u(3)=u(4)
      else if (ble.eq.2) then ! supersonic inflow: fix 2
         r(2)=r(2)
         u(3)=u(3)
      endif
c     right boundary conditions
      if (1.eq.2) then ! automatic detection
      if (u(ig+3).gt.a) then
         bri=2
      else
         bri=1
      endif
      endif ! 1.eq.2
      if (bri.eq.1) then      ! subsonic outflow: fix 1, extrapolate 1
         r(ig+3)=r(ig+2)
         u(ig+3)=u(ig+3)
      else if (bri.eq.2) then ! supersonic outflow: extrapolate 2
         r(ig+3)=r(ig+2)
         u(ig+3)=u(ig+2)
      endif
c     symmetry reductions
      r(1)=r(2)
      u(2)=u(3)
      r(ig+4)=r(ig+3)
      u(ig+4)=u(ig+3)
      return
      end

c     ---------------------------------------
      subroutine lagrange (ig,id,pro,dt,a2,u,x,r,dm,vi1,vi2)
c     Lagrange solver, after p. 295-296 in Richtmyer & Morton
c     and p. 85 in Noh (1987)
      implicit NONE
      integer i,ig,id,pro
      integer iau1,iau2
      real au1,au2,vi1,vi2
      real dt,a2,C02
      real u(ig+4),x(ig+4),r(ig+4),dm(ig+4),p(ig+4)
c     step 0: define pressure which includes artificial pressure
      C02=3.
      do i=3,ig+2
      p(i)=r(i)*(a2 + C02*min(0.,u(i+1)-u(i))**2)
      enddo
c     step 1a: Euler pressure equation, interior mesh
      do i=4,ig+2
      u(i)=u(i)-id*dt*x(i)**(id-1)*(p(i)-p(i-1))/(dm(i)+dm(i-1))*2.
      enddo
c--------------------------
c     step 1b: other forces
      if (pro.eq.7) then
      do i=4,ig+2
      if (x(i).ge.0.) u(i)=u(i) -dt*x(i)/(1.+x(i)**2)**1.5
     $  +dt*vi1*sqrt(abs(u(i+1)-u(i-1))/(x(i+1)-x(i-1))*
     $             2./(r(i-1)+r(i)))
      enddo
      endif
c--------------------------
c     step 2: boundary conditions
      u(3)=   u(4)    ! CAUTION
      u(ig+3)=u(ig+2) ! CAUTION
c     step 3: kinematics including boundaries
      do i=3,ig+3
      x(i)=x(i)+dt*u(i)
      enddo
c     step 4: density, interior mesh
      do i=3,ig+2
      r(i)=dm(i)/(x(i+1)**id-x(i)**id)
      enddo
      return
      end