"The last physics bug fix PR for RRFS.v1"

physics/MP/Thompson/module_mp_thompson.F90

@@ -5184,20 +5184,23 @@ 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
 
       vt = 1.
       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                  &

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

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

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,15 +113,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
+             ! 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)
@@ -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

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

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             )

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