SUBROUTINE VELSET
C     =================
C
C     Determination of the macroscopic velocity as a function of depth
C
C     Input:  
C
C     RSTAR   - stellar radius (in solar radii or in cm)
C     RMAX    - maximum radial extent (in stellar radii)
C     AMLOSS  - mass loss rate ( in solar masses per year)
C     VELMAX  - maximum velocity (= V_infinity) - in km/s
C     BETA    - beta exponent in the beta-law for velocity
C     NDRAD   - Number of layers
C     NRCORE  - Number of core rays
C
C
c      parameter (un=1.,two=2.)
      INCLUDE 'PARAMS.FOR'
      INCLUDE 'MODELP.FOR'
      INCLUDE 'WINCOM.FOR'
      dimension zz(mdepth),vel0(mdepth),rrel(mdepth),
c     *          dvel0(mdepth),vel1(mdepth),hstt(mdepth),
     *          den0(mdepth),vel00(mdepth),ind(mdepth),
     *          densa(mdepth),eleca(mdepth),tempa(mdepth),
     *          rda(mdepth),rrela(mdepth),vel0a(mdepth)
c
      un=1
      two=2.
      read(55,*,err=100,end=100) rstar,rmax,amloss,vinf,beta,
     *           ndrad,nrcore,nfiry,ndf,nda
      rstr=rstar
      if(rstar.lt.1.e5) rstr=rstar*6.9598e10
      amdot=amloss*6.3029e25
      RCORE=RSTR
      XMDOT=amdot
      BETAV=beta
      con=amdot/12.566e5
      conr=con/rstr/rstr
      nrext0=ndrad-nd
      zz(nd+nrext0)=0.
      rd(nd+nrext0)=rstr
      rrel(nd+nrext0)=1.
      do iid=1,nd-1
         id=nd-iid
         zz(id+nrext0)=zz(id+1+nrext0)+2.*(dm(id+1)-dm(id))/
     *                 (dens(id+1)+dens(id))
         rd(id+nrext0)=rstr+zz(id+nrext0)
         rrel(id+nrext0)=rd(id+nrext0)/rstr
      end do
C
      do id=1+nrext0,nd+nrext0
         vel0(id)=con/rd(id)**2/dens(id-nrext0)
         vel00(id)=vel0(id)
         if(vel00(id).gt.vinf) vel00(id)=vinf
      end do
      vin=vel0(nrext0+1)
      r1=rrel(nrext0+1)
C
      if(rrel(1+nrext0).lt.rmax.and.nd.lt.ndrad) then
      rl1=1.-1./rrel(1+nrext0)
      rl2=1.-1./rmax
      drl=(rl2-rl1)/nrext0
      do id=1,nrext0
         rlo=rl2-(id-1)*drl
         rrel(id)=1./(1.-rlo)
         rd(id)=rrel(id)*rstr
      end do
      end if
c
      do id=nd+nrext0-1,nrext0+1,-1
         r0=rrel(id)
         numid=0
         do id1=nd+nrext0-1,nrext0+1,-1
            x=un-r0/rrel(id1)
            if(x.lt.1.e-6) x=1.e-6
            v2=vinf*x**beta 
            ind(id1)=0
            if(v2.ge.vel0(id1)) then
              ind(id1)=id1
              numid=numid+1
            end if
         end do
         if(numid.eq.0) go to 10
         rsum=0.
         isum=0
         do id1=nd+nrext0-1,nrext0+1,-1
            if(ind(id1).gt.0) then
               rsum=rsum+rrel(id1)
               isum=isum+id1
            endif
         end do
         rc=rsum/numid
         idc=isum/numid
         numid0=numid
         r00=r0
      end do
   10 continue
      v1=vel0(idc)
      r0=(r0+r00)*0.5
      if(r0.lt.rc) v2=vinf*(un-r0/rc)**beta 
      write(6,602) numid0,idc,rc,r0,v1,v2
  602 format('numid,idc,rc,r0,v1,v2 ',2i4,4f10.5)
c
      do id=nd+nrext0-1,1,-1
         if(rrel(id).gt.rc.and.rrel(id).gt.r0) 
     *      vel0(id)=vinf*(1.-r0/rrel(id))**beta
      end do
c
      t1=temp(1)
      erel=elec(1)/dens(1)
      do id=nd,1,-1
         temp(id+nrext0)=temp(id)
         den0(id+nrext0)=dens(id)
         elec(id+nrext0)=elec(id)
         do i=1,nlevel
            popul(i,id+nrext0)=popul(i,id)
         end do
         WMM(ID+nrext0)=WMM(id)
         WMY(ID+nrext0)=WMY(id)
         YTOT(ID+nrext0)=YTOT(id)
         do i=1,natom
            relab(i,id+nrext0)=relab(i,id)
            abund(i,id+nrext0)=abund(i,id)
         end do
         do i=1,matom
            abndd(i,id+nrext0)=abndd(i,id)
         end do
      end do
C
      do id=1,nrext0
         TEMP(ID)=T1
         WMM(ID)=WMM(NREXT0+1)
         WMY(ID)=WMY(NREXT0+1)
         YTOT(ID)=YTOT(NREXT0+1)
         do i=1,natom
            relab(i,id)=relab(i,nrext0+1)
            abund(i,id)=abund(i,nrext0+1)
         end do
         do i=1,matom
            abndd(i,id)=abndd(i,nrext0+1)
         end do
      end do
      idstd=idstd+nrext0
c
      VINF=vinf*1.e5
      write(6,600)
      do id=1,nd+nrext0
         if(vel0(id).gt.0.) dens(id)=con/rd(id)**2/vel0(id)
         VEL(ID)=vel0(id)*1.e5
c        velc(id)=vel0(id)/2.997925e5
      end do
c
      do id=nd,1,-1
         id1=id+nrext0
         elec(id1)=elec(id1)*dens(id1)/den0(id1)
         do i=1,nlevel
            popul(i,id1)=popul(i,id1)*dens(id1)/den0(id1)
         end do
      end do
c
      do id=1,nrext0
         elec(id)=elec(nrext0+1)*dens(id)/dens(nrext0+1)
         do i=1,nlevel
            popul(i,id)=popul(i,nrext0+1)*dens(id)/dens(nrext0+1)
         end do
      end do
C
      ND=NDRAD
      if(ndf.eq.0) ndf=nd
      do id=1,nd
         write(6,601) id,dm(id),temp(id),elec(id),dens(id),rd(id),
     *                rrel(id),vel0(id)
         write(96,601) id,dm(id),temp(id),elec(id),dens(id),rd(id),
     *                rrel(id),vel0(id),vel00(id)
      end do
  600 format('    ID    M      TEMP       ELEC      DENS             ',
     *       'R       Rrel     VEL'/)
  601 format(1h ,i3,1pe10.3,0pf8.0,1p3e12.3,0pf10.4,0p2f8.2)
C
C
      if(nda.gt.0) then
         XR1=LOG(DENS(1))
         XR2=LOG(DENS(ND))
         DXR=(XR2-XR1)/FLOAT(NDA-1)
         DO ID=1,NDA
            DENSA(ID)=EXP(XR1+FLOAT(ID-1)*DXR)
         END DO
         CALL INTERP(DENS,TEMP,DENSA,TEMPA,ND,NDA,3,1,1)
         CALL INTERP(DENS,ELEC,DENSA,ELECA,ND,NDA,3,1,1)
         CALL INTERP(DENS,RD,DENSA,RDA,ND,NDA,3,1,1)
         CALL INTERP(DENS,RREl,DENSA,RRELA,ND,NDA,3,1,1)
         CALL INTERP(DENS,VEL0,DENSA,VEL0A,ND,NDA,3,1,1)
         do id=1,nda
         write(6,603) id,tempa(id),eleca(id),densa(id),rda(id),
     *                rrela(id),vel0a(id)
         write(96,603) id,tempa(id),eleca(id),densa(id),rda(id),
     *                rrela(id),vel0a(id)
         end do
      end if
  603 format(1h ,i3,0pf8.0,1p3e12.3,0pf10.4,0p2f8.2)
C
  100 continue
      return
      end