[SCM] Gerris Flow Solver branch, upstream, updated. b3aa46814a06c9cb2912790b23916ffb44f1f203
Stephane Popinet
popinet at users.sf.net
Fri May 15 02:56:04 UTC 2009
The following commit has been merged in the upstream branch:
commit 3a7ab2dba6fddca3de0e81114aa1fa27df15e399
Author: Stephane Popinet <popinet at users.sf.net>
Date: Mon Jul 28 16:03:05 2008 +1000
Iterative Stokes wave solver of John Chaplin
darcs-hash:20080728060305-d4795-3347b6910bd9e7dfbc03937ce1d207cb11e7d219.gz
diff --git a/modules/CW263.f b/modules/CW263.f
new file mode 100644
index 0000000..20332dc
--- /dev/null
+++ b/modules/CW263.f
@@ -0,0 +1,887 @@
+ program CW263
+c
+c Stream function wave theory code:
+c * automatic selection of order
+c * uniform current
+c * returns proportion of limiting height
+c j.r.chaplin at soton.ac.uk
+c
+c Subroutine cw260 solves the wave. Kinematics are then available
+c through subtroutine kmts as illustrated below in the main program.
+c Stores results in CW263.PSI in a format compatible with that used
+c by CW6.FOR (24/3/99)
+c
+ character ans*1
+c
+ 10 write(*,'(a)') ' Water depth (m) = '
+ read*,d
+ write(*,'(a)') ' Period (s) = '
+ read*,t
+ write(*,'(a)') ' Wave height (m) = '
+ read*,h
+ write(*,'(a)') ' Current (m/s) = '
+ read*,u
+ nverb= 1
+c
+ call cw260(d,t,h,u,nverb,n,el)
+c
+ 20 write(*,'(/a)')
+ : ' (H)orizontal, (V)ertical, (S)urface, (N)ew wave, (Q)uit : '
+ read(*,'(a)') ans
+ goto (20,21,22,23,10,29) (index('HhVvSsNnQq',ans)+3)/2
+c
+c Horizontal
+c
+ 21 write(*,'(a)') ' y (m) = '
+ read*,yy
+ npt= 21
+ write(*,3)
+ do 31 i=1,npt
+ xx= el*(i-1)/float(npt-1)
+ tt= 0.0
+ call kmts(n,xx,yy,tt,uu,vv,ut,vt,du,dv,etah)
+ ans= ' '
+ if (yy.gt.etah) ans='*'
+ 31 write(*,'(f9.3,8f8.3,1x,a)') xx,yy,uu,vv,ut,vt,du,dv,etah,ans
+ goto 20
+c
+c Vertical
+c
+ 22 write(*,'(a)') ' x/L = '
+ read*,xl
+ xx= xl*el
+ yy= 0.0
+ tt= 0.0
+ call kmts(n,xx,yy,tt,uu,vv,ut,vt,du,dv,etah)
+ npt= 21
+ write(*,3)
+ do 32 i=1,npt
+ yy= (d+etah)*(npt-i)/float(npt-1)-d
+ call kmts(n,xx,yy,tt,uu,vv,ut,vt,du,dv,etah)
+ 32 write(*,'(f9.3,8f8.3,1x,1a)') xx,yy,uu,vv,ut,vt,du,dv,etah
+ goto 20
+c
+c Free surface
+c
+ 23 npt= 21
+ write(*,3)
+ do 33 i=1,npt
+ xx= el*(i-1)/float(npt-1)
+ yy= h
+ tt= 0.0
+ call kmts(n,xx,yy,tt,uu,vv,ut,vt,du,dv,etah)
+ 33 write(*,'(f9.3,8f8.3)') xx,etah,uu,vv,ut,vt,du,dv,etah
+ goto 20
+c
+ 29 stop
+ 3 format(' x y u v ut ',
+ : ' vt du dv eta')
+ end
+c
+c
+c
+ subroutine cw260(zd,zt,zh,zu,nverb,nfun,zel)
+c
+c Input: zd=depth; zt=period; zh=height; u=current; nverb=verbosity
+c Output: nfun=order; zel=wavelength
+c
+ implicit double precision (a-h,o-z)
+ parameter (nmax=25)
+ real zd,zt,zh,zu,zel
+ double precision k
+ character itl*79,datim*22
+ common /one/ d,t,h,u,k
+ : /two/ eta(nmax),c(nmax),amp(0:nmax)
+c
+ pi= 4*atan(1.0)
+c
+ d= zd
+ t= zt
+ hw= zh
+ u= zu
+c
+c Get required solution order
+c
+ call wavecel(t,d,u,tr,cel)
+ call limit(hw,d,tr,rat,1)
+ dl0= d/(9.81*t*t/(2*pi))
+ a= 0.86/sqrt(dl0)
+ b= 7+2.2*log(dl0)
+ cc= 2.7-3*log(dl0)
+ nve= nint((a+b*rat+cc*rat**2)/2)*2
+ nve= nve+2
+ hb= hw/rat
+c
+ k= 2*pi/(cel*t)
+c
+c Start with height=<Hb/2
+c First step up orders; then heights ...
+c
+ if (nverb.ne.0) write(*,'(/2a)')
+ : ' d T H U order iter ',
+ : 'rms error code L'
+ nfun = 6
+ h= min(hw,0.5*hb)
+ do 10 i=0,nmax
+ 10 amp(i)= 0.0
+ amp(1)= h/2
+ 12 continue
+ call cw261(nfun,iter,fsumsq,ifail)
+ if (ifail.eq.0) then
+ el= 2*pi/k
+ if (nverb.ne.0) write(*,'(a,4f8.3,i5,i7,1p1e12.3,0p,i5,f10.3)')
+ : ' ',d,t,h,u,nfun,iter,fsumsq,ifail,el
+ else
+ if (nverb.ne.0) write(*,'(a,4f8.3,i5,i7,1p1e12.3,0p,i5)')
+ : ' ',d,t,h,u,nfun,iter,fsumsq,ifail
+ stop
+ endif
+ if (nfun.lt.nve) then
+ nfun= nfun+2
+ goto 12
+ endif
+c
+ if (hw.gt.0.5*hb) then
+ fac= 1.1
+ 11 hm= h
+ if (h*fac.gt.hw) then
+ h= hw
+ ilast= 1
+ else
+ h= h*fac
+ ilast= 0
+ endif
+ do 14 i=1,nfun-1
+ 14 amp(i)= (h/hm)*amp(i)
+ call cw261(nfun,iter,fsumsq,ifail)
+ if (ifail.eq.0) then
+ el= 2*pi/k
+ if (nverb.ne.0)
+ : write(*,'(a,4f8.3,i5,i7,1p1e12.3,0p,i5,f10.3)')
+ : ' ',d,t,h,u,nfun,iter,fsumsq,ifail,el
+ else
+ if (nverb.ne.0) write(*,'(a,4f8.3,i5,i7,1p1e12.3,0p,i5)')
+ : ' ',d,t,h,u,nfun,iter,fsumsq,ifail
+ stop
+ endif
+ if (ilast.ne.1) then
+ fac= 0.995*fac
+ goto 11
+ endif
+ endif
+ zel= 2*pi/k
+c
+ open(12,file='cw263.psi')
+ itl= 'Solved by CW263'
+ datim= 'Space for date & time:'
+ ver= 5.01
+ el0= 9.81*zt**2/(2*pi)
+ hl0= zh/el0
+ dl0= zd/el0
+ wl0= zel/el0
+ write(12,101) itl,datim,ver,hl0,dl0,wl0,zh,zd,zt,nfun-1,zu
+ write(12,102) (eta(i),c(i+1),i=1,nfun-1),eta(nfun)
+ write(12,107) (amp(i)/zh,i=0,nfun-1)
+ close(12)
+c
+ 101 format(a/a,f10.2/3f16.10/3f16.10,i5,f16.10)
+ 102 format(1p2e25.16)
+ 107 format(1p1e25.16)
+c
+ return
+ end
+c
+c
+c
+ subroutine cw261(nfun,iter,fsumsq,ifail)
+ implicit double precision (a-h,o-z)
+ parameter (nmax=25)
+ double precision k
+ dimension x(nmax),f(nmax),etas(2*nmax),bb(0:nmax)
+ common /one/ d,t,h,u,k
+ : /two/ eta(nmax),c(nmax),amp(0:nmax)
+c
+ pi= 4*atan(1.0)
+c
+ do 10 i=1,nfun-2
+ th= (i-1)*pi/(nfun-1)
+ x(i)= 0.0
+ do 10 j=1,nfun-1
+ 10 x(i)= x(i)+cos(j*th)*amp(j)
+ x(nfun-1)= k
+ jverb= 0
+ call gaf(nfun,nfun-1,x,f,fsumsq,jverb,iter,ifail)
+ if (ifail.eq.1) return
+ k= x(nfun-1)
+ do 11 i=1,nfun
+ 11 etas(i)= eta(i)
+ do 12 i=1,nfun-1
+ 12 etas(nfun+i)= eta(nfun-i)
+ call four(etas,2*nfun-2,amp,bb,nfun-1)
+ amp(nfun)= 0.0
+ bb(nfun)= 0.0
+ return
+ end
+c
+c
+c
+ subroutine lsfun(nfun,x,f,jac,sq,ifail)
+c
+c Computes function errors and Jacobian
+c x(1:nfun-1) = unknowns, f(1:nfun) = functions
+c
+ implicit double precision (a-h,o-z)
+ parameter (nmax=25)
+ double precision k,jac,ks
+ dimension x(nmax),f(nmax),jac(nmax,nmax)
+ dimension dex(nmax,nmax),th(nmax),s(nmax),
+ : chcs(nmax,nmax),chsn(nmax,nmax),shcs(nmax,nmax),shsn(nmax,nmax),
+ : fu(nmax,nmax),df(nmax,nmax),a(nmax),da(nmax,nmax),
+ : dax(nmax,nmax),dcx(nmax,nmax),dfx(nmax,nmax),
+ : us(nmax),vs(nmax),dux(nmax,nmax),dvx(nmax,nmax)
+ common /one/ d,t,h,u,k
+ : /two/ eta(nmax),c(nmax),amp(0:nmax)
+c
+ pi= 4*atan(1.0d0)
+ g= 9.81d0
+ ifail= 0
+ om= 2*pi/t
+c
+c Get surface elevations
+c
+ do 40 i=1,nfun-2
+ 40 eta(i)= x(i)
+ eta(nfun)= eta(1)-h
+ sm= 0.0
+ do 10 i=2,nfun-2
+ 10 sm= sm+eta(i)
+ eta(nfun-1)= -sm-(eta(1)+eta(nfun))/2
+ do 11 i=2,nfun
+ if (eta(i).gt.eta(i-1)+5*h/nfun) then
+ ifail= 1
+c write(*,'(7f11.3)') (eta(m),m=1,nfun-1)
+ return
+ endif
+ 11 continue
+ k= x(nfun-1)
+c
+c Get dex(i,j) = d(eta(i))/d(eta(j)), 1 <= j <= nfun-2
+c
+ do 18 i=1,nfun
+ dex(i,nfun-1)= 0.0d0
+ do 18 j=1,nfun-2
+ if (i.lt.nfun-1) then
+ if (i.eq.j) then
+ dex(i,j)= 1.0d0
+ else
+ dex(i,j)= 0.0d0
+ endif
+ else if (i.eq.nfun-1) then
+ dex(i,j)= -1.0d0
+ else if (i.eq.nfun) then
+ if (j.eq.1) then
+ dex(i,j)= 1.0d0
+ else
+ dex(i,j)= 0.0d0
+ endif
+ endif
+ 18 continue
+c
+c Set some useful functions
+c
+ do 23 i=1,nfun
+ th(i)= pi*(i-1)/float(nfun-1)
+ s(i)= d+eta(i)
+ ks= k*s(i)
+ do 23 n=1,nfun
+ ch= cosh(n*ks)
+ sh= sinh(n*ks)
+ cs= cos(n*th(i))
+ sn= sin(n*th(i))
+ chcs(i,n)= ch*cs
+ shsn(i,n)= sh*sn
+ chsn(i,n)= ch*sn
+ 23 shcs(i,n)= sh*cs
+c
+c Get normalised stream function coefficients a(j)
+c
+ do 12 i=1,nfun
+ fu(i,1)= k*s(i)
+ df(i,1)= k
+ do 12 j=2,nfun
+ n= j-1
+ fu(i,j)= shcs(i,n)
+ 12 df(i,j)= chcs(i,n)*n*k
+ call trans2(nfun,fu,df,a,da)
+c
+c Get dax(i,j) = d(a(i))/d(eta(j)), 1 <=j <= nfun-2
+c
+ do 41 i=1,nfun
+ do 41 j=1,nfun-2
+ dij= 0.0
+ do 42 n=1,nfun
+ 42 dij= dij+da(i,n)*dex(n,j)
+ 41 dax(i,j)= dij
+c
+c Get dax(i,nfun-1) = d(a(i))/dk
+c
+ do 19 i=1,nfun
+ sm= 0.0
+ do 13 m=1,nfun
+ 13 sm= sm+da(i,m)*s(m)/k
+ 19 dax(i,nfun-1)= sm
+c
+c De-normalise to get correct reverse mean flow
+c
+ a1k2= a(1)*k**2
+ r= (u*k-om)/a1k2
+ do 24 i=1,nfun
+ c(i)= a(i)*r
+ do 24 j=1,nfun-1
+ drx= -(u*k-om)*dax(1,j)*(k/a1k2)**2
+ if (j.eq.nfun-1) then
+ drx= drx+u/a1k2-2*a(1)*k*(u*k-om)/a1k2**2
+ endif
+ 24 dcx(i,j)= r*dax(i,j)+a(i)*drx
+c
+c Get surface velocities and functions
+c
+ do 14 i=1,nfun
+ su= k*c(1)
+ sv= 0.0d0
+ do 15 n=1,nfun-1
+ su= su+n*k*c(n+1)*chcs(i,n)
+ 15 sv= sv+n*k*c(n+1)*shsn(i,n)
+ us(i)= su
+ 14 vs(i)= sv
+ sm= 0.0d0
+ do 16 i=1,nfun
+ f(i)= (us(i)**2+vs(i)**2)/(2*g)+eta(i)
+ 16 sm= sm+f(i)
+ sm= sm/float(nfun)
+ sq= 0.0d0
+ do 17 i=1,nfun
+ f(i)= f(i)-sm
+ 17 sq= sq+f(i)**2
+ sq= sqrt(sq/nfun)/h
+c
+c Get d(u(i))/d(eta(j)) and d(u(i))/dk
+c
+ do 20 i=1,nfun
+ do 21 j=1,nfun-2
+ su= k*dcx(1,j)
+ sv= 0.0d0
+ do 22 n=1,nfun-1
+ su= su+n*k*dcx(n+1,j)*chcs(i,n)+(n*k)**2*c(n+1)*shcs(i,n)*dex(i,j)
+ 22 sv= sv+n*k*dcx(n+1,j)*shsn(i,n)+(n*k)**2*c(n+1)*chsn(i,n)*dex(i,j)
+ dux(i,j)= su
+ 21 dvx(i,j)= sv
+ j= nfun-1
+ su= c(1)+k*dcx(1,j)
+ sv= 0.0d0
+ do 25 n=1,nfun-1
+ su= su+(c(n+1)+k*dcx(n+1,j))*n*chcs(i,n)+
+ : n**2*k*s(i)*c(n+1)*shcs(i,n)
+ 25 sv= sv+(c(n+1)+k*dcx(n+1,j))*n*shsn(i,n)+
+ : n**2*k*s(i)*c(n+1)*chsn(i,n)
+ dux(i,j)= su
+ 20 dvx(i,j)= sv
+c
+c Get derivatives of functions
+c
+ do 44 i=1,nfun
+ do 44 j=1,nfun-1
+ 44 dfx(i,j)= (us(i)*dux(i,j)+vs(i)*dvx(i,j))/g+dex(i,j)
+c
+ do 45 j=1,nfun-1
+ sm= 0.0
+ do 46 i=1,nfun
+ 46 sm= sm+dfx(i,j)
+ sm= sm/nfun
+ do 45 i=1,nfun
+ 45 jac(i,j)= dfx(i,j)-sm
+c
+ return
+ end
+c
+c
+c
+ subroutine trans2(n,f,df,a,da)
+c
+c Gramm-Schmidt orthogonalisation with derivatives
+c
+ implicit double precision (a-h,o-z)
+ parameter (nmax=25)
+ dimension f(nmax,nmax),df(nmax,nmax),a(nmax),da(nmax,nmax)
+ dimension c(nmax,nmax),b(nmax),cr(nmax,nmax),g(nmax,nmax),
+ : dly(nmax),dc(nmax,nmax),br(nmax,nmax),db(nmax)
+ do 10 i=1,n
+ do 10 j=1,n
+ 10 c(i,j)= 0.0d0
+c
+ st= 0.0d0
+ do 54 i=1,n
+ 54 st= st+f(i,1)**2
+ st= dsqrt(st)
+c
+ do 16 i=1,n
+ 16 g(i,1)= f(i,1)/st
+ c(1,1)= st
+ do 18 k=2,n
+ km= k-1
+ do 29 j=1,km
+ st= 0.0d0
+ do 20 i=1,n
+ 20 st= st+g(i,j)*f(i,k)
+ 29 c(j,k)= st
+ c(k,k)= 1.0d0
+ do 21 i=1,n
+ st= 0.0d0
+ do 22 j=1,km
+ 22 st= st+c(j,k)*g(i,j)
+ 21 g(i,k)= f(i,k)-st
+ st= 0.0d0
+ do 50 i=1,n
+ 50 st= st+g(i,k)**2
+ st= dsqrt(st)
+ do 52 i=1,n
+ 52 g(i,k)= g(i,k)/st
+ c(k,k)= st
+ do 200 j=1,km
+ st= 0.0d0
+ do 202 i=1,n
+ 202 st= st+g(i,j)*g(i,k)
+ 200 db(j)= st
+ do 204 i=1,n
+ st= 0.0d0
+ do 206 j=1,km
+ 206 st= st+db(j)*g(i,j)
+ 204 g(i,k)=g(i,k)-st
+ 18 continue
+ do 24 j=1,n
+ st= 0.0d0
+ do 26 i=1,n
+ 26 st= st+g(i,j)
+ 24 b(j)= st
+ st= 0.0d0
+ do 208 i=1,n
+ sb= 0.0d0
+ do 210 j=1,n
+ 210 sb= sb+b(j)*g(i,j)
+ 208 st= st+(sb-1.0d0)**2
+ call trinv(n,c,cr)
+ do 28 i=1,n
+ st= 0.0d0
+ do 30 j=1,n
+ 30 st= st+cr(i,j)*b(j)
+ 28 a(i)= st
+ st= 0.0d0
+ do 72 i=1,n
+ sb= 0.0d0
+ do 74 j=1,n
+ 74 sb= sb+a(j)*f(i,j)
+ 72 st= st+(sb-1.0d0)**2
+ st= dsqrt(st/float(n))
+ do 80 k=1,n
+ do 81 i=1,n
+ do 81 j=1,n
+ 81 dc(i,j)= 0.0d0
+ dly(1)= 2.0d0*f(k,1)*df(k,1)
+ dc(1,1)= 0.5d0*dly(1)/c(1,1)
+ do 82 j=2,n
+ 82 dc(1,j)= (f(k,j)*df(k,1)+df(k,j)*f(k,1))/c(1,1)
+ : -0.5d0*dly(1)*c(1,j)/c(1,1)**2
+ do 84 l=2,n
+ lm= l-1
+ st= 0.0d0
+ do 86 m=1,lm
+ 86 st= st+c(m,l)*dc(m,l)
+ dly(l)= 2.0d0*(f(k,l)*df(k,l)-st)
+ dc(l,l)= 0.5d0*dly(l)/c(l,l)
+ if (l.eq.n) goto 84
+ lp= l+1
+ do 85 m=lp,n
+ st= 0.0d0
+ do 90 j=1,lm
+ 90 st= st+c(j,m)*dc(j,l)+dc(j,m)*c(j,l)
+ 85 dc(l,m)= (f(k,m)*df(k,l)+df(k,m)*f(k,l))/c(l,l)
+ : -st/c(l,l)-0.5d0*dly(l)*c(l,m)/c(l,l)**2
+ 84 continue
+ do 92 l=1,n
+ do 92 m=1,n
+ st= 0.0d0
+ do 96 j=1,n
+ 96 st= st+g(l,j)*dc(j,m)
+ 92 br(l,m)= -st
+ do 94 m=1,n
+ 94 br(k,m)= br(k,m)+df(k,m)
+ do 32 j=1,n
+ sb= 0.0d0
+ do 34 i=1,n
+ st= 0.0d0
+ do 36 l=1,n
+ 36 st= st+br(i,l)*cr(l,j)
+ 34 sb= sb+st
+ 32 db(j)= sb
+ do 38 i=1,n
+ st= db(i)
+ do 40 j=1,n
+ 40 st= st-dc(i,j)*a(j)
+ 38 dly(i)= st
+ do 42 i=1,n
+ st= 0.0d0
+ do 44 j=1,n
+ 44 st= st+cr(i,j)*dly(j)
+ 42 da(i,k)= st
+ 80 continue
+ return
+ end
+c
+c
+c
+ subroutine trinv(n,c,w)
+c
+c 3-diagonal system solver
+c
+ implicit double precision (a-h,o-z)
+ parameter (nmax=25)
+ dimension c(nmax,nmax),w(nmax,nmax)
+ do 10 i=1,n
+ do 10 j=1,n
+ 10 w(i,j)= 0.0d0
+ do 12 jj=1,n
+ j= n+1-jj
+ w(j,j)= 1.0d0
+ if (j.eq.1) goto 12
+ k= j-1
+ do 13 ii=1,k
+ i= k-ii+1
+ sum= 0.0d0
+ do 14 l=1,ii
+ m= n+2-jj-l
+ 14 sum= sum+c(i,m)*w(m,j)/c(m,m)
+ 13 w(i,j)= -sum
+ 12 continue
+ do 16 i=1,n
+ sum= c(i,i)
+ do 16 j=1,n
+ 16 w(i,j)= w(i,j)/sum
+ return
+ end
+c
+c
+c
+ subroutine wavecel(ta,d,u,tr,c)
+c
+c Linear theory C by series approximation for waves on a current u
+c d = still water depth
+c ta = absolute period (fixed reference frame)
+c tr = relative period (reference frame moving with the current)
+c c = celerity (reference frame moving with the current)
+c
+ implicit double precision (a-h,o-z)
+ pi= 4*atan(1.0d0)
+ g= 9.81
+ sigma= 2*pi/ta
+ y= sigma*sigma*d/g
+ a= 1.0/(1.0+y*(0.6667+y*(0.3556+y*(0.1608+y*(0.06321+
+ : y*(0.02744+y*0.01))))))
+ c= dsqrt((d*g)/(y+a))
+ if (abs(u).lt.1.0d-6) then
+ tr= ta
+ return
+ else
+ el= c*ta
+ iter= 0
+ 10 tr= el/(el/ta-u)
+ elp= (g*tr**2/(2*pi))*tanh(2*pi*d/el)
+ del= elp-el
+ el= el+del*0.4
+ if (abs(del/el).gt.1.0d-6) then
+ iter= iter+1
+ if (iter.eq.100) then
+ write(*,'(a)') ' WAVECEL error'
+ stop
+ endif
+ goto 10
+ endif
+ tr= el/(el/ta-u)
+ c= el/tr
+ return
+ endif
+ end
+c
+c
+c
+ subroutine slpds(a,b,n,cc)
+c
+c Linear system solution
+c
+ implicit double precision (a-h,o-z)
+ parameter (nmax=25)
+ dimension a(nmax,nmax),b(nmax),cc(nmax)
+ n1= n-1
+ do 100 k=1,n1
+ c= a(k,k)
+ k1= k+1
+ if (abs(c).lt.1.0d-10) then
+ write(*,'(a,i5)') ' Matrix error 1: ',k
+ stop
+ endif
+ do 4 j=k1,n
+ 4 a(k,j)= a(k,j)/c
+ b(k)= b(k)/c
+ do 10 i=k1,n
+ c= a(i,k)
+ do 5 j=k1,n
+ 5 a(i,j)= a(i,j)-c*a(k,j)
+ 10 b(i)= b(i)-c*b(k)
+ 100 continue
+ if (abs(a(n,n)).lt.1.0d-10) then
+ write(*,'(a,i5)') ' Matrix error 2: ',n
+ stop
+ endif
+ b(n)= b(n)/a(n,n)
+ do 200 l=1,n1
+ k= n-l
+ k1= k+1
+ do 200 j=k1,n
+ 200 b(k)= b(k)-a(k,j)*b(j)
+ do 300 i=1,n
+ 300 cc(i)= b(i)
+ return
+ end
+c
+c
+c
+ subroutine gaf(nf,nv,xc,fvecc,fsumsq,nverb,iter,ifail)
+c
+c Non-linear system error minimisation
+c ifail = 0: OK
+c 1: looks hopeless
+c 2: poor convergence
+c
+ implicit double precision (a-h,o-z)
+ parameter (nmax=25)
+ dimension xc(nmax),fvecc(nmax)
+ dimension fjacc(nmax,nmax),
+ : aa(nmax,nmax),bb(nmax),cc(nmax),xcm(nmax)
+c
+ iter = 0
+ icalls= 0
+ orf= 2.0/nf
+ fsumsm= 100.0
+c
+ 15 do 16 i=1,nv
+ 16 xcm(i)= xc(i)
+ 19 call lsfun(nf,xc,fvecc,fjacc,fsumsq,ifl)
+ icalls= icalls+1
+ if ((fsumsq.gt.fsumsm.and.iter.gt.1).or.ifl.ne.0) then
+ if (orf.lt.0.05) then
+ ifail= 1
+ return
+ endif
+ orf= 0.8*orf
+ iter= max(iter-1,0)
+ do 18 i=1,nv
+ 18 xc(i)= xcm(i)
+ fsumsm= 100.0
+ goto 19
+ endif
+ call monit(nf,fvecc,icalls,nverb)
+c
+ do 10 i=1,nv
+ do 10 j=1,nv
+ ra= 0.0d0
+ do 11 l=1,nf
+ 11 ra= ra+fjacc(l,i)*fjacc(l,j)
+ 10 aa(i,j)= ra
+c
+ do 12 i=1,nv
+ ra= 0.0d0
+ do 13 l=1,nf
+ 13 ra= ra+fvecc(l)*fjacc(l,i)
+ 12 bb(i)= -ra
+c
+ call slpds(aa,bb,nv,cc)
+c
+ dxmax= 0.0d0
+ do 14 i=1,nv
+ dxmax= max(dxmax,abs(cc(i)))
+ 14 xc(i)= xc(i)+orf*cc(i)
+c
+ iter= iter+1
+ fsumsm= fsumsq
+ orf= min(1.0d0,orf*1.1)
+c
+ if (iter.ge.50.and.fsumsq.lt.1.0d-4) then
+ ifail= 2
+ return
+ endif
+ if (iter.ge.50) then
+ ifail= 1
+ return
+ endif
+ if (fsumsq.gt.1.0d-6) goto 15
+ ifail= 0
+ return
+ end
+c
+c
+c
+ subroutine monit(nfun,f,icalls,nverb)
+ implicit double precision (a-h,o-z)
+ double precision k
+ parameter (nmax=25)
+ dimension f(nmax)
+ common /one/ d,t,h,u,k
+ : /two/ eta(nmax),c(nmax),amp(0:nmax)
+c
+c Outputs monitoring data
+c
+ if (nverb.eq.0) return
+ sm= 0.0d0
+ do 10 i=1,nfun
+ 10 sm=sm+f(i)**2
+ sq= sqrt(sm/nfun)/h
+ write(*,'(/i11,1p1e11.3)') icalls,sq
+ write(*,'(1p7e11.3)') (eta(i),i=1,nfun)
+ return
+ end
+c
+c
+c
+ subroutine limit(h,d,t,rat,nverb)
+ double precision h,d,t,rat
+c
+c Estimates H/(limiting height for d and t) = rat
+c
+ dimension dl0(18),hl0(18)
+ data dl0/2,0.578,0.440,0.356,0.293,0.243,0.201,0.166,0.1359,
+ : 0.1100,0.0876,0.0686,0.0524,0.0390,0.0277,0.01879,
+ : 0.01168,0.00638/
+ data hl0/0.1682,0.1665,0.1613,0.1531,0.1423,0.1298,0.1159,
+ : 0.1017,0.0873,0.0735,0.0605,0.0487,0.0380,0.0289,
+ : 0.0208,0.01440,0.00911,0.00501/
+ pi= 4*atan(1.0)
+ el0= 9.81*t**2/(2*pi)
+ ha= h/el0
+ da= d/el0
+ if (da.gt.dl0(1)) then
+ rat= ha/hl0(1)
+ else if (da.lt.dl0(18)) then
+ rat= ha/(0.8*da)
+ else
+ do 10 i=2,18
+ if (dl0(i).lt.da) goto 11
+ 10 continue
+ 11 x1= log(dl0(i))
+ x2= log(dl0(i-1))
+ y1= log(hl0(i))
+ y2= log(hl0(i-1))
+ r= (log(da)-x1)/(x2-x1)
+ hb= exp(y1+r*(y2-y1))
+ rat= ha/hb
+ endif
+c
+ if (nverb.ne.0.or.rat.gt.1.0) then
+ write(*,'(a,f5.3)') ' H/Hb = ',rat
+ endif
+ if (rat.gt.1.0) stop
+c
+ return
+ end
+c
+c
+c
+ subroutine four(f,n,a,b,nb)
+c
+c Fourier analysis
+c
+ implicit double precision (a-h,o-z)
+ dimension a(0:nb),b(0:nb),f(n)
+ pi= 4*atan(1.0d0)
+ rn= 2.0d0/n
+ t= 2*pi/n
+ c= cos(t)
+ s= sin(t)
+ vk= 0.0d0
+ vl= -1.0d0
+ do 10 k=0,nb
+ t= c*vk
+ ck= t-vl
+ vl= vk
+ vk= ck+t
+ t= ck+ck
+ ul= 0.0d0
+ um= f(n)
+ do 12 mm=3,n
+ m= n+2-mm
+ u0= ul
+ ul= um
+ 12 um= t*ul-u0+f(m)
+ a(k)= (ck*um-ul+f(1))*rn
+ 10 b(k)= s*vl*um*rn
+ a(0)= a(0)*0.5d0
+ if (2*nb.ne.n) return
+ a(nb)= a(nb)*0.5d0
+ b(nb)= 0.0d0
+ return
+ end
+c
+c
+c
+ subroutine kmts(nfun,xx,yy,tt,uu,vv,ut,vt,du,dv,etah)
+c
+c Computes
+c horizontal and vertical velocity components (u,v)
+c horizontal and vertical local acceleration components (ut,vt)
+c horizontal and vertical total acceleration components (du,dv)
+c water surface elevation (etah)
+c at t=tt; x=xx; y=yy.
+c If yy>eta kinematics are returned at the free surface
+c
+ implicit double precision (a-h,o-z)
+ real xx,yy,tt,uu,vv,ut,vt,du,dv,etah
+ parameter (nmax=25)
+ double precision k,ks
+ common /one/ d,t,h,u,k
+ : /two/ eta(nmax),c(nmax),amp(0:nmax)
+ pi=4*atan(1.0d0)
+ om= 2*pi/t
+ theta= k*xx-om*tt
+c
+ etah= 0.0
+ do 11 i=1,nfun-1
+ 11 etah= etah+cos(i*theta)*amp(i)
+c
+ ks= k*(d+min(yy,etah))
+ s1= 0.0
+ s2= 0.0
+ s3= 0.0
+ s4= 0.0
+c
+ do 10 i=1,nfun-1
+ ip= i+1
+ ch= cosh(i*ks)
+ sh= sinh(i*ks)
+ cs= cos(i*theta)
+ sn= sin(i*theta)
+ s1= s1+i*ch*cs*c(ip)
+ s2= s2+i*sh*sn*c(ip)
+ s3= s3+i*i*ch*sn*c(ip)
+ 10 s4= s4+i*i*sh*cs*c(ip)
+ uu= u+k*s1
+ vv= k*s2
+ ut= k*om*s3
+ vt= -k*om*s4
+ ux= -k*k*s3
+ vx= k*k*s4
+ uy= vx
+ vy= -ux
+ du= ut+uu*ux+vv*uy
+ dv= vt+uu*vx+vv*vy
+c
+ return
+ end
+
--
Gerris Flow Solver
More information about the debian-science-commits
mailing list