C Subroutine file 'metric.for' to define the coordinate system
C
C By Vlastislav Cerveny, Ludek Klimes, Ivan Psencik
C
C This file consists of the following subroutines:
C     METRB...Block data subroutine defining the common block /METRC/ to
C             store the data describing the coordinate system.
C     METR1...Subroutine designed to store the data into the common
C             block /METRC/.
C     KOOR... Integer function returning the type of the coordinate
C             system.
C     METRIC..Subroutine designed to evaluate the metric tensor and
C             Christoffel symbols at a given point.
C     CARTES..Subroutine designed to transform the model coordinates to
C             Cartesian coordinates and vice versa.  The indexing of
C             coordinate systems should correspond to the subroutine
C             METRIC.
C
C Storage in the memory:
C     The data describing the coordinate system are stored in the common
C     block /METRC/ defined in the following subroutine:
C     ------------------------------------------------------------------
      BLOCK DATA METRB
      INCLUDE 'metric.inc'
      END
C     ------------------------------------------------------------------
C
C Date: 1992, November 2
C Coded by Ludek Klimes
C
C=======================================================================
C
      SUBROUTINE METR1(KOOR)
      INTEGER KOOR
C
C Subroutine METR1 is designed to store the data specifying the
C coordinate system into the common block /METRC/.
C
C Input:
C     KOOR... Specifies the type of the right-handed coordinate system.
C The input parameter is not altered.
C
C No output.
C
C Common block /METRC/:
      INCLUDE 'metric.inc'
C All the storage locations of the common block are defined in this
C subroutine.
C
C Subroutines and external functions required:
      EXTERNAL METRB
C     METRB.. Block data subroutine of this file.
C
C Date: 1993, December 18
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     No auxiliary storage locations.
C
      KOORS=KOOR
      RETURN
      END
C
C=======================================================================
C
      INTEGER FUNCTION KOOR()
C
C Integer function KOOR is designed to return the type of the coordinate
C system.
C
C No input.
C
C Output:
C     KOOR... Specifies the type of the right-handed coordinate system:
C             KOOR.LE.0: Cartesian coordinates.
C             KOOR.EQ.1: Polar spherical coordinates (X1,X2,X3)=
C               (COLATITUDE,LONGITUDE,RADIUS).
C             KOOR.GE.2: Geographic spherical coordinates (X1,X2,X3)=
C               (longitude,latitude,radius).
C
C Common block /METRC/:
      INCLUDE 'metric.inc'
C None of the storage locations of the common block are altered.
C
C No subroutines and external functions required.
C
C Date: 1989, December 18
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     No auxiliary storage locations.
C
      KOOR=KOORS
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE METRIC(COOR,GSQRD,G,GAMMA)
      REAL COOR(3),GSQRD,G(12),GAMMA(18)
C
C This subroutine evaluates the metric tensor and Christoffel
C symbols at a given point.
C
C Input:
C     COOR... Array containing coordinates X1, X2, X3 of the given point
C None of the input parameters are altered.
C
C Output:
C     GSQRD...Square root of the determinant of the covariant metric
C             tensor.
C     G...    Array containing covariant components G11, G12, G22, G13,
C             G23, G33, and contravariant components G11, G12, G22, G13,
C             G23, G33 of the metric tensor at the given point.
C     GAMMA...Array containing Christoffel symbols GAMMA111, GAMMA121,
C             GAMMA221, GAMMA131, GAMMA231, GAMMA331, GAMMA112,
C             GAMMA122, GAMMA222, GAMMA132, GAMMA232, GAMMA332,
C             GAMMA113, GAMMA123, GAMMA223, GAMMA133, GAMMA233,
C             GAMMA333, where the first two indices are subscripts and
C             the third index is a superscript.
C
C Common block /METRC/:
      INCLUDE 'metric.inc'
C None of the storage locations of the common block are altered.
C
C No subroutines and external functions required.
C
C Date: 1991, May 19
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
      INTEGER I
      REAL SMALL,C,S,R
      PARAMETER (SMALL=1.E-12)
C
C     I...    Auxiliary loop variable.
C     SMALL...The lower limit for the distance from the singularities of
C             the coordinate system measured in the coordinate units.
C     C,S,R...Auxiliary storage locations.
C
C.......................................................................
C
      DO 1 I=1,12
        G(I)=0.
    1 CONTINUE
      DO 2 I=1,18
        GAMMA(I)=0.
    2 CONTINUE
C
      IF(KOORS.LE.0) THEN
        GSQRD=1.
        G(1) =1.
        G(3) =1.
        G(6) =1.
        G(7) =1.
        G(9) =1.
        G(12)=1.
      ELSE IF(KOORS.EQ.1) THEN
        C=COS(COOR(1))
        S=SIN(COOR(1))
        R=COOR(3)
        IF(S.EQ.0.) S=SMALL
        IF(R.EQ.0.) R=SMALL
        GSQRD=R*R*ABS(S)
        G(1) =R*R
        G(3) =(R*S)**2
        G(6) =1.
        G(7) =1./G(1)
        G(9) =1./G(3)
        G(12)=1.
        GAMMA(3) =-S*C
        GAMMA(4) =1./R
        GAMMA(8) =C/S
        GAMMA(11)=1./R
        GAMMA(13)=-R
        GAMMA(15)=-R*S*S
      ELSE
        C=COS(COOR(2))
        S=SIN(COOR(2))
        R=COOR(3)
        IF(C.EQ.0.) C=SMALL
        IF(R.EQ.0.) R=SMALL
        GSQRD=R*R*ABS(C)
        G(1)=(R*C)**2
        G(3)=R*R
        G(6)=1.
        G(7)=1./G(1)
        G(9)=1./G(3)
        G(12)=1.
        GAMMA(3)=-S/C
        GAMMA(4)=1./R
        GAMMA(8)=S*C
        GAMMA(11)=1./R
        GAMMA(13)=-R*C*C
        GAMMA(15)=-R
      END IF
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE CARTES(COOR,TO,CART,PDER)
      LOGICAL TO
      REAL COOR(3),CART(3),PDER(9)
C
C This subroutine transforms the model coordinates to Cartesian
C coordinates and vice versa.  This subroutine has to correspond to the
C subroutine METRIC.
C
C Arguments:
C     COOR... Array containing the model coordinates X1,X2,X3 of the
C             given point.
C     TO...   .TRUE. To transform the model coordinates to the Cartesian
C               coordinates.
C               Input: COOR,
C               Output: CART,PDER.
C             .FALSE. To transform the Cartesian coordinates to the
C               model coordinates.
C               Input: CART,
C               Output: COOR,PDER.
C     CART... Array containing the Cartesian coordinates C1, C2, C3 of
C             the given point.
C     PDER... Partial derivatives of the output coordinates with respect
C             to the input coordinates.  I.e. the transformation matrix
C             of contravariant vectors, corresponding to the coordinate
C             transformation.  I.e. the transposed transformation matrix
C             of covariant vectors, corresponding to the inverse
C             transformation.
C
C Subroutines and external functions required:
      EXTERNAL KOOR
      INTEGER  KOOR
C     KOOR...  This file.
C
C Date: 1996, September 30
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER I
      REAL C,S,R
C
C.......................................................................
C
      IF(KOOR().LE.0) THEN
C       Cartesian coordinates:
        IF(TO) THEN
          CART(1)=COOR(1)
          CART(2)=COOR(2)
          CART(3)=COOR(3)
        ELSE
          COOR(1)=CART(1)
          COOR(2)=CART(2)
          COOR(3)=CART(3)
        END IF
        DO 11 I=2,8
          PDER(I)=0.
   11   CONTINUE
        DO 12 I=1,9,4
          PDER(I)=1.
   12   CONTINUE
      ELSE IF(KOOR().EQ.1) THEN
C       Polar spherical coordinates:
        IF(TO) THEN
          R=COOR(3)
          S=R*SIN(COOR(1))
          CART(1)=S*COS(COOR(2))
          CART(2)=S*SIN(COOR(2))
          CART(3)=R*COS(COOR(1))
          PDER(1)= CART(1)*CART(3)/S
          PDER(2)= CART(2)*CART(3)/S
          PDER(3)=-S
          PDER(4)=-CART(2)
          PDER(5)= CART(1)
          PDER(6)= 0.
          PDER(7)= CART(1)/R
          PDER(8)= CART(2)/R
          PDER(9)= CART(3)/R
        ELSE
          S=CART(1)**2+CART(2)**2
          R=SQRT(S+CART(3)**2)
          S=SQRT(S)
          IF(R.NE.0.) THEN
            COOR(1)=ATAN2(S,CART(3))
          ELSE
            COOR(1)=0.
          END IF
          IF(S.NE.0.) THEN
            COOR(2)=ATAN2(CART(2),CART(1))
          ELSE
            COOR(2)=0.
          END IF
          COOR(3)=R
          PDER(1)= CART(1)*CART(3)/S/R
          PDER(2)=-CART(2)/S
          PDER(3)= CART(1)/R
          PDER(4)= CART(2)*CART(3)/S/R
          PDER(5)= CART(1)/S
          PDER(6)= CART(2)/R
          PDER(7)=-S/R
          PDER(8)= 0.
          PDER(9)= CART(3)/R
        END IF
      ELSE
C       Geographic spherical coordinates:
        IF(TO) THEN
          R=COOR(3)
          C=R*COS(COOR(2))
          CART(1)=C*COS(COOR(1))
          CART(2)=C*SIN(COOR(1))
          CART(3)=R*SIN(COOR(2))
          PDER(1)=-CART(2)
          PDER(2)= CART(1)
          PDER(3)= 0.
          PDER(4)=-CART(1)*CART(3)/C
          PDER(5)=-CART(2)*CART(3)/C
          PDER(6)= C
          PDER(7)= CART(1)/R
          PDER(8)= CART(2)/R
          PDER(9)= CART(3)/R
        ELSE
          C=CART(1)**2+CART(2)**2
          R=SQRT(C+CART(3)**2)
          C=SQRT(C)
          IF(R.NE.0.) THEN
            COOR(2)= ATAN2(CART(3),C)
            PDER(3)= CART(1)/R
            PDER(6)= CART(2)/R
            PDER(8)= C/R
            PDER(9)= CART(3)/R
          ELSE
            COOR(2)= 0.
            PDER(3)= 1.
            PDER(6)= 1.
            PDER(8)= 0.
            PDER(9)= 1.
          END IF
          IF(C.NE.0.) THEN
            COOR(1)= ATAN2(CART(2),CART(1))
            PDER(1)=-CART(2)/C
            PDER(2)=-CART(1)*CART(3)/C/R
            PDER(4)= CART(1)/C
            PDER(5)=-CART(2)*CART(3)/C/R
          ELSE
            COOR(1)= 0.
            PDER(1)= 0.
            PDER(2)=-1.
            PDER(4)= 0.
            PDER(5)=-1.
          END IF
          COOR(3)= R
          PDER(7)= 0.
        END IF
      END IF
      RETURN
      END
C
C=======================================================================
C