chem/la_srb.F

   1 !=============================================================================
   2 ! This file contains the following subroutines, related to the calculation
   3 ! of radiation at Lyman-alpha and Schumann-Runge wavelengths:
   4 !     la_srb
   5 !     lymana
   6 !     schum
   7 !     effxs
   8 !     calc_params
   9 !     init_xs
  10 !     sjo2
  11 ! and the following functions
  12 !     chebev
  13 !=============================================================================
  14
  15       module SRB
  16
  17       implicit none
  18
  19       private
  20       public :: la_srb, sjo2, init_srb
  21       public :: nchebev_term, nchebev_wave
  22       public :: chebev_ac, chebev_bc
  23       public :: ila, isrb
  24
  25       INTEGER, parameter :: kla = 2
  26       INTEGER, PARAMETER :: ksrb = 18
  27       integer, parameter :: nla =  kla - 1
  28       integer, parameter :: nsrb = ksrb - 1
  29
  30       integer :: nchebev_term, nchebev_wave
  31
  32       integer :: ila, isrb
  33       REAL(8) :: b(3), c(3), d(3), e(3)
  34       REAL(8), allocatable :: chebev_ac(:,:)
  35       REAL(8), allocatable :: chebev_bc(:,:)
  36
  37       REAL    :: xslod(nsrb)
  38       REAL    :: wlsrb(ksrb)
  39       REAL    :: wlla(kla)
  40
  41       CONTAINS
  42
  43       SUBROUTINE init_srb
  44
  45       b(:) = (/ 6.8431e-01_8,  2.29841e-01_8,  8.65412e-02_8 /)
  46       c(:) = (/ 8.22114e-21_8, 1.77556e-20_8,  8.22112e-21_8 /)
  47       d(:) = (/ 6.0073e-21_8,  4.28569e-21_8,  1.28059e-20_8 /)
  48       e(:) = (/ 8.21666e-21_8, 1.63296e-20_8,  4.85121e-17_8 /)
  49       xslod(:) = (/6.2180730E-21, 5.8473627E-22, 5.6996334E-22, &
  50                    4.5627094E-22, 1.7668250E-22, 1.1178808E-22, &
  51                    1.2040544E-22, 4.0994668E-23, 1.8450616E-23, &
  52                    1.5639540E-23, 8.7961075E-24, 7.6475608E-24, &
  53                    7.6260556E-24, 7.5565696E-24, 7.6334338E-24, &
  54                    7.4371992E-24, 7.3642966E-24 /)
  55       wlla(:)  = (/ 121.4, 121.9/)
  56       wlsrb(:) = (/174.4, 177.0, 178.6, 180.2, 181.8, &
  57                    183.5, 185.2, 186.9, 188.7, 190.5, &
  58                    192.3, 194.2, 196.1, 198.0, 200.0, &
  59                    202.0, 204.1, 205.8/)
  60
  61       END SUBROUTINE init_srb
  62
  63       SUBROUTINE la_srb( nlyr, z, tlev, wmin, &
  64                          vcol, scol, o2_xs, dto2, srb_o2_xs )
  65 !-----------------------------------------------------------------------------
  66 !=  PURPOSE:
  67 !=  Compute equivalent optical depths for O2 absorption, and O2 effective
  68 !=  absorption cross sections, parameterized in the Lyman-alpha and SR bands
  69 !-----------------------------------------------------------------------------
  70 !=  PARAMETERS:
  71 !=  NZ      - INTEGER, number of specified altitude levels in the working (I)
  72 !=            grid
  73 !=  Z       - REAL, specified altitude working grid (km)                  (I)
  74 !=  NW      - INTEGER, number of specified intervals + 1 in working       (I)
  75 !=            wavelength grid
  76 !=  WL      - REAL, vector of lxower limits of wavelength intervals in    (I)
  77 !=            working wavelength grid
  78 !=  CZ      - REAL, number of air molecules per cm^2 at each specified    (I)
  79 !=            altitude layer
  80 !=  ZEN     - REAL, solar zenith angle                                    (I)
  81 !=
  82 !=  O2XS1   - REAL, O2 cross section from rdo2xs                          (I)
  83 !=
  84 !=  DTO2    - REAL, optical depth due to O2 absorption at each specified  (O)
  85 !=            vertical layer at each specified wavelength
  86 !=  O2XS    - REAL, molecular absorption cross section in SR bands at     (O)
  87 !=            each specified altitude and wavelength.  Includes Herzberg
  88 !=            continuum.
  89 !-----------------------------------------------------------------------------
  90
  91       use params_mod, only : o2vmr, largest
  92
  93 !-----------------------------------------------------------------------------
  94 !     ... dummy arguments
  95 !-----------------------------------------------------------------------------
  96       INTEGER, intent(in) :: nlyr
  97       REAL, intent(in) :: wmin
  98       REAL, intent(in) :: z(:)
  99       REAL, intent(in) :: tlev(:)
 100
 101       REAL, intent(in) :: vcol(:)
 102       REAL, intent(in) :: scol(:)
 103       REAL, intent(in) :: o2_xs(:)
 104       REAL, intent(inout) :: dto2(:,:)
 105       REAL, intent(inout) :: srb_o2_xs(:,:)
 106
 107 !-----------------------------------------------------------------------------
 108 !     ... local variables
 109 !-----------------------------------------------------------------------------
 110       REAL :: secchi(nlyr)
 111       REAL :: o2col(nlyr)
 112
 113 !-----------------------------------------------------------------------------
 114 ! Lyman-alpha variables
 115 ! O2 optical depth and equivalent cross section in the Lyman-alpha region
 116 !-----------------------------------------------------------------------------
 117       INTEGER :: nlev
 118       INTEGER :: nlev_srb
 119       INTEGER :: k, iw, wn
 120       REAL    :: dto2la(nlyr,nla), o2xsla(nlyr,nla)
 121
 122 !-----------------------------------------------------------------------------
 123 ! grid on which Koppers' parameterization is defined
 124 ! O2 optical depth and equivalent cross section on Koppers' grid
 125 !-----------------------------------------------------------------------------
 126       REAL    :: dto2k(nlyr,nsrb), o2xsk(nlyr,nsrb)
 127
 128       nlev_srb = size( srb_o2_xs,dim=2 )
 129       nlev = nlyr
 130 !----------------------------------------------------------------------
 131 ! initalize O2 cross sections
 132 !----------------------------------------------------------------------
 133       DO k = 1, nlev_srb
 134         srb_o2_xs(:,k) = o2_xs(:)
 135       END DO
 136
 137       IF( wmin <= wlsrb(nsrb) ) THEN
 138 !----------------------------------------------------------------------
 139 ! Slant O2 column and x-sections.
 140 !----------------------------------------------------------------------
 141         o2col(:nlyr) = o2vmr * scol(:nlyr)
 142 !----------------------------------------------------------------------
 143 ! Effective secant of solar zenith angle.
 144 ! Use 2.0 if no direct sun (value for isotropic radiation)
 145 ! For nz, use value at nz-1
 146 !----------------------------------------------------------------------
 147         WHERE( scol(:nlyr) > .1*largest )
 148           secchi(:nlyr) = 2.
 149         ELSEWHERE
 150           secchi(:nlyr) = scol(:nlyr)/vcol(:nlyr)
 151         ENDWHERE
 152
 153 !---------------------------------------------------------------------
 154 ! Lyman-Alpha parameterization, output values of O2 optical depth
 155 ! and O2 effective (equivalent) cross section
 156 !----------------------------------------------------------------------
 157         CALL lymana( nlyr, o2col, secchi, dto2la, o2xsla )
 158         DO wn = ila, ila + nla - 1
 159           iw = wn - ila + 1
 160           dto2(:nlyr,wn)          = dto2la(:nlyr,iw)
 161           srb_o2_xs(wn,:nlev_srb) = o2xsla(2:nlev_srb+1,iw)
 162         ENDDO
 163
 164 !------------------------------------------------------------------------------
 165 ! Koppers' parameterization of the SR bands, output values of O2
 166 ! optical depth and O2 equivalent cross section
 167 !------------------------------------------------------------------------------
 168         CALL schum( nlyr, o2col, tlev, secchi, dto2k, o2xsk )
 169         DO wn = isrb, isrb + nsrb - 1
 170           iw = wn - isrb + 1
 171           dto2(:nlyr,wn)          = dto2k(:nlyr,iw)
 172           srb_o2_xs(wn,:nlev_srb) = o2xsk(2:nlev_srb+1,iw)
 173         ENDDO
 174       ENDIF
 175
 176       END SUBROUTINE la_srb
 177
 178       SUBROUTINE lymana( nlyr, o2col, secchi, dto2la, o2xsla )
 179 !-----------------------------------------------------------------------------
 180 !=  PURPOSE:
 181 !=  Calculate the effective absorption cross section of O2 in the Lyman-Alpha
 182 !=  bands and an effective O2 optical depth at all altitudes.  Parameterized
 183 !=  after:  Chabrillat, S., and G. Kockarts, Simple parameterization of the
 184 !=  absorption of the solar Lyman-Alpha line, Geophysical Research Letters,
 185 !=  Vol.24, No.21, pp 2659-2662, 1997.
 186 !-----------------------------------------------------------------------------
 187 !=  PARAMETERS:
 188 !=  NZ      - INTEGER, number of specified altitude levels in the working (I)
 189 !=            grid
 190 !=  O2COL   - REAL, slant overhead O2 column (molec/cc) at each specified (I)
 191 !=            altitude
 192 !=  DTO2LA  - REAL, optical depth due to O2 absorption at each specified  (O)
 193 !=            vertical layer
 194 !=  O2XSLA  - REAL, molecular absorption cross section in LA bands        (O)
 195 !-----------------------------------------------------------------------------
 196
 197
 198 !-----------------------------------------------------------------------------
 199 !     ... dummy arguments
 200 !-----------------------------------------------------------------------------
 201       INTEGER, intent(in) :: nlyr
 202       REAL,    intent(in) :: o2col(:)
 203       REAL,    intent(in) :: secchi(:)
 204       REAL, intent(inout) :: dto2la(nlyr,nla), o2xsla(nlyr,nla)
 205
 206 !-----------------------------------------------------------------------------
 207 !     ... local variables
 208 !-----------------------------------------------------------------------------
 209       REAL, parameter    :: xsmin = 1.e-20
 210       REAL(8), parameter :: rmmin = 1.e-100_8
 211
 212       INTEGER :: k, kp1, wn
 213       REAL(8) :: o2_col
 214       REAL(8) :: rm(nlyr), ro2(nlyr)
 215       REAL(8) :: rm_wrk(3), ro2_wrk(3)
 216
 217       do wn = 1,nla
 218         dto2la(:nlyr,wn) = 0.
 219         o2xsla(:nlyr,wn) = 0.
 220       end do
 221 !-----------------------------------------------------------------------------
 222 ! calculate reduction factors at every layer
 223 !-----------------------------------------------------------------------------
 224       rm(:nlyr)  = 0._8
 225       ro2(:nlyr) = 0._8
 226       DO k = 1, nlyr
 227         o2_col = real( o2col(k),8 )
 228         rm_wrk(:)  = b(:) * EXP( -c(:) * o2_col )
 229         ro2_wrk(:) = d(:) * EXP( -e(:) * o2_col )
 230         rm(k)  = sum( rm_wrk )
 231         ro2(k) = sum( ro2_wrk )
 232       ENDDO
 233
 234 !-----------------------------------------------------------------------------
 235 ! calculate effective O2 optical depths and effective O2 cross sections
 236 !-----------------------------------------------------------------------------
 237       DO k = 1, nlyr-1
 238         kp1 = k + 1
 239         IF (rm(k) > rmmin) THEN
 240           IF (ro2(k) > rmmin) THEN
 241             o2xsla(k,1) = REAL( ro2(k)/rm(k) )
 242           ELSE
 243             o2xsla(k,1) = xsmin
 244           ENDIF
 245
 246           IF (rm(kp1) > 0._8) THEN
 247             dto2la(k,1) = LOG( rm(kp1) )/secchi(kp1)  &
 248                         - LOG( rm(k))   /secchi(k)
 249           ELSE
 250             dto2la(k,1) = 1000.
 251           ENDIF
 252         ELSE
 253           dto2la(k,1) = 1000.
 254           o2xsla(k,1) = xsmin
 255         ENDIF
 256       END DO
 257
 258 !-----------------------------------------------------------------------------
 259 ! do top layer separately
 260 !-----------------------------------------------------------------------------
 261       IF( rm(nlyr) > rmmin ) THEN
 262         o2xsla(nlyr,1) = REAL( ro2(nlyr)/rm(nlyr) )
 263       ELSE
 264         o2xsla(nlyr,1) = xsmin
 265       ENDIF
 266
 267       END SUBROUTINE lymana
 268
 269       SUBROUTINE schum( nlyr, o2col, tlev, secchi, dto2, o2xsk )
 270 !-----------------------------------------------------------------------------
 271 !=  PURPOSE:
 272 !=  Calculate the equivalent absorption cross section of O2 in the SR bands.
 273 !=  The algorithm is based on parameterization of G.A. Koppers, and
 274 !=  D.P. Murtagh [ref. Ann.Geophys., 14 68-79, 1996]
 275 !=  Final values do include effects from the Herzberg continuum.
 276 !-----------------------------------------------------------------------------
 277 !=  PARAMETERS:
 278 !=  NZ      - INTEGER, number of specified altitude levels in the working (I)
 279 !=            grid
 280 !=  O2COL   - REAL, slant overhead O2 column (molec/cc) at each specified (I)
 281 !=            altitude
 282 !=  TLEV    - tmeperature at each level
 283 !=  SECCHI  - ratio of slant to vertical o2 columns
 284 !=  DTO2    - REAL, optical depth due to O2 absorption at each specified
 285 !=            vertical layer at each specified wavelength
 286 !=  O2XSK  - REAL, molecular absorption cross section in SR bands at
 287 !=            each specified wavelength.  Includes Herzberg continuum
 288 !-----------------------------------------------------------------------------
 289
 290       use params_mod, only : precis
 291
 292 !-----------------------------------------------------------------------------
 293 !     ... dummy arguments
 294 !-----------------------------------------------------------------------------
 295       INTEGER, intent(in) :: nlyr
 296       REAL,    intent(in) :: o2col(:)
 297       REAL,    intent(in) :: tlev(:), secchi(:)
 298       REAL, intent(inout) :: dto2(nlyr,nsrb), o2xsk(nlyr,nsrb)
 299
 300 !-----------------------------------------------------------------------------
 301 !     ... local variables
 302 !-----------------------------------------------------------------------------
 303       REAL, parameter :: o2col_min = exp( 38. )
 304
 305       INTEGER :: wn, k, ktop, ktop1, kbot, nlyrm1
 306       REAL    :: x
 307       REAL    :: o2col1(nlyr)
 308       REAL    :: xs(nsrb)
 309
 310       nlyrm1 = nlyr - 1
 311 !-----------------------------------------------------------------------------
 312 !     ...Initialize cross sections to values at large optical depth
 313 !-----------------------------------------------------------------------------
 314       DO wn = 1, nsrb
 315         o2xsk(:nlyr,wn) = xslod(wn)
 316       END DO
 317
 318 !-----------------------------------------------------------------------------
 319 !     Calculate cross sections
 320 !     Set smallest O2col = exp(38.) molec cm-2
 321 !     to stay in range of parameterization
 322 !     given by Koppers et al. at top of atm.
 323 !-----------------------------------------------------------------------------
 324       ktop = 2*nlyr
 325       kbot = 0
 326
 327       DO k = 1,nlyr
 328         o2col1(k) = MAX( o2col(k),o2col_min )
 329         x  = LOG( o2col1(k) )
 330         IF (x < 38.0) THEN
 331           ktop1 = k-1
 332           ktop  = MIN(ktop1,ktop)
 333         ELSE IF (x > 56.0) THEN
 334           kbot = k
 335         ELSE
 336           CALL effxs( x, tlev(k), xs )
 337           o2xsk(k,:nsrb) = xs(:nsrb)
 338         ENDIF
 339       END DO
 340
 341 !-----------------------------------------------------------------------------
 342 !  fill in cross section where X is out of range by repeating edge table values
 343 !  Do not allow kbot = nlyr to avoid division by zero in no light case.
 344 !-----------------------------------------------------------------------------
 345       IF( kbot == nlyr) then
 346         kbot = nlyrm1
 347       ENDIF
 348
 349       IF( kbot > 0 ) THEN
 350         DO wn = 1,nsrb
 351           o2xsk(:kbot,wn) = o2xsk(kbot+1,wn)
 352         END DO
 353       ENDIF
 354
 355       IF( ktop < nlyr ) THEN
 356         DO wn = 1,nsrb
 357           o2xsk(ktop+1:nlyr,wn) = o2xsk(ktop,wn)
 358         END DO
 359       ENDIF
 360
 361 !-----------------------------------------------------------------------------
 362 !  Calculate incremental optical depths
 363 !-----------------------------------------------------------------------------
 364       dto2(nlyr,1:nsrb) = 0.0       ! set optical depth to zero at top
 365       DO wn = 1,nsrb
 366 !-----------------------------------------------------------------------------
 367 !     ... calculate an optical depth weighted by density,
 368 !         put in mean value estimate, if in shade
 369 !-----------------------------------------------------------------------------
 370         WHERE (ABS(1. - o2col1(2:nlyr)/o2col1(:nlyrm1)) <= 2.*precis)
 371           dto2(:nlyrm1,wn) = o2xsk(2:nlyr,wn)*o2col1(2:nlyr)/real(nlyrm1)
 372         ELSEWHERE
 373           dto2(:nlyr-1,wn) = ABS( &
 374             (o2xsk(2:nlyr,wn)*o2col1(2:nlyr) - o2xsk(:nlyrm1,wn)*o2col1(:nlyrm1)) &
 375             /(1. + LOG(o2xsk(2:nlyr,wn)/o2xsk(:nlyrm1,wn))  &
 376               / LOG(o2col1(2:nlyr)/o2col1(:nlyrm1))) )
 377 !-----------------------------------------------------------------------------
 378 !     ... change to vertical optical depth
 379 !-----------------------------------------------------------------------------
 380           dto2(:nlyrm1,wn) = 2. * dto2(:nlyrm1,wn)/(secchi(:nlyr-1)+secchi(2:nlyr))
 381         ENDWHERE
 382       END DO
 383
 384       END SUBROUTINE schum
 385
 386       SUBROUTINE EFFXS( x, t, xs )
 387 !-----------------------------------------------------------------------------
 388 !     Subroutine for evaluating the effective cross section
 389 !     of O2 in the Schumann-Runge bands using parameterization
 390 !     of G.A. Koppers, and D.P. Murtagh [ref. Ann.Geophys., 14
 391 !     68-79, 1996]
 392 !
 393 !     method:
 394 !     ln(xs) = A(X)[T-220]+B(X)
 395 !     X = log of slant column of O2
 396 !     A,B calculated from Chebyshev polynomial coeffs
 397 !     AC and BC using NR routine chebev.  Assume interval
 398 !     is 38<ln(NO2)<56.
 399 !-----------------------------------------------------------------------------
 400
 401 !-----------------------------------------------------------------------------
 402 !     ... dummy arguments
 403 !-----------------------------------------------------------------------------
 404       REAL, intent(in)  :: t, x
 405       REAL, intent(out) :: xs(nsrb)
 406
 407 !-----------------------------------------------------------------------------
 408 !     ... local variables
 409 !-----------------------------------------------------------------------------
 410       INTEGER :: i
 411       REAL    :: a(nsrb), b(nsrb)
 412
 413       call calc_params( x, a, b )
 414
 415       xs(:nsrb) = EXP( a(:nsrb)*( t - 220.) + b(:nsrb) )
 416
 417       END SUBROUTINE EFFXS
 418
 419       SUBROUTINE CALC_PARAMS( x, a, b )
 420 !-----------------------------------------------------------------------------
 421 !     calculates coefficients (A,B), used in calculating the
 422 !     effective cross section, for nsrb wavelength intervals
 423 !     as a function of log O2 column density (X)
 424 !     Wavelength intervals are defined in WMO1985
 425 !-----------------------------------------------------------------------------
 426
 427 !-----------------------------------------------------------------------------
 428 !     ... dummy arguments
 429 !-----------------------------------------------------------------------------
 430       REAL, intent(in)  :: x
 431       REAL, intent(out) :: a(nsrb), b(nsrb)
 432
 433 !-----------------------------------------------------------------------------
 434 !     ... local variables
 435 !-----------------------------------------------------------------------------
 436       INTEGER :: wn
 437
 438 !-----------------------------------------------------------------------------
 439 !     call Chebyshev Evaluation routine to calc A and B from
 440 !     set of 20 coeficients for each wavelength
 441 !-----------------------------------------------------------------------------
 442
 443       DO wn = 1,nsrb
 444         a(wn) = chebev( 38.0 , 56.0, chebev_ac(:,wn), nchebev_term, x )
 445         b(wn) = chebev( 38.0 , 56.0, chebev_bc(:,wn), nchebev_term, x )
 446       END DO
 447
 448       END SUBROUTINE CALC_PARAMS
 449
 450       REAL FUNCTION chebev( a, b, c, m, x )
 451 !-------------------------------------------------------------
 452 !     Chebyshev evaluation algorithm
 453 !     See Numerical recipes p193
 454 !-------------------------------------------------------------
 455
 456 !-------------------------------------------------------------
 457 !       ... dummy arguments
 458 !-------------------------------------------------------------
 459       INTEGER, intent(in) :: m
 460       REAL,    intent(in) :: a, b, x
 461       REAL(8), intent(in) :: c(:)
 462
 463 !-------------------------------------------------------------
 464 !       ... local variables
 465 !-------------------------------------------------------------
 466       INTEGER :: j
 467       REAL    :: d, dd, sv, y, y2
 468
 469       IF( (x - a)*(x - b) > 0.) THEN
 470         chebev = 0.0
 471       ELSE
 472         d  = 0.
 473         dd = 0.
 474         y  = (2.*x - a - b)/(b - a)
 475         y2 = 2.*y
 476         DO J = m,2,-1
 477           sv = d
 478           d  = y2*d - dd + real( c(J),4 )
 479           dd = sv
 480         END DO
 481         chebev = y*d - dd + 0.5*real( c(1),4 )
 482       ENDIF
 483
 484       END FUNCTION chebev
 485
 486       SUBROUTINE sjo2( nlyr, nwave, xso2, xsqy )
 487 !-----------------------------------------------------------------------------
 488 !=  PURPOSE:
 489 !=  Update the weighting function (cross section x quantum yield) for O2
 490 !=  photolysis.  The strong spectral variations in the O2 cross sections are
 491 !=  parameterized into a few bands for Lyman-alpha (121.4-121.9 nm, one band)
 492 !=  and Schumann-Runge (174.4-205.8, nsrb bands) regions. The parameterizations
 493 !=  depend on the overhead O2 column, and therefore on altitude and solar
 494 !=  zenith angle, so they need to be updated at each time/zenith step.
 495 !-----------------------------------------------------------------------------
 496 !=  PARAMETERS:
 497 !=  NZ     - INTEGER, number of altitude levels in working altitude grid  (I)
 498 !=  NW     - INTEGER, number of specified intervals + 1 in working        (I)
 499 !=           wavelength grid
 500 !=  XSO2   - REAL, molecular absorption cross section in SR bands at      (I)
 501 !=           each specified altitude and wavelength.  Includes Herzberg
 502 !=            continuum.
 503 !=  NJ     - INTEGER, index of O2 photolysis in array SQ                  (I)
 504 !=  xsqy   - REAL, cross section x quantum yield (cm^2) for each          (O)
 505 !=           photolysis reaction, at each wavelength and each altitude level
 506 !-----------------------------------------------------------------------------
 507
 508
 509 !-----------------------------------------------------------------------------
 510 !     ... dummy arguments
 511 !-----------------------------------------------------------------------------
 512       INTEGER, intent(in)    :: nlyr, nwave
 513       REAL,    intent(in)    :: xso2(:,:)
 514       REAL,    intent(inout) :: xsqy(:,:)
 515
 516 !-----------------------------------------------------------------------------
 517 !     ... local variables
 518 !-----------------------------------------------------------------------------
 519       INTEGER :: k
 520
 521 !-----------------------------------------------------------------------------
 522 ! O2 + hv -> O + O
 523 ! quantum yield assumed to be unity
 524 ! assign cross section values at all wavelengths and at all altitudes
 525 !      qy = 1.
 526 !-----------------------------------------------------------------------------
 527       DO k = 1, nlyr
 528         xsqy(:nwave,k) = xso2(:nwave,k)
 529       END DO
 530
 531       END SUBROUTINE sjo2
 532
 533       end module SRB