      program calc_lines



      implicit none

!     Purpose: to compute the CO2-H2O continuum

      integer nS,nX1,nX2,N
      parameter(nS=224515,nX1=49,nX2=52,N=1600000)

      double precision wn_arr(nS)
      double precision abs_arr(nS)
      double precision EH2O(nS)
      double precision GAMCO2(nS)
      double precision GAMs(nS)
      double precision Beta(nS)
      double precision betas(nS)
      double precision gH2OR(nS)

      double precision wn_X1(nX1)
      double precision abs_X1(nX1)
      double precision wn_X2(nX2)
      double precision abs_X2(nX2)

      double precision wn_f(N)
      double precision abs_f(N)

      double precision temp1, temp2, temp3

      double precision x1_min, x1_max, x2_min, x2_max
      double precision y1_min, y1_max, y2_min, y2_max

      integer k, i, i1, i2

      character*20 file_in
      character*28 file_X1
      character*28 file_X2

      character(LEN=*), parameter :: fmat1 = "(A20, F10.3, F10.3, I7, F7.1)"
      character(LEN=*), parameter :: fmat2 = "(F10.4, E11.3)"

      double precision Na, h, hb, kb, pi_n

      double precision cutoff

      character(LEN=*), parameter :: fmat3 = "(A3,F18.4,A6,E11.3)"

      integer n_K
      parameter(n_K=1000)

      double precision FopartH2O,Fpart,QH2O

      double precision T,T0

      double precision Qcoef1,Qcoef2,Qcoef3,Qcoef4

      double precision Ptot, frH2O, frCO2, dH2O, dCO2

      Qcoef1=-0.37688d+01 ! coefficients for the partition function
      Qcoef2=0.26168d+00
      Qcoef3=0.13497d-02
      Qcoef4=0.66013d-06

      pi_n=3.14159265359 ! constants
      h=6.62e-34
      hb=h/(2.*pi_n)
      kb=1.38e-23
      Na=6.02e23

      T0=296.d0             ! reference temperature 
      FopartH2O=Qcoef1+Qcoef2*T0+Qcoef3*T0*T0+Qcoef4*T0*T0*T0 ! partition function at T0 (296K)

      T=296.d0              ! temperature of the calculation
      Fpart=Qcoef1+Qcoef2*T+Qcoef3*T*T+Qcoef4*T*T*T ! partition function at T

      file_in='01_hit12_Gamache.par'
      file_X1='chi_positif_H2O-CO2_296K.dat'
      file_X2='chi_negatif_H2O-CO2_296K.dat'

      cutoff=1d-6 ! non-zero cutoff
      Ptot=1000./1013.25 ! total pressure in atm
      frH2O=0.01                    ! fraction of H2O
      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)

1001  FORMAT(3X,F12.6,E10.3,10X,F5.4,F5.3,F10.4,F4.2,115X,F4.2)       

      open(1,file=file_in,form='formatted',status='old')
         do k=1,nS
            read(1,1001) wn_arr(k),abs_arr(k),GAMCO2(k),GAMs(k),EH2O(k),Beta(k),betas(k)
            gH2OR(k)=Ptot*(1.d0-frH2O)*GAMCO2(k)*(T0/T)**Beta(k)+Ptot*frH2O*GAMs(k)*(T0/T)**betas(k)
         end do         
      close(1)

      do k=1,nS
         abs_arr(k)=abs_arr(k)*(FopartH2O/Fpart)*  &
         exp(-1.4388*EH2O(k)*(1/T-1/T0))/  &
         (1-exp(-1.4388*wn_arr(k)/T0))*  &
         (1-exp(-1.4388*wn_arr(k)/T))
      end do    

      open(2,file=file_X1,form='formatted',status='old') ! Read Khi+ factors
         do k=1,nX1
            read(2,*) wn_X1(k), abs_X1(k)
         end do         
      close(2)  

      open(3,file=file_X2,form='formatted',status='old') ! Read Khi- factors
         do k=1,nX2
            read(3,*) wn_X2(k), abs_X2(k)
         end do         
      close(3)

      ! Normalize Khi-factors at 1 at sigma-sigma0 = 0.
      do k=2,nX1
         abs_X1(k)=abs_X1(k)/abs_X1(1)
      end do
      abs_X1(1)=1.d0
      do k=2,nX2
         abs_X2(k)=abs_X2(k)/abs_X2(1)
      end do
      abs_X2(1)=1.d0

      ! generate the wavenumber grid
      do i=1,N 
         wn_f(i)=1.d-3+(i-1)*1600./N
         abs_f(i)=0.
         !write(*,*) wn_f(i)
      enddo

      ! loop to calculate the Ma & Tipping 92 terms of equation 40
      do i=1,N 
         do k=1, nS

           temp1=0.
           temp2=0.
           temp3=0.

           ! STEP 1 - INTERPOLATION OF THE KHI+ AND KHI- FACTORS
           if( ((wn_f(i)+wn_arr(k)) .lt. 1274.7) .and. (abs(wn_f(i)+wn_arr(k)) .gt. cutoff)) then
              i2=1
              do while(abs(wn_f(i)+wn_arr(k)) .gt. wn_X2(i2))
                 i2=i2+1
              enddo
              x2_min=wn_X2(i2-1)
              x2_max=wn_X2(i2)
              y2_min=abs_X2(i2-1)
              y2_max=abs_X2(i2)
              temp2=y2_min*exp((log(abs(wn_f(i)+wn_arr(k)))-log(x2_min))*((log(y2_max)-log(y2_min))/  &  
              (log(x2_max)-log(x2_min))))
              !write(*,*), 'case 2 ', x2_min, x2_max, y2_min, y2_max, wn_f(i)+wn_arr(k), temp2
           endif

           if( ((wn_f(i)-wn_arr(k)) .gt. cutoff) .and. ((wn_f(i)-wn_arr(k)) .lt. 1274.7) ) then
               i1=1
               do while((wn_f(i)-wn_arr(k)) .gt. wn_X2(i1))
                  i1=i1+1
               enddo
               x1_min=wn_X2(i1-1)
               x1_max=wn_X2(i1)
               y1_min=abs_X2(i1-1)
               y1_max=abs_X2(i1)
               temp1=y1_min*exp((log((wn_f(i)-wn_arr(k)))-log(x1_min))*((log(y1_max)-log(y1_min))/  &  
               (log(x1_max)-log(x1_min))))
               !write(*,*), 'case 1 ', x1_min, x1_max, y1_min, y1_max, wn_f(i)-wn_arr(k), temp1
           endif

           if( ((-wn_f(i)+wn_arr(k)) .gt. cutoff) .and. ((-wn_f(i)+wn_arr(k)) .lt. 1986.9) ) then
               i1=1
               do while((-wn_f(i)+wn_arr(k)) .gt. wn_X1(i1))
                  i1=i1+1
               enddo
               x1_min=wn_X1(i1-1)
               x1_max=wn_X1(i1)
               y1_min=abs_X1(i1-1)
               y1_max=abs_X1(i1)
               temp1=y1_min*exp((log((-wn_f(i)+wn_arr(k)))-log(x1_min))*((log(y1_max)-log(y1_min))/  &  
               (log(x1_max)-log(x1_min))))
               !write(*,*), 'case 1 ', x1_min, x1_max, y1_min, y1_max, wn_f(i)-wn_arr(k), temp1
           endif

           ! CALCULATE THE EXPRESSION (equation 40 of Ma & Tipping) OF THE ABSORPTION BY H2O LINES
           temp3=abs_arr(k)*wn_f(i)*sinh(hb*wn_f(i)/(2.*kb*T))/(wn_arr(k)*pi_n*sinh(hb*wn_arr(k)/(2.*kb*T)))
           abs_f(i)=abs_f(i)+temp3*( (temp1*gH2OR(k)/((wn_f(i)-wn_arr(k))*(wn_f(i)-wn_arr(k))+gH2OR(k)**2.)) & 
                                    +(temp2*gH2OR(k)/((wn_f(i)+wn_arr(k))*(wn_f(i)+wn_arr(k))+gH2OR(k)**2.)) ) ! 0.1 comes from Ma & Tipping 92 fig 4 at 1cm-1
         enddo
         write(*,*) 'We are at i = ', i, ' for sigma = ', wn_f(i), ' cm-1', ' and abs_f(i) = ', abs_f(i)
      enddo

      open(10,file='h2o-co2_296K_recalc_nov18_0cm_cutoff.dat',form='formatted',status='new')
         do i=1,N
            write(10,*) wn_f(i), abs_f(i)
         end do         
      close(10)


      end program




