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           0 :   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           0 :     if (masterproc) then
      51           0 :        open( newunit=unitn, file=trim(nlfile), status='old' )
      52           0 :        call find_group_name(unitn, 'mee_ion_nl', status=ierr)
      53           0 :        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           0 :        close(unitn)
      60             :     end if
      61             : 
      62             :     ! Broadcast namelist variables
      63           0 :     call mpi_bcast(mee_ion_inline, 1, mpi_logical, masterprocid, mpicom, ierr)
      64           0 :     call mpi_bcast(mee_ion_blc, 1, mpi_real8, masterprocid, mpicom, ierr)
      65           0 :     call mpi_bcast(mee_ion_diagonly, 1, mpi_logical, masterprocid, mpicom, ierr)
      66           0 :     if ( masterproc ) then
      67           0 :        write(iulog,*) subname//':: mee_ion_inline = ', mee_ion_inline
      68           0 :        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           0 :     if (mee_ion_inline) then
      75           0 :        call mee_fluxes_readnl(nlfile)
      76             :     end if
      77             : 
      78           0 :   end subroutine mee_ion_readnl
      79             : 
      80             :   !-----------------------------------------------------------------------------
      81             :   !-----------------------------------------------------------------------------
      82           0 :   subroutine mee_ion_init()
      83           0 :     use cam_history, only: addfld
      84             :     use mee_fluxes,  only: mee_fluxes_init
      85             : 
      86             :     integer :: err
      87             : 
      88           0 :     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           0 :   end subroutine mee_ion_init
      98             : 
      99             :   !-----------------------------------------------------------------------------
     100             :   !-----------------------------------------------------------------------------
     101           0 :   subroutine mee_ion_final()
     102           0 :     use mee_fluxes,      only: mee_fluxes_final
     103             :     use mee_ap_util_mod, only: mee_ap_final
     104             : 
     105           0 :     if (.not.mee_ion_inline) return
     106             : 
     107           0 :     call mee_fluxes_final()
     108           0 :     call mee_ap_final()
     109             : 
     110           0 :   end subroutine mee_ion_final
     111             : 
     112             :   !-----------------------------------------------------------------------------
     113             :   !-----------------------------------------------------------------------------
     114           0 :   subroutine mee_ionpairs(ncol, lchnk, pmid, alt, temp, ionpairs)
     115             : 
     116           0 :     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           0 :     if (.not.mee_ion_inline) then
     136           0 :        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           0 :   end subroutine mee_ionpairs
     161             : 
     162             : 
     163             : end module mee_ionization