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
C 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.
C             File MODEL must be already opened and the position in the
C             file must be just after reading
C             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                                                                 ```