Skip to content

Commit

Permalink
Merge pull request #180 from haiqinli/production/RRFS.v1-mynn
Browse files Browse the repository at this point in the history 
[production/RRFS.v1] Update MYNN PBL & Smoke for RRFS.v1
  • Loading branch information
@dustinswales
dustinswales authored Mar 18, 2024
2 parents 09d7f2f + 0935794 commit 85fabc0
Show file tree
Hide file tree
Showing 9 changed files with 395 additions and 456 deletions.
7 changes: 2 additions & 5 deletions physics/Interstitials/UFS_SCM_NEPTUNE/sgscloud_radpre.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 )

Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down
7 changes: 7 additions & 0 deletions physics/Interstitials/UFS_SCM_NEPTUNE/sgscloud_radpre.meta
View file
Original file line number Diff line number Diff line change
Expand Up @@ -273,6 +273,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
Expand Down
38 changes: 19 additions & 19 deletions physics/PBL/MYNN_EDMF/module_bl_mynn.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -1932,11 +1932,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
Expand All @@ -1956,7 +1956,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
Expand Down Expand Up @@ -2014,14 +2014,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
Expand All @@ -2034,14 +2034,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
Expand Down Expand Up @@ -2117,13 +2117,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

Expand Down Expand Up @@ -3356,8 +3356,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

Expand Down Expand Up @@ -6504,11 +6504,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
Expand Down
95 changes: 49 additions & 46 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
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 85fabc0

Please sign in to comment.