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