C*********************************************************************** C* THE FOLLOWING PARAMETERS ARE READ FROM THE NAMELISTS ON UNIT 5. C* C* LIMITS C* XLATN = NORTHERN LATITUDE BAND C* XLATS = SOUTHERN LATITUDE BAND C* XLONGW = WESTERNMOST LONGTITUDE BAND C* XLONGE = EASTERNMOST LONGTITUDE BAND C* XLONGW < XLONGE C* XMIN = SPACING OF A POINT GRIDDED FILE IN MINUTES C* C* DIRACC C* LGLEN = NUMBER OF LONGTITUDE BANDS IN THE DIRECT ACCESS FILE C* IRECL = NUMBER OF WORDS PER RECORD C* NREC = NUMBER OF RECORDS IN DIRECT ACCESS FILE C* XINC = SPACING IN MINUTES OF DIRECT ACCESS FILE C* TOPLAT = NORTHERNMOST LATITUDE BAND C* BOTLAT = SOUTHERNMOST LATITUDE BAND C* WLONG = WESTERNMOST LONGITUDE BAND C* ELONG = EASTERNMOST LONGITUDE BAND C* C* OUTPUT C* NWR = NUMBER OF WORDS PER RECORD C* NRB = NUMBER OF RECORDS PER BLOCK C* IPRT = EVERY IPRT RECORD IS PRINTED C* VARIABLES C* NPTS - NUMBER OF POINTS C* NWB - NUMBER OF WORDS PER BLOCK C* NWR - NUMBER OF WORDS PER RECORD C* NRB - NUMBER OF RECORDS PER BLOCK C* IOUT - BLOCK COUNTER C************************************************************************* C C GEOID HEIGHT BI-LINEAR INTERPOLATION METHOD DIMENSION PHINS(2), XLAMWE(2), IREC(2,2), XNHT(2,2) DIMENSION OUT(500),W(10) NAMELIST /LIMITS/ XLATN,XLATS,XLONGW,XLONGE,XMIN NAMELIST /DIRACC/ LGLEN,IRECL,NREC,XINC,TOPLAT,BOTLAT,WLONG,ELONG NAMELIST /OUTPUT/ NWR,NRB,IPRT DATA OUT/500*0.0/, PAD/-5401.0/ C C INPUT PROGRAM PARAMETERS C READ(5,LIMITS) READ(5,DIRACC) READ(5,OUTPUT) WRITE(6,40) WRITE(6,41)TOPLAT,BOTLAT,WLONG,ELONG,LGLEN,IRECL,NREC,XINC WRITE(6,42)XLATN,XLATS,XLONGW,XLONGE,NWR,NRB,XMIN IF (XLATN .GT. TOPLAT ) THEN XLATN = TOPLAT PRINT *,' ************* WARNING ************* ' PRINT *,' NORTHERNMOST LATITUDE VALUE OF OUTPUT FILE (XLATN) ' PRINT *,' MUST BE LESS THAN OR EQUAL TO THE TOP LATITUDE OF ' PRINT *,' THE INPUT DIRECT ACCESS FILE (TOPLAT)!!' PRINT *,' CHECK XLATN AND/OR TOPLAT FOR VALIDITY!!!' PRINT *,' PROGRAM WILL CONTINUE WITH:' PRINT *,' XLATN = TOPLAT = ',XLATN END IF IF (XLATS .LT. BOTLAT ) THEN XLATS = BOTLAT PRINT *,' ************* WARNING ************* ' PRINT *,' SOUTHERNMOST LATITUDE VALUE OF OUTPUT FILE (XLATS) ' PRINT *,' MUST BE GREATER THAN OR EQUAL TO THE BOTTOM LATITUDE' PRINT *,' THE INPUT DIRECT ACCESS FILE (BOTLAT)!!' PRINT *,' CHECK XLATS AND/OR BOTLAT FOR VALIDITY!!!' PRINT *,' PROGRAM WILL CONTINUE WITH:' PRINT *,' XLATS = BOTLAT = ',XLATS END IF IF (XLONGW .LT. WLONG ) THEN XLONGW = WLONG PRINT *,' ************* WARNING ************* ' PRINT *,' WESTERNMOST LONGITUDE VALUE OF OUTPUT FILE (XLONGW) ' PRINT *,' MUST BE GREATER THAN OR EQUAL TO THE WEST LONGITUDE' PRINT *,' THE INPUT DIRECT ACCESS FILE (WLONG)!!' PRINT *,' CHECK XLONGW AND/OR WLONG FOR VALIDITY!!!' PRINT *,' PROGRAM WILL CONTINUE WITH:' PRINT *,' XLONGW = WLONG = ',XLONGW END IF IF (XLONGE .GT. ELONG ) THEN XLONGE = ELONG PRINT *,' ************* WARNING ************* ' PRINT *,' EASTERNMOST LONGITUDE VALUE OF OUTPUT FILE (XLONGE) ' PRINT *,' MUST BE LESS THAN OR EQUAL TO THE EAST LONGITUDE OF ' PRINT *,' THE INPUT DIRECT ACCESS FILE (ELONG)!!' PRINT *,' CHECK XLONGE AND/OR ELONG FOR VALIDITY!!!' PRINT *,' PROGRAM WILL CONTINUE WITH:' PRINT *,' XLONGE = ELONG = ',XLONGE END IF C WRITE(6,43)IPRT C NPTS = 0 NWB = NWR * NRB IOUT = 1 C C OPEN(11,ACCESS='DIRECT',RECL=IRECL,RCDS=NREC) C COORDINATES OF GRID IN MINUTES LATN = XLATN LATS = XLATS LONGW = XLONGW LONGE = XLONGE MIN = XMIN TSUP1 = ZCLOK('ON') TCPU1 = ZTIME('ON') DO 10 IPHI = LATN,LATS,-MIN IF ( LONGW .GT. LONGE ) LONGW = LONGW - 21600. PHI = FLOAT( IPHI ) DO 10 ILAM = LONGW,LONGE,MIN IF ( LONGW .LT. 0 .AND. ILAM .LT.0) THEN XLAM = FLOAT( ILAM + 21600) ELSE XLAM = FLOAT( ILAM ) END IF IF ( XLAM .EQ. 21600.0) XLAM = 0.0 CALL BORDER(PHI,XLAM,XINC,PHINS,XLAMWE) IF(PHINS(1).LT.BOTLAT)PHINS(1)=BOTLAT IF(PHINS(2).GT.TOPLAT)PHINS(2)=TOPLAT IF(XLAMWE(2).GT.ELONG)XLAMWE(2)=ELONG DO 20 I=1,2 DO 20 J=1,2 CALL RECORD( PHINS(I),XLAMWE(J),TOPLAT,WLONG,LGLEN, > XINC,IREC(I,J)) IF ( IREC(I,J) .GT. NREC .OR. > IREC(I,J) .LT. 1) THEN PRINT*,' ' PRINT*,' OUT OF BOUNDS RECORD ' PRINT*,' ' WRITE( 6,*) PHI,XLAM,IREC(I,J) GO TO 10 END IF READ(11,REC=IREC(I,J)) (W(IW),IW=1,IRECL) XNHT(I,J) = W(IRECL) 20 CONTINUE C C C INTERPOLATE THE GEOID HEIGHT FOR THE GIVEN POINT. C A0 = XNHT(1,1) A1 = XNHT(1,2) - XNHT(1,1) A2 = XNHT(2,1) - XNHT(1,1) A3 = XNHT(1,1) + XNHT(2,2) - XNHT(1,2) - XNHT(2,1) X = (XLAM - XLAMWE(1)) / (XLAMWE(2) - XLAMWE(1)) Y = (PHI - PHINS(1)) / (PHINS(2) - PHINS(1)) XNHTI = A0 + A1 * X + A2 * Y + A3 * X * Y C IF(XNHTI.LT.-200.) GO TO 10 C IF(XNHTI.LT.-200.) GO TO 10 C OUT(IOUT) = PHI/60. OUT(IOUT + 1) = XLAM/60. OUT(IOUT + 2) = XNHTI IF ( MOD(NPTS,IPRT) .EQ. 0 ) THEN WRITE(6,*) (OUT(IOUT+JJ-1), JJ = 1, NWR ) END IF IOUT = IOUT + NWR IF ( IOUT .GT. NWB) THEN WRITE(12) ( OUT(II) , II = 1, NWB ) IOUT = 1 END IF NPTS = NPTS + 1 10 CONTINUE C DO 30 I = IOUT, NWB OUT(I) = PAD 30 CONTINUE WRITE(12) ( OUT(II) , II = 1, NWB ) TSUP2 = ZCLOK('STOP') TCPU2 = ZTIME('STOP') PRINT *,' THE NUMBER OF POINTS COMPUTED ', NPTS PRINT *,' THE TOTAL SUP TIME IN SECONDS ', TSUP2 PRINT *,' THE TOTAL CPU TIME IN SECONDS ', TCPU2 PRINT *,' THE SUP TIME PER MEAN IN SECONDS',TSUP2/NPTS PRINT *,' THE CPU TIME PER MEAN IN SECONDS',TCPU2/NPTS 40 FORMAT('1','GRIDDED GEOID HEIGHTS COMPUTED FROM THE 30 X 30 MIN. ' *,'FILE',/,' BY BILINEAR INTERPOLATION',/,' (COORDINATES IN MIN.)') 41 FORMAT(///,5X,'INPUT DIRECT ACCESS FILE PARAMETERS:',//, *' TOP LATITUDE = ',F10.3,/, *' BOTTOM LATITUDE = ',F10.3,/, *' WEST LONGITUDE = ',F10.3,/, *' EAST LONGITUDE = ',F10.3,/, *' # OF LONG. BANDS = ',I5,/, *' # OF WORDS/RECORD = ',I5,/, *' # OF RECORDS = ',I5,/, *' GRID SPACING = ',F6.1) 42 FORMAT(///,5X,'OUTPUT FILE PARAMETERS',//, *' TOP LATITUDE = ',F10.3,/, *' BOTTOM LATITUDE = ',F10.3,/, *' WEST LONGITUDE = ',F10.3,/, *' EAST LONGITUDE = ',F10.3,/, *' # WORDS/RECORD = ',I5,/, *' # RECORDS/BLOCK = ',I5,/, *' GRID SPACING = ',F6.1,/) 43 FORMAT('1',4X,'EVERY NTH RECORD OF THE OUTPUT FILE IS PRINTED',/, *5X,'N = ',I5,/,5X,'COORDINATES IN DECIMAL DEGREES',///) 99 STOP END SUBROUTINE RECORD( LAT,LONG,TOPLAT,WLONG,LGLEN,XINC,IREC) IMPLICIT REAL*4 (A-H,O-Z) REAL*4 LAT, LONG IPOSLT = NINT((TOPLAT - LAT)/XINC) IPOSLG = NINT((LONG-WLONG)/XINC) + 1 IREC = IPOSLT * LGLEN + IPOSLG RETURN END SUBROUTINE BORDER(PHI,XLAM,XINC,PHINS,XLAMWE) IMPLICIT REAL*4 (A-H, O-Z) DIMENSION PHINS(2), XLAMWE(2) C LOCATION OF GRID POINTS PHIINT = INT(PHI/XINC) * XINC XLMINT = INT(XLAM/XINC) * XINC IF ( PHI .LT. 0.0 ) THEN XLIM1 = PHIINT + SIGN(XINC,PHI) XLIM2 = PHIINT ELSE XLIM1 = PHIINT XLIM2 = PHIINT + SIGN(XINC,PHI) END IF PHINS(1) = XLIM1 PHINS(2) = XLIM2 XLAMWE(1) = XLMINT XLAMWE(2) = XLAMWE(1) + XINC IF (XLAMWE(1) .EQ. 21600.0) XLAMWE(1) = 0.0 IF (XLAMWE(2) .EQ. 21600.0) XLAMWE(2) = 0.0 RETURN END END OF FILE