LCOV - code coverage report
Current view: | top level - chemistry/mozart - photo_bkgrnd.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 0 | 48 | 0.0 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 0 | 3 | 0.0 %

              Line data    Source code
           1             : !----------------------------------------------------------------------------
           2             : ! For calculating background ionization due to star light, geo-corona radiation
           3             : ! Applicable to high altitudes of WACCM and WACCMX
           4             : ! Module created by Francis Vitt 14 Feb 2013
           5             : ! Background ionization algorithm provided by Stan Solomn
           6             : !----------------------------------------------------------------------------
           7             : module photo_bkgrnd
           8             :
           9             :   use shr_kind_mod,   only : r8 => shr_kind_r8
          10             :   use mo_chem_utls,   only : get_rxt_ndx
          11             :   use ppgrid,         only : pver
          12             :
          13             :   implicit none
          14             :
          15             :   private
          16             :   public :: photo_bkgrnd_calc
          17             :   public :: photo_bkgrnd_init
          18             :
          19             :   integer :: jo_ndx, jo2_ndx, jn2_ndx, jn_ndx, jno_ndx
          20             :
          21             : contains
          22             :
          23             :   !----------------------------------------------------------------------------
          24             :   ! look up corresponding reaction rate indices
          25             :   !----------------------------------------------------------------------------
          26           0 :   subroutine photo_bkgrnd_init()
          27           0 :     jo_ndx  = get_rxt_ndx( 'jeuv_1' ) ! O + hv  -> Op + e
          28           0 :     jo2_ndx = get_rxt_ndx( 'jeuv_5' ) ! O2 + hv -> O2p + e
          29           0 :     jn2_ndx = get_rxt_ndx( 'jeuv_6' ) ! N2 + hv -> N2p + e
          30           0 :     jn_ndx  = get_rxt_ndx( 'jeuv_10') ! N2 + hv -> N + Np + e
          31           0 :     jno_ndx = get_rxt_ndx( 'jno_i' )  ! NO + hv -> NOp + e
          32           0 :   end subroutine photo_bkgrnd_init
          33             :
          34             :   !----------------------------------------------------------------------------
          35             :   ! Adds background ionization rates to WACCM's photolysis rates
          36             :   !----------------------------------------------------------------------------
          37           0 :   subroutine photo_bkgrnd_calc(f107, o_den, o2_den, n2_den, no_den, zint, rates, &
          38           0 :                                qbko1_out, qbko2_out, qbkn2_out, qbkn1_out, qbkno_out )
          39             :
          40             :    ! arguments
          41             :     real(r8), intent(in) :: f107
          42             :     real(r8), intent(in) :: o_den(:)           ! N density (molecules/cm^3)
          43             :     real(r8), intent(in) :: o2_den(:)          ! O2 density (molecules/cm^3)
          44             :     real(r8), intent(in) :: n2_den(:)          ! N2 density (molecules/cm^3)
          45             :     real(r8), intent(in) :: no_den(:)          ! NO density (molecules/cm^3)
          46             :     real(r8), intent(in) :: zint(:)            ! interface height (km)
          47             :
          48             :     real(r8), intent(inout) :: rates(:,:)      ! photolysis rates (sec-1)
          49             :
          50             :     real(r8), intent(out), optional :: qbko1_out(:) ! rate (cm-3 sec-1) of O  + hv ->  Op  + e
          51             :     real(r8), intent(out), optional :: qbko2_out(:) ! rate (cm-3 sec-1) of O2 + hv ->  O2p + e
          52             :     real(r8), intent(out), optional :: qbkn2_out(:) ! rate (cm-3 sec-1) of N2 + hv ->  N2p + e
          53             :     real(r8), intent(out), optional :: qbkn1_out(:) ! rate (cm-3 sec-1) of N  + hv ->  Np  + e
          54             :     real(r8), intent(out), optional :: qbkno_out(:) ! rate (cm-3 sec-1) of N2 + hv ->  NOp + e
          55             :
          56             :     ! local vars
          57             :     real(r8), parameter :: km2cm = 1.e5_r8
          58             :     integer, parameter :: nmaj = 3
          59             :
          60             :     real(r8) :: zmaj(nmaj,pver)
          61             :     real(r8) :: zno(pver)
          62             :     real(r8) :: zvcd(nmaj,pver)
          63             :     real(r8) :: delz(pver)
          64             :     real(r8) :: qbko1(pver)
          65             :     real(r8) :: qbko2(pver)
          66             :     real(r8) :: qbkn2(pver)
          67             :     real(r8) :: qbkn1(pver)
          68             :     real(r8) :: qbkno(pver)
          69             :
          70             :     integer :: k
          71             :
          72           0 :     zmaj(1,:) = o_den(:)
          73           0 :     zmaj(2,:) = o2_den(:)
          74           0 :     zmaj(3,:) = n2_den(:)
          75           0 :     zno(:)    = no_den(:)
          76             :
          77             :     ! thickness of each layer (cm)
          78           0 :     delz(1:pver-1) = km2cm*(zint(1:pver-1) - zint(2:pver))
          79           0 :     delz(pver) = delz(pver-1)
          80             :
          81           0 :     zvcd(:,:) = 0._r8
          82             :
          83             :     ! intergate column above each layer
          84           0 :     do k = 2,pver
          85           0 :        zvcd(1,k) = zvcd(1,k-1) + delz(k) * o_den(k)
          86           0 :        zvcd(2,k) = zvcd(2,k-1) + delz(k) * o2_den(k)
          87           0 :        zvcd(3,k) = zvcd(3,k-1) + delz(k) * n2_den(k)
          88             :     enddo
          89             :
          90             :     ! invoke Stan's background ionization method -- returns rates (cm-3 sec-1)
          91           0 :     call qback (zmaj,zno,zvcd,f107,nmaj,pver,qbko1,qbko2,qbkn2,qbkn1,qbkno)
          92             :
          93             :     ! divide by densities to get photolysis rates (sec-1)
          94           0 :     if (jo_ndx>0)  rates(:,jo_ndx)  = rates(:,jo_ndx)  + qbko1(:)/o_den(:)
          95           0 :     if (jo2_ndx>0) rates(:,jo2_ndx) = rates(:,jo2_ndx) + qbko2(:)/o2_den(:)
          96           0 :     if (jn2_ndx>0) rates(:,jn2_ndx) = rates(:,jn2_ndx) + qbkn2(:)/n2_den(:)
          97           0 :     if (jn_ndx >0) rates(:,jn_ndx)  = rates(:,jn_ndx)  + qbkn1(:)/n2_den(:)
          98           0 :     if (jno_ndx>0) rates(:,jno_ndx) = rates(:,jno_ndx) + qbkno(:)/no_den(:)
          99             :
         100             :     ! for diagnostics
         101           0 :     if (present(qbko1_out)) qbko1_out(:) = qbko1(:)
         102           0 :     if (present(qbko2_out)) qbko2_out(:) = qbko2(:)
         103           0 :     if (present(qbkn2_out)) qbkn2_out(:) = qbkn2(:)
         104           0 :     if (present(qbkn1_out)) qbkn1_out(:) = qbkn1(:)
         105           0 :     if (present(qbkno_out)) qbkno_out(:) = qbkno(:)
         106             :
         107           0 :   endsubroutine photo_bkgrnd_calc
         108             :
         109             : !----------------------------------------------------------------------------
         110             : ! Private Method
         111             : !----------------------------------------------------------------------------
         112             : !
         113             : ! Stan Solomon, 11/88, 11/92
         114             : ! Comment updated 3/05
         115             : ! New version uses updated TIE-GCM and TIME-GCM qinite.F formulation, 1/13
         116             : !
         117             : ! Estimates background ("nighttime") ionization rates.
         118             : ! Four components are used:
         119             : !    Geocoronal Lyman-beta 102.6 nm (ionizes O2 only)
         120             : !    Geocoronal He I 58.4 nm
         121             : !    Geocoronal He II 30.4 nm
         122             : !    Geocoronal Lyman-alpha 121.6 nm (ionizes NO only)
         123             : !
         124             : ! Definitions:
         125             : !
         126             : ! zmaj    major species O, O2, N2 at each altitude
         127             : ! zno     nitric oxide at each altitude
         128             : ! zvcd    vertical column density for each major species above each altitude
         129             : ! photoi  photoionization rates for each state, species, altitude
         130             : ! f107    solar 10.7 cm radio flux activity index
         131             : ! jm      number of altitude levels
         132             : ! nmaj    number of major species (3)
         133             : ! nst     number of states
         134             : !
         135             : ! al      photon flux at 102.6 nm, 58.4 nm, 30.4 nm
         136             : ! flyan   photon flux at 121.6 nm
         137             : ! sa      absorption cross sections for O, O2, N2 at each wavelength
         138             : ! si      ionization cross sections for O, O2, N2 at each wavelength
         139             : ! flyan   photon flux at 121.6 nm
         140             : ! salyao2 absorption cross section for O2 at 121.6 nm
         141             : ! silyano ionization cross section for NO at 121.6 nm
         142             : ! bn2p    branching ratio for N2+ from ionization of N2
         143             : ! bn1p    branching ratio for N+ from ionization of N2
         144             : ! tau     optical depth
         145             : ! qbko1    production rate of O+
         146             : ! qbko2    production rate of O2+
         147             : ! qbkn2    production rate of N2+
         148             : ! qbkn1    production rate of N+
         149             : ! qbkno    production rate of NO+
         150             : !
         151             : ! All units cgs.
         152             : !
         153             : !
         154           0 : subroutine qback (zmaj,zno,zvcd,f107,nmaj,jm,qbko1,qbko2,qbkn2,qbkn1,qbkno)
         155             :
         156             :  ! args:
         157             :   integer,  intent(in)  :: nmaj,jm
         158             :   real(r8), intent(in)  :: f107
         159             :   real(r8), intent(in)  :: zmaj(nmaj,jm), zno(jm), zvcd(nmaj,jm)
         160             :   real(r8), intent(out) :: qbko1(jm),qbko2(jm),qbkn2(jm),qbkn1(jm),qbkno(jm)
         161             :
         162             :  ! local vars:
         163             :   real(r8) :: al(3), sa(3,3), si(3,3)
         164             :   real(r8) :: salyao2, silyano, bn2p, bn1p
         165             :   real(r8) :: flyan
         166             :   real(r8) :: tau
         167             :   integer :: j,l
         168             :
         169             :   data al /1.5e7_r8, 1.5e6_r8, 1.5e6_r8/
         170             :
         171             :   data sa /      0._r8,  1.6e-18_r8,       0._r8, &
         172             :            10.2e-18_r8, 22.0e-18_r8, 23.1e-18_r8, &
         173             :             8.4e-18_r8, 16.0e-18_r8, 11.6e-18_r8/
         174             :
         175             :   data si /      0._r8,  1.0e-18_r8,       0._r8, &
         176             :            10.2e-18_r8, 22.0e-18_r8, 23.1e-18_r8, &
         177             :             8.4e-18_r8, 16.0e-18_r8, 11.6e-18_r8/
         178             :
         179             :   data salyao2/8.0e-21_r8/
         180             :   data silyano/2.0e-18_r8/
         181             :   data bn2p/0.86_r8/
         182             :   data bn1p/0.14_r8/
         183             :
         184             : !
         185             : ! Calculate Lyman-alpha 121.6 nm geocoronal flux as a function of F10.7:
         186             : !
         187           0 :   flyan = 5.E9_r8*(1._r8+0.002_r8*(f107-65._r8))
         188             : !
         189             : ! Loop over altitudes:
         190             : !
         191           0 :   do j=1,jm
         192             : !
         193             : ! Calculate optical depth and ionization rates for major species:
         194             : !
         195           0 :      qbko1(j)=0._r8
         196           0 :      qbko2(j)=0._r8
         197           0 :      qbkn2(j)=0._r8
         198           0 :      qbkn1(j)=0._r8
         199           0 :      do l=1,3
         200           0 :         tau=(sa(1,l)*zvcd(1,j)+sa(2,l)*zvcd(2,j)+sa(3,l)*zvcd(3,j))
         201           0 :         qbko1(j) = qbko1(j) +       al(l)*si(1,l)*zmaj(1,j)*exp(-tau)
         202           0 :         qbko2(j) = qbko2(j) +       al(l)*si(2,l)*zmaj(2,j)*exp(-tau)
         203           0 :         qbkn2(j) = qbkn2(j) + bn2p*(al(l)*si(3,l)*zmaj(3,j)*exp(-tau))
         204           0 :         qbkn1(j) = qbkn1(j) + bn1p*(al(l)*si(3,l)*zmaj(3,j)*exp(-tau))
         205             :      enddo
         206             : !
         207             : ! Calculate optical depth of Ly-alpha, and ionization rate of NO:
         208             : !
         209           0 :      tau = salyao2*zvcd(2,j)
         210           0 :      qbkno(j) = flyan*silyano*zno(j)*exp(-tau)
         211             :
         212             :   enddo
         213             :
         214           0 :   return
         215             : end subroutine qback
         216             :
         217             : end module photo_bkgrnd

Generated by: LCOV version 1.14