module serv_xnl4v5
contains
SUBROUTINE y_gauleg(x1,x2,x,w,n)
!-------------------------------------------------------------------
INTEGER, intent(in) :: n ! Number of intervals
real, intent(in) :: x1 ! lower limit of integration interval
real, intent(in) :: x2 ! upper limit of integration interval
real, intent(out) :: x(n) ! Positions for function evaluations
real, intent(out) :: w(n) ! Weights
!
!-----------------------------------------------------------------------
DOUBLE PRECISION EPS
PARAMETER (EPS=3.d-14)
INTEGER i,j,m
DOUBLE PRECISION p1,p2,p3,pp,xl,xm,z,z1
!-----------------------------------------------------------------------
m=(n+1)/2
xm=0.5d0*(x2+x1)
xl=0.5d0*(x2-x1)
do 12 i=1,m
z=cos(3.141592654d0*(i-.25d0)/(n+.5d0))
1 continue
p1=1.d0
p2=0.d0
do 11 j=1,n
p3=p2
p2=p1
p1=((2.d0*j-1.d0)*z*p2-(j-1.d0)*p3)/j
11 continue
pp=n*(z*p1-p2)/(z*z-1.d0)
z1=z
z=z1-p1/pp
if(abs(z-z1).gt.EPS)goto 1
x(i)=xm-xl*z
x(n+1-i)=xm+xl*z
w(i)=2.d0*xl/((1.d0-z*z)*pp*pp)
w(n+1-i)=w(i)
12 continue
!
return
END subroutine
!-----------------------------------------------------------------------------!
subroutine z_cmpcg(sigma,depth,grav,cg)
!-----------------------------------------------------------------------------!
!
! +-------+ ALKYON Hydraulic Consultancy & Research
! | | Gerbrant van Vledder
! | +---+
! | | +---+
! +---+ | |
! +---+
!
implicit none
!
! 0. Update history
!
! 12/01/2001 Initial version
! 11/04/2001 Check included for the cases tat sigma < 0 or depth <0
! Result is cg = -10
!
! 1. Purpose:
!
! Compute group velocity for a given radian frequency and depth
!
! 2. Method
!
! Linear wave theory
!
! 3. Parameter list:
!
!Type I/O Name Description
!------------------------------------------------------------------------------
real, intent(in) :: sigma ! radian frequency (rad)
real, intent(in) :: depth ! water depth (m)
real, intent(in) :: grav ! gravitational acceleration (m/s^2)
real, intent(out) :: cg ! group velocity (m/s)
!
real k ! wave number
!/A
!!real z_wnumb ! compute wave number
!/Z
!-----------------------------------------------------------------------------
k = z_wnumb(sigma,depth,grav)
!
if(depth <= 0. .or. sigma <= 0.) then
cg = -10.
else
if(depth*k > 30.) then
cg = grav/(2.*sigma)
else
cg = sigma/k*(0.5+depth*k/sinh(2.*depth*k))
end if
end if
!
return
end subroutine
!
!-----------------------------------------------------------------------------!
subroutine z_intp1(x1,y1,x2,y2,n1,n2,ierr) !
!-----------------------------------------------------------------------------!
!
! +-------+ ALKYON Hydraulic Consultancy & Research
! | | Gerbrant van Vledder
! | +---+
! | | +---+
! +---+ | |
! +---+
!
implicit none
!
! 0. Update history
!
! 30/03/1999 Initical version
! 9/04/1999 Check included for monotonicity of x-data
! 11/10/1999 Error messages added and updated
! 18/01/2001 Check include if n1==1
! 24/01/2001 Check for equality of y2 data loosened if n2==1
! 13/09/2001 Documentation updated
!
! 1. Purpose
!
! Interpolate function values
!
! 2. Method
!
! Linear interpolation
! If a requested point falls outside the input domain, then
! the nearest point is used (viz. begin or end point of x1/y1 array
!
! If the input array has only one point. A constant value is assumed
!
! 3. Parameter list
!
! Name I/O Type Description
!
integer, intent(in) :: n1 ! number of data points in x1-y1 arrays
integer, intent(in) :: n2 ! number of data points in x2-y2 arrays
real, intent(in) :: x1(n1) ! x-values of input data
real, intent(in) :: y1(n1) ! y-values of input data
real, intent(in) :: x2(n2) ! x-values of output data
real, intent(out) :: y2(n2) ! y-values of output data
integer, intent(out) :: ierr ! Error indicator
!
! 4. Subroutines used
!
! 5. Error messages
!
! ierr = 0 No errors detected
! = 1 x1-data not monotonic increasing
! = 10 x2-data not monotonic increasing
! = 11 x1- and x2 data not monotonic increasing
! = 2 x1-data not monotonic decreasing
! = 20 x1-data not monotonic decreasing
! = 22 x1- and x2 data not monotonic decreasing
!
! = 2 No variation in x1-data
! = 3 No variation in x2-data is allowed if n2=1
!
! 6. Remarks
!
! It is assumed that the x1- and x2-data are either
! monotonic increasing or decreasing
!
! If a requested x2-value falls outside the range of x1-values
! it is assumed that the corresponding y2-value is equal to
! the nearest boundary value of the y1-values
!
! Example: x1 = [0 1 2 3]
! y1 = [1 2 1 0]
!
! x2 = -1, y2 = 1
! x2 = 5, y2 = 0
!
!------------------------------------------------------------------------------
integer i1,i2 ! counters
!
real ds ! step size
real fac ! factor in linear interpolation
real s1,s2 ! search values
real xmin1,xmax1 ! minimum and maximum of x1-data
real xmin2,xmax2 ! minimum and maximum of x2-data
!
real, parameter :: eps=1.e-20
!------------------------------------------------------------------------------
! initialisation
!
ierr = 0
!
! check number of points of input array
!
if(n1==1) then
y2 = y1(1)
goto 9999
end if
!
! check minimum and maximum data values
!
xmin1 = minval(x1)
xmax1 = maxval(x1)
xmin2 = minval(x2)
xmax2 = maxval(x2)
!
if (abs(xmin1-xmax1) < eps .or. abs(x1(1)-x1(n1)) < eps) then
ierr = 2
goto 9999
end if
!
if ((abs(xmin2-xmax2) < eps .or. abs(x2(1)-x2(n2)) < eps) .and. n2 > 1) then
ierr = 3
goto 9999
end if
!
! check input data for monotonicity
!
if(x1(1) < x1(n1)) then ! data increasing
do i1=1,n1-1
if(x1(i1) > x1(i1+1)) then
ierr=1
write(*,*) 'z_intp1: i1 x1(i1) x1(i1+1):',i1,x1(i1),x1(i1+1)
goto 9999
end if
end do
!
do i2=1,n2-1
if(x2(i2) > x2(i2+1)) then
ierr=ierr+10
write(*,*) 'z_intp1: i2 x2(i2) x2(i2+1):',i2,x2(i2),x2(i2+1)
goto 9999
end if
end do
!
else ! data decreasing
do i1=1,n1-1
if(x1(i1) < x1(i1+1)) then
ierr=2
write(*,*) 'z_intp1: i1 x1(i1) x1(i1+1):',i1,x1(i1),x1(i1+1)
goto 9999
end if
end do
!
do i2=1,n2-1
if(x2(i2) < x2(i2+1)) then
ierr=ierr + 20
write(*,*) 'z_intp1: i2 x2(i2) x2(i2+1):',i2,x2(i2),x2(i2+1)
goto 9999
end if
end do
end if
!
!------------------------------------------------------------------------------
! initialize
!------------------------------------------------------------------------------
if(ierr==0) then
i1 = 1
s1 = x1(i1)
!
do i2 = 1,n2
s2 = x2(i2)
do while (s1 <= s2 .and. i1 < n1)
i1 = i1 + 1
s1 = x1(i1)
end do
!
! special point
! choose lowest s1-value if x2(:) < x1(1)
!
if(i1 ==1) then
y2(i2) = y1(i1)
else
ds = s2 - x1(i1-1)
fac = ds/(x1(i1)-x1(i1-1))
y2(i2) = y1(i1-1) + fac*(y1(i1)-y1(i1-1))
end if
!
! special case at end: choose s2(n2) > s1(n1), choose last value of y1(1)
!
if(i2==n2 .and. s2>s1) y2(n2) = y1(n1)
end do
end if
!
9999 continue
!
return
end subroutine
!
!-----------------------------------------------------------------------------!
subroutine z_polyarea(xpol,ypol,npol,area)
!-----------------------------------------------------------------------------!
!
! +-------+ ALKYON Hydraulic Consultancy & Research
! | | P.O. Box 248
! | +---+ 8300 AE Emmeloord
! | | +---+ Tel: +31 527 620909
! +---+ | | Fax: +31 527 610020
! +---+ http://www.alkyon.nl
!
! Gerbrant van Vledder
!
! 0. Update history
!
! 0.01 12/06/2003 Initial version
!
! 1. Purpose
!
! Computes area of a closed polygon
!
! 2. Method
!
! The area of the polygon
!
! 3. Parameter list
!
! Name I/O Type Description
!
integer, intent(in) :: npol ! Number of points of polygon
real, intent(in) :: xpol(npol) ! x-coodinates of polygon
real, intent(in) :: ypol(npol) ! y-coordinates of polygon
real, intent(out) :: area ! area of polygon
!
! 4. Subroutines used
!
! 5. Error messages
!
! 6. Remarks
!
integer ipol,ipol1 ! counters
real xmin,xmax,ymin,ymax ! minima and maxima of polygon
real xmean,ymean ! mean values
real xa,ya,xb,yb ! temporary variables
real sumx,sumy ! sums
real darea ! piece of area
!-------------------------------------------------------------------------------
if(npol<=1) then
crf = 0.
xz = 0.
yz = 0.
area = 0.
return
end if
!
! compute minimum and maximum coordinates
!
xmin = minval(xpol)
xmax = maxval(xpol)
ymin = minval(ypol)
ymax = maxval(ypol)
!
! compute mean of range of x- and y-coordinates
!
xmean = 0.5*(xmin + xmax)
ymean = 0.5*(ymin + ymax)
!
! compute area and center of gravity
! do loop over all line pieces of polygon
!
area = 0.
sumx = 0.
sumy = 0.
!
do ipol=1,npol
ipol1 = ipol + 1
if(ipol==npol) ipol1 = 1
xa = xpol(ipol)
ya = ypol(ipol)
xb = xpol(ipol1)
yb = ypol(ipol1)
!
darea = 0.5*((xa-xmean)*(yb-ymean) - (xb-xmean)*(ya-ymean))
area = area + darea
sumx = sumx + darea*(xa+xb+xmean)/3.
sumy = sumy + darea*(ya+yb+ymean)/3.
end do
!
return
end subroutine
!
!-----------------------------------------------------------------------------!
subroutine z_steps(x,dx,nx)
!-----------------------------------------------------------------------------!
!
!
! +-------+ ALKYON Hydraulic Consultancy & Research
! | | Gerbrant van Vledder
! | +---+
! | | +---+ Creation date: September 28, 1998
! +---+ | | Last Update: march 19, 2003
! +---+
!
! 0. Update history
!
! 19/03/2003 Input argument defined using intent option
! check included nx > 0
!
! 1. Purpose
!
! Compute bandwidth of spectral discretization
!
implicit none
!
integer, intent(in) :: nx ! Number of elements in array
real, intent(in) :: x(nx) ! Input data array with elements
real, intent(out) :: dx(nx) ! Output array with step sizes
!
integer ix ! counter
!------------------------------------------------------------------------------
if (nx<1) then
return
!
elseif (nx==1) then
dx = 0
else
do ix=2,nx-1
dx(ix) = 0.5 * (x(ix+1) - x(ix-1))
end do
!
if (nx >= 4) then
dx(1) = dx(2)*dx(2)/dx(3)
dx(nx) = dx(nx-1)*dx(nx-1)/dx(nx-2)
else
dx(1) = dx(2)
dx(nx) = dx(nx-1)
end if
end if
!
return
end subroutine
!-----------------------------------------------------------------------------!
real function z_root2(func,x1,x2,xacc,iprint,ierr)
!-----------------------------------------------------------------------------!
!
! +-------+ ALKYON Hydraulic Consultancy & Research
! | |
! | +---+
! | | +---+
! +---+ | |
! +---+
!
! 0. Update history
!
! Version Date Modification
!
! 0.01 29/11/1999 Initial version
! 0.02 07/11/1999 Test added to check boundaries, and reverse if necessary
! Bug fixed in assigning answer
! 0.03 02/09/2002 Maximum number of iterations set to 20, instead of 10
!
! 1. Purpose
!
! Find zero crossing point of function FUNC between the
! initial values on either side of zero crossing
!
! 2. Method
!
! Ridders method of root finding
!
! adapted from routine zridddr
! Numerical Recipes
! The art if scientific computing, second edition, 1992
! W.H. Press, S.A. Teukolsky, W.T. Vetterling and B.P. Flannery
!
! 3. Parameter list
!
! Name I/O Type Description
!
! func i r real function
! x1 i r initial x-value on left/right side of zero-crossing
! x2 i r initial x-value on right/left side of zero-crossing
! xacc i r accuracy, used as |x1(i)-x2(i)|< xacc
! iprint i i Output channel and test level
! ierr o i Error indicator
!
! 4. Subroutines used
!
! Func user supplied real function
!
! 5. Error messages
!
! ierr = 0 No errors occured during iteration process
! 1 Iteration halted in dead end, this combination may NEVER occur
! 2 Maximum number of iterations exceeded
! 3 Solution jumped outside interval
!
! 6. Remarks
!
! It is assumed that the x1- and x2-coordinate lie
! on different sides of the actual zero crossing
!
! The input parameter IPRINT is used to generate test output.
! If IPRINT==0, no test output is created
! > 0, test output is directed to the file connected to unit LUPRINT=IPRINT
! if no file is connected to this unit, no output is written
!
!
implicit none
!
real func ! external function
real, intent (in) :: x1 ! x-value at one side of interval
real, intent (in) :: x2 ! x-value at other side of interval
real, intent (in) :: xacc ! requested accuracy
integer, intent (in) :: iprint ! number of output channel, only used when
integer, intent (out) :: ierr ! error indicator
!
real unused ! default value
real zriddr ! intermediate function value
real xx1,xx2,xx ! local boundaries during iteration
integer maxit ! maximum number of iteration
integer luprint ! unit of test output
logical lopen ! check if a file is opened
parameter (maxit = 20)
external func
!
integer iter ! counter for number of iterations
real fh ! function value FUNC(xh)
real fl ! function value FUNC(xl)
real fm ! function value FUNC(xm)
real fnew ! function value FUNC(xnew)
real s ! temp. function value, used for inverse quadratic interpolation
real xh ! upper (high) boundary of interval
real xl ! lower boundary of interval
real xm ! middle point of interval
real xnew ! new estimate according to Ridders method
!
ierr = 0 ! set error level
unused =-1.11e30 ! set start value
!
xx1 = x1 ! copy boundaries of interval to local variables
xx2 = x2
!
luprint = iprint
!
if(luprint > 0) then
inquire(unit=luprint,opened=lopen)
if(.not.lopen) then
luprint = 0
write(*,'(a,i4)') 'Z_ROOT2: invalid unit number:',iprint
end if
end if
!
! check boundaries on requirement x2 > x1
!
if(xx1 > xx2) then
xx = xx1
xx1 = xx2
xx2 = xx
end if
!
fl = func(xx1)
fh = func(xx2)
!
!if(luprint > 0) write(luprint,'(a,4e13.5)') &
!& 'Z_ROOT2: xx1 xx2 fl fh:',xx1,xx2,fl,fh
!
if((fl > 0. .and. fh < 0.) .or. (fl < 0. .and. fh > 0.))then
xl = xx1
xh = xx2
zriddr = unused
!
do iter=1,maxit
xm = 0.5*(xl+xh)
fm = func(xm)
s = sqrt(fm**2-fl*fh)
if(s == 0.) goto 9000
xnew = xm+(xm-xl)*(sign(1.,fl-fh)*fm/s)
!
! if(luprint>0) write(luprint,'(a,4e13.5)') &
!& 'Z_ROOT2: xm,fm,s,xnew:',xm,fm,s,xnew
!
if (abs(xnew-zriddr) <= xacc) then
! if(luprint>0) write(luprint,'(a)') 'Z_ROOT2: xnew=zriddr'
goto 9000
end if
!
zriddr = xnew
fnew = func(zriddr)
if (fnew == 0.) goto 9000
!
if(sign(fm,fnew) /= fm) then
xl = xm
fl = fm
xh = zriddr
fh = fnew
elseif(sign(fl,fnew) /= fl) then
xh = zriddr
fh = fnew
elseif(sign(fh,fnew) /= fh) then
xl = zriddr
fl = fnew
else
ierr = 1
goto 9000
endif
!
if(abs(xh-xl) <= xacc) goto 9000
!
if(luprint > 0) write(luprint,'(a,i4,5e14.6)') &
& 'Z_ROOT2: iter,x1,x2,|x1-x2|,xacc,z:', iter,xl,xh,abs(xl-xh),xacc,fnew
!
end do
ierr = 2
if(luprint > 0) write(luprint,'(a)') 'Z_ROOT2: -> ierr=2'
goto 9000
else if (fl == 0.) then
zriddr = xx1
else if (fh == 0.) then
zriddr = xx2
else
ierr = 3
goto 9999
! 'root must be bracketed in zriddr'
endif
!
9000 continue
!
z_root2 = zriddr
!
if(luprint > 0) write(luprint,'(a,2i3,5e13.5)') &
& 'Z_ROOT2: ierr,iter,xl,xh,acc,x0,z0:', ierr,iter,xl,xh,xacc,z_root2,func(z_root2)
!
9999 continue
!
return
end function
!
!-----------------------------------------------------------------------------!
subroutine z_upper(str)
!-----------------------------------------------------------------------------!
!
! +-------+ ALKYON Hydraulic Consultancy & Research
! | | Gerbrant van Vledder
! | +---+
! | | +---+ Creation date: July 3,1998
! +---+ | | Last Update:
! +---+
!
! 0. Update history
!
! 1. Purpose
!
! Transform all lower capitals to UPPER CAPITALS in string STR
!
! 2. Method
!
! 3. Parameter list
!
! Name I/O Type Description
!
implicit none
character(len=*), intent(inout) :: str ! Character string to be converted
!
! 4. Subroutines used
!
! 5. Error messages
!
! 6. Remarks
!
integer nlen
integer i,ial,iau,izl
!
nlen = len(str)
!
ial = ichar('a')
iau = ichar('A')
izl = ichar('z')
!
do i=1,nlen
if(ichar(str(i:i)) >= ial.and. ichar(str(i:i)) <= izl) then
str(i:i) = char(ichar(str(i:i))-ial+iau)
end if
end do
!
return
end subroutine
!
!-----------------------------------------------------------------------------!
real function z_wnumb(w,d,grav)
!-----------------------------------------------------------------------------!
!
! +-------+ ALKYON Hydraulic Consultancy & Research
! | | Gerbrant van Vledder
! | +---+
! | | +---+
! +---+ | |
! +---+
!
implicit none
!
! 0. Update history
!
! 01/04/1999 Initial version
! 12/01/2001 grav added as input parameter
! 11/04/2001 Check included for the case w < 0 or d < 0
!
! 1. Purpose:
!
! Compute wave number k for a given radian frequency and water depth
!
! 2. Method
!
! finite depth linear dispersion relation, using a Pade approximation
!
! 3. Parameter list:
!
!Type I/O Name Description
!------------------------------------------------------------------------------
real, intent(in) :: w ! radian frequency (rad)
real, intent(in) :: d ! water depth (m)
real, intent(in) :: grav ! graviational acceleration (m/s^2)
!
! 4. Subroutines used
!
! 5. Error messages
!
! 6. Remarks
!
! The Pade approximation has been described in Hunt, 198.
!
!
real x,xx,y,omega
!
if(d<=0 .or. w<= 0.) then
z_wnumb = -10.
else
omega = w**2/grav
y = omega*d
xx = y*(y+1./(1.+y*(0.66667+y*(0.35550+y*(0.16084+y*(0.06320+y* &
& (0.02174+y*(0.00654+y*(0.00171+y*(0.00039+y*0.00011))))))))))
x = sqrt(xx)
z_wnumb = x/d
end if
!
return
end function
end module