Skip to content

Commit

Permalink
Implicit Explicit Vertical Advection (IEVA) (#1373)
Browse files Browse the repository at this point in the history 
TYPE: new feature

KEYWORDS: IEVA, vertical advection

SOURCE: Louis Wicker (NOAA/NSSL)

DESCRIPTION OF CHANGES:

For grids with large aspect ratios (dx/dz >> 1) that permit explicit convection, the large time step is limited by the 
strongest updraft that occurs during integration.  This results in time step often 20-30% smaller, or requires the use 
of w-filtering, such as latent-heat tendency limiting.  Regions of large vertical velocities are also often very small 
relative to the domain.  The Implicit-Explicit Vertical Advection (IEVA) scheme has been implemented (see Wicker, L. J., 
and W. C. Skamarock, 2020: An Implicit–Explicit Vertical Transport Scheme for Convection-Allowing Models. Mon. Wea. 
Rev., 148, 3893–3910) and that permits a larger time step by partitioning the vertical transport into an explicit piece, 
which uses the normal vertical schemes present in WRF, and a implicit piece which uses implicit transport (which is 
unconditionally stable).  The combined scheme permits a larger time step than has been previously been used and 
reduced w-filtering.
 
The scheme will be useful for CONUS scale CAM (convection allowing model) simulations (dx ~ 2-3 km) when the 
number of vertical levels > 50.  Time steps can increase to as large as 25 s, depending on the problem.  The Wicker 
and Skamarock paper demonstrates IEVA's advantages on the 27 April 2011 Alabama tornado outbreak by comparing 
it to the operational CAM (the High Resolution Rapid Refresh) configuration.  Results are shown that the HRRR 
simulation is stable up to a dt=20 s and only with latent-heat limiting on.  Using the IEVA scheme the dt can be increased 
to 24 s and no latent-heat limiting is needed.  Overall integration efficiency increases ~ 15%, and the IEVA solutions 
are closer to a benchmark run using smaller dt (12 s) than the HRRR simulation.

LIST OF MODIFIED FILES: 
M       Registry/Registry.EM_COMMON
M       dyn_em/Makefile
M       dyn_em/module_advect_em.F
M       dyn_em/module_big_step_utilities_em.F
M       dyn_em/module_em.F
A       dyn_em/module_ieva_em.F
M       dyn_em/solve_em.F
M       run/README.namelist
A       test/em_real/namelist.input.IEVA.4km 

TESTS CONDUCTED:
1. This pull requested code has been tested repeatedly on 27 April for a 24 hour simulation with the parameters set as in the Wicker and Skamarock (2020) as well as a 10 step difference test with the em_quarter_ss case using a single node and 4 nodes.  The differences in the results ~ 10^-5.

2. Jenkins tests are all pass.

RELEASE NOTES: The Implicit-Explicit Vertical Advection (IEVA) scheme has been implemented (see Wicker, L. J., and W. C. Skamarock, 2020: An Implicit–Explicit Vertical Transport Scheme for Convection-Allowing Models. Mon. Wea. Rev., 148, 3893–3910) and that permits a larger time step by partitioning the vertical transport into an explicit piece, which uses the normal vertical schemes present in WRF, and a implicit piece which uses implicit transport (which is unconditionally stable).  The combined scheme permits a larger time step than has been previously been used and reduced w-filtering. The scheme will be useful for CONUS-scale CAM (convection allowing model) simulations (dx ~ 2-3 km) when the number of vertical levels > 50.  In these cases, time steps can increase to as large as 25 s, depending on the problem. Overall integration efficiency increases ~ 15%, and the IEVA solutions are closer to a benchmark run using smaller time step.
  • Loading branch information
@louiswicker
louiswicker authored Mar 22, 2021
1 parent 3f0aae4 commit 4412521
Show file tree
Hide file tree
Showing 9 changed files with 2,206 additions and 542 deletions.
10 changes: 7 additions & 3 deletions Registry/Registry.EM_COMMON
View file
Original file line number Diff line number Diff line change
Expand Up @@ -209,6 +209,8 @@ state real w ikjb dyn_em 2 Z \
state real ww ikj dyn_em 1 Z r "ww" "mu-coupled eta-dot" "Pa s-1"
state real rw ikj dyn_em 1 Z - "rw" "mu-coupled w" "Pa m s-1"
i1 real ww1 ikj dyn_em 1 Z
i1 real wwE ikj dyn_em 1 Z - "wwE" "Explicit vertical velocity" "Pa s-1"
i1 real wwI ikj dyn_em 1 Z - "wwI" "Implicit vertical velocity" "Pa s-1"
state real ww_m ikj dyn_em 1 Z r "ww_m" "time-avg mu-coupled eta-dot" "Pa s-1"
i1 real wwp ikj dyn_em 1 Z
i1 real rw_tend ikj dyn_em 1 Z
Expand Down Expand Up @@ -2772,11 +2774,13 @@ rconfig integer scm_force_flux namelist,scm 1 0 rh
# dynamics option (see package definitions, below)
rconfig integer dyn_opt namelist,dynamics 1 2
rconfig integer rk_ord namelist,dynamics 1 3 irh "rk_order" "" ""
rconfig integer w_damping namelist,dynamics 1 0 irh "w_damping" "" ""
rconfig integer w_damping namelist,dynamics 1 0 irh "w_damping" "" ""
rconfig real w_crit_cfl namelist,dynamics 1 1.2 irh "w_crit_cfl" "W-CFL where w-damping is on" ""
rconfig integer zadvect_implicit namelist,dynamics 1 0 irh "zadvect_implicit" "Turns on IEVA for vertical adv" ""
# diff_opt 1=old diffusion, 2=new
rconfig integer diff_opt namelist,dynamics max_domains -1 irh "diff_opt" "" ""
rconfig integer diff_opt namelist,dynamics max_domains -1 irh "diff_opt" "" ""
# diff_opt_dfi is needed for backwards integration in dfi
rconfig integer diff_opt_dfi namelist,dynamics max_domains 0 irh "diff_opt_dfi" "" ""
rconfig integer diff_opt_dfi namelist,dynamics max_domains 0 irh "diff_opt_dfi" "" ""
# km_opt 1=old coefs, 2=tke, 3=Smagorinksy
rconfig integer km_opt namelist,dynamics max_domains -1 irh "km_opt" "" ""
# km_opt_dfi is needed for backward integration in dfi
Expand Down
19 changes: 10 additions & 9 deletions dyn_em/Makefile
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,26 +7,27 @@ RM = rm -f

MODULES = \
module_advect_em.o \
module_diffusion_em.o \
module_small_step_em.o \
module_ieva_em.o \
module_diffusion_em.o \
module_small_step_em.o \
module_big_step_utilities_em.o \
module_em.o \
module_solvedebug_em.o \
module_bc_em.o \
module_init_utilities.o \
module_wps_io_arw.o \
module_damping_em.o \
module_polarfft.o \
module_damping_em.o \
module_polarfft.o \
module_force_scm.o \
module_first_rk_step_part1.o \
module_first_rk_step_part2.o \
module_avgflx_em.o \
module_sfs_nba.o \
module_sfs_nba.o \
module_convtrans_prep.o \
module_sfs_driver.o \
module_stoch.o \
module_after_all_rk_steps.o \
$(CASE_MODULE)
module_sfs_driver.o \
module_stoch.o \
module_after_all_rk_steps.o \
$(CASE_MODULE)

# possible CASE_MODULE settings
# module_initialize_b_wave.o \
Expand Down
287 changes: 154 additions & 133 deletions dyn_em/module_advect_em.F
View file

Large diffs are not rendered by default.

265 changes: 147 additions & 118 deletions dyn_em/module_big_step_utilities_em.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2501,13 +2501,13 @@ END SUBROUTINE pg_buoy_w
!-------------------------------------------------------------------------------
SUBROUTINE w_damp( rw_tend, max_vert_cfl,max_horiz_cfl, &
u, v, ww, w, mut, c1f, c2f, rdnw, &
rdx, rdy, msfux, msfuy, &
msfvx, msfvy, dt, &
config_flags, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
u, v, ww, w, mut, c1f, c2f, rdnw, &
rdx, rdy, msfux, msfuy, &
msfvx, msfvy, dt, &
config_flags, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
USE module_llxy
IMPLICIT NONE
Expand Down Expand Up @@ -2542,7 +2542,6 @@ SUBROUTINE w_damp( rw_tend, max_vert_cfl,max_horiz_cfl, &
REAL :: vert_cfl, cf_n, cf_d, maxdub, maxdeta
INTEGER :: itf, jtf, i, j, k, maxi, maxj, maxk
INTEGER :: some
CHARACTER*512 :: temp
CHARACTER (LEN=256) :: time_str
Expand All @@ -2551,9 +2550,18 @@ SUBROUTINE w_damp( rw_tend, max_vert_cfl,max_horiz_cfl, &
integer :: total
REAL :: msfuxt , msfxffl
REAL :: w_damp_on = 1.0
REAL :: w_crit_cfl = 2.0
REAL :: w_flag_cfl = 1.2
LOGICAL :: wflags_differ = .false.
LOGICAL :: print_flag = .true.
SAVE :: print_flag
INTEGER :: some1 ! Now have two catagories of CFL information, hence some1 & some2
INTEGER :: some2
!<DESCRIPTION>
!
! w_damp computes a damping term for the vertical velocity when the
! W_damp computes a damping term for the vertical velocity when the
! vertical Courant number is too large. This was found to be preferable to
! decreasing the timestep or increasing the diffusion in real-data applications
! that produced potentially-unstable large vertical velocities because of
Expand All @@ -2564,154 +2572,175 @@ SUBROUTINE w_damp( rw_tend, max_vert_cfl,max_horiz_cfl, &
! horizontal motion. These values are returned via the max_vert_cfl and
! max_horiz_cfl variables. (Added by T. Hutchinson, WSI, 3/5/2007)
!
!-----
!
! W_damp modified to be more flexible with IEVA capability. Two things changed.
! First, the value of the W-CFL where damping is turned on can now be set in the namelist
! using "w_crit_cfl". Second, the code has been modified to report two levels of
! W-CFL information. Previously, both damping and W-CFL information was
! printed at a W-CFL == 1.2. This new capability does not need IEVA on to run.
!
! We keep the old level reporting (for consistency), but use "w_flag_cfl" to flag
! and print information from those points. "w_flag_cfl" is set here in w_damp.
! A second level of print occurs, which is useful for when we run
! IEVA, because many points can run with W-CFL > 1.2. We now print this
! second level separately, as these points might make the model blow up.
! and where we now want Rayleigh damping turned on. The default value "w_crit_cfl"
! which replaced "w_beta" in the old code, is set in the namelist/Registery
! turn on w_damp at 1.2 consistent with the WRFV3/4 (before this) code.
!
! Summary
! -------
! W_damp will write out CFL information when the vertical courant number is > w_flag_cfl,
! and W_damp will also turn on Rayleigh damping if vertical courant number > w_crit_cfl.
!
! Typical settings for non-IEVA use: w_crit_cfl = 1.2
! Typical settings for IEVA use: w_crit_cfl = 2.0
!
! I also commented out a lot of code put in 8-13 years ago for timing assuming that
! this is no longer needed....there is some pretty hacky stuff. HTH.
!
! (Added by L. Wicker, NSSL, 5/4/2020)
!
!</DESCRIPTION>
itf=MIN(ite,ide-1)
jtf=MIN(jte,jde-1)
some = 0
max_vert_cfl = 0.
some1 = 0
some2 = 0
max_vert_cfl = 0.
max_horiz_cfl = 0.
total = 0
w_crit_cfl = config_flags%w_crit_cfl
IF( abs(w_crit_cfl - w_flag_cfl) > 0.1 ) THEN
wflags_differ = .true.
ELSE
wflags_differ = .false.
ENDIF
IF( print_flag ) THEN
write(wrf_err_message,*) '----------------------------------------'
CALL wrf_debug( 0, wrf_err_message )
WRITE(temp,*) 'W-DAMPING BEGINS AT W-COURANT NUMBER = ',w_crit_cfl
CALL wrf_debug ( 0 , TRIM(temp) )
write(wrf_err_message,*) '----------------------------------------'
CALL wrf_debug( 0, wrf_err_message )
print_flag = .false.
ENDIF
IF(config_flags%polar ) then
msfxffl = 1.0/COS(config_flags%fft_filter_lat*degrad)
END IF
IF ( config_flags%w_damping == 1 ) THEN
! Routine has been reorganized to hopefully reduce redundant code while maintaining efficiency
#ifdef OPTIMIZE_CFL_TEST
! 20121025, L. Meadows vector optimization does not include special case for Cassini
IF(config_flags%polar ) then
CALL wrf_error_fatal('module_big_step_utilities_em.F: -DOPTIMIZE_CFL_TEST option does not support global domains')
END IF
IF(config_flags%polar ) then
CALL wrf_error_fatal('module_big_step_utilities_em.F: -DOPTIMIZE_CFL_TEST option does not support global domains')
END IF
#endif
DO j = jts,jtf
DO j = jts,jtf
DO k = 2,kde-1
DO i = its,itf
vert_cfl = abs(ww(i,k,j)/(c1f(k)*mut(i,j)+c2f(k))*rdnw(k)*dt)
DO k = 2, kde-1
DO i = its,itf
#if 1
# ifdef OPTIMIZE_CFL_TEST
! L. Meadows, Intel, MIC optimization, 20121025
msfuxt = msfux(i,j)
vert_cfl = abs(ww(i,k,j)/(c1f(k)*mut(i,j)+c2f(k))*rdnw(k)*dt)
IF ( vert_cfl > max_vert_cfl ) THEN
max_vert_cfl = vert_cfl
ENDIF
msfuxt = msfux(i,j)
IF ( vert_cfl > max_vert_cfl ) THEN
max_vert_cfl = vert_cfl
ENDIF
# else
IF(config_flags%polar ) then
msfuxt = MIN(msfux(i,j), msfxffl)
ELSE
msfuxt = msfux(i,j)
END IF
vert_cfl = abs(ww(i,k,j)/(c1f(k)*mut(i,j)+c2f(k))*rdnw(k)*dt)
IF(config_flags%polar ) THEN
msfuxt = MIN(msfux(i,j), msfxffl)
ELSE
msfuxt = msfux(i,j)
ENDIF
IF ( vert_cfl > max_vert_cfl ) THEN
max_vert_cfl = vert_cfl ; maxi = i ; maxj = j ; maxk = k
maxdub = w(i,k,j) ; maxdeta = -1./rdnw(k)
ENDIF
IF ( vert_cfl > max_vert_cfl ) THEN
max_vert_cfl = vert_cfl ; maxi = i ; maxj = j ; maxk = k
maxdub = w(i,k,j) ; maxdeta = -1./rdnw(k)
ENDIF
# endif
horiz_cfl = max( abs(u(i,k,j) * rdx * msfuxt * dt), &
abs(v(i,k,j) * rdy * msfvy(i,j) * dt) )
if (horiz_cfl > max_horiz_cfl) then
max_horiz_cfl = horiz_cfl
endif
horiz_cfl = max( abs(u(i,k,j) * rdx * msfuxt * dt), &
abs(v(i,k,j) * rdy * msfvy(i,j) * dt) )
IF (horiz_cfl > max_horiz_cfl) THEN
max_horiz_cfl = horiz_cfl
ENDIF
! Dump out more information without affecting performance
if(vert_cfl .gt. w_beta)then
#else
! restructure to get rid of divide
!
! This had been used for efficiency, but with the addition of returning the cfl values,
! the old version (above) was reinstated. (T. Hutchinson, 3/5/2007)
!
cf_n = abs(ww(i,k,j)*rdnw(k)*dt)
cf_d = abs((c1f(k)*mut(i,j)+c2f(k)))
if(cf_n .gt. cf_d*w_beta )then
#endif
#ifndef OPTIMIZE_CFL_TEST
! This internal write is costly on newer Xeon processors because it breaks
! vectorization. (J. Michalakes for L. Meadows at Intel, 12/13/2012)
WRITE(temp,*)i,j,k,' vert_cfl,w,d(eta)=',vert_cfl,w(i,k,j),-1./rdnw(k)
IF (vert_cfl .gt. w_crit_cfl) THEN
some1 = some1 + 1
WRITE(temp,FMT="(3(1x,i5,1x),'W-CRIT_CFL: ',f7.4,2x,'W-CFL: ',f7.4,2x,'dETA: ',f7.4))") &
i,j,k,vert_cfl,w(i,k,j),-1./rdnw(k)
CALL wrf_debug ( 100 , TRIM(temp) )
#endif
if ( vert_cfl > 2. ) some = some + 1
rw_tend(i,k,j) = rw_tend(i,k,j)-sign(1.,w(i,k,j))*w_alpha*(vert_cfl-w_beta)*(c1f(k)*mut(i,j)+c2f(k))
endif
END DO
ENDDO
ENDDO
ELSE
! just print
DO j = jts,jtf
DO k = 2, kde-1
DO i = its,itf
ENDIF
#if 1
# ifdef OPTIMIZE_CFL_TEST
msfuxt = msfux(i,j)
vert_cfl = abs(ww(i,k,j)/(c1f(k)*mut(i,j)+c2f(k))*rdnw(k)*dt)
IF ( vert_cfl > max_vert_cfl ) THEN
max_vert_cfl = vert_cfl
ENDIF
# else
! L. Meadows MIC optimization, 20121025
IF(config_flags%polar ) then
msfuxt = MIN(msfux(i,j), msfxffl)
ELSE
msfuxt = msfux(i,j)
END IF
vert_cfl = abs(ww(i,k,j)/(c1f(k)*mut(i,j)+c2f(k))*rdnw(k)*dt)
IF ( vert_cfl > max_vert_cfl ) THEN
max_vert_cfl = vert_cfl ; maxi = i ; maxj = j ; maxk = k
maxdub = w(i,k,j) ; maxdeta = -1./rdnw(k)
ENDIF
# endif
horiz_cfl = max( abs(u(i,k,j) * rdx * msfuxt * dt), &
abs(v(i,k,j) * rdy * msfvy(i,j) * dt) )
IF ((vert_cfl .gt. w_flag_cfl) .and. wflags_differ) THEN
if (horiz_cfl > max_horiz_cfl) then
max_horiz_cfl = horiz_cfl
endif
if(vert_cfl .gt. w_beta)then
#else
! restructure to get rid of divide
!
! This had been used for efficiency, but with the addition of returning the cfl values,
! the old version (above) was reinstated. (T. Hutchinson, 3/5/2007)
!
cf_n = abs(ww(i,k,j)*rdnw(k)*dt)
cf_d = abs((c1f(k)*mut(i,j)+c2f(k)))
if(cf_n .gt. cf_d*w_beta )then
#endif
#ifndef OPTIMIZE_CFL_TEST
WRITE(temp,*)i,j,k,' vert_cfl,w,d(eta)=',vert_cfl,w(i,k,j),-1./rdnw(k)
some2 = some2 + 1
WRITE(temp,FMT="(3(1x,i5,1x),'W-FLAG_CFL: ',f7.4,2x,'W-CFL: ',f7.4,2x,'dETA: ',f7.4))") &
i,j,k,vert_cfl,w(i,k,j),-1./rdnw(k)
CALL wrf_debug ( 100 , TRIM(temp) )
ENDIF
#endif
if ( vert_cfl > 2. ) some = some + 1
endif
END DO
ENDDO
ENDDO
ENDIF
IF ( some .GT. 0 ) THEN
IF ((vert_cfl .gt. w_crit_cfl) .and. (config_flags%w_damping == 1) ) THEN
rw_tend(i,k,j) = rw_tend(i,k,j)-sign(1.,w(i,k,j))*w_alpha*(vert_cfl-w_crit_cfl)*(c1f(k)*mut(i,j)+c2f(k))
ENDIF
ENDDO ! end i-loop
ENDDO ! end k-loop
ENDDO ! end j-loop
IF ( some1 .GT. 0 ) THEN
CALL get_current_time_string( time_str )
CALL get_current_grid_name( grid_str )
WRITE(temp,*)some, &
' points exceeded cfl=2 in domain '//TRIM(grid_str)//' at time '//TRIM(time_str)//' hours'
WRITE(temp,*) some1, &
' points exceeded W_CRITICAL_CFL in domain '//TRIM(grid_str)//' at time '//TRIM(time_str)//' hours'
CALL wrf_debug ( 0 , TRIM(temp) )
#ifndef OPTIMIZE_CFL_TEST
WRITE(temp,*)'MAX AT i,j,k: ',maxi,maxj,maxk,' vert_cfl,w,d(eta)=',max_vert_cfl, &
maxdub,maxdeta
WRITE(temp,FMT="('Max W: ',3(1x,i5,1x),'W: ',f7.2,2x,'W-CFL: ',f7.2,2x,'dETA: ',f7.2))") &
maxi, maxj, maxk, maxdub, max_vert_cfl, maxdeta
CALL wrf_debug ( 0 , TRIM(temp) )
WRITE(temp,FMT="('Max U/V: ',3(1x,i5,1x),'U: ',f7.2,2x,'U-CFL: ',f7.2,2x,'V: ',f7.2,2x,'V-CFL: ',f7.2))") &
maxi, maxj, maxk, u(maxi,maxk,maxj), dt*u(maxi,maxk,maxj)*rdx, v(maxi,maxk,maxj), dt*v(maxi,maxk,maxj)*rdy
CALL wrf_debug ( 0 , TRIM(temp) )
ENDIF
IF ( some2 .GT. 0 ) THEN
CALL get_current_time_string( time_str )
CALL get_current_grid_name( grid_str )
WRITE(temp,*) some2, &
' points exceeded W_FLAG_CFL in domain '//TRIM(grid_str)//' at time '//TRIM(time_str)//' hours'
CALL wrf_debug ( 0 , TRIM(temp) )
#endif
ENDIF
END SUBROUTINE w_damp
END SUBROUTINE W_DAMP
!-------------------------------------------------------------------------------
Expand Down
Loading

0 comments on commit 4412521

Please sign in to comment.