C Subroutine file 'model.for' to specify a seismic model.
C
C By Vlastislav Cerveny, Ludek Klimes, Ivan Psencik
C
C This file consists of the following external procedures:
C MODELB..Block data subroutine defining common blocks MODELT and
C MODELC to store the data describing the model.
C MODEL1..Subroutine designed to read the input data for the
C description of the model and to store them in the memory.
C Subroutine MODEL1 reads the input data (1)-(8) listed
C below itself, then calls subroutine SRFC1 (included in the
C subroutine file 'srfc.for') to read the input data (9) for
C smooth surfaces, and finaly calls subroutine PARM1
C (included in the subroutine file 'parm.for') to read the
C input data (10) for the parameters of the medium.
C NSRFC...Integer function returning the number of surfaces covering
C structural interfaces.
C BLOCK...Subroutine designed to determine the mutual position of a
C point and a simple and a complex block.
C ISIDE...Auxiliary function to subroutine BLOCK.
C INTERF..Auxiliary subroutine to subroutine BLOCK.
C VELOC...Subroutine transforming the values of a medium parameters
C into velocities and loss factors.
C FPOWER..Subroutine evaluating the value and, possibly, the three
C first and six second partial derivatives of a function, if
C the value and the three first and six second partial
C derivatives of the POWER-th power of the function are
C are known. Particularly, this is an auxiliary subroutine
C to the subroutine VELOC.
C
C Attention:
C The lines denoted by '*V' in the first 2 columns should be removed
C from the subroutines veloc and fpower before compilation.
C If the model variations with respect to the model parameters were
C to be calculated, these lines would remain while replacing each
C '*V' in the first columns by spaces. In such a case, subroutines
C VAR4 and VAR5 of the file 'var.for' would be called to handle the
C variations.
C
C Input data for the specification of the model:
C The data are read in by the list directed input (free format). In
C the list of input data below, each numbered paragraph indicates
C the beginning of a new input operation (new READ statement). If
C the first letter of the symbolic name of the input variable is
C I-N, the corresponding value in input data must be of the type
C INTEGER. Otherwise (except TEXTM), the input parameter is of the
C type REAL.
C (1) TEXTM
C The string containing the name of the model. Only the first 80
C characters of the string are significant.
C (2) KOORS,NEXPV,NEXPQ,IVERT,/
C KOORS...Specifies the type of the right-handed coordinate system
C (see the subroutine file 'metric.for'):
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 NEXPV,NEXPQ... Specify exponents of the power of velocities
C (NEXPV) and loss factors (NEXPQ) in input data. For
C example, unit values of NEXPV and nexpq indicate that the
C parameters of the medium are velocities and loss factors,
C indices equal -1 indicate reciprocal values of these
C quantities, i.e. slownesses and quality factors.
C While using the basic submitted version of the subroutine
C file 'parm.for', the specification of NEXPV=1, NEXPQ=1 is
C highly recommended.
C IVERT...Orientation of the vertical axis:
C IVERT=0: unknown (default),
C IVERT=+1: X1 vertical, pointing upwards,
C IVERT=-1: X1 vertical, pointing downwards,
C IVERT=+2: X2 vertical, pointing upwards,
C IVERT=-2: X2 vertical, pointing downwards,
C IVERT=+3: X3 vertical, pointing upwards,
C IVERT=-3: X3 vertical, pointing downwards,
C Has no influence on the calculations, and need not be
C specified. If it is non-zero, it may be considered for
C plotting purposes.
C /... Obligatory slash at the end of line for future extensions.
C Default: KOORS=0, NEXPV=1, NEXPQ=1, IVERT=0.
C (3) X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX
C Boundaries of the model.
C (4) NSRFC
C Number of smooth surfaces in the model. The surfaces are indexed
C sequentially in any order by positive integers ISRFC from 1 to
C NSRFC.
C (5) NSB
C Number of simple blocks in the model, defined in this data set.
C The defined blocks are indexed sequentially by positive integers
C ISB from 1 to NSB in the same order as they are specified in data
C (6). All material simple blocks in the model must be defined.
C Free-space blocks are usually not defined in this data set,
C however they may be defined here at user's discretion.
C Intersecting simple blocks are allowed but not recommended.
C (6) NSB input operations (READ statements):
C For each simple block with index ISB, the indices of the surfaces
C forming the set F(+) and the indices of the surfaces forming the
C set F(-). The indices of surfaces from F(+) must be positive, the
C indices of surfaces from F(-) must be indicated by negative signs.
C The indices may be specified in an arbitrary order and must be
C terminated by a slash.
C (7) NCB
C Number of material complex blocks in the model. The blocks are
C indexed sequentially by positive integers ICB from 1 to NCB. The
C free-space blocks are not indexed.
C (8) NCB input operations (READ statements):
C For each material complex block, the indices of material simple
C blocks forming the complex block. The indices may be specified
C in an arbitrary order and must be terminated by a slash.
C Each material simple block must appear exactly once within these
C data lines. Simple blocks defined by data (6) but not listed here
C are deemed to be free-space simple blocks.
C (9) The data specifying NSRFC functions F(X1,X2,X3) describing the
C smooth surfaces in the model. The data are read by subroutine
C SRFC1. For their description refer to subroutine SRFC1 (included
C in the subroutine file 'srfc.for').
C (10) The data specifying the distribution of parameters of the model
C in all NCB material complex blocks. The data are read by
C subroutine PARM1. For their description refer to subroutine
C PARM1 (included in the subroutine file 'parm.for').
C For an example refer to the sample input data for the model.
C
C Storage in the memory:
C The input data (1) to (8) describing the structure of the model
C are stored in common blocks /MODELT/ and /MODELC/ defined in the
C following subroutine:
C ------------------------------------------------------------------
BLOCK DATA MODELB
INCLUDE 'model.inc'
END
C ------------------------------------------------------------------
C
C Date: 1996, September 30
C Coded by Ludek Klimes
C
C=======================================================================
C
SUBROUTINE MODEL1(LUN)
INTEGER LUN
C
C Subroutine MODEL1 reads the input data (1)-(8) for the description
C of the model and stores them in common blocks /MODELT/ and /MODELC/.
C Then it calls subroutine SRFC1 (included in the subroutine file
C 'srfc.for') to read the input data (9) for smooth surfaces and to
C store them in the memory. Finaly, it calls subroutine PARM1 (included
C in the subroutine file 'parm.for') to read the input data (10) for the
C parameters of the medium and to store them in the memory.
C
C Input:
C LUN... Logical unit number of the external input device
C containing the input data.
C The input parameter is not altered.
C
C No output.
C
C Common blocks /MODELT/ and /MODELC/:
INCLUDE 'model.inc'
C All the storage locations of the common blocks are defined in this
C subroutine.
C
C Subroutines and external functions required:
EXTERNAL MODELB,METR1,SRFC1,PARM1
C MODELB..Block data subroutine of this file.
C METR1...File 'metric.for'.
C SRFC1 and subsequent routines... File 'srfc.for' and subsequent
C files (like 'val.for' and 'fit.for').
C PARM1 and subsequent routines... File 'parm.for' and subsequent
C files (like 'val.for' and 'fit.for').
C
C Error messages:
C 310... Insufficient memory in /MODELC/:
C Insufficient memory for the input data in common block
C /MODELC/. The dimensions MSB or MCB of arrays KSB or KCB,
C respectively, must be enlarged. See the block data
C subroutine MODELB.
C 311... Block bounded by wrong interface:
C Index of the surface bounding the simple block is greater
C than the specified number of surfaces.
C 312... C. block composed of wrong s. block:
C Index of a simple block composing the complex block is
C greater than the specified number of simple blocks.
C
C Date: 1996, September 30
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C Auxiliary storage locations:
INTEGER I,J,L
C
READ(LUN,*) TEXTM
I=0
NEXPV=1
NEXPQ=1
IVERT=0
READ(LUN,*) I,NEXPV,NEXPQ,IVERT
CALL METR1(I)
READ(LUN,*) BOUNDM
READ(LUN,*) NSRFCS
C
C Simple blocks:
C Number of simple blocks
READ(LUN,*) NSB
C Initializing memory for indices of surfaces bounding simple blocks
L=NSB+1
DO 11 I=L,MSB
KSB(I)=0
11 CONTINUE
C Reading indices of surfaces bounding simple blocks:
DO 14 J=1,NSB
READ(LUN,*) (KSB(I),I=L,MSB)
DO 12 I=L,MSB
IF(IABS(KSB(I)).GT.NSRFCS) THEN
PAUSE 'Error 311 in MODEL: Block bounded by wrong interface'
STOP
ELSE IF(KSB(I).EQ.0) THEN
KSB(J)=I-1
L=I
GO TO 13
END IF
12 CONTINUE
GO TO 99
13 CONTINUE
14 CONTINUE
C
C Complex blocks:
C Number of complex blocks
READ(LUN,*) NCB
C Initializing memory for indices of simple blocks forming c. blocks
L=NCB+1
DO 21 I=L,MCB
KCB(I)=0
21 CONTINUE
C Reading indices of simple blocks forming complex blocks
DO 24 J=1,NCB
READ(LUN,*) (KCB(I),I=L,MCB)
DO 22 I=L,MCB
IF(KCB(I).LT.0.OR.KCB(I).GT.NSB) THEN
PAUSE
* 'Error 312 in MODEL: C. block composed of wrong s. block.'
STOP
ELSE IF(KCB(I).EQ.0) THEN
KCB(J)=I-1
L=I
GO TO 23
END IF
22 CONTINUE
GO TO 99
23 CONTINUE
24 CONTINUE
C
C Smooth surfaces:
CALL SRFC1(LUN,NSRFCS)
C
C Material parameters:
CALL PARM1(LUN,NCB)
RETURN
C
99 CONTINUE
PAUSE 'Error 310 in MODEL1: Insufficient memory in /MODELC/'
STOP
END
C
C=======================================================================
C
INTEGER FUNCTION NSRFC()
C
C Integer function NSRFC is designed to return the number of surfaces
C covering structural interfaces.
C
C No input.
C
C Output:
C NSRFC...Number of surfaces covering structural interfaces.
C
C Common block /MODELC/:
INCLUDE 'model.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 19
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C No auxiliary storage locations.
C
NSRFC=NSRFCS
RETURN
END
C
C=======================================================================
C
SUBROUTINE BLOCK(INSIDE,COOR,ISRF1,ISB1,ISRF2,ISB2,ICB2,F)
LOGICAL INSIDE
REAL COOR(3)
INTEGER ISRF1,ISB1,ISB2,ICB2,ISRF2
REAL F(10)
C
C This subroutine searches for the simple block and the complex block in
C which a given point is situated. This routine may be also used to
C determine the index of a block touching a specified block at a given
C point situated on the boundary of the specified block (the situation
C which may occur when a ray impinges on a boundary of a block).
C Another function of the routine is to determine the index of any of
C the surfaces bounding a block and separating the block from the given
C point.
C
C Input:
C INSIDE..INSIDE=.TRUE.: the given point is checked for its
C location inside simple block ISB1.
C INSIDE=.FALSE.: the given point is expected to be
C situated outside simple block ISB1, and simple block
C ISB2 (different from ISB1) in which the point is
C situated is being searched for. If ISRF1.NE.0, the
C position with respect to surface ISRF1 is not checked,
C and the simple blocks situated at the other side of
C ISRF1 than simple block ISB1 are preferred when
C searching for ISB2.
C COOR... Array containing coordinates X1, X2, X3 of a given point.
C ISRF1...ISRF1.NE.0: index of the surface at which the given point
C is situated. The sign is ignored.
C ISRF1.EQ.0: the given point is not situated at any surface
C ISB1... For ISRF1.NE.0: index of a simple block bounded by the
C surface IABS(ISRF1) - search for the index of a
C neighbouring simple block touching ISB1 at the given
C point.
C For ISRF1.EQ.0: index of an arbitrary simple block or
C ISB1=0.
C None of the input parameters are altered.
C
C Output:
C ISRF2...for the given point not situated inside the block ISB1 or
C at its boundary, ISRF2 has the meaning of the index of one
C of the surfaces bounding the simple block ISB1 and
C separating the given point from the simple block ISB1,
C supplemented by a sign '+' or '-' for the simple block
C ISB1 situated at the positive or negative side of the
C surface, respectively.
C Otherwise zero.
C If ISRF1.NE.0, the position with respect to surface ISRF1
C is not checked, and if INSIDE=.FALSE., ISRF2=0 without any
C checking.
C ISB2... For ISRF1.NE.0, ISRF2.EQ.0: index of a simple block
C neighbouring to the simple block isb1 and separated from
C the simple block ISB1 by surface IABS(ISRF1) at the
C given point.
C For ISRF1.NE.0, ISRF2.NE.0: index of a simple block in
C which the given point is situated. In this case, the
C given point may be situated either inside the simple
C block ISB2 or in the vicinity of its boundary formed by
C the surface IABS(ISRF1). From the two possible simple
C blocks touching the surface IABS(ISRF1) at the given
C point, that being situated at the same side of the
C surface ISRF1 as the simple block ISB1, is selected.
C For ISRF1.EQ.0: index of the simple block in which the
C given point is situated.
C If there is no simple material block of the above
C properties, ISB2=0.
C ICB2... Index of the complex block in which the simple block ISB2
C is situated. ICB2=0 if the simple block ISB2 is not
C situated in any complex block. For ISB2=0: ICB2=0.
C F... For ISRF2.NE.0: array containing the value and partial
C derivatives F, F1, F3, F11, F12, F22, F13, F23, F33 of
C the function describing the surface IABS(ISRF2).
C For ISRF2.EQ.0: undefined.
C
C Common block /MODELC/:
INCLUDE 'model.inc'
C None of the storage locations of the common block are altered.
C
C Subroutines and external functions required:
EXTERNAL ISIDE,INTERF,SRFC2
INTEGER ISIDE
C ISIDE,INTERF... This file.
C SRFC2 and subsequent routines... File 'srfc.for' and subsequent
C files (like 'val.for' and 'fit.for').
C
C Error messages:
C 313... Wrong index of surface:
C Absolute value of the input parameter ISRFC1 (index of the
C surface) is greater than the number NSRFC of the surfaces
C covering structural interfaces.
C 314... Wrong index of simple block:
C Parameter ISB1 (index of the simple block) is either
C negative or greater than the number nsb of simple blocks.
C
C Date: 1995, October 29
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C Auxiliary storage locations:
INTEGER I,J,I1,I2,ISIDE1,ISIDE2,ISRF
REAL FAUX(10)
C
C Checking input values:
IF(ISRF1.LT.-NSRFCS.OR.ISRF1.GT.NSRFCS) THEN
PAUSE 'Error 313 in BLOCK: Wrong index of surface'
STOP
END IF
IF(ISB1.LT.0.OR.ISB1.GT.NSB) THEN
PAUSE 'Error 314 in BLOCK: Wrong index of simple block'
STOP
END IF
C
C Initialization:
ISRF2=0
IF(ISRF1.EQ.0.OR.ISB1.EQ.0) THEN
ISIDE1=2
ELSE
ISIDE1=-ISIDE(ISRF1,ISB1)
END IF
C
C Position of the given point with respect to the given simple block
C ISB1:
IF(INSIDE) THEN
IF(ISB1.NE.0) THEN
CALL INTERF(COOR,ISRF1,ISB1,ISRF2,F)
IF(ISRF1.EQ.0) THEN
IF(ISRF2.EQ.0) THEN
C The point is inside simple block ISB1:
ISB2=ISB1
GO TO 10
END IF
ELSE
IF(ISRF2.NE.0) THEN
C The point being situated at surface ISRF1 bounding simple
C block ISB1 is not situated at the boundary of simple block
C ISB1:
ISIDE1=-ISIDE1
END IF
END IF
END IF
END IF
C
C Search for the simple block in which the given point is
C situated:
I2=MAX0(ISB1-1,NSB-ISB1)
DO 2 J=1,I2
DO 1 I=1,-1,-2
ISB2=ISB1+I*J
IF(ISB2.GT.0.AND.ISB2.LE.NSB) THEN
IF(ISRF1.NE.0) THEN
ISIDE2=ISIDE(ISRF1,ISB2)
END IF
IF(ISRF1.EQ.0.OR.ISIDE1.EQ.ISIDE2) THEN
CALL INTERF(COOR,ISRF1,ISB2,ISRF,FAUX)
IF(ISRF.EQ.0) GO TO 10
END IF
END IF
1 CONTINUE
2 CONTINUE
IF(.NOT.INSIDE) THEN
DO 4 J=1,I2
DO 3 I=1,-1,-2
ISB2=ISB1+I*J
IF(ISB2.GT.0.AND.ISB2.LE.NSB) THEN
CALL INTERF(COOR,ISRF1,ISB2,ISRF,FAUX)
IF(ISRF.EQ.0) GO TO 10
END IF
3 CONTINUE
4 CONTINUE
END IF
C No simple block has been found:
ISB2=0
ICB2=0
RETURN
C
C Determination of the complex block:
10 CONTINUE
DO 12 J=1,NCB
I1=KCB(J-1)+1
I2=KCB(J)
DO 11 I=I1,I2
IF(KCB(I).EQ.ISB2) THEN
ICB2=J
RETURN
END IF
11 CONTINUE
12 CONTINUE
C No complex block:
ICB2=0
C
RETURN
END
C
C=======================================================================
C
INTEGER FUNCTION ISIDE(ISRF,ISB)
INTEGER ISRF,ISB
C
C This is an auxiliary function to the subroutine BLOCK.
C This function determines the mutual position of a surface and a simple
C block.
C
C Input:
C ISRF... Index of the surface. The sign is ignored.
C ISB... Index of the simple block.
C None of the input parameters are altered.
C
C Output:
C ISIDE...ISIDE=-1: The simple block is bounded by the surface and
C is situated on its negative side.
C ISIDE= 1: The simple block is bounded by the surface and
C is situated on its positive side.
C ISIDE= 2: The simple block is not bounded by the surface.
C
C Common block /MODELC/:
INCLUDE 'model.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 19
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C Auxiliary storage locations:
INTEGER IS,LS,MS
C
LS=KSB(ISB-1)+1
MS=KSB(ISB)
C
C Loop for surfaces bounding simple block ISB:
DO 1 IS=LS,MS
IF(IABS(KSB(IS)).EQ.IABS(ISRF)) THEN
ISIDE=ISIGN(1,KSB(IS))
RETURN
END IF
1 CONTINUE
C
ISIDE=2
RETURN
END
C
C=======================================================================
C
SUBROUTINE INTERF(COOR,ISRF1,ISB,ISRF2,F)
REAL COOR(3),F(10)
INTEGER ISRF1,ISB,ISRF2
C
C This is an auxiliary subroutine to the subroutine BLOCK.
C This subroutine determines the position of a given point with respect
C to a given simple block.
C
C Input:
C COOR... Array containing coordinates X1, X2, X3 of a given point.
C ISRF1...ISRF1.NE.0: Index of the surface at which the given point
C is situated. The sign is ignored.
C ISRF1.EQ.0: The given point is not situated at any
C surface.
C ISB... Index of the given simple block.
C None of the input parameters are altered.
C
C Output:
C ISRF2...Index of a surface separating the given point and the
C simple block ISB, supplemented by a sign '+' or '-' for
C the simple block ISB1 situated at the positive or
C negative side of the surface, respectively.
C ISRF2=0 if the given point is situated inside the simple
C block ISB.
C F... For ISRF2.NE.0: Array containing the value and partial
C derivatives F, F1, F3, F11, F12, F22, F13, F23, F33 of
C the function describing the surface IABS(ISRF2) at the
C given point.
C For ISRF2.EQ.0: Undefined.
C
C Common block /MODELC/:
INCLUDE 'model.inc'
C None of the storage locations of the common block are altered.
C
C Subroutines and external functions required:
EXTERNAL SRFC2
C SRFC2 and subsequent routines... File 'srfc.for' and subsequent
C files (like 'val.for' and 'fit.for').
C
C Date: 1989, December 19
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C Auxiliary storage locations:
INTEGER IS,JS,KS,LS,MS
C
LS=KSB(ISB-1)+1
MS=KSB(ISB)
C
C Loop for surfaces bounding simple block ISB:
DO 1 IS=LS,MS
KS=KSB(IS)
JS=IABS(KS)
IF(JS.NE.IABS(ISRF1)) THEN
CALL SRFC2(JS,COOR,F)
IF(F(1)*FLOAT(KS).LT.0.) THEN
ISRF2=KS
RETURN
END IF
END IF
1 CONTINUE
C
ISRF2=0
RETURN
END
C
C=======================================================================
C
SUBROUTINE VELOC(IWAVE,UP,US,QP,QS,VP,VS,VD,QL)
INTEGER IWAVE
REAL UP(10),US(10),QP,QS,VP,VS,VD(10),QL
C
C This subroutine transforms the values of parameters of the medium into
C velocities and loss factors.
C
C Input:
C IWAVE...Type of wave.
C IWAVE.GE.0: P wave,
C IWAVE.LT.0: S wave.
C UP,US...Powers of P and S wave velocities and their first and
C second partial derivatives (the exponent of the powers is
C NEXPV, see 'Input data for the model'), in order U, U1,
C U2, U3, U11, U12, U22, U13, U23, U33.
C QP,QS...Powers of the loss factors of P and S waves (the exponent
C of the powers is NEXPQ, see 'Input data for the model').
C None of the input parameters are altered.
C
C Output:
C VP,VS...P and S wave velocities.
C VD... Velocity and its first and second partial derivatives
C ordered as UP, US, corresponding to the wave specified by
C IWAVE, in order V, V1, V2, V3, V11, V12, V22, V13, V23,
C V33.
C QL... Loss factor corresponding to the wave specified by IWAVE.
C
C Common block /MODELC/:
INCLUDE 'model.inc'
C None of the storage locations of the common block are altered.
C
C Subroutines and external functions required:
EXTERNAL FPOWER
C FPOWER...This file.
C
C Date: 1992, December 31
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C Auxiliary storage location:
REAL POWER,AUX1(1),AUX2(1)
C
POWER=FLOAT(NEXPV)
IF(IWAVE.GE.0) THEN
CALL FPOWER(10,UP,POWER,VD)
*V CALL VAR5(1,1)
VP=VD(1)
CALL FPOWER(1,US,POWER,AUX2)
VS=AUX2(1)
AUX1(1)=QP
ELSE
CALL FPOWER(1,UP,POWER,AUX2)
VP=AUX2(1)
CALL FPOWER(10,US,POWER,VD)
*V CALL VAR5(2,2)
VS=VD(1)
AUX1(1)=QS
END IF
CALL FPOWER(1,AUX1,FLOAT(NEXPQ),AUX2)
QL=AUX2(1)
RETURN
END
C
C=======================================================================
C
SUBROUTINE FPOWER(N,FINP,POWER,FOUT)
INTEGER N
REAL FINP(N),POWER,FOUT(N)
C
C This subroutine evaluates the value and, possibly, the three first and
C six second partial derivatives of a function if the value and the
C three first and six second partial derivatives of the POWER-th power
C of the function are known.
C
C Input:
C N... For N=1: only the function value is evaluated. The
C derivatives are ignored.
C For N=4: the value and the three first partial derivatives
C are evaluated.
C For N=10: the value and the three first and six second
C partial derivatives are evaluated.
C FINP... Array containing the value, the first and second partial
C derivatives of the POWER-th power of the function to be
C evaluated, in the order F, F1, F2, F3, F11, F12, F22, F13,
C F23, F33. For N=1, only the function value is required.
C POWER...The specified function is equal to the POWER-th power of
C the corresponding physical quantity.
C POWER=0: Zero output array FOUT is generated.
C None of the input parameters are altered (except FINP if this
C parameter and FOUT are identical in the calling sequence).
C
C Output:
C FOUT... Array containing the value, the first and second partial
C derivatives of the evaluated function, in the order F, F1,
C F2, F3, F11, F12, F22, F13, F23, F33. This parameter may
C coincide with FINP, in which case FINP is destroyed on
C output. Note that this coincidence is an exception to
C ANSI standard of FORTRAN 77.
C
C No subroutines and external functions required.
C
C Date: 1996, September 30
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C Auxiliary storage locations:
INTEGER I
REAL F,AUX1,AUX2
C
IF(POWER.EQ.0.) THEN
DO 1 I=1,N
FOUT(I)=0.
1 CONTINUE
ELSE IF(0.999.LT.POWER.AND.POWER.LT.1.001) THEN
DO 2 I=1,N
FOUT(I)=FINP(I)
2 CONTINUE
*V CALL VAR4(0,1.)
ELSE
IF(-1.001.LT.POWER.AND.POWER.LT.-0.999) THEN
F=1./FINP(1)
ELSE
F=FINP(1)**(1./POWER)
END IF
FOUT(1)=F
IF(N.GT.1) THEN
AUX1= F/(FINP(1)*POWER)
AUX2= (POWER-1.)/F
FOUT(2)=AUX1*FINP(2)
FOUT(3)=AUX1*FINP(3)
FOUT(4)=AUX1*FINP(4)
IF(N.GT.4) THEN
FOUT(5)=AUX1*FINP(5)-AUX2*FOUT(2)*FOUT(2)
FOUT(6)=AUX1*FINP(6)-AUX2*FOUT(2)*FOUT(3)
FOUT(7)=AUX1*FINP(7)-AUX2*FOUT(3)*FOUT(3)
FOUT(8)=AUX1*FINP(8)-AUX2*FOUT(2)*FOUT(4)
FOUT(9)=AUX1*FINP(9)-AUX2*FOUT(3)*FOUT(4)
FOUT(10)=AUX1*FINP(10)-AUX2*FOUT(4)*FOUT(4)
END IF
*V CALL VAR4(0,AUX1)
*V CALL VAR4(2,-AUX2*FOUT(2))
*V CALL VAR4(3,-AUX2*FOUT(3))
*V CALL VAR4(4,-AUX2*FOUT(4))
END IF
END IF
RETURN
END
C
C=======================================================================
C