c**********************************************************************
c Creation de profils initiaux pour CECILE 1D
c 
c Author: Amaury juillet 2002 (CERFACS)
c modifs: Amaury sept.   2002
c**********************************************************************
      PROGRAM mprofils
      IMPLICIT NONE

      REAL*8 tmax,tmin,tiu,tfu,demix
      real*8 n
      INTEGER i, choix
      REAL tl,dtl,tu,dtu
      real*8 r, dr, x,xcentre,x2,rmin,rmax
      real*8 temp,cte,tbord
      real*8 a, b,c, perturb
      real*8 titi, toto
      dimension r(0:200),x(0:200),temp(0:200)
      PRINT *,'Rayon min(paroi)'
      READ *,rmin
      PRINT *,'Rayon max(paroi)? '
      READ *,rmax
       PRINT *,'nombre de volume (PAIR ou IMPAIR)? '
      READ *,n
      r(0) = rmin
      r(n+1) = rmax
      r(n)   = rmax
      demix=(rmax-rmin)/2.d+0
      demix = demix + rmin
      dr = (rmax-rmin)/n
      do i=1,n
         r(i)=r(0)+i*dr
      enddo
      open(50,file='rayons.in', status='unknown')
      do i=0,n+1
         write(50,*) r(i)
      enddo
      close(50)
      print *, 'FIN RAYON'
!_________________________________________________________
      !recherche des n+1 valeurs ou j'ai les donnees
      x(0)=r(0)
      x(n+1)=r(n+1)
      ! pour les rayons: bord ext des mailles prendre centre de maille
      do i=1,n
         x(i)= (r(i)+r(i-1))/2.d+0
      enddo
!__________________________________________________________      
      print *,'profil para temp N+2'
      PRINT *,'T milieu? '
      READ *,tmin
      print *,'ou xmilieu: ',(rmax-rmin)/2.d+0 
      xcentre = (rmax-rmin)/2.d+0
      PRINT *,'T bord? '
      READ *,tmax
      print *,'ou xmax:', rmax
      x2 = rmax
      a=(tmin-tmax)/(xcentre*(xcentre-x2))
      c=tmax
      temp(0)=tmax
      temp(n+1)=tmax
      print *, 'variations sur profil de temp ?'
      print *, '0: non'
      print *, '1: oui'
      read (*,*) choix
      do i=1,n
         perturb = 0.d+0
         if (choix.eq.1) perturb = sin(10.*n * x(i)) * 20.d+0
         temp(i) = a*x(i)*(x(i)-x2) + c + perturb
         temp(i) = dabs(temp(i))
      enddo

      OPEN(10,file='temperat.in',status='unknown')
      DO i=0,n+1
         write(10,*) temp(i)
      ENDDO
      CLOSE(10)
!_____________________________________________________________________
      !profil parabol
      !! y=Ax(x-1)+B
      titi = 0.d+0                    !0.155d+0
      PRINT *,'cte fmh2o? '
      read (*,*) titi
      print *,'fm_h2o = ',titi 
      OPEN(10,file='fmh2o.in',status='unknown')
      DO i=1,n
         write(10,*) titi
      ENDDO
      CLOSE(10)
      toto = 0.116d+0
      PRINT *,'cte fmco2? '
      toto = 0.d+0
      print *,'fm_co2 = ',toto 
      OPEN(10,file='fmco2.in',status='unknown')
      DO i=1,n
         write(10,*) toto
      ENDDO
      CLOSE(10)
!.............
      cte = 0.d+0
      OPEN(10,file='fmco.in',status='unknown')
      DO i=1,n
         write(10,*) cte
      ENDDO
      CLOSE(10)
!________________________________________________________________________
      OPEN(10,file='n_tirage.in',status='unknown')
      DO i=0,n+1
         write(10,*) 10000  !!!0
      ENDDO
      CLOSE(10)
      print *, '//!\\ METTRE ', n, ' DANS CECILE.IN'
      END
