Skip to content

Add guassian plume rise for fire emissions. namelist parameters for c… #120

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
KaiWang-NOAA merged 4 commits into from
Nov 18, 2025
Merged

Add guassian plume rise for fire emissions. namelist parameters for c… #120

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
7e2824a
Add guassian plume rise for fire emissions. namelist parameters for c…
bbakernoaa Sep 12, 2025
fa4ebd2
Refactor emission modules: remove unused fire fractions, add plume ri…
bbakernoaa Sep 19, 2025
7b1be04
Ensure non-negative profile values in aqm_fires_mod and clean up whit…
bbakernoaa Sep 23, 2025
14c01eb
Merge branch 'NOAA-EMC:develop' into feature/guassian_fire_plume
bbakernoaa Oct 27, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion src/shr/aqm_config_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -812,4 +812,4 @@ subroutine aqm_config_log(config, name, rc)

end subroutine aqm_config_log

end module aqm_config_mod
end module aqm_config_mod
84 changes: 70 additions & 14 deletions src/shr/aqm_emis_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -206,6 +206,8 @@ subroutine aqm_emis_src_create(config, em, rc)
em(item) % layers = 1
em(item) % scalefactor = 1.0
em(item) % topfraction = -1.0
em(item) % fires_surface_frac = 0.01
em(item) % fires_adjacent_frac = 0.15
em(item) % gridded = .true.
em(item) % sync = .false.
em(item) % verbose = .false.
Expand Down Expand Up @@ -636,6 +638,60 @@ subroutine aqm_emis_src_init(model, em, rc)
return ! bail out
end if
end if
if (trim(em % plumerise) /= "none") then
call ESMF_ConfigGetAttribute(config, em % fires_adjacent_frac, &
label=trim(em % name)//"_plume_top_fraction:", default=-1.0, rc=localrc)
if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, &
line=__LINE__, &
file=__FILE__, &
rcToReturn=rc)) &
return ! bail out
if (em % verbose) then
write(msgString,'(g20.8)') em % fires_adjacent_frac
call ESMF_LogWrite(trim(em % logprefix)//": "//pName &
//": plume_top_fraction: "//adjustl(msgString), ESMF_LOGMSG_INFO, rc=localrc)
if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, &
line=__LINE__, &
file=__FILE__, &
rcToReturn=rc)) &
return ! bail out
end if
if (em % fires_adjacent_frac > 1.0) then
call ESMF_LogSetError(ESMF_RC_NOT_VALID, &
msg="plume top fraction must not exceed 1.0", &
line=__LINE__, &
file=__FILE__, &
rcToReturn=rc)
return ! bail out
end if
end if
if (trim(em % plumerise) /= "none") then
call ESMF_ConfigGetAttribute(config, em % fires_surface_frac, &
label=trim(em % name)//"_plume_top_fraction:", default=-1.0, rc=localrc)
if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, &
line=__LINE__, &
file=__FILE__, &
rcToReturn=rc)) &
return ! bail out
if (em % verbose) then
write(msgString,'(g20.8)') em % fires_surface_frac
call ESMF_LogWrite(trim(em % logprefix)//": "//pName &
//": plume_top_fraction: "//adjustl(msgString), ESMF_LOGMSG_INFO, rc=localrc)
if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, &
line=__LINE__, &
file=__FILE__, &
rcToReturn=rc)) &
return ! bail out
end if
if (em % fires_surface_frac > 1.0) then
call ESMF_LogSetError(ESMF_RC_NOT_VALID, &
msg="plume top fraction must not exceed 1.0", &
line=__LINE__, &
file=__FILE__, &
rcToReturn=rc)
return ! bail out
end if
end if
case ("point-source")
em % gridded = .false.
! -- get plumerise type
Expand Down Expand Up @@ -1592,6 +1648,7 @@ subroutine aqm_emis_grd_read(em, spcname, buffer, localDe, rc)
character(len=ESMF_MAXSTR) :: msgString
real(ESMF_KIND_R4), pointer :: fptr(:,:)
type(aqm_state_type), pointer :: stateIn
real(ESMF_KIND_R4) :: contrib

! -- begin
if (present(rc)) rc = AQM_RC_SUCCESS
Expand Down Expand Up @@ -1634,9 +1691,8 @@ subroutine aqm_emis_grd_read(em, spcname, buffer, localDe, rc)
do i = lb(1), ub(1)
k = k + 1
if (abs(fptr(i,j)) < emAccept) then
buffer(k) = buffer(k) &
+ em % factors(item) * fptr(i,j) / stateIn % area(i,j) &
/ stateIn % area(i,j)
contrib = em % factors(item) * fptr(i,j) / real(stateIn % area(i,j), ESMF_KIND_R4) / real(stateIn % area(i,j), ESMF_KIND_R4)
buffer(k) = buffer(k) + max(0.0_ESMF_KIND_R4, contrib)
end if
end do
end do
Expand All @@ -1647,8 +1703,8 @@ subroutine aqm_emis_grd_read(em, spcname, buffer, localDe, rc)
do i = lb(1), ub(1)
k = k + 1
if (abs(fptr(i,j)) < emAccept) then
buffer(k) = buffer(k) &
+ em % factors(item) * fptr(i,j) / stateIn % area(i,j)
contrib = em % factors(item) * fptr(i,j) / real(stateIn % area(i,j), ESMF_KIND_R4)
buffer(k) = buffer(k) + max(0.0_ESMF_KIND_R4, contrib)
end if
end do
end do
Expand All @@ -1659,8 +1715,8 @@ subroutine aqm_emis_grd_read(em, spcname, buffer, localDe, rc)
do i = lb(1), ub(1)
k = k + 1
if (abs(fptr(i,j)) < emAccept) then
buffer(k) = buffer(k) &
+ em % factors(item) * fptr(i,j)
contrib = em % factors(item) * fptr(i,j)
buffer(k) = buffer(k) + max(0.0_ESMF_KIND_R4, contrib)
end if
end do
end do
Expand Down Expand Up @@ -1705,6 +1761,7 @@ subroutine aqm_emis_pts_read(em, spcname, buffer, ip, jp, ijmap, localDe, rc)
character(len=ESMF_MAXSTR) :: msgString
real(ESMF_KIND_R4) :: em_min, em_max
type(aqm_state_type), pointer :: stateIn
real(ESMF_KIND_R4) :: contrib

! -- begin
if (present(rc)) rc = AQM_RC_SUCCESS
Expand Down Expand Up @@ -1737,25 +1794,24 @@ subroutine aqm_emis_pts_read(em, spcname, buffer, ip, jp, ijmap, localDe, rc)
n = em % ijmap(m)
i = em % ip(n)
j = em % jp(n)
buffer(n) = buffer(n) &
+ em % factors(item) * em % rates(item) % values(n) / stateIn % area(i,j) &
/ stateIn % area(i,j)
contrib = em % factors(item) * em % rates(item) % values(n) / real(stateIn % area(i,j), ESMF_KIND_R4) / real(stateIn % area(i,j), ESMF_KIND_R4)
buffer(n) = buffer(n) + max(0.0_ESMF_KIND_R4, contrib)
end do
case (0)
! -- emissions are totals over each grid cell
do m = 1, size(em % ijmap)
n = em % ijmap(m)
i = em % ip(n)
j = em % jp(n)
buffer(n) = buffer(n) &
+ em % factors(item) * em % rates(item) % values(n) / stateIn % area(i,j)
contrib = em % factors(item) * em % rates(item) % values(n) / real(stateIn % area(i,j), ESMF_KIND_R4)
buffer(n) = buffer(n) + max(0.0_ESMF_KIND_R4, contrib)
end do
case (1:)
! -- emissions are already provided as surface densities, no need to normalize
do m = 1, size(em % ijmap)
n = em % ijmap(m)
buffer(n) = buffer(n) &
+ em % factors(item) * em % rates(item) % values(n)
contrib = em % factors(item) * em % rates(item) % values(n)
buffer(n) = buffer(n) + max(0.0_ESMF_KIND_R4, contrib)
end do
case default
! -- this case should never occur
Expand Down
95 changes: 83 additions & 12 deletions src/shr/aqm_fires_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,12 +21,15 @@ subroutine aqm_plume_sofiev(em, frp, profile, rc)

! -- local variables
integer :: localrc
integer :: is, ie, js, je, nl, nx, ny
integer :: c, r, l
integer :: lev0, lev1
integer :: is, ie, js, je, nl, nx, ny
real :: Hp, pblh, th0, th1, dz
real :: tfrac, w
real(AQM_KIND_R8) :: hbl
real :: fixed_surface, remaining_w
real :: dz_local, sigma, fixed_below, fixed_above, central_frac, gauss_sum, total_sum
real, dimension(:), allocatable :: gauss_weights
real(AQM_KIND_R8) :: hbl, dist
real(AQM_KIND_R8), pointer :: phi(:)
type(aqm_state_type), pointer :: state

Expand Down Expand Up @@ -55,6 +58,8 @@ subroutine aqm_plume_sofiev(em, frp, profile, rc)
nx = ie - is + 1
ny = je - js + 1

allocate(gauss_weights(nl))

! -- compute layer empirical weights
if (aqm_rc_test((em % topfraction > 1.0), &
msg="Plume top fraction must be between 0.0 and 1.0", &
Expand Down Expand Up @@ -89,19 +94,85 @@ subroutine aqm_plume_sofiev(em, frp, profile, rc)
! -- distribute linearly between surface and plume top height
lev0 = 1
lev1 = maxloc(phi, 1, mask = phi <= grav * Hp)

! Allocate fires_surface_frac of emissions to surface layer
fixed_surface = em % fires_surface_frac * w
if (fixed_surface > w) fixed_surface = w
profile(c,r,1) = fixed_surface
remaining_w = w - fixed_surface

if (lev1 > lev0) then
dz = 0.0
tfrac = 0.0
do l = lev0, lev1-1
profile(c,r,l) = w * ( phi(l) - dz ) / phi(lev1)
dz = phi(l)
tfrac = tfrac + profile(c,r,l)
end do
profile(c,r,lev1) = 1.0 - tfrac
! Gaussian distribution around plume height
phi1 = phi(lev1)
if (lev1 > 1) then
phi2 = phi(lev1 - 1)
else
profile(c,r,lev0) = 1.0
phi2 = phi(lev1)
end if
dz_local = onebg * (phi1 - phi2)
if (lev1 == 1) then
if (nl > 1) then
dz_local = onebg * (phi(2) - phi(1))
else
dz_local = onebg * phi(1)
end if
else if (lev1 == nl) then
dz_local = onebg * (phi(nl) - phi(nl-1))
end if

sigma = dz_local / 2.0

fixed_below = 0.0
fixed_above = 0.0
central_frac = remaining_w

if (lev1 > 1) then
fixed_below = em % fires_adjacent_frac * remaining_w
central_frac = central_frac - fixed_below
end if
if (lev1 < nl) then
fixed_above = em % fires_adjacent_frac * remaining_w
central_frac = central_frac - fixed_above
end if

gauss_sum = 0.0
do l = 1, nl
dist = real(l - lev1, AQM_KIND_R8)
gauss_weights(l) = exp(-0.5 * (dist / sigma)**2)
gauss_sum = gauss_sum + gauss_weights(l)
end do

if (gauss_sum > 0.0) then
do l = 1, nl
gauss_weights(l) = (gauss_weights(l) / gauss_sum) * central_frac
end do
end if

if (lev1 > 1) then
profile(c,r,lev1-1) = fixed_below
end if
if (lev1 < nl) then
profile(c,r,lev1+1) = fixed_above
end if

do l = 1, nl
profile(c,r,l) = profile(c,r,l) + gauss_weights(l)
end do

! Renormalize to ensure total sums to w
total_sum = sum(profile(c,r,1:nl))
if (abs(total_sum - w) > 1e-10) then
profile(c,r,1:nl) = profile(c,r,1:nl) * (w / total_sum)
end if

! Special case for single layer
if (nl == 1 .and. lev1 == 1) then
profile(c,r,1) = w
end if

! Ensure non-negative profile values
do l = 1, nl
profile(c,r,l) = max(0.0, profile(c,r,l))
end do

end do
end do
Expand Down
100 changes: 51 additions & 49 deletions src/shr/aqm_internal_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,64 +5,66 @@ module aqm_internal_mod
implicit none

type aqm_internal_rate_type
real(ESMF_KIND_R4), dimension(:), pointer :: values => null()
real(ESMF_KIND_R4), dimension(:), pointer :: values => null() ! Array of rate values
end type aqm_internal_rate_type

type aqm_internal_emis_type
character(len=ESMF_MAXSTR) :: name
character(len=ESMF_MAXSTR) :: type
character(len=ESMF_MAXSTR) :: format
character(len=ESMF_MAXPATHLEN) :: path
character(len=ESMF_MAXSTR) :: file
character(len=ESMF_MAXSTR) :: frequency
character(len=ESMF_MAXSTR) :: logprefix
character(len=ESMF_MAXSTR) :: plumerise
character(len=ESMF_MAXSTR) :: specfile
character(len=ESMF_MAXSTR) :: specprofile
character(len=ESMF_MAXSTR) :: latname
character(len=ESMF_MAXSTR) :: lonname
character(len=ESMF_MAXSTR) :: stkdmname
character(len=ESMF_MAXSTR) :: stkhtname
character(len=ESMF_MAXSTR) :: stktkname
character(len=ESMF_MAXSTR) :: stkvename
character(len=6) :: period
logical :: gridded
logical :: sync
logical :: verbose
real :: scalefactor
real :: topfraction
integer(ESMF_KIND_I4) :: layers
integer :: count
integer :: irec
integer :: iofmt
character(len=ESMF_MAXSTR) :: iomode
type(ESMF_GridComp) :: IO
type(ESMF_Alarm) :: alarm
character(len=ESMF_MAXSTR), dimension(:), pointer :: sources => null()
character(len=ESMF_MAXSTR), dimension(:), pointer :: species => null()
character(len=ESMF_MAXSTR), dimension(:), pointer :: units => null()
integer, dimension(:), pointer :: dens_flag => null()
integer, dimension(:), pointer :: ip => null()
integer, dimension(:), pointer :: jp => null()
integer, dimension(:), pointer :: ijmap => null()
real(ESMF_KIND_R4), dimension(:), pointer :: lat => null()
real(ESMF_KIND_R4), dimension(:), pointer :: lon => null()
real(ESMF_KIND_R4), dimension(:), pointer :: stkdm => null()
real(ESMF_KIND_R4), dimension(:), pointer :: stkht => null()
real(ESMF_KIND_R4), dimension(:), pointer :: stktk => null()
real(ESMF_KIND_R4), dimension(:), pointer :: stkve => null()
real(ESMF_KIND_R4), dimension(:), pointer :: factors => null()
type(ESMF_Field), dimension(:), pointer :: fields => null()
type(aqm_internal_rate_type), dimension(:), pointer :: rates => null()
character(len=ESMF_MAXSTR), dimension(:,:), pointer :: table => null()
character(len=ESMF_MAXSTR) :: name ! Name of the emission source or dataset
character(len=ESMF_MAXSTR) :: type ! Type of emission (e.g., area, point, mobile)
character(len=ESMF_MAXSTR) :: format ! Data format (e.g., netCDF, ASCII)
character(len=ESMF_MAXPATHLEN) :: path ! Directory path to emission files
character(len=ESMF_MAXSTR) :: file ! Base filename for emission data
character(len=ESMF_MAXSTR) :: frequency ! Update frequency (e.g., hourly, daily)
character(len=ESMF_MAXSTR) :: logprefix ! Prefix for log file names
character(len=ESMF_MAXSTR) :: plumerise ! Plume rise model configuration
character(len=ESMF_MAXSTR) :: specfile ! Species definition file
character(len=ESMF_MAXSTR) :: specprofile ! Vertical profile specification for species
character(len=ESMF_MAXSTR) :: latname ! Variable name for latitude in input files
character(len=ESMF_MAXSTR) :: lonname ! Variable name for longitude in input files
character(len=ESMF_MAXSTR) :: stkdmname ! Variable name for stack diameter
character(len=ESMF_MAXSTR) :: stkhtname ! Variable name for stack height
character(len=ESMF_MAXSTR) :: stktkname ! Variable name for stack temperature
character(len=ESMF_MAXSTR) :: stkvename ! Variable name for stack velocity
character(len=6) :: period ! Time period identifier (e.g., YYYYMM)
logical :: gridded ! Flag for gridded (true) vs point source (false) emissions
logical :: sync ! Flag to synchronize emission reads with model time
logical :: verbose ! Flag for verbose logging output
real :: fires_surface_frac ! Fraction of fire emissions at surface
real :: fires_adjacent_frac ! Fraction of fire emissions in adjacent cells
real :: scalefactor ! Scaling factor for emission values
real :: topfraction ! Fraction of emissions released at model top
integer(ESMF_KIND_I4) :: layers ! Number of vertical layers for emission distribution
integer :: count ! Count of emission sources or records
integer :: irec ! Current record index for file reading
integer :: iofmt ! I/O format flag (e.g., binary, unformatted)
character(len=ESMF_MAXSTR) :: iomode ! Mode for I/O operations (e.g., read, write)
type(ESMF_GridComp) :: IO ! ESMF Grid Component for I/O handling
type(ESMF_Alarm) :: alarm ! ESMF Alarm for scheduling emission updates
character(len=ESMF_MAXSTR), dimension(:), pointer :: sources => null() ! Array of emission source names
character(len=ESMF_MAXSTR), dimension(:), pointer :: species => null() ! Array of chemical species names
character(len=ESMF_MAXSTR), dimension(:), pointer :: units => null() ! Array of units for emission values
integer, dimension(:), pointer :: dens_flag => null() ! Flags for density-based calculations
integer, dimension(:), pointer :: ip => null() ! I-index (longitude) for point sources
integer, dimension(:), pointer :: jp => null() ! J-index (latitude) for point sources
integer, dimension(:), pointer :: ijmap => null() ! Mapping indices for grid cells
real(ESMF_KIND_R4), dimension(:), pointer :: lat => null() ! Latitude values for sources
real(ESMF_KIND_R4), dimension(:), pointer :: lon => null() ! Longitude values for sources
real(ESMF_KIND_R4), dimension(:), pointer :: stkdm => null() ! Stack diameter values
real(ESMF_KIND_R4), dimension(:), pointer :: stkht => null() ! Stack height values
real(ESMF_KIND_R4), dimension(:), pointer :: stktk => null() ! Stack temperature values
real(ESMF_KIND_R4), dimension(:), pointer :: stkve => null() ! Stack velocity values
real(ESMF_KIND_R4), dimension(:), pointer :: factors => null() ! Scaling or conversion factors
type(ESMF_Field), dimension(:), pointer :: fields => null() ! Array of ESMF Fields for emissions
type(aqm_internal_rate_type), dimension(:), pointer :: rates => null() ! Array of emission rate data structures
character(len=ESMF_MAXSTR), dimension(:,:), pointer :: table => null() ! Lookup table for emission parameters
end type

type aqm_internal_data_type
type(aqm_internal_emis_type), pointer :: emis(:) => null()
type(aqm_internal_emis_type), pointer :: emis(:) => null() ! Array of emission configurations
end type

type aqm_internal_state_type
type(aqm_internal_data_type), pointer :: wrap => null()
type(aqm_internal_data_type), pointer :: wrap => null() ! Wrapper for internal data structures
end type

private
Expand Down