C ********************************************************* C * * C * PROGRAM GIN * C * """"""""""""""""" * C * VERSION 4.0 * C * """""""""""" * C * APRIL 7 1987 * C * """""""""""" * C * * C * SURV. ENGINEERING DEPT. 1986 * C * * C * A SOFTWARE PACKAGE (IN SINGLE PRECISION) * C * FOR THE COMPUTATION OF A GRAVIMETRIC GEOID * C * USING A MODIFIED STOKES SPHEROIDAL KERNEL * C * TO DEGREE 20 , AN INTEGRATION CAP SIZE OF A * C * RADIUS OF 6 DEGREES AND A COMBINATION OF * C * POINT, MEAN 5X5 MIN., AND MEAN 1X1 DEG. * C * (REDUCED BY DEGREE 20) GRAVITY ANOMALIES. * C * """""""""""""""""""""""""""""""""" * C * THIS PROGRAM COMPUTES THE PARTIAL GRAVIMETRIC * C * GEOID FROM THREE CONTRIBUTIONS: * C * A. THE INNERMOST ZONE CONTRIBUTION * C * B. THE INNER ZONE CONTRIBUTION * C * C. THE OUTER ZONE CONTRIBUTION * C * * C ********************************************************* C_______________________________________________________________________ C INPUT : C CARD:1) THIS SECTION DETERMINES THE CORRECTIONS YOU WISH TO C APPLY TO THE GEOIDAL SURFACES.THE FLAGS USED C IN THE PROGRAM ARE EITHER A (1) FOR YES AND A C (0) FOR NO. C A) FIRST COLLUMN -->ATMOSPHERIC EFFECT. C B) THIRD COLLUMN -->INDIRECT EFFECT. C C) FIFTH COLLUMN-->TOPOGRAPHIC EFFECT FOR INNERMOST C AND INNER ZONES. C C CARD:2) IMODE : THE COMPUTATION MODE OPTION: C --- ZERO ( 0 ) IF POINT MODE IS USED C --- ONE ( 1 ) IF GRID MODE IS DESIRED C CARD:3) NUMDAT : THE NUMBER OF EXPECTED EVALUATION POINTS C WHEN POINT MODE IS SELECTED. C CARD:4) PSIO : THE INTEGRATION CAP SIZE C CARD:5) PHISW : C ALASW : THE S-W AND N-E CORNER COORDINATES OF THE C PHINW : AREA TO BE GRIDDED IN THE GRID MODE C ALASE : IN DEGREES C XNTPHI : THE GRID INTERVAL IN LAT. DIR. IN MINUTES C XNTALA : THE GRID INTERVAL IN LON. DIR. IN MINUTES C CARD:6) SWPHIO : THE S.-W. CORNER COORD. OF THE MEAN 1X1 DEG. C SWLAMO : GRAV. ANOMALY FILE USED IN THE COMPUTATIONS C THIS FILE IS SPECIFICALLY TAYLORED TO HAVE C A SIZE OF 50 DEG IN LATITUDE AND 150 DEG C IN LONGITUDE HERE. C OUTPUT : C 1) PHIO,ALAO : THE EVALUATION POINT COORDINATES. C 2) RESO,SIGMAO : THE INNERMOST ZONE CONTRIBUTION AND ITS RMS. C 3) RES1,SIGMA1 : THE INNER ZONE CONTRIBUTION AND ITS RMS. C 4) RES2,SIGMA2 : THE OUTER ZONE CONTRIBUTION AND ITS RMS. C 5) GEOID,GEODAC : THE TOTAL PARTIAL STOKES CONTR. AND ITS RMS. C C----------------------------------------------------------------------- C INITIALIZATION OF ARRAYS C """""""""""""""""""""""""" C----------------------------------------------------------------------- IMPLICIT REAL*4(A-H,O-Z) INTEGER ATM,IND,TINZ DIMENSION PHI(10000),ALA(10000) DIMENSION PHIARY(9000),ALAARY(9000),DG(9000),SDG(9000), * HETARY(9000) DIMENSION PHARY(9000),XLARAY(9000),FARAY(9000),SARAY(9000), * REHARY(9000) DIMENSION X(9000),Y(9000) DIMENSION C(1,6),ALFA(6),COV(6,6) DIMENSION DGI(9216),SDGI(9216),ELVANI(9216),TVANI(9216) DIMENSION PHINZ(576),ALANZ(576),DGINZ(576),SDGINZ(576),HEINZ(576), * TOPIN(576) DIMENSION DGO(7701),SDGO(7701),ELVANO(7701),TVANIO(7701) DIMENSION PHIOZ(576),ALAOZ(576),DGOZ(576),SDGOZ(576),HEOZ(576), * TOPOV(576) DIMENSION PHIA(4),ALI(4),GED(4),SGED(4),HED(4),TED(4) DIMENSION PHIA1(4),ALI1(4),GED1(4),SGED1(4),HED1(4),TED1(4) EXTERNAL FCT C______________________________________________________________________ C READING IN INITIAL FLAGS AND LIMITATIONS C---------------------------------------------------------------------- READ(5,*)ATM,IND,TINZ READ(5,*)IMODE READ(5,*)NUMDAT READ(5,*)PSIO READ(5,*)PHISW,ALASW,PHINW,ALASE,XNTPHI,XNTALA READ(5,*)SWPHIO,SWLAMO C________________________________________________________________________ C INITIALIZATION OF CONSTANTS C------------------------------------------------------------------------ PI=DARCOS(-1.D0) DR=PI/180.0 SCALE=720.0 DISTO=DR*PSIO DISTI=DISTO ITOL1=0 ITOL2=0 FF=-999.9 C_____________________________________________________________________ C IF IMODE = 0 CALL SRB#1 -DATA- TO READ EVALUATION POINTS C FOR POINT MODE C IF IMODE = 1 CALL SRB#2 -AREA- TO COMPUTE EVALUATION POINTS C FOR GRID MODE C----------------------------------------------------------------------- IF(IMODE.EQ.0)THEN CALL DATA(NUMDAT,NR,PHI,ALA) LM=NR ELSE CALL AREA(PHISW,ALASW,PHINW,ALASE,XNTPHI,XNTALA,PHI,ALA,K,IL,IP) LM=K ENDIF C________________________________________________________________________ C C LOAD THE 1X1 DEGREE MEAN GRAVITY ANOMALY(MGA) DATA FROM A C SEQUENTIALLY ORGANIZED FILE AND CREATE : C DGO : ARRAY CONTAINING THE 1X1 DEG. MGAS C SDGO : ARRAY CONTAINING THE ASSOCIATED RMS VALUES OF DGO C ELVANO : ARRAY CONTAINING THE ASSOCIATED HEIGHTS C C NOTE : THE MAXIMUM SIZE OF THESE ARRAYS IS 7701 C------------------------------------------------------------------------ DO 966 NM=1,7701 READ(33) DGO(NM),SDGO(NM),ELVANO(NM) 966 CONTINUE C________________________________________________________________________ C OUTER DO LOOP FOR 1 TO LM EVALUATION POINTS C C LM : NUMBER OF EVALUATION POINTS DETERMINED IN SRB#1 OR SRB#2 C PHIO : LATITUDE OF CURRENT EVALUATION POINT C ALAO : LONGITUDE OF CURRENT EVALUATION POINT C ALA1 : LONGITUDE OF THE NEXT EVALUATION POINT C------------------------------------------------------------------------ DO 5555 I=1,LM PHIO=PHI(I) ALAO=ALA(I) ALA1=ALA(I+1) C------------------------------------------------------------------------ C CALL SRB#3 -CELSZE- TO COMPUTE THE MGA CELL SIZE FOR PHIO C C NOTE: FOR PHIO BELOW 54 DEG. 5X5 MIN. MGA ARE USED C FOR PHIO ABOVE 54 DEG. 5X10 MIN. MGA ARE USED C THERE ARE 9216 CELLS IN THE MGA 8X8 DEG. DATA BLOCK C------------------------------------------------------------------------ CALL CELSZE(9216,PHIO,XKPHII,XKALAI) C________________________________________________________________________ C C MEAN GRAVITY ANOMALY(MGA) CONDITIONS: C IF I=1 LOAD THE FIRST MGA DATA BLOCK C IF ITOL1=0 KEEP THE MGA DATA BLOCK C IF ITOL1=1 LOAD A NEW MGA DATA BLOCK FOR THE NEXT EVALUATION POINT C IF IFFGG=0 EVALUATION POINT LIES WITHIN THE MGA COVERAGE C IF IFFGG.GT.0 EVALUATION POINT LIES OUTSIDE THE MGA COVERAGE C C IF I=1 OR ITOL1=1 CALL C C SRB#4 -CHANG5- TO COMPUTE THE MAX AND MIN. COORDINATES C AND READ THE DATA FOR THE APPROPRIATE C 5X5(10) MIN. MGA DATA BLOCK OF PHIO C TO SET THE FLAG IFFGG C C C IF IFFGG=0; LET ITOL1=0 C C ELSE IF IND=1; PRINT DEFAULT VALUES WITH INDIRECT EFFECT C C ELSE IF IND=0; PRINT DEFAULT VALUES WITHOUT INDIRECT EFFECT C C IF IND=0 OR 1 GO TO NEXT POINT C C IF THE EVALUATION POINT LIES WITHIN THE ONE DEGREE OVERLAP C DEFINED BY PN,(SWLAM+1),(PMAX-1), AND (SWPHII+1) THEN SET ITOL1=1 C C------------------------------------------------------------------------ IF(I.EQ.1.OR.ITOL1.EQ.1)THEN CALL CHANG5(PHIO,ALAO,M,SWPHII,SWLAMI,PMAX5,DMAX5,DGI,SDGI, * ELVANI,TVANI,IFFGG) PN=DMAX5-1.0 ENDIF IF(IFFGG.EQ.0)THEN ITOL1=0 ELSE IF(IND.EQ.1)THEN ITOL1=1 WRITE(6,3334)PHIO,ALAO,FF,FF,FF,FF,FF,FF,FF,FF,FF,FF,FF C WRITE(43,3334)PHIO,ALAO,FF,FF,FF,FF,FF,FF,FF,FF,FF,FF,FF GOTO 5555 ELSE ITOL1=1 WRITE(6,3334)PHIO,ALAO,FF,FF,FF,FF,FF,FF,FF,FF,FF,FF C WRITE(43,3334)PHIO,ALAO,FF,FF,FF,FF,FF,FF,FF,FF,FF,FF GOTO 5555 END IF ENDIF C IF(ALA1.GT.PN.OR.ALA1.LT.(SWLAMI+1.0))ITOL1=1 IF(PHI(I+1).GT.(PMAX5-1.0).OR.PHI(I+1).LT.(SWPHII+1.0))ITOL1=1 C________________________________________________________________________ C C POINT GRAVITY ANOMALY(PGA) DECISIONS: C IF I=1 LOAD THE THE FIRST PGA DATA BLOCK C IF ITOL2=0 KEEP THE PGA DATA BLOCK C IF ITOL2=1 LOAD A NEW PGA DATA BLOCK FOR THE NEXT C EVALUATION POINT C IF IFK=0 EVALUATION POINT LIES WITHIN THE PGA COVERAGE C IF IFK.GT.0 EVALUATION POINT LIES OUTSIDE THE PGA COVERAGE C C IF I=1 OR ITOL2=1 CALL C C SRB#5 -CHANGP- TO COMPUTE THE MAX. AND MIN. COORDINATES AND READ C IN THE DATA FOR THE APPROPRIATE PGA DATA BLOCK. C TO SET THE FLAG IFK C C IF IFK=0; LET ITOL2=0 C C ELSE IF IND=1; PRINT DEFAULT VALUES WITH INDIRECT EFFECT C C ELSE IF IND=0; PRINT DEFAULT VALUES WITHOUT INDIRECT EFFECT C C IF IND=0 OR 1 GO TO NEXT POINT C C IF THE EVALUATION POINT LIES WITHIN THE 5 AND 7.5 MIN OVERLAPS C DEFINED BY XN,ALAO,(PMAX-.0833), AND (PMAX-.1667)+(.0833/2) C LET ITOL2=1 C C------------------------------------------------------------------------ IF(I.EQ.1.OR.ITOL2.EQ.1)THEN CALL CHANGP(PHIO,ALAO,MK,PMAXP,DMAXP,DG,SDG,HETARY,PHIARY, * ALAARY,IFK) XN=DMAXP-0.083 ENDIF IF(IFK.EQ.0)THEN ITOL2=0 ELSE IF(IND.EQ.1)THEN ITOL2=1 WRITE(6,3334)PHIO,ALAO,FF,FF,FF,FF,FF,FF,FF,FF,FF,FF,FF C WRITE(43,3334)PHIO,ALAO,FF,FF,FF,FF,FF,FF,FF,FF,FF,FF,FF GOTO 5555 ELSE ITOL2=1 WRITE(6,3334)PHIO,ALAO,FF,FF,FF,FF,FF,FF,FF,FF,FF,FF C WRITE(43,3334)PHIO,ALAO,FF,FF,FF,FF,FF,FF,FF,FF,FF,FF GOTO 5555 END IF ENDIF C IF(ALA1.GT.XN.OR.ALA1.LT.ALAO.OR. & PHI(I+1).GT.(PMAXP-0.0833).OR. & PHI(I+1).LT.((PMAXP-0.1667)+(0.0833/2.0)))ITOL2=1 C________________________________________________________________________ C INNERMOST ZONE: CALL C """""""""""""" C SRB#6 -INMBND- TO DETERMINE THE BOUNDARIES OF THE INNERMOST ZONE C C SRB#7 -BOUND- TO COMPUTE THE INNER MOST ZONE INTEGRATION LIMITS C C SBR#8 -NUMERC- TO COMPUTE THE VALUES OF THE VECTOR C USED IN THE C INNERMOST ZONE NUMERICAL INTEGRATION C------------------------------------------------------------------------ CALL INMBND(5,PHIO,ALAO,PHICS,PHICN,ALACW,ALACE) CALL BOUND(PHIO,ALAO,PHICS,PHICN,ALACW,ALACE,X1,X2,Y1,Y2) CALL NUMERC(X1,X2,Y1,Y2,C) C________________________________________________________________________ C INNER ZONE: CALL C """"""""""" C SRB#9 -INBND- TO FORM THE BOUNDARIES OF THE INNER ZONE C C SRB#10 -CELNUM- TO DETERMINE THE NUMBER OF CELLS IN THE INNER ZONE C----------------------------------------------------------------------- CALL INBND(1,PHIO,ALAO,PHIIS,PHIIN,ALAIW,ALAIE) CALL CELNUM(8,8,XKPHII,XKALAI,XNPHII,XNALAI) C________________________________________________________________________ C INNERMOST AND INNER ZONE: CALL C """""""""""""""""""""""" C SRB#11 -CHECK- TO DETERMINE THE POSITION OF THE INNERMOST W.R.T. C THE INNER ZONE. C NOTE:A)THE INMST ZONE MUST BE AN INTEGER MULTIPLE OF THE CELL SIZE C B)THE INMST MUST BE WITHIN THE INNER ZONE C C SRB#12 -SCAN- TO FIND THE CELLS OF THE INNERMOST AND INNER ZONE C TO CHECK IF THE CELLS ARE EMPTY(NO MGAS) C NOTE: A)4 CELLS FOR THE INNERMOST ZONE. C B)572 CELLS FOR THE INNER ZONE. C------------------------------------------------------------------------ CALL CHECK(PHIIS,PHIIN,ALAIW,ALAIE,PHICS,PHICN,ALACW,ALACE, * PHIO,ALAO,XKPHII,XKALAI,SWPHII,SWLAMI,IFLG, * PHIISM,PHIINM,ALAIWM,ALAIEM,PHICSM,PHICNM,ALACWM, * ALACEM,SWPHIM,SWLAIM) C CALL SCAN(M,572,SWPHIM,SWLAIM,DGI,SDGI,ELVANI,TVANI,XKPHII,XKALAI, * XNALAI,PHIISM,ALAIWM,PHIINM,ALAIEM,PHICSM,PHICNM,ALACWM,ALACEM, * PHINZ,ALANZ,DGINZ,SDGINZ,HEINZ,TOPIN,NUMBI,ISIGA,ISIGB,PHIA,ALI, * GED,SGED,HED,TED,NUM2) C________________________________________________________________________ C OUTER ZONE: CALL C """""""""" C SRB#13 -OUTER- TO FORM THE BOUNDARIES OF THE OUTER ZONE C C SRB#3 -CELSZE- TO COMPUTE THE 1X1 DEG. MGA CELL SIZE FOR PHIO C C SRB#10 -CELNUM- TO DETERMINE THE NUMBER OF CELLS IN THE OUTER ZONE C C SRB#11 -CHECK- TO DETERMINE THE POSITION OF THE INNER W.R.T. C THE OUTER ZONE. C NOTE: A)THE INNER MUST BE INTEGER MULTIPLES OF THE CELL SIZE C B)THE INNER MUST BE WITHIN THE OUTER ZONE C------------------------------------------------------------------------ CALL OUTER(PSIO,PHIO,ALAO,PHIOS,PHION,ALAOW,ALAOE) CALL CELSZE(7701,PHIO,XKPHIO,XKALAO) CALL CELNUM(51,151,XKPHIO,XKALAO,XNPHIO,XNALAO) CALL CHECK(PHIOS,PHION,ALAOW,ALAOE,PHIIS,PHIIN,ALAIW,ALAIE, * PHIO,ALAO,XKPHIO,XKALAO,SWPHIO,SWLAMO,IFLG, * PHIOSM,PHIONM,ALAOWM,ALAOEM,PHIISM,PHIINM,ALAIWM, * ALAIEM,SWPHOM,SWLAMM) C________________________________________________________________________ C ERROR STATEMENTS FOR THE ERROR FLAG IFLG C------------------------------------------------------------------------ IF(IFLG.EQ.-1)THEN WRITE(6,555) ELSE IF(IFLG.EQ.-2)THEN WRITE(6,444) ELSE IF(IFLG.EQ.-3)THEN WRITE(6,333) ENDIF ENDIF ENDIF C C IFLG FORMAT STATMENTS C 555 FORMAT('0',5X,'BOUNDARIES NOT INTEGER MULTIPLES OF CELL SIZE') 444 FORMAT('0',5X,'INNER BOUNDARY IS NOT WITHIN OUTER BOUNDARY') 333 FORMAT('0',5X,'EVALUATION POINT NOT IN THE INNER BOUNDARY') C___________________________________________________________________ C CALL SRB#12 -SCAN- TO FIND THE CELLS OF THE INNER AND OUTER ZONE C TO CHECK IF THE CELLS ARE EMPTY(NO MGAS) C------------------------------------------------------------------------ CALL SCAN(7701,500,SWPHOM,SWLAMM,DGO,SDGO,ELVANO, * TVANIO,XKPHIO,XKALAO,XNALAO,PHIOSM,ALAOWM, * PHIONM,ALAOEM,PHIISM,PHIINM,ALAIWM,ALAIEM, * PHIOZ,ALAOZ,DGOZ,SDGOZ,HEOZ,TOPOV,NUMBO,ISIGA1,ISIGB1,PHIA1, * ALI1,GED1,SGED1,HED1,TED1,NUM21) C________________________________________________________________________ C IF ISIGA=-4 AND ISIGB=-5 C SET MFLAG=0 TO INDICATE THAT THERE ARE NO C MGAS IN THE 5X5(10) MIN. CELLS OF THE INNER AND INNERMOST ZONE C LOAD 1X1 DEG. MGAS TO REPLACE 5X5(10) MIN. CELLS C IF MFLAG=-9 THEN THERE ARE NO 1X1 DEG. MGAS. C C ELSE IF IND=1 PRINT DEFAULT VALUES WITH INDIRECT EFFECT C C ELSE IF IND=0 PRINT DEFAULT VALUES WITHOUT INDIRECT EFFECT C C IF IND=1 OR 0 GO TO NEXT POINT C------------------------------------------------------------------------ IF(ISIGA.EQ.-4.OR.ISIGB.EQ.-5) THEN MFLAG=0 DO 6661 IOZ=1,NUM21 IF(SGED1(IOZ).EQ.50.0)MFLAG=-9 6661 CONTINUE IF(MFLAG.EQ.-9)THEN IF(IND.EQ.1)THEN WRITE(6,3334) PHIO,ALAO,FF,FF,FF,FF,FF,FF,FF,FF,FF C WRITE(43,3334) PHIO,ALAO,FF,FF,FF,FF,FF,FF,FF,FF,FF GOTO 5555 ELSE WRITE(6,3334) PHIO,ALAO,FF,FF,FF,FF,FF,FF,FF,FF C WRITE(43,3334) PHIO,ALAO,FF,FF,FF,FF,FF,FF,FF,FF GOTO 5555 END IF ENDIF C PHIM=PHIIS+1.0 DLAMM=ALAIW+1.0 C________________________________________________________________________ C IF ISIGA=-4 REPLACE 5X5(10) MIN. WITH A 1X1 DEG. MGA C FOR THE INNERMOST ZONE CELLS C------------------------------------------------------------------------ IF(ISIGA.EQ.-4) THEN DO 5551 INMZ=1,NUM2 IF(GED(INMZ).EQ.9999.0.AND. * PHIA(INMZ).LT.PHIM.AND. * ALI(INMZ).LE.DLAMM) THEN GED(INMZ)=GED1(1) SGED(INMZ)=50.0 HED(INMZ)=HED1(1) TED(INMZ)=0.0 ELSEIF(GED(INMZ).EQ.9999.0.AND. * PHIA(INMZ).LT.PHIM.AND. * ALI(INMZ).GT.DLAMM) THEN GED(INMZ)=GED1(2) SGED(INMZ)=50.0 HED(INMZ)=HED1(2) TED(INMZ)=0.0 ELSEIF(GED(INMZ).EQ.9999.0.AND. * PHIA(INMZ).GT.PHIM.AND. * ALI(INMZ).LE.DLAMM) THEN GED(INMZ)=GED1(3) SGED(INMZ)=50.0 HED(INMZ)=HED1(3) TED(INMZ)=0.0 ELSEIF(GED(INMZ).EQ.9999.0.AND. * PHIA(INMZ).GT.PHIM.AND. * ALI(INMZ).GT.DLAMM) THEN GED(INMZ)=GED1(4) SGED(INMZ)=50.0 HED(INMZ)=HED1(4) TED(INMZ)=0.0 ENDIF 5551 CONTINUE ENDIF C________________________________________________________________________ C IF ISIGA=-5 REPLACE 5X5(10) MIN. WITH A 1X1 DEG. MGA C FOR THE INNER ZONE CELLS C------------------------------------------------------------------------ IF(ISIGB.EQ.-5) THEN DO 5552 INZ=1,NUMBI IF(DGINZ(INZ).EQ.9999.0.AND. * PHINZ(INZ).LT.PHIM.AND. * ALANZ(INZ).LE.DLAMM) THEN DGINZ(INZ)=GED1(1) SDGINZ(INZ)=50.0 HEINZ(INZ)=HED1(1) TOPIN(INZ)=0.0 ELSEIF(DGINZ(INZ).EQ.9999.0.AND. * PHINZ(INZ).LT.PHIM.AND. * ALANZ(INZ).GT.DLAMM) THEN DGINZ(INZ)=GED1(2) SDGINZ(INZ)=50.0 HEINZ(INZ)=HED1(2) TOPIN(INZ)=0.0 ELSEIF(DGINZ(INZ).EQ.9999.0.AND. * PHINZ(INZ).GT.PHIM.AND. * ALANZ(INZ).LE.DLAMM) THEN DGINZ(INZ)=GED1(3) SDGINZ(INZ)=50.0 HEINZ(INZ)=HED1(3) TOPIN(INZ)=0.0 ELSEIF(DGINZ(INZ).EQ.9999.0.AND. * PHINZ(INZ).GT.PHIM.AND. * ALANZ(INZ).GT.DLAMM) THEN DGINZ(INZ)=GED1(4) SDGINZ(INZ)=50.0 HEINZ(INZ)=HED1(4) TOPIN(INZ)=0.0 ENDIF 5552 CONTINUE ENDIF ENDIF C________________________________________________________________________ C SURFACE FITTING OF THE INNERMOST ZONE C """"""""""""""""""""""""""""""""""""" C IF MK=O THERE ARE NO PGAS C USE MGAS INSTEAD OF PGAS FOR THE INNERMOST ZONE C ELSE C USE PGAS FOR THE INNERMOST ZONE BY CALLING C C CALL SRB#14 -FIND- TO FIND THE PGAS OF THE INNERMOST ZONE C C------------------------------------------------------------------------ IF(MK.EQ.0)THEN NUM=NUM2 ELSE CALL FIND(MK,PHIO,ALAO,PHIARY,ALAARY,DG,SDG,HETARY,PHICS,PHICN, * ALACW,ALACE,NUM,PHARY,XLARAY,FARAY,SARAY,REHARY) ENDIF C________________________________________________________________________ C IF ATM=1 AND NUM=0 APPLY THE ATMOSPHERIC CORRECTIONS TO C C 1)THE INNERMOST ZONE USING PGAS C 2)THE INNERMOST ZONE USING 5X5(10) MIN. MGAS C 3)THE INNER ZONE USING 5X5(10) MIN MGAS C C BY CALLING FUNC#7 -ATMATR- TO COMPUTE CORRECTIONS C------------------------------------------------------------------------ IF(ATM.EQ.1)THEN IF(NUM.GT.0)THEN DO 8334 JVV=1,NUM IF(FARAY(JVV).EQ.9999.0)GOTO 8333 CGA=ATMATR(REHARY(JVV)) FARAY(JVV)=FARAY(JVV)+CGA 8334 CONTINUE ENDIF C DO 8333 JV=1,NUM2 IF(GED(JV).EQ.9999.0)GOTO 8333 CGA=ATMATR(HED(JV)) GED(JV)=GED(JV)+CGA 8333 CONTINUE C DO 8335 JW=1,NUMBI IF(DGINZ(JW).EQ.9999.0)GOTO 8334 CGA=ATMATR(HEINZ(JW)) DGINZ(JW)=DGINZ(JW)+CGA 8335 CONTINUE ENDIF C________________________________________________________________________ C IF TINZ=1 AND NUM=0 APPLY THE TOPOGRAPHIC CORRECTIONS TO C C 1)THE INNERMOST ZONE USING PGAS C 2)THE INNERMOST ZONE USING 5X5(10) MIN. MGAS C 3)THE INNER ZONE USING 5X5(10) MIN MGAS C------------------------------------------------------------------------ IF(TINZ.EQ.1)THEN IF(NUM.GT.0)THEN DO 8555 KIV=1,NUM IF(PHARY(KIV).LE.PHIO.AND.XLARAY(KIV).LE.ALAO)THEN FARAY(KIV)=FARAY(KIV)+TED(1) ELSEIF(PHARY(KIV).LT.PHIO.AND.XLARAY(KIV).GT.ALAO)THEN FARAY(KIV)=FARAY(KIV)+TED(2) ELSEIF(PHARY(KIV).GT.PHIO.AND.XLARAY(KIV).LE.ALAO)THEN FARAY(KIV)=FARAY(KIV)+TED(3) ELSEIF(PHARY(KIV).GT.PHIO.AND.XLARAY(KIV).GT.ALAO)THEN FARAY(KIV)=FARAY(KIV)+TED(4) ENDIF 8555 CONTINUE C DO 8336 KLV=1,NUM2 GED(KLV)=GED(KLV)+TED(KLV) 8336 CONTINUE ENDIF C DO 8338 KIW=1,NUMBI DGINZ(KIW)=DGINZ(KIW)+TOPIN(KIW) 8338 CONTINUE ENDIF C_______________________________________________________________________ C CALL FUNC#9 -MEAN- TO COMPUTE THE MEAN HEIGHT C C IF IND=1 CALL FUNC#8 -ENDEFT- TO COMPUTE THE INDIRECT EFFECT C C NOTE: THE INDIRECT EFFECT IS CALCULATED USING THE HEIGHTS OF C THE FOUR 5X5 MIN. OR TWO 5X10 MIN. MGAS CLOSEST TO THE C EVALUATION POINT C---------------------------------------------------------------------- HMEAN=MEAN(NUM2,HED) IF(IND.EQ.1)THEN XNDFTI=ENDEFT(NUM2,PHIO,ALAO,PHIA,ALI,XKALAI,XKPHII,HED * ,HMEAN) ENDIF C________________________________________________________________________ C THIS SECTION PERFORMS THE SURFACE FITTING FOR THE C INNERMOST ZONE BY CALLING C C SRB#15 -COM- TO COMPUTE THE LOCAL CARTESIAN COORDINATES C OF THE INNERMOST ZONE GRAVITY ANOMALIES. C C FUNC#3 -ICOUNT- TO DETERMINE THE NUMBER OF PARAMETERS TO SOLVED C FOR IN THE SURFACE FITTING. C C SRB#16 -SURF- TO CHECK THE GEOMETRY OF THE PGAS AND PERFORM C A LEAST SQUARE'S FIT TO THE INNERMOST ZONE. C C THE NEXT 3 CONDITIONS DECIDE ON WHETHER TO USE PGAS OR C 5X5(10) MIN. MGAS FOR THE INNERMOST ZONE C C A)IF THE NUMBER OF PGA VALUES GREATER THAN FOUR CALL COM C B)IF THE NUMBER OF PGA VALUES IS EQUAL TO ZERO OR LESS THAN C THE NUMBER OF 5X5(10) MGAS USE THE MGA VALUES FOR COM C C)IF PHIO GT 54 DEG. AND THE NUMBER OF PGAS IS LT 4 USE THE C MGA VALUES FOR COM C C THE PROGRAM IS THEN DIVIDED INTO THREE PARTS: C C PART#1) C"""""""" ------>NUM=0,1,2,3,4 C | C TWO POSSIBILITIES C C A)--->PHIO < 54 DEGREES C THE 4(5X5)MGA ARE LOADED AND THE SURFACE C IS FITTED SOLVING FOR NP=3 COEFFICIENTS. C C C B)--->PHIO > 54 DEGREES C THE 2(5X10)MGA ARE LOADED AND THE SURFACE C IS FITTED SOLVING FOR NP=1 COEFFICIENT. C C PART#2) C"""""""" ------>NUM=5,6 C THE DISTRIBUTION OF THE PGAS IS CHECKED WITHIN THE SURFACE C FITTING PROGRAM.IF THE DISTRIBUTION OF THE PGAS IS O.K.THEN THE C SURFACE PROGRAM SOLVES FOR 4 COEFFICIENTS(NP=4). C C IF A PROBLEM EXIST WITH THE GEOMETRY OR THE DISTRIBUTION C OF THE PGAS THE PROGRAM OPTS FOR A LOWER ORDER SURFACE,CHOOSES C 4 PGAS AND SOLVES FOR 3 COEFFICIENTS.(NUM=4,NP=3) C | C TWO POSSIBILITIES C C A)--->PHIO < 54 DEGREES C THE 4(5X5)MGA ARE LOADED AND THE SURFACE C IS FITTED SOLVING FOR NP=3 COEFFICIENTS. C C C B)--->PHIO > 54 DEGREES C THE 2(5X10)MGA ARE LOADED AND THE SURFACE C IS FITTED SOLVING FOR NP=1 COEFFICIENT. C C PART #3 C"""""""" ------>NUM=6> C THE DISTRIBUTION OF THE PGAS ARE CHECK WITHIN THE SURFACE C FITTING PROGRAM.IF THE DISTRIBUTION OF THE PGAS IS O.K.THEN THE C SURFACE PROGRAM SOLVES FOR 6 COEFFICIENTS(NUM=6>,NP=6). C C IF A PROBLEM EXIST WITH THE GEOMETRY OR THE DISTRIBUTION OF C THE PGAS,THE PROGRAM OPTS FOR A LOWER ORDER SURFACE,CHOOSES 5 OR 6 C PGAS AND SOLVES FOR 4 COEFFICIENTS(NUM=5,6 NP=4). C C IF A PROBLEM STILL EXIST WITH THE GEOMETRY OR THE DISTRIBUTION C OF THE PGAS,THE PROGRAM OPTS FOR A LOWER ORDER SURFACE,CHOOSES C 4(5X5) OR THE 2(5X10) MGA AND SOLVES FOR 3 OR 1 COEFFICIENTS. C | (NUM=4,NP=3) C | (NUM=2,NP=1) C TWO POSSIBILITIES C C A)--->PHIO < 54 DEGREES C THE 4(5X5)MGA ARE LOADED AND THE SURFACE C IS FITTED SOLVING FOR NP=3 COEFFICIENTS. C C C B)--->PHIO > 54 DEGREES C THE 2(5X10)MGA ARE LOADED AND THE SURFACE C IS FITTED SOLVING FOR NP=1 COEFFICIENT. C------------------------------------------------------------------------ IF(NUM.GT.4)GO TO 7771 IF(NUM.EQ.0.OR.NUM.LE.NUM2)NUM=NUM2 IF(PHIO.GE.54.0.AND.NUM.LE.4)NUM=NUM2 C DO 9012 KII=1,NUM PHARY(KII)=PHIA(KII) XLARAY(KII)=ALI(KII) FARAY(KII)=GED(KII) SARAY(KII)=SGED(KII) REHARY(KII)=HED(KII) 9012 CONTINUE C 7771 CALL COM(PHIO,ALAO,PHARY,XLARAY,NUM,X,Y,X2,Y2) C________________________________________________________________________ C PART#1 C------------------------------------------------------------------------ NP=ICOUNT(NUM) C IF(NP.EQ.1)THEN CALL SURF(NUM,NP,X,Y,PHARY,XLARAY,FARAY,SARAY,REHARY,COV,ALFA, & DET,PHIO,ALAO,NFLAG) ENDIF C IF(NP.EQ.3)THEN CALL SURF(NUM,NP,X,Y,PHARY,XLARAY,FARAY,SARAY,REHARY,COV,ALFA, & DET,PHIO,ALAO,NFLAG) ENDIF C________________________________________________________________________ C PART#2 C------------------------------------------------------------------------ IF(NP.EQ.4)THEN CALL SURF(NUM,NP,X,Y,PHARY,XLARAY,FARAY,SARAY,REHARY,COV,ALFA, & DET,PHIO,ALAO,NFLAG) IF(NFLAG.EQ.-9)THEN NUM=NUM2 DO 9013 KJI=1,NUM PHARY(KJI)=PHIA(KJI) XLARAY(KJI)=ALI(KJI) FARAY(KJI)=GED(KJI) SARAY(KJI)=SGED(KJI) REHARY(KJI)=HED(KJI) 9013 CONTINUE CALL COM(PHIO,ALAO,PHARY,XLARAY,NUM,X,Y,X2,Y2) NP=ICOUNT(NUM) CALL SURF(NUM,NP,X,Y,PHARY,XLARAY,FARAY,SARAY,REHARY,COV,ALFA, & DET,PHIO,ALAO,NFLAG) ENDIF ENDIF C________________________________________________________________________ C PART#3 C------------------------------------------------------------------------ IF(NP.EQ.6)THEN CALL SURF(NUM,NP,X,Y,PHARY,XLARAY,FARAY,SARAY,REHARY,COV,ALFA, & DET,PHIO,ALAO,NFLAG) IF (NFLAG.EQ.-9)THEN NP=4 CALL SURF(NUM,NP,X,Y,PHARY,XLARAY,FARAY,SARAY,REHARY,COV,ALFA, & DET,PHIO,ALAO,NFLAG) IF (NFLAG.EQ.-9)THEN NUM=NUM2 DO 9014 KJK=1,NUM PHARY(KJK)=PHIA(KJK) XLARAY(KJK)=ALI(KJK) FARAY(KJK)=GED(KJK) SARAY(KJK)=SGED(KJK) REHARY(KJK)=HED(KJK) 9014 CONTINUE CALL COM(PHIO,ALAO,PHARY,XLARAY,NUM,X,Y,X2,Y2) NP=ICOUNT(NUM) CALL SURF(NUM,NP,X,Y,PHARY,XLARAY,FARAY,SARAY,REHARY,COV,ALFA, & DET,PHIO,ALAO,NFLAG) ENDIF ENDIF ENDIF C------------------------------------------------------------------------ C DETERMINING THE CONTRIBUTION FROM THE THREE ZONES BY CALLING C C SRB#17 -INNER- TO CALCULATE THE INNERMOST CONTRIBUTIONS AND C RMS VALUES C SRB#18 -XNTGRT- TO INTEGRATE OVER THE INNER OR OUTER ZONE C TO DETERMINE THE CONTRIBUTION AND ITS RMS VALUES. C------------------------------------------------------------------------ CALL INNER(C,NUM,NP,ALFA,COV,RESO,SIGMAO) CALL XNTGRT(NUMBI,RES1,SIGMA1,FCT,PHINZ,ALANZ,DGINZ,SDGINZ, * HEINZ,PHIO,ALAO,DISTI,XKALAI,XKPHII) CALL XNTGRT(NUMBO,RES2,SIGMA2,FCT,PHIOZ,ALAOZ,DGOZ,SDGOZ, * HEOZ,PHIO,ALAO,DISTO,XKALAO,XKPHIO) C________________________________________________________________________ C THE RESULTS ARE COMBINED BY CALLING C C SRB#19 -TOTAL- TO COMPUTE THE GEOIDAL HEIGHT AND ITS RMS VALUE C FROM THE THREE ZONES. C C IF IND EQ TO 1 PRINT THE RESULTS WITH THE INDIRECT EFFECT C IF IND EQ TO 0 PRINT THE RESULTS WITHOUT THE INDIRECT EFFECT C C LOGICAL UNIT DESCRIPTIONS C 6 : ON DATA SET 107 C 43 : ON DISK C------------------------------------------------------------------------ IF(IND.EQ.1)THEN CALL TOTAL(RESO,RES1,RES2,SIGMAO,SIGMA1,SIGMA2,GEOID,GEODAC) GEOID=GEOID+XNDFTI WRITE(6,3336) PHIO,ALAO,RESO,SIGMAO,RES1,SIGMA1,RES2,SIGMA2, * GEOID,GEODAC,XNDFTI C WRITE(43,3336) PHIO,ALAO,RESO,SIGMAO,RES1,SIGMA1,RES2,SIGMA2, C * GEOID,GEODAC,XNDFTI ELSE CALL TOTAL(RESO,RES1,RES2,SIGMAO,SIGMA1,SIGMA2,GEOID,GEODAC) WRITE(6,3334) PHIO,ALAO,RESO,SIGMAO,RES1,SIGMA1,RES2,SIGMA2, * GEOID,GEODAC C WRITE(43,3334) PHIO,ALAO,RESO,SIGMAO,RES1,SIGMA1,RES2,SIGMA2, C * GEOID,GEODAC C________________________________________________________________________ C FINAL RESULT FORMAT STATEMENTS C------------------------------------------------------------------------ 3336 FORMAT(11F7.2) 3334 FORMAT(10F7.2) ENDIF C________________________________________________________________________ C CONTINUE OUTER LOOP C------------------------------------------------------------------------ 5555 CONTINUE STOP END C END OF MAIN PROGRAM C """"""""""""""""""" C________________________________________________________________________ C------------------------------------------------------------------------ C SRB#1 C """"""" SUBROUTINE DATA(NM,NR,PHI,ALA) IMPLICIT REAL*4(A-H,O-Z) DIMENSION PHI(10000),ALA(10000) C C DATA READS THE LIST OF EVALUATION POINTS MERGED C TO THE BOTTOM OF THE PROGRAM. C C INPUT: C """"" C FROM MAIN PROGRAM C """"""""""""""""" C NUMDAT:NUMBER OF EVALUATION POINTS TO BE READ. C C OUTPUT: C """""" NR:NUMBER OF EVALUATION POINTS TO BE READ. C PHI:ARRAY CONTAINING THE LATITUDES OF THE EVALUATION POINTS C IN DECIMAL DEGREES. C ALA:ARRAY CONTAINING THE LONGITUDES OF THE EVALUATION POINTS. C IN DECIMAL DEGREES. C NR=NM READ(5,*)(PHI(K),ALA(K),K=1,NM) RETURN END C C________________________________________________________________________ C------------------------------------------------------------------------ C SRB#2 C """"""" SUBROUTINE AREA(PHISW,ALASW,PHINW,ALASE,XNTPHI,XNTALA,PHI,ALA,K, * IL,IP) C C AREA GENERATES THE EVALUATION POINTS GIVEN THE BOUNDARIES C OF THE GRID AND THE INTERVAL ALONG THE LAT. ANF LONG. DIRECTIONS. C C INPUT: C """"" C C FROM MAIN PROGRAM C """"""""""""""""" 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 EVALUATION POINTS. C ALA:ARRAY OF LONGITUDES OF EVALUATION POINTS. C K:NUMBER OF EVALUATION POINTS IN THE AREA C IP:NUMBER OF EVALUATION POINTS IN THE LAT. DIRECTION C IL:NUMBER OF EVALUATION POINTS IN THE LONG. DIRECTION C IMPLICIT REAL*4(A-H,O-Z) DIMENSION PHI(10000),ALA(10000) C C COMPUTES THE # OF POINTS IN THE C LAT. AND LONG. DIRECTIONS, AND THE TOTAL AREA XMIN=FLOAT(60) XMINE=(ALASE-ALASW)*XMIN IL=(XMINE/XNTALA)+1 XMINW=(PHINW-PHISW)*XMIN IP=(XMINW/XNTPHI)+1 K=IL*IP C C COMPUTES THE EVALUATION POINTS C 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#3 C """""" SUBROUTINE CELSZE(NO,PHIO,XKPHI,XKALA) C C CELSZE COMPUTES THE CELL LAT. AND LONG. CELL DIMENSIONS FOR THE C THREE TYPES OF MGA CELLS: 5X5 MIN, 5X10 MIN, AND 1X1 DEG.. C C INPUT: C """"" C FROM MAIN PROGRAM: C """"""""""""""""" C NO:THE NUMBER OF 5X5 OR 5X10 MINUTE CELLS IN THE C 8X8 DEGREE MEAN GRAVITY ANOMALY DATA BLOCK. C C FROM DATA OR AREA: C """"""""""""""""" C PHIO:LATITUDE OF THE EVALUATION POINT C C OUTPUT: C """"""" C XKPHI:CELL SIZE IN LATITUDINAL DIRECTION IN MINUTES C XKALA:CELL SIZE IN THE LONGITUDINAL DIRECTION IN C MINUTES. C IMPLICIT REAL*4(A-H,O-Z) C C INITIAL VALUES C FIVE=FLOAT(5) TEN=FLOAT(10) FIFT4=FLOAT(54) SIXTY=FLOAT(60) C C IF A 1X1 DEG MGA CELL IS SPECIFIED C SET THE LAT. AND LONG TO 60 MIN. C IF(NO.EQ.7701)THEN XKPHI=SIXTY XKALA=SIXTY C C ELSE USE THE 5X5(10) MGA CELL SIZE C C IF POINT IS BELOW 54 DEG. SET THE C LAT. AND LONG. TO 5 MIN. C C ELSE IF POINT IS ABOVE 54 DEG. SET C THE LAT. TO 5 MIN. AND LONG. TO 10 MIN. C C ELSE IF(PHIO.LE.FIFT4)THEN XKPHI=FIVE XKALA=FIVE ELSE XKPHI=FIVE XKALA=TEN ENDIF ENDIF RETURN END C C________________________________________________________________________ C------------------------------------------------------------------------ C SRB#4 C ----- C SUBROUTINE CHANG5(PC,DC,M,PSW,DSW,PMAX,DMAX,DG,SDG,H,TOP,IFFGG) C C CHANG5 SETS A FLAG TO CHECK IF THE EVALUATION POINT IS WITHIN C THE 5X5(10) MGA COVERAGE; AND FINDS THE CORRECT 8X8 DEG. DATA BLOCK C BEFORE IT COMPUTES MAX. AND MIN. COORDINATES OF AND READS DATA C FROM THIS BLOCK. C C 5X5(10) MGA FILE DESCRIPTION: C C THE 5X5 MGA SUBFILE IS DIVIDED INTO 75 8X8 DEG. BLOCKS; EACH BLOCK C HAS 9216 RECORDS AND A FOUR DEGREE OVERLAP; EACH RECORD HAS FOUR C FIELDS: DG, SDG,H, AND TOP. C C THE 5X10 MGA SUBFILE IS DIVIDED INTO 125 8X8 DEG. BLOCKS; EACH C BLOCK HAS 4608 RECORDS AND A FOUR DEGREE OVERLAP; EACH RECORD HAS C FOUR FIELDS: DG, SDG,H, AND TOP. C C INPUT: C """"" C C FROM DATA OR AREA: C """"""""""""""""" C PC:LATITUDE OF THE EVALUATION POINT. C DC:LONGITUDE OF THE EVALUATION POINT. C C OUTPUT: C """"""" M: THE NUMBER OF CELLS IN THE DATA 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 PGAS FOR THE 8X8 DEGREE BLOCK. C SDG:ARRAY CONTAINING THE RMS VALUES OF THE 5X5 OR 5X10 C MINITE MEAN PGAS FOR THE 8X8 DEGREE BLOCK. C H:ARRAY CONTAINING THE ASSOCIATED HEIGHTS OF THE C 5X5 OR THE 5X10 MINUTE MEAN PGAS C FOR THE 8X8 DEGREE BLOCK. C TOP:CONTAINS THE VECTOR FOR THE TOPOGRAPHIC CORRECTION. C IMPLICIT REAL*4 (A-H,O-Z) REAL*4 DG(9216) REAL*4 SDG(9216) REAL*4 H(9216) REAL*4 TOP(9216) DIMENSION D1(48),D2(48),D3(48),D4(48) INTEGER*4 REC INTEGER*2 I C C OPENS THE MGA FILE C OPEN(1,STATUS='OLD',ACCESS='DIRECT',FORM='UNFORMATTED',RECL=768) IFFGG=0 C C DETERMINES THE FLAG IFFGG FOR FIVE POSSIBLE CONDITIONS: C C IFFGG=0 POINT LIES WITHIN THE BOUNDARIES OF THE 5X5 MINUTE C MEAN GRAVITY ANOMALY CELL. C IFFGG=1 POINT LIES OUTSIDE THE LATITUDE BOUNDARIES OF THE C 5X5 MINUTE MEAN GRAVITY ANOMALY CELL. C IFFGG=2 POINT LIES OUTSIDE THE LONGITUDINAL BOUNDARIES OF C THE 5X5 MINUTE MEAN GRAVITY ANOMALY CELL. C IFFGG=3 POINT LIES OUTSIDE THE LATITUDE BOUNDARIES OF THE C 5X10 MINUTE MEAN GRAVITY ANOMALY CELL. C IFFGG=4 POINT LIES OUTSIDE THE LONGITUDINAL BOUNDARIES OF C THE 5X10 MINUTE MEAN GRAVITY ANOMALY CELL. C 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 C C INITIALIZES RECORD FIELDS C DO 10 L=1,9216 DG(L)=0.0 SDG(L)=0.0 H(L)=0.0 TOP(L)=0.0 10 CONTINUE C C DETERMINES APPROPIATE 8X8 DEGREE BLOCK WITHIN THE MGA FILE C 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 C C COMPUTES THE MAX. AND MIN. COORDINATES OF THE C 8X8 DEG. BLOCK. C PSW=40.D0+FLOAT((I-1)*4) DSW=214.D0+FLOAT((J-1)*4) PMAX=PSW+8.0 DMAX=DSW+8.0 C C DETERMINES THE POSITION OF THE RECORD WITHIN THE BLOCK C 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 C C READING IN AND REDEFINING THE RECORD FIELDS INTO C THE PROGRAM NOTATION FOR THE 5X5 MGA DATA C DO 20 I=1,192 NO=LLL+I-1 READ(1,REC=NO) (D1(J),D2(J),D3(J),,D4(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) TOP(MN)=D4(JN) 21 CONTINUE 20 CONTINUE GO TO 50 C C READING IN AND REDEFINING THE RECORD FIELDS INTO C THE PROGRAM NOTATION FOR THE 5X10 MGA DATA C 30 DO 40 I=1,96 NO=LLL+I-1 READ(1,REC=NO) (D1(J),D2(J),D3(J),D4(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) TOP(MN)=D4(JN) 31 CONTINUE 40 CONTINUE 50 RETURN END C C________________________________________________________________________ C------------------------------------------------------------------------ C SRB#5 C """""" C SUBROUTINE CHANGP(PC,DC,K,PMAX,DMAX,DG,SDG,H,P,D,IFK) C C CHANGP SETS A FLAG TO CHECK IF THE EVALUATION POINT IS WITHIN C THE PGA COVERAGE; FINDS THE APPROPRIATE PGA SUBFILE BEFORE IT C COMPUTES MAX. AND MIN. COORDINATES OF AND READS DATA FROM THE C CORRECT 10 MIN. X 20 DEG. DATA STRIP. C C PGA FILE DESCRIPTION C C THE PGA FILE CONSISTS OF 20 SUBFILES(PGA1 TO PGA20): 16 10X20 C DEG. SUBFILES AND 4 10X22 DEG. SUBFILES; EACH SUBFILE IS C EQUIVALENT TO ONE BLOCK; EACH BLOCK HAS AN OVERLAP OF 10 MIN., C 61 10 MIN. X 20(22) DEG. DATA STRIPS AND SPACE FOR 9000 RECORDS; C AND EACH RECORD HAS FIVE FIELDS: DG, SDG, H, P, AND D. C C INPUT: C """"" C FROM DATA OR AREA: C """"""""""""""""" C DC: THE LATITUDE OF THE EVALUATION POINT C PC: THE LONGITUDE OF THE EVALUATION POINT C C OUTPUT: C """"""" K: THE NUMBER OF PGAS IN THE STRIP C PMAX: THE MAXIMUM LATITUDE BOUNDARY OF THE STRIP C DMAX: THE MAXIMUM LONGITUDE BOUNDARY OF THE STRIP C DG: ARRAY OF PGAS IN THE STRIP C SDG: ARRAY OF CORRESPNDING PGA STD. DEVS. C H: ARRAY OFCORRESPONDING PGA HEIGHTS(METERS) C P: ARRAY OF CORRESPONDING PGA LATITUDES(DEGREES) C D: ARRAY OF CORRESPONDING PGA LONGITUDES(DEGREES) C IFK: A FLAG USED TO CHECK IF THE POINT IS WITHIN THE C PGA COVERAGE C IMPLICIT REAL*4 (A-H,O-Z) DIMENSION DG(9000),SDG(9000),H(9000),P(9000),D(9000) C C DETERMINING CONDITIONS FOR FLAG IFK: C C IFK=0 :POINT LIES INSIDE THE SPECIFIED BOUNDARIES C IFK=1 :POINT LIES OUTSIDE OF LATITUDE BOUNDARY OF C POINT PGAS. C IFK=2 :POINT LIES OUTSIDE OF LONGITUDE BOUNDARY OF C POINT PGAS. C IFK=0 IF(PC.LE.40.0833.OR.PC.GE.79.9167) THEN IFK=1 GO TO 100 ELSE END IF IF(DC.LE.218.0833.OR.DC.GE.319.9167) THEN IFK=2 GO TO 100 ELSE END IF C C INITIALIZING MGA RECORD FIELDS C DO 20 L=1,9000 DG(L)=0.0 SDG(L)=0.0 H(L)=0.0 P(L)=0.0 D(L)=0.0 20 CONTINUE C C DETERMINING WHICH 10X20 DEG DATA BLOCK TO USE C I=(ABS(PC-45.0)/10.0)+1.0 J=(ABS(DC-228.0)/20.0)+1.0 KCC=45+10*(I-1) LCC=228+20*(J-1) TP=PC-FLOAT(KCC) TD=DC-FLOAT(LCC) C C CHECKS FOR THE 10 MIN. OVERLAP AND THE 10X22 DEG DATA BLOCK C IF(TP.GT.5.041667) THEN I=I+1 KCC=KCC+10 ELSE END IF IF(TD.GT.10.041667)THEN J=J+1 LCC=LCC+20 ELSE END IF IF(J.EQ.5)THEN DMAX=FLOAT(LCC)+12.0 ELSE DMAX=FLOAT(LCC)+10.0 END IF C C FINDS POSITION OF BLOCK WITHIN PGA FILE C LMIN=LCC-10 LF=(I-1)*5+J JK=LF+9 C C FINDS POSITION OF 10 MIN.X 20 DEG. STRIP WITHIN DATA BLOCK C K=0 IS=HFIX(PC) BSC=PC-FLOAT(IS) IMU=HFIX(BSC*60.0) IMI=IMU/5 LC=5*IMI IDD=IMU-LC IF(IDD.LE.2.5) THEN MUN=LC ELSE MUN=LC+5 END IF C C COMPUTES MAX AND MIN COORDINATES OF DATA STRIP C PHS=FLOAT(IS)+FLOAT(MUN)/60.0 PMAX=PHS+0.083333 PMIN=PHS-0.083333 C C READING IN RECORD FIELDS C 50 READ(LF+9,88,END=100)P1,D1,ZG,ZDG,ZH,ZDH 88 FORMAT(F13.5,2X,F13.5,F9.1,F8.3,F10.2,F8.3) C C CHECKING THE VALIDITY OF THE DATA C IF(P1.GT.PMAX.OR.P1.LT.PMIN) THEN GO TO 50 ELSE END IF IF(K.LE.1)THEN GO TO 99 ELSE END IF IF(P1.EQ.P(K).AND.D1.EQ.D(K)) THEN GO TO 50 ELSE END IF C C REDEFINING THE FIELDS C 99 K=K+1 P(K)=P1 D(K)=D1 DG(K)=ZG SDG(K)=ZDG H(K)=ZH C GO TO 50 C 100 REWIND (LF+9) RETURN END C C________________________________________________________________________ C------------------------------------------------------------------------ C SRB#6 C """""" C SUBROUTINE INMBND(NK,PHIO,ALAO,PHICS,PHICN,ALACW,ALACE) C C INMBND DETERMINES THE BOUNDARIES OF THE INNERMOST ZONE IN MINUTES C FOR THE EVALUATION POINT GIVEN IN DECIMAL DEGREES. C C INPUT: C """"" C FROM MAIN PROGRAM C """"""""""""""""" C NK :1/2 THE INNERMOST ZONE SIZE IN MINUTES C C FROM DATA OR AREA: C """"""""""""""""" C PHIO: THE LATITUDE OF THE EVALUATION POINT IN DEGREES C ALAO: THE LONGITUDE OF THE EVALUATION POINT IN DEGREES C C OUTPUT: C """"""" PHICS: THE SOUTHERN BORDER OF THE INNERMOST ZONE C PHICN: THE NORTHERN BORDER OF THE INNERMOST ZONE C ALACW: THE WESTERN BORDER OF THE INNERMOST ZONE C ALACE: THE EASTERN BORDER OF THE INNERMOST ZONE C IMPLICIT REAL*4(A-H,O-Z) XN=FLOAT(NK) XM=FLOAT(60) XPHIO=ROUND(PHIO) XALAO=ROUND(ALAO) PHICS=(XPHIO-XN)/XM PHICN=(XPHIO+XN)/XM ALACW=(XALAO-XN)/XM ALACE=(XALAO+XN)/XM RETURN END C C_______________________________________________________________________ C------------------------------------------------------------------------ C SRB#7 C """""" SUBROUTINE BOUND(PHIO,ALAO,PHICS,PHICN,ALACW,ALACE,X1,X2,Y1,Y2) C C BOUND COMPUTES THE INTEGRATION LIMITS FOR THE INNERMOST IN RADIANS C GIVEN THE EVALUATION POINT IN DEG. AND THE INNERMOST ZONE BORDERS C IN MINUTES. C C INPUT: C """"" C FROM DATA OR AREA C """"""""""""""""" C PHIO:IS THE LATITUDE OF THE EVALUATION POINT IN DEGREES C ALAO:IS THE LONGITUDE OF THE EVALUATION POINT IN DEGREES C C FROM INMBND: C """"""""""" C PHICS:IS THE SOUTHERN BORDER OF THE INNERMOST ZONE C PHICN:IS THE NORTHERN BORDER OF THE INNERMOST ZONE C ALACW:IS THE WESTERN BORDER OF THE INNERMOST ZONE C ALACE:IS THE EASTERN BORDER OF THE INNERMOST ZONE C C OUTPUT: C """"""" X1:IS THE LOWER LIMIT IN THE X-DIRECTION C X2:IS THE UPPER LIMIT IN THE X-DIRECTION C Y1:IS THE LOWER LIMIT IN THE Y-DIRECTION C Y2:IS THE UPPER LIMIT IN THE Y-DIRECTION C IMPLICIT REAL*4(A-H,O-Z) PI=DARCOS(-1.D0) DR=PI/180.0 XLAT=PHIO*DR X1=COS(XLAT)*(ALACW-ALAO)*DR X2=COS(XLAT)*(ALACE-ALAO)*DR Y1=(PHICS-PHIO)*DR Y2=(PHICN-PHIO)*DR RETURN END C C________________________________________________________________________ C------------------------------------------------------------------------ C SRB#8 C """""" SUBROUTINE NUMERC(X1,X2,Y1,Y2,C) C C NUMERC COMPUTES THE VECTOR 'C' USED IN THE NUMERICAL INTEGRATION C OF THE INNERMOST ZONE GIVEN THE INTEGRATION LIMITS IN RADIANS. C C INPUT: C """"" C FROM BOUND C """""""""" C X1: THE LOWER LIMIT IN THE X-DIRECTION C X2: THE UPPER LIMIT IN THE X-DIRECTION C Y1: THE LOWER LIMIT IN THE Y-DIRECTION C Y2: THE UPPER LIMIT IN THE Y-DIRECTION C C OUTPUT: C """"""" C: A VECTOR CONTAINING THE RESULTS OF THE SURFACE C INTEGRATION (REF."THE CANADIAN GEOID",FEB.1986, C SECTION 2.2.3.1 THE INNERMOST ZONE INTEGRATION C PAGE 22.) C IMPLICIT REAL*4(A-H,O-Z) DIMENSION C(1,6) C TT= 0.6378137/1.2312529 C C TO INITIALIZE THE NUMERICAL INTEGRATION VECTOR 'C' C -------------------------------------------------- C DO 10 I=1,6 C(1,I)=0.0 10 CONTINUE C C C COEFF OF A0 C ------------- C R1=SQRT(X2**2+Y2**2) R2=SQRT(X1**2+Y2**2) R3=SQRT(X1**2+Y1**2) R4=SQRT(X2**2+Y1**2) R5=ABS((Y2+R1)/(Y1+R4)) R6=ABS((Y2+R2)/(Y1+R3)) R7=ABS((X2+R1)/(X1+R2)) R8=ABS((X2+R4)/(X1+R3)) R9=X2*LOG(R5) R10=X1*LOG(R6) R11=Y2*LOG(R7) R12=Y1*LOG(R8) XI1=R9-R10+R11-R12 C C COEF OF B0 C ------------ C T1=SQRT(R3/2.0) T2=SQRT(R2/2.0) T3=SQRT(R4/2.0) T4=SQRT(R1/2.0) T5=X1*Y1*(-3.0/2.0+LOG(T1)) T6=X2*Y2*(-3.0/2.0+LOG(T4)) T7=ABS(X1*Y2)*(-3.0/2.0+LOG(T2)) T8=ABS(X2*Y1)*(-3.0/2.0+LOG(T3)) T9=T5+T6+T7+T8 T10=X1**2*ARCOS(X1/R3) T11=X2**2*ARCOS(X2/R1) T12=R3**2*ARSIN(X1/R3)/2.0 T13=R1**2*ARSIN(X2/R1)/2.0 T14=T10+T11+T12+T13 T15=X1**2*ARCOS(ABS(X1)/R2) T16=X2**2*ARCOS(ABS(X2)/R4) T17=R2**2*ARCOS(ABS(Y2)/R2)/2.0 T18=R4**2*ARCOS(ABS(Y1)/R4) T19=T15+T16+T17+T18 XI7=T9+T14+T19 C C COEF OF CO C ---------- C XI8=(X2-X1)*(Y2-Y1) C C COEF OF A1 C ---------- C C1=Y2*(R1-R2) C2=Y1*(R4-R3) C3=X2*R9 C4=X1*R10 XI2=(C1-C2+C3-C4)/2.0 C C C C C COEF OF A2 C ---------- C B1=X2*(R1-R4) B2=X1*(R2-R3) B3=Y2*R11 B4=Y1*R12 XI3=(B1-B2+B3-B4)/2.0 C C COEF OF A3 C ---------- C D1=R1**3 D2=R4**3 D3=R2**3 D4=R3**3 XI4=(D1-D2-D3+D4)/3.0 C C COEF 0F A4 C ------------ C E1=X2*Y2*R1 E2=X2*Y1*R4 E3=X1*Y2*R2 E4=X1*Y1*R3 E5=(E1-E2-E3+E4)/6.0 E6=(X2*C3-X1*C4)/3.0 E7=(Y2*B3-Y1*B4)/6.0 XI5=E5+E6-E7 C C C COEF OF A5 C ----------- C F1=Y2*X2*R1 F2=Y2*X1*R2 F3=Y1*X2*R4 F4=Y1*X1*R3 F5=(F1-F2-F3+F4)/6.0 F6=(Y2*B3-Y1*B4)/3.0 F7=(X2*C3-X1*C4)/6.0 XI6=F5+F6-F7 C C C C(1,1)= TT*(2.0*XI1-3.0*XI7-4.0*XI8) C(1,2)=2.0*TT*XI2 C(1,3)=2.0*TT*XI3 C(1,4)=2.0*TT*XI4 C(1,5)=2.0*TT*XI5 C(1,6)=2.0*TT*XI6 C RETURN END C________________________________________________________________________ C------------------------------------------------------------------------ C SRB#9 C """""" SUBROUTINE INBND(NJ,PHIO,ALAO,PHIIS,PHIIN,ALAIW,ALAIE) C C INBND FORMS THE BOUNDARIES OF THE INNER ZONE GIVEN THE EVALUATION C POINT IN DECIMAL DEGREES. C C INPUT: C """"" FROM THE MAIN PROGRAM C """"""""""""""""""""" C NJ: 1/2 THE INNER ZONE SIZE IN DEG. C C FROM DATA OR AREA: C """"""""""""""""" C PHIO: LATITUDE OF THE EVALUATION POINT IN DEG. C ALAO: THE LONGITUDE OF THE EVALUATION POINT IN DEG. C C OUTPUT: C """"""" PHIIS: SOUTHERN BORDER OF THE INNER ZONE IN DEG. C PHIIN: NORTHERN BORDER OF THE INNER ZONE IN DEG. C ALAIW: WESTERN BORDER OF THE INNER ZONE IN DEG. C ALAIE: EASTERN BORDER OF THE INNER ZONE IN DEG. C IMPLICIT REAL*4(A-H,O-Z) C C SEPARATING POINT INTO DEGREES AND THE DECIMAL PORTION C OF THE DEGREE C XNJ=FLOAT(NJ) IPHI=PHIO IALA=ALAO XIPHI=IPHI XIALA=IALA RP=PHIO-XIPHI RL=ALAO-XIALA C C FORMING BORDERS OF THE INNER ZONE C 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#10 C """"""" SUBROUTINE CELNUM(ILAT,ILONG,XKPHI,XKALA,XNPHI,XNALA) C C CELNUM COMPUTES THE NUMBER OF CELLS IN THE LAT. AND LONG. C DIRECTIONS OF THE DESIRED ZONE. C C INPUT: C """"" FROM MAIN PROGRAM C """"""""""""""""" C ILAT: INTERVAL IN SOUTH TO NORTH DIRECTION IN MIN. C ILONG: INTERVAL IN THE WEST TO EAST DIRECTION IN MIN. C C FROM CELSZE C """"""""""" C XKPHI: CELLSIZE IN LAT. DIRECTION IN MIN. C XKALA: CELLSIZE IN LONG. DIRECTION IN MIN. C C OUTPUT: C """"""" XNPHI: THE NUMBER OF CELLS IN THE LATITUDINAL DIRECTION C XNALA: 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#11 C ------- SUBROUTINE CHECK(PHOUTS,PHOUTN,ALOUTW,ALOUTE,PHNERS,PHNERN, * ALNERW,ALNERE,PHIO,ALAO,XKPHI,XKALA,PHISW,ALASW, * IFLG,SPHOUT,XNPHOT,WALOUT,EALOUT,SPHNER,XNPHNR,WALNER, * EALNER,XPHISW,XALASW) C C C CHECK CONVERTS THE BORDERS OF TWO ZONES FROM DEGREES TO C MINUTES AND THEN SETS A FLAG FOR THREE CONDITIONS THAT C WILL BE MENTIONED LATER. INNER REFERS TO THE ZONE INSIDE THE C OUTER ZONE. C C INPUT: C """""" C FROM INBND OR OUTER C """"""""""""""""""" C PHOUTS: SOUTHERN BORDER OF THE OUTER ZONE C PHOUTN: NORTHERN BORDER OF THE OUTER ZONE C ALOUTW: WESTERN BORDER OF THE OUTER ZONE C ALOUTE: EASTERN BORDER OF THE OUTER ZONE C C FROM INMBND OR INBND C """""""""""""""""""" C PHNERS: SOUTHERN BORDER OF THE INNER ZONE C PHNERN: NORTHERN BORDER OF THE INNER ZONE C ALNERW: WESTERN BORDER OF THE INNER ZONE C ALNERE: EASTERN BORDER OF THE INNER ZONE C C FROM DATA OR AREA C """"""""""""""""" C PHIO: LATITUDE OF THE EVALUATION POINT C ALAO: LONGITUDE OF THE EVALUATION POINT C C FROM CELSZE C """"""""""" C XKPHI: CELL LENGTH IN THE PHI DIRECTION C XKALA: CELL LENGTH IN THE LAMDA DIRECTION C C FROM CHANG5 C """"""""""" C PHISW: SOUTHERN LATITUDE OF THE MGA DATA BLOCK C ALASW: WESTERN LONGITUDE OF THE MGA DATA BLOCK C C NOTE: C """"" THE INPUT PARAMETERS ARE ALL IN DEGREES C C OUTPUT: C """"""" IFLG:IS AN ERROR FLAG WITH FOUR POSSIBLE VALUES C ----------------------------------------------- C #1) IF IFLG=0 EVERYTHING IS OK C #2) IF IFLG=-1 THEN THE DIFFERENCES BETWEEN BOUNDARIES C ARE NOT INTEGER MULTIPLES OF CELL SIZE C #3) IF IFLG=-2 THEN THE BOUNDARIES OF THE INNERMOST ZONE C IS OUTSIDE OF THE INNER ZONE CELL C #4) IF IFLG=-3 THEN THE EVALUATION POINT LIES OUTSIDE C OF THE INNERMOST ZONE CELL C C OUTER ZONE C """"""""""" C SPHOUT: SOUTHERN BORDER OF THE OUTER ZONE C XNPHOT: NORTHERN BORDER OF THE OUTER ZONE C WALOUT: WESTERN BORDER OF THE OUTER ZONE C EALOUT: EASTERN BORDER OF THE OUTER ZONE C C INNER ZONE C """""""""" C SPHNER: SOUTHERN BORDER OF THE INNER ZONE C XNPHNR: NORTHERN BORDER OF THE INNER ZONE C WALNER: WESTERN BORDER OF THE INNER ZONE C EALNER: EASTERN BORDER OF THE INNER ZONE C C DATA STRIP C """"""""""" C XPHISW: SOUTHERN LATITUDE OF THE DATA STRIP C XALASW: WESTERN LONGITUDE OF THE DATA STRIP C IMPLICIT REAL*4(A-H,O-Z) C KPHI=XKPHI KALA=XKALA C C --------------------------------------- C TO CONVERT THE BOUNDARIES INTO MINUTES BY CALLING C FUNC#2 -DEGMIN- WHICH CONVERTS DEGREES TO MINUTES C --------------------------------------- C XPHISW=DEGMIN(PHISW) XALASW=DEGMIN(ALASW) SPHOUT=DEGMIN(PHOUTS) XNPHOT=DEGMIN(PHOUTN) WALOUT=DEGMIN(ALOUTW) EALOUT=DEGMIN(ALOUTE) SPHNER=DEGMIN(PHNERS) XNPHNR=DEGMIN(PHNERN) WALNER=DEGMIN(ALNERW) EALNER=DEGMIN(ALNERE) C C ----------------------------------------------------- C TO FIND DIFFERNCES BETWEEN OUTER AND INNER BOUNDARIES C ----------------------------------------------------- C IDIF1=WALNER-WALOUT IDIF2=EALOUT-EALNER IDIF3=SPHNER-SPHOUT IDIF4=XNPHOT-XNPHNR C C -------------------------------------------------- C TO CHECK WHETHER THE DIFFERNCES BETWEEN BOUNDARIES C ARE INTEGER MULTIPLES OF CELL SIZE C ------------------------------------------------- C IF(MOD(IDIF1,KALA).NE.0.OR. * MOD(IDIF2,KALA).NE.0.OR. * MOD(IDIF3,KPHI).NE.0.OR. * MOD(IDIF4,KPHI).NE.0)THEN IFLG=-1 ELSE C C --------------------------------------------------------- C TO CHECK WHETHER INNER BOUNDARIES WITHIN OUTER BOUNDARIES C --------------------------------------------------------- C IF(PHNERS.LE.PHOUTS.OR. * PHNERN.LE.PHNERS.OR. * PHOUTN.LE.PHNERN.OR. * ALNERW.LE.ALOUTW.OR. * ALNERE.LE.ALNERW.OR. * ALOUTE.LE.ALNERE)THEN IFLG=-2 ELSE C C ---------------------------------------------- C TO CHECK WHETHER POINT OF COMPUTATION WITHIN C INNER ZONE C ------------------------------------------- C IF(PHIO.LE.PHNERS.OR. * PHNERN.LE.PHIO.OR. * ALAO.LE.ALNERW.OR. * ALNERE.LE.ALAO) THEN IFLG=-3 ELSE IFLG=0 ENDIF ENDIF ENDIF RETURN END C C________________________________________________________________________ C------------------------------------------------------------------------ C SRB#12 C """""" SUBROUTINE SCAN(K,M,XPHISW,XALASW,DGARY,SDGARY,HARAY,TOPAY, * XKPHI,XKALA,XNALA,SPHOUT,WALOUT,XNPHOT,EALOUT, * SPHNER,XNPHNR,WALNER,EALNER,RETPHI,RETALA,RETDG,RETSDG, * RETHRY,TTOP,NUMB,ISIGA,ISIGB,PHIA,ALI,GED,SGED,HED,TED,NUM2) C C SCAN FINDS THE CELLS OF THE INNERMOST AND THE INNER ZONE; OR THE C INNER AND THE OUTER ZONE; AND SEPARATES THEM INTO TWO SEPARATE C ARRAYS SCAN ALSO SETS A FLAG TO IF THE INNERMOST, INNER, OR OUTER C ZONES HAVE NO 5X5(10) MIN MGAS. C C INPUT: C"""""" C FROM MAIN PROGRAM C """"""""""""""""" C K: NUMBER OF RECORDS IN THE DATA BLOCK C M: NUMBER OF CELLS IN INNER OR OUTER ZONE C C FROM SUBROUTINE CHECK C """""""""""""""""""""" C XPHISW: SOUTHERN LAT. OF THE DATA STRIP IN MINUTES C XALASW: WESTERN LONG. OF THE DATA STRIP IN MINUTES C C FROM SUBROUTINE CHANG5 C """""""""""""""""""""" C DGARY: AN ARRAY CONTAINING THE 5X5(10) MGAS FOR THE DATA BLOCK C SDGARY: CORRESPONDING RMS OF THE MGAS C HARAY: CORRESPONDING HEIGHTS. C C FROM SUBROUTINE CELSZE C """""""""""""""""""""" C XKPHI: CELL SIZE IN LATITUDINAL DIRECTION C XKALA: CELL SIZE IN LONGITUDINAL DIRECTION C C FROM SUBROUTINE CELNUM C """""""""""""""""""""" C XNALA: NUMBER OF CELLS IN THE LONGITUDINAL DIRECTION C C FROM SUBROUTINE CHECK C """""""""""""""""""""" C SPHOUT: SOUTHERN BORDER OF THE INNER OR OUTER ZONE C WALOUT: WESTERN BORDER OF THE INNER OR OUTER ZONE C XNPHOT: NORTHERN BORDER OF THE INNER OR OUTER ZONE C EALOUT: EASTERN BORDER OF THE INNER OR OUTER ZONE C SPHNER: SOUTHERN BORDER OF INNERMOST OR INNER ZONE C XNPHNR: NORTHERN BORDER OF INNERMOST OR INNER ZONE C WALNER: WESTERN BORDER OF INNERMOST OR INNER ZONE C EALNER: EASTERN BORDER OF INNERMOST OR INNER ZONE C C OUTPUT: C """"""" FOR INNER OR OUTER ZONE C """"""""""""""""""""""" C RETPHI: AN ARRAY OF THE LATITUDES IN THE INNER OR OUTER ZONE C RETALA: AN ARRAY OF THE LONGITUDES IN THE INNER OR OUTER ZONE C RETDG: CORRESPONDING ARRAY OF MGAS C RETSDG: CORRESPONDING ARRAY OF RMS OF THE MGAS C RETHRY: CORRESPONDING ARRAY OF HEIGHTS C NUMB: NUMBER OF CELLS IN THE INNER OR OUTER ZONE C C ISIGB: 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 OR OUTER ZONE. C IF ISIGB=-5 THERE ARE EMPTY CELL OF 5'X 5'OR 5'X 10' C WITHIN THE INNER OR OUTER ZONE. C C INNERMOST OR INNER ZONE C """"""""""""""""""""""" C PHIA: AN ARRAY OF THE LATITUDES IN THE INNERMOST OR INNER ZONE C ALI: CORRESPONDING ARRAY OF LONGITUDES C GED: CORRESPONDING ARRAY OF MGAS C SGED: CORRESPONDING ARRAY THE RMS OF THE MGAS C HED: CORRESPONDING ARRAY OF THE HEIGHTS C NUM2: NUMBER OF MGAS C C ISIGA: AN ERROR FLAG WITH TWO POSSIBLE VALUES C IF ISIGA=0 THERE ARE NO EMPTY CELLS OF 5'X 5'OR 5'X 10' C WITHIN THE INNNERMOST OR INNER ZONE BOUNDARIES. C IF ISIGA=-4 EMPTY CELL OF 5'X 5' OR 5'X 10' IS C LOCATED INSIDE THE INNERMOST OR INNER ZONE BOUNDARIES. C IMPLICIT REAL*4(A-H,O-Z) DIMENSION DGARY(K),SDGARY(K),HARAY(K),RETPHI(M),TOPAY(M), * RETALA(M),RETDG(M),RETSDG(M),RETHRY(M),TTOP(M), * PHIA(4),ALI(4),GED(4),SGED(4),HED(4),TED(4) LOGICAL INSIDE C C INITIAL VALUES C ISIGA=0 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 PHII=SPHOUT+XKPHI*HALF ALAI=WALOUT+XKALA*HALF NUMB=0 NUM2=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 TO CHECK WHETHER PHII & ALAI ARE ARE WITHIN C INNER BOUNDARIES(IE.INNERMOST OR INNER C DEPENDING ON THE CALLING SUBROUTINE) C INSIDE=.FALSE. IF(PHII.GE.SPHNER.AND.PHII.LE.XNPHNR)THEN IF(ALAI2.GE.WALNER.AND.ALAI2.LT.EALNER)THEN INSIDE=.TRUE. END IF END IF C C DEFINING INNERMOST OR INNER ZONE CELLS C IF(INSIDE)THEN NUM2=NUM2+1 GED(NUM2)=DGARY(IPOS) IF(GED(NUM2).EQ.9999.0) THEN ISIGA=-4 C RETURN ENDIF FPHII=PHII/XM FALAI=ALAI2/XM PHIA(NUM2)=FPHII ALI(NUM2)=FALAI SGED(NUM2)=SDGARY(IPOS) HED(NUM2)=HARAY(IPOS) TED(NUM2)=TOPAY(IPOS) ENDIF C C DEFINING INNER OR OUTER ZONE CELLS C IF(.NOT. INSIDE)THEN 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) TTOP(NUMB)=TOPAY(IPOS) ENDIF C C TO SHIFT TO THE NEXT BLOCK C ALAI2=ALAI2+XKALA GOTO 565 300 PHII=PHII+XKPHI GOTO 564 200 CONTINUE C RETURN END C C________________________________________________________________________ C------------------------------------------------------------------------ C SRB#13 C """"""" C SUBROUTINE OUTER(PSIO,PHIO,ALAO,PHIOS,PHION,ALAOW,ALAOE) C C ----------- C OUTER IS USED TO FORM THE BOUNDARIES OF THE OUTER ZONE C GIVEN THE EVALUATION POINT AND THE CAP SIZE. C C INPUT: C """"" FROM MAIN PROGRAM C """"""""""""""""" C PSIO: RADIUS OF THE CAP OF INTEGRATION IN RADIANS C C FROM DATA AND AREA: C """""""""""""""""" C PHIO: LATITUDE OF THE EVALUATION POINT C ALAO: LONGITUDE OF THE EVALUATION POINT C C OUTPUT: C """""""PHIOS: SOUTH BORDER OF THE OUTER ZONE C PHION: NORTH BORDER OF THE OUTER ZONE C ALAOW: WEST BORDER OF THE OUTER ZONE C ALAOE: EAST BORDER OF THE OUTER ZONE C NOTE: C """"" ALL PARAMETERS ARE OUTPUTED IN DEGREES C IMPLICIT REAL*4(A-H,O-Z) C C CONVERTING POINT AND CAP SIZE INTO RADIANS C PI=DARCOS(-1.D0) DR=PI/180.0 XPHIO=PHIO*DR SCOS=COS(XPHIO) IDLAM=PSIO/SCOS DLAM=IDLAM+1 C C SEPARATING POINTS INTO DEGREES AND THE DECIMAL PORTION C OF THE DEGREE C IPHI=PHIO IALA=ALAO XIPHI=IPHI XIALA=IALA RP=PHIO-XIPHI RL=ALAO-XIALA C C FORMING OUTER BOUNDARIES C 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 PHIOS=PHI-PSIO PHION=PHI+PSIO ALAOW=ALA-DLAM ALAOE=ALA+DLAM RETURN END C C________________________________________________________________________ C------------------------------------------------------------------------ C SRB#14 C """"""" SUBROUTINE FIND(INER,PHIO,ALAO,PHIARY,ALAARY,DG,SDG,HETARY, * PHICS,PHICN,ALACW,ALACE,NUM,PHARY,XLARAY,FARAY,SARAY,REHARY) C C FIND: LOCATES THE PGA FILE AND EXTRACTS THE PGAS FOR THE C INNERMOST ZONE. C C INPUT: C """""" C FROM SUBROUTINE DATA OR AREA C """""""""""""""""""""""""""" C PHIO: LATITUDE OF THE EVALUATION POINT C ALAO: LONGITUDE OF THE EVALUATION POINT C C FROM SUBROUTINE CHANGP C """""""""""""""""""""" C INER: NUMBER OF PGAS IN THE DATA STRIP C PHIARY: AN ARRAY OF LATITUDE PGAS IN DATA STRIP C ALAARY: CORRESPONDING ARRAY OF LONGITUDES C DG: CORRESPONDING ARRAY OF PGAS C SDG: CORRESPONDING ARRAY OF THE RMS OF THE ANOMALIES C HETARY: CORRESPONDING ARRAY OF HEIGHTS C C FROM SUBROUTINE INMBD C """"""""""""""""""""" C PHICS: SOUTHERN BORDER OF THE INNERMOST ZONE C PHICN: NORTHERN BORDER OF THE INNERMOST ZONE C ALACW: WESTERN BORDER OF THE INNERMOST ZONE C ALACE: EASTERN BORDER OF THE INNERMOST ZONE C C OUTPUT: C """"""" NUM: NUMBER OF PGAS IN THE INNERMOST ZONE C PHARY: AN ARRAY OF LATITUDES OF THE PGAS IN THE INNERMOST ZONE C XLARAY: CORRESPONDING ARRAY OF LONGITUDES C FARAY: CORRESPONDING ARRAY OF PGAS C SARAY: CORRESPONDING ARRAY OF THE RMS OF THE PGAS C REHARY: CORRESPONDING ARRAY OF HEIGHTS C IMPLICIT REAL*4(A-H,O-Z) DIMENSION DG(INER),SDG(INER),PHIARY(INER),ALAARY(INER), * HETARY(INER) DIMENSION PHARY(INER),XLARAY(INER),FARAY(INER),SARAY(INER), * REHARY(INER) NUM=0 C DO 10 J=1,INER IF(ALAARY(J).GT.ALACW.AND.ALAARY(J).LT.ALACE)THEN NUM=NUM+1 PHARY(NUM)=PHIARY(J) XLARAY(NUM)=ALAARY(J) FARAY(NUM)=DG(J) SARAY(NUM)=SDG(J) REHARY(NUM)=HETARY(J) ENDIF 10 CONTINUE RETURN END C C________________________________________________________________________ C------------------------------------------------------------------------ C SRB#15 C """""" C SUBROUTINE COM(PHIO,ALAO,PHARY,XLARAY,NUM,X,Y,X2,Y2) C C COMPUTES LOCAL CARTESIAN COORDINATES OF THE INNERMOST ZONE C PGAS. C C INPUT: C """"" FROM DATA OR AREA C """"""""""""""""" C PHIO: LATITUDE (IN DEG.) OF THE EVALUATION POINT C ALAO: LONGITUDE (IN DEG.) OF THE EVALUATION POINT C C FROM FIND C """"""""" C PHARY: ARRAY OF LATITUDES OF PGAS WITHIN C THE INNERMOST ZONE(IN DEG.) C XLARAY: CORRESPONDING ARRAY OF LONGITUDES(IN DEG.) C NUM: NUMBER OF PGAS WITHIN THE INNERMOST ZONE C C FROM NUMERC C """"""""""" C X2: UPPER X-LIMIT OF THE INNERMOST ZONE (IN RAD) C Y2: UPPER Y-LIMIT OF THE INNERMOST ZONE (IN RAD) C C OUTPUT: C """"""" X: AN ARRAY OF X-VALUES (IN RAD) OF PGAS IN THE C INNERMOST ZONE C Y: CORRESPONDING ARRAY OF Y-VALUES (IN RAD) C C C C IMPLICIT REAL*4(A-H,O-Z) DIMENSION PHARY(NUM),XLARAY(NUM),X(NUM),Y(NUM) PI=DARCOS(-1.D0) DR=PI/180.0 C R=6731000 XLAT=PHIO*DR DO 40 J=1,NUM PHI=PHARY(J) ALA=XLARAY(J) X(J)=COS(XLAT)*(ALA-ALAO)*DR/X2 Y(J)=(PHI-PHIO)*DR/Y2 40 CONTINUE RETURN END C C________________________________________________________________________ C------------------------------------------------------------------------ C SRB#16 C """"""" SUBROUTINE SURF(NUM,NP,X,Y,PHARY,XLARAY,FARAY,SARAY,REHARY,COV, * ALFA,XDET,PHIO,ALAO,NFLAG) C C SURF IS USED TO CHECK THE GEOMETRY OF THE PGAS C AND PERFORM A LEAST SQUARES FIT TO THE INNERMOST ZONE. C C INPUT: C """"" C FROM FIND C """"""""" C NUM: NUMBER OF PGAS IN THE INNERMOST ZONE C C FROM FUNCTION ICOUNT C """""""""""""""""""" C NP: NUMBER OF PARAMETERS TO BE SOLVED FOR C IN THE LEAST SQUARE ADJUSTMENT C C FROM COM C """"""""" C X: X COORDINATE ARRAY OF PGAS IN INNERMOST ZONE C Y: Y COORDINATE ARRAY OF PGAS IN INNERMOST ZONE C C FROM FIND C """"""""" C PHARY: AN ARRAY OF LATITUDES OF PGAS IN INNERMOST ZONE C XLARAY: THE CORRESPONDING ARRAY OF LONGITUDES C FARAY: THE CORRESPONDING ARRAY OF PGAS C SARAY: THE CORRESPONDING ARRAY OF THE RMS OF THE PGAS C REHARY: THE CORRESPONDING ARRAY OF HEIGHTS C C FROM DATA OR AREA C """"""""""""""""" C ALAO: LONGITUDE OF THE EVALUATION POINT C PHIO: LATITUDE OF THE EVALUATION POINT C C C OUTPUT: C """"""" COV: THE COVARIANCE MATRIX OF THE SURFACE FITTING C SOLUTION C ALFA: THE VECTOR OF CORRECTIONS TO THE SURFACE C FITTED PGAS C XDET: THE DETERMINANT OF COVARIANCE MATRIX C NFLAG:REFERS TO THE GEOMETRY OF THE PGAS BEING C FITTED TO THE SURFACE C IMPLICIT REAL*4(A-H,O-Z) REAL*8 NX(6,6),DET,ATPAD(6,6),ATCAD(6,6),ATPLD(6) REAL*8 CXD(6,6),CXXD(6,6),RD(6) C DIMENSION FARAY(NUM),SARAY(NUM),REHARY(NUM),X(NUM),Y(NUM) DIMENSION PHARY(NUM),XLARAY(NUM),EING(6),A(6,6) DIMENSION ATPA(6,6),ATPL(6),R(6),ALFA(6),COV(6,6),VEC(6), *IW1(6),IW2(6),ATCA(6,6),CX(6,6),CXX(6,6),NP1(6),Z(6,6),WK(6) INTEGER*4 NNP(4) C C TESTING OF THE GEOMETRY OF THE PGAS C """"""""""""""""""""""""""""""""""" NFLAG=0 IF(NP.EQ.1.OR.NP.EQ.3)GO TO 2222 C NNP(1)=0 NNP(2)=0 NNP(3)=0 NNP(4)=0 C DO 222 I=1,NUM NN=0 IF(PHARY(I).LT.PHIO)NN=2 IF(XLARAY(I).LT.ALAO)NN=NN+2 IF(XLARAY(I).GT.ALAO)NN=NN+1 NNP(NN)=NNP(NN)+1 222 CONTINUE C IF (NNP(1).GE.1.AND.NNP(2).GE.1.AND. & NNP(3).GE.1)GO TO 2222 IF (NNP(2).GE.1.AND.NNP(3).GE.1.AND. & NNP(4).GE.1)GO TO 2222 IF (NNP(3).GE.1.AND.NNP(4).GE.1.AND. & NNP(1).GE.1)GO TO 2222 IF (NNP(4).GE.1.AND.NNP(1).GE.1.AND. & NNP(2).GE.1)GO TO 2222 GO TO 2223 C 2222 CONTINUE C C INITIALIZATION OF MATICES C """"""""""""""""""""""""" C DO 29 II=1,6 ALFA(II)=0.0 DO 30 JJ=1,6 COV(II,JJ)=0.0 30 CONTINUE 29 CONTINUE IF(NP.EQ.1)THEN XDET=0.0 C XNUM=NUM XN=0.0 VAR=0.0 DO 52 IJJ=1,NUM XN=XN+FARAY(IJJ) VAR=VAR+SARAY(IJJ)**2 52 CONTINUE ALFA(1)=XN/XNUM COV(1,1)=VAR ELSE C DO 428 LP=1,6 ATPL(LP)=0.0 ATPLD(LP)=0.0D0 DO 429 LD=1,6 A(LP,LD)=0.0 ATPA(LP,LD)=0.0 ATPAD(LP,LD)=0.0D0 ATCA(LP,LD)=0.0 ATCAD(LP,LD)=0.0D0 CXD(LP,LD)=0.0D0 CXX(LP,LD)=0.0 CXXD(LP,LD)=0.0D0 COV(LP,LD)=0.0 429 CONTINUE 428 CONTINUE C DO 135 IK=1,NUM G=FARAY(IK) SG=SARAY(IK) IF(SG.EQ.0.0)THEN SG=1.0 ENDIF P=1.0 C=SG**2 VEC(1)=1.0 VEC(2)=X(IK) VEC(3)=Y(IK) VEC(4)=X(IK)*Y(IK) VEC(5)=X(IK)**2 VEC(6)=Y(IK)**2 C DO 134 IM=1,NP DO 133 JM=1,NP ATPA(IM,JM)=ATPA(IM,JM)+VEC(IM)*VEC(JM) A(IM,JM)=ATPA(IM,JM) ATCA(IM,JM)=ATCA(IM,JM)+VEC(IM)*VEC(JM)*C 133 CONTINUE 134 CONTINUE C C DO 132 KI=1,NP ATPL(KI)=ATPL(KI)+VEC(KI)*G ATPLD(KI)=DBLE(ATPL(KI)) 132 CONTINUE 135 CONTINUE C C TESTING SINGULARITY OF NORMAL EQUATION MATRIX C CALL EIGRS(A,6,10,EING,Z ,6,WK,IER) CUTOF=EING(6-NP+1)/EING(6) IF(CUTOF.LT.0.001)GO TO 2223 C C GOING FROM DOUBLE TO SIGNAL PRECISION DO 141 I=1,NP DO 142 J=1,NP ATPAD(I,J)=DBLE(ATPA(I,J)) ATCAD(I,J)=DBLE(ATCA(I,J)) 142 CONTINUE 141 CONTINUE C C CALL MINVD(ATPAD,6,NP,DET,IW1,IW2) C XDET=DET CALL SPIN(ATPAD,NP,6,DET,IDEXP) XDET=DET*10**IDEXP C CALL MMULD(RD,6,ATPAD,6,ATPLD,6,NP,NP,1) CALL MMULD(CXD,6,ATCAD,6,ATPAD,6,NP,NP,NP) CALL MMULD(CXXD,6,ATPAD,6,CXD,6,NP,NP,NP) DO 151 I=1,NP R(I)=RD(I) DO 152 J=1,NP ATPA(I,J)=ATPAD(I,J) CXX(I,J)=CXXD(I,J) 152 CONTINUE 151 CONTINUE C DO 2760 KF=1,NP ALFA(KF)=R(KF) DO 2761 JF=1,NP COV(KF,JF)=CXX(KF,JF) 2761 CONTINUE 2760 CONTINUE C ENDIF GO TO 2224 C 2223 NFLAG=-9 C 2224 CONTINUE RETURN END C C________________________________________________________________________ C------------------------------------------------------------------------ C SRB#17 C """""" SUBROUTINE INNER(C,NUM,NP,ALFA,COV,RESO,SIGMAO) C C C C INNER IS USED TO CALCULATE THE INNERMOST ZONE CONTRIBUTIONS AND C ITS ASSOCIATED RMS VALUE IN METERS WITH THE VALUES THAT C HAVE BEEN CALCULATED IN THE SURFACE FITTING ROUTINE. C C INPUT: C """""" FROM SUBROUTINE NUMERC C ----------------------- C C:VECTOR CONTAINING THE RESULTS OF THE SURFACE INTEGRATION C VECTOR WITHIN THE INNERMOST ZONE. C C FROM FUNCTION ICOUNT C --------------------- C NUM:NUMBER OF POINTS USED IN THE SURFACE FITTING SOLUTION C C FROM SUBROUTINE SURF C --------------------- C ALFA:THE CORRECTION VECTOR OF THE SURFACE FITTING SOLUTION C COV:THE COVARIANCE MATRIX OF THE SURFACE FITTING SOLUTION C C OUTPUT: C """"""" C RESO:THE INNERMOST ZONE CONTRIBUTION IN METERS. C SIGMAO:THE ASSOCIATED RMS VALUE. C IMPLICIT REAL*4(A-H,O-Z) DIMENSION ALFA(6),COV(6,6),C(1,6),CT(6,1),CN(1,6),CNCT(1,1), * R(1,1) C RESULT=0.0 STDV=0.0 C DO 14 L=1,6 CT(L,1)=C(1,L) 14 CONTINUE C IF(NP.EQ.1)THEN RESULT=RESULT+C(1,1)*ALFA(1) STDV=STDV+(C(1,1)**2*COV(1,1)) RESO=RESULT SIGMAO=SQRT(STDV) ELSE CALL MMULS(R,1,C,1,ALFA,6,1,NP,1) CALL MMULS(CN,1,C,1,COV,6,1,NP,NP) CALL MMULS(CNCT,1,CN,1,CT,6,1,NP,1) RESO=R(1,1) SIGMAO=SQRT(CNCT(1,1)) ENDIF RETURN END C C________________________________________________________________________ C------------------------------------------------------------------------ C SRB#18 C SUBROUTINE XNTGRT(N,RES,SIGMA,FCT,XLATD,XLONGD,GRAV,SDGRAV, * HEGHT,PHIO,ALAO,DIST,XLAI,XHII) C C XNTGRT IS USED TO INTEGRATE OVER THE INNER AND C OUTER ZONES BY THE PARTIAL STOKE'S EQUATION. C C INPUT: C """""" C FROM SUBROUTINE SCAN C """""""""""""""""""" C N: THE NUMBER OF POINTS IN THE ZONE C C FROM FUNCTION FCT C """"""""""""""""" C FCT: AN EXTERNAL KERNEL C C FROM SUBROUTINE SCAN C """""""""""""""""""" C XLATD: AN ARRAY OF LATITUDES OF POINTS IN THE C ZONE C XLONGD: THE CORRESPONDING ARRAY OF LONGITUDES C GRAV: THE CORRESPONDING ARRAY OF POINT GRAVITY C ANOMALIES OR IF THERE ARE NO PGA MGA IS USED C SDGRAV: THE CORRESPONDING ARRAY OF THE RMS OF THE C PGAS C HEGHT: THE CORRESPONDING ARRAY OF HEIGHTS C C FROM AREA OR DATA C """"""""""""""""" C PHIO: THE LATITUDE OF THE EVALUATION POINT C ALAO: THE LONGITUDE OF THE EVALUATION POINT C DIST: THE CAP SIZE IN RAD(INITIALZED AT BEGINNING OF C PROGAM) C FROM SUBROUTINE CELLSIZE C """"""""""""""""""""""""" C XLAI: THE CELL SIZE IN THE LAMDA DIRECTION C XHII: THE CELL SIZE IN THE PHI DIRECTION C NOTE: C """"" THE CELL SIZE DESCRIBE ABOVE IS FOR THE INNER ZONE C C OUTPUT: C """""""RES: THE ZONAL CONTRIBUTION AS A RESULT OF THE C INTERGRATION OVER THE SPECIFIED SPHERICAL DISTANCE. C SIGMA: ASSOCIATED RMS VALUE OF RES. C C IMPLICIT REAL*4(A-H,O-Z) DIMENSION XLATD(N),XLONGD(N),GRAV(N),SDGRAV(N),HEGHT(N) EXTERNAL FCT C TT=6378137.D0/1.2312529D+7 RES=0.D0 SIGMA=0.D0 DO 10 I=1,N PSI=SPHER(PHIO,ALAO,XLATD(I),XLONGD(I)) IF(PSI.GT.DIST)GOTO 20 XKERNL=FCT(PSI) DAREA=ELEMNT(XLATD(I),XLAI,XHII) RES=RES+(XKERNL*DAREA*GRAV(I)) SIGMA=SIGMA+((XKERNL*SDGRAV(I)*DAREA)**2) 20 CONTINUE 10 CONTINUE RES=RES*TT SIGMA=SQRT(SIGMA)*TT RETURN END C C________________________________________________________________________ C------------------------------------------------------------------------ C SRB#19 C """""" SUBROUTINE TOTAL(RESO,RES1,RES2,SIGMAO,SIGMA1,SIGMA2,GEOID, * GEODAC) C TOTAL IS USED TO COMPUTE THE TOTAL GEOID HEIGHT C AND IT'S RMS VALUE FROM THE THREE ZONES. C C INPUT: C """""" RESO: THE RESULT FROM THE INNERMOST ZONE C RES1: THE RESULT FROM THE INNER ZONE C RES2: THE RESULT FROM THE OUTER ZONE C SIGMAO: THE RMS OF THE RESULT FROM THE INNERMOST ZONE C SIGMA1: THE RMS OF THE RESULT FROM THE INNER ZONE C SIGMA2: THE RMS OF THE RESULT FROM THE OUTER ZONE C C OUTPUT: C """""""GEOID: THE GEOID HEIGHT IN METERS C GEODAC: ASSOCIATED RMS VALUE C C IMPLICIT REAL*4(A-H,O-Z) GEOID=RESO+RES1+RES2 GEODAC=SQRT(SIGMAO**2+SIGMA1**2+SIGMA2**2) RETURN END C________________________________________________________________________ C------------------------------------------------------------------------ C FUNCTION #1 C """"""""""" C ROUND CONVERTS A DECIMAL DEGREE VALUE TO C TO A VALUE ROUNDED OFF TO THE NEAREST MINUTE. C C INPUT: C """"" XPHI: DECIMAL DEGREE VALUE C FUNCTION ROUND(XPHI) IMPLICIT REAL*4(A-H,O-Z) C TWO=FLOAT(2) FIV=FLOAT(5) TEN=FLOAT(10) ONETW=FLOAT(120) W=XPHI*ONETW+FIV IR=W/TEN R=IR ROUND=R*TEN/TWO RETURN END C________________________________________________________________________ C------------------------------------------------------------------------ C FUNCTION #2 C """"""""""" C DEGMIN CONVERTS DEGREES(NOT DECIMAL DEGREES) TO MINUTES C C INPUT: C """"" DEG: NUMBER OF DEGREES C FUNCTION DEGMIN(DEG) IMPLICIT REAL*4(A-H,O-Z) F=FLOAT(60) DEGMIN=DEG*F RETURN END C________________________________________________________________________ C------------------------------------------------------------------------ C FUNCTION #3 C """"""""""" C ICOUNT DETERMINES THE NUMBER OF PARAMETERS TO BE SOLVED C FOR IN THE LEAST SQUARE'S SURFACE FITTING ROUTINE. C C INPUT: C """"" NUM: THE NUMBER OF GRAVTIY ANOMALIES IN THE INNERMOST ZONE C FUNCTION ICOUNT(NUM) IMPLICIT REAL*4(A-H,O-Z) IF(NUM.LE.3)THEN ICOUNT=1 ELSE IF(NUM.EQ.4)THEN ICOUNT=3 ELSE IF(NUM.EQ.5.OR.NUM.EQ.6)THEN ICOUNT=4 ELSE ICOUNT=6 ENDIF ENDIF ENDIF RETURN END C________________________________________________________________________ C------------------------------------------------------------------------ C FUNCTION #4 C """"""""""" C SPHER COMPUTES THE SPHERICAL DISTANCE(IN RADIANS) BETWEEN TWO POINTS 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 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 #5 C """"""""""" C FCT COMPUTES THE MODIFIED SPHERICAL STOKES KERNAL C FROM AN ANGLE(IN RADIANS) IN XNTGTR. C C INPUT: C """"" ANG: DISTANCE IN RADIANS C FUNCTION FCT(ANG) IMPLICIT REAL*4 (A-H,O-Z) PI=DARCOS(-1.D0) A1=1.99727 A2=-3.44927 A3=-173.24417 A4=-32.43544 X=LOG(ANG/2.0) Z=ANG**2 Y=Z*X FCT=A1/ANG+A2*X+A3*Y+A4 RETURN END C________________________________________________________________________ C------------------------------------------------------------------------ C FUNCTION #6 C """"""""""" C ELEMNT COMPUTES THE SURFACE ELEMENT IN XNTGTR C C INPUT: C PHI:LATITUDE OF THE POINT(MIN.) C DPHI:CELL SIZE IN THE LAT. DIRECTION(MIN.) C DLAM:CELL SIZE IN THE LONG. DIRECTION(MIN.) 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 C------------------------------------------------------------------------ C FUNCTION #7 C """"""""""" C THIS FUNCTION COMPUTES THE ATMOSPERIC ATTRACTION EFFECT C ON THE ANOMALIES C C FUNCTION ATMATR(H) IMPLICIT REAL*4(A-H,O-Z) DIMENSION DELTAG(10) C DELTAG(1)=.87 DELTAG(2)=.82 DELTAG(3)=.78 DELTAG(4)=.73 DELTAG(5)=.69 DELTAG(6)=.65 DELTAG(7)=.61 DELTAG(8)=.57 DELTAG(9)=.54 DELTAG(10)=.50 XH=H IF(XH.LT.0.0)XH=0.0 IF(XH.GT.4500.0) THEN XH=4500.0 ENDIF I=(XH/500.0)+1 DH=XH-(I-1)*500.0 ATMATR=DELTAG(I)+(DH/500.0)*(DELTAG(I+1)-DELTAG(I)) RETURN END C________________________________________________________________________ C FUNCTION #8 C """"""""""" C THIS FUNCTION COMPUTES THE INDIRECT EFFECT C FUNCTION ENDEFT(NUMBER,PHIO,ALAO,PHIC,ALAC,XKL,XKPH,H,HMEAN) IMPLICIT REAL*4(A-H,O-Z) DIMENSION PHIC(NUMBER),ALAC(NUMBER),H(NUMBER) C PI=DARCOS(-1.D0) GAMA=9.81 F=1000.0**3 G=66.7/F ROW=2.7 R=6.371*F/1000.0 C=(G*ROW*R**2)/2.0 HMEAN2=0.0 HMEAN3=0.0 IF(HMEAN.GT.0.0)HMEAN2=HMEAN**2 IF(HMEAN.GT.0.0)HMEAN3=HMEAN**3 XB=(PI*G*ROW*HMEAN2)/GAMA C1=C/(GAMA*3.0) GT=0.0 DO 12 I=1,NUMBER HDUM3=0.0 IF(H(I).GT.0.0)HDUM3=H(I)**3 EH=HDUM3-HMEAN3 PSI=SPHER(PHIO,ALAO,PHIC(I),ALAC(I)) IF(PSI.LT.0.0001)GOTO 2457 DAREA=ELEMNT(PHIC(I),XKL,XKPH) XL=(R*PSI)**3 GT=GT+((EH*DAREA)/XL) 2457 CONTINUE 12 CONTINUE ENDEFT=-(GT*C1+XB) RETURN END C________________________________________________________________________ C FUNCTION #9 C """"""""""" C THIS FUNCTION DETERMINES THE AVERAGE HEIGHT C IN THE FOUR CELLS OF THE INNERMOST ZONE C FUNCTION MEAN(NUM2,HED) IMPLICIT REAL*4(A-H,O-Z) DIMENSION HED(NUM2) S=0.0 XNUM2=NUM2 DO 9234 JC=1,NUM2 S=S+HED(JC) 9234 CONTINUE MEAN=S/XNUM2 RETURN END