chem/rdxs.F

   1 ! This file contains subroutines related to reading the
   2 ! absorption cross sections of gases that contribute to atmospheric transmission:
   3 ! Some of these subroutines are also called from rxn.f when loading photolysis cross sections
   4 ! for these same gases. It is possible to have different cross sections for
   5 ! transmission and for photolysis, e.g. for ozone, Bass et al. could be used
   6 ! for transmission while Molina and Molina could be used for photolysis.
   7 ! This flexibility can be useful but users should be aware.
   8 ! For xsections that are temperature dependent, caution should be used in passing the proper
   9 ! temperature to the data routines.  Usually, transmission is for layers, TLAY(NZ-1), while
  10 ! photolysis is at levels, T(NZ).
  11 ! The following subroutines are her:
  12 !     rdo3xs
  13 !       o3_mol
  14 !       o3_rei
  15 !       o3_bas
  16 !       o3_wmo
  17 !       o3_jpl
  18 !     rdo2xs
  19 !     rdno2xs
  20 !       no2xs_d
  21 !       no2xs_jpl94
  22 !       no2xs_har
  23 !       no2xs_jpl06a
  24 !       no2xs_jpl06b
  25 !     rdso2xs
  26 !=============================================================================*
  27
  28       module module_xsections
  29
  30       use module_params
  31
  32       IMPLICIT NONE
  33
  34       public :: o3xs, rdo2xs, rdso2xs, no2xs_jpl06a
  35       public :: rdxs_init
  36       private
  37
  38       REAL, allocatable :: rei218(:), rei228(:), rei243(:), rei295(:)
  39       REAL :: v195, v345, v830
  40       REAL, allocatable :: wmo203(:), wmo273(:)
  41       REAL :: v176, v850
  42
  43       REAL, allocatable :: jpl295(:), jpl218(:)
  44       REAL :: v186, v825
  45
  46       REAL, allocatable :: mol226(:), mol263(:), mol298(:)
  47       REAL :: v185, v240, v350
  48
  49       REAL, allocatable :: c0(:), c1(:), c2(:)
  50       REAL vb245, vb342
  51
  52       REAL, allocatable :: no2xs_a(:), no2xs_b(:)
  53
  54       CONTAINS
  55
  56       SUBROUTINE rdxs_init( nw, wl )
  57
  58       integer, intent(in) :: nw
  59       real, intent(in)    :: wl(nw)
  60
  61       integer :: istat, astat
  62
  63       istat = 0
  64       if( .not. allocated( rei218 ) ) then
  65         allocate( rei218(nw-1),rei228(nw-1),rei243(nw-1),rei295(nw-1),stat=astat )
  66         istat = istat + astat
  67       endif
  68       if( .not. allocated( wmo203 ) ) then
  69         allocate( wmo203(nw-1),wmo273(nw-1),stat=astat )
  70         istat = istat + astat
  71       endif
  72       if( .not. allocated( jpl218 ) ) then
  73         allocate( jpl218(nw-1),jpl295(nw-1),stat=astat )
  74         istat = istat + astat
  75       endif
  76       if( .not. allocated( mol226 ) ) then
  77         allocate( mol226(nw-1),mol263(nw-1),mol298(nw-1),stat=astat )
  78         istat = istat + astat
  79       endif
  80       if( .not. allocated( c0 ) ) then
  81         allocate( c0(nw-1),c1(nw-1),c2(nw-1),stat=astat )
  82         istat = istat + astat
  83       endif
  84       if( .not. allocated( no2xs_a ) ) then
  85         allocate( no2xs_a(nw-1),no2xs_b(nw-1),stat=astat )
  86         istat = istat + astat
  87       endif
  88
  89 !_______________________________________________________________________
  90 ! read data from different sources
  91 ! rei = Reims group (Malicet et al., Brion et al.)
  92 ! jpl = JPL 2006 evaluation
  93 ! wmo = WMO 1985 O3 assessment
  94 ! mol = Molina and Molina
  95 ! bas = Bass et al.
  96 !_______________________________________________________________________
  97       CALL o3_rei(nw,wl)
  98       CALL o3_jpl(nw,wl)
  99       CALL o3_wmo(nw,wl)
 100       CALL o3_mol(nw,wl)
 101       CALL o3_bas(nw,wl)
 102
 103       CALL rdno2xs(nw,wl)
 104
 105       END SUBROUTINE rdxs_init
 106
 107       SUBROUTINE o3xs(nz,t,nw,wl, xs)
 108 !-----------------------------------------------------------------------------*
 109 !=  PURPOSE:                                                                 =*
 110 !=  Read ozone molecular absorption cross section.  Re-grid data to match    =*
 111 !=  specified wavelength working grid. Interpolate in temperature as needed  =*
 112 !-----------------------------------------------------------------------------*
 113 !=  PARAMETERS:                                                              =*
 114 !=  MABS   - INTEGER, option for splicing different combinations of       (I)=*
 115 !=           absorption cross secttions                                      =*
 116 !=  NZ     - INTEGER, number of altitude levels or layers                 (I)=*
 117 !=  T      - REAL, temperature of levels or layers                        (I)=*
 118 !=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
 119 !=           wavelength grid                                                 =*
 120 !=  WL     - REAL, vector of lower limits of wavelength intervals in      (I)=*
 121 !=           working wavelength grid. In vacuum, nm                          =*
 122 !=  XS     - REAL, molecular absoprtion cross section (cm^2) of O3 at     (O)=*
 123 !=           each specified wavelength (WMO value at 273)                    =*
 124 !-----------------------------------------------------------------------------*
 125
 126 ! input: (altitude working grid)
 127
 128       INTEGER, intent(in) :: nw
 129       REAL, intent(in)    :: wl(nw)
 130
 131       INTEGER, intent(in) :: nz
 132       REAL, intent(in) :: t(nz)
 133
 134 ! output:
 135 ! ozone absorption cross sections interpolated to
 136 !   working wavelength grid (iw)
 137 !   working altitude grid (iz) for temperature of layer or level (specified in call)
 138 ! Units are cm2 molecule-1 in vacuum
 139
 140       REAL, intent(inout) :: xs(:,:)
 141
 142 ! internal
 143
 144       INTEGER :: iz, iw
 145       REAL    :: factor, factor1
 146       REAL    :: tc(nz)
 147
 148 !***** option 1:
 149 ! assign according to wavelength range:
 150 !  175.439 - 185.185  1985WMO (203, 273 K)
 151 !  185.185 - 195.00   2006JPL_O3 (218, 295 K)
 152 !  195.00  - 345.00   Reims group (218, 228, 243, 295 K)
 153 !  345.00  - 830.00   Reims group (295 K)
 154 !  no extrapolations in temperature allowed
 155
 156       DO iw = 1, nw-1
 157         IF(wl(iw) < v185) THEN
 158           factor = (wmo273(iw) - wmo203(iw))/(273. - 203.)
 159           xs(1:nz,iw) = wmo203(iw) + (t(1:nz) - 203.)*factor
 160           WHERE (t(1:nz) <= 203.)
 161             xs(1:nz,iw) = wmo203(iw)
 162           ELSEWHERE (t(1:nz) >= 273.)
 163             xs(1:nz,iw) = wmo273(iw)
 164           ENDWHERE
 165         ELSEIF(wl(iw) >= v185 .AND. wl(iw) < v195) THEN
 166           factor = (jpl295(iw) - jpl218(iw))/(295. - 218.)
 167           xs(1:nz,iw) = jpl218(iw) + (t(1:nz) - 218.)*factor
 168           WHERE (t(1:nz) <= 218.)
 169             xs(1:nz,iw) = jpl218(iw)
 170           ELSEWHERE (t(1:nz) >= 295.)
 171             xs(1:nz,iw) = jpl295(iw)
 172           ENDWHERE
 173         ELSEIF(wl(iw) >= v195 .AND. wl(iw) < v345) THEN
 174           factor = .1*(rei228(iw) - rei218(iw))
 175           WHERE( t(1:nz) < 218. )
 176             xs(1:nz,iw) = rei218(iw)
 177           ELSEWHERE( t(1:nz) >= 218. .AND. t(1:nz) < 228. )
 178             xs(1:nz,iw) = rei218(iw) + (t(1:nz) - 218.)*factor
 179           ENDWHERE
 180           factor = (rei243(iw) - rei228(iw))/15.
 181           WHERE( t(1:nz) >= 228. .AND. t(1:nz) < 243. )
 182             xs(1:nz,iw) = rei228(iw) + (t(1:nz) - 228.)*factor
 183           ENDWHERE
 184           factor = (rei295(iw) - rei243(iw))/(295. - 243.)
 185           WHERE( t(1:nz) >= 243. .AND. t(1:nz) < 295.)
 186             xs(1:nz,iw) = rei243(iw) + (t(1:nz) - 243.)*factor
 187           ELSEWHERE( t(1:nz) >= 295. )
 188             xs(1:nz,iw) = rei295(iw)
 189           ENDWHERE
 190         ELSEIF(wl(iw) >= v345) THEN
 191           xs(1:nz,iw) = rei295(iw)
 192         ENDIF
 193       END DO
 194
 195       END SUBROUTINE o3xs
 196
 197 !=============================================================================*
 198
 199       SUBROUTINE o3_rei(nw,wl)
 200 !-----------------------------------------------------------------------------*
 201 !=  PURPOSE:                                                                 =*
 202 !=  Read and interpolate the O3 cross section from Reims group               =*
 203 !-----------------------------------------------------------------------------*
 204 !=  PARAMETERS:                                                              =*
 205 !=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
 206 !=           wavelength grid                                                 =*
 207 !=  WL     - REAL, vector of lower limits of wavelength intervals in      (I)=*
 208 !=           working wavelength grid                                         =*
 209 !=  REI218 - REAL, cross section (cm^2) for O3 at 218K                    (O)=*
 210 !=  REI228 - REAL, cross section (cm^2) for O3 at 218K                    (O)=*
 211 !=  REI243 - REAL, cross section (cm^2) for O3 at 218K                    (O)=*
 212 !=  REI295 - REAL, cross section (cm^2) for O3 at 218K                    (O)=*
 213 !=  V195   - REAL, exact wavelength in vacuum for data breaks             (O)=*
 214 !=              e.g. start, stop, or other change                            =*
 215 !=  V345   - REAL, exact wavelength in vacuum for data breaks             (O)=*
 216 !=  V830   - REAL, exact wavelength in vacuum for data breaks             (O)=*
 217 !-----------------------------------------------------------------------------*
 218
 219 !  input
 220
 221       INTEGER, intent(in) :: nw
 222       REAL, intent(in)    :: wl(nw)
 223
 224 !* internal
 225
 226       INTEGER, PARAMETER :: kdata = 70000
 227
 228       INTEGER n1, n2, n3, n4, iw
 229       REAL x1(kdata), x2(kdata), x3(kdata), x4(kdata)
 230       REAL y1(kdata), y2(kdata), y3(kdata), y4(kdata)
 231
 232       INTEGER i
 233       INTEGER ierr
 234
 235 ! used for air-to-vacuum wavelength conversion
 236
 237       REAL ri(kdata)
 238       CHARACTER(len=256) :: emsg
 239
 240 ! data from the Reims group:
 241 !=  For Hartley and Huggins bands, use temperature-dependent values from     =*
 242 !=  Malicet et al., J. Atmos. Chem.  v.21, pp.263-273, 1995.                 =*
 243 !=  over 345.01 - 830.00, use values from Brion, room temperature only
 244
 245       OPEN(UNIT=kin,FILE='DATAE1/O3/1995Malicet_O3.txt',STATUS='old',iostat=ierr)
 246       if( ierr /= 0 ) then
 247         call wrf_error_fatal( 'o3_rei: Failed to open DATAE1/O3/1985Malicet_O3.txt' )
 248       endif
 249       DO i = 1, 2
 250          READ(kin,*,iostat=ierr)
 251          if( ierr /= 0 ) then
 252            call wrf_error_fatal( 'o3_rei: Failed to read DATAE1/O3/1985Malicet_O3.txt' )
 253          endif
 254       ENDDO
 255       n1 = 15001
 256       n2 = 15001
 257       n3 = 15001
 258       n4 = 15001
 259       DO i = 1, n1
 260          READ(kin,*,iostat=ierr) x1(i), y1(i), y2(i), y3(i), y4(i)
 261          if( ierr /= 0 ) then
 262            call wrf_error_fatal( 'o3_rei: Failed to read DATAE1/O3/1985Malicet_O3.txt' )
 263          endif
 264          x2(i) = x1(i)
 265          x3(i) = x1(i)
 266          x4(i) = x1(i)
 267       ENDDO
 268       CLOSE (kin)
 269
 270 !=  over 345.01 - 830.00, use values from Brion, room temperature only
 271 ! skip datum at 345.00 because already read in from 1995Malicet
 272
 273       OPEN(UNIT=kin,FILE='DATAE1/O3/1998Brion_295.txt',STATUS='old')
 274       DO i = 1, 15
 275          READ(kin,*)
 276       ENDDO
 277       DO i = 1, 48515-15
 278          n1 = n1 + 1
 279          READ(kin,*) x1(n1), y1(n1)
 280       ENDDO
 281       CLOSE (kin)
 282
 283       DO i = 1, n1
 284          ri(i) = refrac(x1(i), 2.45E19)
 285       ENDDO
 286       DO i = 1, n1
 287          x1(i) = x1(i) * ri(i)
 288       ENDDO
 289
 290       CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax),0.)
 291       CALL addpnt(x1,y1,kdata,n1,               0.,0.)
 292       CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.)
 293       CALL addpnt(x1,y1,kdata,n1,            1.e+38,0.)
 294       CALL inter2(nw,wl,rei295,n1,x1,y1,ierr)
 295       IF (ierr .NE. 0) THEN
 296          WRITE(emsg,'(''o3_rei: interp err = '',i5,'' in O3 xsect - Reims 295K'')') ierr
 297          call wrf_error_fatal( trim(emsg) )
 298       ENDIF
 299
 300       DO i = 1, n2
 301          ri(i) = refrac(x2(i), 2.45E19)
 302       ENDDO
 303       DO i = 1, n2
 304          x2(i) = x2(i) * ri(i)
 305          x3(i) = x2(i)
 306          x4(i) = x2(i)
 307       ENDDO
 308
 309       CALL addpnt(x2,y2,kdata,n2,x2(1)*(1.-deltax),0.)
 310       CALL addpnt(x2,y2,kdata,n2,               0.,0.)
 311       CALL addpnt(x2,y2,kdata,n2,x2(n2)*(1.+deltax),0.)
 312       CALL addpnt(x2,y2,kdata,n2,            1.e+38,0.)
 313       CALL inter2(nw,wl,rei243,n2,x2,y2,ierr)
 314       IF (ierr .NE. 0) THEN
 315          WRITE(emsg,'(''o3_wmo: interp err = '',i5,'' in O3 xsect - Reims 243K'')') ierr
 316          call wrf_error_fatal( trim(emsg) )
 317       ENDIF
 318
 319       CALL addpnt(x3,y3,kdata,n3,x3(1)*(1.-deltax),0.)
 320       CALL addpnt(x3,y3,kdata,n3,               0.,0.)
 321       CALL addpnt(x3,y3,kdata,n3,x3(n3)*(1.+deltax),0.)
 322       CALL addpnt(x3,y3,kdata,n3,            1.e+38,0.)
 323       CALL inter2(nw,wl,rei228,n3,x3,y3,ierr)
 324       IF (ierr .NE. 0) THEN
 325          WRITE(emsg,'(''o3_wmo: interp err = '',i5,'' in O3 xswct - Reims 228K'')') ierr
 326          call wrf_error_fatal( trim(emsg) )
 327       ENDIF
 328
 329       CALL addpnt(x4,y4,kdata,n4,x4(1)*(1.-deltax),0.)
 330       CALL addpnt(x4,y4,kdata,n4,               0.,0.)
 331       CALL addpnt(x4,y4,kdata,n4,x4(n4)*(1.+deltax),0.)
 332       CALL addpnt(x4,y4,kdata,n4,            1.e+38,0.)
 333       CALL inter2(nw,wl,rei218,n4,x4,y4,ierr)
 334       IF (ierr .NE. 0) THEN
 335          WRITE(emsg,'(''o3_wmo: interp err = '',i5,'' in O3 xswct - Reims 218K'')') ierr
 336          call wrf_error_fatal( trim(emsg) )
 337       ENDIF
 338
 339 ! wavelength breaks must be converted to vacuum:
 340
 341       v195 = 195.00 * refrac(195.00, 2.45E19)
 342       v345 = 345.00 * refrac(345.00, 2.45E19)
 343       v830 = 830.00 * refrac(830.00, 2.45E19)
 344
 345       END SUBROUTINE o3_rei
 346
 347 !=============================================================================*
 348
 349       SUBROUTINE o3_wmo(nw,wl)
 350 !-----------------------------------------------------------------------------*
 351 !=  PURPOSE:                                                                 =*
 352 !=  Read and interpolate the O3 cross section                                =*
 353 !=  data from WMO 85 Ozone Assessment                                        =*
 354 !-----------------------------------------------------------------------------*
 355 !=  PARAMETERS:                                                              =*
 356 !=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
 357 !=           wavelength grid                                                 =*
 358 !=  WL     - REAL, vector of lower limits of wavelength intervals in      (I)=*
 359 !=           working wavelength grid                                         =*
 360 !=  WMO203 - REAL, cross section (cm^2) for O3 at 203K                    (O)=*
 361 !=  WMO273 - REAL, cross section (cm^2) for O3 at 273K                    (O)=*
 362 !=  V176   - REAL, exact wavelength in vacuum for data breaks             (O)=*
 363 !=              e.g. start, stop, or other change                            =*
 364 !=  V850   - REAL, exact wavelength in vacuum for data breaks             (O)=*
 365 !-----------------------------------------------------------------------------*
 366
 367 !  input
 368
 369       INTEGER, intent(in) :: nw
 370       REAL, intent(in)    :: wl(nw)
 371
 372 ! internal
 373
 374       INTEGER, parameter :: kdata = 200
 375
 376       INTEGER n1, n2
 377       REAL x1(kdata), x2(kdata)
 378       REAL y1(kdata), y2(kdata)
 379
 380       INTEGER i, idum, iw
 381       REAL a1, a2, dum
 382       INTEGER ierr
 383
 384 ! used for air-to-vacuum wavelength conversion
 385
 386       REAL ri(kdata)
 387       CHARACTER(len=256) :: emsg
 388
 389 ! output
 390
 391 !----------------------------------------------------------
 392 ! cross sections from WMO 1985 Ozone Assessment
 393 ! from 175.439 to 847.500 nm
 394
 395       OPEN(UNIT=kin,FILE='DATAE1/wmo85',STATUS='old',iostat=ierr)
 396       if( ierr /= 0 ) then
 397         call wrf_error_fatal( 'o3_wmo: Failed to open DATAE1/wmo85' )
 398       endif
 399       DO i = 1, 3
 400          read(kin,*,iostat=ierr)
 401          if( ierr /= 0 ) then
 402            call wrf_error_fatal( 'o3_wmo: Failed to read DATAE1/wmo85' )
 403          endif
 404       ENDDO
 405       n1 = 158
 406       n2 = 158
 407       DO i = 1, n1
 408          READ(kin,*,iostat=ierr) idum, a1, a2, dum, dum, dum, y1(i), y2(i)
 409          if( ierr /= 0 ) then
 410            call wrf_error_fatal( 'o3_wmo: Failed to read DATAE1/wmo85' )
 411          endif
 412          x1(i) = (a1+a2)/2.
 413          x2(i) = (a1+a2)/2.
 414       ENDDO
 415       CLOSE (kin)
 416
 417 ! convert wavelengths to vacuum
 418
 419       DO i = 1, n1
 420          ri(i) = refrac(x1(i), 2.45E19)
 421       ENDDO
 422       DO i = 1, n1
 423          x1(i) = x1(i) * ri(i)
 424          x2(i) = x2(i) * ri(i)
 425       ENDDO
 426
 427       CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax),0.)
 428       CALL addpnt(x1,y1,kdata,n1,               0.,0.)
 429       CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.)
 430       CALL addpnt(x1,y1,kdata,n1,           1.e+38,0.)
 431       CALL inter2(nw,wl,wmo203,n1,x1,y1,ierr)
 432       IF (ierr .NE. 0) THEN
 433          WRITE(emsg,'(''o3_wmo: interp err = '',i5,'' in O3 cross section - WMO - 203K'')') ierr
 434          call wrf_error_fatal( trim(emsg) )
 435       ENDIF
 436
 437       CALL addpnt(x2,y2,kdata,n2,x2(1)*(1.-deltax),0.)
 438       CALL addpnt(x2,y2,kdata,n2,               0.,0.)
 439       CALL addpnt(x2,y2,kdata,n2,x2(n2)*(1.+deltax),0.)
 440       CALL addpnt(x2,y2,kdata,n2,           1.e+38,0.)
 441       CALL inter2(nw,wl,wmo273,n2,x2,y2,ierr)
 442       IF (ierr .NE. 0) THEN
 443          WRITE(emsg,'(''o3_wmo: interp err = '',i5,'' in O3 cross section - WMO - 273K'')') ierr
 444          call wrf_error_fatal( trim(emsg) )
 445       ENDIF
 446
 447 ! wavelength breaks must be converted to vacuum:
 448
 449       a1 = (175.438 + 176.991) / 2.
 450       v176 = a1 * refrac(a1,2.45E19)
 451
 452       a1 = (847.5 + 852.5) / 2.
 453       v850 = a1 * refrac(a1, 2.45E19)
 454
 455       END SUBROUTINE o3_wmo
 456
 457 !=============================================================================*
 458
 459       SUBROUTINE o3_jpl(nw,wl)
 460 !-----------------------------------------------------------------------------*
 461 !=  PURPOSE:                                                                 =*
 462 !=  Read and interpolate the O3 cross section from JPL 2006                  =*
 463 !-----------------------------------------------------------------------------*
 464 !=  PARAMETERS:                                                              =*
 465 !=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
 466 !=           wavelength grid                                                 =*
 467 !=  WL     - REAL, vector of lower limits of wavelength intervals in      (I)=*
 468 !=           working wavelength grid                                         =*
 469 !=  JPL218 - REAL, cross section (cm^2) for O3 at 218K                    (O)=*
 470 !=  JPL295 - REAL, cross section (cm^2) for O3 at 295K                    (O)=*
 471 !=  V186   - REAL, exact wavelength in vacuum for data breaks             (O)=*
 472 !=              e.g. start, stop, or other change                            =*
 473 !=  V825   - REAL, exact wavelength in vacuum for data breaks             (O)=*
 474 !-----------------------------------------------------------------------------*
 475
 476 !  input
 477
 478       INTEGER, intent(in) :: nw
 479       REAL, intent(in)    :: wl(nw)
 480
 481 ! internal
 482
 483       INTEGER, parameter :: kdata = 200
 484
 485       INTEGER n1, n2, iw
 486       REAL x1(kdata), x2(kdata)
 487       REAL y1(kdata), y2(kdata)
 488
 489       INTEGER i
 490       REAL dum
 491       INTEGER ierr
 492       CHARACTER(len=256) :: emsg
 493
 494 ! used for air-to-vacuum wavelength conversion
 495
 496       REAL ri(kdata)
 497
 498 ! output
 499
 500       OPEN(UNIT=kin,FILE='DATAE1/O3/2006JPL_O3.txt',STATUS='old',iostat=ierr)
 501       if( ierr /= 0 ) then
 502         call wrf_error_fatal( 'o3_jpl: Failed to open DATAE1/O3/2006JPL_O3.txt' )
 503       endif
 504       DO i = 1, 2
 505          read(kin,*,iostat=ierr)
 506          if( ierr /= 0 ) then
 507            call wrf_error_fatal( 'o3_jpl: Failed to read DATAE1/O3/2006JPL_O3.txt' )
 508          endif
 509       ENDDO
 510       n1 = 167
 511       n2 = 167
 512       DO i = 1, n1
 513          READ(kin,*,iostat=ierr) dum, dum, x1(i), y1(i), y2(i)
 514          if( ierr /= 0 ) then
 515            call wrf_error_fatal( 'o3_jpl: Failed to read DATAE1/O3/2006JPL_O3.txt' )
 516          endif
 517          y1(i) = y1(i) * 1.e-20
 518          y2(i) = y2(i) * 1.e-20
 519       ENDDO
 520       CLOSE (kin)
 521
 522 ! convert wavelengths to vacuum
 523
 524       DO i = 1, n1
 525          ri(i) = refrac(x1(i), 2.45E19)
 526       ENDDO
 527       DO i = 1, n1
 528          x1(i) = x1(i) * ri(i)
 529          x2(i) = x1(i)
 530       ENDDO
 531
 532       CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax),0.)
 533       CALL addpnt(x1,y1,kdata,n1,               0.,0.)
 534       CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.)
 535       CALL addpnt(x1,y1,kdata,n1,           1.e+38,0.)
 536       CALL inter2(nw,wl,jpl295,n1,x1,y1,ierr)
 537       IF (ierr .NE. 0) THEN
 538          WRITE(emsg,'(''o3_jpl: interp err = '',i5,'' in file O3 cross section - WMO - 295K'')') ierr
 539          call wrf_error_fatal( trim(emsg) )
 540       ENDIF
 541
 542       CALL addpnt(x2,y2,kdata,n2,x2(1)*(1.-deltax),0.)
 543       CALL addpnt(x2,y2,kdata,n2,               0.,0.)
 544       CALL addpnt(x2,y2,kdata,n2,x2(n2)*(1.+deltax),0.)
 545       CALL addpnt(x2,y2,kdata,n2,           1.e+38,0.)
 546       CALL inter2(nw,wl,jpl218,n2,x2,y2,ierr)
 547       IF (ierr .NE. 0) THEN
 548          WRITE(emsg,'(''o3_jpl: interp err = '',i5,'' in file O3 cross section - WMO - 218K'')') ierr
 549          call wrf_error_fatal( trim(emsg) )
 550       ENDIF
 551
 552 ! wavelength breaks must be converted to vacuum:
 553
 554       v186 = 186.051 * refrac(186.051, 2.45E19)
 555       v825 = 825.    * refrac(825.   , 2.45E19)
 556
 557
 558       END SUBROUTINE o3_jpl
 559
 560 !=============================================================================*
 561
 562       SUBROUTINE o3_mol(nw,wl)
 563 !-----------------------------------------------------------------------------*
 564 !=  PURPOSE:                                                                 =*
 565 !=  Read and interpolate the O3 cross section from Molina and Molina 1986    =*
 566 !-----------------------------------------------------------------------------*
 567 !=  PARAMETERS:                                                              =*
 568 !=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
 569 !=           wavelength grid                                                 =*
 570 !=  WL     - REAL, vector of lower limits of wavelength intervals in      (I)=*
 571 !=           working wavelength grid                                         =*
 572 !=  MOL226 - REAL, cross section (cm^2) for O3 at 226 K                   (O)=*
 573 !=  MOL263 - REAL, cross section (cm^2) for O3 at 263 K                   (O)=*
 574 !=  MOL298 - REAL, cross section (cm^2) for O3 at 298 K                   (O)=*
 575 !=  V185   - REAL, exact wavelength in vacuum for data breaks             (O)=*
 576 !=              e.g. start, stop, or other change                            =*
 577 !=  V240   - REAL, exact wavelength in vacuum for data breaks             (O)=*
 578 !=  V350   - REAL, exact wavelength in vacuum for data breaks             (O)=*
 579 !-----------------------------------------------------------------------------*
 580
 581 !  input
 582
 583       INTEGER, intent(in) :: nw
 584       REAL, intent(in)    :: wl(nw)
 585
 586 ! internal
 587
 588       INTEGER i
 589       INTEGER ierr
 590
 591       INTEGER, parameter :: kdata = 335
 592       INTEGER n1, n2, n3, iw
 593       REAL x1(kdata), x2(kdata), x3(kdata)
 594       REAL y1(kdata), y2(kdata), y3(kdata)
 595
 596 ! used for air-to-vacuum wavelength conversion
 597
 598       REAL ri(kdata)
 599       CHARACTER(len=256) :: emsg
 600
 601 ! output
 602
 603       OPEN(UNIT=kin,FILE='DATAE1/O3/1986Molina.txt',STATUS='old',iostat=ierr)
 604       if( ierr /= 0 ) then
 605         call wrf_error_fatal( 'o3_mol: Failed to open DATAE1/O3/1986Molina.txt' )
 606       endif
 607       DO i = 1, 10
 608         READ(kin,*,iostat=ierr)
 609         if( ierr /= 0 ) then
 610           call wrf_error_fatal( 'o3_mol: Failed to read DATAE1/O3/1986Molina.txt' )
 611         endif
 612       ENDDO
 613       n1 = 0
 614       n2 = 0
 615       n3 = 0
 616       DO i = 1, 121-10
 617          n1 = n1 + 1
 618          n3 = n3 + 1
 619          READ(kin,*,iostat=ierr) x1(n1), y1(n1),  y3(n3)
 620          if( ierr /= 0 ) then
 621            call wrf_error_fatal( 'o3_mol: Failed to read DATAE1/O3/1986Molina.txt' )
 622          endif
 623          x3(n3) = x1(n1)
 624       ENDDO
 625       DO i = 1, 341-122
 626          n1 = n1 + 1
 627          n2 = n2 + 1
 628          n3 = n3 + 1
 629          READ(kin,*,iostat=ierr) x1(n1), y1(n1), y2(n2), y3(n3)
 630          if( ierr /= 0 ) then
 631            call wrf_error_fatal( 'o3_mol: Failed to read DATAE1/O3/1986Molina.txt' )
 632          endif
 633          x2(n2) = x1(n1)
 634          x3(n3) = x1(n1)
 635       ENDDO
 636       CLOSE (kin)
 637
 638 ! convert all wavelengths from air to vacuum
 639
 640       DO i = 1, n1
 641          ri(i) = refrac(x1(i), 2.45E19)
 642       ENDDO
 643       DO i = 1, n1
 644          x1(i) = x1(i) * ri(i)
 645       ENDDO
 646
 647       DO i = 1, n2
 648          ri(i) = refrac(x2(i), 2.45E19)
 649       ENDDO
 650       DO i = 1, n2
 651          x2(i) = x2(i) * ri(i)
 652       ENDDO
 653
 654       DO i = 1, n3
 655          ri(i) = refrac(x3(i), 2.45E19)
 656       ENDDO
 657       DO i = 1, n3
 658          x3(i) = x3(i) * ri(i)
 659       ENDDO
 660
 661 ! convert wavelength breaks from air to vacuum
 662
 663       v185 = 185.  * refrac(185. , 2.45E19)
 664       v240 = 240.5 * refrac(240.5, 2.45E19)
 665       v350 = 350.  * refrac(350. , 2.45E19)
 666
 667 ! interpolate to working grid
 668
 669       CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax),0.)
 670       CALL addpnt(x1,y1,kdata,n1,               0.,0.)
 671       CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.)
 672       CALL addpnt(x1,y1,kdata,n1,            1.e+38,0.)
 673       CALL inter2(nw,wl,mol226,n1,x1,y1,ierr)
 674       IF (ierr .NE. 0) THEN
 675          WRITE(emsg,'(''o3_mol: interp err = '',i5,'' in O3 xsect - 226K Molina'')') ierr
 676          call wrf_error_fatal( trim(emsg) )
 677       ENDIF
 678
 679       CALL addpnt(x2,y2,kdata,n2,x2(1)*(1.-deltax),0.)
 680       CALL addpnt(x2,y2,kdata,n2,               0.,0.)
 681       CALL addpnt(x2,y2,kdata,n2,x2(n2)*(1.+deltax),0.)
 682       CALL addpnt(x2,y2,kdata,n2,            1.e+38,0.)
 683       CALL inter2(nw,wl,mol263,n2,x2,y2,ierr)
 684       IF (ierr .NE. 0) THEN
 685          WRITE(emsg,'(''o3_mol: interp err = '',i5,'' in O3 xsect - 263K Molina'')') ierr
 686          call wrf_error_fatal( trim(emsg) )
 687       ENDIF
 688
 689       CALL addpnt(x3,y3,kdata,n3,x3(1)*(1.-deltax),0.)
 690       CALL addpnt(x3,y3,kdata,n3,               0.,0.)
 691       CALL addpnt(x3,y3,kdata,n3,x3(n3)*(1.+deltax),0.)
 692       CALL addpnt(x3,y3,kdata,n3,            1.e+38,0.)
 693       CALL inter2(nw,wl,mol298,n3,x3,y3,ierr)
 694       IF (ierr .NE. 0) THEN
 695          WRITE(emsg,'(''o3_mol: interp err = '',i5,'' in O3 xsect - 298K Molina'')') ierr
 696          call wrf_error_fatal( trim(emsg) )
 697       ENDIF
 698
 699       END SUBROUTINE o3_mol
 700
 701 !=============================================================================*
 702
 703       SUBROUTINE o3_bas(nw,wl)
 704 !-----------------------------------------------------------------------------*
 705 !=  PURPOSE:                                                                 =*
 706 !=  Read and interpolate the O3 cross section from Bass 1985                 =*
 707 !-----------------------------------------------------------------------------*
 708 !=  PARAMETERS:                                                              =*
 709 !=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
 710 !=           wavelength grid                                                 =*
 711 !=  WL     - REAL, vector of lower limits of wavelength intervals in      (I)=*
 712 !=           working wavelength grid                                         =*
 713 !=  c0     - REAL, coefficint for polynomial fit to cross section (cm^2)  (O)=*
 714 !=  c1     - REAL, coefficint for polynomial fit to cross section (cm^2)  (O)=*
 715 !=  c2     - REAL, coefficint for polynomial fit to cross section (cm^2)  (O)=*
 716 !=  Vb245   - REAL, exact wavelength in vacuum for data breaks            (O)=*
 717 !=              e.g. start, stop, or other change                            =*
 718 !=  Vb342   - REAL, exact wavelength in vacuum for data breaks            (O)=*
 719 !-----------------------------------------------------------------------------*
 720
 721 ! input:
 722
 723       INTEGER, intent(in) :: nw
 724       REAL, intent(in)    :: wl(nw)
 725
 726 ! internal:
 727
 728       INTEGER, parameter :: kdata = 2000
 729
 730       INTEGER i, iw
 731       INTEGER ierr
 732
 733       INTEGER n1, n2, n3
 734       REAL x1(kdata), x2(kdata), x3(kdata)
 735       REAL y1(kdata), y2(kdata), y3(kdata)
 736
 737 ! used for air-to-vacuum wavelength conversion
 738
 739       REAL ri(kdata)
 740       CHARACTER(len=256) :: emsg
 741
 742       OPEN(UNIT=kin,FILE='DATAE1/O3/1985Bass_O3.txt',STATUS='old',iostat=ierr)
 743       if( ierr /= 0 ) then
 744         call wrf_error_fatal( 'o3_bas: Failed to open DATAE1/O3/1985Bass_O3.txt' )
 745       endif
 746       DO i = 1, 8
 747          READ(kin,*,iostat=ierr)
 748       ENDDO
 749       n1 = 1915
 750       n2 = 1915
 751       n3 = 1915
 752       DO i = 1, n1
 753         READ(kin,*) x1(i), y1(i), y2(i), y3(i)
 754         if( ierr /= 0 ) then
 755           call wrf_error_fatal( 'o3_bas: Failed to read DATAE1/O3/1985Bass_O3.txt' )
 756         endif
 757       ENDDO
 758       CLOSE (kin)
 759       y1(1:n1) = 1.e-20 * y1(1:n1)
 760       y2(1:n1) = 1.e-20 * y2(1:n1)
 761       y3(1:n1) = 1.e-20 * y3(1:n1)
 762
 763 ! convert all wavelengths from air to vacuum
 764
 765       DO i = 1, n1
 766          ri(i) = refrac(x1(i), 2.45E19)
 767       ENDDO
 768       x1(1:n1) = x1(1:n1) * ri(1:n1)
 769       x2(1:n1) = x1(1:n1)
 770       x3(1:n1) = x1(1:n1)
 771
 772 ! convert wavelength breaks to vacuum
 773
 774       vb245 = 245.018 * refrac(245.018, 2.45E19)
 775       vb342 = 341.981 * refrac(341.981, 2.45E19)
 776
 777 ! interpolate to working grid
 778
 779       CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax),0.)
 780       CALL addpnt(x1,y1,kdata,n1,               0.,0.)
 781       CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),0.)
 782       CALL addpnt(x1,y1,kdata,n1,            1.e+38,0.)
 783       CALL inter2(nw,wl,c0,n1,x1,y1,ierr)
 784       IF (ierr .NE. 0) THEN
 785          WRITE(emsg,'(''o3_bas: interp err = '',i5,'' in O3 xsect - c0 Bass'')') ierr
 786          call wrf_error_fatal( trim(emsg) )
 787       ENDIF
 788
 789       CALL addpnt(x2,y2,kdata,n2,x2(1)*(1.-deltax),0.)
 790       CALL addpnt(x2,y2,kdata,n2,               0.,0.)
 791       CALL addpnt(x2,y2,kdata,n2,x2(n2)*(1.+deltax),0.)
 792       CALL addpnt(x2,y2,kdata,n2,            1.e+38,0.)
 793       CALL inter2(nw,wl,c1,n2,x2,y2,ierr)
 794       IF (ierr .NE. 0) THEN
 795          WRITE(emsg,'(''o3_bas: interp err = '',i5,'' in O3 xsect - c1 Bass'')') ierr
 796          call wrf_error_fatal( trim(emsg) )
 797       ENDIF
 798
 799       CALL addpnt(x3,y3,kdata,n3,x3(1)*(1.-deltax),0.)
 800       CALL addpnt(x3,y3,kdata,n3,               0.,0.)
 801       CALL addpnt(x3,y3,kdata,n3,x3(n3)*(1.+deltax),0.)
 802       CALL addpnt(x3,y3,kdata,n3,            1.e+38,0.)
 803       CALL inter2(nw,wl,c2,n3,x3,y3,ierr)
 804       IF (ierr .NE. 0) THEN
 805          WRITE(emsg,'(''o3_bas: interp err = '',i5,'' in O3 xsect - c2 Bass'')') ierr
 806          call wrf_error_fatal( trim(emsg) )
 807       ENDIF
 808
 809       END SUBROUTINE o3_bas
 810
 811 !=============================================================================*
 812
 813       SUBROUTINE rdo2xs(nw,wl,o2xs1)
 814 !-----------------------------------------------------------------------------*
 815 !=  PURPOSE:                                                                 =*
 816 !=  Compute equivalent O2 cross section, except                              =*
 817 !=  the SR bands and the Lyman-alpha line.                                   =*
 818 !-----------------------------------------------------------------------------*
 819 !=  PARAMETERS:
 820 !=  NW      - INTEGER, number of specified intervals + 1 in working       (I)=*
 821 !=            wavelength grid                                                =*
 822 !=  WL      - REAL, vector of lower limits of wavelength intervals in     (I)=*
 823 !=            working wavelength grid
 824 !=            vertical layer at each specified wavelength                    =*
 825 !=  O2XS1   - REAL, O2 molecular absorption cross section                    =*
 826 !=
 827 !-----------------------------------------------------------------------------*
 828
 829 ! Input
 830
 831       INTEGER, intent(in) :: nw
 832       REAL, intent(in)    :: wl(nw)
 833
 834 ! Output O2 xsect, temporary, will be over-written in Lyman-alpha and
 835 !   Schumann-Runge wavelength bands.
 836
 837       REAL, intent(inout) :: o2xs1(:)
 838
 839 ! Internal
 840
 841       INTEGER, parameter :: kdata = 200
 842       INTEGER :: i, n
 843       INTEGER :: ierr
 844       REAL    :: x, y
 845       REAL    :: x1(kdata), y1(kdata)
 846       CHARACTER(len=256) :: emsg
 847
 848 ! Read O2 absorption cross section data:
 849 !  116.65 to 203.05 nm = from Brasseur and Solomon 1986
 850 !  205 to 240 nm = Yoshino et al. 1988
 851
 852 ! Note that subroutine la_srb.f will over-write values in the spectral regions
 853 !   corresponding to:
 854 ! - Lyman-alpha (LA: 121.4-121.9 nm, Chabrillat and Kockarts parameterization)
 855 ! - Schumann-Runge bands (SRB: 174.4-205.8 nm, Koppers parameteriaztion)
 856
 857       n = 0
 858
 859       OPEN(UNIT=kin,FILE='DATAE1/O2/O2_brasseur.abs',iostat=ierr)
 860       if( ierr /= 0 ) then
 861         call wrf_error_fatal( 'rdso2xs: Failed to open DATAE1/O2/O2_brasseur.abs' )
 862       endif
 863       DO i = 1, 7
 864          READ(kin,*,iostat=ierr)
 865          if( ierr /= 0 ) then
 866            call wrf_error_fatal( 'rdso2xs: Failed to read DATAE1/O2/O2_brasseur.abs' )
 867          endif
 868       ENDDO
 869       DO i = 1, 78
 870          READ(kin,*,iostat=ierr) x, y
 871          if( ierr /= 0 ) then
 872            call wrf_error_fatal( 'rdso2xs: Failed to read DATAE1/O2/O2_brasseur.abs' )
 873          endif
 874          IF (x .LE. 204.) THEN
 875             n = n + 1
 876             x1(n) = x
 877             y1(n) = y
 878          ENDIF
 879       ENDDO
 880       CLOSE(kin)
 881
 882       OPEN(UNIT=kin,FILE='DATAE1/O2/O2_yoshino.abs',STATUS='old',iostat=ierr)
 883       if( ierr /= 0 ) then
 884          call wrf_error_fatal( 'rdso2xs: Failed to open DATAE1/O2/O2_yoshino.abs' )
 885       endif
 886
 887       DO i = 1, 8
 888          READ(kin,*,iostat=ierr)
 889          if( ierr /= 0 ) then
 890            call wrf_error_fatal( 'rdso2xs: Failed to read DATAE1/O2/O2_yoshino.abs' )
 891          endif
 892       ENDDO
 893       DO i = 1, 36
 894          n = n + 1
 895          READ(kin,*,iostat=ierr) x, y
 896          if( ierr /= 0 ) then
 897            call wrf_error_fatal( 'rdso2xs: Failed to read DATAE1/O2/O2_yoshino.abs' )
 898          endif
 899          y1(n) = y*1.E-24
 900          x1(n) = x
 901       END DO
 902       CLOSE (kin)
 903
 904 ! Add termination points and interpolate onto the
 905 !  user grid (set in subroutine gridw):
 906
 907       CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),y1(1))
 908       CALL addpnt(x1,y1,kdata,n,0.               ,y1(1))
 909       CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
 910       CALL addpnt(x1,y1,kdata,n,              1.E+38,0.)
 911       CALL inter2(nw,wl,o2xs1, n,x1,y1, ierr)
 912       IF (ierr .NE. 0) THEN
 913          WRITE(emsg,'(''rdo2xs: interp err = '',i5,'' in O2 -> O + O'')') ierr
 914          call wrf_error_fatal( trim(emsg) )
 915       ENDIF
 916
 917       END SUBROUTINE rdo2xs
 918
 919 !=============================================================================*
 920
 921       SUBROUTINE rdno2xs(nw,wl)
 922 !-----------------------------------------------------------------------------*
 923 !=  PURPOSE:                                                                 =*
 924 !=  Read NO2 molecular absorption cross section.  Re-grid data to match      =*
 925 !=  specified wavelength working grid.                                       =*
 926 !-----------------------------------------------------------------------------*
 927 !=  PARAMETERS:                                                              =*
 928 !=  NW     - INTEGER, number of specified intervals + 1 in working        (I)=*
 929 !=           wavelength grid                                                 =*
 930 !=  WL     - REAL, vector of lower limits of wavelength intervals in      (I)=*
 931 !=           working wavelength grid                                         =*
 932 !=  NO2XS  - REAL, molecular absoprtion cross section (cm^2) of NO2 at    (O)=*
 933 !=           each specified wavelength                                       =*
 934 !-----------------------------------------------------------------------------*
 935
 936 ! input:
 937
 938       INTEGER, intent(in) :: nw
 939       REAL, intent(in)    :: wl(nw)
 940
 941 ! locals:
 942       INTEGER, parameter :: kdata = 100
 943       INTEGER :: iw
 944       INTEGER :: i, n1, n2, ierr
 945       REAL    :: dum1, dum2
 946       REAL    :: x1(kdata), x2(kdata), y1(kdata), y2(kdata)
 947
 948 ! NO2 absorption cross section from JPL2006
 949 ! with interpolation of bin midpoints
 950
 951       OPEN(UNIT=kin,FILE='DATAE1/NO2/NO2_jpl2006.abs',STATUS='old')
 952       DO i = 1, 3
 953          READ(kin,*)
 954       ENDDO
 955       n1 = 81
 956       DO i = 1, n1
 957          READ(kin,*) dum1, dum2, y1(i), y2(i)
 958          x1(i) = 0.5 * (dum1 + dum2)
 959          x2(i) = x1(i)
 960          y1(i) = y1(i)*1.E-20
 961          y2(i) = y2(i)*1.E-20
 962       ENDDO
 963       CLOSE(kin)
 964       n2 = n1
 965
 966       CALL addpnt(x1,y1,kdata,n1,x1(1)*(1.-deltax),0.)
 967       CALL addpnt(x1,y1,kdata,n1,               0.,0.)
 968       CALL addpnt(x1,y1,kdata,n1,x1(n1)*(1.+deltax),   0.)
 969       CALL addpnt(x1,y1,kdata,n1,            1.e+38,   0.)
 970       CALL inter2(nw,wl,no2xs_a,n1,x1,y1,ierr)
 971
 972       CALL addpnt(x2,y2,kdata,n2,x2(1)*(1.-deltax),0.)
 973       CALL addpnt(x2,y2,kdata,n2,               0.,0.)
 974       CALL addpnt(x2,y2,kdata,n2,x2(n2)*(1.+deltax),   0.)
 975       CALL addpnt(x2,y2,kdata,n2,            1.e+38,   0.)
 976       CALL inter2(nw,wl,no2xs_b,n2,x2,y2,ierr)
 977
 978       END SUBROUTINE rdno2xs
 979
 980 !=============================================================================*
 981
 982       SUBROUTINE no2xs_jpl06a(nz,t,nw,wl, no2xs)
 983
 984 ! interpolate NO2 xs from JPL2006
 985
 986 ! input:
 987
 988       INTEGER, intent(in) :: nz
 989       INTEGER, intent(in) :: nw
 990       REAL, intent(in)    :: t(nz)
 991       REAL, intent(in)    :: wl(nw)
 992
 993 ! output:
 994
 995       REAL, intent(inout) :: no2xs(:,:)
 996
 997 ! local
 998
 999       INTEGER :: iw
1000       REAL    :: tfac(nz)
1001
1002       tfac(1:nz) = (t(1:nz) - 220.)/74.
1003       DO iw = 1, nw-1
1004         no2xs(1:nz,iw) = no2xs_a(iw) + (no2xs_b(iw)-no2xs_a(iw))*tfac(1:nz)
1005       ENDDO
1006
1007       END SUBROUTINE no2xs_jpl06a
1008
1009 !=============================================================================*
1010
1011       SUBROUTINE rdso2xs(nw,wl,so2xs)
1012 !-----------------------------------------------------------------------------*
1013 !=  PURPOSE:                                                                 =*
1014 !=  Read SO2 molecular absorption cross section.  Re-grid data to match      =*
1015 !=  specified wavelength working grid.                                       =*
1016 !-----------------------------------------------------------------------------*
1017
1018       INTEGER, parameter :: kdata = 1000
1019
1020 ! input: (altitude working grid)
1021       INTEGER, intent(in) :: nw
1022       REAL, intent(in)    :: wl(nw)
1023
1024 ! output:
1025
1026       REAL, intent(inout) :: so2xs(nw)
1027
1028 ! local:
1029       REAL x1(kdata)
1030       REAL y1(kdata)
1031       REAL yg(kw)
1032       REAL dum
1033       INTEGER ierr
1034       INTEGER i, l, n, idum
1035       CHARACTER(len=40)  :: fil
1036       CHARACTER(len=256) :: emsg
1037 !************ absorption cross sections:
1038 ! SO2 absorption cross sections from J. Quant. Spectrosc. Radiat. Transfer
1039 ! 37, 165-182, 1987, T. J. McGee and J. Burris Jr.
1040 ! Angstrom vs. cm2/molecule, value at 221 K
1041
1042       fil = 'DATA/McGee87'
1043       OPEN(UNIT=kin,FILE='DATAE1/SO2/SO2xs.all',STATUS='old')
1044       DO i = 1,3
1045         read(kin,*)
1046       ENDDO
1047       n = 704
1048       DO i = 1, n
1049         READ(kin,*) x1(i), y1(i)
1050         x1(i) = .1*x1(i)
1051       ENDDO
1052       CLOSE (kin)
1053
1054       CALL addpnt(x1,y1,kdata,n,x1(1)*(1.-deltax),0.)
1055       CALL addpnt(x1,y1,kdata,n,          0.,0.)
1056       CALL addpnt(x1,y1,kdata,n,x1(n)*(1.+deltax),0.)
1057       CALL addpnt(x1,y1,kdata,n,      1.e+38,0.)
1058       CALL inter2(nw,wl,so2xs,n,x1,y1,ierr)
1059       IF (ierr .NE. 0) THEN
1060          WRITE(emsg,'(''rdso2xs: interp err = '',i5,'' in file '',a)') ierr, fil
1061          call wrf_error_fatal( trim(emsg) )
1062       ENDIF
1063
1064       END SUBROUTINE rdso2xs
1065
1066       real FUNCTION refrac(w,airden)
1067
1068       IMPLICIT NONE
1069
1070 ! input vacuum wavelength, nm and air density, molec cm-3
1071
1072       REAL, intent(in) :: w, airden
1073
1074 ! output refractive index for standard air
1075 ! (dry air at 15 deg. C, 101.325 kPa, 0.03% CO2)
1076
1077 ! internal
1078
1079       REAL :: sig,  sigsq, dum
1080
1081 ! from CRC Handbook, originally from Edlen, B., Metrologia, 2, 71, 1966.
1082 ! valid from 200 nm to 2000 nm
1083 ! beyond this range, use constant value
1084
1085       IF (w < 200.) then
1086         dum = 5.e-3
1087       ELSEIF (w > 2000.) then
1088         dum = 5.e-4
1089       ELSE
1090         dum = 1./w
1091       ENDIF
1092       sig = 1.E3*dum
1093       sigsq = sig * sig
1094
1095       dum = 8342.13 + 2406030./(130. - sigsq) + 15997./(38.9 - sigsq)
1096
1097 ! adjust to local air density
1098       dum = dum * airden/(2.69e19 * 273.15/288.15)
1099
1100 ! index of refraction:
1101       refrac = 1. + 1.E-8 * dum
1102
1103       END function refrac
1104
1105       end module module_xsections