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