From cfbc6444cdff1876f422213ed611fed5f22ca3cb Mon Sep 17 00:00:00 2001 From: Julio Bacmeister Date: Sat, 31 Oct 2026 10:12:21 -0600 Subject: [PATCH 1/8] Fix moving-mountain GWD momentum conservation: sign calc_taucd stress with ubm at the source level Passes calc_taucd a corrected copy of ubi with the source-interface element replaced by ubm(src_level), so the momentum flux sign matches the drag gw_drag_prof actually applied. Previously the midpoint-interpolated ubi(src_level+1) could carry the opposite sign, making momentum_fixer double the column momentum change instead of cancelling it. gw_drag_prof and the UBI_MOVMTN diagnostic are unchanged, so the legacy source-level critical-level behavior is preserved bit-for-bit. Co-Authored-By: Claude Opus 5 (1M context) --- .../gravity_wave_drag_moving_mountain.F90 | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/schemes/gravity_wave_drag/gravity_wave_drag_moving_mountain.F90 b/schemes/gravity_wave_drag/gravity_wave_drag_moving_mountain.F90 index 9e74bbffb..47686c6f6 100644 --- a/schemes/gravity_wave_drag/gravity_wave_drag_moving_mountain.F90 +++ b/schemes/gravity_wave_drag/gravity_wave_drag_moving_mountain.F90 @@ -294,8 +294,9 @@ subroutine gravity_wave_drag_moving_mountain_run(ncol, pver, pcnst, & integer, intent(out) :: errflg ! Local variables - integer :: k, m + integer :: k, m, i real(kind_phys) :: xpwp_clubb(ncol, pver + 1) + real(kind_phys) :: ubi_src(ncol, pver + 1) ! Reynolds stress for waves propagating in each cardinal direction. real(kind_phys) :: taucd(ncol, pver + 1, 4) @@ -355,8 +356,13 @@ subroutine gravity_wave_drag_moving_mountain_run(ncol, pver, pcnst, & ttgw, qtgw, egwdffi, gwut, dttdf, dttke, & lapply_effgw_in=gw_apply_tndmax) + ! Adjust ubi to yield correct directions when ubi(src)>0>ubi(src+1) + ubi_src(:ncol,:) = ubi(:ncol,:) + do i=1,ncol + ubi_src(i,src_level(i)+1) = ubm(i ,src_level(i)) + end do ! Project stress into directional components. - taucd = calc_taucd(ncol, band%ngwv, tend_level, tau, phase_speeds, xv, yv, ubi) + taucd = calc_taucd(ncol, band%ngwv, tend_level, tau, phase_speeds, xv, yv, ubi_src) taucd_west(:,:pver+1) = taucd(:,:pver+1,west) taucd_east(:,:pver+1) = taucd(:,:pver+1,east) From d0d2c01432043e65f3466a49e13431362d4c363f Mon Sep 17 00:00:00 2001 From: Julio Bacmeister Date: Sat, 31 Oct 2026 10:12:21 -0600 Subject: [PATCH 2/8] Add roo tilt-dominance source blocking with terrain-following centroid bounds to moving-mountain GWD Adds the 'roo' (near-surface tilt-dominance) criterion to tilt_centroid_levels, with terrain-following centroid bounds scaled by local surface pressure instead of a fixed 800 hPa cutoff, and a roo_block output flagging columns whose source was zeroed for that reason. Threads zm through the source routines. Only gw_movmtn_source 3, 4 and 5 run the centroid search, so the default source type is unaffected. Co-Authored-By: Claude Opus 5 (1M context) --- .../gravity_wave_drag_moving_mountain.F90 | 173 ++++++++++++++---- .../gravity_wave_drag_moving_mountain.meta | 6 + 2 files changed, 140 insertions(+), 39 deletions(-) diff --git a/schemes/gravity_wave_drag/gravity_wave_drag_moving_mountain.F90 b/schemes/gravity_wave_drag/gravity_wave_drag_moving_mountain.F90 index 47686c6f6..978574225 100644 --- a/schemes/gravity_wave_drag/gravity_wave_drag_moving_mountain.F90 +++ b/schemes/gravity_wave_drag/gravity_wave_drag_moving_mountain.F90 @@ -210,7 +210,7 @@ subroutine gravity_wave_drag_moving_mountain_run(ncol, pver, pcnst, & ttgw, qtgw, egwdffi_tot, dttdf, dttke, & tau0, gwut0, & usteer, vsteer, CS, steer_level, xpwp_src, & - tilt, p_steer, p_launch, & + tilt, p_steer, p_launch, roo_block, & taucd_west, taucd_east, taucd_south, taucd_north, & errmsg, errflg) @@ -285,6 +285,7 @@ subroutine gravity_wave_drag_moving_mountain_run(ncol, pver, pcnst, & real(kind_phys), intent(out) :: tilt(:, :) ! Vortex tilting magnitude |zeta|*|dV/dz| [s-2] real(kind_phys), intent(out) :: p_steer(:) ! Pressure at steering level [Pa] real(kind_phys), intent(out) :: p_launch(:) ! Pressure at launch level [Pa] + real(kind_phys), intent(out) :: roo_block(:) ! 1=column blocked by near-surface tilt-dominance (roo) criterion, 0=not blocked real(kind_phys), intent(out) :: taucd_west(:, :) ! Reynolds stress for waves in W direction, interfaces [N m-2] real(kind_phys), intent(out) :: taucd_east(:, :) ! Reynolds stress for waves in E direction, interfaces [N m-2] real(kind_phys), intent(out) :: taucd_south(:, :) ! Reynolds stress for waves in S direction, interfaces [N m-2] @@ -339,7 +340,7 @@ subroutine gravity_wave_drag_moving_mountain_run(ncol, pver, pcnst, & tau, ubm, ubi, xv, yv, & phase_speeds, hdepth, use_gw_movmtn_pbl, rair, gravit, & usteer, vsteer, CS, steer_level, xpwp_src, & - tilt, p_steer, p_launch, & + tilt, p_steer, p_launch, roo_block, & errmsg, errflg) if (errflg /= 0) return @@ -410,7 +411,7 @@ subroutine gw_movmtn_src(ncol, pver, & src_level, tend_level, tau, ubm, ubi, xv, yv, & c, hdepth, use_gw_movmtn_pbl, rair, gravit, & usteer, vsteer, CS, steer_level, xpwp_src, & - tilt, p_steer, p_launch, & + tilt, p_steer, p_launch, roo_block, & errmsg, errflg) !------------------------------------------------------------------------ ! Driver for gravity wave source from obstacle effects produced by @@ -473,6 +474,11 @@ subroutine gw_movmtn_src(ncol, pver, & real(kind_phys), intent(out) :: tilt(:,:) ! Pressure at steering/launch level [Pa]. real(kind_phys), intent(out) :: p_steer(:), p_launch(:) + ! Flag: 1 if this column's source was zeroed because the near-surface + ! tilt-dominance (roo) criterion blocked the centroid-based search (see + ! tilt_centroid_levels); 0 otherwise. Only ever 1 for source_type=3/4/5; + ! source_type=1/2 use fixed levels and never run that search. + real(kind_phys), intent(out) :: roo_block(:) character(len=512), intent(out) :: errmsg integer, intent(out) :: errflg @@ -524,6 +530,7 @@ subroutine gw_movmtn_src(ncol, pver, & q0 = 0._kind_phys p_steer = 0._kind_phys p_launch = 0._kind_phys + roo_block = 0._kind_phys !---------------------------------------------------------------------- ! source_type>=3 deposits xpwp_src directly at the launch level, which @@ -564,16 +571,16 @@ subroutine gw_movmtn_src(ncol, pver, & end do case (3) - call tilt_dyn_src(tilt, u, v, pmid, delp, ncol, pver, alpha_gw_movmtn, & - xpwp_src, Steer_k, Launch_k, p_steer, p_launch, usteer, vsteer) + call tilt_dyn_src(tilt, u, v, pmid, delp, zm, ncol, pver, alpha_gw_movmtn, & + xpwp_src, Steer_k, Launch_k, p_steer, p_launch, usteer, vsteer, roo_block) case (4) - call tilt_precip_6param_src(tilt, prect, u, v, pmid, ncol, pver, alpha_gw_movmtn, & - xpwp_src, Steer_k, Launch_k, p_steer, p_launch, usteer, vsteer) + call tilt_precip_6param_src(tilt, prect, u, v, pmid, zm, ncol, pver, alpha_gw_movmtn, & + xpwp_src, Steer_k, Launch_k, p_steer, p_launch, usteer, vsteer, roo_block) case (5) - call tilt_precip_cx17_src(tilt, prect, u, v, pmid, delp, ncol, pver, alpha_gw_movmtn, & - xpwp_src, Steer_k, Launch_k, p_steer, p_launch, usteer, vsteer) + call tilt_precip_cx17_src(tilt, prect, u, v, pmid, delp, zm, ncol, pver, alpha_gw_movmtn, & + xpwp_src, Steer_k, Launch_k, p_steer, p_launch, usteer, vsteer, roo_block) case default errflg = 1 @@ -756,8 +763,8 @@ end subroutine gw_movmtn_src !========================================================================== - subroutine tilt_centroid_levels(tilt, u, v, pmid, ncol, pver, & - steering_level, launch_level, p_steer, p_launch, usteer, vsteer, valid) + subroutine tilt_centroid_levels(tilt, u, v, pmid, zm, ncol, pver, & + steering_level, launch_level, p_steer, p_launch, usteer, vsteer, valid, roo_block) !------------------------------------------------------------------------ ! Flow-dependent steering and launch levels shared by source_type 3-5. ! @@ -767,33 +774,94 @@ subroutine tilt_centroid_levels(tilt, u, v, pmid, ncol, pver, & ! valid=.false. and levels fall back to safe values so downstream ! code is never handed an out-of-bounds index. Callers must zero ! the source where valid=.false. + ! 3. roo criterion: the centroid search can "succeed" (return a + ! well-defined in-bounds pair of levels) and still be meaningless, + ! when a column's tilt profile is dominated by an isolated, extreme + ! near-surface spike rather than a genuine free-tropospheric signal + ! -- the regime observed over Antarctica, the Arctic Ocean, and + ! Greenland. roo = max(|tilt|, 0-500 m AGL) / max(|tilt|, 500-15000 m + ! AGL); roo > roo_thresh means the near-surface layer so overwhelms + ! the layer above it that the steering/launch search cannot be + ! trusted, so the column is blocked (valid=.false., same safe + ! fallback levels as a failed centroid search) regardless of whether + ! the search itself returned in-bounds levels. roo_block records + ! which columns were blocked for this reason specifically (as + ! opposed to a plain centroid-search failure). ! ! CAM level ordering: k=1 TOA, k=pver surface, pressure increases ! with k. Therefore launch_level < steering_level in index space - ! (launch is higher altitude = lower pressure = smaller k). + ! (launch is higher altitude = lower pressure = smaller k). zm is + ! height AGL (m), decreasing with k, consistent with its use in + ! compute_tilt. !------------------------------------------------------------------------ integer, intent(in) :: ncol, pver real(kind_phys), intent(in) :: tilt(ncol,pver) real(kind_phys), intent(in) :: u(ncol,pver), v(ncol,pver) real(kind_phys), intent(in) :: pmid(ncol,pver) + real(kind_phys), intent(in) :: zm(ncol,pver) integer, intent(out) :: steering_level(ncol), launch_level(ncol) real(kind_phys), intent(out) :: p_steer(ncol), p_launch(ncol) real(kind_phys), intent(out) :: usteer(ncol), vsteer(ncol) logical, intent(out) :: valid(ncol) - - ! Pressure bounds for centroid search [Pa]: 0 (TOA) to 800 hPa (excl. BL). - real(kind_phys), parameter :: p_min_centroid = 0._kind_phys - real(kind_phys), parameter :: p_max_centroid = 80000._kind_phys - + real(kind_phys), intent(out) :: roo_block(ncol) + + ! Pressure bounds for centroid search [Pa]. + ! + ! Top bound: fixed at TOA (unchanged). + ! + ! Bottom (BL-exclusion) bound: was a single fixed 800 hPa, which is + ! terrain-blind. Over elevated terrain (e.g. the Antarctic plateau, + ! Andes, Himalaya/Tibetan Plateau, Rockies -- anywhere local surface + ! pressure is already below ~800 hPa) that fixed cutoff is never + ! reached, so no level is ever excluded and the steering-level + ! centroid is free to lock onto raw near-surface shear/tilt noise + ! (which is exactly where this parameterization's local eta-level + ! vorticity calc and polar grid-metric terms are most contaminated). + ! + ! Fix: exclude the bottom (1 - frac_pbl_excl) of each column's own + ! pressure range instead, i.e. scale the cutoff by local surface + ! pressure (pmid(:,pver), the lowest model layer midpoint, used as + ! a close proxy for true surface pressure -- using the true + ! interface/surface pressure instead would require threading it in + ! from the caller, which is not done here to keep the change local). + real(kind_phys), parameter :: p_min_centroid = 0._kind_phys + real(kind_phys), parameter :: frac_pbl_excl = 0.8_kind_phys ! fraction of local p_sfc + + ! roo criterion bounds: a fixed, terrain-following (AGL, not + ! pressure) comparison, since the regimes where this fires are + ! precisely where a fixed-pressure near-surface layer would span + ! wildly different physical depths depending on surface elevation. + real(kind_phys), parameter :: z_low_roo = 500._kind_phys ! m AGL + real(kind_phys), parameter :: z_mid_top_roo = 15000._kind_phys ! m AGL + real(kind_phys), parameter :: roo_thresh = 5._kind_phys ! try 10 for a stricter trigger + real(kind_phys), parameter :: roo_eps = 1.e-7_kind_phys ! ~typical tilt scale, avoids /0 + + real(kind_phys) :: p_min_bound(ncol), p_max_bound(ncol) + real(kind_phys) :: z_steer(ncol), z_launch(ncol) + real(kind_phys) :: max_low, max_mid, roo + logical :: mask_low(pver), mask_mid(pver) + logical :: blocked_roo integer :: i + p_min_bound = p_min_centroid + p_max_bound = frac_pbl_excl * pmid(:,pver) + call vorticity_centroid_levels(tilt, pmid, ncol, pver, & - p_min_centroid, p_max_centroid, & - steering_level, launch_level) + p_min_bound, p_max_bound, & + steering_level, launch_level, & + z_steer, z_launch) do i = 1, ncol + mask_low = zm(i,:) <= z_low_roo + mask_mid = zm(i,:) > z_low_roo .and. zm(i,:) <= z_mid_top_roo + max_low = maxval(abs(tilt(i,:)), mask=mask_low) + max_mid = maxval(abs(tilt(i,:)), mask=mask_mid) + roo = max_low / (max_mid + roo_eps) + blocked_roo = roo > roo_thresh + roo_block(i) = merge(1._kind_phys, 0._kind_phys, blocked_roo) + valid(i) = steering_level(i) >= 2 .and. launch_level(i) >= 2 .and. & - steering_level(i) > launch_level(i) + steering_level(i) > launch_level(i) .and. .not. blocked_roo if (.not. valid(i)) then steering_level(i) = 2 launch_level(i) = 1 @@ -808,26 +876,29 @@ end subroutine tilt_centroid_levels !========================================================================== - subroutine tilt_dyn_src(tilt, u, v, pmid, delp, ncol, pver, alpha_gw_movmtn, & + subroutine tilt_dyn_src(tilt, u, v, pmid, delp, zm, ncol, pver, alpha_gw_movmtn, & xpwp_src, steering_level, launch_level, p_steer, p_launch, & - usteer, vsteer) + usteer, vsteer, roo_block) !------------------------------------------------------------------------ ! source_type=3: flow-dependent GW source using the tilt field. ! ! Source flux = pressure-weighted mean tilt between launch and steering ! levels, scaled by alpha_gw_movmtn and scale_factor. Zero where the - ! centroid search failed. + ! centroid search failed or the roo criterion blocked it (see + ! tilt_centroid_levels). !------------------------------------------------------------------------ integer, intent(in) :: ncol, pver real(kind_phys), intent(in) :: tilt(ncol,pver) real(kind_phys), intent(in) :: u(ncol,pver), v(ncol,pver) real(kind_phys), intent(in) :: pmid(ncol,pver) real(kind_phys), intent(in) :: delp(ncol,pver) + real(kind_phys), intent(in) :: zm(ncol,pver) real(kind_phys), intent(in) :: alpha_gw_movmtn real(kind_phys), intent(out) :: xpwp_src(ncol) integer, intent(out) :: steering_level(ncol), launch_level(ncol) real(kind_phys), intent(out) :: p_steer(ncol), p_launch(ncol) real(kind_phys), intent(out) :: usteer(ncol), vsteer(ncol) + real(kind_phys), intent(out) :: roo_block(ncol) ! Dimensional scale factor: converts tilt [s-2] to flux [Pa]. real(kind_phys), parameter :: scale_factor = 1.e6_kind_phys @@ -835,8 +906,8 @@ subroutine tilt_dyn_src(tilt, u, v, pmid, delp, ncol, pver, alpha_gw_movmtn, & logical :: valid(ncol) integer :: i, ks, kl - call tilt_centroid_levels(tilt, u, v, pmid, ncol, pver, & - steering_level, launch_level, p_steer, p_launch, usteer, vsteer, valid) + call tilt_centroid_levels(tilt, u, v, pmid, zm, ncol, pver, & + steering_level, launch_level, p_steer, p_launch, usteer, vsteer, valid, roo_block) xpwp_src = 0._kind_phys do i = 1, ncol @@ -853,9 +924,9 @@ end subroutine tilt_dyn_src !========================================================================== - subroutine tilt_precip_6param_src(tilt, prect, u, v, pmid, ncol, pver, & + subroutine tilt_precip_6param_src(tilt, prect, u, v, pmid, zm, ncol, pver, & alpha_gw_movmtn, xpwp_src, steering_level, launch_level, p_steer, p_launch, & - usteer, vsteer) + usteer, vsteer, roo_block) !------------------------------------------------------------------------ ! source_type=4: six-parameter tilt + precipitation source (Aug 2026). ! @@ -879,11 +950,13 @@ subroutine tilt_precip_6param_src(tilt, prect, u, v, pmid, ncol, pver, & real(kind_phys), intent(in) :: prect(ncol) real(kind_phys), intent(in) :: u(ncol,pver), v(ncol,pver) real(kind_phys), intent(in) :: pmid(ncol,pver) + real(kind_phys), intent(in) :: zm(ncol,pver) real(kind_phys), intent(in) :: alpha_gw_movmtn real(kind_phys), intent(out) :: xpwp_src(ncol) integer, intent(out) :: steering_level(ncol), launch_level(ncol) real(kind_phys), intent(out) :: p_steer(ncol), p_launch(ncol) real(kind_phys), intent(out) :: usteer(ncol), vsteer(ncol) + real(kind_phys), intent(out) :: roo_block(ncol) ! Reference scales. real(kind_phys), parameter :: d_ref = 3.0e-7_kind_phys ! [s-2] @@ -901,8 +974,8 @@ subroutine tilt_precip_6param_src(tilt, prect, u, v, pmid, ncol, pver, & real(kind_phys) :: d_val, dyn_term, conv_term integer :: i - call tilt_centroid_levels(tilt, u, v, pmid, ncol, pver, & - steering_level, launch_level, p_steer, p_launch, usteer, vsteer, valid) + call tilt_centroid_levels(tilt, u, v, pmid, zm, ncol, pver, & + steering_level, launch_level, p_steer, p_launch, usteer, vsteer, valid, roo_block) xpwp_src = 0._kind_phys do i = 1, ncol @@ -929,9 +1002,9 @@ end subroutine tilt_precip_6param_src !========================================================================== - subroutine tilt_precip_cx17_src(tilt, prect, u, v, pmid, delp, ncol, pver, & + subroutine tilt_precip_cx17_src(tilt, prect, u, v, pmid, delp, zm, ncol, pver, & alpha_gw_movmtn, xpwp_src, steering_level, launch_level, p_steer, p_launch, & - usteer, vsteer) + usteer, vsteer, roo_block) !------------------------------------------------------------------------ ! source_type=5: PySR equation-discovery source, complexity 17 ! (front_full_both.json). @@ -957,11 +1030,13 @@ subroutine tilt_precip_cx17_src(tilt, prect, u, v, pmid, delp, ncol, pver, & real(kind_phys), intent(in) :: u(ncol,pver), v(ncol,pver) real(kind_phys), intent(in) :: pmid(ncol,pver) real(kind_phys), intent(in) :: delp(ncol,pver) + real(kind_phys), intent(in) :: zm(ncol,pver) real(kind_phys), intent(in) :: alpha_gw_movmtn real(kind_phys), intent(out) :: xpwp_src(ncol) integer, intent(out) :: steering_level(ncol), launch_level(ncol) real(kind_phys), intent(out) :: p_steer(ncol), p_launch(ncol) real(kind_phys), intent(out) :: usteer(ncol), vsteer(ncol) + real(kind_phys), intent(out) :: roo_block(ncol) ! Reference scales (fixed before the search). real(kind_phys), parameter :: tau_ref = 4.68e-3_kind_phys ! [Pa] @@ -976,8 +1051,8 @@ subroutine tilt_precip_cx17_src(tilt, prect, u, v, pmid, delp, ncol, pver, & real(kind_phys) :: p_s, u_s, d_g, d_l integer :: i, ks, kl - call tilt_centroid_levels(tilt, u, v, pmid, ncol, pver, & - steering_level, launch_level, p_steer, p_launch, usteer, vsteer, valid) + call tilt_centroid_levels(tilt, u, v, pmid, zm, ncol, pver, & + steering_level, launch_level, p_steer, p_launch, usteer, vsteer, valid, roo_block) xpwp_src = 0._kind_phys do i = 1, ncol @@ -1004,7 +1079,8 @@ end subroutine tilt_precip_cx17_src subroutine vorticity_centroid_levels(weight, p, ncol, pver, & p_min, p_max, & - steering_level, launch_level) + steering_level, launch_level, & + z_steer, z_launch) !------------------------------------------------------------------------ ! Derive steering and launch level indices from a weight profile ! (typically |vorticity| or tilt) via two successive centroid passes. @@ -1015,15 +1091,21 @@ subroutine vorticity_centroid_levels(weight, p, ncol, pver, & ! above the steering level (smaller k = higher altitude). ! ! p is the vertical coordinate for the centroid; callers pass pmid [Pa]. + ! + ! p_min/p_max are now per-column (ncol) bounds, so the steering-level + ! search window can depend on local surface pressure rather than a + ! single fixed value (see tilt_centroid_levels). !------------------------------------------------------------------------ integer, intent(in) :: ncol, pver real(kind_phys), intent(in) :: weight(ncol,pver) real(kind_phys), intent(in) :: p(ncol,pver) - real(kind_phys), intent(in) :: p_min, p_max + real(kind_phys), intent(in) :: p_min(ncol), p_max(ncol) integer, intent(out) :: steering_level(ncol), launch_level(ncol) + real(kind_phys), intent(out) :: z_steer(ncol), z_launch(ncol) real(kind_phys) :: abswgt(ncol,pver), w_top(ncol,pver) real(kind_phys) :: z_cent(ncol), z_top(ncol) + real(kind_phys) :: unbound_min(ncol), unbound_max(ncol) integer :: i abswgt = abs(weight) @@ -1042,10 +1124,15 @@ subroutine vorticity_centroid_levels(weight, p, ncol, pver, & end if end do + unbound_min = -huge(1._kind_phys) + unbound_max = huge(1._kind_phys) call weighted_centroid(w_top, p, ncol, pver, & - -huge(1._kind_phys), huge(1._kind_phys), & + unbound_min, unbound_max, & z_top, launch_level) + z_steer = z_cent + z_launch = z_top + end subroutine vorticity_centroid_levels !========================================================================== @@ -1056,12 +1143,16 @@ subroutine weighted_centroid(w, z, ncol, pver, z_min, z_max, z_centroid, k_centr ! ! z_centroid = integral( z*w dz ) / integral( w dz ) ! - ! restricted to z_min <= z <= z_max. Integration is trapezoidal. + ! restricted to z_min(i) <= z <= z_max(i). Integration is trapezoidal. ! Sentinel on failure: z_centroid = -1, k_centroid = -1. + ! + ! z_min/z_max are now per-column (ncol) rather than scalar, so callers + ! can set a column-dependent bound (e.g. a fraction of local surface + ! pressure, instead of a fixed value that is terrain-blind). !------------------------------------------------------------------------ integer, intent(in) :: ncol, pver real(kind_phys), intent(in) :: w(ncol,pver), z(ncol,pver) - real(kind_phys), intent(in) :: z_min, z_max + real(kind_phys), intent(in) :: z_min(ncol), z_max(ncol) real(kind_phys), intent(out) :: z_centroid(ncol) integer, intent(out) :: k_centroid(ncol) @@ -1072,7 +1163,7 @@ subroutine weighted_centroid(w, z, ncol, pver, z_min, z_max, z_centroid, k_centr do i = 1, ncol zi = z(i,:) wi = w(i,:) - where (zi < z_min .or. zi > z_max) wi = 0._kind_phys + where (zi < z_min(i) .or. zi > z_max(i)) wi = 0._kind_phys num = 0._kind_phys denom = 0._kind_phys @@ -1083,6 +1174,10 @@ subroutine weighted_centroid(w, z, ncol, pver, z_min, z_max, z_centroid, k_centr if (denom /= 0._kind_phys) then z_centroid(i) = num / denom + ! Block centroid search from returning levels above 10_000Pa (100hPa) + if ( z_centroid(i) < 10000._kind_phys) then + z_centroid(i) = 10000._kind_phys + end if k_centroid(i) = minloc(abs(zi - z_centroid(i)), dim=1) else z_centroid(i) = -1._kind_phys diff --git a/schemes/gravity_wave_drag/gravity_wave_drag_moving_mountain.meta b/schemes/gravity_wave_drag/gravity_wave_drag_moving_mountain.meta index 79fcb5fd7..e8622ac79 100644 --- a/schemes/gravity_wave_drag/gravity_wave_drag_moving_mountain.meta +++ b/schemes/gravity_wave_drag/gravity_wave_drag_moving_mountain.meta @@ -419,6 +419,12 @@ type = real | kind = kind_phys dimensions = (horizontal_loop_extent) intent = out +[ roo_block ] + standard_name = flag_for_near_surface_tilt_dominance_blocking_of_moving_mountain_gravity_wave_drag_source + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = out [ taucd_west ] standard_name = westward_reynolds_stress_due_to_moving_mountain_gravity_wave_drag units = N m-2 From f696a8f178aefd1333087d1aa734cd55dcc05b6e Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Mon, 9 Nov 2026 12:45:07 -0500 Subject: [PATCH 3/8] Clean up moving-mountain roo blocking and mirror ROO_BLOCK_MOVMTN in SIMA diagnostics No answer changes. Removes the unused z_steer/z_launch outputs of vorticity_centroid_levels, names the 100 hPa centroid floor in weighted_centroid, makes the centroid-bound comments self-contained, and documents the momentum fixer sign fix in terms of what calc_taucd and momentum_flux read. gravity_wave_drag_moving_mountain_diagnostics gains ROO_BLOCK_MOVMTN so CAM-SIMA writes the same history field as CAM. Assisted-by: claude-fable:5.1 --- .../gravity_wave_drag_moving_mountain.F90 | 69 ++++++++----------- ..._wave_drag_moving_mountain_diagnostics.F90 | 6 +- ...wave_drag_moving_mountain_diagnostics.meta | 6 ++ 3 files changed, 41 insertions(+), 40 deletions(-) diff --git a/schemes/gravity_wave_drag/gravity_wave_drag_moving_mountain.F90 b/schemes/gravity_wave_drag/gravity_wave_drag_moving_mountain.F90 index 978574225..970911605 100644 --- a/schemes/gravity_wave_drag/gravity_wave_drag_moving_mountain.F90 +++ b/schemes/gravity_wave_drag/gravity_wave_drag_moving_mountain.F90 @@ -357,10 +357,15 @@ subroutine gravity_wave_drag_moving_mountain_run(ncol, pver, pcnst, & ttgw, qtgw, egwdffi, gwut, dttdf, dttke, & lapply_effgw_in=gw_apply_tndmax) - ! Adjust ubi to yield correct directions when ubi(src)>0>ubi(src+1) + ! calc_taucd signs the stress at each interface with (c - ubi). + ! At the source interface src_level+1 the interpolated ubi can have the opposite sign to + ! ubm(src_level), the midpoint wind that sets the drag direction gw_drag_prof applies in + ! the source layer. + ! momentum_flux reads taucd at that interface, so use ubm(src_level) there to keep the + ! momentum fixer consistent with the applied drag. ubi itself is left unchanged. ubi_src(:ncol,:) = ubi(:ncol,:) - do i=1,ncol - ubi_src(i,src_level(i)+1) = ubm(i ,src_level(i)) + do i = 1, ncol + ubi_src(i, src_level(i)+1) = ubm(i, src_level(i)) end do ! Project stress into directional components. taucd = calc_taucd(ncol, band%ngwv, tend_level, tau, phase_speeds, xv, yv, ubi_src) @@ -806,24 +811,16 @@ subroutine tilt_centroid_levels(tilt, u, v, pmid, zm, ncol, pver, & real(kind_phys), intent(out) :: roo_block(ncol) ! Pressure bounds for centroid search [Pa]. - ! - ! Top bound: fixed at TOA (unchanged). - ! - ! Bottom (BL-exclusion) bound: was a single fixed 800 hPa, which is - ! terrain-blind. Over elevated terrain (e.g. the Antarctic plateau, - ! Andes, Himalaya/Tibetan Plateau, Rockies -- anywhere local surface - ! pressure is already below ~800 hPa) that fixed cutoff is never - ! reached, so no level is ever excluded and the steering-level - ! centroid is free to lock onto raw near-surface shear/tilt noise - ! (which is exactly where this parameterization's local eta-level - ! vorticity calc and polar grid-metric terms are most contaminated). - ! - ! Fix: exclude the bottom (1 - frac_pbl_excl) of each column's own - ! pressure range instead, i.e. scale the cutoff by local surface - ! pressure (pmid(:,pver), the lowest model layer midpoint, used as - ! a close proxy for true surface pressure -- using the true - ! interface/surface pressure instead would require threading it in - ! from the caller, which is not done here to keep the change local). + ! Top bound: TOA. + ! Bottom (BL-exclusion) bound: the lowest (1 - frac_pbl_excl) of each column's own + ! pressure range is excluded, so the cutoff scales with local surface pressure. + ! A fixed cutoff such as 800 hPa is terrain-blind: over elevated terrain (Antarctic + ! plateau, Andes, Tibetan Plateau, Rockies) local surface pressure is already below it, + ! no level is excluded, and the steering-level centroid can lock onto raw near-surface + ! shear/tilt noise, exactly where the local eta-level vorticity and polar grid-metric + ! terms are most contaminated. + ! pmid(:,pver), the lowest layer midpoint, stands in for surface pressure so that the + ! surface pressure does not have to be threaded in from the caller. real(kind_phys), parameter :: p_min_centroid = 0._kind_phys real(kind_phys), parameter :: frac_pbl_excl = 0.8_kind_phys ! fraction of local p_sfc @@ -837,7 +834,6 @@ subroutine tilt_centroid_levels(tilt, u, v, pmid, zm, ncol, pver, & real(kind_phys), parameter :: roo_eps = 1.e-7_kind_phys ! ~typical tilt scale, avoids /0 real(kind_phys) :: p_min_bound(ncol), p_max_bound(ncol) - real(kind_phys) :: z_steer(ncol), z_launch(ncol) real(kind_phys) :: max_low, max_mid, roo logical :: mask_low(pver), mask_mid(pver) logical :: blocked_roo @@ -848,8 +844,7 @@ subroutine tilt_centroid_levels(tilt, u, v, pmid, zm, ncol, pver, & call vorticity_centroid_levels(tilt, pmid, ncol, pver, & p_min_bound, p_max_bound, & - steering_level, launch_level, & - z_steer, z_launch) + steering_level, launch_level) do i = 1, ncol mask_low = zm(i,:) <= z_low_roo @@ -1079,8 +1074,7 @@ end subroutine tilt_precip_cx17_src subroutine vorticity_centroid_levels(weight, p, ncol, pver, & p_min, p_max, & - steering_level, launch_level, & - z_steer, z_launch) + steering_level, launch_level) !------------------------------------------------------------------------ ! Derive steering and launch level indices from a weight profile ! (typically |vorticity| or tilt) via two successive centroid passes. @@ -1092,16 +1086,14 @@ subroutine vorticity_centroid_levels(weight, p, ncol, pver, & ! ! p is the vertical coordinate for the centroid; callers pass pmid [Pa]. ! - ! p_min/p_max are now per-column (ncol) bounds, so the steering-level - ! search window can depend on local surface pressure rather than a - ! single fixed value (see tilt_centroid_levels). + ! p_min/p_max are per-column bounds so the steering-level search window can depend + ! on local surface pressure (see tilt_centroid_levels). !------------------------------------------------------------------------ integer, intent(in) :: ncol, pver real(kind_phys), intent(in) :: weight(ncol,pver) real(kind_phys), intent(in) :: p(ncol,pver) real(kind_phys), intent(in) :: p_min(ncol), p_max(ncol) integer, intent(out) :: steering_level(ncol), launch_level(ncol) - real(kind_phys), intent(out) :: z_steer(ncol), z_launch(ncol) real(kind_phys) :: abswgt(ncol,pver), w_top(ncol,pver) real(kind_phys) :: z_cent(ncol), z_top(ncol) @@ -1130,9 +1122,6 @@ subroutine vorticity_centroid_levels(weight, p, ncol, pver, & unbound_min, unbound_max, & z_top, launch_level) - z_steer = z_cent - z_launch = z_top - end subroutine vorticity_centroid_levels !========================================================================== @@ -1146,9 +1135,8 @@ subroutine weighted_centroid(w, z, ncol, pver, z_min, z_max, z_centroid, k_centr ! restricted to z_min(i) <= z <= z_max(i). Integration is trapezoidal. ! Sentinel on failure: z_centroid = -1, k_centroid = -1. ! - ! z_min/z_max are now per-column (ncol) rather than scalar, so callers - ! can set a column-dependent bound (e.g. a fraction of local surface - ! pressure, instead of a fixed value that is terrain-blind). + ! z_min/z_max are per-column so callers can set a column-dependent bound + ! (e.g. a fraction of local surface pressure). !------------------------------------------------------------------------ integer, intent(in) :: ncol, pver real(kind_phys), intent(in) :: w(ncol,pver), z(ncol,pver) @@ -1156,6 +1144,10 @@ subroutine weighted_centroid(w, z, ncol, pver, z_min, z_max, z_centroid, k_centr real(kind_phys), intent(out) :: z_centroid(ncol) integer, intent(out) :: k_centroid(ncol) + ! Floor on the centroid [Pa]: callers pass pressure as z, and neither the steering nor + ! the launch level may be placed above 100 hPa. + real(kind_phys), parameter :: z_centroid_floor = 10000._kind_phys + real(kind_phys) :: wi(pver), zi(pver) real(kind_phys) :: num, denom integer :: i, k @@ -1174,9 +1166,8 @@ subroutine weighted_centroid(w, z, ncol, pver, z_min, z_max, z_centroid, k_centr if (denom /= 0._kind_phys) then z_centroid(i) = num / denom - ! Block centroid search from returning levels above 10_000Pa (100hPa) - if ( z_centroid(i) < 10000._kind_phys) then - z_centroid(i) = 10000._kind_phys + if (z_centroid(i) < z_centroid_floor) then + z_centroid(i) = z_centroid_floor end if k_centroid(i) = minloc(abs(zi - z_centroid(i)), dim=1) else diff --git a/schemes/sima_diagnostics/gravity_wave_drag_moving_mountain_diagnostics.F90 b/schemes/sima_diagnostics/gravity_wave_drag_moving_mountain_diagnostics.F90 index 4f88624a9..1ecf51615 100644 --- a/schemes/sima_diagnostics/gravity_wave_drag_moving_mountain_diagnostics.F90 +++ b/schemes/sima_diagnostics/gravity_wave_drag_moving_mountain_diagnostics.F90 @@ -51,6 +51,8 @@ subroutine gravity_wave_drag_moving_mountain_diagnostics_init(errmsg, errflg) call history_add_field('XPWP_SRC_MOVMTN', 'momentum_flux_source_for_moving_mountain_gravity_wave_drag', horiz_only, 'inst', 'Pa') call history_add_field('PSTEER_MOVMTN', 'air_pressure_at_steering_level_due_to_moving_mountain_gravity_wave_drag', horiz_only, 'inst', 'Pa') call history_add_field('PLAUNCH_MOVMTN', 'air_pressure_at_launch_level_due_to_moving_mountain_gravity_wave_drag', horiz_only, 'inst', 'Pa') + call history_add_field('ROO_BLOCK_MOVMTN', & + 'flag_for_near_surface_tilt_dominance_blocking_of_moving_mountain_gravity_wave_drag_source', horiz_only, 'inst', '1') call history_add_field('MMTAUE', 'eastward_reynolds_stress_due_to_moving_mountain_gravity_wave_drag', 'ilev', 'avg', 'Pa') call history_add_field('MMTAUW', 'westward_reynolds_stress_due_to_moving_mountain_gravity_wave_drag', 'ilev', 'avg', 'Pa') @@ -69,7 +71,7 @@ subroutine gravity_wave_drag_moving_mountain_diagnostics_run( & hdepth, ttend_dp, & ttend_clubb, & upwp_clubb, vpwp_clubb, vorticity, tilt, & - usteer, vsteer, CS, xpwp_src, p_steer, p_launch, & + usteer, vsteer, CS, xpwp_src, p_steer, p_launch, roo_block, & taucd_west, taucd_east, taucd_south, taucd_north, & errmsg, errflg) @@ -105,6 +107,7 @@ subroutine gravity_wave_drag_moving_mountain_diagnostics_run( & real(kind_phys), intent(in) :: xpwp_src(:) ! [Pa] real(kind_phys), intent(in) :: p_steer(:) ! [Pa] real(kind_phys), intent(in) :: p_launch(:) ! [Pa] + real(kind_phys), intent(in) :: roo_block(:) ! 1 = source blocked by roo criterion [1] ! Reynolds stress for waves propagating in each cardinal direction [N m-2] real(kind_phys), intent(in) :: taucd_west(:,:) @@ -142,6 +145,7 @@ subroutine gravity_wave_drag_moving_mountain_diagnostics_run( & call history_out_field('XPWP_SRC_MOVMTN', xpwp_src) call history_out_field('PSTEER_MOVMTN', p_steer) call history_out_field('PLAUNCH_MOVMTN', p_launch) + call history_out_field('ROO_BLOCK_MOVMTN', roo_block) call history_out_field('MMTAUE', taucd_east) call history_out_field('MMTAUW', taucd_west) call history_out_field('MMTAUS', taucd_south) diff --git a/schemes/sima_diagnostics/gravity_wave_drag_moving_mountain_diagnostics.meta b/schemes/sima_diagnostics/gravity_wave_drag_moving_mountain_diagnostics.meta index dda07e574..8e8a136d0 100644 --- a/schemes/sima_diagnostics/gravity_wave_drag_moving_mountain_diagnostics.meta +++ b/schemes/sima_diagnostics/gravity_wave_drag_moving_mountain_diagnostics.meta @@ -165,6 +165,12 @@ type = real | kind = kind_phys dimensions = (horizontal_loop_extent) intent = in +[ roo_block ] + standard_name = flag_for_near_surface_tilt_dominance_blocking_of_moving_mountain_gravity_wave_drag_source + units = 1 + type = real | kind = kind_phys + dimensions = (horizontal_loop_extent) + intent = in [ taucd_west ] standard_name = westward_reynolds_stress_due_to_moving_mountain_gravity_wave_drag units = N m-2 From 58f46eaf41438d3fcc29a5f399bab21b48652143 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Mon, 9 Nov 2026 15:35:00 -0500 Subject: [PATCH 4/8] Comment trim --- .../gravity_wave_drag_moving_mountain.F90 | 109 +++++------------- 1 file changed, 32 insertions(+), 77 deletions(-) diff --git a/schemes/gravity_wave_drag/gravity_wave_drag_moving_mountain.F90 b/schemes/gravity_wave_drag/gravity_wave_drag_moving_mountain.F90 index 970911605..72d5a842f 100644 --- a/schemes/gravity_wave_drag/gravity_wave_drag_moving_mountain.F90 +++ b/schemes/gravity_wave_drag/gravity_wave_drag_moving_mountain.F90 @@ -285,7 +285,7 @@ subroutine gravity_wave_drag_moving_mountain_run(ncol, pver, pcnst, & real(kind_phys), intent(out) :: tilt(:, :) ! Vortex tilting magnitude |zeta|*|dV/dz| [s-2] real(kind_phys), intent(out) :: p_steer(:) ! Pressure at steering level [Pa] real(kind_phys), intent(out) :: p_launch(:) ! Pressure at launch level [Pa] - real(kind_phys), intent(out) :: roo_block(:) ! 1=column blocked by near-surface tilt-dominance (roo) criterion, 0=not blocked + real(kind_phys), intent(out) :: roo_block(:) ! 1 = source blocked by near-surface tilt dominance; 0 otherwise real(kind_phys), intent(out) :: taucd_west(:, :) ! Reynolds stress for waves in W direction, interfaces [N m-2] real(kind_phys), intent(out) :: taucd_east(:, :) ! Reynolds stress for waves in E direction, interfaces [N m-2] real(kind_phys), intent(out) :: taucd_south(:, :) ! Reynolds stress for waves in S direction, interfaces [N m-2] @@ -357,12 +357,9 @@ subroutine gravity_wave_drag_moving_mountain_run(ncol, pver, pcnst, & ttgw, qtgw, egwdffi, gwut, dttdf, dttke, & lapply_effgw_in=gw_apply_tndmax) - ! calc_taucd signs the stress at each interface with (c - ubi). - ! At the source interface src_level+1 the interpolated ubi can have the opposite sign to - ! ubm(src_level), the midpoint wind that sets the drag direction gw_drag_prof applies in - ! the source layer. - ! momentum_flux reads taucd at that interface, so use ubm(src_level) there to keep the - ! momentum fixer consistent with the applied drag. ubi itself is left unchanged. + ! Use the source midpoint wind to sign calc_taucd stress at src_level+1: + ! interpolated ubi can have the opposite sign, making the momentum fixer + ! inconsistent with the drag applied by gw_drag_prof. ubi_src(:ncol,:) = ubi(:ncol,:) do i = 1, ncol ubi_src(i, src_level(i)+1) = ubm(i, src_level(i)) @@ -479,10 +476,8 @@ subroutine gw_movmtn_src(ncol, pver, & real(kind_phys), intent(out) :: tilt(:,:) ! Pressure at steering/launch level [Pa]. real(kind_phys), intent(out) :: p_steer(:), p_launch(:) - ! Flag: 1 if this column's source was zeroed because the near-surface - ! tilt-dominance (roo) criterion blocked the centroid-based search (see - ! tilt_centroid_levels); 0 otherwise. Only ever 1 for source_type=3/4/5; - ! source_type=1/2 use fixed levels and never run that search. + ! 1 = source blocked by near-surface tilt dominance; 0 otherwise. + ! Applies to source_type=3-5 only (see tilt_centroid_levels). real(kind_phys), intent(out) :: roo_block(:) character(len=512), intent(out) :: errmsg integer, intent(out) :: errflg @@ -771,67 +766,42 @@ end subroutine gw_movmtn_src subroutine tilt_centroid_levels(tilt, u, v, pmid, zm, ncol, pver, & steering_level, launch_level, p_steer, p_launch, usteer, vsteer, valid, roo_block) !------------------------------------------------------------------------ - ! Flow-dependent steering and launch levels shared by source_type 3-5. + ! Flow-dependent steering/launch levels from |tilt| pressure centroids + ! (for source_type=3-5). ! - ! 1. Centroid of |tilt| in pressure space gives steering and launch - ! levels (see vorticity_centroid_levels). - ! 2. Sentinel handling: if the centroid search fails (returns -1), - ! valid=.false. and levels fall back to safe values so downstream - ! code is never handed an out-of-bounds index. Callers must zero - ! the source where valid=.false. - ! 3. roo criterion: the centroid search can "succeed" (return a - ! well-defined in-bounds pair of levels) and still be meaningless, - ! when a column's tilt profile is dominated by an isolated, extreme - ! near-surface spike rather than a genuine free-tropospheric signal - ! -- the regime observed over Antarctica, the Arctic Ocean, and - ! Greenland. roo = max(|tilt|, 0-500 m AGL) / max(|tilt|, 500-15000 m - ! AGL); roo > roo_thresh means the near-surface layer so overwhelms - ! the layer above it that the steering/launch search cannot be - ! trusted, so the column is blocked (valid=.false., same safe - ! fallback levels as a failed centroid search) regardless of whether - ! the search itself returned in-bounds levels. roo_block records - ! which columns were blocked for this reason specifically (as - ! opposed to a plain centroid-search failure). + ! roo = max(|tilt|, 0-500 m AGL) / (max(|tilt|, 500-15000 m AGL) + roo_eps). + ! roo > roo_thresh blocks sources dominated by near-surface tilt spikes. + ! roo_block flags this rejection only, not centroid-search failures. + ! Invalid columns get safe fallback indices. ! - ! CAM level ordering: k=1 TOA, k=pver surface, pressure increases - ! with k. Therefore launch_level < steering_level in index space - ! (launch is higher altitude = lower pressure = smaller k). zm is - ! height AGL (m), decreasing with k, consistent with its use in - ! compute_tilt. + ! If the centroid search fails (returns -1), valid=.false. and invalid + ! columns fall back to safe values. Callers must zero their source + ! where valid=.false. !------------------------------------------------------------------------ integer, intent(in) :: ncol, pver real(kind_phys), intent(in) :: tilt(ncol,pver) real(kind_phys), intent(in) :: u(ncol,pver), v(ncol,pver) real(kind_phys), intent(in) :: pmid(ncol,pver) - real(kind_phys), intent(in) :: zm(ncol,pver) - integer, intent(out) :: steering_level(ncol), launch_level(ncol) + real(kind_phys), intent(in) :: zm(ncol,pver) ! [m] + integer, intent(out) :: steering_level(ncol), launch_level(ncol) ! valid launch_level < steering_level. real(kind_phys), intent(out) :: p_steer(ncol), p_launch(ncol) real(kind_phys), intent(out) :: usteer(ncol), vsteer(ncol) logical, intent(out) :: valid(ncol) real(kind_phys), intent(out) :: roo_block(ncol) - ! Pressure bounds for centroid search [Pa]. - ! Top bound: TOA. - ! Bottom (BL-exclusion) bound: the lowest (1 - frac_pbl_excl) of each column's own - ! pressure range is excluded, so the cutoff scales with local surface pressure. - ! A fixed cutoff such as 800 hPa is terrain-blind: over elevated terrain (Antarctic - ! plateau, Andes, Tibetan Plateau, Rockies) local surface pressure is already below it, - ! no level is excluded, and the steering-level centroid can lock onto raw near-surface - ! shear/tilt noise, exactly where the local eta-level vorticity and polar grid-metric - ! terms are most contaminated. - ! pmid(:,pver), the lowest layer midpoint, stands in for surface pressure so that the - ! surface pressure does not have to be threaded in from the caller. + ! Exclude near-surface tilt from the centroid search + ! using a fraction of local surface pressure, + ! approximated by pmid(:,pver). + ! Unlike a fixed pressure cutoff, this also excludes the + ! near-surface layer over elevated terrain. real(kind_phys), parameter :: p_min_centroid = 0._kind_phys real(kind_phys), parameter :: frac_pbl_excl = 0.8_kind_phys ! fraction of local p_sfc - ! roo criterion bounds: a fixed, terrain-following (AGL, not - ! pressure) comparison, since the regimes where this fires are - ! precisely where a fixed-pressure near-surface layer would span - ! wildly different physical depths depending on surface elevation. + ! AGL bounds keep the roo comparison depths consistent across terrain elevations. real(kind_phys), parameter :: z_low_roo = 500._kind_phys ! m AGL real(kind_phys), parameter :: z_mid_top_roo = 15000._kind_phys ! m AGL real(kind_phys), parameter :: roo_thresh = 5._kind_phys ! try 10 for a stricter trigger - real(kind_phys), parameter :: roo_eps = 1.e-7_kind_phys ! ~typical tilt scale, avoids /0 + real(kind_phys), parameter :: roo_eps = 1.e-7_kind_phys ! Tilt scale [s-2]; avoids division by zero real(kind_phys) :: p_min_bound(ncol), p_max_bound(ncol) real(kind_phys) :: max_low, max_mid, roo @@ -875,12 +845,9 @@ subroutine tilt_dyn_src(tilt, u, v, pmid, delp, zm, ncol, pver, alpha_gw_movmtn, xpwp_src, steering_level, launch_level, p_steer, p_launch, & usteer, vsteer, roo_block) !------------------------------------------------------------------------ - ! source_type=3: flow-dependent GW source using the tilt field. - ! - ! Source flux = pressure-weighted mean tilt between launch and steering - ! levels, scaled by alpha_gw_movmtn and scale_factor. Zero where the - ! centroid search failed or the roo criterion blocked it (see - ! tilt_centroid_levels). + ! source_type=3: pressure-weighted mean tilt between launch and steering + ! levels, scaled by alpha_gw_movmtn and scale_factor. Invalid columns from + ! tilt_centroid_levels have zero source flux. !------------------------------------------------------------------------ integer, intent(in) :: ncol, pver real(kind_phys), intent(in) :: tilt(ncol,pver) @@ -1076,18 +1043,9 @@ subroutine vorticity_centroid_levels(weight, p, ncol, pver, & p_min, p_max, & steering_level, launch_level) !------------------------------------------------------------------------ - ! Derive steering and launch level indices from a weight profile - ! (typically |vorticity| or tilt) via two successive centroid passes. - ! - ! Pass 1 — steering level: centroid of weight over the full column - ! within [p_min, p_max] (pressure bounds, Pa). - ! Pass 2 — launch level: centroid of weight restricted to levels - ! above the steering level (smaller k = higher altitude). - ! - ! p is the vertical coordinate for the centroid; callers pass pmid [Pa]. - ! - ! p_min/p_max are per-column bounds so the steering-level search window can depend - ! on local surface pressure (see tilt_centroid_levels). + ! Steering level: centroid of |weight| within per-column pressure bounds + ! [p_min, p_max]. Launch level: centroid above the steering level (smaller k). + ! Callers pass pmid [Pa] as the pressure coordinate p. !------------------------------------------------------------------------ integer, intent(in) :: ncol, pver real(kind_phys), intent(in) :: weight(ncol,pver) @@ -1134,9 +1092,6 @@ subroutine weighted_centroid(w, z, ncol, pver, z_min, z_max, z_centroid, k_centr ! ! restricted to z_min(i) <= z <= z_max(i). Integration is trapezoidal. ! Sentinel on failure: z_centroid = -1, k_centroid = -1. - ! - ! z_min/z_max are per-column so callers can set a column-dependent bound - ! (e.g. a fraction of local surface pressure). !------------------------------------------------------------------------ integer, intent(in) :: ncol, pver real(kind_phys), intent(in) :: w(ncol,pver), z(ncol,pver) @@ -1144,8 +1099,8 @@ subroutine weighted_centroid(w, z, ncol, pver, z_min, z_max, z_centroid, k_centr real(kind_phys), intent(out) :: z_centroid(ncol) integer, intent(out) :: k_centroid(ncol) - ! Floor on the centroid [Pa]: callers pass pressure as z, and neither the steering nor - ! the launch level may be placed above 100 hPa. + ! Callers pass pressure as z. + ! limit steering/launch centroids to >= 100 hPa (centroid floor) real(kind_phys), parameter :: z_centroid_floor = 10000._kind_phys real(kind_phys) :: wi(pver), zi(pver) From c38a742b786c386308359bb31ca476092dd569ec Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Mon, 9 Nov 2026 20:54:11 -0500 Subject: [PATCH 5/8] Address review comments --- .../gravity_wave_drag_moving_mountain.F90 | 68 +++++++++++-------- .../gravity_wave_drag_moving_mountain.meta | 4 +- ..._wave_drag_moving_mountain_diagnostics.F90 | 12 +++- ...wave_drag_moving_mountain_diagnostics.meta | 10 ++- 4 files changed, 59 insertions(+), 35 deletions(-) diff --git a/schemes/gravity_wave_drag/gravity_wave_drag_moving_mountain.F90 b/schemes/gravity_wave_drag/gravity_wave_drag_moving_mountain.F90 index 72d5a842f..af5b240b6 100644 --- a/schemes/gravity_wave_drag/gravity_wave_drag_moving_mountain.F90 +++ b/schemes/gravity_wave_drag/gravity_wave_drag_moving_mountain.F90 @@ -285,7 +285,7 @@ subroutine gravity_wave_drag_moving_mountain_run(ncol, pver, pcnst, & real(kind_phys), intent(out) :: tilt(:, :) ! Vortex tilting magnitude |zeta|*|dV/dz| [s-2] real(kind_phys), intent(out) :: p_steer(:) ! Pressure at steering level [Pa] real(kind_phys), intent(out) :: p_launch(:) ! Pressure at launch level [Pa] - real(kind_phys), intent(out) :: roo_block(:) ! 1 = source blocked by near-surface tilt dominance; 0 otherwise + logical, intent(out) :: roo_block(:) ! Source zeroed by the near-surface tilt-dominance (roo) criterion real(kind_phys), intent(out) :: taucd_west(:, :) ! Reynolds stress for waves in W direction, interfaces [N m-2] real(kind_phys), intent(out) :: taucd_east(:, :) ! Reynolds stress for waves in E direction, interfaces [N m-2] real(kind_phys), intent(out) :: taucd_south(:, :) ! Reynolds stress for waves in S direction, interfaces [N m-2] @@ -357,9 +357,14 @@ subroutine gravity_wave_drag_moving_mountain_run(ncol, pver, pcnst, & ttgw, qtgw, egwdffi, gwut, dttdf, dttke, & lapply_effgw_in=gw_apply_tndmax) - ! Use the source midpoint wind to sign calc_taucd stress at src_level+1: - ! interpolated ubi can have the opposite sign, making the momentum fixer - ! inconsistent with the drag applied by gw_drag_prof. + ! calc_taucd assigns stress to a cardinal direction from the sign of + ! c - ubi at each interface. + ! At the source interface src_level+1 the interpolated ubi can have the + ! opposite sign to ubm(src_level), the wind gw_drag_prof used to apply + ! the source-layer drag, so the momentum fixer would double the column + ! momentum change instead of cancelling it. + ! Pass calc_taucd a copy of ubi with that one element replaced by + ! ubm(src_level). Note: gw_drag_prof itself still uses the unmodified ubi. ubi_src(:ncol,:) = ubi(:ncol,:) do i = 1, ncol ubi_src(i, src_level(i)+1) = ubm(i, src_level(i)) @@ -476,9 +481,10 @@ subroutine gw_movmtn_src(ncol, pver, & real(kind_phys), intent(out) :: tilt(:,:) ! Pressure at steering/launch level [Pa]. real(kind_phys), intent(out) :: p_steer(:), p_launch(:) - ! 1 = source blocked by near-surface tilt dominance; 0 otherwise. - ! Applies to source_type=3-5 only (see tilt_centroid_levels). - real(kind_phys), intent(out) :: roo_block(:) + ! True where the source was zeroed by the near-surface tilt-dominance + ! (roo) criterion in tilt_centroid_levels. + ! Always false for source_type = {1,2} since it does not run that search + logical, intent(out) :: roo_block(:) character(len=512), intent(out) :: errmsg integer, intent(out) :: errflg @@ -530,7 +536,7 @@ subroutine gw_movmtn_src(ncol, pver, & q0 = 0._kind_phys p_steer = 0._kind_phys p_launch = 0._kind_phys - roo_block = 0._kind_phys + roo_block = .false. !---------------------------------------------------------------------- ! source_type>=3 deposits xpwp_src directly at the launch level, which @@ -769,14 +775,21 @@ subroutine tilt_centroid_levels(tilt, u, v, pmid, zm, ncol, pver, & ! Flow-dependent steering/launch levels from |tilt| pressure centroids ! (for source_type=3-5). ! - ! roo = max(|tilt|, 0-500 m AGL) / (max(|tilt|, 500-15000 m AGL) + roo_eps). - ! roo > roo_thresh blocks sources dominated by near-surface tilt spikes. - ! roo_block flags this rejection only, not centroid-search failures. - ! Invalid columns get safe fallback indices. + ! The centroid search can return in-bounds levels that are still meaningless + ! when the column's tilt profile is dominated by an isolated near-surface + ! spike rather than a free-tropospheric signal, as seen over Antarctica, + ! the Arctic Ocean and Greenland. + ! The "roo" criterion rejects such columns using the ratio of the + ! near-surface to the free-tropospheric maximum |tilt|, + ! roo = max(|tilt|, 0-500 m AGL) / (max(|tilt|, 500-15000 m AGL) + roo_eps), + ! and blocks the source where roo > roo_thresh. + ! This is an empirical criterion from moving-mountain source development; + ! it has no published reference. + ! roo_block is true only for this rejection, not for centroid-search failures. ! - ! If the centroid search fails (returns -1), valid=.false. and invalid - ! columns fall back to safe values. Callers must zero their source - ! where valid=.false. + ! If the centroid search fails (returns -1) or roo blocks the column, + ! valid=.false. and the levels fall back to safe in-bounds values. + ! Callers must zero their source where valid=.false. !------------------------------------------------------------------------ integer, intent(in) :: ncol, pver real(kind_phys), intent(in) :: tilt(ncol,pver) @@ -787,13 +800,10 @@ subroutine tilt_centroid_levels(tilt, u, v, pmid, zm, ncol, pver, & real(kind_phys), intent(out) :: p_steer(ncol), p_launch(ncol) real(kind_phys), intent(out) :: usteer(ncol), vsteer(ncol) logical, intent(out) :: valid(ncol) - real(kind_phys), intent(out) :: roo_block(ncol) + logical, intent(out) :: roo_block(ncol) - ! Exclude near-surface tilt from the centroid search - ! using a fraction of local surface pressure, - ! approximated by pmid(:,pver). - ! Unlike a fixed pressure cutoff, this also excludes the - ! near-surface layer over elevated terrain. + ! Pressure bounds for the centroid search [Pa]: TOA down to a fraction + ! of the local surface pressure (see p_max_bound below). real(kind_phys), parameter :: p_min_centroid = 0._kind_phys real(kind_phys), parameter :: frac_pbl_excl = 0.8_kind_phys ! fraction of local p_sfc @@ -806,9 +816,12 @@ subroutine tilt_centroid_levels(tilt, u, v, pmid, zm, ncol, pver, & real(kind_phys) :: p_min_bound(ncol), p_max_bound(ncol) real(kind_phys) :: max_low, max_mid, roo logical :: mask_low(pver), mask_mid(pver) - logical :: blocked_roo integer :: i + ! Exclude near-surface tilt from the centroid search by capping it at + ! a fraction of the local surface pressure, approximated by pmid(:,pver). + ! Unlike a fixed pressure cutoff, this also excludes the near-surface + ! layer over elevated terrain. p_min_bound = p_min_centroid p_max_bound = frac_pbl_excl * pmid(:,pver) @@ -822,11 +835,10 @@ subroutine tilt_centroid_levels(tilt, u, v, pmid, zm, ncol, pver, & max_low = maxval(abs(tilt(i,:)), mask=mask_low) max_mid = maxval(abs(tilt(i,:)), mask=mask_mid) roo = max_low / (max_mid + roo_eps) - blocked_roo = roo > roo_thresh - roo_block(i) = merge(1._kind_phys, 0._kind_phys, blocked_roo) + roo_block(i) = roo > roo_thresh valid(i) = steering_level(i) >= 2 .and. launch_level(i) >= 2 .and. & - steering_level(i) > launch_level(i) .and. .not. blocked_roo + steering_level(i) > launch_level(i) .and. .not. roo_block(i) if (.not. valid(i)) then steering_level(i) = 2 launch_level(i) = 1 @@ -860,7 +872,7 @@ subroutine tilt_dyn_src(tilt, u, v, pmid, delp, zm, ncol, pver, alpha_gw_movmtn, integer, intent(out) :: steering_level(ncol), launch_level(ncol) real(kind_phys), intent(out) :: p_steer(ncol), p_launch(ncol) real(kind_phys), intent(out) :: usteer(ncol), vsteer(ncol) - real(kind_phys), intent(out) :: roo_block(ncol) + logical, intent(out) :: roo_block(ncol) ! Dimensional scale factor: converts tilt [s-2] to flux [Pa]. real(kind_phys), parameter :: scale_factor = 1.e6_kind_phys @@ -918,7 +930,7 @@ subroutine tilt_precip_6param_src(tilt, prect, u, v, pmid, zm, ncol, pver, & integer, intent(out) :: steering_level(ncol), launch_level(ncol) real(kind_phys), intent(out) :: p_steer(ncol), p_launch(ncol) real(kind_phys), intent(out) :: usteer(ncol), vsteer(ncol) - real(kind_phys), intent(out) :: roo_block(ncol) + logical, intent(out) :: roo_block(ncol) ! Reference scales. real(kind_phys), parameter :: d_ref = 3.0e-7_kind_phys ! [s-2] @@ -998,7 +1010,7 @@ subroutine tilt_precip_cx17_src(tilt, prect, u, v, pmid, delp, zm, ncol, pver, & integer, intent(out) :: steering_level(ncol), launch_level(ncol) real(kind_phys), intent(out) :: p_steer(ncol), p_launch(ncol) real(kind_phys), intent(out) :: usteer(ncol), vsteer(ncol) - real(kind_phys), intent(out) :: roo_block(ncol) + logical, intent(out) :: roo_block(ncol) ! Reference scales (fixed before the search). real(kind_phys), parameter :: tau_ref = 4.68e-3_kind_phys ! [Pa] diff --git a/schemes/gravity_wave_drag/gravity_wave_drag_moving_mountain.meta b/schemes/gravity_wave_drag/gravity_wave_drag_moving_mountain.meta index e8622ac79..b9f014a3d 100644 --- a/schemes/gravity_wave_drag/gravity_wave_drag_moving_mountain.meta +++ b/schemes/gravity_wave_drag/gravity_wave_drag_moving_mountain.meta @@ -421,8 +421,8 @@ intent = out [ roo_block ] standard_name = flag_for_near_surface_tilt_dominance_blocking_of_moving_mountain_gravity_wave_drag_source - units = 1 - type = real | kind = kind_phys + units = flag + type = logical dimensions = (horizontal_loop_extent) intent = out [ taucd_west ] diff --git a/schemes/sima_diagnostics/gravity_wave_drag_moving_mountain_diagnostics.F90 b/schemes/sima_diagnostics/gravity_wave_drag_moving_mountain_diagnostics.F90 index 1ecf51615..c5adc3106 100644 --- a/schemes/sima_diagnostics/gravity_wave_drag_moving_mountain_diagnostics.F90 +++ b/schemes/sima_diagnostics/gravity_wave_drag_moving_mountain_diagnostics.F90 @@ -64,7 +64,7 @@ end subroutine gravity_wave_drag_moving_mountain_diagnostics_init !> \section arg_table_gravity_wave_drag_moving_mountain_diagnostics_run Argument Table !! \htmlinclude gravity_wave_drag_moving_mountain_diagnostics_run.html subroutine gravity_wave_drag_moving_mountain_diagnostics_run( & - u, v, & + ncol, u, v, & src_level, tend_level, steer_level, & ubi, ubm, & tau0, gwut0, utgw, vtgw, & @@ -79,6 +79,8 @@ subroutine gravity_wave_drag_moving_mountain_diagnostics_run( & use cam_history, only: history_out_field + integer, intent(in) :: ncol + ! Input winds [m s-1] real(kind_phys), intent(in) :: u(:,:) real(kind_phys), intent(in) :: v(:,:) @@ -107,7 +109,7 @@ subroutine gravity_wave_drag_moving_mountain_diagnostics_run( & real(kind_phys), intent(in) :: xpwp_src(:) ! [Pa] real(kind_phys), intent(in) :: p_steer(:) ! [Pa] real(kind_phys), intent(in) :: p_launch(:) ! [Pa] - real(kind_phys), intent(in) :: roo_block(:) ! 1 = source blocked by roo criterion [1] + logical, intent(in) :: roo_block(:) ! source zeroed by the near-surface tilt-dominance (roo) criterion ! Reynolds stress for waves propagating in each cardinal direction [N m-2] real(kind_phys), intent(in) :: taucd_west(:,:) @@ -118,6 +120,9 @@ subroutine gravity_wave_drag_moving_mountain_diagnostics_run( & character(len=512), intent(out) :: errmsg integer, intent(out) :: errflg + ! history_out_field takes real fields, so roo_block is written as 1 (blocked) or 0. + real(kind_phys) :: roo_block_real(ncol) + errmsg = '' errflg = 0 @@ -145,7 +150,8 @@ subroutine gravity_wave_drag_moving_mountain_diagnostics_run( & call history_out_field('XPWP_SRC_MOVMTN', xpwp_src) call history_out_field('PSTEER_MOVMTN', p_steer) call history_out_field('PLAUNCH_MOVMTN', p_launch) - call history_out_field('ROO_BLOCK_MOVMTN', roo_block) + roo_block_real = merge(1._kind_phys, 0._kind_phys, roo_block) + call history_out_field('ROO_BLOCK_MOVMTN', roo_block_real) call history_out_field('MMTAUE', taucd_east) call history_out_field('MMTAUW', taucd_west) call history_out_field('MMTAUS', taucd_south) diff --git a/schemes/sima_diagnostics/gravity_wave_drag_moving_mountain_diagnostics.meta b/schemes/sima_diagnostics/gravity_wave_drag_moving_mountain_diagnostics.meta index 8e8a136d0..e0105c80f 100644 --- a/schemes/sima_diagnostics/gravity_wave_drag_moving_mountain_diagnostics.meta +++ b/schemes/sima_diagnostics/gravity_wave_drag_moving_mountain_diagnostics.meta @@ -21,6 +21,12 @@ [ccpp-arg-table] name = gravity_wave_drag_moving_mountain_diagnostics_run type = scheme +[ ncol ] + standard_name = horizontal_loop_extent + units = count + type = integer + dimensions = () + intent = in [ u ] standard_name = eastward_wind units = m s-1 @@ -167,8 +173,8 @@ intent = in [ roo_block ] standard_name = flag_for_near_surface_tilt_dominance_blocking_of_moving_mountain_gravity_wave_drag_source - units = 1 - type = real | kind = kind_phys + units = flag + type = logical dimensions = (horizontal_loop_extent) intent = in [ taucd_west ] From 4bb70460f2af89807d546f3403dc55d821f00932 Mon Sep 17 00:00:00 2001 From: Peter Hjort Lauritzen Date: Fri, 13 Nov 2026 13:24:41 -0700 Subject: [PATCH 6/8] Apply ridge GW vertical mixing once, with the summed diffusivity of all ridges gw_rdg_calc solved the implicit vertical diffusion of constituents and dry static energy once per ridge (and residual), each from the same input state, and summed the increments. Each solve alone is stable, but the sum over-mixes the column: a 2dz mode is multiplied by 1 - sum_n(1 - f_n) < -1 when several ridges mix strongly, so it flips sign and grows every step (ESCOMP/CAM issue #1753). Fix: keep each ridge's diffusivity, drag and kinetic-energy heating; drop its own q/DSE diffusion; after the ridge loop, do one implicit solve with the diffusivity summed over all ridges and the residual. All edits are tagged "! SUMK". This is the version used in the WACCM-MA ne30pg3 L135 validation runs. --- .../gravity_wave_drag_ridge.F90 | 45 +++++++++++++++++++ 1 file changed, 45 insertions(+) diff --git a/schemes/gravity_wave_drag/gravity_wave_drag_ridge.F90 b/schemes/gravity_wave_drag/gravity_wave_drag_ridge.F90 index 9798e17b1..8a4cbd695 100644 --- a/schemes/gravity_wave_drag/gravity_wave_drag_ridge.F90 +++ b/schemes/gravity_wave_drag/gravity_wave_drag_ridge.F90 @@ -1571,6 +1571,10 @@ subroutine gw_rdg_calc(band_oro, & use coords_1d, only: coords1d use gw_common, only: GWBand, gw_drag_prof, energy_change + use gw_diffusion, only: gw_diff_tend ! SUMK + use vdiff_lu_solver, only: fin_vol_lu_decomp ! SUMK + use linear_1d_operators, only: TriDiagDecomp ! SUMK + use physconst, only: gravit ! SUMK type(GWBand), intent(in) :: band_oro real(kind_phys), intent(in) :: vramp(:) @@ -1707,6 +1711,9 @@ subroutine gw_rdg_calc(band_oro, & ! Effective gravity wave diffusivity at interfaces. real(kind_phys) :: egwdffi(ncol, pver + 1) + ! SUMK: diffusivity summed over ridges + residual, (g*rho)^2, LU decomposition + real(kind_phys) :: egwdffi_sum(ncol, pver + 1), dpidz_sq(ncol, pver + 1) + type(TriDiagDecomp) :: decomp ! Temperature tendencies from diffusion and kinetic energy. real(kind_phys) :: dttdf(ncol, pver) @@ -1736,6 +1743,7 @@ subroutine gw_rdg_calc(band_oro, & utrdg = 0._kind_phys vtrdg = 0._kind_phys tau_diag = -9999._kind_phys + egwdffi_sum = 0._kind_phys ! SUMK do nn = 1, n_rdg kwvrdg = 0.001_kind_phys/(hwdth(:, nn) + 0.001_kind_phys) @@ -1772,6 +1780,15 @@ subroutine gw_rdg_calc(band_oro, & kwvrdg=kwvrdg, & satfac_in=1._kind_phys, lapply_vdiff=gw_rdg_do_vdiff, tau_diag=tau_diag) + ! SUMK: keep this ridge's K; drop its own diffusion of q and DSE (done once below) + if (gw_rdg_do_vdiff) then + egwdffi_sum = egwdffi_sum + egwdffi + do k = 1, pver + ttgw(:, k) = ttgw(:, k) - dttdf(:, k)*vramp(k) + end do + qtgw = 0._kind_phys + end if + ! Add the tendencies from each ridge to the totals. do k = 1, pver ! diagnostics @@ -1822,6 +1839,15 @@ subroutine gw_rdg_calc(band_oro, & kwvrdg=kwvrdg, & satfac_in=1._kind_phys, lapply_vdiff=gw_rdg_do_vdiff, tau_diag=tau_diag) + ! SUMK: keep the residual's K; its diffusion is done once below + if (gw_rdg_do_vdiff) then + egwdffi_sum = egwdffi_sum + egwdffi + do k = 1, pver + ttgw(:, k) = ttgw(:, k) - dttdf(:, k)*vramp(k) + end do + qtgw = 0._kind_phys + end if + ! Add the tendencies from isotropic residual to the totals. do k = 1, pver ! diagnostics @@ -1848,6 +1874,25 @@ subroutine gw_rdg_calc(band_oro, & end do end if ! end of residual variance calc + ! SUMK: one implicit diffusion of q and DSE with K summed over all ridges and + ! the residual (as gw_ediff already does over the wave spectrum). Summing the + ! per-ridge increments instead gives a 2dz factor 1-sum(1-f_n) < -1 (unstable). + if (gw_rdg_do_vdiff) then + dpidz_sq = (rhoi(:ncol, :pver + 1)*gravit)**2 + decomp = fin_vol_lu_decomp(dt, p%section([1, ncol], [1, pver]), & + coef_q_diff=egwdffi_sum*dpidz_sq) + do m = 1, pcnst + call gw_diff_tend(ncol, pver, pver, 1, q(:, :, m), dt, decomp, qtgw(:, :, m)) + q_tend(:ncol, :, m) = q_tend(:ncol, :, m) + qtgw(:, :, m) + end do + call gw_diff_tend(ncol, pver, pver, 1, dse, dt, decomp, dttdf) + do k = 1, pver + ttrdg(:, k) = ttrdg(:, k) + dttdf(:, k)*vramp(k) + s_tend(:ncol, k) = s_tend(:ncol, k) + dttdf(:, k)*vramp(k) + end do + call decomp%finalize() + end if + ! Calculate energy change for output to CAM's energy checker. call energy_change(dt, p, u, v, u_tend(:ncol, :), & v_tend(:ncol, :), s_tend(:ncol, :), de) From 6c738330b8b3a2e93e2cefbf6a5395f3990a8530 Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Fri, 13 Nov 2026 16:55:44 -0500 Subject: [PATCH 7/8] Cleanup comments --- .../gravity_wave_drag/gravity_wave_drag_ridge.F90 | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/schemes/gravity_wave_drag/gravity_wave_drag_ridge.F90 b/schemes/gravity_wave_drag/gravity_wave_drag_ridge.F90 index 8a4cbd695..bba8f9686 100644 --- a/schemes/gravity_wave_drag/gravity_wave_drag_ridge.F90 +++ b/schemes/gravity_wave_drag/gravity_wave_drag_ridge.F90 @@ -1711,7 +1711,7 @@ subroutine gw_rdg_calc(band_oro, & ! Effective gravity wave diffusivity at interfaces. real(kind_phys) :: egwdffi(ncol, pver + 1) - ! SUMK: diffusivity summed over ridges + residual, (g*rho)^2, LU decomposition + ! diffusivity summed over ridges + residual, (g*rho)^2, LU decomposition real(kind_phys) :: egwdffi_sum(ncol, pver + 1), dpidz_sq(ncol, pver + 1) type(TriDiagDecomp) :: decomp @@ -1780,7 +1780,7 @@ subroutine gw_rdg_calc(band_oro, & kwvrdg=kwvrdg, & satfac_in=1._kind_phys, lapply_vdiff=gw_rdg_do_vdiff, tau_diag=tau_diag) - ! SUMK: keep this ridge's K; drop its own diffusion of q and DSE (done once below) + ! keep this ridge's K but drop own diffusion of q and DSE which is done once below. if (gw_rdg_do_vdiff) then egwdffi_sum = egwdffi_sum + egwdffi do k = 1, pver @@ -1839,7 +1839,7 @@ subroutine gw_rdg_calc(band_oro, & kwvrdg=kwvrdg, & satfac_in=1._kind_phys, lapply_vdiff=gw_rdg_do_vdiff, tau_diag=tau_diag) - ! SUMK: keep the residual's K; its diffusion is done once below + ! keep K of residual, its diffusion is done once below if (gw_rdg_do_vdiff) then egwdffi_sum = egwdffi_sum + egwdffi do k = 1, pver @@ -1874,9 +1874,10 @@ subroutine gw_rdg_calc(band_oro, & end do end if ! end of residual variance calc - ! SUMK: one implicit diffusion of q and DSE with K summed over all ridges and - ! the residual (as gw_ediff already does over the wave spectrum). Summing the - ! per-ridge increments instead gives a 2dz factor 1-sum(1-f_n) < -1 (unstable). + ! do one implicit diffusion of q and DSE with K summed over all ridges + ! and the residual + ! if the per-ridge increments are summed instead it gives an unstable + ! 2dz factor 1-sum(1-f_n) < -1. if (gw_rdg_do_vdiff) then dpidz_sq = (rhoi(:ncol, :pver + 1)*gravit)**2 decomp = fin_vol_lu_decomp(dt, p%section([1, ncol], [1, pver]), & From 2d8e67cfbe0d7f7872d1fafb45bfd6743df0a18c Mon Sep 17 00:00:00 2001 From: Haipeng Lin Date: Fri, 13 Nov 2026 17:42:14 -0500 Subject: [PATCH 8/8] Cleanup roo comment --- .../gravity_wave_drag_moving_mountain.F90 | 32 +++++++++---------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/schemes/gravity_wave_drag/gravity_wave_drag_moving_mountain.F90 b/schemes/gravity_wave_drag/gravity_wave_drag_moving_mountain.F90 index af5b240b6..c1564709d 100644 --- a/schemes/gravity_wave_drag/gravity_wave_drag_moving_mountain.F90 +++ b/schemes/gravity_wave_drag/gravity_wave_drag_moving_mountain.F90 @@ -774,22 +774,6 @@ subroutine tilt_centroid_levels(tilt, u, v, pmid, zm, ncol, pver, & !------------------------------------------------------------------------ ! Flow-dependent steering/launch levels from |tilt| pressure centroids ! (for source_type=3-5). - ! - ! The centroid search can return in-bounds levels that are still meaningless - ! when the column's tilt profile is dominated by an isolated near-surface - ! spike rather than a free-tropospheric signal, as seen over Antarctica, - ! the Arctic Ocean and Greenland. - ! The "roo" criterion rejects such columns using the ratio of the - ! near-surface to the free-tropospheric maximum |tilt|, - ! roo = max(|tilt|, 0-500 m AGL) / (max(|tilt|, 500-15000 m AGL) + roo_eps), - ! and blocks the source where roo > roo_thresh. - ! This is an empirical criterion from moving-mountain source development; - ! it has no published reference. - ! roo_block is true only for this rejection, not for centroid-search failures. - ! - ! If the centroid search fails (returns -1) or roo blocks the column, - ! valid=.false. and the levels fall back to safe in-bounds values. - ! Callers must zero their source where valid=.false. !------------------------------------------------------------------------ integer, intent(in) :: ncol, pver real(kind_phys), intent(in) :: tilt(ncol,pver) @@ -799,7 +783,23 @@ subroutine tilt_centroid_levels(tilt, u, v, pmid, zm, ncol, pver, & integer, intent(out) :: steering_level(ncol), launch_level(ncol) ! valid launch_level < steering_level. real(kind_phys), intent(out) :: p_steer(ncol), p_launch(ncol) real(kind_phys), intent(out) :: usteer(ncol), vsteer(ncol) + + ! If the centroid search fails (returns -1) or "roo" (see below) blocks the column, + ! valid=.false. and the levels fall back to safe in-bounds values. + ! Callers must zero their source where valid=.false. logical, intent(out) :: valid(ncol) + + ! "Roo" criterion: + ! The centroid search can return in-bounds levels that are still meaningless + ! when the column's tilt profile is dominated by an isolated near-surface + ! spike rather than a free-tropospheric signal, as seen over Antarctica, + ! the Arctic Ocean and Greenland. + ! This empirical criterion rejects such columns using the ratio of the + ! near-surface to the free-tropospheric maximum |tilt|, + ! roo = max(|tilt|, 0-500 m AGL) / (max(|tilt|, 500-15000 m AGL) + roo_eps), + ! and blocks the source where roo > roo_thresh. + ! + ! roo_block is true only for this rejection, not for centroid-search failures. logical, intent(out) :: roo_block(ncol) ! Pressure bounds for the centroid search [Pa]: TOA down to a fraction