From 274a69673ed190a0f3c5e0d48e6643b0f0af1580 Mon Sep 17 00:00:00 2001 From: "Haiqin.Li" Date: Tue, 5 Mar 2024 03:32:42 +0000 Subject: [PATCH 01/15] "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 79b9522c5..a0f647083 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) @@ -1932,11 +1932,11 @@ SUBROUTINE mym_length ( & h1=MIN(h1,maxdz) ! 1/2 transition layer depth h2=h1/2.0 ! 1/4 transition layer depth - qkw(kts) = SQRT(MAX(qke(kts),1.0e-10)) + qkw(kts) = SQRT(MAX(qke(kts), qkemin)) DO k = kts+1,kte afk = dz(k)/( dz(k)+dz(k-1) ) abk = 1.0 -afk - qkw(k) = SQRT(MAX(qke(k)*abk+qke(k-1)*afk,1.0e-3)) + qkw(k) = SQRT(MAX(qke(k)*abk+qke(k-1)*afk, qkemin)) END DO elt = 1.0e-5 @@ -1956,7 +1956,7 @@ SUBROUTINE mym_length ( & elt = alp1*elt/vsc vflx = ( vt(kts)+1.0 )*flt +( vq(kts)+tv0 )*flq - vsc = ( gtr*elt*MAX( vflx, 0.0 ) )**(1.0/3.0) + vsc = ( gtr*elt*MAX( vflx, 0.0 ) )**onethird ! ** Strictly, el(i,k=1) is not zero. ** el(kts) = 0.0 @@ -2014,14 +2014,14 @@ SUBROUTINE mym_length ( & h1=MIN(h1,600.) ! 1/2 transition layer depth h2=h1/2.0 ! 1/4 transition layer depth - qtke(kts)=MAX(0.5*qke(kts), 0.01) !tke at full sigma levels + qtke(kts)=MAX(0.5*qke(kts), 0.5*qkemin) !tke at full sigma levels thetaw(kts)=theta(kts) !theta at full-sigma levels - qkw(kts) = SQRT(MAX(qke(kts),1.0e-10)) + qkw(kts) = SQRT(MAX(qke(kts), qkemin)) DO k = kts+1,kte afk = dz(k)/( dz(k)+dz(k-1) ) abk = 1.0 -afk - qkw(k) = SQRT(MAX(qke(k)*abk+qke(k-1)*afk,1.0e-3)) + qkw(k) = SQRT(MAX(qke(k)*abk+qke(k-1)*afk, qkemin)) qtke(k) = 0.5*(qkw(k)**2) ! q -> TKE thetaw(k)= theta(k)*abk + theta(k-1)*afk END DO @@ -2034,14 +2034,14 @@ SUBROUTINE mym_length ( & zwk = zw(k) DO WHILE (zwk .LE. zi2+h1) dzk = 0.5*( dz(k)+dz(k-1) ) - qdz = min(max( qkw(k)-qmin, 0.02 ), 30.0)*dzk + qdz = min(max( qkw(k)-qmin, 0.01 ), 30.0)*dzk elt = elt +qdz*zwk vsc = vsc +qdz k = k+1 zwk = zw(k) END DO - elt = MIN( MAX( alp1*elt/vsc, 10.), 400.) + elt = MIN( MAX( alp1*elt/vsc, 8.), 400.) !avoid use of buoyancy flux functions which are ill-defined at the surface !vflx = ( vt(kts)+1.0 )*flt + ( vq(kts)+tv0 )*flq vflx = fltv @@ -2117,13 +2117,13 @@ SUBROUTINE mym_length ( & h1=MIN(h1,600.) h2=h1*0.5 ! 1/4 transition layer depth - qtke(kts)=MAX(0.5*qke(kts),0.01) !tke at full sigma levels - qkw(kts) = SQRT(MAX(qke(kts),1.0e-4)) + qtke(kts)=MAX(0.5*qke(kts), 0.5*qkemin) !tke at full sigma levels + qkw(kts) = SQRT(MAX(qke(kts), qkemin)) DO k = kts+1,kte afk = dz(k)/( dz(k)+dz(k-1) ) abk = 1.0 -afk - qkw(k) = SQRT(MAX(qke(k)*abk+qke(k-1)*afk,1.0e-3)) + qkw(k) = SQRT(MAX(qke(k)*abk+qke(k-1)*afk, qkemin)) qtke(k) = 0.5*qkw(k)**2 ! qkw -> TKE END DO @@ -3356,8 +3356,8 @@ SUBROUTINE mym_predict (kts,kte, & CALL tridiag2(kte,a,b,c,d,x) DO k=kts,kte -! qke(k)=max(d(k-kts+1), 1.e-4) - qke(k)=max(x(k), 1.e-4) +! qke(k)=max(d(k-kts+1), qkemin) + qke(k)=max(x(k), qkemin) qke(k)=min(qke(k), 150.) ENDDO @@ -6504,11 +6504,11 @@ SUBROUTINE DMP_mf( & do k=kts,kte-1 do I=1,nup edmf_a(K) =edmf_a(K) +UPA(K,i) - edmf_w(K) =edmf_w(K) +rhoz(k)*UPA(K,i)*UPW(K,i) - edmf_qt(K) =edmf_qt(K) +rhoz(k)*UPA(K,i)*UPQT(K,i) - edmf_thl(K)=edmf_thl(K)+rhoz(k)*UPA(K,i)*UPTHL(K,i) - edmf_ent(K)=edmf_ent(K)+rhoz(k)*UPA(K,i)*ENT(K,i) - edmf_qc(K) =edmf_qc(K) +rhoz(k)*UPA(K,i)*UPQC(K,i) + edmf_w(K) =edmf_w(K) +UPA(K,i)*UPW(K,i) + edmf_qt(K) =edmf_qt(K) +UPA(K,i)*UPQT(K,i) + edmf_thl(K)=edmf_thl(K)+UPA(K,i)*UPTHL(K,i) + edmf_ent(K)=edmf_ent(K)+UPA(K,i)*ENT(K,i) + edmf_qc(K) =edmf_qc(K) +UPA(K,i)*UPQC(K,i) enddo enddo do k=kts,kte-1 From 7eab4efc48f27d662eb372294361431aefae7900 Mon Sep 17 00:00:00 2001 From: "Haiqin.Li" Date: Wed, 13 Mar 2024 03:44:24 +0000 Subject: [PATCH 02/15] "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 | 274 ++++++++++++------ physics/smoke_dust/rrfs_smoke_wrapper.meta | 70 +++-- 6 files changed, 369 insertions(+), 434 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 145b23934..f4e533284 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 @@ -8,11 +8,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 @@ -30,8 +31,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 @@ -42,31 +41,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 @@ -92,9 +91,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 @@ -123,12 +124,12 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, ebu_smoke,fhist,min_fplume, & max_fplume, hwp, hwp_ave, wetness, ndvel, ddvel_inout, & 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, ntsmoke, ntdust, ntcoarsepm, ndvel, nlcat real(kind_phys),intent(in) :: dt, julian, g, pi, con_cp, con_rd, con_fv @@ -145,7 +146,7 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, real(kind_phys), dimension(:,:), intent(in) :: 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 @@ -154,6 +155,7 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, real(kind_phys), dimension(:), intent(inout) :: emdust, emseas, emanoc real(kind_phys), dimension(:), intent(inout) :: ebb_smoke_in,coef_bb, frp_output, fhist real(kind_phys), dimension(:,:), intent(inout) :: ebu_smoke + real(kind_phys), dimension(:,:), intent(inout) :: rho_dry real(kind_phys), dimension(:), intent(out ) :: fire_heat_flux_out, frac_grid_burned_out real(kind_phys), dimension(:), intent(inout) :: max_fplume, min_fplume, uspdavg, hpbl_thetav real(kind_phys), dimension(:), intent(inout) :: hwp, peak_hr_out @@ -166,7 +168,6 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, real(kind_phys), dimension(:), intent(out) :: lu_nofire_out,lu_qfire_out integer, dimension(:), intent(out) :: fire_type_out 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 @@ -196,11 +197,12 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, 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 @@ -215,6 +217,7 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, !> -- 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 @@ -229,6 +232,14 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, 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 = '' @@ -244,13 +255,14 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, 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. @@ -298,16 +310,16 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, 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, & @@ -339,20 +351,26 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, 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 @@ -390,11 +408,11 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, ! 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 @@ -406,7 +424,7 @@ 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) = 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 @@ -417,21 +435,24 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, 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 ) @@ -506,9 +527,9 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, !---- 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 @@ -552,6 +573,13 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, 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 @@ -559,17 +587,15 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, 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 @@ -607,16 +633,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, & @@ -629,19 +655,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) :: & @@ -663,8 +690,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 @@ -672,19 +700,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 @@ -692,7 +721,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 @@ -734,6 +765,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 @@ -765,7 +798,6 @@ subroutine rrfs_smoke_prep( & rmol (i,1)=recmol (i) enddo - do k=1,nsoil do j=jts,jte do i=its,ite @@ -790,6 +822,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 @@ -874,52 +917,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 @@ -927,7 +1011,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 @@ -947,7 +1030,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 271d2dd36..2d793875f 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 @@ -191,6 +199,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] @@ -589,7 +604,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 @@ -790,6 +805,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 @@ -901,13 +924,6 @@ type = real kind = kind_phys intent = out -[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 @@ -916,6 +932,14 @@ type = real kind = kind_phys intent = in +[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 [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP From 766360bf7ec501d436bec92bf80f4a449f151114 Mon Sep 17 00:00:00 2001 From: "Haiqin.Li" Date: Wed, 13 Mar 2024 23:36:33 +0000 Subject: [PATCH 03/15] "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 f4e533284..0bff3fbde 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 0935794e948a5501b1d45d1b050b21e23399bfc2 Mon Sep 17 00:00:00 2001 From: "Haiqin.Li" Date: Thu, 14 Mar 2024 14:45:35 +0000 Subject: [PATCH 04/15] "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 936393d5b..9c2e27611 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 a9635efa5..e094c3e12 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/sgscloud_radpre.meta +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/sgscloud_radpre.meta @@ -273,6 +273,13 @@ dimensions = () type = integer intent = in +[conv_cf_opt] + standard_name = option_for_convection_scheme_cloud_fraction_computation + long_name = option for convection scheme cloud fraction computation + units = flag + dimensions = () + type = integer + intent = in [qc_save] standard_name = cloud_condensed_water_mixing_ratio_save long_name = ratio of mass of cloud water to mass of dry air plus vapor (without condensates) before entering a physics scheme From 3dd1dd34d6f94b589107ed43b57a994277e83b58 Mon Sep 17 00:00:00 2001 From: "jordan.schnell" Date: Thu, 28 Mar 2024 16:03:36 +0000 Subject: [PATCH 05/15] Fix improperly assigned fire emissions for ebb_dcycle==1 for retrospectives (NOT operational!) --- physics/smoke_dust/module_add_emiss_burn.F90 | 9 --------- 1 file changed, 9 deletions(-) diff --git a/physics/smoke_dust/module_add_emiss_burn.F90 b/physics/smoke_dust/module_add_emiss_burn.F90 index 50a56c1bc..ff1fb09b6 100755 --- a/physics/smoke_dust/module_add_emiss_burn.F90 +++ b/physics/smoke_dust/module_add_emiss_burn.F90 @@ -66,15 +66,6 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi,ebb_min, & timeq= mod(timeq,timeq_max) coef_con=1._kind_phys/((2._kind_phys*pi)**0.5) -! 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) - enddo - enddo - endif - if (ebb_dcycle==2) then do j=jts,jte From 1b089ec49ff3dc110f244b2fa89a7b2fba39c528 Mon Sep 17 00:00:00 2001 From: Dustin Swales Date: Tue, 23 Apr 2024 20:22:56 +0000 Subject: [PATCH 06/15] Add SW clear-sky downward flux at surface to available diagnostics --- .../UFS_SCM_NEPTUNE/GFS_debug.F90 | 1 + .../Interstitials/UFS_SCM_NEPTUNE/dcyc2t3.f | 18 +++++++++++------- .../Interstitials/UFS_SCM_NEPTUNE/dcyc2t3.meta | 16 ++++++++++++++++ physics/Radiation/RRTMG/rrtmg_sw_post.F90 | 6 ++++-- physics/Radiation/RRTMG/rrtmg_sw_post.meta | 8 ++++++++ 5 files changed, 40 insertions(+), 9 deletions(-) diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_debug.F90 b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_debug.F90 index ed26b795f..7f85da7f6 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_debug.F90 +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_debug.F90 @@ -688,6 +688,7 @@ subroutine GFS_diagtoscreen_run (Model, Statein, Stateout, Sfcprop, Coupling, call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Diag%dlwsfci ', Diag%dlwsfci) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Diag%ulwsfci ', Diag%ulwsfci) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Diag%dswsfci ', Diag%dswsfci) + call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Diag%dswsfcci ', Diag%dswsfcci) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Diag%nswsfci ', Diag%nswsfci) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Diag%uswsfci ', Diag%uswsfci) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Diag%dusfci ', Diag%dusfci) diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/dcyc2t3.f b/physics/Interstitials/UFS_SCM_NEPTUNE/dcyc2t3.f index 36299651d..d8ef0a86d 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/dcyc2t3.f +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/dcyc2t3.f @@ -36,7 +36,7 @@ module dcyc2t3 ! ( solhr,slag,sdec,cdec,sinlat,coslat, ! ! xlon,coszen,tsfc_lnd,tsfc_ice,tsfc_wat, ! ! tf,tsflw,sfcemis_lnd,sfcemis_ice,sfcemis_wat, ! -! sfcdsw,sfcnsw,sfcdlw,sfculw,swh,swhc,hlw,hlwc, ! +! sfcdsw,sfcdswc,sfcnsw,sfcdlw,sfculw,swh,swhc,hlw,hlwc, ! ! sfcnirbmu,sfcnirdfu,sfcvisbmu,sfcvisdfu, ! ! sfcnirbmd,sfcnirdfd,sfcvisbmd,sfcvisdfd, ! ! im, levs, deltim, fhswr, ! @@ -44,7 +44,7 @@ module dcyc2t3 ! input/output: ! ! dtdt,dtdtnp, ! ! outputs: ! -! adjsfcdsw,adjsfcnsw,adjsfcdlw,adjsfculw, ! +! adjsfcdsw,adjsfcdswc,adjsfcnsw,adjsfcdlw, ! ! adjsfculw_lnd,adjsfculw_ice,adjsfculw_wat,xmu,xcosz, ! ! adjnirbmu,adjnirdfu,adjvisbmu,adjvisdfu, ! ! adjdnnbmd,adjdnndfd,adjdnvbmd,adjdnvdfd) ! @@ -68,6 +68,7 @@ module dcyc2t3 ! sfcemis_wat(im) - real, surface emissivity (fraction) o. ocean (k)! ! tsflw (im) - real, sfc air (layer 1) temp in k saved in lw call ! ! sfcdsw (im) - real, total sky sfc downward sw flux ( w/m**2 ) ! +! sfcdswc (im) - real, clear sky sfc downward sw flux ( w/m**2 ) ! ! sfcnsw (im) - real, total sky sfc net sw into ground (w/m**2) ! ! sfcdlw (im) - real, total sky sfc downward lw flux ( w/m**2 ) ! ! sfculw (im) - real, total sky sfc upward lw flux ( w/m**2 ) ! @@ -98,6 +99,7 @@ module dcyc2t3 ! ! ! outputs: ! ! adjsfcdsw(im)- real, time step adjusted sfc dn sw flux (w/m**2) ! +! adjsfcdswc(im)- real, time step adjusted sfc dn sw flux (w/m**2) ! ! adjsfcnsw(im)- real, time step adj sfc net sw into ground (w/m**2)! ! adjsfcdlw(im)- real, time step adjusted sfc dn lw flux (w/m**2) ! ! adjsfculw_lnd(im)- real, sfc upw. lw flux at current time (w/m**2)! @@ -169,7 +171,7 @@ subroutine dcyc2t3_run & & con_g, con_cp, con_pi, con_sbc, & & xlon,coszen,tsfc_lnd,tsfc_ice,tsfc_wat,tf,tsflw,tsfc, & & sfcemis_lnd, sfcemis_ice, sfcemis_wat, & - & sfcdsw,sfcnsw,sfcdlw,swh,swhc,hlw,hlwc, & + & sfcdsw,sfcdswc,sfcnsw,sfcdlw,swh,swhc,hlw,hlwc, & & sfcnirbmu,sfcnirdfu,sfcvisbmu,sfcvisdfu, & & sfcnirbmd,sfcnirdfd,sfcvisbmd,sfcvisdfd, & & im, levs, deltim, fhswr, & @@ -181,7 +183,7 @@ subroutine dcyc2t3_run & ! --- input/output: & dtdt,dtdtnp,htrlw, & ! --- outputs: - & adjsfcdsw,adjsfcnsw,adjsfcdlw,adjsfculw, & + & adjsfcdsw,adjsfcdswc,adjsfcnsw,adjsfculw,adjsfcdlw, & & adjsfculw_lnd,adjsfculw_ice,adjsfculw_wat,xmu,xcosz, & & adjnirbmu,adjnirdfu,adjvisbmu,adjvisdfu, & & adjnirbmd,adjnirdfd,adjvisbmd,adjvisdfd, & @@ -213,8 +215,9 @@ subroutine dcyc2t3_run & & deltim, fhswr, lfnc_k, lfnc_p0 real(kind=kind_phys), dimension(:), intent(in) :: & - & sinlat, coslat, xlon, coszen, tf, tsflw, sfcdlw, & - & sfcdsw, sfcnsw, sfculw, sfculw_med, tsfc, tsfc_radtime + & sinlat, coslat, xlon, coszen, tf, tsflw, sfcdlw, & + & sfcdsw, sfcdswc, sfcnsw, sfculw, sfculw_med, tsfc, & + & tsfc_radtime real(kind=kind_phys), dimension(:), intent(in) :: & & tsfc_lnd, tsfc_ice, tsfc_wat, & @@ -244,7 +247,7 @@ subroutine dcyc2t3_run & real(kind=kind_phys), dimension(:), intent(out) :: & & adjsfcdsw, adjsfcnsw, adjsfcdlw, adjsfculw, xmu, xcosz, & & adjnirbmu, adjnirdfu, adjvisbmu, adjvisdfu, & - & adjnirbmd, adjnirdfd, adjvisbmd, adjvisdfd + & adjnirbmd, adjnirdfd, adjvisbmd, adjvisdfd, adjsfcdswc real(kind=kind_phys), dimension(:), intent(out) :: & & adjsfculw_lnd, adjsfculw_ice, adjsfculw_wat @@ -367,6 +370,7 @@ subroutine dcyc2t3_run & adjsfcnsw(i) = sfcnsw(i) * xmu(i) adjsfcdsw(i) = sfcdsw(i) * xmu(i) + adjsfcdswc(i)= sfcdswc(i) * xmu(i) adjnirbmu(i) = sfcnirbmu(i) * xmu(i) adjnirdfu(i) = sfcnirdfu(i) * xmu(i) diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/dcyc2t3.meta b/physics/Interstitials/UFS_SCM_NEPTUNE/dcyc2t3.meta index 95b3f341b..386b0ee0f 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/dcyc2t3.meta +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/dcyc2t3.meta @@ -191,6 +191,14 @@ type = real kind = kind_phys intent = in +[sfcdswc] + standard_name = surface_downwelling_shortwave_flux_on_radiation_timestep_assuming_clear_sky + long_name = clear sky surface downwelling shortwave flux on radiation time step + units = W m-2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = in [sfcdlw] standard_name = surface_downwelling_longwave_flux_on_radiation_timestep long_name = total sky surface downwelling longwave flux on radiation time step @@ -508,6 +516,14 @@ type = real kind = kind_phys intent = out +[adjsfcdswc] + standard_name = surface_downwelling_shortwave_flux_assuming_clear_sky + long_name = surface downwelling shortwave flux at current time assuming clear sky + units = W m-2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out [adjsfcnsw] standard_name = surface_net_downwelling_shortwave_flux long_name = surface net downwelling shortwave flux at current time diff --git a/physics/Radiation/RRTMG/rrtmg_sw_post.F90 b/physics/Radiation/RRTMG/rrtmg_sw_post.F90 index f39cba71c..8fa9494b8 100644 --- a/physics/Radiation/RRTMG/rrtmg_sw_post.F90 +++ b/physics/Radiation/RRTMG/rrtmg_sw_post.F90 @@ -14,7 +14,7 @@ subroutine rrtmg_sw_post_run (im, levr, levs, ltp, nday, lm, kd, lsswr, & swhtr, sfcalb1, sfcalb2, sfcalb3, sfcalb4, htswc, htsw0, & nirbmdi, nirdfdi, visbmdi, visdfdi, nirbmui, nirdfui, visbmui,& visdfui, sfcdsw, sfcnsw, htrsw, swhc, scmpsw, sfcfsw, topfsw, & - errmsg, errflg) + sfcdswc, errmsg, errflg) use machine, only: kind_phys use module_radsw_parameters, only: topfsw_type, sfcfsw_type, & @@ -33,7 +33,8 @@ subroutine rrtmg_sw_post_run (im, levr, levs, ltp, nday, lm, kd, lsswr, & visbmdi, visdfdi, & nirbmui, nirdfui, & visbmui, visdfui, & - sfcdsw, sfcnsw + sfcdsw, sfcnsw, & + sfcdswc real(kind=kind_phys), dimension(:,:), intent(inout) :: htrsw, swhc type(cmpfsw_type), dimension(:), intent(inout) :: scmpsw @@ -122,6 +123,7 @@ subroutine rrtmg_sw_post_run (im, levr, levs, ltp, nday, lm, kd, lsswr, & do i=1,im sfcnsw(i) = sfcfsw(i)%dnfxc - sfcfsw(i)%upfxc sfcdsw(i) = sfcfsw(i)%dnfxc + sfcdswc(i)= sfcfsw(i)%dnfx0 enddo endif ! end_if_lsswr diff --git a/physics/Radiation/RRTMG/rrtmg_sw_post.meta b/physics/Radiation/RRTMG/rrtmg_sw_post.meta index 9914051ce..0a48b04d7 100644 --- a/physics/Radiation/RRTMG/rrtmg_sw_post.meta +++ b/physics/Radiation/RRTMG/rrtmg_sw_post.meta @@ -198,6 +198,14 @@ type = real kind = kind_phys intent = inout +[sfcdswc] + standard_name = surface_downwelling_shortwave_flux_on_radiation_timestep_assuming_clear_sky + long_name = clear sky sfc downward sw flux + units = W m-2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout [htrsw] standard_name = tendency_of_air_temperature_due_to_shortwave_heating_on_radiation_timestep long_name = total sky sw heating rate From 7d1a68987b520e4421359cec7687949dc7a75a45 Mon Sep 17 00:00:00 2001 From: "Haiqin.Li" Date: Wed, 26 Jun 2024 04:23:18 +0000 Subject: [PATCH 07/15] "update deep soil temperature of RUC LSM" --- physics/SFC_Models/Land/RUC/module_sf_ruclsm.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) 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 From 7a28e0bfaf4b055662f18f3352ca38a7d27beb68 Mon Sep 17 00:00:00 2001 From: "Haiqin.Li" Date: Tue, 23 Jul 2024 21:11:40 +0000 Subject: [PATCH 08/15] "Update to scale the RAVE fire emission and duirnal cycle of agriculute fire" --- physics/smoke_dust/module_add_emiss_burn.F90 | 46 ++++++++++---------- physics/smoke_dust/rrfs_smoke_config.F90 | 1 + physics/smoke_dust/rrfs_smoke_wrapper.F90 | 34 ++++++++++----- physics/smoke_dust/rrfs_smoke_wrapper.meta | 14 ++++-- 4 files changed, 58 insertions(+), 37 deletions(-) diff --git a/physics/smoke_dust/module_add_emiss_burn.F90 b/physics/smoke_dust/module_add_emiss_burn.F90 index ff1fb09b6..cf942a4c6 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,17 @@ 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 + real(kind_phys), parameter :: avgx1=0., sigmx1=2.2, C1=0.2 ! Ag fires + real(kind_phys), parameter :: avgx2=0.5, sigmx2=0.8, C2=1.1 ! Grass fires, slash burns timeq= gmt*3600._kind_phys + real(time_int,4) timeq= mod(timeq,timeq_max) @@ -70,34 +75,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)-2.0)) !One hour delay is due to the latency of the RAVE files, hours; one more hour subtracted to have fire_end_hr in the range of 0-24 instead of 0-25 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) ! 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 (3,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 +115,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 - 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) + ! 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 +148,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..750b5c701 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)>260. .AND. xlat(i,j)>25. .AND. xlat(i,j)<41.) then + fire_type(i,j) = 2 ! slash burn and wildfires in the east, eastern temperate forest ecosystem + else if (lu_sfire(i,j)>0.8) then + fire_type(i,j) = 3 ! savanna and grassland fires + 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 ) @@ -941,7 +951,7 @@ subroutine rrfs_smoke_prep( & !---- 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 + precip_factor = 2.5_kind_phys + real(hour_int, kind=kind_phys)*2.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 diff --git a/physics/smoke_dust/rrfs_smoke_wrapper.meta b/physics/smoke_dust/rrfs_smoke_wrapper.meta index 2d793875f..d5afb1831 100755 --- a/physics/smoke_dust/rrfs_smoke_wrapper.meta +++ b/physics/smoke_dust/rrfs_smoke_wrapper.meta @@ -93,9 +93,9 @@ type = integer intent = in [hwp_method_in] - standard_name = do_smoke_forecast - long_name = index for rrfs smoke forecast - units = index + standard_name = control_for_HWP_equation + long_name = control for HWP equation + units = 1 dimensions = () type = integer intent = in @@ -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 = 1 + dimensions = () + type = real + kind = kind_phys + intent = in [plume_alpha_in] standard_name = alpha_for_plumerise_scheme long_name = alpha paramter for plumerise scheme From ceb35f311e2dd5634c3ed64a9c7ba529ab5502cf Mon Sep 17 00:00:00 2001 From: Jili Dong Date: Thu, 22 Aug 2024 15:13:01 +0000 Subject: [PATCH 09/15] add a new parameter to control if G-F cold starts or not --- physics/CONV/Grell_Freitas/cu_gf_driver.F90 | 6 +++--- physics/CONV/Grell_Freitas/cu_gf_driver.meta | 7 +++++++ physics/CONV/Grell_Freitas/cu_gf_driver_pre.F90 | 5 +++-- physics/CONV/Grell_Freitas/cu_gf_driver_pre.meta | 7 +++++++ 4 files changed, 20 insertions(+), 5 deletions(-) diff --git a/physics/CONV/Grell_Freitas/cu_gf_driver.F90 b/physics/CONV/Grell_Freitas/cu_gf_driver.F90 index 31161f818..efef60293 100644 --- a/physics/CONV/Grell_Freitas/cu_gf_driver.F90 +++ b/physics/CONV/Grell_Freitas/cu_gf_driver.F90 @@ -56,7 +56,7 @@ end subroutine cu_gf_driver_init !! \htmlinclude cu_gf_driver_run.html !! !>\section gen_gf_driver Grell-Freitas Cumulus Scheme Driver General Algorithm - subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,& + subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart, gf_coldstart, & cactiv,cactiv_m,g,cp,xlv,r_v,forcet,forceqv_spechum,phil,raincv, & qv_spechum,t,cld1d,us,vs,t2di,w,qv2di_spechum,p2di,psuri, & hbot,htop,kcnv,xland,hfx2,qfx2,aod_gf,cliw,clcw, & @@ -97,7 +97,7 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,& integer :: its,ite, jts,jte, kts,kte integer, intent(in ) :: im,km,ntracer,nchem,kdt integer, intent(in ) :: ichoice_in,ichoicem_in,ichoice_s_in - logical, intent(in ) :: flag_init, flag_restart, do_mynnedmf + logical, intent(in ) :: flag_init, flag_restart, do_mynnedmf, gf_coldstart logical, intent(in ) :: flag_for_scnv_generic_tend,flag_for_dcnv_generic_tend real (kind=kind_phys), intent(in) :: g,cp,xlv,r_v logical, intent(in ) :: ldiag3d @@ -427,7 +427,7 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,& ccn_m(i) = 0. ! set aod and ccn - if (flag_init .and. .not.flag_restart) then + if ((flag_init .and. .not.flag_restart) .or. gf_coldstart) then aod_gf(i)=aodc0 else if((cactiv(i).eq.0) .and. (cactiv_m(i).eq.0))then diff --git a/physics/CONV/Grell_Freitas/cu_gf_driver.meta b/physics/CONV/Grell_Freitas/cu_gf_driver.meta index d0b661fd8..57b4e900b 100644 --- a/physics/CONV/Grell_Freitas/cu_gf_driver.meta +++ b/physics/CONV/Grell_Freitas/cu_gf_driver.meta @@ -577,6 +577,13 @@ dimensions = () type = logical intent = in +[gf_coldstart] + standard_name = flag_for_cold_start_gf + long_name = flag to cold start G-F + units = flag + dimensions = () + type = logical + intent = in [ichoice_in] standard_name = identifier_for_c3_or_gf_deep_convection_closure long_name = flag for C3 or GF deep convection closure diff --git a/physics/CONV/Grell_Freitas/cu_gf_driver_pre.F90 b/physics/CONV/Grell_Freitas/cu_gf_driver_pre.F90 index 7ff66be21..f48daa112 100644 --- a/physics/CONV/Grell_Freitas/cu_gf_driver_pre.F90 +++ b/physics/CONV/Grell_Freitas/cu_gf_driver_pre.F90 @@ -15,7 +15,7 @@ module cu_gf_driver_pre !> \section arg_table_cu_gf_driver_pre_run Argument Table !! \htmlinclude cu_gf_driver_pre_run.html !! - subroutine cu_gf_driver_pre_run (flag_init, flag_restart, kdt, fhour, dtp, t, q, prevst, prevsq, & + subroutine cu_gf_driver_pre_run (flag_init, flag_restart, gf_coldstart, kdt, fhour, dtp, t, q, prevst, prevsq, & forcet, forceq, cactiv, cactiv_m, conv_act, conv_act_m, & rrfs_sd, ntsmoke, ntdust, ntcoarsepm, chem3d, gq0, & errmsg, errflg) @@ -26,6 +26,7 @@ subroutine cu_gf_driver_pre_run (flag_init, flag_restart, kdt, fhour, dtp, t, q, logical, intent(in) :: flag_init logical, intent(in) :: flag_restart + logical, intent(in) :: gf_coldstart logical, intent(in) :: rrfs_sd integer, intent(in) :: kdt real(kind_phys), intent(in) :: fhour @@ -58,7 +59,7 @@ subroutine cu_gf_driver_pre_run (flag_init, flag_restart, kdt, fhour, dtp, t, q, ! For restart runs, can assume that prevst and prevsq ! are read from the restart files beforehand, same ! for conv_act. - if(flag_init .and. .not.flag_restart) then + if((flag_init .and. .not.flag_restart) .or. gf_coldstart) then !$acc kernels forcet(:,:)=0.0 forceq(:,:)=0.0 diff --git a/physics/CONV/Grell_Freitas/cu_gf_driver_pre.meta b/physics/CONV/Grell_Freitas/cu_gf_driver_pre.meta index ff22f1583..b86a378b6 100644 --- a/physics/CONV/Grell_Freitas/cu_gf_driver_pre.meta +++ b/physics/CONV/Grell_Freitas/cu_gf_driver_pre.meta @@ -122,6 +122,13 @@ type = real kind = kind_phys intent = in +[gf_coldstart] + standard_name = flag_for_cold_start_gf + long_name = flag to cold start G-F + units = flag + dimensions = () + type = logical + intent = in [rrfs_sd] standard_name = do_smoke_coupling long_name = flag controlling rrfs_sd collection From 7223b8c7c083202cc6ff533408392a9140a5a579 Mon Sep 17 00:00:00 2001 From: Jili Dong Date: Fri, 13 Sep 2024 20:25:20 +0000 Subject: [PATCH 10/15] saSAS sigmab initialization modification --- physics/CONV/SAMF/samfdeepcnv.f | 9 +++++---- physics/CONV/SAMF/samfdeepcnv.meta | 7 +++++++ physics/CONV/progsigma_calc.f90 | 28 +++++++++++++++++----------- 3 files changed, 29 insertions(+), 15 deletions(-) diff --git a/physics/CONV/SAMF/samfdeepcnv.f b/physics/CONV/SAMF/samfdeepcnv.f index 5853254c0..00553b351 100644 --- a/physics/CONV/SAMF/samfdeepcnv.f +++ b/physics/CONV/SAMF/samfdeepcnv.f @@ -12,10 +12,11 @@ module samfdeepcnv contains subroutine samfdeepcnv_init(imfdeepcnv,imfdeepcnv_samf, & - & errmsg, errflg) + & sigmab_coldstart,errmsg, errflg) integer, intent(in) :: imfdeepcnv integer, intent(in) :: imfdeepcnv_samf + logical, intent(in) :: sigmab_coldstart character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -84,7 +85,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & & clam,c0s,c1,betal,betas,evef,pgcon,asolfac, & & do_ca, ca_closure, ca_entr, ca_trigger, nthresh,ca_deep, & & rainevap,sigmain,sigmaout,betadcu,betamcu,betascu, & - & maxMF, do_mynnedmf,errmsg,errflg) + & maxMF, do_mynnedmf,sigmab_coldstart,errmsg,errflg) ! use machine , only : kind_phys use funcphys , only : fpvs @@ -100,7 +101,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & & prslp(:,:), garea(:), hpbl(:), dot(:,:), phil(:,:) real(kind=kind_phys), dimension(:), intent(in) :: fscav logical, intent(in) :: first_time_step,restart,hwrf_samfdeep, & - & progsigma,do_mynnedmf + & progsigma,do_mynnedmf,sigmab_coldstart real(kind=kind_phys), intent(in) :: nthresh,betadcu,betamcu, & & betascu real(kind=kind_phys), intent(in) :: ca_deep(:) @@ -2915,7 +2916,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & if(progsigma)then !Initial computations, dynamic q-tendency - if(first_time_step .and. .not.restart)then + if(first_time_step .and. (.not.restart .or. sigmab_coldstart))then do k = 1,km do i = 1,im qadv(i,k)=0. diff --git a/physics/CONV/SAMF/samfdeepcnv.meta b/physics/CONV/SAMF/samfdeepcnv.meta index 2dbd4407c..f015de39e 100644 --- a/physics/CONV/SAMF/samfdeepcnv.meta +++ b/physics/CONV/SAMF/samfdeepcnv.meta @@ -488,6 +488,13 @@ dimensions = () type = logical intent = in +[sigmab_coldstart] + standard_name = flag_to_cold_start_for_sigmab_init + long_name = flag to cold start for sigmab initialization + units = flag + dimensions = () + type = logical + intent = in [qlcn] standard_name = mass_fraction_of_convective_cloud_liquid_water long_name = mass fraction of convective cloud liquid water diff --git a/physics/CONV/progsigma_calc.f90 b/physics/CONV/progsigma_calc.f90 index 469df49f6..5ab471002 100644 --- a/physics/CONV/progsigma_calc.f90 +++ b/physics/CONV/progsigma_calc.f90 @@ -54,7 +54,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_shallow,& real(kind=kind_phys) :: gcvalmx,epsilon,ZZ,cvg,mcon,buy2, & fdqb,dtdyn,dxlim,rmulacvg,tem, & - DEN,dp1,invdelt + DEN,dp1,invdelt,sigmind_new !Parameters gcvalmx = 0.1 @@ -63,6 +63,12 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_shallow,& km1=km-1 invdelt = 1./delt + if (flag_init) then + sigmind_new=0.0 + else + sigmind_new=sigmind + end if + !Initialization 2D do k = 1,km do i = 1,im @@ -168,13 +174,13 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_shallow,& enddo !sigmab - if(flag_init .and. .not. flag_restart)then - do i = 1,im - if(cnvflg(i))then - sigmab(i)=0.03 - endif - enddo - else +! if(flag_init .and. .not. flag_restart)then +! do i = 1,im +! if(cnvflg(i))then +! sigmab(i)=0.03 +! endif +! enddo +! else do i = 1,im if(cnvflg(i))then DEN=MIN(termC(i)+termB(i),1.E8) @@ -186,11 +192,11 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_shallow,& sigmab(i)=(ZZ*(termA(i)+cvg))/(DEN+(1.0-ZZ)) if(sigmab(i)>0.)then sigmab(i)=MIN(sigmab(i),0.95) - sigmab(i)=MAX(sigmab(i),0.01) + sigmab(i)=MAX(sigmab(i),sigmind_new) endif endif!cnvflg enddo - endif +! endif do k=1,km do i=1,im @@ -219,7 +225,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_shallow,& do i= 1, im if(cnvflg(i)) then sigmab(i)=sigmab(i)/betadcu - sigmab(i)=MAX(sigmind,sigmab(i)) + sigmab(i)=MAX(sigmind_new,sigmab(i)) endif enddo endif From 8302de4afc9b693974d26d98f6d15cbcda65e93a Mon Sep 17 00:00:00 2001 From: Jili Dong Date: Mon, 16 Sep 2024 13:40:09 +0000 Subject: [PATCH 11/15] format fix --- physics/CONV/SAMF/samfdeepcnv.f | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/physics/CONV/SAMF/samfdeepcnv.f b/physics/CONV/SAMF/samfdeepcnv.f index 00553b351..beeafcd14 100644 --- a/physics/CONV/SAMF/samfdeepcnv.f +++ b/physics/CONV/SAMF/samfdeepcnv.f @@ -12,11 +12,10 @@ module samfdeepcnv contains subroutine samfdeepcnv_init(imfdeepcnv,imfdeepcnv_samf, & - & sigmab_coldstart,errmsg, errflg) + & errmsg, errflg) integer, intent(in) :: imfdeepcnv integer, intent(in) :: imfdeepcnv_samf - logical, intent(in) :: sigmab_coldstart character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -2916,7 +2915,8 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & if(progsigma)then !Initial computations, dynamic q-tendency - if(first_time_step .and. (.not.restart .or. sigmab_coldstart))then + if(first_time_step .and. (.not.restart + & .or. sigmab_coldstart))then do k = 1,km do i = 1,im qadv(i,k)=0. From a2f5f8337be64f941121424bd6ec14a90976f946 Mon Sep 17 00:00:00 2001 From: Jili Dong Date: Thu, 19 Sep 2024 18:32:35 +0000 Subject: [PATCH 12/15] code cleanup --- physics/CONV/progsigma_calc.f90 | 7 ------- 1 file changed, 7 deletions(-) diff --git a/physics/CONV/progsigma_calc.f90 b/physics/CONV/progsigma_calc.f90 index 5ab471002..2b38ffe51 100644 --- a/physics/CONV/progsigma_calc.f90 +++ b/physics/CONV/progsigma_calc.f90 @@ -174,13 +174,6 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_shallow,& enddo !sigmab -! if(flag_init .and. .not. flag_restart)then -! do i = 1,im -! if(cnvflg(i))then -! sigmab(i)=0.03 -! endif -! enddo -! else do i = 1,im if(cnvflg(i))then DEN=MIN(termC(i)+termB(i),1.E8) From b2d503256ababff56b0ba1b03bb7ee4e8c5ccf83 Mon Sep 17 00:00:00 2001 From: Jili Dong Date: Tue, 24 Sep 2024 17:40:24 +0000 Subject: [PATCH 13/15] code clean up #2 --- physics/CONV/progsigma_calc.f90 | 1 - 1 file changed, 1 deletion(-) diff --git a/physics/CONV/progsigma_calc.f90 b/physics/CONV/progsigma_calc.f90 index 2b38ffe51..f1415e89f 100644 --- a/physics/CONV/progsigma_calc.f90 +++ b/physics/CONV/progsigma_calc.f90 @@ -189,7 +189,6 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_shallow,& endif endif!cnvflg enddo -! endif do k=1,km do i=1,im From 1e5a6d3af61eae38a6e431d47548e7a33e54893f Mon Sep 17 00:00:00 2001 From: Grant Firl Date: Fri, 7 Feb 2025 10:06:44 -0500 Subject: [PATCH 14/15] remove changes not in ufs/dev PR#225 --- .../Interstitials/UFS_SCM_NEPTUNE/GFS_debug.F90 | 1 - .../Interstitials/UFS_SCM_NEPTUNE/dcyc2t3.meta | 16 ---------------- physics/Radiation/RRTMG/rrtmg_sw_post.F90 | 6 ++---- physics/Radiation/RRTMG/rrtmg_sw_post.meta | 8 -------- 4 files changed, 2 insertions(+), 29 deletions(-) diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_debug.F90 b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_debug.F90 index 37557f533..cdd3d8e2b 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_debug.F90 +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_debug.F90 @@ -694,7 +694,6 @@ subroutine GFS_diagtoscreen_run (Model, Statein, Stateout, Sfcprop, Coupling, call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Diag%dlwsfci ', Diag%dlwsfci) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Diag%ulwsfci ', Diag%ulwsfci) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Diag%dswsfci ', Diag%dswsfci) - call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Diag%dswsfcci ', Diag%dswsfcci) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Diag%nswsfci ', Diag%nswsfci) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Diag%uswsfci ', Diag%uswsfci) call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Diag%dusfci ', Diag%dusfci) diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/dcyc2t3.meta b/physics/Interstitials/UFS_SCM_NEPTUNE/dcyc2t3.meta index bf6fb1a47..b2187f0c5 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/dcyc2t3.meta +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/dcyc2t3.meta @@ -191,14 +191,6 @@ type = real kind = kind_phys intent = in -[sfcdswc] - standard_name = surface_downwelling_shortwave_flux_on_radiation_timestep_assuming_clear_sky - long_name = clear sky surface downwelling shortwave flux on radiation time step - units = W m-2 - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = in [sfcdlw] standard_name = surface_downwelling_longwave_flux_on_radiation_timestep long_name = total sky surface downwelling longwave flux on radiation time step @@ -523,14 +515,6 @@ type = real kind = kind_phys intent = out -[adjsfcdswc] - standard_name = surface_downwelling_shortwave_flux_assuming_clear_sky - long_name = surface downwelling shortwave flux at current time assuming clear sky - units = W m-2 - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = out [adjsfcnsw] standard_name = surface_net_downwelling_shortwave_flux long_name = surface net downwelling shortwave flux at current time diff --git a/physics/Radiation/RRTMG/rrtmg_sw_post.F90 b/physics/Radiation/RRTMG/rrtmg_sw_post.F90 index 8fa9494b8..f39cba71c 100644 --- a/physics/Radiation/RRTMG/rrtmg_sw_post.F90 +++ b/physics/Radiation/RRTMG/rrtmg_sw_post.F90 @@ -14,7 +14,7 @@ subroutine rrtmg_sw_post_run (im, levr, levs, ltp, nday, lm, kd, lsswr, & swhtr, sfcalb1, sfcalb2, sfcalb3, sfcalb4, htswc, htsw0, & nirbmdi, nirdfdi, visbmdi, visdfdi, nirbmui, nirdfui, visbmui,& visdfui, sfcdsw, sfcnsw, htrsw, swhc, scmpsw, sfcfsw, topfsw, & - sfcdswc, errmsg, errflg) + errmsg, errflg) use machine, only: kind_phys use module_radsw_parameters, only: topfsw_type, sfcfsw_type, & @@ -33,8 +33,7 @@ subroutine rrtmg_sw_post_run (im, levr, levs, ltp, nday, lm, kd, lsswr, & visbmdi, visdfdi, & nirbmui, nirdfui, & visbmui, visdfui, & - sfcdsw, sfcnsw, & - sfcdswc + sfcdsw, sfcnsw real(kind=kind_phys), dimension(:,:), intent(inout) :: htrsw, swhc type(cmpfsw_type), dimension(:), intent(inout) :: scmpsw @@ -123,7 +122,6 @@ subroutine rrtmg_sw_post_run (im, levr, levs, ltp, nday, lm, kd, lsswr, & do i=1,im sfcnsw(i) = sfcfsw(i)%dnfxc - sfcfsw(i)%upfxc sfcdsw(i) = sfcfsw(i)%dnfxc - sfcdswc(i)= sfcfsw(i)%dnfx0 enddo endif ! end_if_lsswr diff --git a/physics/Radiation/RRTMG/rrtmg_sw_post.meta b/physics/Radiation/RRTMG/rrtmg_sw_post.meta index 0a48b04d7..9914051ce 100644 --- a/physics/Radiation/RRTMG/rrtmg_sw_post.meta +++ b/physics/Radiation/RRTMG/rrtmg_sw_post.meta @@ -198,14 +198,6 @@ type = real kind = kind_phys intent = inout -[sfcdswc] - standard_name = surface_downwelling_shortwave_flux_on_radiation_timestep_assuming_clear_sky - long_name = clear sky sfc downward sw flux - units = W m-2 - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = inout [htrsw] standard_name = tendency_of_air_temperature_due_to_shortwave_heating_on_radiation_timestep long_name = total sky sw heating rate From 26783c859477fd13519c814a28b40b0824b644fe Mon Sep 17 00:00:00 2001 From: Grant Firl Date: Fri, 14 Feb 2025 10:43:05 -0500 Subject: [PATCH 15/15] fix line leading to reproducibility problems in progsigma_calc.f90 --- physics/CONV/progsigma_calc.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/physics/CONV/progsigma_calc.f90 b/physics/CONV/progsigma_calc.f90 index f1415e89f..a6f216a1c 100644 --- a/physics/CONV/progsigma_calc.f90 +++ b/physics/CONV/progsigma_calc.f90 @@ -63,7 +63,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart,flag_shallow,& km1=km-1 invdelt = 1./delt - if (flag_init) then + if(flag_init .and. .not. flag_restart) then sigmind_new=0.0 else sigmind_new=sigmind