Skip to content

Commit

Permalink
"smoke updates for RRFS.v1"
Browse files Browse the repository at this point in the history
  • Loading branch information
@haiqinli @grantfirl
haiqinli authored and grantfirl committed Oct 29, 2024
1 parent 228c550 commit 728ecf1
Show file tree
Hide file tree
Showing 6 changed files with 370 additions and 435 deletions.
97 changes: 50 additions & 47 deletions physics/smoke_dust/module_add_emiss_burn.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,11 +6,12 @@ module module_add_emiss_burn
use machine , only : kind_phys
use rrfs_smoke_config
CONTAINS
subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi, &
subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi,ebb_min, &
chem,julday,gmt,xlat,xlong, &
fire_end_hr, peak_hr,time_int, &
coef_bb_dc, fire_hist, hwp, hwp_prevd, &
swdown,ebb_dcycle, ebu_in, ebu,fire_type,&
q_vap, add_fire_moist_flux, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte,mpiid )
Expand All @@ -27,102 +28,99 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi, &
INTENT(INOUT ) :: chem ! shall we set num_chem=1 here?

real(kind_phys), DIMENSION( ims:ime, kms:kme, jms:jme ), &
INTENT(INOUT ) :: ebu
INTENT(INOUT ) :: ebu, q_vap ! SRB: added q_vap

real(kind_phys), DIMENSION(ims:ime,jms:jme), INTENT(IN) :: xlat,xlong, swdown
real(kind_phys), DIMENSION(ims:ime,jms:jme), INTENT(IN) :: hwp, peak_hr, fire_end_hr, ebu_in !RAR: Shall we make fire_end integer?
real(kind_phys), DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: coef_bb_dc ! RAR:
real(kind_phys), DIMENSION(ims:ime,jms:jme), INTENT(IN) :: hwp_prevd
real(kind_phys), DIMENSION(ims:ime,jms:jme), INTENT(IN) :: xlat,xlong, swdown
real(kind_phys), DIMENSION(ims:ime,jms:jme), INTENT(IN) :: hwp, peak_hr, fire_end_hr, ebu_in !RAR: Shall we make fire_end integer?
real(kind_phys), DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: coef_bb_dc ! RAR:
real(kind_phys), DIMENSION(ims:ime,jms:jme), INTENT(IN) :: hwp_prevd

real(kind_phys), DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(IN) :: dz8w,rho_phy !,rel_hum
real(kind_phys), INTENT(IN) :: dtstep, gmt
real(kind_phys), INTENT(IN) :: time_int,pi ! RAR: time in seconds since start of simulation
real(kind_phys), INTENT(IN) :: time_int, pi, ebb_min ! RAR: time in seconds since start of simulation
INTEGER, DIMENSION(ims:ime,jms:jme), INTENT(IN) :: fire_type
integer, INTENT(IN) :: ebb_dcycle ! RAR: this is going to be namelist dependent, ebb_dcycle=means
real(kind_phys), DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: fire_hist
!>--local
logical, intent(in) :: add_fire_moist_flux
integer :: i,j,k,n,m
integer :: icall=0
integer :: icall !=0
real(kind_phys) :: conv_rho, conv, dm_smoke, dc_hwp, dc_gp, dc_fn !daero_num_wfa, daero_num_ifa !, lu_sum1_5, lu_sum12_14

INTEGER, PARAMETER :: kfire_max=51 ! max vertical level for BB plume rise

real(kind_phys) :: timeq, fire_age, age_hr, dt1,dt2,dtm, coef_con ! For BB emis. diurnal cycle calculation
real(kind_phys), PARAMETER :: ef_h2o=324.22 ! Emission factor for water vapor
! Constants for the fire diurnal cycle calculation ! JLS - needs to be
! defined below due to intent(in) of pi
real(kind_phys) :: coef_con
real(kind_phys) :: timeq, fire_age, age_hr, dt1,dt2,dtm ! For BB emis. diurnal cycle calculation

! For Gaussian diurnal cycle
real(kind_phys), PARAMETER :: sc_factor=1. ! to scale up the wildfire emissions, TBD later
real(kind_phys), PARAMETER :: rinti=2.1813936e-8, ax2=3400., const2=130., &
coef2=10.6712963e-4, cx2=7200., timeq_max=3600.*24.
!>-- Fire parameters
!>-- Fire parameters: Fores west, Forest east, Shrubland, Savannas, Grassland, Cropland
real(kind_phys), dimension(1:5), parameter :: avg_fire_dur = (/8.9, 4.2, 3.3, 3.0, 1.4/)
real(kind_phys), dimension(1:5), parameter :: sigma_fire_dur = (/8.7, 6.0, 5.5, 5.2, 2.4/)

timeq= gmt*3600._kind_phys + real(time_int,4)
timeq= mod(timeq,timeq_max)
coef_con=1._kind_phys/((2._kind_phys*pi)**0.5)


! RAR: Grasslands (29% of ther western HRRR CONUS domain) probably also need to
! be added below, check this later
! RAR: In the HRRR CONUS domain (western part) crop 11%, 2% cropland/natural
! vegetation and 0.4% urban of pixels
!.OR. lu_index(i,j)==14) then ! Croplands(12), Urban and Built-Up(13),
!cropland/natural vegetation (14) mosaic in MODI-RUC vegetation classes
! Peak hours for the fire activity depending on the latitude
! if (xlong(i,j)<-130.) then max_ti= 24.041288* 3600. !
! peak at 24 UTC, fires in Alaska
! elseif (xlong(i,j)<-100.) then max_ti= 22.041288* 3600.
! ! peak at 22 UTC, fires in the western US
! elseif (xlong(i,j)<-70.) then ! peak at 20 UTC, fires in
! the eastern US, max_ti= 20.041288* 3600.
! else max_ti= 18.041288* 3600.
! endif
! RAR: for option #1 ebb and frp are ingested for 24 hours. No modification is
! applied!
! RAR: for option #1 ebb and frp are ingested for 24 hours. No modification is applied!
if (ebb_dcycle==1) then
do k=kts,kte
do i=its,ite
ebu(i,k,1)=ebu_in(i,1) ! RAR:
ebu(i,k,1)=ebu_in(i,1)
enddo
enddo
endif

if (ebb_dcycle==2) then

! Constants for the fire diurnal cycle calculation
coef_con=1._kind_phys/((2._kind_phys*pi)**0.5_kind_phys)

do j=jts,jte
do i=its,ite
fire_age= time_int + (fire_end_hr(i,j)-1._kind_phys)*3600._kind_phys !One hour delay is due to the latency of the RAVE files
fire_age= MAX(0._kind_phys,fire_age)
fire_age= time_int/3600._kind_phys + (fire_end_hr(i,j)-1._kind_phys) !One hour delay is due to the latency of the RAVE files
fire_age= MAX(0.1_kind_phys,fire_age) ! in hours

SELECT CASE ( fire_type(i,j) ) !Ag, urban fires, bare land etc.
CASE (1)
! these fires will have exponentially decreasing diurnal cycle,
coef_bb_dc(i,j) = coef_con*1._kind_phys/(sigma_fire_dur(1) *fire_age) * &
exp(- ( log(fire_age) - avg_fire_dur(1))**2 /(2._kind_phys*sigma_fire_dur(1)**2 ))
coef_bb_dc(i,j) = coef_con*1._kind_phys/(sigma_fire_dur(5) *fire_age) * &
exp(- ( log(fire_age) - avg_fire_dur(5))**2 /(2._kind_phys*sigma_fire_dur(5)**2 ))

IF ( dbg_opt .AND. time_int<5000.) then
WRITE(6,*) 'i,j,peak_hr(i,j) ',i,j,peak_hr(i,j)
WRITE(6,*) 'coef_bb_dc(i,j) ',coef_bb_dc(i,j)
END IF

CASE (2) ! Savanna and grassland fires
coef_bb_dc(i,j) = coef_con*1._kind_phys/(sigma_fire_dur(4) *fire_age) * &
exp(- ( log(fire_age) - avg_fire_dur(4))**2 /(2._kind_phys*sigma_fire_dur(4)**2 ))

IF ( dbg_opt .AND. time_int<5000.) then
WRITE(6,*) 'i,j,peak_hr(i,j) ',i,j,peak_hr(i,j)
WRITE(6,*) 'coef_bb_dc(i,j) ',coef_bb_dc(i,j)
END IF



CASE (3)
age_hr= fire_age/3600._kind_phys
!age_hr= fire_age/3600._kind_phys

IF (swdown(i,j)<.1 .AND. age_hr> 12. .AND. fire_hist(i,j)>0.75) THEN
IF (swdown(i,j)<.1 .AND. fire_age> 12. .AND. fire_hist(i,j)>0.75) THEN
fire_hist(i,j)= 0.75_kind_phys
ENDIF
IF (swdown(i,j)<.1 .AND. age_hr> 24. .AND. fire_hist(i,j)>0.5) THEN
IF (swdown(i,j)<.1 .AND. fire_age> 24. .AND. fire_hist(i,j)>0.5) THEN
fire_hist(i,j)= 0.5_kind_phys
ENDIF
IF (swdown(i,j)<.1 .AND. age_hr> 48. .AND. fire_hist(i,j)>0.25) THEN
IF (swdown(i,j)<.1 .AND. fire_age> 48. .AND. fire_hist(i,j)>0.25) THEN
fire_hist(i,j)= 0.25_kind_phys
ENDIF

! this is based on hwp, hourly or instantenous TBD
dc_hwp= hwp(i,j)/ MAX(5._kind_phys,hwp_prevd(i,j))
dc_hwp= hwp(i,j)/ MAX(10._kind_phys,hwp_prevd(i,j))
dc_hwp= MAX(0._kind_phys,dc_hwp)
dc_hwp= MIN(25._kind_phys,dc_hwp)
dc_hwp= MIN(20._kind_phys,dc_hwp)

! RAR: Gaussian profile for wildfires
dt1= abs(timeq - peak_hr(i,j))
Expand All @@ -131,8 +129,8 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi, &
dc_gp = rinti*( ax2 * exp(- dtm**2/(2._kind_phys*cx2**2) ) + const2 - coef2*timeq )
dc_gp = MAX(0._kind_phys,dc_gp)

dc_fn = MIN(dc_hwp/dc_gp,3._kind_phys)
coef_bb_dc(i,j) = fire_hist(i,j)* dc_hwp
!dc_fn = MIN(dc_hwp/dc_gp,3._kind_phys)
coef_bb_dc(i,j) = fire_hist(i,j)* dc_hwp

IF ( dbg_opt .AND. time_int<5000.) then
WRITE(6,*) 'i,j,fire_hist(i,j),peak_hr(i,j) ', i,j,fire_hist(i,j),peak_hr(i,j)
Expand All @@ -152,7 +150,7 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi, &
do j=jts,jte
do i=its,ite
do k=kts,kfire_max
if (ebu(i,k,j)<0.001_kind_phys) cycle
if (ebu(i,k,j)<ebb_min) cycle

if (ebb_dcycle==1) then
conv= dtstep/(rho_phy(i,k,j)* dz8w(i,k,j))
Expand All @@ -163,6 +161,12 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi, &
chem(i,k,j,p_smoke) = chem(i,k,j,p_smoke) + dm_smoke
chem(i,k,j,p_smoke) = MIN(MAX(chem(i,k,j,p_smoke),0._kind_phys),5.e+3_kind_phys)

! SRB: Modifying Water Vapor content based on Emissions
if (add_fire_moist_flux) then
q_vap(i,k,j) = q_vap(i,k,j) + (dm_smoke * ef_h2o * 1.e-9) ! kg/kg:used 1.e-9 as dm_smoke is in ug/kg
q_vap(i,k,j) = MIN(MAX(q_vap(i,k,j),0._kind_phys),1.e+3_kind_phys)
endif

if ( dbg_opt .and. (k==kts .OR. k==kfire_max) .and. (icall .le. n_dbg_lines) ) then
WRITE(1000+mpiid,*) 'add_emiss_burn:xlat,xlong,curr_secs,fire_type,fire_hist,peak_hr', xlat(i,j),xlong(i,j),int(time_int),fire_type(i,j),fire_hist(i,j),peak_hr(i,j)
WRITE(1000+mpiid,*) 'add_emiss_burn:xlat,xlong,curr_secs,coef_bb_dc,ebu',xlat(i,j),xlong(i,j),int(time_int),coef_bb_dc(i,j),ebu(i,k,j)
Expand All @@ -172,7 +176,6 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi, &
enddo
enddo


END subroutine add_emis_burn

END module module_add_emiss_burn
Expand Down
Loading

0 comments on commit 728ecf1

Please sign in to comment.