Skip to content

Commit

Permalink
Merge pull request #211 from haiqinli/production/RRFS.v1-final
Browse files Browse the repository at this point in the history 
The final physics update PR to RRFS.v1
  • Loading branch information
@dustinswales
dustinswales authored Aug 3, 2024
2 parents 4128e06 + 7a28e0b commit 5cb99bf
Show file tree
Hide file tree
Showing 4 changed files with 58 additions and 37 deletions.
46 changes: 24 additions & 22 deletions physics/smoke_dust/module_add_emiss_burn.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi,ebb_min, &
coef_bb_dc, fire_hist, hwp, hwp_prevd, &
swdown,ebb_dcycle, ebu_in, ebu,fire_type,&
q_vap, add_fire_moist_flux, &
sc_factor, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte,mpiid )
Expand All @@ -34,7 +35,6 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi,ebb_min, &
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, ebb_min ! RAR: time in seconds since start of simulation
Expand All @@ -55,12 +55,17 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi,ebb_min, &
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), INTENT(IN) :: sc_factor ! to scale up the wildfire emissions, Jordan please make this a namelist option
real(kind_phys), PARAMETER :: rinti=2.1813936e-8, ax2=3400., const2=130., &
coef2=10.6712963e-4, cx2=7200., timeq_max=3600.*24.
!>-- 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/)
! For fire diurnal cycle calculation
!real(kind_phys), parameter :: avgx1=-2.0, sigmx1=0.7, C1=0.083 ! Ag fires
!real(kind_phys), parameter :: avgx2=-0.1, sigmx2=0.8, C2=0.55 ! Grass fires, slash burns
real(kind_phys), parameter :: avgx1=0., sigmx1=2.2, C1=0.2 ! Ag fires
real(kind_phys), parameter :: avgx2=0.5, sigmx2=0.8, C2=1.1 ! Grass fires, slash burns

timeq= gmt*3600._kind_phys + real(time_int,4)
timeq= mod(timeq,timeq_max)
Expand All @@ -70,34 +75,31 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi,ebb_min, &

do j=jts,jte
do i=its,ite
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
fire_age= MAX(0.01_kind_phys,time_int/3600. + (fire_end_hr(i,j)-2.0)) !One hour delay is due to the latency of the RAVE files, hours; one more hour subtracted to have fire_end_hr in the range of 0-24 instead of 0-25

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(5) *fire_age) * &
exp(- ( log(fire_age) - avg_fire_dur(5))**2 /(2._kind_phys*sigma_fire_dur(5)**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 ))
coef_bb_dc(i,j)= C1/(sigmx1* fire_age)* exp(- (log(fire_age) - avgx1)**2 /(2.*sigmx1**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 ))
CASE (2) ! Savanna and grassland fires, or fires in the eastern US
! 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 ))
coef_bb_dc(i,j)= C2/(sigmx2* fire_age)* exp(- (log(fire_age) - avgx2)**2 /(2.*sigmx2**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

CASE (3,4) ! wildfires
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
Expand All @@ -113,15 +115,15 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi,ebb_min, &
dc_hwp= MAX(0._kind_phys,dc_hwp)
dc_hwp= MIN(20._kind_phys,dc_hwp)

! RAR: Gaussian profile for wildfires
dt1= abs(timeq - peak_hr(i,j))
dt2= timeq_max - peak_hr(i,j) + timeq ! peak hour is always <86400.
dtm= MIN(dt1,dt2)
dc_gp = rinti*( ax2 * exp(- dtm**2/(2._kind_phys*cx2**2) ) + const2 - coef2*timeq )
dc_gp = MAX(0._kind_phys,dc_gp)
! RAR: Gaussian profile for wildfires, to be used later
!dt1= abs(timeq - peak_hr(i,j))
!dt2= timeq_max - peak_hr(i,j) + timeq ! peak hour is always <86400.
!dtm= MIN(dt1,dt2)
!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
coef_bb_dc(i,j) = sc_factor* fire_hist(i,j)* dc_hwp ! RAR: scaling factor is applied to the forest fires only, except the eastern US

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 @@ -146,7 +148,7 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi,ebb_min, &
if (ebb_dcycle==1) then
conv= dtstep/(rho_phy(i,k,j)* dz8w(i,k,j))
elseif (ebb_dcycle==2) then
conv= sc_factor*coef_bb_dc(i,j)*dtstep/(rho_phy(i,k,j)* dz8w(i,k,j))
conv= coef_bb_dc(i,j)*dtstep/(rho_phy(i,k,j)* dz8w(i,k,j))
endif
dm_smoke= conv*ebu(i,k,j)
chem(i,k,j,p_smoke) = chem(i,k,j,p_smoke) + dm_smoke
Expand Down
1 change: 1 addition & 0 deletions physics/smoke_dust/rrfs_smoke_config.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,6 +43,7 @@ module rrfs_smoke_config
logical :: extended_sd_diags = .false.
real(kind_phys) :: wetdep_ls_alpha = .5 ! scavenging factor
real(kind_phys) :: plume_alpha = 0.05
real(kind_phys) :: sc_factor = 1.0

! --
integer, parameter :: CHEM_OPT_GOCART= 1
Expand Down
34 changes: 22 additions & 12 deletions physics/smoke_dust/rrfs_smoke_wrapper.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,8 @@ module rrfs_smoke_wrapper
num_moist, num_chem, num_emis_seas, num_emis_dust, &
p_qv, p_atm_shum, p_atm_cldq,plume_wind_eff, &
p_smoke, p_dust_1, p_coarse_pm, epsilc, &
n_dbg_lines, add_fire_moist_flux, plume_alpha
n_dbg_lines, add_fire_moist_flux, plume_alpha, &
sc_factor
use dust_data_mod, only : dust_alpha, dust_gamma, dust_moist_opt, &
dust_moist_correction, dust_drylimit_factor
use seas_mod, only : gocart_seasalt_driver
Expand Down Expand Up @@ -46,7 +47,8 @@ subroutine rrfs_smoke_wrapper_init( seas_opt_in,
rrfs_sd, do_plumerise_in, plumerisefire_frq_in, & ! smoke namelist
plume_wind_eff_in,add_fire_heat_flux_in, & ! smoke namelist
addsmoke_flag_in, ebb_dcycle_in, hwp_method_in, & ! smoke namelist
add_fire_moist_flux_in, plume_alpha_in, & ! smoke namelist
add_fire_moist_flux_in, & ! smoke namelist
sc_factor_in, plume_alpha_in, & ! smoke namelist
dust_opt_in, dust_alpha_in, dust_gamma_in, & ! dust namelist
dust_moist_opt_in, & ! dust namelist
dust_moist_correction_in, dust_drylimit_factor_in, & ! dust namelist
Expand All @@ -58,6 +60,7 @@ subroutine rrfs_smoke_wrapper_init( seas_opt_in,
real(kind_phys), intent(in) :: dust_alpha_in, dust_gamma_in, wetdep_ls_alpha_in, plume_alpha_in
real(kind_phys), intent(in) :: dust_moist_correction_in
real(kind_phys), intent(in) :: dust_drylimit_factor_in
real(kind_phys), intent(in) :: sc_factor_in
integer, intent(in) :: dust_opt_in,dust_moist_opt_in, wetdep_ls_opt_in, pm_settling_in, seas_opt_in
integer, intent(in) :: drydep_opt_in
logical, intent(in) :: aero_ind_fdb_in,dbg_opt_in, extended_sd_diags_in, add_fire_heat_flux_in, add_fire_moist_flux_in
Expand Down Expand Up @@ -96,6 +99,7 @@ subroutine rrfs_smoke_wrapper_init( seas_opt_in,
add_fire_heat_flux = add_fire_heat_flux_in
add_fire_moist_flux = add_fire_moist_flux_in
plume_alpha = plume_alpha_in
sc_factor = sc_factor_in
!>-Feedback
aero_ind_fdb = aero_ind_fdb_in
!>-Other
Expand Down Expand Up @@ -360,17 +364,18 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate,
! cropland, urban, cropland/natural mosaic, barren and sparsely
! vegetated and non-vegetation areas:
lu_qfire(i,j) = lu_nofire(i,j) + vegfrac(i,12,j) + vegfrac(i,13,j) + vegfrac(i,14,j) + vegfrac(i,16,j)
! Savannas and grassland fires, these fires last longer than the Ag
! fires:
lu_sfire(i,j) = lu_nofire(i,j) + vegfrac(i,8,j) + vegfrac(i,9,j) + vegfrac(i,10,j)
! Savannas and grassland fires, these fires last longer than the Ag fires:
lu_sfire(i,j) = lu_qfire(i,j) + vegfrac(i,8,j) + vegfrac(i,9,j) + vegfrac(i,10,j)
if (lu_nofire(i,j)>0.95) then ! no fires
fire_type(i,j) = 0
else if (lu_qfire(i,j)>0.9) then ! Ag. and urban fires
fire_type(i,j) = 1
else if (lu_sfire(i,j)>0.9) then ! savanna and grassland fires
fire_type(i,j) = 2
else
fire_type(i,j) = 3 ! wildfires, new approach is necessary for the controlled burns in the forest areas
else if (xlong(i,j)>260. .AND. xlat(i,j)>25. .AND. xlat(i,j)<41.) then
fire_type(i,j) = 2 ! slash burn and wildfires in the east, eastern temperate forest ecosystem
else if (lu_sfire(i,j)>0.8) then
fire_type(i,j) = 3 ! savanna and grassland fires
else
fire_type(i,j) = 4 ! potential wildfires
end if
end if
end do
Expand Down Expand Up @@ -424,7 +429,11 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate,
! Apply the diurnal cycle coefficient to frp_inst ()
do j=jts,jte
do i=its,ite
frp_inst(i,j) = MIN(frp_in(i,j)*coef_bb_dc(i,j),frp_max)
IF ( fire_type(i,j) .eq. 4 ) THEN ! only apply scaling factor to wildfires
frp_inst(i,j) = MIN(sc_factor*frp_in(i,j)*coef_bb_dc(i,j),frp_max)
ELSE
frp_inst(i,j) = MIN(frp_in(i,j)*coef_bb_dc(i,j),frp_max)
ENDIF
enddo
enddo

Expand Down Expand Up @@ -452,7 +461,8 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate,
fire_end_hr, peak_hr,curr_secs, &
coef_bb_dc,fire_hist,hwp_local,hwp_day_avg, &
swdown,ebb_dcycle,ebu_in,ebu,fire_type, &
moist(:,:,:,p_qv), add_fire_moist_flux, &
moist(:,:,:,p_qv), add_fire_moist_flux, &
sc_factor, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte , mpiid )
Expand Down Expand Up @@ -941,7 +951,7 @@ subroutine rrfs_smoke_prep( &

!---- Calculate HWP based on selected method
hwp_local = 0._kind_phys
precip_factor = 5._kind_phys + real(hour_int)*5._kind_phys/24._kind_phys
precip_factor = 2.5_kind_phys + real(hour_int, kind=kind_phys)*2.5_kind_phys/24._kind_phys
! total precip is only in the SMOKE_RRFS_DATA if ebb_dcycle == 2 and should be
! filled here before calculating HWP
! !!WARNING!! IF EBB_DYCLE != 2 and HWP_METHOD = 1 | 3, HWP will not take into account totprcp_24hrs
Expand Down
14 changes: 11 additions & 3 deletions physics/smoke_dust/rrfs_smoke_wrapper.meta
View file
Original file line number Diff line number Diff line change
Expand Up @@ -93,9 +93,9 @@
type = integer
intent = in
[hwp_method_in]
standard_name = do_smoke_forecast
long_name = index for rrfs smoke forecast
units = index
standard_name = control_for_HWP_equation
long_name = control for HWP equation
units = 1
dimensions = ()
type = integer
intent = in
Expand All @@ -106,6 +106,14 @@
dimensions = ()
type = logical
intent = in
[sc_factor_in]
standard_name = scale_factor_for_wildfire_emissions
long_name = scale factor for wildfire emissions
units = 1
dimensions = ()
type = real
kind = kind_phys
intent = in
[plume_alpha_in]
standard_name = alpha_for_plumerise_scheme
long_name = alpha paramter for plumerise scheme
Expand Down

0 comments on commit 5cb99bf

Please sign in to comment.