PBLH Corrections by hgopalan · Pull Request #3486 · erf-model/ERF

hgopalan · 2026-07-23T22:38:33Z

MRF and YSU depend on planetary boundary layer height to compute the turbulent viscosity. The variable becomes prognostic instead of diagnostic. So we have to modify the calculation and filling up the surface layer. Few other changes:

Fixed the low wind speed limits for MRF and YSU. A geostrophic wind speed of 0.1 m/s will be still stable.
Added PBLH smoothing between neighbouring cells and improved PBLH diagnosis for unstable ABL
Added QNSE stability condition for very stable ABL for MRF/YSUNew. This prevents the collapse of the very stable ABL under low geostrophic winds and high cooling rate.

…SUNew PBL schemes Co-authored-by: hgopalan <1108371+hgopalan@users.noreply.github.com>

Co-authored-by: hgopalan <1108371+hgopalan@users.noreply.github.com>

Add Vogelezang & Holtslag (1996) shear-correction and PBLH spatial smoothing for MRF and YSUNew PBL schemes

hgopalan · 2026-08-04T20:56:55Z

#3541 will be resolved here.

AMLattanzi · 2026-08-13T16:22:46Z

+
+        // Swap arrays for next iteration (or final result)
+        if (pass < passes - 1) {
+            pblh_fab.copy(pblh_temp, 0, 0, 1);


The copies at L482 and 485 need a template parameter for where they should be executed --- e.g., pblh_fab.copy<RunOn::Device>(pblh_temp, 0, 0, 1);

This should fix the GPU compilation failures shown in the workflow.

AMLattanzi · 2026-08-13T16:24:07Z

                  const int hgam_comp = (stagdir == 0) ? EddyDiff::HGAMU_v : EddyDiff::HGAMV_v;
                  const Real gam_hi = myhalf * (mu_turb(i,j,klo,hgam_comp) + mu_turb(i,j,klo+1,hgam_comp));
-                  RHS_a(i,j,klo) += Fact * gfac * rhoAlpha_hi * gam_hi / met_h_zeta_hi;
+                  RHS_a(i,j,klo) -= Fact * gfac * rhoAlpha_hi * gam_hi / met_h_zeta_hi;


All of the met_h_zeta divisions need removing. This includes the state and momentum diffusion.

AMLattanzi · 2026-08-13T16:24:46Z

                  const int hgam_comp = (stagdir == 0) ? EddyDiff::HGAMU_v : EddyDiff::HGAMV_v;
                  const Real gam_hi = myhalf * (mu_turb(i,j,klo,hgam_comp) + mu_turb(i,j,klo+1,hgam_comp));
-                  RHS_a(i,j,klo) += Fact * gfac * dz_inv * rhoAlpha_hi * gam_hi;
+                  RHS_a(i,j,klo) += Fact * gfac  * rhoAlpha_hi * gam_hi;


Need to flip the sign to -=

AMLattanzi

See comments at respective lines

AMLattanzi · 2026-08-17T17:22:52Z

                  const Real gam_hi  = myhalf * (gam_k + gam_kp1);
                  const Real gam_lo  = myhalf * (gam_k + gam_km1);
-                  RHS_a(i,j,k) += Fact * gfac * (rhoAlpha_hi * gam_hi / met_h_zeta_hi - rhoAlpha_lo * gam_lo / met_h_zeta_lo);
+                  RHS_a(i,j,k) -= Fact * gfac * (rhoAlpha_hi * gam_hi / met_h_zeta_hi - rhoAlpha_lo * gam_lo / met_h_zeta_lo);


The met_h_zeta still needs removing. There is no vertical derivative that needs correcting here.

The neutral log law on the tangential wind cannot shed heat from a hot face in calm air, which is why the day canonical's roof reached 340 K. Two switches, both off by default so every earlier result is unchanged: erf.ibseb.convective_velocity = deardorff puts a convective velocity scale into the wind the wall function sees, U_eff^2 = U_tan^2 + (beta w*)^2 (Beljaars), with w* = (g/theta H/(rho c_p) depth)^(1/3) from the previous step's flux out of the face. The depth is the mixed layer above a roof (z_i - z_face, floored at the building height) and the building height for a wall; z_i follows the day by the bulk Richardson diagnostic on the horizontal-mean profile (erf.ibseb.z_i_mode = bulk_ri), from the surface layer's pblh at the column (pblh), or fixed (z_i). erf.ibseb.stability_correction now iterates Dyer's similarity functions on the roofs to convergence on the face's own Obukhov length, under- relaxed as the surface layer's iteration is in erf-model erf-model#3486, seeded from the ground surface layer's 2D field at the column (erf.ibseb.obukhov_seed). Walls stay on the log law. The face's L stays its own because a roof in a separation zone or a sunlit wall can be in the opposite regime from the ground. The face dump gains w_star, Olen, z_i and h_bld; the summary gains w_star_max. Regtest SEB/Phase8_WallFunction: the cube in calm air under a strong sun; the neutral law sheds 0.2 W/m2 from a 342 K roof and the scale 270 W/m2; w*, the depth, u* and H follow the formulas to 1e-9; the roofs' L is negative and consistent with u* and theta* and the corrected log law to 1e-7; the bulk Richardson depth on a capped sounding is the first cell above the inversion. The default path is bit-identical to the phase 6 reference. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>

Conflicts in the MRF and YSUNew kernels, where development's PBLH work (erf-model#3486) meets this branch's immersed-boundary awareness (local surface ksurf and effective surface scales us_eff/ts_eff/qs_eff/ol_eff): - MRF passes 2 and 4: keep development's convective/shear w* blend, evaluated on the column's local surface (theta_v at ksurf) with the effective u* and theta*, as the branch does for every other surface term. - MRF K-profile: development's ghost-cell guard and notes, with the branch's k < ksurf (no diffusivity inside a solid) condition. - YSUNew K-profile: development's rho guard, plus the branch's ksurf and effective Obukhov length. - VH96 shear correction (new on development, off by default): use the effective u*. On columns without immersed cells it equals u*, so runs without erf.pbl_ib_aware are unchanged. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

…ore a garbage PBLH); restore #3961's HOL_abs bases (#3972) * PBL: keep MRF, YSUNew and YSU loops inside their MFIter tiles ComputeDiffusivityMRF and ComputeDiffusivityYSUNew build the corrected PBL height in a per-tile planar FArrayBox (growntilebox(1,1,0) under TileNoZ), then copied it into pblh_mf with a ParallelFor over mfi.validbox(). With the CPU default tile size (1024000 8 8) any box 16 or more cells wide in y is tiled, so the copy read past the tile's array. Debug and bound-check builds abort at step 1 with "(0,8,0,0) is out of bound (-1:16,-1:7,0:0,0:0)". Release builds silently store out-of-bounds values for every tile but the last in SurfaceLayer's pblh (the 2D plotfile pblh, the Beljaars w* term, the dust layer). Loop over mfi.tilebox() instead. It spans the full column under TileNoZ, and pblh_mf has no ghost cells, so the tiles together fill every cell before set_pblh copies it. The legacy ComputeDiffusivityYSU iterated with TilingIfNotGPU(), which also tiles in z (tile size 8), so its full-column AMREX_ALWAYS_ASSERT aborted at step 1 whenever nz >= 16. Use TileNoZ() like the other column schemes. Add CTest ABL_MRF_Tiling, ABL_YSUNew_Tiling and ABL_YSU_Tiling (Tests/RunTilingParity.cmake). The same 32x32x32 deck on four 16x16x32 boxes runs with tile size 1024000 8 8 and with tiling off, and the 3D and 2D plotfiles must agree. With the fix every field, pblh included, is bit-identical. On the old code the Release MRF and YSUNew runs give a pblh difference of 16 m (relative 1.0), the bound-check MRF run aborts in the tiled case, and YSU fails its assert in the tiled case. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * MRF: restore the |HOL| stability-function bases from #3961 #3961 wrote the unstable-arm bases in ComputeDiffusivityMRF as 1 + 16*|HOL|, which equals 1 - 16*HOL on that arm (HOL <= 0 there) and is at least 1 for every HOL, so the optimiser cannot hoist a pow of a negative base above the selection and raise FE_INVALID under amrex.fpe_trap_invalid. #3486 put back pow(max(1 - 16*HOL, 0.01), ...) in the Pass 4 wstar recompute and in the K-profile phiM and phit, leaving the HOL_abs declared for Pass 4 unused. That is an -Wunused-variable warning, which fails the CI warning check. Use HOL_abs again at those three sites. On the unstable arm max(1 - 16*HOL, 0.01) is 1 - 16*HOL (it is >= 1), so the results are unchanged; the tiling regression deck gives bit-identical plotfiles. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Tiling parity tests: register only with MPI off Windows, harden the script Review fixes for the tiling parity harness: - Tests/CTestList.cmake: register ABL_MRF_Tiling, ABL_YSUNew_Tiling and ABL_YSU_Tiling only under ERF_ENABLE_MPI AND NOT WIN32. The Windows job configures with MPI off and runs ctest -L regression, so MPIEXEC arrived empty, execute_process dropped it and tried to run "1". On Windows the erf_exec and fcompare paths are globs that only sh -c expands, and execute_process does not. - Tests/CTestList.cmake: note that ABL_YSU_Tiling covers the full-column assert only. Legacy YSU never calls set_pblh, so the 2D pblh holds the SurfaceLayer's initial value in both runs. - Tests/RunTilingParity.cmake: treat an empty -D value as missing, since add_test_tiling_parity always passes every argument. - Tests/RunTilingParity.cmake: remove the step-suffixed plotfiles (<run>_plt*) before each run. The old REMOVE_RECURSE named tiled_plt/tiled_plt2d, which never exist, so stale step-10 plotfiles could be compared. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> * Tiling parity tests: make the compared PBL height vary, and require it to Auditing the review findings on this PR showed a fifth weakness in the harness: amrex_fextrema on the untiled MRF and YSUNew plotfiles gave pblh (2D) and Lturb (3D) = 16 in every column. That is the floor 0.5*dz (dz is 32 m, not the 200 m the description implied): with a 15 m/s wind and the inversion above the domain top the bulk-Richardson search never crossed Ribcr. A parity check on a field that is one constant everywhere cannot see a tile bug that permutes columns, only one that reads garbage. - Tests/test_files/ABL_MRF_Tiling: put a 6 K inversion at 150-250 m, lower the wind to 5 m/s and use local MOST averaging (erf.most.average_policy = 1). The Rib crossing now lies inside the perturbed layer, so pblh ranges over 168-204 m (MRF) and 180-201 m (YSUNew), and u_star, t_star and Olen differ between columns too. - Tests/RunTilingParity.cmake: new VARYING_3D / VARYING_2D lists. For each named field the script runs amrex_fextrema on the untiled plotfile and fails unless min < max, so an agreement between two copies of a constant is reported as a failure, not a pass. On the step-0 2D plotfile, where pblh is still the 1e150 initial value, the guard stops with "pblh is uniform (1e+150)". - Tests/CTestList.cmake: pass FEXTREMA (built next to fcompare by the AMReX plotfile tools) and the lists. MRF and YSUNew check Lturb, Kmv, pblh and u_star. Legacy YSU never writes the surface-layer pblh, so its run drops pblh from erf.plot2d_vars_1 instead of comparing the constant, and checks Lturb, Kmv and u_star. Its cooling is -0.02 instead of -0.05: with the 5 m/s wind the stronger cooling stops the MOST iteration converging. With the old PBL sources (489a2aa) the tiled MRF and YSUNew runs still abort on the bound check and YSU still fails its full-column assert; with the fix all three pass with every field bit-identical. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> --------- Co-authored-by: Claude Opus 5 <noreply@anthropic.com> Co-authored-by: Ann Almgren <asalmgren@lbl.gov> Co-authored-by: Aaron M. Lattanzi <103702284+AMLattanzi@users.noreply.github.com>

hgopalan and others added 21 commits July 23, 2026 16:36

PBLH Corrections

e3da8ff

Space

b4d6f80

Space

55eb91e

Merge branch 'erf-model:development' into ERF-PBLH

1f752de

YSUNew

d9cf705

YSUNew

fe56af5

Merge branch 'development' into ERF-PBLH

14a1e1b

Merge branch 'development' into ERF-PBLH

874f8f6

Implement Vogelezang & Holtslag (1996) shear-correction for MRF and Y…

8dec6e9

…SUNew PBL schemes Co-authored-by: hgopalan <1108371+hgopalan@users.noreply.github.com>

Add PBLH spatial smoothing feature for both MRF and YSUNew PBL schemes

4b14fce

Co-authored-by: hgopalan <1108371+hgopalan@users.noreply.github.com>

Merge pull request #277 from hgopalan/copilot/erf-pblh

9142853

Add Vogelezang & Holtslag (1996) shear-correction and PBLH spatial smoothing for MRF and YSUNew PBL schemes

Smoothing should be model independent

a41ff68

QNSE

9c123a6

Fix QNSE turbChoice

d215faa

QNSE branch proper

afc19b3

Merge branch 'erf-model:development' into ERF-PBLH

9e586e2

QNSE Test Case

e87bc4a

Fix

e6ec425

Fix zero wind speed MOST

3f40378

Revert unintended submodule update

aa2a2b7

YSUNew

f0d1680

hgopalan mentioned this pull request Aug 4, 2026

Clarify and complete YSU momentum countergradient integration #3541

Closed

hgopalan added 4 commits August 4, 2026 15:36

Scale with PBLH

32ee34b

Merge branch 'development' into ERF-PBLH

308a63e

Merge branch 'development' into ERF-PBLH

2480285

Counter Gradient Comments

cbb5413

hgopalan mentioned this pull request Aug 5, 2026

Fix terrain YSU countergradient scaling #3543

Closed

hgopalan added 2 commits August 5, 2026 11:44

Merge branch 'development' into ERF-PBLH

89acd21

Merge branch 'erf-model:development' into ERF-PBLH

e9cdc45

Fix Implicit Diffusion Signs

097374c

AMLattanzi reviewed Aug 13, 2026

View reviewed changes

AMLattanzi self-requested a review August 13, 2026 16:29

AMLattanzi requested changes Aug 13, 2026

View reviewed changes

hgopalan added 5 commits August 13, 2026 11:54

Merge branch 'development' into ERF-PBLH

fb1fb5e

Update signs

0b54b09

Fix the denominator

b1fd669

Merge branch 'development' into ERF-PBLH

d39cebe

Merge branch 'development' into ERF-PBLH

30696b3

AMLattanzi reviewed Aug 17, 2026

View reviewed changes

AMLattanzi requested changes Aug 17, 2026

View reviewed changes

hgopalan added 3 commits August 17, 2026 12:19

Merge branch 'development' into ERF-PBLH

a19ecd1

Added RegTest and fixed a compile bug

19f645b

Added RegTest and fixed a compile bug

82d5c8c

AMLattanzi approved these changes Aug 17, 2026

View reviewed changes

hgopalan and others added 4 commits August 18, 2026 11:18

Merge branch 'erf-model:development' into ERF-PBLH

372a4f3

Merge branch 'development' into ERF-PBLH

407fe55

Merge branch 'erf-model:development' into ERF-PBLH

01c9551

Merge branch 'development' into ERF-PBLH

5a3e7a0

Merge branch 'development' into ERF-PBLH

a121cc5

asalmgren merged commit f4fce41 into erf-model:development Sep 9, 2026
19 of 47 checks passed

hgopalan deleted the ERF-PBLH branch September 11, 2026 21:43

Navigation Menu

Uh oh!

PBLH Corrections - #3486

PBLH Corrections#3486
asalmgren merged 45 commits into
erf-model:developmentfrom
hgopalan:ERF-PBLH

hgopalan commented Jul 23, 2026 •

edited

Loading

Uh oh!

hgopalan commented Aug 4, 2026

Uh oh!

AMLattanzi Aug 13, 2026 •

edited

Loading

Uh oh!

AMLattanzi Aug 13, 2026 •

edited

Loading

Uh oh!

hgopalan Aug 13, 2026

Uh oh!

AMLattanzi Aug 13, 2026

Uh oh!

hgopalan Aug 13, 2026

Uh oh!

AMLattanzi left a comment

Uh oh!

AMLattanzi Aug 17, 2026

Uh oh!

hgopalan Aug 17, 2026

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

Uh oh!

Conversation

hgopalan commented Jul 23, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

hgopalan commented Aug 4, 2026

Uh oh!

AMLattanzi Aug 13, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

AMLattanzi Aug 13, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

hgopalan Aug 13, 2026

Choose a reason for hiding this comment

Uh oh!

AMLattanzi Aug 13, 2026

Choose a reason for hiding this comment

Uh oh!

hgopalan Aug 13, 2026

Choose a reason for hiding this comment

Uh oh!

AMLattanzi left a comment

Choose a reason for hiding this comment

Uh oh!

AMLattanzi Aug 17, 2026

Choose a reason for hiding this comment

Uh oh!

hgopalan Aug 17, 2026

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

hgopalan commented Jul 23, 2026 •

edited

Loading

AMLattanzi Aug 13, 2026 •

edited

Loading

AMLattanzi Aug 13, 2026 •

edited

Loading