From af6762b141a2c94e2b2864bd9ebf1901e7927701 Mon Sep 17 00:00:00 2001 From: phamkh Date: Sun, 9 Nov 2025 20:28:45 -0700 Subject: [PATCH 1/2] Fixing Kp=0 bug for Gallagher --- src/base/earthhelper.F90 | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/base/earthhelper.F90 b/src/base/earthhelper.F90 index 26b89af6..eaed95e8 100644 --- a/src/base/earthhelper.F90 +++ b/src/base/earthhelper.F90 @@ -96,6 +96,11 @@ module earthhelper D = 0.0 + xfn = 0.0_rp + yfn = 0.0_rp + xfd = 0.0_rp + yfd = 0.0_rp + if (present(kpO)) then kp = nint(kpO) ! Gallagher model takes integer Kp else From 83f503863510016741ef167616643960c87a85c9 Mon Sep 17 00:00:00 2001 From: phamkh Date: Mon, 10 Nov 2025 06:59:02 -0700 Subject: [PATCH 2/2] Fixing Out of bound Gallagher bug --- src/base/earthhelper.F90 | 29 ++++++++++++++++++----------- 1 file changed, 18 insertions(+), 11 deletions(-) diff --git a/src/base/earthhelper.F90 b/src/base/earthhelper.F90 index eaed95e8..92bfe5d0 100644 --- a/src/base/earthhelper.F90 +++ b/src/base/earthhelper.F90 @@ -61,13 +61,24 @@ module earthhelper real(rp), intent(in) :: kp !Can only be 1,5 - if ( (kp>=1.0) .and. (kp<=5.0) ) then - kpDefault = kp - else if (kp>5.0) then - kpDefault = 5.0 - endif + kpDefault = GallagherKp(kp) end subroutine SetKp0 + + ! Routine to enforce kp to Gallagher's valid range [1,5] + pure function GallagherKp(kp) result(validKp) + real(rp), intent(in) :: kp + real(rp) :: validKp + + if (kp < 1.0_rp) then + validKp = 1.0_rp + else if (kp > 5.0_rp) then + validKp = 5.0_rp + else + validKp = kp + endif + end function GallagherKp + !Return 2D gallagher density afa r,phi (rad) function GallagherRP(r,phi,kpO) result(D) real(rp), intent(in) :: r,phi @@ -96,13 +107,9 @@ module earthhelper D = 0.0 - xfn = 0.0_rp - yfn = 0.0_rp - xfd = 0.0_rp - yfd = 0.0_rp - if (present(kpO)) then - kp = nint(kpO) ! Gallagher model takes integer Kp + !Can only be 1,5 + kp = nint(GallagherKp(kpO)) else kp = nint(kpDefault) endif