C Subroutine file 'model.for' to specify a seismic model. C C Date: 1997, October 25 C Coded by Ludek Klimes C C....................................................................... C C This file consists of the following external procedures: 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 finally 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 MODEL1 C NSRFC...Integer function returning the number of surfaces covering C structural interfaces. C NSRFC C BLOCK...Subroutine designed to determine the mutual position of a C point and a simple and a complex block. C BLOCK C ISIDE...Auxiliary function to subroutine BLOCK. C ISIDE C INTERF..Auxiliary subroutine to subroutine BLOCK. C INTERF C VELOC...Subroutine transforming the values of a medium parameters C into velocities and loss factors. C VELOC 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 FPOWER C C Note: C The lines denoted by '*V' in the first two columns of file C 'model.for' in subroutines VELOC and FPOWER are designed to C calculate the model variations with respect to the model C parameters. C File 'modelv.for', intended for the model inversion, is created C from 'model.for' by replacing each '*V' in the first two columns C by spaces using program 'clean.for'. Subroutines VAR4 and VAR5 C of file 'var.for' may then be called to handle the variations. C C....................................................................... C C C Input data MODEL 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 coordinate system: C KOORS.LE.0: Cartesian coordinates (default). C KOORS.EQ.1: polar spherical coordinates in radians, C (X1,X2,X3)=(colatitude,longitude,radius). C KOORS.GE.2: geographic spherical coordinates in radians, C (X1,X2,X3)=(longitude,latitude,radius). C If the coordinate system is right-handed (recommended), C all vectorial products are evaluated using the right-hand C rule, otherwise using the left-hand rule. C KOORS is passed to the subroutines of file 'metric.for' C by means of invocation of subroutine METR1 and presently C represents the only input data for the coordinate system. C Note that possible future additional data for the C coordinate system should be read by subroutine METR1 and C should be located between input data (2) and (3). C metric.for C NEXPV,NEXPQ... Specify exponents of the power of velocities C (NEXPV) and loss factors (NEXPQ) reported by the C subroutines evaluating isotropic material parameters. 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 When using the basic submitted version of the subroutine C file 'parm.for', the default values of NEXPV=1, NEXPQ=1 C are highly recommended. Other values make sense only if a C user is submitting his own subroutines to evaluate C isotropic material parameters which, e.g., evaluate the C slowness instead of the velocity. In such a case, C switching NEXPV from 1 to -1 may avoid the modification C of the user's subroutines. 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 It is recommended to define only surfaces covering structural C interfaces (model surfaces) in this data set. Auxiliary surfaces C related to particular source-receiver configurations, numerical C procedures, etc., should preferably be defined in other data sets. 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). Intersecting simple blocks are allowed but not recommended. C All material simple blocks in the model must be defined. C Free-space blocks need not be (and usually are not) defined in C this data set. Free-space blocks may be defined at user's C discretion in order to enhance the possibility to check for the C model consistency. Free-space blocks may be defined either here C or in the additional data file designed just for the consistency C check by program MODCHK. C Program MODCHK 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 srfc.for: Input data 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 parm.for: Input data C For an example refer to the sample input data for the model. C Example of data set MODEL C C....................................................................... 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 include file 'model.inc'. C model.inc C C======================================================================= 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. Finally, it calls subroutine PARM1 C (included in the subroutine file 'parm.for') to read the input data C (10) for the 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 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 METR1,SRFC1,PARM1 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 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 C 311 PAUSE 'Error 311 in MODEL: Block bounded by wrong interface' STOP C Index of the surface bounding the simple block is greater C than the specified number of surfaces. 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 C 312 PAUSE * 'Error 312 in MODEL: C. block composed of wrong s. block.' STOP C Complex block composed of wrong simple block: C Index of a simple block composing the complex block is C greater than the specified number of simple blocks. 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 C 310 PAUSE 'Error 310 in MODEL1: Insufficient memory in /MODELC/' STOP C Insufficient memory for the input data in common block /MODELC/. C The dimensions MSB or MCB of arrays KSB or KCB, respectively, C must be enlarged. C Refer to include file model.inc END C C======================================================================= 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 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 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 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 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 C 313 PAUSE 'Error 313 in BLOCK: Wrong index of surface' STOP 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. END IF IF(ISB1.LT.0.OR.ISB1.GT.NSB) THEN C 314 PAUSE 'Error 314 in BLOCK: Wrong index of simple block' STOP C Parameter ISB1 (index of the simple block) is either C negative or greater than the number nsb of simple blocks. 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 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 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 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 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 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 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) 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) 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 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 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 CALL VAR4(0,AUX1) CALL VAR4(2,-AUX2*FOUT(2)) CALL VAR4(3,-AUX2*FOUT(3)) CALL VAR4(4,-AUX2*FOUT(4)) END IF END IF RETURN END C C======================================================================= C