diff --git a/src/shr/aqm_config_mod.F90 b/src/shr/aqm_config_mod.F90 index 8344ac1..9a5b442 100644 --- a/src/shr/aqm_config_mod.F90 +++ b/src/shr/aqm_config_mod.F90 @@ -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 \ No newline at end of file diff --git a/src/shr/aqm_emis_mod.F90 b/src/shr/aqm_emis_mod.F90 index 24bda78..8cdb72a 100644 --- a/src/shr/aqm_emis_mod.F90 +++ b/src/shr/aqm_emis_mod.F90 @@ -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. @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -1737,9 +1794,8 @@ 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 @@ -1747,15 +1803,15 @@ 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) + 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 diff --git a/src/shr/aqm_fires_mod.F90 b/src/shr/aqm_fires_mod.F90 index 405d74c..e0f88d4 100644 --- a/src/shr/aqm_fires_mod.F90 +++ b/src/shr/aqm_fires_mod.F90 @@ -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 @@ -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", & @@ -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 diff --git a/src/shr/aqm_internal_mod.F90 b/src/shr/aqm_internal_mod.F90 index 3b4a110..512a986 100644 --- a/src/shr/aqm_internal_mod.F90 +++ b/src/shr/aqm_internal_mod.F90 @@ -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