"MYNN, GF, RUC LSM and smoke plumerise updates for RRFSv1 code freeze"

physics/CONV/Grell_Freitas/cu_gf_deep.F90

@@ -425,9 +425,9 @@ subroutine cu_gf_deep_run(        &
      integer :: turn,pmin_lev(its:ite),start_level(its:ite),ktopkeep(its:ite)
      real(kind=kind_phys),    dimension (its:ite,kts:kte) :: dtempdz
      integer, dimension (its:ite,kts:kte) ::  k_inv_layers 
-     real(kind=kind_phys),    dimension (its:ite) :: c0    ! HCB
+     real(kind=kind_phys),    dimension (its:ite) :: c0, rrfs_factor  ! HCB
      real(kind=kind_phys),    dimension (its:ite,kts:kte) :: c0t3d    ! hli for smoke/dust wet scavenging
-!$acc declare create(pmin_lev,start_level,ktopkeep,dtempdz,k_inv_layers,c0,c0t3d)
+!$acc declare create(pmin_lev,start_level,ktopkeep,dtempdz,k_inv_layers,c0,rrfs_factor,c0t3d)
 
 ! rainevap from sas
      real(kind=kind_phys) zuh2(40)
@@ -486,6 +486,7 @@ subroutine cu_gf_deep_run(        &
 ! Set cloud water to rain water conversion rate (c0)
 !$acc kernels
       c0(:)=0.004
+      rrfs_factor(:)=1.
       do i=its,itf
          xland1(i)=int(xland(i)+.0001) ! 1.
          if(xland(i).gt.1.5 .or. xland(i).lt.0.5)then
@@ -495,6 +496,7 @@ subroutine cu_gf_deep_run(        &
          if(imid.eq.1)then
            c0(i)=0.002
          endif
+         if(kdt.le.(4500./dtime))rrfs_factor(i)=1.-(float(kdt)/(4500./dtime)-1.)**2
       enddo
 !$acc end kernels
 
@@ -591,7 +593,6 @@ 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
@@ -2029,7 +2030,7 @@ subroutine cu_gf_deep_run(        &
             zuo,pre,pwo_ens,xmb,ktop,                                    &
             edto,pwdo,'deep',ierr2,ierr3,                                &
             po_cup,pr_ens,maxens3,                                       &
-            sig,closure_n,xland1,xmbm_in,xmbs_in,                        &
+            sig,closure_n,xland1,xmbm_in,xmbs_in,rrfs_factor,            &
             ichoice,imid,ipr,itf,ktf,                                    &
             its,ite, kts,kte,                                            &
             dicycle,xf_dicycle )
@@ -4056,7 +4057,7 @@ subroutine cup_output_ens_3d(xff_mid,xf_ens,ierr,dellat,dellaq,dellaqc,  &
               zu,pre,pw,xmb,ktop,                                           &
               edt,pwd,name,ierr2,ierr3,p_cup,pr_ens,                        &
               maxens3,                                                      &
-              sig,closure_n,xland1,xmbm_in,xmbs_in,                         &
+              sig,closure_n,xland1,xmbm_in,xmbs_in,rrfs_factor,             &
               ichoice,imid,ipr,itf,ktf,                                     &
               its,ite, kts,kte,                                             &
               dicycle,xf_dicycle )
@@ -4118,7 +4119,7 @@ subroutine cup_output_ens_3d(xff_mid,xf_ens,ierr,dellat,dellaq,dellaqc,  &
         ,intent (inout)                   ::                           &
         ierr,ierr2,ierr3
      integer, intent(in) :: dicycle
-     real(kind=kind_phys),    intent(in), dimension (its:ite) :: xf_dicycle
+     real(kind=kind_phys),    intent(in), dimension (its:ite) :: xf_dicycle, rrfs_factor
 !$acc declare copyin(zu,pwd,p_cup,sig,xmbm_in,xmbs_in,edt,xff_mid,dellat,dellaqc,dellaq,pw,ktop,xland1,xf_dicycle)
 !$acc declare copy(xf_ens,pr_ens,outtem,outq,outqc,pre,xmb,closure_n,ierr,ierr2,ierr3)
 !
@@ -4198,6 +4199,7 @@ subroutine cup_output_ens_3d(xff_mid,xf_ens,ierr,dellat,dellaq,dellaqc,  &
          clos_wei=16./max(1.,closure_n(i))
          xmb_ave(i)=min(xmb_ave(i),100.)
          xmb(i)=clos_wei*sig(i)*xmb_ave(i)
+         if(dx(i)<dx_thresh) xmb(i)=rrfs_factor(i)*xmb(i)
 
            if(xmb(i) < 1.e-16)then
               ierr(i)=19

physics/Interstitials/UFS_SCM_NEPTUNE/GFS_MP_generic_post.F90

@@ -21,11 +21,11 @@ module GFS_MP_generic_post
       subroutine GFS_MP_generic_post_run(                                                                                 &
         im, levs, kdt, nrcm, nncl, ntcw, ntrac, imp_physics, imp_physics_gfdl, imp_physics_thompson, imp_physics_nssl,    &
         imp_physics_mg, imp_physics_fer_hires, cal_pre, cplflx, cplchm, cpllnd, progsigma, con_g, rhowater, rainmin, dtf, &
-        frain, rainc, rain1, rann, xlat, xlon, gt0, gq0, prsl, prsi, phii, tsfc, ice, phil, htop, refl_10cm,              & 
+        frain, rainc, rain1, rann, xlat, xlon, gt0, gq0, prsl, prsi, phii, tsfc, ice, phil, htop, refl_10cm,maxupmf,xland,&
         imfshalcnv,imfshalcnv_gf,imfdeepcnv,imfdeepcnv_gf,imfdeepcnv_samf, con_t0c, snow, graupel, save_t, save_q,        &
         rain0, ice0, snow0, graupel0, del, rain, domr_diag, domzr_diag, domip_diag, doms_diag, tprcp, srflag, sr, cnvprcp,&
         totprcp, totice, totsnw, totgrp, cnvprcpb, totprcpb, toticeb, totsnwb, totgrpb, rain_cpl, rainc_cpl, snow_cpl,    &
-        pwat, frzr, frzrb, frozr, frozrb, tsnowp, tsnowpb, rhonewsn1, exticeden,                                          & 
+        pwat, frzr, frzrb, frozr, frozrb, tsnowp, tsnowpb, rhonewsn1, exticeden,                                          &
         drain_cpl, dsnow_cpl, lsm, lsm_ruc, lsm_noahmp, raincprv, rainncprv, iceprv, snowprv,                             &
         graupelprv, draincprv, drainncprv, diceprv, dsnowprv, dgraupelprv, dtp, dfi_radar_max_intervals,                  &
         dtend, dtidx, index_of_temperature, index_of_process_mp,ldiag3d, qdiag3d,dqdt_qmicro, lssav, num_dfi_radar,       &
@@ -42,15 +42,15 @@ subroutine GFS_MP_generic_post_run(
       logical, intent(in) :: cal_pre, lssav, ldiag3d, qdiag3d, cplflx, cplchm, cpllnd, progsigma, exticeden
       integer, intent(in) :: index_of_temperature,index_of_process_mp,use_lake_model(:)
       integer, intent(in) :: imfshalcnv,imfshalcnv_gf,imfdeepcnv,imfdeepcnv_gf,imfdeepcnv_samf
-      integer, dimension (:), intent(in) :: htop
+      integer, dimension (:), intent(in) :: htop, xland
       integer                                                :: dfi_radar_max_intervals
       real(kind=kind_phys),                    intent(in)    :: fh_dfi_radar(:), fhour, con_t0c
       real(kind=kind_phys),                    intent(in)    :: radar_tten_limits(:)
       integer                                                :: ix_dfi_radar(:)
       real(kind=kind_phys), dimension(:,:),    intent(inout) :: gt0,refl_10cm
 
       real(kind=kind_phys),                    intent(in)    :: dtf, frain, con_g, rainmin, rhowater
-      real(kind=kind_phys), dimension(:),      intent(in)    :: rain1, xlat, xlon, tsfc
+      real(kind=kind_phys), dimension(:),      intent(in)    :: rain1, xlat, xlon, tsfc, maxupmf
       real(kind=kind_phys), dimension(:),      intent(inout) :: ice, snow, graupel, rainc
       real(kind=kind_phys), dimension(:),      intent(in)    :: rain0, ice0, snow0, graupel0
       real(kind=kind_phys), dimension(:,:),    intent(in)    :: rann
@@ -171,6 +171,9 @@ subroutine GFS_MP_generic_post_run(
                  fctz = 10.**(factor(i)*delz)
                endif
                cuprate = rainc(i) * 3.6e6 / dtp  ! cu precip rate (mm/h)
+               if (imfdeepcnv==imfdeepcnv_gf .and. xland(i)==0)then
+                 if( maxupmf(i).lt.0.1 .or. cuprate.lt.0.05) cuprate=0.
+               endif
                ze_conv = 300.0 * cuprate**1.4
                ze_conv = fctz * ze_conv
                ze_mp = 10._kind_phys ** (0.1 * refl_10cm(i,k))

physics/Interstitials/UFS_SCM_NEPTUNE/GFS_MP_generic_post.meta

@@ -277,6 +277,21 @@
   type = real
   kind = kind_phys
   intent = inout
+[maxupmf]
+  standard_name = maximum_convective_updraft_mass_flux
+  long_name = maximum convective updraft mass flux within a column
+  units = m s-1
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = in
+[xland]
+  standard_name = sea_land_ice_mask
+  long_name = landmask: sea/land/ice=0/1/2
+  units = flag
+  dimensions = (horizontal_loop_extent)
+  type = integer
+  intent = in
 [imfshalcnv]
   standard_name = control_for_shallow_convection_scheme
   long_name = flag for mass-flux shallow convection scheme

physics/PBL/MYNN_EDMF/module_bl_mynn.F90

@@ -2000,7 +2000,7 @@ SUBROUTINE  mym_length (                     &
         ugrid = sqrt(u1(kts)**2 + v1(kts)**2)
         uonset= 15. 
         wt_u = (1.0 - min(max(ugrid - uonset, 0.0)/30.0, 0.5)) 
-        cns  = 2.7 !was 3.5
+        cns  = 3.5
         alp1 = 0.23
         alp2 = 0.3
         alp3 = 2.5 * wt_u !taper off bouyancy enhancement in shear-driven pbls
@@ -2034,7 +2034,7 @@ 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.03 ), 30.0)*dzk
+           qdz = min(max( qkw(k)-qmin, 0.02 ), 30.0)*dzk
            elt = elt +qdz*zwk
            vsc = vsc +qdz
            k   = k+1
@@ -5031,7 +5031,7 @@ SUBROUTINE mynn_tendencies(kts,kte,i,       &
     IF (FLAG_QI) THEN
       DO k=kts,kte
          Dth(k)=(thl(k) + xlvcp/exner(k)*sqc2(k)          &
-           &            + xlscp/exner(k)*(sqi2(k)+sqs(k)) &
+           &            + xlscp/exner(k)*(sqi2(k))        & !+sqs(k)) &
            &            - th(k))/delt
          !Use form from Tripoli and Cotton (1981) with their
          !suggested min temperature to improve accuracy:
@@ -6052,7 +6052,7 @@ SUBROUTINE DMP_mf(                            &
     if ((landsea-1.5).LT.0) then  !land
        acfac = .5*tanh((fltv2 - 0.02)/0.05) + .5
     else                          !water
-       acfac = .5*tanh((fltv2 - 0.01)/0.03) + .5
+       acfac = .5*tanh((fltv2 - 0.015)/0.04) + .5
     endif
     !add a windspeed-dependent adjustment to acfac that tapers off
     !the mass-flux scheme linearly above sfc wind speeds of 10 m/s.

physics/SFC_Models/Land/RUC/lsm_ruc.F90

@@ -1304,8 +1304,7 @@ subroutine lsm_ruc_run                                            & ! inputs
 
         ! --- ... accumulated total runoff and surface runoff
         runoff(i)  = runoff(i)  + (drain(i)+runof(i)) * delt  ! accum total kg m-2
-        !srunoff(i) = srunoff(i) + runof(i) * delt             ! accum surface kg m-2
-        srunoff(i) = acrunoff(i,j)        ! accum surface kg m-2
+        srunoff(i) = srunoff(i) + runof(i) * delt             ! accum surface kg m-2
 
         ! --- ... accumulated frozen precipitation (accumulation in lsmruc)
         snowfallac_lnd(i) = snfallac_lnd(i,j) ! accum kg m-2

physics/SFC_Models/Land/RUC/module_sf_ruclsm.F90

@@ -1740,7 +1740,7 @@ SUBROUTINE SFCTMP (debug_print, delt,ktau,conflx,i,j,         & !--- input varia
          !-- will reduce warm bias in western Canada
          !-- and US West coast, where max snow albedo is low (0.3-0.5).
            !print *,'ALB increase to 0.7',alb_snow,snhei,snhei_crit,albsn,i,j
-           !ALBsn = 0.7_kind_phys
+           ALBsn = 0.7_kind_phys
          endif
 
          Emiss= emissn
@@ -1753,7 +1753,7 @@ SUBROUTINE SFCTMP (debug_print, delt,ktau,conflx,i,j,         & !--- input varia
          !-- will reduce warm bias in western Canada
          !-- and US West coast, where max snow albedo is low (0.3-0.5).
            !print *,'ALB increase to 0.7',alb_snow,snhei,snhei_crit,albsn,i,j
-           !ALBsn = 0.7_kind_phys
+           ALBsn = 0.7_kind_phys
            !print *,'NO mosaic ALB increase to 0.7',alb_snow,snhei,snhei_crit,alb,i,j
          endif
 

physics/smoke_dust/module_smoke_plumerise.F90

@@ -169,12 +169,20 @@ subroutine plumerise(m1,m2,m3,ia,iz,ja,jz,                      &
             WRITE(1000+mpiid,*) 'inside plumerise: xlat,xlong,curr_secs,imm,FRP,burnt_area ', lat, long, int(curr_secs), imm, FRP,burnt_area
         END IF
 
+       IF (frp_inst<frp_threshold) THEN
+         k1=1
+         k2=2
+         !exit
+         return
+       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
 
+
        !------  generates the plume rise    ------
-       call makeplume (coms,kmt,ztopmax(imm),ixx,imm)
+       call makeplume (coms,kmt,ztopmax(imm),ixx,imm,mpiid)
 
        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)
@@ -562,7 +570,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)  
+SUBROUTINE MAKEPLUME (coms,kmt,ztopmax,ixx,imm,mpiid)
 !
 ! *********************************************************************
 !
@@ -621,22 +629,23 @@ SUBROUTINE MAKEPLUME (coms,kmt,ztopmax,ixx,imm)
 !
 !
 !**********************************************************************
-!**********************************************************************               
-!use module_zero_plumegen_coms 
-implicit none 
-!logical :: endspace  
+!**********************************************************************
+!use module_zero_plumegen_coms
+implicit none
+!logical :: endspace
 type(plumegen_coms), pointer :: coms
 character (len=10) :: varn
 integer ::  izprint, iconv,  itime, k, kk, kkmax, deltak,ilastprint,kmt &
            ,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,
-character (len=2) :: cixx 
+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()
     REAL(kind=kind_phys) :: DELZ_THRESOLD = 100. 
 
-    INTEGER     :: imm
+    INTEGER     :: imm, dtknt
 
 !  real(kind=kind_phys), external:: esat_pr!
 !
@@ -654,6 +663,7 @@ SUBROUTINE MAKEPLUME (coms,kmt,ztopmax,ixx,imm)
 coms%viscosity = 500.!- coms%viscosity constant (original value: 0.001)
 
 nrectotal=150
+dtknt = 0
 !
 !*************** PROBLEM SETUP AND INITIAL CONDITIONS *****************
 coms%mintime = 1  
@@ -697,9 +707,13 @@ SUBROUTINE MAKEPLUME (coms,kmt,ztopmax,ixx,imm)
 !sam 81  format('nm1=',I0,' from kmt=',I0,' kkmax=',I0,' deltak=',I0)
 !sam     write(0,81) coms%nm1,kmt,kkmax,deltak
 !-- set timestep
-    !coms%dt = (coms%zm(2)-coms%zm(1)) / (coms%tstpf * wmax)  
-    coms%dt = min(5.,(coms%zm(2)-coms%zm(1)) / (coms%tstpf * wmax))
-
+    !coms%dt = (coms%zm(2)-coms%zm(1)) / (coms%tstpf * wmax) i
+    coms%dt = max(0.01,min(5.,(coms%zm(2)-coms%zm(1)) / (coms%tstpf * wmax)))
+    dtknt = dtknt + 1
+!    if (coms%dt .ne. 5.)then
+!    WRITE(1000+mpiid,*) 'dtknt,zm(2),zm(1) ', dtknt,coms%zm(2),coms%zm(1)
+!    WRITE(1000+mpiid,*) 'coms%tstpf,wmax,dt =', coms%tstpf,wmax,coms%dt
+!    endif
 !-- elapsed time, sec
     coms%time = coms%time+coms%dt 
 !-- elapsed time, minutes                                      

physics/smoke_dust/rrfs_smoke_wrapper.F90

@@ -123,7 +123,7 @@ 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,                     &
+                   fire_heat_flux_out, frac_grid_burned_out, kpbl,oro, totprcp,            &
                    uspdavg, hpbl_thetav, mpicomm, mpirank, mpiroot, errmsg,errflg          )
 
     implicit none
@@ -145,7 +145,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
@@ -329,7 +329,7 @@ subroutine rrfs_smoke_wrapper_run(im, kte, kme, ktau, dt, garea, land, jdate,
       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)
+        ebu(i,k,1)=ebu_smoke(i,k) / MAX(1.E-4,coef_bb_dc(i,1))
       enddo
       enddo
     ENDIF
@@ -734,6 +734,9 @@ subroutine rrfs_smoke_prep(                                               &
     moist          = 0._kind_phys  
     chem           = 0._kind_phys
     z_at_w         = 0._kind_phys
+    if ( ebb_dcycle == 1 ) then
+       coef_bb_dc  = 1._kind_phys
+    endif
 
     do i=its,ite
      u10  (i,1)=u10m (i)

physics/smoke_dust/rrfs_smoke_wrapper.meta

@@ -916,6 +916,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