Added thorn TGRtensor. This is a utility thorn from the TGR code. It

contains generic tensor operations, some useful constants, code for
numeric derivatives, interfaces to Lapack and BLAS, and interfaces for
Cactus routines.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinUtils/TGRtensor/trunk@2 b716e942-a2de-43ad-8f52-f3dfe468e4e7

25 files changed, 2430 insertions, 0 deletions

diff --git a/COPYING b/COPYING
new file mode 100644
index 0000000..5b6e7c6
--- /dev/null
+++ b/COPYING
@@ -0,0 +1,340 @@
+		    GNU GENERAL PUBLIC LICENSE
+		       Version 2, June 1991
+
+ Copyright (C) 1989, 1991 Free Software Foundation, Inc.
+                       59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+			    Preamble
+
+  The licenses for most software are designed to take away your
+freedom to share and change it.  By contrast, the GNU General Public
+License is intended to guarantee your freedom to share and change free
+software--to make sure the software is free for all its users.  This
+General Public License applies to most of the Free Software
+Foundation's software and to any other program whose authors commit to
+using it.  (Some other Free Software Foundation software is covered by
+the GNU Library General Public License instead.)  You can apply it to
+your programs, too.
+
+  When we speak of free software, we are referring to freedom, not
+price.  Our General Public Licenses are designed to make sure that you
+have the freedom to distribute copies of free software (and charge for
+this service if you wish), that you receive source code or can get it
+if you want it, that you can change the software or use pieces of it
+in new free programs; and that you know you can do these things.
+
+  To protect your rights, we need to make restrictions that forbid
+anyone to deny you these rights or to ask you to surrender the rights.
+These restrictions translate to certain responsibilities for you if you
+distribute copies of the software, or if you modify it.
+
+  For example, if you distribute copies of such a program, whether
+gratis or for a fee, you must give the recipients all the rights that
+you have.  You must make sure that they, too, receive or can get the
+source code.  And you must show them these terms so they know their
+rights.
+
+  We protect your rights with two steps: (1) copyright the software, and
+(2) offer you this license which gives you legal permission to copy,
+distribute and/or modify the software.
+
+  Also, for each author's protection and ours, we want to make certain
+that everyone understands that there is no warranty for this free
+software.  If the software is modified by someone else and passed on, we
+want its recipients to know that what they have is not the original, so
+that any problems introduced by others will not reflect on the original
+authors' reputations.
+
+  Finally, any free program is threatened constantly by software
+patents.  We wish to avoid the danger that redistributors of a free
+program will individually obtain patent licenses, in effect making the
+program proprietary.  To prevent this, we have made it clear that any
+patent must be licensed for everyone's free use or not licensed at all.
+
+  The precise terms and conditions for copying, distribution and
+modification follow.
+
+		    GNU GENERAL PUBLIC LICENSE
+   TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
+
+  0. This License applies to any program or other work which contains
+a notice placed by the copyright holder saying it may be distributed
+under the terms of this General Public License.  The "Program", below,
+refers to any such program or work, and a "work based on the Program"
+means either the Program or any derivative work under copyright law:
+that is to say, a work containing the Program or a portion of it,
+either verbatim or with modifications and/or translated into another
+language.  (Hereinafter, translation is included without limitation in
+the term "modification".)  Each licensee is addressed as "you".
+
+Activities other than copying, distribution and modification are not
+covered by this License; they are outside its scope.  The act of
+running the Program is not restricted, and the output from the Program
+is covered only if its contents constitute a work based on the
+Program (independent of having been made by running the Program).
+Whether that is true depends on what the Program does.
+
+  1. You may copy and distribute verbatim copies of the Program's
+source code as you receive it, in any medium, provided that you
+conspicuously and appropriately publish on each copy an appropriate
+copyright notice and disclaimer of warranty; keep intact all the
+notices that refer to this License and to the absence of any warranty;
+and give any other recipients of the Program a copy of this License
+along with the Program.
+
+You may charge a fee for the physical act of transferring a copy, and
+you may at your option offer warranty protection in exchange for a fee.
+
+  2. You may modify your copy or copies of the Program or any portion
+of it, thus forming a work based on the Program, and copy and
+distribute such modifications or work under the terms of Section 1
+above, provided that you also meet all of these conditions:
+
+    a) You must cause the modified files to carry prominent notices
+    stating that you changed the files and the date of any change.
+
+    b) You must cause any work that you distribute or publish, that in
+    whole or in part contains or is derived from the Program or any
+    part thereof, to be licensed as a whole at no charge to all third
+    parties under the terms of this License.
+
+    c) If the modified program normally reads commands interactively
+    when run, you must cause it, when started running for such
+    interactive use in the most ordinary way, to print or display an
+    announcement including an appropriate copyright notice and a
+    notice that there is no warranty (or else, saying that you provide
+    a warranty) and that users may redistribute the program under
+    these conditions, and telling the user how to view a copy of this
+    License.  (Exception: if the Program itself is interactive but
+    does not normally print such an announcement, your work based on
+    the Program is not required to print an announcement.)
+
+These requirements apply to the modified work as a whole.  If
+identifiable sections of that work are not derived from the Program,
+and can be reasonably considered independent and separate works in
+themselves, then this License, and its terms, do not apply to those
+sections when you distribute them as separate works.  But when you
+distribute the same sections as part of a whole which is a work based
+on the Program, the distribution of the whole must be on the terms of
+this License, whose permissions for other licensees extend to the
+entire whole, and thus to each and every part regardless of who wrote it.
+
+Thus, it is not the intent of this section to claim rights or contest
+your rights to work written entirely by you; rather, the intent is to
+exercise the right to control the distribution of derivative or
+collective works based on the Program.
+
+In addition, mere aggregation of another work not based on the Program
+with the Program (or with a work based on the Program) on a volume of
+a storage or distribution medium does not bring the other work under
+the scope of this License.
+
+  3. You may copy and distribute the Program (or a work based on it,
+under Section 2) in object code or executable form under the terms of
+Sections 1 and 2 above provided that you also do one of the following:
+
+    a) Accompany it with the complete corresponding machine-readable
+    source code, which must be distributed under the terms of Sections
+    1 and 2 above on a medium customarily used for software interchange; or,
+
+    b) Accompany it with a written offer, valid for at least three
+    years, to give any third party, for a charge no more than your
+    cost of physically performing source distribution, a complete
+    machine-readable copy of the corresponding source code, to be
+    distributed under the terms of Sections 1 and 2 above on a medium
+    customarily used for software interchange; or,
+
+    c) Accompany it with the information you received as to the offer
+    to distribute corresponding source code.  (This alternative is
+    allowed only for noncommercial distribution and only if you
+    received the program in object code or executable form with such
+    an offer, in accord with Subsection b above.)
+
+The source code for a work means the preferred form of the work for
+making modifications to it.  For an executable work, complete source
+code means all the source code for all modules it contains, plus any
+associated interface definition files, plus the scripts used to
+control compilation and installation of the executable.  However, as a
+special exception, the source code distributed need not include
+anything that is normally distributed (in either source or binary
+form) with the major components (compiler, kernel, and so on) of the
+operating system on which the executable runs, unless that component
+itself accompanies the executable.
+
+If distribution of executable or object code is made by offering
+access to copy from a designated place, then offering equivalent
+access to copy the source code from the same place counts as
+distribution of the source code, even though third parties are not
+compelled to copy the source along with the object code.
+
+  4. You may not copy, modify, sublicense, or distribute the Program
+except as expressly provided under this License.  Any attempt
+otherwise to copy, modify, sublicense or distribute the Program is
+void, and will automatically terminate your rights under this License.
+However, parties who have received copies, or rights, from you under
+this License will not have their licenses terminated so long as such
+parties remain in full compliance.
+
+  5. You are not required to accept this License, since you have not
+signed it.  However, nothing else grants you permission to modify or
+distribute the Program or its derivative works.  These actions are
+prohibited by law if you do not accept this License.  Therefore, by
+modifying or distributing the Program (or any work based on the
+Program), you indicate your acceptance of this License to do so, and
+all its terms and conditions for copying, distributing or modifying
+the Program or works based on it.
+
+  6. Each time you redistribute the Program (or any work based on the
+Program), the recipient automatically receives a license from the
+original licensor to copy, distribute or modify the Program subject to
+these terms and conditions.  You may not impose any further
+restrictions on the recipients' exercise of the rights granted herein.
+You are not responsible for enforcing compliance by third parties to
+this License.
+
+  7. If, as a consequence of a court judgment or allegation of patent
+infringement or for any other reason (not limited to patent issues),
+conditions are imposed on you (whether by court order, agreement or
+otherwise) that contradict the conditions of this License, they do not
+excuse you from the conditions of this License.  If you cannot
+distribute so as to satisfy simultaneously your obligations under this
+License and any other pertinent obligations, then as a consequence you
+may not distribute the Program at all.  For example, if a patent
+license would not permit royalty-free redistribution of the Program by
+all those who receive copies directly or indirectly through you, then
+the only way you could satisfy both it and this License would be to
+refrain entirely from distribution of the Program.
+
+If any portion of this section is held invalid or unenforceable under
+any particular circumstance, the balance of the section is intended to
+apply and the section as a whole is intended to apply in other
+circumstances.
+
+It is not the purpose of this section to induce you to infringe any
+patents or other property right claims or to contest validity of any
+such claims; this section has the sole purpose of protecting the
+integrity of the free software distribution system, which is
+implemented by public license practices.  Many people have made
+generous contributions to the wide range of software distributed
+through that system in reliance on consistent application of that
+system; it is up to the author/donor to decide if he or she is willing
+to distribute software through any other system and a licensee cannot
+impose that choice.
+
+This section is intended to make thoroughly clear what is believed to
+be a consequence of the rest of this License.
+
+  8. If the distribution and/or use of the Program is restricted in
+certain countries either by patents or by copyrighted interfaces, the
+original copyright holder who places the Program under this License
+may add an explicit geographical distribution limitation excluding
+those countries, so that distribution is permitted only in or among
+countries not thus excluded.  In such case, this License incorporates
+the limitation as if written in the body of this License.
+
+  9. The Free Software Foundation may publish revised and/or new versions
+of the General Public License from time to time.  Such new versions will
+be similar in spirit to the present version, but may differ in detail to
+address new problems or concerns.
+
+Each version is given a distinguishing version number.  If the Program
+specifies a version number of this License which applies to it and "any
+later version", you have the option of following the terms and conditions
+either of that version or of any later version published by the Free
+Software Foundation.  If the Program does not specify a version number of
+this License, you may choose any version ever published by the Free Software
+Foundation.
+
+  10. If you wish to incorporate parts of the Program into other free
+programs whose distribution conditions are different, write to the author
+to ask for permission.  For software which is copyrighted by the Free
+Software Foundation, write to the Free Software Foundation; we sometimes
+make exceptions for this.  Our decision will be guided by the two goals
+of preserving the free status of all derivatives of our free software and
+of promoting the sharing and reuse of software generally.
+
+			    NO WARRANTY
+
+  11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY
+FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW.  EXCEPT WHEN
+OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES
+PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED
+OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.  THE ENTIRE RISK AS
+TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU.  SHOULD THE
+PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING,
+REPAIR OR CORRECTION.
+
+  12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR
+REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES,
+INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING
+OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED
+TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY
+YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER
+PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE
+POSSIBILITY OF SUCH DAMAGES.
+
+		     END OF TERMS AND CONDITIONS
+
+	    How to Apply These Terms to Your New Programs
+
+  If you develop a new program, and you want it to be of the greatest
+possible use to the public, the best way to achieve this is to make it
+free software which everyone can redistribute and change under these terms.
+
+  To do so, attach the following notices to the program.  It is safest
+to attach them to the start of each source file to most effectively
+convey the exclusion of warranty; and each file should have at least
+the "copyright" line and a pointer to where the full notice is found.
+
+    <one line to give the program's name and a brief idea of what it does.>
+    Copyright (C) <year>  <name of author>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program; if not, write to the Free Software
+    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+
+
+Also add information on how to contact you by electronic and paper mail.
+
+If the program is interactive, make it output a short notice like this
+when it starts in an interactive mode:
+
+    Gnomovision version 69, Copyright (C) year name of author
+    Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
+    This is free software, and you are welcome to redistribute it
+    under certain conditions; type `show c' for details.
+
+The hypothetical commands `show w' and `show c' should show the appropriate
+parts of the General Public License.  Of course, the commands you use may
+be called something other than `show w' and `show c'; they could even be
+mouse-clicks or menu items--whatever suits your program.
+
+You should also get your employer (if you work as a programmer) or your
+school, if any, to sign a "copyright disclaimer" for the program, if
+necessary.  Here is a sample; alter the names:
+
+  Yoyodyne, Inc., hereby disclaims all copyright interest in the program
+  `Gnomovision' (which makes passes at compilers) written by James Hacker.
+
+  <signature of Ty Coon>, 1 April 1989
+  Ty Coon, President of Vice
+
+This General Public License does not permit incorporating your program into
+proprietary programs.  If your program is a subroutine library, you may
+consider it more useful to permit linking proprietary applications with the
+library.  If this is what you want to do, use the GNU Library General
+Public License instead of this License.
diff --git a/README b/README
new file mode 100644
index 0000000..30b1a69
--- /dev/null
+++ b/README
@@ -0,0 +1,8 @@
+Cactus Code Thorn TGRtensor
+Authors    : Erik Schnetter <schnetter@uni-tuebingen.de>
+CVS info   : $Header$
+--------------------------------------------------------------------------
+
+Purpose of the thorn:
+
+It contains generic tensor operations and other generic stuff.
diff --git a/interface.ccl b/interface.ccl
new file mode 100644
index 0000000..5eb7d12
--- /dev/null
+++ b/interface.ccl
@@ -0,0 +1,4 @@
+# Interface definition for thorn TGRtensor
+# $Header$
+
+IMPLEMENTS: TGRtensor
diff --git a/param.ccl b/param.ccl
new file mode 100644
index 0000000..e5004ff
--- /dev/null
+++ b/param.ccl
@@ -0,0 +1,2 @@
+# Parameter definitions for thorn TGRtensor
+# $Header$
diff --git a/schedule.ccl b/schedule.ccl
new file mode 100644
index 0000000..b6b05ea
--- /dev/null
+++ b/schedule.ccl
@@ -0,0 +1,2 @@
+# Schedule definitions for thorn TGRtensor
+# $Header$
diff --git a/src/cactus.F90 b/src/cactus.F90
new file mode 100644
index 0000000..5d0c278
--- /dev/null
+++ b/src/cactus.F90
@@ -0,0 +1,389 @@
+! $Header$
+
+#include "cctk.h"
+
+module cactus
+  implicit none
+  public
+  
+  interface
+     
+     ! from Cactus:
+     
+     subroutine cctk_barrier (ierr, cctkgh)
+       implicit none
+       integer      ierr
+       CCTK_POINTER cctkgh
+     end subroutine cctk_barrier
+     
+     subroutine cctk_coordrange (ierr, cctkgh, lower, upper, dir, name, &
+          &                      systemname)
+       implicit none
+       integer ierr
+       CCTK_POINTER  cctkgh
+       CCTK_REAL     lower
+       CCTK_REAL     upper
+       integer       dir
+       character*(*) name
+       character*(*) systemname
+     end subroutine cctk_coordrange
+     
+     subroutine cctk_coordsystemhandle (handle, coordsystem)
+       implicit none
+       integer       handle
+       character*(*) coordsystem
+     end subroutine cctk_coordsystemhandle
+     
+     function cctk_equals (arg1, arg2)
+       implicit none
+       integer       cctk_equals
+       CCTK_POINTER  arg1
+       character*(*) arg2
+     end function cctk_equals
+     
+     subroutine cctk_fortranstring (clength, cstring, fortstring)
+       implicit none
+       CCTK_INT      clength    ! intent(out)
+       CCTK_POINTER  cstring
+       character*(*) fortstring
+     end subroutine cctk_fortranstring
+     
+     subroutine cctk_groupbboxgn (ierr, cctkgh, nbbox, bbox, groupname)
+       implicit none
+       integer       ierr
+       CCTK_POINTER  cctkgh
+       integer       nbbox
+       integer       bbox(nbbox)
+       character*(*) groupname
+     end subroutine cctk_groupbboxgn
+     
+     subroutine cctk_groupgshgn (ierr, cctkgh, ngsh, gsh, groupname)
+       implicit none
+       integer       ierr
+       CCTK_POINTER  cctkgh
+       integer       ngsh
+       integer       gsh(ngsh)
+       character*(*) groupname
+     end subroutine cctk_groupgshgn
+     
+     subroutine cctk_grouplbndgn (ierr, cctkgh, nlbnd, lbnd, groupname)
+       implicit none
+       integer       ierr
+       CCTK_POINTER  cctkgh
+       integer       nlbnd
+       integer       lbnd(nlbnd)
+       character*(*) groupname
+     end subroutine cctk_grouplbndgn
+     
+     subroutine cctk_grouplshgn (ierr, cctkgh, nlsh, lsh, groupname)
+       implicit none
+       integer       ierr
+       CCTK_POINTER  cctkgh
+       integer       nlsh
+       integer       lsh(nlsh)
+       character*(*) groupname
+     end subroutine cctk_grouplshgn
+     
+     subroutine cctk_groupnghostzonesgn &
+          (ierr, cctkgh, nnghostzones, nghostzones, groupname)
+       implicit none
+       integer       ierr
+       CCTK_POINTER  cctkgh
+       integer       nnghostzones
+       integer       nghostzones(nnghostzones)
+       character*(*) groupname
+     end subroutine cctk_groupnghostzonesgn
+     
+     subroutine cctk_groupubndgn (ierr, cctkgh, nubnd, ubnd, groupname)
+       implicit none
+       integer       ierr
+       CCTK_POINTER  cctkgh
+       integer       nubnd
+       integer       ubnd(nubnd)
+       character*(*) groupname
+     end subroutine cctk_groupubndgn
+     
+     subroutine cctk_info (thorn, message)
+       implicit none
+       character*(*) thorn
+       character*(*) message
+     end subroutine cctk_info
+     
+     subroutine cctk_interphandle (handle, interp)
+       implicit none
+       integer       handle
+       character*(*) interp
+     end subroutine cctk_interphandle
+     
+     function cctk_isthornactive (name)
+       implicit none
+       integer       cctk_isthornactive
+       character*(*) name
+     end function cctk_isthornactive
+     
+     function cctk_myproc (cctkgh)
+       implicit none
+       integer      cctk_myproc
+       CCTK_POINTER cctkgh
+     end function cctk_myproc
+     
+     function cctk_nprocs (cctkgh)
+       implicit none
+       integer      cctk_nprocs
+       CCTK_POINTER cctkgh
+     end function cctk_nprocs
+     
+     subroutine cctk_paramwarn (thorn, message)
+       implicit none
+       character*(*) thorn
+       character*(*) message
+     end subroutine cctk_paramwarn
+     
+     subroutine cctk_reductionarrayhandle (handle, reduction)
+       implicit none
+       integer       handle
+       character*(*) reduction
+     end subroutine cctk_reductionarrayhandle
+     
+     subroutine cctk_registerbanner (ierr, banner)
+       implicit none
+       integer       ierr
+       character*(*) banner
+     end subroutine cctk_registerbanner
+     
+     subroutine cctk_syncgroup (ierr, cctkgh, groupname)
+       implicit none
+       integer       ierr
+       CCTK_POINTER  cctkgh
+       character*(*) groupname
+     end subroutine cctk_syncgroup
+     
+     subroutine cctk_syncgroupi (ierr, cctkgh, group)
+       implicit none
+       integer      ierr
+       CCTK_POINTER cctkgh
+       integer      group
+     end subroutine cctk_syncgroupi
+     
+     subroutine cctk_syncgroupwithvar (ierr, cctkgh, varname)
+       implicit none
+       integer       ierr
+       CCTK_POINTER  cctkgh
+       character*(*) varname
+     end subroutine cctk_syncgroupwithvar
+     
+     subroutine cctk_syncgroupwithvari (ierr, cctkgh, var)
+       implicit none
+       integer      ierr
+       CCTK_POINTER cctkgh
+       integer      var
+     end subroutine cctk_syncgroupwithvari
+     
+     subroutine cctk_varindex (vi, vname)
+       implicit none
+       integer       vi
+       character*(*) vname
+     end subroutine cctk_varindex
+     
+     subroutine cctk_vartypei (vartype, varindex)
+       implicit none
+       integer vartype
+       integer varindex
+     end subroutine cctk_vartypei
+     
+     subroutine cctk_warn (level, line, file, thorn, message)
+       implicit none
+       integer       level
+       integer       line
+       character*(*) file
+       character*(*) thorn
+       character*(*) message
+     end subroutine cctk_warn
+     
+     
+     
+     ! from thorn CactusBase/Boundary:
+     subroutine bndcopygn (ierr, cctkgh, sw, togroupname, fromgroupname)
+       implicit none
+       integer       ierr
+       CCTK_POINTER  cctkgh
+       integer       sw(3)
+       character*(*) togroupname
+       character*(*) fromgroupname
+     end subroutine bndcopygn
+     
+     subroutine bndflatgn (ierr, cctkgh, sw, groupname)
+       implicit none
+       integer       ierr
+       CCTK_POINTER  cctkgh
+       integer       sw(3)
+       character*(*) groupname
+     end subroutine bndflatgn
+     
+     subroutine bndradiativegn (ierr, cctkgh, sw, var0, v0, &
+          &                     togroupname, fromgroupname)
+       implicit none
+       integer       ierr
+       CCTK_POINTER  cctkgh
+       integer       sw(3)
+       CCTK_REAL     var0
+       CCTK_REAL     v0
+       character*(*) togroupname
+       character*(*) fromgroupname
+     end subroutine bndradiativegn
+     
+     subroutine bndradiativevn (ierr, cctkgh, sw, var0, v0, &
+          &                     tovarname, fromvarname)
+       implicit none
+       integer       ierr
+       CCTK_POINTER  cctkgh
+       integer       sw(3)
+       CCTK_REAL     var0
+       CCTK_REAL     v0
+       character*(*) tovarname
+       character*(*) fromvarname
+     end subroutine bndradiativevn
+     
+     subroutine bndrobinvn (ierr, cctkgh, sw, finf, npow, varname)
+       implicit none
+       integer       ierr
+       CCTK_POINTER  cctkgh
+       integer       sw(3)
+       CCTK_REAL     finf
+       integer       npow
+       character*(*) varname
+     end subroutine bndrobinvn
+     
+     subroutine bndscalargn (ierr, cctkgh, sw, var0, groupname)
+       implicit none
+       integer       ierr
+       CCTK_POINTER  cctkgh
+       integer       sw(3)
+       CCTK_REAL     var0
+       character*(*) groupname
+     end subroutine bndscalargn
+     
+     subroutine bndscalarvn (ierr, cctkgh, sw, var0, varname)
+       implicit none
+       integer       ierr
+       CCTK_POINTER  cctkgh
+       integer       sw(3)
+       CCTK_REAL     var0
+       character*(*) varname
+     end subroutine bndscalarvn
+     
+     
+     
+     ! from thorn CactusBase/CartGrid3D:
+     
+     subroutine cartsymgn (ierr, cctkgh, groupname)
+       implicit none
+       integer       ierr
+       CCTK_POINTER  cctkgh
+       character*(*) groupname
+     end subroutine cartsymgn
+     
+     subroutine setcartsymvn (ierr, cctkgh, sw, varname)
+       implicit none
+       integer       ierr
+       CCTK_POINTER  cctkgh
+       integer       sw(3)
+       character*(*) varname
+     end subroutine setcartsymvn
+     
+     
+     
+     ! from thorn AlphaThorns/Cart3D:
+     
+     subroutine cart3dsymvi (ierr, cctkgh, vi)
+       implicit none
+       integer       ierr
+       CCTK_POINTER  cctkgh
+       integer       vi
+     end subroutine cart3dsymvi
+     
+     subroutine cart3dsymvn (ierr, cctkgh, vn)
+       implicit none
+       integer       ierr
+       CCTK_POINTER  cctkgh
+       character*(*) vn
+     end subroutine cart3dsymvn
+     
+     subroutine cart3dsymgi (ierr, cctkgh, gi)
+       implicit none
+       integer       ierr
+       CCTK_POINTER  cctkgh
+       integer       gi
+     end subroutine cart3dsymgi
+     
+     subroutine cart3dsymgn (ierr, cctkgh, gn)
+       implicit none
+       integer       ierr
+       CCTK_POINTER  cctkgh
+       character*(*) gn
+     end subroutine cart3dsymgn
+     
+     subroutine cart3dsettensortypevi (ierr, cctkgh, nvars, vi, typename)
+       implicit none
+       integer       ierr
+       CCTK_POINTER  cctkgh
+       integer       nvars
+       integer       vi(nvars)
+       character*(*) typename
+     end subroutine cart3dsettensortypevi
+     
+     subroutine cart3dsettensortypevn (ierr, cctkgh, vn, typename)
+       implicit none
+       integer       ierr
+       CCTK_POINTER  cctkgh
+       character*(*) vn
+       character*(*) typename
+     end subroutine cart3dsettensortypevn
+     
+     subroutine cart3dgetsymmetriesvi (ierr, cctkgh, sym, vi)
+       implicit none
+       integer      ierr
+       CCTK_POINTER cctkgh
+       integer      sym(6)
+       integer      vi
+     end subroutine cart3dgetsymmetriesvi
+     
+     subroutine cart3dgetsymmetriesvn (ierr, cctkgh, sym, vn)
+       implicit none
+       integer       ierr
+       CCTK_POINTER  cctkgh
+       integer       sym(6)
+       character*(*) vn
+     end subroutine cart3dgetsymmetriesvn
+     
+     subroutine cart3dgetsymmetryboundaries (ierr, cctkgh, symbnd)
+       implicit none
+       integer      ierr
+       CCTK_POINTER cctkgh
+       integer      symbnd(6)
+     end subroutine cart3dgetsymmetryboundaries
+     
+     subroutine cart3dgetsymmetric (ierr, cctkgh, symmetric)
+       implicit none
+       integer      ierr
+       CCTK_POINTER cctkgh
+       integer      symmetric(6)
+     end subroutine cart3dgetsymmetric
+     
+     subroutine cart3dgetstaggered (ierr, cctkgh, staggered)
+       implicit none
+       integer      ierr
+       CCTK_POINTER cctkgh
+       integer      staggered(6)
+     end subroutine cart3dgetstaggered
+     
+     subroutine cart3dgetperiodic (ierr, cctkgh, periodic)
+       implicit none
+       integer      ierr
+       CCTK_POINTER cctkgh
+       integer      periodic(6)
+     end subroutine cart3dgetperiodic
+     
+end interface
+  
+end module cactus
diff --git a/src/constants.F90 b/src/constants.F90
new file mode 100644
index 0000000..d9b948d
--- /dev/null
+++ b/src/constants.F90
@@ -0,0 +1,52 @@
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+module constants
+  implicit none
+  DECLARE_CCTK_PARAMETERS
+  private
+  
+  public delta2, delta3, delta4
+  public eta4
+  public epsilon2, epsilon3
+  
+  public pi
+  
+  
+  
+  CCTK_REAL, parameter :: zero = 0
+  
+  integer, parameter :: rk = kind(zero)
+  
+  
+  CCTK_REAL, parameter :: delta2(2,2) &
+       = reshape((/ 1,0,   0,1 /), (/2,2/))
+  
+  CCTK_REAL, parameter :: delta3(3,3) &
+       = reshape((/ 1,0,0,   0,1,0,   0,0,1 /), (/3,3/))
+  
+  CCTK_REAL, parameter :: delta4(0:3,0:3) &
+       = reshape((/ 1,0,0,0,   0,1,0,0,   0,0,1,0,   0,0,0,1 /), (/4,4/))
+  
+  
+  
+  CCTK_REAL, parameter :: eta4(0:3,0:3) &
+       = reshape((/ -1,0,0,0,   0,1,0,0,   0,0,1,0,   0,0,0,1 /), (/4,4/))
+  
+  
+  
+  CCTK_REAL, parameter :: epsilon2(2,2) &
+       = reshape ((/ 0,+1,   -1,0 /), (/2,2/))
+  
+  CCTK_REAL, parameter :: epsilon3(3,3,3) &
+       = reshape ((/ 0,0,0,    0,0,+1,   0,-1,0, &
+       &             0,0,-1,   0,0,0,    +1,0,0, &
+       &             0,+1,0,   -1,0,0,   0,0,0   /), (/3,3,3/))
+  
+  
+  
+  CCTK_REAL, parameter :: pi = 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068_rk
+  
+end module constants
diff --git a/src/conversion.F90 b/src/conversion.F90
new file mode 100644
index 0000000..d6e30e9
--- /dev/null
+++ b/src/conversion.F90
@@ -0,0 +1,185 @@
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+module conversion
+  implicit none
+  DECLARE_CCTK_PARAMETERS
+  private
+  
+  public make_cylindrical2cartesian
+  public make_spherical2cartesian
+  public make_cartesian2spherical
+  
+contains
+  
+  subroutine make_cylindrical2cartesian (xx, tt)
+    CCTK_REAL, intent(in)  :: xx(3)
+    CCTK_REAL, intent(out) :: tt(3,3)
+    CCTK_REAL, parameter :: eps = 1.0d-20
+    CCTK_REAL :: x,y,z
+    CCTK_REAL :: rho
+    
+    ! cylindrical coordinates: rho phi z
+    ! Cartesian   coordinates: x   y   z
+    !    rho^2    = x^2 + y^2
+    !    tan(phi) = y/x
+    
+    ! A[cart](p)_i = TT(p)_i^j A[cyl](p)_j
+    !    where p = p[cart]
+    
+    ! Transformation:
+    !    x = rho cos(phi)
+    !    y = rho sin(phi)
+    !    z = z
+    !    rho = sqrt(x^2+y^2)
+    !    phi = atan(y/x)
+    !    z   = z
+    
+    ! Jacobian:
+    !    drho/dx = x/rho
+    !    drho/dy = y/rho
+    !    drho/dz = 0
+    !    dphi/dx = -y/rho^2
+    !    dphi/dy = x/rho^2
+    !    dphi/dz = 0
+    !    dz  /dx = 0
+    !    dz  /dy = 0
+    !    dz  /dz = 1
+    
+    x = xx(1)
+    y = xx(2)
+    z = xx(3)
+    
+    rho = sqrt(x**2+y**2)
+    
+    tt(1,1) = x/(rho+eps)
+    tt(2,1) = y/(rho+eps)
+    tt(3,1) = 0
+    tt(1,2) = -y/(rho**2+eps)
+    tt(2,2) = x/(rho**2+eps)
+    tt(3,2) = 0
+    tt(1,3) = 0
+    tt(2,3) = 0
+    tt(3,3) = 1
+    
+  end subroutine make_cylindrical2cartesian
+  
+  subroutine make_spherical2cartesian (xx, tt)
+    CCTK_REAL, intent(in)  :: xx(3)
+    CCTK_REAL, intent(out) :: tt(3,3)
+    CCTK_REAL, parameter :: eps = 1.0d-20
+    CCTK_REAL :: x,y,z
+    CCTK_REAL :: rho,r
+    
+    ! Cartesian coordinates: x y z
+    ! spherical coordinates: r theta phi
+    !    r^2        = x^2 + y^2 + z^2
+    !    cos(theta) = z/r
+    !    tan(phi)   = y/x
+    
+    ! A[cart](p)_i = TT(p)_i^j A[spher](p)_j
+    !    where p = p[cart]
+    
+    ! Definition:
+    !    rho   = r sin(theta)
+    !    z/rho = cos(theta) / sin(theta)
+    
+    ! Transformation:
+    !    x = r sin(theta) cos(phi)
+    !    y = r sin(theta) sin(phi)
+    !    z = r cos(theta)
+    !    r     = sqrt(x^2+y^2+z^2)
+    !    theta = acos(z/r)
+    !    phi   = atan(y/x)
+    
+    ! Jacobian:
+    !    dr    /dx =  x/r
+    !    dr    /dy =  y/r
+    !    dr    /dz =  z/r
+    !    dtheta/dx =  xz/r^2rho           ! 0
+    !    dtheta/dy =  yz/r^2rho           ! 0
+    !    dtheta/dz = -rho/r^2             ! 1/rho
+    !    dphi  /dx = -y/rho^2
+    !    dphi  /dy =  x/rho^2
+    !    dphi  /dz =  0
+    
+    x = xx(1)
+    y = xx(2)
+    z = xx(3)
+    
+    rho = sqrt(x**2+y**2)
+    r = sqrt(x**2+y**2+z**2)
+    
+    tt(1,1) =  x/(r+eps)
+    tt(2,1) =  y/(r+eps)
+    tt(3,1) =  z/(r+eps)
+    tt(1,2) =  x*z/(r**2*rho+eps)
+    tt(2,2) =  y*z/(r**2*rho+eps)
+    tt(3,2) = -rho**2/(r**2*rho+eps)
+    tt(1,3) = -y/(rho**2+eps)
+    tt(2,3) =  x/(rho**2+eps)
+    tt(3,3) =  0
+    
+  end subroutine make_spherical2cartesian
+  
+  subroutine make_cartesian2spherical (xx, tt)
+    CCTK_REAL, intent(in)  :: xx(3)
+    CCTK_REAL, intent(out) :: tt(3,3)
+    CCTK_REAL, parameter :: eps = 1.0d-20
+    CCTK_REAL :: x,y,z
+    CCTK_REAL :: rho,r
+    
+    ! Cartesian coordinates: x y z
+    ! spherical coordinates: r theta phi
+    !    r^2        = x^2 + y^2 + z^2
+    !    cos(theta) = z/r
+    !    tan(phi)   = y/x
+    
+    ! A[spher](p)_i = TT(p)_i^j A[cart](p)_j
+    !    where p = p[cart]
+    
+    ! Definition:
+    !    rho   = r sin(theta)
+    !    z/rho = cos(theta) / sin(theta)
+    
+    ! Transformation:
+    !    x = r sin(theta) cos(phi)
+    !    y = r sin(theta) sin(phi)
+    !    z = r cos(theta)
+    !    r     = sqrt(x^2+y^2+z^2)
+    !    theta = acos(z/r)
+    !    phi   = atan(y/x)
+    
+    ! Jacobian:
+    !    dx/dr     =     sin(theta) cos(phi) =  x/r
+    !    dx/dtheta =   r cos(theta) cos(phi) =  x z/rho
+    !    dx/dphi   = - r sin(theta) sin(phi) = -y
+    !    dy/dr     =     sin(theta) sin(phi) =  y/r
+    !    dy/dtheta =   r cos(theta) sin(phi) =  y z/rho
+    !    dy/dphi   =   r sin(theta) cos(phi) =  x
+    !    dz/dr     =     cos(theta)          =  z/r
+    !    dz/dtheta = - r sin(theta)          = -rho
+    !    dz/dphi   =   0                     =  0
+    
+    x = xx(1)
+    y = xx(2)
+    z = xx(3)
+    
+    rho = sqrt(x**2+y**2)
+    r = sqrt(x**2+y**2+z**2)
+    
+    tt(1,1) =  x/(r+eps)
+    tt(2,1) =  x * z/(rho+eps)
+    tt(3,1) = -y
+    tt(1,2) =  y/(r+eps)
+    tt(2,2) =  y * z/(rho+eps)
+    tt(3,2) =  x
+    tt(1,3) =  z/(r+eps)
+    tt(2,3) = -rho
+    tt(3,3) =  0
+    
+  end subroutine make_cartesian2spherical
+  
+end module conversion
diff --git a/src/covderivs.F90 b/src/covderivs.F90
new file mode 100644
index 0000000..50e3624
--- /dev/null
+++ b/src/covderivs.F90
@@ -0,0 +1,201 @@
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+module covderivs
+  implicit none
+  DECLARE_CCTK_PARAMETERS
+  private
+  
+  public calc_scalargrad
+  public calc_covectorgrad
+  public calc_vectorgrad
+  public calc_tensorgrad
+  public calc_scalargradgrad
+  public calc_covectorgradgrad
+  public calc_longitudinal
+  public calc_gradlongitudinal
+  
+contains
+  
+  subroutine calc_scalargrad (f, df, gamma, gf)
+    CCTK_REAL, intent(in)  :: f
+    CCTK_REAL, intent(in)  :: df(3)
+    CCTK_REAL, intent(in)  :: gamma(3,3,3)
+    CCTK_REAL, intent(out) :: gf(3)
+    integer :: i
+    ! f;i = f,i
+    do i=1,3
+       gf(i) = df(i)
+    end do
+  end subroutine calc_scalargrad
+  
+  subroutine calc_covectorgrad (f, df, gamma, gf)
+    CCTK_REAL, intent(in)  :: f(3)
+    CCTK_REAL, intent(in)  :: df(3,3)
+    CCTK_REAL, intent(in)  :: gamma(3,3,3)
+    CCTK_REAL, intent(out) :: gf(3,3)
+    integer :: i,j,k
+    ! f_i;j = f_i,j - Gamma^k_ij f_k
+    do i=1,3
+       do j=1,3
+          gf(i,j) = df(i,j)
+          do k=1,3
+             gf(i,j) = gf(i,j) - gamma(k,i,j) * f(k)
+          end do
+       end do
+    end do
+  end subroutine calc_covectorgrad
+  
+  subroutine calc_vectorgrad (f, df, gamma, gf)
+    CCTK_REAL, intent(in)  :: f(3)
+    CCTK_REAL, intent(in)  :: df(3,3)
+    CCTK_REAL, intent(in)  :: gamma(3,3,3)
+    CCTK_REAL, intent(out) :: gf(3,3)
+    integer :: i,j,k
+    ! f^i;j = f^i,j + Gamma^i_kj f^k
+    do i=1,3
+       do j=1,3
+          gf(i,j) = df(i,j)
+          do k=1,3
+             gf(i,j) = gf(i,j) + gamma(i,k,j) * f(k)
+          end do
+       end do
+    end do
+  end subroutine calc_vectorgrad
+  
+  subroutine calc_tensorgrad (f, df, gamma, gf)
+    CCTK_REAL, intent(in)  :: f(3,3)
+    CCTK_REAL, intent(in)  :: df(3,3,3)
+    CCTK_REAL, intent(in)  :: gamma(3,3,3)
+    CCTK_REAL, intent(out) :: gf(3,3,3)
+    integer :: i,j,k,l
+    ! f_ij;k = f_ij,k - Gamma^l_ik f_lj - Gamma^l_jk f_il
+    do i=1,3
+       do j=1,3
+          do k=1,3
+             gf(i,j,k) = df(i,j,k)
+             do l=1,3
+                gf(i,j,k) = gf(i,j,k) - gamma(l,i,k) * f(l,j) &
+                     &                - gamma(l,j,k) * f(i,l)
+             end do
+          end do
+       end do
+    end do
+  end subroutine calc_tensorgrad
+  
+  
+  
+  subroutine calc_scalargradgrad (f, df, ddf, gamma, ggf)
+    CCTK_REAL, intent(in)  :: f
+    CCTK_REAL, intent(in)  :: df(3)
+    CCTK_REAL, intent(in)  :: ddf(3,3)
+    CCTK_REAL, intent(in)  :: gamma(3,3,3)
+    CCTK_REAL, intent(out) :: ggf(3,3)
+    integer :: i,j,k
+    ! f;ij = f,ij - Gamma^k_ij f,k
+    do i=1,3
+       do j=1,3
+          ggf(i,j) = ddf(i,j)
+          do k=1,3
+             ggf(i,j) = ggf(i,j) - gamma(k,i,j) * df(k)
+          end do
+       end do
+    end do
+  end subroutine calc_scalargradgrad
+  
+  subroutine calc_covectorgradgrad (f, df, ddf, gamma, dgamma, ggf)
+    CCTK_REAL, intent(in)  :: f(3)
+    CCTK_REAL, intent(in)  :: df(3,3)
+    CCTK_REAL, intent(in)  :: ddf(3,3,3)
+    CCTK_REAL, intent(in)  :: gamma(3,3,3), dgamma(3,3,3,3)
+    CCTK_REAL, intent(out) :: ggf(3,3,3)
+    CCTK_REAL :: gf(3,3)
+    integer :: i,j,k,l
+    ! f_i;j  = f_i,j - Gamma^l_ij f_l
+    ! f_i;jk = f_i,jk - Gamma^l_ij,k f_l - Gamma^l_ij f_l,k
+    !          - Gamma^l_ik f_l;j - Gamma^l_jk f_i;l
+    call calc_covectorgrad (f, df, gamma, gf)
+    do i=1,3
+       do j=1,3
+          do k=1,3
+             ggf(i,j,k) = ddf(i,j,k)
+             do l=1,3
+                ggf(i,j,k) = ggf(i,j,k) &
+                     - dgamma(l,i,j,k) * f(l) - gamma(l,i,j) * df(l,k) &
+                     - gamma(l,i,k) * gf(l,j) - gamma(l,j,k) * gf(i,l)
+             end do
+          end do
+       end do
+    end do
+  end subroutine calc_covectorgradgrad
+  
+  subroutine calc_vectorgradgrad (f, df, ddf, gamma, dgamma, ggf)
+    CCTK_REAL, intent(in)  :: f(3)
+    CCTK_REAL, intent(in)  :: df(3,3)
+    CCTK_REAL, intent(in)  :: ddf(3,3,3)
+    CCTK_REAL, intent(in)  :: gamma(3,3,3), dgamma(3,3,3,3)
+    CCTK_REAL, intent(out) :: ggf(3,3,3)
+    CCTK_REAL :: gf(3,3)
+    integer :: i,j,k,l
+    ! f^i;j  = f^i,j + Gamma^i_lj f^l
+    ! f^i;jk = f^i,jk + Gamma^i_lj,k f^l + Gamma^i_lj f^l,k
+    !          - Gamma^i_lk f^l;j - Gamma^l_jk f^i;l
+    call calc_vectorgrad (f, df, gamma, gf)
+    do i=1,3
+       do j=1,3
+          do k=1,3
+             ggf(i,j,k) = ddf(i,j,k)
+             do l=1,3
+                ggf(i,j,k) = ggf(i,j,k) &
+                     + dgamma(i,l,j,k) * f(l) + gamma(i,l,j) * df(l,k) &
+                     - gamma(i,l,k) * gf(l,j) - gamma(l,j,k) * gf(i,l)
+             end do
+          end do
+       end do
+    end do
+  end subroutine calc_vectorgradgrad
+  
+  
+  
+  subroutine calc_longitudinal (gf, gg, gu, lf)
+    CCTK_REAL, intent(in)  :: gf(3,3)
+    CCTK_REAL, intent(in)  :: gg(3,3), gu(3,3)
+    CCTK_REAL, intent(out) :: lf(3,3)
+    integer :: i,j,k,l
+    ! (Lf)_ij = D_i f_j + D_j f_i - 2/3 g_ij D^l f_l
+    do i=1,3
+       do j=1,3
+          lf(i,j) = gf(i,j) + gf(j,i)
+          do k=1,3
+             do l=1,3
+                lf(i,j) = lf(i,j) - (2.d0/3.d0) * gg(i,j) * gu(k,l) * gf(k,l)
+             end do
+          end do
+       end do
+    end do
+  end subroutine calc_longitudinal
+  
+  subroutine calc_gradlongitudinal (ggf, gg, gu, glf)
+    CCTK_REAL, intent(in)  :: ggf(3,3,3)
+    CCTK_REAL, intent(in)  :: gg(3,3), gu(3,3)
+    CCTK_REAL, intent(out) :: glf(3,3,3)
+    integer :: i,j,k,l,m
+    ! D_k (Lf)_ij = D_k (D_i f_j + D_j f_i - 2/3 g_ij D^l f_l)
+    do i=1,3
+       do j=1,3
+          do k=1,3
+             glf(i,j,k) = ggf(j,i,k) + ggf(i,j,k)
+             do l=1,3
+                do m=1,3
+                   glf(i,j,k) = glf(i,j,k) &
+                        - (2.d0/3.d0) * gg(i,j) * gu(l,m) * ggf(l,m,k)
+                end do
+             end do
+          end do
+       end do
+    end do
+  end subroutine calc_gradlongitudinal
+  
+end module covderivs
diff --git a/src/covderivs2.F90 b/src/covderivs2.F90
new file mode 100644
index 0000000..15e0143
--- /dev/null
+++ b/src/covderivs2.F90
@@ -0,0 +1,204 @@
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+module covderivs2
+  implicit none
+  DECLARE_CCTK_PARAMETERS
+  private
+  
+  public calc_2scalargrad
+  public calc_2covectorgrad
+  public calc_2vectorgrad
+  public calc_2tensorgrad
+  
+  public calc_2scalargradgrad
+  public calc_2covectorgradgrad
+  public calc_2vectorgradgrad
+  
+  public calc_2longitudinal
+  public calc_2gradlongitudinal
+  
+contains
+  
+  subroutine calc_2scalargrad (f, df, gamma, gf)
+    CCTK_REAL, intent(in)  :: f
+    CCTK_REAL, intent(in)  :: df(2)
+    CCTK_REAL, intent(in)  :: gamma(2,2,2)
+    CCTK_REAL, intent(out) :: gf(2)
+    integer :: i
+    ! f;i = f,i
+    do i=1,2
+       gf(i) = df(i)
+    end do
+  end subroutine calc_2scalargrad
+  
+  subroutine calc_2covectorgrad (f, df, gamma, gf)
+    CCTK_REAL, intent(in)  :: f(2)
+    CCTK_REAL, intent(in)  :: df(2,2)
+    CCTK_REAL, intent(in)  :: gamma(2,2,2)
+    CCTK_REAL, intent(out) :: gf(2,2)
+    integer :: i,j,k
+    ! f_i;j = f_i,j - Gamma^k_ij f_k
+    do i=1,2
+       do j=1,2
+          gf(i,j) = df(i,j)
+          do k=1,2
+             gf(i,j) = gf(i,j) - gamma(k,i,j) * f(k)
+          end do
+       end do
+    end do
+  end subroutine calc_2covectorgrad
+  
+  subroutine calc_2vectorgrad (f, df, gamma, gf)
+    CCTK_REAL, intent(in)  :: f(2)
+    CCTK_REAL, intent(in)  :: df(2,2)
+    CCTK_REAL, intent(in)  :: gamma(2,2,2)
+    CCTK_REAL, intent(out) :: gf(2,2)
+    integer :: i,j,k
+    ! f^i;j = f^i,j + Gamma^i_kj f^k
+    do i=1,2
+       do j=1,2
+          gf(i,j) = df(i,j)
+          do k=1,2
+             gf(i,j) = gf(i,j) + gamma(i,k,j) * f(k)
+          end do
+       end do
+    end do
+  end subroutine calc_2vectorgrad
+  
+  subroutine calc_2tensorgrad (f, df, gamma, gf)
+    CCTK_REAL, intent(in)  :: f(2,2)
+    CCTK_REAL, intent(in)  :: df(2,2,2)
+    CCTK_REAL, intent(in)  :: gamma(2,2,2)
+    CCTK_REAL, intent(out) :: gf(2,2,2)
+    integer :: i,j,k,l
+    ! f_ij;k = f_ij,k - Gamma^l_ik f_lj - Gamma^l_jk f_il
+    do i=1,2
+       do j=1,2
+          do k=1,2
+             gf(i,j,k) = df(i,j,k)
+             do l=1,2
+                gf(i,j,k) = gf(i,j,k) - gamma(l,i,k) * f(l,j) &
+                     &                - gamma(l,j,k) * f(i,l)
+             end do
+          end do
+       end do
+    end do
+  end subroutine calc_2tensorgrad
+  
+  
+  
+  subroutine calc_2scalargradgrad (f, df, ddf, gamma, ggf)
+    CCTK_REAL, intent(in)  :: f
+    CCTK_REAL, intent(in)  :: df(2)
+    CCTK_REAL, intent(in)  :: ddf(2,2)
+    CCTK_REAL, intent(in)  :: gamma(2,2,2)
+    CCTK_REAL, intent(out) :: ggf(2,2)
+    integer :: i,j,k
+    ! f;ij = f,ij - Gamma^k_ij f,k
+    do i=1,2
+       do j=1,2
+          ggf(i,j) = ddf(i,j)
+          do k=1,2
+             ggf(i,j) = ggf(i,j) - gamma(k,i,j) * df(k)
+          end do
+       end do
+    end do
+  end subroutine calc_2scalargradgrad
+  
+  subroutine calc_2covectorgradgrad (f, df, ddf, gamma, dgamma, ggf)
+    CCTK_REAL, intent(in)  :: f(2)
+    CCTK_REAL, intent(in)  :: df(2,2)
+    CCTK_REAL, intent(in)  :: ddf(2,2,2)
+    CCTK_REAL, intent(in)  :: gamma(2,2,2), dgamma(2,2,2,2)
+    CCTK_REAL, intent(out) :: ggf(2,2,2)
+    CCTK_REAL :: gf(2,2)
+    integer :: i,j,k,l
+    ! f_i;j  = f_i,j - Gamma^l_ij f_l
+    ! f_i;jk = f_i,jk - Gamma^l_ij,k f_l - Gamma^l_ij f_l,k
+    !          - Gamma^l_ik f_l;j - Gamma^l_jk f_i;l
+    call calc_2covectorgrad (f, df, gamma, gf)
+    do i=1,2
+       do j=1,2
+          do k=1,2
+             ggf(i,j,k) = ddf(i,j,k)
+             do l=1,2
+                ggf(i,j,k) = ggf(i,j,k) &
+                     - dgamma(l,i,j,k) * f(l) - gamma(l,i,j) * df(l,k) &
+                     - gamma(l,i,k) * gf(l,j) - gamma(l,j,k) * gf(i,l)
+             end do
+          end do
+       end do
+    end do
+  end subroutine calc_2covectorgradgrad
+  
+  subroutine calc_2vectorgradgrad (f, df, ddf, gamma, dgamma, ggf)
+    CCTK_REAL, intent(in)  :: f(2)
+    CCTK_REAL, intent(in)  :: df(2,2)
+    CCTK_REAL, intent(in)  :: ddf(2,2,2)
+    CCTK_REAL, intent(in)  :: gamma(2,2,2), dgamma(2,2,2,2)
+    CCTK_REAL, intent(out) :: ggf(2,2,2)
+    CCTK_REAL :: gf(2,2)
+    integer :: i,j,k,l
+    ! f^i;j  = f^i,j + Gamma^i_lj f^l
+    ! f^i;jk = f^i,jk + Gamma^i_lj,k f^l + Gamma^i_lj f^l,k
+    !          - Gamma^i_lk f^l;j - Gamma^l_jk f^i;l
+    call calc_2vectorgrad (f, df, gamma, gf)
+    do i=1,2
+       do j=1,2
+          do k=1,2
+             ggf(i,j,k) = ddf(i,j,k)
+             do l=1,2
+                ggf(i,j,k) = ggf(i,j,k) &
+                     + dgamma(i,l,j,k) * f(l) + gamma(i,l,j) * df(l,k) &
+                     - gamma(i,l,k) * gf(l,j) - gamma(l,j,k) * gf(i,l)
+             end do
+          end do
+       end do
+    end do
+  end subroutine calc_2vectorgradgrad
+  
+  
+  
+  subroutine calc_2longitudinal (gf, gg, gu, lf)
+    CCTK_REAL, intent(in)  :: gf(2,2)
+    CCTK_REAL, intent(in)  :: gg(2,2), gu(2,2)
+    CCTK_REAL, intent(out) :: lf(2,2)
+    integer :: i,j,k,l
+    ! (Lf)_ij = D_i f_j + D_j f_i - 2/3 g_ij D^l f_l
+    do i=1,2
+       do j=1,2
+          lf(i,j) = gf(i,j) + gf(j,i)
+          do k=1,2
+             do l=1,2
+                lf(i,j) = lf(i,j) - (2.d0/3.d0) * gg(i,j) * gu(k,l) * gf(k,l)
+             end do
+          end do
+       end do
+    end do
+  end subroutine calc_2longitudinal
+  
+  subroutine calc_2gradlongitudinal (ggf, gg, gu, glf)
+    CCTK_REAL, intent(in)  :: ggf(2,2,2)
+    CCTK_REAL, intent(in)  :: gg(2,2), gu(2,2)
+    CCTK_REAL, intent(out) :: glf(2,2,2)
+    integer :: i,j,k,l,m
+    ! D_k (Lf)_ij = D_k (D_i f_j + D_j f_i - 2/3 g_ij D^l f_l)
+    do i=1,2
+       do j=1,2
+          do k=1,2
+             glf(i,j,k) = ggf(j,i,k) + ggf(i,j,k)
+             do l=1,2
+                do m=1,2
+                   glf(i,j,k) = glf(i,j,k) &
+                        - (2.d0/3.d0) * gg(i,j) * gu(l,m) * ggf(l,m,k)
+                end do
+             end do
+          end do
+       end do
+    end do
+  end subroutine calc_2gradlongitudinal
+  
+end module covderivs2
diff --git a/src/depend.pl b/src/depend.pl
new file mode 100644
index 0000000..b41b8f6
--- /dev/null
+++ b/src/depend.pl
@@ -0,0 +1,46 @@
+#!/usr/bin/perl -w
+# $Header$
+
+# Create dependencies for Fortran 90 "use" and "include" statements
+
+$src = $ARGV[0];
+$srcfile = $ARGV[1];
+$dest = $ARGV[2];
+$srcdir = $ARGV[3];
+$moddir = $ARGV[4];
+$srcsuffix = $ARGV[5];
+$modsuffix = $ARGV[6];
+@otherdirs = @ARGV[7..$#ARGV];
+
+print "# $src\n";
+print "$dest:";
+
+open (FILE, $src);
+while (<FILE>) {
+  if (/^\s*include\s*['"]([^'"]+)['"]/i) {
+    print " $srcdir$1";
+  } elsif (/^\s*use\s+(\w+)/i) {
+    $found = 0;
+    if (! $found) {
+      if (-e "$srcdir$1$srcsuffix") {
+	$found = 1;
+	print " $moddir$1$modsuffix";
+      }
+    }
+    if (! $found) {
+      foreach $dir (@otherdirs) {
+	if (-e "$dir$1$modsuffix") {
+	  $found = 1;
+	  print " $dir$1$modsuffix";
+	  last;
+	}
+      }
+    }
+    if (! $found) {
+      die "\nWhile tracing depencencies:\nFortran module $1 (referenced in file $srcfile) not found.\nAborting.\n\n";
+    }
+  }
+}
+close (FILE);
+
+print "\n";
diff --git a/src/derivs.F90 b/src/derivs.F90
new file mode 100644
index 0000000..1d425c8
--- /dev/null
+++ b/src/derivs.F90
@@ -0,0 +1,46 @@
+! $Header$
+
+#ifndef TGR_INCLUDED
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+module derivs
+  implicit none
+  private
+  public get_derivs
+  public get_derivs2
+contains
+#endif
+  subroutine get_derivs (a, f, pos, off, dx)
+    CCTK_REAL, intent(in)  :: a(*)
+    CCTK_REAL, intent(out) :: f(3)
+    integer,   intent(in)  :: pos, off(3)
+    CCTK_REAL, intent(in)  :: dx(3)
+    integer :: i
+    do i=1,3
+       f(i) = (a(pos+off(i)) - a(pos-off(i))) / (2*dx(i))
+    end do
+  end subroutine get_derivs
+  subroutine get_derivs2 (a, f, pos, off, dx)
+    CCTK_REAL, intent(in)  :: a(*)
+    CCTK_REAL, intent(out) :: f(3,3)
+    integer,   intent(in)  :: pos, off(3)
+    CCTK_REAL, intent(in)  :: dx(3)
+    integer :: i
+    do i=1,3
+       f(i,i) = (a(pos+off(i)) - 2*a(pos) + a(pos-off(i))) / dx(i)**2
+    end do
+    f(1,2) = (  a(pos-off(1)-off(2)) - a(pos+off(1)-off(2)) &
+         &    - a(pos-off(1)+off(2)) + a(pos+off(1)+off(2))) / (4*dx(1)*dx(2))
+    f(2,1) = f(1,2)
+    f(1,3) = (  a(pos-off(1)-off(3)) - a(pos+off(1)-off(3)) &
+         &    - a(pos-off(1)+off(3)) + a(pos+off(1)+off(3))) / (4*dx(1)*dx(3))
+    f(3,1) = f(1,3)
+    f(2,3) = (  a(pos-off(2)-off(3)) - a(pos+off(2)-off(3)) &
+         &    - a(pos-off(2)+off(3)) + a(pos+off(2)+off(3))) / (4*dx(2)*dx(3))
+    f(3,2) = f(2,3)
+  end subroutine get_derivs2
+#ifndef TGR_INCLUDED
+end module derivs
+#endif
diff --git a/src/derivs2.F90 b/src/derivs2.F90
new file mode 100644
index 0000000..4ca51bf
--- /dev/null
+++ b/src/derivs2.F90
@@ -0,0 +1,40 @@
+! $Header$
+
+#ifndef TGR_INCLUDED
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+module derivs2
+  implicit none
+  private
+  public get_2derivs
+  public get_2derivs2
+contains
+#endif
+  subroutine get_2derivs (a, f, pos, off, dx)
+    CCTK_REAL, intent(in)  :: a(*)
+    CCTK_REAL, intent(out) :: f(2)
+    integer,   intent(in)  :: pos, off(2)
+    CCTK_REAL, intent(in)  :: dx(2)
+    integer :: i
+    do i=1,2
+       f(i) = (a(pos+off(i)) - a(pos-off(i))) / (2*dx(i))
+    end do
+  end subroutine get_2derivs
+  subroutine get_2derivs2 (a, f, pos, off, dx)
+    CCTK_REAL, intent(in)  :: a(*)
+    CCTK_REAL, intent(out) :: f(2,2)
+    integer,   intent(in)  :: pos, off(2)
+    CCTK_REAL, intent(in)  :: dx(2)
+    integer :: i
+    do i=1,2
+       f(i,i) = (a(pos+off(i)) - 2*a(pos) + a(pos-off(i))) / dx(i)**2
+    end do
+    f(1,2) = (  a(pos-off(1)-off(2)) - a(pos+off(1)-off(2)) &
+         &    - a(pos-off(1)+off(2)) + a(pos+off(1)+off(2))) / (4*dx(1)*dx(2))
+    f(2,1) = f(1,2)
+  end subroutine get_2derivs2
+#ifndef TGR_INCLUDED
+end module derivs2
+#endif
diff --git a/src/hyperslab.F90 b/src/hyperslab.F90
new file mode 100644
index 0000000..e540cb2
--- /dev/null
+++ b/src/hyperslab.F90
@@ -0,0 +1,33 @@
+! $Header$
+
+#include "cctk.h"
+
+module hyperslab
+  implicit none
+  public
+  
+  interface
+     
+     subroutine hyperslab_fillhyperslab (ierr, cctkgh, target_proc, &
+          vindex, vtimelevel, &
+          hdim, global_startpoint, directions, lengths, downsample, &
+          nhdata, hdata, hsize)
+       implicit none
+       integer      ierr
+       CCTK_POINTER cctkgh
+       integer      target_proc
+       integer      vindex
+       integer      vtimelevel
+       integer      hdim
+       integer      global_startpoint(*) ! vdim
+       integer      directions(*) ! vdim
+       integer      lengths(hdim)
+       integer      downsample(hdim)
+       integer      nhdata
+       CCTK_REAL    hdata(*)
+       integer      hsize(hdim)
+     end subroutine hyperslab_fillhyperslab
+     
+  end interface
+  
+end module hyperslab
diff --git a/src/lapack.F90 b/src/lapack.F90
new file mode 100644
index 0000000..5bcd70b
--- /dev/null
+++ b/src/lapack.F90
@@ -0,0 +1,51 @@
+! $Header$
+
+#include "cctk.h"
+
+module lapack
+  implicit none
+  public
+  
+  interface geev
+     subroutine sgeev (jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, &
+          &            ldvr, work, lwork, info)
+       character jobvl, jobvr
+       integer   info, lda, ldvl, ldvr, lwork, n
+       real      a(lda,n), vl(ldvl,n), vr(ldvr,n), &
+            &    wi(n), work(lwork), wr(n)
+     end subroutine sgeev
+     subroutine dgeev (jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, &
+          &            ldvr, work, lwork, info)
+       character        jobvl, jobvr
+       integer          info, lda, ldvl, ldvr, lwork, n
+       double precision a(lda,n), vl(ldvl,n), vr(ldvr,n), &
+            &           wi(n), work(lwork), wr(n)
+     end subroutine dgeev
+  end interface
+
+  interface gesv
+     subroutine sgesv (n, nrhs, a, lda, ipiv, b, ldb, info)
+       implicit none
+       integer n
+       integer nrhs
+       integer lda
+       real    a(lda,n)
+       integer ipiv(n)
+       integer ldb
+       real    b(ldb,nrhs)
+       integer info
+     end subroutine sgesv
+     subroutine dgesv (n, nrhs, a, lda, ipiv, b, ldb, info)
+       implicit none
+       integer          n
+       integer          nrhs
+       integer          lda
+       double precision a(lda,n)
+       integer          ipiv(n)
+       integer          ldb
+       double precision b(ldb,nrhs)
+       integer          info
+     end subroutine dgesv
+  end interface
+  
+end module lapack
diff --git a/src/make.code.defn b/src/make.code.defn
new file mode 100644
index 0000000..4b157e4
--- /dev/null
+++ b/src/make.code.defn
@@ -0,0 +1,25 @@
+# Main make.code.defn file for thorn TGRtensor   -*-Makefile-*-
+# $Header$
+
+# Source files in this directory
+SRCS =	cactus.F90				\
+	constants.F90				\
+	conversion.F90				\
+	covderivs.F90				\
+	covderivs2.F90				\
+	derivs.F90				\
+	derivs2.F90				\
+	hyperslab.F90				\
+	hyperslab.c				\
+	lapack.F90				\
+	matexp.F90				\
+	matinv.F90				\
+	rotation.F90				\
+	tensor.F90				\
+	tensor2.F90				\
+	pointwise.F90				\
+	pointwise2.F90
+
+# Subdirectories containing source files
+SUBDIRS = 
+
diff --git a/src/make.code.deps b/src/make.code.deps
new file mode 100644
index 0000000..e8ea1b9
--- /dev/null
+++ b/src/make.code.deps
@@ -0,0 +1,11 @@
+# Main make.code.deps file for thorn TGRtensor   -*-Makefile-*-
+# $Header$
+
+USESTHORNS = 
+
+# Automatically create dependencies for Fortran modules and Fortran includes
+define F90_DEPENDENCIES
+$(F_DEPEND) $(INC_DIRS:%=-I%) $(EXTRA_DEFINES:%=-D%) -DFCODE $< $(F_DEPEND_OUT)
+$(DEPENDENCY_FIXER)
+dir=`pwd`; $(CPP) $(INC_DIRS:%=-I%) $(EXTRA_DEFINES:%=-D%) -DFCODE $< | $(PERL) $(SRCDIR)/depend.pl - $< $(basename $(notdir $<)).F90.o $(SRCDIR)/ ./ .F90 .F90.o $(USESTHORNS:%=$$dir/../%/) >> $@ || { rm $@; exit 1; }
+endef
diff --git a/src/make.configuration.deps b/src/make.configuration.deps
new file mode 100644
index 0000000..6a85f9c
--- /dev/null
+++ b/src/make.configuration.deps
@@ -0,0 +1,6 @@
+# Main make.configuration.deps file for thorn TGRtensor   -*-Makefile-*-
+# $Header$
+
+USESTHORNS = 
+
+$(CCTK_LIBDIR)$(DIRSEP)libTGRtensor.a: $(USESTHORNS:%=$(CCTK_LIBDIR)$(DIRSEP)lib%.a)
diff --git a/src/matexp.F90 b/src/matexp.F90
new file mode 100644
index 0000000..d2ecc77
--- /dev/null
+++ b/src/matexp.F90
@@ -0,0 +1,38 @@
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+module matexp
+  use constants
+  implicit none
+  DECLARE_CCTK_PARAMETERS
+  private
+  
+  public calc_exp3
+  
+contains
+  
+  subroutine calc_exp3 (h3, g3)
+    CCTK_REAL, intent(in)  :: h3(3,3)
+    CCTK_REAL, intent(out) :: g3(3,3)
+    CCTK_REAL :: nfact
+    CCTK_REAL :: tmp(3,3)
+    integer   :: n
+    integer   :: i, j
+    
+    g3 = delta3
+    
+    nfact = 1
+    tmp = delta3
+    do n = 1, 18
+       nfact = nfact * n
+       tmp = matmul (h3, tmp)
+       
+       ! exp(x) = sum_n x^n / n!
+       g3 = g3 + tmp / nfact
+    end do
+    
+  end subroutine calc_exp3
+  
+end module matexp
diff --git a/src/matinv.F90 b/src/matinv.F90
new file mode 100644
index 0000000..ab6c515
--- /dev/null
+++ b/src/matinv.F90
@@ -0,0 +1,54 @@
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+module matinv
+  use constants
+  use lapack
+  implicit none
+  DECLARE_CCTK_PARAMETERS
+  private
+  
+  public calc_inv3
+  public calc_inv4
+  
+contains
+  
+  subroutine calc_inv3 (g3, gu3)
+    CCTK_REAL, intent(in)  :: g3(3,3)
+    CCTK_REAL, intent(out) :: gu3(3,3)
+    CCTK_REAL :: tmp(3,3)
+    integer   :: ipiv(3)
+    integer   :: info
+    character :: msg*1000
+    
+    tmp = g3
+    gu3 = delta3
+    call gesv (3, 3, tmp, 3, ipiv, gu3, 3, info)
+    
+    if (info /= 0) then
+       write (msg, '("Error in call to GESV, info=",i2)') info
+       call CCTK_WARN (0, trim(msg))
+    end if
+  end subroutine calc_inv3
+  
+  subroutine calc_inv4 (g4, gu4)
+    CCTK_REAL, intent(in)  :: g4(0:3,0:3)
+    CCTK_REAL, intent(out) :: gu4(0:3,0:3)
+    CCTK_REAL :: tmp(0:3,0:3)
+    integer   :: ipiv(0:3)
+    integer   :: info
+    character :: msg*1000
+    
+    tmp = g4
+    gu4 = delta4
+    call gesv (4, 4, tmp, 4, ipiv, gu4, 4, info)
+    
+    if (info /= 0) then
+       write (msg, '("Error in call to GESV, info=",i2)') info
+       call CCTK_WARN (0, trim(msg))
+    end if
+  end subroutine calc_inv4
+  
+end module matinv
diff --git a/src/pointwise.F90 b/src/pointwise.F90
new file mode 100644
index 0000000..827ad6e
--- /dev/null
+++ b/src/pointwise.F90
@@ -0,0 +1,252 @@
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+module pointwise
+!!$  use derivs
+  implicit none
+  DECLARE_CCTK_PARAMETERS
+  private
+  public calc_position
+  public calc_offsets
+  public get_scalar
+  public get_vector
+  public get_tensor
+  public get_connections
+  public set_scalar
+  public set_vector
+  public set_tensor
+  public set_connections
+  public get_scalarderivs
+  public get_vectorderivs
+  public get_tensorderivs
+  public get_scalarderivs2
+  public get_vectorderivs2
+  public get_tensorderivs2
+contains
+#define TGR_INCLUDED
+#include "derivs.F90"
+#undef TGR_INCLUDED
+  subroutine calc_position (shape, i,j,k, pos)
+    integer, intent(in)  :: shape(3)
+    integer, intent(in)  :: i,j,k
+    integer, intent(out) :: pos
+    pos = 1 + i-1 + shape(1) * (j-1 + shape(2) * (k-1))
+  end subroutine calc_position
+  subroutine calc_offsets (shape, off)
+    integer, intent(in)  :: shape(3)
+    integer, intent(out) :: off(3)
+    off(1) = 1
+    off(2) = shape(1)
+    off(3) = shape(1) * shape(2)
+  end subroutine calc_offsets
+  subroutine get_scalar (a, f, pos)
+    CCTK_REAL, intent(in)  :: a(*)
+    CCTK_REAL, intent(out) :: f
+    integer,   intent(in)  :: pos
+    f = a(pos)
+  end subroutine get_scalar
+  subroutine get_vector (ax,ay,az, f, pos)
+    CCTK_REAL, intent(in)  :: ax(*),ay(*),az(*)
+    CCTK_REAL, intent(out) :: f(3)
+    integer,   intent(in)  :: pos
+    f(1) = ax(pos)
+    f(2) = ay(pos)
+    f(3) = az(pos)
+  end subroutine get_vector
+  subroutine get_tensor (axx,axy,axz,ayy,ayz,azz, f, pos)
+    CCTK_REAL, intent(in)  :: axx(*),axy(*),axz(*),ayy(*),ayz(*),azz(*)
+    CCTK_REAL, intent(out) :: f(3,3)
+    integer,   intent(in)  :: pos
+    f(1,1) = axx(pos)
+    f(1,2) = axy(pos)
+    f(1,3) = axz(pos)
+    f(2,1) = axy(pos)
+    f(2,2) = ayy(pos)
+    f(2,3) = ayz(pos)
+    f(3,1) = axz(pos)
+    f(3,2) = ayz(pos)
+    f(3,3) = azz(pos)
+  end subroutine get_tensor
+  subroutine get_connections (&
+       gammaxxx,gammaxxy,gammaxxz,gammaxyy,gammaxyz,gammaxzz, &
+       gammayxx,gammayxy,gammayxz,gammayyy,gammayyz,gammayzz, &
+       gammazxx,gammazxy,gammazxz,gammazyy,gammazyz,gammazzz, gamma, pos)
+    CCTK_REAL, intent(in)  :: gammaxxx(*),gammaxxy(*),gammaxxz(*)
+    CCTK_REAL, intent(in)  :: gammaxyy(*),gammaxyz(*),gammaxzz(*)
+    CCTK_REAL, intent(in)  :: gammayxx(*),gammayxy(*),gammayxz(*)
+    CCTK_REAL, intent(in)  :: gammayyy(*),gammayyz(*),gammayzz(*)
+    CCTK_REAL, intent(in)  :: gammazxx(*),gammazxy(*),gammazxz(*)
+    CCTK_REAL, intent(in)  :: gammazyy(*),gammazyz(*),gammazzz(*)
+    CCTK_REAL, intent(out) :: gamma(3,3,3)
+    integer,   intent(in)  :: pos
+    gamma(1,1,1) = gammaxxx(pos)
+    gamma(1,1,2) = gammaxxy(pos)
+    gamma(1,1,3) = gammaxxz(pos)
+    gamma(1,2,1) = gammaxxy(pos)
+    gamma(1,2,2) = gammaxyy(pos)
+    gamma(1,2,3) = gammaxyz(pos)
+    gamma(1,3,1) = gammaxxz(pos)
+    gamma(1,3,2) = gammaxyz(pos)
+    gamma(1,3,3) = gammaxzz(pos)
+    gamma(2,1,1) = gammayxx(pos)
+    gamma(2,1,2) = gammayxy(pos)
+    gamma(2,1,3) = gammayxz(pos)
+    gamma(2,2,1) = gammayxy(pos)
+    gamma(2,2,2) = gammayyy(pos)
+    gamma(2,2,3) = gammayyz(pos)
+    gamma(2,3,1) = gammayxz(pos)
+    gamma(2,3,2) = gammayyz(pos)
+    gamma(2,3,3) = gammayzz(pos)
+    gamma(3,1,1) = gammazxx(pos)
+    gamma(3,1,2) = gammazxy(pos)
+    gamma(3,1,3) = gammazxz(pos)
+    gamma(3,2,1) = gammazxy(pos)
+    gamma(3,2,2) = gammazyy(pos)
+    gamma(3,2,3) = gammazyz(pos)
+    gamma(3,3,1) = gammazxz(pos)
+    gamma(3,3,2) = gammazyz(pos)
+    gamma(3,3,3) = gammazzz(pos)
+  end subroutine get_connections
+  subroutine set_scalar (f, a, pos)
+    CCTK_REAL, intent(in) :: f
+    CCTK_REAL             :: a(*)
+    integer,   intent(in) :: pos
+    a(pos) = f
+  end subroutine set_scalar
+  subroutine set_vector (f, ax,ay,az, pos)
+    CCTK_REAL, intent(in) :: f(3)
+    CCTK_REAL             :: ax(*),ay(*),az(*)
+    integer,   intent(in) :: pos
+    ax(pos) = f(1)
+    ay(pos) = f(2)
+    az(pos) = f(3)
+  end subroutine set_vector
+  subroutine set_tensor (f, axx,axy,axz,ayy,ayz,azz, pos)
+    CCTK_REAL, intent(in) :: f(3,3)
+    CCTK_REAL             :: axx(*),axy(*),axz(*),ayy(*),ayz(*),azz(*)
+    integer,   intent(in) :: pos
+    axx(pos) = f(1,1)
+    axy(pos) = f(1,2)
+    axz(pos) = f(1,3)
+    ayy(pos) = f(2,2)
+    ayz(pos) = f(2,3)
+    azz(pos) = f(3,3)
+  end subroutine set_tensor
+  subroutine set_connections (gamma, &
+       gammaxxx,gammaxxy,gammaxxz,gammaxyy,gammaxyz,gammaxzz, &
+       gammayxx,gammayxy,gammayxz,gammayyy,gammayyz,gammayzz, &
+       gammazxx,gammazxy,gammazxz,gammazyy,gammazyz,gammazzz, pos)
+    CCTK_REAL, intent(in) :: gamma(3,3,3)
+    CCTK_REAL             :: gammaxxx(*),gammaxxy(*),gammaxxz(*)
+    CCTK_REAL             :: gammaxyy(*),gammaxyz(*),gammaxzz(*)
+    CCTK_REAL             :: gammayxx(*),gammayxy(*),gammayxz(*)
+    CCTK_REAL             :: gammayyy(*),gammayyz(*),gammayzz(*)
+    CCTK_REAL             :: gammazxx(*),gammazxy(*),gammazxz(*)
+    CCTK_REAL             :: gammazyy(*),gammazyz(*),gammazzz(*)
+    integer,   intent(in)  :: pos
+    gammaxxx(pos) = gamma(1,1,1)
+    gammaxxy(pos) = gamma(1,1,2)
+    gammaxxz(pos) = gamma(1,1,3)
+    gammaxyy(pos) = gamma(1,2,2)
+    gammaxyz(pos) = gamma(1,2,3)
+    gammaxzz(pos) = gamma(1,3,3)
+    gammayxx(pos) = gamma(2,1,1)
+    gammayxy(pos) = gamma(2,1,2)
+    gammayxz(pos) = gamma(2,1,3)
+    gammayyy(pos) = gamma(2,2,2)
+    gammayyz(pos) = gamma(2,2,3)
+    gammayzz(pos) = gamma(2,3,3)
+    gammazxx(pos) = gamma(3,1,1)
+    gammazxy(pos) = gamma(3,1,2)
+    gammazxz(pos) = gamma(3,1,3)
+    gammazyy(pos) = gamma(3,2,2)
+    gammazyz(pos) = gamma(3,2,3)
+    gammazzz(pos) = gamma(3,3,3)
+  end subroutine set_connections
+  subroutine get_scalarderivs (a, f, pos, off, dx)
+    CCTK_REAL, intent(in)  :: a(*)
+    CCTK_REAL, intent(out) :: f(3)
+    integer,   intent(in)  :: pos, off(3)
+    CCTK_REAL, intent(in)  :: dx(3)
+    call get_derivs (a, f, pos, off, dx)
+  end subroutine get_scalarderivs
+  subroutine get_vectorderivs (ax,ay,az, f, pos, off, dx)
+    CCTK_REAL, intent(in)  :: ax(*),ay(*),az(*)
+    CCTK_REAL, intent(out) :: f(3,3)
+    integer,   intent(in)  :: pos, off(3)
+    CCTK_REAL, intent(in)  :: dx(3)
+    CCTK_REAL :: fx(3),fy(3),fz(3)
+    call get_derivs (ax, fx, pos, off, dx)
+    call get_derivs (ay, fy, pos, off, dx)
+    call get_derivs (az, fz, pos, off, dx)
+    f(1,:) = fx
+    f(2,:) = fy
+    f(3,:) = fz
+  end subroutine get_vectorderivs
+  subroutine get_tensorderivs (axx,axy,axz,ayy,ayz,azz, f, pos, off, dx)
+    CCTK_REAL, intent(in)  :: axx(*),axy(*),axz(*),ayy(*),ayz(*),azz(*)
+    CCTK_REAL, intent(out) :: f(3,3,3)
+    integer,   intent(in)  :: pos, off(3)
+    CCTK_REAL, intent(in)  :: dx(3)
+    CCTK_REAL :: fxx(3),fxy(3),fxz(3),fyy(3),fyz(3),fzz(3)
+    call get_derivs (axx, fxx, pos, off, dx)
+    call get_derivs (axy, fxy, pos, off, dx)
+    call get_derivs (axz, fxz, pos, off, dx)
+    call get_derivs (ayy, fyy, pos, off, dx)
+    call get_derivs (ayz, fyz, pos, off, dx)
+    call get_derivs (azz, fzz, pos, off, dx)
+    f(1,1,:) = fxx
+    f(1,2,:) = fxy
+    f(1,3,:) = fxz
+    f(2,1,:) = fxy
+    f(2,2,:) = fyy
+    f(2,3,:) = fyz
+    f(3,1,:) = fxz
+    f(3,2,:) = fyz
+    f(3,3,:) = fzz
+  end subroutine get_tensorderivs
+  subroutine get_scalarderivs2 (a, f, pos, off, dx)
+    CCTK_REAL, intent(in)  :: a(*)
+    CCTK_REAL, intent(out) :: f(3,3)
+    integer,   intent(in)  :: pos, off(3)
+    CCTK_REAL, intent(in)  :: dx(3)
+    call get_derivs2 (a, f, pos, off, dx)
+  end subroutine get_scalarderivs2
+  subroutine get_vectorderivs2 (ax,ay,az, f, pos, off, dx)
+    CCTK_REAL, intent(in)  :: ax(*),ay(*),az(*)
+    CCTK_REAL, intent(out) :: f(3,3,3)
+    integer,   intent(in)  :: pos, off(3)
+    CCTK_REAL, intent(in)  :: dx(3)
+    CCTK_REAL :: fx(3,3),fy(3,3),fz(3,3)
+    call get_derivs2 (ax, fx, pos, off, dx)
+    call get_derivs2 (ay, fy, pos, off, dx)
+    call get_derivs2 (az, fz, pos, off, dx)
+    f(1,:,:) = fx
+    f(2,:,:) = fy
+    f(3,:,:) = fz
+  end subroutine get_vectorderivs2
+  subroutine get_tensorderivs2 (axx,axy,axz,ayy,ayz,azz, f, pos, off, dx)
+    CCTK_REAL, intent(in)  :: axx(*),axy(*),axz(*),ayy(*),ayz(*),azz(*)
+    CCTK_REAL, intent(out) :: f(3,3,3,3)
+    integer,   intent(in)  :: pos, off(3)
+    CCTK_REAL, intent(in)  :: dx(3)
+    CCTK_REAL :: fxx(3,3),fxy(3,3),fxz(3,3),fyy(3,3),fyz(3,3),fzz(3,3)
+    call get_derivs2 (axx, fxx, pos, off, dx)
+    call get_derivs2 (axy, fxy, pos, off, dx)
+    call get_derivs2 (axz, fxz, pos, off, dx)
+    call get_derivs2 (ayy, fyy, pos, off, dx)
+    call get_derivs2 (ayz, fyz, pos, off, dx)
+    call get_derivs2 (azz, fzz, pos, off, dx)
+    f(1,1,:,:) = fxx
+    f(1,2,:,:) = fxy
+    f(1,3,:,:) = fxz
+    f(2,1,:,:) = fxy
+    f(2,2,:,:) = fyy
+    f(2,3,:,:) = fyz
+    f(3,1,:,:) = fxz
+    f(3,2,:,:) = fyz
+    f(3,3,:,:) = fzz
+  end subroutine get_tensorderivs2
+end module pointwise
diff --git a/src/pointwise2.F90 b/src/pointwise2.F90
new file mode 100644
index 0000000..1759b28
--- /dev/null
+++ b/src/pointwise2.F90
@@ -0,0 +1,148 @@
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+module pointwise2
+!!$  use derivs2
+  implicit none
+  DECLARE_CCTK_PARAMETERS
+  private
+  public calc_2position
+  public calc_2offsets
+  public get_2scalar
+  public get_2vector
+  public get_2tensor
+  public set_2scalar
+  public set_2vector
+  public set_2tensor
+  public get_2scalarderivs
+  public get_2vectorderivs
+  public get_2tensorderivs
+  public get_2scalarderivs2
+  public get_2vectorderivs2
+  public get_2tensorderivs2
+contains
+#define TGR_INCLUDED
+#include "derivs2.F90"
+#undef TGR_INCLUDED
+  subroutine calc_2position (shape, i,j, pos)
+    integer, intent(in)  :: shape(2)
+    integer, intent(in)  :: i,j
+    integer, intent(out) :: pos
+    pos = 1 + i-1 + shape(1) * (j-1)
+  end subroutine calc_2position
+  subroutine calc_2offsets (shape, off)
+    integer, intent(in)  :: shape(2)
+    integer, intent(out) :: off(2)
+    off(1) = 1
+    off(2) = shape(1)
+  end subroutine calc_2offsets
+  subroutine get_2scalar (a, f, pos)
+    CCTK_REAL, intent(in)  :: a(*)
+    CCTK_REAL, intent(out) :: f
+    integer,   intent(in)  :: pos
+    f = a(pos)
+  end subroutine get_2scalar
+  subroutine get_2vector (ax,ay, f, pos)
+    CCTK_REAL, intent(in)  :: ax(*),ay(*)
+    CCTK_REAL, intent(out) :: f(2)
+    integer,   intent(in)  :: pos
+    f(1) = ax(pos)
+    f(2) = ay(pos)
+  end subroutine get_2vector
+  subroutine get_2tensor (axx,axy,ayy, f, pos)
+    CCTK_REAL, intent(in)  :: axx(*),axy(*),ayy(*)
+    CCTK_REAL, intent(out) :: f(2,2)
+    integer,   intent(in)  :: pos
+    f(1,1) = axx(pos)
+    f(1,2) = axy(pos)
+    f(2,1) = axy(pos)
+    f(2,2) = ayy(pos)
+  end subroutine get_2tensor
+  subroutine set_2scalar (f, a, pos)
+    CCTK_REAL, intent(in) :: f
+    CCTK_REAL             :: a(*)
+    integer,   intent(in) :: pos
+    a(pos) = f
+  end subroutine set_2scalar
+  subroutine set_2vector (f, ax,ay, pos)
+    CCTK_REAL, intent(in) :: f(2)
+    CCTK_REAL             :: ax(*),ay(*)
+    integer,   intent(in) :: pos
+    ax(pos) = f(1)
+    ay(pos) = f(2)
+  end subroutine set_2vector
+  subroutine set_2tensor (f, axx,axy,ayy, pos)
+    CCTK_REAL, intent(in) :: f(2,2)
+    CCTK_REAL             :: axx(*),axy(*),ayy(*)
+    integer,   intent(in) :: pos
+    axx(pos) = f(1,1)
+    axy(pos) = f(1,2)
+    ayy(pos) = f(2,2)
+  end subroutine set_2tensor
+  subroutine get_2scalarderivs (a, f, pos, off, dx)
+    CCTK_REAL, intent(in)  :: a(*)
+    CCTK_REAL, intent(out) :: f(2)
+    integer,   intent(in)  :: pos, off(2)
+    CCTK_REAL, intent(in)  :: dx(2)
+    call get_2derivs (a, f, pos, off, dx)
+  end subroutine get_2scalarderivs
+  subroutine get_2vectorderivs (ax,ay, f, pos, off, dx)
+    CCTK_REAL, intent(in)  :: ax(*),ay(*)
+    CCTK_REAL, intent(out) :: f(2,2)
+    integer,   intent(in)  :: pos, off(2)
+    CCTK_REAL, intent(in)  :: dx(2)
+    CCTK_REAL :: fx(2),fy(2)
+    call get_2derivs (ax, fx, pos, off, dx)
+    call get_2derivs (ay, fy, pos, off, dx)
+    f(1,:) = fx
+    f(2,:) = fy
+  end subroutine get_2vectorderivs
+  subroutine get_2tensorderivs (axx,axy,ayy, f, pos, off, dx)
+    CCTK_REAL, intent(in)  :: axx(*),axy(*),ayy(*)
+    CCTK_REAL, intent(out) :: f(2,2,2)
+    integer,   intent(in)  :: pos, off(2)
+    CCTK_REAL, intent(in)  :: dx(2)
+    CCTK_REAL :: fxx(2),fxy(2),fyy(2)
+    call get_2derivs (axx, fxx, pos, off, dx)
+    call get_2derivs (axy, fxy, pos, off, dx)
+    call get_2derivs (ayy, fyy, pos, off, dx)
+    f(1,1,:) = fxx
+    f(1,2,:) = fxy
+    f(2,1,:) = fxy
+    f(2,2,:) = fyy
+  end subroutine get_2tensorderivs
+  subroutine get_2scalarderivs2 (a, f, pos, off, dx)
+    CCTK_REAL, intent(in)  :: a(*)
+    CCTK_REAL, intent(out) :: f(2,2)
+    integer,   intent(in)  :: pos, off(2)
+    CCTK_REAL, intent(in)  :: dx(2)
+    call get_2derivs2 (a, f, pos, off, dx)
+  end subroutine get_2scalarderivs2
+  subroutine get_2vectorderivs2 (ax,ay, f, pos, off, dx)
+    CCTK_REAL, intent(in)  :: ax(*),ay(*)
+    CCTK_REAL, intent(out) :: f(2,2,2)
+    integer,   intent(in)  :: pos, off(2)
+    CCTK_REAL, intent(in)  :: dx(2)
+    CCTK_REAL :: fx(2,2),fy(2,2)
+    call get_2derivs2 (ax, fx, pos, off, dx)
+    call get_2derivs2 (ay, fy, pos, off, dx)
+    f(1,:,:) = fx
+    f(2,:,:) = fy
+  end subroutine get_2vectorderivs2
+  subroutine get_2tensorderivs2 (axx,axy,ayy, f, pos, off, dx)
+    CCTK_REAL, intent(in)  :: axx(*),axy(*),ayy(*)
+    CCTK_REAL, intent(out) :: f(2,2,2,2)
+    integer,   intent(in)  :: pos, off(2)
+    CCTK_REAL, intent(in)  :: dx(2)
+    CCTK_REAL :: fxx(2,2),fxy(2,2),fyy(2,2)
+    call get_2derivs2 (axx, fxx, pos, off, dx)
+    call get_2derivs2 (axy, fxy, pos, off, dx)
+    call get_2derivs2 (ayy, fyy, pos, off, dx)
+    f(1,1,:,:) = fxx
+    f(1,2,:,:) = fxy
+    f(2,1,:,:) = fxy
+    f(2,2,:,:) = fyy
+  end subroutine get_2tensorderivs2
+end module pointwise2
diff --git a/src/rotation.F90 b/src/rotation.F90
new file mode 100644
index 0000000..d76973e
--- /dev/null
+++ b/src/rotation.F90
@@ -0,0 +1,59 @@
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+module rotation
+  use constants
+  implicit none
+  DECLARE_CCTK_PARAMETERS
+  private
+  
+  public make_euler
+  
+contains
+  
+  subroutine make_euler1 (phi, tt)
+    CCTK_REAL, intent(in)  :: phi
+    CCTK_REAL, intent(out) :: tt(0:3,0:3)
+    ! y^a = T^a_b x^b
+    tt = delta4
+    tt(1,1) =  cos(phi)
+    tt(1,2) = -sin(phi)
+    tt(2,1) =  sin(phi)
+    tt(2,2) =  cos(phi)
+  end subroutine make_euler1
+  
+  subroutine make_euler2 (theta, tt)
+    CCTK_REAL, intent(in)  :: theta
+    CCTK_REAL, intent(out) :: tt(0:3,0:3)
+    ! y^a = T^a_b x^b
+    tt = delta4
+    tt(2,2) =  cos(theta)
+    tt(2,3) = -sin(theta)
+    tt(3,2) =  sin(theta)
+    tt(3,3) =  cos(theta)
+  end subroutine make_euler2
+  
+  subroutine make_euler3 (psi, tt)
+    CCTK_REAL, intent(in)  :: psi
+    CCTK_REAL, intent(out) :: tt(0:3,0:3)
+    ! y^a = T^a_b x^b
+    tt = delta4
+    tt(1,1) =  cos(psi)
+    tt(1,2) = -sin(psi)
+    tt(2,1) =  sin(psi)
+    tt(2,2) =  cos(psi)
+  end subroutine make_euler3
+  
+  subroutine make_euler (phi, theta, psi, tt)
+    CCTK_REAL, intent(in)  :: phi, theta, psi
+    CCTK_REAL, intent(out) :: tt(0:3,0:3)
+    CCTK_REAL :: tt1(0:3,0:3), tt2(0:3,0:3), tt3(0:3,0:3)
+    call make_euler1 (phi, tt1)
+    call make_euler2 (theta, tt2)
+    call make_euler3 (psi, tt3)
+    tt = matmul(tt3, matmul(tt2, tt1))
+  end subroutine make_euler
+  
+end module rotation
diff --git a/src/tensor.F90 b/src/tensor.F90
new file mode 100644
index 0000000..0fef831
--- /dev/null
+++ b/src/tensor.F90
@@ -0,0 +1,127 @@
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+module tensor
+  implicit none
+  DECLARE_CCTK_PARAMETERS
+  private
+  
+  public calc_trace
+  
+  public calc_det
+  public calc_detderiv
+  public calc_detdot
+  
+  public calc_inv
+  public calc_invderiv
+  public calc_invdot
+  
+contains
+  
+  subroutine calc_trace (kk, gu, trk)
+    CCTK_REAL, intent(in)  :: kk(3,3)
+    CCTK_REAL, intent(in)  :: gu(3,3)
+    CCTK_REAL, intent(out) :: trk
+    integer :: i,j
+    trk = 0
+    do i=1,3
+       do j=1,3
+          trk = trk + gu(i,j) * kk(i,j)
+       end do
+    end do
+  end subroutine calc_trace
+  
+  
+  
+  subroutine calc_det (gg, dtg)
+    CCTK_REAL, intent(in)  :: gg(3,3)
+    CCTK_REAL, intent(out) :: dtg
+    dtg =      gg(1,1) * gg(2,2) * gg(3,3) &
+         & + 2*gg(1,2) * gg(1,3) * gg(2,3) &
+         & -   gg(1,1) * gg(2,3)**2        &
+         & -   gg(2,2) * gg(1,3)**2        &
+         & -   gg(3,3) * gg(1,2)**2
+  end subroutine calc_det
+  
+  subroutine calc_pdet (gg, pgg, pdtg)
+    CCTK_REAL, intent(in)  :: gg(3,3), pgg(3,3)
+    CCTK_REAL, intent(out) :: pdtg
+    pdtg =     pgg(1,1) *   gg(2,2) *  gg(3,3) &
+         & +    gg(1,1) *  pgg(2,2) *  gg(3,3) &
+         & +    gg(1,1) *   gg(2,2) * pgg(3,3) &
+         & + 2*pgg(1,2) *   gg(1,3) *  gg(2,3) &
+         & + 2* gg(1,2) *  pgg(1,3) *  gg(2,3) &
+         & + 2* gg(1,2) *   gg(1,3) * pgg(2,3) &
+         & -   pgg(1,1) *   gg(2,3)**2         &
+         & -    gg(1,1) * 2*gg(2,3) * pgg(2,3) &
+         & -   pgg(2,2) *   gg(1,3)**2         &
+         & -    gg(2,2) * 2*gg(1,3) * pgg(1,3) &
+         & -   pgg(3,3) *   gg(1,2)**2         &
+         & -    gg(3,3) * 2*gg(1,2) * pgg(1,2)
+  end subroutine calc_pdet
+  
+  subroutine calc_detderiv (gg, dgg, ddtg)
+    CCTK_REAL, intent(in)  :: gg(3,3), dgg(3,3,3)
+    CCTK_REAL, intent(out) :: ddtg(3)
+    integer :: i
+    do i=1,3
+       call calc_pdet (gg, dgg(:,:,i), ddtg(i))
+    end do
+  end subroutine calc_detderiv
+  
+  subroutine calc_detdot (gg, gg_dot, dtg_dot)
+    CCTK_REAL, intent(in)  :: gg(3,3), gg_dot(3,3)
+    CCTK_REAL, intent(out) :: dtg_dot
+    call calc_pdet (gg, gg_dot, dtg_dot)
+  end subroutine calc_detdot
+  
+  
+  
+  subroutine calc_inv (gg, dtg, gu)
+    CCTK_REAL, intent(in)  :: gg(3,3), dtg
+    CCTK_REAL, intent(out) :: gu(3,3)
+    gu(1,1) = (gg(2,2) * gg(3,3) - gg(2,3) ** 2     ) / dtg
+    gu(2,2) = (gg(1,1) * gg(3,3) - gg(1,3) ** 2     ) / dtg
+    gu(3,3) = (gg(1,1) * gg(2,2) - gg(1,2) ** 2     ) / dtg
+    gu(1,2) = (gg(1,3) * gg(2,3) - gg(1,2) * gg(3,3)) / dtg
+    gu(1,3) = (gg(1,2) * gg(2,3) - gg(1,3) * gg(2,2)) / dtg
+    gu(2,3) = (gg(1,3) * gg(1,2) - gg(2,3) * gg(1,1)) / dtg
+    gu(2,1) = gu(1,2)
+    gu(3,1) = gu(1,3)
+    gu(3,2) = gu(2,3)
+  end subroutine calc_inv
+  
+  subroutine calc_pinv (gu, pgg, pgu)
+    CCTK_REAL, intent(in)  :: gu(3,3), pgg(3,3)
+    CCTK_REAL, intent(out) :: pgu(3,3)
+    integer :: i,j,k,l
+    do i=1,3
+       do j=1,3
+          pgu(i,j) = 0
+          do k=1,3
+             do l=1,3
+                pgu(i,j) = pgu(i,j) - gu(i,k) * gu(j,l) * pgg(k,l)
+             end do
+          end do
+       end do
+    end do
+  end subroutine calc_pinv
+  
+  subroutine calc_invderiv (gu, dgg, dgu)
+    CCTK_REAL, intent(in)  :: gu(3,3), dgg(3,3,3)
+    CCTK_REAL, intent(out) :: dgu(3,3,3)
+    integer :: i
+    do i=1,3
+       call calc_pinv (gu, dgg(:,:,i), dgu(:,:,i))
+    end do
+  end subroutine calc_invderiv
+  
+  subroutine calc_invdot (gu, gg_dot, gu_dot)
+    CCTK_REAL, intent(in)  :: gu(3,3), gg_dot(3,3)
+    CCTK_REAL, intent(out) :: gu_dot(3,3)
+    call calc_pinv (gu, gg_dot, gu_dot)
+  end subroutine calc_invdot
+  
+end module tensor
diff --git a/src/tensor2.F90 b/src/tensor2.F90
new file mode 100644
index 0000000..b02cd06
--- /dev/null
+++ b/src/tensor2.F90
@@ -0,0 +1,107 @@
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+module tensor2
+  implicit none
+  DECLARE_CCTK_PARAMETERS
+  private
+  
+  public calc_2trace
+  
+  public calc_2det
+  public calc_2detderiv
+  public calc_2detdot
+  
+  public calc_2inv
+  public calc_2invderiv
+  public calc_2invdot
+  
+contains
+  
+  subroutine calc_2trace (kk, gu, trk)
+    CCTK_REAL, intent(in)  :: kk(2,2)
+    CCTK_REAL, intent(in)  :: gu(2,2)
+    CCTK_REAL, intent(out) :: trk
+    integer :: i,j
+    trk = 0
+    do i=1,2
+       do j=1,2
+          trk = trk + gu(i,j) * kk(i,j)
+       end do
+    end do
+  end subroutine calc_2trace
+  
+  
+  
+  subroutine calc_2det (gg, dtg)
+    CCTK_REAL, intent(in)  :: gg(2,2)
+    CCTK_REAL, intent(out) :: dtg
+    dtg =  gg(1,1) * gg(2,2) - gg(1,2)**2
+  end subroutine calc_2det
+  
+  subroutine calc_2pdet (gg, pgg, pdtg)
+    CCTK_REAL, intent(in)  :: gg(2,2), pgg(2,2)
+    CCTK_REAL, intent(out) :: pdtg
+    pdtg = pgg(1,1) * gg(2,2) + gg(1,1) * pgg(2,2) - 2 * gg(1,2) * pgg(1,2)
+  end subroutine calc_2pdet
+  
+  subroutine calc_2detderiv (gg, dgg, ddtg)
+    CCTK_REAL, intent(in)  :: gg(2,2), dgg(2,2,2)
+    CCTK_REAL, intent(out) :: ddtg(2)
+    integer :: i
+    do i=1,2
+       call calc_2pdet (gg, dgg(:,:,i), ddtg(i))
+    end do
+  end subroutine calc_2detderiv
+  
+  subroutine calc_2detdot (gg, gg_dot, dtg_dot)
+    CCTK_REAL, intent(in)  :: gg(2,2), gg_dot(2,2)
+    CCTK_REAL, intent(out) :: dtg_dot
+    call calc_2pdet (gg, gg_dot, dtg_dot)
+  end subroutine calc_2detdot
+  
+  
+  
+  subroutine calc_2inv (gg, dtg, gu)
+    CCTK_REAL, intent(in)  :: gg(2,2), dtg
+    CCTK_REAL, intent(out) :: gu(2,2)
+    gu(1,1) = gg(2,2) / dtg
+    gu(2,2) = gg(1,1) / dtg
+    gu(1,2) = gg(1,2) / dtg
+    gu(2,1) = gu(1,2)
+  end subroutine calc_2inv
+  
+  subroutine calc_2pinv (gu, pgg, pgu)
+    CCTK_REAL, intent(in)  :: gu(2,2), pgg(2,2)
+    CCTK_REAL, intent(out) :: pgu(2,2)
+    integer :: i,j,k,l
+    do i=1,2
+       do j=1,2
+          pgu(i,j) = 0
+          do k=1,2
+             do l=1,2
+                pgu(i,j) = pgu(i,j) - gu(i,k) * gu(j,l) * pgg(k,l)
+             end do
+          end do
+       end do
+    end do
+  end subroutine calc_2pinv
+  
+  subroutine calc_2invderiv (gu, dgg, dgu)
+    CCTK_REAL, intent(in)  :: gu(2,2), dgg(2,2,2)
+    CCTK_REAL, intent(out) :: dgu(2,2,2)
+    integer :: i
+    do i=1,2
+       call calc_2pinv (gu, dgg(:,:,i), dgu(:,:,i))
+    end do
+  end subroutine calc_2invderiv
+  
+  subroutine calc_2invdot (gu, gg_dot, gu_dot)
+    CCTK_REAL, intent(in)  :: gu(2,2), gg_dot(2,2)
+    CCTK_REAL, intent(out) :: gu_dot(2,2)
+    call calc_2pinv (gu, gg_dot, gu_dot)
+  end subroutine calc_2invdot
+  
+end module tensor2