LCOV - code coverage report
Current view: | top level - chemistry/mozart - mee_ap_util_mod.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 0 | 104 | 0.0 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 0 | 9 | 0.0 %

              Line data    Source code
           1             : !------------------------------------------------------------------------------
           2             : ! Utility routines for medium energy electron ionization base on Ap
           3             : ! geomagnetic activity index or observed electron fluxes
           4             : !------------------------------------------------------------------------------
           5             : module mee_ap_util_mod
           6             :   use shr_kind_mod, only: r8 => shr_kind_r8
           7             :   use shr_const_mod, only: pi => shr_const_pi
           8             :   use mee_fluxes, only: mee_fluxes_nenergy, mee_fluxes_energy, mee_fluxes_denergy
           9             :   use mee_fluxes, only: mee_fluxes_active, mee_fluxes_extract
          10             :   use cam_abortutils, only : endrun
          11             :
          12             :   implicit none
          13             :
          14             :   private
          15             :   public :: mee_ap_iprs
          16             :   public :: mee_ap_init
          17             :   public :: mee_ap_final
          18             :   public :: mee_ap_error
          19             :   public :: mee_ap_noerror
          20             :
          21             :   integer, parameter :: mee_ap_error=-1
          22             :   integer, parameter :: mee_ap_noerror=0
          23             :
          24             :   integer :: nbins=0
          25             :
          26             :   real(r8),pointer :: energies(:) => null()  ! energy bins
          27             :   real(r8),pointer :: denergies(:) => null() ! width of each energy bin
          28             :   real(r8) :: solid_angle_factor = -huge(1._r8)
          29             :   real(r8),allocatable :: fang_coefs(:,:)
          30             :
          31             : contains
          32             :
          33             :   !-----------------------------------------------------------------------------
          34             :   ! Sets up energy spectrum grid for medium energy electrons in earth's
          35             :   ! radiation belt
          36             :   !-----------------------------------------------------------------------------
          37           0 :   subroutine mee_ap_init(loss_cone_angle, status)
          38             :
          39             :     real(r8), intent(in) :: loss_cone_angle ! Bounce Loss Cone angle (degrees)
          40             :     integer, intent(out) :: status          ! error status
          41             :
          42             :     integer :: ierr
          43             :     character(len=*), parameter :: subname = 'mee_ap_init: '
          44             :
          45           0 :     status = mee_ap_noerror
          46           0 :     if ( loss_cone_angle<0._r8 .or. loss_cone_angle>90._r8 ) then
          47           0 :        status = mee_ap_error
          48           0 :        return
          49             :     endif
          50             :
          51             :     ! The area of the BLC in sr is 2pi(1-cos(BLC))
          52           0 :     solid_angle_factor = 2._r8*pi*(1._r8-cos(loss_cone_angle*pi/180._r8))
          53             :
          54           0 :     if (mee_fluxes_active) then
          55           0 :        nbins=mee_fluxes_nenergy
          56           0 :        energies=>mee_fluxes_energy
          57           0 :        denergies=>mee_fluxes_denergy
          58             :     else
          59           0 :        nbins=100
          60           0 :        allocate(energies(nbins), stat=ierr)
          61           0 :        if (ierr/=0) call endrun(subname//'not able to allocate energies')
          62           0 :        allocate(denergies(nbins), stat=ierr)
          63           0 :        if (ierr/=0) call endrun(subname//'not able to allocate denergies')
          64           0 :        call gen_energy_grid(nbins, energies, denergies)
          65             :     endif
          66             :
          67           0 :     call init_fang_coefs()
          68             :
          69             :   end subroutine mee_ap_init
          70             :
          71             :   !-----------------------------------------------------------------------------
          72             :   !-----------------------------------------------------------------------------
          73           0 :   subroutine mee_ap_final()
          74             :
          75           0 :     deallocate(fang_coefs)
          76             :
          77           0 :     if (.not.mee_fluxes_active) then
          78           0 :        deallocate(energies)
          79           0 :        deallocate(denergies)
          80             :     endif
          81           0 :     nullify(energies)
          82           0 :     nullify(denergies)
          83             :
          84           0 :  end subroutine mee_ap_final
          85             :
          86             :   !-----------------------------------------------------------------------------
          87             :   ! Computes ion pair production rates base on Ap geomagnetic activity index
          88             :   ! or observed electron fluxes
          89             :   !-----------------------------------------------------------------------------
          90           0 :   subroutine mee_ap_iprs( ncols, nlyrs, airden, scaleh, Ap, ionpairs, status, maglat, lshell)
          91             :
          92             :     integer ,intent(in) :: ncols, nlyrs
          93             :     real(r8),intent(in) :: airden(ncols,nlyrs)     ! g/cm3
          94             :     real(r8),intent(in) :: scaleh(ncols,nlyrs)     ! cm
          95             :     real(r8),intent(in) :: Ap                      ! geomagnetic activity index
          96             :     real(r8),intent(out) :: ionpairs(ncols,nlyrs)  ! pairs/cm3/sec
          97             :     integer, intent(out) :: status                 ! error status
          98             :     real(r8),intent(in), optional :: maglat(ncols) ! magnetic latitude (radians)
          99             :     real(r8),intent(in), optional :: lshell(ncols) ! magnetosphere L-Shells
         100             :
         101             :     integer :: i,k
         102           0 :     real(r8) :: flux_sd(nbins)
         103           0 :     logical  :: valid(nbins)
         104           0 :     real(r8) :: flux(nbins)
         105           0 :     real(r8) :: ipr(nbins,nlyrs)
         106           0 :     real(r8) :: l_shells(ncols)
         107             :
         108           0 :     status = 0
         109             :
         110           0 :     if (present(lshell)) then
         111           0 :        l_shells(:) = lshell(:)
         112           0 :     elseif (present(maglat)) then
         113             :        ! get L-shell values corresponeding to the column geo-mag latitudes
         114           0 :        l_shells = maglat2lshell( maglat )
         115             :     else
         116           0 :        ionpairs(:,:) = -huge(1._r8)
         117           0 :        status = mee_ap_error
         118             :        return
         119             :     endif
         120             :
         121           0 :     ionpairs(:,:) = 0._r8
         122             :
         123           0 :     do i = 1,ncols
         124             :
         125           0 :        if (  l_shells(i) >= 2._r8 .and. l_shells(i) <= 10._r8 ) then
         126             :
         127           0 :           valid(:) = .false.
         128           0 :           flux_sd(:) = 0._r8
         129             :
         130           0 :           if (mee_fluxes_active) then
         131             :              ! use prescribed top-of-atmosphere fluxes
         132             :              call mee_fluxes_extract( l_shells(i), flux_sd, valid )
         133             :           end if
         134             :
         135           0 :           where ( (.not.valid(:)) .and. (energies(:)>=30._r8) .and. (energies(:)<=1000._r8) )
         136             :              ! calculate the top of the atmosphere energetic electron energy spectrum
         137             :              ! van de Kamp is per steradian (electrons / (cm2 sr s keV))
         138             :              flux_sd(:) = FluxSpectrum(energies(:), l_shells(i), Ap)
         139             :           end where
         140             :
         141             :           ! assume flux is isotropic inside a nominal bounce loss cone (BLC) angle.
         142             :           ! The area of the BLC in sr is 2pi(1-cosd(BLC))
         143           0 :           flux(:) = solid_angle_factor*flux_sd(:)
         144             :
         145             :           ! calculate the IPR as a function of height and energy
         146           0 :           ipr(:,:) = iprmono(energies, flux, airden(i,:), scaleh(i,:))
         147             :
         148             :           ! integrate across the energy range to get total IPR
         149           0 :           do k=1,nlyrs
         150           0 :              ionpairs(i,k) = sum(ipr(:,k)*denergies(:))
         151             :           end do
         152             :
         153             :        end if
         154             :
         155             :     end do
         156             :
         157           0 :   end subroutine mee_ap_iprs
         158             :
         159             :   !------------------------------------------------------------------------------
         160             :   ! Electron fluxes for a specific L-shell and Ap
         161             :   !
         162             :   ! Based on:
         163             :   !
         164             :   ! van de Kamp, M., Seppala, A., Clilverd, M. A., Rodger, C. J., Verronen, P. T., and Whittaker, I. C. (2016),
         165             :   ! A model providing long‐term datasets of energetic electron precipitation during geomagnetic storms,
         166             :   ! J. Geophys. Res. Atmos., 121, 12,520– 12,540,
         167             :   ! [doi:10.1002/2015JD024212](https://agupubs.onlinelibrary.wiley.com/doi/full/10.1002/2015JD024212)
         168             :   !------------------------------------------------------------------------------
         169           0 :   pure function FluxSpectrum( en, lshell, Ap ) result(flux)
         170             :
         171             :     real(r8), intent(in) :: en(:)  ! electron energy bins (keV)
         172             :     real(r8), intent(in) :: lshell ! magnetosphere L-Shell number
         173             :     real(r8), intent(in) :: Ap     ! geomagnetic activity index
         174             :
         175             :     real(r8) :: flux(size(en))
         176             :
         177             :     real(r8) :: lpp
         178             :     real(r8) :: Spp
         179             :     real(r8) :: A
         180             :     real(r8) :: b
         181             :     real(r8) :: c
         182             :     real(r8) :: s
         183             :     real(r8) :: d
         184             :     real(r8) :: F30
         185             :     real(r8) :: E
         186             :     real(r8) :: bk
         187             :     real(r8) :: sk
         188             :     real(r8) :: k
         189             :     real(r8) :: x
         190             :
         191           0 :     lpp = -0.7430_r8*log(Ap) + 6.5257_r8
         192           0 :     Spp = lshell - lpp
         193             :
         194             :     ! vdK2016 eqn.(8)
         195             :
         196           0 :     A = 8.2091_r8*Ap**(0.16255_r8)
         197           0 :     b = 1.3754_r8*Ap**(0.33042_r8)
         198           0 :     c = 0.13334_r8*Ap**(0.42616_r8)
         199           0 :     s = 2.2833_r8*Ap**(-0.22990_r8)
         200           0 :     d = 2.7563e-4_r8*Ap**(2.6116_r8)
         201             :
         202           0 :     F30 = exp(A) / (exp(-b*(Spp-s)) + exp(c*(Spp-s)) + d)
         203             :
         204             :     ! vdK2016 eqn.(9)
         205             :
         206           0 :     E = 3.3777_r8*Ap**(-1.7038_r8) + 0.15_r8
         207           0 :     bk = 3.7632_r8*Ap**(-0.16034_r8)
         208           0 :     sk = 12.184_r8*Ap**(-0.30111_r8)
         209             :
         210           0 :     k = -1.0_r8 / (E*exp(-bk*Spp) + 0.30450_r8*cosh(0.20098_r8*(Spp-sk))) - 1._r8
         211             :
         212           0 :     x=k+1
         213           0 :     c = F30*x / ((1e3_r8**x) - (30._r8**x))
         214           0 :     flux(:) = c * (en(:)**k)
         215             :
         216           0 :   end function FluxSpectrum
         217             :
         218             :   !------------------------------------------------------------------------------
         219             :   ! Calculate how energy from top of atmosphere is deposited in rest of atmosphere
         220             :   !
         221             :   ! The function takes the energy spectrum at the top of the atmosphere and
         222             :   ! calculates how that energy is deposited in the atmosphere using the parameterization
         223             :   ! described in [Fang et al., (2010)](https://opensky.ucar.edu/islandora/object/articles:10653)
         224             :   !
         225             :   ! Fang, X., C. E. Randall, D. Lummerzheim, W. Wang, G. Lu, S. C. Solomon,
         226             :   ! and R. A. Frahm (2010), Parameterization of monoenergetic electron impact
         227             :   ! ionization, Geophys. Res. Lett., 37, L22106, [doi:10.1029/2010GL045406.]
         228             :   ! (https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2010GL045406)
         229             :   !
         230             :   ! Application of the new parameterization requires the following steps:
         231             :   !
         232             :   ! 1. Fang coefficients using equation (5) and Table 1. are calculated at initialization time
         233             :   ! 2. Calculate the y values throughout the atmosphere using equation (1).
         234             :   ! 3. Calculate the normalized energy dissipation f values using equation (4).
         235             :   ! 4. Obtain the altitude profile of qtot by substituting the f values into equation (3).
         236             :   !
         237             :   !------------------------------------------------------------------------------
         238           0 :   function iprmono(e, flux, rho, scaleh) result(ipr)
         239             :     real(r8), intent(in) :: e(:)      ! electron energy bins (keV)
         240             :     real(r8), intent(in) :: flux(:)   ! top of atmos electron fluxes (electrons / (cm2 sr s keV))
         241             :     real(r8), intent(in) :: rho(:)    ! density of atmos (g/cm3)
         242             :     real(r8), intent(in) :: scaleh(:) ! scale  height (cm)
         243             :
         244             :     real(r8) :: ipr(size(e),size(rho))
         245             :     integer :: nenergies, n
         246             :     integer :: nlayers, k
         247             :
         248           0 :     real(r8) :: y(size(rho))
         249           0 :     real(r8) :: f(size(rho))
         250             :     real(r8) :: Qmono
         251             :
         252             :     ! assign constants
         253             :     real(r8), parameter :: epsilon = 0.035_r8 ! keV energy loss per ion pair produced
         254             :
         255           0 :     ipr = 0._r8
         256           0 :     nenergies = size(e)
         257           0 :     nlayers = size(rho)
         258             :
         259           0 :     do n = 1,nenergies
         260             :
         261             :        ! step 1. (eq. 1)
         262           0 :        y(:) = (2/e(n))*(rho(:)*scaleh(:)/6.0e-6_r8)**0.7_r8
         263             :
         264           0 :        do k = 1,nlayers
         265           0 :           f(k) = fang(y(k), n)
         266             :        end do
         267             :
         268             :        ! calculate ipr (qtot) using eq. (3) for a specified flux at ea. energy
         269           0 :        Qmono = flux(n)*e(n) ! (keV cm−2 s−1)
         270           0 :        ipr(n,:) = f(:)*Qmono/(epsilon*scaleh(:))
         271             :     end do
         272             :
         273             :   contains
         274             :
         275           0 :     pure function fang(y,n) result(f)
         276             :       real(r8), intent(in) :: y
         277             :       integer,  intent(in) :: n ! energy ndx
         278             :
         279             :       real(r8) :: f
         280             :       ! Input:
         281             :       ! y - normalized atmospheric column mass as a function of vertical location (z)
         282             :       ! Emono - is incident electron energy (keV)
         283             :       ! Output:
         284             :       ! f - quanity calculated by eqn. (4)
         285             :
         286             :
         287             :       ! eq. (4) - Normalized energy deposition
         288           0 :       f = fang_coefs(1,n) * (y**fang_coefs(2,n)) * exp(-fang_coefs(3,n)*y**fang_coefs(4,n)) &
         289           0 :         + fang_coefs(5,n) * (y**fang_coefs(6,n)) * exp(-fang_coefs(7,n)*y**fang_coefs(8,n))
         290             :
         291           0 :     end function fang
         292             :
         293             :   end function iprmono
         294             :
         295             :   !------------------------------------------------------------------------------
         296             :   ! initializes the coeffs used in the fang function in iprmono
         297             :   !------------------------------------------------------------------------------
         298           0 :   subroutine init_fang_coefs
         299             :
         300             :     integer :: n,i, ierr
         301             :
         302             :     real(r8) :: lne, lne2, lne3
         303             :     ! Table 1. of Fang et al., 2010
         304             :     ! (https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2010GL045406)
         305             :     real(r8), parameter :: p1(8,4) = reshape( &
         306             :          (/ 1.24616E+0_r8,  1.45903E+0_r8, -2.42269E-1_r8,  5.95459E-2_r8, &
         307             :             2.23976E+0_r8, -4.22918E-7_r8,  1.36458E-2_r8,  2.53332E-3_r8, &
         308             :             1.41754E+0_r8,  1.44597E-1_r8,  1.70433E-2_r8,  6.39717E-4_r8, &
         309             :             2.48775E-1_r8, -1.50890E-1_r8,  6.30894E-9_r8,  1.23707E-3_r8, &
         310             :            -4.65119E-1_r8, -1.05081E-1_r8, -8.95701E-2_r8,  1.22450E-2_r8, &
         311             :             3.86019E-1_r8,  1.75430E-3_r8, -7.42960E-4_r8,  4.60881E-4_r8, &
         312             :            -6.45454E-1_r8,  8.49555E-4_r8, -4.28581E-2_r8, -2.99302E-3_r8, &
         313             :             9.48930E-1_r8,  1.97385E-1_r8, -2.50660E-3_r8, -2.06938E-3_r8 /), &
         314             :          shape=(/8,4/),order=(/2,1/))
         315             :
         316           0 :     allocate(fang_coefs(8,nbins), stat=ierr)
         317           0 :     if (ierr/=0) call endrun('init_fang_coefs: not able to allocate fang_coefs')
         318             :
         319           0 :     do n = 1,nbins
         320             :        ! terms in eq. (5)
         321           0 :        lne = log(energies(n))
         322           0 :        lne2 = lne*lne
         323           0 :        lne3 = lne*lne2
         324             :
         325             :        ! step 2. calculate the C array in (5)
         326           0 :        do i = 1,8
         327           0 :           fang_coefs(i,n) = exp(p1(i,1) + p1(i,2)*lne + p1(i,3)*lne2 + p1(i,4)*lne3)
         328             :        end do
         329             :
         330             :     end do
         331           0 :   end subroutine init_fang_coefs
         332             :
         333             :   !------------------------------------------------------------------------------
         334             :   ! Generate a grid of energy bins for the flux spectrum.
         335             :   ! The energy range of the spectrum is 30-1000 keV,
         336             :   ! with nbins of logarithmically spaced grid points.
         337             :   !------------------------------------------------------------------------------
         338           0 :   subroutine gen_energy_grid(nbins, energies, deltas)
         339             :     integer, intent(in) :: nbins
         340             :     real(r8),intent(out) :: energies(nbins)
         341             :     real(r8),intent(out) :: deltas(nbins)
         342             :
         343             :     integer :: i
         344             :     real(r8) :: low,med,hig
         345             :
         346             :     ! for energy bins ranging from 30 keV to 1000 keV
         347             :     real(r8), parameter :: e1 = log10(30._r8)
         348             :     real(r8), parameter :: e2 = log10(1000._r8)
         349             :
         350           0 :     do i = 1,nbins
         351           0 :        low = e1 + (e2-e1)*(i-1.0_r8)/nbins
         352           0 :        med = e1 + (e2-e1)*(i-0.5_r8)/nbins
         353           0 :        hig = e1 + (e2-e1)*(i)/nbins
         354             :
         355           0 :        energies(i) = 10**med
         356           0 :        deltas(i) = (10**hig)-(10**low)
         357             :     end do
         358             :
         359           0 :   end subroutine gen_energy_grid
         360             :
         361             :   !------------------------------------------------------------------------------
         362             :   ! returns L-Shell number for a given magnetic latitude (radians)
         363             :   !------------------------------------------------------------------------------
         364           0 :   pure elemental function maglat2lshell( rmaglat ) result(lshell)
         365             :     real(r8), intent(in) :: rmaglat ! mag latitude in radians
         366             :     real(r8) :: lshell  ! magnetosphere L-Shell number
         367             :
         368             :     ! lshell = r/(cos(rmaglat)**2) (https://en.wikipedia.org/wiki/L-shell)
         369             :     ! where r is the radial distance (in planetary radii) to a point on the line.
         370             :     ! r = 1.01 corresponds to an altitude in the lower mesosphere (~64 km)
         371             :     ! where medium-energy electrons typically deposit their energy.
         372           0 :     lshell = 1.01_r8/(cos(rmaglat)**2)
         373             :
         374           0 :   end function maglat2lshell
         375             :
         376             : end module mee_ap_util_mod

Generated by: LCOV version 1.14