subroutine read_model(m,alt,pres,temp,Tsol,Tesp,
     &     dens_mix,adia,grav,zeta,capa)
      implicit none
c
c     Purpose: to read atmospheric data contained in file "modelPT.VIRA" for Venus' atmosphere
c
c     Output:
c       + m: number of atmospheric levels
c       + alt: altitude at each level (m)
c       + pres: pressure at each level (Pa)
c       + temp: temperature at each level (K)
c       + Tsol: temperature of the ground (K)
c       + Tesp: temperature of space (K)
c       + dens_mix: density of the gas mixture, at each level (kg/m³)
c       + adia: adiabatic lapse rate at each level (K/m)
c       + grav: gravity acceleration at each level (m/s²)
c       + zeta: non-ideality parameter (P.V=zeta.n.R.T) at eahc level 
c       + capa: calorific capacity at each level (J/kg/K)
c
      include 'max.inc'
      include 'formats.inc'
      integer i,m,ios,strlen,dim,n
      double precision tmp1(1:Nmax),tmp2(1:Nmax)
      double precision tmp3(1:Nmax),tmp4(1:Nmax)
      double precision tmp5(1:Nmax),tmp6(1:Nmax)
      double precision tmp7(1:Nmax),tmp8(1:Nmax)
      double precision alt(1:Nmax)
      double precision pres(1:Nmax)
      double precision adia(1:Nmax)
      double precision temp(1:Nmax)
      double precision Tsol
      double precision Tesp
      double precision dens_mix(1:Nmax)
      double precision grav(1:Nmax)
      double precision zeta(1:Nmax)
      double precision capa(1:Nmax)
      double precision raf
      character*(Nchar_mx) file_model,file_Tsmooth
      character*(Nchar_mx) strtmp1,strtmp2

c     Récupération des données d'altitude et de température
      file_model='./modelPT.VIRA'
      open(11,file=file_model(1:strlen(file_model)),
     &        status='old',iostat=ios)
      if (ios.ne.0) then        ! file not found
         write(*,*) 'Error from routine read_model:'
         write(*,*) 'Data file could not be found:'
         write(*,*) file_model(1:strlen(file_model))
         stop
      endif
      write(*,*) 'Reading atmospheric profile data from file:'
      write(*,*) file_model(1:strlen(file_model))
      do i=1,2
         read(11,*)
      enddo
      read(11,33) strtmp1,dim,strtmp2,m
c     Debug
c      write(*,*) 'dim=',dim,' m=',m
c     Debug
      do i=1,4
         read(11,*)
      enddo
c     m= number of atmospheric levels (interfaces)
      do i=1,m
         read(11,*) tmp1(i),tmp2(i),tmp3(i),tmp4(i),tmp5(i),
     &        tmp6(i),tmp7(i),tmp8(i)
      enddo
      close(11)

      do i=1,m
         alt(i)=tmp1(m-i+1)*1.0D+3    ! altitudes (m)
         temp(i)=tmp2(m-i+1)          ! temperatures (K)
         pres(i)=tmp3(m-i+1)*1.0D+5 ! pressures (bar -> Pa)
         dens_mix(i)=tmp4(m-i+1)          ! densities of the mix (kg/m³)
         adia(i)=tmp5(m-i+1)/1.0D+3   ! adiabatic lapse rate (K/m)
         grav(i)=tmp6(m-i+1)          ! gravity accelerations (m/s²)
         zeta(i)=tmp7(m-i+1)          ! non-ideality parameters
         capa(i)=tmp8(m-i+1)          ! calorific capacities (J/kg/K)
      enddo

c     Overwrite temperature profile by a smoothed temperature profile
c     (no discontinuity of the second-order derivative)
      file_Tsmooth='./Tsmooth.txt'
      open(15,file=file_Tsmooth(1:strlen(file_Tsmooth)),
     &        status='old',iostat=ios)
      if (ios.ne.0) then        ! file not found
         write(*,*) 'Error from routine read_model:'
         write(*,*) 'Data file could not be found:'
         write(*,*) file_Tsmooth(1:strlen(file_Tsmooth))
         stop
      endif
      write(*,*) 'Using smooth temperature profile from file:'
      write(*,*) file_Tsmooth(1:strlen(file_Tsmooth))
      do i=1,m
         read(15,*) raf,temp(i)
      enddo
      close(15)

c     Ground temperature is taken in the gas, at ground level (no discontinuity)
c     Space temperature is given at 3K
      Tsol=temp(1)
      Tesp=3.0D+0

      return
      end