*" PACKAGE GTSTSPCT(NEW) 時空間パワースペクトル、位相 *" *" PACKAGE GTSTSPCT 時空間パワースペクトル、位相 *" [HIS] 93/11/25 堀之内 第一版（未完成部分多々あり） *" 93/12/08 堀之内 今のところ"空間"軸(SDIM)が x方向であると仮定 *" して軸ファイルがつくられるようになっている。 *" オプション -TIMSQ は正常に働かない。 *" PACKAGE GTSTSPCT(NEW) 時空間パワースペクトル、位相 *" 96/06/02 オプション -steady で定常成分(振動数＝0)も出力する *" 96/06/02 配列詰め替えによりベクトル化率を高める *********************************************************************** PROGRAM GTSTSP * #ifdef SYS_IBMS INCLUDE (GZSIZE) INCLUDE (GTSINC) #else #include "gzsize.F" #include "gtsinc.F" #endif * *" [PARAM] 必要な領域は約 (2*NDIM4+3*IJKDIM)*4 バイト. 確保可能か注意 ** PARAMETER (NDIM4 = 30E5 ) !" 読み込む全データ数の最大値 (+α) ** PARAMETER (NDIM4 = 4.5E6 ) !" 読み込む全データ数の最大値 (+α) ** PARAMETER (NDIM4 = 1E7 ) !" 読み込む全データ数の最大値 (+α) * REAL GDAT4( NDIM4 ) * CHARACTER HHEAD ( NDC )*(NCC) CHARACTER HHEADS( NDC )*(NCC) CHARACTER HHEADW( NDC )*(NCC) REAL GDATA ( IJKDIM ) * CHARACTER HFILE( 1 ) *(NFILN) DATA HFILE / '$GTTMPDIR/gtool.out' / DATA IFILE / 50 / DATA JFILE / 60 / DATA JFILEP/ 61 / * CHARACTER OUT( 1 ) *(NFILN) DATA OUT / '$GTTMPDIR/gtool.out' / CHARACTER OUTP( 1 ) *(NFILN) DATA OUTP/'$GTTMPDIR/gtoolph.out' / LOGICAL APND DATA APND / .FALSE. / * CHARACTER ITEM *(NCC) CHARACTER UNIT *(NCC) CHARACTER TITLE *(NCC*2) CHARACTER DSET *(NCC) DATA ITEM, UNIT, TITLE, DSET /4*' '/ LOGICAL GRESET DATA GRESET / .TRUE. / LOGICAL HELP DATA HELP / .FALSE. / LOGICAL NOPH !" T -> 位相を出力しない DATA NOPH / .FALSE. / LOGICAL STEADY !" T:定常成分(ω＝0)も出力する DATA STEADY / .FALSE. / LOGICAL MIS00 !" T: s=ω=0成分を欠損値にする(強制的にSTEADY=Tにする) DATA MIS00 / .FALSE. / INTEGER SDIM !" Fourier変換される空間次元は何次元目か PARAMETER (SDIM=1) ** LOGICAL TIMSQ !" T -> 入力データは既に時系列化されている ***" 但し、時系列データは、１レコードの最終次元が時間軸であると仮定する ** DATA TIMSQ / .FALSE. / INTEGER SEQ !" 元が空間３次元データの場合、できた4次元データを * !" 第 SEQ 次元以外の3次元データの列として出力する DATA SEQ / 4 / !" デフォルトは、周波数毎のデータの列 LOGICAL RSTRD !" T -> 空間フーリエの前にリニアトレンドを引く DATA RSTRD / .FALSE. / LOGICAL RTTRD !" T -> 時間フーリエの前にリニアトレンドを引く DATA RTTRD / .FALSE. / LOGICAL LOG !" T-> スペクトルのスケーリングタイプを log にする DATA LOG / .FALSE. / LOGICAL NOAX !" T-> 軸ファイルを書かない DATA NOAX / .FALSE. / CHARACTER UTIM *4 !" 共通パラメター 'UTIM': 変えない方が無難 DATA UTIM / 'DAY ' / * NAMELIST /OPTION/ OUT, OUTP, APND, & ITEM, UNIT, TITLE, DSET, GRESET, & HELP, HFILE, NOPH, STEADY, MIS00, & SEQ, RSTRD, RTTRD, LOG, NOAX, UTIM * CHARACTER EDIT *(NCC) CHARACTER ETTL *(NCC) DATA EDIT / 'STPWSP' / DATA ETTL / 's-t power spectr' / * CHARACTER EDITP *(NCC) CHARACTER ETTLP *(NCC) DATA EDITP / 'STPH' / DATA ETTLP / 's-t phase' / * CHARACTER AITMF *(NCC) !" 時間軸ファイル名 DATA AITMF / '@FREQ' / !" その頭の部分 CHARACTER AITMS *(NCC) !" 空間軸ファイル名 DATA AITMS / '@WVNWE' / !" その頭の部分 INTEGER ITIMD !" 入力の何次元目が時間軸か。 DATA ITIMD / 4 / INTEGER ISN DATA ISN / 1 / CHARACTER HFTIME *(NCC) CHARACTER HFTIM1 *(NCC) LOGICAL LFIRST DATA LFIRST /.TRUE./ LOGICAL LINIT DATA LINIT /.TRUE./ CHARACTER HAITM1*(NCC),HAITM2*(NCC),HAITM3*(NCC) CHARACTER HAXSEQ*(NCC) CHARACTER HAXSDM*(NCC) INTEGER IXDIM,IYDIM,IZDIM,IWDIM !" ４次元データの各要素数 INTEGER IRDIM,ITDIM !" 変換に使われる空間、時間の要素数 INTEGER IQDIM,ISDIM !" 時系列は(IQ,IR,IS,IT)DIM の４次元として扱われる INTEGER ITIMDO !" 出力の何次元目が時間軸か LOGICAL LSEQ !" T -> 入力が空間３次元データの(時系)列である。 DATA LSEQ /.FALSE./ * * [WORK] REAL WX(IJKDIM),WT(IJKDIM) !" work(consts) for FFT REAL WY(NDIM4) !" work for FFT REAL SP(NDIM4/2),PH(NDIM4/2) !" power spectrum, phase EQUIVALENCE (WY(1),SP) EQUIVALENCE (WY(NDIM4/2+1),PH) * INTEGER IWT(5),IWX(5) * CALL OPTARG ( 91, 'OPTION', 'HFILE', NOPT, NFILE ) READ (91,OPTION,IOSTAT=IOS) CLOSE(91) IF ( IOS.NE.0 .OR. HELP ) THEN WRITE(6,OPTION) STOP ENDIF * CALL GTOPEN CALL GTSIZE ( HHEAD, IJKDIM ) * CALL GURNTF ( HFILE( 1 ), OUT(1) , '$GTTMPDIR/gtool.in' ) CALL GURNTF ( HFILE( 1 ), OUTP(1) , '$GTTMPDIR/gtool.in' ) * CALL GFROPN ( IFILE , HFILE( 1 ) ) CALL GFOOPN ( JFILE , OUT(1) , APND ) CALL GUNENV( OUT(1),'.',.FALSE. ) IL=LENC(OUT(1)) WRITE (6,*) 'output='//OUT(1)(1:IL) IF(.NOT.NOPH) THEN CALL GFOOPN ( JFILEP, OUTP(1) , APND ) CALL GUNENV( OUTP(1),'.',.FALSE. ) IL=LENC(OUTP(1)) WRITE (6,*) 'output='//OUTP(1)(1:IL) ENDIF * CALL GTCSET ( 'UTIM', UTIM ) * * < read data > * IANUM = 0 1100 CONTINUE CALL GHCOPY( HHEAD,HHEADS) CALL GFREAD O ( HHEAD , GDATA , IEOD , I IFILE , 1 ) * IF ( IEOD .EQ.0 ) THEN IANUM = IANUM + 1 IF (LFIRST) THEN LFIRST = .FALSE. CALL GHCGET( HHEAD, 'TIME' , HFTIM1 ) CALL GHCGET( HHEAD, 'AITM1', HAITM1 ) CALL GHCGET( HHEAD, 'AITM2', HAITM2 ) CALL GHCGET( HHEAD, 'AITM3', HAITM3 ) IF ( HAITM1 .EQ. ' ' ) & CALL MSGDMP( 'E',' MAIN ','ZERO SPACE DIMENSION' ) CALL GUQSIZ I ( HHEAD , O IXSTR , IXEND , IXDIM , O IYSTR , IYEND , IYDIM , O IZSTR , IZEND , IZDIM ) IXYZDM=IXDIM*IYDIM*IZDIM IYZDIM=IYDIM*IZDIM * IF(HAITM2.EQ.' ') THEN ITIMDO = 2 ELSE IF(HAITM3.EQ.' ') THEN ITIMDO = 3 ELSE LSEQ = .TRUE. IF (SEQ.EQ.4) THEN ITIMDO = 4 ELSE ITIMDO = 3 ENDIF ENDIF ENDIF IF ( IANUM*IXDIM*IYDIM*IZDIM .GT. NDIM4 ) & CALL MSGDMP('E',' MAIN ','TOO SMALL "NDIM4"') CALL XYZIRK(GDAT4(1+(IANUM-1)*IXYZDM),GDATA,IYZDIM,IXDIM) IF (MOD(IANUM,2).EQ.0) CALL GHCGET( HHEAD, 'TIME' , HFTIME ) GOTO 1100 ENDIF CALL GHCOPY( HHEADS,HHEAD) CALL GHCOPY( HHEADS,HHEADW) * CALL GHCSET(HHEAD, 'TIME', HFTIME) * * < set dimension > * IWDIM = IANUM * IQDIM = IYZDIM IRDIM = IXDIM ISDIM = 1 ITDIM = 2*(IWDIM/2) !" 偶数に限る. IFDIM = ITDIM/2+1 !" 周波数軸の個数(計算) * IF(MIS00) THEN STEADY=.TRUE. ENDIF IF(STEADY) THEN IFDIMO= ITDIM/2+1 !" 周波数軸の個数(実際の出力) : ω=0 付き ELSE IFDIMO= ITDIM/2 !" 周波数軸の個数(実際の出力) : ω=0 なし ENDIF * IF(ITDIM.NE.IWDIM) THEN WRITE(6,*) '** MESSAGE ** USE ONLY THE FIRST',ITDIM,' DATA ', + ' OF',IWDIM,' DATA' ENDIF * * < Fourier trans, spectrum,phase > * IF (RSTRD) THEN !" 空間リニアトレンドを引く CALL SUBTRD M ( GDAT4, I IQDIM, IRDIM, IRDIM, ISDIM*ITDIM ) ENDIF CALL AFFTRU !" 空間 Fourier M ( GDAT4, W WY, C WX, IWX, I IQDIM, IRDIM, IRDIM, ISDIM*ITDIM, ISN, LINIT ) * IF (RTTRD) THEN !" 時間リニアトレンドを引く CALL SUBTRD M ( GDAT4, I IQDIM*IRDIM*ISDIM, ITDIM, ITDIM, 1 ) ENDIF CALL AFFTRU !" 時間 Fourier M ( GDAT4, W WY, C WT, IWT, I IQDIM*IRDIM*ISDIM, ITDIM, ITDIM, 1, ISN, LINIT ) LINIT = .FALSE. * CALL FST2S2 !" space-time fourier --> space-time power spectrum O ( SP , I GDAT4 , IQDIM, IRDIM, ISDIM, ITDIM ) * IF (.NOT.NOPH) THEN IF ( IQDIM*(IRDIM+1)*ISDIM*IFDIM .GT. NDIM4/2 ) & CALL MSGDMP('E',' MAIN ','TOO Small "NDIM4"') CALL FST2P2 !" space-time fourier --> space-time fourier phase O ( PH , I GDAT4 , IQDIM, IRDIM, ISDIM, ITDIM ) ENDIF * * < header set > * IF ( ITEM .NE. ' ' ) THEN CALL GHCSET( HHEAD , 'ITEM', ITEM ) ENDIF IF ( UNIT .NE. ' ' ) THEN CALL GHCSET( HHEAD , 'UNIT', UNIT ) ELSE CALL GHCGET( HHEAD , 'UNIT', UNIT ) IL = LENZ(UNIT) IF (IL.NE.0) THEN UNIT = '('//UNIT(1:IL)//')|2"' CALL GHCSET( HHEAD , 'UNIT', UNIT ) ENDIF ENDIF IF ( TITLE .NE. ' ' ) THEN CALL GHCSTS( HHEAD , 'TITL', TITLE ) ELSE CALL GHCGTS( HHEAD , 'TITL', TITLE ) IL = LENC(TITLE) TITLE = TITLE(1:IL)//' s-t pw.sp' CALL GHCSTS( HHEAD , 'TITL', TITLE ) ENDIF IF ( DSET .NE. ' ' ) THEN CALL GHCSET( HHEAD , 'DSET', DSET ) ENDIF IF ( GRESET ) THEN CALL GHRSGP( HHEAD ) ENDIF * CALL GHEADD( HHEAD , EDIT , ETTL ) * IL = LENC(AITMS) CALL CHVAL( '(I8)', REAL(IRDIM/2) , AITMS(IL+1:NCC) ) * CALL GHCGET( HHEAD, 'AITM1', HAXSDM) ** IF (HAXSDM(1:4).EQ.'CSIG') THEN ** AITMS(5:6) = 'VT' ** HAXSDM(1:4) = 'CLNS' ** CALL GHCSET( HHEADW, 'AITM1', HAXSDM) ** ENDIF CALL GHCSET( HHEAD, 'AITM1', AITMS) CALL GHPSET( HHEAD, 'ASTR1', 1 ) CALL GHPSET( HHEAD, 'AEND1', IRDIM+1) * IL = LENC(AITMF) CALL CHVAL( '(I8)', REAL(IFDIMO) , AITMF(IL+1:NCC) ) CALL FFTTIM I ( TIMSQ, ITDIM, I HHEAD, HFTIM1, I ITIMDO, O FTIM ) IL = LENC(AITMF) AITMF(IL+1:IL+1) = '-' CALL CHVAL( '(I8)', REAL(FTIM) , AITMF(IL+2:NCC) ) * CALL HEDAXF M ( HHEAD, O HAXSEQ, I AITMF, IFDIMO, ITIMDO, SEQ ,LSEQ ) * CALL GHCOPY( HHEAD,HHEADS) IF ( TITLE .NE. ' ' ) THEN CALL GHCSTS( HHEADS , 'TITL', TITLE ) ELSE CALL GHCGTS( HHEADS , 'TITL', TITLE ) IL = LENC(TITLE) TITLE = TITLE(1:IL)//' s-t phase' CALL GHCSTS( HHEADS , 'TITL', TITLE ) ENDIF CALL GHEADD( HHEADS , EDITP , ETTLP ) CALL GHCSET( HHEADS ,'UNIT', 'RAD') * IF (LOG) THEN CALL GHPSET( HHEAD , 'STYP' , 2 ) ENDIF * * < output > * IF(MIS00) THEN CALL GHPGET ( HHEAD, 'MISS', VMISS ) CALL SMIS00(SP,VMISS,IYZDIM,IXDIM+1,IFDIMO) IF(.NOT.NOPH) THEN CALL SMIS00(PH,VMISS,IYZDIM,IXDIM+1,IFDIMO) ENDIF ENDIF * IF (LSEQ) THEN IXDIMO = IXDIM + 1 IYDIMO = IYDIM IZDIMO = IZDIM IF(SEQ.EQ.1) NSEQ=IXDIMO IF(SEQ.EQ.2) THEN NSEQ=IYDIMO KSEQ=IYSTR ENDIF IF(SEQ.EQ.3) THEN NSEQ=IZDIMO KSEQ=IZSTR ENDIF IF(SEQ.EQ.4) NSEQ=IFDIMO CALL XYZRCV(GDAT4,SP,IXDIM+1,IYZDIM,IFDIM,IFDIMO) DO 100 ISEQ = 1,NSEQ CALL COPY43(SP,GDAT4,SEQ,ISEQ,IXDIMO,IYDIMO,IZDIMO,IFDIMO) IF (SEQ.EQ.4) THEN CALL EADDPR M ( HHEAD , I ISEQ, FTIM, AITMF ) ELSE IF (SEQ.EQ.SDIM) THEN CALL EADDSX M ( HHEAD, I ISEQ, IRDIM , AITMS ) ELSE CALL EADDSQ M ( HHEAD, I HAXSEQ, SEQ, (KSEQ-1)+ISEQ ) ENDIF CALL GFWRIT I ( HHEAD , SP , I JFILE , 1 , 0 ) 100 CONTINUE IF (.NOT.NOPH) THEN CALL XYZRCV(GDAT4,PH,IXDIM+1,IYZDIM,IFDIM,IFDIMO) DO 200 ISEQ = 1,NSEQ CALL COPY43(PH,GDAT4,SEQ,ISEQ,IXDIMO,IYDIMO,IZDIMO,IFDIMO) IF (SEQ.EQ.4) THEN CALL EADDPR M ( HHEADS , I ISEQ, FTIM, AITMF ) ELSE IF (SEQ.EQ.SDIM) THEN CALL EADDSX M ( HHEADS, I ISEQ, IRDIM , AITMS ) ELSE CALL EADDSQ M ( HHEADS, I HAXSEQ, SEQ, ISEQ ) ENDIF CALL GFWRIT I ( HHEADS , PH , I JFILEP, 1 , 0 ) 200 CONTINUE ENDIF ELSE CALL XYZRCV(GDAT4,SP,IXDIM+1,IYZDIM,IFDIM,IFDIMO) CALL GFWRIT I ( HHEAD ,GDAT4 , I JFILE , 1 , 0 ) IF (.NOT.NOPH) THEN CALL XYZRCV(GDAT4,PH,IXDIM+1,IYZDIM,IFDIM,IFDIMO) CALL GFWRIT I ( HHEADS ,GDAT4 , I JFILEP , 1 , 0 ) ENDIF ENDIF * * < make axis file > * IF (.NOT.NOAX) THEN CALL JIKFRQ I ( HHEAD, AITMF, FTIM, I IFDIMO, ITDIM ) * IF (HAXSDM(1:4).EQ.'GLON') THEN CALL JIKSPC I ( HHEAD, AITMS, I IRDIM, SDIM ) * ELSE IF(HAXSDM(1:4).EQ.'CLNS') THEN * CALL JIKVRT * I ( HHEAD, HHEADW, HAXSDM, * I IRDIM, SDIM ) * ENDIF ENDIF * STOP END ********************************************************************* SUBROUTINE COPY(X,Y,N) REAL X(N),Y(N) DO 10 I=1,N X(I) = Y(I) 10 CONTINUE RETURN END ********************************************************************* * copy 4dim array --> 3dim array SUBROUTINE COPY43(X,Y,KD,ID,N1,N2,N3,N4) REAL X(*),Y(*) IF (KD.EQ.1) THEN CALL COP431(X,Y,ID,N1,N2,N3,N4) ELSE IF (KD.EQ.2) THEN CALL COP432(X,Y,ID,N1,N2,N3,N4) ELSE IF (KD.EQ.3) THEN CALL COP433(X,Y,ID,N1,N2,N3,N4) ELSE IF (KD.EQ.4) THEN CALL COP434(X,Y,ID,N1,N2,N3,N4) ELSE CALL MSGDMP('E','COPY43','INVALID "KD"') ENDIF RETURN END ********************************************************************* * copy 4dim array --> 3dim array SUBROUTINE COP431(X,Y,ID,N1,N2,N3,N4) REAL X(N2,N3,N4),Y(N1,N2,N3,N4) DO 10 I4 = 1,N4 DO 10 I3 = 1,N3 DO 10 I2 = 1,N2 X(I2,I3,I4) = Y(ID,I2,I3,I4) 10 CONTINUE RETURN END ********************************************************************* * copy 4dim array --> 3dim array SUBROUTINE COP432(X,Y,ID,N1,N2,N3,N4) REAL X(N1,N3,N4),Y(N1,N2,N3,N4) DO 10 I4 = 1,N4 DO 10 I3 = 1,N3 DO 10 I1 = 1,N1 X(I1,I3,I4) = Y(I1,ID,I3,I4) 10 CONTINUE RETURN END ********************************************************************* * copy 4dim array --> 3dim array SUBROUTINE COP433(X,Y,ID,N1,N2,N3,N4) REAL X(N1,N2,N4),Y(N1,N2,N3,N4) DO 10 I4 = 1,N4 DO 10 I2 = 1,N2 DO 10 I1 = 1,N1 X(I1,I2,I4) = Y(I1,I2,ID,I4) 10 CONTINUE RETURN END ********************************************************************* * copy 4dim array --> 3dim array SUBROUTINE COP434(X,Y,ID,N1,N2,N3,N4) REAL X(N1,N2,N3),Y(N1,N2,N3,N4) DO 10 I3 = 1,N3 DO 10 I2 = 1,N2 DO 10 I1 = 1,N1 X(I1,I2,I3) = Y(I1,I2,I3,ID) 10 CONTINUE RETURN END ********************************************************************* SUBROUTINE JIKFRQ I ( HHEAD, AITMF, FTIM, I IFDIM, ITDIM ) * #ifdef SYS_IBMS INCLUDE (GZSIZE) !" NCC, NDC, NED #else #include "gzsize.F" !" NCC, NDC, NED #endif CHARACTER HHEAD ( NDC )*(NCC) CHARACTER AITMF *(NCC) CHARACTER HUTIM *4 !"共通パラメタ UTIM (default:'DAY') CHARACTER HUNIT *6 PARAMETER ( NAX=1024 ) REAL AX ( NAX ) REAL DAX ( NAX ) CHARACTER HAX (NDC)*(NCC) * CALL GTCGET ( 'UTIM', HUTIM ) IL = LENC(HUTIM) HUNIT = '1/'//HUTIM(1:IL) CALL GHPGET ( HHEAD, 'MISS', VMISS ) * DF = 1. / FTIM DO 10 I = 1,IFDIM IF(IFDIM.NE.ITDIM/2) THEN AX(I) = DF*(I-1) ELSE AX(I) = DF*I ENDIF DAX(I) = 1. 10 CONTINUE * CALL GHCSET( HAX, 'AITM1', AITMF ) CALL GHPSET( HAX, 'AEND1', IFDIM ) * CALL GHPACA O ( HAX, I 'IAXLOC', AITMF, I 'frequency', HUNIT , I IFDIM, I 'UR4' , VMISS , I VMISS , VMISS , VMISS , VMISS , 1 ) * CALL GUSIAX I ( HAX , AX , DAX , I 1 ) * CALL GUWIAX I ( HAX ) * RETURN END ********************************************************************* SUBROUTINE JIKSPC I ( HHEAD, AITMS, I NS , IAX ) * #ifdef SYS_IBMS INCLUDE (GZSIZE) !" NCC, NDC, NED #else #include "gzsize.F" !" NCC, NDC, NED #endif CHARACTER HHEAD ( NDC )*(NCC) CHARACTER AITMS *(NCC) CHARACTER HAITM *(NCC) PARAMETER ( NAX=1024 ) REAL AX ( NAX ) REAL DAX ( NAX ) CHARACTER HAX (NDC)*(NCC) * CALL GHPGET ( HHEAD, 'MISS', VMISS ) CALL GHPSET( HAX, 'AEND1', NS+1 ) * IF (IAX.EQ.1) THEN CALL GHCGET( HHEAD, 'AITM1', HAITM ) ELSE IF (IAX.EQ.2) THEN CALL GHCGET( HHEAD, 'AITM2', HAITM ) ELSE IF (IAX.EQ.3) THEN CALL GHCGET( HHEAD, 'AITM3', HAITM ) ENDIF * DO 10 I = -NS/2,NS/2 AX(I+NS/2+1) = 1.*I DAX(I+NS/2+1) = 1. 10 CONTINUE * CALL GHPACA O ( HAX, I 'IAXLOC', AITMS, I 'wave number', ' ', I NS+1, I 'UR4' , VMISS , I VMISS , VMISS , VMISS , VMISS , 1 ) * CALL GUSIAX I ( HAX , AX , DAX , I 1 ) * CALL GUWIAX I ( HAX ) * RETURN END ********************************************************************* SUBROUTINE JIKVRT I ( HHEAD, HHEADW, HAXSDM, I NS , IAX ) * #ifdef SYS_IBMS INCLUDE (GZSIZE) !" NCC, NDC, NED #else #include "gzsize.F" !" NCC, NDC, NED #endif CHARACTER HHEAD ( NDC )*(NCC) CHARACTER HAITM *(NCC) CHARACTER HAXSDM *(NCC) PARAMETER ( NAX=1024 ) REAL WGT ( NAX ) REAL AX ( NAX ) REAL DAX ( NAX ) CHARACTER HAX (NDC)*(NCC) CHARACTER HHEADO( NDC )*(NCC) LOGICAL LREQ * CALL GHPGET ( HHEAD, 'MISS', VMISS ) CALL GHPSET( HAX, 'AEND1', NS+1 ) * IF (IAX.EQ.1) THEN CALL GHCGET( HHEAD, 'AITM1', HAITM ) ELSE IF (IAX.EQ.2) THEN CALL GHCGET( HHEAD, 'AITM2', HAITM ) ELSE IF (IAX.EQ.3) THEN CALL GHCGET( HHEAD, 'AITM3', HAITM ) ENDIF * CALL GTSIZE ( HHEADO , NAX ) CALL GUQAXV I ( HHEADW, IAX, 'WGT' , O HHEADO,WGT, IEOD ) * W = WGT(1) S = WGT(1) DO 5 I = 2,NS W2 = WGT(I) IF(.NOT.LREQ(W,W2)) & CALL MSGDMP('W','JIKVRT','AXIS IN NOT EQUALLY DIVDED') W = W2 S = S + WGT(I) 5 CONTINUE * DS = -1 / S * DO 10 I = -NS/2,NS/2 AX(I+NS/2+1) = DS*I DAX(I+NS/2+1) = 1. 10 CONTINUE * CALL GHPACA O ( HAX, I 'IAXLOC', HAITM, I 'vertical wn', ' ', I NS+1, I 'UR4' , VMISS , I VMISS , VMISS , VMISS , VMISS , 1 ) * CALL GUSIAX I ( HAX , AX , DAX , I 1 ) * CALL GUWIAX I ( HAX ) * RETURN END *********************************************************************** * GET TIME FOR TIME FFT * SUBROUTINE FFTTIM I ( LTIMSQ, ITDIM, I HHEAD, HFTIM1, I ITIMD, O FTIM ) * #ifdef SYS_IBMS INCLUDE (GZSIZE) !" NCC, NDC, NED #else #include "gzsize.F" !" NCC, NDC, NED #endif * [INPUT] LOGICAL LTIMSQ !" T -> 入力データは既に時系列化されている。 INTEGER ITDIM !" 変換に使われる時間の要素数 CHARACTER HHEAD ( NDC )*(NCC) !" (for .NOT.LTMSQ) CHARACTER HFTIM1 *(NCC) !" 始めの時刻 (for .NOT.LTMSQ) INTEGER ITIMD !" 時間軸は何次元目か (for LTIMSQ) * [OUTPUT] REAL FTIM !" フーリエ変換に使われた時間 (単位 HUTIM) * [WORK] CHARACTER HUTIM *4 !"共通パラメタ UTIM (default:'DAY') となる CHARACTER HUTIMF *4 !" (for LTIMSQ) CHARACTER HAW *5 !" (for LTIMSQ) CHARACTER HHEADW( NDC )*(NCC) PARAMETER ( MAXAX = 1024 ) REAL AXIS(MAXAX) * IF (.NOT.LTIMSQ) THEN CALL GUQTIM ( HHEAD , RTF , HUTIM ) !" 'DAY'?単位で時刻を読む CALL GHCOPY( HHEAD,HHEADW) CALL GHCSET(HHEADW, 'TIME', HFTIM1) CALL GUQTIM ( HHEADW , RTI , HUTIM ) !" 'DAY'?単位で時刻を読む FTIM = ( RTF - RTI ) * ITDIM / (ITDIM-1) ELSE IF (ITDIM.GT.MAXAX) & CALL MSGDMP('E','FFTTIM','TOO LARGE TIME SEQUENCE') CALL GUQAXV I ( HHEAD, ITIMD, 'LOC' , O HHEADW,AXIS , IEOD ) HAW = 'ASTR#' WRITE(HAW(5:5),'(I1)') ITIMD CALL GHPGET(HHEAD,HAW,IASTR) RTI = AXIS(IASTR) RTF = AXIS(IASTR+ITDIM) FTIM = ( RTF - RTI ) * (ITDIM+1) / ITDIM CALL GHCGET ( HHEADW, 'UNIT', HUTIMF ) CALL GTCGET ( 'UTIM', HUTIM ) IF ( HUTIMF(1:1) .NE. HUTIM(1:1) ) THEN IT = FTIM CALL GUCT2R I ( IT , HUTIMF, HUTIM , 1 , O FTIM ) ENDIF ENDIF * RETURN END *********************************************************************** SUBROUTINE HEDAXF M ( HHEAD, O HAXSEQ, I AITMF, IFDIM, ITIMDO, ISEQ ,LSEQ ) * #ifdef SYS_IBMS INCLUDE (GZSIZE) !" NCC, NDC, NED #else #include "gzsize.F" !" NCC, NDC, NED #endif * [MODIFY] CHARACTER HHEAD ( NDC )*(NCC) * [OUTPUT] CHARACTER HAXSEQ *(NCC) * [INPUT] CHARACTER AITMF *(NCC) !" 時間軸ファイル名 INTEGER ITIMDO !" 出力の何次元目が時間軸か INTEGER ISEQ !" 元が空間３次元データの場合、できた4次元データを * !" 第 ISEQ 次元以外の3次元データの列として出力する。 LOGICAL LSEQ !" T -> 入力が空間３次元データの(時系)列である。 * [WORK] CHARACTER*(NCC) CAITM,CASTR,CAEND CHARACTER*(NCC) HAITM,HASTR,HAEND * CAITM = 'AITM#' CASTR = 'ASTR#' CAEND = 'AEND#' * IF (ISEQ.NE.4) THEN WRITE(CAITM(5:5),'(I1)') ISEQ CALL GHCGET( HHEAD, CAITM, HAXSEQ ) ENDIF * IF (LSEQ .AND. ISEQ.LT.3) THEN DO 10 KS = ISEQ+1,3 WRITE(CAITM(5:5),'(I1)') KS WRITE(CASTR(5:5),'(I1)') KS WRITE(CAEND(5:5),'(I1)') KS CALL GHCGET( HHEAD, CAITM, HAITM ) CALL GHPGET( HHEAD, CASTR, HASTR ) CALL GHPGET( HHEAD, CAEND, HAEND ) WRITE(CAITM(5:5),'(I1)') KS-1 WRITE(CASTR(5:5),'(I1)') KS-1 WRITE(CAEND(5:5),'(I1)') KS-1 CALL GHCSET( HHEAD, CAITM, HAITM ) CALL GHPSET( HHEAD, CASTR, HASTR ) CALL GHPSET( HHEAD, CAEND, HAEND ) 10 CONTINUE ENDIF * WRITE(CAITM(5:5),'(I1)') ITIMDO WRITE(CASTR(5:5),'(I1)') ITIMDO WRITE(CAEND(5:5),'(I1)') ITIMDO IF (ITIMDO.NE.4) THEN CALL GHCSET( HHEAD, CAITM, AITMF ) CALL GHPSET( HHEAD, CASTR, 1 ) CALL GHPSET( HHEAD, CAEND, IFDIM) ENDIF * RETURN END ********************************************************************* *" 編集記述子追加 (周期) SUBROUTINE EADDPR M ( HHEAD, I IFRQ, FTIM, AITMF ) * #ifdef SYS_IBMS INCLUDE (GZSIZE) !" NCC, NDC, NED #else #include "gzsize.F" !" NCC, NDC, NED #endif CHARACTER HHEAD ( NDC )*(NCC) * INTEGER IFRQ !" 何番目か REAL FTIM !" フーリエ変換に使われた時間 (単位 'UTIM') CHARACTER AITMF *(NCC) !" 時間軸ファイル名 * CHARACTER EDIT *(NCC) CHARACTER ETTL *(NCC) CHARACTER HUTIM *4 !"共通パラメタ UTIM (default:'DAY') となる CHARACTER HNUM *(NCC) REAL PERIOD INTEGER IENUM LOGICAL OFIRST DATA OFIRST / .TRUE. / SAVE IENUM , OFIRST * DF = 1. / FTIM FRQ = DF*IFRQ PERIOD = 1/FRQ * EDIT = AITMF * CALL GTCGET ( 'UTIM', HUTIM ) CALL CHVAL( 'B' , PERIOD , HNUM ) IL = LENC(HNUM) ETTL = 'PERIOD'//HNUM(1:IL)//HUTIM * IF (OFIRST) THEN CALL GHQENM ( HHEAD, IENUM ) OFIRST = .FALSE. ENDIF CALL GHESET( HHEAD , EDIT , ETTL , IENUM+1 ) * RETURN END ********************************************************************* *" 編集記述子追加 (zonal w.n.) SUBROUTINE EADDSX M ( HHEAD, I ISX, IRDIM , AITMS ) * #ifdef SYS_IBMS INCLUDE (GZSIZE) !" NCC, NDC, NED #else #include "gzsize.F" !" NCC, NDC, NED #endif CHARACTER HHEAD ( NDC )*(NCC) * INTEGER ISX !" 何番目か INTEGER IRDIM !" フーリエ変換に使われたデータ数 CHARACTER AITMS *(NCC) !" 軸ファイル名 * CHARACTER EDIT *(NCC) CHARACTER ETTL *(NCC) CHARACTER HNUM *(NCC) REAL WN INTEGER IENUM LOGICAL OFIRST DATA OFIRST / .TRUE. / SAVE IENUM , OFIRST * WN = -IRDIM/2 + (ISX-1) * EDIT = AITMS * CALL CHVAL( '(I8)' , WN , HNUM ) ETTL = 'WVN'//HNUM * IF (OFIRST) THEN CALL GHQENM ( HHEAD, IENUM ) OFIRST = .FALSE. ENDIF CALL GHESET( HHEAD , EDIT , ETTL , IENUM+1 ) * RETURN END *********************************************************************** *" 編集記述子追加 SUBROUTINE EADDSQ M ( HHEAD, I HAXSEQ, KD, IS ) * #ifdef SYS_IBMS INCLUDE (GZSIZE) !" NCC, NDC, NED #else #include "gzsize.F" !" NCC, NDC, NED #endif CHARACTER HHEAD ( NDC )*(NCC) * CHARACTER HAXSEQ *(NCC) INTEGER KD !" 何次元目か INTEGER IS !" 何番目か * CHARACTER HEDIT *(NCC) !" 編集略記号 CHARACTER HETTL *(NCC) !" 編集タイトル CHARACTER HED *2 LOGICAL OFIRST DATA OFIRST / .TRUE. / INTEGER IENUM SAVE IENUM , OFIRST * IF(KD.EQ.1) THEN HED = 'XS' ELSE IF(KD.EQ.2) THEN HED = 'YS' ELSE IF(KD.EQ.3) THEN HED = 'ZS' ELSE HED = 'SS' ENDIF * CALL GUEAXN I ( HAXSEQ, HED , '-select', I IS , IS , O HEDIT , HETTL ) IF (OFIRST) THEN CALL GHQENM ( HHEAD, IENUM ) OFIRST = .FALSE. ENDIF CALL GHESET( HHEAD , HEDIT , HETTL , IENUM+1 ) * RETURN END ************************************************************************ SUBROUTINE XYZIRK(ROUT,RIN,L,M) REAL ROUT(L,M) REAL RIN (M,L) DO 100 J=1,M DO 100 I=1,L ROUT(I,J)=RIN(J,I) 100 CONTINUE RETURN END ************************************************************************ SUBROUTINE XYZRCV(ROUT,RIN,L,M,NI,NO) REAL ROUT(L,M,NO) REAL RIN (M,L,NI) IF(NI.EQ.NO) THEN DO 100 K=1,NI DO 100 I=1,L DO 100 J=1,M ROUT(I,J,K)=RIN(J,I,K) 100 CONTINUE ELSE IF(NI.EQ.NO+1) THEN DO 200 K=1,NO DO 200 I=1,L DO 200 J=1,M ROUT(I,J,K)=RIN(J,I,K+1) 200 CONTINUE ELSE CALL MSGDMP('E','XYZRCV','NI.NE.NO .AND. NI.NE.NO+1') ENDIF RETURN END ************************************************************************ SUBROUTINE SMIS00(R,VMISS,L,M,N) REAL R(L,M,N) DO 10 I=1,L R(I,M/2+1,1)=VMISS 10 CONTINUE RETURN END