C**** LPDIF.FOR DIFerence LinePlot1 minus LinePlot2 2000/02/13 C**** C**** Compile into executable module: FCE90 LPDIF.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=8) CHARACTER*80 FILEIN(2),TITLEI(2),XAXISI(2),YAXISI(2) CHARACTER*72 TITLEO CHARACTER*256 DATA,LABELI(2) INTEGER*4 LINM(2),IINM(2), JMIN(0:LMAX),JMAX(0:LMAX) LOGICAL*4 QY(LMAX) REAL*4 SY(LMAX), YOUT(IMAX,LMAX) COMMON /INCB/ YIN(IMAX,0:LMAX,2) COMMON /OUTCB/ XOUT(IMAX) EQUIVALENCE (YOUT,YIN(1,1,1)) C**** IF(IARGC().ne.2) GO TO 800 C**** C**** Read in LinePlot input files C**** DO 150 K=1,2 CALL GETARG (K,FILEIN(K)) OPEN (1,FILE=FILEIN(K),STATUS='OLD',ERR=810) READ (1,910) TITLEI(K) READ (1,910) XAXISI(K) READ (1,910) YAXISI(K) READ (1,910) LABELI(K) C**** Parse column labels LINM(K) = 1 + LMAX CALL PARSE (LABELI(K),256,LINM(K),JMIN,JMAX) LINM(K) = LINM(K) - 1 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) LOUTM = MIN(LINM(1),LINM(2)) C**** C**** Subtract LinePlot2 from LinePlot1 C**** DO 320 L=1,LOUTM IOUTM = IINM(1) DO 310 I=1,IOUTM 310 XOUT(I) = YIN(I,0,1) 320 CALL ADTERM (IMAX, IOUTM ,XOUT ,YIN(1,L,1), * IINM(2),YIN(1,0,2),YIN(1,L,2), -1.) C**** C**** Kill leading and trailing data lines with all undefined Y values C**** 500 DO 510 L=1,LOUTM 510 IF(YOUT(1,L).ne.-999999.) GO TO 540 C**** All Y values are undefined for leading line IF(IOUTM.le.1) GO TO 851 IOUTM = IOUTM - 1 DO 520 I=1,IOUTM XOUT(I) = XOUT(I+1) DO 520 L=1,LOUTM 520 YOUT(I,L) = YOUT(I+1,L) GO TO 500 C**** 540 DO 550 L=1,LOUTM 550 IF(YOUT(IOUTM,L).ne.-999999.) GO TO 600 C**** All Y values are undefined for trailing line IF(IOUTM.le.1) GO TO 851 IOUTM = IOUTM - 1 GO TO 540 C**** C**** Write lineplot output file C**** 600 NTOK = 1 CALL PARSE (TITLEI(1)(51:80),30,NTOK,JMIN(1),JMAX(1)) CALL PARSE (TITLEI(2)(51:80),30,NTOK,JMIN(2),JMAX(2)) TITLEO = TITLEI(1)(1:50) // TITLEI(1)(50+JMIN(1):50+JMAX(1)) // * '-' // TITLEI(2)(50+JMIN(2):50+JMAX(2)) CALL WRITLP (6, IMAX,IOUTM,LOUTM, * TITLEO,XAXISI(1),YAXISI(1),LABELI(1), XOUT,YOUT) GO TO 999 C**** 800 WRITE (0,*) ' Usage: LPDIF LinePlot1 LinePlot2 2000/02/13' WRITE (0,*) ' console receives LinePlot1 minus LinePlot2' GO TO 999 810 WRITE (0,*) ' Error opening lineplot file: ',FILEIN(K)(1:48) STOP 810 851 WRITE (0,*) ' Lineplot output file was not created. ', * ' It would have had no valid Y values.' STOP 851 C**** 910 FORMAT (A) 920 FORMAT (1X,A6,1X,A72) 921 FORMAT (' Files Xaxis ',8(' Line',I1,2X) / * ' ----- ----- ',8(' -----',2X)) 922 FORMAT (' File',A1,1X,9(2X,A6)) 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).gt.XJ(J)+.0005) GO TO 400 IF(XI(I).ge.XJ(J)-.0005) GO TO 300 C**** C**** XI(I) < XJ(J) C**** 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 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**** WRITE (IUNIT,900) TRIM(TITLE) WRITE (IUNIT,900) TRIM(XAXIS) WRITE (IUNIT,900) TRIM(YAXIS) WRITE (IUNIT,900) TRIM(LABEL) 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,8A10) 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