LCOV - code coverage report
Current view: | top level - chemistry/mozart - mo_setext.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 67 | 130 | 51.5 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 2 | 2 | 100.0 %

              Line data    Source code
           1             : module mo_setext
           2             :
           3             :   use shr_kind_mod, only : r8 => shr_kind_r8
           4             :   use shr_const_mod,only : pi => shr_const_pi
           5             :   use cam_logfile,  only : iulog
           6             :
           7             :   private
           8             :   public :: setext_inti, setext, has_ions
           9             :
          10             :   save
          11             :
          12             :   integer :: co_ndx, no_ndx, xno_ndx, o_ndx
          13             :   integer :: op_ndx, o2p_ndx, np_ndx, n2p_ndx, n2d_ndx, n_ndx, e_ndx, oh_ndx
          14             :   logical :: has_ions = .false.
          15             :   logical :: has_dregion_ions = .false.
          16             :   integer :: aoa_nh_ndx
          17             :
          18             : contains
          19             :
          20     1490712 :   subroutine setext_inti
          21             :     !--------------------------------------------------------
          22             :     !   ... Initialize the external forcing module
          23             :     !--------------------------------------------------------
          24             :
          25             :     use mo_chem_utls, only : get_extfrc_ndx, get_spc_ndx
          26             :     use cam_history,  only : addfld
          27             :     use spmd_utils,   only : masterproc
          28             :     use mee_ionization,only : mee_ion_init
          29             :
          30             :     implicit none
          31             :
          32        1536 :     co_ndx    = get_extfrc_ndx( 'CO' )
          33        1536 :     no_ndx    = get_extfrc_ndx( 'NO' )
          34        1536 :     xno_ndx   = get_extfrc_ndx( 'XNO' )
          35        1536 :     aoa_nh_ndx = get_extfrc_ndx('AOA_NH')
          36             :
          37        1536 :     op_ndx   = get_extfrc_ndx( 'Op' )
          38        1536 :     o2p_ndx  = get_extfrc_ndx( 'O2p' )
          39        1536 :     np_ndx   = get_extfrc_ndx( 'Np' )
          40        1536 :     n2p_ndx  = get_extfrc_ndx( 'N2p' )
          41        1536 :     n2d_ndx  = get_extfrc_ndx( 'N2D' )
          42        1536 :     n_ndx    = get_extfrc_ndx( 'N' )
          43        1536 :     e_ndx    = get_extfrc_ndx( 'e' )
          44        1536 :     oh_ndx   = get_extfrc_ndx( 'OH' )
          45        1536 :     o_ndx    = get_extfrc_ndx( 'O' )
          46             :
          47        1536 :     has_ions = op_ndx > 0 .and. o2p_ndx > 0 .and. np_ndx > 0 .and. n2p_ndx > 0 .and. e_ndx > 0
          48        1536 :     has_dregion_ions = has_ions .and. (get_spc_ndx( 'O2m' )>0)
          49             :
          50        1536 :     if (masterproc) then
          51           2 :        write(iulog,*) ' '
          52           2 :        write(iulog,*) 'setext_inti: diagnostics: co_ndx, no_ndx, xno_ndx'
          53           2 :        write(iulog,'(10i5)') co_ndx, no_ndx, xno_ndx
          54             :     endif
          55             :
          56        3072 :     call addfld( 'NO_Lightning', (/ 'lev' /), 'A','molec/cm3/s', 'lightning NO source' )
          57        3072 :     call addfld( 'NO_Aircraft',  (/ 'lev' /), 'A', 'molec/cm3/s', 'aircraft NO source' )
          58        3072 :     call addfld( 'CO_Aircraft',  (/ 'lev' /), 'A', 'molec/cm3/s', 'aircraft CO source' )
          59             :
          60        3072 :     call addfld( 'EPP_ionpairs', (/ 'lev' /), 'A', 'pairs/cm3/s', 'EPP ionization forcing' )
          61        3072 :     call addfld( 'GCR_ionpairs', (/ 'lev' /), 'A', 'pairs/cm3/s', 'GCR ionization forcing' )
          62             :
          63        1536 :     if (.not.has_dregion_ions) then
          64        1536 :        if ( n2d_ndx > 0 .and. n_ndx>0 ) then
          65           0 :           call addfld( 'N4S_EPP', (/ 'lev' /), 'I', 'molec/cm3/s', 'solar proton event N(4S) source' )
          66           0 :           call addfld( 'N2D_EPP', (/ 'lev' /), 'I', 'molec/cm3/s', 'solar proton event N(2D) source' )
          67        1536 :        elseif ( no_ndx > 0 .and. n_ndx>0 ) then
          68           0 :           call addfld( 'N4S_EPP', (/ 'lev' /), 'I', 'molec/cm3/s', 'solar proton event N(4S) source' )
          69           0 :           call addfld( 'NO_EPP',  (/ 'lev' /), 'I',  'molec/cm3/s', 'solar proton event NO source' )
          70             :        endif
          71        1536 :        if ( oh_ndx > 0 ) then
          72           0 :           call addfld( 'OH_EPP',  (/ 'lev' /), 'I',  'molec/cm3/s', 'solar proton event OH source' )
          73             :        endif
          74             :     endif
          75        1536 :     if ( has_ions ) then
          76           0 :        call addfld( 'P_Op',    (/ 'lev' /), 'I',   '/s', 'production o+' )
          77           0 :        call addfld( 'P_O2p',   (/ 'lev' /), 'I',  '/s', 'production o2+' )
          78           0 :        call addfld( 'P_N2p',   (/ 'lev' /), 'I',  '/s', 'production n2+' )
          79           0 :        call addfld( 'P_Np',    (/ 'lev' /), 'I',   '/s', 'production n+' )
          80           0 :        call addfld( 'P_IONS',  (/ 'lev' /), 'I', '/s', 'total ion production' )
          81             :     endif
          82             :
          83        1536 :     if ( aoa_nh_ndx > 0 ) then
          84           0 :        call addfld('AOA_NH_XFRC', (/ 'lev' /), 'A', 'molec/cm3/s',  'external forcing for AOA_NH' )
          85             :     endif
          86             :
          87        1536 :     call mee_ion_init()
          88             :
          89        1536 :   end subroutine setext_inti
          90             :
          91     1489176 :   subroutine setext( extfrc, zint_abs, zint_rel, cldtop, &
          92     1489176 :        zmid, lchnk, tfld, o2mmr, ommr, &
          93     1489176 :        pmid, mbar, rlats, calday, ncol, rlons, pbuf )
          94             :     !--------------------------------------------------------
          95             :     !     ... for this latitude slice:
          96             :     !         - form the production from datasets
          97             :     !         - form the nox (xnox) production from lighing
          98             :     !         - form the nox (xnox) production from airplanes
          99             :     !         - form the co production from airplanes
         100             :     !--------------------------------------------------------
         101             :
         102        1536 :     use cam_history,  only : outfld
         103             :     use ppgrid,       only : pver, pcols
         104             :     use mo_airplane,  only : airpl_set
         105             :     use chem_mods,    only : extcnt
         106             :     use mo_lightning, only : prod_no
         107             :
         108             :     use mo_extfrc,    only : extfrc_set
         109             :     use hco_cc_emissions, only : hco_set_extfrc
         110             :     use chem_mods,    only : extcnt
         111             :     use tracer_srcs,  only : num_tracer_srcs, tracer_src_flds, get_srcs_data
         112             :     use mo_chem_utls, only : get_extfrc_ndx
         113             :
         114             :     use mo_aurora,      only : aurora
         115             :     use gcr_ionization, only : gcr_ionization_ionpairs
         116             :     use epp_ionization, only : epp_ionization_ionpairs
         117             :     use mee_ionization, only : mee_ionpairs
         118             :     use spehox,         only : hox_prod_factor
         119             :
         120             :     use physics_buffer, only : physics_buffer_desc
         121             :
         122             :     use phys_control,   only : use_hemco                   ! Use Harmonized Emissions Component (HEMCO)
         123             :
         124             :     implicit none
         125             :
         126             :     !--------------------------------------------------------
         127             :     !     ... dummy arguments
         128             :     !--------------------------------------------------------
         129             :     !--------------------------------------------------------
         130             :     integer,  intent(in)  ::   lchnk                       ! chunk id
         131             :     integer,  intent(in)  ::   ncol                        ! columns in chunk
         132             :     real(r8), intent(in)  ::   zint_abs(ncol,pver+1)           ! interface geopot height ( km )
         133             :     real(r8), intent(in)  ::   zint_rel(ncol,pver+1)           ! interface geopot height ( km )
         134             :     real(r8), intent(in)  ::   cldtop(ncol)                ! cloud top index
         135             :     real(r8), intent(out) ::   extfrc(ncol,pver,extcnt)    ! the "extraneous" forcing
         136             :
         137             :     real(r8), intent(in)  ::   calday                      ! calendar day of year
         138             :     real(r8), intent(in)  ::   rlats(ncol)                 ! column latitudes (radians)
         139             :     real(r8), intent(in)  ::   rlons(ncol)                 ! column longitudes (radians)
         140             :     real(r8), intent(in)  ::   zmid(ncol,pver)             ! midpoint geopot height ( km )
         141             :     real(r8), intent(in)  ::   pmid(pcols,pver)            ! midpoint pressure (Pa)
         142             :     real(r8), intent(in)  ::   tfld(pcols,pver)            ! midpoint temperature (K)
         143             :     real(r8), intent(in)  ::   mbar(ncol,pver)             ! mean molecular mass (g/mole)
         144             :     real(r8), intent(in)  ::   o2mmr(ncol,pver)            ! o2 concentration (kg/kg)
         145             :     real(r8), intent(in)  ::   ommr(ncol,pver)             ! o concentration (kg/kg)
         146             :
         147             :     type(physics_buffer_desc),pointer :: pbuf(:)
         148             :
         149             :     !--------------------------------------------------------
         150             :     !     ... local variables
         151             :     !--------------------------------------------------------
         152             :     integer :: i, k, cldind
         153     2978352 :     real(r8) :: srcs_offline( ncol, pver )
         154             :     integer :: ndx
         155             :
         156     2978352 :     real(r8), dimension(ncol,pver) :: no_lgt
         157             :
         158     2978352 :     real(r8) :: gcr_ipr(ncol,pver)        ! ion pairs production rate associated with galactic comsic rays
         159     2978352 :     real(r8) :: epp_ipr(ncol,pver)        ! ion pairs production rate associated with energetic particles
         160     2978352 :     real(r8) :: epp_hox(ncol,pver)        ! HOx production rate associated with energetic particles
         161             :
         162             :     real(r8), parameter :: rad2deg = 180._r8/pi                ! radians to degrees conversion factor
         163             :     real(r8) :: xlat
         164             :
         165     2978352 :     real(r8) :: mee_ap_ipr(ncol,pver) ! ion pairs production rate from Ap formulation
         166             :
         167  2314006344 :     call mee_ionpairs(ncol, lchnk, pmid, zmid*1.e3_r8, tfld, mee_ap_ipr)
         168             :
         169 25455558960 :     extfrc(:,:,:) = 0._r8
         170             :
         171  2314006344 :     no_lgt(:,:) = 0._r8
         172             :
         173     1489176 :     if(use_hemco) then
         174             :        !--------------------------------------------------------
         175             :        !     ... set frcing from datasets (HEMCO)
         176             :        !--------------------------------------------------------
         177           0 :        call hco_set_extfrc( lchnk, zint_rel, extfrc, ncol, pbuf )
         178             :     else
         179             :        !--------------------------------------------------------
         180             :        !     ... set frcing from datasets
         181             :        !--------------------------------------------------------
         182     1489176 :        call extfrc_set( lchnk, zint_rel, extfrc, ncol )
         183             :     endif
         184             :
         185             :     !--------------------------------------------------------
         186             :     !     ... set nox production from lighting
         187             :     !         note: from ground to cloud top production is c shaped
         188             :     !--------------------------------------------------------
         189     1489176 :     if ( no_ndx > 0 ) then
         190           0 :        do i = 1,ncol
         191           0 :           cldind = nint( cldtop(i) )
         192           0 :           if( cldind < pver .and. cldind > 0  ) then
         193           0 :              extfrc(i,cldind:pver,no_ndx) = extfrc(i,cldind:pver,no_ndx) &
         194           0 :                   + prod_no(i,cldind:pver,lchnk)
         195           0 :              no_lgt(i,cldind:pver) = prod_no(i,cldind:pver,lchnk)
         196             :           end if
         197             :        end do
         198             :     endif
         199     1489176 :     if ( xno_ndx > 0 ) then
         200           0 :        do i = 1,ncol
         201           0 :           cldind = nint( cldtop(i) )
         202           0 :           if( cldind < pver .and. cldind > 0  ) then
         203           0 :              extfrc(i,cldind:pver,xno_ndx) = extfrc(i,cldind:pver,xno_ndx) &
         204           0 :                   + prod_no(i,cldind:pver,lchnk)
         205             :           end if
         206             :        end do
         207             :     endif
         208             :
         209     1489176 :     call outfld( 'NO_Lightning', no_lgt(:ncol,:), ncol, lchnk )
         210             :
         211     1489176 :     call airpl_set( lchnk, ncol, no_ndx, co_ndx, xno_ndx, cldtop, zint_abs, extfrc)
         212             :
         213     1489176 :     do i = 1,num_tracer_srcs
         214             :
         215           0 :        ndx =  get_extfrc_ndx( tracer_src_flds(i) )
         216           0 :        call get_srcs_data( tracer_src_flds(i), srcs_offline,  ncol, lchnk, pbuf )
         217     1489176 :        do k = 1,pver
         218           0 :           extfrc(:ncol,k,ndx) = extfrc(:ncol,k,ndx) + srcs_offline(:ncol,k)
         219             :        enddo
         220             :
         221             :     enddo
         222             :
         223             : !
         224             : ! CCMI : external forcing for AOA_NH
         225             : ! set to 100 ppb per year
         226             : !
         227     1489176 :     if ( aoa_nh_ndx > 0 ) then
         228           0 :        extfrc(:ncol,:,aoa_nh_ndx) = 1.e-7_r8/(86400._r8*365._r8)
         229           0 :        do i=1,ncol
         230           0 :          xlat = rlats(i)*rad2deg              ! convert to degrees
         231           0 :          if ( xlat >= 30._r8 .and. xlat <= 50._r8 ) then
         232           0 :             extfrc(i,pver,aoa_nh_ndx) = 0._r8
         233             :          end if
         234             :        end do
         235           0 :        call outfld('AOA_NH_XFRC',extfrc(:ncol,:,aoa_nh_ndx), ncol, lchnk )
         236             :     end if
         237             :
         238     1489176 :     if ( has_ions ) then
         239             :        !---------------------------------------------------------------------
         240             :        !     ... set ion auroral production
         241             :        !---------------------------------------------------------------------
         242             :
         243             :        call aurora( tfld, o2mmr, ommr, mbar, rlats, &
         244             :             extfrc(:,:,o2p_ndx), extfrc(:,:,op_ndx), extfrc(:,:,n2p_ndx), extfrc(:,:,np_ndx), pmid, &
         245           0 :             lchnk, calday, ncol, rlons, pbuf )
         246             :        !---------------------------------------------------------------------
         247             :        !     ... set n(2d) and n(4s) auroral production
         248             :        !         Stan Solomon HAO
         249             :        !         include production of N by secondary auroral hot electrons (e_s*):
         250             :        !         this is not a "real" reaction; instead, the production is parameterized in terms
         251             :        !         of the production rate of N2+ by primary electrons, QN2P (which is in the model),
         252             :        !         as follows:
         253             :        !---------------------------------------------------------------------
         254           0 :        do k = 1,pver
         255           0 :           extfrc(:,k,n2d_ndx) = 1.57_r8*.6_r8*extfrc(:,k,n2p_ndx)
         256           0 :           extfrc(:,k,n_ndx) = 1.57_r8*.4_r8*extfrc(:,k,n2p_ndx)
         257             :        end do
         258             :
         259             :     endif
         260             :
         261             :     !---------------------------------------------------------------------
         262             :     !     ... set SPE NOx and HOx production
         263             :     !     Jackman et al., JGR, 2005
         264             :     !     production of 1.25 Nitrogen atoms/ion pair with branching ratios
         265             :     !     of 0.55 N(4S) and 0.7 N(2D).
         266             :     !---------------------------------------------------------------------
         267             :     !---------------------------------------------------------------------
         268             :     ! ion pairs production due to Galactic Cosmic Rays
         269             :     !---------------------------------------------------------------------
         270     1489176 :     call gcr_ionization_ionpairs( ncol, lchnk, gcr_ipr )
         271     1489176 :     call outfld( 'GCR_ionpairs', gcr_ipr, ncol, lchnk )
         272             :
         273             :     !---------------------------------------------------------------------
         274             :     ! ion pairs production due to Energetic Particle Precipitation
         275             :     !---------------------------------------------------------------------
         276     1489176 :     call epp_ionization_ionpairs( ncol, lchnk, pmid, tfld, epp_ipr )
         277     1489176 :     call outfld( 'EPP_ionpairs', epp_ipr, ncol, lchnk )
         278             :
         279  2314006344 :     epp_ipr(:ncol,:pver) = epp_ipr(:ncol,:) + gcr_ipr(:ncol,:) + mee_ap_ipr(:ncol,:)
         280             :
         281     1489176 :     if (has_dregion_ions) then
         282             :        ! D-region ion chemistry is active ...
         283             :        ! N2p production
         284           0 :        extfrc(:ncol,:pver,n2p_ndx) = extfrc(:ncol,:pver,n2p_ndx) + 0.585_r8 * epp_ipr(:ncol,:pver)
         285             :        ! O2p production
         286           0 :        extfrc(:ncol,:pver,o2p_ndx) = extfrc(:ncol,:pver,o2p_ndx) + 0.15_r8 * epp_ipr(:ncol,:pver)
         287             :        ! Np
         288           0 :        extfrc(:ncol,:pver,np_ndx)  = extfrc(:ncol,:pver,np_ndx)  + 0.185_r8 * epp_ipr(:ncol,:pver)
         289             :        ! Op
         290           0 :        extfrc(:ncol,:pver,op_ndx)  = extfrc(:ncol,:pver,op_ndx)  + 0.076_r8 * epp_ipr(:ncol,:pver)
         291             :        ! N2D/N4S branching
         292             :        ! new initial rates
         293           0 :        extfrc(:ncol,:pver,n2d_ndx) = extfrc(:ncol,:pver,n2d_ndx) + 0.583_r8 * epp_ipr(:ncol,:pver)
         294           0 :        extfrc(:ncol,:pver,n_ndx)   = extfrc(:ncol,:pver,n_ndx)   + 0.502_r8 * epp_ipr(:ncol,:pver)
         295             :        ! O
         296           0 :        extfrc(:ncol,:pver,o_ndx)   = extfrc(:ncol,:pver,o_ndx)   + 1.074_r8 * epp_ipr(:ncol,:pver)
         297             :
         298             :     else
         299             :        ! D-region ion chemistry is NOT active
         300     1489176 :        if ( n2d_ndx>0 .and. n_ndx>0 ) then
         301           0 :           extfrc(:ncol,:pver,n2d_ndx) = extfrc(:ncol,:pver,n2d_ndx) +  0.7_r8*epp_ipr(:ncol,:pver)
         302           0 :           extfrc(:ncol,:pver,  n_ndx) = extfrc(:ncol,:pver,  n_ndx) + 0.55_r8*epp_ipr(:ncol,:pver)
         303           0 :           call outfld( 'N2D_EPP', 0.7_r8*epp_ipr(:ncol,:), ncol, lchnk ) ! N(2D) produciton (molec/cm3/s)
         304           0 :           call outfld( 'N4S_EPP',0.55_r8*epp_ipr(:ncol,:), ncol, lchnk ) ! N(4S) produciton (molec/cm3/s)
         305     1489176 :        elseif ( no_ndx>0 .and. n_ndx>0 ) then
         306             :           ! for mechanisms that do not include N2D -- the EPP produce NO
         307           0 :           extfrc(:ncol,:pver, no_ndx) = extfrc(:ncol,:pver, no_ndx) +  0.7_r8*epp_ipr(:ncol,:pver)
         308           0 :           extfrc(:ncol,:pver,  n_ndx) = extfrc(:ncol,:pver,  n_ndx) + 0.55_r8*epp_ipr(:ncol,:pver)
         309           0 :           call outfld( 'NO_EPP',  0.7_r8*epp_ipr(:ncol,:), ncol, lchnk ) ! NO produciton (molec/cm3/s)
         310           0 :           call outfld( 'N4S_EPP',0.55_r8*epp_ipr(:ncol,:), ncol, lchnk ) ! N(4S) produciton (molec/cm3/s)
         311             :        endif
         312     1489176 :        if ( oh_ndx > 0 ) then
         313           0 :           do i = 1,ncol
         314           0 :              epp_hox(i,:pver) = epp_ipr(i,:pver)*hox_prod_factor( epp_ipr(i,:pver), zmid(i,:pver) )
         315             :           end do
         316           0 :           extfrc(:ncol,:pver, oh_ndx) = extfrc(:ncol,:pver, oh_ndx) + epp_hox(:ncol,:pver)
         317           0 :           call outfld( 'OH_EPP' ,  epp_hox(:ncol,:), ncol, lchnk ) ! HOX produciton (molec/cm3/s)
         318             :        endif
         319             :     endif
         320             :
         321     1489176 :     if ( has_ions ) then
         322             :        !---------------------------------------------------------------------
         323             :        !     ... set total electron production
         324             :        !---------------------------------------------------------------------
         325           0 :        do k = 1,pver
         326           0 :           extfrc(:,k,e_ndx) = extfrc(:,k,op_ndx) + extfrc(:,k,o2p_ndx) &
         327           0 :                + extfrc(:,k,np_ndx) + extfrc(:,k,n2p_ndx)
         328             :        end do
         329           0 :        call outfld( 'P_Op',  extfrc(:,:,op_ndx), ncol, lchnk )
         330           0 :        call outfld( 'P_O2p', extfrc(:,:,o2p_ndx), ncol, lchnk )
         331           0 :        call outfld( 'P_Np',  extfrc(:,:,np_ndx), ncol, lchnk )
         332           0 :        call outfld( 'P_N2p', extfrc(:,:,n2p_ndx), ncol, lchnk )
         333           0 :        call outfld( 'P_IONS',extfrc(:,:,e_ndx), ncol, lchnk )
         334             :     end if
         335             :
         336     1489176 :   end subroutine setext
         337             :
         338             : end module mo_setext

Generated by: LCOV version 1.14