C**** MERGE.FOR MERGE several LinePlot files 2006/01/26 C**** C**** Compile into executable module: FCE MERGE.FOR C**** C**** Indicies: I = numbered data rows in each LinePlot file C**** L = numbered Y columns in each LinePlot file C**** K = numbered LinePlot input files C**** PARAMETER (IMAX=1800,LMAX=16,KMAX=8) CHARACTER*80 FILEIN(KMAX),TITLEI(KMAX),XAXISI(KMAX),YAXISI(KMAX) CHARACTER*80 TITLEO,XAXISO,YAXISO CHARACTER*256 DATA CHARACTER*32 LABIN(0:LMAX,KMAX), LABOUT(0:LMAX) INTEGER*4 LINM(KMAX),IINM(KMAX),IOUTM(LMAX), * JMIN(0:LMAX),JMAX(0:LMAX) COMMON /INCB/ YIN(IMAX,0:LMAX,KMAX) COMMON /OUTCB/ XOUT(IMAX,LMAX),YOUT(IMAX,LMAX) C**** KFILES = IARGC() IF(KFILES.le.0) GO TO 800 IF(KFILES.gt.KMAX) KFILES = KMAX C**** C**** Read in LinePlot input files C**** DO 150 K=1,KFILES CALL GETARG (K,FILEIN(K)) OPEN (1,FILE=FILEIN(K),STATUS='OLD',ERR=801) READ (1,910) TITLEI(K) READ (1,910) XAXISI(K) READ (1,910) YAXISI(K) READ (1,910) DATA C**** Parse column labels LINM(K) = 1 + LMAX CALL PARSE (DATA,256,LINM(K),JMIN,JMAX) LINM(K) = LINM(K) - 1 DO 110 L=0,LINM(K) 110 LABIN(L,K) = DATA(JMIN(L):JMAX(L)) C**** Read in data, replace missing data with -999999. DO 130 I=1,IMAX READ (1,910,END=140) DATA NTOK = 1 + LINM(K) CALL PARSE (DATA,256,NTOK,JMIN,JMAX) IF(NTOK.le.0) GO TO 140 NTOK = NTOK - 1 DO 120 L=0,NTOK READ (DATA(JMIN(L):JMAX(L)),*,IOSTAT=IOSQ) YIN(I,L,K) 120 IF(IOSQ.ne.0) YIN(I,L,K) = -999999. IF(YIN(I,0,K).eq.-999999.) GO TO 140 DO 130 L=NTOK+1,LINM(K) 130 YIN(I,L,K) = -999999. 140 IINM(K) = I - 1 150 CLOSE (1) C**** C**** Write column labels to console C**** DO 210 K=1,KFILES WRITE (6,920) 'TITLE:',TITLEI(K) WRITE (6,920) 'Xaxis:',XAXISI(K) WRITE (6,920) 'Yaxis:',YAXISI(K) 210 WRITE (6,910) WRITE (6,921) (L,L=1,8) DO 220 K=1,KFILES 220 WRITE (6,922) CHAR(Z'40'+K),(LABIN(L,K),L=0,LINM(K)) C**** C**** Calculate data and read in labels for output columns C**** WRITE (6,930) DO 410 L=1,LMAX 310 WRITE (6,931) L READ (5,910) DATA NMAX = LEN_TRIM(DATA) IF(NMAX.le.0) GO TO 420 IOUTM(L) = 0 ADDCON = 0. 320 IF(NMAX.le.0) GO TO 380 C**** Decode current term NF = SCAN(DATA(1:NMAX),'ABCDEFGHabcdefgh',.true.) NS = SCAN(DATA(1:NMAX),'+-',.true.) NT = SCAN(DATA(1:NMAX),'*',.true.) IF(NT.gt.NS) GO TO 340 IF(NF.gt.NS) GO TO 330 C**** Current term is a signed constant READ (DATA(NS+1:NMAX),*,ERR=350) FAC IF(DATA(NS:NS).eq.'-') FAC = -FAC ADDCON = ADDCON + FAC NMAX = NS - 1 GO TO 320 C**** Current term is a signed data record 330 KIN = ICHAR(DATA(NF:NF)) - Z'40' IF(KIN.ge.Z'20') KIN = KIN - Z'20' IF(KIN.gt.KFILES) GO TO 360 READ (DATA(NF+1:NMAX),*,ERR=370) LIN IF(LIN.le.0 .or. LIN.gt.LINM(KIN)) GO TO 370 FAC = 1. IF(DATA(NS:NS).eq.'-') FAC = -1. CALL ADTERM (IMAX, IOUTM(L),XOUT(1,L),YOUT(1,L), * IINM(KIN),YIN(1,0,KIN),YIN(1,LIN,KIN), FAC) NMAX = NS - 1 GO TO 320 C**** Current term is a signed number times a data record 340 KIN = ICHAR(DATA(NF:NF)) - Z'40' IF(KIN.ge.Z'20') KIN = KIN - Z'20' IF(KIN.gt.KFILES) GO TO 360 READ (DATA(NF+1:NMAX),*,ERR=370) LIN IF(LIN.le.0 .or. LIN.gt.LINM(KIN)) GO TO 370 READ (DATA(NS+1:NT-1),*,ERR=350) FAC IF(DATA(NS:NS).eq.'-') FAC = -FAC CALL ADTERM (IMAX, IOUTM(L),XOUT(1,L),YOUT(1,L), * IINM(KIN),YIN(1,0,KIN),YIN(1,LIN,KIN), FAC) NMAX = NS - 1 GO TO 320 C**** Error reading multiplicative factor 350 WRITE (6,*) ' Unable to decipher multiplicative factor: ', * DATA(NS+1:NMAX) GO TO 310 C**** File letter exceeds those received from command line 360 WRITE (6,*) ' File letter exceeds files received from command ', * 'line: ',DATA(NF:NF) GO TO 310 C**** Error reading column number integer or column number is invalid 370 WRITE (6,*) ' Column number exceeds input columns or is ', * 'undecipherable: ',DATA(NF+1:NMAX) GO TO 310 C**** Column has no X coordinates: most likely a constant number 380 IF(IOUTM(L).gt.0) GO TO 390 WRITE (6,*) ' Column has no X coordinates.' WRITE (6,*) *' To create a column with constant number, try: A1 - A1 + number' GO TO 310 C**** Add ADDCON to YOUT values 390 DO 400 I=1,IOUTM(L) 400 IF(YOUT(I,L).ne.-999999.) YOUT(I,L) = YOUT(I,L) + ADDCON C**** Determine output column label IF(L.eq.1) TITLEO = TITLEI(KIN) IF(L.eq.1) XAXISO = XAXISI(KIN) IF(L.eq.1) YAXISO = YAXISI(KIN) IF(L.eq.1) LABOUT(0) = LABIN(0,KIN) WRITE (6,940) LABIN(LIN,KIN) READ (5,910) LABOUT(L) IF(LABOUT(L).eq.' ') LABOUT(L) = LABIN(LIN,KIN) 410 continue 420 LOUTM = L - 1 IF(LOUTM.le.0) GO TO 842 C**** C**** Interpolate each output column so that all have same X coodinates C**** IF(LOUTM.le.1) GO TO 600 C**** Interpolate columns (XOUT(L),YOUT(1)) into (XOUT(1),YOUT(1)) DO 510 L=2,LOUTM 510 CALL INTERP (IMAX,IOUTM(1),XOUT(1,1),YOUT(1,1), * IOUTM(L),XOUT(1,L)) C**** Interpolate columns (XOUT(1),YOUT(L)) into (XOUT(L),YOUT(L)) DO 520 L=2,LOUTM 520 CALL INTERP (IMAX,IOUTM(L),XOUT(1,L),YOUT(1,L), * IOUTM(1),XOUT(1,1)) C**** C**** Kill leading and trailing data rows with all undefined Y values C**** 600 DO 610 L=1,LOUTM 610 IF(YOUT(1,L).ne.-999999.) GO TO 640 C**** All Y values are undefined for leading row IF(IOUTM(1).le.1) GO TO 861 IOUTM(1) = IOUTM(1) - 1 DO 620 I=1,IOUTM(1) XOUT(I,1) = XOUT(I+1,1) DO 620 L=1,LOUTM 620 YOUT(I,L) = YOUT(I+1,L) GO TO 600 C**** 640 DO 650 L=1,LOUTM 650 IF(YOUT(IOUTM(1),L).ne.-999999.) GO TO 700 C**** All Y values are undefined for trailing line IF(IOUTM(1).le.1) GO TO 861 IOUTM(1) = IOUTM(1) - 1 GO TO 640 C**** C**** Write LinePlot output file C**** 700 OPEN (2,FILE='MERGE.LP') C**** Calculate column labels DO 710 L=1,LOUTM JMAX(L) = LEN_TRIM(LABOUT(L)) DO 710 I=2,JMAX(L)-1 710 IF(LABOUT(L)(I:I).eq.' ') LABOUT(L)(I:I) = '_' WRITE (DATA,971) LABOUT(0),(LABOUT(L)(1:JMAX(L)),L=1,LOUTM) C**** CALL WRITLP (2, IMAX,IOUTM(1),LOUTM, TITLEO,XAXISO,YAXISO,DATA, * XOUT,YOUT) CLOSE (2) GO TO 999 C**** 800 WRITE (0,*) ' Usage: MERGE LinePlot1 [LinePlot2 ... LinePlot8] ', * ' 2006/01/26' WRITE (0,*) ' Create new LinePlot file from linear combinations', * ' of old LinePlot columns' GO TO 999 801 WRITE (0,*) ' Error opening LinePlot file: ',FILEIN(K)(1:48) STOP 801 842 WRITE (0,*) ' Lineplot output file was not created. ', * ' No output columns were defined.' STOP 842 861 WRITE (0,*) ' Lineplot output file was not created. ', * ' It would have had no valid Y values.' STOP 861 C**** 910 FORMAT (A) 920 FORMAT (1X,A6,1X,A72) 921 FORMAT (' Files Xaxis',8(' Col',I2) / * ' ----- -----',8(' -----')) 922 FORMAT (' File',A1,1X,17(2X,A5)) 930 FORMAT (/ *' Define output columns as linear combinations of input columns.'/ *' For example ... Colum1: 2*A2 - 273.16 + B1 - 2.*C1') 931 FORMAT (/ ' Colum',I1,': ',$) 940 FORMAT (' Default column label: ',A32 / * ' Enter column label: ',$) 971 FORMAT (A8,16(2X,A)) 999 END SUBROUTINE ADTERM (IMAX, IM,XI,YI, JM,XJ,YJ, FAC) C**** C**** ADTERM adds the term YJ*FAC to YI, augments the X coordinates, C**** and performs linear interpolation for newly defined Y coordinates C**** C**** Input: IMAX = dimension size of XI, YI, XJ and YJ C**** IM = number of (XI,YI) pairs C**** XI,YI = arrays of X and Y coordinates C**** JM = number of (XJ,YJ) pairs C**** XJ,YJ = arrays of X and Y coordinates C**** FAC = constant factor that will multiply YJ C**** C**** Output: IM = number of augmented (XI,JI) pairs C**** XI = augmented X coordinates, superset of XI and XJ C**** YI = YI + YJ*FAC C**** REAL*4 XI(IMAX),YI(IMAX), XJ(IMAX),YJ(IMAX) C**** IF(JM.gt.0) GO TO 10 WRITE (0,*) ' Number of (XJ,YJ) pairs received by ADTERM is 0.' RETURN 10 IF(IM.gt.0) GO TO 30 C**** C**** XI and YI arrays are empty, copy (XJ,YJ*FAC) to (XI,YI) C**** IM = JM DO 20 I=1,IM XI(I) = XJ(I) YI(I) = YJ(I)*FAC 20 IF(YJ(I).eq.-999999.) YI(I) = -999999. RETURN C**** C**** XI and YI arrays are not empty C**** 30 I = 1 J = 1 100 IF(XI(I)-XJ(J)) 200,300,400 C**** C**** XI(I) < XJ(J) C**** 200 IF(J.le.1) GO TO 240 C**** XJ(J-1) < XI(I) < XJ(J): insert YI+YJ*FAC into YI IF(YI(I).eq.-999999. .or. YJ(J-1).eq.-999999. .or. * YJ(J).eq.-999999.) then YI(I) = -999999. else YI(I) = YI(I) + ((XJ(J)-XI(I))*YJ(J-1) + * (XI(I)-XJ(J-1))*YJ(J))*FAC / (XJ(J)-XJ(J-1)) endif I = I + 1 IF(I.le.IM) GO TO 100 C**** XI(IM) < XJ(J): insert (XJ,-999999) pairs into (XI,YI) arrays 220 DO 230 JJ=J,JM IF(IM.ge.IMAX) GO TO 800 IM = IM + 1 XI(IM) = XJ(JJ) 230 YI(IM) = -999999. RETURN C**** XI(I) < XJ(1): insert (XI,-999999) into (XI,YI) 240 YI(I) = -999999. I = I + 1 IF(I.le.IM) GO TO 100 GO TO 220 C**** C**** XI(I) = XJ(J): insert (XJ,YJ*FAC) into (XI,YI) arrays C**** 300 IF(YI(I).eq.-999999. .or. YJ(J).eq.-999999.) then YI(I) = -999999. else YI(I) = YI(I) + YJ(J)*FAC endif I = I + 1 J = J + 1 IF(I.gt.IM) GO TO 220 IF(J.gt.JM) GO TO 420 GO TO 100 C**** C**** XI(I) > XJ(J): move (XI,YI) pairs up one I C**** 400 IF(IM.ge.IMAX) GO TO 800 DO 410 II=IM,I,-1 XI(II+1) = XI(II) 410 YI(II+1) = YI(II) IM = IM + 1 IF(I.le.1) GO TO 440 C**** XI(I-1) < XJ(J) < XI(I): insert (XJ,YI+YJ*FAC) into (XI,YI) IF(YI(I).eq.-999999. .or. YI(I-1).eq.-999999. .or. * YJ(J).eq.-999999.) then YI(I) = -999999. else YI(I) = YJ(J)*FAC + ((XI(I)-XJ(J))*YI(I-1) + * (XJ(J)-XJ(I-1))*YI(I)) / (XI(I)-XI(I-1)) endif XI(I) = XJ(J) I = I + 1 J = J + 1 IF(J.le.JM) GO TO 100 C**** XI(I) > XJ(JM): reset YI to -999999 for all original I's 420 DO 430 II=I,IM 430 YI(II) = -999999. RETURN C**** XI(1) > XJ(J) 440 XI(1) = XJ(J) YI(1) = -999999. I = I + 1 J = J + 1 IF(J.le.JM) GO TO 100 GO TO 420 C**** 800 WRITE (0,*) ' Number of X coordinates in ADTERM exceeds ', * 'internal dimensions of:',IMAX RETURN END SUBROUTINE INTERP (IMAX, IM,XI,YI, JM,XJ) C**** C**** INTERP adds new coordinates (XJ,YI) into existing pairs (XI,YI) C**** performing linear interpolation for newly defined YI coordinates C**** C**** Input: IMAX = dimension size of XI, YI, XJ and YJ C**** IM = number of (XI,YI) pairs C**** XI,YI = arrays of X and Y coordinates C**** JM = number of XJ coordinates C**** XJ = array of X coordinates C**** C**** Output: IM = number of augmented (XI,JI) pairs C**** XI = augmented X coordinates, superset of XI and XJ C**** YI = YI (possibly via linear interpolation) C**** REAL*4 XI(IMAX),YI(IMAX), XJ(IMAX) C**** IF(IM.gt.0) GO TO 10 WRITE (0,*) ' Number of (XI,YI) pairs received by INTERP is 0.' RETURN 10 IF(JM.gt.0) GO TO 20 WRITE (0,*) ' Number of (XJ,YJ) pairs received by INTERP is 0.' RETURN C**** C**** XI and YI arrays are not empty C**** 20 I = 1 J = 1 100 IF(XI(I)-XJ(J)) 200,300,400 C**** C**** XI(I) < XJ(J) C**** C**** XI(I) < XJ(J): get next I 200 I = I + 1 IF(I.le.IM) GO TO 100 C**** XI(IM) < XJ(J): insert (XJ,-999999) pairs into (XI,YI) arrays 220 DO 230 JJ=J,JM IF(IM.ge.IMAX) GO TO 800 IM = IM + 1 XI(IM) = XJ(JJ) 230 YI(IM) = -999999. RETURN C**** C**** XI(I) = XJ(J): get next I and J C**** 300 I = I + 1 J = J + 1 IF(I.gt.IM) GO TO 220 IF(J.gt.JM) GO TO 420 GO TO 100 C**** C**** XI(I) > XJ(J): move (XI,YI) pairs up one I C**** 400 IF(IM.ge.IMAX) GO TO 800 DO 410 II=IM,I,-1 XI(II+1) = XI(II) 410 YI(II+1) = YI(II) IM = IM + 1 IF(I.le.1) GO TO 440 C**** XI(I-1) < XJ(J) < XI(I): insert (XJ,interpolated YI) into (XI,YI) IF(YI(I).eq.-999999. .or. YI(I-1).eq.-999999.) then YI(I) = -999999. else YI(I) = ((XI(I)-XJ(J))*YI(I-1) + * (XJ(J)-XI(I-1))*YI(I)) / (XI(I)-XI(I-1)) endif XI(I) = XJ(J) I = I + 1 J = J + 1 IF(J.le.JM) GO TO 100 C**** XI(I) > XJ(JM): no more XJs to insert 420 RETURN C**** XI(1) > XJ(J) 440 XI(1) = XJ(J) YI(1) = -999999. I = I + 1 J = J + 1 IF(J.le.JM) GO TO 100 GO TO 420 C**** 800 WRITE (0,*) ' Number of X coordinates in ADTERM exceeds ', * 'internal dimensions of:',IMAX RETURN END SUBROUTINE WRITLP (IUNIT, KMAX,KM,LM, TITLE,XAXIS,YAXIS,LABEL,X,Y) C**** C**** WRITLP writes the data lines of a LinePlot file to unit IUNIT, C**** using a reasonable format. Y values that are -999999. are C**** replaced with * in the output file. C**** C**** Input: KMAX = row dimension of X and Y arrays C**** KM = number of used rows in X and Y arrays C**** LM = number of Y columns C**** TITLE = title character string C**** XAXIS = X axis label character string C**** YAXIS = Y axis label character string C**** LABEL = column labels character string C**** X = one dimensional array of X coordinates C**** Y = two dimensional array of Y values C**** C**** Output: IUNIT = unit to send LinePlot data to, exclude titles C**** CHARACTER TITLE*(*),XAXIS*(*),YAXIS*(*),LABEL*(*), C(0:16)*10 REAL*4 X(KMAX),Y(KMAX,LM) C**** C**** Write 4 header lines of LinePlot file C**** ITITLE = LEN_TRIM (TITLE) IXAXIS = LEN_TRIM (XAXIS) IYAXIS = LEN_TRIM (YAXIS) ILABEL = LEN_TRIM (LABEL) WRITE (IUNIT,900) TITLE(1:ITITLE) WRITE (IUNIT,900) XAXIS(1:IXAXIS) WRITE (IUNIT,900) YAXIS(1:IYAXIS) WRITE (IUNIT,900) LABEL(1:ILABEL) C**** C**** Calculate maximum absolute Y value C**** ABSMAX = 0. DO 1 L=1,LM DO 1 K=1,KM IF(Y(K,L).eq.-999999.) GO TO 1 IF(ABSMAX.lt.ABS(Y(K,L))) ABSMAX = ABS(Y(K,L)) 1 continue KPOWER = LOG10(ABSMAX) + 2 IF(KPOWER.lt.1) KPOWER = 1 IF(KPOWER.gt.8) KPOWER = 8 C**** C**** Write data lines C**** GO TO (700,600,500,400,300,200,100),KPOWER C**** 1000000 < ABSMAX: output numbers have no decimal places DO 20 K=1,KM DO 10 L=1,LM IF(Y(K,L).ne.-999999.) WRITE (C(L),901) NINT(Y(K,L)) 10 IF(Y(K,L).eq.-999999.) C(L) = ' *' 20 WRITE (IUNIT,902) X(K),(C(L),L=1,LM) RETURN C**** 100000 < ABSMAX < 1000000: output numbers have 1 decimal place 100 DO 120 K=1,KM DO 110 L=1,LM IF(Y(K,L).ne.-999999.) WRITE (C(L),911) Y(K,L) 110 IF(Y(K,L).eq.-999999.) C(L) = ' * ' 120 WRITE (IUNIT,902) X(K),(C(L),L=1,LM) RETURN C**** 10000 < ABSMAX < 100000: output numbers have 2 decimal places 200 DO 220 K=1,KM DO 210 L=1,LM IF(Y(K,L).ne.-999999.) WRITE (C(L),921) Y(K,L) 210 IF(Y(K,L).eq.-999999.) C(L) = ' * ' 220 WRITE (IUNIT,902) X(K),(C(L),L=1,LM) RETURN C**** 1000 < ABSMAX < 10000: output numbers have 3 decimal places 300 DO 320 K=1,KM DO 310 L=1,LM IF(Y(K,L).ne.-999999.) WRITE (C(L),931) Y(K,L) 310 IF(Y(K,L).eq.-999999.) C(L) = ' * ' 320 WRITE (IUNIT,902) X(K),(C(L),L=1,LM) RETURN C**** 100 < ABSMAX < 1000: output numbers have 4 decimal places 400 DO 420 K=1,KM DO 410 L=1,LM IF(Y(K,L).ne.-999999.) WRITE (C(L),941) Y(K,L) 410 IF(Y(K,L).eq.-999999.) C(L) = ' * ' 420 WRITE (IUNIT,902) X(K),(C(L),L=1,LM) RETURN C**** 10 < ABSMAX < 100: output numbers have 5 decimal places 500 DO 520 K=1,KM DO 510 L=1,LM IF(Y(K,L).ne.-999999.) WRITE (C(L),951) Y(K,L) 510 IF(Y(K,L).eq.-999999.) C(L) = ' * ' 520 WRITE (IUNIT,902) X(K),(C(L),L=1,LM) RETURN C**** 1 < ABSMAX < 10: output numbers have 6 decimal places 600 DO 620 K=1,KM DO 610 L=1,LM IF(Y(K,L).ne.-999999.) WRITE (C(L),961) Y(K,L) 610 IF(Y(K,L).eq.-999999.) C(L) = ' * ' 620 WRITE (IUNIT,902) X(K),(C(L),L=1,LM) RETURN C**** ABSMAX < 1: output numbers have 7 decimal places 700 DO 720 K=1,KM DO 710 L=1,LM IF(Y(K,L).ne.-999999.) WRITE (C(L),961) Y(K,L) 710 IF(Y(K,L).eq.-999999.) C(L) = ' * ' 720 WRITE (IUNIT,902) X(K),(C(L),L=1,LM) RETURN C**** 900 FORMAT (A) 901 FORMAT (1X,I9) 902 FORMAT (F8.3,16A10) 911 FORMAT (1X,F9.1) 921 FORMAT (1X,F9.2) 931 FORMAT (1X,F9.3) 941 FORMAT (1X,F9.4) 951 FORMAT (1X,F9.5) 961 FORMAT (1X,F9.6) 971 FORMAT (1X,F9.7) END SUBROUTINE PARSE (S,LENS,NTOK,JMIN,JMAX) C**** C**** PARSE parses a character string S of length LENS into its first C**** NTOK tokens. JMIN(n) and JMAX(n) are the first and last C**** characters of S for the n-th token. C**** INTEGER*4 JMIN(*),JMAX(*) CHARACTER S*(*) N = 0 K = 0 10 IF(N.ge.NTOK) GO TO 70 N = N + 1 C**** Skip over blank characters between tokens 20 IF(K.ge.LENS) GO TO 50 K = K + 1 IF(S(K:K).eq.' ') GO TO 20 JMIN(N) = K C**** Contiguous nonblank characters are part of the token 30 IF(K.ge.LENS) GO TO 40 K = K + 1 IF(S(K:K).ne.' ') GO TO 30 JMAX(N) = K-1 GO TO 10 40 JMAX(N) = K N = N + 1 C**** Remaining tokens were not included in the string 50 DO 60 N0=N,NTOK JMIN(N0) = 0 60 JMAX(N0) = 0 N = N - 1 70 NTOK = N RETURN END