c This program computes RT CIA spectra of CO2 at low density limit C at temperatures between 200 and 800K. C C ============================================================================ C Copyright (C) 1997 Marcin Gruszka C ============================================================================ C C Copyright Notice: C You may use this program for your scientific applications, C but please do not distribute it yourself. C Direct all your inquires to the author: e-mail aborysow@stella.nbi.dk C Final request: If you publish your work which benefited from this C program, please acknowledge using this program and quote C the original paper describing the procedure: C M. Gruszka and A. Borysow, C Roto-translational collision-induced absorption of CO2 C for the Atmosphere of Venus at Frequencies from 0 to 250 cm$^{-1}$ C and at temperature from 200K to 800K, C Icarus, vol. 129, pp. 172-177, (1997). C C The paper (Icarus) describing this computer model, is based on original paper by C M. Gruszka and A. Borysow, C "Computer Simulation of the Far Infrared Collision-induced absorption spectra C of gaseous CO2", Mol. Phys., vol. 93, pp. 1007-1016, 1998. C ============================================================================ C C NOTE: The program gives data at spacing corresponding to c Freq(max)-Freq(min)/(npoint-1), c Since reasons are unknown (author's secret), it is recommended to c request one more point that neccesary and get the right freq. step c ================================================================== c c The following input parametres are required (*): c ------------------------------------------------ c (*) some input parameters are Fortran REAL numbers c and should be entered with a decimal point, c not like INTEGER! c c 1. Temperature: (in K) (REAL) c -the model is restricted to the temperature range from c 200 K to 800 K. c c 2. Frequency limits: (in cm^-1) (REAL,REAL) c -the model has been designed to reproduce the MD and the c experimental data up to 250 cm^-1. There is no restriction c on requesting a higher value, but the accuracy may decrease c significantly. c c 3. Number of points in the output file: (less than 1000) (INTEGER) c -the frequency grid is controled by the number of points c in the output file. c (i.e. if the upper limit = 250.0 and the number of points = 50.0, c the spectrum will be computed every 5 cm^-1 ) c c 4. The output in automaticaly directed to output.dat file. c c STRUCTURE OF THE PROGRAM: c -------------------------------------------------------------- c The MAIN program consists of three modules: c c 1. getdat - provides the user interface, c in this subroutine the Temperature, Frequency Limits and c the number of points are read from the keyboard. c c 2. getspc - the main subroutine, where the spectrum is computed. c The three input parameters (temp, range and n; all REAL numbers) c are transmited by value and the absorption coefficient abcoef(..) c is returned. The default lenght of abcoef(..) is 1000 and if c n < 1000, the empty places are set to 0. If a denser grid is c required please change the n appropriaetly in the data declaration c of the MAIN program and inside the getspc. c c 3. putdat - generates the output file. c The frequency (in cm^-1) and c the absorption coefficient (in cm^-1 amagat^-2) are written c to the 'output.dat' file in the followinf format: c format(f10.3,e20.7) c c=========================================================================== c c c BEGINING OF THE MAIN SEGMENT c ------------------------------------------------------- program model implicit none real*8 temp !temperature real*8 range(2) !frequency range real*8 n !number of points in the output real*8 abcoef(1000) !absorption coefficient c ------------------------------------------------------- call getdat(temp,range,n) !reads input parameters call getspc(temp,range,n,abcoef) !computes the spectrum call putdat(temp,range,n,abcoef) !creates the output file end c ------------------------------------------------------- c END OF THE MAIN SEGMENT c c c subroutine getdat(temp,range,n) !reads input parameters implicit none real*8 temp !temperature real*8 range(2) !frequency range real*8 n !number of points in the output c ------------------------------------------------------- real*8 incrmt integer res, no, i 10 continue i=0 write(*,*) 'Input the temperature (200-800 in K),(REAL): ' read(*,*) temp if ((temp .lt. 200.0) .or. (temp .gt. 800.0)) then i=1 write(*,*) 'Temperature outside the 200-800 range!!!' endif if (i .eq. 1) goto 10 c 11 continue i=1 write(*,*) 'Input the frequency range (0-250 cm^-1),' write(*,*) 'the lower and the upper limit (in cm^-1),(REAL,REAL): ' read(*,*) range(1),range(2) if ( range(2) .gt. 250.0 ) then write(*,*) 'The accuracy of the absorption coefficient' write(*,*) 'above 250 cm^-1 may be unpredictable.' write(*,*) 'Do you want to continue? (1-Yes, 0-No) (INTEGER)' read(*,*) i endif if (i .eq. 0) goto 11 c 12 continue i=0 write(*,*) 'Input the number of points in the spectrum. (INTEGER):' read(*,*) no if ((no .lt. 1) .or. (no .gt. 1000)) then i=1 write(*,*) 'Input number between 1 and 1000' write(*,*) 'If you need a denser grid, see point 2' write(*,*) 'in the STRUCTURE OF THE PROGRAM section of the program' endif if (i .eq. 1) goto 12 c n=dfloat(no) incrmt=(range(2)-range(1))/(n-1.0d0) write(*,*) 'OK. Thanks.' write(*,*) 'You entered:' write(*,900) 'Temperature: ', temp 900 format(a13,f12.4) write(*,901) 'Frequency range: ',range(1),' - ',range(2),' cm^-1' 901 format(a17,f12.4,a3,f12.4,a6) write(*,902) 'Frequency increment: ',incrmt 902 format(a21,f12.6) write(*,*) 'Do you want to: proceed (1), or edit the input (0): ? ' read(*,*) res if (res .eq. 0) goto 10 write(*,*) 'The result will be written to "output.dat" file.' write(*,*) 'OK. Computation begins......' return end c c c subroutine getspc(temp,range,n,abcoef) !computes the spectrum implicit none integer nop parameter (nop=1000) real*8 temp !temperature real*8 range(2) !frequency range real*8 n !number of points in the output real*8 abcoef(nop) !absorption coefficient c ------------------------------------------------------- c The A, B and C parameters as given in the paper: integer ndeg parameter (ndeg=3) real*8 ah(ndeg),bh(ndeg),at(ndeg),bt(ndeg),gam(ndeg) c tau^{L}_{1}: data ah /0.1586914382D+13,-0.9344296879D+01,0.6943966881D+00/ c tau^{L}_{2}: data bh /0.1285676961D-12,0.9420973263D+01,-0.7855988401D+00/ c tau^{H}_{1}: data at /0.3312598766D-09,0.7285659464D+01,-0.6732642658D+00/ c tau^{H}_{2}: data bt /0.1960966173D+09,-0.6834613750D+01,0.5516825232D+00/ c gamma_{1}: data gam /0.1059151675D+17,-0.1048630307D+02, 0.7321430968D+00/ c ************************************************************** c The S - shape function parameters: real*8 w1,w2 data w1,w2 /50.0d0, 100.0d0/ c -------------------------------------------------------------- c local variables: c ---------------------------------- real*8 incrmt integer p1,p2,i real*8 a1,b1,a2,b2,gamma real*8 bc1(nop),bc2(nop),bcbc(nop) real*8 scon,mtot,spunit,gm0con real*8 icm, frq, lntemp real*8 xk1,x1,x2 external xk1 incrmt=(range(2)-range(1))/(n-1.0d0) p1=idint(dnint(w1/incrmt)) p2=idint(dnint(w2/incrmt)) lntemp=dlog(temp) a1=1.0d0/(ah(1)*dexp(ah(2)*lntemp+ah(3)*lntemp*lntemp)) b1=bh(1)*dexp(bh(2)*lntemp+bh(3)*lntemp*lntemp) a2=1.0d0/(at(1)*dexp(at(2)*lntemp+at(3)*lntemp*lntemp)) b2=bt(1)*dexp(bt(2)*lntemp+bt(3)*lntemp*lntemp) gamma=gam(1)*dexp(gam(2)*lntemp+gam(3)*lntemp*lntemp) icm=0.1885d0 do 100 i=1,idint(n) frq=dble(i-1)*incrmt+range(1) x1=b1*dsqrt(a1*a1+icm*icm*frq*frq) bc1(i)=dexp(a1*b1)*a1*XK1(x1)/(a1*a1+icm*icm*frq*frq) x2=b2*dsqrt(a2*a2+icm*icm*frq*frq) bc2(i)=dexp(a2*b2)*a2*XK1(x2)/(a2*a2+icm*icm*frq*frq) 100 continue do 200 i=1,idint(n) if (i .le. p1) then bcbc(i)=bc1(i) else if (i .gt. p2) then bcbc(i)=bc2(i) else bcbc(i)=dexp( (1.0d0-dble(i-p1)/dble(p2-p1))*dlog(bc1(i))+ & dble(i-p1)/dble(p2-p1)*dlog(bc2(i)) ) endif 200 continue spunit=1.296917d55 gm0con=1.259009d-6 scon=spunit/temp mtot=gm0con*temp*gamma*1.0d-56 do 300 i=1,idint(n) frq=dble(i-1)*incrmt+range(1) abcoef(i)=scon*mtot*bcbc(i)*frq*frq 300 continue return end c c c subroutine putdat(temp,range,n,abcoef) !creates the output file implicit none real*8 temp !temperature real*8 range(2) !frequency range real*8 n !number of points in the output real*8 abcoef(1000) !absorption coefficient c ----------------------------------------- integer i real*8 incrmt, frq character*10 outfil outfil='output.dat' incrmt=(range(2)-range(1))/(n-1) open(unit=1, file=outfil, form='formatted', status='unknown') write(1,*) 'CIA CO_2-CO_2. Temp=',temp write(1,*) 'Frequency Absorption Coefficient' write(1,*) ' cm^-1 cm^-1 amagat^-2' do 100 i=1,idint(n) frq=dble(i-1)*incrmt+range(1) write(1,900) frq, abcoef(i) 100 continue close(1) 900 format(f12.5,e20.7) write(*,*) 'Done.' return end c c c FUNCTION XK1(X) C MODIFIED BESSEL FUNCTION K1(X) TIMES X C PRECISION IS BETTER THAN 2.2e-7 EVERYWHERE. C ABRAMOWITZ AND S,TEGUN, P.379; TABLES P.417. implicit double precision (a-h,o-z) IF(X-2.) 10,10,20 10 T=(X/3.75)**2 FI1=X*((((((.00032411*T+.00301532)*T+.02658733)*T+.15084934) 1 *T+.51498869)*T+.87890594)*T+.5) T=(X/2.)**2 P=(((((-.00004686*T-.00110404)*T-.01919402)*T-.18156897)*T- 1 .67278579)*T+.15443144)*T+1. XK1=X*dLOG(X/2)*FI1+P RETURN 20 T=2./X P=(((((-.00068245*T+.00325614)*T-.00780353)*T+.01504268)*T- 1 .03655620)*T+.23498619)*T+1.25331414 X=dMIN1(X,330.d0) XK1=dSQRT(X)*dEXP(-X)*P RETURN END