C C This is the subroutine of the three-points sixth-order combined compact C difference (CCD) scheme. It is used for reducing numerical errors C especially the errors in horizontal pressure gradient computation. C C C The references for the high-order CCD schemes are listed as follows: C C Chu, P.C., and C.W. Fan, 1997: Sixth-order difference scheme for sigma C coordinate ocean models. 27, 2064-2071. C C Chu, P.C., and C.W. Fan, 1998a: A three-point combined compact difference C scheme. Journal of Computational Physics, 140, 370-399. C C Chu, P.C., and C.W. Fan, 1998b: Improvement of estuarine and coastal C modeling using high-order difference schemes. Estuarine and Coastal C Modeling, 5, American Society of Engineering, 28-34. C C Chu, P.C., and C.W. Fan, 1999: A non-uniform three-point combined C compact difference scheme. Journal of Computational Physics, 148, C 663-674. C C Chu, P.C., and C.W. Fan, 2000: A staggered three-point combined C compact difference scheme. Mathematical and Computer Modeling, in press. C C C If you have any question or need reprints, please contact C Peter Chu (chu@nps.navy.mil) or Chenwu Fan (fan@nps.navy.mil) at C the Naval Postgraduate School. C C----------------------------------------------------------------------------- C C Additional notes to POM users: C C 1. Simply replace BAROPG; only other required change is to add C parameter ijm = max(im,jm) to comblk97.h C 2. On 64bit machines it is adviced to comment out "double precision" C in subroutines getd2f & getd2fm. C 3. The code does not vectorized well on a vector machine like the Cray T90 C and thus may require more than 10 times the cpu of the standard pom97.f code. C On an SGI machine additional computational cost is only around 60%. C T.E C----------------------------------------------------------------------------- SUBROUTINE BAROPG(DRHOX,DRHOY) INCLUDE 'comblk97.h' c DIMENSION DRHOX(IM,JM,KB),DRHOY(IM,JM,KB), & drhos(kbm1,im,jm),aa(3,ijm),r1d(im+1,jm,kbm1), & ttt(ijm),t1d(ijm+1),t2d(ijm),tint(ijm), & xy(ijm),rh1d(ijm),f1m(jm), & rh(im,jm,kbm1),dh(im,jm),d1d(im+1,jm) C --- included in comblk --- C grav=9.8 do i=1,ijm xy(i)=i enddo DO 299 K=1,KB DO 299 J=1,JM DO 299 I=1,IM 299 RHO(I,J,K)=RHO(I,J,K)-RMEAN(I,J,K) C C X COMPONENT OF BAROCLINIC PRESSURE GRADIENT C c c inteperate rho to u location: c do j=2,jmm1 xl=1.0 xr=im do k=1,kbm1 do i=1,im ttt(i)=rho(i,j,k) enddo call getd2fm(im,ttt,t2d,xy,xl,xr,0.0,0.0,fsm(1,j)) call apfdcgnprdx(im+1,aa,ttt,t1d,t2d,fsm(1,j)) call apfdcgnprdx2(im,aa,ttt,t2d,t1d,fsm(1,j)) call apfdcgnprdx(im+1,aa,ttt,r1d(1,j,k),t2d,fsm(1,j)) call apfdcgnprdx2(im,aa,ttt,t2d,r1d(1,j,k),fsm(1,j)) call fmidn(im,ttt,r1d(1,j,k),t2d,rh(1,j,k),fsm(1,j)) enddo do i=1,im ttt(i)=DT(i,j) enddo call getd2fm(im,ttt,t2d,xy,xl,xr,0.0,0.0,fsm(1,j)) call apfdcgnprdx(im+1,aa,ttt,t1d,t2d,fsm(1,j)) call apfdcgnprdx2(im,aa,ttt,t2d,t1d,fsm(1,j)) call apfdcgnprdx(im+1,aa,ttt,d1d(1,j),t2d,fsm(1,j)) call apfdcgnprdx2(im,aa,ttt,t2d,d1d(1,j),fsm(1,j)) call fmidn(im,ttt,d1d(1,j),t2d,dh(1,j),fsm(1,j)) enddo c c using 4th compact differencing calculate dr/dsig c do i=2,imm1 do j=2,jmm1 do k=1,kbm1 ttt(k)=rh(i,j,k) enddo call getd2f(kbm1,ttt,t2d,zz,0.0,-1.0,0.0,0.0) call getdf(kbm1,ttt,drhos(1,i,j),t2d,zz) enddo enddo do j=2,jmm1 do i=2,imm1 do k=1,kbm1 ttt(k)=grav*(r1d(i,j,k)*dh(i,j) & -zz(k)*d1d(i,j)*drhos(k,i,j)) enddo call getd2f(kbm1,ttt,t2d,zz,0.0,-1.0,0.0,0.0) call getdf(kbm1,ttt,t1d,t2d,zz) call getint(kbm1,zz,ttt,t1d,t2d,zz,tint,kbm1) do k=1,kbm1 drhox(i,j,k)=-dh(i,j)*tint(k)*DUM(i,j)* & (dy(i,j)+dy(i-1,j))/2.0*RAMP enddo enddo enddo C C X COMPONENT OF BAROCLINIC PRESSURE GRADIENT C c DO 300 J=2,JMM1 c DO 300 I=2,IMM1 c 300 DRHOX(I,J,1)=.5E0*GRAV*(-ZZ(1))*(DT(I,J)+DT(I-1,J)) c 2 *(RHO(I,J,1)-RHO(I-1,J,1)) c DO 310 K=2,KBM1 c DO 310 J=2,JMM1 c DO 310 I=2,IMM1 c 310 DRHOX(I,J,K)=DRHOX(I,J,K-1) c 1 +GRAV*.25E0*(ZZ(K-1)-ZZ(K))*(DT(I,J)+DT(I-1,J)) c 2 *(RHO(I,J,K)-RHO(I-1,J,K)+RHO(I,J,K-1)-RHO(I-1,J,K-1)) c 3 +GRAV*.25E0*(ZZ(K-1)+ZZ(K))*(DT(I,J)-DT(I-1,J)) c 4 *(RHO(I,J,K)+RHO(I-1,J,K)-RHO(I,J,K-1)-RHO(I-1,J,K-1)) C c DO 360 K=1,KBM1 c DO 360 J=2,JMM1 c DO 360 I=2,IMM1 c 360 DRHOX(I,J,K)=.25E0*(DT(I,J)+DT(I-1,J))*DRHOX(I,J,K)*DUM(I,J) c 1 *(DY(I,J)+DY(I-1,J)) c c c Y COMPONENT OF BAROCLINIC PRESSURE GRADIENT C c c inteperate rho to v location: c do i=2,imm1 xl=1.0 xr=jm do j=1,jm f1m(j)=fsm(i,j) enddo do k=1,kbm1 do j=1,jm ttt(j)=rho(i,j,k) enddo call getd2fm(jm,ttt,t2d,xy,xl,xr,0.0,0.0,f1m) call apfdcgnprdx(jm+1,aa,ttt,t1d,t2d,f1m) call apfdcgnprdx2(jm,aa,ttt,t2d,t1d,f1m) call apfdcgnprdx(jm+1,aa,ttt,t1d,t2d,f1m) call apfdcgnprdx2(jm,aa,ttt,t2d,t1d,f1m) call fmidn(jm,ttt,t1d,t2d,rh1d,f1m) do j=1,jm rh(i,j,k)=rh1d(j) r1d(i,j,k)=t1d(j) enddo enddo do j=1,jm ttt(j)=DT(i,j) enddo call getd2fm(jm,ttt,t2d,xy,xl,xr,0.0,0.0,f1m) call apfdcgnprdx(jm+1,aa,ttt,t1d,t2d,f1m) call apfdcgnprdx2(jm,aa,ttt,t2d,t1d,f1m) call apfdcgnprdx(jm+1,aa,ttt,t1d,t2d,f1m) call apfdcgnprdx2(jm,aa,ttt,t2d,t1d,f1m) call fmidn(jm,ttt,t1d,t2d,rh1d,f1m) do j=1,jm dh(i,j)=rh1d(j) d1d(i,j)=t1d(j) enddo enddo c c using 4th compact differencing calculate dr/dsig c do i=2,imm1 do j=2,jmm1 do k=1,kbm1 ttt(k)=rh(i,j,k) enddo call getd2f(kbm1,ttt,t2d,zz,0.0,-1.0,0.0,0.0) call getdf(kbm1,ttt,drhos(1,i,j),t2d,zz) enddo enddo do j=2,jmm1 do i=2,imm1 do k=1,kbm1 ttt(k)=grav*(r1d(i,j,k)*dh(i,j) & -zz(k)*d1d(i,j)*drhos(k,i,j)) enddo call getd2f(kbm1,ttt,t2d,zz,0.0,-1.0,0.0,0.0) call getdf(kbm1,ttt,t1d,t2d,zz) call getint(kbm1,zz,ttt,t1d,t2d,zz,tint,kbm1) do k=1,kbm1 drhoy(i,j,k)=-dh(i,j)*tint(k)*DVM(i,j)* & (dx(i,j)+dx(i,j-1))/2.0*RAMP enddo enddo enddo c DO 500 J=2,JMM1 c DO 500 I=2,IMM1 c 500 DRHOY(I,J,1)=.5E0*GRAV*(-ZZ(1))*(DT(I,J)+DT(I,J-1)) c 1 *(RHO(I,J,1)-RHO(I,J-1,1)) c DO 510 K=2,KBM1 c DO 510 J=2,JMM1 c DO 510 I=2,IMM1 c 510 DRHOY(I,J,K)=DRHOY(I,J,K-1) c 1 +GRAV*.25E0*(ZZ(K-1)-ZZ(K))*(DT(I,J)+DT(I,J-1)) c 2 *(RHO(I,J,K)-RHO(I,J-1,K)+RHO(I,J,K-1)-RHO(I,J-1,K-1)) c 3 +GRAV*.25E0*(ZZ(K-1)+ZZ(K))*(DT(I,J)-DT(I,J-1)) c 4 *(RHO(I,J,K)+RHO(I,J-1,K)-RHO(I,J,K-1)-RHO(I,J-1,K-1)) C c DO 560 K=1,KBM1 c DO 560 J=2,JMM1 c DO 560 I=2,IMM1 c 560 DRHOY(I,J,K)=.25E0*(DT(I,J)+DT(I,J-1))*DRHOY(I,J,K)*DVM(I,J) c 1 *(DX(I,J)+DX(I,J-1)) c DO 561 K=1,KB c DO 561 J=2,JMM1 c DO 561 I=2,IMM1 c DRHOX(I,J,K)=RAMP*DRHOX(I,J,K) c 561 DRHOY(I,J,K)=RAMP*DRHOY(I,J,K) C DO 571 K=1,KB DO 571 J=1,JM DO 571 I=1,IM 571 RHO(I,J,K)=RHO(I,J,K)+RMEAN(I,J,K) C RETURN END C c The follow software is provided by Chenwu Fan and "as is" without c warranty of any kind. c c Chenwu Fan (fan@nps.navy.mil) 8/20/99 subroutine triapd(m,n,a,b) dimension a(3,m) real b(n,m) integer i,j,m,n do j=2,m w=a(1,j)/a(2,j-1) a(2,j)=a(2,j)-w*a(3,j-1) do i=1,n b(i,j)=b(i,j)-w*b(i,j-1) enddo enddo do i=1,n b(i,m)=b(i,m)/a(2,m) do j=m-1,1,-1 b(i,j)=(b(i,j)-a(3,j)*b(i,j+1))/a(2,j) enddo enddo return end c c subroutine fmidn is used to calculate the middle way value. c c 297 137 c f(0)= ---(f(1/2)+f(-1/2)) - ---(f(3/2)+f(-3/2)) + c 320 320 c c 324 207 c --- (f'(1)-f'(-1)) - --- (f''(1/2)+f''(-1/2)) c 320 320 c c 171 1 c + --- --- f''''''''(0)*h^8 c 32 8! c subroutine fmidn (l2,f,df,d2f,fh,fsm) integer l1,l2,i dimension f(l2),fh(l2),df(l2+1),d2f(l2),fsm(l2) dimension b(4),bb(6) integer nb(4,2),nbb(6,2) c data nb/ 297, -137, 324, -207, & 320, 320, 320, 320/ data nbb/ 13, 3, 127, 1, 5, -1, & 16, 16, 432, 54, 144, 32/ do i=1,4 b(i)=real(nb(i,1))/real(nb(i,2)) enddo do i=1,6 bb(i)=real(nbb(i,1))/real(nbb(i,2)) enddo if(fsm(1)+fsm(2) .lt. 0.5) then fh(2)=(f(1)+f(2))/2.0 else fh(2)=f(1)*bb(1)+f(2)*bb(2)+df(1)*bb(3) & +df(3)*bb(4)+d2f(1)*bb(5)+d2f(2)*bb(6) endif do i=3,l2-1 if(fsm(i-1)+fsm(i) .lt. 0.5) then fh(i)=(f(i-1)+f(i))/2.0 else fh(i)=b(1)*(f(i-1)+f(i))+b(2)*(f(i+1)+f(i-2))+ & b(3)*(df(i+1)-df(i-1))+b(4)*(d2f(i-1)+d2f(i)) endif enddo if(fsm(l2-1)+fsm(l2) .lt. 0.5) then fh(l2)=(f(l2-1)+f(l2))/2.0 else fh(l2)=f(l2)*bb(1)+f(l2-1)*bb(2)-df(l2+1)*bb(3) & -df(l2-1)*bb(4)+d2f(l2)*bb(5)+d2f(l2-1)*bb(6) endif return end c subroutine triprd(m,n,a,b) dimension a(3,m) real b(n,m) integer i,j,m,n do j=2,m-2 w=a(1,j)/a(2,j-1) a(1,j)=-w*a(1,j-1) a(2,j)=a(2,j)-w*a(3,j-1) do i=1,n b(i,j)=b(i,j)-w*b(i,j-1) enddo w=a(3,m)/a(2,j-1) a(3,m)=-w*a(3,j-1) a(2,m)=a(2,m)-w*a(1,j-1) do i=1,n b(i,m)=b(i,m)-w*b(i,j-1) enddo enddo j=m-2 w=a(1,m-1)/a(2,j) a(2,m-1)=a(2,m-1)-w*a(3,j) a(3,m-1)=a(3,m-1)-w*a(1,j) do i=1,n b(i,m-1)=b(i,m-1)-w*b(i,j) enddo w=a(3,m)/a(2,j) a(1,m)=a(1,m)-w*a(3,j) a(2,m)=a(2,m)-w*a(1,j) do i=1,n b(i,m)=b(i,m)-w*b(i,j) enddo j=m-1 w=a(1,m)/a(2,j) a(2,m)=a(2,m)-w*a(3,j) do i=1,n b(i,m)=b(i,m)-w*b(i,j) b(i,m)=b(i,m)/a(2,m) b(i,m-1)=(b(i,m-1)-a(3,m-1)*b(i,m))/a(2,m-1) do j=m-2,1,-1 b(i,j)=(b(i,j)-a(1,j)*b(i,m)-a(3,j)*b(i,j+1))/a(2,j) enddo enddo return end subroutine apfdcgnprdx(n,a,f,df,d2f,fsm) real a(3,n) dimension f(n-1),df(n),d2f(n-1),b(3), & aa(3,20),fsm(*) real ff(3),ddf(2) integer b1(3),b2(3),bb1(6),bb2(6) data b1/-7,127,-17/ data b2/240,120,240/ data one/1.0d0/ do i=1,3 b(i)=real(b1(i))/real(b2(i)) enddo df(1)=f(2)-f(1) df(n)=f(n-1)-f(n-2) do i=2,n-1 if(fsm(i)+fsm(i-1) .lt. 0.9) then a(1,i)=0.0 a(2,i)=1.0 a(3,i)=0.0 df(i)=f(i)-f(i-1) else a(1,i)=b(1) a(2,i)=b(2) a(3,i)=b(1) df(i)=(f(i)-f(i-1))+b(3)*(d2f(i)-d2f(i-1)) endif enddo call triapd(n-2,1,a(1,2),df(2)) return end subroutine apfdcgnprdx2(n,a,f,d2f,d1f5,fsm) real a(3,n) dimension f(n),d2f(n),d1f5(n+1),b(4),fsm(*) integer b1(4),b2(4) data b1/-5,47,-17,24/ data b2/84,42,7,7/ data one,two/1.0d0,2.0d0/ do i=1,4 b(i)=real(b1(i))/real(b2(i)) enddo a(1,1)=0.0 a(2,1)=1.0 a(3,1)=0.0 if(fsm(1)+fsm(2) .lt. 0.9) d2f(1)=0.0 do i=2,n-1 if(fsm(i-1)+fsm(i)+fsm(i+1) .lt. 0.9) then a(1,i)=0.0 a(2,i)=1.0 a(3,i)=0.0 d2f(i)=0.0 else a(1,i)=b(1) a(2,i)=b(2) a(3,i)=b(1) d2f(i)=b(3)*(f(i-1)-two*f(i)+f(i+1))+ & b(4)*(d1f5(i+1)-d1f5(i)) endif enddo a(1,n)=0.0 a(2,n)=1.0 a(3,n)=0.0 if(fsm(n)+fsm(n-1) .lt. 0.9) d2f(n)=0.0 call triapd(n,1,a,d2f) return end subroutine getdf(m,f,df,d2f,x) dimension f(*),df(*),d2f(*),x(*) h=x(2)-x(1) s=(x(3)-x(2))/h df(1)=((2.0-2.0*s**2-s**3)*(f(1)-f(2))+(f(3)-f(2))/s) & /(h*(s**3+2.0*s**2-1.0)) & +((2.0*s**2+s-1.0)*d2f(1)+(1.0+s)**2*d2f(2)) & /(6.0*(1.0-s-s**2))*h do i=2,m-1 h=x(i+1)-x(i) s=(x(i)-x(i-1))/h df(i)=((f(i-1)-f(i))/s+(2.*s**2+s**3)*(f(i+1)-f(i)))/ & (s**3+2.*s**2-1.0)/h - h*((3.*s+2.*s**2)*d2f(i) & +s**2*d2f(i+1))/(6.0*(s**2+s-1.0)) enddo h=x(m)-x(m-1) s=(x(m-1)-x(m-2))/h df(m)=((f(m-2)-f(m-1))/s+(2.0-2.0*s**2-s**3)*(f(m)-f(m-1))) & /(h*(1.0-2.0*s**2-s**3)) & +((1.0+s)**2*d2f(m-1)+(2.0*s**2+s-1.0)*d2f(m)) & /(6.0*(s**2+s-1.0))*h return end subroutine getd2fm(m,f,d2f,x,xl,xr,d2l,d2r,fsm) parameter (nn=200) implicit double precision (a-h, o-z) real*4 f(*),x(*),d2f(*),xl,xr,d2l,d2r,fsm(*) dimension a(4,nn),dx(nn),s(nn) do i=2,m dx(i)=x(i)-x(i-1) enddo do i=2,m-1 s(i)=dx(i+1)/dx(i) enddo if(fsm(1)+fsm(2) .lt. 0.9) then a(1,1)=0.0 a(2,1)=1.0 a(3,1)=0.0 a(4,1)=0.0 else tt=(xl-x(2))/dx(2) ss=s(2) s2=ss*ss t2=tt*tt sss=1.0-ss-s2 a(1,1) = 0.0 a(2,1)=(2.0*(1.0-ss)*t2+(1.0-ss+s2)*tt)/sss a(3,1)=(2.0*(2.0+ss)*t2+(5.0+ss-s2)*tt+sss)/sss a(4,1)=d2l-12.0*(t2+tt)*(f(1)-f(2)+(f(3)-f(2))/ss) & /(-(1.0+ss)*sss*dx(2)**2) endif do k=2,m-1 if(fsm(m-1)+fsm(k)+fsm(k+1) .lt. 0.9) then a(1,k)=0.0 a(2,k)=1.0 a(3,k)=0.0 a(4,k)=0.0 else a(1,k)= 1.0+s(k)-s(k)**2 a(2,k)=(1.0+4.0*s(k)+4.0*s(k)**2+s(k)**3)/s(k) a(3,k)=(s(k)**2+s(k)-1.0)/s(k) a(4,k)=12.0/dx(k)**2*((f(k+1)-f(k))/s(k)-(f(k)-f(k-1))) endif enddo if(fsm(m)+fsm(m-1) .lt. 0.9) then a(1,m)=0.0 a(2,m)=1.0 a(3,m)=0.0 a(4,m)=0.0 else tt=(xr-x(m-1))/dx(m) t2=tt*tt ss=dx(m-1)/dx(m) s2=ss*ss sss=s2+ss-1.0 a(1,m)=(sss-2.0*(2.0+ss)*t2+(5.0+ss-s2)*tt)/sss a(2,m)=(2.0*(ss-1.0)*t2+(1.0-ss+s2)*tt)/sss a(3,m) = 0.0 a(4,m)=d2r-12.0*(t2-tt)*((f(m-2)-f(m-1))/ss+f(m)-f(m-1)) & /((ss**3+2.0*s2-1.0)*dx(m)**2) endif do k=2,m w=a(1,k)/a(2,k-1) a(2,k)=a(2,k)-w*a(3,k-1) a(4,k)=a(4,k)-w*a(4,k-1) enddo a(4,m)=a(4,m)/a(2,m) do k=m-1,1,-1 a(4,k)=(a(4,k)-a(3,k)*a(4,k+1))/a(2,k) enddo do k=1,m d2f(k)=a(4,k) enddo return end subroutine getd2f(m,f,d2f,x,xl,xr,d2l,d2r) parameter (nn=200) implicit double precision (a-h, o-z) real*4 f(*),x(*),d2f(*),xl,xr,d2l,d2r dimension a(4,nn),dx(nn),s(nn) do i=2,m dx(i)=x(i)-x(i-1) enddo do i=2,m-1 s(i)=dx(i+1)/dx(i) enddo tt=(xl-x(2))/dx(2) ss=s(2) s2=ss*ss t2=tt*tt sss=1.0-ss-s2 a(1,1) = 0.0 a(2,1)=(2.0*(1.0-ss)*t2+(1.0-ss+s2)*tt)/sss a(3,1)=(2.0*(2.0+ss)*t2+(5.0+ss-s2)*tt+sss)/sss a(4,1)=d2l-12.0*(t2+tt)*(f(1)-f(2)+(f(3)-f(2))/ss) & /(-(1.0+ss)*sss*dx(2)**2) do k=2,m-1 a(1,k)= 1.0+s(k)-s(k)**2 a(2,k)=(1.0+4.0*s(k)+4.0*s(k)**2+s(k)**3)/s(k) a(3,k)=(s(k)**2+s(k)-1.0)/s(k) a(4,k)=12.0/dx(k)**2*((f(k+1)-f(k))/s(k)-(f(k)-f(k-1))) enddo tt=(xr-x(m-1))/dx(m) t2=tt*tt ss=dx(m-1)/dx(m) s2=ss*ss sss=s2+ss-1.0 a(1,m)=(sss-2.0*(2.0+ss)*t2+(5.0+ss-s2)*tt)/sss a(2,m)=(2.0*(ss-1.0)*t2+(1.0-ss+s2)*tt)/sss a(3,m) = 0.0 a(4,m)=d2r-12.0*(t2-tt)*((f(m-2)-f(m-1))/ss+f(m)-f(m-1)) & /((ss**3+2.0*s2-1.0)*dx(m)**2) do k=2,m w=a(1,k)/a(2,k-1) a(2,k)=a(2,k)-w*a(3,k-1) a(4,k)=a(4,k)-w*a(4,k-1) enddo a(4,m)=a(4,m)/a(2,m) do k=m-1,1,-1 a(4,k)=(a(4,k)-a(3,k)*a(4,k+1))/a(2,k) enddo do k=1,m d2f(k)=a(4,k) enddo return end function herm(t,f,df,d2f,dx,m) c c f'' f'' c f' f' c f f c |------- dx --------| c o-------------------o c i i+1 c c m: -1 integrate c 0 value c 1 first derivative c 2 second dericative c c herm=dx**(-m)*(para(i,1)*x(i)*f(1) + para(i,2)*x(i)*f(2) c + (para(i,3)*x(i)*df(1) + para(i,4)*x(i)*df(2))*dx c + (para(i,5)*x(i)*d2f(1)+para(i,6)*x(i)*d2f(2))*dx**2) c dimension para(6,6),x(6) dimension f(*),df(*),d2f(*) data para/1.0, 0.0, 0.0, -10.0, 15.0, -6.0, & 0.0, 0.0, 0.0, 10.0, -15.0, 6.0, & 0.0, 1.0, 0.0, -6.0, 8.0, -3.0, & 0.0, 0.0, 0.0, -4.0, 7.0, -3.0, & 0.0, 0.0, 0.5, -1.5, 1.5, -0.5, & 0.0, 0.0, 0.0, 0.5, -1.0, 0.5/ if(m.eq.-1) then do i=1,6 x(i)=t**i/i enddo elseif(m.eq.1) then x(1)=0.0 do i=2,6 x(i)=(i-1.0)*t**(i-2) enddo elseif(m.eq.2) then x(1)=0.0 x(2)=0.0 do i=3,6 x(i)=(i-1.0)*(i-2.0)*t**(i-3) enddo else do i=1,6 x(i)=t**(i-1) enddo endif herm=0.0 do i=1,6 herm=herm+x(i)*(para(i,1)*f(1)+para(i,2)*f(2) & +(para(i,3)*df(1) + para(i,4)*df(2))*dx & +(para(i,5)*d2f(1)+ para(i,6)*d2f(2))*dx**2) enddo herm=herm/dx**m return end subroutine getint(m,x,f,df,d2f,xx,fint,n) dimension x(*),f(*),df(*),d2f(*),xx(*),fint(*), & fi(2000) fi(1)=-herm(-x(1)/(x(2)-x(1)),f,df,d2f,x(2)-x(1),-1) do i=1,m-1 fi(i+1)= fi(i) & + herm(1.0,f(i),df(i),d2f(i),x(i+1)-x(i),-1) enddo do i=1,n k=inx(m,xx(i),x) fint(i)=fi(k)+herm((xx(i)-x(k))/(x(k+1)-x(k)), & f(k),df(k),d2f(k),x(k+1)-x(k),-1) enddo return end function inx(m,t,x) dimension x(*) idr=1 if(x(2).lt.x(1)) idr=-1 inx=1 ih=m do while(ih.gt.inx+1) k=(inx+ih)/2 if(idr*x(k) .gt. idr*t) then ih=k else inx=k endif enddo return end C