Fix to allow use_theta_m=1 with vertical nesting and simplify rebalan…

…ce options (#1408)

TYPE: bug fix

KEYWORDS: vertical nesting, vert_refine_method, use_theta_m, rebalance

SOURCE: Robert Arthur and Katie Lundquist, LLNL

DESCRIPTION OF CHANGES:
Problem:
1. When vertical nesting is used with the moist theta option (use_theta_m=1), large errors occur in the vertically-
refined domain. This is the main problem addressed in this PR.
2. Another minor issue is that rebalance option 2 is redundant, and is therefore unnecessary.

Solution:
1. Several use_theta_m flags were added to the rebalancing routines in start_em.F and module_dm.F, which are 
called when vertical nesting is used. Thus, the rebalanced pressure, inverse density, and geopotential fields are 
calculated correctly and consistently with the rest of the code, regardless of whether use_theta_m=0 or 1.
2. Rebalance option 2 has been removed. Now, rebalance must be set to 1 when vertical nesting is turned on.

LIST OF MODIFIED FILES: 
Registry/Registry.EM_COMMON
dyn_em/module_initialize_real.F
dyn_em/start_em.F
external/RSL_LITE/module_dm.F
share/module_check_a_mundo.F

TESTS CONDUCTED: 
1. A 4-domain nested case was completed with vertical refinement on d04. All domains were initialized from a 
wrfinput file and eta levels were not provided. Pertinent namelist options are:

input_from_file        = .true., .true., .true., .true.,
e_vert                      =  45,     45,     45,     55,
vert_refine_method = 0,        0,       0,       2,
rebalance                = 1,
hybrid_opt               = 0,
use_theta_m            = 0 or 1, depending on the case

Results are summarized in the plot below, which shows U at the first grid point above the surface on d04 at 
various output times.
![vnest_fix_use_theta_m](https://user-images.githubusercontent.com/46732079/108439600-7304d980-7206-11eb-99e5-269b06f04b1e.png)
With the original code, vertical nesting works with use_theta_m=0, but not use_theta_m=1. With the updated code, 
results look reasonable for use_theta_m=1 and are roughly the same as for use_theta_m=0. Similar patterns are also 
seen for other variables, such as V and MU, although they are not shown here.

2. Several similar tests were also conducted, as above but with
   * d04 initialized via interpolation from d03 during runtime, rather than from wrfinput_d04
   *  eta levels specified for each domain in the namelist

These tests provided qualitatively similar results to the above, showing that the fix works for a range of possible 
domain setups.

3. Finally, a test was completed with use_theta_m=0, but using the updated code. This showed bit-for-bit identical 
results (using diffwrf) to the same case using the original code, confirming that previous vertical nesting functionality 
with use_theta_m=0 is maintained.

RELEASE NOTE: Vertical refinemnet now works with the moist theta option use_theta_m=1. The redundant "rebalance" option 2 was removed. The rebalance option must now be set to 1 when vertical refinement is turned on.

Registry/Registry.EM_COMMON

@@ -2108,7 +2108,7 @@ rconfig   integer num_metgrid_soil_levels namelist,domains	1             4
 rconfig   real    p_top_requested         namelist,domains      1              5000    irh    "p_top_requested" "Pa"      ""
 rconfig   logical interp_theta            namelist,domains	1              .false. irh    "interp_theta"  "inside real, vertically interpolate theta (T) or temperature (F)"  ""
 rconfig   integer interp_type             namelist,domains	1              2       irh    "interp_type"  "1=interp in pressure, 2=interp in LOG pressure"  ""
-rconfig   integer rebalance               namelist,domains      1              0       irh    "rebalance"  "0=no; 1=yes, always; 2=yes, but only when doing vertical nesting"
+rconfig   integer rebalance               namelist,domains      1              0       irh    "rebalance"  "0=no; 1=yes (must be 1 if vertical nesting is used)"
 rconfig   integer vert_refine_method      namelist,domains      max_domains    0       irh    "vert_refine_method"  "0=no vertical nesting, 1=integer refinement, 2=use specified eta levels or compute_eta routine" ""
 rconfig   integer vert_refine_fact        namelist,domains      1              1       irh    "vertical refinment factor for ndown, not used for concurrent vertical grid nesting"  "" ""
 rconfig   integer extrap_type             namelist,domains	1              2       irh    "extrap_type"  "1= use 2 lowest levels, 2=constant"  ""

dyn_em/module_initialize_real.F

@@ -3967,9 +3967,8 @@ SUBROUTINE init_domain_rk ( grid &
       !  Compute pressure similarly to how computed within model, with final Qv.
 
       !  Do a re-balance or not?  0 = NOPE
-
-      IF ( ( config_flags%rebalance .EQ. 0 ) .OR. &
-           ( ( config_flags%rebalance .EQ. 2 ) .AND.  ( config_flags%vert_refine_method .NE. 2 ) ) ) THEN
+      !  Note that rebalance must be 1 for vertical nesting
+      IF ( config_flags%rebalance .EQ. 0 ) THEN
 
           DO j = jts, min(jde-1,jte)
              DO k=kts,kte-1

dyn_em/start_em.F

@@ -691,8 +691,8 @@ SUBROUTINE start_domain_em ( grid, allowed_to_read &
       ENDDO
       ENDDO
       ENDDO
-      IF ( ( config_flags%rebalance .EQ. 1 ) .OR. &
-           ( ( config_flags%rebalance .EQ. 2 ) .AND. ( config_flags%vert_refine_method .EQ. 2 ) ) ) THEN
+      !Note that rebalance must be set to 1 for vertical nesting
+      IF ( config_flags%rebalance .EQ. 1 ) THEN
         !  Integrate base geopotential, starting at terrain elevation.  This assures that
         !  the base state is in exact hydrostatic balance with respect to the model equations.
         !  This field is on full levels.
@@ -773,9 +773,9 @@ SUBROUTINE start_domain_em ( grid, allowed_to_read &
   !  it is an easy decision.
   ELSE
      ! We request rebalancing for vertical grid nesting, or when the user asks for rebalancing.
+     ! (Note that rebalance now must be set to 1 for vertical nesting)
      ! Rebalance recomputes 1/rho, p, ph_2, ph0, p_hyd
-     IF ( ( config_flags%rebalance .EQ. 1 ) .OR. &
-        ( ( config_flags%rebalance .EQ. 2 ) .AND. ( config_flags%vert_refine_method .EQ. 2 ) ) ) THEN
+     IF ( config_flags%rebalance .EQ. 1 ) THEN
        call rebalance_driver_cycl (grid )
        DO j = jts,min(jte,jde-1)
        DO k = kts,kte
@@ -2308,7 +2308,11 @@ SUBROUTINE rebalance_cycl ( grid  &
                qvf1 = qtot*qvf2
 
                grid%p(i,kk,j) = - 0.5*((grid%c1f(k)*grid%Mu_2(i,j))+qvf1*(grid%c1f(k)*grid%Mub(i,j)+grid%c2f(k)))/grid%rdnw(kk)/qvf2
-               qvf = 1.+rvovrd*moist(i,kk,j,P_QV)
+               IF ( grid%use_theta_m == 0) THEN
+                  qvf = 1.+rvovrd*moist(i,kk,j,P_QV)
+               ELSE
+                  qvf = 1.
+               ENDIF
                grid%alt(i,kk,j) = (r_d/p1000mb)*(grid%t_2(i,kk,j)+t0)*qvf*         &
                       (((grid%p(i,kk,j)+grid%pb(i,kk,j))/p1000mb)**cvpm)
                grid%al(i,kk,j) = grid%alt(i,kk,j) - grid%alb(i,kk,j)
@@ -2329,7 +2333,11 @@ SUBROUTINE rebalance_cycl ( grid  &
                qvf1 = qtot*qvf2
                grid%p(i,kk,j) = grid%p(i,kk+1,j) - ((grid%c1f(k)*grid%Mu_2(i,j)) +       &
                                qvf1*(grid%c1f(k)*grid%Mub(i,j)+grid%c2f(k)))/qvf2/grid%rdn(kk+1)
-               qvf = 1. + rvovrd*moist(i,kk,j,P_QV)
+               IF ( grid%use_theta_m == 0) THEN
+                  qvf = 1. + rvovrd*moist(i,kk,j,P_QV)
+               ELSE
+                  qvf = 1.
+               ENDIF
                grid%alt(i,kk,j) = (r_d/p1000mb)*(grid%t_2(i,kk,j)+t0)*qvf* &
                            (((grid%p(i,kk,j)+grid%pb(i,kk,j))/p1000mb)**cvpm)
                grid%al(i,kk,j) = grid%alt(i,kk,j) - grid%alb(i,kk,j)
@@ -2377,7 +2385,11 @@ SUBROUTINE rebalance_cycl ( grid  &
                qvf2 = 1./(1.+qvf1)
                qvf1 = qvf1*qvf2
                grid%p(i,kk,j) = - 0.5*((grid%c1f(k)*grid%Mu_2(i,j))+qvf1*(grid%c1f(k)*grid%Mub(i,j)+grid%c2f(k)))/grid%rdnw(kk)/qvf2
-               qvf = 1. + rvovrd*moist(i,kk,j,P_QV)
+               IF ( grid%use_theta_m == 0) THEN
+                  qvf = 1. + rvovrd*moist(i,kk,j,P_QV)
+               ELSE
+                  qvf = 1.
+               ENDIF
                grid%alt(i,kk,j) = (r_d/p1000mb)*(grid%t_2(i,kk,j)+t0)*qvf    &
                            *(((grid%p(i,kk,j)+grid%pb(i,kk,j))/p1000mb)**cvpm)
                grid%al(i,kk,j) = grid%alt(i,kk,j) - grid%alb(i,kk,j)
@@ -2390,7 +2402,11 @@ SUBROUTINE rebalance_cycl ( grid  &
                qvf2 = 1./(1.+qvf1)
                qvf1 = qvf1*qvf2
                grid%p(i,kk,j) = grid%p(i,kk+1,j) - ((grid%c1f(k)*grid%Mu_2(i,j)) + qvf1*(grid%c1f(k)*grid%Mub(i,j)+grid%c2f(k)))/qvf2/grid%rdn(kk+1)
-               qvf = 1. + rvovrd*moist(i,kk,j,P_QV)
+               IF ( grid%use_theta_m == 0) THEN
+                  qvf = 1. + rvovrd*moist(i,kk,j,P_QV)
+               ELSE
+                  qvf = 1.
+               ENDIF
                grid%alt(i,kk,j) = (r_d/p1000mb)*(grid%t_2(i,kk,j)+t0)*qvf* &
                            (((grid%p(i,kk,j)+grid%pb(i,kk,j))/p1000mb)**cvpm)
                grid%al(i,kk,j) = grid%alt(i,kk,j) - grid%alb(i,kk,j)

external/RSL_LITE/module_dm.F

@@ -4310,7 +4310,11 @@ SUBROUTINE force_domain_em_part2 ( grid, ngrid, pgrid, config_flags    &
                qvf1 = qvf1*qvf2
       
                p(i,kk,j) = - 0.5*((ngrid%c1f(k)*grid%Mu_2(i,j))+qvf1*(ngrid%c1f(k)*grid%Mub(i,j)+ngrid%c2f(k)))/ngrid%rdnw(kk)/qvf2
-               qvf = 1. + rvovrd*moist(i,kk,j,P_QV)
+               IF ( config_flags%use_theta_m == 0) THEN
+                  qvf = 1. + rvovrd*moist(i,kk,j,P_QV)
+               ELSE
+                  qvf = 1.
+               ENDIF
                al(i,kk,j) = (r_d/p1000mb)*(grid%t_2(i,kk,j)+t0)*qvf* &
                            (((p(i,kk,j)+grid%pb(i,kk,j))/p1000mb)**cvpm) - grid%alb(i,kk,j)
       
@@ -4323,7 +4327,11 @@ SUBROUTINE force_domain_em_part2 ( grid, ngrid, pgrid, config_flags    &
                   qvf2 = 1./(1.+qvf1)
                   qvf1 = qvf1*qvf2
                   p(i,kk,j) = p(i,kk+1,j) - ((ngrid%c1f(k)*grid%Mu_2(i,j)) + qvf1*(ngrid%c1f(k)*grid%Mub(i,j)+ngrid%c2f(k)))/qvf2/ngrid%rdn(kk+1)
-                  qvf = 1. + rvovrd*moist(i,kk,j,P_QV)
+                  IF ( config_flags%use_theta_m == 0) THEN
+                     qvf = 1. + rvovrd*moist(i,kk,j,P_QV)
+                  ELSE
+                     qvf = 1.
+                  ENDIF
                   al(i,kk,j) = (r_d/p1000mb)*(grid%t_2(i,kk,j)+t0)*qvf* &
                               (((p(i,kk,j)+grid%pb(i,kk,j))/p1000mb)**cvpm) - grid%alb(i,kk,j)
                END DO

share/module_check_a_mundo.F

@@ -2322,6 +2322,26 @@ END FUNCTION bep_bem_nbui_max
         count_fatal_error = count_fatal_error + 1
       END IF
 
+!-----------------------------------------------------------------------
+!  Consistency checks for vertical refinement:
+!  rebalance must be set to 1 
+!-----------------------------------------------------------------------
+      oops = 0 
+      DO i = 2, model_config_rec % max_dom
+        IF ( .NOT. model_config_rec % grid_allowed(i) ) CYCLE
+        IF (model_config_rec%vert_refine_method(i) .NE. 0) THEN
+          IF ( model_config_rec%rebalance .NE. 1 ) THEN
+            oops = oops + 1
+          END IF
+        END IF
+      END DO
+
+      IF ( oops .GT. 0 ) THEN
+        wrf_err_message = '--- ERROR: vert_refine_method=2 only works with rebalance=1 '
+        CALL wrf_debug ( 0, TRIM( wrf_err_message ) )
+        count_fatal_error = count_fatal_error + 1
+      END IF
+
 !-----------------------------------------------------------------------
 ! This WRF version does not support trajectories on a global domain
 !-----------------------------------------------------------------------