Commit ddf18bee authored by thomas.forbriger's avatar thomas.forbriger

libseife [TASK]: handle case of nstep=1 in tide removal explicitly

parent e73b9e98
c this is <seife_tides.f>
c------------------------------------------------------------------------------
c ($Id$)
c
c Copyright 1984 by Erhard Wielandt
c This code was part of seife.f. A current version of seife.f can be obtained
......@@ -26,6 +25,7 @@ c extracted from libseife.f
c
c REVISIONS and CHANGES
c 25/10/2000 V1.0 Thomas Forbriger
c 22/06/2020 V1.1 do computation for nstep=1 explicitly
c
c==============================================================================
c
......@@ -50,6 +50,11 @@ c to the total length of the record.
step=nstep
tstep2=(step-one)*dth
tint=tstep2+dth
if (nstep.eq.1) then
step=1.d0
tstep2=0.d0
tint=dth
endif
c determine the number of frequencies required for a good fit
dur=n*dt/3600.d0
nfreq=6
......@@ -66,6 +71,55 @@ c determine partial amplitudes by least-squares fit
rs(i)=0.d0
do 1 k=1,nco
1 a(i,k)=0.d0
if (nstep.eq.1) then
c ======================================================================
c do not average
c ==============
c correction for averaging over nstep samples
do j=1,nfreq
cor(j)=1.d0
enddo
c set up system of linear equations
c(1)=one
do j=1,n
sx=zero
do jj=j,j
sx=sx+x(jj)
enddo
t=(j-1)*dt
do k=1,nfreq
c(2*k) =dcos(omeg(k)*t)
c(2*k+1)=dsin(omeg(k)*t)
enddo
do i=1,nco
rs(i)=rs(i)+sx*c(i)
do k=1,nco
a(i,k)=a(i,k)+c(i)*c(k)
enddo
enddo
enddo
c solve for partial amplitudes
call seife_gauss(a,nco,ndim,rs,f)
do j=1,n
t=(j-1)*dt
c remove average and tides
xgez1=f(1)
do k=1,nfreq
k2=2*k
k3=k2+1
c(k2)=dcos(omeg(k)*t)
c(k3)=dsin(omeg(k)*t)
xgez1=xgez1+f(k2)*c(k2)+f(k3)*c(k3)
enddo
x(j)=x(j)-xgez1
enddo
write(6, *) f(1),(omega(j),f(2*j),f(2*j+1), j=1,nfreq)
write(msg, '("remove tides with nfreq: ",i10)')nfreq
nil=.false.
else
c ======================================================================
c average over nstep samples
c ==========================
c correction for averaging over nstep samples
do 102 j=1,nfreq
102 cor(j)=step*dsin(omeg(j)*dth)/dsin(step*omeg(j)*dth)
......@@ -111,9 +165,10 @@ c remove average and tides
do 3 jj=j,j+nstep-1
tdif=(jj-j)*dt-tstep2
3 x(jj)=x(jj)-xgez1-xgez2*tdif-xgez3*tdif*tdif
write(6, *) f(1),(omega(j),f(2*j),f(2*j+1), j=1,nfreq)
write(6, *) f(1),(omega(j),f(2*j),f(2*j+1), j=1,nfreq)
write(msg, '("remove tides with nfreq: ",i10)')nfreq
nil=.false.
endif
return
end
c
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment