From c67cb9ea59c63b3029820d07546d1e20ed46faea Mon Sep 17 00:00:00 2001 From: "haiqin.li" Date: Fri, 24 May 2024 17:38:54 +0000 Subject: [PATCH] "The last physics bug fix PR for RRFS.v1" --- physics/MP/Thompson/module_mp_thompson.F90 | 13 ++++--- physics/PBL/MYNN_EDMF/module_bl_mynn.F90 | 18 +++++----- .../SFC_Models/Land/RUC/module_sf_ruclsm.F90 | 8 ++--- physics/smoke_dust/module_add_emiss_burn.F90 | 34 +++++++++---------- physics/smoke_dust/rrfs_smoke_config.F90 | 1 + physics/smoke_dust/rrfs_smoke_wrapper.F90 | 28 ++++++++++----- physics/smoke_dust/rrfs_smoke_wrapper.meta | 12 +++++-- 7 files changed, 69 insertions(+), 45 deletions(-) diff --git a/physics/MP/Thompson/module_mp_thompson.F90 b/physics/MP/Thompson/module_mp_thompson.F90 index 11f2b27e2..327da88d5 100644 --- a/physics/MP/Thompson/module_mp_thompson.F90 +++ b/physics/MP/Thompson/module_mp_thompson.F90 @@ -5184,9 +5184,9 @@ end subroutine table_Efsw real function Eff_aero(D, Da, visc,rhoa,Temp,species) implicit none - real:: D, Da, visc, rhoa, Temp - character(LEN=1):: species - real:: aval, Cc, diff, Re, Sc, St, St2, vt, Eff + real(wp):: D, Da, visc, rhoa, Temp + character(LEN=1),intent(in):: species + real:: aval, Cc, diff, Re, Sc, St, St2, vt, Eff, rho_p real(wp), parameter:: boltzman = 1.3806503E-23 real(wp), parameter:: meanPath = 0.0256E-6 @@ -5194,10 +5194,13 @@ real function Eff_aero(D, Da, visc,rhoa,Temp,species) if (species .eq. 'r') then vt = -0.1021 + 4.932E3*D - 0.9551E6*D*D & + 0.07934E9*D*D*D - 0.002362E12*D*D*D*D + rho_p = rho_w elseif (species .eq. 's') then vt = av_s*D**bv_s + rho_p = rho_s elseif (species .eq. 'g') then vt = av_g*D**bv_g + rho_p = rho_g endif Cc = 1. + 2.*meanPath/Da *(1.257+0.4*exp(-0.55*Da/meanPath)) @@ -5206,8 +5209,8 @@ real function Eff_aero(D, Da, visc,rhoa,Temp,species) Re = 0.5*rhoa*D*vt/visc Sc = visc/(rhoa*diff) - St = Da*Da*vt*1000./(9.*visc*D) - aval = 1.+LOG(1.+Re) + St = (rho_p - rhoa)*Da*Da*vt*CC/(9.*visc*D) !Da*Da*vt*1000./(9.*visc*D) + aval = LOG(1.+Re) St2 = (1.2 + 1./12.*aval)/(1.+aval) Eff = 4./(Re*Sc) * (1. + 0.4*SQRT(Re)*Sc**0.3333 & diff --git a/physics/PBL/MYNN_EDMF/module_bl_mynn.F90 b/physics/PBL/MYNN_EDMF/module_bl_mynn.F90 index a0f647083..7c22c0387 100644 --- a/physics/PBL/MYNN_EDMF/module_bl_mynn.F90 +++ b/physics/PBL/MYNN_EDMF/module_bl_mynn.F90 @@ -1286,7 +1286,7 @@ SUBROUTINE mynn_bl_driver( & &s_aw1,s_awchem1, & &emis_ant_no(i), & &frp(i), rrfs_sd, & - &enh_mix, smoke_dbg ) + &enh_mix, smoke_dbg, znt(i) ) else call mynn_mix_chem(kts,kte,i, & &delt, dz1, pblh(i), & @@ -1298,7 +1298,7 @@ SUBROUTINE mynn_bl_driver( & &s_aw1,s_awchem1, & &zero, & &zero, rrfs_sd, & - &enh_mix, smoke_dbg ) + &enh_mix, smoke_dbg, znt(i) ) endif do ic = 1,nchem do k = kts,kte @@ -5229,7 +5229,7 @@ SUBROUTINE mynn_mix_chem(kts,kte,i, & dfh, & s_aw, s_awchem, & emis_ant_no, frp, rrfs_sd, & - enh_mix, smoke_dbg ) + enh_mix, smoke_dbg, znt ) !------------------------------------------------------------------- integer, intent(in) :: kts,kte,i @@ -5242,7 +5242,7 @@ SUBROUTINE mynn_mix_chem(kts,kte,i, & real(kind_phys), dimension( kts:kte, nchem ), intent(inout) :: chem1 real(kind_phys), dimension( kts:kte+1,nchem), intent(in) :: s_awchem real(kind_phys), dimension( ndvel ), intent(in) :: vd1 - real(kind_phys), intent(in) :: emis_ant_no,frp + real(kind_phys), intent(in) :: emis_ant_no,frp, znt logical, intent(in) :: rrfs_sd,enh_mix,smoke_dbg !local vars @@ -5260,8 +5260,8 @@ SUBROUTINE mynn_mix_chem(kts,kte,i, & real(kind_phys), dimension(kts:kte) :: rhoinv real(kind_phys), dimension(kts:kte+1) :: rhoz,khdz real(kind_phys), parameter :: NO_threshold = 10.0 ! For anthropogenic sources - real(kind_phys), parameter :: frp_threshold = 10.0 ! RAR 02/11/22: I increased the frp threshold to enhance mixing over big fires - real(kind_phys), parameter :: pblh_threshold = 100.0 + real(kind_phys), parameter :: frp_threshold = 20.0 ! RAR 02/11/22: I increased the frp threshold to enhance mixing over big fires + real(kind_phys), parameter :: pblh_threshold = 150.0 dztop=.5*(dz(kte)+dz(kte-1)) @@ -5302,9 +5302,11 @@ SUBROUTINE mynn_mix_chem(kts,kte,i, & khdz(k) = MAX(1.1*khdz(k),sqrt((emis_ant_no / NO_threshold)) / dz(k) * rhoz(k)) ! JLS 12/21/21 ! khdz(k) = MAX(khdz(k),khdz_back) ENDIF - IF ( frp > frp_threshold ) THEN + IF ( frp > frp_threshold .and. znt .ge. 0.6) THEN ! SRB: znt>0.6 correlates to forests [05/23/24] kmaxfire = ceiling(log(frp)) - khdz(k) = MAX(1.1*khdz(k), (1. - k/(kmaxfire*2.)) * ((log(frp))**2.- 2.*log(frp)) / dz(k)*rhoz(k)) ! JLS 12/21/21 + IF ( k .ge. 3 .and. k .le. kmaxfire ) THEN ! SRB: k=3 is ~60m above ground, avg evergreen forest canopy height [05/23/24] + khdz(k) = MAX(1.1*khdz(k), (1. - k/(real(kmaxfire)*2.)) * ((log(frp))**2.- 2.*log(frp)) / dz(k)*rhoz(k)) ! JLS 12/21/21 + ENDIF ! khdz(k) = MAX(khdz(k),khdz_back) ENDIF ENDIF diff --git a/physics/SFC_Models/Land/RUC/module_sf_ruclsm.F90 b/physics/SFC_Models/Land/RUC/module_sf_ruclsm.F90 index f1647ef81..9677e7bf1 100644 --- a/physics/SFC_Models/Land/RUC/module_sf_ruclsm.F90 +++ b/physics/SFC_Models/Land/RUC/module_sf_ruclsm.F90 @@ -1077,8 +1077,6 @@ SUBROUTINE LSMRUC(xlat,xlon, & tso(i,k,j) = tso1d(k) enddo - tso(i,nzs,j) = tbot(i,j) - do k=1,nzs smfr3d(i,k,j) = smfrkeep(k) keepfr3dflag(i,k,j) = keepfr (k) @@ -1104,8 +1102,10 @@ SUBROUTINE LSMRUC(xlat,xlon, & if(snow(i,j)==zero) EMISSL(i,j) = EMISBCK(i,j) EMISS (I,J) = EMISSL(I,J) ! SNOW is in [mm], SNWE is in [m]; CANWAT is in mm, CANWATR is in m - SNOW (i,j) = SNWE*1000._kind_phys - SNOWH (I,J) = SNHEI + !-- 17 may 2024 - cap snow for points at high elevations where all year round skin temperatures are close to 0 C + !-- Snow density for these points will be 3000/7.5=400 [kg/m^3] + SNOW (i,j) = min(3._kind_phys,SNWE)*1000._kind_phys ! cap to be < 3 m + SNOWH (I,J) = min(7.5_kind_phys,SNHEI) ! cap to be < 7.5 m CANWAT (I,J) = CANWATR*1000._kind_phys if (debug_print) then diff --git a/physics/smoke_dust/module_add_emiss_burn.F90 b/physics/smoke_dust/module_add_emiss_burn.F90 index ff1fb09b6..7466a0a8c 100755 --- a/physics/smoke_dust/module_add_emiss_burn.F90 +++ b/physics/smoke_dust/module_add_emiss_burn.F90 @@ -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 ) @@ -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 @@ -55,12 +55,15 @@ 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 timeq= gmt*3600._kind_phys + real(time_int,4) timeq= mod(timeq,timeq_max) @@ -70,34 +73,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)-1.0)) !One hour delay is due to the latency of the RAVE files, 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(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:3) ! 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 (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 @@ -113,7 +113,7 @@ 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 + ! 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) @@ -121,7 +121,7 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi,ebb_min, & 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) @@ -146,7 +146,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 diff --git a/physics/smoke_dust/rrfs_smoke_config.F90 b/physics/smoke_dust/rrfs_smoke_config.F90 index 3df0c5303..67469d943 100755 --- a/physics/smoke_dust/rrfs_smoke_config.F90 +++ b/physics/smoke_dust/rrfs_smoke_config.F90 @@ -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 diff --git a/physics/smoke_dust/rrfs_smoke_wrapper.F90 b/physics/smoke_dust/rrfs_smoke_wrapper.F90 index 0bff3fbde..eed8d9116 100755 --- a/physics/smoke_dust/rrfs_smoke_wrapper.F90 +++ b/physics/smoke_dust/rrfs_smoke_wrapper.F90 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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)>-100. .AND. xlat(i,j)>25. .AND. xlat(i,j)<41.) then + fire_type(i,j) = 3 ! slash burn and wildfires in the east, eastern temperate forest ecosystem + else + fire_type(i,j) = 4 ! potential wildfires end if end if end do @@ -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 @@ -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 ) diff --git a/physics/smoke_dust/rrfs_smoke_wrapper.meta b/physics/smoke_dust/rrfs_smoke_wrapper.meta index 2d793875f..14ec275ae 100755 --- a/physics/smoke_dust/rrfs_smoke_wrapper.meta +++ b/physics/smoke_dust/rrfs_smoke_wrapper.meta @@ -93,8 +93,8 @@ type = integer intent = in [hwp_method_in] - standard_name = do_smoke_forecast - long_name = index for rrfs smoke forecast + standard_name = control_for_HWP_equation + long_name = control for HWP equation units = index dimensions = () type = integer @@ -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 = none + dimensions = () + type = real + kind = kind_phys + intent = in [plume_alpha_in] standard_name = alpha_for_plumerise_scheme long_name = alpha paramter for plumerise scheme