C 
C Subroutine file 'metric.for' to define the coordinate system
C
C Date: 2006, March 4
C Coded by Ludek Klimes
C
C.......................................................................
C
C This file consists of the following subroutines:
C     METR1...Subroutine designed to store the data into the common
C             block /METRC/.
C             METR1
C     KOOR... Integer function returning the type of the coordinate
C             system.
C             KOOR
C     METRIC..Subroutine designed to evaluate the metric tensor and
C             Christoffel symbols at a given point.
C             METRIC
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             CARTES
C
C.......................................................................
C
C                                                    
C Additional input data for the coordinate system:
C The data depend on the value of KOORS.
C No data for KOORS=0, KOORS=1 or KOORS=2.
C For KOORS=3 or KOORS=4 following line:
C (1) CIRCLE,RADIUS,/
C     CIRCLE..Number angular units per a circle for spherical
C             coordinates.
C             For example, CIRCLE=6.2831853 means radians,
C             CIRCLE=360 means degrees.
C             Default: CIRCLE=6.2831853
C     RADIUS..Reference Earth radius defining elevation X3 in spherical
C             coordinates.
C             The radius in spherical coordinates is RADIUS+X3.
C             Default: RADIUS=0.
C
C.......................................................................
C
C Storage in the memory:
C     The data describing the coordinate system are stored in the common
C     block /METRC/ defined in the include file 'metric.inc'.
C     metric.inc
C
C=======================================================================
C
C     
C
      SUBROUTINE METR1(KOOR,LUN)
      INTEGER KOOR,LUN
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     LUN ... Number of the logical unit connected to the
C             file MODEL. File MODEL must
C             be already opened and the position in the file must be
C             just after reading data (2).
C             LUN is required only if KOOR.EQ.3 or KOOR.EQ.4.
C The input parameters are not altered.
C
C No output.
C
C Common block /METRC/:
      INCLUDE 'metric.inc'
C     metric.inc
C All the storage locations of the common block are defined in this
C subroutine.
C
C No subroutines and external functions required.
C
C Date: 2006, March 4
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage location:
      REAL PI2
      PARAMETER (PI2=6.2831853)
C
      KOORS=KOOR
      RADIAN=1.
      RADIUS=0.
      IF(KOORS.EQ.3.OR.KOORS.EQ.4) THEN
        KOORS=KOORS-2
        RADIAN=PI2
        READ(LUN,*) RADIAN,RADIUS
        RADIAN=PI2/RADIAN
      END IF
      RETURN
      END
C
C=======================================================================
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     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
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     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: 2006, February 23
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(RADIAN*COOR(1))
        S=SIN(RADIAN*COOR(1))
        R=RADIUS+COOR(3)
        IF(S.EQ.0.) S=SMALL
        IF(R.EQ.0.) R=SMALL
        R=RADIAN*R
        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) =-RADIAN*S*C
        GAMMA(4) = RADIAN/R
        GAMMA(8) = RADIAN*C/S
        GAMMA(11)= GAMMA(4)
        GAMMA(13)=-RADIAN*R
        GAMMA(15)=-RADIAN*R*S*S
      ELSE
        C=COS(RADIAN*COOR(2))
        S=SIN(RADIAN*COOR(2))
        R=RADIUS+COOR(3)
        IF(C.EQ.0.) C=SMALL
        IF(R.EQ.0.) R=SMALL
        R=RADIAN*R
        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(2)= -RADIAN*S/C
        GAMMA(4)=  RADIAN/R
        GAMMA(7)=  RADIAN*S*C
        GAMMA(11)= GAMMA(4)
        GAMMA(13)=-RADIAN*R*C*C
        GAMMA(15)=-RADIAN*R
      END IF
      RETURN
      END
C
C=======================================================================
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 Common block /METRC/:
      INCLUDE 'metric.inc'
C     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: 2006, February 21
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER I
      REAL C1,S1,C2,S2,C,S,R
C
C.......................................................................
C
      IF(KOORS.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(KOORS.EQ.1) THEN
C       Polar spherical coordinates:
        IF(TO) THEN
          C1=COS(RADIAN*COOR(1))
          S1=SIN(RADIAN*COOR(1))
          C2=COS(RADIAN*COOR(2))
          S2=SIN(RADIAN*COOR(2))
          R=RADIUS+COOR(3)
          S=R*S1
          CART(1)=S*C2
          CART(2)=S*S2
          CART(3)=R*C1
          PDER(1)= RADIAN*CART(3)*C2
          PDER(2)= RADIAN*CART(3)*S2
          PDER(3)=-RADIAN*S
          PDER(4)=-RADIAN*CART(2)
          PDER(5)= RADIAN*CART(1)
          PDER(6)= 0.
          PDER(7)= S1*C2
          PDER(8)= S1*S2
          PDER(9)= C1
        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))/RADIAN
            PDER(3)= CART(1)/R
            PDER(6)= CART(2)/R
            PDER(7)=-S/R                 /R/RADIAN
            PDER(9)= CART(3)/R
          ELSE
            COOR(1)= 0.
            PDER(3)= 1.
            PDER(6)= 0.
            PDER(7)=-1.
            PDER(9)= 0.
          END IF
          IF(S.NE.0.) THEN
            COOR(2)= ATAN2(CART(2),CART(1))/RADIAN
            PDER(1)= CART(1)*CART(3)/S/R /R/RADIAN
            PDER(2)=-CART(2)/S           /S/RADIAN
            PDER(4)= CART(2)*CART(3)/S/R /R/RADIAN
            PDER(5)= CART(1)/S           /S/RADIAN
          ELSE
            COOR(2)= 0.
            PDER(1)= 0.
            PDER(2)= 0.
            PDER(4)= 0.
            PDER(5)= 1.
          END IF
          COOR(3)=R-RADIUS
          PDER(8)= 0.
        END IF
      ELSE
C       Geographic spherical coordinates:
        IF(TO) THEN
          C1=COS(RADIAN*COOR(1))
          S1=SIN(RADIAN*COOR(1))
          C2=COS(RADIAN*COOR(2))
          S2=SIN(RADIAN*COOR(2))
          R=RADIUS+COOR(3)
          C=R*C2
          CART(1)=C*C1
          CART(2)=C*S1
          CART(3)=R*S2
          PDER(1)=-RADIAN*CART(2)
          PDER(2)= RADIAN*CART(1)
          PDER(3)= 0.
          PDER(4)=-RADIAN*CART(3)*C1
          PDER(5)=-RADIAN*CART(3)*S1
          PDER(6)= RADIAN*C
          PDER(7)= C2*C1
          PDER(8)= C2*S1
          PDER(9)= S2
        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)/RADIAN
            PDER(3)= CART(1)/R
            PDER(6)= CART(2)/R
            PDER(8)= C/R                 /R/RADIAN
            PDER(9)= CART(3)/R
          ELSE
            COOR(2)= 0.
            PDER(3)= 1.
            PDER(6)= 0.
            PDER(8)= 1.
            PDER(9)= 0.
          END IF
          IF(C.NE.0.) THEN
            COOR(1)= ATAN2(CART(2),CART(1))/RADIAN
            PDER(1)=-CART(2)/C           /C/RADIAN
            PDER(2)=-CART(1)*CART(3)/C/R /R/RADIAN
            PDER(4)= CART(1)/C           /C/RADIAN
            PDER(5)=-CART(2)*CART(3)/C/R /R/RADIAN
          ELSE
            COOR(1)= 0.
            PDER(1)= 0.
            PDER(2)= 0.
            PDER(4)= 1.
            PDER(5)= 0.
          END IF
          COOR(3)= R-RADIUS
          PDER(7)= 0.
        END IF
      END IF
      RETURN
      END
C
C=======================================================================
C