program lod c A program to evaluate the length of day and its time derivative c after removal of atmospheric signal and running averaging. c Code derived from ufm modelling code c Questions / comments --> holme@liv.ac.uk (Richard Holme) c--------------------------------------------------------------------- c c uses code by: Jeremy Bloxham, Andrew Jackson, and Carl de Boor c c--------------------------------------------------------------------- c for details of B-splines see: c Carl de Boor "A Practical Guide to Splines" c Springer-Verlag, 1978. c--------------------------------------------------------------------- c c nspl is number of B-splines c c tknts is an array of knot points c c--------------------------------------------------------------------- c c jord = 4 (order of B-splines) c c--------------------------------------------------------------------- c c CALLS: getlod - the driving subroutine to interogate the c spline coefficients c c interv - calculates which knot lies immediately left c of the current time point c c bspline - calculates B-splines at current time point c c bspline1 - calculates first time derivative of c B-splines at current time point c c b0, b1 - functions required by bspline1 c c---------------------------------------------------------------------- implicit undefined (a-z) c note that the parameters in this subroutine are different from c the rest of the code, because this bit is self-contained integer nspl parameter (nspl=100) integer irough,ns,k real*8 timein,dlodt,ddlodt,tstart,tend real*8 dlod(nspl),tknts(nspl+4) character*10 splinefile c choose the rough or smooth fit write(6,*) 'Do you want the rough or smooth fit?' write(6,*) 'Enter 1 for rough.' read(5,*) irough if (irough .eq. 1) then splinefile = 'spl_rough' else splinefile = 'spl_smooth' endif open(1,file=splinefile) write(6,*) 'Using the spline file ', splinefile write(6,*) c----- c input current model read(1,*) tstart,tend read(1,*) ns,(tknts(k),k=1,ns+4) read(1,*) (dlod(k),k=1,ns) close(1) write(6,*) 'Do you want a single estimate, or a time series?' write(6,*) 'Enter 1 for a single estimate.' read(5,*) irough write(6,*) if (irough .eq. 1) then write(6,*) 'Enter the calculation time:' write(6,*) 'Allowed values between',tstart,tend read(5,*) timein write(6,*) c if timein is out of the range of splines, we cant do this if (timein .lt. tstart .or. timein .gt. tend) then write(6,*) 'Cant calculate for time',timein stop endif call getlod(timein,dlodt,ddlodt,nspl,ns,tknts,dlod) write(6,"(a24,f9.3)") 'Results for time: ',timein write(6,"(a24,f9.3)") 'Excess LOD (ms): ',dlodt write(6,"(a24,f9.3)") 'Rate of change (ms/yr): ',ddlodt else open(unit=1,file='lodout') do k=0,10000 timein = tstart + k/12.0 if (timein .ge. tend) goto 1000 call getlod(timein,dlodt,ddlodt,nspl,ns,tknts,dlod) write(1,"(f10.3,x,f8.4,x,f8.4)") timein,dlodt,ddlodt end do 1000 close(1) write(6,*) 'Output writen to lodout.' write(6,*) 'Format: time (years), dlod, rate of change' endif end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine getlod(timein,dlodt,ddlodt,ns,nspl,tknts,dlod) c calculate excess lod and rate of change at timein implicit none integer jord,ns,k,nleft,nspl real*8 dlodt,ddlodt,timein,spl real*8 spl(nspl),tknts(nspl+4),dlod(nspl) data jord/4/ c----- c calculate length of day excess at time timein c do k=1,nspl spl(k)=0. enddo call interv(tknts,timein,ns,nleft) call bspline(tknts,timein,nspl,jord,nleft,spl(nleft-3)) dlodt = 0.0 do k=1,4 dlodt = dlodt + spl(k+nleft-4)*dlod(k+nleft-4) enddo c----- c calculate rate of change of LOD per year c do k=1,nspl spl(k)=0. enddo call bspline1(tknts,timein,nspl,nleft,spl) ddlodt=0.0 do k=1,4 ddlodt = ddlodt + spl(k+nleft-4)*dlod(k+nleft-4) enddo return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine interv(tknts,time,nspl,nleft) implicit real*8 (a-h,o-z) dimension tknts(nspl+1) c interval starts at tknts(4) and ends at tknts(nspl+1) c----- c calculate nleft: c tknts(nleft) < tknts(nleft+1) c tknts(nleft) <= time <= tknts(nleft+1) c c check we are in-range if(time.lt.tknts(4).or.time.gt.tknts(nspl+1)) return do 200 n=5,nspl+1 if(time.le.tknts(n)) then nleft=n-1 goto 210 endif 200 continue 210 continue return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine bspline(tknts,t,nspl,jorder,nleft,spl) implicit real*8(a-h,o-z) c calculate splines of order jorder where 1 <= jorder <= 4 dimension tknts(nspl+4) dimension spl(4) dimension deltal(4),deltar(4) spl(1)=1.0 do 200 j=1,jorder-1 deltar(j) = tknts(nleft+j) - t deltal(j) = t - tknts(nleft+1-j) saved=0.0 do 100 i=1,j term = spl(i)/(deltar(i)+deltal(j+1-i)) spl(i) = saved + deltar(i)*term saved = deltal(j+1-i)*term 100 continue spl(j+1) = saved 200 continue return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine bspline1(tknts,t,nspl,nl,spl) implicit real*8(a-h,o-z) dimension tknts(nspl+4) dimension spl(nspl),spl1(500) data jord/3/ do 100 is=1,nspl spl(is)=0.0 100 continue call bspline(tknts,t,nspl,jord,nl,spl1(nl-2)) spl(nl) = b0(tknts,nl)*spl1(nl) spl(nl-1) = b0(tknts,nl-1)*spl1(nl-1) + b1(tknts,nl-1)*spl1(nl) spl(nl-2) = b0(tknts,nl-2)*spl1(nl-2) + b1(tknts,nl-2)*spl1(nl-1) spl(nl-3) = b1(tknts,nl-3)*spl1(nl-2) return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc function b0(tknts,i) implicit real*8(a-h,o-z) dimension tknts(500) b0 = 3.0/( tknts(i+3) - tknts(i) ) return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc function b1(tknts,i) implicit real*8(a-h,o-z) dimension tknts(500) b1 = -3.0/( tknts(i+4) - tknts(i+1) ) return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc