C Programme qui calcule la transmission d'un mélange de CO2+H2O
      implicit none
      integer nsig,nsigmx,i,nraimx,nrCO2,nrH2O,isig,ir,irai
      Parameter(nsigmx=5000000,nraimx=5000000)
      double precision T0,A,pi,c,cte,cte1
      double precision T,Ptot,frH2O,frCO2
      double precision FopartH2O,FopartCO2,QH2O,QCO2,amasseH2O,amasseCO2
     1,Fpart
      double precision sgmin,sgmax,sigmoy,pas
      double precision ai,bi,ci,qi,af,bf,df,ef,ff,cf
C Donnees exp
      integer iexp,Nexp,Ntot
      character fichier*21
      double precision Temp,PH2O,aL,dCO2,dH2O
      double precision sg,TransI,TransI0
      double precision sexp(nsigmx),Trans(Nsigmx),Trans0(Nsigmx)
C Fonction d'appareil
      Integer nconv,iconv
      double precision sum,sum2,OPD
      double precision SigC(-20000:20000),Fapp(-20000:20000)
C Donnees HITRAN
      integer ISO0,J,m,mabs
      double precision sigma,ST0,GAMair,E,beta,del
      character band*23,branch*1
      double precision anw(nraimx),gCO2w(nraimx)
C Données pour CA de H2O pur
      double precision SS,AA,SS2,AA2,SS3,AA3
      double precision SCA(Nraimx),ACA(Nraimx),CAH2O(Nraimx),
     1TCAH2O(Nraimx),CAH2OCO2(Nraimx),TCAH2OCO2(Nraimx),
     1ACA2(Nraimx),SCA2(Nraimx),aCA3(Nraimx),aCA4(Nraimx),
     1sCA3(Nraimx) 
C Données retenues pour CO2
      double precision gDCO2
      double precision SgCO2R(nraimx),SCO2R(nraimx),gCO2R(nraimx)
     1,ECO2(nraimx),bCO2R(nraimx),dCO2R(nraimx)
      double precision SgCO2L(nraimx),SCO2L(nraimx),gCO2wL(nraimx)
     1,gCO2sL(nraimx),gCO2L(Nraimx)
C Données retenues pour H2O
      double precision gDH2O,betas,gamCO2,gamS
      double precision SgH2OR(Nraimx),SH2OR(nraimx),gH2OR(nraimx),
     1EH2O(nraimx),bH2OR(nraimx),dH2OR(nraimx)
      double precision SgH2OL(nraimx),gH2OL(nraimx),SH2OL(nraimx)
C Données pour facteur khi
      double precision sig1,sig2,sig3,BB1,BB2,BB3
C
      integer indc,indloc
      double precision AAkhi,akhi,dsig,x,y,wr,wi,loschmidt
	double precision aa0(Nsigmx)
      double precision sig(Nsigmx),AbsCO2H2O(Nsigmx),absCO2CO2(Nsigmx)
     1,AbsH2O(nsigmx),AbsV(nsigmx),AbsVl(nsigmx),AbsCO2H2Olor(Nsigmx)
     1,absCO2CO2lor(Nsigmx),AV(Nsigmx)
      double precision sumH2O,sumCO2,sumCACO2
      
C
      T0=296.d0             ! reference temperature 
      A=1.4388D0				
      pi=4.d0*datan(1.d0)
      FopartH2O=QH2O(T0)    ! partition function of H2O at T0
      FopartCO2=QCO2(T0)
      amasseH2O=18.d-03     ! mass of H2O
      amasseCO2=44.d-03
      loschmidt = 2.686777d19
      
	aL=1.5175d4        ! longueur de la cellule
      sgmin=10.d0         ! first wavenumber 
      sgmax=1600.d0	 ! last wavenumber		
      sigmoy=(sgmin+sgmax)/2		  
      pas=10.d+0                   ! pas de calcul
      nsig=(sgmax-sgmin)/pas+1
      if(nsig.gt.nsigmx)then
          write(6,*)'pb nsigmx trop petit'
          stop
      end if            
       
C      
C Lecture du fichier des données P,T      
C      
      open(5,file='Liste_donnees_tst.txt')
      read(5,*)
C501   format(a22,6x,f5.2,12x,f5.1,4x,f5.1,6x,f6.2,6x,f6.4,3x,f6.4) 
501   format(a21,3x,f5.1,3x,f6.1,3x,f7.5)
C150   read(5,501,end=15)fichier,Temp,Ptoti,Ptotf,PH2O,dH2O,dCO2   ! Ptoti, Ptotf: pressions totales au début et à la fin de la manip, dH2O et dCO2 sont en amagat, les pressions sont en Torr
150   read(5,501,end=15)fichier,Temp,Ptot,frH2O ! frH2O = fraction of H2O
      
      Ptot=Ptot/1013.25
      T=Temp
      frCO2=1.d0-frH2O              ! fraction of CO2
      dH2O=Ptot*frH2O*(273.15/T)    ! density = (p/p0)*(T0/T) in amagat
      dCO2=Ptot*frCO2*(273.15/T)
      write(6,*)fichier,T,Ptot,frCO2



C*************************************************************
C Lecture de données de H2O
C      
        
        write(6,*)'lecture de la base de donnes de raies pour H2O'
        i=0
1001  FORMAT(2X,I1,F12.6,E10.3,10X,F5.4,F5.3,F10.4,F4.2,F8.6,50X,57X,
     1  f4.2)               
        OPEN(10,file='01_hit12_Gamache.par',status='old')
1013  READ(10,1001,END=2013)ISO0,Sigma,ST0,GAMCO2,GAMs,E,Beta,del,
     1 betas
        IF(Sigma.LT.(sgmin-25.d0)) GOTO 1013
        IF(Sigma.GT.(sgmax+25.d0)) GOTO 1013
        i=i+1
          if(i.gt.nraimx)then
          write(6,*)'pb nraimx trop petit'
          stop
          end if
        SgH2OR(i)=Sigma
        SH2OR(i)=ST0
        gH2OR(i)=Ptot*(1.d0-frH2O)*GAMCO2*(T0/T)**Beta+
     1       Ptot*frH2O*GAMs*(T0/T)**Betas
        EH2O(i)=E
        bH2OR(i)=Beta
        dH2OR(i)=del
        goto 1013
2013    close(10)
        nrH2O=i
        write(6,*)'nbre de raies de H2O',nrH2O
C      
      Fpart=QH2O(T)
      c=Ptot*frH2O*0.1013d0/(1.38d-23*T)	 
      gDH2O=dsqrt((2.*dlog(2.d0)*8.314*T)/
     +       (AmasseH2O*(3.d+08**2)))*sigmoy
      cte=dsqrt(dlog(2.d0))/gDH2O
      
      do 1538 irai=1,nrH2O
          SH2OL(irai)=c*
     /       SH2OR(irai)*
     /       (FopartH2O/Fpart)*
     /       exp(-1.4388*EH2O(irai)*(1/T-1/T0))/
     /       (1-exp(-1.4388*SgH2OR(irai)/T0))*
     /       (1-exp(-1.4388*SgH2OR(irai)/T))
          SgH2OL(irai)=sgH2OR(irai)+dH2OR(irai)*Ptot
          gH2OL(irai)= gH2OR(irai)
1538  continue  
      write(6,*)'H2O'

C      
C Calcul de l'absorption ******************************************************************************
C
C
C due aux raies permises de H2O
      cte=dsqrt(dlog(2.d0))/gDH2O
      cte1=cte/dsqrt(pi)
C Determination de l'absorption à 25 cm-1 pour chaque raie
	do irai=1,nrH2O
	X=25.d0*cte
	Y=gH2OL(irai)*cte
	call CPF0(x,y,wr,wi)
	AA0(irai)=cte1*(WR*SH2OL(irai))
	enddo
C Calcul de l'absoprtion
      do 201 isig=1,Nsig
          AbsH2O(isig)=0.d0
	      sig(isig)=sgmin+(isig-1)*pas
            do 101 irai=1,nrH2O
		      dsig=sgH2OR(irai)-sig(isig)
	          
                if (dabs(dsig).le.25.d0) then
                    aa=AA0(irai)
		    AbsH2O(isig)=AbsH2O(isig)+aa
                endif
101        continue
!	write(17,*)sig(isig),AbsH2O(isig)
201   continue
       write(6,*)'fin du calcul de labsprtion de H2O'

106   format(f12.7,3(3x,e14.7))            
2002  format(f16.8,50(2x,e12.7))      
      open(16,file=fichier//'_plinthe_25cm-1.dat')!,type='unknown')       ! Fichier pour écrire le résultat
      do isig=1,Nsig
         write(16,*)sig(isig),AbsH2O(isig)
      enddo
      close(16) 
      !stop
      
      go to 150
15    close(5)         
      stop
      end   


      Subroutine CPF0(X,Y,WR,WI)
C*********************************************************************
C "CPF": Complex Probability Function
C .........................................................
C         .       Subroutine to Compute the Complex       .
C         .        Probability Function W(z=X+iY)         .
C         .     W(z)=exp(-z**2)*Erfc(-i*z) with Y>=0      .
C         .    Which Appears when Convoluting a Complex   .
C         .     Lorentzian Profile by a Gaussian Shape    .
C         .................................................
C
C             WR : Real Part of W(z)
C             WI : Imaginary Part of W(z)
C
C This Routine was Taken from the Paper by J. Humlicek, which 
C is Available in Page 309 of Volume 21 of the 1979 Issue of
C the Journal of Quantitative Spectroscopy and Radiative Transfer
C Please Refer to this Paper for More Information
C*********************************************************************
C      
      Implicit None
      Integer I
      doubleprecision X,Y,WR,WI
      doubleprecision T,U,S,Y1,Y2,Y3,R,R2,D,D1,D2,D3,D4
C      
      Dimension T(6),U(6),S(6)
      Data T/.314240376d0,.947788391d0,1.59768264d0,2.27950708d0
     ,        ,3.02063703d0,3.8897249d0/
      Data U/1.01172805d0,-.75197147d0,1.2557727d-2,1.00220082d-2
     ,        ,-2.42068135d-4,5.00848061d-7/
      Data S/1.393237d0,.231152406d0,-.155351466d0,6.21836624d-3
     ,        ,9.19082986d-5,-6.27525958d-7/
      WR=0.d0
      WI=0.d0
      Y1=Y+1.5d0
      Y2=Y1*Y1
      If( (Y.GT.0.85d0) .OR. (DABS(X).LT.(18.1d0*Y+1.65d0)) )GoTo 2
C
C       Region II
C
      If( DABS(X).LT.12.d0 )WR=DEXP(-X*X)
      Y3=Y+3.d0
      DO 1 I=1,6
      R=X-T(I)
      R2=R*R
      D=1.d0/(R2+Y2)
      D1=Y1*D
      D2=R*D
      WR=WR+Y*(U(I)*(R*D2-1.5d0*D1)+S(I)*Y3*D2)/(R2+2.25d0)
      R=X+T(I)
      R2=R*R
      D=1.d0/(R2+Y2)
      D3=Y1*D
      D4=R*D
      WR=WR+Y*(U(I)*(R*D4-1.5d0*D3)-S(I)*Y3*D4)/(R2+2.25d0)
      WI=WI+U(I)*(D2+D4)+S(I)*(D1-D3)
 1    Continue  
      Return
C
C       Region I
C
 2    Continue
      DO 3 I=1,6
      R=X-T(I)
      D=1.d0/(R*R+Y2)
      D1=Y1*D
      D2=R*D
      R=X+T(I)
      D=1.d0/(R*R+Y2)
      D3=Y1*D
      D4=R*D
      WR=WR+U(I)*(D1+D3)-S(I)*(D2-D4)
      WI=WI+U(I)*(D2+D4)+S(I)*(D1-D3)
 3    Continue  
      Return
      End
 

  
	          doubleprecision Function QH2O(T)
C Compute the H2O partition function       
       doubleprecision T,Qcoef(4)
       Data Qcoef/-0.37688d+01,0.26168d+00,
     ,             0.13497d-02,0.66013d-06/
       QH2O=Qcoef(1)+
     +     Qcoef(2)*T+
     +     Qcoef(3)*T*T+
     +     Qcoef(4)*t*t*t
       Return
       End
				doubleprecision Function QCO2(T)
C Compute the CO2 partition function       
       doubleprecision T,Qcoef(4)
       Data Qcoef/-0.21995d+01,0.96751d+00,
     ,             -0.80827d-03,0.28040d-05/
       QCO2=Qcoef(1)+
     +     Qcoef(2)*T+
     +     Qcoef(3)*T*T+
     +     Qcoef(4)*t*t*t
       Return
       End
