From 5df79c53745fde5d6c3340a2979b1429cd5892c1 Mon Sep 17 00:00:00 2001 From: Henrik Rydberg Date: Sat, 8 Oct 2011 20:30:28 +0200 Subject: Initial import of htcd system 1.0 Signed-off-by: Henrik Rydberg --- src/labat/grabgr.f | 32 ++++++++++++++++++++++++++++++++ 1 file changed, 32 insertions(+) create mode 100644 src/labat/grabgr.f (limited to 'src/labat/grabgr.f') diff --git a/src/labat/grabgr.f b/src/labat/grabgr.f new file mode 100644 index 0000000..1ed039c --- /dev/null +++ b/src/labat/grabgr.f @@ -0,0 +1,32 @@ + subroutine grabgr(r,gradup,graddn,i,imin,imax,done,gagd,gagup, + j gagdn) +c calculates grad(abs(grad(density))) +c input r : radial coordinate (array) +c input gradup : grad(density), spin up (array) +c input graddn : grad(density), spin down (array) +c input i : counter, i.e. current position = r(i) +c input imin : min value of i +c input imax : max value of i +c in/output done : true if abs(grad(density)) is splined +c output gagd : grad(abs(grad(density))) +c output gagup,gagdn : grad(abs(grad(density))), spin up and down + implicit logical (a-z) + logical done + double precision r,gradup,graddn,gagd,gagup,gagdn,temp,gag2,deriv + integer i,imin,imax,j + common/ggagag/temp,gag2 + dimension r(561),gradup(561),graddn(561) + dimension temp(2,561),gag2(2,561) + if (.not.done) then + do 10 j = imin,imax + temp(1,j) = dabs(gradup(j)) + temp(2,j) = dabs(graddn(j)) + 10 continue + call spline(r,imin+1,imax,temp,gag2) + done = .true. + endif + gagup = deriv(r,temp,gag2,1,i,imin,imax) + gagdn = deriv(r,temp,gag2,-1,i,imin,imax) + gagd = gagup + gagdn + return + end -- cgit v1.2.3