Merge branch 'production/RRFS.v1' into thompson_refactor_RRFS

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/CONV/Grell_Freitas/cu_gf_driver.F90

@@ -880,6 +880,13 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,&
                  cutenm(i)=0.
               endif   ! pret > 0
 
+              maxupmf(i)=0.
+              if(forcing2(i,6).gt.0.)then
+                maxupmf(i)=maxval(xmb(i)*zu(i,kts:ktf)/forcing2(i,6))
+              endif
+              if (xland(i)==0)then ! cu precip rate (mm/h)
+                 if((maxupmf(i).lt.0.1) .or. (pret(i)*3600.lt.0.05)) pret(i)=0.
+              endif
               if(pret(i).gt.0.)then
                  cuten(i)=1.
                  cutenm(i)=0.
@@ -996,10 +1003,6 @@ subroutine cu_gf_driver_run(ntracer,garea,im,km,dt,flag_init,flag_restart,&
             gdc(i,15,10)=qfx(i)
             gdc(i,16,10)=pret(i)*3600.
 
-            maxupmf(i)=0.
-            if(forcing2(i,6).gt.0.)then
-              maxupmf(i)=maxval(xmb(i)*zu(i,kts:ktf)/forcing2(i,6))
-            endif
 
             if(ktop(i).gt.2 .and.pret(i).gt.0.)dt_mf(i,ktop(i)-1)=ud_mf(i,ktop(i))
             endif

physics/GWD/drag_suite.F90

@@ -1363,7 +1363,8 @@ subroutine drag_suite_run(                                           &
             DO k=kts,km
                wsp=SQRT(uwnd1(i,k)**2 + vwnd1(i,k)**2)
                ! alpha*beta*Cmd*Ccorr*2.109 = 12.*1.*0.005*0.6*2.109 = 0.0759
-               var_temp = 0.0759*EXP(-(zl(i,k)/H_efold)**1.5)*a2*          &
+               ! Change alpha to 35 -- 0.0759 becomes 0.2214
+               var_temp = 0.2214*EXP(-(zl(i,k)/H_efold)**1.5)*a2*          &
                                  zl(i,k)**(-1.2)*ss_taper(i) ! this is greater than zero
                !  Note:  This is a semi-implicit treatment of the time differencing
                !  per Beljaars et al. (2004, QJRMS)

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

@@ -109,12 +109,6 @@ subroutine plumerise(m1,m2,m3,ia,iz,ja,jz,                      &
 ! print *,' Plumerise_scalar 1',ncall
   coms => get_thread_coms()
 
-IF (frp_inst<frp_threshold) THEN
-   k1=1
-   k2=2
-   !return
-END IF
-
 ! print *,' Plumerise_scalar 2',m1
   j=1
   i=1
@@ -169,12 +163,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 +564,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 +623,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 +657,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 +701,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

@@ -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
     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)