Skip to content

Commit

Permalink
Merge branch 'ufs/dev' into rrfsv1-to-ufs/dev5
Browse files Browse the repository at this point in the history
grantfirl committed Dec 6, 2024
2 parents c69a7c3 + dab57fc commit 5d7c8f1
Showing 9 changed files with 395 additions and 437 deletions.
7 changes: 2 additions & 5 deletions physics/Interstitials/UFS_SCM_NEPTUNE/sgscloud_radpre.F90
Original file line number Diff line number Diff line change
@@ -55,7 +55,7 @@ subroutine sgscloud_radpre_run( &
nlay, plyr, xlat, dz,de_lgth, &
cldsa,mtopa,mbota, &
imp_physics, imp_physics_gfdl,&
imp_physics_fa, &
imp_physics_fa, conv_cf_opt, &
iovr, &
errmsg, errflg )

@@ -75,7 +75,7 @@ subroutine sgscloud_radpre_run( &
real(kind=kind_phys) :: gfac
integer, intent(in) :: im, levs, imfdeepcnv, imfdeepcnv_gf, &
& nlay, imfdeepcnv_sas, imfdeepcnv_c3, imp_physics, &
& imp_physics_gfdl, imp_physics_fa
& imp_physics_gfdl, imp_physics_fa, conv_cf_opt
logical, intent(in) :: flag_init, flag_restart, do_mynnedmf

real(kind=kind_phys), dimension(:,:), intent(inout) :: qc, qi
@@ -120,9 +120,6 @@ subroutine sgscloud_radpre_run( &
real :: a, f, sigq, qmq, qt, xl, th, thl, rsl, cpm, cb_cf
real(kind=kind_phys) :: tlk

!Option to convective cloud fraction
integer, parameter :: conv_cf_opt = 0 !0: C-B, 1: X-R

! Initialize CCPP error handling variables
errmsg = ''
errflg = 0
7 changes: 7 additions & 0 deletions physics/Interstitials/UFS_SCM_NEPTUNE/sgscloud_radpre.meta
Original file line number Diff line number Diff line change
@@ -275,6 +275,13 @@
dimensions = ()
type = integer
intent = in
[conv_cf_opt]
standard_name = option_for_convection_scheme_cloud_fraction_computation
long_name = option for convection scheme cloud fraction computation
units = flag
dimensions = ()
type = integer
intent = in
[qc_save]
standard_name = cloud_condensed_water_mixing_ratio_save
long_name = ratio of mass of cloud water to mass of dry air plus vapor (without condensates) before entering a physics scheme
38 changes: 19 additions & 19 deletions physics/PBL/MYNN_EDMF/module_bl_mynn.F90
Original file line number Diff line number Diff line change
@@ -302,7 +302,7 @@ MODULE module_bl_mynn
! Note that the following mixing-length constants are now specified in mym_length
! &cns=3.5, alp1=0.23, alp2=0.3, alp3=3.0, alp4=10.0, alp5=0.2

real(kind_phys), parameter :: gpw=5./3., qcgmin=1.e-8, qkemin=1.e-12
real(kind_phys), parameter :: qkemin=1.e-4
real(kind_phys), parameter :: tliq = 269. !all hydrometeors are liquid when T > tliq

! Constants for cloud PDF (mym_condensation)
@@ -1937,11 +1937,11 @@ SUBROUTINE mym_length ( &
h1=MIN(h1,maxdz) ! 1/2 transition layer depth
h2=h1/2.0 ! 1/4 transition layer depth

qkw(kts) = SQRT(MAX(qke(kts),1.0e-10))
qkw(kts) = SQRT(MAX(qke(kts), qkemin))
DO k = kts+1,kte
afk = dz(k)/( dz(k)+dz(k-1) )
abk = 1.0 -afk
qkw(k) = SQRT(MAX(qke(k)*abk+qke(k-1)*afk,1.0e-3))
qkw(k) = SQRT(MAX(qke(k)*abk+qke(k-1)*afk, qkemin))
END DO

elt = 1.0e-5
@@ -1961,7 +1961,7 @@ SUBROUTINE mym_length ( &

elt = alp1*elt/vsc
vflx = ( vt(kts)+1.0 )*flt +( vq(kts)+tv0 )*flq
vsc = ( gtr*elt*MAX( vflx, 0.0 ) )**(1.0/3.0)
vsc = ( gtr*elt*MAX( vflx, 0.0 ) )**onethird

! ** Strictly, el(i,k=1) is not zero. **
el(kts) = 0.0
@@ -2019,14 +2019,14 @@ SUBROUTINE mym_length ( &
h1=MIN(h1,600.) ! 1/2 transition layer depth
h2=h1/2.0 ! 1/4 transition layer depth

qtke(kts)=MAX(0.5*qke(kts), 0.01) !tke at full sigma levels
qtke(kts)=MAX(0.5*qke(kts), 0.5*qkemin) !tke at full sigma levels
thetaw(kts)=theta(kts) !theta at full-sigma levels
qkw(kts) = SQRT(MAX(qke(kts),1.0e-10))
qkw(kts) = SQRT(MAX(qke(kts), qkemin))

DO k = kts+1,kte
afk = dz(k)/( dz(k)+dz(k-1) )
abk = 1.0 -afk
qkw(k) = SQRT(MAX(qke(k)*abk+qke(k-1)*afk,1.0e-3))
qkw(k) = SQRT(MAX(qke(k)*abk+qke(k-1)*afk, qkemin))
qtke(k) = 0.5*(qkw(k)**2) ! q -> TKE
thetaw(k)= theta(k)*abk + theta(k-1)*afk
END DO
@@ -2039,14 +2039,14 @@ SUBROUTINE mym_length ( &
zwk = zw(k)
DO WHILE (zwk .LE. zi2+h1)
dzk = 0.5*( dz(k)+dz(k-1) )
qdz = min(max( qkw(k)-qmin, 0.02 ), 30.0)*dzk
qdz = min(max( qkw(k)-qmin, 0.01 ), 30.0)*dzk
elt = elt +qdz*zwk
vsc = vsc +qdz
k = k+1
zwk = zw(k)
END DO

elt = MIN( MAX( alp1*elt/vsc, 10.), 400.)
elt = MIN( MAX( alp1*elt/vsc, 8.), 400.)
!avoid use of buoyancy flux functions which are ill-defined at the surface
!vflx = ( vt(kts)+1.0 )*flt + ( vq(kts)+tv0 )*flq
vflx = fltv
@@ -2122,13 +2122,13 @@ SUBROUTINE mym_length ( &
h1=MIN(h1,600.)
h2=h1*0.5 ! 1/4 transition layer depth

qtke(kts)=MAX(0.5*qke(kts),0.01) !tke at full sigma levels
qkw(kts) = SQRT(MAX(qke(kts),1.0e-4))
qtke(kts)=MAX(0.5*qke(kts), 0.5*qkemin) !tke at full sigma levels
qkw(kts) = SQRT(MAX(qke(kts), qkemin))

DO k = kts+1,kte
afk = dz(k)/( dz(k)+dz(k-1) )
abk = 1.0 -afk
qkw(k) = SQRT(MAX(qke(k)*abk+qke(k-1)*afk,1.0e-3))
qkw(k) = SQRT(MAX(qke(k)*abk+qke(k-1)*afk, qkemin))
qtke(k) = 0.5*qkw(k)**2 ! qkw -> TKE
END DO

@@ -3361,8 +3361,8 @@ SUBROUTINE mym_predict (kts,kte, &
CALL tridiag2(kte,a,b,c,d,x)

DO k=kts,kte
! qke(k)=max(d(k-kts+1), 1.e-4)
qke(k)=max(x(k), 1.e-4)
! qke(k)=max(d(k-kts+1), qkemin)
qke(k)=max(x(k), qkemin)
qke(k)=min(qke(k), 150.)
ENDDO

@@ -6509,11 +6509,11 @@ SUBROUTINE DMP_mf( &
do k=kts,kte-1
do I=1,nup
edmf_a(K) =edmf_a(K) +UPA(K,i)
edmf_w(K) =edmf_w(K) +rhoz(k)*UPA(K,i)*UPW(K,i)
edmf_qt(K) =edmf_qt(K) +rhoz(k)*UPA(K,i)*UPQT(K,i)
edmf_thl(K)=edmf_thl(K)+rhoz(k)*UPA(K,i)*UPTHL(K,i)
edmf_ent(K)=edmf_ent(K)+rhoz(k)*UPA(K,i)*ENT(K,i)
edmf_qc(K) =edmf_qc(K) +rhoz(k)*UPA(K,i)*UPQC(K,i)
edmf_w(K) =edmf_w(K) +UPA(K,i)*UPW(K,i)
edmf_qt(K) =edmf_qt(K) +UPA(K,i)*UPQT(K,i)
edmf_thl(K)=edmf_thl(K)+UPA(K,i)*UPTHL(K,i)
edmf_ent(K)=edmf_ent(K)+UPA(K,i)*ENT(K,i)
edmf_qc(K) =edmf_qc(K) +UPA(K,i)*UPQC(K,i)
enddo
enddo
do k=kts,kte-1
74 changes: 47 additions & 27 deletions physics/smoke_dust/module_add_emiss_burn.F90
Original file line number Diff line number Diff line change
@@ -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 )
@@ -27,38 +28,43 @@ 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
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
@@ -78,42 +84,51 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi, &
! 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))
@@ -122,8 +137,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)
@@ -143,7 +158,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))
@@ -154,6 +169,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)
@@ -163,7 +184,6 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi, &
enddo
enddo


END subroutine add_emis_burn

END module module_add_emiss_burn
98 changes: 45 additions & 53 deletions physics/smoke_dust/module_plumerise.F90
Original file line number Diff line number Diff line change
@@ -21,29 +21,28 @@ subroutine ebu_driver ( flam_frac,ebu_in,ebu, &
z_at_w,z,g,con_cp,con_rd, & ! scale_fire_emiss is part of config_flags
frp_inst, k_min, k_max, & ! RAR:
wind_eff_opt, &
kpbl_thetav, & ! SRB: added kpbl_thetav
kpbl_thetav, kpbl, & ! SRB: added kpbl_thetav and kpbl
curr_secs,xlat, xlong , uspdavg2d, &
hpbl2d, mpiid, alpha, &
frp_min, frp_wthreshold,zpbl_lim,uspd_lim, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte, errmsg, errflg,curr_secs, &
xlat, xlong , uspdavg2, hpbl_thetav2, mpiid)
its,ite, jts,jte, kts,kte, &
errmsg, errflg )

use rrfs_smoke_config
!use plume_data_mod
USE module_zero_plumegen_coms
USE module_smoke_plumerise
IMPLICIT NONE

REAL(kind_phys), PARAMETER :: frp_threshold= 1.e+7 ! Minimum FRP (Watts) to distribute smoke in PBL

REAL(kind_phys), PARAMETER :: zpbl_threshold = 2000. ! SRB: Minimum PBL depth to have plume rise
REAL(kind_phys), PARAMETER :: uspd_threshold = 5. ! SRB: Wind speed averaged across PBL depth to control smoke release levels
REAL(kind_phys), PARAMETER :: frp_threshold500 = 500.e+6 ! SRB: Minimum FRP (Watts) to have plume rise
REAL(kind_phys), intent(in) :: frp_min, frp_wthreshold, zpbl_lim, uspd_lim

real(kind=kind_phys), DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: frp_inst ! RAR: FRP array

real(kind_phys), DIMENSION(ims:ime,jms:jme), INTENT(IN) :: xlat,xlong ! SRB

real(kind_phys), DIMENSION(ims:ime, jms:jme), INTENT(IN) :: kpbl_thetav ! SRB
integer, dimension(ims:ime, jms:jme), intent(in) :: kpbl, kpbl_thetav

character(*), intent(inout) :: errmsg
integer, intent(inout) :: errflg
@@ -52,6 +51,7 @@ subroutine ebu_driver ( flam_frac,ebu_in,ebu, &
its,ite, jts,jte, kts,kte
real(kind_phys) :: curr_secs
INTEGER, INTENT(IN ) :: wind_eff_opt
REAL(kind_phys), INTENT(IN) :: alpha ! SRB: Enrainment constant for plumerise scheme
real(kind=kind_phys), DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT ) :: ebu
real(kind=kind_phys), INTENT(IN ) :: g, con_cp, con_rd
real(kind=kind_phys), DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: ebu_in
@@ -61,19 +61,19 @@ subroutine ebu_driver ( flam_frac,ebu_in,ebu, &

! Local variables...
INTEGER :: nv, i, j, k, kp1, kp2
INTEGER :: icall=0
INTEGER :: icall
INTEGER, DIMENSION(ims:ime, jms:jme), INTENT (OUT) :: k_min, k_max ! Min and max ver. levels for BB injection spread
REAL, DIMENSION(ims:ime, jms:jme), INTENT (OUT) :: uspdavg2, hpbl_thetav2 ! SRB
REAL, DIMENSION(ims:ime, jms:jme), INTENT (IN) :: uspdavg2d, hpbl2d ! SRB
real(kind_phys), dimension (kte) :: u_in ,v_in ,w_in ,theta_in ,pi_in, rho_phyin ,qv_in ,zmid, z_lev, uspd ! SRB
real(kind=kind_phys) :: dz_plume, cpor, con_rocp, uspdavg ! SRB
real(kind=kind_phys) :: dz_plume, cpor, con_rocp ! SRB

! MPI variables
INTEGER, INTENT(IN) :: mpiid

cpor =con_cp/con_rd
con_rocp=con_rd/con_cp

if (mod(int(curr_secs),1800) .eq. 0) then
if ( dbg_opt .and. (mod(int(curr_secs),1800) .eq. 0) ) then
icall = 0
endif

@@ -94,18 +94,6 @@ subroutine ebu_driver ( flam_frac,ebu_in,ebu, &
enddo
!enddo

! For now the flammable fraction is constant, based on the namelist. The next
! step to use LU index and meteorology to parameterize it
do j=jts,jte
do i=its,ite
flam_frac(i,j)= 0.
if (frp_inst(i,j) > frp_threshold) then
flam_frac(i,j)= 0.9
end if
enddo
enddo


! RAR: new FRP based approach
! Haiqin: do_plumerise is added to the namelist options
check_pl: IF (do_plumerise) THEN ! if the namelist option is set for plumerise
@@ -122,11 +110,10 @@ subroutine ebu_driver ( flam_frac,ebu_in,ebu, &
z_lev(k)= z_at_w(i,k,j)-z_at_w(i,kts,j)
rho_phyin(k)= rho_phy(i,k,j)
theta_in(k)= theta_phy(i,k,j)
uspd(k)= wind_phy(i,k,j) ! SRB
!uspd(k)= wind_phy(i,k,j) ! SRB
enddo


IF (dbg_opt .and. (icall .le. n_dbg_lines) .and. (frp_inst(i,j) .ge. frp_threshold) ) then
IF (dbg_opt .and. (icall .le. n_dbg_lines) .and. (frp_inst(i,j) .ge. frp_min) ) then
WRITE(1000+mpiid,*) 'module_plumerise_before:xlat,xlong,curr_secs,ebu(kts),frp_inst',xlat(i,j), xlong(i,j), int(curr_secs),ebu(i,kts,j),frp_inst(i,j)
WRITE(1000+mpiid,*) 'module_plumerise_before:xlat,xlong,curr_secs,u(10),v(10),w(10),qv(10)',xlat(i,j), xlong(i,j),int(curr_secs), u_in(10),v_in(10),w_in(kte),qv_in(10)
END IF
@@ -139,46 +126,51 @@ subroutine ebu_driver ( flam_frac,ebu_in,ebu, &
frp_inst(i,j), k_min(i,j), &
k_max(i,j), dbg_opt, g, con_cp, &
con_rd, cpor, errmsg, errflg, &
icall, mpiid, xlat(i,j), xlong(i,j), curr_secs )
icall, mpiid, xlat(i,j), xlong(i,j), &
curr_secs, alpha, frp_min )
if(errflg/=0) return

kp1= k_min(i,j)
kp2= k_max(i,j)
dz_plume= z_at_w(i,kp2,j) - z_at_w(i,kp1,j)

! SRB: Adding condition for overwriting plumerise levels
uspdavg=SUM(uspd(kts:kpbl_thetav(i,j)))/kpbl_thetav(i,j) !Average wind speed within the boundary layer
!uspdavg=SUM(uspd(kts:kpbl(i)))/kpbl(i) !Average wind speed within the boundary layer

! SRB: Adding output
uspdavg2(i,j) = uspdavg
hpbl_thetav2(i,j) = z_lev(kpbl_thetav(i,j))

IF ((frp_inst(i,j) .gt. frp_threshold) .AND. (frp_inst(i,j) .le. frp_threshold500) .AND. &
(z_lev(kpbl_thetav(i,j)) .gt. zpbl_threshold) .AND. (wind_eff_opt .eq. 1)) THEN
kp1=1
IF (uspdavg .ge. uspd_threshold) THEN ! Too windy
kp2=kpbl_thetav(i,j)/3
ELSE
kp2=kpbl_thetav(i,j)
END IF
dz_plume= z_at_w(i,kp2,j) - z_at_w(i,kp1,j)
do k=kp1,kp2-1
ebu(i,k,j)= ebu_in(i,j)* (z_at_w(i,k+1,j)-z_at_w(i,k,j))/dz_plume
enddo
!uspdavg2(i,j) = uspdavg
!hpbl_thetav2(i,j) = z_lev(kpbl(i))

IF (frp_inst(i,j) .le. frp_min) THEN
!kp1=1
!kp2=2
flam_frac(i,j)= 0.
ELSE IF ( (frp_inst(i,j) .le. frp_wthreshold) .AND. ( uspdavg2d(i,1) .ge. uspd_lim ) .AND. &
( hpbl2d(i,1) .gt. zpbl_lim) .AND. (wind_eff_opt .eq. 1)) THEN
kp1=2
kp2=MAX(3,NINT(real(kpbl(i,j))/3._kind_phys))
flam_frac(i,j)=0.85
ELSE
do k=kp1,kp2-1
ebu(i,k,j)= flam_frac(i,j)* ebu_in(i,j)* (z_at_w(i,k+1,j)-z_at_w(i,k,j))/dz_plume
enddo
ebu(i,kts,j)= (1.-flam_frac(i,j))* ebu_in(i,j)
flam_frac(i,j)=0.9 ! kp1,2 come from the plumerise scheme
END IF
! SRB: End modification

IF ( dbg_opt .and. (icall .le. n_dbg_lines) .and. (frp_inst(i,j) .ge. frp_threshold) ) then
! RAR: emission distribution
dz_plume= z_at_w(i,kp2,j) - z_at_w(i,kp1,j)
do k=kp1,kp2-1
ebu(i,k,j)=flam_frac(i,j)*ebu_in(i,j)*(z_at_w(i,k+1,j)-z_at_w(i,k,j))/dz_plume
enddo
ebu(i,kts,j)= (1.-flam_frac(i,j))* ebu_in(i,j)

! For output diagnostic
k_min(i,j) = kp1
k_max(i,j) = kp2

IF ( dbg_opt .and. (icall .le. n_dbg_lines) .and. (frp_inst(i,j) .ge. frp_min) ) then
WRITE(1000+mpiid,*) 'mod_plumerise_after:xlat,xlong,curr_secs,k_min(i,j), k_max(i,j) ',xlat(i,j),xlong(i,j),int(curr_secs),kp1,kp2
WRITE(1000+mpiid,*) 'mod_plumerise_after:xlat,xlong,curr_secs,ebu(kts),frp_inst',xlat(i,j),xlong(i,j),int(curr_secs),ebu(i,kts,j),frp_inst(i,j)
WRITE(1000+mpiid,*) 'mod_plumerise_after:xlat,xlong,curr_secs,u(10),v(10),w(10),qv(10)',xlat(i,j),xlong(i,j),int(curr_secs),u_in(10),v_in(10),w_in(kte),qv_in(10)
WRITE(1000+mpiid,*) 'mod_plumerise_after:xlat,xlong,curr_secs,uspdavg,kpbl_thetav',xlat(i,j),xlong(i,j),int(curr_secs),uspdavg,kpbl_thetav(i,j)
IF ( frp_inst(i,j) .ge. 3.e+9 ) then
!WRITE(1000+mpiid,*) 'mod_plumerise_after:xlat,xlong,curr_secs,uspdavg,kpbl_thetav,kpbl',xlat(i,j),xlong(i,j),int(curr_secs),uspdavg,kpbl_thetav(i,j),kpbl(i)
IF ( frp_inst(i,j) .ge. 2.e+10 ) then
WRITE(1000+mpiid,*) 'mod_plumerise_after:High FRP at : xlat,xlong,curr_secs,frp_inst',xlat(i,j),xlong(i,j),int(curr_secs),frp_inst(i,j)
END IF
icall = icall + 1
252 changes: 40 additions & 212 deletions physics/smoke_dust/module_smoke_plumerise.F90

Large diffs are not rendered by default.

12 changes: 9 additions & 3 deletions physics/smoke_dust/rrfs_smoke_config.F90
Original file line number Diff line number Diff line change
@@ -18,25 +18,31 @@ module rrfs_smoke_config
!-- aerosol module configurations
integer :: chem_opt = 1
integer :: kemit = 1
integer :: dust_opt = 5
integer :: dust_opt = 1
real(kind=kind_phys) :: dust_drylimit_factor = 1.0
real(kind=kind_phys) :: dust_moist_correction = 1.0
real(kind=kind_phys) :: dust_alpha = 0.
real(kind=kind_phys) :: dust_gamma = 0.
integer :: seas_opt = 0 ! turn off by default
logical :: do_plumerise = .true.
integer :: addsmoke_flag = 1
integer :: smoke_forecast = 1
integer :: plumerisefire_frq=60
integer :: n_dbg_lines = 3
integer :: wetdep_ls_opt = 1
integer :: drydep_opt = 1
integer :: pm_settling = 1
integer :: nfire_types = 5
integer :: ebb_dcycle = 2 ! 1: read in ebb_smoke(i,24), 2: daily
integer :: hwp_method = 2
logical :: dbg_opt = .true.
logical :: aero_ind_fdb = .false.
logical :: add_fire_heat_flux= .false.
logical :: add_fire_moist_flux= .false.
logical :: do_rrfs_sd = .true.
! integer :: wind_eff_opt = 1
integer :: plume_wind_eff = 1
logical :: extended_sd_diags = .false.
real(kind_phys) :: wetdep_ls_alpha = .5 ! scavenging factor
real(kind_phys) :: plume_alpha = 0.05

! --
integer, parameter :: CHEM_OPT_GOCART= 1
273 changes: 178 additions & 95 deletions physics/smoke_dust/rrfs_smoke_wrapper.F90

Large diffs are not rendered by default.

71 changes: 48 additions & 23 deletions physics/smoke_dust/rrfs_smoke_wrapper.meta
Original file line number Diff line number Diff line change
@@ -50,13 +50,6 @@
dimensions = ()
type = logical
intent = in
[add_fire_heat_flux_in]
standard_name = flag_for_fire_heat_flux
long_name = flag to add fire heat flux to LSM
units = flag
dimensions = ()
type = logical
intent = in
[do_plumerise_in]
standard_name = do_smoke_plumerise
long_name = rrfs smoke plumerise option
@@ -71,20 +64,20 @@
dimensions = ()
type = integer
intent = in
[n_dbg_lines_in]
standard_name = smoke_debug_lines
long_name = rrfs smoke add smoke option
units = index
dimensions = ()
type = integer
intent = in
[plume_wind_eff_in]
standard_name = option_for_wind_effects_on_smoke_plumerise
long_name = wind effect plumerise option
units = index
dimensions = ()
type = integer
intent = in
[add_fire_heat_flux_in]
standard_name = flag_for_fire_heat_flux
long_name = flag to add fire heat flux to LSM
units = flag
dimensions = ()
type = logical
intent = in
[addsmoke_flag_in]
standard_name = control_for_smoke_biomass_burning_emissions
long_name = rrfs smoke add smoke option
@@ -99,13 +92,28 @@
dimensions = ()
type = integer
intent = in
[smoke_forecast_in]
[hwp_method_in]
standard_name = do_smoke_forecast
long_name = index for rrfs smoke forecast
units = index
dimensions = ()
type = integer
intent = in
[add_fire_moist_flux_in]
standard_name = flag_for_fire_moisture_flux
long_name = flag to add fire moisture flux
units = flag
dimensions = ()
type = logical
intent = in
[plume_alpha_in]
standard_name = alpha_for_plumerise_scheme
long_name = alpha paramter for plumerise scheme
units = none
dimensions = ()
type = real
kind = kind_phys
intent = in
[dust_opt_in]
standard_name = control_for_smoke_dust
long_name = rrfs smoke dust chem option
@@ -188,6 +196,13 @@
dimensions = ()
type = integer
intent = out
[n_dbg_lines_in]
standard_name = smoke_debug_lines
long_name = rrfs smoke add smoke option
units = index
dimensions = ()
type = integer
intent = in

#####################################################################
[ccpp-arg-table]
@@ -598,7 +613,7 @@
standard_name = emission_smoke_prvd_RRFS
long_name = emission fire RRFS daily
units = various
dimensions = (horizontal_loop_extent,4)
dimensions = (horizontal_loop_extent,5)
type = real
kind = kind_phys
intent = in
@@ -821,6 +836,15 @@
dimensions = ()
type = integer
intent = in
[rho_dry]
standard_name = dry_air_density
long_name = dry air density
units = kg m-3
dimensions = (horizontal_loop_extent,vertical_layer_dimension)
type = real
kind = kind_phys
intent = inout
optional = True
[uspdavg]
standard_name = mean_wind_speed_in_boundary_layer
long_name = average wind speed within the boundary layer
@@ -945,13 +969,6 @@
kind = kind_phys
intent = out
optional = True
[kpbl]
standard_name = vertical_index_at_top_of_atmosphere_boundary_layer
long_name = PBL top model level index
units = index
dimensions = (horizontal_loop_extent)
type = integer
intent = in
[oro]
standard_name = height_above_mean_sea_level
long_name = height_above_mean_sea_level
@@ -969,6 +986,14 @@
kind = kind_phys
intent = in
optional = True
[totprcp]
standard_name = accumulated_lwe_thickness_of_precipitation_amount
long_name = accumulated total precipitation
units = m
dimensions = (horizontal_loop_extent)
type = real
kind = kind_phys
intent = in
[cpl_fire]
standard_name = do_fire_coupling
long_name = flag controlling fire_behavior collection (default off)

0 comments on commit 5d7c8f1

Please sign in to comment.