Skip to content

Conversation

@vitorlavor
Copy link

As pointed out in Issue #1055, $c_{2m}$ and $c_{2h}$ from Eqs. 28a, b (HF08) are being being set to 0.5 - several tests were performed for multiple neighbourhoods and there were no apparent issues. This PR reverts the fixed values and uses the $c_{2m}$ and $c_{2h}$ calculated directly from Eqs. 28a, b (HF08)

@vitorlavor vitorlavor requested a review from sunt05 January 23, 2026 12:35
@vitorlavor vitorlavor linked an issue Jan 23, 2026 that may be closed by this pull request
Copy link

@sunt05 sunt05 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for working on Issue #1055. Removing the hardcoded c2 = 0.5 makes sense if the physics-based formula produces better results.

Clarification needed: Is the goal to always use the formula, or to use it except when numerically unsafe? If the latter, some guard logic should remain (see inline comments).

Test update required: This change affects model outputs but does not update the test baseline. Please:

# 1. Rebuild and run tests
source .venv/bin/activate
make dev
pytest test/core/test_sample_output.py -v

# 2. Verify outputs are physically reasonable (no NaN/Inf, sensible magnitudes)

# 3. Regenerate baseline
python scripts/gen_sample_output.py

# 4. Commit updated baseline with your changes
git add test/fixtures/data_test/sample_output.csv.gz
git add src/suews/src/suews_phys_rslprof.f95

c2 = 0.5
!c2 = 0.5

c2 = (kappa*(3.-(2.*beta**2.*Lc/phim_zh*dphi)))/(2.*beta*phim_zh - kappa)
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The old IF (phi_hatmZh >= 1.) block was not just forcing c2 = 0.5 — it also short-circuited the formula evaluation. By removing it, we expose this calculation to edge cases:

  • Denominator (2.*beta*phim_zh - kappa) could be zero or near-zero
  • Large c2 values will overflow EXP(c2/2.)
  • Even when (1.-phi_hatmZh) = 0, multiplying by NaN/Inf still produces NaN

Consider keeping a guard:

IF (phi_hatmZh >= 1.) THEN
   cm = 0.
ELSE
   c2 = (kappa*(3.-(2.*beta**2.*Lc/phim_zh*dphi)))/(2.*beta*phim_zh - kappa)
   cm = (1.-phi_hatmZh)*EXP(c2/2.)
END IF

! force c2h to 0.5 for better stability. TS 14 Jul 2020
! TODO: a more proper threshold needs to be determined
c2h = 0.5
!c2h = 0.5
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same numerical stability concern as c2/cm above — the denominator (2.*beta*phih_zh - kappa*Scc) can cause division errors, and EXP(c2h/2.) can overflow.

c2 = (kappa*(3.-(2.*beta**2.*Lc/phim_zh*dphi)))/(2.*beta*phim_zh - kappa)
END IF
! Removing fixed c2m value #Issue 1055
!IF (phi_hatmZh >= 1.) THEN
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please delete this commented-out block. Git preserves history, so there is no need to keep old logic as comments — it adds noise and maintenance burden.

! ! if very unstable this might cause some high values of psihat_z
! c2 = (kappa*(3.-(2.*beta**2.*Lc/phim_zh*dphi)))/(2.*beta*phim_zh - kappa)
!END IF
! force c2 to 0.5 for better stability. TS 14 Jul 2020
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These comments now contradict the code behaviour. Please remove or update them.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

Complementary revision for displacement height #302

2 participants