diff --git a/physics/CONV/Grell_Freitas/cu_gf_deep.F90 b/physics/CONV/Grell_Freitas/cu_gf_deep.F90 index a1bca36c9..8a2c73600 100644 --- a/physics/CONV/Grell_Freitas/cu_gf_deep.F90 +++ b/physics/CONV/Grell_Freitas/cu_gf_deep.F90 @@ -142,13 +142,13 @@ subroutine cu_gf_deep_run( & !! betwee -1 and +1 ,do_capsuppress,cap_suppress_j & ! ,k22 & ! - ,jmin,tropics) ! + ,jmin,kdt,tropics) ! implicit none integer & ,intent (in ) :: & - nranflag,itf,ktf,its,ite, kts,kte,ipr,imid + nranflag,itf,ktf,its,ite, kts,kte,ipr,imid,kdt integer, intent (in ) :: & ichoice,nchem real(kind=kind_phys), dimension (its:ite,4) & @@ -591,6 +591,7 @@ subroutine cu_gf_deep_run( & sig(i)=(1.-frh)**2 !frh_out(i) = frh if(forcing(i,7).eq.0.)sig(i)=1. + if(kdt.le.(3600./dtime))sig(i)=1. frh_out(i) = frh*sig(i) enddo !$acc end kernels diff --git a/physics/CONV/Grell_Freitas/cu_gf_driver.F90 b/physics/CONV/Grell_Freitas/cu_gf_driver.F90 index 92f8760b0..54a23ca74 100644 --- a/physics/CONV/Grell_Freitas/cu_gf_driver.F90 +++ b/physics/CONV/Grell_Freitas/cu_gf_driver.F90 @@ -68,7 +68,7 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,& dfi_radar_max_intervals,ldiag3d,qci_conv,do_cap_suppress, & maxupmf,maxMF,do_mynnedmf,ichoice_in,ichoicem_in,ichoice_s_in, & spp_cu_deep,spp_wts_cu_deep,nchem,chem3d,fscav,wetdpc_deep, & - do_smoke_transport,errmsg,errflg) + do_smoke_transport,kdt,errmsg,errflg) !------------------------------------------------------------- implicit none integer, parameter :: maxiens=1 @@ -95,7 +95,7 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,& integer :: ishallow_g3 ! depend on imfshalcnv !------------------------------------------------------------- integer :: its,ite, jts,jte, kts,kte - integer, intent(in ) :: im,km,ntracer, nchem + 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_for_scnv_generic_tend,flag_for_dcnv_generic_tend @@ -766,7 +766,7 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,& ! betwee -1 and +1 ,do_cap_suppress_here,cap_suppress_j & ,k22m & - ,jminm,tropics) + ,jminm,kdt,tropics) !$acc kernels do i=its,itf do k=kts,ktf @@ -853,7 +853,7 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,& ! betwee -1 and +1 ,do_cap_suppress_here,cap_suppress_j & ,k22 & - ,jmin,tropics) + ,jmin,kdt,tropics) jpr=0 ipr=0 !$acc kernels diff --git a/physics/CONV/Grell_Freitas/cu_gf_driver.meta b/physics/CONV/Grell_Freitas/cu_gf_driver.meta index fe9b4c375..d0b661fd8 100644 --- a/physics/CONV/Grell_Freitas/cu_gf_driver.meta +++ b/physics/CONV/Grell_Freitas/cu_gf_driver.meta @@ -651,6 +651,13 @@ type = real kind = kind_phys intent = inout +[kdt] + standard_name = index_of_timestep + long_name = current forecast iteration + units = index + dimensions = () + type = integer + intent = in [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP diff --git a/physics/PBL/MYNN_EDMF/module_bl_mynn.F90 b/physics/PBL/MYNN_EDMF/module_bl_mynn.F90 index 6840f80bf..cc7a47ce6 100644 --- a/physics/PBL/MYNN_EDMF/module_bl_mynn.F90 +++ b/physics/PBL/MYNN_EDMF/module_bl_mynn.F90 @@ -2001,9 +2001,9 @@ SUBROUTINE mym_length ( & uonset= 15. wt_u = (1.0 - min(max(ugrid - uonset, 0.0)/30.0, 0.5)) cns = 2.7 !was 3.5 - alp1 = 0.22 + alp1 = 0.23 alp2 = 0.3 - alp3 = 2.0 * wt_u !taper off bouyancy enhancement in shear-driven pbls + alp3 = 2.5 * wt_u !taper off bouyancy enhancement in shear-driven pbls alp4 = 5.0 alp5 = 0.3 alp6 = 50. @@ -2059,12 +2059,12 @@ SUBROUTINE mym_length ( & ! ** Length scale limited by the buoyancy effect ** IF ( dtv(k) .GT. 0.0 ) THEN - bv = max( sqrt( gtr*dtv(k) ), 0.001) + bv = max( sqrt( gtr*dtv(k) ), 0.0001) elb = MAX(alp2*qkw(k), & & alp6*edmf_a1(k-1)*edmf_w1(k-1)) / bv & & *( 1.0 + alp3*SQRT( vsc/(bv*elt) ) ) elb = MIN(elb, zwk) - elf = 0.80 * qkw(k)/bv + elf = 1.0 * qkw(k)/bv elBLavg(k) = MAX(elBLavg(k), alp6*edmf_a1(k-1)*edmf_w1(k-1)/bv) ELSE elb = 1.0e10 @@ -2084,8 +2084,10 @@ SUBROUTINE mym_length ( & !add blending to use BouLac mixing length in free atmos; !defined relative to the PBLH (zi) + transition layer (h1) !el(k) = MIN(elb/( elb/elt+elb/els+1.0 ),elf) - !try squared-blending - el(k) = SQRT( els**2/(1. + (els**2/elt**2) +(els**2/elb**2))) + !try squared-blending - but take out elb (makes it underdiffusive) + !el(k) = SQRT( els**2/(1. + (els**2/elt**2) +(els**2/elb**2))) + el(k) = sqrt( els**2/(1. + (els**2/elt**2))) + el(k) = min(el(k), elb) el(k) = MIN (el(k), elf) el(k) = el(k)*(1.-wt) + alp5*elBLavg(k)*wt @@ -3633,13 +3635,13 @@ SUBROUTINE mym_condensation (kts,kte, & real(kind_phys):: qsl,esat,qsat,dqsl,cld0,q1k,qlk,eq1,qll, & &q2p,pt,rac,qt,t,xl,rsl,cpm,Fng,qww,alpha,beta,bb, & - &ls,wt,qpct,cld_factor,fac_damp,liq_frac,ql_ice,ql_water, & + &ls,wt,wt2,qpct,cld_factor,fac_damp,liq_frac,ql_ice,ql_water, & &qmq,qsat_tk,q1_rh,rh_hack,dzm1,zsl,maxqc real(kind_phys), parameter :: qpct_sfc=0.025 real(kind_phys), parameter :: qpct_pbl=0.030 real(kind_phys), parameter :: qpct_trp=0.040 real(kind_phys), parameter :: rhcrit =0.83 !for cloudpdf = 2 - real(kind_phys), parameter :: rhmax =1.01 !for cloudpdf = 2 + real(kind_phys), parameter :: rhmax =1.02 !for cloudpdf = 2 integer :: i,j,k real(kind_phys):: erf @@ -3864,25 +3866,18 @@ SUBROUTINE mym_condensation (kts,kte, & !Add condition for falling/settling into low-RH layers, so at least !some cloud fraction is applied for all qc, qs, and qi. rh_hack= rh(k) + wt2 = min(max( zagl - pblh2, 0.0 )/300., 1.0) !ensure adequate RH & q1 when qi is at least 1e-9 (above the PBLH) - if (qi(k)>1.e-9 .and. zagl .gt. pblh2) then - rh_hack =min(rhmax, rhcrit + 0.07*(9.0 + log10(qi(k)))) + if ((qi(k)+qs(k))>1.e-9 .and. (zagl .gt. pblh2)) then + rh_hack =min(rhmax, rhcrit + wt2*0.045*(9.0 + log10(qi(k)+qs(k)))) rh(k) =max(rh(k), rh_hack) !add rh-based q1 q1_rh =-3. + 3.*(rh(k)-rhcrit)/(1.-rhcrit) q1(k) =max(q1_rh, q1(k) ) endif - !ensure adequate RH & q1 when qc is at least 1e-6 - if (qc(k)>1.e-6) then - rh_hack =min(rhmax, rhcrit + 0.09*(6.0 + log10(qc(k)))) - rh(k) =max(rh(k), rh_hack) - !add rh-based q1 - q1_rh =-3. + 3.*(rh(k)-rhcrit)/(1.-rhcrit) - q1(k) =max(q1_rh, q1(k) ) - endif - !ensure adequate RH & q1 when qs is at least 1e-8 (above the PBLH) - if (qs(k)>1.e-8 .and. zagl .gt. pblh2) then - rh_hack =min(rhmax, rhcrit + 0.07*(8.0 + log10(qs(k)))) + !ensure adequate rh & q1 when qc is at least 1e-6 (above the PBLH) + if (qc(k)>1.e-6 .and. (zagl .gt. pblh2)) then + rh_hack =min(rhmax, rhcrit + wt2*0.08*(6.0 + log10(qc(k)))) rh(k) =max(rh(k), rh_hack) !add rh-based q1 q1_rh =-3. + 3.*(rh(k)-rhcrit)/(1.-rhcrit) @@ -3994,7 +3989,7 @@ SUBROUTINE mym_condensation (kts,kte, & fac_damp = min(zagl * 0.0025, 1.0) !cld_factor = 1.0 + fac_damp*MAX(0.0, ( RH(k) - 0.75 ) / 0.26 )**1.9 !HRRRv4 !cld_factor = 1.0 + fac_damp*min((max(0.0, ( RH(k) - 0.92 )) / 0.25 )**2, 0.3) - cld_factor = 1.0 + fac_damp*min((max(0.0, ( RH(k) - 0.92 )) / 0.145)**2, 0.35) + cld_factor = 1.0 + fac_damp*min((max(0.0, ( RH(k) - 0.92 )) / 0.145)**2, 0.37) cldfra_bl1D(K) = min( 1., cld_factor*cldfra_bl1D(K) ) enddo @@ -4181,38 +4176,33 @@ SUBROUTINE mynn_tendencies(kts,kte,i, & k=kts -!original approach (drag in b-vector): -! a(1)=0. -! b(1)=1. + dtz(k)*(dfm(k+1)+ust**2/wspd) - 0.5*dtz(k)*s_aw(k+1)*onoff -! c(1)=-dtz(k)*dfm(k+1) - 0.5*dtz(k)*s_aw(k+1)*onoff -! d(1)=u(k) + dtz(k)*uoce*ust**2/wspd - dtz(k)*s_awu(k+1)*onoff + & -! sub_u(k)*delt + det_u(k)*delt - !rho-weighted (drag in b-vector): a(k)= -dtz(k)*kmdz(k)*rhoinv(k) - b(k)=1.+dtz(k)*(kmdz(k+1)+rhosfc*ust**2/wspd)*rhoinv(k) & - & - 0.5*dtz(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff + b(k)=1.+dtz(k)*(kmdz(k+1)+rhosfc*ust**2/wspd)*rhoinv(k) & + & - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff & + & - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff c(k)= -dtz(k)*kmdz(k+1)*rhoinv(k) & - & - 0.5*dtz(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff - d(k)=u(k) + dtz(k)*uoce*ust**2/wspd - dtz(k)*s_awu(k+1)*onoff - & - & dtz(k)*rhoinv(k)*sd_awu(k+1)*onoff + sub_u(k)*delt + det_u(k)*delt - -!rho-weighted with drag term moved out of b-array -! a(k)= -dtz(k)*kmdz(k)*rhoinv(k) -! b(k)=1.+dtz(k)*(kmdz(k+1))*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff -! c(k)= -dtz(k)*kmdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff -! d(k)=u(k)*(1.-ust**2/wspd*dtz(k)*rhosfc/rho(k)) + dtz(k)*uoce*ust**2/wspd - & -! !!!d(k)=u(k)*(1.-ust**2/wspd*dtz(k)) + dtz(k)*uoce*ust**2/wspd - & -! & dtz(k)*rhoinv(k)*s_awu(k+1)*onoff - dtz(k)*rhoinv(k)*sd_awu(k+1)*onoff + sub_u(k)*delt + det_u(k)*delt - - DO k=kts+1,kte-1 - a(k)= -dtz(k)*kmdz(k)*rhoinv(k) + 0.5*dtz(k)*rhoinv(k)*s_aw(k)*onoff + 0.5*dtz(k)*rhoinv(k)*sd_aw(k)*onoff - b(k)=1.+dtz(k)*(kmdz(k)+kmdz(k+1))*rhoinv(k) + & - & 0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1))*onoff + 0.5*dtz(k)*rhoinv(k)*(sd_aw(k)-sd_aw(k+1))*onoff - c(k)= -dtz(k)*kmdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff - d(k)=u(k) + dtz(k)*rhoinv(k)*(s_awu(k)-s_awu(k+1))*onoff + dtz(k)*rhoinv(k)*(sd_awu(k)-sd_awu(k+1))*onoff + & - & sub_u(k)*delt + det_u(k)*delt - ENDDO + & - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff & + & - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff + d(k)=u(k) + dtz(k)*uoce*ust**2/wspd & + & - dtz(k)*rhoinv(k)*s_awu(k+1)*onoff & + & + dtz(k)*rhoinv(k)*sd_awu(k+1)*onoff & + & + sub_u(k)*delt + det_u(k)*delt + + do k=kts+1,kte-1 + a(k)= -dtz(k)*kmdz(k)*rhoinv(k) & + & + 0.5*dtz(k)*rhoinv(k)*s_aw(k)*onoff & + & + 0.5*dtz(k)*rhoinv(k)*sd_aw(k)*onoff + b(k)=1.+ dtz(k)*(kmdz(k)+kmdz(k+1))*rhoinv(k) & + & + 0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1))*onoff & + & + 0.5*dtz(k)*rhoinv(k)*(sd_aw(k)-sd_aw(k+1))*onoff + c(k)= - dtz(k)*kmdz(k+1)*rhoinv(k) & + & - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff & + & - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff + d(k)=u(k) + dtz(k)*rhoinv(k)*(s_awu(k)-s_awu(k+1))*onoff & + & - dtz(k)*rhoinv(k)*(sd_awu(k)-sd_awu(k+1))*onoff & + & + sub_u(k)*delt + det_u(k)*delt + enddo !! no flux at the top ! a(kte)=-1. @@ -4247,37 +4237,33 @@ SUBROUTINE mynn_tendencies(kts,kte,i, & k=kts -!original approach (drag in b-vector): -! a(1)=0. -! b(1)=1. + dtz(k)*(dfm(k+1)+ust**2/wspd) - 0.5*dtz(k)*s_aw(k+1)*onoff -! c(1)= - dtz(k)*dfm(k+1) - 0.5*dtz(k)*s_aw(k+1)*onoff -! d(1)=v(k) + dtz(k)*voce*ust**2/wspd - dtz(k)*s_awv(k+1)*onoff + & -! sub_v(k)*delt + det_v(k)*delt - !rho-weighted (drag in b-vector): a(k)= -dtz(k)*kmdz(k)*rhoinv(k) - b(k)=1.+dtz(k)*(kmdz(k+1) + rhosfc*ust**2/wspd)*rhoinv(k) & - & - 0.5*dtz(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff - c(k)= -dtz(k)*kmdz(k+1)*rhoinv(k) - 0.5*dtz(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff - d(k)=v(k) + dtz(k)*voce*ust**2/wspd - dtz(k)*s_awv(k+1)*onoff - dtz(k)*rhoinv(k)*sd_awv(k+1)*onoff + & - & sub_v(k)*delt + det_v(k)*delt - -!rho-weighted with drag term moved out of b-array -! a(k)= -dtz(k)*kmdz(k)*rhoinv(k) -! b(k)=1.+dtz(k)*(kmdz(k+1))*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff -! c(k)= -dtz(k)*kmdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff -! d(k)=v(k)*(1.-ust**2/wspd*dtz(k)*rhosfc/rho(k)) + dtz(k)*voce*ust**2/wspd - & -! !!!d(k)=v(k)*(1.-ust**2/wspd*dtz(k)) + dtz(k)*voce*ust**2/wspd - & -! & dtz(k)*rhoinv(k)*s_awv(k+1)*onoff - dtz(k)*rhoinv(k)*sd_awv(k+1)*onoff + sub_v(k)*delt + det_v(k)*delt - - DO k=kts+1,kte-1 - a(k)= -dtz(k)*kmdz(k)*rhoinv(k) + 0.5*dtz(k)*rhoinv(k)*s_aw(k)*onoff + 0.5*dtz(k)*rhoinv(k)*sd_aw(k)*onoff - b(k)=1.+dtz(k)*(kmdz(k)+kmdz(k+1))*rhoinv(k) + & - & 0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1))*onoff + 0.5*dtz(k)*rhoinv(k)*(sd_aw(k)-sd_aw(k+1))*onoff - c(k)= -dtz(k)*kmdz(k+1)*rhoinv(k) - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff - d(k)=v(k) + dtz(k)*rhoinv(k)*(s_awv(k)-s_awv(k+1))*onoff + dtz(k)*rhoinv(k)*(sd_awv(k)-sd_awv(k+1))*onoff + & - & sub_v(k)*delt + det_v(k)*delt - ENDDO + b(k)=1.+dtz(k)*(kmdz(k+1) + rhosfc*ust**2/wspd)*rhoinv(k) & + & - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff & + & - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff + c(k)= -dtz(k)*kmdz(k+1)*rhoinv(k) & + & - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff & + & - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff + d(k)=v(k) + dtz(k)*voce*ust**2/wspd & + & - dtz(k)*rhoinv(k)*s_awv(k+1)*onoff & + & + dtz(k)*rhoinv(k)*sd_awv(k+1)*onoff & + & + sub_v(k)*delt + det_v(k)*delt + + do k=kts+1,kte-1 + a(k)= -dtz(k)*kmdz(k)*rhoinv(k) & + & + 0.5*dtz(k)*rhoinv(k)*s_aw(k)*onoff & + & + 0.5*dtz(k)*rhoinv(k)*sd_aw(k)*onoff + b(k)=1.+dtz(k)*(kmdz(k)+kmdz(k+1))*rhoinv(k) & + & + 0.5*dtz(k)*rhoinv(k)*(s_aw(k)-s_aw(k+1))*onoff & + & + 0.5*dtz(k)*rhoinv(k)*(sd_aw(k)-sd_aw(k+1))*onoff + c(k)= -dtz(k)*kmdz(k+1)*rhoinv(k) & + & - 0.5*dtz(k)*rhoinv(k)*s_aw(k+1)*onoff & + & - 0.5*dtz(k)*rhoinv(k)*sd_aw(k+1)*onoff + d(k)=v(k) + dtz(k)*rhoinv(k)*(s_awv(k)-s_awv(k+1))*onoff & + & - dtz(k)*rhoinv(k)*(sd_awv(k)-sd_awv(k+1))*onoff & + & + sub_v(k)*delt + det_v(k)*delt + enddo !! no flux at the top ! a(kte)=-1. diff --git a/physics/SFC_Models/Land/RUC/module_sf_ruclsm.F90 b/physics/SFC_Models/Land/RUC/module_sf_ruclsm.F90 index b15592052..2d01f96c9 100644 --- a/physics/SFC_Models/Land/RUC/module_sf_ruclsm.F90 +++ b/physics/SFC_Models/Land/RUC/module_sf_ruclsm.F90 @@ -1687,7 +1687,7 @@ SUBROUTINE SFCTMP (debug_print, delt,ktau,conflx,i,j, & !--- input varia endif if(newsn > zero ) then - SNOWFRACnewsn=MIN(one,SNHEI/SNHEI_CRIT_newsn) + SNOWFRACnewsn=MIN(one,snowfallac*1.e-3_kind_phys/SNHEI_CRIT_newsn) endif !-- due to steep slopes and blown snow, limit snow fraction in the @@ -1700,7 +1700,7 @@ SUBROUTINE SFCTMP (debug_print, delt,ktau,conflx,i,j, & !--- input varia if(snowfrac < 0.75_kind_phys) snow_mosaic = one KEEP_SNOW_ALBEDO = zero - IF (NEWSN > zero .and. snowfracnewsn > 0.99_kind_phys .and. rhosnfall < 450._kind_phys) THEN + IF (snowfracnewsn > 0.99_kind_phys .and. rhosnfall < 450._kind_phys) THEN ! new snow KEEP_SNOW_ALBEDO = one ! turn off separate treatment of snow covered and snow-free portions of the grid cell @@ -1735,7 +1735,7 @@ SUBROUTINE SFCTMP (debug_print, delt,ktau,conflx,i,j, & !--- input varia ! hwlps with these biases.. if( snow_mosaic == one) then ALBsn=alb_snow - if(newsn > zero .and. KEEP_SNOW_ALBEDO > 0.9_kind_phys .and. albsn < 0.4_kind_phys) then + if(KEEP_SNOW_ALBEDO > 0.9_kind_phys .and. albsn < 0.4_kind_phys) then !-- Albedo correction with fresh snow and deep snow pack !-- will reduce warm bias in western Canada !-- and US West coast, where max snow albedo is low (0.3-0.5). diff --git a/physics/smoke_dust/dep_dry_mod_emerson.F90 b/physics/smoke_dust/dep_dry_mod_emerson.F90 index 76fdc2411..e69d6bc3f 100755 --- a/physics/smoke_dust/dep_dry_mod_emerson.F90 +++ b/physics/smoke_dust/dep_dry_mod_emerson.F90 @@ -9,7 +9,7 @@ module dep_dry_emerson_mod use machine , only : kind_phys use dep_data_mod ! JLS - use rrfs_smoke_config, only : num_chem, p_smoke, p_dust_1, p_coarse_pm + use rrfs_smoke_config, only : num_chem, p_smoke, p_dust_1, p_coarse_pm, n_dbg_lines implicit none @@ -23,7 +23,7 @@ subroutine dry_dep_driver_emerson(rmol,ustar,znt,ndvel,ddvel,vgrav, & settling_flag,drydep_flux,settling_flux,dbg_opt, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & - its,ite, jts,jte, kts,kte ) + its,ite, jts,jte, kts,kte, curr_secs, mpiid, xlat, xlong ) ! ! compute dry deposition velocity for aerosol particles ! Based on Emerson et al. (2020), PNAS, @@ -37,6 +37,9 @@ subroutine dry_dep_driver_emerson(rmol,ustar,znt,ndvel,ddvel,vgrav, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte + + REAL(kind_phys) :: curr_secs + REAL(kind_phys), DIMENSION( ims:ime , jms:jme ) , & INTENT(IN) :: ustar, rmol, znt, snowh REAL(kind_phys), DIMENSION( ims:ime , kms:kme , jms:jme ), & @@ -80,6 +83,9 @@ subroutine dry_dep_driver_emerson(rmol,ustar,znt,ndvel,ddvel,vgrav, & real(kind_phys), dimension( kts:kte, ndvel ) :: cblk_col, vg_col integer, dimension(ndvel) :: ndt_settl integer :: i, j, k, ntdt, nv + integer :: icall=0 + integer, INTENT(IN) :: mpiid + real(kind_phys), DIMENSION(ims:ime,jms:jme), INTENT(IN) :: xlat,xlong ! chem pointers (p_*) are not sequentially numbered, need to define so nv loops work integer, dimension(ndvel) :: chem_pointers !> -- Gas constant @@ -87,11 +93,15 @@ subroutine dry_dep_driver_emerson(rmol,ustar,znt,ndvel,ddvel,vgrav, & chem_pointers(1) = p_smoke chem_pointers(2) = p_dust_1 chem_pointers(3) = p_coarse_pm - + growth_fac = 1.0 conver=1.e-9 converi=1.e9 + if (mod(int(curr_secs),1800) .eq. 0) then + icall = 0 + endif + do j = jts, jte do i = its, ite aer_res(i,j) = 0.0 @@ -116,7 +126,7 @@ subroutine dry_dep_driver_emerson(rmol,ustar,znt,ndvel,ddvel,vgrav, & aer_res(i,j) = max(aer_res(i,j)/100._kind_phys,0._kind_phys) ! Air kinematic viscosity (cm^2/s) airkinvisc = ( 1.8325e-4 * ( 416.16 / ( t_phy(i,k,j) + 120.0 ) ) * & - ( ( t_phy(i,k,j) / 296.16 )**1.5 ) ) / ( rho_phy(i,k,j) / 28.966e3 ) ! Convert density to mol/cm^3 + ( ( t_phy(i,k,j) / 296.16 )**1.5 ) ) / ( rho_phy(i,k,j) / 1.e3 ) ! Convert density to mol/cm^3 ! Air molecular freepath (cm) ! Check against XLM from above freepath = 7.39758e-4 * airkinvisc / sqrt( t_phy(i,k,j) ) do nv = 1, ndvel @@ -141,11 +151,11 @@ subroutine dry_dep_driver_emerson(rmol,ustar,znt,ndvel,ddvel,vgrav, & amu_corrected = amu / Cc ! Gravitational Settling vg = aerodens * dp * dp * gravity * 100. * Cc / & ! Convert gravity to cm/s^2 - ( 18. * airkinvisc * ( rho_phy(i,k,j) / 28.966e3 ) ) ! Convert density to mol/cm^3 + ( 18. * airkinvisc * ( rho_phy(i,k,j) / 1.e3 ) ) ! Convert density to mol/cm^3 ! -- Rest of loop for the surface when deposition velocity needs to be cacluated if ( k == kts ) then ! Brownian Diffusion - DDp = ( boltzmann * t_phy(i,k,j) ) * Cc / (3. * pi * airkinvisc * ( rho_phy(i,k,j) / 28.966e3 ) * dp) ! Convert density to mol/cm^3 + DDp = ( boltzmann * t_phy(i,k,j) ) * Cc / (3. * pi * airkinvisc * ( rho_phy(i,k,j) / 1.e3 ) * dp) ! Convert density to mol/cm^3 ! Schmit number Sc = airkinvisc / DDp ! Brownian Diffusion @@ -169,23 +179,23 @@ subroutine dry_dep_driver_emerson(rmol,ustar,znt,ndvel,ddvel,vgrav, & eps0 = eps0_grs end if ! Set if snow greater than 1 cm -! if ( snowh(i,j) .gt. 0.01 ) then ! snow -! A = A_wat -! eps0 = eps0_wat -! endif ! Interception Ein = Cin * ( dp / A )**vv ! Surface resistance Rs = 1. / ( ( ustar(i,j) * 100.) * ( Eb + Eim + Ein) * eps0 ) ! Convert ustar to cm/s ! Compute final ddvel = aer_res + RS, set max at max_dep_vel in dep_data_mod.F[ m/s] ! The /100. term converts from cm/s to m/s, required for MYNN. - ddvel(i,j,nv) = min( (1. / (aer_res(i,j) + Rs ))/100., max_dep_vel) - if ( dbg_opt ) then - WRITE(6,*) 'dry_dep_mod_emerson: i,j,nv',i,j,nv - WRITE(6,*) 'dry_dep_mod_emerson: deposition velocity (m/s) ',ddvel(i,j,nv) + if ( settling_flag == 1 ) then + ddvel(i,j,nv) = max(min( ( vg + 1./(aer_res(i,j)+Rs) )/100., max_dep_vel),0._kind_phys) + else + ddvel(i,j,nv) = max(min( ( 1./(aer_res(i,j)+Rs) )/100., max_dep_vel),0._kind_phys) endif - drydep_flux(i,j,nv) = chem(i,kts,j,chem_pointers(nv))*p_phy(i,kts,j) / & - (RSI*t_phy(i,kts,j))*ddvel(i,j,nv)*dt*1.E-6 + if ( dbg_opt .and. (icall .le. n_dbg_lines) ) then + WRITE(1000+mpiid,*) 'dry_dep_mod_emer:xlat,xlong,curr_secs,nv',xlat(i,j),xlong(i,j),int(curr_secs),nv + WRITE(1000+mpiid,*) 'dry_dep_mod_emer:xlat,xlong,curr_secs,deposition velocity (m/s)',xlat(i,j),xlong(i,j),int(curr_secs),ddvel(i,j,nv) + icall = icall + 1 + endif + drydep_flux(i,j,nv) = chem(i,kts,j,chem_pointers(nv))*rho_phy(i,k,j)*ddvel(i,j,nv)/100.0*dt endif ! k == kts vgrav(i,k,j,nv) = vg ! Fill column variables @@ -220,25 +230,15 @@ subroutine dry_dep_driver_emerson(rmol,ustar,znt,ndvel,ddvel,vgrav, & IF (ndt_settl(nv) > 12) ndt_settl(nv) = 12 dt_settl(nv) = REAL(ntdt,kind=kind_phys) /REAL(ndt_settl(nv),kind=kind_phys) enddo - do nv = 1, ndvel - chem_before(nv) = 0._kind_phys - do k = kts, kte - chem_before(nv) = chem_before(nv) + (cblk_col(k,nv) * rho_phy(i,k,j) * delz(i,k,j) ) ! ug/m2 - enddo - enddo ! Perform gravitational settling if desired if ( settling_flag == 1 ) then call particle_settling(cblk_col,rho_col,delz_col,vg_col,dt_settl,ndt_settl,ndvel,kts,kte) endif ! Put cblk back into chem array do nv= 1, ndvel - chem_after(nv) = 0._kind_phys - settling_flux(i,j,nv) = 0._kind_phys do k = kts, kte chem(i,k,j,chem_pointers(nv)) = cblk_col(k,nv) - chem_after(nv) = chem_after(nv) + (cblk_col(k,nv) * rho_phy(i,k,j) * delz(i,k,j) ) ! ug/m2 enddo ! k - settling_flux(i,j,nv) = settling_flux(i,j,nv) + (chem_before(nv) - chem_after(nv)) ! ug/m2 enddo ! nv end do ! j end do ! i diff --git a/physics/smoke_dust/dust_fengsha_mod.F90 b/physics/smoke_dust/dust_fengsha_mod.F90 index 54e66712d..6ec8f8d4a 100755 --- a/physics/smoke_dust/dust_fengsha_mod.F90 +++ b/physics/smoke_dust/dust_fengsha_mod.F90 @@ -21,8 +21,8 @@ module dust_fengsha_mod contains subroutine gocart_dust_fengsha_driver(dt, & - chem,rho_phy,smois,p8w,ssm, & - isltyp,snowh,xland,area,g,emis_dust, & + chem,rho_phy,smois,stemp,p8w,ssm, & + isltyp,snowh,xland,area,g,emis_dust, & ust,znt,clay,sand,rdrag,uthr, & num_emis_dust,num_chem,num_soil_layers, & ids,ide, jds,jde, kds,kde, & @@ -54,7 +54,7 @@ subroutine gocart_dust_fengsha_driver(dt, & REAL(kind_phys), DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN) :: rho_phy REAL(kind_phys), DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), INTENT(INOUT) :: chem REAL(kind_phys), DIMENSION( ims:ime, 1, jms:jme,num_emis_dust),OPTIONAL, INTENT(INOUT) :: emis_dust - REAL(kind_phys), DIMENSION( ims:ime, num_soil_layers, jms:jme ), INTENT(IN) :: smois + REAL(kind_phys), DIMENSION( ims:ime, num_soil_layers, jms:jme ), INTENT(IN) :: smois, stemp !0d input variables REAL(kind_phys), INTENT(IN) :: dt ! time step @@ -146,6 +146,11 @@ subroutine gocart_dust_fengsha_driver(dt, & ilwi = 0 endif + ! Don't emit over frozen soil + if (stemp(i,1,j) < 268.0) then ! -5C + ilwi = 0 + endif + ! Do not allow areas with bedrock, lava, or land-ice to loft IF (isltyp(i,j) .eq. 15 .or. isltyp(i,j) .eq. 16. .or. & diff --git a/physics/smoke_dust/module_add_emiss_burn.F90 b/physics/smoke_dust/module_add_emiss_burn.F90 index 95005b973..80d91bb0e 100755 --- a/physics/smoke_dust/module_add_emiss_burn.F90 +++ b/physics/smoke_dust/module_add_emiss_burn.F90 @@ -9,11 +9,11 @@ module module_add_emiss_burn subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi, & chem,julday,gmt,xlat,xlong, & fire_end_hr, peak_hr,time_int, & - coef_bb_dc, fhist, hwp, hwp_prevd, & + coef_bb_dc, fire_hist, hwp, hwp_prevd, & swdown,ebb_dcycle, ebu_in, ebu,fire_type,& ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & - its,ite, jts,jte, kts,kte ) + its,ite, jts,jte, kts,kte,mpiid ) IMPLICIT NONE @@ -22,6 +22,7 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte + INTEGER, INTENT(IN) :: mpiid real(kind_phys), DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), & INTENT(INOUT ) :: chem ! shall we set num_chem=1 here? @@ -29,7 +30,7 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi, & INTENT(INOUT ) :: ebu 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 !, vfrac + 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 @@ -38,17 +39,17 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi, & real(kind_phys), INTENT(IN) :: time_int,pi ! 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) :: fhist + real(kind_phys), DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: fire_hist !>--local integer :: i,j,k,n,m + integer :: icall=0 real(kind_phys) :: conv_rho, conv, dm_smoke, dc_hwp, dc_gp, dc_fn !daero_num_wfa, daero_num_ifa !, lu_sum1_5, lu_sum12_14 INTEGER, PARAMETER :: kfire_max=51 ! max vertical level for BB plume rise - ! real, parameter :: cx = 2.184936 * 3600., rinti = 2.1813936e-8 , ax = 2000.6038 - ! bx_bburn = 20.041288 * 3600., RAR: this depends on the vegetation class, location (local time) etc. - real(kind_phys) :: timeq, fire_age, age_hr, dt1,dt2,dtm ! For BB emis. diurnal cycle calculation - ! For Gaussian diurnal cycle + real(kind_phys) :: timeq, fire_age, age_hr, dt1,dt2,dtm, coef_con ! 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. @@ -88,35 +89,40 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi, & 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))*3600. + 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) SELECT CASE ( fire_type(i,j) ) !Ag, urban fires, bare land etc. CASE (1) ! these fires will have exponentially decreasing diurnal cycle, - ! We assume 1hr latency in ingesting the sat. data - coef_bb_dc(i,j) = 1._kind_phys/((2*pi)**0.5_kind_phys * sigma_fire_dur(1) *fire_age) * & - exp(- ( log(fire_age) - avg_fire_dur(1))**2 /(2*sigma_fire_dur(1)**2 )) + 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 )) + + 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 - IF (swdown(i,j)<.1 .AND. age_hr> 12. .AND. fhist(i,j)>0.75) THEN - fhist(i,j)= 0.75_kind_phys + IF (swdown(i,j)<.1 .AND. age_hr> 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. fhist(i,j)>0.5) THEN - fhist(i,j)= 0.5_kind_phys + IF (swdown(i,j)<.1 .AND. age_hr> 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. fhist(i,j)>0.25) THEN - fhist(i,j)= 0.25_kind_phys + IF (swdown(i,j)<.1 .AND. age_hr> 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= ebu_in(i,j)* hwp(i,j)/ MAX(1._kind_phys,hwp_prevd(i,j)) + dc_hwp= hwp(i,j)/ MAX(5._kind_phys,hwp_prevd(i,j)) dc_hwp= MAX(0._kind_phys,dc_hwp) - - !coef_bb_dc(i,j)= sc_factor* fhist(i,j)* rate_ebb2(i,j)* (1. + log( - !hwp_(i,j)/ hwp_day_avg(i,j))) + dc_hwp= MIN(25._kind_phys,dc_hwp) ! RAR: Gaussian profile for wildfires dt1= abs(timeq - peak_hr(i,j)) @@ -125,17 +131,29 @@ 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 = MAX(dc_hwp/dc_gp,3._kind_phys) - coef_bb_dc(i,j) = fhist(i,j)* dc_fn + 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) + WRITE(6,*) 'dc_gp,dc_hwp,dc_fn,coef_bb_dc(i,j) ',dc_gp,dc_hwp,dc_fn,coef_bb_dc(i,j) + END IF + CASE DEFAULT END SELECT enddo enddo endif + if (mod(int(time_int),1800) .eq. 0) then + icall = 0 + endif + do j=jts,jte do i=its,ite do k=kts,kfire_max + if (ebu(i,k,j)<0.001_kind_phys) cycle + if (ebb_dcycle==1) then conv= dtstep/(rho_phy(i,k,j)* dz8w(i,k,j)) elseif (ebb_dcycle==2) then @@ -143,14 +161,14 @@ subroutine add_emis_burn(dtstep,dz8w,rho_phy,pi, & endif dm_smoke= conv*ebu(i,k,j) chem(i,k,j,p_smoke) = chem(i,k,j,p_smoke) + dm_smoke - chem(i,k,j,p_smoke) = MIN(chem(i,k,j,p_smoke),5.e+3) + chem(i,k,j,p_smoke) = MIN(MAX(chem(i,k,j,p_smoke),0._kind_phys),5.e+3_kind_phys) - if ( dbg_opt .and. (k==kts .OR. k==kfire_max) ) then - WRITE(6,*) 'add_emiss_burn: i,j,k ',i,j,k - WRITE(6,*) 'add_emiss_burn: rho_phy(i,k,j),dz8w(i,k,j),conv',rho_phy(i,k,j),dz8w(i,k,j),conv - WRITE(6,*) 'add_emiss_burn: ebu(i,k,j),dm_smoke ', ebu(i,k,j),dm_smoke - endif + if ( dbg_opt .and. (k==kts .OR. k==kfire_max) .and. (icall .le. n_dbg_lines) ) then + WRITE(1000+mpiid,*) 'add_emiss_burn:xlat,xlong,curr_secs,fire_type,fire_hist,peak_hr', xlat(i,j),xlong(i,j),int(time_int),fire_type(i,j),fire_hist(i,j),peak_hr(i,j) + WRITE(1000+mpiid,*) 'add_emiss_burn:xlat,xlong,curr_secs,coef_bb_dc,ebu',xlat(i,j),xlong(i,j),int(time_int),coef_bb_dc(i,j),ebu(i,k,j) + endif enddo + icall = icall + 1 enddo enddo diff --git a/physics/smoke_dust/module_plumerise.F90 b/physics/smoke_dust/module_plumerise.F90 index 8a1d6ab25..5f7ef2a0e 100755 --- a/physics/smoke_dust/module_plumerise.F90 +++ b/physics/smoke_dust/module_plumerise.F90 @@ -24,10 +24,11 @@ subroutine ebu_driver ( flam_frac,ebu_in,ebu, & kpbl_thetav, & ! SRB: added kpbl_thetav ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & - its,ite, jts,jte, kts,kte, errmsg, errflg) + its,ite, jts,jte, kts,kte, errmsg, errflg,curr_secs, & + xlat, xlong , uspdavg2, hpbl_thetav2, mpiid) use rrfs_smoke_config - use plume_data_mod + !use plume_data_mod USE module_zero_plumegen_coms USE module_smoke_plumerise IMPLICIT NONE @@ -40,6 +41,8 @@ subroutine ebu_driver ( flam_frac,ebu_in,ebu, & real(kind=kind_phys), DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: frp_inst ! RAR: FRP array + real(kind_phys), DIMENSION(ims:ime,jms:jme), INTENT(IN) :: xlat,xlong ! SRB + real(kind_phys), DIMENSION(ims:ime, jms:jme), INTENT(IN) :: kpbl_thetav ! SRB character(*), intent(inout) :: errmsg @@ -47,6 +50,7 @@ subroutine ebu_driver ( flam_frac,ebu_in,ebu, & INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte + real(kind_phys) :: curr_secs INTEGER, INTENT(IN ) :: wind_eff_opt real(kind=kind_phys), DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT ) :: ebu real(kind=kind_phys), INTENT(IN ) :: g, con_cp, con_rd @@ -57,23 +61,32 @@ subroutine ebu_driver ( flam_frac,ebu_in,ebu, & ! Local variables... INTEGER :: nv, i, j, k, kp1, kp2 + INTEGER :: icall=0 INTEGER, DIMENSION(ims:ime, jms:jme), INTENT (OUT) :: k_min, k_max ! Min and max ver. levels for BB injection spread + REAL, DIMENSION(ims:ime, jms:jme), INTENT (OUT) :: uspdavg2, hpbl_thetav2 ! SRB real(kind_phys), dimension (kte) :: u_in ,v_in ,w_in ,theta_in ,pi_in, rho_phyin ,qv_in ,zmid, z_lev, uspd ! SRB real(kind=kind_phys) :: dz_plume, cpor, con_rocp, uspdavg ! SRB +! MPI variables + INTEGER, INTENT(IN) :: mpiid + cpor =con_cp/con_rd con_rocp=con_rd/con_cp - IF ( dbg_opt ) then - WRITE(*,*) 'module_plumerise: its,ite,jts,jte ', its,ite,jts,jte - WRITE(*,*) 'module_plumerise: ims,ime,jms,jme ', ims,ime,jms,jme - WRITE(*,*) 'module_plumerise: maxval(ebu(:,kts,:)) ', maxval(ebu(:,kts,:)) + if (mod(int(curr_secs),1800) .eq. 0) then + icall = 0 + endif + + IF ( dbg_opt .and. icall .le. n_dbg_lines) then + WRITE(1000+mpiid,*) 'module_plumerise: its,ite,jts,jte ', its,ite,jts,jte + WRITE(1000+mpiid,*) 'module_plumerise: ims,ime,jms,jme ', ims,ime,jms,jme + WRITE(1000+mpiid,*) 'module_plumerise: maxval(ebu(:,kts,:)) ', maxval(ebu(:,kts,:)) END IF ! RAR: setting to zero the ebu emissions at the levels k>1, this is necessary when the plumerise is called, so the emissions at k>1 are updated !do nv=1,num_ebu do j=jts,jte - do k=kts+1,kte + do k=kts,kte do i=its,ite ebu(i,k,j)=0. enddo @@ -112,12 +125,10 @@ subroutine ebu_driver ( flam_frac,ebu_in,ebu, & uspd(k)= wind_phy(i,k,j) ! SRB enddo - IF (dbg_opt) then - WRITE(*,*) 'module_plumerise: i,j ',i,j - WRITE(*,*) 'module_plumerise: frp_inst(i,j) ',frp_inst(i,j) - WRITE(*,*) 'module_plumerise: ebu(i,kts,j) ',ebu(i,kts,j) - WRITE(*,*) 'module_plumerise: u_in(10),v_in(10),w_in(kte),qv_in(10),pi_in(10) ',u_in(10),v_in(10),w_in(kte),qv_in(10),pi_in(10) - WRITE(*,*) 'module_plumerise: zmid(kte),z_lev(kte),rho_phyin(kte),theta_in(kte) ',zmid(kte),z_lev(kte),rho_phyin(kte),theta_in(kte) + + IF (dbg_opt .and. (icall .le. n_dbg_lines) .and. (frp_inst(i,j) .ge. frp_threshold) ) 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 ! RAR: the plume rise calculation step: @@ -127,7 +138,8 @@ subroutine ebu_driver ( flam_frac,ebu_in,ebu, & wind_eff_opt, & frp_inst(i,j), k_min(i,j), & k_max(i,j), dbg_opt, g, con_cp, & - con_rd, cpor, errmsg, errflg ) + con_rd, cpor, errmsg, errflg, & + icall, mpiid, xlat(i,j), xlong(i,j), curr_secs ) if(errflg/=0) return kp1= k_min(i,j) @@ -136,9 +148,13 @@ subroutine ebu_driver ( flam_frac,ebu_in,ebu, & ! 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 + +! 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_at_w(i,kpbl_thetav(i,j),j) .gt. zpbl_threshold) .AND. (wind_eff_opt .eq. 1)) THEN + (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 @@ -157,11 +173,18 @@ subroutine ebu_driver ( flam_frac,ebu_in,ebu, & END IF ! SRB: End modification - IF ( dbg_opt ) then - WRITE(*,*) 'module_plumerise: i,j ',i,j - WRITE(*,*) 'module_plumerise: k_min(i,j), k_max(i,j) ',kp1, kp2 ! SRB: replaced k_min, k_max with kp1, kp2 + IF ( dbg_opt .and. (icall .le. n_dbg_lines) .and. (frp_inst(i,j) .ge. frp_threshold) ) 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: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 END IF -! endif check_frp +! endif check_frp +! icall = icall + 1 enddo enddo diff --git a/physics/smoke_dust/module_smoke_plumerise.F90 b/physics/smoke_dust/module_smoke_plumerise.F90 index 0fca91de4..aa45890f4 100755 --- a/physics/smoke_dust/module_smoke_plumerise.F90 +++ b/physics/smoke_dust/module_smoke_plumerise.F90 @@ -14,10 +14,11 @@ module module_smoke_plumerise use machine , only : kind_phys - use plume_data_mod, only : num_frp_plume, p_frp_hr, p_frp_std + !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 @@ -28,12 +29,13 @@ 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 ) + cpor, errmsg, errflg, icall, mpiid, lat, long, curr_secs ) implicit none LOGICAL, INTENT (IN) :: dbg_opt - INTEGER, INTENT (IN) :: wind_eff_opt + INTEGER, INTENT (IN) :: wind_eff_opt, mpiid + real(kind_phys), INTENT(IN) :: lat,long, curr_secs ! SRB ! INTEGER, PARAMETER :: ihr_frp=1, istd_frp=2!, imean_fsize=3, istd_fsize=4 ! RAR: @@ -70,7 +72,7 @@ subroutine plumerise(m1,m2,m3,ia,iz,ja,jz, & ! real(kind=kind_phys), dimension(nveg_agreg) :: firesize,mean_fct INTEGER :: wind_eff - + INTEGER, INTENT(IN) :: icall type(plumegen_coms), pointer :: coms ! Set wind effect from namelist @@ -162,19 +164,11 @@ subroutine plumerise(m1,m2,m3,ia,iz,ja,jz, & endif burnt_area= max(1.0e4,burnt_area) - IF (dbg_opt) THEN - WRITE(*,*) 'plumerise: m1 ', m1 - WRITE(*,*) 'plumerise: imm, FRP,burnt_area ', imm, FRP,burnt_area - ! WRITE(*,*) 'convert_smold_to_flam ',convert_smold_to_flam - WRITE(*,*) 'plumerise: zcon ', coms%zcon - WRITE(*,*) 'plumerise: zzcon ', coms%zzcon - END IF - - IF (dbg_opt) then - WRITE(*,*) 'plumerise: imm ', imm - WRITE(*,*) 'plumerise: burnt_area ',burnt_area - END IF - + IF ( dbg_opt .and. (icall .le. n_dbg_lines) .and. (frp_inst .ge. frp_threshold) ) 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 + END IF + !- get fire properties (burned area, plume radius, heating rates ...) call get_fire_properties(coms,imm,iveg_ag,burnt_area,FRP,errmsg,errflg) if(errflg/=0) return @@ -182,8 +176,8 @@ subroutine plumerise(m1,m2,m3,ia,iz,ja,jz, & !------ generates the plume rise ------ call makeplume (coms,kmt,ztopmax(imm),ixx,imm) - IF (dbg_opt) then - WRITE(*,*) 'plumerise after makeplume: imm,kmt,ztopmax(imm) ',imm,kmt,ztopmax(imm) + IF ( dbg_opt .and. (icall .le. n_dbg_lines) .and. (frp_inst .ge. frp_threshold) ) then + WRITE(1000+mpiid,*) 'inside plumerise after makeplume:xlat,xlong,curr_secs,imm,kmt,ztopmax(imm) ', lat, long, int(curr_secs), imm,kmt, ztopmax(imm) END IF enddo lp_minmax @@ -203,12 +197,12 @@ subroutine plumerise(m1,m2,m3,ia,iz,ja,jz, & ! 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 + !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 + !END IF ! enddo lp_veg ! sub-grid vegetation, currently it's aggregated diff --git a/physics/smoke_dust/module_wetdep_ls.F90 b/physics/smoke_dust/module_wetdep_ls.F90 index 8ba8f67d9..2ef07e38c 100755 --- a/physics/smoke_dust/module_wetdep_ls.F90 +++ b/physics/smoke_dust/module_wetdep_ls.F90 @@ -31,11 +31,18 @@ subroutine wetdep_ls(dt,var,rain,moist, real(kind_phys) :: dvar,factor,clsum integer :: nv,i,j,k,km,kb,kbeg !real(kind_phys), parameter :: alpha = .5 ! scavenging factor + integer, save :: print_alpha = 0 wetdpr_smoke =0. wetdpr_dust =0. wetdpr_coarsepm=0. + !if ( print_alpha == 0 ) then + ! write(*,*) 'wetdep_ls, alpha = ',alpha + ! print_alpha = print_alpha + 1 + !endif + + do nv=1,nchem do i=its,ite do j=jts,jte @@ -76,11 +83,11 @@ subroutine wetdep_ls(dt,var,rain,moist, dvar=alpha*factor/(1+factor)*var(i,k,j,nv) ! Accumulate diags if (nv .eq. p_smoke ) then - wetdpr_smoke(i,j) = wetdpr_smoke(i,j) + dvar * rho(i,k,j) * dt * 1.E-6 + wetdpr_smoke(i,j) = wetdpr_smoke(i,j) + dvar * rho(i,k,j) / dt elseif (nv .eq. p_dust_1 ) then - wetdpr_dust(i,j) = wetdpr_dust(i,j) + dvar * rho(i,k,j) * dt * 1.E-6 + wetdpr_dust(i,j) = wetdpr_dust(i,j) + dvar * rho(i,k,j) / dt elseif (nv .eq. p_coarse_pm ) then - wetdpr_coarsepm(i,j) = wetdpr_coarsepm(i,j) + dvar * rho(i,k,j) * dt * 1.E-6 + wetdpr_coarsepm(i,j) = wetdpr_coarsepm(i,j) + dvar * rho(i,k,j) / dt endif var(i,k,j,nv)=max(1.e-16,var(i,k,j,nv)-dvar) endif diff --git a/physics/smoke_dust/plume_data_mod.F90 b/physics/smoke_dust/plume_data_mod.F90 deleted file mode 100755 index 3f0bcdecd..000000000 --- a/physics/smoke_dust/plume_data_mod.F90 +++ /dev/null @@ -1,51 +0,0 @@ -!>\file plume_data_mod.F90 -!! This file contains data for the fire plume rise module. - -module plume_data_mod - - use machine , only : kind_phys - - implicit none - - ! -- FRP parameters - integer, dimension(0:20), parameter :: & - catb = (/ & - 0, & - 2, 1, 2, 1, & !floresta tropical 2 and 4 / extra trop fores 1,3,5 - 2, 3, 3, 3, 3, & !cerrado/woody savanna :6 a 9 - 4, 4, 4, 4, 4, 0, 4, 0, 0, 0, 0 & !pastagem/lavouras: 10 ... - /) - - real(kind=kind_phys), dimension(0:4), parameter :: & - flaming = (/ & - 0.00, & ! - 0.45, & ! % biomass burned at flaming phase : tropical forest igbp 2 and 4 - 0.45, & ! % biomass burned at flaming phase : extratropical forest igbp 1 , 3 and 5 - 0.75, & ! % biomass burned at flaming phase : cerrado/woody savanna igbp 6 to 9 - 0.00 & ! % biomass burned at flaming phase : pastagem/lavoura: igbp 10 a 17 - /) - - real(kind=kind_phys), dimension(0:20), parameter :: & - msize= (/ & - 0.00021, & !0near water,1Evergreen needleleaf,2EvergreenBroadleaf,!3Deciduous Needleleaf,4Deciduous Broadleaf - 0.00021, 0.00021, 0.00021, 0.00021, & !5Mixed forest,6Closed shrublands,7Open shrublands,8Woody savannas,9Savannas, - 0.00023, 0.00022, 0.00022, 0.00022, 0.00029, &! 10Grassland,11Permanent wetlands,12cropland,13'Urban and Built-Up' - 0.00029, 0.00021, 0.00026, 0.00021, 0.00026, &!14cropland/natural vegetation mosaic,15Snow and ice,16Barren or sparsely vegetated - 0.00021, 0.00021, 0.00021, 0.00021, 0.00021, 0.00021 & !17Water,18Wooded Tundra,19Mixed Tundra,20Bare Ground Tundra - /) - - ! -- FRP buffer indices - integer, parameter :: p_frp_hr = 1 - integer, parameter :: p_frp_std = 2 - integer, parameter :: num_frp_plume = 2 - - ! -- plumerise parameters - integer, parameter :: tropical_forest = 1 - integer, parameter :: boreal_forest = 2 - integer, parameter :: savannah = 3 - integer, parameter :: grassland = 4 - integer, parameter :: nveg_agreg = 4 - - public - -end module plume_data_mod diff --git a/physics/smoke_dust/rrfs_smoke_config.F90 b/physics/smoke_dust/rrfs_smoke_config.F90 index d7478986b..c20d6e2db 100755 --- a/physics/smoke_dust/rrfs_smoke_config.F90 +++ b/physics/smoke_dust/rrfs_smoke_config.F90 @@ -24,6 +24,7 @@ module rrfs_smoke_config integer :: addsmoke_flag = 1 integer :: smoke_forecast = 1 integer :: plumerisefire_frq=60 + integer :: n_dbg_lines = 3 integer :: wetdep_ls_opt = 1 integer :: drydep_opt = 1 integer :: pm_settling = 1 diff --git a/physics/smoke_dust/rrfs_smoke_postpbl.meta b/physics/smoke_dust/rrfs_smoke_postpbl.meta index 50fbb4e03..8d7481ec4 100755 --- a/physics/smoke_dust/rrfs_smoke_postpbl.meta +++ b/physics/smoke_dust/rrfs_smoke_postpbl.meta @@ -1,8 +1,7 @@ [ccpp-table-properties] name = rrfs_smoke_postpbl type = scheme - dependencies = ../hooks/machine.F,dep_dry_simple_mod.F90,module_wetdep_ls.F90,dust_data_mod.F90,dust_fengsha_mod.F90,module_add_emiss_burn.F90,module_plumerise.F90,module_smoke_plumerise.F90,module_zero_plumegen_coms.F90,plume_data_mod.F90,rrfs_smoke_config.F90,seas_data_mod.F90,seas_mod.F90,seas_ngac_mod.F90,dep_dry_mod_emerson.F90,dep_data_mod.F90 - + dependencies = dep_dry_simple_mod.F90,module_wetdep_ls.F90,dust_data_mod.F90,dust_fengsha_mod.F90,module_add_emiss_burn.F90,module_plumerise.F90,module_smoke_plumerise.F90,module_zero_plumegen_coms.F90,rrfs_smoke_config.F90,seas_data_mod.F90,seas_mod.F90,seas_ngac_mod.F90,dep_dry_mod_emerson.F90,dep_data_mod.F90 ######################################################################## [ccpp-arg-table] name = rrfs_smoke_postpbl_run diff --git a/physics/smoke_dust/rrfs_smoke_wrapper.F90 b/physics/smoke_dust/rrfs_smoke_wrapper.F90 index eb7f83af6..3842cba54 100755 --- a/physics/smoke_dust/rrfs_smoke_wrapper.F90 +++ b/physics/smoke_dust/rrfs_smoke_wrapper.F90 @@ -12,7 +12,7 @@ module rrfs_smoke_wrapper 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 + p_smoke, p_dust_1, p_coarse_pm, epsilc, n_dbg_lines 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 @@ -52,7 +52,7 @@ subroutine rrfs_smoke_wrapper_init( seas_opt_in, 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 ) + errmsg, errflg, n_dbg_lines_in ) !>-- Namelist @@ -62,7 +62,7 @@ subroutine rrfs_smoke_wrapper_init( seas_opt_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 + integer, intent(in) :: smoke_forecast_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 @@ -100,6 +100,7 @@ subroutine rrfs_smoke_wrapper_init( seas_opt_in, !>-Other extended_sd_diags = extended_sd_diags_in dbg_opt = dbg_opt_in + n_dbg_lines = n_dbg_lines_in end subroutine rrfs_smoke_wrapper_init @@ -111,17 +112,19 @@ end subroutine rrfs_smoke_wrapper_init subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, & u10m, v10m, ustar, rlat, rlon, tskin, pb2d, t2m, dpt2m, & pr3d, ph3d,phl3d, prl3d, tk3d, us3d, vs3d, spechum, w, & - nsoil, smc, vegtype_dom, vegtype_frac, soiltyp, nlcat, & + nsoil, smc, tslb, vegtype_dom, vegtype_frac, soiltyp, nlcat, & dswsfc, zorl, snow, julian,recmol, & idat, rain_cpl, rainc_cpl, hf2d, g, pi, con_cp, con_rd, con_fv, & dust12m_in, emi_ant_in, smoke_RRFS, smoke2d_RRFS, & ntrac, qgrs, gq0, chem3d, tile_num, & ntsmoke, ntdust, ntcoarsepm, imp_physics, imp_physics_thompson, & nwfa, nifa, emanoc, emdust, emseas, drydep_flux_out, wetdpr, & - ebb_smoke_in, frp_output, coef_bb, ebu_smoke,fhist,min_fplume, & - max_fplume, hwp, hwp_ave, wetness, ndvel, ddvel_inout,fire_in, & + ebb_smoke_in, frp_output, coef_bb, fire_type_out, & + 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, & - errmsg,errflg ) + uspdavg, hpbl_thetav, mpicomm, mpirank, mpiroot, errmsg,errflg ) implicit none @@ -135,7 +138,7 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, integer, parameter :: its=1,jts=1,jte=1, kts=1 integer, dimension(:), intent(in) :: land, vegtype_dom, soiltyp - real(kind_phys), dimension(:,:), intent(in) :: smc + real(kind_phys), dimension(:,:), intent(in) :: smc, tslb real(kind_phys), dimension(:,:,:), intent(in) :: dust12m_in real(kind_phys), dimension(:,:,:), intent(in) :: smoke_RRFS real(kind_phys), dimension(:,:), intent(in) :: smoke2d_RRFS @@ -151,16 +154,17 @@ 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) :: fire_in - real(kind_phys), dimension(:), intent(out) :: fire_heat_flux_out, frac_grid_burned_out - real(kind_phys), dimension(:), intent(inout) :: max_fplume, min_fplume - real(kind_phys), dimension(:), intent( out) :: hwp + 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 real(kind_phys), dimension(:), intent(inout) :: hwp_ave real(kind_phys), dimension(:,:), intent(inout) :: nwfa, nifa real(kind_phys), dimension(:,:), intent(inout) :: ddvel_inout real(kind_phys), dimension(:,:), intent(inout) :: drydep_flux_out real(kind_phys), dimension(:,:), intent(inout) :: wetdpr real(kind_phys), dimension(:), intent(in) :: wetness + 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 @@ -187,16 +191,17 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, real(kind_phys), dimension(ims:im, jms:jme) :: ssm, rdrag, uthr, snowh ! fengsha dust real(kind_phys), dimension(ims:im, jms:jme) :: rmol, swdown, znt, clayf, sandf real(kind_phys), dimension(ims:im, nlcat, jms:jme) :: vegfrac - real(kind_phys), dimension(ims:im, nsoil, jms:jme) :: smois + real(kind_phys), dimension(ims:im, nsoil, jms:jme) :: smois, stemp real(kind_phys), dimension(ims:im, 1:1, jms:jme, 1:num_emis_dust) :: emis_dust 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 + fire_end_hr, hwp_day_avg, kpbl_thetav,& + uspdavg2, hpbl_thetav2 integer, dimension(ims:im, jms:jme ) :: min_fplume2, max_fplume2, fire_type - logical :: call_plume, reset_hwp_ave + 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 real(kind_phys), dimension(ims:im, kms:kme, jms:jme, ndvel) :: vgrav @@ -218,6 +223,13 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, real(kind_phys) :: factor, factor2, factor3 integer :: nbegin, nv integer :: i, j, k, kp, n +! MPI variables + integer :: mpiid + integer, intent(in) :: mpicomm + integer, intent(in) :: mpirank + integer, intent(in) :: mpiroot + + mpiid = mpirank errmsg = '' errflg = 0 @@ -232,6 +244,8 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, min_fplume2 = 0 max_fplume2 = 0 + uspdavg2 = 0. + hpbl_thetav2 = 0. emis_seas = 0. emis_dust = 0. peak_hr = 0. @@ -260,12 +274,17 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, ! -- compute incremental convective and large-scale rainfall do i=its,ite rnav(i,1)=max((rain_cpl(i)-rainc_cpl(i))*1000., 0.) ! meter to mm +! coef_bb initializes as clear_val (from GFS_typedefs.F90) +! at ktau = 1, coef_bb_dc is set = 1.0 coef_bb_dc(i,1) = coef_bb(i) +! fhist initializes as 1. (from GFS_typedefs.F90) fire_hist (i,1) = fhist (i) + peak_hr (i,1) = peak_hr_out(i) enddo ! Is this a reset timestep (00:00 + dt)? reset_hwp_ave = mod(int(curr_secs-dt),3600) == 0 + avg_hwp_ave = mod(int(curr_secs),3600) == 0 ! plumerise frequency in minutes set up by the namelist input call_plume = (do_plumerise .and. (plumerisefire_frq > 0)) @@ -276,7 +295,8 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, ktau,current_month, current_hour, gmt, con_rd, con_fv, con_cp, & u10m,v10m,ustar,land,garea,rlat,rlon,tskin, & pr3d,ph3d,phl3d,tk3d,prl3d,us3d,vs3d,spechum,w, & - nsoil,smc,vegtype_dom,soiltyp,nlcat,vegtype_frac,dswsfc,zorl, & + 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, & u10,v10,ust,tsk,xland,xlat,xlong,dxy, & @@ -287,20 +307,14 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, num_chem,num_moist, & ntsmoke, ntdust,ntcoarsepm, & moist,chem,ebu_in,kpbl_thetav,ebb_smoke_in, & - fhist,frp_in, hwp_day_avg, fire_end_hr, & - emis_anoc,smois,ivgtyp,isltyp,vegfrac,rmol,swdown,znt,hfx,pbl, & - snowh,clayf,rdrag,sandf,ssm,uthr,oro, hwp_local, & + fire_hist,frp_in, hwp_day_avg, 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, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) - do j=jts,jte - do i=its,ite - peak_hr(i,j)= fire_in(i,1) - enddo - enddo - IF (ktau==1) THEN ebu = 0. do j=jts,jte @@ -311,6 +325,13 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, enddo enddo enddo + ELSE + do k=kts,kte + do i=its,ite + ! ebu is divided by coef_bb_dc since it is applied in the output + ebu(i,k,1)=ebu_smoke(i,k) / coef_bb_dc(i,1) + enddo + enddo ENDIF !RAR: change this to the fractional LU type; fire_type: 0- no fires, 1- Ag @@ -320,6 +341,7 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, do i=its,ite if (ebu_in(i,j)<0.01) then fire_type(i,j) = 0 + lu_nofire(i,j) = 1.0 else ! 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) @@ -350,9 +372,10 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, endif !-- compute dust (opt=5) - if (dust_opt==5) then - call gocart_dust_fengsha_driver(dt,chem,rho_phy,smois,p8w,ssm, & - isltyp,snowh,xland,dxy,g,emis_dust,ust,znt, & + if (dust_opt==1) then + call gocart_dust_fengsha_driver(dt,chem,rho_phy, & + smois,stemp,p8w,ssm, & + isltyp,snowh,xland,dxy,g,emis_dust,ust,znt, & clayf,sandf,rdrag,uthr, & num_emis_dust,num_chem,nsoil, & ids,ide, jds,jde, kds,kde, & @@ -367,6 +390,7 @@ 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 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) / & @@ -396,7 +420,8 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, kpbl_thetav, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & - its,ite, jts,jte, kts,kte, errmsg, errflg ) + its,ite, jts,jte, kts,kte, errmsg, errflg, curr_secs, & + xlat, xlong, uspdavg2, hpbl_thetav2, mpiid ) if(errflg/=0) return end if @@ -409,7 +434,7 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, swdown,ebb_dcycle,ebu_in,ebu,fire_type, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & - its,ite, jts,jte, kts,kte ) + its,ite, jts,jte, kts,kte , mpiid ) endif !>-- compute coarsepm setting if using simple dry dep option and @@ -431,7 +456,7 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, pm_settling,drydep_flux_local,settling_flux,dbg_opt, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & - its,ite, jts,jte, kts,kte ) + its,ite, jts,jte, kts,kte, curr_secs, mpiid, xlat, xlong ) do nv=1,ndvel do i=its,ite ddvel_inout(i,nv)=ddvel(i,1,nv) @@ -470,10 +495,12 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, endif endif -! Smoke emisisons diagnostic +! Smoke emisisons diagnostic, RAR: let's multiply by coef_bb_dc before output +! Since ebu_smoke includes coef_bb_dc, we need to divide by coef_bb_dc when it +! comes back into the wrapper. do k=kts,kte do i=its,ite - ebu_smoke(i,k)=ebu(i,k,1) + ebu_smoke(i,k)=ebu(i,k,1) * coef_bb_dc(i,1) enddo enddo @@ -485,15 +512,21 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, do i=its,ite hwp(i)=hwp_local(i,1) hwp_ave(i) = hwp_ave(i) + hwp(i)*dt + if ( ktau == 1) then + hwp_ave(i) = hwp_ave(i) / dt + elseif ( avg_hwp_ave ) then + hwp_ave(i) = hwp_ave(i) / 3600._kind_phys + endif enddo - + + !---- diagnostic output of dry deposition & gravitational settling fluxes if ( drydep_opt == 1 .and. (extended_sd_diags .or. dbg_opt) ) then do nv = 1, ndvel do i=its,ite drydep_flux_out(i,nv) = drydep_flux_out(i,nv) + & - drydep_flux_local(i,1,nv) + & - settling_flux(i,1,nv) + drydep_flux_local(i,1,nv) !+ & + !settling_flux(i,1,nv) enddo enddo endif @@ -520,6 +553,7 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, !------------------------------------- !-- to output for diagnostics 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 emanoc (i) = emis_anoc (i) ! anthropogenic organic carbon: ug/m2/s emdust (i) = emis_dust(i,1,1,1) + emis_dust(i,1,1,2) + & @@ -529,10 +563,13 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate, 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 - fire_in(i,1) = peak_hr(i,1) + peak_hr_out(i) = peak_hr(i,1) enddo !-- to provide real aerosol emission for Thompson MP @@ -568,7 +605,7 @@ subroutine rrfs_smoke_prep( & ktau,current_month,current_hour,gmt,con_rd,con_fv,con_cp, & u10m,v10m,ustar,land,garea,rlat,rlon,ts2d, & pr3d,ph3d,phl3d,tk3d,prl3d,us3d,vs3d,spechum,w, & - nsoil,smc,vegtype_dom,soiltyp,nlcat,vegtype_frac,dswsfc,zorl, & + 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, & u10,v10,ust,tsk,xland,xlat,xlong,dxy, & @@ -579,9 +616,9 @@ subroutine rrfs_smoke_prep( & num_chem, num_moist, & ntsmoke, ntdust, ntcoarsepm, & moist,chem,ebu_in,kpbl_thetav,ebb_smoke_in, & - fhist,frp_in, hwp_day_avg, fire_end_hr, & - emis_anoc,smois,ivgtyp,isltyp,vegfrac,rmol,swdown,znt,hfx,pbl, & - snowh,clayf,rdrag,sandf,ssm,uthr,oro,hwp_local, & + fire_hist,frp_in, hwp_day_avg, 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, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & @@ -599,7 +636,7 @@ subroutine rrfs_smoke_prep( & u10m, v10m, ustar, garea, rlat, rlon, ts2d, dswsfc, & zorl, snow_cpl, pb2d, hf2d, oro, t2m, dpt2m, wetness, recmol real(kind=kind_phys), dimension(ims:ime, nlcat), intent(in) :: vegtype_frac - real(kind=kind_phys), dimension(ims:ime, nsoil), intent(in) :: smc + 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 @@ -633,8 +670,8 @@ subroutine rrfs_smoke_prep( & real(kind_phys), dimension(ims:ime, kms:kme, jms:jme, num_chem), intent(out) :: chem 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 - real(kind_phys), dimension(ims:ime,jms:jme), intent(inout) :: frp_in, fire_end_hr, fhist, coef_bb_dc + 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), intent(inout) :: emis_anoc,ebb_smoke_in real(kind_phys), parameter :: conv_frp = 1.e+06_kind_phys ! FRP conversion factor, MW to W @@ -730,6 +767,7 @@ subroutine rrfs_smoke_prep( & do j=jts,jte do i=its,ite smois(i,k,j)=smc(i,k) + stemp(i,k,j)=tslb(i,k) enddo enddo enddo @@ -776,13 +814,9 @@ subroutine rrfs_smoke_prep( & rri(i,k,j)=1./rho_phy(i,k,j) vvel(i,k,j)=-w(i,kkp)*rri(i,k,j)/g moist(i,k,j,:)=0. - moist(i,k,j,1)=gq0(i,kkp,p_atm_shum) - if (t_phy(i,k,j) > 265.) then - moist(i,k,j,2)=gq0(i,kkp,p_atm_cldq) - if (moist(i,k,j,2) < 1.e-8) moist(i,k,j,2)=0. - else - moist(i,k,j,2)=0. - endif + moist(i,k,j,1)=gq0(i,kkp,1) + moist(i,k,j,2)=gq0(i,kkp,2) + if (moist(i,k,j,2) < 1.e-8) moist(i,k,j,2)=0. !-- zmid(i,k,j)=phl3d(i,kkp)/g enddo @@ -862,11 +896,11 @@ subroutine rrfs_smoke_prep( & 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 - fire_end_hr(i,j) = 0.0 - hwp_day_avg(i,j) = 0.0 + ! 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 @@ -890,7 +924,8 @@ subroutine rrfs_smoke_prep( & if (ktau==1) then do j=jts,jte do i=its,ite - fhist (i,j) = 1. + ! 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 peak_hr(i,j)= 0.0* 3600. ! peak at 24 UTC, fires in Alaska @@ -909,7 +944,7 @@ subroutine rrfs_smoke_prep( & enddo endif - ! We will add a namelist variable, real :: flam_frac_global + ! 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 af61ac05e..271d2dd36 100755 --- a/physics/smoke_dust/rrfs_smoke_wrapper.meta +++ b/physics/smoke_dust/rrfs_smoke_wrapper.meta @@ -1,7 +1,7 @@ [ccpp-table-properties] name = rrfs_smoke_wrapper type = scheme - dependencies = ../hooks/machine.F,dep_dry_simple_mod.F90,module_wetdep_ls.F90,dust_data_mod.F90,dust_fengsha_mod.F90,module_add_emiss_burn.F90,module_plumerise.F90,module_smoke_plumerise.F90,module_zero_plumegen_coms.F90,plume_data_mod.F90,rrfs_smoke_config.F90,seas_data_mod.F90,seas_mod.F90,seas_ngac_mod.F90,coarsepm_settling_mod.F90, dep_dry_mod_emerson.F90,dep_data_mod.F90 + dependencies = dep_dry_simple_mod.F90,module_wetdep_ls.F90,dust_data_mod.F90,dust_fengsha_mod.F90,module_add_emiss_burn.F90,module_plumerise.F90,module_smoke_plumerise.F90,module_zero_plumegen_coms.F90,rrfs_smoke_config.F90,seas_data_mod.F90,seas_mod.F90,seas_ngac_mod.F90,coarsepm_settling_mod.F90, dep_dry_mod_emerson.F90,dep_data_mod.F90 ######################################################################## [ccpp-arg-table] @@ -71,6 +71,13 @@ 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 @@ -406,6 +413,14 @@ type = real kind = kind_phys intent = inout +[tslb] + standard_name = soil_temperature_for_land_surface_model + long_name = soil temperature for land surface model + units = K + dimensions = (horizontal_loop_extent,vertical_dimension_of_soil_internal_to_land_surface_scheme) + type = real + kind = kind_phys + intent = in [vegtype_dom] standard_name = vegetation_type_classification long_name = vegetation type at each grid cell @@ -715,6 +730,13 @@ type = real kind = kind_phys intent = inout +[fire_type_out] + standard_name = fire_type + long_name = type of fire + units = 1 + dimensions = (horizontal_loop_extent) + type = integer + intent = out [ebu_smoke] standard_name = ebu_smoke long_name = buffer of vertical fire emission @@ -747,6 +769,43 @@ type = real kind = kind_phys intent = inout +[mpicomm] + standard_name = mpi_communicator + long_name = MPI communicator + units = index + dimensions = () + type = integer + intent = in +[mpirank] + standard_name = mpi_rank + long_name = current MPI-rank + units = index + dimensions = () + type = integer + intent = in +[mpiroot] + standard_name = mpi_root + long_name = master MPI-rank + units = index + dimensions = () + type = integer + intent = in +[uspdavg] + standard_name = mean_wind_speed_in_boundary_layer + long_name = average wind speed within the boundary layer + units = m s-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout +[hpbl_thetav] + standard_name = atmosphere_boundary_layer_thickness_from_modified_parcel + long_name = pbl height based on modified parcel method + units = m + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout [drydep_flux_out] standard_name = dry_deposition_flux long_name = rrfs dry deposition flux @@ -802,14 +861,30 @@ type = real kind = kind_phys intent = inout -[fire_in] - standard_name = smoke_fire_auxiliary_input - long_name = smoke fire auxiliary input variables - units = various - dimensions = (horizontal_loop_extent,fire_auxiliary_data_extent) - type = real - kind = kind_phys - intent = inout +[peak_hr_out] + standard_name = peak_hr_fire + long_name = time_of_peak_fire_emissions + units = s + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out +[lu_nofire_out] + standard_name = sum_of_land_use_fractions_for_no_fire_pixels + long_name = land use of no fire pixels for type + units = 1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out +[lu_qfire_out] + standard_name = sum_of_land_use_fractions_for_cropland_fire_pixels + long_name = land use of fire pixels for type + units = 1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out [fire_heat_flux_out] standard_name = surface_fire_heat_flux long_name = heat flux of fire at the surface