| 1 | c This program computes RT CIA spectra of CO2 at low density limit |
|---|
| 2 | C at temperatures between 200 and 800K. |
|---|
| 3 | C |
|---|
| 4 | C ============================================================================ |
|---|
| 5 | C Copyright (C) 1997 Marcin Gruszka |
|---|
| 6 | C ============================================================================ |
|---|
| 7 | C |
|---|
| 8 | C Copyright Notice: |
|---|
| 9 | C You may use this program for your scientific applications, |
|---|
| 10 | C but please do not distribute it yourself. |
|---|
| 11 | C Direct all your inquires to the author: e-mail aborysow@stella.nbi.dk |
|---|
| 12 | C Final request: If you publish your work which benefited from this |
|---|
| 13 | C program, please acknowledge using this program and quote |
|---|
| 14 | C the original paper describing the procedure: |
|---|
| 15 | C M. Gruszka and A. Borysow, |
|---|
| 16 | C Roto-translational collision-induced absorption of CO2 |
|---|
| 17 | C for the Atmosphere of Venus at Frequencies from 0 to 250 cm$^{-1}$ |
|---|
| 18 | C and at temperature from 200K to 800K, |
|---|
| 19 | C Icarus, vol. 129, pp. 172-177, (1997). |
|---|
| 20 | |
|---|
| 21 | C |
|---|
| 22 | C The paper (Icarus) describing this computer model, is based on original paper by |
|---|
| 23 | C M. Gruszka and A. Borysow, |
|---|
| 24 | C "Computer Simulation of the Far Infrared Collision-induced absorption spectra |
|---|
| 25 | C of gaseous CO2", Mol. Phys., vol. 93, pp. 1007-1016, 1998. |
|---|
| 26 | C ============================================================================ |
|---|
| 27 | C |
|---|
| 28 | C NOTE: The program gives data at spacing corresponding to |
|---|
| 29 | c Freq(max)-Freq(min)/(npoint-1), |
|---|
| 30 | c Since reasons are unknown (author's secret), it is recommended to |
|---|
| 31 | c request one more point that neccesary and get the right freq. step |
|---|
| 32 | c ================================================================== |
|---|
| 33 | |
|---|
| 34 | c |
|---|
| 35 | c The following input parametres are required (*): |
|---|
| 36 | c ------------------------------------------------ |
|---|
| 37 | c (*) some input parameters are Fortran REAL numbers |
|---|
| 38 | c and should be entered with a decimal point, |
|---|
| 39 | c not like INTEGER! |
|---|
| 40 | c |
|---|
| 41 | c 1. Temperature: (in K) (REAL) |
|---|
| 42 | c -the model is restricted to the temperature range from |
|---|
| 43 | c 200 K to 800 K. |
|---|
| 44 | c |
|---|
| 45 | c 2. Frequency limits: (in cm^-1) (REAL,REAL) |
|---|
| 46 | c -the model has been designed to reproduce the MD and the |
|---|
| 47 | c experimental data up to 250 cm^-1. There is no restriction |
|---|
| 48 | c on requesting a higher value, but the accuracy may decrease |
|---|
| 49 | c significantly. |
|---|
| 50 | c |
|---|
| 51 | c 3. Number of points in the output file: (less than Npmx) (INTEGER) |
|---|
| 52 | c -the frequency grid is controled by the number of points |
|---|
| 53 | c in the output file. |
|---|
| 54 | c (i.e. if the upper limit = 250.0 and the number of points = 50.0, |
|---|
| 55 | c the spectrum will be computed every 5 cm^-1 ) |
|---|
| 56 | c |
|---|
| 57 | c 4. The output in automaticaly directed to output.dat file. |
|---|
| 58 | c |
|---|
| 59 | c STRUCTURE OF THE PROGRAM: |
|---|
| 60 | c -------------------------------------------------------------- |
|---|
| 61 | c The MAIN program consists of three modules: |
|---|
| 62 | c |
|---|
| 63 | c 1. getdat - provides the user interface, |
|---|
| 64 | c in this subroutine the Temperature, Frequency Limits and |
|---|
| 65 | c the number of points are read from the keyboard. |
|---|
| 66 | c |
|---|
| 67 | c 2. getspc - the main subroutine, where the spectrum is computed. |
|---|
| 68 | c The three input parameters (temp, range and n; all REAL numbers) |
|---|
| 69 | c are transmited by value and the absorption coefficient abcoef(..) |
|---|
| 70 | c is returned. The default lenght of abcoef(..) is Npmx and if |
|---|
| 71 | c n < Npmx, the empty places are set to 0. If a denser grid is |
|---|
| 72 | c required please change the n appropriaetly in the data declaration |
|---|
| 73 | c of the MAIN program and inside the getspc. |
|---|
| 74 | c |
|---|
| 75 | c 3. putdat - generates the output file. |
|---|
| 76 | c The frequency (in cm^-1) and |
|---|
| 77 | c the absorption coefficient (in cm^-1 amagat^-2) are written |
|---|
| 78 | c to the 'output.dat' file in the followinf format: |
|---|
| 79 | c format(f10.3,e20.7) |
|---|
| 80 | c |
|---|
| 81 | c=========================================================================== |
|---|
| 82 | c |
|---|
| 83 | c |
|---|
| 84 | c --------------------------------------------------------------------------------------------------------------- |
|---|
| 85 | c 2008/10/15 |
|---|
| 86 | c |
|---|
| 87 | c Modifications of the code by V. Eymet as follows: |
|---|
| 88 | c + only the "getspc" routine remains |
|---|
| 89 | c + inputs: temp (temperature in K), wn (wavenumber in cm^-1); output: abcoef (absorption coefficient, in cm^-1 amagat^-2) |
|---|
| 90 | c --------------------------------------------------------------------------------------------------------------- |
|---|
| 91 | c |
|---|
| 92 | c |
|---|
| 93 | c |
|---|
| 94 | c --------------------------------------------------------------------------------------------------------------- |
|---|
| 95 | c 2009/2/28 |
|---|
| 96 | c |
|---|
| 97 | c Modifications of the code by R. Wordsworth as follows: |
|---|
| 98 | c + "baranov" routine added to calculate CIA dimer spectrum between 1100 and 1600 cm^-1 |
|---|
| 99 | c + this uses recorded data rather than an analytical function, so bilinear interpolation is needed |
|---|
| 100 | c + inputs and outputs are identical to those of getspc |
|---|
| 101 | c --------------------------------------------------------------------------------------------------------------- |
|---|
| 102 | |
|---|
| 103 | |
|---|
| 104 | |
|---|
| 105 | subroutine getspc(temp,wn,abcoef) ! computes the spectrum |
|---|
| 106 | implicit none |
|---|
| 107 | real*8 temp !temperature |
|---|
| 108 | real*8 wn !wavenumber |
|---|
| 109 | real*8 abcoef !absorption coefficient |
|---|
| 110 | c ------------------------------------------------------- |
|---|
| 111 | c The A, B and C parameters as given in the paper: |
|---|
| 112 | integer ndeg |
|---|
| 113 | parameter (ndeg=3) |
|---|
| 114 | real*8 ah(ndeg),bh(ndeg),at(ndeg),bt(ndeg),gam(ndeg) |
|---|
| 115 | c tau^{L}_{1}: |
|---|
| 116 | data ah /0.1586914382D+13,-0.9344296879D+01,0.6943966881D+00/ |
|---|
| 117 | c tau^{L}_{2}: |
|---|
| 118 | data bh /0.1285676961D-12,0.9420973263D+01,-0.7855988401D+00/ |
|---|
| 119 | c tau^{H}_{1}: |
|---|
| 120 | data at /0.3312598766D-09,0.7285659464D+01,-0.6732642658D+00/ |
|---|
| 121 | c tau^{H}_{2}: |
|---|
| 122 | data bt /0.1960966173D+09,-0.6834613750D+01,0.5516825232D+00/ |
|---|
| 123 | c gamma_{1}: |
|---|
| 124 | data gam /0.1059151675D+17,-0.1048630307D+02, 0.7321430968D+00/ |
|---|
| 125 | c ************************************************************** |
|---|
| 126 | c The S - shape function parameters: |
|---|
| 127 | real*8 w1,w2 |
|---|
| 128 | data w1,w2 /50.0d0, 100.0d0/ |
|---|
| 129 | c -------------------------------------------------------------- |
|---|
| 130 | c local variables: |
|---|
| 131 | c ---------------------------------- |
|---|
| 132 | real*8 incrmt |
|---|
| 133 | integer p1,p2,n |
|---|
| 134 | real*8 a1,b1,a2,b2,gamma |
|---|
| 135 | real*8 bc1,bc2,bcbc |
|---|
| 136 | real*8 scon,mtot,spunit,gm0con |
|---|
| 137 | real*8 icm, frq, lntemp |
|---|
| 138 | real*8 xk1,x1,x2 |
|---|
| 139 | external xk1 |
|---|
| 140 | incrmt=250.0D+0 |
|---|
| 141 | n=1 |
|---|
| 142 | p1=idint(dnint(w1/incrmt)) |
|---|
| 143 | p2=idint(dnint(w2/incrmt)) |
|---|
| 144 | lntemp=dlog(temp) |
|---|
| 145 | a1=1.0d0/(ah(1)*dexp(ah(2)*lntemp+ah(3)*lntemp*lntemp)) |
|---|
| 146 | b1=bh(1)*dexp(bh(2)*lntemp+bh(3)*lntemp*lntemp) |
|---|
| 147 | a2=1.0d0/(at(1)*dexp(at(2)*lntemp+at(3)*lntemp*lntemp)) |
|---|
| 148 | b2=bt(1)*dexp(bt(2)*lntemp+bt(3)*lntemp*lntemp) |
|---|
| 149 | gamma=gam(1)*dexp(gam(2)*lntemp+gam(3)*lntemp*lntemp) |
|---|
| 150 | icm=0.1885d0 |
|---|
| 151 | frq=wn |
|---|
| 152 | x1=b1*dsqrt(a1*a1+icm*icm*frq*frq) |
|---|
| 153 | bc1=dexp(a1*b1)*a1*XK1(x1)/(a1*a1+icm*icm*frq*frq) |
|---|
| 154 | x2=b2*dsqrt(a2*a2+icm*icm*frq*frq) |
|---|
| 155 | bc2=dexp(a2*b2)*a2*XK1(x2)/(a2*a2+icm*icm*frq*frq) |
|---|
| 156 | if (p1.ge.1) then |
|---|
| 157 | bcbc=bc1 |
|---|
| 158 | else if (p2.lt.1) then |
|---|
| 159 | bcbc=bc2 |
|---|
| 160 | else |
|---|
| 161 | bcbc=dexp( (1.0d0-dble(1-p1)/dble(p2-p1))*dlog(bc1)+ |
|---|
| 162 | & dble(1-p1)/dble(p2-p1)*dlog(bc2) ) |
|---|
| 163 | endif |
|---|
| 164 | spunit=1.296917d55 |
|---|
| 165 | gm0con=1.259009d-6 |
|---|
| 166 | scon=spunit/temp |
|---|
| 167 | mtot=gm0con*temp*gamma*1.0d-56 |
|---|
| 168 | frq=wn |
|---|
| 169 | abcoef=scon*mtot*bcbc*frq*frq |
|---|
| 170 | return |
|---|
| 171 | end |
|---|
| 172 | c |
|---|
| 173 | c |
|---|
| 174 | FUNCTION XK1(X) |
|---|
| 175 | C MODIFIED BESSEL FUNCTION K1(X) TIMES X |
|---|
| 176 | C PRECISION IS BETTER THAN 2.2e-7 EVERYWHERE. |
|---|
| 177 | C ABRAMOWITZ AND S,TEGUN, P.379; TABLES P.417. |
|---|
| 178 | implicit double precision (a-h,o-z) |
|---|
| 179 | ! IF(X-2.) 10,10,20 |
|---|
| 180 | IF(X-2..le.0.0)then |
|---|
| 181 | ! 10 T=(X/3.75)**2 |
|---|
| 182 | T=(X/3.75)**2 |
|---|
| 183 | FI1=X*((((((.00032411*T+.00301532)*T+.02658733)*T+.15084934) |
|---|
| 184 | 1 *T+.51498869)*T+.87890594)*T+.5) |
|---|
| 185 | T=(X/2.)**2 |
|---|
| 186 | P=(((((-.00004686*T-.00110404)*T-.01919402)*T-.18156897)*T- |
|---|
| 187 | 1 .67278579)*T+.15443144)*T+1. |
|---|
| 188 | XK1=X*dLOG(X/2)*FI1+P |
|---|
| 189 | ! RETURN |
|---|
| 190 | else |
|---|
| 191 | ! 20 T=2./X |
|---|
| 192 | T=2./X |
|---|
| 193 | P=(((((-.00068245*T+.00325614)*T-.00780353)*T+.01504268)*T- |
|---|
| 194 | 1 .03655620)*T+.23498619)*T+1.25331414 |
|---|
| 195 | X=dMIN1(X,330.d0) |
|---|
| 196 | XK1=dSQRT(X)*dEXP(-X)*P |
|---|
| 197 | endif |
|---|
| 198 | RETURN |
|---|
| 199 | END |
|---|
| 200 | |
|---|
| 201 | |
|---|
| 202 | |
|---|
| 203 | |
|---|
| 204 | c ---------------------------------------------------------------------------------------------- |
|---|
| 205 | c Added by RDW for inclusion of Baranov (2004) data |
|---|
| 206 | subroutine baranov(temp,wn,abcoef) |
|---|
| 207 | ! computes the dimer spectrum |
|---|
| 208 | |
|---|
| 209 | implicit none |
|---|
| 210 | real*8 temp !temperature |
|---|
| 211 | real*8 wn !wavenumber |
|---|
| 212 | real*8 abcoef !absorption coefficient |
|---|
| 213 | |
|---|
| 214 | integer nS,nT |
|---|
| 215 | ! parameter(nS=1713) |
|---|
| 216 | parameter(nS=3000) |
|---|
| 217 | parameter(nT=9) |
|---|
| 218 | |
|---|
| 219 | real*8 wn_arr(nS) |
|---|
| 220 | real*8 temp_arr(nT) |
|---|
| 221 | real*8 dim_arr(nS,nT) |
|---|
| 222 | |
|---|
| 223 | integer pop |
|---|
| 224 | logical firstcall |
|---|
| 225 | |
|---|
| 226 | save wn_arr, temp_arr, dim_arr |
|---|
| 227 | |
|---|
| 228 | character*100 dt_file |
|---|
| 229 | integer strlen,ios |
|---|
| 230 | character*100 label |
|---|
| 231 | label='subroutine baranov' |
|---|
| 232 | |
|---|
| 233 | dt_file='/san/home/rdword/kspectrum/line_data/CO2_dimer_data_NEW' |
|---|
| 234 | ! dt_file='./data/CO2_dimer_data' |
|---|
| 235 | open(33,file=dt_file(1:strlen(dt_file)), |
|---|
| 236 | & form='unformatted', |
|---|
| 237 | & status='old',iostat=ios) |
|---|
| 238 | if (ios.ne.0) then ! file not found |
|---|
| 239 | write(*,*) 'Error from ',label(1:strlen(label)),' :' |
|---|
| 240 | write(*,*) 'Data file could not be found:' |
|---|
| 241 | write(*,*) dt_file(1:strlen(dt_file)) |
|---|
| 242 | stop |
|---|
| 243 | else |
|---|
| 244 | read(33) wn_arr |
|---|
| 245 | read(33) temp_arr |
|---|
| 246 | read(33) dim_arr |
|---|
| 247 | endif |
|---|
| 248 | close(33) |
|---|
| 249 | |
|---|
| 250 | |
|---|
| 251 | ! print*,wn_arr(1)! |
|---|
| 252 | ! print*,wn_arr(3000) |
|---|
| 253 | ! print*,temp_arr |
|---|
| 254 | ! print*,dim_arr(1500,1) |
|---|
| 255 | ! stop |
|---|
| 256 | |
|---|
| 257 | |
|---|
| 258 | call bilinear(wn_arr,temp_arr,nS,nT,dim_arr,wn,temp,abcoef) |
|---|
| 259 | |
|---|
| 260 | 111 continue |
|---|
| 261 | return |
|---|
| 262 | end |
|---|
| 263 | c ---------------------------------------------------------------------------------------------- |
|---|
| 264 | |
|---|
| 265 | c ---------------------------------------------------------------------------------------------- |
|---|
| 266 | c Added by RDW for inclusion of Baranov (2004) data |
|---|
| 267 | subroutine bilinear(x_arr,y_arr,nX,nY,f2d_arr,x,y,f) |
|---|
| 268 | |
|---|
| 269 | implicit none |
|---|
| 270 | |
|---|
| 271 | integer nX,nY,i,j,a,b |
|---|
| 272 | |
|---|
| 273 | real*8 x,y,x1,x2,y1,y2 |
|---|
| 274 | real*8 f,f11,f12,f21,f22,fA,fB |
|---|
| 275 | real*8 x_arr(nX) |
|---|
| 276 | real*8 y_arr(nY) |
|---|
| 277 | real*8 f2d_arr(nX,nY) |
|---|
| 278 | |
|---|
| 279 | integer strlen |
|---|
| 280 | character*100 label |
|---|
| 281 | label='subroutine bilinear' |
|---|
| 282 | |
|---|
| 283 | c 1st check we're within the wavenumber range |
|---|
| 284 | if ((x.lt.x_arr(2)).or.(x.gt.x_arr(nX-2))) then |
|---|
| 285 | f=0.0D+0 |
|---|
| 286 | c print*,'x=',x,x_Arr(2),x_arr(nX-2) |
|---|
| 287 | c stop('Outside CIA dimer wavenumber range!\n') |
|---|
| 288 | goto 112 |
|---|
| 289 | else |
|---|
| 290 | |
|---|
| 291 | c in the x (wavenumber) direction 1st |
|---|
| 292 | i=1 |
|---|
| 293 | 10 if (i.lt.(nX+1)) then ! what passes for a 'while' loop in this sorry language... |
|---|
| 294 | if (x_arr(i).gt.x) then |
|---|
| 295 | x1=x_arr(i-1) |
|---|
| 296 | x2=x_arr(i) |
|---|
| 297 | a=i-1 |
|---|
| 298 | i=9999 |
|---|
| 299 | endif |
|---|
| 300 | i=i+1 |
|---|
| 301 | goto 10 |
|---|
| 302 | endif |
|---|
| 303 | endif |
|---|
| 304 | |
|---|
| 305 | if ((y.lt.y_arr(1)).or.(y.gt.y_arr(nY))) then |
|---|
| 306 | write(*,*) 'Error from ',label(1:strlen(label)),' :' |
|---|
| 307 | write(*,*) 'Outside CIA dimer temperature range!' |
|---|
| 308 | if(y.lt.y_arr(1))then |
|---|
| 309 | y=y_arr(1)+0.01 |
|---|
| 310 | endif |
|---|
| 311 | if(y.gt.y_arr(nY))then |
|---|
| 312 | y=y_arr(nY)-0.01 |
|---|
| 313 | endif |
|---|
| 314 | ! stop |
|---|
| 315 | else |
|---|
| 316 | c in the y (temperature) direction 2nd |
|---|
| 317 | j=1 |
|---|
| 318 | 20 if (j.lt.(nY+1)) then |
|---|
| 319 | if (y_arr(j).gt.y) then |
|---|
| 320 | y1=y_arr(j-1) |
|---|
| 321 | y2=y_arr(j) |
|---|
| 322 | b=j-1 |
|---|
| 323 | j=9999 |
|---|
| 324 | endif |
|---|
| 325 | j=j+1 |
|---|
| 326 | goto 20 |
|---|
| 327 | endif |
|---|
| 328 | endif |
|---|
| 329 | |
|---|
| 330 | f11=f2d_arr(a,b) |
|---|
| 331 | f21=f2d_arr(a+1,b) |
|---|
| 332 | f12=f2d_arr(a,b+1) |
|---|
| 333 | f22=f2d_arr(a+1,b+1) |
|---|
| 334 | |
|---|
| 335 | c 1st in x-direction |
|---|
| 336 | fA=f11*(x2-x)/(x2-x1)+f21*(x-x1)/(x2-x1) |
|---|
| 337 | fB=f12*(x2-x)/(x2-x1)+f22*(x-x1)/(x2-x1) |
|---|
| 338 | |
|---|
| 339 | c then in y-direction |
|---|
| 340 | f=fA*(y2-y)/(y2-y1)+fB*(y-y1)/(y2-y1) |
|---|
| 341 | |
|---|
| 342 | 112 continue |
|---|
| 343 | return |
|---|
| 344 | end |
|---|
| 345 | c ---------------------------------------------------------------------------------------------- |
|---|
| 346 | |
|---|
