Avoid double staggering of ww in vertical advection of geopotential

Registry/Registry.EM_COMMON

@@ -2662,6 +2662,7 @@ rconfig   integer  chem_adv_opt           namelist,dynamics	max_domains    1
 rconfig   integer  tracer_adv_opt         namelist,dynamics     max_domains    1       rh    "tracer_adv_opt"        "positive-definite RK3 transport switch"      ""
 rconfig   integer  scalar_adv_opt         namelist,dynamics	max_domains    1       rh    "scalar_adv_opt"        "positive-definite RK3 transport switch"      ""
 rconfig   integer  tke_adv_opt            namelist,dynamics	max_domains    1       rh    "tke_adv_opt"           "positive-definite RK3 transport switch"      ""
+rconfig   integer  phi_adv_z              namelist,dynamics	max_domains    1       h     "phi_adv_z"             "Vertical advection option for geopotential; 1: original 2: avoid double staggering of omega"   ""
 # switches for selectively deactivating 2nd and 6th order horizontal filters for specific scalar variable classes
 rconfig   logical  moist_mix2_off         namelist,dynamics	max_domains    .false. rh    "moist_mix2_off"        "de-activate 2nd-order horizontal mixing for moisture"      ""
 rconfig   logical  chem_mix2_off          namelist,dynamics	max_domains    .false. rh    "chem_mix2_off"         "de-activate 2nd-order horizontal mixing for chem species"  ""

dyn_em/module_big_step_utilities_em.F

@@ -1454,23 +1454,42 @@ SUBROUTINE rhs_ph( ph_tend, u, v, ww,               &
 
 ! advective form for the geopotential equation
 
+   !  RHS term 3 is: - omega partial d/dnu(phi)
+
    DO j = jts, jtf
+      IF (config_flags%phi_adv_z == 2) THEN
+         !  First get staggered partial d/dnu(phi) and then multiply with omega
+         !  to avoid double staggering of omega
 
-     DO k = 2, kte
-     DO i = its, itf
-          wdwn(i,k) = .5*(ww(i,k,j)+ww(i,k-1,j))*rdnw(k-1)               &
-                        *(ph(i,k,j)-ph(i,k-1,j)+phb(i,k,j)-phb(i,k-1,j))
-     ENDDO
-     ENDDO
+         DO k = 2, kte
+         DO i = its, itf
+              wdwn(i,k) = rdnw(k-1)*(ph(i,k,j)-ph(i,k-1,j)+phb(i,k,j)-phb(i,k-1,j))
+         ENDDO
+         ENDDO
+
+         DO k = 2, kte-1
+         DO i = its, itf
+              ph_tend(i,k,j) = ph_tend(i,k,j)                           &
+                                - ww(i,k,j)*(fnm(k)*wdwn(i,k+1)+fnp(k)*wdwn(i,k))
+         ENDDO
+         ENDDO
+      ELSE
+         !  First destagger omega and multiply with partial d/dnu(phi), then stagger the product
 
-!  RHS term 3 is: - omega partial d/dnu(phi)
+         DO k = 2, kte
+         DO i = its, itf
+              wdwn(i,k) = .5*(ww(i,k,j)+ww(i,k-1,j))*rdnw(k-1)               &
+                            *(ph(i,k,j)-ph(i,k-1,j)+phb(i,k,j)-phb(i,k-1,j))
+         ENDDO
+         ENDDO
 
-     DO k = 2, kte-1
-     DO i = its, itf
-           ph_tend(i,k,j) = ph_tend(i,k,j)                           &
-                             - (fnm(k)*wdwn(i,k+1)+fnp(k)*wdwn(i,k))
-     ENDDO
-     ENDDO
+         DO k = 2, kte-1
+         DO i = its, itf
+               ph_tend(i,k,j) = ph_tend(i,k,j)                           &
+                                 - (fnm(k)*wdwn(i,k+1)+fnp(k)*wdwn(i,k))
+         ENDDO
+         ENDDO
+     ENDIF
 
    ENDDO
 

dyn_em/module_small_step_em.F

@@ -1315,20 +1315,45 @@ SUBROUTINE advance_w( w, rw_tend, ww, w_save, u, v,  &
                     +(1.-epssm)*t_2ave(i,k,j))
       t_2ave(i,k,j)=(t_2ave(i,k,j) + (c1h(k)*Muave(i,j))*t0) &
                     /((c1h(k)*Muts(i,j)+c2h(k))*(t0+t_1(i,k,j)))
-      wdwn(i,k+1)=.5*(ww(i,k+1,j)+ww(i,k,j))*rdnw(k)    &
-           *(ph_1(i,k+1,j)-ph_1(i,k,j)+phb(i,k+1,j)-phb(i,k,j))
       rhs(i,k+1) = dts*(ph_tend(i,k+1,j) + .5*g*(1.-epssm)*w(i,k+1,j))
     ENDDO
     ENDDO
 !   building up RHS of phi equation
 !   ph_tend contains terms 1 and 2; now adding 3 and 4 in stages:
 !   here rhs = delta t [ph_tend + ~g*w/2 - ~ww * partial d phi/dz]
-    DO k=2,k_end
-    DO i=i_start, i_end
-       rhs(i,k) = rhs(i,k)-dts*( fnm(k)*wdwn(i,k+1)  &
-                                +fnp(k)*wdwn(i,k  ) )
-    ENDDO
-    ENDDO
+    IF (config_flags%phi_adv_z == 2) THEN
+      !  First get staggered partial d/dnu(phi) and then multiply with omega
+      !  to avoid double staggering of omega
+
+      DO k=1, k_end
+      DO i=i_start, i_end
+         wdwn(i,k+1)=rdnw(k)*(ph_1(i,k+1,j)-ph_1(i,k,j)+phb(i,k+1,j)-phb(i,k,j))
+      ENDDO
+      ENDDO
+
+      DO k=2,k_end
+      DO i=i_start, i_end
+         rhs(i,k) = rhs(i,k)-dts*ww(i,k,j)*( fnm(k)*wdwn(i,k+1)  &
+                                            +fnp(k)*wdwn(i,k  ) )
+      ENDDO
+      ENDDO
+    ELSE
+      !  First destagger omega and multiply with partial d/dnu(phi), then stagger the product
+
+      DO k=1, k_end
+      DO i=i_start, i_end
+         wdwn(i,k+1)=.5*(ww(i,k+1,j)+ww(i,k,j))*rdnw(k)    &
+              *(ph_1(i,k+1,j)-ph_1(i,k,j)+phb(i,k+1,j)-phb(i,k,j))
+      ENDDO
+      ENDDO
+
+      DO k=2,k_end
+      DO i=i_start, i_end
+         rhs(i,k) = rhs(i,k)-dts*( fnm(k)*wdwn(i,k+1)  &
+                                  +fnp(k)*wdwn(i,k  ) )
+      ENDDO
+      ENDDO
+    ENDIF
 !  NOTE:  phi'' is not coupled with the map-scale factor  (1/m),
 !         but it's tendency is, so must multiply by msft here
 !   Comments on map scale factors:

run/README.namelist

@@ -1517,6 +1517,7 @@ The following are for observation nudging:
  chem_adv_opt (max_dom)              = 1        ! for chem variables
  tracer_adv_opt (max_dom)            = 1        ! for tracer variables (WRF-Chem activated)
  tke_adv_opt (max_dom)               = 1        ! for tke
+ phi_adv_z (max_dom)                 = 1        ! vertical advection option for geopotential; 1: original (default) 2: avoid double staggering of omega
  moist_mix2_off (max_dom)            = .false.  ! if set to T, deactivate 2nd-order horizontal mixing for moisture. default is F.
  chem_mix2_off (max_dom)             = .false.  ! if set to T, deactivate 2nd-order horizontal mixing for chem species. default is F.
  tracer_mix2_off (max_dom)           = .false.  ! if set to T, deactivate 2nd-order horizontal mixing for tracers. default is F.