Skip to content
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.

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

The final physics update PR to RRFS.v1 #211

Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
46 changes: 24 additions & 22 deletions physics/smoke_dust/module_add_emiss_burn.F90
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
grantfirl marked this conversation as resolved.
Show resolved Hide resolved
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
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
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
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
Loading