//ALTINT JOB NOTIFY=6461 /*SERVICE -4 /*JOBPARM T=10,L=99,R=4096,PRINT=ALL //S1 EXEC FORTVCLG,REGION=4096K,PARM.FORT='OPTIMIZE(3)' //FORT.SYSIN DD * C----------------------------------------------------------------------- C*********************************************************************** C * C PROGRAM NAME : ALTINT * C FUNCTION : TWO-DIMENSIONAL NUMERICAL INTEGRATION OVER * C RESIDUAL SEA SURFACE HEIGHTS USING DIFFERENT * C ALTIMETRIC TRUNCATION KERNELS * C COMPILER : VS FORTRAN VERSION 2 * C AUTHOR : SPIROS PAGIATAKIS * C HISTORY : FEBRUARY 4, 1986 - VERSION 1.0 * C REFERENCE : U.N.B. TECHNICAL REPORT # * C*********************************************************************** C IMPLICIT REAL*8(A-H,O-Z) REAL *8 KER(100),KERNEL(200) CHARACTER *5 MODE DIMENSION ALAT0(16), ALON0(16), RING(200), RRING(200), X(500), * B(500), C(500), D(500) D0 = 3.14159265D0 / 180.D0 C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C INPUT: ALAT0 = NORTH LATITUDE OF POINT OF INTEREST IN DECIMAL DEGREES C ALON0 = EAST LONGITUDE OF POINT OF INTEREST IN DECIMAL DEGREES C TRING = THICKNESS OF RINGS IN DECIMAL MINUTES OF ARC C PSIMAX = MAXIMUM PSI FOR INTEGRATION IN DECIMAL DEGREES C RHO = RADIUS OF STOKES' RING IN DECIMAL DEGREES C PHIMIN = LATITUDE OF SOUTH-WEST CORNER OF AVAILABLE GRID IN C DECIMAL DEGREES C ALAMIN = EAST LONGITUDE OF SOUTH-WEST CORNER OF AVAILABLE GRID C IN DECIMAL DEGREES C PHIMAX = LATITUDE OF NORTH-EAST CORNER OF AVAILABLE GRID C IN DECIMAL DEGREES C ALAMAX = EAST LONGITUDE OF NORTH-EAST CORNER OF AVAILABLE GRID C IN DECIMAL DEGREES C MODE = POINT OR GRID COMPUTATION (MODE=POINT OR MODE=GRID) C W1 = LATITUDE OF SOUTH-WEST CORNER OF GRID OF INTEREST C W2 = EAST LONGITUDE OF SOUTH-WEST CORNER OF GRID OF INTEREST C W3 = LATITUDE OF NORTH-EAST CORNER OF GRID OF INTEREST C W4 = EAST LONGITUDE OF NORTH-EAST CORNER OF GRID OF INTEREST C STEP = GRID STEP IN DECIMAL MINUTES OF ARC C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C READ MODE OF COMPUTATION C READ (5,1000) MODE C C READ NUMBER OF POINTS OF INTEREST C READ (5,1001) NPNTS C C READ LATITUDE AND LONGITUDE OF POINTS OF INTEREST AND THICKNESS C OF RINGS C DO 10 I = 1, NPNTS READ (5,1002) ALAT0(I), ALON0(I) IF(ALON0(I) .GT. 180) ALON0(I) = 180.D0 - ALON0(I) 10 CONTINUE READ (5,1003) TRING C C READ MAXIMUM PSI FOR INTEGRATION, RHO OF CONSTANT RING C READ (5,1004) PSIMAX, RHO WRITE(6,1009) PSIMAX, RHO, TRING TRING = TRING / 60.D0 C C READ NUMBER OF KERNEL VALUES , AND KERNEL C READ (5,1005) NK DO 20 I = 1, NK READ (5,1006) X(I), KERNEL(I) 20 CONTINUE C C READ AREA BOUNDARIES OF AVAILABLE GRIDDED DATA C READ (5,1007) PHIMIN, ALAMIN, PHIMAX, ALAMAX C C READ WINDOW OF GRID REQUESTED AND STEP C READ (5,1008) W1, W2, W3, W4, STEP STEP = STEP / 60.D0 C C CHECK IF DATA ON DISK IS AVAILABLE FOR POINTS OF INTEREST C DO 30 I = 1, NPNTS CALL PSI0(I, ALAT0, ALON0, ALAT0, ALAMIN, PSI) IF(PSI .LT. PSIMAX) CALL ERROR (3) CALL PSI0(I, ALAT0, ALON0, PHIMAX, ALON0, PSI) IF(PSI .LT. PSIMAX) CALL ERROR(3) CALL PSI0(I, ALAT0, ALON0, ALAT0, ALAMAX, PSI) IF(PSI .LT. PSIMAX) CALL ERROR(3) CALL PSI0(I, ALAT0, ALON0, PHIMIN, ALON0, PSI) IF(PSI .LT. PSIMAX) CALL ERROR(3) 30 CONTINUE C C CHECK IF REQUESTED GRID FALLS WITHIN AVAILABLE DATA GRID C W11 = W1 - PSIMAX - 1.0D0 W22 = W2 - PSIMAX - 1.0D0 W33 = W3 + PSIMAX + 1.0D0 W44 = W4 + PSIMAX + 1.0D0 IF(W11 .LT. PHIMIN) CALL ERROR(1) IF(W22 .LT. ALAMIN) CALL ERROR(1) IF(W33 .GT. PHIMAX) CALL ERROR(1) IF(W44 .GT. ALAMAX) CALL ERROR(1) C C CHECK IF DIMENSIONS OF GRID ARE INTEGER MULTIPLES OF STEP C D1 = DMOD ((W3-W1), STEP) D2 = DMOD ((W4-W2), STEP) C C COMPUTE OUTER AND MEAN RADIUS OF INNERMOST RING C DR = RHO - TRING / 2.D0 DM = DMOD(DR,TRING) DM = DABS(DM - TRING) IF(DM .LT. 0.000001D0) GO TO 40 IRING = DR / TRING RING(1) = DR - TRING * IRING C WRITE(6,*) 'RING(1)= ',RING(1) RRING(1) = RING(1) / 2.D0 GO TO 50 40 RING(1) = TRING RRING(1) = TRING / 2.D0 C C COMPUTE OUTER AND MEAN RADII OF INTERMEDIATE RINGS C IDENTIFY RING THAT IS IDENTICAL WITH STOKES' RING C 50 NRINGS = PSIMAX / TRING + 2.D0 C DO 60 I = 2, NRINGS RING(I) = RING(1) + (I-1) * TRING IF(RING(I) .GE. PSIMAX) GO TO 70 RRING(I) = RING(I) - TRING/2.D0 DRING = DABS(RRING(I) - RHO) IF(DRING .LT. 0.000001D0) KK = I 60 CONTINUE C C COMPUTE RADIUS OF THE OUTERMOST RING C 70 RING(I) = PSIMAX RRING(I) = RING(I-1) + (PSIMAX - RING(I-1)) / 2.D0 NRINGS = I C C COMPUTE CUBIC SPLINE FOR KERNEL C NI = NK / 3 + 1 C DO 80 II = 1,NI B(II) = 0.D0 C(II) = 0.D0 D(II) = 0.D0 80 CONTINUE CALL SPLIND (NK, X, KERNEL, B, C, D) C C COMPUTE KERNEL VALUES FOR EACH OF THE RINGS C U = RRING(1) SEV = SEVALD(NK, U, X, KERNEL, B, C, D) KER(1) = SEV*DSIN(U*D0)*2.D0*U*D0/(DCOS(U*D0)-DCOS(RHO*D0)) DO 90 I = 2, NRINGS-1 IF(I .EQ. KK) GO TO 90 U = RRING(I) SEV = SEVALD(NK, U, X, KERNEL, B, C, D) KER(I) = SEV*DSIN(U*D0)*TRING*D0/(DCOS(U*D0)-DCOS(RHO*D0)) 90 CONTINUE U = RRING(NRINGS) SEV = SEVALD(NK, U, X, KERNEL, B, C, D) KER(NRINGS) = SEV*DSIN(U*D0)*(PSIMAX - RING(NRINGS-1))*D0 * /(DCOS(U*D0)-DCOS(RHO*D0)) C C OPEN DIRECT ACCESS FILE ON DISK (THIS IS AN IBM DIRECT ACCESS FILE) C OPEN (1,STATUS='OLD',ACCESS='DIRECT',FORM='UNFORMATTED', * RECL=4088) C IF(MODE .EQ. 'POINT') CALL POINTM(NPNTS,NRINGS,RING, * PHIMIN,ALAMIN,PSIMAX,ALAT0,ALON0,KER,KK) C IF(MODE .EQ. 'GRID ') CALL GRIDM (NRINGS,RING, * PSIMAX,W1,W2,W3,W4,KER,KK,STEP,PHIMIN,ALAMIN) C 1000 FORMAT (5A) 1001 FORMAT (I4) 1002 FORMAT (F9.5,1X,F9.5) 1003 FORMAT (F8.4) 1004 FORMAT (F9.5,1X,F9.5) 1005 FORMAT (I4) 1006 FORMAT (F13.7,F17.10) 1007 FORMAT (4(F9.5,1X)) 1008 FORMAT (5(F9.5,1X)) 1009 FORMAT (1X,'MAXIMUM PSI = ',F10.5,/,1X,'RHO = ', * F10.5,/,1X,'RING THICKNESS = ',F10.5, * ///,2X,'POINT',2X,'LATITUDE',2X,'LONGITUDE',1X, * ' DN(M)',1X,'SIGMA(M)') STOP END C*********************************************************************** SUBROUTINE PSI0 (I,LATA, LONA, LAT, LON, PSI) C----------------------------------------------------------------------- C SUBROUTINE PSI0 DETERMINES THE SPHERICAL ANGLE PSI, DEFINED BY THE C THE POINT OF INTEREST AND THE DUMMY POINT. C C INPUT C ===== C LATA = LATITUDE OF THE POINT OF INTEREST (NORTH OR SOUTH) C LONA = LONGITUDE OF THE POINT OF INTEREST (EAST OR WEST) C LAT = LATITUDE OF THE DUMMY POINT C LON = LONGITUDE OF THE DUMMY POINT C C OUTPUT C ====== C PSI = GEOCENTRIC SPHERICAL ANGLE IN DEGREES. C C NOTE: C ======= C ALL THE INPUT ANGLES AS WELL AS THE VARIABLES ARE IN DEGREES. C C NORTH LATITUDE AND EAST LONGITUDE ARE CONSIDERED POSITIVE WHILE C SOUTH LATITUDE AND WEST LONGITUDE ARE CONSIDERED NEGATIVE. C---------------------------------------------------------------------- IMPLICIT REAL *8 (A-H,O-Z) REAL *8 LATA(16), LONA(16), LAT, LON D0 = 3.14159265D0 / 180.D0 P1 = DSIN(LATA(I)*D0) * DSIN(LAT*D0) P2 = DCOS(LATA(I)*D0) * DCOS(LAT*D0) * DCOS((LON-LONA(I))*D0) PSI = DARCOS(P1 + P2) / D0 RETURN END C********************************************************************** SUBROUTINE ERROR(IER) INTEGER IER,IPR DATA IPR /6/ C---------------------------------------------------------------------- C FUNCTION: ERROR DETECTS WHETHER ERROR HAS OCCURED C C ARGUMENT: IER = ERROR INDEX C----------------------------------------------------------------------- IF(IER .EQ. 1) GO TO 1 IF(IER .EQ. 2) GO TO 2 IF(IER .EQ. 3) GO TO 3 IF(IER .EQ. 4) GO TO 4 IF(IER .EQ. 5) GO TO 5 IF(IER .EQ. 6) GO TO 6 1 WRITE (6,1001) STOP 2 WRITE (6,1002) STOP 3 WRITE (6,1003) STOP 4 WRITE (6,1004) STOP 5 WRITE (6,1005) STOP 6 WRITE (6,1006) STOP 1001 FORMAT (' ',///,1X,'*** ERROR ***',//,1X,'REQUESTED GRID IS OUT OF * RANGE OF AVAILABLE ALTIMETRY DATA') 1002 FORMAT (' ',///,1X,'*** ERROR ***',//,1X,'REQUESTED DIMENSIONS OF * GRID ARE NOT INTEGER MULTIPLES OF STEP') 1003 FORMAT (' ',///,1X,'*** ERROR ***',//,1X,'ALTIMETRY DATA NOT AVAIL *ABLE. CHECK POSITION OF POINTS OF INTEREST') 1004 FORMAT (' ',///,1X,'*** ERROR ***',//,1X,'THERE ARE NO POINTS IN * ONE OR MORE RINGS') 1005 FORMAT (' ',///,1X,'*** ERROR ***',//,1X,'FILE NUMBER NOT AVAILABL *E ON DISK') 1006 FORMAT (' ',///,'*** ERROR ***') END C*********************************************************************** SUBROUTINE POINTM (NPNTS,NRINGS,RING,PHIMIN,ALAMIN, * PSIMAX,ALAT0,ALON0,KER,KK) IMPLICIT REAL*8(A-H,O-Z) REAL *8 KER(100),NUM(200) REAL *4 AN(212,541),SIGMAN(212,541) DIMENSION ALAT0(16), ALON0(16), RING(200), * NNN(200),VAR(200), SAN(200), STD(200) D0 = 3.14159265D0 / 180.D0 C C READ RECORDS FROM DISK C DO 130 MM = 1, NPNTS ALMIN = 180.D0 - ALON0(MM) - PSIMAX - 3.0D0 ALMAX = 180.D0 - ALON0(MM) + PSIMAX + 3.0D0 N = 6.D0 * (ALAT0(MM) - PSIMAX - PHIMIN) NREC = 12.D0 * PSIMAX + 8.D0 NREC = NREC + N DO 10 I = N, NREC READ (1,REC=I) (AN(I,J), SIGMAN(I,J), J = 1,541) IF(AN(I,J) .GT. 100.0) AN(I,J) = 0.0 10 CONTINUE C C INITIALIZE VARIABLES C COMPUTE N(RHO,ALFA) AND N(PSI,ALFA) C K = 0 SUM = 0.D0 DO 20 I = 1, NRINGS SAN(I) = 0.D0 VAR(I) = 0.D0 NUM(I) = 0.D0 NNN(I) = 0 20 CONTINUE C DO 90 I = N, NREC PHI = PHIMIN + (I - 1) / 6.D0 DO 80 J = 1, 541 ALAM = ALAMIN + (J - 1) / 6.D0 IF(ALAM .LT. ALMIN .OR. ALAM .GT. ALMAX) GO TO 80 IF(ALAM .GT. 180.D0) ALAM = 180.D0 - ALAM P1 = DSIN(ALAT0(MM)*D0) * DSIN(PHI*D0) P2 = DCOS(ALAT0(MM)*D0) * DCOS(PHI*D0) * * DCOS((ALAM-ALON0(MM))*D0) PSI = DARCOS(P1+P2) / D0 IF(PSI .GT. PSIMAX) GO TO 80 C 30 DO 40 K = 2, NRINGS K1 = NRINGS - K + 1 K2 = K1 + 1 IF(PSI .GT. RING(K1) .AND. PSI .LE. RING(K2)) GO TO 50 40 CONTINUE C GO TO 60 50 NUM(K2) = NUM(K2) + DCOS(PHI*D0) NNN(K2) = NNN(K2) + 1 SAN(K2) = SAN(K2) + DBLE(AN(I,J))*DCOS(PHI*D0) VAR(K2) = VAR(K2) + DBLE(SIGMAN(I,J) * SIGMAN(I,J)) GO TO 80 60 IF(PSI .LE. RING(1)) GO TO 70 GO TO 80 70 NUM(1) = NUM(1) + DCOS(PHI*D0) NNN(1) = NNN(1) + 1 SAN(1) = SAN(1) + DBLE(AN(I,J))*DCOS(PHI*D0) VAR(1) = VAR(1) + DBLE(SIGMAN(I,J) * SIGMAN(I,J)) 80 CONTINUE 90 CONTINUE C DO 100 IJ = 1 ,NRINGS IF(NUM(IJ) .LE. 0.000001D0) CALL ERROR(4) WRITE(6,*) 'RING',IJ,'# POINTS', NNN(IJ),'SUM N',SAN(IJ) 100 CONTINUE C SUM = SAN(KK) / NUM(KK) STDSUM = DSQRT(VAR(KK)) / NUM(KK) SAN(1) = SAN(1) / NUM(1) STD(1) = DSQRT(VAR(1)) / NUM(1) SAN(1) = (SAN(1)-SUM) * KER(1) STD(1) = KER(1) * DSQRT(STD(1)*STD(1) + STDSUM*STDSUM) C DO 110 L = 2, NRINGS IF(L .EQ. KK) GO TO 110 SAN(L) = SAN(L) / NUM(L) C WRITE(6,*) 'RING',L,'AVERAGE N', SAN(L) SAN(L) = (SAN(L)-SUM) * KER(L) STD(L) = KER(L) * DSQRT(STD(L)*STD(L) + STDSUM*STDSUM) 110 CONTINUE C DN = 0.D0 VARDN = 0.D0 DO 120 M = 1, NRINGS IF(M .EQ. KK) GO TO 120 DN = DN + SAN(M) C VARDN = VARDN + (SEV*SEV) * (STD(M)*STD(M)) 120 CONTINUE C DN = 0.5D0 * DN DNT=SAN(1)-DN STDDN = 0.5D0 * DSQRT(VARDN) IF(ALON0(MM) .LT. 0) ALON0(MM) = 180.D0 - ALON0(MM) C C PRINT RESULTS C WRITE(6,1001) MM, ALAT0(MM), ALON0(MM), DNT, STDDN 130 CONTINUE 1001 FORMAT (' ',2X,I3,1X,F10.5,1X,F10.5,1X,F6.3,1X,F6.3) STOP END C*********************************************************************** SUBROUTINE GRIDM(NRINGS,RING,PSIMAX,W1,W2,W3,W4, * KER,KK,STEP,PHIMIN,ALAMIN) IMPLICIT REAL*8(A-H,O-Z) REAL *8 KER(100),NUM(200) REAL *4 AN(211,511),SIGMAN(211,511) DIMENSION RING(200), * NNN(200),VAR(200), SAN(200), STD(200) D0 = 3.141592653589793D0 / 180.D0 C C READ RECORDS FROM DISK C N = 6.D0 * (W1 - PSIMAX - PHIMIN) NREC = 6.D0 * (2.D0 * PSIMAX + W3 - W1) + 8.D0 NREC = NREC + N C DO 10 I = N, NREC READ (1,REC=I) (AN(I,J), SIGMAN(I,J), J = 1,511) IF(ABS(AN(I,J)).GT.100.0) AN(I,J)= 1.D12 10 CONTINUE C C INITIALIZE VARIABLES C COMPUTE N(RHO,ALFA) AND N(PSI,ALFA) C GLAT = W1 GLON = W2 15 ALMIN = GLON -(PSIMAX/DCOS(GLAT*D0)+.2D-3) ALMAX = GLON +(PSIMAX/DCOS(GLAT*D0)+.2D-3) CLM1 = GLAT - PSIMAX CLM2 = GLAT + PSIMAX IF(GLON .GE. 180.D0) GLON = 180.D0 - GLON C DO 20 I = 1, NRINGS SAN(I) = 0.D0 VAR(I) = 0.D0 NUM(I) = 0.D0 NNN(I) = 0 20 CONTINUE C DO 90 I = N, NREC PHI = PHIMIN + (I - 1) / 6.D0 IF(PHI .LT. CLM1 .OR. PHI .GT. CLM2) GO TO 90 DO 80 J = 1, 511 ALAM = ALAMIN + (J - 1) / 6.D0 IF(ALAM .LT. ALMIN .OR. ALAM .GT. ALMAX) GO TO 80 IF(ALAM .GT. 180.D0) ALAM = 180.D0 - ALAM P1 = DSIN(GLAT*D0) * DSIN(PHI*D0) P2 = DCOS(GLAT*D0) * DCOS(PHI*D0) * * DCOS((ALAM-GLON)*D0) PSI = DARCOS(P1+P2) / D0 IF(PSI .GT. PSIMAX) GO TO 80 DO 40 K = 2, NRINGS K1 = NRINGS - K + 1 K2 = K1 + 1 IF(PSI .GT. RING(K1) .AND. PSI .LT. RING(K2)) GO TO 50 40 CONTINUE GO TO 60 50 NUM(K2) = NUM(K2) + DCOS(PHI*D0) NNN(K2) = NNN(K2) + 1 SAN(K2) = SAN(K2) + DBLE(AN(I,J))*DCOS(PHI*D0) VAR(K2) = VAR(K2) + DBLE(SIGMAN(I,J) * SIGMAN(I,J)) GO TO 80 60 IF(PSI .LE. RING(1)) GO TO 70 GO TO 80 70 NUM(1) = NUM(1) + DCOS(PHI*D0) NNN(1) = NNN(1) + 1 SAN(1) = SAN(1) + DBLE(AN(I,J))*DCOS(PHI*D0) VAR(1) = VAR(1) + DBLE(SIGMAN(I,J) * SIGMAN(I,J)) 80 CONTINUE 90 CONTINUE C DO 100 IJ = 1 ,NRINGS IF(NUM(IJ) .LE. 0.000001D0) CALL ERROR(4) 100 CONTINUE C SUM = SAN(KK) / NUM(KK) STDSUM = DSQRT(VAR(KK)) / NUM(KK) SAN(1) = SAN(1) / NUM(1) C STD(1) = DSQRT(VAR(1)) / NUM(1) SSS = SAN(1) SAN(1) = (SAN(1)-SUM) * KER(1) STD(1) = KER(1) * DSQRT(STD(1)*STD(1) + STDSUM*STDSUM) C DO 110 L = 2, NRINGS IF(L .EQ. KK) GO TO 110 SAN(L) = SAN(L) / NUM(L) STD(L) = DSQRT(VAR(L)) / NUM(L) SAN(L) = (SAN(L)-SUM) * KER(L) STD(L) = KER(L) * DSQRT(STD(L)*STD(L) + STDSUM*STDSUM) 110 CONTINUE C DN = 0.D0 VARDN = 0.D0 DO 120 M = 1, NRINGS IF(M .EQ. KK) GO TO 120 DN = DN + SAN(M) VARDN = VARDN + (STD(M)*STD(M)) 120 CONTINUE C DN = 0.5D0 * DN DNT=SSS-DN STDDN = 0.5D0 * DSQRT(VARDN) IF(GLON .LT. 0) GLON = 180.D0 - GLON C C PRINT RESULTS C WRITE(6,1001) GLAT, GLON, DNT, STDDN GLON = GLON + STEP IF(GLON .GT. W4) GO TO 130 GO TO 15 130 GLAT = GLAT + STEP IF(GLAT .GT. W3) GO TO 140 GLON = W2 GO TO 15 1001 FORMAT (' ',6X,F10.5,1X,F10.5,1X,F10.3,1X,F6.3) 140 STOP END C********************************************************************** SUBROUTINE SPLIND(N,X,Y,B,C,D) INTEGER N REAL*8 X(N),Y(N),B(N),C(N),D(N) INTEGER NM1,IB,I REAL*8 T NM1=N-1 IF(N.LT.2) RETURN IF (N.LT.3) GO TO 50 D(1)=X(2)-X(1) C(2)=(Y(2)-Y(1))/D(1) DO 10 I=2,NM1 D(I)=X(I+1)-X(I) B(I)=2.*(D(I-1)+D(I)) C(I+1)=(Y(I+1)-Y(I))/D(I) C(I)=C(I+1)-C(I) 10 CONTINUE B(1)=-D(1) B(N)=-D(N-1) C(1)=0. C(N)=0. IF (N.EQ.3) GO TO 15 C(1)=C(3)/(X(4)-X(2))-C(2)/(X(3)-X(1)) C(N)=C(N-1)/(X(N)-X(N-2))-C(N-2)/(X(N-1)-X(N-3)) C(1)=C(1)*D(1)**2/(X(4)-X(1)) C(N)=-C(N)*D(N-1)**2/(X(N)-X(N-3)) 15 DO 20 I=2,N T=D(I-1)/B(I-1) B(I)=B(I)-T*D(I-1) C(I)=C(I)-T*C(I-1) 20 CONTINUE C(N)=C(N)/B(N) DO 30 IB=1,NM1 I=N-IB C(I)=(C(I)-D(I)*C(I+1))/B(I) 30 CONTINUE B(N)=(Y(N)-Y(NM1))/D(NM1)+D(NM1)*(C(NM1)+2.*C(N)) DO 40 I=1,NM1 B(I)=(Y(I+1)-Y(I))/D(I)-D(I)*(C(I+1)+2.*C(I)) D(I)=(C(I+1)-C(I))/D(I) C(I)=3.*C(I) 40 CONTINUE C(N)=3.*C(N) D(N)=D(N-1) RETURN 50 I=1 B(I)=(Y(2)-Y(1))/(X(2)-X(1)) C(1)=0. D(1)=0. B(2)=B(1) C(2)=0. D(2)=0. RETURN END //GO.FT01F001 DD DSN=A.M66774.SW35265.NE70350.REDUCED.DAFILE, // DISP=SHR,LABEL=(,,,IN) //GO.SYSIN DD * GRID 0004 050.00000 330.00000 050.50000 330.00000 051.00000 330.00000 051.00000 331.00000 10.0000 002.98497 000.50000 0035 .0500000 -.6352274155 .1500000 -.6097598817 .2500000 -.5937992761 .3500000 -.6137125445 .4500000 -.6524700664 .5500000 -.6641680044 .6500000 -.6191783843 .7500000 -.5284651491 .8500000 -.4251898481 .9500000 -.3313138102 1.0500000 -.2463315798 1.1500000 -.1616584493 1.2500000 -.0757554642 1.3500000 .0067780342 1.4500000 .0840416263 1.5500000 .1587937783 1.6500000 .2329758779 1.7500000 .3047420898 1.8500000 .3716928187 1.9500000 .4341187881 2.0500000 .4934516996 2.1500000 .5492422862 2.2500000 .5995406127 2.3500000 .6435180810 2.4500000 .6819575638 2.5500000 .7152043921 2.6500000 .7421347971 2.7500000 .7616109284 2.8500000 .7738343657 2.9500000 .7795058968 3.0500000 .7783813695 3.1500000 .7695661770 3.2500000 .7528953127 3.3500000 .7291113778 3.4500000 .6986867250 035.00000 265.00000 070.00000 350.00000 053.00000 300.00000 054.00000 301.00000 010.00000 //