Skip to content

Commit

Permalink
Avoid double staggering of ww in vertical advection of geopotential (w…
Browse files Browse the repository at this point in the history
…rf-model#1338)

TYPE: enhancement

KEYWORDS: geopotential,vertical advection

SOURCE: Matthias Göbel (University of Innsbruck)

DESCRIPTION OF CHANGES:
The geopotential equation in WRF is formulated in the following (advective) form where staggering/destaggering is 
marked with overbars:

<img src="https://latex.codecogs.com/svg.latex?\partial_{t}\phi^{\prime}&plus;\mu_{d}^{-1}[m_{x}m_{y}\left(\overline{\overline{U}^\eta\:\delta_{x}\phi}^x&plus;\overline{\overline{V}^\eta\:\delta_{y}\phi}^y\right)&plus;m_{y}\overline{\overline{\Omega}^\eta\:\delta_{\eta}\phi}^\eta-m_{y}gW]=0" title="http://latex.codecogs.com/svg.latex?\partial_{t}\phi^{\prime}+\mu_{d}^{-1}[m_{x}m_{y}\left(\overline{\overline{U}^\eta\:\delta_{x}\phi}^x+\overline{\overline{V}^\eta\:\delta_{y}\phi}^y\right)+m_{y}\overline{\overline{\Omega}^\eta\:\delta_{\eta}\phi}^\eta-m_{y}gW]=0" />

The vertical advection term thus employs a double staggering of Omega. First, Omega is destaggered to mass levels 
and multiplied with the phi gradient. Then the product is staggered back to w-levels. This formulation is analogous to the horizontal terms, but the double averaging of Omega may introduce additional inaccuracies and is not necessary. 

The phi gradient can be first staggered to w-levels and then multiplied with the unaveraged Omega:
 <img src="https://latex.codecogs.com/svg.latex?\overline{\overline{\Omega}^\eta\:\delta_{\eta}\phi}^\eta\rightarrow\Omega\:\overline{\delta_{\eta}\phi}^\eta"  />
For consistency and to avoid numerical instabilities, the same thing has to be done in the acoustic step.

The new namelist option `phi_adv_z` allows switching between the two formulations: 
`phi_adv_z = 1`: old formulation (default)
`phi_adv_z = 2`: new formulation

LIST OF MODIFIED FILES: 
Registry/Registry.EM_COMMON
dyn_em/module_big_step_utilities_em.F
dyn_em/module_small_step_em.F
run/README.namelist 

TESTS CONDUCTED: 
Idealized test cases em_les and em_hill2d_x. 
The vertical velocity is affected by the modification. Updrafts and downdrafts in the em_hill2d_x case become ~10% 
stronger with the new formulation:
![W_xz](https://user-images.githubusercontent.com/17001470/107921009-243d1280-6f6e-11eb-97d0-a897d2d78701.png)
In the LES case stronger updrafts and downdrafts occur more frequently:
![hist_w](https://user-images.githubusercontent.com/17001470/108021629-3c209f00-701f-11eb-80c2-5640201690a0.png)

I came across this issue when trying to transform advective components of temperature, humidity, and momentum 
into Cartesian form in a consistent way, using w for the vertical advection. Without the proposed modifications I 
could not close the budget, because the vertical advection of phi was not consistent with the vertical advection 
in the other budget equations.


RELEASE NOTE: A new namelist option is included that allows a user to avoiding double stagger averaging of Omega in the vertical advection of geopotential. The new namelist option phi_adv_z allows switching between the two formulations: phi_adv_z = 1 (the old formulation, which remains the default), and phi_adv_z = 2 (the new formulation).
  • Loading branch information
matzegoebel authored and vlakshmanan-scala committed Apr 4, 2024
1 parent 11a6680 commit a1f4f73
Show file tree
Hide file tree
Showing 4 changed files with 67 additions and 21 deletions.
1 change: 1 addition & 0 deletions Registry/Registry.EM_COMMON
Original file line number Diff line number Diff line change
Expand Up @@ -2738,6 +2738,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" ""
Expand Down
45 changes: 32 additions & 13 deletions dyn_em/module_big_step_utilities_em.F
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
41 changes: 33 additions & 8 deletions dyn_em/module_small_step_em.F
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down
1 change: 1 addition & 0 deletions run/README.namelist
Original file line number Diff line number Diff line change
Expand Up @@ -1543,6 +1543,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.
Expand Down

0 comments on commit a1f4f73

Please sign in to comment.