Initial import of htcd system 1.0

Signed-off-by: Henrik Rydberg <rydberg@euromail.se>

1 files changed, 32 insertions, 0 deletions

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