LCOV - code coverage report
Current view: | top level - chemistry/mozart - mee_ionization.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 27 | 52 | 51.9 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 4 | 4 | 100.0 %

              Line data    Source code
           1             : !--------------------------------------------------------------------------
           2             : ! CAM interface layer for inline computation of atmosphere ionization rates
           3             : ! due to medium energy electrons in the magnetosphere radiation belts impacting
           4             : ! the atmosphere.  Fluxes of electrons incident on the upper atmosphere can
           5             : ! be computed based on Ap or read from file.
           6             : !--------------------------------------------------------------------------
           7             : module mee_ionization
           8             :   use shr_kind_mod, only: r8 => shr_kind_r8
           9             :   use solar_parms_data, only: Ap=>solar_parms_ap ! geomag activity index
          10             :   use mo_apex, only: alatm ! mag latitude at each column (radians)
          11             :   use ppgrid, only: pcols, pver
          12             :   use cam_logfile,  only: iulog
          13             :   use spmd_utils,   only: masterproc
          14             :   use cam_abortutils, only: endrun
          15             :   use mee_ap_util_mod,only: mee_ap_init, mee_ap_error, mee_ap_iprs
          16             :
          17             :   implicit none
          18             :
          19             :   private
          20             :   public :: mee_ion_readnl
          21             :   public :: mee_ion_init
          22             :   public :: mee_ion_final
          23             :   public :: mee_ionpairs
          24             :
          25             :   logical :: mee_ion_inline = .false.
          26             :   logical :: mee_ion_diagonly = .false.
          27             :   real(r8) :: mee_ion_blc = -huge(1._r8) ! bounce cone angle (degrees)
          28             :
          29             : contains
          30             :
          31             :   !-----------------------------------------------------------------------------
          32             :   ! reads namelist options
          33             :   !-----------------------------------------------------------------------------
          34     1490712 :   subroutine mee_ion_readnl(nlfile)
          35             :
          36             :     use namelist_utils, only: find_group_name
          37             :     use spmd_utils,     only: mpicom, mpi_logical, mpi_real8, masterprocid
          38             :     use mee_fluxes,     only: mee_fluxes_readnl
          39             :
          40             :     character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input
          41             :
          42             :     ! Local variables
          43             :     integer :: unitn, ierr
          44             :     character(len=*), parameter :: subname = 'mee_ion_readnl'
          45             :
          46             :     namelist /mee_ion_nl/ mee_ion_inline, mee_ion_blc, mee_ion_diagonly
          47             :
          48             :
          49             :     ! Read namelist
          50        1536 :     if (masterproc) then
          51           2 :        open( newunit=unitn, file=trim(nlfile), status='old' )
          52           2 :        call find_group_name(unitn, 'mee_ion_nl', status=ierr)
          53           2 :        if (ierr == 0) then
          54           0 :           read(unitn, mee_ion_nl, iostat=ierr)
          55           0 :           if (ierr /= 0) then
          56           0 :              call endrun(subname // ':: ERROR reading namelist')
          57             :           end if
          58             :        end if
          59           2 :        close(unitn)
          60             :     end if
          61             :
          62             :     ! Broadcast namelist variables
          63        1536 :     call mpi_bcast(mee_ion_inline, 1, mpi_logical, masterprocid, mpicom, ierr)
          64        1536 :     call mpi_bcast(mee_ion_blc, 1, mpi_real8, masterprocid, mpicom, ierr)
          65        1536 :     call mpi_bcast(mee_ion_diagonly, 1, mpi_logical, masterprocid, mpicom, ierr)
          66        1536 :     if ( masterproc ) then
          67           2 :        write(iulog,*) subname//':: mee_ion_inline = ', mee_ion_inline
          68           2 :        if (mee_ion_inline) then
          69           0 :           write(iulog,*) subname//':: mee_ion_blc = ', mee_ion_blc
          70           0 :           write(iulog,*) subname//':: mee_ion_diagonly = ', mee_ion_diagonly
          71             :        endif
          72             :     endif
          73             :
          74        1536 :     if (mee_ion_inline) then
          75           0 :        call mee_fluxes_readnl(nlfile)
          76             :     end if
          77             :
          78        1536 :   end subroutine mee_ion_readnl
          79             :
          80             :   !-----------------------------------------------------------------------------
          81             :   !-----------------------------------------------------------------------------
          82        1536 :   subroutine mee_ion_init()
          83        1536 :     use cam_history, only: addfld
          84             :     use mee_fluxes,  only: mee_fluxes_init
          85             :
          86             :     integer :: err
          87             :
          88        1536 :     if (.not.mee_ion_inline) return
          89             :
          90           0 :     call mee_fluxes_init()
          91           0 :     call mee_ap_init(mee_ion_blc,err)
          92           0 :     if (err==mee_ap_error) then
          93           0 :        call endrun('mee_ion_init: not able to initialize Ap based MEE ionization')
          94             :     endif
          95             :
          96           0 :     call addfld( 'APMEEionprs', (/ 'lev' /), 'A', 'pairs/cm3/sec', 'Ap generated MEE ionization rate' )
          97        1536 :   end subroutine mee_ion_init
          98             :
          99             :   !-----------------------------------------------------------------------------
         100             :   !-----------------------------------------------------------------------------
         101        1536 :   subroutine mee_ion_final()
         102        1536 :     use mee_fluxes,      only: mee_fluxes_final
         103             :     use mee_ap_util_mod, only: mee_ap_final
         104             :
         105        1536 :     if (.not.mee_ion_inline) return
         106             :
         107           0 :     call mee_fluxes_final()
         108           0 :     call mee_ap_final()
         109             :
         110        1536 :   end subroutine mee_ion_final
         111             :
         112             :   !-----------------------------------------------------------------------------
         113             :   !-----------------------------------------------------------------------------
         114     1489176 :   subroutine mee_ionpairs(ncol, lchnk, pmid, alt, temp, ionpairs)
         115             :
         116        1536 :     use air_composition, only: mbarv  ! kg/kmole
         117             :     use physconst,       only: gravit
         118             :     use air_composition, only: rairv  ! composition dependent gas constant (J/K/kg)
         119             :     use physconst,       only: boltz  ! Boltzman's constant (J/K/molecule)
         120             :     use physconst,       only: avogad ! Avogadro's number (molecules/kmole)
         121             :     use physconst,       only: rearth ! radius of earth (m)
         122             :     use cam_history,     only: outfld
         123             :
         124             :     integer,  intent(in) :: ncol,lchnk
         125             :     real(r8), intent(in) :: pmid(:,:)
         126             :     real(r8), intent(in) :: alt(:,:) ! meters
         127             :     real(r8), intent(in) :: temp(:,:)
         128             :     real(r8), intent(out) :: ionpairs(:,:) ! ion pairs /cm3/sec
         129             :
         130             :     real(r8) :: rho(pcols,pver)
         131             :     real(r8) :: scaleh(pcols,pver)
         132             :     real(r8) :: grvty(pcols,pver)
         133             :     integer :: err
         134             :
         135     1489176 :     if (.not.mee_ion_inline) then
         136  2314006344 :        ionpairs(:,:) = 0._r8
         137             :        return
         138             :     end if
         139             :
         140           0 :     rho(:ncol,:) = pmid(:ncol,:)/(rairv(:ncol,:,lchnk)*temp(:ncol,:)) ! kg/m3
         141           0 :     rho(:ncol,:) = rho(:ncol,:)*1.0e-3_r8 ! kg/m3 --> g/cm3
         142             :
         143           0 :     grvty(:ncol,:) = gravit * ( (rearth/(rearth+alt(:ncol,:)))**2 )
         144             :
         145           0 :     scaleh(:ncol,:) = avogad * boltz*temp(:ncol,:)/(mbarv(:ncol,:,lchnk)*grvty(:ncol,:)) ! m
         146           0 :     scaleh(:ncol,:) = scaleh(:ncol,:) * 1.0e2_r8 ! m -> cm
         147             :
         148           0 :     call mee_ap_iprs(ncol, pver, rho(:ncol,:), scaleh(:ncol,:), Ap, ionpairs(:ncol,:), &
         149           0 :                      status=err, maglat=alatm(:ncol,lchnk))
         150           0 :     if (err==mee_ap_error) then
         151           0 :        call endrun('mee_ionpairs: error in Ap based MEE ionization calculation')
         152             :     end if
         153             :
         154           0 :     call outfld( 'APMEEionprs', ionpairs(:ncol,:), ncol, lchnk )
         155             :
         156           0 :     if (mee_ion_diagonly) then
         157           0 :        ionpairs(:,:) = 0._r8
         158             :     end if
         159             :
         160     1489176 :   end subroutine mee_ionpairs
         161             :
         162             :
         163             : end module mee_ionization

Generated by: LCOV version 1.14