From 61605d3478252be8b413389cc6d7f9a01efad4b1 Mon Sep 17 00:00:00 2001 From: Michele Guerini Rocco Date: Mon, 18 Sep 2023 17:27:36 +0200 Subject: [PATCH] fix density calculation in analytical mode MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit In the analytical mode the density does not have a tail, so tail%end is not initialised. The density is exactly zero at ψ≥1. --- src/coreprofiles.f90 | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/coreprofiles.f90 b/src/coreprofiles.f90 index 3e0b7ba..014b7d6 100644 --- a/src/coreprofiles.f90 +++ b/src/coreprofiles.f90 @@ -65,9 +65,8 @@ contains dens = 0 ddens = 0 - ! Outside the tail end both density and its - ! derivatives are identically zero - if (psin >= tail%end .or. psin < 0) return + ! Bail out when ψ is not available + if (psin < 0) return if (iprof == 0) then ! Use the analytical model @@ -81,6 +80,10 @@ contains else ! Use the numerical data + ! Past the tail end both density and its + ! derivatives are identically zero + if (psin >= tail%end) return + if (psin < tail%start) then ! Use the interpolating spline when in range