c - Program geopot97 c - Version 0.4d, 4/14/2004 c - by Dru Smith, National Geodetic Survey c - original programmer of GPTDR : Goad/Tscherning c - All versions prior to 0.4d contain known errors c - Changes: 0.4d: c - Fixed a typo in "function 6" where the flattening c - was being computed incorrectly if the user c - input a & b values. c - Changes: 0.4c: c - Fixed the error I inadvertantly introduced in version 0.4a c - in the subroutine "norm". c - Changes: 0.4b: c - Fixed a "goto 34" line to "goto 134" in Function # 3 c - Changes: 0.4a: c - Modified subroutine "norm" for the computation of the c - J2,J4, etc coefficients to a simpler form. The slight c - change in code yields a very slight change in the c - numbers (around the 10th decimal place) c - Changes: 0.4 c Removed a strict definition of gravity anomaly, and instead c gave the user the choice: c Whenever gravity anomalies are to be calculated, give the c user 2 options: c c Spherical approximation (compatable with Stokes' eqn and c spherical harmonics) c 1) delta_g = -dT/dr - 2*T/r - delta_W*2/r c Eval Pt: P-Q P P P P-Q P c c Direct evaluation to 1st order terms (incompatable with c Stokes' eqn and spherical harmonics) c 2) delta_g =-dT/dh +T*(dgamma/dh)/gamma +delta_W*(dgamma/dh)/gamma cEvalPt: P-Q P P Q Q P-Q Q Q c The "Eval Pt" is the evaluation point, where P is the point c on the geoid and Q the point on the ellipsoid c - Changes: 0.3: c Removed the following errors found while debugging: c 1) sin(lat) calls were replaced with sin(lat*d2r) in general c where they occurred (dtdr computations mostly) c 2) In a type (9) run, I was constantly re-setting the c evaluation height (h) to zero, rather than just on the c first time through. c 3) In a type (8) run, the gravity anomaly was calculated c 4) Changed a typo in type (8) run -- units of N changed c from saying "m**2/s**2" to "m" c 5) Changed the gravity anomaly calculation BACK to that c of geopot96 (that is, the 2nd term is now -2T/r rather c than (dgamma/dh)/gamma ). This allows for a consistancy c between spherical harmonic N and delta-g. c 6) Added the degree zero terms to gravity anomalies for c the non-equivalence of potential (W0<>U0). (The degree c zero term due to non-equivalence of the masses is already c accounted for in the degree zero term of disturbing potential). c 7) In a type #11 run, the total potential call was fixed c to use "fomega" rather than the erroneous use of "domega" c c Added the following: c 1) Added the degree zero component of gravity anomalies that c is induced by W0.ne.U0 (2nd half of equation 2-186 in H/M) c 2) Added a gravity/normal gravity printout to a type #7 run. c - Changes: 0.2b co-ordinates the non-variable nature of the c - call to INITAL and the subroutine header itself. Also removed c - references to non-used variables in INITAL. c - Changes: 0.2a removes variable "OM2" from subroutine c - GPTDR. This variable was a left-over from the days when c - GPTDR used "omega" and did the centrifugal potential part. c - It *shouldn't* have been a problem, since OM2 is undefined, c - and thus *probably* set to zero. But safest just to remove c - the code, right? c - Changes: 0.2 adds the option of inputting an c - (ellipsoidal or orthometric) height DTED and getting out surface c - gravity values on the DTED grid, at the Earth's surface c - Also, the subroutine "getgeom" was added to remove a redundancy c - between numerous functions. Also, subroutine getmsl was added c - to allow one to modify the DTED based on it's divergence from c - best global mean sea level. The bias for NAVD 88 is hard-coded c - in the parameter statement c - Changes: 0.1a added a warning about the mean tide system c - being incompatable with Laplace's equation. ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Program to use geopotential coefficients for point and c c grid calculation of various gravimetric quantities c c c c Dru Smith, 1/13/97 c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c History (geopot97): c c -- Program started by Tscherning and Goad to use Clenshaw c c summation for potential, undulation, anomaly, and c c gravity gradient computation. c c -- Modified by Dennis Milbert into grd360 and geo360 c c (grid and point computations) c c -- Modified by Dru Smith, 1996, into geopot96, a point or c c grid method, for either geoid or ellipsoid quantities. c c Added upward derivatives of height anomaly. c c Added ability to use any reference ellipsoid. c c -- Modified by Dru Smith, 1997, into geopot97 to properly c c account for tides, and clean up the code which was c c cluttered and clunky. c c -- geopot97 work begun 1/13/97 c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c FUNCTIONS PERFORMED (Assuming you've input the coefficients): c (geopot97) c c (1) Given a spherical lat/lon, and a value of GRAVITY potential, c find the radial distance to the surface of that potential. c (2) Given a geometric ellipsoid and an ellipsoidal lat/lon, and c a value of GRAVITY potential, find the ellipsoidal height to the c surface of that potential. c (3) Same as (1), but do the computations for grid points c (4) Same as (2), but do the computations for grid points c c (5) Given a spherical lat/lon/distance, compute the potential c and gravity at that point. c (6) Given a geometric ellipsoid, and an ellipsoidal lat/lon/height, c compute the potential and gravity at that point. c c (7) Given an ellipsoidal normal field, and a lat/lon/height, c compute potential, gravity, disturbing potential, gravity c anomaly, height anomaly, gravity disturbance, deflection c of the vertical, and 3-D gravity gradients at that point. c c (8) Given an ellipsoidal normal field, and either a best fitting c geoid or gravity potential value on the geoid, compute the c geoid undulation and gravity anomaly at that point. c c (9) Given an ellipsoidal normal field, and lat/lon boundaries, c compute the geoid undulation and/or gravity anomaly inside c the masses (theoretically incorrect / mathematically consistant) c c (10) Given an ellipsoidal normal field, and lat/lon boundaries, c compute the 0th, 1st and 2nd order derivatives of the c height anomaly on the ellipsoid. c c - Added for version 0.2: c (11) Given a DTED, compute the surface c gravity value from the spherical harmonic model, on a grid c at the Earth's surface c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c INPUT (General) (geopot97): c 1) The name of a set of fully normalized spherical harmonic c coefficients of the Earth's external GRAVITATIONAL potential. c (The name leads to a hard-coded identification of the c proper scale factors, a and GM, of the coefficient set) c c 2) The coefficient set named in (1) c c 3) The tide system (mean, zero or non) to which the c coefficient set refers (if known). c c 4) A standard value of omega (often taken as 7292115d-11 rad/sec) c c INPUT (Functions 2,4,6): c 1) Geometric parameters of an ellipsoid (a/f or a/b. c *NOT* a and J2, as the conversion to "f" requires c GM and omega, which are non-geometric). c c INPUT (Functions 7,8,9,10): c 1) Any four parameters defining the ellipsoidal normal field c of the Earth's NORMAL external GRAVITY potential. c 2) The tide system of the normal field, if known. c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c IMPORTANT RULES (geopot97): c 1) Never use: N(p,q) = T(p)/gamma(q) c Always use: N(p,q) = T(p)/gamma(q) - [W(p) - U(q)]/gamma(q) c c 2) The Earth spins; The ellipsoid spins; Their spins cancel c c 3) Before using any approximation, ask yourself this: Would I c bet my house on the accuracy of that approximation? ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c - Assumptions(geopot96 and geopot97): c - 1) The program will do calculations INSIDE THE MASSES, but c - the results are NOT guaranteed. Only values outside of c - the masses have any grounds in theoretical correctness. c - 2) The TOTAL CENTRIFUGAL and NORMAL CENTRIFUGAL potentials c - are identical, thus the "DISTURBING" potential will have c - no dependence on OMEGA (being both the TOTAL SPIN and NORMAL SPIN c - rate of the Earth). c - 3) The summation of the degree 1 through NMAX terms (and c - the 1st and 2nd order derivatives of those sums (wrt the c - lat, lon and height) are calculated through a function c - called "GPTDR", written by Tscherning/Goad, which uses c - Clenshaw summations and either un-normalized or c - quasi-normalized coefficients (*NOT* FULLY Normalized, c - though that is what is assumed to come INTO this c - program -- they will be converted to QUASI-normalized c - coefficients later). ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c - Generalizations(geopot96 and geopot97): c - 1) It is generalized that the "GM" and "a" values associated c - with the TOTAL field expansion will not necessarily be c - the same as the "GM" and "a" values associated with the c - NORMAL field. As such, a DISTURBING potential will have c - a NON-ZERO value for the degree zero term. Subroutine c - GPTDR was not set up for such a contingency, and was c - therefore modified to allow for such. ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc implicit real*8 (a-h,o-z) parameter(max = 720,njmax=15,b88=-0.314) c - max = highest possible degree/order of the coefficients to c use in calculations c (max+1)*(max+1)-1 = number of coefficients stored, knowing that all c even zonal S terms are zero (and thus not stored), and c the degree zero term of C is stored in the 0th place. c - njmax => njmax*2 = highest even zonal computed for the c normal field. Thus, njmax=10 means going up to J(20) c - b88 = bias in NAVD88 (needed for type 11 run.) The number c is how far the H(88)=0 level is above the H(true)=0 level c Thus a negative b88 value puts it below. real*8 root((max+1)*2) real*8 smlp1(max+2),cmlp1(max+2),c(0:(max+1)*(max+1)-1) real*8 cj(0:(max+1)*(max+1)-1) real*8 su(10*(max+1)) c real*8 typot1(3),tionpot1(3),cenpot1(3) real*8 typot2(3,3),tionpot2(3,3),cenpot2(3,3) c - tionpot1(3) = g1(3) in GPTDR c - tionpot2(3,3) = g2(3,3) in GPTDR real*8 j2norm,j2bf,jnorm(njmax) c - Fully normalize Legendre Polynomials c real*8 pleg(njmax*2+1) c - Output data real*8 glamn,glomn,dla,dlo integer*4 nla,nlo,ikind c - Functions 3 and 4: logical lfile,lscreen c - Function 3: real*4 raddist(50000) c - Function 4: real*4 outht(50000) c - Function 6: real*4 geoid(50000),ganom(50000) logical lgeoid,lganom c - Function 7: real*4 t0(50000),t1(50000),t2(50000) logical lt0,lt1,lt2 c - Function 11: real*4 elev(50000),outg(50000) c - To force the summation in gptdr: logical force character*1 yesno character*4 tidcoef,tidsum,tidbf character*80 name,fname logical iform common/block1/pi,d2r pi = 2.d0*dasin(1.d0) d2r = pi/180.d0 force = .false. c - Before doing ANYTHING, initialze variables that help us keep c track of our trips to GPTDR call inital c - Initialize the cj vector (quasi-normalized NORMAL coefficients, c - scaled to work with the gm and a of the coefficients). do 771 ii=0,(max+1)*(max+1)-1 cj(ii) = 0.d0 771 continue c - Establish which coefficients, and defining parameters write(6,*) ' ' write(6,*) ' Program geopot97, version 0.4d, 4/14/2004' write(6,*) ' ' write(6,*) ' ' write(6,*) ' ' write(6,*) ' ' write(6,*) ' ***************** WARNING ******************* ' write(6,*) ' Before proceeding, a word of warning is needed' write(6,*) ' regarding the mean tide system. Changing between' write(6,*) ' tide systems is often performed by modifying ' write(6,*) ' the C(2,0) harmonic coefficient and re-inserting' write(6,*) ' it into the spherical harmonic expansion. ' write(6,*) ' The *MEAN* tide system, however, is incompatable' write(6,*)' with Laplaces equation, and therefore can not ' write(6,*) ' be represented by the simple spherical harmonic' write(6,*) ' expansion of the potential. No current fix of ' write(6,*) ' this problem exists. Take heed when using any' write(6,*) ' data given in the mean tide system. Thank you.' write(6,*) ' ***************** WARNING ******************* ' write(6,*) ' ' write(6,*) ' ' write(6,*) ' ' 3 write(6,1) 1 format(1x,'Name of the coefficient set? :',$) read(5,'(a)')name if(name.eq.'egm96' .or. name.eq.'EGM96')then gmcoef = 3.986004415d14 acoef = 6378136.3d0 nmaxcoef = 360 tidcoef = 'non-' elseif(name.eq.'osu91a' .or. name.eq.'OSU91A')then gmcoef = 3.98600436d14 acoef = 6378137.d0 nmaxcoef = 360 tidcoef = '????' else write(6,2) 2 format(1x,' That coefficient set is unknown. What', * 'to do?',/,1x,'1 - Re-enter name',/, * '2 - Proceed and prompt for defining parameters',/, * '3 - Abort program ',/,$) read(5,*)ians if(ians.le.0 .or. ians.ge.3)stop if(ians.eq.1) goto 3 if(ians.eq.2)then write(6,4) 4 format(1x,'GM of the coefficients? ',$) read(5,*)gmcoef write(6,5) 5 format(1x,'a of the coefficients? ',$) read(5,*)acoef write(6,6) 6 format(1x,'maximum degree of the coefficients? ',$) read(5,*)nmaxcoef write(6,7) 7 format(1x,'tidal system of the coefficients?',/, * ' 1 - mean',/,'2 - zero',/,'3 - non-tidal',/,$) read(5,*)ians if(ians.le.0 .or. ians.ge.4)stop if(ians.eq.1)tidcoef = 'mean' if(ians.eq.2)tidcoef = 'zero' if(ians.eq.3)tidcoef = 'non-' endif endif c - Establish format and name of file and open it write(6,23) 23 format(1x,'Is the coefficient file formatted? :',$) read(5,'(a)')yesno iform = .true. if(yesno.eq.'n' .or. yesno.eq.'N')iform = .false. write(6,24) 24 format(1x,'Name of the coefficient file? :',$) read(5,'(a)')fname if(iform)then open(1,file=fname,form='formatted',status='old') else open(1,file=fname,form='unformatted',status='old') endif c - Establish tide system of incoming coefficients write(6,10)tidcoef 10 format(/1x,'Your geopotential model is assumed to be in the ', * 'system of: ',a4,'tide',/,' Change this (y/n)? ',$) read(5,'(a)')yesno if(yesno.eq.'y' .or. yesno.eq.'Y')then write(6,11) 11 format(1x,'New tide system: ',/,'1 - mean',/,'2 - zero',/, * '3 - non-tidal',/,$) read(5,*)ians if(ians.le.0 .or. ians.ge.4)stop if(ians.eq.1)tidcoef = 'mean' if(ians.eq.2)tidcoef = 'zero' if(ians.eq.3)tidcoef = 'non-' endif c - Print out defining parameters write(6,13) 13 format(/,1x,60('*'),/) write(6,8) 8 format(1x,'Your coeffiecients will be read in with the' *' following defining parameters:') write(6,112)gmcoef 112 format(1x,' GM = ',d19.13) write(6,113)acoef 113 format(1x,' a = ',f15.7) write(6,114)nmaxcoef 114 format(1x,' maximum degree = ',i4) write(6,12)tidcoef 12 format(1x,' tidal system = ',a4,' tidal') write(6,36) 36 format(1x,' normalization = fully-normalized') c - Establish Omega write(6,13) write(6,14) c - fomega = "FULL" omega (that is the true rate) c - domega = "DISTURBED" omega (that is, zero) fomega = 7292115d-11 domega = 0.d0 14 format(1x,'Establish spin rate of Earth (which must'/ * 1x,'also be the same rate as any reference '/ * 1x,'ellipsoid you use, for this program) :'/ * 1x,'Omega = 7292115d-11 rad/sec. Ok (y/n)? ',$) read(5,'(a)')yesno if(yesno.eq.'n' .or. yesno.eq.'N')then write(6,15) 15 format(1x,'New value of Omega (rad/sec) : ',$) read(5,*) fomega endif c - Pick which function we'll be doing: write(6,13) write(6,16) write(6,120) write(6,17) 16 format(1x,' Perform which of the following functions?'// *' 1 = Given a spherical lat/lon and a value of GRAVITY'/ *' POTENTIAL, find the radial distance to the surface'/ *' of that potential.'// *' 2 = Given a geometric ellipsoid, an ellipsoidal '/ *' latitude and longitude and a value of GRAVITY'/ *' POTENTIAL, find the ellipsoidal height to the '/ *' surface of that potential.'// *' 3 = Same as function #1, but done on a lat/lon grid'// *' 4 = Same as function #2, but done on a lat/lon grid'// *' 5 = Given a spherical lat/lon/distance, compute the'/ *' GRAVITY POTENTIAL, and GRAVITATIONAL POTENTIAL at'/ *' that location.'// *' 6 = Given a geometric ellipsoid, and an ellipsoidal'/ *' latitude, longitude and height, find the GRAVITY,'/ *' GRAVITY POTENTIAL and GRAVITATIONAL POTENTIAL at'/ *' that location.'/) 120 format( *' 7 = Given and ellipsoidal normal gravity field, and an'/ *' ellipsoidal latitude, longitude, and height, compute'/ *' the gravity potential, gravitational potential, and'/ *' gravity *AND ALSO* the gravity anomaly, height anomaly,'/ *' gravity disturbance, deflections of the vertical, and'/ *' 3-D gradients of gravity at that location.'// *' 8 = Given an ellipsoidal normal gravity field, and an'/ *' ellipsoidal latitude, longitude, and a defined'/ *' geoid (either through a best fit ellipsoid or a gravity'/ *' potential value), give the geoid undulation and '/ *' gravity anomaly.'/) 17 format( *' 9 = Given an ellipsoidal normal gravity field, and'/ *' boundaries of ellipsoidal latitude and longitude,'/ *' compute the geoid undulations and/or gravity anomalies'/ *' on a grid, in the masses.'// *' 10= Given an ellipsoidal normal gravity field, and'/ *' boundaries of ellipsoidal latitude and longitude,'/ *' compute the height anomaly, as well as the 1st and'/ *' 2nd order upward derivatives of the height anomaly,'/ *' all on the ellipsoid.'// *' 11= Given a DTED (H or h), compute a grid of gravity values'/ *' at the Surface of the Earth'/) read(5,*)ifun if(ifun.lt.1 .or. ifun.gt.11)stop c - Determine maximum degree we'll be using write(6,18)nmaxcoef 18 format(/1x,'The maximum degree in the coefficient set is:',i4,/ * 1x,'What is the maximum degree you want used? ',$) read(5,*)nmax c - Read in fully-normalized coefficients (assume 2,0 is first one) c - C(0,0) = c(0) c - S(0,0) = not stored c - C(1,0) = c(1) c - S(1,0) = not stored c - C(1,1) = c(2) c - S(1,1) = c(3) c - C(2,0) = c(4) <- this is why isto starts at isto=4 c - S(2,0) = not stored c - C(2,1) = c(5) c - S(2,1) = c(6) c - C(2,2) = c(7) c - S(2,2) = c(8) etc.... c(0) = 1.d0 ! <- may change if we work with disturbing potentials c(1) = 0.d0 c(2) = 0.d0 c(3) = 0.d0 isto = 4 c - Formatted if(iform)then do 19 n=2,nmax do 20 m = 0,n read(1,*)nnn,mmm,ccc,sss c(isto) = ccc if(m.ne.0)then isto = isto + 1 c(isto) = sss endif isto = isto + 1 if(nnn.ne.n .or. mmm.ne.m)stop 20 continue 19 continue c - Unformatted else do 21 n=2,nmax do 22 m = 0,n read(1)ccc,sss c(isto) = ccc if(m.ne.0)then isto = isto + 1 c(isto) = sss endif isto = isto + 1 22 continue 21 continue endif c - Convert Fully Normalized Coefficients into c - quasi-normalized coefficients, for efficiency in GPTDR c - and also compute a table of square roots for quick c - look-up. call setcm(nmax,root,smlp1,cmlp1,c,su,max) c - Find the tide system of the output call gettid(tidsum,tidcoef) ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c - FUNCTION # 1 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc if(ifun.eq.1)then 78 write(6,25) 25 format(/1x,' Input spherical latitude (dec. deg): ',$) read(5,*)xlat write(6,26) 26 format(1x,' Input spherical longitude (dec. deg): ',$) read(5,*)xlon write(6,27) 27 format(1x,' Input GRAVITY potential value (m**2/s**2): ',$) read(5,*)xpot c - Take care of the tides if(tidsum.ne.tidcoef)call fixc20(tidcoef,tidsum,c(4), * gmcoef,acoef) c - iterate until we find the surface to 0.01 mm r = 6371000.d0 34 isw = 0 call pot(xlat,xlon,r,fomega,tionpot,tionpot1,tionpot2, * cenpot,cenpot1,cenpot2,typot,typot1,typot2,isw,nmax, * root,smlp1,cmlp1,c,su,max,gmcoef,acoef,force) dv = typot - xpot c write(6,*) r,typot,dv if(dabs(dv).lt.0.000098)then write(6,60)xlat 60 format(1x,'Spherical Latitude (dec. deg.) = ',f20.8) write(6,61)xlon 61 format(1x,'Spherical Longitude(dec. deg.) = ',f20.8) write(6,65)typot 65 format(1x,'Gravity Pot. (m**2/s**2) = ',f20.8) write(6,62)r 62 format(1x,'Radial Distance (m) = ',f20.8) write(6,63)tionpot 63 format(1x,'Gravitational Pot. (m**2/s**2) = ',f20.8) write(6,64)cenpot 64 format(1x,'Centrifugal Pot. (m**2/s**2) = ',f20.8) write(6,72)tidsum 72 format(1x,'Tide system = ',15x,a4,' tide') else r = r + dv/9.8 goto 34 endif write(6,77) 77 format(1x,'Go again (y/n)? ',$) read(5,'(a)')yesno if(yesno.eq.'y' .or. yesno.eq.'Y')goto 78 stop ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c - FUNCTION # 2 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc elseif(ifun.eq.2)then c - Get the geometric ellipsoid parameters call getgeom(anorm,finvnorm) 80 write(6,31) 31 format(/1x,'Input ellipsoidal latitude (dec. deg.): ',$) read(5,*)geodlat write(6,32) 32 format(1x,'Input ellipsoidal longitude (dec. deg.): ',$) read(5,*)geodlon write(6,33) 33 format(1x,' Input GRAVITY potential value (m**2/s**2): ',$) read(5,*)xpot c - Take care of the tides if(tidsum.ne.tidcoef)call fixc20(tidcoef,tidsum,c(4), * gmcoef,acoef) h = 0 35 call ell2sph(geodlat,geodlon,h,anorm,finvnorm,xlat,xlon,r) isw = 0 call pot(xlat,xlon,r,fomega,tionpot,tionpot1,tionpot2, * cenpot,cenpot1,cenpot2,typot,typot1,typot2,isw,nmax, * root,smlp1,cmlp1,c,su,max,gmcoef,acoef,force) dv = typot - xpot if(dabs(dv).lt.0.000098)then write(6,66)geodlat 66 format(1x,'Ellipsoidal Latitude (dec. deg.) = ',f20.8) write(6,67)geodlon 67 format(1x,'Ellipsoidal Longitude(dec. deg.) = ',f20.8) write(6,68)typot 68 format(1x,'Gravity Pot. (m**2/s**2) = ',f20.8) write(6,69)h 69 format(1x,'Ellipsoidal Height (m) = ',f20.8) write(6,70)tionpot 70 format(1x,'Gravitational Pot. (m**2/s**2) = ',f20.8) write(6,71)cenpot 71 format(1x,'Centrifugal Pot. (m**2/s**2) = ',f20.8) write(6,73)tidsum 73 format(1x,'Tide system = ',15x,a4,' tide') else h = h + dv/9.8 goto 35 endif write(6,79) 79 format(1x,'Go again (y/n)? ',$) read(5,'(a)')yesno if(yesno.eq.'y' .or. yesno.eq.'Y')goto 80 stop ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c - FUNCTION # 3 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc elseif(ifun.eq.3)then c - Get the grid and spacings in spherical coordinates call getgrid2(glamn,glomn,dla,dlo,nla,nlo,ikind) c - Output to file or screen lfile = .false. lscreen = .false. write(6,152) 152 format(1x,'Output to a file(y/n)? :',$) read(5,'(a)')yesno if(yesno.eq.'y' .or. yesno.eq.'Y')then lfile = .true. write(6,153) 153 format(1x,' Output file name? :',$) read(5,'(a)')fname open(31,file=fname,status='new',form='unformatted') write(31)glamn,glomn,dla,dlo,nla,nlo,ikind endif write(6,154) 154 format(1x,'Output to screen(y/n)? :',$) read(5,'(a)')yesno if(yesno.eq.'y' .or. yesno.eq.'Y')then lscreen = .true. endif c - Get the gravity potential value write(6,127) 127 format(1x,' Input GRAVITY potential value (m**2/s**2): ',$) read(5,*)xpot c - Take care of the tides if(tidsum.ne.tidcoef)call fixc20(tidcoef,tidsum,c(4), * gmcoef,acoef) c - Set initial radial vector r = 6371000.d0 c - Loop over the grid, iterating at each point to 0.01 mm do 150 irow = 1,nla xlat = glamn + (irow-1)*dla do 151 icol = 1,nlo xlon = glomn + (icol-1)*dlo 134 isw = 0 call pot(xlat,xlon,r,fomega,tionpot,tionpot1,tionpot2, * cenpot,cenpot1,cenpot2,typot,typot1,typot2,isw,nmax, * root,smlp1,cmlp1,c,su,max,gmcoef,acoef,force) dv = typot - xpot c write(6,*) r,typot,dv if(dabs(dv).lt.0.000098)then if(lscreen)then write(6,260)xlat 260 format(1x,'Spherical Latitude (dec. deg.) = ',f20.8) write(6,261)xlon 261 format(1x,'Spherical Longitude(dec. deg.) = ',f20.8) write(6,265)typot 265 format(1x,'Gravity Pot. (m**2/s**2) = ',f20.8) write(6,262)r 262 format(1x,'Radial Distance (m) = ',f20.8) write(6,263)tionpot 263 format(1x,'Gravitational Pot. (m**2/s**2) = ',f20.8) write(6,264)cenpot 264 format(1x,'Centrifugal Pot. (m**2/s**2) = ',f20.8) write(6,272)tidsum 272 format(1x,'Tide system = ',15x,a4,' tide') endif if(lfile)then raddist(icol) = r endif else r = r + dv/9.8 goto 134 endif 151 continue if(lfile)write(31)(raddist(j),j=1,nlo) 150 continue stop ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c - FUNCTION # 4 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc elseif(ifun.eq.4)then c - Get the geometric ellipsoid parameters call getgeom(anorm,finvnorm) c - Get the grid and spacings in geodetic coordinates call getgrid(glamn,glomn,dla,dlo,nla,nlo,ikind) c - Output to file or screen lfile = .false. lscreen = .false. write(6,155) 155 format(1x,'Output to a file(y/n)? :',$) read(5,'(a)')yesno if(yesno.eq.'y' .or. yesno.eq.'Y')then lfile = .true. write(6,156) 156 format(1x,' Output file name? :',$) read(5,'(a)')fname open(32,file=fname,status='new',form='unformatted') write(32)glamn,glomn,dla,dlo,nla,nlo,ikind endif write(6,157) 157 format(1x,'Output to screen(y/n)? :',$) read(5,'(a)')yesno if(yesno.eq.'y' .or. yesno.eq.'Y')then lscreen = .true. endif c - Get the gravity potential value write(6,158) 158 format(1x,' Input GRAVITY potential value (m**2/s**2): ',$) read(5,*)xpot c - Take care of the tides if(tidsum.ne.tidcoef)call fixc20(tidcoef,tidsum,c(4), * gmcoef,acoef) c - Set initial height value h = 0.d0 c - Loop over the grid, iterating at each point to 0.01 mm do 160 irow = 1,nla geodlat = glamn + (irow-1)*dla do 161 icol = 1,nlo write(6,*) ' irow, icol = ',irow,icol geodlon = glomn + (icol-1)*dlo 162 call ell2sph(geodlat,geodlon,h,anorm,finvnorm,xlat,xlon,r) isw = 0 call pot(xlat,xlon,r,fomega,tionpot,tionpot1,tionpot2, * cenpot,cenpot1,cenpot2,typot,typot1,typot2,isw,nmax, * root,smlp1,cmlp1,c,su,max,gmcoef,acoef,force) dv = typot - xpot c write(6,*) r,typot,dv if(dabs(dv).lt.0.000098)then if(lscreen)then write(6,360)geodlat 360 format(1x,'Geodetic Latitude (dec. deg.) = ',f20.8) write(6,361)geodlon 361 format(1x,'Geodetic Longitude(dec. deg.) = ',f20.8) write(6,365)typot 365 format(1x,'Gravity Pot. (m**2/s**2) = ',f20.8) write(6,362)h 362 format(1x,'Ellipsoidal Height (m) = ',f20.8) write(6,363)tionpot 363 format(1x,'Gravitational Pot. (m**2/s**2) = ',f20.8) write(6,364)cenpot 364 format(1x,'Centrifugal Pot. (m**2/s**2) = ',f20.8) write(6,372)tidsum 372 format(1x,'Tide system = ',15x,a4,' tide') endif if(lfile)then outht(icol) = h endif else h = h + dv/9.8 goto 162 endif 161 continue if(lfile)write(32)(outht(j),j=1,nlo) 160 continue stop ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c - FUNCTION # 5 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc elseif(ifun.eq.5) then 82 write(6,38) 38 format(1x,' Input spherical latitude (dec. deg): ',$) read(5,*)xlat write(6,39) 39 format(1x,' Input spherical longitude (dec. deg): ',$) read(5,*)xlon write(6,40) 40 format(1x,' Input radial distance (meters): ',$) read(5,*)r c - Take care of the tides if(tidsum.ne.tidcoef)call fixc20(tidcoef,tidsum,c(4), * gmcoef,acoef) isw = 0 call pot(xlat,xlon,r,fomega,tionpot,tionpot1,tionpot2, * cenpot,cenpot1,cenpot2,typot,typot1,typot2,isw,nmax, * root,smlp1,cmlp1,c,su,max,gmcoef,acoef,force) write(6,41)xlat 41 format(1x,'Spherical Latitude (dec. deg.) = ',f20.8) write(6,42)xlon 42 format(1x,'Spherical Longitude(dec. deg.) = ',f20.8) write(6,43)r 43 format(1x,'Radial Distance (m) = ',f20.8) write(6,44)tionpot 44 format(1x,'Gravitational Pot. (m**2/s**2) = ',f20.8) write(6,45)cenpot 45 format(1x,'Centrifugal Pot. (m**2/s**2) = ',f20.8) write(6,46)typot 46 format(1x,'Gravity Pot. (m**2/s**2) = ',f20.8) write(6,75)tidsum 75 format(1x,'Tide system = ',15x,a4,' tide') write(6,81) 81 format(1x,'Go again (y/n)? ',$) read(5,'(a)')yesno if(yesno.eq.'y' .or. yesno.eq.'Y')goto 82 stop ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c - FUNCTION # 6 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc elseif(ifun.eq.6) then c - Get the geometric ellipsoid parameters 84 call getgeom(anorm,finvnorm) c - Get the ellipsoidal coordinates write(6,51) 51 format(1x,'Enter ellipsoidal latitude (dec. deg.): ',$) read(5,*)geodlat write(6,52) 52 format(1x,'Enter ellipsoidal longitude (dec. deg.): ',$) read(5,*)geodlon write(6,53) 53 format(1x,' Input ellipsoidal height (meters): ',$) read(5,*)h c - Take care of the tides if(tidsum.ne.tidcoef)call fixc20(tidcoef,tidsum,c(4), * gmcoef,acoef) call ell2sph(geodlat,geodlon,h,anorm,finvnorm,xlat,xlon,r) isw = 0 call pot(xlat,xlon,r,fomega,tionpot,tionpot1,tionpot2, * cenpot,cenpot1,cenpot2,typot,typot1,typot2,isw,nmax, * root,smlp1,cmlp1,c,su,max,gmcoef,acoef,force) write(6,54)geodlat 54 format(1x,'Ellipsoidal Latitude (dec. deg.) = ',f20.8) write(6,55)geodlon 55 format(1x,'Ellipsoidal Longitude(dec. deg.) = ',f20.8) write(6,56)h 56 format(1x,'Ellipsoidal Height (m) = ',f20.8) write(6,57)tionpot 57 format(1x,'Gravitational Pot. (m**2/s**2) = ',f20.8) write(6,58)cenpot 58 format(1x,'Centrifugal Pot. (m**2/s**2) = ',f20.8) write(6,59)typot 59 format(1x,'Gravity Pot. (m**2/s**2) = ',f20.8) write(6,76)tidsum 76 format(1x,'Tide system = ',15x,a4,' tide') write(6,83) 83 format(1x,'Go again (y/n)? ',$) read(5,'(a)')yesno if(yesno.eq.'y' .or. yesno.eq.'Y')goto 84 stop ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c - FUNCTION # 7 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc elseif(ifun.eq.7) then c - NOTE: In THIS particular function, there is no need to c - define the geoid through "bestfit", as everything c - will be relative to the assumption that Wsurf=Utelluroid, c - since a dW correction at altitude is both complicated c - and undeveloped. c call bestfit(gmcoef,fomega,abf,j2bf,finvbf,w0,tidbf,tidcoef) c - get the 4 parameters and the tide system of the normal field, c - BUT return the normal field in the tidsum tide system call getref(anorm,bnorm,finvnorm,j2norm,gmnorm,fomega,tidsum) c - Now, anorm,bnorm,finvnorm,j2norm,gmnorm,fomega are in the c - tidsum system and need no further tinkering :) c - calculate normal values from the normal field call norm(anorm,j2norm,gmnorm,fomega,bnorm,jnorm,njmax, * u0norm,gammae,gammap) c - The "jnorm" values are NOT fully-normalized here. They c - will be fully normalized in subroutine dist. c - Now, correct the C(2,0) coefficient for tides, if necessary if(tidsum.ne.tidcoef)then call fixc20(tidcoef,tidsum,c(4),gmcoef,acoef) endif c - At this point, everything is in the "tidsum" tide system. c - get the location of the computation point 91 write(6,92) 92 format(/1x,'Input ellipsoidal latitude (dec. deg.): ',$) read(5,*)geodlat write(6,93) 93 format(1x,'Input ellipsoidal longitude (dec. deg.): ',$) read(5,*)geodlon write(6,94) 94 format(1x,'Input ellipsoidal height (meters): ',$) read(5,*)h c - calculate normal gravity at the point c - as of 1/29/97, the following subroutines are either c - corrupt or (more likely) the actual formulation c - is bad. Bad meaning disagreements to 0.5 mgals or so c - (function gam(*) is ok, and gamh(*) replaces the need c - for the other two. gamma = gam(geodlat,anorm,bnorm,gammae,gammap) cbad dgamdh = dgdh(geodlat,anorm,finvnorm,gamma,fomega) cbad d2gamdh2 = d2gdh2(gamma,anorm) cbad gammah = gamma + h*dgamdh + (0.5d0)*h*h*d2gamdh2 gammah = gamh(anorm,bnorm,gmnorm,fomega,h,geodlat,gamma) c - Convert to spherical coordinates call ell2sph(geodlat,geodlon,h,anorm,finvnorm,xlat,xlon,r) c------------------------ c - STEP #1: c - Compute the gravity/gravitational/centrifugal total potentials c do 802 i=0,20 c write(6,*) ' c(i) = ',c(i) c write(6,*) ' cj(i)= ',cj(i) c 802 continue isw = 1 call pot(xlat,xlon,r,fomega,tionpot,tionpot1,tionpot2, * cenpot,cenpot1,cenpot2,typot,typot1,typot2,isw,nmax, * root,smlp1,cmlp1,c,su,max,gmcoef,acoef,force) c - Store the 3 potentials for printing xp1 = tionpot xp2 = typot xp3 = cenpot c - Store the 3 derivs for computing gamma from disturbings later xg1 = typot1(1) xg2 = typot1(2) xg3 = typot1(3) gravityh = dsqrt(typot1(1)**2+typot1(2)**2+typot1(3)**2) c write(6,*) ' BACK FROM GPTDR using total field....' c write(6,*) ' typot = ',typot c write(6,*) ' typot1(1) = ',typot1(1) c write(6,*) ' typot1(2) = ',typot1(2) c write(6,*) ' typot1(3) = ',typot1(3) c write(6,*) ' tionpot = ',tionpot c write(6,*) ' tionpot1(1)= ',tionpot1(1) c write(6,*) ' tionpot1(2)= ',tionpot1(2) c write(6,*) ' tionpot1(3)= ',tionpot1(3) c write(6,*) ' cenpot = ',cenpot c write(6,*) ' cenpot1(1)= ',cenpot1(1) c write(6,*) ' cenpot1(2)= ',cenpot1(2) c write(6,*) ' cenpot1(3)= ',cenpot1(3) c write(6,*) ' gravityh = ',gravityh c write(6,*) ' *************************************' c--------------------------- c - STEP #1a: c - Compute the NORMAL values c - Create the cj vector, containing the quasi-normalized j coefs c - scaled to work with gmcoef and acoef call mknorm(cj,max,nmax,gmcoef,acoef,jnorm,njmax,gmnorm,anorm) c - Compute potential from the normal field alone c do 801 i=0,20 c write(6,*) ' c(i) = ',c(i) c write(6,*) ' cj(i)= ',cj(i) c 801 continue c isw = 1 c force = .true. c call pot(xlat,xlon,r,fomega,tionpot,tionpot1,tionpot2, c * cenpot,cenpot1,cenpot2,typot,typot1,typot2,isw,nmax, c * root,smlp1,cmlp1,cj,su,max,gmcoef,acoef,force) c force = .false. c xnp1 = tionpot c xnp2 = typot c xnp3 = cenpot c c gammah2= dsqrt(typot1(1)**2+typot1(2)**2+typot1(3)**2) c c write(6,*) ' BACK FROM GPTDR using normal field....' c write(6,*) ' typot = ',typot c write(6,*) ' typot1(1) = ',typot1(1) c write(6,*) ' typot1(2) = ',typot1(2) c write(6,*) ' typot1(3) = ',typot1(3) c c write(6,*) ' tionpot = ',tionpot c write(6,*) ' tionpot1(1)= ',tionpot1(1) c write(6,*) ' tionpot1(2)= ',tionpot1(2) c write(6,*) ' tionpot1(3)= ',tionpot1(3) c c write(6,*) ' cenpot = ',cenpot c write(6,*) ' cenpot1(1)= ',cenpot1(1) c write(6,*) ' cenpot1(2)= ',cenpot1(2) c write(6,*) ' cenpot1(3)= ',cenpot1(3) c c write(6,*) ' gammah2= ',gammah2 c write(6,*) ' *************************************' c c write(6,*) ' dist tionpot = ',xp1-xnp1 c write(6,*) ' dist typot = ',xp2-xnp2 c write(6,*) ' dist cenpot = ',xp3-xnp3 c write(6,*) ' *************************************' c - Now try it with a different method...not full clenshaw c nn = njmax*2+1 c call legpol(xlat,nn,pleg) c call sumnorm(xlat,xlon,r,gmnorm,anorm,jnorm,njmax,pleg, c * xnpn,xnpn1,xnpn2,xnpn3,xypn,xypn1,xypn2,xypn3, c * xcpn,xcpn1,xcpn2,xcpn3,fomega) c c write(6,*) ' BACK FROM SUMNORM...' c write(6,*) ' sumnorm: normal tionpot: ',xnpn c write(6,*) ' sumnorm: normal tionpot/N: ',xnpn1 c write(6,*) ' sumnorm: normal tionpot/E: ',xnpn2 c write(6,*) ' sumnorm: normal tionpot/U: ',xnpn3 c c write(6,*) ' sumnorm: normal typot: ',xypn c write(6,*) ' sumnorm: normal typot/N: ',xypn1 c write(6,*) ' sumnorm: normal typot/E: ',xypn2 c write(6,*) ' sumnorm: normal typot/U: ',xypn3 c c write(6,*) ' sumnorm: normal cenpot: ',xcpn c write(6,*) ' sumnorm: normal cenpot/N: ',xcpn1 c write(6,*) ' sumnorm: normal cenpot/E: ',xcpn2 c write(6,*) ' sumnorm: normal cenpot/U: ',xcpn3 c c gammasu = dsqrt(xypn1**2+xypn2**2+xypn3**2) c write(6,*) ' sumnorm gammasu = ',gammasu c write(6,*) ' *************************************' c - Lastly try a closed formula for gamma at height: c xxm = fomega*fomega*anorm*anorm*bnorm/gmnorm c fnorm = (anorm-bnorm)/anorm c gammah3 = gamma*(1.d0 - (2.d0/anorm)*(1.d0 + fnorm c * + xxm - 2*fnorm*dsin(geodlat*d2r)**2)*h c * + (3.d0/anorm/anorm)*h*h) c write(6,*) ' gammah3 (closed form) = ',gammah3 c write(6,*) ' *************************************' c------------------------------ c - STEP #2: c - Compute the gravity/gravitational/centrifugal disturbing potentials c - Calculate the disturbing coefficients, set up such c - that they work with gmcoef and acoef (which means c - you may need to scale the j2 coefficients from c - working with gmnorm and anorm) call dist(c,max,nmax,gmcoef,acoef,jnorm,njmax,gmnorm,anorm) c - the "c" vector now contains the fully normalized disturbing c - coefficients c - Do the summation c - All "*pot" values are DISTURBING potential values!!! isw = 2 c - Don't have to "force", as we've changed the order, c - (that is, isw), and thus gptdr has a built in re-summation call pot(xlat,xlon,r,domega,tionpot,tionpot1,tionpot2, * cenpot,cenpot1,cenpot2,typot,typot1,typot2,isw,nmax, * root,smlp1,cmlp1,c,su,max,gmcoef,acoef,force) c - STEP #2a: c - Compute Normal values by just subtracting the disturbing from total xnp1 = xp1 - tionpot xnp2 = xp2 - typot xnp3 = xp3 - cenpot c - Compute a gamma value to check against the closed formula gammachk = dsqrt( (xg1-typot1(1))**2+(xg2-typot1(2))**2 * + (xg3-typot1(3))**2 ) dd = (gammah - gammachk) * 1.d5 if(dabs(dd).gt.0.01)then write(6,117) 117 format(1x,'***WARNING*** The closed formula for gamma',/, * 'disagrees with that implied by summing the',/, * 'Normal field coefficients by more than ',/, * '0.01 mgals! ***') write(6,*) ' gammah (gamh subrtn) = ',gammah write(6,*) ' gammachk(g - Delta-g) = ',gammachk write(6,*) ' dif, mgals = ',dd endif c-------------------------------------- c - STEP #3: c - Calculate gravimetric quantities of interest c - Height Anomaly (Assuming W(surf) = U(telluroid) htanom = typot/gammah c - Deflections of the Vertical (assuming d(gamma)/ds = 0) c - Convert deflections from radians to arcseconds defns = ((-typot1(1)/gamma)/d2r)*3600.d0 defew = ((-typot1(2)/gamma)/d2r)*3600.d0 c - Get dr/dh, needed for gravity disturbance: xw = dsqrt(1.d0 - e2norm*dsin(geodlat*d2r)**2) xn = anorm/xw drdh=((xn+h)-xn*e2norm*dsin(geodlat*d2r)**2)/r c - Gravity Disturbance (two ways; convert both to mgals) grdist1 = (gravityh - gammah)*1.d5 c write(6,*) ' typot1(1)*1.d5 = ',typot1(1)*1.d5 c write(6,*) ' typot1(2)*1.d5 = ',typot1(2)*1.d5 c write(6,*) ' typot1(3)*1.d5 = ',typot1(3)*1.d5 c write(6,*) ' drdh = ',drdh c xxx = dsqrt(typot1(1)**2+typot1(2)**2+typot1(3)**2) c write(6,*) ' xxx = ',xxx*1.d5 c grdist2 = (-typot1(3)*drdh)*1.d5 c - Gravity Anomaly (Assuming W(surf) = U(telluroid)): c - Modified for version 0.4: call getdgtyp(idgtyp) if(idgtyp.eq.1)then dtdr = typot1(3) anom = (-dtdr - 2.d0*typot/r)*1.d5 else if(h.ne.0)then write(6,*) ' I cannot calculate dgamma/dh' write(6,*) ' at altitudes, so I will give' write(6,*) ' the spherical approximation' write(6,*) ' to the gravity anomaly' dtdr = typot1(3) anom = (-dtdr - 2.d0*typot/r)*1.d5 else dtdh = typot1(3) * drdh dgamdh = dgdh(geodlat,anorm,finvnorm,gamma,fomega) anom = (-dtdh + typot*dgamdh/gamma)*1.d5 endif endif c - Gravity Gradients in mgals/km: x11 = typot2(1,1)*1.d8 x12 = typot2(1,2)*1.d8 x13 = typot2(1,3)*1.d8 x21 = typot2(2,1)*1.d8 x22 = typot2(2,2)*1.d8 x23 = typot2(2,3)*1.d8 x31 = typot2(3,1)*1.d8 x32 = typot2(3,2)*1.d8 x33 = typot2(3,3)*1.d8 c - Print it all out neatly: write(6,103)name,tidsum 103 format(1x,60('*'),/, * 1x,'Coefficient Set: ',a9, * 8x,'Tide System: ',a4,'tide') write(6,104)gmnorm,anorm,j2norm,fomega 104 format(1x,'Normal Field: GM = ',d20.10,5x, * 'a = ',d20.10,/, * 1x,' J2 = ',d20.10,5x, * 'omega = ',d20.10,/) write(6,105)geodlat,geodlon,h 105 format(1x,'Geodetic Latitude(deg) : ',f19.11,/, * 1x,'Geodetic Longitude(deg): ',f19.11,/, * 1x,'Ellipsoidal Height(m) : ',f19.11,/) write(6,106)xp1,xp3,xp2 106 format(1x,'Total Gravitational Potential (m**2/s**2): ', *f20.10,/, * 1x,'Total Centrifugal Potential (m**2/s**2): ', *f20.10,/, * 1x,'Total Gravity Potential (m**2/s**2): ', *f20.10,/) write(6,115)xnp1,xnp3,xnp2 115 format(1x,'Normal Gravitational Potential (m**2/s**2): ', *f20.10,/, * 1x,'Normal Centrifugal Potential (m**2/s**2): ', *f20.10,/, * 1x,'Normal Gravity Potential (m**2/s**2): ', *f20.10,/) write(6,116)tionpot,cenpot,typot 116 format(1x,'Dist. Gravitational Potential (m**2/s**2): ', *f20.10,/, * 1x,'Dist. Centrifugal Potential (m**2/s**2): ', *f20.10,/, * 1x,'Dist. Gravity Potential (m**2/s**2): ', *f20.10,//) write(6,107)htanom 107 format(1x,'Height Anomaly (m, between Wsurf=Utell) : ',f15.5) write(6,108)anom,defns,defew 108 format(1x,'Gravity Anomaly (mgals, btwn Wsurf=Utell): ',f10.5 * ,/, * 1x,'Deflection N/S (arcsec): ',f10.5,5x, * 'Deflection E/W (arcsec): ',f10.5,/) write(6,708)gravityh,gammachk 708 format(1x,'Total Gravity at h (m/s**2) = ',f15.12,/, * 1x,'Normal Gravity at h (m/s**2) = ',f15.12) write(6,109)grdist1 109 format(1x, * 'Gravity Disturbance (abs(g)-abs(gamma), mgals ): ',f10.5,/) write(6,110)x11,x12,x13,x22,x23,x33 110 format(1x,'Gravity Gradients (mgals/km):',/, * 11x,' dNorth ',3x,' dEast ',3x,' dUp ',/, * 1x,'dNorth ',f10.7,3x,f10.7,3x,f10.7,/, * 1x,'dEast ',10x,3x,f10.7,3x,f10.7,/, * 1x,'dUp ',10x,3x,10x,3x,f10.7,//) write(6,111) 111 format(1x,'Go again? :',$) read(5,'(a)')yesno c - If we go again, go back to TOTAL coefficients for 1st calc if(yesno.eq.'y' .or. yesno.eq.'Y')then call undist(c,max,nmax,gmcoef,acoef,jnorm,njmax,gmnorm * ,anorm) goto 91 endif stop ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c - FUNCTION # 8 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc elseif(ifun.eq.8) then call bestfit(gmcoef,fomega,abf,j2bf,finvbf,w0,tidbf,tidcoef) call getref(anorm,bnorm,finvnorm,j2norm,gmnorm,fomega,tidsum) fnorm = 1.d0/finvnorm e2norm = 2.d0*fnorm - fnorm*fnorm call norm(anorm,j2norm,gmnorm,fomega,bnorm,jnorm,njmax, * u0norm,gammae,gammap) c - Now, correct the C(2,0) coefficient for tides, if necessary if(tidsum.ne.tidcoef)then call fixc20(tidcoef,tidsum,c(4),gmcoef,acoef) endif c - At this point, everything is in the "tidsum" tide system. c - Calculate the disturbing coefficients, set up such c - that they work with gmcoef and acoef (which means c - you may need to scale the j2 coefficients from c - working with gmnorm and anorm) call dist(c,max,nmax,gmcoef,acoef,jnorm,njmax,gmnorm,anorm) c - Get the computation point 480 write(6,431) 431 format(/1x,'Input ellipsoidal latitude (dec. deg.): ',$) read(5,*)geodlat write(6,432) 432 format(1x,'Input ellipsoidal longitude (dec. deg.): ',$) read(5,*)geodlon c - Get normal gravity on the normal field ellipsoid gamma = gam(geodlat,anorm,bnorm,gammae,gammap) dgamdh = dgdh(geodlat,anorm,finvnorm,gamma,fomega) c - Iterate to find the undulation h = 0.d0 435 call ell2sph(geodlat,geodlon,h,anorm,finvnorm,xlat,xlon,r) isw = 1 call pot(xlat,xlon,r,domega,tionpot,tionpot1,tionpot2, * cenpot,cenpot1,cenpot2,typot,typot1,typot2,isw,nmax, * root,smlp1,cmlp1,c,su,max,gmcoef,acoef,force) und = (typot/gamma) - (w0 - u0norm)/gamma dv = h - und0 if(dabs(dv).lt.0.001)then c - Compute gravity anomaly after we're done iterating xw = dsqrt(1.d0 - e2norm*dsin(geodlat*d2r)**2) xn = anorm/xw drdh=((xn+h)-xn*e2norm*dsin(geodlat*d2r)**2)/r dtdr = typot1(3) dtdh = dtdr * drdh c - v 0.3: Added the zero order term (eqn 2-186, H/M) c - v 0.4: Gave a choice on anomaly computation call getdgtyp(idgtyp) if(idgtyp.eq.1)then anom = ( -dtdr - 2.d0*typot/r )*(1.d5) * + (2*(w0-u0norm)/r)*(1.d5) else anom = ( -dtdh + (dgamdh*typot)/gamma )*(1.d5) * - (dgamdh*(w0-u0norm)/gamma)*(1.d5) endif c write(6,*) ' drdh = ',drdh c write(6,*) ' typot1(3) = ',typot1(3) c write(6,*) ' dgamdh = ',dgamdh c write(6,*) ' typot = ',typot c write(6,*) ' gamma = ',gamma c write(6,*) ' -dtdh (mgals) = ',-dtdh*1d5 c write(6,*) ' dgamdh*typot/gamma (mgals) = ', c * (dgamdh*typot/gamma)*(1d5) write(6,470)w0,u0norm 470 format(1x, * 'Total Gravity Potential on the Geoid (Wo) = ', * f20.8,/,1x, * 'Normal Gravity Potential on the Ellipsoid (Uo) = ', * f20.8) write(6,466)geodlat 466 format(1x,'Ellipsoidal Latitude (dec. deg.) = ',f20.8) write(6,467)geodlon 467 format(1x,'Ellipsoidal Longitude(dec. deg.) = ',f20.8) write(6,468)und 468 format(1x,'Geoid Undulation. (m) = ',f20.8) write(6,469)anom 469 format(1x,'Gravity Anomaly (mgals) = ',f20.8) write(6,473)tidsum 473 format(1x,'Tide system = ',15x,a4,' tide') else h = und goto 435 endif write(6,479) 479 format(1x,'Go again (y/n)? ',$) read(5,'(a)')yesno if(yesno.eq.'y' .or. yesno.eq.'Y')goto 480 stop ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c - FUNCTION # 9 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c - See function #7 for a lot of similar comments elseif(ifun.eq.9) then call bestfit(gmcoef,fomega,abf,j2bf,finvbf,w0,tidbf,tidcoef) call getref(anorm,bnorm,finvnorm,j2norm,gmnorm,fomega,tidsum) fnorm = 1.d0/finvnorm e2norm = 2.d0*fnorm - fnorm*fnorm call norm(anorm,j2norm,gmnorm,fomega,bnorm,jnorm,njmax, * u0norm,gammae,gammap) c - version 0.4: call getdgtyp(idgtyp) c - Now, correct the C(2,0) coefficient for tides, if necessary if(tidsum.ne.tidcoef)then call fixc20(tidcoef,tidsum,c(4),gmcoef,acoef) endif c - At this point, everything is in the "tidsum" tide system. c - Calculate the disturbing coefficients, set up such c - that they work with gmcoef and acoef (which means c - you may need to scale the j2 coefficients from c - working with gmnorm and anorm) c write(6,*) ' total:' c do 9901 iii=0,20 c write(6,*) ' iii,c(iii,0) = ',iii,c(iii*iii) c if(mod(iii,2).eq.0)then c if(iii.eq.0)write(6,*) ' iii,j(iii) = ',iii,1.d0 c if(iii.ne.0)write(6,*) ' iii,j(iii) = ',iii,jnorm(iii/2) c endif c9901 continue call dist(c,max,nmax,gmcoef,acoef,jnorm,njmax,gmnorm,anorm) c write(6,*) ' disturbing:' c do 9902 iii=0,20 c write(6,*) ' iii,c(iii,0) = ',iii,c(iii*iii) c9902 continue c - the "c" vector now contains the fully normalized disturbing c - coefficients c - Get the grid boundaries/spacings call getgrid(glamn,glomn,dla,dlo,nla,nlo,ikind) c - Prompt for file names lgeoid = .false. lganom = .false. write(6,87) 87 format(1x,'Write undulations to a file? :',$) read(5,'(a)')yesno if(yesno.eq.'y' .or. yesno.eq.'Y')then write(6,88) 88 format(1x,'Geoid Undulation file name? :',$) read(5,'(a)')fname open(11,file=fname,status='new',form='unformatted') lgeoid = .true. write(11)glamn,glomn,dla,dlo,nla,nlo,ikind endif write(6,89) 89 format(1x,'Write anomalies to a file? :',$) read(5,'(a)')yesno if(yesno.eq.'y' .or. yesno.eq.'Y')then write(6,90) 90 format(1x,'Gravity Anomaly file name? :',$) read(5,'(a)')fname open(12,file=fname,status='new',form='unformatted') lganom = .true. write(12)glamn,glomn,dla,dlo,nla,nlo,ikind endif c - Loop lat/lon, iterating on h as necessary do 85 irow=1,nla glat = glamn + (irow-1)*dla write(6,*) ' irow,glat = ',irow,glat c - Get the normal gravity at this point and it's deriv wrt h gamma = gam(glat,anorm,bnorm,gammae,gammap) dgamdh = dgdh(glat,anorm,finvnorm,gamma,fomega) do 86 icol = 1,nlo glon = glomn + (icol-1)*dlo c write(6,*) ' icol,glon = ',icol,glon c - Pseudo-Iterate as necessary c - That is, as N changes by more then 10 meters, c - change the evaluation height (h) to N. c - Until that point, +/- 10 meter changes in the c - evaluation point translate to 0.001 m changes c - in the geoid. c - Set the first evaluation height: if(irow.eq.1 .and. icol.eq.1)then c - version 0.3 : moved these 2 lines inside h = 0.d0 call ell2sph(glat,glon,h,anorm,finvnorm,xlat,xlon,r) isw = 0 call pot(xlat,xlon,r,domega,tionpot,tionpot1,tionpot2, * cenpot,cenpot1,cenpot2,typot,typot1,typot2,isw,nmax, * root,smlp1,cmlp1,c,su,max,gmcoef,acoef,force) und = (typot/gamma) - ((w0 - u0norm)/gamma) h = und c - Flags c write(6,*) ' first hit ...' c write(6,*) ' h = ',h c write(6,*) ' gamma = ',gamma c write(6,*) ' typot = ',typot c write(6,*) ' w0 = ',w0 c write(6,*) ' u0 = ',u0norm c write(6,*) ' N1 = ',typot/gamma c write(6,*) ' N2 = ',(w0-u0norm)/gamma c write(6,*) ' N = ',und c write(6,*) ' -------------------' endif c - Do the summation at each point c - All "*pot" values are DISTURBING potential values!!! c - version 0.3: Added this line: call ell2sph(glat,glon,h,anorm,finvnorm,xlat,xlon,r) isw = 1 call pot(xlat,xlon,r,domega,tionpot,tionpot1,tionpot2, * cenpot,cenpot1,cenpot2,typot,typot1,typot2,isw,nmax, * root,smlp1,cmlp1,c,su,max,gmcoef,acoef,force) und = (typot/gamma) - ((w0 - u0norm)/gamma) c - Derivative of Disturbing Potential wrt h c - dT/dh = dT/dNorth * dNorth/dh c + dT/dEast * dEast/dh c + dT/dr * dr/dh c - dEast/dh = 0 and dNorth/dh is assumed zero, though c it is not (but is only about 2mm/m, see Derivation notebook) xw = dsqrt(1.d0 - e2norm*dsin(glat*d2r)**2) xn = anorm/xw drdh=((xn+h)-xn*e2norm*dsin(glat*d2r)**2)/r dtdr = typot1(3) dtdh = dtdr * drdh c - Gravity Anomaly c - H/M, p 86, eqn 2-147c (converted to mgals): c - This formula is more accurate than 2-151c, but is less c - compatable with N=T/gamma + Stokes in the spherical harmonic domain c anom = ( -dtdh + (dgamdh*typot)/gamma )*(1.d5) c - v 0.3: Added the zero order term (eqn 2-186, H/M) c - v 0.3: After studying the problem, it seems I need to keep c gravity anomalies computed with 2nd term of: -2T/r c - A geopot96 type anomaly, eqn 2-151c (converted to mgals): c - This formula is less accurate than 2-147c, but is compatable c - with N=T/gamma + Stokes in the spherical harmonic domain c anom = -dtdh * 1.d5 - 1.d5*2.d0*typot/r c - v 0.3: Added the zero order term (eqn 2-186, H/M) c - v0.4: Allowed a choice for anomaly compuation: if(idgtyp.eq.1)then anom = (-dtdr - 2.d0*typot/r)*1.d5 * + (2*(w0-u0norm)/r)*(1.d5) else anom = ( -dtdh + (dgamdh*typot)/gamma )*(1.d5) * - (dgamdh*(w0-u0norm)/gamma)*(1.d5) endif c write(6,*) ' -------------' c write(6,*) ' for geopot97: ' c write(6,*) ' drdh = ',drdh c write(6,*) ' typot1(3) = ',typot1(3) c write(6,*) ' dgamdh = ',dgamdh c write(6,*) ' typot = ',typot c write(6,*) ' gamma = ',gamma c write(6,*) ' -dtdh (mgals) = ',-dtdh*1d5 c write(6,*) ' dgamdh*typot/gamma (mgals) = ', c * (dgamdh*typot/gamma)*(1d5) c write(6,*) ' ' c write(6,*) ' for geopot96: ' c write(6,*) ' tp (typot) = ',typot c write(6,*) ' r = ',r c write(6,*) ' -2*tp/r (mgals) = ',-1.d5*2.d0*typot/r c write(6,*) ' ' c write(6,*) ' Dif of 2nd terms(mgals) = ', c * 1.d5*( ((dgamdh*typot)/gamma) - (-2.d0*typot/r)) c - From geopot96: c em = (7292115d-11)**2 * anorm**2 * bnorm / gmcoef c glat2 = glat+glat c sglat2 = dsin(glat2*d2r) c eh = em*sglat2*typot1(1) c anom96 = (typot1(3) + 2.d0*typot/r)*(1.d5) + eh*1.d5 c anom96p= (typot1(3)*drdh + 2.d0*typot/r)*(1.d5) c - Flags c write(6,*) ' h = ',h c write(6,*) ' gamma = ',gamma c write(6,*) ' typot = ',typot c write(6,*) ' w0 = ',w0 c write(6,*) ' u0 = ',u0norm c write(6,*) ' N1 = ',typot/gamma c write(6,*) ' N2 = ',(w0-u0norm)/gamma c write(6,*) ' N = ',und c write(6,*) ' ----------' c write(6,*) ' xn = ',xn c write(6,*) ' dtdr = ',dtdr c write(6,*) ' drdh = ',drdh c write(6,*) ' dtdh = ',dtdh c write(6,*) ' dgamdh=',dgamdh c write(6,*) ' Dg1 = ',-dtdh * (1.d5) c write(6,*) ' Dg2 = ',(1.d5)*(dgamdh*typot)/gamma c write(6,*) ' Dg= ',anom c write(6,*) ' -------------------' c write(6,*) ' typot1(1) = ',typot1(1) c write(6,*) ' sglat2 = ',sglat2 c write(6,*) ' em = ',em c write(6,*) ' Dg1 96 = ',typot1(3)*1.d5 c write(6,*) ' Dg2 96 = ',1.d5*2.d0*typot/r c write(6,*) ' Dg3 96 = ',eh*1.d5 c write(6,*) ' Dg 96(old) = ',anom96 c write(6,*) ' -----------------------' c write(6,*) ' Dg 96(new) = ',anom96p c - v 0.3: Moved these lines to here: dd = dabs(h-und) if(dd.gt.10.d0)then ! <- 10 meters works well h = und endif c - Output to files and loop back if(lgeoid)geoid(icol) = und if(lganom)ganom(icol) = anom 86 continue if(lgeoid)write(11)(geoid(j),j=1,nlo) if(lganom)write(12)(ganom(j),j=1,nlo) 85 continue stop ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c - FUNCTION # 10 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c - See function #7 for a lot of similar comments elseif(ifun.eq.10) then call bestfit(gmcoef,fomega,abf,j2bf,finvbf,w0,tidbf,tidcoef) call getref(anorm,bnorm,finvnorm,j2norm,gmnorm,fomega,tidsum) fnorm = 1.d0/finvnorm e2norm = 2.d0*fnorm - fnorm*fnorm call norm(anorm,j2norm,gmnorm,fomega,bnorm,jnorm,njmax, * u0norm,gammae,gammap) c - Now, correct the C(2,0) coefficient for tides, if necessary if(tidsum.ne.tidcoef)then call fixc20(tidcoef,tidsum,c(4),gmcoef,acoef) endif c - At this point, everything is in the "tidsum" tide system. c - Calculate the disturbing coefficients, set up such c - that they work with gmcoef and acoef (which means c - you may need to scale the j2 coefficients from c - working with gmnorm and anorm) write(6,*) ' total:' do 9903 iii=0,20 write(6,*) ' iii,c(iii,0) = ',iii,c(iii*iii) if(mod(iii,2).eq.0)then if(iii.eq.0)write(6,*) ' iii,j(iii) = ',iii,1.d0 if(iii.ne.0)write(6,*) ' iii,j(iii) = ',iii,jnorm(iii/2) endif 9903 continue call dist(c,max,nmax,gmcoef,acoef,jnorm,njmax,gmnorm,anorm) c write(6,*) ' disturbing:' c do 9904 iii=0,20 c write(6,*) ' iii,c(iii,0) = ',iii,c(iii*iii) c9904 continue c - the "c" vector now contains the fully normalized disturbing c - coefficients c - Get the grid boundaries/spacings call getgrid(glamn,glomn,dla,dlo,nla,nlo,ikind) c - Prompt for file names lt0 = .false. lt1 = .false. lt2 = .false. write(6,95) 95 format(1x,'Write height anomalies to a file? :',$) read(5,'(a)')yesno if(yesno.eq.'y' .or. yesno.eq.'Y')then write(6,96) 96 format(1x,'Ht Anom on Ellipsoid file name? :',$) read(5,'(a)')fname open(21,file=fname,status='new',form='unformatted') lt0 = .true. write(21)glamn,glomn,dla,dlo,nla,nlo,ikind endif write(6,97) 97 format(1x,'Write 1st upward derivs (dXi/dh) to a file? :',$) read(5,'(a)')yesno if(yesno.eq.'y' .or. yesno.eq.'Y')then write(6,98) 98 format(1x,'1st deriv file name? :',$) read(5,'(a)')fname open(22,file=fname,status='new',form='unformatted') lt1 = .true. write(22)glamn,glomn,dla,dlo,nla,nlo,ikind endif write(6,99) 99 format(1x,'Write 2nd upward derivs (d2Xi/dh2) to a file? :',$) read(5,'(a)')yesno if(yesno.eq.'y' .or. yesno.eq.'Y')then write(6,100) 100 format(1x,'2nd deriv file name? :',$) read(5,'(a)')fname open(23,file=fname,status='new',form='unformatted') lt2 = .true. write(23)glamn,glomn,dla,dlo,nla,nlo,ikind endif c - Loop lat/lon, computing everything on the ellipsoid do 101 irow=1,nla glat = glamn + (irow-1)*dla write(6,*) ' irow,glat = ',irow,glat c - Get the normal gravity at this point and it's derivs wrt h gamma = gam(glat,anorm,bnorm,gammae,gammap) dgamdh = dgdh(glat,anorm,finvnorm,gamma,fomega) d2gamdh2=d2gdh2(gamma,anorm) do 102 icol = 1,nlo glon = glomn + (icol-1)*dlo write(6,*) ' icol,glon = ',icol,glon c - h set to zero...we evaluate everything on the ellipsoid surface h = 0.d0 call ell2sph(glat,glon,h,anorm,finvnorm,xlat,xlon,r) c - Do the summation at each point c - All "*pot" values are DISTURBING potential values!!! isw = 2 call pot(xlat,xlon,r,domega,tionpot,tionpot1,tionpot2, * cenpot,cenpot1,cenpot2,typot,typot1,typot2,isw,nmax, * root,smlp1,cmlp1,c,su,max,gmcoef,acoef,force) c - Height anomaly on ellipsoid: c - We want to have a height anomaly and its derivatives such c - that we can upward continue them, get a height anomaly at c - the surface, and downward correct (using Bouguer) to get c - an undulation between W=W0 and U=U0norm, where W0<>U0norm. c - As such, we have to include the mis-matched potential c - right in the height anomaly (but being nearly constant, it c - won't enter into the height derivatives). This means that c - the height anomaly at the surface will *NOT* be the separation c - of W=U, but will rather have this constant term in there. htanom = (typot/gamma) - ((w0 - u0norm)/gamma) c - Derivative of Disturbing Potential wrt h c - dT/dh = dT/dNorth * dNorth/dh c + dT/dEast * dEast/dh c + dT/dr * dr/dh c - dEast/dh = 0 and dNorth/dh is assumed zero, though c it is not (but is only about 2mm/m, see Derivation notebook) xw = dsqrt(1.d0 - e2norm*dsin(glat*d2r)**2) xn = anorm/xw drdh=((xn+h)-xn*e2norm*dsin(glat*d2r)**2)/r dtdr = typot1(3) dtdh = dtdr * drdh c - 1st derivative of ht anomaly wrt h t1a = dtdh/gamma t1b = -typot*dgamdh/gamma/gamma c - 2nd derivative of ht anomaly wrt h d2tdr2 = typot2(3,3) d2tdh2 = d2tdr2 * (drdh)*(drdh) t2a = d2tdh2/gamma t2b = -typot*(d2gamdh2/gamma/gamma) t2c = -typot*(-2.d0*dgamdh*dgamdh/gamma/gamma/gamma) write(6,*) ' h = ',h write(6,*) ' htanom = ',htanom write(6,*) ' 1st der= ',t1a+t1b write(6,*) ' 2nd der= ',t2a+t2b+t2c c - Output to files and loop back if(lt0)t0(icol) = htanom if(lt1)t1(icol) = t1a+t1b if(lt2)t2(icol) = t2a+t2b+t2c 102 continue if(lt0)write(21)(t0(j),j=1,nlo) if(lt1)write(22)(t1(j),j=1,nlo) if(lt2)write(23)(t2(j),j=1,nlo) 101 continue stop ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c - FUNCTION # 11 - New as of version 0.2 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc elseif(ifun.eq.11) then c - !!!!!!!NOTE!!!!!!! c - By their very nature, DTED's are notoriously inaccurate, c - at least MUCH less accurate than the 10-20 cm level associated c - with the different tide systems. As such, tide systems c - will be ignored in this section. c - !!!!!!!NOTE!!!!!!! c - Get the DTED info write(6,163) 163 format * (/1x,'Input the name of the height DTED: ',$) read(5,'(a)')fname open(41,file=fname,status='old',form='unformatted') 171 write(6,170) 170 format(/1x,'1 = ellipsoidal',/,1x,'2 = orthometric',/, * 1x,'What sort of DTED? : ',$) read(5,*)idtyp if(idtyp.lt.1 .or. idtyp.gt.2)goto 171 write(6,164) 164 format(1x,'Does the header lie? : ',$) read(5,'(a)')yesno xlatlie = 0.d0 xlonlie = 0.d0 if(yesno.eq.'y' .or. yesno.eq.'Y')then write(6,165) 165 format(1x,'Points in DTED are too far south ', * 'by how many seconds?: ',$) read(5,*)xlatlie write(6,166) 166 format(1x,'Points in DTED are too far west ', * 'by how many seconds?: ',$) read(5,*)xlonlie endif read(41)glamn,glomn,dla,dlo,nla,nlo,ikind glamn = glamn + xlatlie/3600.d0 glomn = glomn + xlonlie/3600.d0 glamx = glamn + (nla-1)*dla glomx = glomn + (nlo-1)*dlo c - Get the output info write(6,*) write(6,*) 'kind=',ikind write(6,*) 'LAT min=',glamn write(6,*) ' del=',dla,' # lat=',nla write(6,*) ' max=',glamn+(nla-1)*dla write(6,*) 'LON min=',glomn write(6,*) ' del=',dlo,' # lon=',nlo write(6,*) ' max=',glomn+(nlo-1)*dlo write(6,*) ' ' write(6,*) ' The output can be on a decimated version of the' write(6,*) ' above grid spacing. So...' write(6,167) 167 format(1x,'Output on 1 of every "n" points in lat: ',$) read(5,*) ny write(6,168) 168 format(1x,'Output on 1 of every "n" points in lon: ',$) read(5,*) nx dla2 = ny*dla dlo2 = nx*dlo nla2 = (nla-1)/ny+1 nlo2=(nlo-1)/nx+1 write(6,181) 181 format(1x,'Skip zero-elevation points?: ',$) read(5,'(a)')yesno iskz = 0 if(yesno.eq.'y' .or. yesno.eq.'Y')iskz = 1 172 write(6,169) 169 format(/1x,'Output types: ',/, * 1x,'1 = Gridded ".b" format',/, * 1x,'2 = ASCII lat/lon/H/g file',/, * 1x,'Which type do you want? : ',$) read(5,*)iouttyp if(iouttyp.lt.1 .or. iouttyp.gt.2)goto 172 write(6,173) 173 format(/1x,'Output file name? :',$) read(5,'(a)')fname if(iouttyp.eq.1)then open(42,file=fname,status='new',form='unformatted') else open(42,file=fname,status='new',form='formatted') endif c - Find the geometric ellipsoid (for idtyp=1,iouttyp=1) or c - the normal field (for idtyp=1,iouttyp=2 and idtyp=2) if(idtyp.eq.1 .and. iouttyp.eq.1)then call getgeom(anorm,finvnorm) else write(6,178) 178 format(/1x,'Because of the combination of input/output' * 1x,'you have requested, I need to approximate the' * 1x,'geoid undulation, and for that I need to ask' * 1x,'some questions. Because of the inaccuracy' * 1x,'of DTEDs, and the small magnitude of tide' * 1x,'system signals, the tide system questions you' * 1x,'will get are irrelevant. Just use the tide' * 1x,'system of the coefficients as an answer to any' * 1x,'of these tide questions. Let us proceed....',/) c - Define the Geoid call bestfit(gmcoef,fomega,abf,j2bf,finvbf,w0,tidbf,tidcoef) c - Get the bias to add to the DTED to make H refer to best global MSL dorth = 0.d0 call getmsl(dorth,b88) c - Get the normal field call getref(anorm,bnorm,finvnorm,j2norm,gmnorm,fomega,tidsum) fnorm = 1.d0/finvnorm e2norm = 2.d0*fnorm - fnorm*fnorm c - Get the normal field parameters call norm(anorm,j2norm,gmnorm,fomega,bnorm,jnorm,njmax, * u0norm,gammae,gammap) c - Now, correct the C(2,0) coefficient for tides, if necessary if(tidsum.ne.tidcoef)then call fixc20(tidcoef,tidsum,c(4),gmcoef,acoef) endif c - At this point, everything is in the "tidsum" tide system. endif c - Do the calculations do 174 irow = 1,nla read(41)(elev(j),j=1,nlo) iii = irow-1 if(mod(iii,nx).ne.0)goto 174 ! <- Decimates by row glat = glamn + dla*(irow-1) write(6,*) ' Latitude of this row = ',glat do 175 icol = 1,nlo2 jorig = (icol-1)*ny+1 ! <- Decimates by column glon = glomn + dlo2*(icol-1) ell = 0.d0 orth = 0.d0 if(idtyp.eq.1)ell = elev(jorig) if(idtyp.eq.2)then if(elev(jorig).ne.0.d0)orth = elev(jorig)+dorth if(elev(jorig).eq.0.d0)orth = elev(jorig) endif c - Skip zero points if requested if(idtyp.eq.2 .and. iskz.eq.1 .and. orth.eq.0.d0)goto 175 c - If: c - Given H, want phi/lam/H/g => Need N first to get h for g c - Given H, want gridded => Need N first to get h for g c - Given h, want phi/lam/H/g => Need N first to get H for list c - Given h, want gridded => Don't need N if(idtyp.eq.1 .and. iouttyp.eq.1)goto 176 c - Get Spherical coordinates of point on ellipsoid: call ell2sph(glat,glon,0.d0,anorm,finvnorm,xlat,xlon,r) c - Now get gamma at the ellipsoid gamma = gam(glat,anorm,bnorm,gammae,gammap) c do 8102 izz = 0,15 c write(6,*) ' i,c(i)(orig) = ',i,c(izz) c8102 continue c - Disturb the coefficients call dist(c,max,nmax,gmcoef,acoef,jnorm,njmax,gmnorm,anorm) c - Then call pot to get T(h=0) for a quick'n'dirty N computation isw = 0 call pot(xlat,xlon,r,domega,tionpot,tionpot1,tionpot2, * cenpot,cenpot1,cenpot2,typot,typot1,typot2,isw,nmax, * root,smlp1,cmlp1,c,su,max,gmcoef,acoef,force) und0 = tionpot/gamma - (w0 - u0norm)/gamma c - Then use the quick'n'dirty N to connect h to H if(idtyp.eq.1)then orth = ell - und0 ! <-no dorth, cuz we force orht to msl else if(orth.ne.0.d0)then ell = (orth+dorth) + und0 ! <-dorth cuz we're given orth else ell = (orth ) + und0 ! <- oceans have no bias endif endif c - Lastly, un-disturb the coefficients call undist(c,max,nmax,gmcoef,acoef,jnorm,njmax,gmnorm * ,anorm) c do 8101 izz = 0,15 c write(6,*) ' i,c(i)(undist) = ',i,c(izz) c8101 continue c - Now call pot to get the gravity value at the surface "h" c - But first get the spherical coordinates of the surface pt call ell2sph(glat,glon,ell,anorm,finvnorm,xlat,xlon,r) 176 isw=1 call pot(xlat,xlon,r,fomega,tionpot,tionpot1,tionpot2, * cenpot,cenpot1,cenpot2,typot,typot1,typot2,isw,nmax, * root,smlp1,cmlp1,c,su,max,gmcoef,acoef,force) gravityh = dsqrt(typot1(1)**2+typot1(2)**2+typot1(3)**2) c write(6,*) ' *********************************************' c write(6,*) ' xlat,xlon,r = ',xlat,xlon,r c write(6,*) ' glat,glon = ',glat,glon c write(6,*) ' h = ',ell c write(6,*) ' H = ',ell-und0-dorth c write(6,*) ' gravityh(mgals) = ',gravityh*1.d5 c write(6,177) glat,glon,orth,gravityh*1d5 if(iouttyp.eq.1)outg(icol)=gravityh*1d5 c - Loop and Output (in mgals) if(iouttyp.eq.2)then write(42,177)glat,glon,orth,gravityh*1d5 177 format(f13.7,1x,f13.7,1x,f8.3,1x,f10.3) endif 175 continue if(iouttyp.eq.1)then write(42)(outg(j),j=1,nlo2) endif 174 continue endif end c c c subroutine pot(xlat,xlon,r,omega,tionpot,tionpot1,tionpot2, * cenpot,cenpot1,cenpot2,typot,typot1,typot2,isw,nmax, * root,smlp1,cmlp1,c,su,max,gmsum,asum,force) c c - subroutine to compute quantities based on either the FULL c - potential, or the DISTURBING potential. For the case of the c - disturbing potential, it is assumed that scaling of the c - coefficients has occurred if differing "GM" or "a" values c - are being sent in, and that the coefficients (and also omega!!) c - are DISTURBING values, not FULL values!!! c c - tionpot = gravitaTIONal POTential (full or disturbing, based on coefs) c - typot = graviTY POTential (full or disturbing, based on coefs & omega) c - cenpot = CENtrifugal POTential (full or disturbing, based on omega) c c - tionpot1,typot1,cenpot1 = 1st order derivatives (N,E,U) c - tionpot2,typot2,cenpot2 = 2nd order derivatives c - (NN,NE,NU ; EN,EE,EU ; UN,UE,UU) c c - isw = iSWitch, to tell the summing routine whether we're c looking for: c 0 = just the potential c 1 = potential and 1st derirvatives (gravity) c 2 = potential, 1st and 2nd derivatives (gravity gradients) implicit real*8(a-h,o-z) real*8 root((max+1)*2) real*8 smlp1(max+2),cmlp1(max+2),c(0:(max+1)*(max+1)-1) real*8 su(10*(max+1)) real*8 tionpot1(3),cenpot1(3),typot1(3) real*8 tionpot2(3,3),cenpot2(3,3),typot2(3,3) logical force c write(6,*) ' inside pot' call centpot(xlat,xlon,r,omega,cenpot,cenpot1,cenpot2,isw) c - cenpot = centrifugal potential c - cenpot1 = 1st derivatives, 3 directions c - cenpot2 = 2nd derivatives, 3x3 directions call addemup(xlat,xlon,r,tionpot,tionpot1,tionpot2,isw, * nmax,su, * root,smlp1,cmlp1,c,max,gmsum,asum,force) c - tionpot = gravitational potential c - tionpot1 = 1st derivatives, 3 directions c - tionpot2 = 2nd derivatives, 3x3 directions typot = cenpot + tionpot typot1(1) = cenpot1(1) + tionpot1(1) ! North typot1(2) = cenpot1(2) + tionpot1(2) ! East typot1(3) = cenpot1(3) + tionpot1(3) ! radial(up) typot2(1,1) = cenpot2(1,1) + tionpot2(1,1) ! N,N typot2(1,2) = cenpot2(1,2) + tionpot2(1,2) ! N,E = E,N typot2(1,3) = cenpot2(1,3) + tionpot2(1,3) ! N,r = r,N typot2(2,1) = cenpot2(2,1) + tionpot2(2,1) ! E,N = N,E typot2(2,2) = cenpot2(2,2) + tionpot2(2,2) ! E,E typot2(2,3) = cenpot2(2,3) + tionpot2(2,3) ! E,r = r,E typot2(3,1) = cenpot2(3,1) + tionpot2(3,1) ! r,N = N,r typot2(3,2) = cenpot2(3,2) + tionpot2(3,2) ! r,E = E,r typot2(3,3) = cenpot2(3,3) + tionpot2(3,3) ! r,r return end c c c subroutine centpot(xlat,xlon,r,omega,cp,cp1,cp2,isw) c - subroutine to compute centrifugal potential, and its c - 1st and 2nd order derivatives. implicit real*8(a-h,o-z) real*8 cp1(3),cp2(3,3) common/block1/pi,d2r c write(6,*) ' inside centpot' c - Clear anything out cp = 0.d0 do 1 i=1,3 cp1(i) = 0.d0 do 2 j=1,3 cp2(i,j) = 0.d0 2 continue 1 continue c = dcos(xlat*d2r) cp = 0.5d0 * omega*omega * r*r *c*c if(isw.eq.0)return s = dsin(xlat*d2r) cp1(1) = -omega*omega*r*c*s cp1(2) = 0.d0 cp1(3) = omega*omega*r*c*c if(isw.eq.1)return c2 = dcos(2*xlat*d2r) cp2(1,1) = -omega*omega*c2 cp2(1,2) = 0.d0 cp2(1,3) = -2*omega*omega*c*s cp2(2,1) = cp2(1,2) cp2(2,2) = 0.d0 cp2(2,3) = 0.d0 cp2(3,1) = cp2(1,3) cp2(3,2) = cp2(2,3) cp2(3,3) = omega*omega*c*c return end c c c subroutine addemup(xlat,xlon,r,tionpot,tionpot1,tionpot2,isw, * nmax,su, * root,smlp1,cmlp1,c,max,gmsum,asum,force) c - subroutine to prepare to go into GPTDR, the Clenshaw c - summation subroutine. c c - tionpot = Gravitational Potential c - tionpot1 = 1st derivs of tionpot c - tionpot2 = 2nd derivs of tionpot c c - isw = switch (0,1 or 2) telling whether you want 0th, 1st c - or 2nd order derivatives added up c implicit real*8(a-h,o-z) real*8 root((max+1)*2) real*8 smlp1(max+2),cmlp1(max+2),c(0:(max+1)*(max+1)-1) real*8 su(10*(max+1)) real*8 tionpot1(3),tionpot2(3,3) real*8 po(6) integer*4 order logical force common/block1/pi,d2r c write(6,*)' inside addemup' c - Clear anything out tionpot = 0.d0 do 1 i=1,3 tionpot1(i) = 0.d0 do 2 j=1,3 tionpot2(i,j) = 0.d0 2 continue 1 continue po(1) = r * dcos(xlat*d2r) po(2) = r po(3) = dcos( (90.d0 - xlat) * d2r) po(4) = dsin( (90.d0 - xlat) * d2r) po(5) = dsin( xlon * d2r) po(6) = dcos( xlon * d2r) c - call gptdr with negative nmax to indicate the use of c - quasi-normalized coefficients order = isw negn = -nmax tionpot = gptdr(po,negn,order,su,tionpot1,tionpot2, * root,smlp1,cmlp1,c,max,gmsum,asum,force) return end c c c SUBROUTINE SETCM(CAPN,root,smlp1,cmlp1,c,su,max) *** THIS ROUTINE SETS THE SQUARE ROOT TABLE IN 'root' *** changes the coefficients from fully normalized to quasi-normalized c - CAPN = actual maximum degree to be used c - max = maximal possible degree read in c - Fully normalized and quasi-normalized C(0,0) is the same IMPLICIT REAL*8(A-H,O-Z),INTEGER*4(I-N) INTEGER*4 CAPN real*8 root((max+1)*2) real*8 smlp1(max+2),cmlp1(max+2),c(0:(max+1)*(max+1)-1) real*8 su(10*(max+1)) c write(6,*) ' inside setcm' DO 50 I=1,2*(max + 1) X=I ROOT(I)=DSQRT(X) 50 continue SQ2=DSQRT(2.D0) ! <- to change the m=0 conversion to m<>0 conversion c - Don't loop from 0, since deg 0 fully norm and quasi norm are same DO 200 N=1,CAPN N2=N+N c S21=DSQRT(N2+1.D0) ! <- stupid to re-calculate! S21 = root(n2+1) ! <- use pre-calculated values!!! K=N**2 ! <- neat...c(n**2) = C(n,m=0) *** D IS THE QUASI-NORMALIZATION FACTOR FOR ZONAL TERMS D=S21 C(K)=C(K)*D *** GG IS THE QUASI-NORMALIZATION FACTOR FOR NON-ZONAL TERMS GG=D*SQ2 DO 100 J=1,N KJ2=J+J+K C(KJ2-1)=C(KJ2-1)*GG ! <- the C(n,m<>0) terms C(KJ2)=C(KJ2)*GG ! <- the S(n,m<>0) terms 100 CONTINUE 200 CONTINUE RETURN END c c ------------------------------------------------------------------- c double precision function gptdr(po,nmax,order,su,g1,g2, * root,smlp1,cmlp1,c,max, * gmsum,asum,force) c - Version 1/14/97, by Dru Smith, National Geodetic Survey c - Uses the core of the original function, but cleans it up c - for readability and moved all non-clenshaw functions c - OUTSIDE this function. c - Added logical function "force" which, if true, forces c - a re-summation (useful when you've switched the coefficients) c - Rather than have commons/parameters, I just pass all the c - variable addresses: g1,g2,root,smlp1,cmlp1,c,max c - I also pass ONLY ONE value of GM (gmsum) and c - ONLY ONE value of "a" (asum), assuming that all of the c - coefficients are scaled. I also do NOT pass omega, as it has c - no bearing on the gravitational potential (which is this c - subroutine) and do NOT pass J2 (or other parts of the reference c - normal field) as they also have nothing to do with the c - Clenshaw summation of gravitational potential! c c - IN short, all I need (variable-wise) is: c 1) The "c" array holding C(n,m) and S(n,m) in a set c order, and already quasi-normalized. These coefficients c represent either full or disturbing gravitational c potential, and are scaled to work with "gmsum" and "asum" c 2) The value of GM used in the summation, gmsum c 3) The value of a used in the summation, asum c 4) The positional information, in the "po" array c 5) The maximum degree to take the summation, nmax c 6) The highest order derivatives to compute, order. c - What comes out: c 1) GPTDR itself returns the gravitational potential (full or c disturbing. c 2) g1(3) returns the N,E and up derivatives of GPTDR c 3) g2(3,3) returns the NN,NE,NU,EN,EE,EU,UN,UE,UU derivatives c of GPTDR c - THAT IS ALL THIS SUBROUTINE HAS BEEN MODIFIED TO DO! c - Original comments: *** GI REG.NO. 81013 AUTHOR -C.C.TSCHERNING, DANISH GEODETIC INSTITUTE *** JULY 1981 IN ALGOL REF.(2) *** -C.C.GOAD, NOAA/NOS/NATIONAL GEODETIC SURVEY *** MAY 1982 TRANSLATED TO FORTRAN *** REFERENCES: *** (1) TSCHERNING, C.C.:ON THE CHAIN-RULE METHOD FOR COMPUTING *** POTENTIAL DERIVATIVES. MANUSCRIPTA GEODAETICA, VOL.1, *** PP. 125-141, 1976 *** (2) TSCHERNING, C.C., AND PODER, K.: SOME APPLICATIONS OF CLENSHAW *** SUMMATION, PRESENTED AT VIII SYMPOSIUM ON MATHEMATICAL GEODESY, *** COMO, ITALY, SEPT 7-9, 1981 *** THE PROCEDURE COMPUTES THE VALUE AND UP TO THE SECOND-ORDER *** DERIVATIVES OF THE POTENTIAL OF THE EARTH (W) OR OF ITS *** CORRESPONDING ANOMALOUS POTENTIAL(T). *** THE POTENTIAL IS REPRESENTED BY A SERIES OF SOLID SPHERICAL *** HARMONICS, WITH UN-NORMALIZED OR QUASI-NORMALIZED COEFFICIENTS. *** THE CHAIN-RULE IS USED ALONG WITH THE CLENSHAW ALGORITHM. *** THE ARRAY C MUST HOLD THE COEFFICIENTS C(1)=C(1,0),C(2)=C(1,1), *** C(3)=S(1,1), ETC. UP TO C((N+1)**2-1=S(N,N). C(0,0) IS STORED IN C0 *** WHICH IMPLICITLY ACTS AS C(0) (SEE THE COMMON BLOCK CM). *** PARAMETERS: *** (A) INPUT VALUES: *** NMAX *** THE ABSOLUTE VALUE OF NMAX IS EQUAL TO THE MAXIMAL DEGREE AND *** ORDER OF THE SERIES. NEGATIVE NMAX INDICATES THAT THE COEFFICIENTS *** ARE QUASI-NORMALIZED. *** ORDER *** ORDER OF DERIVATIVES *** 0 FOR POTENTIAL ONLY *** 1 FOR POTENTIAL AND FIRST DERIVATIVES *** 2 FOR POTENTIAL, FIRST DERIVATIVES, AND SECOND DERIVATIVES *** PO *** ARRAY HOLDING POSITION INFORMATION. PO(6) *** PO(1)=P, THE DISTANCE FROM THE Z (ROTATION) AXIS, *** PO(2)=R, THE DISTANCE FROM THE ORIGIN, *** PO(3),PO(4) COS AND SIN OF GEOCENTRIC POLAR ANGLE(COLATITUDE), *** PO(5),PO(6) SIN AND COS OF THE LONGITUDE. *** C *** C((ABS(NMAX)+1)**2-1) ARRAY OF C'S AND S'S DESCRIBED ABOVE c - GEOPOT97: Ignore the following historical comments. c - das *** CM3=GM *** CM2=A THE SEMI-MAJOR AXIS OF THE REFERENCE ELLIPSOID *** CM1=THE ANGULAR VELOCITY (=0,WHEN DEALING WITH T) *** C0=1.D0 FOR W AND =0.D0 FOR T c - ACTUALLY!!!!! -- CM1,2,3 are truly those values assigned c - to the UNSCALED coeffients. However, c - the coefficients have been SCALED, so c - that the proper use of the scale factors c - is to use OMEGA, AE, and GM!!!!!!!!! c - These values come in under the "normal" c - common block, and are hard-coded into c - this subroutine, replacing CM1,CM2,CM3 c - ALSO - C0 is not necessarily 1 or 0. It c - is (CM3/GM)-1, for T, and 1 for W. c - This is because of the c - possible difference between CM3 and GM c - das c - GEOPOT97: Til here for ignoring *** ROOT(K)=SQRT(K), 0.LE.K.LE.2(ABS(N)+1)-1 WHEN NMAX.LT.0 *** (B) RETURN VALUES: *** G1 AND G2 *** THE RESULT IS STORED IN G1 AND G2 AS FOLLOWS: *** G1(1)=DW/DX, G1(2)=DW/DY, G1(3)=DW/DZ *** G2(1,1)=DDW/DXX, G2(1,2)=G2(2,1)=DDW/DXDY, *** G2(1,3)=G2(3,1)=DDW/DXDZ, G2(2,2)=DDW/DYY, *** G2(2,3)=G2(3,2)=DDW/DYDZ AND G2(3,3)=DDW/DZZ *** WHERE W MAY BE INTERCHANGED WITH T AND *** VARIABLES X, Y, Z ARE THE CARTESIAN COORDINATES *** IN A LOCAL (FIXED) FRAME WITH ORIGIN IN THE POINT *** OF EVALUATION, X POSITIVE NORTH, Y POSITIVE EAST, *** AND Z POSITIVE IN THE DIRECTION OF THE RADIUS *** VECTOR, (CF. REF.(1),EQ (4) AND (5)). *** THE VALUES OF W OR T WILL BE RETURNED IN GPTDR. *** (C) PASSED AND RETURNED VALUES: *** SU *** ARRAY OF DIMENSION K*(N+1), WHERE K=2 FOR NO DERIVATIVES, *** =6 FOR 0-TH AND FIRST DERIVATIVES, =10 FOR 0-TH, FIRST AND *** SECOND DERIVATIVES. HERE ARE STORED THE PARTIAL SUMS, CF. *** REF.(2), EQ. (29), OF P(N,M)*(A/R)**(N+1-M)/P(M,M)*(C(N,M) OR *** S(N,M)) FROM N=M TO N=N, AND THE DERIVATIVES OF THESE SUMS. *** THIS MAKES IT UNNECESSARY TO RECALCULATE THESE QUANTITIES, IF *** THE PROCEDURE IS CALLED SUBSEQUENTLY WITH THE SAME VALUE OF T *** AND R, AND THE SAME ORDER. ********************************************************************** IMPLICIT REAL*8 (A-H,O-Z),INTEGER*4(I-N) INTEGER*4 CAPN,ORDER,CAPN21,OLDORD LOGICAL QUASI,DERIV1,DERIV2,POLE LOGICAL FIRST,NEW,OLD,NPOLE logical force REAL*8 M21,M21T,M21U,M21U0 real*8 root((max+1)*2) real*8 smlp1(max+2),cmlp1(max+2),c(0:(max+1)*(max+1)-1) real*8 su(10*(max+1)) real*8 po(6),g1(3),g2(3,3) c - These determine if we've been here before (efficiency!) COMMON/GPTCM/OLDT,OLDR,IZ,FIRST,OLDORD,I1,I2,I3,I4, * I5,I6,I7,I8,I9,NMAXSV c - das: Removed this (can't have "adjustables" in equivalences) c - das: and entirely removed reference to SML and CML c EQUIVALENCE(SML(1),SMLP1(2)),(CML(1),CMLP1(2)) c write(6,*) ' inside gptdr' c do 8001 i=0,20 c write(6,*) ' inside gptdr, i,c(i) = ',i,c(i) c8001 continue c - convert to the local variable to avoid a lengthy replacement CM2 = asum CM3 = gmsum IF(NMAXSV.NE.NMAX)FIRST=.FALSE. NMAXSV=NMAX if(force) first = .false. ! <-seems backwards, but logically ok IF(FIRST)GO TO 100 FIRST=.TRUE. OLDT=2.D0 ! < - picked to force disagreement with -1<=T<=1 J=IABS(NMAX) I=J+1 I1=I+1 I2=I1+I I3=I2+I I4=I3+I I5=I4+I I6=I5+I I7=I6+I I8=I7+I I9=I8+I 100 CAPN=NMAX *** DISTANCE FROM ROTATION AXIS P=PO(1) *** DISTANCE FROM ORIGIN R=PO(2) *** COSINE OF COLATITUDE T=PO(3) *** SINE OF COLATITUDE U=PO(4) *** SINE OF LONGITUDE SL=PO(5) *** COSINE OF LONGITUDE CL=PO(6) T2=T+T POLE=DABS(U).LE.1.D-9 NEW=DABS(OLDR-R).GT.1.D-3 .OR. DABS(OLDT-T).GT.1.D-9 .OR. * OLDORD.NE.ORDER .OR. POLE OLD=.NOT.NEW NPOLE=.NOT.POLE IF(OLD)GO TO 200 OLDR=R OLDT=T OLDORD=ORDER 200 QUASI=.FALSE. IF(CAPN.LT.0)QUASI=.TRUE. IF(QUASI)CAPN=-CAPN *** COMPUTE AE/R S=CM2/R S2=S**2 CMLP1(1)=1.D0 *** CML(0)=1.D0 SMLP1(1)=0.D0 *** SML(0)=0.D0 DERIV1=.FALSE. IF(ORDER.GT.0)DERIV1=.TRUE. DERIV2=.FALSE. IF(ORDER.GT.1) DERIV2=.TRUE. c - das: SMLP1(M) = sin( (m-1)*longitude) c - das: CMLP1(M) = cos( (m-1)*longitude) c - das: SML(M) and CML(M) removed entirely *** SML(M) AND CML(M) ARE THE SINE AND COSINE OF M*LONGITUDE c SML(1)=SL c CML(1)=CL SMLP1(2)=SL CMLP1(2)=CL c M1=1 c DO 300 M=2,CAPN c SML(M)=SML(M1)*CL+CML(M1)*SL c CML(M)=CML(M1)*CL-SML(M1)*SL c 300 M1=M M1=2 DO 300 M=3,CAPN SMLP1(M)=SMLP1(M1)*CL+CMLP1(M1)*SL CMLP1(M)=CMLP1(M1)*CL-SMLP1(M1)*SL 300 M1=M CAPN21=CAPN+CAPN+1 VM=0.D0 VXM=0.D0 VYM=0.D0 VZM=0.D0 SQNM1=1.D0 SQNPM1=1.D0 IF(.NOT.DERIV2) GO TO 400 VXXM=0.D0 VYYM=0.D0 VZZM=0.D0 VXYM=0.D0 VXZM=0.D0 VYZM=0.D0 400 KM=(CAPN+1)**2 MAX2=CAPN21 *** WE NOW USE THE CLENSHAW ALGORITHM, CF. REF.(2), EQ(8-13), *** MODIFIED IN AN OBVIOUS WAY FOLLOWING REF.(1). ITWO=2 c write(6,*) ' po(1) = ',po(1) c write(6,*) ' po(2) = ',po(2) c write(6,*) ' po(3) = ',po(3) c write(6,*) ' po(4) = ',po(4) c write(6,*) ' po(5) = ',po(5) c write(6,*) ' po(6) = ',po(6) c do 8801 i=0,20 c write(6,*) ' i,c(i) = ',i,c(i) c8801 continue DO 1700 IM=IZ,CAPN M=CAPN-IM MPLUS1=M+1 IF(M.EQ.0)ITWO=1 KM=KM-ITWO K=KM N21=CAPN21 VS=0.D0 VC=0.D0 VS1=0.D0 VC1=0.D0 VXS1=0.D0 VXC1=0.D0 VZS=0.D0 VZC=0.D0 VZS1=0.D0 VZC1=0.D0 VXC=0.D0 VXS=0.D0 IF(.NOT.DERIV2)GO TO 500 VXXC=0.D0 VXXS=0.D0 VXXC1=0.D0 VXXS1=0.D0 VZZC=0.D0 VZZS=0.D0 VZZC1=0.D0 VZZS1=0.D0 VXZC=0.D0 VXZS=0.D0 VXZC1=0.D0 VXZS1=0.D0 500 CM=CMLP1(MPLUS1) SM=SMLP1(MPLUS1) NM1=CAPN-M+2 N1=CAPN+1 NPM1=CAPN+M+2 IF(DERIV2)M2=M*M IF(OLD)GO TO 1300 N=CAPN+1 DO 1000 IN=M,CAPN N=N-1 NM2=NM1 NM1=NM1-1 NPM1=NPM1-1 *** REF.(2) EQ.(40) IF(.NOT.QUASI) GO TO 600 *** REF.(2) EQ(30B) SQNM2=SQNM1 SQNM1=ROOT(NM1) SQNPM2=SQNPM1 SQNPM1=ROOT(NPM1) SQ1=SQNM1*SQNPM1 A1=S*N21/SQ1 B2=-S2*SQ1/(SQNM2*SQNPM2) GO TO 700 *** REF.(2), EQ.(30) 600 A1=(S*N21)/NM1 B2=-(S2*NPM1)/NM2 700 A1T=A1*T A1U=A1*U N21=N21-2 CK=C(K) CK1=C(K+1) K=K-N21 *** REF.(2), EQ(33) V2=VC1 VC1=VC VC=VC1*A1T+V2*B2+CK V2=VS1 VS1=VS VS=VS1*A1T+V2*B2+CK1 IF(.NOT.DERIV1) GO TO 1000 CKZ=CK*N1 CK1Z=CK1*N1 *** REF.(2), EQ(10) V2=VXC1 VXC1=VXC VXC=VXC1*A1T+VC1*A1U+V2*B2 V2=VXS1 VXS1=VXS VXS=VXS1*A1T+VS1*A1U+V2*B2 V2=VZC1 VZC1=VZC VZC=VZC1*A1T+V2*B2-CKZ V2=VZS1 VZS1=VZS VZS=VZS1*A1T+V2*B2-CK1Z N1=N IF(.NOT.DERIV2) GO TO 1000 N2=N+2 *** REF.(2), EQ(41) V2=VZZC1 VZZC1=VZZC VZZC=VZZC1*A1T+V2*B2+N2*CKZ V2=VZZS1 VZZS1=VZZS VZZS=VZZS1*A1T+V2*B2+N2*CK1Z IF(NPOLE) GO TO 800 *** REF.(2), EQ(12) SECOND-ORDER DERIVATIVE WRT LATITUDE V2=VXXC1 VXXC1=VXXC VXXC=A1T*(VXXC1-VC1)+(A1U+A1U)*VXC1+V2*B2 V2=VXXS1 VXXS1=VXXS VXXS=A1T*(VXXS1-VS1)+(A1U+A1U)*VXS1+V2*B2 *** REF.(2) EQ(10,40) DERIVATIVE WRT R AND LATITUDE 800 V2=VXZC1 VXZC1=VXZC VXZC=VXZC1*A1T+VZC1*A1U+V2*B2 V2=VXZS1 VXZS1=VXZS VXZS=VXZS1*A1T+VZS1*A1U+V2*B2 1000 CONTINUE SU(M+1)=VC SU(M+I1)=VS IF(.NOT.DERIV1) GO TO 1500 SU(M+I2)=VXC SU(M+I3)=VXS SU(M+I4)=VZC SU(M+I5)=VZS IF(.NOT.DERIV2) GO TO 1500 SU(M+I6)=VZZC SU(M+I7)=VZZS SU(M+I8)=VXZC SU(M+I9)=VXZS GO TO 1500 1300 VC=SU(M+1) VS=SU(M+I1) IF(.NOT.QUASI) GO TO 1400 SQNPM1=ROOT(MAX2) SQNPM2=ROOT(MAX2+1) 1400 NPM1=MAX2 MAX2=MAX2-2 IF(.NOT.DERIV1) GO TO 1500 c - 1st order derivative stuff VXC=SU(M+I2) VXS=SU(M+I3) VZC=SU(M+I4) VZS=SU(M+I5) IF(.NOT.DERIV2)GO TO 1500 c - 2nd order derivative stuff VZZC=SU(M+I6) VZZS=SU(M+I7) VXZC=SU(M+I8) VXZS=SU(M+I9) *** THE COMPUTATION OF DERIVATIVES IN DIRECTION OF POS LONGITUDE,Y, *** INVOLVES THE DIVISION BY U=COS(LATITUDE). THIS DIVISION IS *** PERFORMED IMPLICITLY, BY STOPPING THE CLENSHAW SUMMATION AT M=1. *** THIS IS DONE BY PUTTING U0=1.0. THIS TRICK PERMITS THE USE OF THE *** PROCEDURE AT POLES, EXCEPT FOR THE SECOND-ORDER DERIVATIVE WRT *** LONGITUDE. BUT HERE WE USE THE VALIDITY OF THE LAPLACE EQUATION *** AND COMPUTE THE SECOND-ORDER DERIVATIVES WRT X AND Z 1500 U0=U IF(M.EQ.0)U0=1.D0 *** REF.(2) EQ.(35) AUX=NPM1 IF(QUASI)AUX=SQNPM1/SQNPM2 M21=S*AUX M21U=M21*U IF(.NOT.DERIV1) GO TO 1700 M21T=M21*T M21U0=M21*U0 IF(.NOT.DERIV2) GO TO 1600 c - 2nd order derivatives: VZZM=VZZC*CM+VZZS*SM+M21U*VZZM IF(M.GT.0)VXYM=M*(VXS*CM-VXC*SM)+M21U*VXYM-M21T*VYM VXZM=VXZC*CM+VXZS*SM-M21T*VZM+M21U*VXZM VYZM=(VZS*CM-VZC*SM)*M+M21U0*VYZM IF(POLE) VXXM=VXXC*CM+VXXS*SM+M21*(U*(VXXM-VM)-T2*VXM) IF(NPOLE)VYYM=-(VC*CM+VS*SM)*M2+M21U0*VYYM c - 1st order derivatives: 1600 VXM=VXC*CM+VXS*SM-M21T*VM+M21U*VXM VYM=M*(VS*CM-VC*SM)+M21U0*VYM VZM=VZC*CM+VZS*SM+M21U*VZM c - 0th order 1700 VM=VC*CM+VS*SM+M21U*VM c write(6,*) ' VM = ',VM cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c - ASK YOURSELF THIS -> DOES VM represent a sum from 0 to NMAX????? c - Answer, 3/14/97: The sum over loop 1700 is from c - IZ to CAPN; IZ is initialized to ZERO in subroutine inital c - and never changes. CAPN is equal to NMAX (actually it is c - equal to -NEGN where NEGN = -NMAX, so CAPN = --NMAX = NMAX. c - So the answer is : YES! cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc *** COMPUTE GM/R S=CM3/R GPTDR=S*VM IF(.NOT.DERIV1) RETURN *** COMPUTE GM/R**2 S=S/R G1(1)=S*VXM G1(2)=S*VYM c - das (see above, remember CM3 currently is GM) c - changed for v0.2a of geopot97: G1(3)=VZM*S c G1(3)=VZM*S+U**2*OM2*R c G1(3)=VZM*S+U**2*OM2*R - (CM3/R/R)*C0 c - das c - das c - changed back for geopot97 (but specifically for v0.2a): IF(.NOT.DERIV2) RETURN c IF(.NOT.DERIV2) then c CM1 = CM1dum c CM2 = CM2dum c CM3 = CM3dum c RETURN c endif c - das *** COMPUTE GM/R**3 S=S/R *** HERE THE LAPLACE EQUATION IS USED IF(NPOLE) GO TO 1900 VXXM=VXXM+VZM VYYM=-(VXXM+VZZM) GO TO 2000 1900 VYYM=VZM+(VYYM-T*VXM)/U VXXM=-(VZZM+VYYM) c - changed for v0.2a of geopot97: c2000 G2(1,1)=VXXM*S+OM2*T**2 2000 G2(1,1)=VXXM*S *** man. geod. patch G2(1,2)=S*VXYM*m21 G2(2,1)=G2(1,2) c - changed for v0.2a of geopot97: G2(1,3)=S*(VXZM-VXM) c G2(1,3)=S*(VXZM-VXM)-U*T*OM2 G2(3,1)=G2(1,3) c - changed for v0.2a of geopot97: G2(2,2)=VYYM*S c G2(2,2)=VYYM*S+OM2 G2(2,3)=S*(VYZM-VYM) G2(3,2)=G2(2,3) c - das (see above, remember that CM3 currently is GM) c - changed for v0.2a of geopot97: G2(3,3)=S*VZZM c G2(3,3)=S*VZZM+U**2*OM2 c G2(3,3)=S*VZZM+U**2*OM2 + (2*CM3/R/R/R)*C0 c - das c - das c - switched back for geopot97 (but specifically for v0.2a): c CM1 = CM1dum c CM2 = CM2dum c CM3 = CM3dum c - das RETURN END c c c subroutine norm(a,j2,gm,omega,b,j,nmax,u0,gammae,gammap) c - This subroutine takes, as input, the normal gravity field values c - of a,j2,gm,omega and b (b is gotten from "getb" first), and c - computes the values of the normal gravity field: c - j(nmax) = coefficients of the spherical harmonic expansion of c - the normal gravitational potential c - nmax is equal to the highest degree desired divided by 2 c - (thus if j2 through j20 are desired, set n = 10, c - and receive j2 = j(1), j4 = j(2), ..., j20 = j(10) c - (see H/M, eqn after 2-91) c - (note that j2 is input, not calculated) c - u0 = value of normal gravity potential on the ellipsoid c - gammae = normal gravity, on the ellipsoid, at the equator c - gammap = normal gravity, on the ellipsoid, at the pole implicit real*8(a-z) integer n,nmax real*8 j(nmax) c write(6,*) ' inside norm' bige = dsqrt(a**2 - b**2) e = bige/a ! Added for version 0.4a ep = bige/b c - H/M (exact formula), 2-61, p. 67 u0 = (gm/bige)*datan2(bige,b) + omega**2*a**2/3 c write(6,*) ' u0 = ',u0 q0 = 0.5d0*((1.d0+3*b**2/bige**2)*datan2(bige,b) * - 3*b/bige) c - H/M, exact formula, 2-67 (applied at u=b), p. 68 qp0 = 3 * (1.d0 + b**2/bige**2) * *(1.d0 - datan2(bige,b)*b/bige) - 1.d0 c x1 = omega**2 * a**2 / 3 / q0 c x2 = gm/bige c j(1) = j2 c do 1 n=2,nmax c xn = (-1)**n * bige**(2*n+1) / (2*n+1) c j(n) = (x1*(2*n)/(2*n+3) - x2)*xn / a**(2*n) / gm c write(6,*) 2*n,j(n) c 1 continue c - Change for version 0.4c: c - Added the definition of j(1)!!!!!!! j(1) = j2 c - Change for version 0.4a: c - Replaced 10/2/1997. What was I thinking above??? do 101 n=2,nmax j(n) = ((-1)**(n+1))*((3*e**(2*n))/(2*n+1)/(2*n+3)) * *(1-n+5*n*j(1)/e/e) write(6,*) 2*n,j(n) 101 continue c - H/M, exact formula, 2-70, p. 69 m = omega**2 * a**2 * b / gm c - H/M, exact formula, 2-73, 2-74, p. 69 gammae = (gm/a/b)*(1.d0 - m - m*ep*qp0/6/q0) gammap = (gm/a/a)*(1.d0 + m*ep*qp0/3/q0) c write(6,*) ' gammae = ',gammae c write(6,*) ' gammap = ',gammap return end c c c SUBROUTINE INITAL *** INITIALIZE COMMON BLOCKS IMPLICIT REAL*8 (A-H,O-Z),INTEGER*4(I-N) LOGICAL FIRST INTEGER*4 OLDORD COMMON/GPTCM/OLDT,OLDR,IZ,FIRST,OLDORD,I1,I2,I3,I4, * I5,I6,I7,I8,I9,NMAXSV common/block1/pi,d2r c write(6,*) ' inside inital' c - Initialize block "GPTCM" IZ=0 FIRST=.FALSE. OLDT=0.D0 OLDR=0.D0 OLDORD=0 I1=0 I2=0 I3=0 I4=0 I5=0 I6=0 I7=0 I8=0 I9=0 c - Initialize block "block1" pi = 2.d0 * dasin(1.d0) d2r = pi/180.d0 RETURN END c c c subroutine ell2sph(glat,glon,h,a,finv,slat,slon,r) c - subroutine to convert ellipsoidal to spherical coordinates c - glat,glon,h = ellipsoidal (input) coordinates c - slat,slon,r = spherical (output) coordinates c - a,finv = semi-major axis and inverse flattening of ellipsoid implicit real*8(a-h,o-z) c write(6,*) ' inside ell2sph' c - Convert to X,Y,Z first call ell2xyz(glat,glon,h,a,finv,x,y,z) c - Convert X,Y,Z to spherical call xyz2sph(x,y,z,slat,slon,r) c - Added because somehow the longitude got a sign mixed up slon = glon return end c c c subroutine ell2xyz(glat,glon,h,a,finv,x,y,z) c - subroutine to convert ellipsoidal to cartesian coordinates c - glat,glon,h = ellipsoidal (input) coordinates c - x,y,z = cartesian (output) coordinates c - a,finv = semi-major axis and inverse flattening of ellipsoid implicit real*8(a-h,o-z) common/block1/pi,d2r f = 1.d0/finv e2 = 2.d0*f - f*f sphi = dsin(glat*d2r) xw = dsqrt(1.d0-e2*sphi*sphi) xn = a / xw cphi = dcos(glat*d2r) slam = dsin(glon*d2r) clam = dcos(glon*d2r) x = (xn + h)*cphi*clam y = (xn + h)*cphi*slam z = (xn*(1.d0-e2) + h)*sphi return end c c c subroutine xyz2sph(x,y,z,slat,slon,r) c - subroutine to convert cartesian coordinates to spherical implicit real*8(a-h,o-z) common/block1/pi,d2r c - x,y,z = cartesian coordinates c - slat,slon,r = spherical coordinates p = dsqrt(x*x + y*y) r = dsqrt(p*p + z*z) slon = datan2(y,x)/d2r slat = datan2(z,p)/d2r return end c c c subroutine bestfit(gmcoef,omega,abf,j2bf,finvbf,w0,tidbf,tidcoef) c - subroutine to "define" the geoid by 1 of 4 methods: c 1) Define the best fitting ellipsoid to the geoid, in the c same tide system as the coefficients, by giving abf and finvbf c 2) Define the best fitting ellipsoid to the geoid, in the c same tide system as the coefficients, by giving abf and j2bf c (j2 can be given, since we have gmcoef and omega) c 3) Define the best fitting ellipsoid to the geoid, in the c same tide system as the coefficients, by giving abf and bbf c 4) Define the gravity potential value on the geoid, Wo implicit real*8(a-h,o-z) real*8 j2bf character*4 tidbf,tidcoef write(6,1) 1 format(/1x,'To properly use the coefficients in combination',/, * 1x,'with a reference field, you *must* define the',/, * 1x,'geoid in one of 4 ways:',/, * 1x,' 1) Define a best-fitting ellipsoid to the geoid',/, * 1x,' using a and f, and knowing the tide system.',/, * 1x,' 2) Define a best-fitting ellipsoid to the geoid',/, * 1x,' using a and J2, and knowing the tide system.',/, * 1x,' 3) Define a best-fitting ellipsoid to the geoid',/, * 1x,' using a and b, and knowing the tide system.',/, * 1x,' 4) Define the gravity potential on the geoid, Wo') write(6,2) 2 format(/1x,'Which way do you wish to define the geoid? ',$) read(5,*)ians if(ians.le.0 .or. ians.ge.5)stop c - Get a/f, a/b or a/j2 if(ians.eq.1 .or. ians.eq.2 .or. ians.eq.3)then write(6,3) 3 format(1x,'Input the best equatorial radius (a): ',$) read(5,*)abf if(ians.eq.1)then write(6,4) 4 format(1x,'Input the best inverse flattening (1/f): ',$) read(5,*)finvbf bbf = abf*(1.d0 - (1.d0/finvbf)) call getj2(abf,gmcoef,omega,bbf,j2bf) elseif(ians.eq.2)then write(6,7) 7 format(1x,'Input the best J2: ',$) read(5,*)j2bf call getb(abf,j2bf,gmcoef,omega,bbf) finvbf = abf / (abf - bbf) elseif(ians.eq.3)then write(6,9) 9 format(1x,'Input the best semi-major axis (b): ',$) read(5,*)bbf finvbf = abf / (abf - bbf) call getj2(abf,gmcoef,omega,bbf,j2bf) endif c - Get tide system of best-fit ellipsoid write(6,5) 5 format(1x,'Which tide system is this best ellipsoid in?',/, * 1x,' 1) = Mean',/, * 1x,' 2) = Zero',/, * 1x,' 3) = Non-') read(5,*)ians2 if(ians2.le.0 .or. ians2.ge.4)stop if(ians2.eq.1)tidbf = 'mean' if(ians2.eq.2)tidbf = 'zero' if(ians2.eq.3)tidbf = 'non-' c - Deal with possibly differing tide systems c - Steps: 1) Transform abf and bbf into the tide system of c the coefficients (abft, bbft) c 2) Calculate Uo from abft,bbft,gmcoef,omega c 3) Set Wo = Uo c - Step 1 is only performed if tidbf is not tidcoef c - Side note to myself...the transformation should NOT c - be necessary, as Uo is identical for all three tide c - systems, given the same GM/omega and the properly c - transformed a/b values! But, for now, leave it c - in. if(tidbf.ne.tidcoef)then write(6,6) 6 format(1x, *'Your best-fit ellipsoid is in a different tide system than'/ *'the coefficients. Transforming your ellipsoid to the tide'/ *'system of the coefficients...') call fixab(abf,bbf,tidbf,tidcoef,abft,bbft) else abft = abf bbft = bbf endif call getu0(abft,gmcoef,omega,bbft,u0) W0 = u0 return c - Get Wo directly, rather than a/f, a/b or a/j2 elseif(ians.eq.4)then write(6,8) 8 format(1x,'Input the Wo value for the geoid: ',$) read(5,*) W0 c - The "tide system of Wo" is irrelevant. Consider this: c - Given GM/omega/a/b in the mean tide system, you get a c - Uo value. Then if you transform to the zero or the c - non-tidal system (changing a and b), you get almost c - EXACTLY the same Uo value!!! (To within .001 m**s/s**2) c - As such, the tide system is NOT tied to the potential c - value, but just to the SHAPE of the ellipsoid/geoid. endif return end c c c subroutine gettid(tidsum,tidcoef) c - Subroutine to define the tide system for the summation c - and thus the output (irregardless of what tide system c - the coefficients/normal field/best fitting ellipsoid are in). character*4 tidsum,tidcoef write(6,1)tidcoef 1 format(/1x,'The tide system for any output quantities can '/ * 1x,'differ from that of the coefficients and/or the '/ * 1x,'normal field (if used). Therefore, what tide'/ * 1x,'system shall all values be computed in?'/ * 1x,'1: mean',/,1x'2: zero',/,1x,'3: non-',/, * 1x,'(The coefficients are in the ',a4,' tide system)') read(5,*)ians if(ians.le.0 .or. ians.ge.4)stop if(ians.eq.1)tidsum = 'mean' if(ians.eq.2)tidsum = 'zero' if(ians.eq.3)tidsum = 'non-' return end c c c subroutine fixc20(tidcoef,tidsum,c20,gm,a) c - subroutine to modify the C(2,0) term because the tide c - system that C(2,0) is in differs from the tide system c - requested for output character*4 tidcoef,tidsum real*8 kto,kfrom write(6,3)tidcoef write(6,4)tidsum 3 format(/,1x,50('*'),/, * 1x,'The tide system of your coefficients is: ',a4, *' tide') 4 format(1x,'You want the potential computed in the : ',a4, *' tide system,',/,1x,'requiring a transformation of C(2,0)') if(tidsum.eq.'non-')then write(6,1) 1 format(1x,'To transform *to* the non-tidal system, you'/ * 1x,'must input a value of the k Love number'/ * 1x,'(often taken to be k = 0.3) : ',$) read(5,*)kto endif if(tidcoef.eq.'non-')then write(6,2) 2 format(1x,'To transform *from* the non-tidal system, you'/ * 1x,'must input the value of k which was used to'/ * 1x,'put the coefficients in the non-tidal system'/ * 1x,'initially (often k = 0.3, but DO NOT GUESS!): ',$) read(5,*)kfrom endif if(tidcoef.ne.'non-' .and. tidsum.ne.'non-')then write(6,5) 5 format(1x,'Since you are neither going to nor from the',/, * 1x,'non-tidal system, the transformation is simple',/, * 1x,'and does not require the k Love number') endif if(tidsum.eq.'non-')then if(tidcoef.eq.'mean')then conv = 1.d0 + kto elseif(tidcoef.eq.'zero')then conv = kto endif elseif(tidsum.eq.'zero')then if(tidcoef.eq.'non-')then conv = -kfrom elseif(tidcoef.eq.'mean')then conv = 1.d0 endif elseif(tidsum.eq.'mean')then if(tidcoef.eq.'non-')then conv = -(1.d0 + kfrom) elseif(tidcoef.eq.'zero')then conv = -1.d0 endif endif c s5 = 1.d0/dsqrt(5.d0) c tot = -0.198d0 * s5 * gamma * r*r*r / gm / a / a write(6,6) 6 format(/1x,'For now, the conversion factor is HARD-CODED',/, *'as -1.39844d-8 x 1 or 1+k or k. If this is bad, ',/, *'you must get in the code and re-compile.') tot = -1.39844d-8 c20 = c20 - conv * tot return end c c c subroutine getu0(a,gm,omega,b,u0) c - subroutine to compute Uo from a,GM,omega and b implicit real*8(a-h,o-z) bige = dsqrt(a**2 - b**2) c - H/M (exact formula), 2-61, p. 67 u0 = (gm/bige)*datan2(bige,b) + (omega**2*a**2)/3 return end c c c subroutine getb(a,j2,gm,omega,b) c - This subroutine takes, as input, the normal gravity field c - variables of a, j2, gm and omega, and iterates to an c - "exact" solution of "b" (where the iteration stops c - when the difference between the b-implied "j2" value c - and the true "j2" value is less than "tol", in absolute c - value. "tol" is hard-coded as 1.d-15 c - implicit real*8(a-z) step = 100.d0 tol = 1.d-15 fp = 1.5d0*j2 + 0.5d0 * omega**2 * a / 9.8d0 bp = a * (1.d0 - fp) c1 = -(1.d0/3/a**2/gm) c2 = 2.d0*(omega*a)**2 / 15.d0 c3 = gm j2old = j2 1 bigep = dsqrt(a**2 - bp**2) x1 = 1.d0 + 3.d0 * bp**2 / bigep**2 x2 = x1 * datan2(bigep,bp) x3 = x2 - 3.d0*bp/bigep q0p = 0.5d0 * x3 j2p = c1 * (c2*bigep**3/q0p - c3*bigep**2) c write(6,*) bp,j2p dj2 = j2old - j2p c write(6,*) ' j2old = ',j2old c write(6,*) ' j2p = ',j2p c write(6,*) ' dj2 = ',dj2 c write(6,*) ' ---------' if(dabs(dj2).lt.tol)goto 95 dj2 = j2 - j2p if(dj2.lt.0)then idum = +1 if(idum.ne.idum2)step = step / 2 bp = bp + step else idum = -1 if(idum.ne.idum2)step = step / 2 bp = bp - step endif j2old = j2p idum2 = idum goto 1 95 b = bp c write(6,*) ' b = ',bp c write(6,*) ' j2calc = ',j2p c write(6,*) ' j2true = ',j2 return end c c c subroutine getj2(a,gm,omega,b,j2) implicit real*8(a-z) c - This subroutine calcultes, non-iteratively, the value c - of J2 implied by 4 parameters of the normal field c - (a,gm,omega and b). bigE = dsqrt(a**2 - b**2) c write(6,*) ' bigE =',bigE ep = bigE/b epinv = 1.d0/ep e2 = bigE**2/a**2 c write(6,*) ' e2 = ',e2 q0 = 0.5d0*( (1.d0+3*epinv**2)*datan2(bigE,b) - 3*epinv) m = omega**2 * a**2 * b / gm c write(6,*) ' m = ',m j2 = (e2/3.d0)*(1.d0 - (2.d0/15.d0)*m*ep/q0) c write(6,*) ' j2 = ',j2 return end c c c subroutine fixab(abf,bbf,tidbf,tidcoef,abft,bbft) c - Subroutine to convert abf/bbf (a and b of the c - ellipsoid in the "tidbf" tide system, to abft/bbft c - (a and b of the ellipsoid in the "tidcoef" tide system). c - That is: transform the best fit ellipsoid to the tide c - system of the coefficients implicit real*8(a-h,o-z) character*4 tidbf,tidcoef real*8 kto,kfrom write(6,3)tidcoef write(6,4)tidbf 3 format(1x,'The tide system of your coefficients is: ',a4, *' tide') 4 format(1x,'You gave the ellipsoid in the: ',a4, *' tide system,',/,1x,'requiring a transformation of a & b') if(tidbf.eq.'non-')then write(6,1) 1 format(1x,'To transform *to* the non-tidal system, you'/ * 1x,'must input a value of the k Love number'/ * 1x,'(often taken to be k = 0.3) : ',$) read(5,*)kto endif if(tidcoef.eq.'non-')then write(6,2) 2 format(1x,'To transform *from* the non-tidal system, you'/ * 1x,'must input the value of k which was used to'/ * 1x,'put the coefficients in the non-tidal system'/ * 1x,'initially (often k = 0.3, but DO NOT GUESS!): ',$) read(5,*)kfrom endif if(tidcoef.ne.'non-' .and. tidbf.ne.'non-')then write(6,5) 5 format(1x,'Since you are neither going to nor from the',/, * 1x,'non-tidal system, the transformation is simple',/, * 1x,'and does not require the k Love number') endif if(tidbf.eq.'non-')then if(tidcoef.eq.'mean')then conv = -(1.d0 + kto) elseif(tidcoef.eq.'zero')then conv = -kto endif elseif(tidbf.eq.'zero')then if(tidcoef.eq.'non-')then conv = kfrom elseif(tidcoef.eq.'mean')then conv = -1.d0 endif elseif(tidbf.eq.'mean')then if(tidcoef.eq.'non-')then conv = 1.d0 + kfrom elseif(tidcoef.eq.'zero')then conv = 1.d0 endif endif abft = abf - (0.099d0)*conv bbft = bbf + (0.198d0)*conv return end c c c subroutine getref(a,b,finv,j2,gm,omega,tidsum) c - Subroutine to input the normal gravity field implicit real*8(a-h,o-z) real*8 j2 character*4 tidsum,tidnorm character*1 yesno 101 write(6,1) 1 format(/1x,'You are now about to input the parameters of the ', *'ellipsoidal normal gravity',/, *1x,'field. This is NOT the so-called best-fitting ellipsoid, ', *'but is rather the ',/, *1x,'actual ellipsoid to which the undulations will refer. ', *'The method to be used is',/, *1x,'method #3 (Notebook DRU-3, p. 58), which has NO geometric ', *'transformations.',/, *1x,'So, let us proceed.....',//) 3 write(6,2) 2 format(1x,'How shall the normal field be described?',/, *1x,'1 - a,J2,GM,omega',/,1x,'2 - a,f,GM,omega',/, *1x,'3 - a,b,GM,omega') read(5,*)ians if(ians.le.0 .or. ians.ge.4)goto 3 write(6,4) 4 format(1x,'Input the equatorial radius (a) = ',$) read(5,*) a write(6,5) 5 format(1x,'Input the GM value = ',$) read(5,*) gm write(6,6)omega 6 format(1x,'To remind you, Omega is defined as:',e13.7) c - Input either j2,finv or b, but compute the other two as well if(ians.eq.1)then write(6,7) 7 format(1x,'Input the J2 value = ',$) read(5,*)j2 call getb(a,j2,gm,omega,b) finv = a/(a-b) elseif(ians.eq.2)then write(6,8) 8 format(1x,'Input the INVERSE flattening (1/f) = ',$) read(5,*)finv b = a - (a/finv) call getj2(a,gm,omega,b,j2) elseif(ians.eq.3)then write(6,9) 9 format(1x,'Input the semi-minor axis (b) = ',$) read(5,*)b finv = a/(a-b) call getj2(a,gm,omega,b,j2) endif c - get the tide system of the normal field (which we will NOT c - use. We will transform the normal field into the tidsum c - system (which is also what the total coefficients will be c - transformed into, in another part of the program, to ensure c - that disturbing potential coefficients are computed from c - total/normal coefficients in a consistant system) 11 write(6,10) 10 format(/1x,'What tide system is this normal field in?',/, * 1x,'1: mean',/,1x,'2: zero',/,1x,'3: non-') read(5,*)ians if(ians.le.0 .or. ians.ge.4)goto 11 if(ians.eq.1)tidnorm = 'mean' if(ians.eq.2)tidnorm = 'zero' if(ians.eq.3)tidnorm = 'non-' c - get the final ref field in the "tidsum" system if(tidnorm.eq.tidsum)return write(6,12) 12 format(/1x,'You have input a normal field in a tide',/, *1x,'system different than you have requested for output.',/, *1x,'In order to have consistancy in computing the ',/, *1x,'disturbing coefficients, we will now transform your',/, *1x,'normal field to the tide system of the output.'//) write(6,13) 13 format(1x,20('!'),'IMPORTANT',20('!')) write(6,14) 14 format(1x,'Note that this *WILL* change the SHAPE of the ',/, *1x,'ellipsoid. I hope you thought about that. If not,',/, *1x,'you can go back and try entering the normal field',/, *1x,'again.') write(6,15) 15 format(1x,'y = YES, please let me re-enter the normal field.', */,1x,'n = NO, I know what I am doing. Proceed, darnit!') read(5,'(a)')yesno if(yesno.eq.'y' .or. yesno.eq.'Y')goto 101 write(6,16) 16 format(1x,'Ok, you asked for it....') call fixab2(a,b,tidnorm,tidsum,afixed,bfixed) a = afixed b = bfixed finv = a/(a-b) call getj2(a,gm,omega,b,j2) return end c c c subroutine dist(c,max,nmax,gmcoef,acoef,jnorm,njmax,gmnorm,anorm) c - subroutine to compute disturbing potential coefficients, based c - on total coefs, normal coefs, and the GM/a values associated c - with them both. It is assumed that total/normal coefs are in c - an identical tide system, whatever that system may be. c - Also, the TOTAL (c) coefficients, are coming in as c - quasi-normlized. The NORMAL (jnorm) coefficients are coming c - in as unnormalized, and of opposite sign as c. We need c - to make the jnorm coefficients of the right sign, and c - quasinormalized, before subtracting. implicit real*8(a-h,o-z) real*8 c(0:(max+1)*(max+1)-1) real*8 jnorm(njmax),j0 c - Do C(0,0) - scale*J(0) first: c - J(0) is the same, whether conventional, fully normalized c - or quasi-normalized (= 1.d0) j0 = 1.d0 scale = (gmnorm/gmcoef)*1.d0 ! the 1.d0 = (anorm/acoef)**0 c(0) = c(0) - j0*scale c - Now do the even zonals from 2 to 2*njmax, transforming the c - jnorm values into fully-normalized values to be consistant c - with the c(*) values. do 1 i=1,njmax n = i*2 scale = (gmnorm/gmcoef) * (anorm/acoef)**n c - Converts from J2 (conventional) to J2(fully normalized) xnorm1 = -1.d0/dsqrt(2.d0*n + 1.d0) c - Converts from J2 (fully normalized) to J2 (quasi-normalized) xnorm2 = dsqrt(2.d0*n + 1.d0) c - Thus, J2(quasi-normalized) = J2(conventional)*(-1.d0) c (xnorm1*xnorm2 = -1.d0) xnorm = xnorm1 * xnorm2 ! = -1.d0 isto = n*n ! <- Just happens to be this way c(isto) = c(isto) - scale*(xnorm)*jnorm(i) 1 continue return end c c c subroutine undist(c,max,nmax,gmcoef,acoef,jnorm,njmax,gmnorm * ,anorm) c - subroutine to UNDISTURB the coefficients (basically c - un-does the dist(*) subroutine implicit real*8(a-h,o-z) real*8 c(0:(max+1)*(max+1)-1) real*8 jnorm(njmax),j0 c - Do C(0,0) first: c - J(0) is the same, whether conventional, fully normalized c - or quasi-normalized (= 1.d0) j0 = 1.d0 scale = (gmnorm/gmcoef)*1.d0 ! the 1.d0 = (anorm/acoef)**0 c(0) = c(0) + j0*scale c - Now do the even zonals from 2 to 2*njmax, transforming the c - jnorm values into fully-normalized values to be consistant c - with the c(*) values. do 1 i=1,njmax n = i*2 scale = (gmnorm/gmcoef) * (anorm/acoef)**n c - Converts J2 (conventional) to J2(fully normalized) xnorm1 = -1.d0/dsqrt(2.d0*n + 1.d0) c - Converts J2 (fully normalized) to J2(quasi-normalized) xnorm2 = dsqrt(2.d0*n + 1.d0) c - Thus, J2(quasi-normalized) = J2(conventional)*(-1.d0) c (xnorm1*xnorm2 = -1.d0) xnorm = xnorm1 * xnorm2 ! = -1.d0 isto = n*n ! < - Just happens to be this way c(isto) = c(isto) + scale*(xnorm)*jnorm(i) 1 continue return end c c c subroutine getgrid(glamn,glomn,dla,dlo,nla,nlo,ikind) implicit real*8(a-h,o-z) write(6,1) 1 format(1x, *'Input south geodetic latitude boundary (degrees): ',$) read(5,*)glamn write(6,2) 2 format(1x, *'Input north geodetic latitude boundary (degrees): ',$) read(5,*)glamx write(6,3) 3 format(1x,'Input latitude grid spacing (minutes): ',$) read(5,*)dlamin write(6,4) 4 format(1x, *'Input west geodetic longitude boundary (degrees): ',$) read(5,*)glomn write(6,5) 5 format(1x, *'Input east geodetic longitude boundary (degrees): ',$) read(5,*)glomx write(6,6) 6 format(1x,'Input longitude grid spacing (minutes): ',$) read(5,*)dlomin dla = dlamin/60.d0 dlo = dlomin/60.d0 nla = nint((glamx-glamn)/dla + 1.d0) nlo = nint((glomx-glomn)/dlo + 1.d0) ikind = 1 return end c c c subroutine getgrid2(glamn,glomn,dla,dlo,nla,nlo,ikind) implicit real*8(a-h,o-z) write(6,1) 1 format(1x, *'Input south spherical latitude boundary (degrees): ',$) read(5,*)glamn write(6,2) 2 format(1x, *'Input north spherical latitude boundary (degrees): ',$) read(5,*)glamx write(6,3) 3 format(1x,'Input latitude grid spacing (minutes): ',$) read(5,*)dlamin write(6,4) 4 format(1x, *'Input west spherical longitude boundary (degrees): ',$) read(5,*)glomn write(6,5) 5 format(1x, *'Input east spherical longitude boundary (degrees): ',$) read(5,*)glomx write(6,6) 6 format(1x,'Input longitude grid spacing (minutes): ',$) read(5,*)dlomin dla = dlamin/60.d0 dlo = dlomin/60.d0 nla = nint((glamx-glamn)/dla + 1.d0) nlo = nint((glomx-glomn)/dlo + 1.d0) ikind = 1 return end c c c double precision function gam(glat,a,b,ge,gp) implicit real*8(a-h,o-z) pi = 2.d0 * dasin(1.d0) d2r = pi/180.d0 c = dcos(glat*d2r) s = dsin(glat*d2r) c2 = c*c s2 = s*s gam = (a*ge*c2 + b*gp*s2 )/ dsqrt(a*a*c2 + b*b*s2) return end c c c double precision function dgdh(glat,a,finv,gamma,omega) implicit real*8(a-h,o-z) pi = 2.d0 * dasin(1.d0) d2r = pi/180.d0 f = 1.d0/finv e2 = 2.d0*f - f*f s2 = dsin(glat*d2r) xw = dsqrt(1.d0 - e2*s2) xm = a*(1.d0 - e2)/xw/xw/xw xn = a/xw xj = ( (1.d0/xm) + (1.d0/xn) ) / 2.d0 c - H/M, p 70, eqn 2-79 dgdh = -2.d0*gamma*xj - 2.d0*omega*omega return end c c c double precision function d2gdh2(gamma,a) implicit real*8(a-h,o-z) d2gdh2 = 6.d0*gamma/a/a return end c c c subroutine mknorm(cj,max,nmax,gmcoef,acoef,jnorm,njmax,gmnorm * ,anorm) implicit real*8(a-h,o-z) real*8 cj(0:(max+1)*(max+1)-1) real*8 jnorm(njmax),j0 c - Do CJ(0,0) = scale*J(0) first: c - J(0) is the same, whether conventional, fully normalized c - or quasi-normalized (= 1.d0) j0 = 1.d0 scale = (gmnorm/gmcoef)*1.d0 ! the 1.d0 = (anorm/acoef)**0 cj(0) = j0*scale c - Now do the even zonals from 2 to 2*njmax, transforming the c - jnorm values into fully-normalized values to be consistant c - with the c(*) values. do 1 i=1,njmax n = i*2 scale = (gmnorm/gmcoef) * (anorm/acoef)**n c - Converts from J2 (conventional) to J2(fully normalized) xnorm1 = -1.d0/dsqrt(2.d0*n + 1.d0) c - Converts from J2 (fully normalized) to J2 (quasi-normalized) xnorm2 = dsqrt(2.d0*n + 1.d0) c - Thus, J2(quasi-normalized) = J2(conventional)*(-1.d0) c (xnorm1*xnorm2 = -1.d0) xnorm = xnorm1 * xnorm2 ! = -1.d0 isto = n*n ! <- Just happens to be this way cj(isto) = scale*(xnorm)*jnorm(i) 1 continue return end c c c subroutine legpol(slat,nmax,p) c - Subroutine to recursively (and stabily) compute c - Legendre polynomials. NOT Fully normalized at the end. c - P(1) = P_sub_zero c - P(2) = P_sub_one c - etc... implicit real*8(a-h,o-z) real*8 p(nmax) pi = 2.d0*dasin(1.d0) d2r = pi/180.d0 c theta = (90.d0 - slat)*d2r c x = dcos(theta) x = dsin(slat*d2r) do 3 i=1,nmax p(i) = 0.d0 3 continue p(1) = 1.d0 p(2) = x do 1 n=2,nmax-1 c write(6,*) ' n,p(n),p(n-1) = ',n,p(n),p(n-1) p(n+1) = 2.d0*x*p(n) - p(n-1) - ( (x*p(n)-p(n-1)) / dble(n)) 1 continue c do 2 n=1,nmax c nn = n-1 c p(n) = dsqrt(2.d0*nn + 1.d0) * p(n) c write(6,*) slat,n-1,p(n) c 2 continue return end c c c double precision function pnm(n,m,xlat) implicit real*8(a-h,o-z) common/block1/pi,d2r c t = dcos(theta*d2r) t = dsin(xlat*d2r) xm2 = dble(m)/2.d0 if(mod((n-m),2).eq.0)then ir = (n-m)/2 else ir = (n-m-1)/2 endif cons = ((1.d0-t*t)**xm2) / (2.d0**n) xsum = 0.d0 do 1 k=0,ir xnum = fact(2*n - 2*k) xden = fact(k)*fact(n-k)*fact(n-m-2*k) sca = (-1.d0)**k * t**(n-m-2*k) xsum = xsum + sca*xnum/xden 1 continue pnm = cons*xsum write(6,*) ' n,m,xlat,pnm = ',n,m,xlat,pnm return end c c c double precision function fact(n) fact = 1.d0 if(n.eq.0)return do 1 i=1,n fact = fact * i 1 continue return end c c c subroutine sumnorm(xlat,xlon,r,gmnorm,anorm,jnorm,njmax,pleg, * V2,V2N,V2E,V2U,U,UN,UE,UU,PHI,PHIN,PHIE,PHIU,omega) implicit real*8(a-h,o-z) c - Subroutine to sum up NORMAL gravitaitonal/gravity/centrifugal c - potentials, and their 1st derivatives wrt N,E,U, without c - resorting to Clenshaw summation. c - jnorm contains CONVENTIONAL j values, which we will stick c - with, and not go to fully nor quasi normalized values. c - Note: d/dphi Pn(t) = Pn1(t) c - Note: d/dNorth = (1/r) * d/dphi real*8 jnorm(njmax),j0,pleg(njmax*2+1) common/block1/pi,d2r j0 = 1.d0 cphi = dcos(xlat*d2r) sphi = dsin(xlat*d2r) V2 = 0.d0 V2N = 0.d0 V2E = 0.d0 V2U = 0.d0 V2 = 0.d0 PHIN = 0.d0 PHIE = 0.d0 PHIU = 0.d0 U = 0.d0 UN = 0.d0 UE = 0.d0 UU = 0.d0 scale = (gmnorm/r) sum = 0.d0 do 1 i=1,njmax n = 2*i xxx = (anorm/r)**n * jnorm(i) V2 = V2 + xxx * pleg(n+1) V2N = V2N + xxx * pnm(n,1,xlat) V2U = V2U + xxx * (dble(-n-1)/r) * pleg(n+1) 1 continue V2 = scale * ( 1.d0 - V2) V2N = (1.d0/r) * scale * ( 0.d0 - V2N) V2U = scale * ((-1.d0/r) - V2U) PHI = 0.5d0*omega*omega*r*r*cphi*cphi PHIN = (1.d0/r) * (-omega*omega*r*r*cphi*sphi) PHIU = omega*omega*r*cphi*cphi U = V2 + PHI UN = V2N + PHIN UU = V2U + PHIU return end c c c double precision function gamh(a,b,gm,omega,h,xlat,gamma) c - Function to compute, through a closed form, the c - normal gravity vector magnitude at an elevation h. c - This closed form was tested 1/29/97 by D. Smith c - against both Clenshaw summation, and pure Legendre c - summation and found to agree to about 0.01 mgals at c - various elevations (up to 3000 meters). The Clenshaw c - and pure Legendre summations agreed to 16 decimal places. c - This function based on equation 2-123 of H/M, p. 79. implicit real*8(a-h,o-z) real*8 m common/block1/pi,d2r f = (a-b)/a m = omega*omega*a*a*b/gm gamh = gamma*(1.d0 - (2.d0/a)*(1.d0 + f + m * - 2*f*dsin(xlat*d2r)**2)*h * + (3.d0/a/a)*h*h) return end c c c subroutine fixab2(anorm,bnorm,tidnorm,tidsum,afixed,bfixed) c - Subroutine to convert anorm/bnorm (a and b of the c - ellipsoid in the "tidnorm" tide system, to afixed/bfixed c - (a and b of the ellipsoid in the "tidsum" tide system). c - That is: transform the normal field ellipsoid to the tide c - system requested for output. implicit real*8(a-h,o-z) character*4 tidnorm,tidsum real*8 kto,kfrom write(6,3)tidsum write(6,4)tidnorm 3 format(1x,'The tide system you requested for output is: ',a4, *' tide') 4 format(1x,'You gave the normal field in the: ',a4, *' tide system,',/,1x,'requiring a transformation of a & b') if(tidnorm.eq.'non-')then write(6,1) 1 format(1x,'To transform