C ********************************************************* C * * C * PROGRAM TOPO * C * """"""""""""""""" * C * VERSION 1.0 * C * """""""""""" * C * * C * COMPUTES TOPOGRAPHIC CORRECTIONS * C * """""""""""""""""""""""""""""""""" * C * * C * UNIVERSITY OF NEW BRUNSWICK * C * SURVEY ENGINEERING DEPT. 1986 * C * * C * THIS PROGRAM COMPUTES THE TOPOGRAPHIC * C * CORRECTIONS TO THE 5X5 OR 5X10 MEAN GRAVITY ANOMALIES * C * USED TO COMPUTE THE UNB GRAVIMETRIC GEOID. * C * * C * DATE:JULY 17 1986 * C * WRITTEN BY: CHRISTER STORHED * C * JOHN MANTHA * C * * C ********************************************************* C------------------------------------------------------------------------ C INITIALIZATION OF ARRAYS C """""""""""""""""""""""""" IMPLICIT REAL*4(A-H,O-Z) DIMENSION PHI(10000),ALA(10000) DIMENSION DGI(9216),SDGI(9216),ELVANI(9216) DIMENSION PHINZ(576),ALANZ(576),DGINZ(576),SDGINZ(576),HEINZ(576) C C------------------------------------------------------------------------ C________________________________________________________________________ C C CARD 1: C """"""" C USED IN SUBROUTINE AREA C ----------------------- C READS IN THE VALUES OF THE DESIRED AREA TO BE GRIDDED ALONG C WITH THE INTERVAL TO BE USED TO GENERATE THE DESIRED GRID. C C PHISW:LATITUDE OF THE SOUTH WEST CORNER OF GRID TO BE COMPUTED C IN DEGREES. C ALASW:LONGITUDE OF THE SOUTH WEST CORNER OF GRID TO BE COMPUTED C IN DEGREES. C PHINW:LATITUDE OF THE NORTH WEST CORNER OF GRID TO BE COMPUTED C IN DEGREES. C ALASE:LONGITUDE OF THE SOUTH EAST CORNER OF GRID TO BE COMPUTED C IN DEGREES. C XNTPHI:INTERVAL ALONG THE LATITUDINAL DIRECTION IN MINUTES C XNTALA:INTERVAL ALONG THE LONGITUDINAL DIRECTION IN MINUTES C READ(5,*)PHISW,ALASW,PHINW,ALASE,XNTPHI,XNTALA C________________________________________________________________________ C------------------------------------------------------------------------ C INITIALIZATION OF CONSTANTS ITOL1=0 FF=0.0 C C_____________________________________________________________________ C GRID MODE C """"""""" C C THE PROGRAM WILL USE THE GRID MODE AND CALCUALTE THE COORDINATES C OF THE GRID INTERSECTION LINES THROUGH THE CALL TO THE SUBROUTINE C AREA. C C THE AREA: C SUBROUTINE GENERATES THE COMPUTATION POINTS FROM MESH (AREA) C GIVEN THE BOUNDERIES OF THE GRID AND THE INTERVAL ALONG EACH C DIRECTION. C C INPUT: C """"""PHISW:LATITUDE OF THE SOUTH WEST CORNER OF GRID TO BE COMPUTED C IN DEGREES. C ALASW:LONGITUDE OF THE SOUTH WEST CORNER OF GRID TO BE COMPUTED C IN DEGREES. C PHINW:LATITUDE OF THE NORTH WEST CORNER OF GRID TO BE COMPUTED C IN DEGREES. C ALASE:LONGITUDE OF THE SOUTH EAST CORNER OF GRID TO BE COMPUTED C IN DEGREES. C XNTPHI:INTERVAL ALONG THE LATITUDINAL DIRECTION IN MINUTES C XNTALA:INTERVAL ALONG THE LONGITUDINAL DIRECTION IN MINUTES C C OUTPUT: C """"""" PHI:ARRAY OF LATITUDES OF COMPUTATION POINTS IN DEGREES. C ALA:ARRAY OF LONGITUDES OF COMPUTATION POINTS IN DEGREES. C K:NUMBER OF EVALUATION POINTS IN THE AREA. C IP:NUMBER OF COMPUTAION POINTS IN THE LATITUDINAL C DIRECTION. C IL:NUMBER OF COMPUTATION POINTS IN THE LONGITUDINAL C DIRECTION. C C THE MAXIMUM NUMBER OF POINTS THAT CAN BE STORED C IN THE ARRAYS PHI AND ALA ARE 10,000. C CSRB#1 CALL AREA(PHISW,ALASW,PHINW,ALASE,XNTPHI,XNTALA,PHI,ALA,K,IL,IP) LM=K C C________________________________________________________________________ C------------------------------------------------------------------------ C EXTERNAL LOOPING PROCESS FOR THE # COMPUTATION POINTS C """""""""""""""""""""""""""""""""""""""""""""""""""""" C THIS DO LOOP:STARTS THE EXTERNAL LOOPING PROCESS FOR THE NUMBER C OF COMPUTATION POINTS INVOLVED IN COMPUTING THE C TOPOGRAPHICAL EFFECT FOR THE INNER ZONE. C C LM:NUMBER OF COMPUTATION POINTS TO BE COMPUTED DETERMINED IN THE C SUBROUTINE AREA. C DO 5555 I=1,LM C C THE COMPUTATION POINT IS THEN INITIALIZED FROM THE ARRAYS C FORMED IN SUBROUTINE AREA. C C INPUT: PHI, ALA, ALA(I+1) C """""" C OUTPUT: C """"""" PHIO=LATITUDE OF COMPUTATION POINT IN DEGREES. C ALAO=LONGITUDE OF COMPUTATION POINT IN DEGREES. C ALA1=LONGITUDE OF THE NEXT COMPUTATION POINT IN DEGREES. C PHIO=PHI(I) ALAO=ALA(I) ALA1=ALA(I+1) C C------------------------------------------------------------------------ C LOADING OF THE 5X5 OR 5X10 MEAN GRAVITY ANOMALY VALUES C """""""""""""""""""""""""""""""""""""""""""""""""""""" C DETERMINING OF THE CELLSIZE C """"""""""""""""""""""""""""" C THE CELSZE: C SUBROUTINE DEPENDING ON THE LATITUDE OF THE COMPUTATION C POINT DECIDES UPON THE CELL SIZE TO BE USED,IN BOTH THE LATITUDINAL C AND LONGITUDINAL DIRECTIONS WITHIN THE 5'X5' OR 5'X10' MINUTE MEAN C GRAVITY ANOMALY DATA STRIPS. C THERE ARE 200 8X8 DEGREE MEAN GRAVITY ANOMALY DATA BLOCKS. C C NOTE: FOR PHIO < 54 DEGREES THE PROGRAM USES THE 5X5 MINUTE MEAN C """"" GRAVITY ANOMALY CELLS. C THERE ARE 9216 5X5 MGA IN EACH BLOCK. C C FOR PHIO > 54 DEGREES THE PROGRAM USES THE 5X10 MINUTE MEAN C GRAVITY ANOMALY CELLS. C THERE ARE 4608 5X10 MGA IN EACH BLOCK. C C INPUT: C """""" C PHIO:LATITUDE OF THE COMPUTATION POINT. C OUTPUT: C """"""" C XKPHII:CELL SIZE IN THE LATITUDINAL DIRECTION IN MINUTES. C XKALAI:CELL SIZE IN THE LONGITUDINAL DIRECTION IN MINUTES. C C SRB#2 CALL CELSZE(PHIO,XKPHII,XKALAI) C C TESTING OF THE LOADED FILE C""""""""""""""""""""""""""" C I=1:LOADS THE 8X8 DEGREE BLOCK OF 5X5 OR 5X10 MIN. MEAN GRAVITY C ANOMALIES OF THE FIRST COMPUTATION POINT FOR THE FIRST RUN. C NOTE: C """" ITOL1 IS USED FOR ALL OTHER COMPUTATION POINTS WITHIN THE LOOP. C C ITOL1=0:KEEP THE 8X8 DEGREE BLOCK OF 5X5 OR 5X10 MIN. MEAN C GRAVITY ANOMALIE CELLS. C ITOL1=1:LOAD A NEW 8X8 DEGREE BLOCK OF 5X5 OR 5X10 MIN. MEAN C GRAVITY ANOMALIE CELLS. C IF(I.EQ.1.OR.ITOL1.EQ.1)THEN C C IF ITOL=1---> THEN SUBROUTINE CHANG5 IS CALLED. C """"""""""""""""""""""""""""""""" C THE CHANG5: C SUBROUTINE FETCHES THE RIGHT 8X8 DEGREE BLOCK OF THE C 5X5 OR 5X10 MINUTE MEAN GRAVITY ANOMALY DATA. C INPUT: C """""" PHIO:LATITUDE OF THE COMPUTATION POINT. C ALAO:LONGITUDE OF THE COMPUTATION POINT. C OUTPUT: C """"""" M:IS THE CELL NUMBER OF THE COMPUTATION BLOCK. C SWPHII:LATITUDE OF THE SOUTH WEST CORNER OF THE C COMPUTATION BLOCK. C SWLAMI:LONGITUDE OF THE SOUTH WEST CORNER OF THE C COMPUTATION BLOCK. C PMAX5:LATITUDE OF THE NORTH EAST CORNER OF THE C COMPUTATION BLOCK. C DMAX5:LONGITUDE OF THE NORTH EAST CORNER OF THE C COMPUTATION BLOCK. C DGI:ARRAY CONTAINING THE 5X5 OR 5X10 MINUTE MEAN C GRAVITY ANOMALIES FOR THE 8X8 DEGREE BLOCK. C SDGI:ARRAY CONTAINING THE RMS VALUES OF THE 5X5 OR 5X10 C MINITE MEAN GRAVITY ANOMALIES FOR THE 8X8 DEGREE BLOCK. C ELVANI:ARRAY CONTAINING THE ASSOCIATED HEIGHTS OF THE C 5X5 OR THE 5X10 MINUTE MEAN GRAVITY ANOMALIES C FOR THE 8X8 DEGREE BLOCK. C C NOTE: IFFGG CAN BE EQUAL TO FIVE VALUES C """"" IFFGG=0 POINT LIES WITHIN THE BOUNDERIES OF THE 5X5 MINUTE C MEAN GRAVITY ANOMALY CELL. C IFFGG=1 POINT LIES OUTSIDE THE LATITUDE BOUNDERIES OF THE C 5X5 MINUTE MEAN GRAVITY ANOMALY CELL. C IFFGG=2 POINT LIES OUTSIDE THE LONGITUDINAL BOUNDERIES OF C THE 5X5 MINUTE MEAN GRAVITY ANOMALY CELL. C IFFGG=3 POINT LIES OUTSIDE THE LATITUDE BOUNDERIES OF THE C 5X10 MINUTE MEAN GRAVITY ANOMALY CELL. C IFFGG=4 POINT LIES OUTSIDE THE LONGITUDINAL BOUNDERIES OF C THE 5X10 MINUTE MEAN GRAVITY ANOMALY CELL. CSRB#3 C CALL CHANG5(PHIO,ALAO,K1,SWPHII,SWLAMI,PMAX5,DMAX5,DGI,SDGI, * ELVANI,IFFGG) WRITE(6,*)'CHANG5' C C FUNCTION #1 C THE DEGMIN FUNCTION CONVERTS FROM DEGREES TO MINUTES C SWPH1=SWPHII SWLAM1=SWLAMI SWPHII=DEGMIN(SWPHII) SWLAMI=DEGMIN(SWLAMI) C C C PN: C IS MAINLY USED IN THE GRID MODE . THE POINTS ARE COMPUTED C SEQUENTIALLY ALONG THE LATITUDES,WHEN THE COMPUTATION POINT C IS WITHIN 1 DEGREE FROM THE END OF THE 8X8 BLOCK,THEN A NEW C 5X5 OR 5X10 MINUTE MEAN GRAVITY ANOMALY BLOCK IS LOADED. C PN=DMAX5-1.0 ENDIF C IF--->IFFGG=0 C THEN-->ITOL1=0 C KEEP THE 5X5 OR 5X10 MINUTE MEAN GRAVITY ANOMALIES. C C ELSE-->IFFGG=1,2,3,4 C A)POINT LIES OUTSIDE OF BOUNDERY C B)SET ITOL1=1 C C)CHANGE THE EXISTING 5X5 OR 5X10 MINUTE MEAN GRAVITY C ANOMALY FILE. C D)EXIT PROGRAM AND C WRITE:PHI,ALA,AND DEFAULT VALUES AND GO PICK UP NEW C COMPUTATION POINT C IF(IFFGG.EQ.0)THEN ITOL1=0 C ELSE C ITOL1=1 WRITE(6,3334)PHIO,ALAO,FF GOTO 3011 ENDIF C C TEST TO SEE IF NEXT POINT LIES WITHIN 8X8 DEGREE BLOCK C """"""""""""""""""""""""""""""""""""""""""""""""""""""" C AN INNER BLOCK IS FORMED 1 DEGREE FROM THE BOUNDERY C OF THE 8X8 DEGREE COMPUTATION BLOCK.THE NEXT COMPUTATION C POINT IS THEN TESTED TO SEE IF A NEW 8X8 BLOCK SHOULD BE C LOADED. C C VARIABLES TESTED: C """"""""""""""""" C ALA1 :LONGITUDE OF THE NEXT COMPUTATION POINT. C PHI(I+1):LATITUDE OF THE NEXT COMPUTATION POINT. C PN :LONGITUDE OF THE EASTERLY BOUNDERY OF THE C COPMPUTATION BLOCK. C PMAX5 :LATITUDE OF THE NORTH WESTERLY BOUNDERY OF C THE COMPUTATION BLOCK. C (SWLAMI+1.0):LONGITUDE OF THE SOUTH WESTERLY CORNER OF C THE COMPUTATION BLOCK. C (SWPHII+1.0):LATITUDE OF THE SOUTH WESTERLY CORNER OF C THE COMPUTATION BLOCK. C C IF THE NEXT COMPUTATION POINT LIES OUTSIDE OF INNER BLOCK C THEN SET ITOL1=1.(LOAD NEW 5X5 OR 5X10 MINUTE MEAN GRAVITY C ANOMALY FILE) C ITOL1=0 IF(ALA1.GT.PN.OR.ALA1.LT.(SWLAM1+1.0))ITOL1=1 WRITE(6,*)ITOL1 WRITE(6,*)ALA1,PN,SWLAM1 C IF(PHI(I+1).GT.(PMAX5-1.0).OR.PHI(I+1).LT.(SWPH1+1.0))ITOL1=1 WRITE(6,*)ITOL1 WRITE(6,*)PHI(I+1),PMAX5,SWPH1 C C------------------------------------------------------------------------ C------------------------------------------------------------------------ C INNER ZONE C """""""""""" C C THE INBND:SUBROUTINE IS USED TO FORM THE BOUNDERIES OF THE C INNER ZONE C INPUT: C """""" OJ:IS THE 1/2 BLOCK SIZE OF THE INNER ZONE TO BE DETERMINED C PHIO:IS THE LATITUDE OF THE COMPUTATION POINT C ALAO:IS THE LONGITUDE OF THE COMPUTATION POINT C C OUTPUT: C """"""" PHIIS:IS THE SOUTHERN BORDER OF THE INNERMOST ZONE C PHIIN:IS THE NORTHERN BORDER OF THE INNERMOST ZONE C ALAIW:IS THE WESTERN BORDER OF THE INNERMOST ZONE C ALAIE:IS THE EASTERN BORDER OF THE INNERMOST ZONE C C NOTE: C """"" THE INPUT PARAMETERS ARE ALL IN DEGREES C THE OUTPUT PARAMETERS ARE ALL IN DEGREES C C C SRB#4 CALL INBND(0.5,PHIO,ALAO,PHIIS,PHIIN,ALAIW,ALAIE) WRITE(6,*)'THE BOUNDS ARE' WRITE(6,*) PHIIS,PHIIN,ALAIW,ALAIE C C NUMBER OF CELLS C """"""""""""""""" C THE CELNUM SUBROUTINE COMPUTES THE NUMBER OF CELLS IN THE C LATITUDINAL AND LONGITUDINAL DIRECTIONS RESPECTIVELY. C C INPUT: C """""" ILAT:IS THE INTERVAL IN SOUTH TO NORTH DIRECTION C ILONG:IS THE INTERVAL IN THE WEST TO EAST DIRECTION C C FROM SUBROUTINE CELSZE C """""""""""""""""""""" C XKPHII:CELLSIZE IN LATITUDINAL DIRECTION C XKALAI:CELLSIZE IN LONGITUDINAL DIRECTION C C NOTE: THE INPUT PARAMETERS ARE IN MINUTES C """"" C C OUTPUT: C """"""" XNPHI:IS THE NUMBER OF CELLS IN THE LATITUDINAL DIRECTION C XNALA:IS THE NUMBER OF CELLS IN THE LONGITUDINAL DIRECTION C SRB#5 CALL CELNUM(8,8,XKPHII,XKALAI,XNPHII,XNALAI) WRITE(6,*)'THE CELLSIZE IN THE LAT AND LON DIRECTION' WRITE(6,*) XKPHII,XKALAI,XNPHII,XNALAI C C________________________________________________________________________ C------------------------------------------------------------------------ C C SCAN AND RETAIN C """""""""""""""" C THE SCAN SUBROUTINE IS USED TO SCAN THE INNER ZONE C AND RETAIN THE EXACT NUMBER OF POINTS C ALONG WITH THE RESPECTIVE ARRAYS OF: C #1)LATITUDE OF POINT C #2)LONGITUDE OF POINT C #3)GRAVITY ANOMALIES OF POINT C #4)RMS OF GRAVITY ANOMALIES OF POINT C NOTE: C """"" POINT USED IN THE CALCULATION C C INPUT: C"""""" C DATA STRIP C """"""""""" C M:IS THE CELL NUMBER OF THE DATA POINTS IN THE STRIP C USED IN THE COMPUTATION C 576:IS THE MAXIMUM CELL NUMBER OF POINTS IN THE DATA STRIP C C FROM SUBROUTINE CHANG5 C """""""""""""""""""""" C SWPHII:IS THE SOUTHERN LAT. OF THE DATA STRIP IN MINUTES C SWLAMI:IS THE WESTERN LONG. OF THE DATA STRIP IN MINUTES C DGI:IS AN ARRAY CONTAINING THE 5X5 OR 5X10 MINUTE C MEAN GRAVITY ANOMALIES FOR THE 8X8 DEGREE DATA STRIP. C SDGI:IS THE CORRESPONDING RMS OF THE ANOMALIES. C ELVANI:IS THE CORRESPONDING HEIGHTS. C C FROM SUBROUTINE CELSZE C """""""""""""""""""""" C XKPHII:IS THE CELL SIZE IN LATITUDINAL DIRECTION IN MINUTES C XKALAI:IS THE CELL SIZE IN LONGITUDINAL DIRECTION IN MINUTES C C FROM SUBROUTINE CELNUM C """""""""""""""""""""" C XNALAI:IS THE NUMBER OF CELLS IN THE LONGITUDINAL DIRECTION C C FROM SUBROUTINE INBND C """""""""""""""""""""" C PHIIS:IS THE SOUTHERN BORDER OF THE INNER ZONE CELL C ALAIW:IS THE WESTERN BORDER OF THE INNER ZONE CELL C PHIIN:IS THE NORTHERN BORDER OF THE INNER ZONE CELL C ALAIE:IS THE EASTERN BORDER OF THE INNER ZONE CELL C NOTE: PARAMETERS INPUTED AS DEGREES C """""" C C OUTPUT: C """"""" C PHINZ:IS AN ARRAY OF THE LATITUDES IN THE STRIP C ALANZ:IS AN ARRAY OF THE LONGITUDES IN THE STRIP C DGINZ:IS THE CORRESPONDING ARRAY OF GRAVITY ANOMALIES C SDGINZ:IS THE CORRESPONDING ARRAY OF RMS OF THE ANOMALIES C HEINZ:IS THE CORRESPONDING ARRAY OF HEIGHTS C NUMBI:IS THE NUMBER OF POINTS IN THE STRIP C C INNER ZONE C ----------- C ISIGB:IS AN ERROR FLAG WITH TWO POSSIBLE VALUES C IF ISIGB=0 THERE ARE NO EMPTY CELLS OF 5'X 5' OR 5'X 10' C WITHIN THE INNER ZONE. C IF ISIGB=-5 THERE ARE EMPTY CELL OF 5'X 5'OR 5'X 10' C WITHIN THE INNER ZONE. C C SRB#6 CALL TSCAN(576,K1,SWPHII,SWLAMI,DGI,SDGI,ELVANI,PHIO,ALAO, * XKPHII,XKALAI,XNALAI,PHIIS,ALAIW,PHIIN,ALAIE, * PHINZ,ALANZ,DGINZ,SDGINZ,HEINZ,NUMBI,ISIGB) C C________________________________________________________________________ C C THE PA SUBROUTINE SELECTS THE RIGHT HEIGHT AND GRAVITY C ANOMALY OF THE COMPUTATION POINT. C C INPUT: C """""" C PHIO: IS THE LATITUDE OF THE COMPUTATION POINT C ALAO: IS THE LONGITUDE OF THE COMPUTATION POINT C PHINZ: IS AN ARRAY OF THE LATITUDES OF THE POINTS C OF THE GRID USED IN THE COMPUTATION C ALANZ: IS THE CORRESPONDING ARRAY OF LONGITUDES C DGINZ: IS THE CORRESPONDING ARRAY OF GRAVITY ANOMALIES C HEINZ: IS THE CORRESPONDING ARRAY OF HEIGHTS C NUMBI: IS THE NUMBER OF POINTS IN THE ABOVE ARRAYS C C OUTPUT: C """"""" C GA: IS THE GRAVITY ANOMALY OF THE POINT C HA: IS THE HEIGHT OF THE POINT C C SRB#7 CALL PA(PHIO,ALAO,PHINZ,ALANZ,DGINZ,HEINZ,NUMBI,GA,HA) C C------------------------------------------------------------------------ C TOPOGRAPHIC EFFECT C """"""""""""""""""" C FUNCTION #2 C C THE TOPATR FUNCTION IS USED TO COMPUTE THE CORRECTIONS C TO THE TOPOGRAPHICAL EFFECT ON THE GEOID. C C INPUT: C """""" C NUMBI: IS THE NUMBER OF POINTS USED IN THE COMPUTATION C PHIO: IS THE LATITUDE OF THE COMPUTATION POINT C ALAO: IS THE LONGITUDE OF THE COMPUTATION POINT C PHINZ: IS AN ARRAY OF LATITUDES OF THE POINTS USED IN C THE COMPUTATION C ALANZ: IS THE CORRESPONDING ARRAY OF LONGITUDES C DGINZ: IS THE CORRESPONDING ARRAY OF GRAVITY ANOMALIES C XKALAI: IS THE CELLSIZE IN THE LONGITUDINAL DIRECTION C XKPHII: IS THE CELLSIZE IN THE LATITUDINAL DIRECTION C HEINZ: IS AN ARRAY OF THE HEIGHTS USED IN THE COMPUTATION C GA: IS THE GRAVITY ANOMALY OF THE COMPUTATION POINT C HA: IS THE HEIGHT OF THE COMPUTATION POINT C C OUTPUT: IS THE CORRECTION TO THE TOPOGRAPHICAL EFFECT ON THE GEOID C""""""" C TOPIN=TOPATR(NUMBI,PHIO,ALAO,PHINZ,ALANZ,DGINZ,XKALAI, * XKPHII,HEINZ,GA,HA) TGA=GA+TOPIN WRITE(6,*)'THE COMPUTATION POINT IS' WRITE(6,*) PHIO,ALAO,GA,TOPIN,TGA WRITE(6,*) 'THE 5X5 MGA USED ARE' DO 77 MI=1,NUMBI WRITE(6,*)PHINZ(MI),ALANZ(MI),DGINZ(MI),HEINZ(MI) 77 CONTINUE C________________________________________________________________________ C------------------------------------------------------------------------ C THE BELOW CONTINUES SENDS YOU TO PICK UP A NEW COMPUTATION POINT. 3334 FORMAT(3F7.2) C 3011 CONTINUE C 5555 CONTINUE STOP END C END OF MAIN PROGRAM C """"""""""""""""""" C________________________________________________________________________ C------------------------------------------------------------------------ C SRB#1 C """"""" SUBROUTINE AREA(PHISW,ALASW,PHINW,ALASE,XNTPHI,XNTALA,PHI,ALA,K, * IL,IP) IMPLICIT REAL*4(A-H,O-Z) DIMENSION PHI(10000),ALA(10000) C C C THE AREA: C SUBROUTINE GENERATES THE COMPUTATION POINTS FROM MESH (AREA) C GIVEN THE BOUNDERIES OF THE GRID AND THE INTERVAL ALONG EACH C DIRECTION. C C INPUT: C """"""PHISW:LATITUDE OF THE SOUTH WEST CORNER OF GRID TO BE COMPUTED C ALASW:LONGITUDE OF THE SOUTH WEST CORNER OF GRID TO BE COMPUTED C PHINW:LATITUDE OF THE NORTH WEST CORNER OF GRID TO BE COMPUTED C ALASE:LONGITUDE OF THE SOUTH EAST CORNER OF GRID TO BE COMPUTED C XNTPHI:INTERVAL ALONG THE LATITUDINAL DIRECTION IN MINUTES C XNTALA:INTERVAL ALONG THE LONGITUDINAL DIRECTION IN MINUTES C C OUTPUT: C """"""" PHI:ARRAY OF LATITUDES OF COMPUTATION POINTS. C ALA:ARRAY OF LONGITUDES OF COMPUTATION POINTS. C K:NUMBER OF EVALUATION POINTS IN THE AREA C IP:NUMBER OF COMPUTAION POINTS IN THE LATITUDINAL DIRECTION C IL:NUMBER OF COMPUTATION POINTS IN THE LONGITUDINAL DIRECTION C C THE MAXIMUM NUMBER OF POINTS THAT CAN BE STORED C IN THE ARRAYS PHI AND ALA ARE 10,000. C XMIN=FLOAT(60) XMINE=(ALASE-ALASW)*XMIN IL=(XMINE/XNTALA)+1 XMINW=(PHINW-PHISW)*XMIN IP=(XMINW/XNTPHI)+1 K=IL*IP DO 10 I=1,IP DO 11 J=1,IL PHI((I-1)*IL+J)=PHISW+FLOAT(I-1)*XNTPHI/XMIN ALA((I-1)*IL+J)=ALASW+FLOAT(J-1)*XNTALA/XMIN 11 CONTINUE 10 CONTINUE RETURN END C C________________________________________________________________________ C------------------------------------------------------------------------ C SRB#2 C """""" SUBROUTINE CELSZE(PHIO,XKPHI,XKALA) IMPLICIT REAL*4(A-H,O-Z) C C THE CELSZE: C SUBROUTINE DEPENDING ON THE LATITUDE OF THE COMPUTATION C POINT DECIDES UPON THE CELL SIZE TO BE USED,IN BOTH THE LATITUDINAL C AND LONGITUDINAL DIRECTIONS WITHIN THE (5X5,5X10) MINUTE MEAN GRAVITY C ANOMALY DATA STRIPS. C THERE ARE 200 8X8 DEGREE MEAN GRAVITY ANOMALY DATA BLOCKS.EACH C BLOCK OVERLAPS WITH THE RIGHT AND UPPER ADJACENT BLOCKS BY 4 DEGREES. C THE RECORDS IN EACH BLOCK ARE SORTED IN THE ORDER OF INCREASING C LATITUDES,AND INCREASING LONGITUDES FOR POINTS OF EQUAL LATITUDE. C C NOTE: FOR PHIO < 54 DEGREES THE PROGRAM USES THE 5X5 MINUTE MEAN C """"" GRAVITY ANOMALY CELLS. C THERE ARE 9216 5X5 MGA IN EACH BLOCK. C C FOR PHIO > 54 DEGREES THE PROGRAM USES THE 5X10 MINUTE MEAN C GRAVITY ANOMALY CELLS. C THERE ARE 4608 5X10 MGA IN EACH BLOCK. C C C C INPUT: C """""" C PHIO:LATITUDE OF THE COMPUTATION POINT C OUTPUT: C """"""" C XKPHI:CELL SIZE IN LATITUDINAL DIRECTION IN MINUTES C XKALA:CELL SIZE IN THE LONGITUDINAL DIRECTION IN C MINUTES. C FIVE=FLOAT(5) TEN=FLOAT(10) FIFT4=FLOAT(54) IF(PHIO.LE.FIFT4)THEN XKPHI=FIVE XKALA=FIVE ELSE XKPHI=FIVE XKALA=TEN ENDIF RETURN END C C________________________________________________________________________ C------------------------------------------------------------------------ C SRB#3 C ----- C SUBROUTINE CHANG5(PC,DC,M,PSW,DSW,PMAX,DMAX,DG,SDG,H,IFFGG) C C C THE CHANG5:SUBROUTINE FETCHES THE RIGHT 8X8 DEGREE BLOCK OF THE C 5X5 OR 5X10 MINUTE MEAN GRAVITY ANOMALY DATA. C R.G. CHANG (JAN. 1986) C READS 5'X5' MEAN GRAVITY ANOMALIES FROM A.M12129.GMA C INPUT: C """""" PC:LATITUDE OF THE COMPUTATION POINT. C DC:LONGITUDE OF THE COMPUTATION POINT. C OUTPUT: C """"""" M:IS THE CELL NUMBER OF THE COMPUTATION BLOCK. C PSW:LATITUDE OF THE SOUTH WEST CORNER OF THE C COMPUTATION BLOCK. C DSW:LONGITUDE OF THE SOUTH WEST CORNER OF THE C COMPUTATION BLOCK. C PMAX:LATITUDE OF THE NORTH EAST CORNER OF THE C COMPUTATION BLOCK. C DMAX:LONGITUDE OF THE NORTH EAST CORNER OF THE C COMPUTATION BLOCK. C DG:ARRAY CONTAINING THE 5X5 OR 5X10 MINUTE MEAN C GRAVITY ANOMALIES FOR THE 8X8 DEGREE BLOCK. C SDG:ARRAY CONTAINING THE RMS VALUES OF THE 5X5 OR 5X10 C MINITE MEAN GRAVITY ANOMALIES FOR THE 8X8 DEGREE BLOCK. C H:ARRAY CONTAINING THE ASSOCIATED HEIGHTS OF THE C 5X5 OR THE 5X10 MINUTE MEAN GRAVITY ANOMALIES C FOR THE 8X8 DEGREE BLOCK. C C NOTE: IFFGG CAN BE EQUAL TO FIVE VALUES C """"" IFFGG=0 POINT LIES WITHIN THE BOUNDERIES OF THE 5X5 MINUTE C MEAN GRAVITY ANOMALY CELL. C IFFGG=1 POINT LIES OUTSIDE THE LATITUDE BOUNDERIES OF THE C 5X5 MINUTE MEAN GRAVITY ANOMALY CELL. C IFFGG=2 POINT LIES OUTSIDE THE LONGITUDINAL BOUNDERIES OF C THE 5X5 MINUTE MEAN GRAVITY ANOMALY CELL. C IFFGG=3 POINT LIES OUTSIDE THE LATITUDE BOUNDERIES OF THE C 5X10 MINUTE MEAN GRAVITY ANOMALY CELL. C IFFGG=4 POINT LIES OUTSIDE THE LONGITUDINAL BOUNDERIES OF C THE 5X10 MINUTE MEAN GRAVITY ANOMALY CELL. C IMPLICIT REAL*4 (A-H,O-Z) REAL*4 DG(9216) REAL*4 SDG(9216) REAL*4 H(9216) DIMENSION D1(48),D2(48),D3(48) INTEGER*4 REC INTEGER*2 I OPEN(1,STATUS='OLD',ACCESS='DIRECT',FORM='UNFORMATTED',RECL=576) IFFGG=0 IF(PC.LT.41.04165.OR.PC.GT.74.95835)THEN IFFGG=1 GO TO 50 ELSE END IF IF(DC.LT.215.04165.OR.DC.GT.316.95835) THEN IFFGG=2 GO TO 50 ELSE END IF IF(PC.GT.54.0.AND.DC.LT.214.0833) THEN IFFGG=3 GO TO 50 ELSE END IF IF(PC.GT.54.0.AND.DC.GT.316.9167) THEN IFFGG=4 GO TO 50 ELSE END IF DO 10 L=1,9216 DG(L)=0.0 SDG(L)=0.0 H(L)=0.0 10 CONTINUE I=ABS((PC-44.D0))/4.0+0.50 J=ABS((DC-218.D0))/4.0+0.50 KCC=44+4*(I-1) LCC=218+4*(J-1) TP=PC-FLOAT(KCC) TD=DC-FLOAT(LCC) IF(TP.GT.2.0) THEN I=I+1 ELSE END IF IF(TD.GT.2.0)THEN J=J+1 ELSE END IF IF(I.GE.8) THEN I=8 ELSE END IF IF(J.GE.25) THEN J=25 ELSE END IF PSW=40.D0+FLOAT((I-1)*4) DSW=214.D0+FLOAT((J-1)*4) PMAX=PSW+8.0 DMAX=DSW+8.0 IF(I.LE.3) THEN M=9216 KP=9216*(I-1)*25+9216*(J-1) LLL=(KP/48)+1 ELSE M=4608 KP=4608*(I-4)*25+4608*(J-1)+691200 LLL=(KP/48)+1 END IF IF(LLL.GT.14400) THEN GO TO 30 ELSE END IF DO 20 I=1,192 NO=LLL+I-1 READ(1,REC=NO) (D1(J),D2(J),D3(J),J=1,48) L=1+48*(I-1) K=L+47 DO 21 MN=L,K JN=MN-(I-1)*48 DG(MN)=D1(JN) SDG(MN)=D2(JN) H(MN)=D3(JN) 21 CONTINUE 20 CONTINUE GO TO 50 30 DO 40 I=1,96 NO=LLL+I-1 READ(1,REC=NO) (D1(J),D2(J),D3(J),J=1,48) L=1+48*(I-1) K=L+47 DO 31 MN=L,K JN=MN-(I-1)*48 DG(MN)=D1(JN) SDG(MN)=D2(JN) H(MN)=D3(JN) 31 CONTINUE 40 CONTINUE 50 RETURN END C C________________________________________________________________________ C------------------------------------------------------------------------ C SRB#4 C """""" SUBROUTINE INBND(OJ,PHIO,ALAO,PHIIS,PHIIN,ALAIW,ALAIE) C C THE INBND:SUBROUTINE IS USED TO FORM THE BOUNDERIES OF THE C INNER ZONE C INPUT: C """""" NJ:IS THE 1/2 BLOCK SIZE OF THE INNER ZONE TO BE DETERMINED C PHIO:IS THE LATITUDE OF THE COMPUTATION POINT C ALAO:IS THE LONGITUDE OF THE COMPUTATION POINT C C OUTPUT: C """"""" PHIIS:IS THE SOUTHERN BORDER OF THE INNERMOST ZONE C PHIIN:IS THE NORTHERN BORDER OF THE INNERMOST ZONE C ALAIW:IS THE WESTERN BORDER OF THE INNERMOST ZONE C ALAIE:IS THE EASTERN BORDER OF THE INNERMOST ZONE C C NOTE: C """"" THE INPUT PARAMETERS ARE ALL IN DEGREES C THE OUTPUT PARAMETERS ARE ALL IN DEGREES C IMPLICIT REAL*4(A-H,O-Z) C XNJ=OJ IPHI=PHIO IALA=ALAO XIPHI=IPHI XIALA=IALA RP=PHIO-XIPHI RL=ALAO-XIALA IF(RP.LT.0.5)THEN PHI=IPHI ELSE PHI=IPHI+1 ENDIF IF(RL.LT.0.5)THEN ALA=IALA ELSE ALA=IALA+1 ENDIF C PHIIS=PHI-XNJ PHIIN=PHI+XNJ ALAIW=ALA-XNJ ALAIE=ALA+XNJ RETURN END C C________________________________________________________________________ C------------------------------------------------------------------------ C SRB#5 C """"" SUBROUTINE CELNUM(ILAT,ILONG,XKPHI,XKALA,XNPHI,XNALA) C C THE CELNUM SUBROUTINE COMPUTES THE NUMBER OF CELLS IN THE C LATITUDINAL AND LONGITUDINAL DIRECTIONS RESPECTIVELY. C C INPUT: C """""" ILAT:IS THE INTERVAL IN SOUTH TO NORTH DIRECTION C ILONG:IS THE INTERVAL IN THE WEST TO EAST DIRECTION C C FROM SUBROUTINE CELSZE C """""""""""""""""""""" C XKPHI:CELLSIZE IN LATITUDINAL DIRECTION C XKALA:CELLSIZE IN LONGITUDINAL DIRECTION C C NOTE: THE INPUT PARAMETERS ARE IN MINUTES C """"" C C OUTPUT: C """"""" XNPHI:IS THE NUMBER OF CELLS IN THE LATITUDINAL DIRECTION C XNALA:IS THE NUMBER OF CELLS IN THE LONGITUDINAL DIRECTION C IMPLICIT REAL*4(A-H,O-Z) C XN=FLOAT(60) XILAT=FLOAT(ILAT) XILONG=FLOAT(ILONG) XNPHI=(XILAT*XN)/XKPHI XNALA=(XILONG*XN)/XKALA RETURN END C C C________________________________________________________________________ C------------------------------------------------------------------------ C SRB#6 C """""" C SUBROUTINE TSCAN(M,K,XPHISW,XALASW,DGARY,SDGARY,HARAY,PHIO,ALAO, * XKPHI,XKALA,XNALA,SPHOUT,WALOUT,XNPHOT,EALOUT, * RETPHI,RETALA,RETDG,RETSDG,RETHRY,NUMB,ISIGB) C C THE TSCAN:SUBROUTINE IS USED TO SCAN THE INNER ZONE AND RETAIN THE C EXACT NUMBER OF POINTS ALONG WITH THE RESPECTIVE C ARRAYS OF: C #1)LATITUDE OF POINT C #2)LONGITUDE OF POINT C #3)GRAVITY ANOMALIES OF POINT C #4)RMS OF GRAVITY ANOMALIES OF POINT C NOTE: C """"" POINT USED IN THE CALCULATION C C INPUT: C"""""" C DATA STRIP C """"""""""" C M:IS THE CELL NUMBER OF THE DATA STRIP USED IN THE COMPUTATION C K:IS THE NUMBER OF POINTS IN THE DATA STRIP C C FROM SUBROUTINE CHANG5 C """""""""""""""""""""" C XPHISW:IS THE SOUTHERN LAT. OF THE DATA STRIP IN DEGREES C XALASW:IS THE WESTERN LONG. OF THE DATA STRIP IN DEGREES C C FROM SUBROUTINE CHANG5 C """""""""""""""""""""" C DGARY:IS AN ARRAY CONTAINING THE 5X5 OR 5X10 MINUTE C MEAN GRAVITY ANOMALIES FOR THE 8X8 DEGREE DATA STRIP. C SDGARY:IS THE CORRESPONDING RMS OF THE ANOMALIES. C HARAY:IS THE CORRESPONDING HEIGHTS. C C FROM SUBROUTINE CELNUM C """""""""""""""""""""" C XKPHI:IS THE CELL SIZE IN LATITUDINAL DIRECTION C XKALA:IS THE CELL SIZE IN LONGITUDINAL DIRECTION C XNALA:IS THE NUMBER OF CELLS IN THE LONGITUDINAL DIRECTION C NOTE: PARAMETERS INPUTED AS DEGREES C """"" C C FROM SUBROUTINE INBND C """""""""""""""""""""" C SPHOUT:IS THE SOUTHERN BORDER OF THE INNER ZONE CELL C WALOUT:IS THE WESTERN BORDER OF THE INNER ZONE CELL C XNPHOT:IS THE NORTHERN BORDER OF THE INNER ZONE CELL C EALOUT:IS THE EASTERN BORDER OF THE INNER ZONE CELL C NOTE: PARAMETERS INPUTED AS DEGREES C """""" C C OUTPUT: C """"""" C RETPHI:IS AN ARRAY OF THE LATITUDES IN THE STRIP C RETALA:IS AN ARRAY OF THE LONGITUDES IN THE STRIP C RETDG:IS THE CORRESPONDING ARRAY OF GRAVITY ANOMALIES C RETSDG:IS THE CORRESPONDING ARRAY OF RMS OF THE ANOMALIES C RETHRY:IS THE CORRESPONDING ARRAY OF HEIGHTS C NUMB:IS THE NUMBER OF POINTS IN THE STRIP C C C INNER ZONE C ----------- C ISIGB:IS AN ERROR FLAG WITH TWO POSSIBLE VALUES C IF ISIGB=0 THERE ARE NO EMPTY CELLS OF 5'X 5' OR 5'X 10' C WITHIN THE INNER ZONE. C IF ISIGB=-5 THERE ARE EMPTY CELL OF 5'X 5'OR 5'X 10' C WITHIN THE INNER ZONE. C C IMPLICIT REAL*4(A-H,O-Z) DIMENSION DGARY(K),SDGARY(K),HARAY(K),RETPHI(M), * RETALA(M),RETDG(M),RETSDG(M),RETHRY(M) C SPHOUT=DEGMIN(SPHOUT) XNPHOT=DEGMIN(XNPHOT) WALOUT=DEGMIN(WALOUT) EALOUT=DEGMIN(EALOUT) C ISIGB=0 NALA=XNALA XM=FLOAT(60) X=FLOAT(1) Y=FLOAT(2) HALF=X/Y C C IPOS: THE POSITION OF POINT OF INTEREST IN DATA FILE C ------------------------------------------------------ C PHII=SPHOUT ALAI=WALOUT NUMB=0 564 IF(PHII.GE.XNPHOT)GOTO 200 N1=(PHII-XPHISW)/XKPHI N2=(ALAI-XALASW)/XKALA IPOS=N1*NALA+N2 ALAI2=ALAI 565 IF(ALAI2.GE.EALOUT)GOTO 300 IPOS=IPOS+1 C C ------------------------------------------ C TO CHECK WHETHER PHII&ALAI ARE ARE WITHIN C INNER BOUNDARIES. C ------------------------------------------ C C NUMB=NUMB+1 RETDG(NUMB)=DGARY(IPOS) IF(RETDG(NUMB).EQ.9999.0) THEN ISIGB=-5 C RETURN ENDIF FPHII=PHII/XM FALAI=ALAI2/XM RETPHI(NUMB)=FPHII RETALA(NUMB)=FALAI RETSDG(NUMB)=SDGARY(IPOS) RETHRY(NUMB)=HARAY(IPOS) C C C ------------------------- C TO SHIFT TO THE NEXT BLOCK C -------------------------- C C ALAI2=ALAI2+XKALA GOTO 565 300 PHII=PHII+XKPHI GOTO 564 200 CONTINUE RETURN END C C________________________________________________________________________ C------------------------------------------------------------------------ C SUBROUTINE #7 C """"""""""""" C THE PA SUBROUTINE SELECTS THE RIGHT HEIGHT AND GRAVITY ANOMALY C FOR THE COMPUTATION POINT C C INPUT: C """""" C PHIO: IS THE LATITUDE OF THE COMPUTATION POINT C ALAO: IS THE LONGITUDE OF THE COMPUTATION POINT C PHINZ: IS AN ARRAY OF LATITUDES OF THE POINTS USED C IN THE COMPUTATION C ALANZ: IS THE CORRESPONDING ARRAY OF LONGITUDES C DGINZ: IS THE CORRESPONDING ARRAY OF GRAVITY ANOMALIES C HEINZ: IS THE CORRESPONDING ARRAY OF HEIGHTS C NUMBER: IS THE NUMBER OF POINTS IN THE ABOVE ARRAYS C C OUTPUT: C """"""" C GA: IS THE GRAVITY ANOMALY OF THE COMPUTATION POINT C HA: IS THE HEIGHT OF THE COMPUTATION POINT C C SUBROUTINE PA(PHIO,ALAO,PHINZ,ALANZ,DGINZ,HEINZ,NUMBER,GA,HA) IMPLICIT REAL*4 (A-H,O-Z) DIMENSION PHINZ(NUMBER),ALANZ(NUMBER),DGINZ(NUMBER),HEINZ(NUMBER) C DO 777 I=1,NUMBER IF(PHIO.EQ.PHINZ(I).AND.ALAO.EQ.ALANZ(I))THEN C IF(PHIO.GT.PHINZ(I-1).AND.PHIO.LT.PHINZ(I+1) C * .AND.ALAO.GT.ALANZ(I-1).AND.ALAO.LT.ALANZ(I+1))THEN GA=DGINZ(I) HA=HEINZ(I) GO TO 7777 ENDIF 777 CONTINUE 7777 CONTINUE RETURN END C________________________________________________________________________ C------------------------------------------------------------------------ C FUNCTION #1 C """"""""""" C THE DEGMIN: FUNCTION IS USED TO CONVERT DEGREES TO MINUTES C C IT IS USED IN SRB #6 'TSCAN AND IN 'MAIN' C ----------------------------------------- C C INPUT: C """"" DEG:IS THE NUMBER OF DEGREES C C OUTPUT: IS DEG IN MINUTES C """""" C FUNCTION DEGMIN(DEG) IMPLICIT REAL*4(A-H,O-Z) C F=FLOAT(60) C DEGMIN=DEG*F RETURN END C________________________________________________________________________ C------------------------------------------------------------------------ C FUNCTION#2 C """"""""""" C THE TOPATR: FUNCTION COMPUTES THE CORRECTION TO THE C TOPOGRAPHICAL EFFECT ON THE GEOID C C IT IS USED IN 'MAIN' C ____________________ C C C INPUT: C """""" C NUMBER: IS THE NUMBER OF DUMMY POINTS IN THE ARRAYS PHIC, C ALAC,DGC AND H C PHIO: IS THE LATITUDE OF THE COMPUTATION POINT C ALAO: IS THE LONGITUDE OF THE COMPUTATION POINT C PHIC: IS AN ARRAY CONTAINING THE LATITUDES OF THE DUMMY POINTS C ALAC: IS THE CORRESPONDING ARRAY OF LONGITUDES C DGC:IS THE CORRESPONDING ARRAY OF GRAVITY ANOMALIES C XKL: IS THE CELLSIZE IN THE LONGITUDINAL DIRECTION C XKPH: IS THE CELLSIZE IN THE LATITUDINAL DIRECTION C H: IS AN ARRAY OF HEIGHTS OF THE DUMMY POINTS C GA: IS THE GRAVITY ANOMALY OF THE COMPUTATION POINT C HA: IS THE HEIGHT OF THE COMPUTATION POINT C NOTE: GA IS USED AS A FLAG WHICH TELLS YOU IF THERE IS AN EMPTY CELL C C OUTPUT: THE CORRECTION TO THE TOPOGRAPHICAL EFFECT ON THE GEOID C """"""" C C FUNCTION TOPATR(NUMBER,PHIO,ALAO,PHIC,ALAC,DGC,XKL,XKPH,H,GA,HA) IMPLICIT REAL*4(A-H,O-Z) DIMENSION PHIC(NUMBER),ALAC(NUMBER),DGC(NUMBER),H(NUMBER) C C=3654894.0 R=6371000.0 GT=0.0 IF(GA.EQ.9999.0)GOTO 2553 HA2=0.0 IF(HA.GT.0.0)HA2=HA**2 DO 12 I=1,NUMBER IF(DGC(I).EQ.9999.0)GOTO 2452 HDUM2=0.0 IF(H(I).GT.0.0)HDUM2=H(I)**2 EH=HDUM2-HA2 PSI=SPHER(PHIO,ALAO,PHIC(I),ALAC(I)) IF(PSI.LT.0.0001)GOTO 2452 DAREA=ELEMNT(PHIC(I),XKL,XKPH) XL=(R*PSI)**3 GT=GT+((EH*DAREA)/XL) 2452 CONTINUE 12 CONTINUE 2553 CONTINUE TOPATR=GT*C*100000.0 RETURN END C________________________________________________________________________ C------------------------------------------------------------------------ C FUNCTION #3 C """"""""""" C THE SPHER: FUNCTION COMPUTES THE SPHERICAL DISTANCE C BETWEEN TWO POINTS HAVING LATITUDE AND C LONGITUDE A,B & C,D RESPECTIVELY C C IT IS USED IN FUNCTION #2 'TOPATR' C ---------------------------------- C C INPUT: C """"" A:LATITUDE OF FIRST POINT IN DEGREES C B:LONGITUDE OF FIRST POINT IN DEGREES C C:LATITUDE OF OTHER POINT IN DEGREES C D:LONGITUDE OF OTHER POINT IN DEGREES C C OUTPUT: SPHERICAL DISTANCE BETWEEN THE TWO POINTS C """""" IN RADIANS C C FUNCTION SPHER(A,B,C,D) IMPLICIT REAL*4 (A-H,O-Z) PI=DARCOS(-1.D0) PR=PI/180.0 XA=A*PR XB=B*PR XC=C*PR XD=D*PR ASIN=SIN(XA) CSIN=SIN(XC) ACOS=COS(XA) CCOS=COS(XC) BD=COS(ABS(XB-XD)) SPHER=ARCOS(ASIN*CSIN+ACOS*CCOS*BD) RETURN END C________________________________________________________________________ C------------------------------------------------------------------------ C FUNCTION #4 C """"""""""" C THE ELEMNT: FUNCTION IS USED TO COMPUTE THE C SURFACE ELEMENT C C IT IS USED IN FUNCTION #2 'TOPATR' C ---------------------------------- C C INPUT: C PHI:LATITUDE OF THE POINT C DPHI:CELL SIZE IN THE LATITUDINAL DIRECTION C DLAM:CELL SIZE IN THE LONGITUDINAL DIRECTION C NOTE: THE INPUT ARE ALL IN MINUTES C C OUTPUT: THE SURFACE ELEMENT IN RADIANS C FUNCTION ELEMNT(PHI,DPHI,DLAM) IMPLICIT REAL*4 (A-H,O-Z) PI=DARCOS(-1.D0) F1=PI/180.0 F2=F1/60.0 XDPHI=DPHI*F2 XDLAM=DLAM*F2 XPHI=PHI*F1 XCOS=COS(XPHI) ELEMNT=XCOS*XDPHI*XDLAM RETURN END