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
      INTEGER KOORS
      COMMON/METRC/KOORS
      SAVE /METRC/
      END
C     ------------------------------------------------------------------
C     KOORS...SPECIFIES THE TYPE OF THE RIGHT-HANDED COORDINATE SYSTEM:
C             KOORS.LE.0: CARTESIAN COORDINATES.
C             KOORS.EQ.1: POLAR SPHERICAL COORDINATES (X1,X2,X3)=
C               (COLATITUDE,LONGITUDE,RADIUS).
C             KOORS.GE.2: GEOGRAPHIC SPHERICAL COORDINATES (X1,X2,X3)=
C               (LONGITUDE,LATITUDE,RADIUS).
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/:
      INTEGER KOORS
      COMMON/METRC/KOORS
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/:
      INTEGER KOORS
      COMMON/METRC/KOORS
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/:
      INTEGER KOORS
      COMMON/METRC/KOORS
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: 1991, JUNE 12
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
        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
        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
        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)
          ELSE
            COOR(2)=0.
          END IF
          IF(S.NE.0.) THEN
            COOR(1)=ATAN2(CART(2),CART(1))
          ELSE
            COOR(1)=0.
          END IF
          COOR(3)=R
          PDER(1)=-CART(2)/C
          PDER(2)=-CART(1)*CART(3)/C/R
          PDER(3)= CART(1)/R
          PDER(4)= CART(1)/C
          PDER(5)=-CART(2)*CART(3)/C/R
          PDER(6)= CART(2)/R
          PDER(7)= 0.
          PDER(8)= C/R
          PDER(9)= CART(3)/R
        END IF
      END IF
      RETURN
      END
C
C=======================================================================
C
