From 15669394a867e09b1cb5aefd01aec7ea46f70854 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Matthias=20G=C3=B6bel?= Date: Sun, 13 Dec 2020 19:50:08 +0100 Subject: [PATCH 1/2] avoid double staggering of ww in vertical advection of geopotential --- dyn_em/module_big_step_utilities_em.F | 5 ++--- dyn_em/module_small_step_em.F | 7 +++---- 2 files changed, 5 insertions(+), 7 deletions(-) diff --git a/dyn_em/module_big_step_utilities_em.F b/dyn_em/module_big_step_utilities_em.F index eceac268da..20cb6738ca 100644 --- a/dyn_em/module_big_step_utilities_em.F +++ b/dyn_em/module_big_step_utilities_em.F @@ -1458,8 +1458,7 @@ SUBROUTINE rhs_ph( ph_tend, u, v, ww, & 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)) + 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 @@ -1468,7 +1467,7 @@ SUBROUTINE rhs_ph( ph_tend, u, v, ww, & 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)) + - ww(i,k,j)*(fnm(k)*wdwn(i,k+1)+fnp(k)*wdwn(i,k)) ENDDO ENDDO diff --git a/dyn_em/module_small_step_em.F b/dyn_em/module_small_step_em.F index 6802822f76..bc9a1f2771 100644 --- a/dyn_em/module_small_step_em.F +++ b/dyn_em/module_small_step_em.F @@ -1315,8 +1315,7 @@ 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)) + 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)) rhs(i,k+1) = dts*(ph_tend(i,k+1,j) + .5*g*(1.-epssm)*w(i,k+1,j)) ENDDO ENDDO @@ -1325,8 +1324,8 @@ SUBROUTINE advance_w( w, rw_tend, ww, w_save, u, v, & ! 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 ) ) + rhs(i,k) = rhs(i,k)-dts*ww(i,k,j)*( fnm(k)*wdwn(i,k+1) & + +fnp(k)*wdwn(i,k ) ) ENDDO ENDDO ! NOTE: phi'' is not coupled with the map-scale factor (1/m), From 9daec2f6b46b728115b16423461f8839de57a22c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Matthias=20G=C3=B6bel?= Date: Fri, 12 Feb 2021 07:34:26 +0100 Subject: [PATCH 2/2] include namelist switch to control the vertical advection of phi --- Registry/Registry.EM_COMMON | 1 + dyn_em/module_big_step_utilities_em.F | 44 +++++++++++++++++++-------- dyn_em/module_small_step_em.F | 40 +++++++++++++++++++----- run/README.namelist | 1 + 4 files changed, 67 insertions(+), 19 deletions(-) diff --git a/Registry/Registry.EM_COMMON b/Registry/Registry.EM_COMMON index a6c2d8f771..98133bea96 100644 --- a/Registry/Registry.EM_COMMON +++ b/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" "" diff --git a/dyn_em/module_big_step_utilities_em.F b/dyn_em/module_big_step_utilities_em.F index 20cb6738ca..2e5bfe4104 100644 --- a/dyn_em/module_big_step_utilities_em.F +++ b/dyn_em/module_big_step_utilities_em.F @@ -1454,22 +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) = 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) & - - ww(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 diff --git a/dyn_em/module_small_step_em.F b/dyn_em/module_small_step_em.F index bc9a1f2771..6b455a1ada 100644 --- a/dyn_em/module_small_step_em.F +++ b/dyn_em/module_small_step_em.F @@ -1315,19 +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)=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*ww(i,k,j)*( 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: diff --git a/run/README.namelist b/run/README.namelist index 397245df1e..5d472521d3 100644 --- a/run/README.namelist +++ b/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.