program SPLUN c c This program performs a cubic interpolation from c nk unequally spaced levels to ni levels. c The interpolation uses a cubic spline with tension. c The tension parameter (tausq) is set to 1.0 for this c example (data intgrity and curvatures are given the c same weight). c parameter (nk = 46,n = nk,nk1 = nk-1,ni = 82) parameter (n2 = nk-2 ,nest = 1000) real t(nk),sigh(nk),delx(nk1),xq(nk),q(nk),xw(nk) real xdum(nest),xout(ni),xint(ni),workint(ni),xt(nk) real da(n),xx(n),f1(nest),work(nest) c c Declare double precision variables c real*8 a(n),b(n),c(n),w(n),g(n),y(n),h(n),DERIVA real*8 ap(n2),bp(n2),cp(n2),up(n2),wp(n2),fp(n2),yp(n2) c c Open output files c open(4,file ='temp.out ',status='unknown') open(7,file ='spechum.out',status='unknown') c c Set the tension coefficient. this coefficient varies between c 0 and 1 inclusive. Values close to 1 give the data quality c high weight whereas values close to 0 give more weight to the c curvature . c tausq=1.0 c c Sample data used for this example. c c pressure levels. c data sigh/1011.,1000.,955.,950.,900.,850.,816.,800.,782.,750., & 722.,700.,687.,661.,655.,650.,641.,600.,594.,578.,566.,550., & 545.,535.,528.,518.,509.,500.,493.,478.,469.,450.,444.,400., & 390.,383.,370.,362.,350.,311.,300.,258.,250.,236.,223.,200./ c c temperature c data t/296.2,295.9,295.2,295.0,291.6,288.8,286.4,285.6,284.5, & 283.1,282.0,280.6,279.2,279.1,278.9,278.5,277.8,274.3,273.8, & 272.6,272.0,271.3,270.7,269.3,268.5,268.5,267.5,267.2,266.8, & 264.7,263.6,261.8,261.1,255.3,254.7,253.7,251.6,251.2,249.2, & 242.2,240.4,231.8,230.1,226.8,223.3,216.9/ c c specific humidity. c data q/.1994,.1914,.1812,.1816,.1730,.1472,.1465,.1424,.1366, & .1313,.1278,.1214,.1070,.0901,.0884,.0863,.0831,.0774,.0762, & .0819,.0735,.0472,.0421,.0445,.0606,.0539,.0435,.0168,.0190, & .0281,.0600,.0398,.0349,.0323,.0164,.0134,.0174,.0123,.0101, & .0049,.0078,.0037,.0026,.0013,.0011,.0009/ c c invert levels ( This step is not necessary). c do 6400 k = 1, nk kk = nk-k+1 xt(k) = t(kk) xq(k) = q(kk) 6400 xw(k) = sigh(kk) write(6,1001) write(6,1002) do 6402 k = 1, nk sigh(k)= xw(k) t(k) = xt(k) q(k) = xq(k) write(6,1003)k,sigh(k),t(k),q(k) 1001 format(15x,'original ',//) 1002 format(8x,'level',2x,'pres ',5x,'temp ',3x,'sp hu',//) 1003 format(6x,i4,2f10.2,f7.4) 6402 continue c c Define a work array for interpolation . c Variable ivar defines the variable to c interpolate. c do 6403 ivar = 1, 2 do 6404 k = 1,nk if (ivar.eq.1) then da(k) = t(k) else da(k) = q(k) endif 6404 continue c c Define the increment array of the independent c variable . c do 6406 k = 1, nk1 delx(k)= sigh(k+1) - sigh(k) 6406 continue c c Compute the values of the multiplyers by solving the c matrix equation.Define the diagonal vectors. c do 6408 i = 3, n2 ap(i) = 36.*tausq*(1./(delx(i)*delx(i-1))) 6408 continue do 6410 i = 2, n2 bp(i) = 2.*delx(i)-36.*tausq*(1./(delx(i-1)*delx(i)) & +1./(delx(i)*delx(i+1))+2./(delx(i)*delx(i))) 6410 continue do 6412 i = 1,n2 cp(i) = 4.*(delx(i)+delx(i+1))+72.*tausq*(1./(delx(i)* & delx(i+1))+1./(delx(i+1)*delx(i+1))+1./(delx(i)*delx(i))) 6412 continue c c Define the forcing vector. c do 6414 i = 1, n2 x1 = (da(i+1) - da(i))/delx(i) x2 = (da(i+2) - da(i+1))/delx(i+1) c x1 = DERIVA (da,delx,i) c x2 = DERIVA (da,delx,i+1) fp(i) = -(x2-x1)*6.0*tausq 6414 continue c call FIVEDI (ap,bp,cp,up,wp,fp,yp,n2) c c Compute the solution ordinates c h(1) = da(1)+6.*yp(1)/delx(1) h(2) = da(2)-(6./delx(1)+6./delx(2)) & *yp(1) + 6.*yp(2)/delx(2) do 6416 i = 3, n2 h(i) = da(i) + 6.*yp(i-2)/delx(i-1)-(6./delx(i-1) & + 6./delx(i))*yp(i-1) + 6.*yp(i)/delx(i) 6416 continue h(n-1) = da(n-1) + 6.*yp(n2-1)/delx(n-2) & -(6./delx(n2) + 6./delx(n-1))*yp(n2) h(n) = da(n) + 6.*yp(n2)/delx(n-1) c c compute the values of the spline coefficients by solving the c tridiagonal matrix equation. Define the diagonal vectors. c do 6418 i = 2,n-1 a(i) = delx(i-1) b(i) = 2.*(delx(i-1)+delx(i)) c(i) = delx(i) 6418 continue b(1) = 2.0*delx(1) b(n) = 2.0*delx(n-1) c(1) = delx(1) a(n) = delx(n-1) c c Compute the forcing vector. c g(1) = yp(1)*delx(1) g(2) = 2.0*yp(1)*(delx(1)+delx(2)) + yp(2)*delx(2) do 6420 i = 3, n2 g(i) = yp(i-2)*delx(i-1)+2.0*yp(i-1)* & (delx(i-1)+delx(i))+yp(i)*delx(i) 6420 continue g(n-1) = (yp(n2-1)*delx(n2))+(2.0*yp(n2)* & (delx(n2)+delx(n-1))) g(n) = yp(n2)*delx(n-1) do 6422 i = 1, n g(i) = -2.0*g(i)/tausq 6422 continue c c Call subroutine tridi to solve the tridiagonal matrix. c call TRIDI (a,b,c,g,y,w,n) c c Define the cubic spline polynomial and generate the c values of dependent variable.Define the sampling locations. c del = 1. icount = 1 do 6424 j = 1, n xx(j) = sigh(j) 6424 continue f1(1) = h(1) work(icount) = f1(1) xdum(icount) = xx(1) do 6426 j = 2, n alower = xx(j-1) upper = xx(j) npts = (upper-alower)/del +1 do 6428 i = 2,npts+1 x = alower + del*(i-1) icount = icount + 1 h1 = (y(j-1)/6.)*((xx(j)-x)**3/delx & (j-1)-delx(j-1)*(xx(j)-x)) h2 = (y(j)/6.)*((x-xx(j-1))**3/delx & (j-1)-delx(j-1)*(x-xx(j-1))) h3 = h(j-1)*(xx(j)-x)/delx(j-1) h4 = h(j)*(x-xx(j-1))/delx(j-1) f1(i) = h1+h2+h3+h4 c c Save output by interval . c xdum(icount) = x work(icount) = f1(i) 6428 continue 6426 continue c c Select the values of the function at the desired points. c delc = 10. do 6430 k = 1, ni c c Define the array at which function is to be interpolated. c xout(k)= xx(1) + (k-1)*delc do 6432 j = 1,nest if (xout(k).eq.xdum(j)) then xint(k) = xout(k) workint(k)= work(j) endif 6432 continue 6430 continue c c Write output to units 4(temperature) and 7(spec. hum). c do 6434 k = 1, ni if (ivar.eq.1) then write (4,1004)xint(k),workint(k) else write (7,1005)xint(k),workint(k) endif 6434 continue 1004 format(2f10.2) 1005 format(2f10.5) c c Display output on the screen. c write(6,1006) if (ivar.eq.1) then write(6,1007) else write(6,1010) endif write(6,1009)1,xdum(1),workint(1),42,xint(42),workint(42) do 6436 k = 2, ni/2 6436 write(6,1009)k,xint(k),workint(k),k+41,xint(k+41),workint(k+41) 1006 format(/,15x,'interpolated ',//) 1007 format(2x,2('level',2x,'pres ',8x,'temp ',3x),//) 1008 format(6x,i4,2f10.2) 1009 format(1x,i2,2x,f10.5,2x,f10.5,2x,i2,2x,f10.5,2x,f10.5) 1010 format(2x,2('level',2x,'pres ',8x,'sp hu',3x),//) c 6403 continue stop end subroutine TRIDI (a,b,c,g,y,work,n) c c This subroutine solves a tridiagonal c matrix system. c a,b,c : are the diagonal vectors of the c tridiagonal coefficient matrix. c g : is the known vector. c y : is the solution vector. c n : is the dimension of the solution. c real*8 a(1),b(1),c(1),g(1),y(1),work(1) if (b(1).eq.0) then print *, ' b(1) = 0)' stop endif denom = b(1) y(1) = g(1)/denom do 6410 i = 2, n work(i-1) = c(i-1)/denom denom = b(i)-a(i)*work(i-1) if (denom.eq.0) then print *, 'denominator vanishes' stop endif 6410 y(i) = (g(i)-a(i)*y(i-1))/denom c c Initiate downward recursion c do 6411 i = n-1, 1,-1 6411 y(i) = y(i)-work(i)*y(i+1) return end function DERIVA (y,x,i) real*8 deriv real x(1),y(1) deriv = (y(i+1)-y(i))/x(i) return end subroutine FIVEDI (a,b,c,u,w,f,y,n) c c Subroutine FIVEDI solve a pentadiagonal matrix system. c a,b,c : are the (repeated) diagonal vectors c of the pentadiagonal matrix. c f : is the known vector. c u,w : are work vectors. c y : is the solution vector. c n : is the solution dimension c real*8 a(1),b(1),c(1),u(1),w(1),f(1),y(1),denom c c Check for singularity c if (c(1).eq.0.0) then print *,'c(1) = 0,singular matrix' stop endif c c setup sequences;initiate upward recursions c u(1) = b(2)/c(1) w(1) = a(3)/c(1) y(1) = f(1)/c(1) denom = c(2)-b(2)*u(1) c c check for singularity c if (denom.eq.0.0) then print *,'denom = 0,singular at step 2' stop endif u(2) = ( b(3)-b(2)*w(1))/denom w(2) = a(4)/denom y(2) = ( f(2)-b(2)*y(1))/denom do 6320 j = 3,n denom = c(j)-b(j)*u(j-1)+a(j)*(u(j-2)*u(j-1)-w(j-2)) c c check for singularity c if (denom.eq.0.0) then print *, 'denom = 0 in loop ' stop endif c c check sequence upper limit c if (j.le.n-2) then w(j) = a(j+2)/denom end if if (j.le.n-1) then u(j) = (b(j+1)-w(j-1)*(b(j)-a(j)*u(j-2)))/denom end if 6320 y(j) = (f(j)-b(j)*y(j-1)+a(j)*(u(j-2) & *y(j-1)-y(j-2)))/denom c c Initiate donward recursion for solution c y(n-1) = y(n-1)-u(n-1)*y(n) do 6321 j = n-2,1,-1 6321 y(j) = y(j)-u(j)*y(j+1)-w(j)*y(j+2) return end