[cig-commits] r12416 - mc/2D/ConMan/trunk/src
wei at geodynamics.org
wei at geodynamics.org
Tue Jul 15 08:54:34 PDT 2008
Author: wei
Date: 2008-07-15 08:54:31 -0700 (Tue, 15 Jul 2008)
New Revision: 12416
Added:
mc/2D/ConMan/trunk/src/fftsubs.f
mc/2D/ConMan/trunk/src/geoid.F
Modified:
mc/2D/ConMan/trunk/src/ConMan.F
mc/2D/ConMan/trunk/src/Makefile
mc/2D/ConMan/trunk/src/common.h
mc/2D/ConMan/trunk/src/eg2.F
mc/2D/ConMan/trunk/src/stress.F
Log:
Fixed the stress.f routine problem. It was converge slowly(first order rather than second order). Scott added the geoid routine, now the runfiles will need a ninth filename(geoid.xxxxx).
Modified: mc/2D/ConMan/trunk/src/ConMan.F
===================================================================
--- mc/2D/ConMan/trunk/src/ConMan.F 2008-07-14 19:28:54 UTC (rev 12415)
+++ mc/2D/ConMan/trunk/src/ConMan.F 2008-07-15 15:54:31 UTC (rev 12416)
@@ -25,7 +25,7 @@
c
c.... remove above card for single-precision operation
c
- character*80 name1,name2,name3,name4,name6,name7,name8
+ character*80 name1,name2,name3,name4,name6,name7,name8,name9
common /temp1 / x(50000)
include 'common.h'
c
@@ -39,6 +39,7 @@
read(5,1001) name6
read(5,1001) name7
read(5,1001) name8
+ read(5,1001) name9
1001 format( a80 )
open (unit=iin, file=name1 , status='old')
open (unit=igeom, file=name2 , status='old')
@@ -48,6 +49,7 @@
open (unit=itsout, file=name6 , status='unknown' )
open (unit=itout, file=name7 , status='unknown')
open (unit=imout, file=name8 , status='unknown' )
+ open (unit=igeoid, file=name9 , status='unknown' )
c
c.... Input phase
c
@@ -99,7 +101,7 @@
c
c.... io channels
common /io / iin , igeom, iout , itsout , itout ,
- & imout , irsin , irsout
+ & imout , irsin , irsout, igeoid
c
c.... useful constants
common /const / zero , pt25, pt33, pt5, pt66, one,
@@ -113,9 +115,9 @@
c & / 1,130000000, 2, 1, 0 /
c
data iin , igeom , iout , itsout , itout , imout ,
- & irsin , irsout
+ & irsin , irsout , igeoid
& / 10, 11 , 12, 13, 14, 15,
- & 16, 17 /
+ & 16, 17 , 18/
c
data zero, pt25,
& pt33, pt5,
@@ -131,4 +133,4 @@
& 1.60000000000000d+1, 1.00000000000000d-7 /
c
c----------------------------------------------------------------------
- end
+ end
\ No newline at end of file
Modified: mc/2D/ConMan/trunk/src/Makefile
===================================================================
--- mc/2D/ConMan/trunk/src/Makefile 2008-07-14 19:28:54 UTC (rev 12415)
+++ mc/2D/ConMan/trunk/src/Makefile 2008-07-15 15:54:31 UTC (rev 12416)
@@ -9,7 +9,7 @@
vadd.o mydate.o mytime.o projct.o prj4q.o print_reg.o prtstr_reg.o \
batchelor.o
-SRC_HOME=/home/scott/ConMan02
+SRC_HOME=$HOME/conman-svn
IMPLICIT=f_tlhs.o f_trhsimp.o
PICARD=f_tlhs.o f_trhsimp.o
Modified: mc/2D/ConMan/trunk/src/common.h
===================================================================
--- mc/2D/ConMan/trunk/src/common.h 2008-07-14 19:28:54 UTC (rev 12415)
+++ mc/2D/ConMan/trunk/src/common.h 2008-07-15 15:54:31 UTC (rev 12416)
@@ -31,7 +31,7 @@
& neqt , neqv , itflag
c
common /io / iin, igeom, iout , itsout , itout , imout,
- & irsin , irsout
+ & irsin , irsout, igeoid
common /ioc / name5
c
common /title / ititle(80)
@@ -46,4 +46,4 @@
c
common /fswtch/ FACSTW
common /timer1/ cpu ,cpunew, cpuold, icode
- parameter (ireg = 0)
+ parameter (ireg = 0)
\ No newline at end of file
Modified: mc/2D/ConMan/trunk/src/eg2.F
===================================================================
--- mc/2D/ConMan/trunk/src/eg2.F 2008-07-14 19:28:54 UTC (rev 12415)
+++ mc/2D/ConMan/trunk/src/eg2.F 2008-07-15 15:54:31 UTC (rev 12416)
@@ -134,6 +134,7 @@
if (task .eq. 'prt_str ') then
call stres ( x , v , t )
call prtstr( x , stress )
+ call geoid ( x, t, stress)
c
c trhs is being used as a temp variable for the pressure at the node
c
@@ -165,4 +166,4 @@
call error ('EG2 ', task, 0)
return
c
- end
+ end
\ No newline at end of file
Added: mc/2D/ConMan/trunk/src/fftsubs.f
===================================================================
--- mc/2D/ConMan/trunk/src/fftsubs.f (rev 0)
+++ mc/2D/ConMan/trunk/src/fftsubs.f 2008-07-15 15:54:31 UTC (rev 12416)
@@ -0,0 +1,1243 @@
+c SUBROUTINE FACTOR(N)
+C FACTOR DETERMINES THE SMALLEST EVEN NUMBER .GE. N THAT CAN BE FACTORED INTO
+C PRODUCTS OF PRIMES IN THE ARRAY NP. THE RESULT IS RETURNED AS N.
+c ONLY THE FIRST FOUR PRIMES ARE USED. THIS ROUTINE IS CALLED BEFORE USING FFTL.
+c DIMENSION NP(4)
+c DATA NP/2,3,5,7/
+c NCHK=N/2
+c IF(N.NE.(NCHK*2)) N=N+1
+c 5 NF=N
+c DO 15 I=1,4
+c M=NP(I)
+c 10 L=MOD(NF,M)
+c IF(L.NE.0) GO TO 15
+c NF=NF/M
+c GO TO 10
+c 15 CONTINUE
+c IF(NF.EQ.1) RETURN
+c N=N+2
+c GO TO 5
+c END
+ SUBROUTINE FACDWN(N)
+C FACDWN DETERMINES THE LARGEST EVEN NUMBER .LE. N THAT CAN BE FACTORED INTO
+C PRODUCTS OF PRIMES IN THE ARRAY NP. THE RESULT IS RETURNED AS N.
+C ONLY THE FIRST FOUR PRIMES ARE USED. THIS ROUTINE IS CALLED BEFORE USING FFTL.
+ DIMENSION NP(4)
+ DATA NP/2,3,5,7/
+ NCHK=N/2
+ IF(N.NE.(NCHK*2)) N=N-1
+ 5 NF=N
+ DO 15 I=1,4
+ M=NP(I)
+ 10 L=MOD(NF,M)
+ IF(L.NE.0) GO TO 15
+ NF=NF/M
+ GO TO 10
+ 15 CONTINUE
+ IF(NF.EQ.1) RETURN
+ N=N-2
+ GO TO 5
+ END
+ SUBROUTINE FFTL(X,N,NDIR,IERR)
+C IF IABS(NDIR)=1 FFTL FOURIER TRANSFORMS THE N POINTS REAL TIME SERIES IN ARRAY
+C X. THE RESULT OVERWRITES X STORED AS (N+2)/2 COMPLEX NUMBERS (NON-NEGATIVE
+C FREQUENCIES ONLY).
+C IF IABS(NDIR)=2 FFTL FOURIER TRANSFORMS THE (N+2)/2 COMPLEX FOURIER COEFFICIENTS
+C (NON-NEGATIVE FREQUENCIES ONLY) IN ARRAY X (ASSUMING THE SERIES IS HERMITIAN).
+C THE RESULTING N POINT REAL TIME SERIES OVERWRITES X.
+C NDIR>0 THE FORWARD TRANSFORM USES THE SIGN CONVENTION EXP(I*W*T).
+C NDIR<0 THE FORWARD TRANSFORM USES THE SIGN CONVENTION EXP(-I*W*T).
+C THE FORWARD TRANSFORM IS NORMALIZED SUCH THAT A SINE WAVE OF UNIT AMPLITUDE
+C IS TRANSFORMED INTO DELTA FUNCTIONS OF UNIT AMPLITUDE. THE BACKWARD TRANSFORM
+C IS NORMALIZED SUCH THAT TRANSFORMING FORWARD AND THEN BACK RECOVERS THE
+C ORIGINAL SERIES. IERR IS NORMALLY ZERO. IF IERR.EQ.1 THEN FFT HAS NOT BEEN
+C ABLE TO FACTOR THE SERIES. HOWEVER, X HAS BEEN SCRAMBLED BY REALTR. NOTE
+C THAT IF N IS ODD THE LAST POINT WILL NOT BE USED IN THE TRANSFORM.
+C -RPB
+ DIMENSION X(*)
+ N2=N/2
+ IDIR=IABS(NDIR)
+ GO TO (1,2),IDIR
+C DO FORWARD TRANSFORM (IE. TIME TO FREQUENCY).
+ 1 CALL FFT(X,X(2),N2,N2,N2,2,IERR)
+ CALL REALTR(X,X(2),N2,2)
+ N1=2*N2+2
+ SCALE=1./N
+ IF(NDIR.GT.0) GO TO 3
+ DO 5 I=4,N,2
+ 5 X(I)=-X(I)
+ GO TO 3
+C DO BACKWARD TRANSFORM (IE. FREQUENCY TO TIME).
+ 2 IF(NDIR.GT.0) GO TO 6
+ DO 7 I=4,N,2
+ 7 X(I)=-X(I)
+ 6 X(2)=0.
+ X(2*N2+2)=0.
+ CALL REALTR(X,X(2),N2,-2)
+ CALL FFT(X,X(2),N2,N2,N2,-2,IERR)
+ N1=2*N2
+ SCALE=.5
+C NORMALIZE THE TRANSFORM.
+ 3 DO 4 I=1,N1
+ 4 X(I)=SCALE*X(I)
+ RETURN
+ END
+ SUBROUTINE REALTR(A,B,N,ISN)
+C IF ISN=1, THIS SUBROUTINE COMPLETES THE FOURIER TRANSFORM OF 2*N REAL
+C DATA VALUES, WHERE THE ORIGINAL DATA VALUES ARE STORED ALTERNATELY IN
+C ARRAYS A AND B, AND ARE FIRST TRANSFORMED BY A COMPLEX FOURIER TRANSFORM
+c OF DIMENSION N. THE COSINE COEFFICIENTS ARE IN A(1),A(2),...A(N+1) AND
+C THE SINE COEFFICIENTS ARE IN B(1),B(2),...B(N+1). A TYPICAL CALLING
+C SEQUENCE IS
+C CALL FFT (A,B,N,1)
+C CALL REALTR (A,B,N,1)
+C THE RESULTS SHOULD BE MULTIPLIED BY 0.5/N TO GIVE THE USUAL SCALING
+C IF ISN=1, THE INVERSE TRANSFORMATION IS DONE, THE FIRST STEP IN EVALUATING
+C A REAL FOURIER SERIES. A TYPICAL CALLING SEQUENCE IS
+C CALL REALTR (A,B,N,-1)
+C CALL FFT (A,B,N,-1)
+C THE RESULTS SHOULD BE MULTIPLIED BY 0.5 TO GIVE THE USUAL SCALING, AND
+C THE TIME DOMAIN RESULTS ALTERNATE IN ARRAYS A AND B, I.E. A(1),B(1), A(2),B(2),
+C ...A(N),B(N). THE DATA MAY ALTERNATIVELY BE STORED IN A SINGLE COMPLEX
+C ARRAY A, THEN THE MAGNITUDE OF ISN CHANGED TO TWO TO GIVE THE CORRECT
+C INDEXING INCREMENT AND A(2) USED TO PASS THE INITIAL ADDRESS FOR THE
+C SEQUENCE OF IMAGINARY VALUES,E.G.
+C CALL FFT(A,A(2),N,2)
+C CALL REALTR(A,A(2),N,2)
+C IN THIS CASE, THE COSINE AND SINE COEFFICIENTS ALTERNATE IN A.
+ DIMENSION A(*),B(*)
+ IF(N .LE. 1) RETURN
+ INC=ISN
+ IF(INC.LT.0) INC=-INC
+ NK = N * INC + 2
+ NH = NK / 2
+ SD = 3.1415926535898E0/(2.E0*N)
+ CD = 2.E0 * SIN(SD)**2
+ SD = SIN(SD+SD)
+ SN = 0.E0
+ IF(ISN .GT. 0) GO TO 10
+ CN = -1.E0
+ SD = -SD
+ GO TO 20
+ 10 CN = 1.E0
+ A(NK-1) = A(1)
+ B(NK-1) = B(1)
+ 20 DO 30 J=1,NH,INC
+ K = NK - J
+ AA = A(J) + A(K)
+ AB = A(J) - A(K)
+ BA = B(J) + B(K)
+ BB = B(J) - B(K)
+ XX = CN * BA + SN * AB
+ YY = SN * BA - CN * AB
+ B(K) = YY - BB
+ B(J) = YY + BB
+ A(K) = AA - XX
+ A(J) = AA + XX
+ AA = CN - (CD * CN + SD * SN)
+ SN = (SD * CN - CD * SN) + SN
+ 30 CN = AA
+ RETURN
+ END
+ SUBROUTINE FFT(A,B,NTOT,N,NSPAN,ISN,IERR)
+C FFT IS A SINGLE PRECISION VERSION OF SINGLETON'S FFT PROCEDURE.
+C MULTIVARIATE COMPLEX FOURIER TRANSFORM, COMPUTED IN PLACE USING
+C MIXED-RADIX FAST FOURIER TRANSFORM ALGORITHM. BY R. C. SINGLETON,
+C STANFORD RESEARCH INSTITUTE, SEPT. 1968
+C ARRAYS A AND B ORIGINALLY HOLD THE REAL AND IMAGINARY COMPONENTS OF THE
+C DATA, AND RETURN THE REAL AND IMAGINARY COMPONENTS OF THE RESULTING FOURIER
+C COEFFICIENTS. MULTIVARIATE DATA IS INDEXED ACCORDING TO THE FORTRAN ARRAY
+C ELEMENT SUCCESSOR FUNCTION, WITHOUT LIMIT ON THE NUMBER OF IMPLIED MULTIPLE
+C SUBSCRIPTS. THE SUBROUTINE IS CALLED ONCE FOR EACH VARIATE. THE CALLS FOR A
+C MULTIVARIATE TRANSFORM MAY BE IN ANY ORDER.
+C NTOT IS THE TOTAL NUMBER OF COMPLEX DATA VALUES.
+C N IS THE DIMENSION OF THE CURRENT VARIABLE.
+C NSPAN/N IS THE SPACING OF CONSECUTIVE DATA VALUES FOR THE CURRENT VARIABLE.
+C IERR IS AN ERROR RETURN INDICATOR. IT IS NORMALLY ZERO, BUT IS SET TO 1 IF
+C THE NUMBER OF TERMS CANNOT BE FACTORED IN THE SPACE AVAILABLE.
+C THE SIGN OF ISN DETERMINES THE SIGN OF THE COMPLEX EXPONENTIAL, AND THE
+C MAGNITUDE OF ISN IS NORMALLY ONE.
+C A TRI-VARIATE TRANSFORM WITH A(N1,N2,N3), B(N1,N2,N3) IS COMPUTED BY
+C CALL FFT(A,B,N1*N2*N3,N1,N1,1)
+C CALL FFT(A,B,N1*N2*N3,N2,N1*N2,1)
+C CALL FFT(A,B,N1*N2*N3,N3,N1*N2*N3,1)
+C FOR A SINGLE-VARIATE TRANSFORM,
+C NTOT = N = NSPAN = (NUMBER OF COMPLEX DATA VALUES), E.G.
+C CALL FFT(A,B,N,N,N,1)
+C WITH MOST FORTRAN COMPILERS THE DATA CAN ALTERNATIVELY BE STORED IN A SINGLE
+C COMPLEX ARRAY A, THEN THE MAGNITUDE OF ISN CHANGED TO TWO TO GIVE THE CORRECT
+C INDEXING INCREMENT AND A(2) USED TO PASS THE INITIAL ADDRESS FOR THE SEQUENCE
+C OF IMAGINARY VALUES, E.G.
+C CALL FFT(A,A(2),NTOT,N,NSPAN,2)
+C ARRAYS AT(MAXF), CK(MAXF), BT(MAXF), SK(MAXF), AND NP(MAXP) ARE USED FOR
+C TEMPORARY STORAGE. IF THE AVAILABLE STORAGE IS INSUFFICIENT, THE PROGRAM IS
+C TERMINATED BY THE ERROR RETURN OPTION
+C MAXF MUST BE .GE. THE MAXIMUM PRIME FACTOR OF N.
+C MAXP MUST BE .GT. THE NUMBER OF PRIME FACTORS OF N.
+C IN ADDITION, IF THE SQUARE-FREE PORTION K OF N HAS TWO OR MORE PRIME FACTORS,
+C THEN MAXP MUST BE .GE. K-1. ARRAY STORAGE IN NFAC FOR A MAXIMUM OF 11 PRIME
+C FACTORS OF N. IF N HAS MORE THAN ONE SQUARE-FREE FACTOR, THE PRODUCT OF THE
+C SQUARE-FREE FACTORS MUST BE .LE. 210. (2**5*7=210)
+ DIMENSION A(*),B(*)
+ DIMENSION NFAC(11),NP(209)
+C ARRAY STORAGE FOR MAXIMUM PRIME FACTOR OF 23
+ DIMENSION AT(23),CK(23),BT(23),SK(23)
+ EQUIVALENCE (I,II)
+C THE FOLLOWING TWO CONSTANTS SHOULD AGREE WITH THE ARRAY DIMENSIONS.
+ MAXF=23
+ MAXP=209
+ IERR=0
+ IF(N .LT. 2) RETURN
+ INC=ISN
+ C72=0.30901699437494E0
+ S72=0.95105651629515E0
+ S120=0.86602540378443E0
+ RAD=6.2831853071796E0
+ IF(ISN .GE. 0) GO TO 10
+ S72=-S72
+ S120=-S120
+ RAD=-RAD
+ INC=-INC
+ 10 NT=INC*NTOT
+ KS=INC*NSPAN
+ KSPAN=KS
+ NN=NT-INC
+ JC=KS/N
+ RADF=RAD*(JC*0.5E0)
+ I=0
+ JF=0
+C DETERMINE THE FACTORS OF N
+ M=0
+ K=N
+ GO TO 20
+ 15 M=M+1
+ NFAC(M)=4
+ K=K/16
+ 20 IF(K-(K/16)*16 .EQ. 0) GO TO 15
+ J=3
+ JJ=9
+ GO TO 30
+ 25 M=M+1
+ NFAC(M)=J
+ K=K/JJ
+ 30 IF(MOD(K,JJ) .EQ. 0) GO TO 25
+ J=J+2
+ JJ=J**2
+ IF(JJ .LE. K) GO TO 30
+ IF(K .GT. 4) GO TO 40
+ KT=M
+ NFAC(M+1)=K
+ IF(K .NE. 1) M=M+1
+ GO TO 80
+ 40 IF(K-(K/4)*4 .NE. 0) GO TO 50
+ M=M+1
+ NFAC(M)=2
+ K=K/4
+ 50 KT=M
+ J=2
+ 60 IF(MOD(K,J) .NE. 0) GO TO 70
+ M=M+1
+ NFAC(M)=J
+ K=K/J
+ 70 J=((J+1)/2)*2+1
+ IF(J .LE. K) GO TO 60
+ 80 IF(KT .EQ. 0) GO TO 100
+ J=KT
+ 90 M=M+1
+ NFAC(M)=NFAC(J)
+ J=J-1
+ IF(J .NE. 0) GO TO 90
+C COMPUTE FOURIER TRANSFORM
+ 100 SD=RADF/(KSPAN)
+ CD=2.E0*SIN(SD)**2
+ SD=SIN(SD+SD)
+ KK=1
+ I=I+1
+ IF(NFAC(I) .NE. 2) GO TO 400
+C TRANSFORM FOR FACTOR OF 2 (INCLUDING ROTATION FACTOR)
+ KSPAN=KSPAN/2
+ K1=KSPAN+2
+ 210 K2=KK+KSPAN
+ AK=A(K2)
+ BK=B(K2)
+ A(K2)=A(KK)-AK
+ B(K2)=B(KK)-BK
+ A(KK)=A(KK)+AK
+ B(KK)=B(KK)+BK
+ KK=K2+KSPAN
+ IF(KK .LE. NN) GO TO 210
+ KK=KK-NN
+ IF(KK .LE. JC) GO TO 210
+ IF(KK .GT. KSPAN) GO TO 800
+ 220 C1=1.E0-CD
+ S1=SD
+ 230 K2=KK+KSPAN
+ AK=A(KK)-A(K2)
+ BK=B(KK)-B(K2)
+ A(KK)=A(KK)+A(K2)
+ B(KK)=B(KK)+B(K2)
+ A(K2)=C1*AK-S1*BK
+ B(K2)=S1*AK+C1*BK
+ KK=K2+KSPAN
+ IF(KK .LT. NT) GO TO 230
+ K2=KK-NT
+ C1=-C1
+ KK=K1-K2
+ IF(KK .GT. K2) GO TO 230
+ AK=CD*C1+SD*S1
+ S1=(SD*C1-CD*S1)+S1
+ C1=C1-AK
+ KK=KK+JC
+ IF(KK .LT. K2) GO TO 230
+ K1=K1+INC+INC
+ KK=(K1-KSPAN)/2+JC
+ IF(KK .LE. JC+JC) GO TO 220
+ GO TO 100
+C TRANSFORM FOR FACTOR OF 3 (OPTIONAL CODE)
+ 320 K1=KK+KSPAN
+ K2=K1+KSPAN
+ AK=A(KK)
+ BK=B(KK)
+ AJ=A(K1)+A(K2)
+ BJ=B(K1)+B(K2)
+ A(KK)=AK+AJ
+ B(KK)=BK+BJ
+ AK=-0.5*AJ+AK
+ BK=-0.5*BJ+BK
+ AJ=(A(K1)-A(K2))*S120
+ BJ=(B(K1)-B(K2))*S120
+ A(K1)=AK-BJ
+ B(K1)=BK+AJ
+ A(K2)=AK+BJ
+ B(K2)=BK-AJ
+ KK=K2+KSPAN
+ IF(KK .LT. NN) GO TO 320
+ KK=KK-NN
+ IF(KK .LE. KSPAN) GO TO 320
+ GO TO 700
+C TRANSFORM FOR FACTOR OF 4
+ 400 IF(NFAC(I) .NE. 4) GO TO 600
+ KSPNN=KSPAN
+ KSPAN=KSPAN/4
+ 410 C1=1.E0
+ S1=0
+ 420 K1=KK+KSPAN
+ K2=K1+KSPAN
+ K3=K2+KSPAN
+ AKP=A(KK)+A(K2)
+ AKM=A(KK)-A(K2)
+ AJP=A(K1)+A(K3)
+ AJM=A(K1)-A(K3)
+ A(KK)=AKP+AJP
+ AJP=AKP-AJP
+ BKP=B(KK)+B(K2)
+ BKM=B(KK)-B(K2)
+ BJP=B(K1)+B(K3)
+ BJM=B(K1)-B(K3)
+ B(KK)=BKP+BJP
+ BJP=BKP-BJP
+ IF(ISN .LT. 0) GO TO 450
+ AKP=AKM-BJM
+ AKM=AKM+BJM
+ BKP=BKM+AJM
+ BKM=BKM-AJM
+ IF(S1 .EQ. 0) GO TO 460
+ 430 A(K1)=AKP*C1-BKP*S1
+ B(K1)=AKP*S1+BKP*C1
+ A(K2)=AJP*C2-BJP*S2
+ B(K2)=AJP*S2+BJP*C2
+ A(K3)=AKM*C3-BKM*S3
+ B(K3)=AKM*S3+BKM*C3
+ KK=K3+KSPAN
+ IF(KK .LE. NT) GO TO 420
+ 440 C2=CD*C1+SD*S1
+ S1=(SD*C1-CD*S1)+S1
+ C1=C1-C2
+ C2=C1**2-S1**2
+ S2=2.E0*C1*S1
+ C3=C2*C1-S2*S1
+ S3=C2*S1+S2*C1
+ KK=KK-NT+JC
+ IF(KK .LE. KSPAN) GO TO 420
+ KK=KK-KSPAN+INC
+ IF(KK .LE. JC) GO TO 410
+ IF(KSPAN .EQ. JC) GO TO 800
+ GO TO 100
+ 450 AKP=AKM+BJM
+ AKM=AKM-BJM
+ BKP=BKM-AJM
+ BKM=BKM+AJM
+ IF(S1 .NE. 0) GO TO 430
+ 460 A(K1)=AKP
+ B(K1)=BKP
+ A(K2)=AJP
+ B(K2)=BJP
+ A(K3)=AKM
+ B(K3)=BKM
+ KK=K3+KSPAN
+ IF(KK .LE. NT) GO TO 420
+ GO TO 440
+C TRANSFORM FOR FACTOR OF 5 (OPTIONAL CODE)
+ 510 C2=C72**2-S72**2
+ S2=2.E0*C72*S72
+ 520 K1=KK+KSPAN
+ K2=K1+KSPAN
+ K3=K2+KSPAN
+ K4=K3+KSPAN
+ AKP=A(K1)+A(K4)
+ AKM=A(K1)-A(K4)
+ BKP=B(K1)+B(K4)
+ BKM=B(K1)-B(K4)
+ AJP=A(K2)+A(K3)
+ AJM=A(K2)-A(K3)
+ BJP=B(K2)+B(K3)
+ BJM=B(K2)-B(K3)
+ AA=A(KK)
+ BB=B(KK)
+ A(KK)=AA+AKP+AJP
+ B(KK)=BB+BKP+BJP
+ AK=AKP*C72+AJP*C2+AA
+ BK=BKP*C72+BJP*C2+BB
+ AJ=AKM*S72+AJM*S2
+ BJ=BKM*S72+BJM*S2
+ A(K1)=AK-BJ
+ A(K4)=AK+BJ
+ B(K1)=BK+AJ
+ B(K4)=BK-AJ
+ AK=AKP*C2+AJP*C72+AA
+ BK=BKP*C2+BJP*C72+BB
+ AJ=AKM*S2-AJM*S72
+ BJ=BKM*S2-BJM*S72
+ A(K2)=AK-BJ
+ A(K3)=AK+BJ
+ B(K2)=BK+AJ
+ B(K3)=BK-AJ
+ KK=K4+KSPAN
+ IF(KK .LT. NN) GO TO 520
+ KK=KK-NN
+ IF(KK .LE. KSPAN) GO TO 520
+ GO TO 700
+C TRANSFORM FOR ODD FACTORS
+ 600 K=NFAC(I)
+ KSPNN=KSPAN
+ KSPAN=KSPAN/K
+ IF(K .EQ. 3) GO TO 320
+ IF(K .EQ. 5) GO TO 510
+ IF(K .EQ. JF) GO TO 640
+ JF=K
+ S1=RAD/(K)
+ C1=COS(S1)
+ S1=SIN(S1)
+ IF(JF .GT. MAXF) GO TO 998
+ CK(JF)=1.E0
+ SK(JF)=0.E0
+ J=1
+ 630 CK(J)=CK(K)*C1+SK(K)*S1
+ SK(J)=CK(K)*S1-SK(K)*C1
+ K=K-1
+ CK(K)=CK(J)
+ SK(K)=-SK(J)
+ J=J+1
+ IF(J .LT. K) GO TO 630
+ 640 K1=KK
+ K2=KK+KSPNN
+ AA=A(KK)
+ BB=B(KK)
+ AK=AA
+ BK=BB
+ J=1
+ K1=K1+KSPAN
+ 650 K2=K2-KSPAN
+ J=J+1
+ AT(J)=A(K1)+A(K2)
+ AK=AT(J)+AK
+ BT(J)=B(K1)+B(K2)
+ BK=BT(J)+BK
+ J=J+1
+ AT(J)=A(K1)-A(K2)
+ BT(J)=B(K1)-B(K2)
+ K1=K1+KSPAN
+ IF(K1 .LT. K2) GO TO 650
+ A(KK)=AK
+ B(KK)=BK
+ K1=KK
+ K2=KK+KSPNN
+ J=1
+ 660 K1=K1+KSPAN
+ K2=K2-KSPAN
+ JJ=J
+ AK=AA
+ BK=BB
+ AJ=0.E0
+ BJ=0.E0
+ K=1
+ 670 K=K+1
+ AK=AT(K)*CK(JJ)+AK
+ BK=BT(K)*CK(JJ)+BK
+ K=K+1
+ AJ=AT(K)*SK(JJ)+AJ
+ BJ=BT(K)*SK(JJ)+BJ
+ JJ=JJ+J
+ IF(JJ .GT. JF) JJ=JJ-JF
+ IF(K .LT. JF) GO TO 670
+ K=JF-J
+ A(K1)=AK-BJ
+ B(K1)=BK+AJ
+ A(K2)=AK+BJ
+ B(K2)=BK-AJ
+ J=J+1
+ IF(J .LT. K) GO TO 660
+ KK=KK+KSPNN
+ IF(KK .LE. NN) GO TO 640
+ KK=KK-NN
+ IF(KK .LE. KSPAN) GO TO 640
+C MULTIPLY BY ROTATION FACTOR (EXCEPT FOR FACTORS OF 2 AND 4)
+ 700 IF(I .EQ. M) GO TO 800
+ KK=JC+1
+ 710 C2=1.E0-CD
+ S1=SD
+ 720 C1=C2
+ S2=S1
+ KK=KK+KSPAN
+ 730 AK=A(KK)
+ A(KK)=C2*AK-S2*B(KK)
+ B(KK)=S2*AK+C2*B(KK)
+ KK=KK+KSPNN
+ IF(KK .LE. NT) GO TO 730
+ AK=S1*S2
+ S2=S1*C2+C1*S2
+ C2=C1*C2-AK
+ KK=KK-NT+KSPAN
+ IF(KK .LE. KSPNN) GO TO 730
+ C2=C1-(CD*C1+SD*S1)
+ S1=S1+(SD*C1-CD*S1)
+ KK=KK-KSPNN+JC
+ IF(KK .LE. KSPAN) GO TO 720
+ KK=KK-KSPAN+JC+INC
+ IF(KK .LE. JC+JC) GO TO 710
+ GO TO 100
+C PERMUTE THE RESULTS TO NORMAL ORDER---DONE IN TWO STAGES
+C PERMUTATION FOR SQUARE FACTORS OF N
+ 800 NP(1)=KS
+ IF(KT .EQ. 0) GO TO 890
+ K=KT+KT+1
+ IF(M .LT. K) K=K-1
+ J=1
+ NP(K+1)=JC
+ 810 NP(J+1)=NP(J)/NFAC(J)
+ NP(K)=NP(K+1)*NFAC(J)
+ J=J+1
+ K=K-1
+ IF(J .LT. K) GO TO 810
+ K3=NP(K+1)
+ KSPAN=NP(2)
+ KK=JC+1
+ K2=KSPAN+1
+ J=1
+ IF(N .NE. NTOT) GO TO 850
+C PERMUTATION FOR SINGLE-VARIATE TRANSFORM (OPTIONAL CODE)
+ 820 AK=A(KK)
+ A(KK)=A(K2)
+ A(K2)=AK
+ BK=B(KK)
+ B(KK)=B(K2)
+ B(K2)=BK
+ KK=KK+INC
+ K2=KSPAN+K2
+ IF(K2 .LT. KS) GO TO 820
+ 830 K2=K2-NP(J)
+ J=J+1
+ K2=NP(J+1)+K2
+ IF(K2 .GT. NP(J)) GO TO 830
+ J=1
+ 840 IF(KK .LT. K2) GO TO 820
+ KK=KK+INC
+ K2=KSPAN+K2
+ IF(K2 .LT. KS) GO TO 840
+ IF(KK .LT. KS) GO TO 830
+ JC=K3
+ GO TO 890
+C PERMUTATION FOR MULTIVARIATE TRANSFORM
+ 850 K=KK+JC
+ 860 AK=A(KK)
+ A(KK)=A(K2)
+ A(K2)=AK
+ BK=B(KK)
+ B(KK)=B(K2)
+ B(K2)=BK
+ KK=KK+INC
+ K2=K2+INC
+ IF(KK .LT. K) GO TO 860
+ KK=KK+KS-JC
+ K2=K2+KS-JC
+ IF(KK .LT. NT) GO TO 850
+ K2=K2-NT+KSPAN
+ KK=KK-NT+JC
+ IF(K2 .LT. KS) GO TO 850
+ 870 K2=K2-NP(J)
+ J=J+1
+ K2=NP(J+1)+K2
+ IF(K2 .GT. NP(J)) GO TO 870
+ J=1
+ 880 IF(KK .LT. K2) GO TO 850
+ KK=KK+JC
+ K2=KSPAN+K2
+ IF(K2 .LT. KS) GO TO 880
+ IF(KK .LT. KS) GO TO 870
+ JC=K3
+ 890 IF(2*KT+1 .GE. M) RETURN
+ KSPNN=NP(KT+1)
+C PERMUTATION FOR SQUARE-FREE FACTORS OF N
+ J=M-KT
+ NFAC(J+1)=1
+ 900 NFAC(J)=NFAC(J)*NFAC(J+1)
+ J=J-1
+ IF(J .NE. KT) GO TO 900
+ KT=KT+1
+ NN=NFAC(KT)-1
+ IF(NN .GT. MAXP) GO TO 998
+ JJ=0
+ J=0
+ GO TO 906
+ 902 JJ=JJ-K2
+ K2=KK
+ K=K+1
+ KK=NFAC(K)
+ 904 JJ=KK+JJ
+ IF(JJ .GE. K2) GO TO 902
+ NP(J)=JJ
+ 906 K2=NFAC(KT)
+ K=KT+1
+ KK=NFAC(K)
+ J=J+1
+ IF(J .LE. NN) GO TO 904
+C DETERMINE THE PERMUTATION CYCLES OF LENGTH GREATER THAN 1
+ J=0
+ GO TO 914
+ 910 K=KK
+ KK=NP(K)
+ NP(K)=-KK
+ IF(KK .NE. J) GO TO 910
+ K3=KK
+ 914 J=J+1
+ KK=NP(J)
+ IF(KK .LT. 0) GO TO 914
+ IF(KK .NE. J) GO TO 910
+ NP(J)=-J
+ IF(J .NE. NN) GO TO 914
+ MAXF=INC*MAXF
+C REORDER A AND B, FOLLOWING THE PERMUTATION CYCLES
+ GO TO 950
+ 924 J=J-1
+ IF(NP(J) .LT. 0) GO TO 924
+ JJ=JC
+ 926 KSPAN=JJ
+ IF(JJ .GT. MAXF) KSPAN=MAXF
+ JJ=JJ-KSPAN
+ K=NP(J)
+ KK=JC*K+II+JJ
+ K1=KK+KSPAN
+ K2=0
+ 928 K2=K2+1
+ AT(K2)=A(K1)
+ BT(K2)=B(K1)
+ K1=K1-INC
+ IF(K1 .NE. KK) GO TO 928
+ 932 K1=KK+KSPAN
+ K2=K1-JC*(K+NP(K))
+ K=-NP(K)
+ 936 A(K1)=A(K2)
+ B(K1)=B(K2)
+ K1=K1-INC
+ K2=K2-INC
+ IF(K1 .NE. KK) GO TO 936
+ KK=K2
+ IF(K .NE. J) GO TO 932
+ K1=KK+KSPAN
+ K2=0
+ 940 K2=K2+1
+ A(K1)=AT(K2)
+ B(K1)=BT(K2)
+ K1=K1-INC
+ IF(K1 .NE. KK) GO TO 940
+ IF(JJ .NE. 0) GO TO 926
+ IF(J .NE. 1) GO TO 924
+ 950 J=K3+1
+ NT=NT-KSPNN
+ II=NT-INC+1
+ IF(NT .GE. 0) GO TO 924
+ RETURN
+C ERROR FINISH, INSUFFICIENT ARRAY STORAGE
+ 998 IERR=1
+ RETURN
+ END
+ SUBROUTINE FFTLDP(X,N,NDIR,DT,IERR)
+C FFTLDP IS A DOUBLE PRECISION VERSION OF FFTL. SEE NOTES TO FFTL FOR USAGE
+C THE SCALING IN FFTLDP DIFFERS FROM THAT IN FFTL :
+C THE FORWARD TRANSFORM IS SPECTRUM(W)=DT*SUM(0,N-1)X(J)*EXP(-I*W*J*DT)
+C THE SPACING IN W IS 2*PI/(N*DT)
+ IMPLICIT REAL*8(A-H,O-Z)
+ DIMENSION X(*)
+ N2=N/2
+ IDIR=NDIR
+ IF(IDIR.LT.0) IDIR=-IDIR
+ GO TO (1,2),IDIR
+C DO FORWARD TRANSFORM (IE. TIME TO FREQUENCY).
+ 1 CALL FFTDP(X,X(2),N2,N2,N2,2,IERR)
+ CALL RLTRDP(X,X(2),N2,2)
+ N1=2*N2+2
+ SCALE=0.5D0*DT
+ IF(NDIR.GT.0) GO TO 3
+ DO 5 I=4,N,2
+ 5 X(I)=-X(I)
+ GO TO 3
+C DO BACKWARD TRANSFORM (IE. FREQUENCY TO TIME).
+ 2 IF(NDIR.GT.0) GO TO 6
+ DO 7 I=4,N,2
+ 7 X(I)=-X(I)
+ 6 X(2)=0.D0
+ N1=2*N2
+ X(N1+2)=0.D0
+ CALL RLTRDP(X,X(2),N2,-2)
+ CALL FFTDP(X,X(2),N2,N2,N2,-2,IERR)
+ SCALE=1.D0/(N*DT)
+C NORMALIZE THE TRANSFORM.
+ 3 DO 4 I=1,N1
+ 4 X(I)=SCALE*X(I)
+ RETURN
+ END
+ SUBROUTINE RLTRDP(A,B,N,ISN)
+C RLTRDP IS A DOUBLE PRECISION VERSION OF REALTR. SEE REALTR FOR NOTES
+C ON USAGE.
+ IMPLICIT REAL*8(A-H,O-Z)
+ DIMENSION A(*),B(*)
+ IF(N .LE. 1) RETURN
+ INC=ISN
+ IF(INC.LT.0) INC=-INC
+ NK = N * INC + 2
+ NH = NK / 2
+ SD = 3.1415926535898D0/(2.D0*N)
+ CD = 2.D0 * DSIN(SD)**2
+ SD = DSIN(SD+SD)
+ SN = 0.D0
+ IF(ISN .GT. 0) GO TO 10
+ CN = -1.D0
+ SD = -SD
+ GO TO 20
+ 10 CN = 1.D0
+ A(NK-1) = A(1)
+ B(NK-1) = B(1)
+ 20 DO 30 J=1,NH,INC
+ K = NK - J
+ AA = A(J) + A(K)
+ AB = A(J) - A(K)
+ BA = B(J) + B(K)
+ BB = B(J) - B(K)
+ XX = CN * BA + SN * AB
+ YY = SN * BA - CN * AB
+ B(K) = YY - BB
+ B(J) = YY + BB
+ A(K) = AA - XX
+ A(J) = AA + XX
+ AA = CN - (CD * CN + SD * SN)
+ SN = (SD * CN - CD * SN) + SN
+ 30 CN = AA
+ RETURN
+ END
+ SUBROUTINE FFTDP(A,B,NTOT,N,NSPAN,ISN,IERR)
+C FFTDP IS A DOUBLE PRECISION VERSION OF SINGLETON'S FFT PROCEDURE. SEE NOTES
+C TO FFT FOR USAGE
+ IMPLICIT REAL*8(A-H),REAL*8(O-Z)
+ DIMENSION A(*),B(*)
+ DIMENSION NFAC(11),NP(209)
+C ARRAY STORAGE FOR MAXIMUM PRIME FACTOR OF 23
+ DIMENSION AT(23),CK(23),BT(23),SK(23)
+ EQUIVALENCE (I,II)
+C THE FOLLOWING TWO CONSTANTS SHOULD AGREE WITH THE ARRAY DIMENSIONS.
+ MAXF=23
+ MAXP=209
+ IERR=0
+ IF(N .LT. 2) RETURN
+ INC=ISN
+ C72=0.30901699437494D0
+ S72=0.95105651629515D0
+ S120=0.86602540378443D0
+ RAD=6.2831853071796D0
+ IF(ISN .GE. 0) GO TO 10
+ S72=-S72
+ S120=-S120
+ RAD=-RAD
+ INC=-INC
+ 10 NT=INC*NTOT
+ KS=INC*NSPAN
+ KSPAN=KS
+ NN=NT-INC
+ JC=KS/N
+ RADF=RAD*(JC*0.5D0)
+ I=0
+ JF=0
+C DETERMINE THE FACTORS OF N
+ M=0
+ K=N
+ GO TO 20
+ 15 M=M+1
+ NFAC(M)=4
+ K=K/16
+ 20 IF(K-(K/16)*16 .EQ. 0) GO TO 15
+ J=3
+ JJ=9
+ GO TO 30
+ 25 M=M+1
+ NFAC(M)=J
+ K=K/JJ
+ 30 IF(MOD(K,JJ) .EQ. 0) GO TO 25
+ J=J+2
+ JJ=J**2
+ IF(JJ .LE. K) GO TO 30
+ IF(K .GT. 4) GO TO 40
+ KT=M
+ NFAC(M+1)=K
+ IF(K .NE. 1) M=M+1
+ GO TO 80
+ 40 IF(K-(K/4)*4 .NE. 0) GO TO 50
+ M=M+1
+ NFAC(M)=2
+ K=K/4
+ 50 KT=M
+ J=2
+ 60 IF(MOD(K,J) .NE. 0) GO TO 70
+ M=M+1
+ NFAC(M)=J
+ K=K/J
+ 70 J=((J+1)/2)*2+1
+ IF(J .LE. K) GO TO 60
+ 80 IF(KT .EQ. 0) GO TO 100
+ J=KT
+ 90 M=M+1
+ NFAC(M)=NFAC(J)
+ J=J-1
+ IF(J .NE. 0) GO TO 90
+C COMPUTE FOURIER TRANSFORM
+ 100 SD=RADF/(KSPAN)
+ CD=2.D0*DSIN(SD)**2
+ SD=DSIN(SD+SD)
+ KK=1
+ I=I+1
+ IF(NFAC(I) .NE. 2) GO TO 400
+C TRANSFORM FOR FACTOR OF 2 (INCLUDING ROTATION FACTOR)
+ KSPAN=KSPAN/2
+ K1=KSPAN+2
+ 210 K2=KK+KSPAN
+ AK=A(K2)
+ BK=B(K2)
+ A(K2)=A(KK)-AK
+ B(K2)=B(KK)-BK
+ A(KK)=A(KK)+AK
+ B(KK)=B(KK)+BK
+ KK=K2+KSPAN
+ IF(KK .LE. NN) GO TO 210
+ KK=KK-NN
+ IF(KK .LE. JC) GO TO 210
+ IF(KK .GT. KSPAN) GO TO 800
+ 220 C1=1.D0-CD
+ S1=SD
+ 230 K2=KK+KSPAN
+ AK=A(KK)-A(K2)
+ BK=B(KK)-B(K2)
+ A(KK)=A(KK)+A(K2)
+ B(KK)=B(KK)+B(K2)
+ A(K2)=C1*AK-S1*BK
+ B(K2)=S1*AK+C1*BK
+ KK=K2+KSPAN
+ IF(KK .LT. NT) GO TO 230
+ K2=KK-NT
+ C1=-C1
+ KK=K1-K2
+ IF(KK .GT. K2) GO TO 230
+ AK=CD*C1+SD*S1
+ S1=(SD*C1-CD*S1)+S1
+ C1=C1-AK
+ KK=KK+JC
+ IF(KK .LT. K2) GO TO 230
+ K1=K1+INC+INC
+ KK=(K1-KSPAN)/2+JC
+ IF(KK .LE. JC+JC) GO TO 220
+ GO TO 100
+C TRANSFORM FOR FACTOR OF 3 (OPTIONAL CODE)
+ 320 K1=KK+KSPAN
+ K2=K1+KSPAN
+ AK=A(KK)
+ BK=B(KK)
+ AJ=A(K1)+A(K2)
+ BJ=B(K1)+B(K2)
+ A(KK)=AK+AJ
+ B(KK)=BK+BJ
+ AK=-0.5*AJ+AK
+ BK=-0.5*BJ+BK
+ AJ=(A(K1)-A(K2))*S120
+ BJ=(B(K1)-B(K2))*S120
+ A(K1)=AK-BJ
+ B(K1)=BK+AJ
+ A(K2)=AK+BJ
+ B(K2)=BK-AJ
+ KK=K2+KSPAN
+ IF(KK .LT. NN) GO TO 320
+ KK=KK-NN
+ IF(KK .LE. KSPAN) GO TO 320
+ GO TO 700
+C TRANSFORM FOR FACTOR OF 4
+ 400 IF(NFAC(I) .NE. 4) GO TO 600
+ KSPNN=KSPAN
+ KSPAN=KSPAN/4
+ 410 C1=1.D0
+ S1=0
+ 420 K1=KK+KSPAN
+ K2=K1+KSPAN
+ K3=K2+KSPAN
+ AKP=A(KK)+A(K2)
+ AKM=A(KK)-A(K2)
+ AJP=A(K1)+A(K3)
+ AJM=A(K1)-A(K3)
+ A(KK)=AKP+AJP
+ AJP=AKP-AJP
+ BKP=B(KK)+B(K2)
+ BKM=B(KK)-B(K2)
+ BJP=B(K1)+B(K3)
+ BJM=B(K1)-B(K3)
+ B(KK)=BKP+BJP
+ BJP=BKP-BJP
+ IF(ISN .LT. 0) GO TO 450
+ AKP=AKM-BJM
+ AKM=AKM+BJM
+ BKP=BKM+AJM
+ BKM=BKM-AJM
+ IF(S1 .EQ. 0) GO TO 460
+ 430 A(K1)=AKP*C1-BKP*S1
+ B(K1)=AKP*S1+BKP*C1
+ A(K2)=AJP*C2-BJP*S2
+ B(K2)=AJP*S2+BJP*C2
+ A(K3)=AKM*C3-BKM*S3
+ B(K3)=AKM*S3+BKM*C3
+ KK=K3+KSPAN
+ IF(KK .LE. NT) GO TO 420
+ 440 C2=CD*C1+SD*S1
+ S1=(SD*C1-CD*S1)+S1
+ C1=C1-C2
+ C2=C1**2-S1**2
+ S2=2.D0*C1*S1
+ C3=C2*C1-S2*S1
+ S3=C2*S1+S2*C1
+ KK=KK-NT+JC
+ IF(KK .LE. KSPAN) GO TO 420
+ KK=KK-KSPAN+INC
+ IF(KK .LE. JC) GO TO 410
+ IF(KSPAN .EQ. JC) GO TO 800
+ GO TO 100
+ 450 AKP=AKM+BJM
+ AKM=AKM-BJM
+ BKP=BKM-AJM
+ BKM=BKM+AJM
+ IF(S1 .NE. 0) GO TO 430
+ 460 A(K1)=AKP
+ B(K1)=BKP
+ A(K2)=AJP
+ B(K2)=BJP
+ A(K3)=AKM
+ B(K3)=BKM
+ KK=K3+KSPAN
+ IF(KK .LE. NT) GO TO 420
+ GO TO 440
+C TRANSFORM FOR FACTOR OF 5 (OPTIONAL CODE)
+ 510 C2=C72**2-S72**2
+ S2=2.D0*C72*S72
+ 520 K1=KK+KSPAN
+ K2=K1+KSPAN
+ K3=K2+KSPAN
+ K4=K3+KSPAN
+ AKP=A(K1)+A(K4)
+ AKM=A(K1)-A(K4)
+ BKP=B(K1)+B(K4)
+ BKM=B(K1)-B(K4)
+ AJP=A(K2)+A(K3)
+ AJM=A(K2)-A(K3)
+ BJP=B(K2)+B(K3)
+ BJM=B(K2)-B(K3)
+ AA=A(KK)
+ BB=B(KK)
+ A(KK)=AA+AKP+AJP
+ B(KK)=BB+BKP+BJP
+ AK=AKP*C72+AJP*C2+AA
+ BK=BKP*C72+BJP*C2+BB
+ AJ=AKM*S72+AJM*S2
+ BJ=BKM*S72+BJM*S2
+ A(K1)=AK-BJ
+ A(K4)=AK+BJ
+ B(K1)=BK+AJ
+ B(K4)=BK-AJ
+ AK=AKP*C2+AJP*C72+AA
+ BK=BKP*C2+BJP*C72+BB
+ AJ=AKM*S2-AJM*S72
+ BJ=BKM*S2-BJM*S72
+ A(K2)=AK-BJ
+ A(K3)=AK+BJ
+ B(K2)=BK+AJ
+ B(K3)=BK-AJ
+ KK=K4+KSPAN
+ IF(KK .LT. NN) GO TO 520
+ KK=KK-NN
+ IF(KK .LE. KSPAN) GO TO 520
+ GO TO 700
+C TRANSFORM FOR ODD FACTORS
+ 600 K=NFAC(I)
+ KSPNN=KSPAN
+ KSPAN=KSPAN/K
+ IF(K .EQ. 3) GO TO 320
+ IF(K .EQ. 5) GO TO 510
+ IF(K .EQ. JF) GO TO 640
+ JF=K
+ S1=RAD/(K)
+ C1=DCOS(S1)
+ S1=DSIN(S1)
+ IF(JF .GT. MAXF) GO TO 998
+ CK(JF)=1.D0
+ SK(JF)=0.D0
+ J=1
+ 630 CK(J)=CK(K)*C1+SK(K)*S1
+ SK(J)=CK(K)*S1-SK(K)*C1
+ K=K-1
+ CK(K)=CK(J)
+ SK(K)=-SK(J)
+ J=J+1
+ IF(J .LT. K) GO TO 630
+ 640 K1=KK
+ K2=KK+KSPNN
+ AA=A(KK)
+ BB=B(KK)
+ AK=AA
+ BK=BB
+ J=1
+ K1=K1+KSPAN
+ 650 K2=K2-KSPAN
+ J=J+1
+ AT(J)=A(K1)+A(K2)
+ AK=AT(J)+AK
+ BT(J)=B(K1)+B(K2)
+ BK=BT(J)+BK
+ J=J+1
+ AT(J)=A(K1)-A(K2)
+ BT(J)=B(K1)-B(K2)
+ K1=K1+KSPAN
+ IF(K1 .LT. K2) GO TO 650
+ A(KK)=AK
+ B(KK)=BK
+ K1=KK
+ K2=KK+KSPNN
+ J=1
+ 660 K1=K1+KSPAN
+ K2=K2-KSPAN
+ JJ=J
+ AK=AA
+ BK=BB
+ AJ=0.D0
+ BJ=0.D0
+ K=1
+ 670 K=K+1
+ AK=AT(K)*CK(JJ)+AK
+ BK=BT(K)*CK(JJ)+BK
+ K=K+1
+ AJ=AT(K)*SK(JJ)+AJ
+ BJ=BT(K)*SK(JJ)+BJ
+ JJ=JJ+J
+ IF(JJ .GT. JF) JJ=JJ-JF
+ IF(K .LT. JF) GO TO 670
+ K=JF-J
+ A(K1)=AK-BJ
+ B(K1)=BK+AJ
+ A(K2)=AK+BJ
+ B(K2)=BK-AJ
+ J=J+1
+ IF(J .LT. K) GO TO 660
+ KK=KK+KSPNN
+ IF(KK .LE. NN) GO TO 640
+ KK=KK-NN
+ IF(KK .LE. KSPAN) GO TO 640
+C MULTIPLY BY ROTATION FACTOR (EXCEPT FOR FACTORS OF 2 AND 4)
+ 700 IF(I .EQ. M) GO TO 800
+ KK=JC+1
+ 710 C2=1.D0-CD
+ S1=SD
+ 720 C1=C2
+ S2=S1
+ KK=KK+KSPAN
+ 730 AK=A(KK)
+ A(KK)=C2*AK-S2*B(KK)
+ B(KK)=S2*AK+C2*B(KK)
+ KK=KK+KSPNN
+ IF(KK .LE. NT) GO TO 730
+ AK=S1*S2
+ S2=S1*C2+C1*S2
+ C2=C1*C2-AK
+ KK=KK-NT+KSPAN
+ IF(KK .LE. KSPNN) GO TO 730
+ C2=C1-(CD*C1+SD*S1)
+ S1=S1+(SD*C1-CD*S1)
+ KK=KK-KSPNN+JC
+ IF(KK .LE. KSPAN) GO TO 720
+ KK=KK-KSPAN+JC+INC
+ IF(KK .LE. JC+JC) GO TO 710
+ GO TO 100
+C PERMUTE THE RESULTS TO NORMAL ORDER---DONE IN TWO STAGES
+C PERMUTATION FOR SQUARE FACTORS OF N
+ 800 NP(1)=KS
+ IF(KT .EQ. 0) GO TO 890
+ K=KT+KT+1
+ IF(M .LT. K) K=K-1
+ J=1
+ NP(K+1)=JC
+ 810 NP(J+1)=NP(J)/NFAC(J)
+ NP(K)=NP(K+1)*NFAC(J)
+ J=J+1
+ K=K-1
+ IF(J .LT. K) GO TO 810
+ K3=NP(K+1)
+ KSPAN=NP(2)
+ KK=JC+1
+ K2=KSPAN+1
+ J=1
+ IF(N .NE. NTOT) GO TO 850
+C PERMUTATION FOR SINGLE-VARIATE TRANSFORM (OPTIONAL CODE)
+ 820 AK=A(KK)
+ A(KK)=A(K2)
+ A(K2)=AK
+ BK=B(KK)
+ B(KK)=B(K2)
+ B(K2)=BK
+ KK=KK+INC
+ K2=KSPAN+K2
+ IF(K2 .LT. KS) GO TO 820
+ 830 K2=K2-NP(J)
+ J=J+1
+ K2=NP(J+1)+K2
+ IF(K2 .GT. NP(J)) GO TO 830
+ J=1
+ 840 IF(KK .LT. K2) GO TO 820
+ KK=KK+INC
+ K2=KSPAN+K2
+ IF(K2 .LT. KS) GO TO 840
+ IF(KK .LT. KS) GO TO 830
+ JC=K3
+ GO TO 890
+C PERMUTATION FOR MULTIVARIATE TRANSFORM
+ 850 K=KK+JC
+ 860 AK=A(KK)
+ A(KK)=A(K2)
+ A(K2)=AK
+ BK=B(KK)
+ B(KK)=B(K2)
+ B(K2)=BK
+ KK=KK+INC
+ K2=K2+INC
+ IF(KK .LT. K) GO TO 860
+ KK=KK+KS-JC
+ K2=K2+KS-JC
+ IF(KK .LT. NT) GO TO 850
+ K2=K2-NT+KSPAN
+ KK=KK-NT+JC
+ IF(K2 .LT. KS) GO TO 850
+ 870 K2=K2-NP(J)
+ J=J+1
+ K2=NP(J+1)+K2
+ IF(K2 .GT. NP(J)) GO TO 870
+ J=1
+ 880 IF(KK .LT. K2) GO TO 850
+ KK=KK+JC
+ K2=KSPAN+K2
+ IF(K2 .LT. KS) GO TO 880
+ IF(KK .LT. KS) GO TO 870
+ JC=K3
+ 890 IF(2*KT+1 .GE. M) RETURN
+ KSPNN=NP(KT+1)
+C PERMUTATION FOR SQUARE-FREE FACTORS OF N
+ J=M-KT
+ NFAC(J+1)=1
+ 900 NFAC(J)=NFAC(J)*NFAC(J+1)
+ J=J-1
+ IF(J .NE. KT) GO TO 900
+ KT=KT+1
+ NN=NFAC(KT)-1
+ IF(NN .GT. MAXP) GO TO 998
+ JJ=0
+ J=0
+ GO TO 906
+ 902 JJ=JJ-K2
+ K2=KK
+ K=K+1
+ KK=NFAC(K)
+ 904 JJ=KK+JJ
+ IF(JJ .GE. K2) GO TO 902
+ NP(J)=JJ
+ 906 K2=NFAC(KT)
+ K=KT+1
+ KK=NFAC(K)
+ J=J+1
+ IF(J .LE. NN) GO TO 904
+C DETERMINE THE PERMUTATION CYCLES OF LENGTH GREATER THAN 1
+ J=0
+ GO TO 914
+ 910 K=KK
+ KK=NP(K)
+ NP(K)=-KK
+ IF(KK .NE. J) GO TO 910
+ K3=KK
+ 914 J=J+1
+ KK=NP(J)
+ IF(KK .LT. 0) GO TO 914
+ IF(KK .NE. J) GO TO 910
+ NP(J)=-J
+ IF(J .NE. NN) GO TO 914
+ MAXF=INC*MAXF
+C REORDER A AND B, FOLLOWING THE PERMUTATION CYCLES
+ GO TO 950
+ 924 J=J-1
+ IF(NP(J) .LT. 0) GO TO 924
+ JJ=JC
+ 926 KSPAN=JJ
+ IF(JJ .GT. MAXF) KSPAN=MAXF
+ JJ=JJ-KSPAN
+ K=NP(J)
+ KK=JC*K+II+JJ
+ K1=KK+KSPAN
+ K2=0
+ 928 K2=K2+1
+ AT(K2)=A(K1)
+ BT(K2)=B(K1)
+ K1=K1-INC
+ IF(K1 .NE. KK) GO TO 928
+ 932 K1=KK+KSPAN
+ K2=K1-JC*(K+NP(K))
+ K=-NP(K)
+ 936 A(K1)=A(K2)
+ B(K1)=B(K2)
+ K1=K1-INC
+ K2=K2-INC
+ IF(K1 .NE. KK) GO TO 936
+ KK=K2
+ IF(K .NE. J) GO TO 932
+ K1=KK+KSPAN
+ K2=0
+ 940 K2=K2+1
+ A(K1)=AT(K2)
+ B(K1)=BT(K2)
+ K1=K1-INC
+ IF(K1 .NE. KK) GO TO 940
+ IF(JJ .NE. 0) GO TO 926
+ IF(J .NE. 1) GO TO 924
+ 950 J=K3+1
+ NT=NT-KSPNN
+ II=NT-INC+1
+ IF(NT .GE. 0) GO TO 924
+ RETURN
+C ERROR FINISH, INSUFFICIENT ARRAY STORAGE
+ 998 IERR=1
+ RETURN
+ END
\ No newline at end of file
Added: mc/2D/ConMan/trunk/src/geoid.F
===================================================================
--- mc/2D/ConMan/trunk/src/geoid.F (rev 0)
+++ mc/2D/ConMan/trunk/src/geoid.F 2008-07-15 15:54:31 UTC (rev 12416)
@@ -0,0 +1,200 @@
+ subroutine geoid (x, t, stress)
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c c
+c subroutine to calculate GEOID from buoyancy and stress field c
+c from ConMan in the 2-D Cartesian case. GEOID uses the c
+c runfile to get the geom file (for coordinates) the c
+c temperature & velocity file for the buoyancies and the c
+c stress file to get the topography at the surface and bottom. c
+c Written by Scott King Feb 9, 1990 at Caltech. c
+c c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+c
+c parameter for array dimensions, greater than the number of nodes
+c
+c
+ implicit double precision (a-h,o-z)
+ PARAMETER (NWAVE=2000)
+c
+#include<common.h>
+#include<materials.h>
+c
+ double precision x(2,*), t(*), stress(6,*)
+ real z1, z2, z3, aratio, pi
+ real tlay(NWAVE), tintz(NWAVE), wave(NWAVE,NWAVE)
+ real geoid1(NWAVE),rtzzt(NWAVE),rtzzb(NWAVE), geoid2(NWAVE)
+c
+c constants for dimensionalization (geoid and topography in meters)
+c the values used for the Blankenback et al. Benchmark 1a
+c (Ra=1.0e4, constant viscosity)
+c
+c d = 1.0e6
+c g = 1.0e1
+c rho = 4.0e3
+c evisc = 2.5e19 * rho * 1.0e4 / ra(1)
+c dif = 1.0e-6
+c bigG = 6.673e-11
+c alfa = 2.5e-5
+c deltaT = 1.0e3
+c
+ d = 3.000e6
+ g = 1.0e1
+ rho = 4.0e3
+ evisc = 4.32e21
+ dif = 1.0e-6
+ bigG = 6.673e-11
+ alfa = 2.0e-5
+ deltaT = 2.000e3
+ pi = 3.14159272e0
+ aratio = x(1,numnp)
+c
+c if the problem uses reflecting boundary conditions, this is correct
+c we also need to reflect the solution about the right boundary
+c to get the correct symmetry.
+c
+ npoints = nelx + 1
+c
+c otherwise discard the last point - it is really the first (not tested)
+c nb: if npoints is odd, the tranform will ignore the last point
+c
+ if (nwrap .gt. 0) npoints = nelx
+ if (nwrap .gt. 0) ntrans = nelx
+ if (nwrap .eq. 0) then
+ do 101 i=1,npoints
+ rtzzt(i) = -stress(5,(nelz+1)*i)
+ rtzzb(i) = -stress(5,1+(nelz+1)*(i-1))
+ rtzzt(2*npoints-i) = rtzzt(i)
+ rtzzb(2*npoints-i) = rtzzb(i)
+101 continue
+ ntrans = 2*npoints - 2
+ endif
+c
+c calculate fourier decomposition of the top and bottom
+c nb: if npoints is odd, the tranform will ignore the last point
+c
+ call fftl(rtzzt, ntrans, 1, ierror)
+c
+c bottom: (assumes iz is the fastest varying index)
+c
+ call fftl(rtzzb, ntrans, 1, ierror)
+c
+c calculate fourier decomposition at each layer
+c
+ do 600 iz = 1, nelz+1
+ do 500 ix = 1, npoints
+ tlay(ix) = t(iz + (nelz+1) * (ix-1))
+500 continue
+ if (nwrap .eq. 0) then
+ do 501 i=1,npoints
+ tlay(2*npoints-i) = tlay(i)
+501 continue
+ endif
+ call fftl(tlay, ntrans, 1, ierror)
+c
+c move twave
+c
+ do 550 ik=1,ntrans
+ wave(iz,ik) = tlay(ik)
+550 continue
+600 continue
+c -kz
+c integrate with depth e
+c this assumes a z fastest node ordering!
+c
+ numwav = ntrans/2 + 1
+ do 900 ik = 1, numwav
+c
+c integral over shape by Simpson's Rule
+c
+c note that "jk" is the array index while "ik" is the wavenumber
+c (there are two coeficients for each wavenumber
+c
+ jk = 2*ik
+ tintz(jk-1) = 0.0e0
+ tintz(jk) = 0.0e0
+c
+c use a 3 point 3rd order rule for the integral
+c
+c do 800 iz = 2, nelz, 2
+ do 800 iz = 2, nelz+1
+ z3 = x(2,iz+1)
+ z2 = x(2,iz)
+ z1 = x(2,iz-1)
+c tintz(jk-1) = tintz(jk-1) +
+c & ( exp(-1.0*pi*float(ik-1)*(1.0-z3)/aratio)*wave(iz+1,jk-1)+
+c & 4.0*exp(-1.0*pi*float(ik-1)*(1.0-z2)/aratio)*wave(iz,jk-1)+
+c & exp(-1.0*pi*float(ik-1)*(1.0-z1)/aratio)*wave(iz-1,jk-1))*
+c & (z2-z1)/3.0
+c tintz(jk) = tintz(jk) +
+c & ( exp(-1.0*pi*float(ik-1)*(1.0-z3)/aratio)*wave(iz+1,jk)+
+c & 4.0*exp(-1.0*pi*float(ik-1)*(1.0-z2)/aratio)*wave(iz,jk)+
+c & exp(-1.0*pi*float(ik-1)*(1.0-z1)/aratio)*wave(iz-1,jk))*
+c & (z2-z1)/3.0
+c two point rule
+ tintz(jk-1) = tintz(jk-1) +
+ & ( exp(-1.0*pi*float(ik-1)*(one-z2)/aratio)*wave(iz,jk-1)+
+ & exp(-1.0*pi*float(ik-1)*(one-z1)/aratio)*wave(iz-1,jk-1))*
+ & (z2-z1)/2.0
+ tintz(jk) = tintz(jk) +
+ & ( exp(-1.0*pi*float(ik-1)*(one-z2)/aratio)*wave(iz,jk)+
+ & exp(-1.0*pi*float(ik-1)*(one-z1)/aratio)*wave(iz-1,jk))*
+ & (z2-z1)/2.0
+800 continue
+c
+900 continue
+c
+c zero out the k=0 harminics to remove the mean
+c
+ rtzzt(1) = 0.0
+ rtzzt(2) = 0.0
+ rtzzb(1) = 0.0
+ rtzzb(2) = 0.0
+ tintz(1) = 0.0
+ tintz(2) = 0.0
+ geoid1(1) = 0.0
+ geoid1(2) = 0.0
+c
+c spit out the result
+c
+ do 1100 ik = 2, numwav
+ jk = 2*ik
+ geoid1(jk-1) = (1.0/float(ik-1))*(rtzzt(jk-1)*evisc*dif/(d*d*g)
+ & -rtzzb(jk-1)*exp(-1.0*pi*float(ik-1)/aratio)*evisc*dif/(d*d*g)
+ & -tintz(jk-1)*rho*alfa*deltaT*d )
+ geoid1(jk) = (1.0/float(ik-1))*(rtzzt(jk)*evisc*dif/(d*d*g)
+ & -rtzzb(jk)*exp(-1.0*pi*float(ik-1)/aratio)*evisc*dif/(d*d*g)
+ & -tintz(jk)*rho*alfa*deltaT*d )
+ geoid2(jk-1)=-(1.0/float(ik-1))*
+ & tintz(jk-1)*rho*alfa*deltaT*d
+ geoid2(jk) =-(1.0/float(ik-1))*
+ & tintz(jk)*rho*afha*deltaT*d
+1100 continue
+c
+1006 format(1x,i5,1x,2(1pe9.3,1x),5(1pe12.5,1x))
+1007 format(1x,i5,1x,5(1pe12.5,1x))
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+c
+c Transforms Wave Number Data Back to real world for use with plot
+c
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+ call fftl(rtzzt, ntrans, 2, ierror)
+ call fftl(rtzzb, ntrans, 2, ierror)
+ call fftl(geoid1, ntrans, 2, ierror)
+ call fftl(geoid2, ntrans, 2, ierror)
+c
+c warning if npoints is odd, the last point is garbage
+c
+ pk = 2 * bigG*d/g
+ do 200 i=1,npoints
+ xm = aratio*float(i-1)/(npoints-1)
+ write(igeoid,1001) xm, rtzzt(i)*evisc*dif / (rho*d*d*g) ,
+ & pk*geoid2(i), pk * geoid1(i)
+200 continue
+c call myflush(igeoid)
+c
+1001 format(4e15.7)
+ return
+ end
\ No newline at end of file
Modified: mc/2D/ConMan/trunk/src/stress.F
===================================================================
--- mc/2D/ConMan/trunk/src/stress.F 2008-07-14 19:28:54 UTC (rev 12415)
+++ mc/2D/ConMan/trunk/src/stress.F 2008-07-15 15:54:31 UTC (rev 12416)
@@ -11,10 +11,10 @@
#include<materials.h>
c
dimension x(nsd,*) , v(ndof,*) , t(*) ,
- & shdx(4,5) , shdy(4,5) , det(5), tl(4)
+ & shdx(4,5) , shdy(4,5) , det(5)
c
c
- common /temp1 / xl(2,4), vl(8) , strtmp(6) , rhsl(8) ,
+ common /temp1 / xl(2,4), vl(8) , strtmp(6) ,
& volume , smass , tmass(4) , evisc(5)
c
real*8 maxsecinv, secinv
@@ -23,6 +23,7 @@
do 50 n = 1, 6
do 45 i = 1, numnp
stress(n,i) = zero
+ pmass(i) = zero
45 continue
50 continue
@@ -35,12 +36,6 @@
xl(2, n) = x(2,ien(ivel,n) )
vl(2*n-1) = v(1,ien(ivel,n) )
vl(2*n ) = v(2,ien(ivel,n) )
-c xl(1, n) = x(1,ien(ivel,n) ) * 3.0e6
-c xl(2, n) = x(2,ien(ivel,n) ) * 3.0e6
-c vl(2*n-1) = v(1,ien(ivel,n) ) * 1.0e-6 / 3.0e6
-c vl(2*n ) = v(2,ien(ivel,n) ) * 1.0e-6 / 3.0e6
- tl(n ) = t(ien(ivel,n) )
-c*3.16e7/6e5
enddo
c
c.... call the global shape function
@@ -50,110 +45,77 @@
c
c find velocity derivative at element centers
c
-c normal stress x
- strtmp(1) = 1.0*(shdx(1,5)*vl(1) + shdx(2,5)*vl(3)
- & + shdx(3,5)*vl(5) + shdx(4,5)*vl(7))
-c normal stress y
- strtmp(2) = 1.0*(shdy(1,5)*vl(2) + shdy(2,5)*vl(4)
- & + shdy(3,5)*vl(6) + shdy(4,5)*vl(8))
-c shear stress
- strtmp(3) = 0.5*(shdx(1,5)*vl(2) + shdx(2,5)*vl(4)
- & + shdx(3,5)*vl(6) + shdx(4,5)*vl(8)
- & + shdy(1,5)*vl(1) + shdy(2,5)*vl(3)
- & + shdy(3,5)*vl(5) + shdy(4,5)*vl(7))
-c pressure
- strtmp(4) = ( strtmp(1) + strtmp(2) )
-c second invariant
- strtmp(5) = abs( strtmp(2)*strtmp(2)
- & + strtmp(1)*strtmp(1)
- & + strtmp(3)*strtmp(3) )
-
-c
-c..... scale by material parameters use last time steps viscosity!
-c
- if (strtmp(5) .gt. zero ) then
- strtmp(5) = sqrt(0.5*strtmp(5))
- endif
-
- if (strtmp(4) .gt. maxsecinv) then
- maxsecinv = strtmp(4)
- endif
-c strtmp(5) = strtmp(5)*1.0e21/tcon(3,mat(ivel))
- xn=3.5
-c strtmp(4) = (strtmp(5)**((1-xn)/xn))
-c
c scale by material parameters
c
call rheol ( x , v , t , ivel , evisc)
-c evisc(5) = 1.0e0
- strtmp(1) = two * evisc(5) * strtmp(1)
- strtmp(2) = two * evisc(5) * strtmp(2)
- strtmp(3) = evisc(5) * strtmp(3)
-c strtmp(4) =-evisc(5)*alam(mat(ivel))* strtmp(4)
- strtmp(4) =-evisc(5)* strtmp(4)
- strtmp(5) = strtmp(2) - strtmp(4)
- strtmp(6) = evisc(5)
c
-c calculte element "weight"
+ do n=1,6
+ strtmp(n) = 0
+ enddo
c
- rhsl(1) = det(1) * shl(1,1) + det(2) * shl(1,2)
- & + det(3) * shl(1,3) + det(4) * shl(1,4)
+c for temperature depenedent problems, you need stresses at
+c the integration points, not at the center, as is in some
+c old versions of this routine. SDK 7/15/08
c
- rhsl(2) = det(1) * shl(2,1) + det(2) * shl(2,2)
- & + det(3) * shl(2,3) + det(4) * shl(2,4)
+ do intp=1,4
+c normal stress x
+ strtmp(1) = strtmp(1) + two*evisc(intp)*
+ & (shdx(1,intp)*vl(1) + shdx(2,intp)*vl(3)
+ & + shdx(3,intp)*vl(5) + shdx(4,intp)*vl(7))*det(intp)
+c normal stress y
+ strtmp(2) = strtmp(2) + two*evisc(intp)*
+ & (shdy(1,intp)*vl(2) + shdy(2,intp)*vl(4)
+ & + shdy(3,intp)*vl(6) + shdy(4,intp)*vl(8))*det(intp)
+c shear stress
+ strtmp(3) = strtmp(3) + evisc(intp)*
+ & + 0.5*(shdx(1,intp)*vl(2) + shdx(2,intp)*vl(4)
+ & + shdx(3,intp)*vl(6) + shdx(4,intp)*vl(8)
+ & + shdy(1,intp)*vl(1) + shdy(2,intp)*vl(3)
+ & + shdy(3,intp)*vl(5) + shdy(4,intp)*vl(7))*det(intp)
+c pressure
+ strtmp(4) = strtmp(4) +
+ & (shdx(1,intp)*vl(1) + shdx(2,intp)*vl(3)
+ & + shdx(3,intp)*vl(5) + shdx(4,intp)*vl(7))*det(intp)
+ & +(shdy(1,intp)*vl(2) + shdy(2,intp)*vl(4)
+ & + shdy(3,intp)*vl(6) + shdy(4,intp)*vl(8))*det(intp)
c
- rhsl(3) = det(1) * shl(3,1) + det(2) * shl(3,2)
- & + det(3) * shl(3,3) + det(4) * shl(3,4)
+c end intp loop
c
- rhsl(4) = det(1) * shl(4,1) + det(2) * shl(4,2)
- & + det(3) * shl(4,3) + det(4) * shl(4,4)
-c
-c assemble element stress contribution to the node
-c
- do n=1,6
- rhsl(5) = rhsl(1) * strtmp(n)
- rhsl(6) = rhsl(2) * strtmp(n)
- rhsl(7) = rhsl(3) * strtmp(n)
- rhsl(8) = rhsl(4) * strtmp(n)
- stress(n,ien(ivel,1)) = stress(n,ien(ivel,1)) + rhsl(5)
- stress(n,ien(ivel,2)) = stress(n,ien(ivel,2)) + rhsl(6)
- stress(n,ien(ivel,3)) = stress(n,ien(ivel,3)) + rhsl(7)
- stress(n,ien(ivel,4)) = stress(n,ien(ivel,4)) + rhsl(8)
enddo
+ volume=det(1)+det(2)+ det(3) + det(4)
+ do n=1,6
+ strtmp(n) = strtmp(n)/volume
+ enddo
+ strtmp(4) =-evisc(5)*alam(mat(ivel))*strtmp(4)
+c strtmp(4) =-evisc(5)* strtmp(4)
+ strtmp(5) = strtmp(2) - strtmp(4)
+ strtmp(6) = evisc(5)
c
-c calculate pmass done above in fluxke left here for completeness
+c assemble element stress contributions to nodes
c
- tmass(1) = shl(1,1) * det(1) + shl(1,2) * det(2)
- & + shl(1,3) * det(3) + shl(1,4) * det(4)
+ do node = 1,4
c
- tmass(2) = shl(2,1) * det(1) + shl(2,2) * det(2)
- & + shl(2,3) * det(3) + shl(2,4) * det(4)
+c calculte element "weight"
c
- tmass(3) = shl(3,1) * det(1) + shl(3,2) * det(2)
- & + shl(3,3) * det(3) + shl(3,4) * det(4)
-c
- tmass(4) = shl(4,1) * det(1) + shl(4,2) * det(2)
- & + shl(4,3) * det(3) + shl(4,4) * det(4)
-c
+ tmass(node) = det(1)*shl(node,1) + det(2)*shl(node,2)
+ & + det(3)*shl(node,3) + det(4)*shl(node,4)
+ do n=1,6
+ stress(n,ien(ivel,node)) = stress(n,ien(ivel,node))
+ & + tmass(node)*strtmp(n)
+ enddo
+ enddo
smass=tmass(1) + tmass(2) + tmass(3) + tmass(4)
-c
volume=det(1)+det(2)+ det(3) + det(4)
c
- tmass(1) = tmass(1) * volume / smass
- tmass(2) = tmass(2) * volume / smass
- tmass(3) = tmass(3) * volume / smass
- tmass(4) = tmass(4) * volume / smass
+ do node = 1,4
+ tmass(node) = tmass(node) * volume / smass
+ pmass(ien(ivel,node)) = pmass(ien(ivel,node))+tmass(node)
+ enddo
c
- pmass( ien(ivel,1) ) = pmass( ien(ivel,1) ) + tmass(1)
- pmass( ien(ivel,2) ) = pmass( ien(ivel,2) ) + tmass(2)
- pmass( ien(ivel,3) ) = pmass( ien(ivel,3) ) + tmass(3)
- pmass( ien(ivel,4) ) = pmass( ien(ivel,4) ) + tmass(4)
-c
c end loop over elements
c
1000 continue
c
-c
c
do n=1,6
do i=1, numnp
@@ -167,7 +129,7 @@
do n=1,6
do i=1,nodebn
c Hughes for testing purposes
-c stress(n,nb(1,i))=two*stress(n,nb(1,i))-stress(n,nb(2,i))
+ stress(n,nb(1,i))=two*stress(n,nb(1,i))-stress(n,nb(2,i))
c234567890123456789012345678901234567890123456789012345678901234567890
enddo
enddo
@@ -202,4 +164,4 @@
1500 format(' node x1 x2 txx tzz',
& ' txz P ')
2000 format(1x,i7,1x,2(1pe13.7,1x),6(1pe12.5,1x))
- end
+ end
\ No newline at end of file
More information about the cig-commits
mailing list