LCOV - code coverage report
Current view: | top level - chemistry/mozart - fire_emissions.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 11 | 96 | 11.5 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 3 | 3 | 100.0 %

              Line data    Source code
           1             : !================================================================================
           2             : ! manages mapping of CLM generated wild fire emissions to chemical constituents
           3             : !================================================================================
           4             : module fire_emissions
           5             :
           6             :   use shr_kind_mod,      only : r8 => shr_kind_r8
           7             :   use shr_fire_emis_mod, only : shr_fire_emis_mechcomps, shr_fire_emis_mechcomps_n, shr_fire_emis_elevated
           8             :   use srf_field_check,   only : active_Fall_flxfire
           9             :   use shr_const_mod,     only : pi => SHR_CONST_PI
          10             :   use shr_const_mod,     only : avogad => SHR_CONST_AVOGAD ! Avogadro's number ~ molecules/kmole
          11             :   use cam_abortutils,    only : endrun
          12             :   use cam_history,       only : addfld, horiz_only, outfld, fieldname_len
          13             :   use cam_logfile,       only : iulog
          14             :   use ppgrid,            only : pver, pverp
          15             :   use constituents,      only : cnst_get_ind
          16             :   use rad_constituents,  only : rad_cnst_get_aer_props, rad_cnst_num_name
          17             :   use mo_chem_utls,      only : get_spc_ndx, get_extfrc_ndx
          18             :   use chem_mods,         only : adv_mass ! g/mole
          19             :   use infnan,            only : nan, assignment(=)
          20             :
          21             :   implicit none
          22             :   private
          23             :   save
          24             :
          25             :   public :: fire_emissions_init
          26             :   public :: fire_emissions_srf
          27             :   public :: fire_emissions_vrt
          28             :
          29             :   ! for surface emissions
          30             :   integer, allocatable :: fire_emis_indices_map(:)
          31             :
          32             :   ! for vertically distributed forcings
          33             :   integer,  allocatable :: frc_spc_map(:)
          34             :   integer,  allocatable :: chm_spc_map(:)
          35             :   integer,  allocatable :: frc_num_map(:)
          36             :   real(r8), allocatable :: spc_mass_factor(:)
          37             :   real(r8), allocatable :: num_mass_factor(:)
          38             :   character(len=fieldname_len), allocatable :: fire_frc_name(:)
          39             :   character(len=fieldname_len), allocatable :: fire_numfrc_name(:)
          40             :   character(len=fieldname_len), allocatable :: fire_sflx_name(:)
          41             :   character(len=fieldname_len), allocatable :: fire_vflx_name(:)
          42             :
          43             : !================================================================================
          44             : contains
          45             : !================================================================================
          46             :
          47             :   !------------------------------------------------------------------------------
          48             :   !------------------------------------------------------------------------------
          49        1536 :   subroutine fire_emissions_init()
          50             :
          51             :     ! local vars
          52             :     integer :: n, ii
          53             :
          54             :     integer :: frc_ndx, spc_ndx
          55             :     integer :: mode, spec
          56             :     character(len=16) :: name
          57             :     character(len=32) :: num_name
          58             :
          59             :     real(r8), parameter :: demis_acc = 0.134e-6_r8 ! meters
          60             :     ! volume-mean emissions diameter of primary BC/OM aerosols, see :
          61             :     ! Liu et al, Toward a minimal representation of aerosols in climate models:
          62             :     ! Description and evaluation in the Community Atmosphere Model CAM5.
          63             :     ! Geosci. Model Dev., 5, 709–739, doi:10.5194/gmd-5-709-2012
          64             :     ! and Table S1 in Supplement: http://www.geosci-model-dev.net/5/709/2012/gmd-5-709-2012-supplement.pdf
          65             :
          66             :     real(r8), parameter :: x_numfact = 1.e-6_r8 * avogad * 6.0_r8 / (pi*(demis_acc**3))   ! 1.e-6 converts m-3 to cm-3.
          67             :     real(r8) :: specdens  ! kg/m3
          68             :     logical :: found
          69             :     character(len=12) :: units
          70             :
          71        1536 :     if (shr_fire_emis_mechcomps_n<1) return
          72             :
          73           0 :     if (shr_fire_emis_elevated) then ! initialize elevated forcings
          74             :
          75           0 :        allocate( frc_spc_map(shr_fire_emis_mechcomps_n) )
          76           0 :        allocate( chm_spc_map(shr_fire_emis_mechcomps_n) )
          77           0 :        allocate( frc_num_map(shr_fire_emis_mechcomps_n) )
          78           0 :        allocate( spc_mass_factor(shr_fire_emis_mechcomps_n) )
          79           0 :        allocate( num_mass_factor(shr_fire_emis_mechcomps_n) )
          80           0 :        allocate( fire_frc_name(shr_fire_emis_mechcomps_n) )
          81           0 :        allocate( fire_numfrc_name(shr_fire_emis_mechcomps_n) )
          82           0 :        allocate( fire_sflx_name(shr_fire_emis_mechcomps_n) )
          83           0 :        allocate( fire_vflx_name(shr_fire_emis_mechcomps_n) )
          84             :
          85           0 :        frc_spc_map(:) = -1
          86           0 :        frc_num_map(:) = -1
          87           0 :        spc_mass_factor(:) = nan
          88           0 :        num_mass_factor(:) = nan
          89             :
          90           0 :        call addfld ('Fire_ZTOP', horiz_only, 'A', 'm', 'top of vertical fire emissions' )
          91             :
          92           0 :        do n=1,shr_fire_emis_mechcomps_n
          93             :
          94           0 :           name = shr_fire_emis_mechcomps(n)%name
          95           0 :           fire_frc_name(n)  = 'FireFrc_'//trim(name)
          96           0 :           fire_sflx_name(n) = 'FireSFLX_'//trim(name)
          97           0 :           fire_vflx_name(n) = 'FireVFLX_'//trim(name)
          98             :
          99           0 :           call addfld (fire_frc_name(n), (/'lev'/), 'A','molecules/cm^3/s', 'vertical fire emissions for '//trim(name))
         100           0 :           call addfld (fire_sflx_name(n),horiz_only,'A','kg/m^2/s', 'surface fire emissions for '//trim(name))
         101           0 :           call addfld (fire_vflx_name(n),horiz_only,'A','kg/m^2/s', 'vertically integrated fire emissions for '//trim(name))
         102             :
         103           0 :           frc_ndx = get_extfrc_ndx( name )
         104           0 :           spc_ndx = get_spc_ndx( name )
         105             :
         106           0 :           if (frc_ndx>0 .and. spc_ndx>0) then
         107           0 :              frc_spc_map(n) = frc_ndx
         108           0 :              chm_spc_map(n) = spc_ndx
         109             :           else
         110           0 :              write(iulog,*) 'fire_emissions_init: not able to map '//trim(name)//' to chem species/forcing '
         111           0 :              write(iulog,*) 'fire_emissions_init:  ... frc_ndx = ',frc_ndx
         112           0 :              write(iulog,*) 'fire_emissions_init:  ... spc_ndx = ',spc_ndx
         113           0 :              call endrun('fire_emissions_init: not able to map '//trim(name)//' to chem species/forcing ')
         114             :           endif
         115             :
         116           0 :           spc_mass_factor(n) = 1.e-6_r8 * avogad / adv_mass(spc_ndx) ! 1.e-6 converts m-3 to cm-3.
         117             :           ! (molecules/kmole) / (g/mole) --> molecules/kg
         118             :
         119             :           ! for MAM need to include cooresponding forcings of number densities
         120             :
         121           0 :           found = rad_cnst_num_name(0, name, num_name, mode_out=mode, spec_out=spec )
         122             :
         123           0 :           if ( found ) then
         124             :
         125           0 :              frc_ndx = get_extfrc_ndx( num_name )
         126             :
         127           0 :              call rad_cnst_get_aer_props(0, mode, spec, density_aer=specdens)
         128           0 :              frc_num_map(n) = frc_ndx
         129           0 :              num_mass_factor(n) = x_numfact / specdens
         130             :
         131           0 :              fire_numfrc_name(n) = 'FireFrc_'//trim(name)//'_'//trim(num_name)
         132           0 :              call addfld (fire_numfrc_name(n),(/'lev'/), 'A', 'molecules/cm^3/s', &
         133           0 :                   'vertical fire emissions for '//trim(num_name)//' due to component '//trim(name))
         134             :
         135             :           endif
         136             :
         137             :        enddo
         138             :
         139             :     else ! initialize surface emissions
         140             :
         141           0 :        allocate( fire_emis_indices_map(shr_fire_emis_mechcomps_n) )
         142             :
         143           0 :        do n=1,shr_fire_emis_mechcomps_n
         144           0 :           call cnst_get_ind (shr_fire_emis_mechcomps(n)%name,  fire_emis_indices_map(n), abort=.false. )
         145             :
         146           0 :           ii = get_spc_ndx(shr_fire_emis_mechcomps(n)%name)
         147           0 :           if (ii<1) then
         148             :              call endrun('gas_phase_chemdr_inti: Fire emissions compound not in chemistry mechanism : '&
         149           0 :                   //trim(shr_fire_emis_mechcomps(n)%name))
         150             :           endif
         151             :
         152           0 :           if (index(shr_fire_emis_mechcomps(n)%name,'num_')>0) then
         153           0 :              units = '1/m2/sec'
         154             :           else
         155           0 :              units = 'kg/m2/sec'
         156             :           endif
         157             :
         158             :           ! Fire emis history fields
         159             :           call addfld( 'FireSF_'//trim(shr_fire_emis_mechcomps(n)%name),horiz_only,'A',units,&
         160           0 :                trim(shr_fire_emis_mechcomps(n)%name)//' Fire emissions flux')
         161             :
         162             :        enddo
         163             :
         164             :     endif
         165             :
         166             :   end subroutine fire_emissions_init
         167             :
         168             :   !------------------------------------------------------------------------------
         169             :   ! sets surface emissions
         170             :   !------------------------------------------------------------------------------
         171     1489176 :   subroutine fire_emissions_srf( lchnk, ncol, fireflx, sflx )
         172             :
         173             :     ! dummy args
         174             :     integer,  intent(in) :: lchnk, ncol
         175             :     real(r8), pointer, intent(in) :: fireflx(:,:)
         176             :     real(r8), intent(inout) :: sflx(:,:)
         177             :
         178             :     ! local vars
         179             :     integer :: i, n
         180             :
         181             :     ! fire surface emissions if not elevated forcing
         182     1489176 :     if ((.not.shr_fire_emis_elevated) .and. active_Fall_flxfire .and. shr_fire_emis_mechcomps_n>0 ) then
         183             :
         184             :        ! set Fire Emis fluxes ( add to other emis )
         185           0 :        do i =1,ncol
         186           0 :           do n = 1,shr_fire_emis_mechcomps_n
         187           0 :              sflx(i,fire_emis_indices_map(n)) &
         188           0 :                   = sflx(i,fire_emis_indices_map(n)) + fireflx(i,n)
         189             :           enddo
         190             :        end do
         191             :
         192             :        ! output fire emis fluxes to history
         193           0 :        do n = 1,shr_fire_emis_mechcomps_n
         194           0 :           call outfld('FireSF_'//trim(shr_fire_emis_mechcomps(n)%name), fireflx(:ncol,n), ncol, lchnk)
         195             :        enddo
         196             :
         197             :     endif
         198             :
         199     1489176 :   end subroutine fire_emissions_srf
         200             :
         201             :   !------------------------------------------------------------------------------
         202             :   ! sets vertical emissions (forcings)
         203             :   ! vertically distributes wild fire emissions
         204             :   !------------------------------------------------------------------------------
         205     1489176 :   subroutine fire_emissions_vrt( ncol, lchnk, zint, fire_sflx, fire_ztop, frcing )
         206             :
         207             :    ! args
         208             :     integer,          intent(in) :: ncol,lchnk
         209             :     real(r8),         intent(in) :: zint(:,:)      ! interface geopot above surface (km)
         210             :     real(r8),pointer, intent(in) :: fire_sflx(:,:) ! fire surface emissions (kg/m2/sec)
         211             :     real(r8),pointer, intent(in) :: fire_ztop(:)   ! top of vert distribution of fire surface emissions (m)
         212             :     real(r8),      intent(inout) :: frcing(:,:,:)  ! insitu forcings (molecules/cm3/sec)
         213             :
         214             :     ! local vars
         215     2978352 :     real(r8) :: vertical_fire(ncol,pver), ztop
         216             :     integer  :: n, i,k
         217     2978352 :     real(r8) :: fire_frc(ncol,pver)
         218     1489176 :     real(r8) :: sflx(ncol)
         219             :
         220     1489176 :     if (.not.shr_fire_emis_elevated) return
         221     1489176 :     if (shr_fire_emis_mechcomps_n<1) return
         222             :
         223             :     !   define vertical_fire from Dentener  units /m
         224           0 :     do k=1,pver
         225           0 :        do i=1,ncol
         226           0 :           ztop = fire_ztop(i)*1.e-3_r8 ! convert m to km
         227           0 :           if(zint(i,k)<ztop)then
         228           0 :              vertical_fire(i,k)=1.e-3_r8/(ztop-zint(i,pverp))
         229           0 :           elseif(zint(i,k)>ztop.and.zint(i,k+1)<ztop)then
         230           0 :              vertical_fire(i,k)=1.e-3_r8*(ztop-zint(i,k+1))/(zint(i,k)-zint(i,k+1))/(ztop-zint(i,pverp))
         231             :           else
         232           0 :              vertical_fire(i,k)=0._r8
         233             :           endif
         234             :        end do
         235             :     end do
         236             :
         237             :     ! extend the surface emissions in the vertical ....
         238             :
         239           0 :     do n=1,shr_fire_emis_mechcomps_n
         240             :
         241           0 :        fire_frc(:,:) = 0._r8
         242             :
         243           0 :        do k = 1,pver !        (kg/m2/sec)        (/m)                    (molecules/kg)*1.e-6   --> molecules/cm3/s
         244           0 :           fire_frc(:ncol,k) = fire_sflx(:ncol,n) * vertical_fire(:ncol,k) * spc_mass_factor(n)  ! molecules/cm3/s
         245             :        enddo
         246           0 :        call outfld( fire_frc_name(n), fire_frc, ncol, lchnk )
         247           0 :        frcing(:ncol,:,frc_spc_map(n)) = frcing(:ncol,:,frc_spc_map(n)) + fire_frc(:ncol,:)
         248             :
         249             :        ! for debugging ...
         250             :        ! vertical intergration of the forcing should get back the surface flux
         251           0 :        sflx(:) = 0._r8
         252           0 :        do k = 1,pver
         253           0 :           sflx(:ncol) = sflx(:ncol) + 1.e5_r8*(zint(:ncol,k)-zint(:ncol,k+1))*fire_frc(:ncol,k) ! molecules/cm3/s --> molecules/cm2/sec
         254             :        enddo
         255           0 :        sflx(:ncol) = sflx(:ncol) * 1.e4_r8 * adv_mass(chm_spc_map(n))/avogad ! molecules/cm2/sec --> kg/m2/sec
         256             :                                                                              ! / avogad     --> kmoles/cm2/sec
         257             :                                                                              ! * adv_mass   -->     kg/cm2/sec
         258             :                                                                              ! * 1.e4       -->     kg/ m2/sec
         259           0 :        call outfld( fire_vflx_name(n),      sflx(:ncol  ), ncol, lchnk )
         260           0 :        call outfld( fire_sflx_name(n), fire_sflx(:ncol,n), ncol, lchnk )
         261             :
         262             :        ! for MAM need to include corresponding forcings of number densities
         263           0 :        if (frc_num_map(n)>0) then
         264           0 :           do k = 1,pver
         265           0 :              fire_frc(:ncol,k) = fire_sflx(:ncol,n) * vertical_fire(:ncol,k) * num_mass_factor(n)  ! molecules/cm3/s
         266             :           enddo
         267           0 :           call outfld( fire_numfrc_name(n), fire_frc, ncol, lchnk )
         268           0 :           frcing(:ncol,:,frc_num_map(n)) = frcing(:ncol,:,frc_num_map(n)) + fire_frc(:ncol,:)
         269             :        endif
         270             :
         271             :     enddo
         272             :
         273           0 :     call outfld( 'Fire_ZTOP', fire_ztop(:ncol), ncol, lchnk )
         274             :
         275           0 :   end subroutine fire_emissions_vrt
         276             :
         277             : end module fire_emissions

Generated by: LCOV version 1.14