c-g1sres.f c--removes tides and pressure from 1 sec files c--modified from g5sres.f - 21 April 97 c-------------------------------------------------------------------------- program G1SRES implicit real*8(a-h,o-y),integer*4(i-n),logical(z) parameter(ndim=180000) dimension g(ndim),gsyn(ndim),p(ndim), & p10(9000),itime(10000),itt(6),pp(6),y(10) real*4 t0,t1,secnds character*40 fname1,fname2,fname3,fname11,fname(400) character*20 fmt,fmtr,dunitg,tunit character*5 fday character*3 fext,fexr character*1 ans common /flags/ir,id,ifcn common // g,gsyn,p c c--initialise with values for Cantley data np,delt/86400,1.0/, & dunitg,tunit,fmt/'microgal','sec','(6f12.6)'/, & dlat,dlon/45.5850,-75.8071/,nhour,dut/0,0.0166667/, & amin/1.d-6/,min,nsec/2*0/,fmtr/'(5g15.7)'/, & gcal,pcal/-63.94d0,60.82105d0/,tdelay/1.54/ c c lenfm = lstr(fmt,-20) c--common block parameters for calculating the synthetic tides ir = 0 ifcn = 1 id = 0 c c--files opened : c 1 = input gravity data (input) c 2 = input pressure data (input) c 3 = output for tides (and/or pressure) residuals c 4 = qinwen.dat (input) c 8 = temporary output c 11 = synthetic tide (output) c 12 = canalyse.dat (local delta factors for Cantley) c c--choose which days to be filtered write(*,'(/a)') ' Program G1SRES ' write(*,'(a)') ' -------------------------------------------' c c--choose which days to be processed write(*,'(a$)') ' Enter first day to be treated --> ' read(*,*) iday1 write(*,'(a$)') ' ..... last ................... --> ' read(*,*) iday2 c c--tides out? write(*,'(a$)') ' Synthetic tides output [y/n] --> ' read(*,'(a1)') ans zstide = .false. if(ans.eq.'y') zstide = .true. c c--include pressure? write(*,'(a$)') ' Include pressure [y/n] --> ' read(*,'(a1)') ans zpress = .false. if(ans.eq.'y') zpress = .true. c c--local tidal delta factors? write(*,'(a$)') ' Local delta factors [y/n] --> ' read(*,'(a1)') ans zdelta = .false. if(ans.eq.'y'.or.ans.eq.'Y') zdelta = .true. c--if local tidal factors, assume rigid earth, no fcn and no time delay! if(zdelta) then ir = 1 ifcn = 0 tdelay = 0.d0 endif c--tidal accuracy write(*,'(a)') ' Doodson cut-off' write(*,'(a)') ' 1 amin = .000001 sums 3070 waves' write(*,'(a)') ' 2 amin = .00001 sums 1079 waves' write(*,'(a)') ' 3 amin = .0001 sums 383 waves' write(*,'(a)') ' 4 other ' write(*,'(a$)') ' choose --> ' read(*,*) iamin if(iamin.eq.1) then amin = 0.000001 elseif(iamin.eq.2) then amin = 0.00001 elseif(iamin.eq.3) then amin = 0.0001 else write(*,'(a$)') ' enter cutoff value --> ' read(*,*) amin endif c c--overwrite existing files? write(*,'(a$)') ' Overwrite existing files [y/n] --> ' read(*,'(a1)') ans zwrite = .false. if(ans.eq.'y') zwrite = .true. write(*,'(a)') ' -------------------------------------------' c c--set the filename extensions fext = 'gss' if(zdelta) fext = 'gls' fexr = 'rts ' if(zdelta) fexr = 'tls' if(zpress) fexr = 'tps' if(zpress.and.zdelta) fexr = 'pls' c c--try to find files for the days specified kk = 0 kday = iday1 c c--convert kday to a character string 'fday ' 5 call intch(kday,fday) c--then transfer fday to the input filename call fnday(fname1,'g',fday,'s14') c--check to see if this file exists inquire(file = fname1,exist = zexist) if(zexist) goto 8 c--if not, try the .s13 file call fnday(fname1,'g',fday,'s13') inquire(file = fname1,exist = zexist) if(zexist) goto 8 c--if not, try the .s12 file call fnday(fname1,'g',fday,'s12') inquire(file = fname1,exist = zexist) if(zexist) goto 8 c--if not, use the original gravity file call fnday(fname1,'g',fday,'s11') inquire(file = fname1,exist =zexist) if(zexist) goto 8 write(*,'(a,i5)') ' No input file found for day ',kday goto 10 c c--copy the filename to the array 8 kk = kk+1 fname(kk) = fname1 c write(*,*) kday,fday,fname1 c c--get another day 10 kday = kday+1 if(kday.le.iday2) goto 5 c c--there are kmax days to be processed kmax = kk if(kmax.eq.0) stop k = 1 c c--select a day from the list 15 fname3 = fname(k) t0 = secnds(0.0) c c--give new filename extension and extract fday (number of the day) call fnext2(fname3,fday,fexr) c c--check for overwriting file if(.not.zwrite) then inquire(file=fname3,exist=zexist) if(zexist) then write(*,*) ' skipping file: ',fname3 goto 800 endif endif c c--read the gravity data, note the data structure write(*,*) ' reading: ',fname(k) open(1,file = fname(k)) ja = 0 20 ja = ja+1 read(1,*) itime(ja),y do 25 i=1,10 jy = 10*(ja-1)+i g(jy) = y(i) 25 continue if(itime(ja).lt.(np-10)) goto 20 close(1) c c--get the correct kday call chint(fday,kday) c c--convert kday to year and # days since beginning of year kkday = kday-1000*(kday/1000) iyr = (kday-kkday)/1000+1900 xday = float(kkday) c c--get the date in the right format c--adjust ut1 for instrument phase lag, converting sec to hr ut1 = tdelay/3.6d2 call jdate2(iyr,mth,nday,nhour,min,nsec,xday,ut1,2) write(*,333) kday,iyr,mth,nday,ut1 333 format(' tides : ',i5,' year ',i4,' month ',i2, & ' day ',i3,' time ',f5.3) c c--compute the local synthetic tide using GWAVE call gwaves(dlat,dlon,iyr,mth,nday,ut1,np,dut,amin,zdelta) c c--save synthetic tides if(zstide) then fname11 = fname(k) call fnext2(fname11,fday,fext) write(*,*) ' writing: ',fname11 open(11,file=fname11) if(zdelta) then write(11,'(2a)') 'local synthetic tide for day ',fday else write(11,'(2a)') ' synthetic tide for day ',fday endif write(11,'(i10,g12.5)') np,delt write(11,'(3a20)') fmt,dunitg,tunit write(11,fmt) (gsyn(i),i = 1,np) close(11) endif c c--find and subtract mean values from gravity and tides call meanv(g,np,gav,1) call meanv(gsyn,np,gsav,1) if(.not.zpress) goto 206 c c--look for a pressure file corresponding to the current day fname2 = fname(k) c--get pressure filename, assuming a manually patched file call fnday(fname2,'p',fday,'st3') inquire(file = fname2,exist = zexist) if(zexist) goto 212 c--if not, try a gap-patched file call fnday(fname2,'p',fday,'st2') inquire(file = fname2,exist = zexist) if(zexist) goto 212 c--if not, use the original pressure file call fnday(fname2,'p',fday,'st1') inquire(file = fname2,exist = zexist) if(zexist) goto 212 write(*,*) ' pressure file not found for day ',fday stop c c--read in pressure data, note data structure 212 write(*,*) ' reading: ',fname2 open(2,file = fname2) ja = 0 202 ja = ja+1 read(2,*) itp,(itt(j),pp(j),j=1,6) do 203 i = 1,6 jy = 6*(ja-1)+i p10(jy) = pp(i) 203 continue if(itp.lt.(np-60)) goto 202 close(2) c c--interpolate the 10 sec values linearly to 1 sec do 205 i = 1,np if(i+10.lt.np) then ji = 1+(i/10) ki = mod(i,10) p0 = p10(ji) p(i) = p0+float(ki)*(p10(ji+1)-p0)/10.0 else p(i) = p10(np/10) endif 205 continue c c--subtract mean value of pressure and convert to mbar call meanv(p,np,pav,1) do 305 i = 1,np p(i) = pcal*p(i) 305 continue c c--PART 1 - Remove tides from observed gravity using simple subtraction c--as in the program GTRES 206 continue do 320 i = 1,np g(i) = g(i)*gcal-gsyn(i) 320 continue c c--write tide corrected gravity if(.not.zpress) then write(*,*) ' writing: ',fname3 open(3,file = fname3) write(3,'(2a)') ' residual gravity for ',fday write(3,'(i10,g12.5)') np,delt write(3,'(3a20)') fmtr,dunitg,tunit write(3,fmtr) (g(i),i=1,np) close(3) goto 800 endif c c--PART 2 - Remove pressure from residual gravity c c--correct gravity for pressure do 340 i = 1,np g(i) = g(i)+0.3d0*p(i) 340 continue c c--write the tide and pressure corrected gravity t1 = secnds(t0) write(*,*) ' writing: ',fname3,' time =',t1,' sec' open(3,file = fname3) write(3,'(2a)') ' residual gravity for ',fday write(3,'(i10,g12.5)') np,delt write(3,'(3a20)') fmtr,dunitg,tunit write(3,fmtr) (g(i),i=1,np) close(3) c c--get another day 800 k = k+1 if(k.le.kmax) goto 15 c 999 stop end c------------------------------------------------------------------------------ subroutine meanv(y,n,ymean,kk) c c--find and possibly remove mean value from a data set c--kk = 0 find mean c--kk = 1 find and remove mean c implicit real*8(a-h,o-z),integer*4(i-n) dimension y(1) c ymean = 0.d0 do 10 i = 1,n ymean = ymean + y(i) 10 continue ymean = ymean/float(n) if(kk.ne.1) return c do 20 i = 1,n y(i) = y(i)-ymean 20 continue c return end c---------------------------------------------------------------------------- subroutine fnday(fname,p,fday,ext) c c--creates a filename fname from prefix 'p' plus c--root 'fday' plus extension 'ext' c implicit real*8(a-h,o-y),integer*4(i-n),logical(z) character*40 fname,fnamec character*5 fday,fdayc character*3 ext,extc character*1 p,dum(40),dum3(3),dum5(5) equivalence (fnamec,dum),(extc,dum3),(fdayc,dum5) c--length of prefix is 1, fday is 5 and extension 3, total 10 fdayc = fday extc = ext dum(1) = p do 10 i = 2,6 dum(i) = dum5(i-1) 10 continue dum(7) = '.' dum(8) = dum3(1) dum(9) = dum3(2) dum(10) = dum3(3) do 20 i = 11,40 dum(i) = ' ' 20 continue fname = fnamec return end c--------------------------------------------------------------------- subroutine fnext2(fname2,fday,ext) c c changes the name of file fname2 c gives extension '.ext', where 'ext' is provided on input c implicit real*8(a-h,o-y),integer*4(i-n),logical(z) character*40 fname2,fnamec character*5 fday,fdayc character*3 ext,extc character*1 dum(40),dum3(3),dum5(5) equivalence (fnamec,dum),(extc,dum3),(fdayc,dum5) c c original length of filename is 'len0' c find length of extension 'len3' by searching last 4 characters c fnamec = fname2 len0 = lstr(fnamec,-40) len3 = 0 ltmp = 0 do 5 i = 1,4 l = len0-i+1 if(l.eq.0) goto 6 if(dum(l).eq.'.') then len3 = ltmp+1 goto 6 endif ltmp = ltmp+1 5 continue c c set fday to first segment name minus the 'g' c 6 do 10 i = 1,5 dum5(i) = dum(i+1) 10 continue fday = fdayc c c add extension '.ext' c 20 extc = ext len2 = len0-len3 dum(len2+4) = dum3(3) dum(len2+3) = dum3(2) dum(len2+2) = dum3(1) dum(len2+1) = '.' fname2 = fnamec c return end c------------------------------------------------------------------ subroutine intch(iday,fday) c c--convert integer 'iday' to character string 'fday' c implicit real*8(a-h,o-y),integer*4(i-n),logical(z) character*5 fday,fdayc character*1 dum(5) equivalence (fdayc,dum) c n = iday do 10 i = 1,5 k = n-10*(n/10) dum(6-i) = char(k+48) n = n/10 10 continue fday = fdayc return end c----------------------------------------------------------------------------- subroutine chint(fday,kday) c c--converts a 5-character string to an integer c implicit real*8(a-h,o-y),integer*4(i-n),logical(z) character*5 fday,fdayc character*1 dum(5) equivalence (fdayc,dum) c fdayc = fday kday = 0 do 10 i = 1,5 k = ichar(dum(6-i))-48 kday = kday+ k*(10**(i-1)) 10 continue c return end c------------------------------------------------------------------------ subroutine jdate2(year,month,day,hour,min,sec,xday,ut1,k) c c converts between date (year, month, day), time (hour, min, sec) c and days in year (xday) and hours (ut1) c k = 1 for date --> ut1 : input year, month, day, hour, min, sec c k = 2 for ut1 --> date : input year, xday, ut1 c implicit real*8(a-h,o-z),integer*4(i-n) integer year,month,day,hour,min,sec,nday(12) logical leapyr data nday/31,59,90,120,151,181,212,243,273,304,334,365/ c c-- check for a leap year leapyr=.false. c--is year a century? if(mod(year,100).eq.0) then c--has to be divisible by 400 if(mod(year,400).eq.0) leapyr=.true. else c--otherwise has to be divisible by 4 if(mod(year,4).eq.0) leapyr=.true. endif if(k.eq.1) then c--date --> ut1 c--calculate # days since beginning of year xday = day if(month.gt.1) xday = xday+nday(month-1) c--add a day for a leap year if(leapyr.and.month.ge.2) xday = xday+1 c--find ut1 in decimal hours ut1 = float(hour)+(float(min)+float(sec)/60.)/60. else c--ut1 --> date c--initialise day iday = aint(xday) day = iday c--check each month do 10 month=1,12 mday = nday(month) c--add a day for a leap year if(leapyr.and.month.ge.2) mday = mday+1 c write(*,*) month,mday,iday if(iday.le.mday) goto 20 day = iday-mday 10 continue c--get time to nearest second 20 hour = aint(ut1) remaind = (ut1-hour)*60. min = aint(remaind) remaind = (remaind-min)*60. sec = anint(remaind) endif c return end c---------------------------------------------------------------------------- function lstr(name,n) c c function to return the length of a string. On input the string c 'name' has a nominal length n; if n is negative, blank characters c at the beginning of 'name' are first truncated. c character*1 name(1) c c--zero length case first if(n.eq.0) then lstr = 0 return endif k = 0 ld = n c c--ld is negative, we truncate leading blanks if(ld.lt.0) then ld = -ld j = 0 10 j = j+1 if(name(j).eq.' ') goto 10 nb = j-1 c c--there are nb blanks to be erased in front of 'name' if(nb.gt.0) then do 20 i = nb+1,ld name(i-nb) = name(i) name(i) = ' ' 20 continue c c--reduce length of string by number of blanks ld = ld-nb endif endif c c--ld is non-zero, find length of string c--first set true length equal to nominal length k = ld+1 c--and decrease k until a non-blank character appears or beginning c--of string is reached 40 k = k-1 if(k.gt.0.and.name(k).eq.' ') goto 40 lstr = k c return end c---------------------------------------------------------------------------- subroutine gwaves(lat,lon,iyr,mth,nday,ut1,nobs,dut,amin,zdelta) c c--subroutine for gwave driver, note similarity of input to gtide c--additional (lower case) code by D. Crossley, 3 April 1990, modified c--from Merriam (94/02/14) c c inputs : lat = station latitude c long = station longitude c iyr = calendar year c mth = month of year c nday = day of month c ut1 = time in decimal hours since beginning of day c nobs = number of readings required c dut = time interval (mins) between each value c amin = minimum cutoff for Doodson amplitude c zdelta = .true. if local tide c---------------------------------------------------------------------------- c C ************************************************************************ C GWAVE IS A GRAVITY TIDE PREDICTION PROGRAM BASED ON THE XI QINWEN SERIES C OF 3070 TIDAL CONSTITUENTS. IT WAS WRITTEN BY J. MERRIAM AT THE UNIV. OF C SASKATCHEWAN TO REMOVE THE GRAVITY TIDE SIGNAL FROM SUPERCONDUCTING C GRAVITY METERS. ITS PREDICTIONS ARE ACCURATE TO ABOUT THE 40 NANOGAL C LEVEL IN THE DIRECT ATTRACTION OF THE SUN AND MOON. C C THE RESPONSE OF THE EARTH IS COMPUTED WITH LATITUDE AND FREQUENCY C DEPENDENT GRAVIMETRIC FACTORS FOLLOWING WAHR (1981), SO THAT THE GRAVITY C PREDICTIONS ARE ON THE ELLIPSOID AT GEODETIC LATITUDE AND LONGITUDE, AND C NORMAL TO THE ELLIPSOID. C C ANY QUESTIONS ARISING FROM THE USE OF THIS PROGRAM SHOULD BE C DIRECTED TO "JAMES@SKERTH.USASK.CA". C ************************************************************************ C C JAN 15 1990 C C ************************************************************************ C implicit real*8(a-h,o-y),integer*4(i-n),logical(z) dimension deltan(47),xkappa(47),omegan(47) real*8 lat,lon,long common /toto/deltan,xkappa,omegan c C ********************************************************************** C THE DOODSON AMPLITUDE IS A RELATIVE MEASURE OF THE AMPLITUDE OF THE C CONSTITUENTS. IT IS THE LAST COLUMN IN THE DATA FILE GWAVE.DAT/1000000. C RESPONDING TO THE NEXT PROMPT WITH .000001 WILL RESULT IN A SUM OVER C ALL 3070 WAVES AND A NOMINAL PRECISION OF A NANOGAL. INCREASING THIS C BY A FACTOR OF 10, WILL DEGRADE THE NOMINAL PRECISION BY A FACTOR OF C ABOUT 10, BUT WILL RESULT IN A FASTER RUN TIME, AS FEWER CONSTITUENTS C WILL BE SUMMED. C ********************************************************************** c c--initialise, reading the appropriate file if local delta factors if(zdelta) call lire pi = 4.d0*datan(1.d0) colat = 90.d0-lat LONG = lon*PI/180.D0 COLAT = COLAT*PI/180.D0 c--the next line converts colatitude from geographic to geocentric COLAT = COLAT+0.003352811*DSIN(2.D0*COLAT) c C *********************************************************************** C CALCULATE TD1, THE TIME IN DAYS FROM 12 HRS JAN 1 2000 TO O HRS ON THE C THE FIRST DAY FOR WHICH A PREDICTION IS DESIRED, ALLOWING FOR LEAP YEARS, C AND FOUR CENTURY LEAP YEARS. C *********************************************************************** TD1=INT(365.25*IYR)+INT(30.6001*(MTH+1))+INT(IYR/400)-INT(IYR/100) IF(MTH.LE.2) TD1=INT(365.25d0*(IYR-1))+ & INT(30.6001d0*(MTH+13))+INT(IYR/400)-INT(IYR/100) TD1=TD1+1720996.5D0-2451545.0D0+nDAY CALL TIDE(TD1,LONG,COLAT,AMIN,NOBS,DUT,UT1,zdelta) c return END c------------------------------------------------------------------------------ SUBROUTINE TIDE(TD1,LONG,COLAT,AMIN,NOBS,DUT,UT1,zdelta) c parameter(ndim=180000) implicit real*8(a-h,o-y),integer*4(i-n),logical(z) dimension gg(ndim),gt(ndim),pp(ndim), & FREQ(3070),AMP(3070),PHASE(3070) REAL*8 LONG c CHARACTER*7 DOODNO common /flags/ir,id,ifcn common // gg,gt,pp c c--initialise C DOOD IS DOODSON'S CONSTANT, AS REDEFINED BY XI QINWEN data DOOD/2.6335838D+04/ pi = 4.d0*datan(1.d0) degrad = pi/180.d0 c C CALCULATE TJ, THE TIME IN JULIAN CENTURIES FROM J2000.0 TJ=TD1/36525.D0 c C CALCULATE ETUT, THE DIFFERENCE BETWEEN ET. AND UT, WHICH MUST BE C ADDED TO TJ TO PRODUCE A PREDICTION FOR UT. ETUT IS CALCULATED AS C A LINEAR RATE OF 1 SEC/YR FROM JAN 1 1991. ETUT=(TJ*100.d0+9.d0)+32.1840d0+26.d0 IF(ETUT.EQ.0.D0) WRITE(*,*)'WARNING: THIS IS ET NOT UT!!' TJ=TJ+ETUT/(3600.D0*24.D0*36525.D0) c C CALCULATE THE RIGHT ASCENSION OF THE FICTITIOUS MEAN SUN IN UT R=18.6973746d0+2400.0513369d0*TJ+0.0000258622d0*TJ**2 & -1.7222D-09*TJ**3 R=R*15.D0 c C CALCULATE THE FUNDAMENTAL ARGUMENTS OF THE TIDAL CONSTITUENTS, C INCLUDING ADDITIVE TERMS DUE TO PLANETARY INTERACTIONS. C C XN - THE NEGATIVE OF THE LONGITUDE OF THE MEAN ASCENDING LUNAR NODE. C S - THE GEOMETRIC MEAN LONGITUDE OF THE MOON. C P1 - THE MEAN LONGITUDE OF SOLAR PERIGEE. C P - THE MEAN LONGITUDE OF LUNAR PERIGEE. C H - THE GEOMETRIC MEAN LONGITUDE OF THE SUN C XN=-125.04452D0+1934.13626D0*TJ-0.00207D0*TJ**2-0.000002D0*TJ**3 XN=XN-0.02666D0*DSIN(-XN*PI/180.D0) &-0.00433D0*DSIN((-XN+272.75-2.30*TJ)*PI/180.D0) &-0.00042D0*DSIN((188.83-1934.14D0*(TJ+1))*DEGRAD) &-0.00017D0*DSIN((51.19+20.20*(TJ+1))*DEGRAD) &-0.00005D0*DSIN((193.44+20.20*(TJ+1))*DEGRAD) S=218.31643D0+481267.88128D0*TJ-0.00161D0*TJ**2+0.000005D0*TJ**3 C ADDITIVE TERMS (NO PLANETARY TERMS ARE INCLUDED WHICH MIGHT RESULT IN C A THIRTY NANOGAL ERROR &+0.00396D0*DSIN((60.57-132.87*TJ)*PI/180.D0) &+0.00202D0*DSIN(-XN*PI/180.D0) &+0.000233*DSIN((51.19+20.20*(TJ+1))*DEGRAD) &+0.000086*DSIN((84.11+16.2*(TJ+1))*DEGRAD) &+0.000078*DSIN((174.23-1936.44*(TJ+1))*DEGRAD) &+0.000065*DSIN((304.33-150.679D0*(TJ+1))*DEGRAD) &+0.000035*DSIN((235.96-1034.10D0*(TJ+1))*DEGRAD) &+0.000030*DSIN((145.27-287.55D0*(TJ+1))*DEGRAD) &+0.000011*DSIN((230.91-6.08*(TJ+1))*DEGRAD) &+0.000011*DSIN((211.50+119.0*(TJ+1))*DEGRAD) &+0.000021*DSIN((76.93+916.72D0*(TJ+1))*DEGRAD) &+0.000008*DSIN((231.63+4311.54D0*(TJ+1))*DEGRAD) &+0.000009*DSIN((275.53-88.34*(TJ+1))*DEGRAD) &+0.000015*DSIN((276.67-723.02D0*(TJ+1))*DEGRAD) &+0.000003*DSIN((113.15+619.69D0*(TJ+1))*DEGRAD) &+0.000004*DSIN((269.69-507.65D0*(TJ+1))*DEGRAD) &+0.000004*DSIN((188.12-658.33D0*(TJ+1))*DEGRAD) &+0.000007*DSIN((46.75-890.44D0*(TJ+1))*DEGRAD) &+0.000005*DSIN((10.59+345.70D0*(TJ+1))*DEGRAD) P1=282.93835D0+1.71946D0*TJ+0.00046D0*TJ**2+0.000003D0*TJ**2 P=83.35345D0+4069.01388D0*TJ-0.01031D0*TJ**2-0.00001D0*TJ**3 &-0.00058D0*DSIN((71.40+20.20*TJ)*PI/180.D0) &-0.00058D0*DSIN(-XN*PI/180.D0) &-0.00023D0*DSIN((174.23-1936.44*(TJ+1))*DEGRAD) &-0.00016D0*DSIN((304.33-150.68*(TJ+1))*DEGRAD) &-0.00003D0*DSIN((193.44-132.87*(TJ+1))*DEGRAD) &-0.00003D0*DSIN((211.50+119.00*(TJ+1))*DEGRAD) H=280.46607D0+36000.76980D0*TJ+0.00030D0*TJ**2 C ADDITIVE TERMS IN H H=H+0.00178D0*DSIN((251.39+20.20*TJ)*PI/180.D0) &+(0.0005183D0-0.0000044D0*TJ)*DSIN((207.51+150.27*TJ)*PI/180.D0) &+0.0000738*DSIN((31.8+119.0*(TJ+1))*DEGRAD) &+0.0000561202*DSIN((315.6+893.3*(TJ+1))*DEGRAD) c write(2,'(5g15.7)') s,h,p,p1,xn C PLANETARY TERMS IN H C (ONLY THOSE RECOMMENDED BY XI ARE INCLUDED. OMITTED TERMS CAN SUM TO C 17 SECONDS OF ARC, OR 30 NANOGAL) c &-0.00200D0*DSIN((67.20+32964.47*TJ)*PI/180.D0) c &-0.00154D0*DSIN((16.85-45036.89*TJ)*PI/180.D0) c &+0.00134D0*DSIN((81.51+22518.44*TJ)*PI/180.D0) C THIS IS A SECULAR TERM IN ECCENTRICITY c &-0.00479D0*(TJ+1)*(1.E0+0.003*(TJ+1))*DSIN((H-P1)*PI/180.D0) C THIS IS THE LUNAR INEQUALITY c &+0.00179D0*DSIN((S-H)*PI/180.E0) C CALCULATE XNUT, THE NUTATION, AND ADD TO ALL THE FUNDAMENTAL ARGUMENTS. c XNUT=-0.00478*DSIN(-XN*PI/180.E0)-0.00037*DSIN(2.D0*H*PI/180.D0) c S=S+XNUT c H=H+XNUT c P=P+XNUT c P1=P1+XNUT c XN=XN-XNUT c EQEQ=XNUT*DCOS(23.43939*PI/180.D0) c R=R+EQEQ C ********************************************************************** C READ IN THE DATA IN GWAVE.DAT. THESE ARE; C C ITAU - AN INTEGER MULTIPLE OF MEAN LUNAR TIME. C NS,NH,NP,NN,NP1 - INTEGER MULTIPLES OF THE ASSOCIATED FUNDAMENTAL C ARGUMENTS. C N,M - DEGREE AND ORDER OF THE CONSTITUENT. C OMEGA - THE FREQUENCY OF THE CONSTITUENET IN DEGREES PER SOLAR HOUR. C IAMP - AN INTEGER BETWEEN +/-1 AND 908175 SPECIFYING THE RELATIVE C AMPLITUDE OF THE CONSTITUENTS. FOR EXAMPLE, THE ENTRY FOR THE C M2 WAVE (THE 465 TH ENTRY IN THE FILE) IS C C 2 0 0 0 0 0 255.555 22 28.98410424 56 75 154 465 908175 C C THE DOODSON NUMBER (255.555) IS NOT READ IN, AS IT IS NOT NEEDED, AND C NEITHER ARE THE FOUR PENULTIMATE COLUMNS AS THEY ONLY INDICATE THE C ROW NUMBER OF THE CONSTITUENT IN A TABLE WITH VARIOUS AMPLITUDE CUTOFFS. C ************************************************************************* c OPEN(UNIT=4,file='qinwen.dat') II=0 DO 1 I=1,3070 READ(4,64) ITAU,NS,NH,NP,NN,NP1,N,M,OMEGA,IAMP c C COMPARE THE DOODSON AMPLITUDE (IAMP*1.D-06) WITH THE CUTOFF, AMIN C AND DO NOT SUM THE CONSTITUENT IF ITS AMPLITUDE IS BELOW THE CUTOFF. c AM=IAMP*1.D-06 IF(DABS(AM).LT.AMIN)GO TO 1 II=II+1 FREQ(II)=OMEGA*PI/180.D0 c C CALCULATE THE COEFFICIENT OF THE GEOGRAPHIC DEPENDANCE OF TIDAL GRAVITY C FOR EACH DEGREE (N) AND ORDER (M). (SEE NOTE IN SUBROUTINE GFACT) c GO TO (1,3,4,5)N 3 IF(M.EQ.0)XL=-1.D0 IF(M.EQ.1)XL=2.D0/3.D0 IF(M.EQ.2)XL=1.D0/3.D0 GO TO 7 4 IF(M.EQ.0)XL=-1.11803*2.D0 IF(M.EQ.1)XL=-0.72618*2.D0/3.D0 IF(M.EQ.2)XL=2.59808/15.D0 IF(M.EQ.3)XL=1.D0/15 GO TO 7 5 IF(M.EQ.0)XL=0.12500*8.D0 IF(M.EQ.1)XL=-0.47346*4.D0/5.D0 IF(M.EQ.2)XL=-0.77778*2.D0/15.D0 IF(M.EQ.3)XL=3.07920/105.D0 IF(M.EQ.4)XL=1.D0/105.D0 7 CONTINUE c C CALL GFACT TO GET THE GRAVIMETRIC FACTOR. c c CALL GFACT(DELTA,N,M,OMEGA,COLAT) c c--new gfact returns also the phase delay phi due to the FCN c call gfact(delta,deltae,phi,n,m,omega,colat) if(zdelta) call inter(omega,deltap,phn) c write(2,'(4g15.7)') omega,delta,deltap,phn c C CALCULATE THE AMPLITUDE OF EACH SUMMED CONSTITUENT (IN GALS), AND C ITS PHASE AT 0 HRS UT ON THE FIRST DAY FOR WHICH A PREDICTION IS DESIRED. AMP(II)=DELTA*DFLOAT(IAMP)*XL*DOOD*DFLOAT(N)/6.378137D+08 c--multiply amplitude by delta factor from inter if(zdelta) amp(ii) = amp(ii)*deltap PHA=NH*H+R*ITAU+(NS-ITAU)*S+NP*P+NN*XN+NP1*P1+M*LONG*180.D0/PI c--add phase shift due to FCN if(ifcn.eq.1) pha=pha+phi*180.d0/pi IX=(N+M)/2 X=DFLOAT(IX) Y=DFLOAT(N+M)/2.D0 IF(X.NE.Y)PHA=PHA-90.D0 PHA=PHA*PI/180.D0 C C IF THE AMPLITUDE OF THE WAVE IS NEGATIVE, CHANGE ITS C SIGN AND SUBTRACT PI FROM THE PHASE. C IF(AMP(II).LT.0.D0) THEN AMP(II)=-AMP(II) PHA=PHA-PI ENDIF PHASE(II)=DATAN2(DSIN(PHA),DCOS(PHA)) c--add phase due to local kappa factors if(zdelta) phase(ii)=phase(ii)+phn C C--CALL OCEAN TO ADJUST AMPLITUDE AND PHASE FOR TIDAL LOADING c IF(LOAD.EQ.'Y'.OR.LOAD.EQ.'y')THEN c CALL OCEAN(FREQ(II),AMP(II),DELTAE,PHASE(II),NOPT) c ENDIF c write(2,'(i10,3g15.7)') ii,xl,amp(ii),phase(ii) 1 CONTINUE close(4) c WRITE(*,65) II c 65 FORMAT(I5,' WAVES INCLUDED IN SUMMATION') 64 FORMAT(6I2,10X,2I1,F12.8,20X,I8) c C CALCULATE THE TIME (IN HRS) FROM 0 HRS ON THE FIRST DAY, TO THE TIME OF EACH C PREDICTION, AND ADVANCE THE PHASE OF EACH CONSTITUENT TO THIS TIME. c DO 2 I=1,NOBS DATE=UT1+DUT*FLOAT(I-1)/60.D0 G=0.D0 DO 6 K=1,II ARG=(FREQ(K)*DATE+PHASE(K)) C SUM OVER ALL CONSTITUENTS WITH AMPLITUDE ABOVE THE CUTOFF G=G-AMP(K)*DCOS(ARG) 6 CONTINUE gt(i) = g 2 CONTINUE c RETURN END c---------------------------------------------------------------------------- subroutine gfact(delta,deltae,phi,n,m,omega,theta) c c--new version of GFACT from Jim Merriam adapted to conform with c--existing versions of GWAVE, GTIDE by D. Crossley 93/06/02 c c SUBROUTINE GFACT(DELTAR,DELTAE,PNM,PHI,N,M,OMEGA,THETA) c C *********************************************************************** C GFACT RETURNS DELTAR, DELTAE, AND PNM. PNM IS THE DEGREE N ORDER M C ASSOCIATED LEGENDRE FUNCTION, DELTAR IS THE RIGID ELLIPSOID C GRAVIMETRIC FACTOR AND DELTAE IS THE ELASTIC ELLIPSOID GRAVIMETRIC C FACTOR. C C TIDAL GRAVITY ON RIGID EQUATORIAL SPHERE IS PNM*[] C WHERE []IS A QUANTITY CALCULATED BY GTIDE AND GWAVE. C C TIDAL GRAVITY ON RIGID ELLIPSOID = DELTAR*PNM*[] C C TIDAL GRAVITY ON ELASTIC ELLIPSOID = DELTAE*DELTAR*PNM*[] C C DELTAE INCLUDES A COMPLEX FCN C ************************************************************************* C C JAN 6 1993 C C ************************************************************************* C C THIS PROGRAM WAS WRITTEN BY J. MERRIAM AT THE UNIVERSITY OF SASKATCHEWAN C C PROBLEMS WITH THIS PROGRAM SHOULD BE ADDRESSED TO C C MERRIAM@GEOID.USASK.CA C C ************************************************************************* C IMPLICIT REAL*8(A-H,L,K,O-Z) COMPLEX*16 FX1 DIMENSION DELTAO(6) common /flags/ir,id,ifcn c DATA DELTAO/1.1563,1.1533,1.1589,1.0719,1.035,1.022/, & F/0.003352811/,O1/13.943036D0/,FCN/15.07590D0/,Q/7300/, & AR/-0.000631/ c C DELTAO(1-6) ARE DEGREE 2,0 2,1 2,2 3 4 5 GRAVIMETRIC FACTOR C FCN - FREQUENCY OF THE FREE CORE NUTATION (DEGREES/SOLAR HOUR). C 15.073727 DEG/HR = 1.004915 CY/DAY (THEORETICAL) C 15.07595 DEG/HR = 1.005063 CY/DAY (RICHTER) C 15.07571 DEG/HR = 1.005047 CY/DAY (NEUBERG) C 15.0758 DEG/HR = 1.005053 CY/DAY (MERRIAM) C THETA - COLATITUDE C PHI=0.D0 IF(N.EQ.2.AND.M.EQ.0)THEN DELTAE=DELTAO(1)-0.0016*DSQRT(9.D0/80.D0)* & (35.D0*DCOS(THETA)**4-30.D0*DCOS(THETA)**2+3.D0)/ & (3.D0*DCOS(THETA)**2-1.D0)+0.0053* & DSQRT(4.D0/5.D0)*(1.D0/(3.D0*DCOS(THETA)**2-1.D0)) ELSEIF(N.EQ.2.AND.M.EQ.1)THEN c--choose whether to make the FCN resonance correction if(ifcn.eq.1) then FX1=(DCMPLX(OMEGA-O1,0.D0)/DCMPLX(FCN-OMEGA,FCN/(2*Q))) else fx1 = dcmplx(0.d0,0.d0) endif FX=CDABS(FX1)*DSIGN(1.D0,DBLE(FX1)) DELTAE=DELTAO(2)+AR*FX+(-0.0015D0-1.06D-05*FX)* & DSQRT(3.D0/8.D0)*(7.D0*DCOS(THETA)**2-3.D0) PHI=DATAN2(DIMAG(FX1)*AR, & DELTAO(2)+DBLE(FX1)*AR+(-.0015D0-1.06D-05*FX)* & DSQRT(3.D0/8.D0)*(7.D0*DCOS(THETA)**2-3.D0)) ELSEIF(N.EQ.2.AND.M.EQ.2)THEN DELTAE=DELTAO(3)-0.0010D0*DSQRT(3.D0/4.D0)* & (7.D0*DCOS(THETA)**2-1.D0) ELSEIF(N.EQ.3.AND.M.EQ.3)THEN DELTAE=DELTAO(4)-0.0034*DSQRT(11.D0/16.D0)* & (9.D0*DCOS(THETA)**2-1.D0) ELSEIF(N.EQ.3.AND.M.NE.3)THEN DELTAE=DELTAO(4) ELSEIF(N.EQ.4)THEN DELTAE=DELTAO(5) ELSEIF(N.EQ.5)THEN DELTAE=DELTAO(6) ENDIF C XN=N XM=M ANM=(XN*(XN-XM+1)*(XN+1+XM)/((2*XN+1)*(2*XN+3)) &-(XN+1)*(XN-XM)*(XN+XM)/((2*XN+1)*(2*XN-1)))/XN XNM=0.D0 IF(N-2.LT.M)GO TO 3 XNM=(XN+XM)*(XN-1+XM)*(XN+1)/((2*XN+1)*(2*XN-1))* &PLGNDR(N-2,M,DCOS(THETA))/XN 3 CONTINUE PNP=(XN-XM+1)*(XN-XM+2)/((2*XN+1)*(2*XN+3))* &PLGNDR(N+2,M,DCOS(THETA)) PNM=PLGNDR(N,M,DCOS(THETA)) DELTAR=1.D0-2*F*(ANM-XNM/PNM+PNP/PNM)-(XN-1.D0)*F*DCOS(THETA)**2 c c--compute delta as if for a rigid ellipsoid and add elasticity as needed c delta = deltar*pnm if(ir.eq.0) delta = delta*deltae c write(2,'(5f15.7)') omega,deltar,deltae,phi,pnm c RETURN END c------------------------------------------------------------------------------- FUNCTION PLGNDR(N,M,X) IMPLICIT REAL*8(A-H,O-Z),integer*4(i-n) C C ADAPTED FROM NUMERICAL RECIPIES C C FUNCTION PLGNDR COMPUTES THE VALUE OF THE ASSOCIATED LEGENDRE C FUNCTION, AS DEFINED BY 12.65 AND 12.84 OF ARFKEN, AT ANY DEGREE C N AND ORDER M. C C THE SIGN OF THE ODD ORDER FUNCTIONS HAS BEEN REVERSED TO COMPLY C WITH ARFKEN. C IF(M.LT.0.OR.M.GT.N.OR.ABS(X).GT.1.D0)PAUSE 'BAD ARGS IN PLGNDR' PMM=1.D0 IF(M.GT.0)THEN SOMX2=DSQRT((1.D0-X)*(1.D0+X)) FACT=1.D0 DO 11 I=1,M PMM=PMM*FACT*SOMX2 FACT=FACT+2 11 CONTINUE ENDIF IF(N.EQ.M)THEN PLGNDR=PMM ELSE PMMP1=X*(2*M+1)*PMM IF(N.EQ.M+1) THEN PLGNDR=PMMP1 ELSE DO 12 LL=M+2,N PLL=(X*(2*LL-1)*PMMP1-(LL+M-1)*PMM)/(LL-M) PMM=PMMP1 PMMP1=PLL 12 CONTINUE PLGNDR=PLL ENDIF ENDIF RETURN END c-------------------------------------------------------------------------- subroutine lire c c--sous programme de lectude du fichier analyse.dat pour le calcul c--de la maree synthetique a l'aide de gwaves, mais ou l'on impose c--sois-meme les fact. delta et kappa c implicit real*8(a-h,o-y),integer*4(i-n),logical(z) dimension deltan(47),xkappa(47),omegan(47) common /toto/deltan,xkappa,omegan c c--initialise pi = 4.d0*datan(1.d0) open(12,file='canalyse.dat') c c--read coefficients and convert phase to radians do 10 j = 1,47 read(12,*) deltan(j),xkappa(j),omegan(j) c print*,j,deltan(j),xkappa(j),omegan(j) xkappa(j) = xkappa(j)*pi/180.d0 10 continue close(12) c return end c----------------------------------------------------------------------------- subroutine inter(omega,deln,phn) c c--programme pour l'interpolation lineaire de delta et kappa c implicit real*8(a-h,o-y),integer*4(i-n),logical(z) dimension deltan(47),xkappa(47),omegan(47) common /toto/deltan,xkappa,omegan c do 1 j=1,46 jj = j if(omegan(jj).gt.omega) goto 2 1 continue c 2 jj = jj-1 x1 = omegan(jj) x2 = omegan(jj+1) difx = x2-x1 deln = deltan(jj)+(deltan(jj+1)-deltan(jj))*(omega-x1)/difx phn = xkappa(jj)+(xkappa(jj+1)-xkappa(jj))*(omega-x1)/difx c print*,jj,x1,omega,x2,deltan(jj),deltan(jj+1) c return end