module viscous_m
  real ::  prk, jrk, bly  !  Prandtl, Rich number, and shear layer depth
  real ::  zdeb=5.        !  Integration bound for ODE integrator
  real ::  eps            !  smallest semi-grid 
  real ::  hstart,tol
contains
function viscous(z,vis,lfirst)
use nr, only: odeint, bsstep,rkqs
implicit none
real :: viscous(12)
real, intent(in):: z
real, save :: start(12)
real, save :: zbef
interface
   function vis(z)
   real, intent(in)::z
   real, dimension(12) :: vis
   end function vis
end interface
logical, intent(inout) :: lfirst
integer, parameter :: ntab=501
real, save :: ztab(ntab)
real, save :: vtab(12,ntab)
real :: dzt,h1,hmin
integer :: ll,kk
integer , save :: ltab
integer, parameter :: neq=12,lenwrk=32*neq,method=2
logical ::   errass, mesage
real :: work(lenwrk) 
!external :: FI

!h1=0.001;hmin=0.;zdeb=5.;tol=1.E-5  !  This always work!!!
h1=0.001;hmin=0.;tol=1.E-6
if(z.ge.zdeb) then ; viscous=vis(z) ; else ;
if(lfirst)then
start=vis(zdeb)
call odeint(start,zdeb,z,tol,h1,hmin,fi,rkqs); viscous=start; lfirst=.false.
else 
if(z.eq.zbef)then;viscous=start;else
call odeint(start,zbef,z,tol,h1,hmin,fi,rkqs); viscous=start; endif
endif
endif
zbef=z
end function viscous

function vis1(z)
implicit none
real, intent(in)::z
real, dimension(12) :: vis1
complex:: cmu
complex:: ci=(0.,1.),vitw,uzz,csqci
csqci=(1.,1.)/sqrt(2.)
if(jrk.ge.0.25) then; cmu=-sqrt(jrk-0.25);else;cmu=-ci*sqrt(0.25-jrk);endif
vitw=sqrt(z)*exp(ci*cmu*log(z))  ! w
if(jrk.lt.0.25)vitw=csqci*sqrt(z)*exp(ci*cmu*log(z))  ! w
vis1( 9)= real(vitw);vis1(10)=aimag(vitw)
vitw=(0.5+ci*cmu)*vitw/z                   !w_z and u
vis1( 3)=-aimag(vitw);vis1(4)=real(vitw)
vitw=(-0.5+ci*cmu)*vitw/z                     !w_zz and chi
vis1( 1)=-aimag(vitw);vis1(2)=real(vitw)
vitw=(-1.5+ci*cmu)*vitw/z                     !w_zzz and p plus bit of xsi
vis1(11)=real(vitw)-z*vis1(3)-vis1(10);vis1(12)=aimag(vitw)-z*vis1(4)+vis1(9)
vis1(5)=z*aimag(vitw);vis1(6)=-z*real(vitw)
vitw=(-2.5+ci*cmu)*vitw/z                     !w_zzzz and b
vis1(7)=real(vitw)-z*vis1(1); vis1(8)=aimag(vitw)-z*vis1(2)
vitw=(-3.5+ci*cmu)*vitw/z                     !w_zzzzz and xsi
vis1(5)=real(vitw)-vis1(1)+vis1(5)
vis1(6)=aimag(vitw)-vis1(2)+vis1(6)
end function vis1
function vis2(z)
implicit none
real, intent(in)::z
real, dimension(12) :: vis2
complex :: cmu
complex :: ci=(0.,1.),vitw,uzz,csqci
csqci=(1.,-1.)/sqrt(2.)
if(jrk.ge.0.25) then; cmu=+sqrt(jrk-0.25);else;cmu=+ci*sqrt(0.25-jrk);endif
vitw=sqrt(z)*exp(ci*cmu*log(z))  ! w
if(jrk.lt.0.25)vitw=csqci*sqrt(z)*exp(ci*cmu*log(z))  ! w
vis2( 9)= real(vitw);vis2(10)=aimag(vitw)
vitw=(0.5+ci*cmu)*vitw/z                   !w_z and u
vis2( 3)=-aimag(vitw);vis2(4)=real(vitw)
vitw=(-0.5+ci*cmu)*vitw/z                     !w_zz and chi
vis2( 1)=-aimag(vitw);vis2(2)=real(vitw)
vitw=(-1.5+ci*cmu)*vitw/z                     !w_zzz and p plus bit of xsi
vis2(11)=real(vitw)-z*vis2(3)-vis2(10);vis2(12)=aimag(vitw)-z*vis2(4)+vis2(9)
vis2(5)=z*aimag(vitw);vis2(6)=-z*real(vitw)
vitw=(-2.5+ci*cmu)*vitw/z                     !w_zzzz and b
vis2(7)=real(vitw)-z*vis2(1); vis2(8)=aimag(vitw)-z*vis2(2)
vitw=(-3.5+ci*cmu)*vitw/z                     !w_zzzzz and xsi
vis2(5)=real(vitw)-vis2(1)+vis2(5)
vis2(6)=aimag(vitw)-vis2(2)+vis2(6)
end function vis2
function vis3(z)
implicit none
real, intent(in)::z
real, dimension(12) :: vis3
complex :: ci=(0.,1.),vitw,csqci
csqci=(1.,1.)/sqrt(2.)
vitw=csqci*exp(-1.25*log(z)-2./3.*csqci*(z**1.5))  ! w
vis3( 9)= real(vitw);vis3(10)=aimag(vitw)
vitw=-csqci*sqrt(z)*vitw             !w_z and u
vis3( 3)=-aimag(vitw);vis3(4)=real(vitw)
vitw=-csqci*sqrt(z)*vitw                      !w_zz and chi
vis3( 1)=-aimag(vitw);vis3(2)=real(vitw)
vitw=-csqci*sqrt(z)*vitw                      !w_zzz and p plus bit of xsi
vis3(11)=real(vitw)-z*vis3(3)-vis3(10);vis3(12)=aimag(vitw)-z*vis3(4)+vis3(9)
vis3(5)=z*aimag(vitw);vis3(6)=-z*real(vitw)
vitw=-csqci*sqrt(z)*vitw                      !w_zzzz and b
vis3(7)=real(vitw)-z*vis3(1); vis3(8)=aimag(vitw)-z*vis3(2)
vitw=-csqci*sqrt(z)*vitw                      !w_zzzzz and xsi
vis3(5)=real(vitw)-vis3(1)+vis3(5)
vis3(6)=aimag(vitw)-vis3(2)+vis3(6)
end function vis3
function vis4(z)
implicit none
real, intent(in)::z
real, dimension(12) :: vis4
complex :: ci=(0.,1.),vitw,csqci
csqci=(1.,1.)/sqrt(2.)*sqrt(prk)
vitw=exp(-2.25*log(z)-2./3.*csqci*(z**1.5))  ! w
vis4( 9)= real(vitw);vis4(10)=aimag(vitw)
vitw=-csqci*sqrt(z)*vitw             !w_z and u
vis4( 3)=-aimag(vitw);vis4(4)=real(vitw)
vitw=-csqci*sqrt(z)*vitw                      !w_zz and chi
vis4( 1)=-aimag(vitw);vis4(2)=real(vitw)
vitw=-csqci*sqrt(z)*vitw                      !w_zzz and p plus bit of xsi
vis4(11)=real(vitw)-z*vis4(3)-vis4(10);vis4(12)=aimag(vitw)-z*vis4(4)+vis4(9)
vis4(5)=z*aimag(vitw);vis4(6)=-z*real(vitw)
vitw=-csqci*sqrt(z)*vitw                      !w_zzzz and b
vis4(7)=real(vitw)-z*vis4(1); vis4(8)=aimag(vitw)-z*vis4(2)
vitw=-csqci*sqrt(z)*vitw                      !w_zzzzz and xsi
vis4(5)=real(vitw)-vis4(1)+vis4(5)
vis4(6)=aimag(vitw)-vis4(2)+vis4(6)
end function vis4
function bab1(xk,z)
implicit none
real, intent(in)::xk
complex, intent(in)::z
real, dimension(12) :: bab1
complex:: cmu,csqci
complex:: ci=(0.,1.),vitw,vitu,buoy,chi,pres
csqci=(1.,1.)/sqrt(2.)
if(abs(xk*real(sqrt(z*conjg(z)))).lt.eps)then; bab1=0.;return; endif
if(jrk.ge.0.25) then; cmu=-sqrt(jrk-0.25);else;cmu=-ci*sqrt(0.25-jrk);endif
vitw=sqrt(z)*exp(ci*cmu*log(z))  ! w
bab1( 9)= real(vitw);bab1(10)=aimag(vitw)
vitu=ci*(0.5+ci*cmu)*vitw/z/xk                 !w_z and u
bab1( 3)=real(vitu);bab1(4)=aimag(vitu)
chi=(-0.5+ci*cmu)*vitu/z                       !w_zz and chi
bab1( 1)=real(chi);bab1(2)=aimag(chi)
pres=-z*vitu+ci/xk*vitw    !Pressure field
bab1(11)=real(pres);bab1(12)=aimag(pres)
buoy=ci/z/xk*jrk*vitw                    !Buoyancy field
bab1(7)=real(buoy);bab1(8)=aimag(buoy)
chi=0.                                       !chi used again for b_z
bab1(5)=real(chi); bab1(6)=aimag(chi)
end function bab1
function bab2(xk,z)
implicit none
real, intent(in)::xk
complex, intent(in)::z
real, dimension(12) :: bab2
complex:: cmu,csqci
complex:: ci=(0.,1.),vitw,vitu,buoy,chi,pres
csqci=(1.,-1.)/sqrt(2.)
if(abs(xk*real(sqrt(z*conjg(z)))).lt.eps)then; bab2=0.;return; endif
if(jrk.ge.0.25) then; cmu=+sqrt(jrk-0.25);else;cmu=+ci*sqrt(0.25-jrk);endif
vitw=sqrt(z)*exp(ci*cmu*log(z))  ! w
bab2( 9)= real(vitw);bab2(10)=aimag(vitw)
vitu=ci*(0.5+ci*cmu)*vitw/z/xk                 !w_z and u
bab2( 3)=real(vitu);bab2(4)=aimag(vitu)
chi=(-0.5+ci*cmu)*vitu/z                       !w_zz and chi
bab2( 1)=real(chi);bab2(2)=aimag(chi)
pres=-z*vitu+ci/xk*vitw    !Pressure field
bab2(11)=real(pres);bab2(12)=aimag(pres)
buoy=ci/z/xk*jrk*vitw                    !Buoyancy field
bab2(7)=real(buoy);bab2(8)=aimag(buoy)
chi=0.                                       !chi used again for b_z
bab2(5)=real(chi); bab2(6)=aimag(chi)
end function bab2

function hanksol(xk,z)
implicit none
real, intent(in)::xk,z
real, dimension(12) :: hanksol
real :: pi=acos(-1.)
complex:: cmu
complex:: ci=(0.,1.),vitw,vitu,buoy,chi,pres,vitp,vitm
!csqci=(1.,1.)/sqrt(2.)
if(abs(xk*z).lt.eps)then; hanksol=0.;return; endif
if(jrk.ge.0.25) then; cmu=+sqrt(jrk-0.25);else;cmu=+ci*sqrt(0.25-jrk);endif
!vitw:
vitw=ci*sqrt(pi*xk*z/2.)*exp(-cmu*pi/2.)*hankel(1,ci*cmu,ci*xk*z,ci*0.)
hanksol( 9)= real(vitw);hanksol(10)=aimag(vitw)
!vitu:
vitp=ci*sqrt(pi*(xk*z+eps)/2.)*exp(-cmu*pi/2.)*hankel(1,ci*cmu,ci*xk*z,ci*eps)
vitm=ci*sqrt(pi*(xk*z-eps)/2.)*exp(-cmu*pi/2.)*hankel(1,ci*cmu,ci*xk*z,-ci*eps)
vitu=ci*(vitp-vitm)/2./eps
hanksol( 3)=real(vitu);hanksol(4)=aimag(vitu)
!chi at at 0
chi=ci*xk*(vitp-2*vitw+vitm)/eps**2;
hanksol(1)=real(chi);hanksol(2)=aimag(chi)
! Pressure:
pres=-z*vitu+ci/xk*vitw  
hanksol(11)=real(pres);hanksol(12)=aimag(pres)
!Buoyancy:
buoy=ci/z/xk*jrk*vitw                    !Buoyancy field
hanksol(7)=real(buoy);hanksol(8)=aimag(buoy)
!b_z (xsi) at 0.
hanksol(5)=0.; hanksol(6)=0.
end function hanksol

function hankel(order,nu,zin,eps)
!  use nr, only: cgamma
! Uses formula 9.1.10 and 9.2.7 from AS tested for z and nu imaginary.
   implicit none
   complex:: cgamma
   complex :: hankel
   integer, intent(in) :: order
   complex, intent(in) :: nu,zin,eps
   complex :: zz
   integer, parameter :: nn=100
   integer :: n
   real*16 :: pi
   complex*16 :: facnu,faczz,bessp,bessm,ci
   complex*16 :: chi,mu
   complex*16 :: facp,facq
   ci=(0.,1._16);pi=acos(-1._16);hankel=0.;zz=zin+eps
!  if(real(sqrt(zz*conjg(zz))).ge.1.e-5)then
   if(real(sqrt(zin*conjg(zin))).le.4.*max(1.,real(nu*conjg(nu))**(0.3)))then
     faczz=1.;facnu=1.
     do n=1,nn
     facnu=-facnu/float(n)/(float(n)+nu)*(zz**2/4);faczz=faczz+facnu
     enddo
     bessp=faczz*(zz/2.)**nu/cgamma(nu+1)
     faczz=1.;facnu=1.
     do n=1,nn
     facnu=-facnu/float(n)/(float(n)-nu)*(zz**2/4);faczz=faczz+facnu
     enddo
     bessm=faczz*(zz/2.)**(-nu)/cgamma(-nu+1)
     if(order.eq.1)hankel=exp(-ci*pi*nu)*bessp-bessm
     if(order.eq.2)hankel=bessm-bessp*exp(+ci*pi*nu)
     if(order.eq.1.or.order.eq.2)then;hankel=hankel*ci/sin(nu*pi)
     else;write(*,*)'order Hankel incorrect';stop;endif
   else
     chi=zz-(nu/2.+0.25)*pi;mu=4.*nu**2
     if(order.eq.1)then;facp=1.;hankel=facp;do n=1,20
        facp=+ci*facp*(mu-(2.*n-1.)**2)/float(n)/8./zz
        hankel=hankel+facp;enddo
     endif
     if(order.eq.2)then;facp=1.;hankel=facp;do n=1,20
        facp=-ci*facp*(mu-(2.*n-1.)**2)/float(n)/8./zz
        hankel=hankel+facp;enddo
     endif
     if(order.eq.1) hankel=hankel*exp(+ci*chi)
     if(order.eq.2) hankel=hankel*exp(-ci*chi)
     if(order.eq.1.or.order.eq.2)then;hankel=hankel*sqrt(2./pi/zz)
     else;write(*,*)'order Hankel incorrect ';stop;endif
   endif
   !endif
   end function hankel

subroutine FI(z,y,yp)
implicit none
real, intent(in)::z
real, dimension(:),intent(in)::y
real, dimension(:), intent(out) :: yp
   yp( 1)=-z*y(4)+y( 9)-y(12)        !chi_z = i(z-zk)u+w+ip
   yp( 2)= z*y(3)+y(10)+y(11)
   yp( 3)= y(1)                            !u_z=chi
   yp( 4)= y(2)
   yp( 5)= prk*(-z*y(8)+jrk*y( 9)) ! xsi_z = Pr(i(z-zk/Pr)b+Jw)
   yp( 6)= prk*(+z*y(7)+jrk*y(10))
   yp( 7)= y(5)                            ! b_z=xsi
   yp( 8)= y(6)
   yp( 9)= y(4)                            !  w_z = -iu
   yp(10)=-y(3)
   yp(11)= y(7)                            ! p_z=b
   yp(12)= y(8)
end subroutine FI

function hypsol(aa,bb,cc,b1,b2,xk,z)
use nr, only: hypgeo
implicit none
real, intent(in)::xk,z
complex, intent(in)::aa,bb,cc,b1,b2
real, dimension(12) :: hypsol
complex :: mu,cm,rr,rm

real :: pi=acos(-1.),ss,cs
complex:: wup(3),wdo(3),wmi(3)
integer :: is
complex:: ci=(0.,1.),vitw,vitu,buoy,chi,pres,vitp,vitm,xsi
!csqci=(1.,1.)/sqrt(2.)
if(z.lt.eps.or.abs(xk).le.eps)then; hypsol=0.;return; endif
if(jrk.ge.0.25) then; mu=+sqrt(jrk-0.25);else;mu=+ci*sqrt(0.25-jrk);endif
cm=cc-aa-bb
ss=tanh(z/bly);rm=ss**2
do is=1,3 
 ss=tanh((z+float(is-2)*eps)/bly) ; rr=ss**2 ; cs=cosh((z+float(is-2)*eps)/bly)
 !if(real(rr).ge.0.5)  wup(is)=2**(-cm)*ss**(0.5+ci*mu)*cs**(-cm)&
 if(real(rm).ge.0.05)  wmi(is)=2**(-cm)*ss**(0.5+ci*mu)*cs**(-cm)&
                           *hypgeo(cc-bb,cc-aa,cc-aa-bb+1.,1.-rr) 
 !if(real(rr).le.0.5)   wdo(is)=cs**(+cm)* &
 if(real(rm).le.0.05)   wmi(is)=cs**(+cm)* &
                     (b1*ss**(0.5-ci*mu)*hypgeo(aa-cc+1,bb-cc+1,2.-cc,rr) &
                     +b2*ss**(0.5+ci*mu)*hypgeo(aa,bb,cc,rr) ) 
! if(real(rr).gt.0.5) wdo(is)=wup(is) ; if(real(rr).lt.0.50) wup(is)=wdo(is)
!if(real(rr).ge.0.5) wmi(is)=wup(is) ; if(real(rr).le.0.5)wmi(is)=wdo(is)
!wmi(is)=((1.+tanh(10.*(rr-0.5)))*wup(is)+(1.-tanh(10.*(rr-0.5)))*wdo(is))/2. 
!wmi(is)=wdo(is) 
enddo

vitw=wmi(2) ; hypsol( 9)= real(vitw);hypsol(10)=aimag(vitw)
!vitu:
vitp=ci/xk*(wmi(3)-wmi(2))/eps ; vitm=ci/xk*(wmi(2)-wmi(1))/eps
vitu=(vitm+vitp)/2.
hypsol( 3)=real(vitu);hypsol(4)=aimag(vitu)
!chi at at 0
chi=(vitp-vitm)/eps; hypsol(1)=real(chi);hypsol(2)=aimag(chi)
! Pressure:
ss=tanh(z/bly) ; pres=-bly*ss*vitu+ci/xk*vitw*(1.-ss**2)
hypsol(11)=real(pres);hypsol(12)=aimag(pres)
!Buoyancy:
buoy=ci/bly/ss/xk*jrk*vitw                    !Buoyancy field
hypsol(7)=real(buoy);hypsol(8)=aimag(buoy)
!b_z (xsi) at 0.
xsi=(-buoy*(1.-ss**2)+jrk*vitu)/ss/bly
hypsol(5)=real(xsi); hypsol(6)=aimag(xsi)
end function hypsol

end module viscous_m



