program tides c tide signal, U. Noell April 1993 c this program uses TIDEPT to compute a single component of the earth c tide signal on an elastic earth. c besides vital variables,format statements 108,109 should be modified dimension xn(120000),xe(120000),xnv(120000),xev(120000) dimension xz(120000),xzv(120000) double precision ts,dt,t,tim69 integer unit,year,month,day character prolog(5)*50,text*50 logical ex inquire(file='tides.par',exist=ex) if(ex.eq..true.) then open(1,file='tides.par') unit=1 else unit=5 endif open(unit=2,file='tides.out') write(*,*) write(*,*) '>> If file tides.par exists, then input data are read' write(*,*) '>> from it. Otherwise, input is interactive.' write(*,*) '>> Input can also be piped in with tides < inputfile' write(*,*) write(*,*) 'latitude?' read(unit,*) glat write(*,*) 'latitude=',glat write(2,*) 'latitude=',glat write(prolog(1),*) 'latitude=',glat write(*,*) 'longitude?' read(unit,*) glong write(*,*) 'longitude=',glong write(2,*) 'longitude=',glong write(prolog(2),*) 'longitude=',glong write(*,*) 'year,month,day (GMT)? ' read(unit,*) year,month,day write(*,*) 'year,month,day (GMT) =',year,month,day write(2,*) 'year,month,day (GMT) =',year,month,day write(prolog(3),*) 'year,month,day (GMT)=',year,month,day write(*,*) 'how many days? (sampled every quarter hour):' read(unit,*) ntage dt=0.25d0 ndata=min(120000,nint(24./dt*ntage)) write(*,*) ntage," days, ",ndata," data points. dt=15 minutes." write(2,*) ntage," days, ",ndata," data points. dt=15 minutes." write(prolog(4),*) & ntage," days, ",ndata," data points. dt=15 minutes." write(*,*) "time in hours GMT, gravity in microgals" write(2,*) "time in hours GMT, gravity in microgals" write(prolog(5),*) "time in hours GMT, gravity in microgals" write(2,*) &' time , z-component , n-component , e-component ' C glat= 48.7708 C glong= 9.19333 c east longitude in degrees azim=90.000 c azimuth in degrees C wink=2.74017 c strainmeter:east from north c tiltmeter:north from east,i.e. 0.0 degrees is east,90 degrees is c north tilt C opt=2 C opt=3 C opt=1 lorder=0 c gravity=1,tilt=2,strain=3 c dt=1.00d0 c timestep in hours ts=tim69(year,month,day) C ts= tim69(1991,7,25) c first point,time in hours from jan 01,1969,0.0hrs GMT c scale= 1.00e+06 c scale= 0.206264806e+09 scale= 1.00e+06 c if scale =1.0e+6 then gravity in microgals,tilt in microrads,strain c in microstrain.for nanorads and nanostrains resp. use 1.0e+09 c for tilt in msecs arc use 0.206264806e+09 C ndata = 32 t=ts do 1 i=1,ndata call tidept(t,xz(i),glat,glong,azim,scale,lorder,1) call tidept(t,xn(i),glat,glong,azim,scale,lorder,2) azm2 = 0.00 call tidept(t,xe(i),glat,glong,azm2,scale,lorder,2) C azm2 = 45.00 C call tidept(t,str1,glat,glong,azm2,scale,lorder,opt) C azm2 = -45.00 C call tidept(t,str2,glat,glong,azm2,scale,lorder,opt) C shear=(str1-str2)/2.00 C dil=(x(i)+str)*(0.66666) C senk=str*cos(wink)**2+x(i)*sin(wink)**2-shear*sin(2.*wink) C write(8,112) i,dil,senk t=t+dt 1 continue do 11 i=1,ndata xzv(i)=-xz(i) xnv(i)=-981.*xn(i) xev(i)=-981.*xe(i) write(2,112) i*dt,xzv(i),xnv(i),xev(i) 11 continue close(2) si=900. tmin=15. tsec=0. text='gravity tide EAST from program TIDES ' call outform('tide.e',text,prolog,ndata,si,tmin,tsec,xev) text='gravity tide NORTH from program TIDES ' call outform('tide.n',text,prolog,ndata,si,tmin,tsec,xnv) text='gravity tide VERTICAL from program TIDES ' call outform('tide.v',text,prolog,ndata,si,tmin,tsec,xzv) 900 format(' i east-strain north-strain shear-strain ') 901 format(' i nord-tilt east-tilt ') 902 format(' i dilatation normal strain ') 112 format(4(f12.4,2x)) c generate plot-parameter file open(4,file='winplot.par') write(4,*) '20, 3, 24, 15, 1,',ndata, ', 0.8' write(4,'(a)') 'tide.e' write(4,'(a)') 'tide.n' write(4,'(a)') 'tide.v' write(4,'(a)') 'synthetic tides from program TIDES' write(*,*) 'plot-parameter file winplot.par was generated' write(*,*) & 'The three time series were also written into file tides.out' stop end subroutine outform(name,text,prolog,n,dt,tmin,tsec,x) c formatted output dimension x(120000) character name*6,text*50,iform*20,prolog(5)*50 xmax=0. do 2 j=1,n 2 xmax=max(xmax,abs(x(j))) nvor=int(log10(max(xmax,2.))+1.) ndec=max(0,10-nvor) write(iform,3) ndec 3 format('(5f13.',i1,')') write(6,5) n,trim(name) 5 format(' writing ',i7,' samples into file ',a6,' in SEIFE format') open(8,file=trim(name)) write(8,'(a)') text do i=1,5 write(8,'(a5,a50)') '% ',prolog(i) enddo write(8,1) n,iform,dt,tmin,tsec 1 format(i10,a20,f10.3,2f10.3) write(8,iform) (x(j),j=1,n) close(8) return end Double Precision Function tim69(myr,mon,mday) c function to return the number of days from the beginning of 1969 c to a given later date. c intended to be used to find the input variable for the program TIDEPT c and its relatives,which need time in hours from January 1,1969, c 0.0 hours G.M.T. c this version goes wrong after 2100 and before 1800 c double precision st69,hours dimension iam(12) data iam/0,31,59,90,120,151,181,212,243,273,304,334/ lday=mday-1 kyr=myr-1969 lyr=kyr if(lyr.lt.0) lyr=lyr-3 st69=kyr*365+iam(mon)+lday+lyr/4 if(myr.eq.myr/4*4.and.mon.gt.2) st69=st69+1.0d0 if(myr.le.1900.and.mon.le.2) st69=st69+1.0d0 hours = st69*24.00d0 tim69=hours return end SUBROUTINE TIDEPT (TIME,RESULT,GLAT,GLONG,ORIENT,SCALE, * LORDER,OPT) C BERGER,EVANS,MCKENZIE C MODIFIED BY YOUNG 16.1.78 TO RETURN STRAIN VALUE AT ONE TIME C ONLY. SEE SUBROUTINE EARTHTD FOR DETAILS OF VARIABLES. C STRAINMETER AZIMUTH IS POSITIVE EAST FROM NORTH REAL LAMDA, LL, H,K,L,KR,MU,MUT,MUL,MUTT,MULL,MUTL,ELL, 1 ETT,ETL DOUBLE PRECISION SAW,CAW,RBARR DOUBLE PRECISION CZ,SZ,Z,X,THFPSI,CNM,SNM,CW,SW,CNU,SNU,NU DOUBLE PRECISION TIME,TT,PSIG,PSIM,HR,AW DOUBLE PRECISION PI20,HALFPI,TWOPI,AZT,HS,T2,PS,T,PSIS,HM,PM,NM, 1 ES,SIG,LS,LM,TANP DOUBLE PRECISION DATANZ INTEGER OPT DIMENSION H(6),K(6),L(6),KR(7),P(7),PP(7),PPP(7) H(1)=0.612 H(2)=0.290 H(3)=0.175 H(4)=0.129 H(5)=0.107 H(6)=0.095 K(1)=0.302 K(2)=0.093 K(3)=0.042 K(4)=0.025 K(5)=0.017 K(6)=0.013 L(1)=0.083 L(2)=0.014 L(3)=0.010 L(4)=0.008 L(5)=0.007 L(6)=0.005 JS=2 JE=3 IF(LORDER.EQ.0) GO TO5 JE=LORDER JS=JE 5 CONTINUE GTSUN=0. TTSUN=0. STNSUN=0. LAMDA=GLONG THETA=90.0-GLAT AZS=-ORIENT AZT=-ORIENT PI20=62.83185307169D0 HALFPI=1.570796326793D0 TWOPI=6.283185307169D0 RE=6.37122E+8 G=982.02 AZT=AZT*0.0174532925199D0 AZS=AZS*0.0174532925199D0 CAZT=COS(AZT) SAZT=SIN(AZT) CAZS=COS(AZS) SAZS=SIN(AZS) TT=TIME+604848.0D0 THETA=THETA*0.0174532925199D0 LAMDA=LAMDA*0.0174532925199D0 CTHETA=COS(THETA) STHETA=SIN(THETA) HR=DMOD(TIME,24.0D0) T=(TT+12.0D0)/876600.0D0 T2=T*T HS=4.881627982482D0+628.3319508731D0*T HS=HS+0.523598775578D-5*T2 HS=DMOD(DMOD(HS,PI20)+PI20,PI20) PS=4.908229466993D0+T*0.0300052641669D0 PS=PS+0.790246300201D-5*T2 ES=0.01675104D0- 0.00004180D0*T-0.000000126D0*T2 PSIG=0.2617993877971D0*(HR-12.0D0)+HS IF(LORDER.NE.0.AND.LORDER.NE.2) GO TO 82 LS=HS+2.0D0*ES*DSIN(HS-PS)+1.25D0*ES**2*DSIN(2.0D0*(HS-PS)) CZ=0.397980654630D0*DSIN(LS) SZ=DSQRT(1.0D0-CZ**2) Z=DATANZ(SZ,CZ) X=0.5D0*(LS+Z-HALFPI) THFPSI=1.52386101015D0*(DSIN(X)/DCOS(X)) PSIS=2.0D0*DATAN(THFPSI) PSIS=DMOD(PSIS+TWOPI,TWOPI) RBARR=1.D0+ES*DCOS(HS-PS)+ES**2*DCOS(2.D0*(HS-PS)) LL=PSIS-PSIG SLL=SIN(LAMDA-LL) CLL=COS(LAMDA-LL) CALFA=CTHETA*CZ+STHETA*SZ*CLL XI=4.2635E-5*RBARR CC=2.120E14*XI KR(2)=CC*XI*XI KR(3)=KR(2)*XI KR(4)=KR(3)*XI MU=CALFA P(2)=0.5*(3.*MU*MU-1.) P(3)=0.5*MU*(5.*MU*MU-3.) P(4)=0.25*(7.*MU*P(3)-3.*P(2)) PP(2)=3.*MU PP(3)=1.5*(5.*MU*MU-1.) PP(4)=0.25*(7.*(P(3)+MU*PP(3))-3.*PP(2)) PPP(2)=3. PPP(3)=15.*MU PPP(4)=7.5*(7.*MU*MU-1.) MUT=-STHETA*CZ+CTHETA*SZ*CLL MUTT=-MU MUL=-STHETA*SZ*SLL MULL=-STHETA*SZ*CLL MUTL=-CTHETA*SZ*SLL J=2 GO TO(50,60,70),OPT 50 GTSUN=0. DEL=1.+2.*H(J-1)/J-(J+1)*K(J-1)/J GTSUN=GTSUN+DEL*J*KR(J)*P(J)*G/RE GO TO 80 60 TTSUN=0. DIM=1.+K(J-1)-H(J-1) TTSUN=TTSUN+(DIM*KR(J)*PP(J)/RE)*(MUT*SAZT+MUL*CAZT/STHETA) GO TO 80 70 ELL=0. ETL=0. ETT=0. M=J-1 ELL=ELL+KR(J)*(L(M)*(PP(J)*MUT*CTHETA/STHETA+(PPP(J)* 1MUL*MUL+PP(J)*MULL)/(STHETA*STHETA))+H(M)*P(J)) ETT=ETT+KR(J)*(L(M)*(PPP(J)*MUT*MUT+PP(J)*MUTT)+ 1 H(M)*P(J)) ETL=ETL+2.*(KR(J)*L(M)*(PPP(J)*MUL*MUT+PP(J)*MUTL- 1 CTHETA*PP(J)*MUL/STHETA))/STHETA STNSUN=(ELL*SAZS*SAZS+ETL*SAZS*CAZS+ETT*CAZS*CAZS)/RE 80 CONTINUE 82 HM=4.720008893870D0+8399.709274530D0*T HM=HM+0.345575191895D-4*T2 HM=DMOD(DMOD(HM,PI20)+PI20,PI20) PM=5.835152597865D0+71.01804120839D0*T PM=PM-0.180108282526D-3*T2 NM=4.523601611611D0-33.75714629423D0*T NM=NM+0.362640633469D-4*T2 CNM=DCOS(NM) SNM=DSIN(NM) CW=0.91369D0-0.03569D0*CNM SW=DSQRT(1.0D0-CW**2) SNU=0.08968D0*SNM/SW CNU=DSQRT(1.0D0-SNU**2) NU=DATAN(SNU/CNU) SAW=0.39798D0*SNM/SW CAW=CNM*CNU+0.91739D0*SNM*SNU AW=2.0D0*DATAN(SAW/(1.0D0+CAW)) AW=DMOD(AW+TWOPI,TWOPI) SIG=HM-NM+AW LM=SIG+0.109801D0*DSIN(HM-PM)+0.003768D0*DSIN(2.0D0*(HM-PM)) 1 -0.224412D0*ES*DSIN(HS-PS)+0.02045D0*DSIN(HM-2.0D0*HS+PM) 2 +0.010809D0*DSIN(2.0D0*(HM-HS))+0.000653D0*DSIN(3.0D0*HM-2.0D0*HS 3 -PM) +0.000451D0*DSIN(2.0D0*HM-3.0D0*HS+PS) CZ=DSIN(LM)*SW WM=DATANZ(SW,CW) SZ=DSQRT(1.0D0-CZ**2) ZM=DATANZ(SZ,CZ) X=0.5D0*(LM+ZM-HALFPI) TANP=(DSIN(X)/DCOS(X))*DSIN(0.5D0*(HALFPI+WM))/DSIN(0.5D0*(HALFPI- 1 WM)) PSIM=2.0D0*DATAN(TANP) PSIM=DMOD(PSIM+TWOPI,TWOPI) RBARR=1.0D0+0.05467D0*DCOS(HM-PM)+0.003011D0*DCOS(2.0D0*(HM-PM)) 1 -0.000141D0*DCOS(HS-PS)+0.009271D0*DCOS(HM-2.0D0*HS+PM) 2+0.007759D0*DCOS(2.0D0*(HM-HS))+0.000633D0*DCOS(3.D0*HM-2.D0*HS-PM 3) +0.000328D0*DCOS(2.D0*HM-3.D0*HS+PS) LL=NU+PSIM-PSIG SLL=SIN(LAMDA-LL) CLL=COS(LAMDA-LL) CALFA=CTHETA*CZ+STHETA*SZ*CLL XI=1.65933E-2*RBARR CC=7.834E6*XI KR(2)=CC*XI*XI J=3 KR(J)=KR(J-1)*XI MU=CALFA P(1)=MU P(2)=0.5*(3.*MU*MU-1.) J=2 P(J+1)=((2*J+1)*MU*P(J)-J*P(J-1))/(J+1) PP(1)=1. PP(2)=3*MU PPP(1)=0. PPP(2)=3. J=2 PP(J+1)=((2*J+1)*(P(J)+MU*PP(J))-J*PP(J-1))/(J+1) PPP(J+1)=((2*J+1)*(2*PP(J)+MU*PPP(J))-J*PPP(J-1))/(J+1) MUT=-STHETA*CZ+CTHETA*SZ*CLL MUTT=-MU MUL=-STHETA*SZ*SLL MULL=-STHETA*SZ*CLL MUTL=-CTHETA*SZ*SLL GO TO (100,110,120),OPT 100 GTMUN=0. DO 105 J=JS,JE DEL=1.+2*H(J-1)/J-(J+1)*K(J-1)/J 105 GTMUN=GTMUN+DEL*J*KR(J)*P(J)*G/RE GO TO 130 110 TTMUN=0. DO 115 J=JS,JE DIM=1.+K(J-1)-H(J-1) 115 TTMUN=TTMUN+(DIM*KR(J)*PP(J)/RE)*(MUT*SAZT+MUL*CAZT/STHETA) GO TO 130 120 ELL=0. ETL=0. ETT=0. DO 125 J=JS,JE M=J-1 ELL=ELL+KR(J)*(L(M)*(PP(J)*MUT*CTHETA/STHETA+(PPP(J)* 1 MUL*MUL+PP(J)*MULL)/(STHETA*STHETA))+H(M)*P(J)) ETT=ETT+KR(J)*(L(M)*(PPP(J)*MUT*MUT+PP(J)*MUTT)+H(M)*P(J)) ETL=ETL+2*(KR(J)*L(M)*(PPP(J)*MUL*MUT+PP(J)*MUTL- 1 CTHETA*PP(J)*MUL/STHETA))/STHETA 125 STNMUN=(ELL*SAZS*SAZS+ETL*SAZS*CAZS+ETT*CAZS*CAZS)/RE 130 CONTINUE GO TO (140,150,160),OPT 140 RESULT=(GTSUN+GTMUN)*SCALE GO TO 170 150 RESULT=(TTSUN+TTMUN)*SCALE GO TO 170 160 RESULT=(STNSUN+STNMUN)*SCALE 170 CONTINUE 180 CONTINUE RETURN END c datanz returns values of phase in rads between plus and minus pi c has been tested in all quadrants c if both y and x are equal to zero,datanz returns pi/2 double precision function datanz(y,x) double precision x,y,halfpi halfpi= 2.00d0*datan(1.000d0) if(x.ne.0.0d0) go to 1 if(y.ge.0.00d0) datanz=halfpi if(y.lt.0.00d0) datanz=-halfpi return 1 datanz=datan(y/x) if(x.lt.0.0d0) datanz=datanz+2.00d0*halfpi if(y.lt.0.0d0.and.x.lt.0.0d0) datanz=datanz-4.00d0*halfpi return end