LCOV - code coverage report
Current view: top level - chemistry/mozart - epp_ionization.F90 (source / functions) Hit Total Coverage
Test: coverage.info Lines: 27 212 12.7 %
Date: 2024-12-17 22:39:59 Functions: 3 14 21.4 %

          Line data    Source code
       1             : !-------------------------------------------------------------------------------
       2             : ! Energetic Particle Precipitation (EPP) forcings module
       3             : !  Manages ionization of the atmosphere due to energetic particles, which consists of
       4             : !  solar protons events (SPE), galactic cosmic rays(GCR), medium energy electrons (MEE)
       5             : !-------------------------------------------------------------------------------
       6             : module epp_ionization
       7             :   use shr_kind_mod,   only : r8 => shr_kind_r8, cs => shr_kind_cs, cl=> shr_kind_cl
       8             :   use spmd_utils,     only : masterproc
       9             :   use cam_abortutils, only : endrun
      10             :   use cam_logfile,    only : iulog
      11             :   use ppgrid,         only : pcols, pver, begchunk, endchunk
      12             :   use phys_grid,      only : get_ncols_p
      13             :   use pio,            only : var_desc_t, file_desc_t
      14             :   use pio,            only : pio_get_var, pio_inq_varid, pio_get_att
      15             :   use pio,            only : pio_inq_varndims, pio_inq_vardimid, pio_inq_dimname, pio_inq_dimlen
      16             :   use pio,            only : PIO_NOWRITE
      17             :   use cam_pio_utils,  only : cam_pio_openfile
      18             :   use ioFileMod,      only : getfil
      19             :   use input_data_utils, only : time_coordinate
      20             : 
      21             :   implicit none
      22             :   private
      23             : 
      24             :   public :: epp_ionization_readnl  ! read namelist variables
      25             :   public :: epp_ionization_init    ! initialization
      26             :   public :: epp_ionization_adv     ! read and time/space interpolate the data
      27             :   public :: epp_ionization_ionpairs! ion pairs production rates
      28             :   public :: epp_ionization_setmag  ! update geomagnetic coordinates mapping
      29             :   public :: epp_ionization_active
      30             : 
      31             :   character(len=cl) :: epp_all_filepath = 'NONE'
      32             :   character(len=cs) :: epp_all_varname  = 'epp_ion_rates'
      33             :   character(len=cl) :: epp_mee_filepath = 'NONE'
      34             :   character(len=cs) :: epp_mee_varname  = 'iprm'
      35             :   character(len=cl) :: epp_spe_filepath = 'NONE'
      36             :   character(len=cs) :: epp_spe_varname  = 'iprp'
      37             :   character(len=cl) :: epp_gcr_filepath = 'NONE'
      38             :   character(len=cs) :: epp_gcr_varname  = 'iprg'
      39             : 
      40             :   logical, protected :: epp_ionization_active = .false.
      41             : 
      42             :   type input_obj_t
      43             :     type(file_desc_t) :: fid
      44             :     type(var_desc_t)  :: vid
      45             :     character(len=32) :: units
      46             :     integer :: nlevs = 0
      47             :     integer :: nglats = 0
      48             :     real(r8), allocatable :: press(:)
      49             :     real(r8), allocatable :: glats(:)
      50             :     real(r8), allocatable :: gwght(:,:)      ! (pcol, begchunk:endchunk)
      51             :     integer,  allocatable :: glatn(:,:)      ! (pcol, begchunk:endchunk)
      52             :     real(r8), allocatable :: indata(:,:,:,:) ! (pcol,nlevs,begchunk:endchunk,2) inputs at indexm and indexp
      53             :     type(time_coordinate) :: time_coord
      54             :   endtype input_obj_t
      55             : 
      56             :   type(input_obj_t), pointer :: epp_in => null()
      57             :   type(input_obj_t), pointer :: spe_in => null()
      58             :   type(input_obj_t), pointer :: mee_in => null()
      59             :   type(input_obj_t), pointer :: gcr_in => null()
      60             : 
      61             : contains
      62             : 
      63             :   !-----------------------------------------------------------------------------
      64             :   !-----------------------------------------------------------------------------
      65        1536 :   subroutine epp_ionization_readnl(nlfile)
      66             : 
      67             :     use namelist_utils, only: find_group_name
      68             :     use units,          only: getunit, freeunit
      69             :     use spmd_utils,     only: mpicom, mpi_character, masterprocid
      70             : 
      71             :     character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input
      72             : 
      73             :     ! Local variables
      74             :     integer :: unitn, ierr
      75             :     character(len=*), parameter :: subname = 'epp_ionization_readnl'
      76             : 
      77             :     namelist /epp_ionization_nl/ epp_all_filepath, epp_all_varname, &
      78             :          epp_mee_filepath, epp_mee_varname, epp_spe_filepath, epp_spe_varname, epp_gcr_filepath, epp_gcr_varname
      79             : 
      80             :     ! Read namelist
      81        1536 :     if (masterproc) then
      82           2 :        unitn = getunit()
      83           2 :        open( unitn, file=trim(nlfile), status='old' )
      84           2 :        call find_group_name(unitn, 'epp_ionization_nl', status=ierr)
      85           2 :        if (ierr == 0) then
      86           0 :           read(unitn, epp_ionization_nl, iostat=ierr)
      87           0 :           if (ierr /= 0) then
      88           0 :              call endrun(subname // ':: ERROR reading namelist')
      89             :           end if
      90             :        end if
      91           2 :        close(unitn)
      92           2 :        call freeunit(unitn)
      93             :     end if
      94             : 
      95             :     ! Broadcast namelist variables
      96        1536 :     call mpi_bcast(epp_all_filepath, len(epp_all_filepath), mpi_character, masterprocid, mpicom, ierr)
      97        1536 :     call mpi_bcast(epp_mee_filepath, len(epp_mee_filepath), mpi_character, masterprocid, mpicom, ierr)
      98        1536 :     call mpi_bcast(epp_spe_filepath, len(epp_spe_filepath), mpi_character, masterprocid, mpicom, ierr)
      99        1536 :     call mpi_bcast(epp_gcr_filepath, len(epp_gcr_filepath), mpi_character, masterprocid, mpicom, ierr)
     100             : 
     101        1536 :     call mpi_bcast(epp_all_varname, len(epp_all_varname), mpi_character, masterprocid, mpicom, ierr)
     102        1536 :     call mpi_bcast(epp_mee_varname, len(epp_mee_varname), mpi_character, masterprocid, mpicom, ierr)
     103        1536 :     call mpi_bcast(epp_spe_varname, len(epp_spe_varname), mpi_character, masterprocid, mpicom, ierr)
     104        1536 :     call mpi_bcast(epp_gcr_varname, len(epp_gcr_varname), mpi_character, masterprocid, mpicom, ierr)
     105             : 
     106        1536 :     epp_ionization_active = epp_all_filepath /= 'NONE'
     107        1536 :     epp_ionization_active = epp_mee_filepath /= 'NONE' .or. epp_ionization_active
     108        1536 :     epp_ionization_active = epp_spe_filepath /= 'NONE' .or. epp_ionization_active
     109        1536 :     epp_ionization_active = epp_gcr_filepath /= 'NONE' .or. epp_ionization_active
     110             : 
     111        1536 :     if ( epp_ionization_active .and. masterproc ) then
     112           0 :        write(iulog,*) subname//':: epp_all_filepath = '//trim(epp_all_filepath)
     113           0 :        write(iulog,*) subname//':: epp_mee_filepath = '//trim(epp_mee_filepath)
     114           0 :        write(iulog,*) subname//':: epp_spe_filepath = '//trim(epp_spe_filepath)
     115           0 :        write(iulog,*) subname//':: epp_gcr_filepath = '//trim(epp_gcr_filepath)
     116             :     endif
     117             : 
     118        1536 :   end subroutine epp_ionization_readnl
     119             : 
     120             :   !-----------------------------------------------------------------------------
     121             :   !-----------------------------------------------------------------------------
     122           0 :   subroutine epp_ionization_init()
     123             :     use cam_history, only : addfld
     124             : 
     125             :     character(len=32) :: fldunits
     126           0 :     fldunits = ''
     127             : 
     128           0 :     if (epp_all_filepath /= 'NONE') then
     129           0 :        epp_in => create_input_obj(epp_all_filepath,epp_all_varname)
     130           0 :        fldunits = trim(epp_in%units)
     131             :     else
     132           0 :        if (epp_mee_filepath /= 'NONE') then
     133           0 :           mee_in => create_input_obj(epp_mee_filepath,epp_mee_varname)
     134           0 :           fldunits = trim(mee_in%units)
     135             :        endif
     136           0 :        if (epp_spe_filepath /= 'NONE') then
     137           0 :           spe_in => create_input_obj(epp_spe_filepath,epp_spe_varname)
     138           0 :           fldunits = trim(spe_in%units)
     139             :        endif
     140           0 :        if (epp_gcr_filepath /= 'NONE') then
     141           0 :           gcr_in => create_input_obj(epp_gcr_filepath,epp_gcr_varname)
     142           0 :           fldunits = trim(gcr_in%units)
     143             :        endif
     144             :     endif
     145           0 :     call addfld( 'EPPions', (/ 'lev' /), 'A', fldunits, 'EPP ionization data' )
     146             : 
     147           0 :   end subroutine epp_ionization_init
     148             : 
     149             :   !-----------------------------------------------------------------------------
     150             :   !-----------------------------------------------------------------------------
     151           0 :   subroutine epp_ionization_setmag( maglat )
     152             :     real(r8), intent(in) :: maglat(pcols,begchunk:endchunk)
     153             : 
     154           0 :     if (.not.epp_ionization_active) return
     155             : 
     156           0 :     if ( associated(epp_in) ) then
     157           0 :        call set_wghts(maglat,epp_in)
     158             :     else
     159           0 :        if ( associated(mee_in) ) then
     160           0 :           call set_wghts(maglat,mee_in)
     161             :        endif
     162           0 :        if ( associated(spe_in) ) then
     163           0 :           call set_wghts(maglat,spe_in)
     164             :        endif
     165           0 :        if ( associated(gcr_in) ) then
     166           0 :           call set_wghts(maglat,gcr_in)
     167             :        endif
     168             :     endif
     169             : 
     170           0 :   end subroutine epp_ionization_setmag
     171             : 
     172             :   !-----------------------------------------------------------------------------
     173             :   !-----------------------------------------------------------------------------
     174      370944 :   subroutine epp_ionization_adv
     175             : 
     176      370944 :     if (.not.epp_ionization_active) return
     177             : 
     178           0 :     if ( associated(epp_in) ) then
     179           0 :        call update_input(epp_in)
     180             :     else
     181           0 :        if ( associated(spe_in) ) then
     182           0 :           call update_input(spe_in)
     183             :        endif
     184           0 :        if ( associated(gcr_in) ) then
     185           0 :           call update_input(gcr_in)
     186             :        endif
     187           0 :        if ( associated(mee_in) ) then
     188           0 :           call update_input(mee_in)
     189             :        endif
     190             :     endif
     191             : 
     192             :   end subroutine epp_ionization_adv
     193             : 
     194             :   !-----------------------------------------------------------------------------
     195             :   !-----------------------------------------------------------------------------
     196     1489176 :   subroutine epp_ionization_ionpairs( ncol, lchnk, pmid, temp, ionpairs )
     197             : 
     198             :     integer,  intent(in) :: ncol, lchnk
     199             :     real(r8), intent(in) :: pmid(:,:), temp(:,:)
     200             :     real(r8), intent(out) :: ionpairs(:,:) ! ion pair production rate
     201             : 
     202  2314006344 :     ionpairs = 0._r8
     203     1489176 :     if (.not.epp_ionization_active) return
     204             : 
     205           0 :     if ( associated(epp_in) ) then
     206           0 :        ionpairs(:ncol,:) = ionpairs(:ncol,:) + interp_ionpairs( ncol, lchnk, pmid, temp, epp_in )
     207             :     else
     208           0 :        if ( associated(spe_in) ) then
     209           0 :           ionpairs(:ncol,:) = ionpairs(:ncol,:) + interp_ionpairs( ncol, lchnk, pmid, temp, spe_in )
     210             :        endif
     211           0 :        if ( associated(gcr_in) ) then
     212           0 :           ionpairs(:ncol,:) = ionpairs(:ncol,:) + interp_ionpairs( ncol, lchnk, pmid, temp, gcr_in )
     213             :        endif
     214           0 :        if ( associated(mee_in) ) then
     215           0 :           ionpairs(:ncol,:) = ionpairs(:ncol,:) + interp_ionpairs( ncol, lchnk, pmid, temp, mee_in )
     216             :        endif
     217             :     endif
     218             : 
     219             :   end subroutine epp_ionization_ionpairs
     220             : 
     221             :   ! private methods
     222             :   !-----------------------------------------------------------------------------
     223             :   !-----------------------------------------------------------------------------
     224           0 :   subroutine update_input( input )
     225             :     type(input_obj_t), pointer :: input
     226             : 
     227           0 :     if ( input%time_coord%read_more() ) then
     228           0 :        call input%time_coord%advance()
     229           0 :        call read_next_data( input )
     230             :     else
     231           0 :        call input%time_coord%advance()
     232             :     endif
     233             : 
     234           0 :   end subroutine update_input
     235             : 
     236             :   !-----------------------------------------------------------------------------
     237             :   !-----------------------------------------------------------------------------
     238           0 :   subroutine read_next_data( input )
     239             :     type(input_obj_t), pointer :: input
     240             : 
     241             :     ! read data corresponding surrounding time indices
     242           0 :     if ( input%nglats > 0 ) then
     243           0 :       call read_2d_profile( input )
     244             :     else
     245           0 :       call read_1d_profile( input )
     246             :     endif
     247             : 
     248           0 :   end subroutine read_next_data
     249             : 
     250             :   !-----------------------------------------------------------------------------
     251             :   !-----------------------------------------------------------------------------
     252           0 :   function interp_ionpairs( ncol, lchnk, pmid, temp, input ) result( ionpairs )
     253             :     use interpolate_data, only : lininterp, lininterp_init, lininterp_finish, extrap_method_zero, interp_type
     254             :     use air_composition,  only : rairv
     255             :     use cam_history,      only : outfld
     256             : 
     257             :     integer,  intent(in) :: ncol, lchnk
     258             :     real(r8), intent(in) :: pmid(:,:) ! Pa
     259             :     real(r8), intent(in) :: temp(:,:) ! K
     260             :     type(input_obj_t), pointer :: input
     261             :     real(r8) :: ionpairs(ncol,pver)
     262             : 
     263             :     real(r8) :: fctr1, fctr2
     264           0 :     real(r8) :: wrk(ncol,input%nlevs)
     265           0 :     real(r8) :: ions_diags(ncol,pver) ! for diagnostics
     266             :     integer :: i
     267             :     type (interp_type) :: interp_wgts
     268             : 
     269           0 :     if (input%time_coord%time_interp) then
     270             :        ! time interpolate
     271           0 :        fctr1 = input%time_coord%wghts(1)
     272           0 :        fctr2 = input%time_coord%wghts(2)
     273           0 :        wrk(:ncol,:) = fctr1*input%indata(:ncol,:,lchnk,1) + fctr2*input%indata(:ncol,:,lchnk,2)
     274             :     else
     275           0 :        wrk(:ncol,:) = input%indata(:ncol,:,lchnk,1)
     276             :     endif
     277             : 
     278             :     ! vertical interpolate ...
     279             :     ! interpolate to model levels
     280           0 :     do i = 1,ncol
     281             : 
     282             :        ! interpolate over log pressure
     283           0 :        call lininterp_init( log(input%press(:input%nlevs)*1.e2_r8), input%nlevs, &
     284           0 :                             log(pmid(i,:pver)), pver, extrap_method_zero, interp_wgts )
     285           0 :        call lininterp( wrk(i,:input%nlevs), input%nlevs, &
     286           0 :                        ionpairs(i,:pver), pver, interp_wgts )
     287           0 :        call lininterp_finish(interp_wgts)
     288             : 
     289           0 :        ions_diags(i,:pver) = ionpairs(i,:pver)
     290             : 
     291           0 :        if ( index(trim(input%units), 'g^-1') > 0 ) then
     292             :           ! convert to ionpairs/cm3/sec
     293           0 :           ionpairs(i,:pver) = ionpairs(i,:pver) *(1.e-3_r8*pmid(i,:pver)/(rairv(i,:pver,lchnk)*temp(i,:pver)))
     294             :        endif
     295             :     enddo
     296             : 
     297           0 :     call outfld( 'EPPions', ions_diags(:ncol,:), ncol, lchnk )
     298             : 
     299           0 :   end function interp_ionpairs
     300             : 
     301             :   !-----------------------------------------------------------------------------
     302             :   ! read 2D profile (geomag-lat vs press) and transfer to geographic grid
     303             :   !-----------------------------------------------------------------------------
     304           0 :   subroutine read_2d_profile( input )
     305             : 
     306             :     type(input_obj_t), pointer :: input
     307             : 
     308             :     ! local vars
     309           0 :     real(r8) :: wrk2d( input%nglats, input%nlevs, 2 )
     310             :     integer  :: t, c, i, ntimes, ncols, ierr
     311             :     real(r8) :: wght1, wght2
     312             :     integer  :: gndx1, gndx2
     313             :     integer  :: cnt(3), strt(3)
     314             : 
     315           0 :     if (input%time_coord%time_interp) then
     316             :        ntimes = 2
     317             :     else
     318           0 :        ntimes = 1
     319             :     endif
     320             : 
     321           0 :     cnt(1) = input%nglats
     322           0 :     cnt(2) = input%nlevs
     323           0 :     cnt(3) = ntimes
     324             : 
     325           0 :     strt(:) = 1
     326           0 :     strt(3) = input%time_coord%indxs(1)
     327             : 
     328           0 :     ierr = pio_get_var( input%fid, input%vid, strt, cnt, wrk2d )
     329             : 
     330           0 :     do t = 1,ntimes
     331           0 :        do c=begchunk,endchunk
     332           0 :           ncols = get_ncols_p(c)
     333           0 :           do i = 1,ncols
     334           0 :              gndx1 = input%glatn(i,c)
     335           0 :              if (gndx1>0) then
     336           0 :                 wght1 = input%gwght(i,c)
     337           0 :                 gndx2 = gndx1+1
     338           0 :                 if (gndx2.le.input%nglats) then
     339           0 :                   wght2 = 1._r8-wght1
     340           0 :                   input%indata(i,:,c,t) = wght1*wrk2d(gndx1,:,t) &
     341           0 :                                         + wght2*wrk2d(gndx2,:,t)
     342             :                 else
     343           0 :                   input%indata(i,:,c,t) = wght1*wrk2d(gndx1,:,t)
     344             :                 endif
     345             :              else
     346           0 :                 input%indata(i,:,c,t) = 0._r8
     347             :              endif
     348             :           end do
     349             :        end do
     350             :     end do
     351             : 
     352           0 :   end subroutine read_2d_profile
     353             : 
     354             :   !-----------------------------------------------------------------------------
     355             :   ! read 1D vertical profile and transfer to geographic grid poleward of 60 degrees geomag-lat
     356             :   !-----------------------------------------------------------------------------
     357           0 :   subroutine read_1d_profile( input )
     358             : 
     359             :     type(input_obj_t), pointer :: input
     360             : 
     361             :     ! local vars
     362           0 :     real(r8) :: wrk( input%nlevs, 2 )
     363             :     integer  :: t, c, i, ntimes, ncols, ierr
     364             :     integer  :: cnt(2), strt(2)
     365             : 
     366           0 :     if (input%time_coord%time_interp) then
     367             :        ntimes = 2
     368             :     else
     369           0 :        ntimes = 1
     370             :     endif
     371             : 
     372           0 :     cnt(1) = input%nlevs
     373           0 :     cnt(2) = ntimes
     374             : 
     375           0 :     strt(:) = 1
     376           0 :     strt(2) = input%time_coord%indxs(1)
     377             : 
     378           0 :     ierr = pio_get_var( input%fid, input%vid, strt, cnt, wrk )
     379             : 
     380           0 :     do t = 1,ntimes
     381           0 :        do c=begchunk,endchunk
     382           0 :           ncols = get_ncols_p(c)
     383           0 :           do i = 1,ncols
     384           0 :              input%indata(i,:,c,t) = input%gwght(i,c)*wrk(:,t)
     385             :           end do
     386             :        end do
     387             :     end do
     388             : 
     389           0 :   end subroutine read_1d_profile
     390             : 
     391             :   !-----------------------------------------------------------------------------
     392             :   !-----------------------------------------------------------------------------
     393           0 :   function create_input_obj( path, varname ) result(in_obj)
     394             :     use infnan,  only : nan, assignment(=)
     395             : 
     396             :     character(*), intent(in) :: path
     397             :     character(*), intent(in) :: varname
     398             :     type(input_obj_t), pointer :: in_obj
     399             : 
     400             :     character(len=cl) :: filen
     401             :     character(len=cl) :: data_units
     402             :     character(len=cs) :: dimname
     403             :     integer :: i, ierr
     404           0 :     integer, allocatable :: dimids(:)
     405             :     integer :: pres_did, pres_vid, glat_did, glat_vid, ndims
     406             : 
     407           0 :     if (path .eq. 'NONE') return
     408             : 
     409           0 :     allocate(in_obj)
     410             : 
     411           0 :     call in_obj%time_coord%initialize( path )
     412             : 
     413           0 :     call getfil( path, filen, 0 )
     414           0 :     call cam_pio_openfile( in_obj%fid, filen, PIO_NOWRITE )
     415             : 
     416           0 :     ierr = pio_inq_varid( in_obj%fid, varname, in_obj%vid )
     417             : 
     418           0 :     ierr = pio_get_att( in_obj%fid, in_obj%vid, 'units', data_units)
     419           0 :     in_obj%units = trim(data_units(1:32))
     420             : 
     421           0 :     ierr = pio_inq_varndims( in_obj%fid, in_obj%vid, ndims )
     422           0 :     allocate( dimids(ndims) )
     423             : 
     424           0 :     ierr = pio_inq_vardimid( in_obj%fid, in_obj%vid, dimids)
     425           0 :     pres_did = -1
     426           0 :     glat_did = -1
     427           0 :     do i = 1,ndims
     428           0 :        ierr = pio_inq_dimname( in_obj%fid, dimids(i), dimname )
     429           0 :        select case( trim(dimname(1:4)) )
     430             :          case ( 'pres', 'lev', 'plev' )
     431           0 :            pres_did = dimids(i)
     432           0 :            ierr = pio_inq_varid( in_obj%fid, dimname, pres_vid)
     433             :          case ( 'glat' )
     434           0 :            glat_did =  dimids(i)
     435           0 :            ierr = pio_inq_varid( in_obj%fid, dimname, glat_vid)
     436             :          case default
     437             :        end select
     438             :     end do
     439             : 
     440           0 :     deallocate( dimids )
     441             : 
     442           0 :     if (pres_did>0) then
     443           0 :        ierr = pio_inq_dimlen( in_obj%fid, pres_did, in_obj%nlevs )
     444           0 :        allocate( in_obj%press(in_obj%nlevs) )
     445           0 :        ierr = pio_get_var( in_obj%fid, pres_vid, in_obj%press )
     446             :     endif
     447           0 :     if (glat_did>0) then
     448           0 :        ierr = pio_inq_dimlen( in_obj%fid, glat_did, in_obj%nglats )
     449           0 :        allocate( in_obj%glats(in_obj%nglats) )
     450           0 :        ierr = pio_get_var( in_obj%fid, glat_vid, in_obj%glats )
     451           0 :        allocate( in_obj%glatn(pcols,begchunk:endchunk) )
     452             :     endif
     453             : 
     454           0 :     allocate( in_obj%gwght(pcols,begchunk:endchunk) )
     455             : 
     456           0 :     if (in_obj%time_coord%time_interp) then
     457           0 :        allocate( in_obj%indata(pcols,in_obj%nlevs,begchunk:endchunk,2) )
     458             :     else
     459           0 :        allocate( in_obj%indata(pcols,in_obj%nlevs,begchunk:endchunk,1) )
     460             :     endif
     461           0 :     in_obj%indata = nan
     462             : 
     463           0 :   end function create_input_obj
     464             : 
     465             :   !-----------------------------------------------------------------------
     466             :   !-----------------------------------------------------------------------
     467           0 :   subroutine set_wghts( maglat, input )
     468             : 
     469             :     real(r8), intent(in) :: maglat(pcols,begchunk:endchunk)
     470             :     type(input_obj_t), pointer :: input
     471             : 
     472             :     integer :: i, c, ncols, imag
     473             : 
     474           0 :     if (input%nglats>1) then ! read in general EPP 2D ionpairs production rates
     475           0 :        do c = begchunk,endchunk
     476           0 :           ncols = get_ncols_p(c)
     477           0 :           col_loop: do i = 1,ncols
     478           0 :              if ( maglat(i,c) .lt.  input%glats(1) ) then
     479           0 :                 input%glatn(i,c) = 1
     480           0 :                 input%gwght(i,c) = 1._r8
     481           0 :              elseif ( maglat(i,c) .gt. input%glats(input%nglats) ) then
     482           0 :                 input%glatn(i,c) = input%nglats
     483           0 :                 input%gwght(i,c) = 1._r8
     484             :              else
     485           0 :                 mag_loop: do imag = 1,input%nglats-1
     486           0 :                    if ( maglat(i,c) .ge. input%glats(imag) .and. &
     487           0 :                         maglat(i,c) .lt. input%glats(imag+1) ) then
     488           0 :                       input%gwght(i,c) = (input%glats(imag+1)-maglat(i,c) ) &
     489           0 :                                        / (input%glats(imag+1)-input%glats(imag))
     490           0 :                       input%glatn(i,c) = imag
     491           0 :                       exit mag_loop
     492             :                    endif
     493             :                 enddo mag_loop
     494             :              endif
     495             :           enddo col_loop
     496             :        enddo
     497             :     else ! read in 1D SPE ionpairs profile ...
     498           0 :        do c = begchunk,endchunk
     499           0 :           ncols = get_ncols_p(c)
     500           0 :           do i = 1,ncols
     501           0 :              if ( abs(maglat(i,c)) .ge. 60._r8 ) then ! poleward of 60 degrees
     502           0 :                 input%gwght(i,c) = 1._r8
     503             :              else
     504           0 :                 input%gwght(i,c) = 0._r8
     505             :              endif
     506             :           enddo
     507             :        enddo
     508             :     endif
     509             : 
     510           0 :     call read_next_data( input ) ! update the inputs when wghts are updated
     511             : 
     512           0 :   end subroutine set_wghts
     513             : 
     514           0 : end module epp_ionization

Generated by: LCOV version 1.14