From 228c550f43a5924a908aea59e2355a9ef1cce914 Mon Sep 17 00:00:00 2001 From: "Haiqin.Li" Date: Tue, 5 Mar 2024 03:32:42 +0000 Subject: [PATCH 1/5] "Update MYNN PBL for RRFS.v1" --- physics/PBL/MYNN_EDMF/module_bl_mynn.F90 | 38 ++++++++++++------------ 1 file changed, 19 insertions(+), 19 deletions(-) diff --git a/physics/PBL/MYNN_EDMF/module_bl_mynn.F90 b/physics/PBL/MYNN_EDMF/module_bl_mynn.F90 index 40995e593..c60a1a017 100644 --- a/physics/PBL/MYNN_EDMF/module_bl_mynn.F90 +++ b/physics/PBL/MYNN_EDMF/module_bl_mynn.F90 @@ -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 From 728ecf11b7b3d449c61c59e8f39cb2745c0371b1 Mon Sep 17 00:00:00 2001 From: "Haiqin.Li" Date: Wed, 13 Mar 2024 03:44:24 +0000 Subject: [PATCH 2/5] "smoke updates for RRFS.v1" --- physics/smoke_dust/module_add_emiss_burn.F90 | 97 +++--- physics/smoke_dust/module_plumerise.F90 | 98 +++---- physics/smoke_dust/module_smoke_plumerise.F90 | 252 +++------------- physics/smoke_dust/rrfs_smoke_config.F90 | 12 +- physics/smoke_dust/rrfs_smoke_wrapper.F90 | 275 ++++++++++++------ physics/smoke_dust/rrfs_smoke_wrapper.meta | 71 +++-- 6 files changed, 370 insertions(+), 435 deletions(-) diff --git a/physics/smoke_dust/module_add_emiss_burn.F90 b/physics/smoke_dust/module_add_emiss_burn.F90 index 80d91bb0e..889e42b43 100755 --- a/physics/smoke_dust/module_add_emiss_burn.F90 +++ b/physics/smoke_dust/module_add_emiss_burn.F90 @@ -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,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 + 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)) @@ -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) @@ -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) 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 diff --git a/physics/smoke_dust/module_smoke_plumerise.F90 b/physics/smoke_dust/module_smoke_plumerise.F90 index 13016d929..9c784a608 100755 --- a/physics/smoke_dust/module_smoke_plumerise.F90 +++ b/physics/smoke_dust/module_smoke_plumerise.F90 @@ -1,43 +1,28 @@ !>\file module_smoke_plumerise.F90 !! This file contains the fire plume rise module. - -!------------------------------------------------------------------------- -!- 12 April 2016 -!- Implementing the fire radiative power (FRP) methodology for biomass burning -!- emissions and convective energy estimation. -!- Saulo Freitas, Gabriel Pereira (INPE/UFJS, Brazil) -!- Ravan Ahmadov, Georg Grell (NOAA, USA) -!- The flag "plumerise_flag" defines the method: -!- =1 => original method -!- =2 => FRP based -!------------------------------------------------------------------------- module module_smoke_plumerise use machine , only : kind_phys - !use plume_data_mod, only : num_frp_plume, p_frp_hr, p_frp_std - !tropical_forest, boreal_forest, savannah, grassland, & - ! wind_eff USE module_zero_plumegen_coms USE rrfs_smoke_config, only : n_dbg_lines - !real(kind=kind_phys),parameter :: rgas=r_d - !real(kind=kind_phys),parameter :: cpor=cp/r_d -CONTAINS + CONTAINS ! RAR: subroutine plumerise(m1,m2,m3,ia,iz,ja,jz, & up,vp,wp,theta,pp,dn0,rv,zt_rams,zm_rams, & wind_eff_opt, & frp_inst,k1,k2, dbg_opt, g, cp, rgas, & - cpor, errmsg, errflg, icall, mpiid, lat, long, curr_secs ) + cpor, errmsg, errflg, icall, mpiid, & + lat, long, curr_secs, alpha, frp_min ) implicit none LOGICAL, INTENT (IN) :: dbg_opt INTEGER, INTENT (IN) :: wind_eff_opt, mpiid - real(kind_phys), INTENT(IN) :: lat,long, curr_secs ! SRB + real(kind_phys), INTENT(IN) :: lat,long, curr_secs, alpha ! SRB -! INTEGER, PARAMETER :: ihr_frp=1, istd_frp=2!, imean_fsize=3, istd_fsize=4 ! RAR: + REAL(kind_phys), INTENT(IN) :: frp_min ! integer, intent(in) :: PLUMERISE_flag real(kind=kind_phys) :: frp_inst ! This is the instantenous FRP, at a given time step @@ -45,15 +30,11 @@ subroutine plumerise(m1,m2,m3,ia,iz,ja,jz, & integer :: ng,m1,m2,m3,ia,iz,ja,jz,ibcon,mynum,i,j,k,imm,ixx,ispc !,nspecies - INTEGER, INTENT (OUT) :: k1,k2 character(*), intent(inout) :: errmsg integer, intent(inout) :: errflg -! integer :: ncall = 0 integer :: kmt -! real(kind=kind_phys),dimension(m1,nspecies), intent(inout) :: eburn_out -! real(kind=kind_phys),dimension(nspecies), intent(in) :: eburn_in real(kind=kind_phys), dimension(m1,m2,m3) :: up, vp, wp,theta,pp,dn0,rv real(kind=kind_phys), dimension(m1) :: zt_rams,zm_rams @@ -61,16 +42,6 @@ subroutine plumerise(m1,m2,m3,ia,iz,ja,jz, & real(kind=kind_phys), dimension(2) :: ztopmax real(kind=kind_phys) :: q_smold_kgm2 - REAL(kind_phys), PARAMETER :: frp_threshold= 1.e+7 ! Minimum FRP (Watts) to have plume rise - -! From plumerise1.F routine - integer, parameter :: iveg_ag=1 -! integer, parameter :: tropical_forest = 1 -! integer, parameter :: boreal_forest = 2 -! integer, parameter :: savannah = 3 -! integer, parameter :: grassland = 4 -! real(kind=kind_phys), dimension(nveg_agreg) :: firesize,mean_fct - INTEGER :: wind_eff INTEGER, INTENT(IN) :: icall type(plumegen_coms), pointer :: coms @@ -78,37 +49,10 @@ subroutine plumerise(m1,m2,m3,ia,iz,ja,jz, & ! Set wind effect from namelist wind_eff = wind_eff_opt -! integer:: iloop - !REAL(kind=kind_phys), INTENT (IN) :: convert_smold_to_flam - - !Fator de conversao de unidades - !!fcu=1. !=> kg [gas/part] /kg [ar] - !!fcu =1.e+12 !=> ng [gas/part] /kg [ar] - !!real(kind=kind_phys),parameter :: fcu =1.e+6 !=> mg [gas/part] /kg [ar] - !---------------------------------------------------------------------- - ! indexacao para o array "plume(k,i,j)" - ! k - ! 1 => area media (m^2) dos focos em biomas floresta dentro do gribox i,j - ! 2 => area media (m^2) dos focos em biomas savana dentro do gribox i,j - ! 3 => area media (m^2) dos focos em biomas pastagem dentro do gribox i,j - ! 4 => desvio padrao da area media (m^2) dos focos : floresta - ! 5 => desvio padrao da area media (m^2) dos focos : savana - ! 6 => desvio padrao da area media (m^2) dos focos : pastagem - ! 7 a 9 => sem uso - !10(=k_CO_smold) => parte da emissao total de CO correspondente a fase smoldering - !11, 12 e 13 => este array guarda a relacao entre - ! qCO( flaming, floresta) e a quantidade total emitida - ! na fase smoldering, isto e; - ! qCO( flaming, floresta) = plume(11,i,j)*plume(10,i,j) - ! qCO( flaming, savana ) = plume(12,i,j)*plume(10,i,j) - ! qCO( flaming, pastagem) = plume(13,i,j)*plume(10,i,j) - !20(=k_PM25_smold),21,22 e 23 o mesmo para PM25 - ! - !24-n1 => sem uso !---------------------------------------------------------------------- ! print *,' Plumerise_scalar 1',ncall coms => get_thread_coms() - + ! print *,' Plumerise_scalar 2',m1 j=1 i=1 @@ -131,12 +75,6 @@ subroutine plumerise(m1,m2,m3,ia,iz,ja,jz, & coms%zcon (k)=zt_rams(k) ! termod-point height coms%zzcon (k)=zm_rams(k) ! W-point height enddo - -! do ispc=2,nspecies - ! eburn_out(1,ispc) = eburn_in(ispc) ! eburn_in is the emissions at the 1st level -! eburn_out(2:m1,ispc)= 0. ! RAR: k>1 are used from eburn_out -! enddo - !- get envinronmental state (temp, water vapor mix ratio, ...) call get_env_condition(coms,1,m1,kmt,wind_eff,g,cp,rgas,cpor,errmsg,errflg) if(errflg/=0) return @@ -146,39 +84,36 @@ subroutine plumerise(m1,m2,m3,ia,iz,ja,jz, & ! iloop=1 ! IF (PLUMERISE_flag == 1) iloop=nveg_agreg - !lp_veg: do iveg_ag=1,iloop - FRP = max(1000.,frp_inst) + !frp_inst = max(1000.,frp_inst) - !- loop over the minimum and maximum heat fluxes/FRP + !- loop over the minimum and maximum heat fluxes/frp_inst lp_minmax: do imm=1,2 if(imm==1 ) then - burnt_area = 0.7* 0.0006* FRP ! 0.00021* FRP ! - 0.5*plume_fre(istd_fsize)) + burnt_area = 0.7* 0.0006* frp_inst ! 0.00021* frp_inst ! - 0.5*plume_fre(istd_fsize)) elseif(imm==2 ) then - burnt_area = 1.3* 0.0006* FRP ! RAR: Based on Laura's paper I increased the fire size *3. This should depend on the fuel type and meteorology/HWP + burnt_area = 1.3* 0.0006* frp_inst ! RAR: Based on Laura's paper I increased the fire size *3. This should depend on the fuel type and meteorology/HWP endif burnt_area= max(1.0e4,burnt_area) - IF ( dbg_opt .and. (icall .le. n_dbg_lines) .and. (frp_inst .ge. frp_threshold) ) THEN + IF ( dbg_opt .and. (icall .le. n_dbg_lines) .and. (frp_inst .ge. frp_min) ) THEN WRITE(1000+mpiid,*) 'inside plumerise: xlat,xlong,curr_secs, m1 ', lat,long, int(curr_secs), m1 - WRITE(1000+mpiid,*) 'inside plumerise: xlat,xlong,curr_secs,imm,FRP,burnt_area ', lat, long, int(curr_secs), imm, FRP,burnt_area + WRITE(1000+mpiid,*) 'inside plumerise: xlat,xlong,curr_secs,imm,frp_inst,burnt_area ', lat, long, int(curr_secs), imm, frp_inst,burnt_area END IF - IF (frp_inst=k1+1 - - !- emission during flaming phase is evenly distributed between levels k1 and k2 - !do k=k1,k2 - ! do ispc= 2,nspecies - ! eburn_out(k,ispc)= dzi* eburn_in(ispc) - ! enddo - !enddo - - !IF (dbg_opt) then - ! WRITE(*,*) 'plumerise after set_flam_vert: nkp,k1,k2, ', nkp,k1,k2 - ! WRITE(*,*) 'plumerise after set_flam_vert: dzi ', dzi - !WRITE(*,*) 'plumerise after set_flam_vert: eburn_in(2) ', eburn_in(2) - !WRITE(*,*) 'plumerise after set_flam_vert: eburn_out(:,2) ',eburn_out(:,2) - !END IF - -! enddo lp_veg ! sub-grid vegetation, currently it's aggregated - end subroutine plumerise !------------------------------------------------------------------------- @@ -285,22 +199,6 @@ subroutine get_env_condition(coms,k1,k2,kmt,wind_eff,g,cp,rgas,cpor,errmsg,errfl !-ewe - env wind effect if(wind_eff < 1) coms%vel_e(1:kmt) = 0. -!-use este para gerar o RAMS.out -! ------- print environment state -!print*,'k,coms%zt(k),coms%pe(k),coms%te(k)-273.15,coms%qvenv(k)*1000' -!do k=1,kmt -! write(*,100) k,coms%zt(k),coms%pe(k),coms%te(k)-273.15,coms%qvenv(k)*1000. -! 100 format(1x,I5,4f20.12) -!enddo -!stop 333 - - -!--------- nao eh necessario este calculo -!do k=1,kmt -! call thetae(coms%pe(k),coms%te(k),coms%qvenv(k),coms%thee(k)) -!enddo - - !--------- converte press de Pa para kPa para uso modelo de plumerise do k=1,kmt coms%pe(k) = coms%pe(k)*1.e-3 @@ -383,62 +281,15 @@ SUBROUTINE set_flam_vert(ztopmax,k1,k2,nkp,zzcon) !,W_VMD,VMD) !print*,'2: ztopmax k=',ztopmax(2), k2 k2= k1+1 ! RAR: I added k1+1 ENDIF - - !- version 2 - !- vertical mass distribution - !- -! w_thresold = 1. -! DO imm=1,2 - -! VMD(1:nkp,imm)= 0. -! xxx=0. -! k_initial= 0 -! k_final = 0 - - !- define range of the upper detrainemnt layer -! do ko=nkp-10,2,-1 - -! if(w_vmd(ko,imm) < w_thresold) cycle - -! if(k_final==0) k_final=ko - -! if(w_vmd(ko,imm)-1. > w_vmd(ko-1,imm)) then -! k_initial=ko -! exit -! endif - -! enddo - !- if there is a non zero depth layer, make the mass vertical distribution -! if(k_final > 0 .and. k_initial > 0) then - -! k_initial=int((k_final+k_initial)*0.5) - - !- parabolic vertical distribution between k_initial and k_final -! kk4 = k_final-k_initial+2 -! do ko=1,kk4-1 -! kl=ko+k_initial-1 -! VMD(kl,imm) = 6.* float(ko)/float(kk4)**2 * (1. - float(ko)/float(kk4)) -! enddo -! if(sum(VMD(1:NKP,imm)) .ne. 1.) then -! xxx= ( 1.- sum(VMD(1:NKP,imm)) )/float(k_final-k_initial+1) -! do ko=k_initial,k_final -! VMD(ko,imm) = VMD(ko,imm)+ xxx !- values between 0 and 1. -! enddo - ! print*,'new mass=',sum(mass)*100.,xxx - !pause -! endif -! endif !k_final > 0 .and. k_initial > - -! ENDDO - + END SUBROUTINE set_flam_vert !------------------------------------------------------------------------- -subroutine get_fire_properties(coms,imm,iveg_ag,burnt_area,FRP,errmsg,errflg) +subroutine get_fire_properties(coms,imm,burnt_area,FRP,errmsg,errflg) !use module_zero_plumegen_coms implicit none type(plumegen_coms), pointer :: coms -integer :: moist, i, icount,imm,iveg_ag !,plumerise_flag +integer :: moist, i, icount,imm real(kind=kind_phys):: bfract, effload, heat, hinc ,burnt_area,heat_fluxW,FRP !real(kind=kind_phys), dimension(2,4) :: heat_flux integer, intent(inout) :: errflg @@ -448,13 +299,8 @@ subroutine get_fire_properties(coms,imm,iveg_ag,burnt_area,FRP,errmsg,errflg) !real(kind=kind_phys), parameter :: beta = 5.0 !ref.: Wooster et al., 2005 REAL(kind=kind_phys), parameter :: beta = 0.88 !ref.: Paugam et al., 2015 -! coms%area = burnt_area! area of burn, m^2 -!IF ( PLUMERISE_flag == 1) THEN -! !fluxo de calor para o bioma -! heat_fluxW = heat_flux(imm,iveg_ag) * 1000. ! converte para W/m^2 - !ELSEIF ( PLUMERISE_flag == 2) THEN ! "beta" factor converts FRP to convective energy heat_fluxW = beta*(FRP/coms%area)/0.55 ! in W/m^2 @@ -475,25 +321,9 @@ subroutine get_fire_properties(coms,imm,iveg_ag,burnt_area,FRP,errmsg,errflg) !heat = 15.5e6 !joules/kg - cerrado heat = 19.3e6 !joules/kg - floresta em alta floresta (mt) !coms%alpha = 0.1 !- entrainment constant -coms%alpha = 0.05 !- entrainment constant +!coms%alpha = 0.05 !- entrainment constant -!-------------------- printout ---------------------------------------- - -!!WRITE ( * , * ) ' SURFACE =', COMS%ZSURF, 'M', ' LCL =', COMS%ZBASE, 'M' -! -!PRINT*,'=======================================================' -!print * , ' FIRE BOUNDARY CONDITION :' -!print * , ' DURATION OF BURN, MINUTES =',COMS%MDUR -!print * , ' AREA OF BURN, HA =',COMS%AREA*1.e-4 -!print * , ' HEAT FLUX, kW/m^2 =',heat_fluxW*1.e-3 -!print * , ' TOTAL LOADING, KG/M**2 =',COMS%BLOAD -!print * , ' FUEL MOISTURE, % =',MOIST !average fuel moisture,percent dry -!print * , ' MODEL TIME, MIN. =',COMS%MAXTIME -! -! -! ! ******************** fix up inputs ********************************* -! !IF (MOD (COMS%MAXTIME, 2) .NE.0) COMS%MAXTIME = COMS%MAXTIME+1 !make coms%maxtime even @@ -503,12 +333,10 @@ subroutine get_fire_properties(coms,imm,iveg_ag,burnt_area,FRP,errmsg,errflg) COMS%FMOIST = MOIST / 100. !- fuel moisture fraction ! -! ! calculate the energy flux and water content at lboundary. ! fills heating() on a minute basis. could ask for a file at this po ! in the program. whatever is input has to be adjusted to a one ! minute timescale. -! DO I = 1, ntime !- make sure of energy release COMS%HEATING (I) = 0.0001 !- avoid possible divide by 0 @@ -564,7 +392,7 @@ subroutine get_fire_properties(coms,imm,iveg_ag,burnt_area,FRP,errmsg,errflg) end subroutine get_fire_properties !------------------------------------------------------------------------------- ! -SUBROUTINE MAKEPLUME (coms,kmt,ztopmax,ixx,imm,mpiid) +SUBROUTINE MAKEPLUME (coms,kmt,ztopmax,ixx,imm,mpiid, alpha) ! ! ********************************************************************* ! @@ -621,7 +449,6 @@ SUBROUTINE MAKEPLUME (coms,kmt,ztopmax,ixx,imm,mpiid) ! ALPHA = ENTRAINMENT CONSTANT ! MAXTIME = TERMINATION TIME (MIN) ! -! !********************************************************************** !********************************************************************** !use module_zero_plumegen_coms @@ -633,7 +460,7 @@ SUBROUTINE MAKEPLUME (coms,kmt,ztopmax,ixx,imm,mpiid) ,ixx,nrectotal,i_micro,n_sub_step real(kind=kind_phys) :: vc, g, r, cp, eps, & tmelt, heatsubl, heatfus, heatcond, tfreeze, & - ztopmax, wmax, rmaxtime, es, esat, heat,dt_save !ESAT_PR, + ztopmax, wmax, rmaxtime, es, esat, heat,dt_save, alpha !ESAT_PR, character (len=2) :: cixx integer, intent(in) :: mpiid ! Set threshold to be the same as dz=100., the constant grid spacing of plume grid model(meters) found in set_grid() @@ -672,7 +499,7 @@ SUBROUTINE MAKEPLUME (coms,kmt,ztopmax,ixx,imm,mpiid) COMS%L = 1 ! COMS%L initialization !--- initialization -CALL INITIAL(coms,kmt) +CALL INITIAL(coms,kmt,alpha) !--- initial print fields: izprint = 0 ! if = 0 => no printout @@ -720,7 +547,7 @@ SUBROUTINE MAKEPLUME (coms,kmt,ztopmax,ixx,imm,mpiid) !-- bounday conditions (k=1) COMS%L=1 - call lbound(coms) + call lbound(coms, alpha) !-- dynamics for the level k>1 !-- W advection @@ -734,10 +561,10 @@ SUBROUTINE MAKEPLUME (coms,kmt,ztopmax,ixx,imm,mpiid) !call scl_advectc_plumerise2(coms,'SC',COMS%NM1) !-- scalars entrainment, adiabatic - call scl_misc(coms,COMS%NM1) + call scl_misc(coms,COMS%NM1, alpha) !-- scalars dinamic entrainment - call scl_dyn_entrain(COMS%NM1,nkp,coms%wbar,coms%w,coms%adiabat,coms%alpha,coms%radius,coms%tt,coms%t,coms%te,coms%qvt,coms%qv,coms%qvenv,coms%qct,coms%qc,coms%qht,coms%qh,coms%qit,coms%qi,& + call scl_dyn_entrain(COMS%NM1,nkp,coms%wbar,coms%w,coms%adiabat,alpha,coms%radius,coms%tt,coms%t,coms%te,coms%qvt,coms%qv,coms%qvenv,coms%qct,coms%qc,coms%qht,coms%qh,coms%qit,coms%qi,& coms%vel_e,coms%vel_p,coms%vel_t,coms%rad_p,coms%rad_t) !-- gravity wave damping using Rayleigh friction layer fot COMS%T @@ -794,7 +621,7 @@ SUBROUTINE MAKEPLUME (coms,kmt,ztopmax,ixx,imm,mpiid) call buoyancy_plumerise(COMS%NM1, COMS%T, COMS%TE, COMS%QV, COMS%QVENV, COMS%QH, COMS%QI, COMS%QC, COMS%WT, COMS%SCR1) !-- Entrainment - call entrainment(coms,COMS%NM1,COMS%W,COMS%WT,COMS%RADIUS,COMS%ALPHA) + call entrainment(coms,COMS%NM1,COMS%W,COMS%WT,COMS%RADIUS,alpha) !-- update W call update_plumerise(coms,coms%nm1,'W') @@ -912,7 +739,7 @@ SUBROUTINE BURN(COMS, EFLUX, WATER) END SUBROUTINE BURN !------------------------------------------------------------------------------- ! -SUBROUTINE LBOUND (coms) +SUBROUTINE LBOUND (coms, alpha) ! ! ********** BOUNDARY CONDITIONS AT ZSURF FOR PLUME AND CLOUD ******** ! @@ -934,7 +761,7 @@ SUBROUTINE LBOUND (coms) real(kind=kind_phys), parameter :: tfreeze = 269.3, pi = 3.14159, e1 = 1./3., e2 = 5./3. real(kind=kind_phys) :: es, esat, eflux, water, pres, c1, c2, f, zv, denscor, xwater !,ESAT_PR ! real(kind=kind_phys), external:: esat_pr! - +REAL(kind=kind_phys) , INTENT(IN) :: alpha ! COMS%QH (1) = COMS%QH (2) !soak up hydrometeors COMS%QI (1) = COMS%QI (2) @@ -947,9 +774,9 @@ SUBROUTINE LBOUND (coms) ! PRES = COMS%PE (1) * 1000. !need pressure in N/m**2 - C1 = 5. / (6. * COMS%ALPHA) !alpha is entrainment constant + C1 = 5. / (6. * alpha) !alpha is entrainment constant - C2 = 0.9 * COMS%ALPHA + C2 = 0.9 * alpha F = EFLUX / (PRES * CP * PI) @@ -1016,7 +843,7 @@ SUBROUTINE LBOUND (coms) END SUBROUTINE LBOUND !------------------------------------------------------------------------------- ! -SUBROUTINE INITIAL (coms,kmt) +SUBROUTINE INITIAL (coms,kmt,alpha) ! ! ************* SETS UP INITIAL CONDITIONS FOR THE PROBLEM ************ !use module_zero_plumegen_coms @@ -1025,6 +852,7 @@ SUBROUTINE INITIAL (coms,kmt) real(kind=kind_phys), parameter :: tfreeze = 269.3 integer :: isub, k, n1, n2, n3, lbuoy, itmp, isubm1 ,kmt real(kind=kind_phys) :: xn1, xi, es, esat!,ESAT_PR +REAL(kind=kind_phys) , INTENT(IN) :: alpha ! COMS%N=kmt ! initialize temperature structure,to the end of equal spaced sounding, @@ -1058,13 +886,13 @@ SUBROUTINE INITIAL (coms,kmt) ! Initialize the entrainment radius, Turner-style plume coms%radius(1) = coms%rsurf do k=2,COMS%N - coms%radius(k) = coms%radius(k-1)+(6./5.)*coms%alpha*(coms%zt(k)-coms%zt(k-1)) + coms%radius(k) = coms%radius(k-1)+(6./5.)*alpha*(coms%zt(k)-coms%zt(k-1)) enddo ! Initialize the entrainment radius, Turner-style plume coms%radius(1) = coms%rsurf coms%rad_p(1) = coms%rsurf DO k=2,COMS%N - coms%radius(k) = coms%radius(k-1)+(6./5.)*coms%alpha*(coms%zt(k)-coms%zt(k-1)) + coms%radius(k) = coms%radius(k-1)+(6./5.)*alpha*(coms%zt(k)-coms%zt(k-1)) coms%rad_p(k) = coms%radius(k) ENDDO @@ -1080,7 +908,7 @@ SUBROUTINE INITIAL (coms,kmt) !ENDDO !stop 333 - CALL LBOUND(COMS) + CALL LBOUND(COMS, alpha) RETURN END SUBROUTINE INITIAL @@ -1537,21 +1365,21 @@ end subroutine tend0_plumerise ! **************************************************************** -subroutine scl_misc(coms,m1) +subroutine scl_misc(coms,m1,alpha) !use module_zero_plumegen_coms implicit none type(plumegen_coms), pointer :: coms real(kind=kind_phys), parameter :: g = 9.81, cp=1004. integer m1,k real(kind=kind_phys) dmdtm - +REAL(kind=kind_phys) , INTENT(IN) :: alpha do k=2,m1-1 COMS%WBAR = 0.5*(COMS%W(k)+COMS%W(k-1)) !-- dry adiabat COMS%ADIABAT = - COMS%WBAR * G / CP ! !-- entrainment - DMDTM = 2. * COMS%ALPHA * ABS (COMS%WBAR) / COMS%RADIUS (k) != (1/M)DM/COMS%DT + DMDTM = 2. * alpha * ABS (COMS%WBAR) / COMS%RADIUS (k) != (1/M)DM/COMS%DT !-- tendency temperature = adv + adiab + entrainment COMS%TT(k) = COMS%TT(K) + COMS%ADIABAT - DMDTM * ( COMS%T (k) - COMS%TE (k) ) diff --git a/physics/smoke_dust/rrfs_smoke_config.F90 b/physics/smoke_dust/rrfs_smoke_config.F90 index c20d6e2db..3df0c5303 100755 --- a/physics/smoke_dust/rrfs_smoke_config.F90 +++ b/physics/smoke_dust/rrfs_smoke_config.F90 @@ -18,11 +18,14 @@ 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 @@ -30,13 +33,16 @@ module rrfs_smoke_config 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 diff --git a/physics/smoke_dust/rrfs_smoke_wrapper.F90 b/physics/smoke_dust/rrfs_smoke_wrapper.F90 index 504014e6a..87817fb7d 100755 --- a/physics/smoke_dust/rrfs_smoke_wrapper.F90 +++ b/physics/smoke_dust/rrfs_smoke_wrapper.F90 @@ -1,6 +1,6 @@ -!>\file rrfs_smoke_wrapper.F90 +!>\file rrfs_smoke_wrapper.F90hwp_method !! This file is CCPP driver of RRFS Smoke and Dust -!! Haiqin.Li@noaa.gov 02/2021 +!! Haiqin.Li@noaa.gov 03/2024 module rrfs_smoke_wrapper @@ -9,11 +9,12 @@ module rrfs_smoke_wrapper use rrfs_smoke_config, only : kemit, dust_opt, seas_opt, do_plumerise, & addsmoke_flag, plumerisefire_frq, wetdep_ls_opt, & drydep_opt, pm_settling, aero_ind_fdb, ebb_dcycle, & - dbg_opt,smoke_forecast,wetdep_ls_alpha,do_rrfs_sd, & + dbg_opt,hwp_method,wetdep_ls_alpha,do_rrfs_sd, & ebb_dcycle, extended_sd_diags,add_fire_heat_flux, & num_moist, num_chem, num_emis_seas, num_emis_dust, & - p_qv, p_atm_shum, p_atm_cldq, & - p_smoke, p_dust_1, p_coarse_pm, epsilc, n_dbg_lines + 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 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 @@ -31,8 +32,6 @@ module rrfs_smoke_wrapper public :: rrfs_smoke_wrapper_run, rrfs_smoke_wrapper_init - integer :: plume_wind_eff - contains !>\defgroup rrfs_smoke_wrapper rrfs-sd emission driver Module @@ -43,31 +42,31 @@ module rrfs_smoke_wrapper !! \htmlinclude rrfs_smoke_wrapper_init.html !! subroutine rrfs_smoke_wrapper_init( seas_opt_in, & ! sea salt namelist - drydep_opt_in, pm_settling_in, & ! Dry Dep namelist - wetdep_ls_opt_in,wetdep_ls_alpha_in, & ! Wet dep namelist + drydep_opt_in, pm_settling_in, & ! dry dep namelist + wetdep_ls_opt_in,wetdep_ls_alpha_in, & ! wet dep namelist 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, smoke_forecast_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 - aero_ind_fdb_in, & ! Feedback namelist - extended_sd_diags_in,dbg_opt_in, & ! Other namelist + addsmoke_flag_in, ebb_dcycle_in, hwp_method_in, & ! smoke namelist + add_fire_moist_flux_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 + aero_ind_fdb_in, & ! feedback namelist + extended_sd_diags_in,dbg_opt_in, & ! other namelist errmsg, errflg, n_dbg_lines_in ) - - + !>-- Namelist - real(kind_phys), intent(in) :: dust_alpha_in, dust_gamma_in, wetdep_ls_alpha_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 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 - integer, intent(in) :: smoke_forecast_in, plume_wind_eff_in, plumerisefire_frq_in, n_dbg_lines_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 + integer, intent(in) :: hwp_method_in, plume_wind_eff_in, plumerisefire_frq_in, n_dbg_lines_in integer, intent(in) :: addsmoke_flag_in, ebb_dcycle_in logical, intent(in) :: do_plumerise_in, rrfs_sd character(len=*),intent(out):: errmsg - integer, intent(out) :: errflg + integer, intent(out) :: errflg errmsg = '' errflg = 0 @@ -93,9 +92,11 @@ subroutine rrfs_smoke_wrapper_init( seas_opt_in, do_plumerise = do_plumerise_in plumerisefire_frq = plumerisefire_frq_in addsmoke_flag = addsmoke_flag_in - smoke_forecast = smoke_forecast_in + hwp_method = hwp_method_in plume_wind_eff = plume_wind_eff_in add_fire_heat_flux = add_fire_heat_flux_in + add_fire_moist_flux = add_fire_moist_flux_in + plume_alpha = plume_alpha_in !>-Feedback aero_ind_fdb = aero_ind_fdb_in !>-Other @@ -126,11 +127,12 @@ subroutine rrfs_smoke_wrapper_run(im, flag_init, kte, kme, ktau, dt, garea, land max_fplume, hwp, hwp_ave, wetness, ndvel, ddvel_inout, & smoke_fire, cpl_fire, & peak_hr_out,lu_nofire_out,lu_qfire_out, & - fire_heat_flux_out, frac_grid_burned_out, kpbl,oro, & - uspdavg, hpbl_thetav, mpicomm, mpirank, mpiroot, errmsg,errflg ) + fire_heat_flux_out, frac_grid_burned_out,oro,totprcp, & + uspdavg, hpbl_thetav, rho_dry, & + mpicomm, mpirank, mpiroot, errmsg,errflg ) + implicit none - integer, intent(in) :: im,kte,kme,ktau,nsoil,tile_num,jdate(8),idat(8) integer, intent(in) :: ntrac, ntfsmoke, ntsmoke, ntdust, ntcoarsepm, ndvel, nlcat logical, intent(in) :: flag_init @@ -148,7 +150,7 @@ subroutine rrfs_smoke_wrapper_run(im, flag_init, kte, kme, ktau, dt, garea, land real(kind_phys), dimension(:,:), intent(in), optional :: emi_ant_in real(kind_phys), dimension(:), intent(in) :: u10m, v10m, ustar, dswsfc, & recmol, garea, rlat,rlon, tskin, pb2d, zorl, snow, & - rain_cpl, rainc_cpl, hf2d, t2m, dpt2m + rain_cpl, rainc_cpl, hf2d, t2m, dpt2m, totprcp real(kind_phys), dimension(:,:), intent(in) :: vegtype_frac real(kind_phys), dimension(:,:), intent(in) :: ph3d, pr3d real(kind_phys), dimension(:,:), intent(in) :: phl3d, prl3d, tk3d, us3d, vs3d, spechum, w @@ -157,6 +159,7 @@ subroutine rrfs_smoke_wrapper_run(im, flag_init, kte, kme, ktau, dt, garea, land real(kind_phys), dimension(:), intent(inout), optional :: emdust, emseas, emanoc real(kind_phys), dimension(:), intent(inout), optional :: ebb_smoke_in,coef_bb, frp_output, fhist real(kind_phys), dimension(:,:), intent(inout), optional :: ebu_smoke + real(kind_phys), dimension(:,:), intent(inout), optional :: rho_dry real(kind_phys), dimension(:), intent(out ), optional :: fire_heat_flux_out, frac_grid_burned_out real(kind_phys), dimension(:), intent(inout), optional :: max_fplume, min_fplume, uspdavg, hpbl_thetav real(kind_phys), dimension(:), intent(inout), optional :: hwp, peak_hr_out @@ -171,7 +174,6 @@ subroutine rrfs_smoke_wrapper_run(im, flag_init, kte, kme, ktau, dt, garea, land real(kind_phys), dimension(:), intent(in), optional :: smoke_fire logical, intent(in) :: cpl_fire integer, intent(in) :: imp_physics, imp_physics_thompson - integer, dimension(:), intent(in) :: kpbl real(kind_phys), dimension(:), intent(in) :: oro character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -201,11 +203,12 @@ subroutine rrfs_smoke_wrapper_run(im, flag_init, kte, kme, ktau, dt, garea, land integer, dimension(ims:im, jms:jme) :: isltyp, ivgtyp !>- plume variables ! -- buffers - real(kind_phys), dimension(ims:im, jms:jme ) :: coef_bb_dc, flam_frac, frp_in, & - fire_hist, peak_hr, lu_nofire, lu_qfire, ebu_in, & - fire_end_hr, hwp_day_avg, kpbl_thetav,& - uspdavg2, hpbl_thetav2 - integer, dimension(ims:im, jms:jme ) :: min_fplume2, max_fplume2, fire_type + real(kind_phys), dimension(ims:im, jms:jme ) :: coef_bb_dc, flam_frac, frp_in, & + fire_hist, peak_hr, lu_nofire, lu_qfire, lu_sfire, & + ebu_in, fire_end_hr, hwp_day_avg, & + uspdavg2d, hpbl2d, totprcp_24hrs + integer, dimension(ims:im, jms:jme ) :: min_fplume2, max_fplume2, fire_type, & + kpbl,kpbl_thetav logical :: call_plume, reset_hwp_ave, avg_hwp_ave !>- optical variables real(kind_phys), dimension(ims:im, jms:jme, ndvel) :: ddvel, settling_flux, drydep_flux_local @@ -220,6 +223,7 @@ subroutine rrfs_smoke_wrapper_run(im, flag_init, kte, kme, ktau, dt, garea, land !> -- aerosol density (kg/m3) real(kind_phys), parameter :: density_dust= 2.6e+3, density_sulfate=1.8e+3 real(kind_phys), parameter :: density_oc = 1.4e+3, density_seasalt=2.2e+3 + real(kind_phys), parameter :: conv_frpi = 1.e-06_kind_phys ! FRP conversion factor, MW to W real(kind_phys), dimension(im) :: daero_emis_wfa, daero_emis_ifa !> -- other real(kind_phys), dimension(im) :: wdgust, snoweq @@ -234,6 +238,14 @@ subroutine rrfs_smoke_wrapper_run(im, flag_init, kte, kme, ktau, dt, garea, land integer, intent(in) :: mpirank integer, intent(in) :: mpiroot +! FRP Thresholds + REAL(kind_phys), PARAMETER :: frp_min = 1.e+7 ! Minimum FRP (Watts) to distribute smoke in PBL, 10MW + REAL(kind_phys), PARAMETER :: frp_max = 2.e+10 ! Maximum FRP over 3km Pixel, 20,000 MW + REAL(kind_phys), PARAMETER :: zpbl_threshold = 2.e+3 ! Minimum PBL depth to have plume rise + REAL(kind_phys), PARAMETER :: uspd_threshold = 5. ! Wind speed averaged across PBL depth to control smoke release levels + REAL(kind_phys), PARAMETER :: frp_wthreshold = 1.e+9 ! Minimum FRP (Watts) to have plume rise in windy conditions + REAL(kind_phys), PARAMETER :: ebb_min = 1.e-3 ! Minimum smoke emissions (ug/m2/s) + mpiid = mpirank errmsg = '' @@ -262,13 +274,14 @@ subroutine rrfs_smoke_wrapper_run(im, flag_init, kte, kme, ktau, dt, garea, land min_fplume2 = 0 max_fplume2 = 0 - uspdavg2 = 0. - hpbl_thetav2 = 0. + uspdavg2d = 0. + hpbl2d = 0. emis_seas = 0. emis_dust = 0. peak_hr = 0. fire_type = 0 lu_qfire = 0. + lu_sfire = 0. lu_nofire = 0. flam_frac = 0. daero_emis_wfa = 0. @@ -316,16 +329,16 @@ subroutine rrfs_smoke_wrapper_run(im, flag_init, kte, kme, ktau, dt, garea, land nsoil,smc,tslb,vegtype_dom,soiltyp, & nlcat,vegtype_frac,dswsfc,zorl, & snow,dust12m_in,emi_ant_in,smoke_RRFS,smoke2d_RRFS,coef_bb_dc, & - hf2d, pb2d, g, pi, hour_int, peak_hr, & + hf2d, pb2d, g, pi, hour_int, peak_hr,uspdavg2d, & u10,v10,ust,tsk,xland,xlat,xlong,dxy, & rri,t_phy,u_phy,v_phy,p_phy,pi_phy,wind_phy,theta_phy, & rho_phy,dz8w,p8w,t8w,recmol, & - z_at_w,vvel,zmid, & - ntrac,gq0, & + z_at_w,vvel,zmid,hpbl2d, & + ntrac,gq0,totprcp, & num_chem,num_moist, & ntsmoke, ntdust,ntcoarsepm, & moist,chem,ebu_in,kpbl_thetav,ebb_smoke_in, & - fire_hist,frp_in, hwp_day_avg, fire_end_hr, & + fire_hist,frp_in, hwp_day_avg, totprcp_24hrs, fire_end_hr, & emis_anoc,smois,stemp,ivgtyp,isltyp,vegfrac,rmol,swdown,znt, & hfx,pbl,snowh,clayf,rdrag,sandf,ssm,uthr,oro, hwp_local, & t2m,dpt2m,wetness,kpbl, & @@ -357,20 +370,26 @@ subroutine rrfs_smoke_wrapper_run(im, flag_init, kte, kme, ktau, dt, garea, land if (ebb_dcycle==2) then do j=jts,jte do i=its,ite - if (ebu_in(i,j)<0.01) then + if (ebu_in(i,j)0.95) then - fire_type(i,j) = 0 - else if (lu_qfire(i,j)>0.95) then - fire_type(i,j) = 1 + ! Permanent wetlands, snow/ice, water, barren tundra: + lu_nofire(i,j)= vegfrac(i,11,j) + vegfrac(i,15,j) + vegfrac(i,17,j) + vegfrac(i,20,j) + ! 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) + 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 ! RAR: need to add another criteria for fire_type=2, i.e. prescribed fires + fire_type(i,j) = 3 ! wildfires, new approach is necessary for the controlled burns in the forest areas end if end if end do @@ -408,11 +427,11 @@ subroutine rrfs_smoke_wrapper_run(im, flag_init, kte, kme, ktau, dt, garea, land ! Every hour (per namelist) the ebu_driver is called to calculate ebu, but ! the plumerise is controlled by the namelist option of plumerise_flag if (add_fire_heat_flux) then - WRITE(1000+mpiid,*) 'Entered add_fire_heat_flux at timestep:',ktau + !WRITE(1000+mpiid,*) 'Entered add_fire_heat_flux at timestep:',ktau do i = its,ite if ( coef_bb_dc(i,1)*frp_in(i,1) .ge. 1.E7 ) then fire_heat_flux_out(i) = min(max(0.,0.88*coef_bb_dc(i,1)*frp_in(i,1) / & - 0.55/dxy(i,1)) ,5000.) ! JLS - W m-2 [0 - 10,000] + 0.55/dxy(i,1)) ,5000.) ! W m-2 [0 - 10,000] frac_grid_burned_out(i) = min(max(0., 1.3*0.0006*coef_bb_dc(i,1)*frp_in(i,1)/dxy(i,1) ),1.) else fire_heat_flux_out(i) = 0.0 @@ -424,7 +443,7 @@ subroutine rrfs_smoke_wrapper_run(im, flag_init, kte, kme, ktau, dt, garea, land ! Apply the diurnal cycle coefficient to frp_inst () do j=jts,jte do i=its,ite - frp_inst(i,j) = frp_in(i,j)*coef_bb_dc(i,j) + frp_inst(i,j) = MIN(frp_in(i,j)*coef_bb_dc(i,j),frp_max) enddo enddo @@ -435,21 +454,24 @@ subroutine rrfs_smoke_wrapper_run(im, flag_init, kte, kme, ktau, dt, garea, land z_at_w,zmid,g,con_cp,con_rd, & frp_inst, min_fplume2, max_fplume2, & plume_wind_eff, & - kpbl_thetav, & + kpbl_thetav,kpbl,curr_secs, & + xlat, xlong, uspdavg2d, hpbl2d, mpiid,plume_alpha, & + frp_min, frp_wthreshold, & + zpbl_threshold, uspd_threshold, & 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 ) if(errflg/=0) return end if ! -- add biomass burning emissions at every timestep if (addsmoke_flag == 1) then - call add_emis_burn(dt,dz8w,rho_phy,pi, & + call add_emis_burn(dt,dz8w,rho_phy,pi,ebb_min, & chem,julday,gmt,xlat,xlong, & 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, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte , mpiid ) @@ -524,9 +546,9 @@ subroutine rrfs_smoke_wrapper_run(im, flag_init, kte, kme, ktau, dt, garea, land !---- diagnostic output of hourly wildfire potential (07/2021) if (ktau == 1 .or. reset_hwp_ave) then - hwp_ave = 0. + hwp_ave = 0._kind_phys endif - hwp = 0. + hwp = 0._kind_phys do i=its,ite hwp(i)=hwp_local(i,1) hwp_ave(i) = hwp_ave(i) + hwp(i)*dt @@ -570,6 +592,13 @@ subroutine rrfs_smoke_wrapper_run(im, flag_init, kte, kme, ktau, dt, garea, land enddo !------------------------------------- !-- to output for diagnostics + + do k=kts,kte + do i=its,ite + rho_dry(i,k) = real(rho_phy(i,k,1)) + enddo + enddo + do i = 1, im ! RAR: let's remove the seas and ant. OC emseas (i) = emis_seas(i,1,1,1)*1.e+9 ! size bin 1 sea salt emission: ug/m2/s @@ -577,17 +606,15 @@ subroutine rrfs_smoke_wrapper_run(im, flag_init, kte, kme, ktau, dt, garea, land emdust (i) = emis_dust(i,1,1,1) + emis_dust(i,1,1,2) + & emis_dust(i,1,1,3) + emis_dust(i,1,1,4) ! dust emission: ug/m2/s coef_bb (i) = coef_bb_dc(i,1) - frp_output (i) = coef_bb_dc(i,1)*frp_in(i,1) + frp_output (i) = coef_bb_dc(i,1)*frp_in(i,1)*conv_frpi ! to get FRP output in MW fhist (i) = fire_hist (i,1) min_fplume (i) = real(min_fplume2(i,1)) max_fplume (i) = real(max_fplume2(i,1)) fire_type_out(i)=fire_type(i,1) lu_nofire_out(i)=lu_nofire(i,1) lu_qfire_out (i)=lu_qfire(i,1) - enddo - - do i = 1, im - peak_hr_out(i) = peak_hr(i,1) + uspdavg (i) = uspdavg2d(i,1) + peak_hr_out(i) = peak_hr(i,1) enddo !-- to provide real aerosol emission for Thompson MP @@ -625,16 +652,16 @@ subroutine rrfs_smoke_prep( & pr3d,ph3d,phl3d,tk3d,prl3d,us3d,vs3d,spechum,w, & nsoil,smc,tslb,vegtype_dom,soiltyp,nlcat,vegtype_frac,dswsfc,zorl, & snow_cpl,dust12m_in,emi_ant_in,smoke_RRFS,smoke2d_RRFS,coef_bb_dc, & - hf2d, pb2d, g, pi, hour_int, peak_hr, & + hf2d, pb2d, g, pi, hour_int, peak_hr,uspdavg2d, & u10,v10,ust,tsk,xland,xlat,xlong,dxy, & rri,t_phy,u_phy,v_phy,p_phy,pi_phy,wind_phy,theta_phy, & rho_phy,dz8w,p8w,t8w,recmol, & - z_at_w,vvel,zmid, & - ntrac,gq0, & + z_at_w,vvel,zmid,hpbl2d, & + ntrac,gq0,totprcp, & num_chem, num_moist, & ntsmoke, ntdust, ntcoarsepm, & moist,chem,ebu_in,kpbl_thetav,ebb_smoke_in, & - fire_hist,frp_in, hwp_day_avg, fire_end_hr, & + fire_hist,frp_in, hwp_day_avg, totprcp_24hrs, fire_end_hr, & emis_anoc,smois,stemp,ivgtyp,isltyp,vegfrac,rmol,swdown, & znt,hfx,pbl,snowh,clayf,rdrag,sandf,ssm,uthr,oro,hwp_local, & t2m,dpt2m,wetness,kpbl, & @@ -647,19 +674,20 @@ subroutine rrfs_smoke_prep( & !FV3 input variables integer, intent(in) :: nsoil, ktau - integer, dimension(ims:ime), intent(in) :: land, vegtype_dom, soiltyp, kpbl + integer, dimension(ims:ime), intent(in) :: land, vegtype_dom, soiltyp integer, intent(in) :: ntrac real(kind=kind_phys), intent(in) :: g, pi, gmt, con_rd, con_fv, con_cp real(kind=kind_phys), dimension(ims:ime), intent(in) :: & - u10m, v10m, ustar, garea, rlat, rlon, ts2d, dswsfc, & - zorl, snow_cpl, pb2d, hf2d, oro, t2m, dpt2m, wetness, recmol + u10m, v10m, ustar, garea, rlat, rlon, ts2d, dswsfc, & + zorl, snow_cpl, pb2d, hf2d, oro, t2m, dpt2m, wetness, recmol, & + totprcp real(kind=kind_phys), dimension(ims:ime, nlcat), intent(in) :: vegtype_frac real(kind=kind_phys), dimension(ims:ime, nsoil), intent(in) :: smc,tslb real(kind=kind_phys), dimension(ims:ime, 12, 5), intent(in) :: dust12m_in real(kind=kind_phys), dimension(ims:ime, 24, 2), intent(in) :: smoke_RRFS ! This is a place holder for ebb_dcycle == 2, currently set to hold a single ! value, which is the previous day's average of hwp, frp, ebb, fire_end - real(kind=kind_phys), dimension(ims:ime, 4), intent(in) :: smoke2d_RRFS + real(kind=kind_phys), dimension(ims:ime, 5), intent(in) :: smoke2d_RRFS real(kind=kind_phys), dimension(ims:ime, 1), intent(in) :: emi_ant_in real(kind=kind_phys), dimension(ims:ime, kms:kme), intent(in) :: pr3d,ph3d real(kind=kind_phys), dimension(ims:ime, kts:kte), intent(in) :: & @@ -681,8 +709,9 @@ subroutine rrfs_smoke_prep( & rri, t_phy, u_phy, v_phy, p_phy, rho_phy, dz8w, p8w, t8w, vvel, & zmid, pi_phy, theta_phy, wind_phy real(kind_phys), dimension(ims:ime, jms:jme), intent(out) :: & - u10, v10, ust, tsk, xland, xlat, xlong, dxy, rmol, swdown, znt, & - pbl, hfx, snowh, clayf, rdrag, sandf, ssm, uthr, hwp_local + u10, v10, ust, tsk, xland, xlat, xlong, dxy, rmol, swdown, znt, & + pbl, hfx, snowh, clayf, rdrag, sandf, ssm, uthr, hwp_local, uspdavg2d, & + hpbl2d real(kind_phys), dimension(ims:ime, nlcat, jms:jme), intent(out) :: vegfrac real(kind_phys), dimension(ims:ime, kms:kme, jms:jme, num_moist), intent(out) :: moist real(kind_phys), dimension(ims:ime, kms:kme, jms:jme, num_chem), intent(out) :: chem @@ -690,19 +719,20 @@ subroutine rrfs_smoke_prep( & real(kind_phys), dimension(ims:ime, kms:kme, jms:jme), intent(out) :: z_at_w real(kind_phys), dimension(ims:ime, nsoil, jms:jme), intent(out) :: smois,stemp real(kind_phys), dimension(ims:ime,jms:jme), intent(inout) :: frp_in, fire_end_hr, fire_hist, coef_bb_dc - real(kind_phys), dimension(ims:ime,jms:jme), intent(inout) :: hwp_day_avg, peak_hr + real(kind_phys), dimension(ims:ime,jms:jme), intent(inout) :: hwp_day_avg, totprcp_24hrs, peak_hr real(kind_phys), dimension(ims:ime), intent(inout) :: emis_anoc,ebb_smoke_in real(kind_phys), parameter :: conv_frp = 1.e+06_kind_phys ! FRP conversion factor, MW to W real(kind_phys), parameter :: frpc = 1._kind_phys ! FRP conversion factor (Regional) ! -- local variables integer i,ip,j,k,k1,kp,kk,kkp,nv,l,ll,n,nl - real(kind_phys) :: SFCWIND,WIND,DELWIND,DZ,wdgust,snoweq,THETA + real(kind_phys) :: SFCWIND,SFCWIND2,WIND,DELWIND,DZ,wdgust,snoweq,THETA real(kind_phys), dimension(ims:ime, kms:kme, jms:jme) :: THETAV real(kind_phys), dimension(ims:ime, jms:jme) :: windgustpot - real(kind_phys), dimension(ims:ime, jms:jme), intent(out) :: kpbl_thetav + integer, dimension(ims:ime, jms:jme),intent(out) :: kpbl, kpbl_thetav real(kind_phys), parameter :: delta_theta4gust = 0.5 real(kind=kind_phys),parameter :: p1000mb = 100000. + real(kind_phys) :: precip_factor ! -- initialize fire emissions ebu_in = 0._kind_phys @@ -710,7 +740,9 @@ subroutine rrfs_smoke_prep( & emis_anoc = 0._kind_phys frp_in = 0._kind_phys hwp_day_avg = 0._kind_phys + totprcp_24hrs = 0._kind_phys fire_end_hr = 0._kind_phys + uspdavg2d = 0._kind_phys ! -- initialize output arrays isltyp = 0._kind_phys @@ -752,6 +784,8 @@ subroutine rrfs_smoke_prep( & moist = 0._kind_phys chem = 0._kind_phys z_at_w = 0._kind_phys + kpbl = 1 + kpbl_thetav = 1 if ( ebb_dcycle == 1 ) then coef_bb_dc = 1._kind_phys endif @@ -783,7 +817,6 @@ subroutine rrfs_smoke_prep( & rmol (i,1)=recmol (i) enddo - do k=1,nsoil do j=jts,jte do i=its,ite @@ -808,6 +841,17 @@ subroutine rrfs_smoke_prep( & enddo enddo + do j=jts,jte + do i=its,ite + do k=kts+1,kte + if(z_at_w(i,k,j).gt.pbl(i,j))then + kpbl(i,j)=max(2,k) + exit + endif + enddo + enddo + enddo + do j=jts,jte do k=kts,kte+1 do i=its,ite @@ -892,52 +936,93 @@ subroutine rrfs_smoke_prep( & enddo enddo -!---- Calculate wind gust potential and HWP +!---- Calculate wind gust potential and average boundary layer wind do i = its,ite SFCWIND = sqrt(u10m(i)**2+v10m(i)**2) windgustpot(i,1) = SFCWIND - if (kpbl_thetav(i,1)+1 .ge. kts+1 ) then - do k=kts+1,int(kpbl_thetav(i,1))+1 + uspdavg2d(i,1) = SFCWIND + if (kpbl(i,1)+1 .ge. kts+1 ) then + do k=kts+1,kpbl(i,1)+1 ! Use kpbl from MYNN WIND = sqrt(us3d(i,k)**2+vs3d(i,k)**2) + uspdavg2d(i,1) = uspdavg2d(i,1) + WIND DELWIND = WIND - SFCWIND DZ = zmid(i,k,1) - oro(i) DELWIND = DELWIND*(1.0-MIN(0.5,DZ/2000.)) windgustpot(i,1) = max(windgustpot(i,1),SFCWIND+DELWIND) enddo endif + uspdavg2d(i,1) = uspdavg2d(i,1) / real(kpbl(i,1)) +! JLS - we have pbl height from MYNN (=pbl), should hpbl2d be renamed to +! pbl_thetav - then whcih should be used for HWP and whcih should be passed to +! plumerise? + hpbl2d(i,1) = z_at_w(i,kpbl(i,1),1) - z_at_w(i,kts,1) ! From MYNN enddo - hwp_local = 0. + +!---- 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 + ! 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 + if ( ebb_dcycle == 2 ) then + do i=its, ite + do j=jts, jte + totprcp_24hrs (i,j) = smoke2d_RRFS(i,5) + enddo + enddo + endif do i=its,ite - wdgust=max(windgustpot(i,1),3.) - snoweq=max((25.-snow_cpl(i))/25.,0.) - hwp_local(i,1)=0.177*wdgust**0.97*max(t2m(i)-dpt2m(i),15.)**1.03*((1.-wetness(i))**0.4)*snoweq ! Eric update 11/2023 + SFCWIND2=max(sqrt(u10m(i)**2+v10m(i)**2),3._kind_phys) + SELECT CASE (hwp_method) + CASE (1) ! Operational method - includes accumulated precip + hwp_local(i,1)=0.022_kind_phys*MAX(precip_factor-(totprcp(i)+totprcp_24hrs(i,1))*1.e+3_kind_phys,0._kind_phys)/precip_factor * & + ((1._kind_phys-wetness(i))**0.51_kind_phys) * & + (SFCWIND2*hpbl2d(i,1))**0.57 * & + MIN(25.0_kind_phys,MAX(15._kind_phys,t2m(i)-dpt2m(i)))**0.74 * & + MIN(3._kind_phys, 1._kind_phys + dswsfc(i)/250._kind_phys)**0.18 !+ 28.67_kind_phys ! Eric update 01/2024 + CASE (2) ! Pre-release of RRFSv1 method - using wind gust calculated via UPP Method + wdgust =max(windgustpot(i,1),3._kind_phys) + snoweq =max((25._kind_phys-snow_cpl(i))/25._kind_phys,0._kind_phys) + hwp_local(i,1)=0.177_kind_phys*wdgust**0.97_kind_phys*max(t2m(i)-dpt2m(i),15._kind_phys)**1.03_kind_phys * & + ((1._kind_phys-wetness(i))**0.4_kind_phys)*snoweq ! Eric update 11/2023 + CASE (3) ! Modified operational method - vent coef calculated using average PBL wind speed + hwp_local(i,1)=0.022_kind_phys*MAX(precip_factor-(totprcp(i)+totprcp_24hrs(i,1))*1.e+3_kind_phys,0._kind_phys)/precip_factor * & + ((1._kind_phys-wetness(i))**0.51_kind_phys) * & + (uspdavg2d(i,1)*hpbl2d(i,1))**0.57 * & + MIN(25.0_kind_phys,MAX(15._kind_phys,t2m(i)-dpt2m(i)))**0.74 * & + MIN(3._kind_phys, 1._kind_phys + dswsfc(i)/250._kind_phys)**0.18 !+ 28.67_kind_phys ! Eric update 01/2024 + CASE (4) ! Modified Pre-release of RRFSv1 methood - using wind gust calculated as scaled surface wind + wdgust =max(1.69*sqrt(u10m(i)**2+v10m(i)**2), 3._kind_phys) + snoweq =max((25._kind_phys-snow_cpl(i))/25._kind_phys,0._kind_phys) + hwp_local(i,1)=0.177_kind_phys*wdgust**0.97_kind_phys*max(t2m(i)-dpt2m(i),15._kind_phys)**1.03_kind_phys * & + ((1._kind_phys-wetness(i))**0.4_kind_phys)*snoweq ! Eric update 11/2023 + CASE DEFAULT + END SELECT + enddo -! Set paramters for ebb_dcycle option + + ! Set paramters for ebb_dcycle option if (ebb_dcycle == 1 ) then if (hour_int .le. 24) then do j=jts,jte do i=its,ite ebu_in (i,j) = smoke_RRFS(i,hour_int+1,1) ! smoke frp_in (i,j) = smoke_RRFS(i,hour_int+1,2)*conv_frp ! frp - ! These 2 arrays aren't needed for this option - ! fire_end_hr(i,j) = 0.0 - ! hwp_day_avg(i,j) = 0.0 ebb_smoke_in (i) = ebu_in(i,j) enddo enddo endif endif - ! RAR: here we need to initialize various arrays in order to apply HWP to - ! diurnal cycle + ! Here we need to initialize various arrays in order to apply HWP to diurnal cycle ! if ebb_dcycle/=2 then those arrays=0, we need to read in temporal if (ebb_dcycle == 2) then do i=its, ite do j=jts, jte - ebu_in (i,j) = smoke2d_RRFS(i,1)!/86400. - frp_in (i,j) = smoke2d_RRFS(i,2)*conv_frp - fire_end_hr (i,j) = smoke2d_RRFS(i,3) - hwp_day_avg (i,j) = smoke2d_RRFS(i,4) - ebb_smoke_in(i ) = ebu_in(i,j) + ebu_in (i,j) = smoke2d_RRFS(i,1)!/86400. + frp_in (i,j) = smoke2d_RRFS(i,2)*conv_frp + fire_end_hr (i,j) = smoke2d_RRFS(i,3) + hwp_day_avg (i,j) = smoke2d_RRFS(i,4) + ebb_smoke_in (i ) = ebu_in(i,j) enddo enddo end if @@ -945,7 +1030,6 @@ subroutine rrfs_smoke_prep( & if (ktau==1) then do j=jts,jte do i=its,ite - ! GFS_typedefs.F90 initializes this = 1, but should be OK to duplicate, RAR?? fire_hist (i,j) = 1. coef_bb_dc (i,j) = 1. if (xlong(i,j)<230.) then @@ -965,7 +1049,6 @@ subroutine rrfs_smoke_prep( & enddo endif - ! We will add a namelist variable, real :: flam_frac_global, RAR?? do k=kms,kte do i=ims,ime chem(i,k,jts,p_smoke )=max(epsilc,gq0(i,k,ntsmoke )) diff --git a/physics/smoke_dust/rrfs_smoke_wrapper.meta b/physics/smoke_dust/rrfs_smoke_wrapper.meta index 739a43d70..d61ce5d7e 100755 --- a/physics/smoke_dust/rrfs_smoke_wrapper.meta +++ b/physics/smoke_dust/rrfs_smoke_wrapper.meta @@ -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,13 +64,6 @@ 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 @@ -85,6 +71,13 @@ 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,14 @@ 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 [uspdavg] standard_name = mean_wind_speed_in_boundary_layer long_name = average wind speed within the boundary layer @@ -944,14 +967,6 @@ type = real 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 +984,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) From 5048dbbbb2ccab5923ff1295fc32b47887fc5828 Mon Sep 17 00:00:00 2001 From: "Haiqin.Li" Date: Wed, 13 Mar 2024 23:36:33 +0000 Subject: [PATCH 3/5] "update to address UFS reviewer's comments" --- physics/smoke_dust/module_add_emiss_burn.F90 | 2 +- physics/smoke_dust/rrfs_smoke_wrapper.F90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/physics/smoke_dust/module_add_emiss_burn.F90 b/physics/smoke_dust/module_add_emiss_burn.F90 index 889e42b43..50a56c1bc 100755 --- a/physics/smoke_dust/module_add_emiss_burn.F90 +++ b/physics/smoke_dust/module_add_emiss_burn.F90 @@ -44,7 +44,7 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi,ebb_min, & !>--local logical, intent(in) :: add_fire_moist_flux integer :: i,j,k,n,m - integer :: icall !=0 + 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 diff --git a/physics/smoke_dust/rrfs_smoke_wrapper.F90 b/physics/smoke_dust/rrfs_smoke_wrapper.F90 index 87817fb7d..32e1c6768 100755 --- a/physics/smoke_dust/rrfs_smoke_wrapper.F90 +++ b/physics/smoke_dust/rrfs_smoke_wrapper.F90 @@ -1,4 +1,4 @@ -!>\file rrfs_smoke_wrapper.F90hwp_method +!>\file rrfs_smoke_wrapper.F90 !! This file is CCPP driver of RRFS Smoke and Dust !! Haiqin.Li@noaa.gov 03/2024 From dac2f23faaa9d19609c9521ee6ca8d670b485002 Mon Sep 17 00:00:00 2001 From: "Haiqin.Li" Date: Thu, 14 Mar 2024 14:45:35 +0000 Subject: [PATCH 4/5] "merge physics PR #186 from Jili" --- physics/Interstitials/UFS_SCM_NEPTUNE/sgscloud_radpre.F90 | 7 ++----- physics/Interstitials/UFS_SCM_NEPTUNE/sgscloud_radpre.meta | 7 +++++++ 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/sgscloud_radpre.F90 b/physics/Interstitials/UFS_SCM_NEPTUNE/sgscloud_radpre.F90 index 95d172d4f..b84cb612a 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/sgscloud_radpre.F90 +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/sgscloud_radpre.F90 @@ -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 diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/sgscloud_radpre.meta b/physics/Interstitials/UFS_SCM_NEPTUNE/sgscloud_radpre.meta index 8e25428cc..813469cd2 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/sgscloud_radpre.meta +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/sgscloud_radpre.meta @@ -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 From ebc3fb310429ba6017c4d27b297f80c8072bc110 Mon Sep 17 00:00:00 2001 From: Grant Firl Date: Fri, 8 Nov 2024 21:17:35 -0500 Subject: [PATCH 5/5] add optional to two variables in rrfs_smoke_wrapper --- physics/smoke_dust/rrfs_smoke_wrapper.meta | 2 ++ 1 file changed, 2 insertions(+) diff --git a/physics/smoke_dust/rrfs_smoke_wrapper.meta b/physics/smoke_dust/rrfs_smoke_wrapper.meta index d61ce5d7e..152686947 100755 --- a/physics/smoke_dust/rrfs_smoke_wrapper.meta +++ b/physics/smoke_dust/rrfs_smoke_wrapper.meta @@ -844,6 +844,7 @@ 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 @@ -967,6 +968,7 @@ type = real kind = kind_phys intent = out + optional = True [oro] standard_name = height_above_mean_sea_level long_name = height_above_mean_sea_level