//OPT3OSU JOB NOTIFY=6461 /*SETUP SLOT=P0147 VOLUME=NL6250 NOWRITE /*SERVICE DEFERRED /*JOBPARM T=80,L=999,R=4096,PRINT=ALL //STEP1 EXEC FORTVCLG,REGION=4096K,PARM.FORT='OPTIMIZE(3)' //FORT.SYSIN DD * C----------------------------------------------------------------------- C*********************************************************************** C * C PROGRAM NAME : GRIDALT * C FUNCTION : GENERATION OF A REGULAR GRID OF ADJUSTED * C SEA SURFACE HEIGHTS FROM SEASAT ADJUSTED * C ALTIMETRY DATA. * C COMPILER : VS FORTRAN VERSION 2 * C AUTHOR : NICK CHRISTOU * C HISTORY : JULY 13, 1985 - VERSION 1.0 * C AUGUST 20, 1985 - VERSION 2.0 (OPTIMIZED) * C REFERENCE : U.N.B. TECHNICAL REPORT # * C * C EXTERNALS : INTPLT , INPUT , LODATA , SPIN . * C*********************************************************************** C----------------------------------------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION GPHI(300),FVALUE(600),SDEVF(600) DIMENSION XLAT(40000),XLON(40000),XASSH(40000),XSD(40000) C COMMON /STRIP/XLAT,XLON,XASSH,XSD COMMON GLON(600),PHIMIN,PHIMAX,DLAMIN,DLAMAX,HWMAX1,HWMAX2, @ HWINT1,HWINT2,HWMIN1,HWMIN2,IPOINT,NOBS,JJ C READ(5,1000) PHIMIN,DLAMIN,PHIMAX,DLAMAX,INCPHI,INCLAM 1000 FORMAT(4F9.4,2I3) C----------------------------------------------------------------------- C PHIMIN : C PHIMAX : REPRESENT THE BOUNDARIES OF THE AREA FOR WHICH C DLAMIN : THE REGULAR GRID WILL BE GENERATED (IN DEC. DEG.) C DLAMAX : C C INCPHI : GRID INTERVAL IN LATITUDE (IN MINUTES) C INCLAM : GRID INTERVAL IN LONGITUDE (IN MINUTES) C----------------------------------------------------------------------- READ(5,1001) HWMAX1,HWMAX2,HWINT1,HWINT2,HWMIN1,HWMIN2 1001 FORMAT(6F5.2) C----------------------------------------------------------------------- C HWMAX1 : MAXIMUM C HWINT1 : ARE THE INTERMEDIATE HALF-WIDTHS OF THE WORKING WINDOWS C HWMIN1 : MINIMUM IN THE LATITUDE DIRECTION C C HWMAX2 : MAXIMUM C HWINT2 : ARE THE INTERMEDIATE HALF-WIDTHS OF THE WORKING WINDOWS C HWMIN2 : MINIMUM IN THE LONGITUDE DIRECTION C----------------------------------------------------------------------- C \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ C GENERATE NOW THE DESIRED GRID . C ////////////////////////////////////////////////////////////// C----------------------------------------------------------------------- DELPHI = PHIMAX - PHIMIN DELLAM = DLAMAX - DLAMIN C DPHI = INCPHI/60.0D0 DLAMDA = INCLAM/60.0D0 C II = DELPHI/DPHI + 1.0D0 JJ = DELLAM/DLAMDA + 1.0D0 C----------------------------------------------------------------------- C II : NUMBER OF GRID POINTS IN LATITUDE (DEFINES DIM. OF GPHI) C JJ : -"- -"- -"- LONGITUDE (DEFINES DIM. OF GLON) C----------------------------------------------------------------------- GPHI(1) = PHIMIN GLON(1) = DLAMIN C DO 11 I=1,(II-1) GPHI(I+1)=GPHI(1)+DPHI*I 11 CONTINUE DO 22 J=1,(JJ-1) GLON(J+1)=GLON(1)+DLAMDA*J 22 CONTINUE C----------------------------------------------------------------------- C THE GENERATED GRID IS CONTAINED IN THE ARRAYS GPHI(I) AND GLON(I) C======================================================================= C NEXT STEP C======================================================================= C C GENERATION OF THE FUNCTIONAL VALUES ON THE REGULAR GRID C ( HERE THE FUNCTION IS THE SEA SURFACE HEIGHT ) C----------------------------------------------------------------------- DO 88 K=1,II GLAT=GPHI(K) CALL INTPLT(GLAT,FVALUE,SDEVF) WRITE(6,9000) (GLAT,GLON(L),FVALUE(L),SDEVF(L),L=1,JJ) 88 CONTINUE 9000 FORMAT(F8.4,F10.4,F8.2,F6.2) STOP END SUBROUTINE INTPLT(GLAT,FVALUE,SDEVF) C----------------------------------------------------------------------- C*********************************************************************** IMPLICIT REAL*8 (A-H,O-Z) INTEGER NP(4) C DIMENSION XLAT(40000),XLON(40000),XASSH(40000),XSD(40000) DIMENSION YLAT(6000),YLON(6000),YVAL(6000),YDEV(6000) DIMENSION ZLAT(6000),ZLON(6000),ZVAL(6000),ZDEV(6000) C DIMENSION X(6000),Y(6000),Z(6000),W(6000) DIMENSION FF(6000,4),DD(10,10),CC(4),XX(4) C DIMENSION FVALUE(600),SDEVF(600) C COMMON /STRIP/ XLAT,XLON,XASSH,XSD COMMON GLON(600),PHIMIN,PHIMAX,DLAMIN,DLAMAX,HWMAX1,HWMAX2, @ HWINT1,HWINT2,HWMIN1,HWMIN2,IPOINT,NOBS,JJ C C PI = 3.141592653589793D0 C----------------------------------------------------------------------- C HWMAX1 : MAXIMUM C HWINT1 : ARE THE INTERMEDIATE HALF-WIDTHS OF THE WORKING WINDOWS C HWMIN1 : MINIMUM IN THE LATITUDE DIRECTION C C HWMAX2 : MAXIMUM C HWINT2 : ARE THE INTERMEDIATE HALF-WIDTHS OF THE WORKING WINDOWS C HWMIN2 : MINIMUM IN THE LONGITUDE DIRECTION C C======================================================================= C C FORM THE WORKING STRIP OF LATITUDE BASED ON THE LATITUDE OF C THE COMPUTATION POINT : I.E. , GLAT . C C DMIN : LOWER BOUNDARY OF THE LATITUDE STRIP C DMAX : UPPER -"- -"- -"- -"- . C----------------------------------------------------------------------- DMIN = GLAT - HWMAX1 DMAX = GLAT + HWMAX1 C----------------------------------------------------------------------- C READ ALL THE ALTIMETRY DATA POINTS WHICH FALL INSIDE THE C ABOVE SPECIFIED LATITUDE STRIP AND STORE THEM IN ARRAYS : C C XLAT(I) : LATITUDE C XLON(I) : LONGITUDE C : CONTAINS THE, ,OF EACH DATA POINT C XASSH(I): ADJ. SEA SURF. HEIGHT C XSD(I) : ST. DEV. OF ASSH C C THE TOTAL NUMBER OF POINTS FALLEN INSIDE THE STRIP IS: *** IPOINT *** C----------------------------------------------------------------------- CALL INPUT(GLAT,DMIN,DMAX) C----------------------------------------------------------------------- C FORM NOW THE MAXIMUM WORKING WINDOW AROUND EACH COMPUTATION POINT C BASED ON THE LONGITUDE OF THE GRID NODE : I.E. , GLON(KK) . C C DLMIN : THE LEFT BOUNDARY OF THE MAXIMUM WORKING WINDOW (KK). C DLMAX : THE RIGHT -"- -"- -"- -"- -"- (KK) . C------------------------------------------------------------------------ DO 66 KK=1,JJ DLMIN = GLON(KK) - HWMAX2 DLMAX = GLON(KK) + HWMAX2 C----------------------------------------------------------------------- C LOAD NOW THOSE DATA POINTS FROM THE LATITUDE STRIP WHICH BELONG C TO THE MAXIMUM WINDOW (KK) , AND PERFORM THE PRELIMINARY TEST C ON THE NUMBER OF POINTS , NO , FALLEN INSIDE THIS WINDOW : C ( I.E., NO >= 8 ) C----------------------------------------------------------------------- ICOUNT=0 DO 77 I=1,IPOINT IF(XLON(I).GT.DLMAX.OR.XLON(I).LT.DLMIN) GO TO 77 ICOUNT = ICOUNT+1 YLAT(ICOUNT) = XLAT(I) YLON(ICOUNT) = XLON(I) YVAL(ICOUNT) = XASSH(I) YDEV(ICOUNT) = XSD(I) 77 CONTINUE NO = ICOUNT IF(NO.LT.8) GO TO 7777 C----------------------------------------------------------------------- C NOTE 1 : IF THE NUMBER OF POINTS , NO , FOUND IN THE MAXIMUM WINDOW C IS LESS THAN 8 , THEN THE PROGRAMME IS DIRECTED TO C STATEMENT 7777 WHICH ASSIGNS A DEFAULT VALUE OF HEIGHT C AND STANDARD DEVIATION TO THIS GRID NODE . C C NOTE 2 : THE DATA POINTS BELONGING TO THE MAXIMUM WINDOW ARE C CONTAINED IN THE ARRAYS : C YLAT(I) , YLON(I) , YVAL(I) , YDEV(I) . C----------------------------------------------------------------------- C C======================================================================= C GENERATE NOW THE MINIMUM WORKING WINDOW TO BE USED FOR THE C INTERPOLATION OF THE HEIGHT FUNCTION ON THE (KK) GRID NODE. C C FMIN : THE LOWER (LATITUDE) BOUNDARY OF THE MINIMUM WINDOW C FMAX : THE UPPER -"- -"- -"- -"- -"- C FLMIN : THE LEFT (LONGITUDE) BOUNDARY OF THE MINIMUM WINDOW C FLMAX : THE RIGHT -"- -"- -"- -"- -"- C======================================================================= FMIN = GLAT - HWMIN1 FMAX = GLAT + HWMIN1 FLMIN = GLON(KK) - HWMIN2 FLMAX = GLON(KK) + HWMIN2 C LCOUNT=0 DO 111 I=1,NO IF(YLAT(I).GT.FMAX.OR.YLAT(I).LT.FMIN) GO TO 111 IF(YLON(I).GT.FLMAX.OR.YLON(I).LT.FLMIN) GO TO 111 LCOUNT = LCOUNT+1 ZLAT(LCOUNT)=YLAT(I) ZLON(LCOUNT)=YLON(I) ZVAL(LCOUNT)=YVAL(I) ZDEV(LCOUNT)=YDEV(I) 111 CONTINUE NOP=LCOUNT IF(NOP.LT.8) GO TO 3333 C----------------------------------------------------------------------- C NOTE 3 : IF THE NUMBER OF POINTS , NOP , IN THE MINIMUM WINDOW C IS LESS THAN 8 , THEN THE PROGRAMME IS DIRECTED TO C STATEMENT 3333 WHICH ASSIGNS THE NEW WINDOW BOUNDARIES C DEFINED BY HWINT1 AND HWINT2 , AND REPEATS THE SAME TEST. C C NOTE 4 : THE DATA POINTS BELONGING IN THE MINIMUM WINDOW ARE C CONTAINED IN THE ARRAYS : C ZLAT(I) , ZLON(I) , ZVAL(I) , ZDEV(I) . C C NOTE 5 : THE NEXT TEST PERFORMED IS RELATED TO DATA-POINT C DISTRIBUTION INSIDE THE MINIMUM WINDOW . C C NP(1) , NP(2) , NP(3) , NP(4) : ------------------------ C THEY ARE ARRAY ELEMENTS | | | C | NP(2) | NP(1) | C REPRESENTING THE NUMBER OF POINTS | | | C ------------------------ C IN EACH QUQDRANT OF THE MINIMUM | | | C | NP(4) | NP(3) | C WINDOW IN THE FOLLOWING MANNER ; | | | C ------------------------ C THE TOTAL NUMBER OF DATA POINTS IN THIS WINDOW IS : **** NOP **** C------------------------------------------------------------------------ NP(1)=0 NP(2)=0 NP(3)=0 NP(4)=0 DO 222 I=1,NOP N=0 IF(ZLAT(I).LT.GLAT) N=2 IF(ZLON(I).LT.GLON(KK)) N=N+2 IF(ZLON(I).GT.GLON(KK)) N=N+1 NP(N)=NP(N)+1 222 CONTINUE IF(NP(1).GE.1.AND.NP(2).GE.1.AND. & NP(3).GE.1.AND.NP(4).GE.1) GO TO 999 C----------------------------------------------------------------------- C NOTE 6 : IF THE NUMBER OF DATA POINTS INSIDE THE MINIMUM WINDOW C IS GREATER THAN 8 AND IN ADDITION THERE IS AT LEAST 1 C DATA POINT IN EACH QUADRANT OF THE WINDOW THEN THE PROGRAMME C IS DIRECTED TO STATEMENT 999 WHICH IS THE FIRST C STATEMENT OF THE INTERPOLATION ALGORITHM . C C NOTE 7 : IF EITHER OF THE ABOVE CONDITIONS IS NOT SATISFIED THEN C THE PROGRAMME IS AGAIN DIRECTED TO STATEMENT 3333 WHICH C WILL ASSIGN THE NEXT WINDOW BOUNDARIES ( IN THIS CASE C INTERMEDIATE SIZE WINDOW ) , AND REPEAT ONCE MORE THE SAME C OPERATIONS AS BEFORE . C----------------------------------------------------------------------- 3333 FMIN = GLAT - HWINT1 FMAX = GLAT + HWINT1 FLMIN = GLON(KK) - HWINT2 FLMAX = GLON(KK) + HWINT2 C JCOUNT=0 DO 333 I=1,NO IF(YLAT(I).GT.FMAX.OR.YLAT(I).LT.FMIN) GO TO 333 IF(YLON(I).GT.FLMAX.OR.YLON(I).LT.FLMIN) GO TO 333 JCOUNT=JCOUNT+1 ZLAT(JCOUNT)=YLAT(I) ZLON(JCOUNT)=YLON(I) ZVAL(JCOUNT)=YVAL(I) ZDEV(JCOUNT)=YDEV(I) 333 CONTINUE NOP=JCOUNT IF(NOP.LT.8) GO TO 4444 C----------------------------------------------------------------------- C NOTE 8 : IF THE NUMBER OF POINTS , NOP , IN THE INTERMEDIATE C WINDOW IS LESS THAN 8 , THEN THE PROGRAMME IS DIRECTED C TO STATEMENT 4444 WHICH IS GOING TO CHECK ONLY ABOUT C THE DISTRIBUTION OF DATA POINTS IN THE MAXIMUM WINDOW, C SINCE THE PRELIMINARY TEST ON THE TOTAL NUMBER OF POINTS C IN THE MAX. WINDOW HAS ALREADY BEEN DONE . C C NOTE 9 : THE ARRAYS : ZLAT(I) , ZLON(I) , ZVAL(I) , ZDEV(I) C CONTAIN NOW THE DATA POINTS FOUND IN THE INTERMEDIATE C WORKING WINDOW . C C NOTE 10 : THE NEXT TEST PERFORMED IS RELATED TO THE DATA-POINT C DISTRIBUTION WITHIN THE INTERMEDIATE WINDOW . C NP(1) , NP(2) , NP(3) , NP(4) REPRESENT NOW THE NUMBER OF C POINTS FOUND IN EACH QUADRANT OF THE INTERMEDIATE WINDOW C AND THEY HAVE THE SAME DESIGNATION AS IN THE FIGURE OF C NOTE # 5 . C C THE TOTAL NUMBER OF POINTS IS DESIGNATED BY : **** NOP **** C----------------------------------------------------------------------- NP(1)=0 NP(2)=0 NP(3)=0 NP(4)=0 DO 444 I=1,NOP N=0 IF(ZLAT(I).LT.GLAT) N=2 IF(ZLON(I).LT.GLON(KK)) N=N+2 IF(ZLON(I).GT.GLON(KK)) N=N+1 NP(N)=NP(N)+1 444 CONTINUE IF(NP(1).GE.1.AND.NP(2).GE.1.AND. & NP(3).GE.1.AND.NP(4).GE.1) GO TO 999 C----------------------------------------------------------------------- C NOTE 11 : IF THE NUMBER OF DATA POINTS IN THE INTERMEDIATE WINDOW C IS LARGER THAN 8 AND IN ADDITION THERE IS AT LEAST 1 C DATA POINT IN EACH QUADRANT THEN THE PROGRAMME IS C DIRECTED TO STATEMENT 999 WHICH IS THE FIRST STATEMENT C OF THE INTERPOLATION ALGORITHM . C C NOTE 12 : IF BOTH THE ABOVE CONDITIONS ARE NOT SATISFIED ,THEN THE C PROGRAMME WILL BE DIRECTED TO STATEMENT 4444 WHICH WILL C TEST THE DATA DISTRIBUTION INSIDE THE MAXIMUM WINDOW . C C NOTE 13 : AT THIS STAGE THE ELEMENTS OF THE Y-DESIGNATED ARRAYS ARE C TRANSFERRED TO THE Z-DESIGNATED ARRAYS . C----------------------------------------------------------------------- 4444 KCOUNT=0 DO 555 I=1,NO KCOUNT=KCOUNT+1 ZLAT(KCOUNT)=YLAT(I) ZLON(KCOUNT)=YLON(I) ZVAL(KCOUNT)=YVAL(I) ZDEV(KCOUNT)=YDEV(I) 555 CONTINUE NOP=KCOUNT C NP(1)=0 NP(2)=0 NP(3)=0 NP(4)=0 DO 666 I=1,NOP N=0 IF(ZLAT(I).LT.GLAT) N=2 IF(ZLON(I).LT.GLON(KK)) N=N+2 IF(ZLON(I).GT.GLON(KK)) N=N+1 NP(N)=NP(N)+1 666 CONTINUE IF(NP(1).GE.1.AND.NP(2).GE.1.AND. & NP(3).GE.1.AND.NP(4).GE.1) GO TO 999 C GO TO 7777 C----------------------------------------------------------------------- C NOTE 14 : IF THE DISTRIBUTION OF DATA POINTS INSIDE THE MAXIMUM C WINDOW IS SATISFACTORY (I.E., ONE DATA POINT IN EACH C QUADRANT) THEN THE PROGRAMME CARRIES ON WITH THE C INTERPOLATION OF THE HEIGHT FUNCTION ON THE GRID NODE . C C NOTE 15 : IF THE DISTRIBUTION OF DATA POINTS IN THE MAXIMUM C WINDOW FAILS THE TEST THEN THE PROGRAMME QUITS TRYING C AND IS DIRECTED TO STATEMENT 7777 WHICH ASSIGNS THE C CHOSEN DEFAULT VALUES FOR HEIGHT AND STANDARD DEVIATION C TO THAT GRID POINT . C----------------------------------------------------------------------- 999 DO 998 I=1,NOP X(I)=69.041*(ZLAT(I)-GLAT) YY=69.041*DCOS((GLAT*PI/180.0D0)) Y(I)=YY*(ZLON(I)-GLON(KK)) Z(I)=ZVAL(I) W(I)=1.0D0/(ZDEV(I)*ZDEV(I)) 998 CONTINUE DO 5 I=1,NOP FF(I,1)=1.0D0 FF(I,2)=X(I) FF(I,3)=Y(I) 5 FF(I,4)=Y(I)*X(I) C 5 CONTINUE DO 6 I=1,4 DO 6 J=1,4 DD(I,J)=0.0D0 DO 6 M=1,NOP 6 DD(I,J)=DD(I,J)+FF(M,I)*FF(M,J)*W(M) C 6 CONTINUE CALL SPIN(DD,4,10,DET,IDEXP) QQ=DD(1,1) DO 7 I=1,4 CC(I)=0.0D0 DO 7 J=1,NOP 7 CC(I)=CC(I)+FF(J,I)*W(J)*Z(J) DO 8 I=1,4 XX(I)=0.0D0 DO 8 J=1,4 8 XX(I)=XX(I)+DD(I,J)*CC(J) C FVALUE(KK)=XX(1) C PVV=0.0D0 DO 9 I=1,NOP POLYNO=0.0D0 DO 10 J=1,4 10 POLYNO=POLYNO+FF(I,J)*XX(J) 9 PVV=PVV+W(I)*(Z(I)-POLYNO)*(Z(I)-POLYNO) SIGMA0=DSQRT(PVV/DFLOAT(NOP-4)) C SDEVF(KK)=SIGMA0*DSQRT(QQ) C GO TO 66 C 7777 FVALUE(KK)=0.0D0 SDEVF(KK)=50.0D0 66 CONTINUE RETURN END SUBROUTINE INPUT(GLAT,DMIN,DMAX) C----------------------------------------------------------------------- C*********************************************************************** IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION XLAT(40000),XLON(40000),XASSH(40000),XSD(40000) C COMMON /STRIP/ XLAT,XLON,XASSH,XSD COMMON GLON(600),PHIMIN,PHIMAX,DLAMIN,DLAMAX,HWMAX1,HWMAX2, @ HWINT1,HWINT2,HWMIN1,HWMIN2,IPOINT,NOBS,JJ C----------------------------------------------------------------------- C NOTE 1 : IF THE INPUT GRID-LATITUDE (GLAT) IS THE FIRST LATITUDE C LINE ALONG WHICH THE INTERPOLATION OF THE ( HEIGHT ) C FUNCTION H(PHI,LAMDA) WILL BE PERFORMED , THEN THE C PROGRAMME IS DIRECTED TO STATEMENT 120 WHICH LOADS THE C DATA POINTS (FROM A TAPE FILE) THAT FALL INSIDE THE C FIRST FORMED LATITUDE STRIP. C THE BOUNDARIES OF THE STRIP ARE INDICATED BY THE VARIABLES: C DMIN , DMAX . C C NOTE 2 : THE DATA POINTS LOADED FROM THE TAPE FILE ARE STORED IN C THE ARRAYS : C XLAT(I) , XLON(I) , XASSH(I) , XSD(I) ; C THEY CONTAIN **** IPOINT ***** NUMBER OF ELEMENTS . C C NOTE 3 : WHEN THE PROGRAMME HAS COMPLETED THE PROCESSING OF THE C FIRST LATITUDE-GRID-LINE , THE CONTROL OF LOADING DATA C WILL BE TRANSFERRED TO STATEMENT 121 . C AT THIS STAGE, A NEW LATITUDE STRIP WILL HAVE BEEN FORMED . C THE OLD DATA POINTS THAT ARE COMMON IN THE NEW LATITUDE C STRIP WILL BE RETAINED , HOWEVER THEY WILL RESUME NEW C POSITIONS AS ELEMENTS OF THE X-DESIGNATED ARRAYS . C----------------------------------------------------------------------- IF(GLAT.EQ.PHIMIN) GO TO 120 GO TO 121 120 NOBS=0 GO TO 122 121 INUM=0 DO 123 I=1,IPOINT IF(XLAT(I).LT.DMIN) GO TO 123 INUM=INUM+1 XLAT(INUM)=XLAT(I) XLON(INUM)=XLON(I) XASSH(INUM)=XASSH(I) XSD(INUM)=XSD(I) 123 CONTINUE NOBS = INUM C----------------------------------------------------------------------- C NOTE 4 : THE DATA POINTS COMMON TO THE OLD AND NEW LATITUDE STRIPS C HAVE NOW BEEN RELOCATED IN THE X-DESIGNATED ARRAYS . C THEY CONTAIN NOW **** NOBS **** NUMBER OF ELEMENTS: C NOBS < OR = IPOINT(OLD) C C NOTE 5 : THE NEXT STEP IS TO UPDATE THE X-ARRAYS WITH ADDITIONAL C DATA POINTS (IF ANY) THAT MAY LIE BETWEEN THE UPPER C BOUNDARIES OF THE OLD AND NEW LATITUDE STRIPS. C THE TOTAL NUMBER OF OBSERVATIONS WILL THEN BECOME : C IPOINT(NEW) > OR = NOBS C----------------------------------------------------------------------- 122 CALL LODATA(DMIN,DMAX) RETURN END SUBROUTINE LODATA(DMIN,DMAX) C----------------------------------------------------------------------- C*********************************************************************** IMPLICIT REAL*8 (A-H,O-Z) REAL*4 XXLAT,XXLON INTEGER*4 REVNO,MJD CHARACTER*7 XXSSH C DIMENSION XLAT(40000),XLON(40000),XASSH(40000),XSD(40000) C COMMON /STRIP/ XLAT,XLON,XASSH,XSD COMMON GLON(600),PHIMIN,PHIMAX,DLAMIN,DLAMAX,HWMAX1,HWMAX2, @ HWINT1,HWINT2,HWMIN1,HWMIN2,IPOINT,NOBS,JJ C----------------------------------------------------------------------- C NOTE 1 : WHENEVER THE AREA FOR WHICH THE REGULAR GRID OF THE HEIGHT C FUNCTION WILL BE CONSTRUCTED IS SMALLER THAN THE GENERAL C AREA OF DATA COVERAGE, THEN IT IS DESIRABLE TO CONFINE OUR C INTEREST TO A SUB-AREA CONTAINING ONLY THE RELEVANT DATA . C ONLY THE LONGITUDINAL BOUNDARIES ARE NEEDED TO BE SET , C SINCE THE LATITUDINAL ONES ARE CONTROLLED THROUGH THE C VARIABLES DMIN AND DMAX . C C THE OSU SEASAT ADJUSTED ALTIMETRY DATA SET COVERS THE AREA : C 35.0 N =< LATITUDE =< 72.1667 N C 260.0 E =< LONGITUDE =< 350.1667 E C C NOTE 2 : DLAM1 : REPRESENT THE LONGITUDINAL BOUNDARIES C DLAM2 : OF THE SUB-AREA OF RELEVANT DATA POINTS. C C IF THIS IS THE CASE THEN REMOVE THE COMMENT CARD FROM THE !!!!! C STATEMENTS BELOW AND THE COMMENT CARD FROM THE IF STATEMENT !!!!! C C NOTE 3 : THE DATA POINTS ARE READ FROM A TAPE FILE . C THE DATA ON THE TAPE ARE SORTED IN INCREASING LONGITUDE C WITHIN INCREASING LATITUDE . C C NOTE 4 : DATA POINTS WITH ZERO STAND. DEV. FOR HEIGHT ARE REJECTED. C C NOTE 5 DATA POINTS WITH ASSH ABSOLUTE VALUE >=70.0M ARE REJECTED. C C NOTE 6 : THE RELEVANT DATA POINTS ARE STORED IN THE ARRAYS : C XLAT(I) , XLON(I) , XASSH(I) , XSD(I) . C----------------------------------------------------------------------- C DLAM1 = DLAMIN - HWMAX2 C DLAM2 = DLAMAX + HWMAX2 NCOUNT=NOBS 2221 READ(20,1002,END=2222) XXLAT,XXLON,SD,IT,GEM9,ASSH,SN,XTIDE,XSSH 1002 FORMAT(2F10.5,F5.2,I10,2F7.2,I10,2F8.2) C IF(XXLAT.LT.DMIN.OR.SD.EQ.0.0E0) GO TO 2221 IF(DABS(ASSH).GT.70.0D0) GO TO 2221 C C IF(XXLON.LT.DLAM1.OR.XXLON.GT.DLAM2) GO TO 2221 C IF(XXLAT.GT.DMAX) GO TO 2222 C NCOUNT=NCOUNT+1 XLAT(NCOUNT)=XXLAT XLON(NCOUNT)=XXLON XASSH(NCOUNT)=ASSH XSD(NCOUNT)=SD GO TO 2221 2222 CONTINUE IPOINT=NCOUNT RETURN END C********************************************************************** SUBROUTINE SPIN(Q,N,MM,DET,IDEXP) C----------------------------------------------------------------------- C FUNCTION : C SUBROUTINE SPIN IS A MATRIX INVERSION ROUTINE FOR SYMMETRIC C POSITIVE-DEFINITE MATRICES. THE MATRIX INVERTED IS THE UPPER C N BY N PORTION OF THE MATRIX Q WHICH IS DIMENSIONED MM BY MM C IN THE CALLING ROUTINE. C AUTHOR : R.R.STEEVES C HISTORY : SEPTEMBER 1979 C----------------------------------------------------------------------- C INPUT: C Q - THE MATRIX DIMENSIONED MM BY MM WHICH CONTAINS THE C MATRIX TO BE INVERTED. C C N - THE DIMENSION OF THE ACTUAL PART( UPPER LEFT CORNER) C OF Q WHICH IS TO BE INVERTED. ( N MAY BE EQUAL BUT MUST C NOT BE LARGER THAN MM) . C MM- DIMENSIONED SIZE OF Q IN THE CALLING ROUTINE. C C C OUTPUT: C C Q - THE UPPER LEFT N BY N PORTION CONTAINS THE INVERSE OF C THE INPUT UPPER LEFT N BY N PORTION. C C DET - THE NON-EXPONENT PORTION OF THE DETERMINANT OF THE C INPUT N BY N (UPPER LEFT PORTION OF Q) MATRIX. SEE C IDEXP BELOW. C C IDEXP - THE EXPONENT (OF 10) PART OF THE DETERMINANT DESCRIBED C UNDER DET ABOVE. THUS THE DETERMINANT IS RETURNED IN C TWO PARTS CORRESPONDING TO C DETERMINANT = DET * 10 ** IDEXP . C THIS IS DONE TO AVOID UNDER OR OVERFLOW IN THE C COMPUTATION OF THE DETERMINANT. TO PRINT THE DETERM- C INANT THE USER SHOULD PRINT BOTH NUMBERS AS FOLLOWS; C (FOR EXAMPLE) C PRINT 10,DET,IDEXP C 10 FORMAT(' ','DETERMINANT= ',F7.4,'D',I4) C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C REAL*8 Q(MM,MM),DET,SUM,DSQRT,DABS,RPART,APART DET=0.D0 DO 4 J=1,N DO 4 I=1,J IF(I.EQ.1) GO TO 2 20 M=I-1 SUM=0.0D0 DO 1 K=1,M 1 SUM=SUM+Q(K,I)*Q(K,J) Q(I,J)=Q(I,J)-SUM 2 IF(I.EQ.J) GO TO 3 Q(I,J)=Q(I,J)/Q(I,I) GO TO 4 3 CONTINUE DET=DET+DLOG10(Q(I,I)) Q(I,I)=DSQRT(Q(I,I)) 4 CONTINUE IDEXP=DET RPART=DET-IDEXP APART=DABS(RPART) IF(APART.LT.1.D-20) DET=1.D0 IF(APART.LT.1.D-20) GO TO 10 DET=10**RPART 10 CONTINUE DO 7 J=1,N DO 7 I=1,J IF(I.LT.J) GO TO 5 Q(J,J)=1.0D0/Q(J,J) GO TO 7 5 SUM=0.0D0 M=J-1 DO 6 K=I,M 6 SUM=SUM-Q(I,K)*Q(K,J) Q(I,J)=SUM/Q(J,J) 7 CONTINUE DO 9 J=1,N DO 9 I=1,J SUM=0.0D0 DO 8 K=J,N 8 SUM=SUM+Q(I,K)*Q(J,K) Q(I,J)=SUM IF(I.EQ.J) GO TO 9 Q(J,I)=SUM 9 CONTINUE RETURN END //GO.SYSIN DD * 35.0000 260.0000 72.0000 350.0000 10 10 2.00 2.00 1.00 1.00 0.50 0.50 //GO.FT20F001 DD UNIT=TAPE6250,VOL=SER=P0147, // LABEL=(1,NL),DISP=OLD, // DCB=(RECFM=FB,LRECL=75,BLKSIZE=32700) //