! v0.2 13 October 2008
!*** Left over from long ago: some comments indicate some code
!*** is "temporary", not sure why.  

subroutine advect5(dqdt,q)
implicit real*8 (a-h,o-z)
include 'eza.size'
include 'eza.par'
include 'eza.common'
dimension dqdt(in,jn),q(in,jn)
common/contra/uc(in,jn),wc(in,jn)

do  j=1,jn
  do  i=2,in-1
    if (uc(i,j).gt.0.) then
      
      if (i.gt.3.and.i.lt.in-1) then
        q1=(-3.*q(i+2,j)+30.*q(i+1,j)+20.*q(i,j) &
          -60.*q(i-1,j)+15.*q(i-2,j)-2.*q(i-3,j))/120.
      else if (i.gt.2.) then
        q1=(2*q(i+1,j)+3*q(i,j)-6*q(i-1,j)+q(i-2,j))/12.
      else
        !          q1=.5*(q(i,j)-q(i-1,j))
        q1=.25*(q(i+1,j)-q(i-1,j))
      end if
      
    else
      if (i.gt.2.and.i.lt.in-2) then
        q1=(2.*q(i+3,j)-15.*q(i+2,j)+60.*q(i+1,j) &
          -20.*q(i,j)-30.*q(i-1,j)+3.*q(i-2,j))/120.
      else if (i.lt.in-1) then
        q1=(-q(i+2,j)+6*q(i+1,j)-3*q(i,j)-2*q(i-1,j))/12.
      else
        !          q1=.5*(q(i+1,j)-q(i,j))
        q1=.25*(q(i+1,j)-q(i-1,j))
      end if
    end if
    
    dqdt(i,j)=dqdt(i,j)-uc(i,j)*q1
  enddo
enddo


do j=2,jn-1
  do i=1,in
    if (wc(i,j).gt.0.) then
      if (j.gt.3.and.j.lt.jn-2) then
        q2=(-3.*q(i,j+2)+30.*q(i,j+1)+20.*q(i,j) &
          -60.*q(i,j-1)+15.*q(i,j-2) &
          -2.*q(i,j-3))/120.
      else if (j.gt.2.) then
        q2=(2*q(i,j+1)+3*q(i,j)-6*q(i,j-1)+q(i,j-2))/12.
      else
        !          q2=.5*(q(i,j)-q(i,j-1))
        q2=.25*(q(i,j+1)-q(i,j-1))
      end if
    else
      if (j.gt.2.and.j.lt.jn-2) then
        q2=(2.*q(i,j+3)-15.*q(i,j+2) &
          +60.*q(i,j+1)-20.*q(i,j)-30.*q(i,j-1) &
          +3.*q(i,j-2))/120.
      else if (j.lt.jn-1) then
        q2=(-q(i,j+2)+6*q(i,j+1)-3*q(i,j)-2*q(i,j-1))/12.
      else
        !          q2=.5*(q(i,j+1)-q(i,j))
        q2=.25*(q(i,j+1)-q(i,j-1))
      end if
    end if
    dqdt(i,j)=dqdt(i,j)-wc(i,j)*q2
  enddo
enddo

return
end

subroutine rterm(dqdt,q)
implicit real*8 (a-h,o-z)
include 'eza.size'
include 'eza.par'
include 'eza.common'
dimension dqdt(in,jn),q(in,jn)
do j=1,jn
  dqdt(1,j)=0.
  do i=2,in
    dqdt(i,j)=dqdt(i,j)-viscu(i,j)*q(i,j)/(xa(i)*xa(i))
  end do
end do
return
end

subroutine advect3(dqdt,q)
implicit real*8 (a-h,o-z)
include 'eza.size'
include 'eza.par'
include 'eza.common'
dimension dqdt(in,jn),q(in,jn)
common/contra/uc(in,jn),wc(in,jn)

do  j=1,jn
  do  i=2,in-1
    if (uc(i,j).gt.0.and.i.ne.2) then
      q1=(2*q(i+1,j)+3*q(i,j)-6*q(i-1,j)+q(i-2,j))/12.
    else if (uc(i,j).lt.0.and.i.ne.im) then
      q1=(-q(i+2,j)+6*q(i+1,j)-3*q(i,j)-2*q(i-1,j))/12.
    else
      q1=.25*(q(i+1,j)-q(i-1,j))
    end if
    dqdt(i,j)=dqdt(i,j)-uc(i,j)*q1
  enddo
enddo

do  j=2,jn-1
  do  i=1,in
    if (wc(i,j).gt.0.and.j.ne.2) then
      q2=(2*q(i,j+1)+3*q(i,j)-6*q(i,j-1)+q(i,j-2))/12.
    else if (wc(i,j).lt.0.and.j.ne.jm) then
      q2=(-q(i,j+2)+6*q(i,j+1)-3*q(i,j)-2*q(i,j-1))/12.
    else
      q2=.25*(q(i,j+1)-q(i,j-1))
    end if
    dqdt(i,j)=dqdt(i,j)-wc(i,j)*q2
  enddo
enddo

return
end

subroutine maketrans
implicit real*8 (a-h,o-z)
include 'eza.size'
include 'eza.par'
include 'eza.common'

do i=1,in
  iru(i)=i
  ilu(i)=i-1
  if(i.eq.in) iru(i)=ilu(i)
  if(i.eq.1) ilu(i)=iru(i)
enddo

do j=1,jn
  jtu(j)=j
  jbu(j)=j-1
  if(j.eq.jn) jtu(j)=jbu(j)
  if(j.eq.1) jbu(j)=jtu(j)
enddo

!  do i=1,in
!  xa(i)=(i-1.)*rwid/(in-1.)
!  end do

!  do j=1,jn
!  za(j)=(j-1.)*hite/(jn-1.)
!  end do

do j=1,jn
  gz=(j-1.)/jm
  za(j)=hite*(stz*gz+(2.-stz)*gz**3)/(1.+gz**2)
end do

do i=1,in
  gx=(i-1.)/im
  xa(i)=rwid*(str*gx+(2.-str)*gx**3)/(1.+gx**2)
end do

if (axi.eq.0.) then
  
  do i=1,in
    xb(i)=1.
  end do
  do i=1,im
    xc(i)=1.
  end do
  
else
  
  do i=2,im
    xb(i)=.5*(xa(i-1)+xa(i+1))
  end do
  xb(1)=.5*(xa(1)+xa(2))
  xb(in)=.5*(xa(im)+xa(in))
  do i=1,im
    xc(i)=.5*(xa(i)+xa(i+1))
  end do
  
end if

do j=1,jn
  rha(j)=1.
end do
do j=1,jn
  rhb(j)=.5*(rha(j-1)+rha(j+1))
end do
rhb(1)=.5*(rha(1)+rha(2))
rhb(jn)=.5*(rha(jn)+rha(jm))
do j=1,jm
  rhc(j)=.5*(rha(j)+rha(j+1))
end do


! u-grid:
do i=2,in-1
  a1(i)=1./(.25*(xa(i+1)-xa(i-1)))
end do
a1(1)=1./(.5*(xa(2)-xa(1)))
a1(in)=1./(.5*(xa(in)-xa(in-1)))

do j=2,jn-1
  a2(j)=1./(.25*(za(j+1)-za(j-1)))
end do
a2(1)=1./(.5*(za(2)-za(1)))
a2(jn)=1./(.5*(za(jn)-za(jn-1)))

! p-grid:
do i=1,im
  a1p(i)=1./(.5*(xa(i+1)-xa(i)))
end do

do j=1,jm
  a2p(j)=1./(.5*(za(j+1)-za(j)))
end do

return
end

subroutine velocity
implicit real*8 (a-h,o-z)
include 'eza.size'
include 'eza.par'
include 'eza.common'

do  j=1,jn
  do i=2,in-1
    wf(i,j)=.25*a1(i)*(stm(i+1,j)-stm(i-1,j))
  end do
  wf(1,j)=.5*a1(1)*(stm(2,j)-stm(1,j))
  wf(in,j)=.5*a1(in)*(stm(in,j)-stm(im,j))*slips
end do

do i=1,in
  do j=2,jn-1
    uf(i,j)=-.25*a2(j)*(stm(i,j+1)-stm(i,j-1))
  end do
  uf(i,1)=-.5*a2(1)*(stm(i,2)-stm(i,1))*slipb
  uf(i,jn)=-.5*a2(jn)*(stm(i,jn)-stm(i,jm))*slipt
end do

if (axi.ne.0.) then
  do j=1,jn
    do i=1,in
      w(i,j)=wf(i,j)/(xb(i)*rha(j))
      u(i,j)=uf(i,j)/(xb(i)*rha(j))
    end do
  end do
else
  do j=1,jn
    do i=1,in
      w(i,j)=wf(i,j)/rha(j)
      u(i,j)=uf(i,j)/rha(j)
    end do
  end do
  
end if

return
end

subroutine hop(k,na,nb,nc,dt,dbdt,b,in,jn)
!*** adams-bashforth  forward time-step
implicit real*8 (a-h,o-z)
dimension dbdt(in,jn,3),b(in,jn)
common/adams/ab(3,3)
do  j=1,jn
  do  i=1,in
    b(i,j)=b(i,j)+dt*(ab(1,k)*dbdt(i,j,na)+ab(2,k)*dbdt(i,j,nb)+ &
      ab(3,k)*dbdt(i,j,nc))
  enddo
enddo
return
end

subroutine strat(dbdt)
implicit real*8 (a-h,o-z)
include 'eza.size'
include 'eza.par'
include 'eza.common'
dimension dbdt(in,jn)
do i=1,in
  do j=2,jn-1
    dbdt(i,j)=dbdt(i,j)+w(i,j)/3.1415
  end do
end do
return
end

subroutine torque(dvordt)
implicit real*8 (a-h,o-z)
include 'eza.size'
include 'eza.par'
include 'eza.common'
dimension dvordt(in,jn)

coef=1.
if (ray.ne.0.) coef=3.141592*ray

do  i=2,in-1
  do  j=1,jn
    dvordt(i,j)=dvordt(i,j) &
      +coef*0.25*a1(i)*(b(i+1,j)-b(i-1,j))
  enddo
enddo
return
end


subroutine diffuse(dqdt,q)
implicit real*8 (a-h,o-z)
include 'eza.size'
include 'eza.par'
include 'eza.common'
dimension dqdt(in,jn),q(in,jn)
common/spare/f1(im,jm),f3(im,jm)

do i=1,im
  ilc=i
  irc=i+1
  do j=1,jm
    jbc=j
    jtc=j+1
    f1(i,j)=xc(i)*rhc(j)*viscp(i,j)*a1p(i)*0.25 &
      *(q(irc,jtc)-q(ilc,jtc)+q(irc,jbc)-q(ilc,jbc))
    f3(i,j)=rhc(j)*viscp(i,j)*a2p(j)*0.25 &
      *(q(irc,jtc)-q(irc,jbc)+q(ilc,jtc)-q(ilc,jbc))
  enddo
enddo

do i=2,im
  do j=1,jn
    dqdt(i,j)=dqdt(i,j)+a1(i)*0.25*( &
      +(f1(iru(i),jtu(j))-f1(ilu(i),jtu(j))) &
      +(f1(iru(i),jbu(j))-f1(ilu(i),jbu(j))))/(xb(i)*rhb(j))
  enddo
enddo

do i=1,in
  do j=2,jm
    dqdt(i,j)=dqdt(i,j)+a2(j)*0.25*( &
      +(f3(iru(i),jtu(j))-f3(iru(i),jbu(j))) &
      +(f3(ilu(i),jtu(j))-f3(ilu(i),jbu(j))))/rhb(j)
  enddo
enddo


do i=1,in
  dqdt(i,1)=dqdt(i,1)+ &
    a2(1)*0.5*(f3(iru(i),jtu(1))+f3(ilu(i),jtu(1)))/rhb(1)
  dqdt(i,jn)=dqdt(i,jn)- &
    a2(jn)*0.5*(f3(iru(i),jbu(jn))+f3(ilu(i),jbu(jn)))/rhb(jn)
end do

do j=1,jn
  dqdt(1,j)=dqdt(1,j)+ &
    a1(1)*0.5*(f1(iru(1),jtu(j))+f1(iru(1),jbu(j)))/ &
    (xb(1)*rhb(j))
  dqdt(in,j)=dqdt(in,j)- &
    a1(in)*0.5*(f1(ilu(in),jtu(j))+f1(ilu(in),jbu(j)))/ &
    (xb(in)*rhb(j))
end do

return
end


subroutine stuff(q,in,jn,c)
implicit real*8 (a-h,o-z)
dimension q(in,jn)
do  j=1,jn
  do  i=1,in
    q(i,j)=c
  enddo
enddo
return
end


subroutine prepare
implicit real*8 (a-h,o-z)
include 'eza.size'
include 'eza.par'
include 'eza.common'
real*4  am(in2),bm(in2),cm(in2),an(jn2),bn(jn2),cn(jn2)
real*4 work(nwork),y(in2,jn2)
common/blks/am,bm,cm,an,bn,cn,work,y

do  j=2,jn-1
  cn(j-1)=.25*a2p(j)*a2(j)
  an(j-1)=.25*a2p(j-1)*a2(j)
  bn(j-1)=-an(j-1)-cn(j-1)
enddo

do  i=2,in-1
  cm(i-1)=.25*a1p(i)*a1(i)
  am(i-1)=.25*a1p(i-1)*a1(i)
  bm(i-1)=-am(i-1)-cm(i-1)
enddo

do  j=2,jn-1
  drhdz=.25*a2(j)*(rha(j+1)-rha(j-1))/rha(j)
  cn(j-1)=cn(j-1)-.25*a2(j)*drhdz
  an(j-1)=an(j-1)+.25*a2(j)*drhdz
end do


if (axi.ne.0.) then
  do  i=2,in-1
    cm(i-1)=cm(i-1)-.25*a1(i)/xa(i)
    am(i-1)=am(i-1)+.25*a1(i)/xa(i)
  end do
end if

if (alf.eq.0) then
  call blktri(0,1,jn2,an,bn,cn, &
    1,in2,am,bm,cm,in2,y,ierr,work)
  write (*,*) ' declared nwork=',nwork,' needed nwork=',work(1), &
    ' blktri ierr=',ierr
endif

return
end

subroutine poiblk(v,s,xa,rha)
implicit real*8 (a-h,o-z)
include 'eza.size'
include 'eza.par'
dimension v(in,jn),s(in,jn),xa(in),rha(in)
real*4  am(in2),bm(in2),cm(in2),an(jn2),bn(jn2),cn(jn2)
real*4 work(nwork),y(in2,jn2)
common/blks/am,bm,cm,an,bn,cn,work,y


if (axi.ne.0.) then
  do j=1,jn2
    do i=1,in2
      y(i,j)=v(i+1,j+1)*xa(i+1)*rha(j+1)
    end do
  end do
else
  do j=1,jn2
    do i=1,in2
      y(i,j)=v(i+1,j+1)*rha(j+1)
    end do
  end do
end if

do j=1,jn2
  y(1,j)=y(1,j)-am(1)*s(1,j+1)
  y(in2,j)=y(in2,j)-cm(in2)*s(in,j+1)
end do

do i=1,in2
  y(i,1)=y(i,1)-an(1)*s(i+1,1)
  y(i,jn2)=y(i,jn2)-cn(jn2)*s(i+1,jn)
end do

call blktri(1,1,jn2,an,bn,cn, &
  1,in2,am,bm,cm,in2,y,ierr,work)

do j=1,jn2
  do i=1,in2
    s(i+1,j+1)=y(i,j)
  end do
end do

return
end


subroutine change(iread)
implicit real*8 (a-h,o-z)
include 'eza.size'
include 'eza.par'
include 'eza.common'
open (unit=62,file=trim(path)//'/eza.input',status='unknown')
if (iread.eq.1) read (62,f)
321   write (*,*) ' current parameters are:'
write (*,f)
write (*,*)' do you want to change any parameters?'
write (*,*) " for 'yes, and rewrite eza.input': 2 , 'yes, but do &
&not rewrite eza.input:' 1, 'no:' 0"
read (*,*) ichn
if (ichn.ne.0) then
  write (*,*) &
    ' enter changes with: (space)$f example=1.0, etc=2. $end'
  read(*,f)
  rewind 62
  if (ichn.eq.2) write(62,f)
  goto 321
end if
close (unit=62)

return
end

function diffmax
implicit real*8 (a-h,o-z)
include 'eza.size'
include 'eza.par'
include 'eza.common'
common/spare/coefs(in,jn)
do i=1,in
  do j=1,jn
    coefs(i,j)=viscu(i,j)*.25*(a1(i)*a1(i)+a2(j)*a2(j))
  end do
end do
diffmax=rm(coefs)
return
end

subroutine abcoef(p,q)
implicit real*8 (a-h,o-z)
common/adams/ab(3,3)
ab(3,1)=0.
ab(2,1)=0.
ab(1,1)=1.
ab(3,2)=0.
ab(2,2)=-0.5/p
ab(1,2)=1-ab(2,2)
ab(3,3)=(2+3*p)/(6*q*(q-p))
ab(2,3)=-(1+2*q*ab(3,3))/(2*p)
ab(1,3)=1.-ab(2,3)-ab(3,3)
return
end

function rm(t)
include 'eza.size'
implicit real*8 (a-h,o-z)
dimension t(in,jn),h(jn)
do j=1,jn
  h(j)=t(1,j)
end do
do i=2,in
  do j=1,jn
    h(j)=dmax1(h(j),dabs(t(i,j)))
  enddo
enddo
rm=h(1)
do j=2,jn
  rm=dmax1(rm,h(j))
enddo
return
end

subroutine utop(uu,pp)
implicit real*8 (a-h,o-z)
include 'eza.size'
dimension pp(im,jm),uu(in,jn)
do i=1,im
  ilc=i
  irc=i+1
  do j=1,jm
    jbc=j
    jtc=j+1
    pp(i,j)=0.25*(uu(ilc,jbc)+uu(ilc,jtc) &
      +uu(irc,jtc)+uu(irc,jbc))
  enddo
enddo
return
end

subroutine ptou(pp,uu)
implicit real*8 (a-h,o-z)
include 'eza.size'
include 'eza.common'
dimension pp(im,jm),uu(in,jn)
do  i=1,in
  do  j=1,jn
    uu(i,j)=0.25*(pp(ilu(i),jbu(j))+pp(ilu(i),jtu(j)) &
      +pp(iru(i),jtu(j))+pp(iru(i),jbu(j)))
    
  enddo
enddo
return
end


subroutine adtime(rmuc,rmwc)
implicit real*8 (a-h,o-z)
include 'eza.size'
include 'eza.par'
include 'eza.common'
common/contra/uc(in,jn),wc(in,jn)
do  j=1,jn
  do  i=1,in
    uc(i,j)=a1(i)*u(i,j)
    wc(i,j)=a2(j)*w(i,j)
  enddo
enddo
rmuc=rm(uc)
rmwc=rm(wc)
return
end

subroutine marat(dqdt,q)
implicit real*8 (a-h,o-z)
include 'eza.size'
include 'eza.par'
include 'eza.common'
dimension dqdt(in,jn),q(in,jn)
common/spare/f1(im,jn),f2(in,jm),uu(im,jn),ww(in,jm)

do i=1,im
  do j=1,jn
    uu(i,j)=.5*(uf(i,j)+uf(i+1,j))
    f1(i,j)=uu(i,j)*.5*(q(i,j)+q(i+1,j))
  end do
end do

do j=1,jm
  do i=1,in
    ww(i,j)=.5*(wf(i,j)+wf(i,j+1))
    f2(i,j)=ww(i,j)*.5*(q(i,j)+q(i,j+1))
  end do
end do

if (order.lt.-2.) then
  
  do j=1,jn
    do i=1,im
      if (uu(i,j).ge.0.0.and.i.ne.1) then
        f1(i,j)=uu(i,j)*(2*q(i+1,j)+5*q(i,j)-q(i-1,j))/6.
      else if (uu(i,j).le.0.0.and.i.ne.im) then
        f1(i,j)=uu(i,j)*(-q(i+2,j)+5*q(i+1,j)+2*q(i,j))/6.
      end if
    end do
  end do
  
  do i=1,in
    do j=1,jm
      if (ww(i,j).ge.0.0.and.j.ne.1) then
        f2(i,j)=ww(i,j)*(2*q(i,j+1)+5*q(i,j)-q(i,j-1))/6.
      else if (ww(i,j).le.0.0.and.j.ne.jm) then
        f2(i,j)=ww(i,j)*(-q(i,j+2)+5*q(i,j+1)+2*q(i,j))/6.
      end if
    end do
  end do
end if

do j=1,jn
  do i=2,im
    dqdt(i,j)=dqdt(i,j)+.5*a1(i)*(f1(i-1,j)-f1(i,j))/(xb(i)*rhb(j))
  end do
  dqdt(1,j)=dqdt(1,j)-a1(1)*f1(1,j)/(xb(1)*rhb(j))
  dqdt(in,j)=dqdt(in,j)+a1(in)*f1(im,j)/(xb(in)*rhb(j))
end do

do i=1,in
  do j=2,jm
    dqdt(i,j)=dqdt(i,j)+.5*a2(j)*(f2(i,j-1)-f2(i,j))/(xb(i)*rhb(j))
  end do
  dqdt(i,1)=dqdt(i,1)-a2(1)*f2(i,1)/(xb(i)*rhb(1))
  dqdt(i,jn)=dqdt(i,jn)+a2(jn)*f2(i,jm)/(xb(i)*rhb(jn))
end do

return
end

subroutine axi1(etar,vr)
implicit real*8 (a-h,o-z)
include 'eza.size'
include 'eza.par'
include 'eza.common'
dimension vr(in,jn),etar(in,jn)
if (axi.ne.0.) then
  do j=1,jn
    do i=2,in
      etar(i,j)=vor(i,j)/xa(i)
      vr(i,j)=v(i,j)*xa(i)
    end do
    etar(1,j)=etar(2,j)
    vr(1,j)=0.
  end do
else
  do j=1,jn
    do i=1,in
      etar(i,j)=vor(i,j)
      vr(i,j)=v(i,j)
    end do
  end do
end if
return
end

subroutine axi2(dvordt,dvdt)
implicit real*8 (a-h,o-z)
include 'eza.size'
include 'eza.par'
include 'eza.common'
dimension dvordt(in,jn),dvdt(in,jn)
do j=1,jn
  dvordt(1,j)=0.
  dvdt(1,j)=0.
  do i=2,in
    dvordt(i,j)=dvordt(i,j)*xa(i)
    dvdt(i,j)=dvdt(i,j)/xa(i)
  end do
end do
return
end


subroutine arakawa(dqdt,q)
implicit real*8 (a-h,o-z)
include 'eza.size'
include 'eza.par'
include 'eza.common'
dimension dqdt(in,jn),q(in,jn)
common/spare/f1(im,jn),f2(in,jm),g1(im,jm),g2(im,jm)

do i=1,im
  do j=2,jm
    f1(i,j)=-(stm(i,j+1)-stm(i,j-1)+stm(i+1,j+1)-stm(i+1,j-1))* &
      (q(i+1,j)+q(i,j))
  end do
  f1(i,1)=-2*(stm(i,2 )+stm(i+1,2 ))*(q(i+1,1 )+q(i,1 ))
  f1(i,jn)= 2*(stm(i,jm)+stm(i+1,jm))*(q(i+1,jn)+q(i,jn))
end do


do j=1,jm
  do i=2,im
    f2(i,j)=(stm(i+1,j)-stm(i-1,j)+stm(i+1,j+1)-stm(i-1,j+1))* &
      (q(i,j+1)+q(i,j))
  end do
  f2(1,j)=2*(stm(2,j)+stm(2,j+1))*(q(1,j+1)+q(1,j))
  f2(in,j)=-2*(stm(im,j)+stm(im,j+1))*(q(in,j+1)+q(in,j))
end do


do i=1,im
  do j=1,jm
    g1(i,j)=(stm(i+1,j)-stm(i,j+1))*(q(i,j)+q(i+1,j+1))
    g2(i,j)=(stm(i+1,j+1)-stm(i,j))*(q(i+1,j+1)+q(i,j))
  end do
end do

do i=2,im
  do j=2,jm
    dqdt(i,j)=dqdt(i,j)+(.25/12.)*a1(i)*a2(j)*( &
      f1(i-1,j)-f1(i,j) &
      +f2(i,j-1)-f2(i,j) &
      +g1(i-1,j-1)-g1(i,j) &
      +g2(i,j-1)-g2(i-1,j))
  end do
end do

do i=2,in-1
  j=1
  dqdt(i,j)=dqdt(i,j)+(.25/12.)*a1(i)*a2(j)*( &
    f1(i-1,j)-f1(i,j) &
    -2*f2(i,j) &
    -2*g1(i,j) &
    -2*g2(i-1,j))
  j=jn
  dqdt(i,j)=dqdt(i,j)+(.25/12.)*a1(i)*a2(j)*( &
    f1(i-1,j)-f1(i,j) &
    +2*f2(i,j-1) &
    +2*g1(i-1,j-1) &
    +2*g2(i,j-1))
end do

do j=2,jm
  i=1
  dqdt(i,j)=dqdt(i,j)+(.25/12.)*a1(i)*a2(j)*( &
    -2*f1(i,j) &
    +f2(i,j-1)-f2(i,j) &
    -2*g1(i,j) &
    +2*g2(i,j-1))
  i=in
  dqdt(i,j)=dqdt(i,j)+(.25/12.)*a1(i)*a2(j)*( &
    2*f1(i-1,j) &
    +f2(i,j-1)-f2(i,j) &
    +2*g1(i-1,j-1) &
    -2*g2(i-1,j))
end do

i=1
j=1
dqdt(i,j)=dqdt(i,j)+(.25/12.)*a1(i)*a2(j)*( &
  -2*f1(i,j) &
  -2*f2(i,j) &
  -2*g1(i,j))

i=1
j=jn
dqdt(i,j)=dqdt(i,j)+(.25/12.)*a1(i)*a2(j)*( &
  -2*f1(i,j) &
  +2*f2(i,j-1) &
  +2*g2(i,j-1))

i=in
j=1
dqdt(i,j)=dqdt(i,j)+(.25/12.)*a1(i)*a2(j)*( &
  2*f1(i-1,j) &
  -2*f2(i,j) &
  -2*g2(i-1,j))

i=in
j=jn
dqdt(i,j)=dqdt(i,j)+(.25/12.)*a1(i)*a2(j)*( &
  2*f1(i-1,j) &
  +2*f2(i,j-1) &
  +2*g1(i-1,j-1))

do i=1,in
  do j=1,jn
    dqdt(i,j)=dqdt(i,j)/(rhb(j)*xb(i))
  end do
end do

return
end

subroutine iterate(v,s,xa,rha)
implicit real*8 (a-h,o-z)
include 'eza.size'
include 'eza.par'
dimension v(in,jn),s(in,jn),xa(in),rha(in)
real*4  am(in2),bm(in2),cm(in2),an(jn2),bn(jn2),cn(jn2)
real*4 work(nwork),y(in2,jn2)
common/blks/am,bm,cm,an,bn,cn,work,y
common/spare/r(in,jn),d(in,jn)

if (axi.ne.0.) then
  do j=2,jm
    do i=2,im
      r(i,j)=v(i,j)*xa(i)*rha(j)
    end do
  end do
else
  do j=2,jm
    do i=2,im
      r(i,j)=v(i,j)*rha(j)
    end do
  end do
end if

call stuff(d,in,jn,0.d0)

do piter=1.d0,pit
  
  do j=2,jm
    do i=2,im
      d(i,j)=(r(i,j)-am(i-1)*s(i-1,j)-cm(i-1)*s(i+1,j) &
        -an(j-1)*s(i,j-1) &
        -cn(j-1)*s(i,j+1))/(bm(i-1)+bn(j-1))
    end do
  end do
  
  do j=2,jm
    do i=2,im
      s(i,j)=(r(i,j)-am(i-1)*d(i-1,j)-cm(i-1)*d(i+1,j) &
        -an(j-1)*d(i,j-1) &
        -cn(j-1)*d(i,j+1))/(bm(i-1)+bn(j-1))
    end do
  end do
  
end do

return
end


subroutine coriolis(dvordt,dvdt)
implicit real*8 (a-h,o-z)
include 'eza.size'
include 'eza.par'
include 'eza.common'

dimension dvordt(in,jn),dvdt(in,jn)

if (axi.gt.0.) then
  do i=2,im
    do j=1,jn
      !  dvdt(i,j)=dvdt(i,j)-u(i,j)*v(i,j)/xa(i)-swirl*u(i,j)
      dvdt(i,j)=dvdt(i,j)-swirl*u(i,j)
    end do
    do j=2,jm
      dvordt(i,j)=dvordt(i,j)-0.25*a2(j)*((v(i,j+1)*v(i,j+1)- &
        v(i,j-1)*v(i,j-1))/xa(i)+swirl*(v(i,j+1)-v(i,j-1)))
    end do
  end do
else
  do i=2,im
    do j=1,jn
      dvdt(i,j)=dvdt(i,j)-swirl*u(i,j)
    end do
    do j=2,jm
      dvordt(i,j)=dvordt(i,j) &
        -0.25*a2(j)*(swirl*(v(i,j+1)-v(i,j-1)))
    end do
  end do
end if
return
end

subroutine afdiff(dqdt,qq)
implicit real*8 (a-h,o-z)
include 'eza.size'
include 'eza.par'
include 'eza.common'
dimension dqdt(in,jn),qq(in,jn)
common/spare/f1(im,jm),f3(im,jm),q(in,jn),drdt(in,jn), &
  x3c(im),x3b(in)

do j=1,jn
  do i=2,in
    q(i,j)=qq(i,j)/xa(i)
    drdt(i,j)=0.
  end do
  q(1,j)=q(2,j)
end do
do i=1,im
  x3c(i)=xc(i)**3
end do
do i=1,in
  x3b(i)=xb(i)**3
end do

do i=1,im
  ilc=i
  irc=i+1
  do j=1,jm
    jbc=j
    jtc=j+1
    f1(i,j)=x3c(i)*rhc(j)*viscp(i,j)*a1p(i)*0.25 &
      *(q(irc,jtc)-q(ilc,jtc)+q(irc,jbc)-q(ilc,jbc))
    f3(i,j)=rhc(j)*viscp(i,j)*a2p(j)*0.25 &
      *(q(irc,jtc)-q(irc,jbc)+q(ilc,jtc)-q(ilc,jbc))
  enddo
enddo

do i=2,im
  do j=1,jn
    drdt(i,j)=a1(i)*0.25*( &
      +(f1(iru(i),jtu(j))-f1(ilu(i),jtu(j))) &
      +(f1(iru(i),jbu(j))-f1(ilu(i),jbu(j))))/(x3b(i)*rhb(j))
  enddo
enddo

do i=1,in
  do j=2,jm
    drdt(i,j)=drdt(i,j)+a2(j)*0.25*( &
      +(f3(iru(i),jtu(j))-f3(iru(i),jbu(j))) &
      +(f3(ilu(i),jtu(j))-f3(ilu(i),jbu(j))))/rhb(j)
  enddo
enddo


do i=1,in
  drdt(i,1)=drdt(i,1) &
    +a2(1)*0.5*(f3(iru(i),jtu(1))+f3(ilu(i),jtu(1)))/rhb(1)
  drdt(i,jn)=drdt(i,jn) &
    -a2(jn)*0.5*(f3(iru(i),jbu(jn))+f3(ilu(i),jbu(jn)))/rhb(jn)
end do

do j=1,jn
  drdt(1,j)=drdt(1,j) &
    +a1(1)*0.5*(f1(iru(1),jtu(j))+f1(iru(1),jbu(j)))/ &
    (x3b(1)*rhb(j))
  drdt(in,j)=drdt(in,j) &
    -a1(in)*0.5*(f1(ilu(in),jtu(j))+f1(ilu(in),jbu(j)))/ &
    (x3b(in)*rhb(j))
end do

do j=1,jn
  do i=2,in
    dqdt(i,j)=dqdt(i,j)+drdt(i,j)*xb(i)
  end do
end do

return
end

subroutine  volume(q,volint,avg)
implicit real*8 (a-h,o-z)
include 'eza.size'
include 'eza.par'
include 'eza.common'
dimension q(in,jn)
sumtt=0.
sumvv=0.
do i=1,in
  dy=.5*(za(2)-za(1))*rhb(1)
  sumt=dy*q(i,1)
  sumv=dy
  dy=.5*(za(jn)-za(jm))*rhb(jn)
  sumt=sumt+dy*q(i,jn)
  sumv=sumv+dy
  do j=2,jn-1
    dy=.5*(za(j+1)-za(j-1))*rhb(j)
    sumt=sumt+dy*q(i,j)
    sumv=sumv+dy
  end do
  if (i.eq.1) then
    dx=.5*(xa(2)-xa(1))*xb(i)
  else if (i.eq.in) then
    dx=.5*(xa(in)-xa(im))*xb(i)
  else
    dx=.5*(xa(i+1)-xa(i-1))*xb(i)
  end if
  sumtt=sumt*dx+sumtt
  sumvv=sumv*dx+sumvv
end do
volint=sumtt
avg=sumtt/sumvv
return
end


subroutine  boundary
implicit real*8 (a-h,o-z)
include 'eza.size'
include 'eza.par'
include 'eza.common'
if (slipb.eq.0.) then
  do i=2,in
    vor(i,1)=.125*a2(1)*a2(1)*(8.*stm(i,2)-stm(i,3)-7*stm(i,1)) &
      *.5/xa(i)+.5*vor(i,1)
    v(i,1)=0.
  end do
end if

if (slipt.eq.0.) then
  do i=2,in
    vor(i,jn)=.125*a2(jn)*a2(jn)* &
      (8.*stm(i,jn-1)-stm(i,jn-2)-7*stm(i,jn)) &
      *.5/xa(i)+.5*vor(i,jn)
    v(i,jn)=0.
  end do
end if

if (slips.eq.0.) then
  do j=2,jm
    vor(in,j)=.125*a1(in)*a1(in)* &
      (8*stm(in-1,j)-stm(in-2,j)-7*stm(in,j)) &
      *.5/xa(in)+.5*vor(in,j)
    v(in,j)=0.
  end do
end if

! temporary???
if (order.eq.0.and.alf.eq.0.and.slipt.eq.1) then
  do i=1,in
    v(i,jn)=v(i,jn-1)
  end do
end if

return
end

subroutine incmp(dudt,dwdt,divm)
implicit real*8 (a-h,o-z)
include 'eza.size'
include 'eza.par'
include 'eza.common'
dimension dudt(in,jn),dwdt(in,jn),dqdt(in,jn),drdt(in,jn)
common/spare/ur(in,jn),wr(in,jn),alpha(im,jm),dum(im,jm)

rdt=.25/dt

do j=1,jn
  do i=1,in
    ur(i,j)=u(i,j)*xb(i)*rha(j)
    wr(i,j)=w(i,j)*xb(i)*rha(j)
  end do
end do

divm=0.
do i=1,im
  il=i
  ir=i+1
  do j=1,jm
    jb=j
    jt=j+1
    alpha(i,j)=rlx/(a1p(i)*a1p(i)+a2p(j)*a2p(j))/ &
      (.5*xb(il)+.5*xb(ir))
    diver= &
      (0.25*a1p(i)*(ur(ir,jt)+ur(ir,jb)-ur(il,jt)-ur(il,jb)) &
      +0.25*a2p(j)*(wr(ir,jt)-wr(ir,jb)+wr(il,jt)-wr(il,jb)))
    
    divm=dmax1(divm,dabs(diver))
    dum(i,j)=rdt*diver
  end do
end do

if (time.eq.0) then
  pitt=1000.
  write (*,*) 'first time in pressure solver, using pit=1000'
else
  pitt=pit
end if

do sn=1.,pitt
  
  call pforce(dqdt,dudt,drdt,dwdt)
  
  do  j=1,jn
    do  i=1,in
      dqdt(i,j)=dqdt(i,j)*xb(i)*rha(j)
      drdt(i,j)=drdt(i,j)*xb(i)*rha(j)
    end do
  end do
  
  do i=1,im
    il=i
    ir=i+1
    do j=1,jm
      jb=j
      jt=j+1
      pr(i,j)=pr(i,j)-alpha(i,j)*(dum(i,j)+ &
        0.25*a1p(i)*(dqdt(ir,jt)+dqdt(ir,jb) &
        -dqdt(il,jt)-dqdt(il,jb)) &
        +0.25*a2p(j)*(drdt(ir,jt)-drdt(ir,jb) &
        +drdt(il,jt)-drdt(il,jb)))
    end do
  end do
end do 

call pforce(dudt,dudt,dwdt,dwdt)

return
end



subroutine pforce(dqdt,dudt,drdt,dwdt)
implicit real*8 (a-h,o-z)
include 'eza.size'
include 'eza.par'
include 'eza.common'
dimension dwdt(in,jn),dudt(in,jn),dqdt(in,jn),drdt(in,jn)
common/pb/dpdz(im,2)
dimension dpd1(jn),dpd2(jn)

!---- compute pressure gradient terms ----

do i=1,in
  do j=1,jn
    dpd1(j)=a1(i)*0.25*( &
      +(pr(iru(i),jtu(j))-pr(ilu(i),jtu(j))) &
      +(pr(iru(i),jbu(j))-pr(ilu(i),jbu(j))))
    dpd2(j)=a2(j)*0.25*( &
      +(pr(iru(i),jtu(j))-pr(iru(i),jbu(j))) &
      +(pr(ilu(i),jtu(j))-pr(ilu(i),jbu(j))))
  enddo
  
  !  dpd1(1)=dpd1(1)
  !     &   -a1(i)*0.5*(dpdz(iru(i),1)-dpdz(ilu(i),1))
  !  dpd1(jn)=dpd1(jn)
  !     &   +a1(i)*0.5*(dpdz(iru(i),2)-dpdz(ilu(i),2))
  
  do j=1,jn
    dqdt(i,j)=dudt(i,j)-dpd1(j)
  enddo
  do j=2,jn-1
    drdt(i,j)=dwdt(i,j)-dpd2(j)
  end do
end do

return
end


subroutine buoyant(dwdt,buoyvar) !ezan
implicit real*8 (a-h,o-z)
include 'eza.size'
include 'eza.par'
include 'eza.common'
dimension dwdt(in,jn)
dimension buoyvar(in,jn)

coef=1.
if (ray.ne.0.) coef=3.141592*ray

do i=1,in
  do j=1,jn
    dwdt(i,j)=dwdt(i,j)+coef*buoyvar(i,j)
  enddo
enddo
return
end


subroutine coriprim(dudt,dvdt)
implicit real*8 (a-h,o-z)
include 'eza.size'
include 'eza.par'
include 'eza.common'
dimension dudt(in,jn),dvdt(in,jn)
if (axi.gt.0.) then
  do i=2,im
    do j=1,jn
      dvdt(i,j)=dvdt(i,j)-swirl*u(i,j)
      dudt(i,j)=dudt(i,j)+swirl*v(i,j)+v(i,j)*v(i,j)/xa(i)
    end do
  end do
else
  do i=2,im
    do j=1,jn
      dvdt(i,j)=dvdt(i,j)-swirl*u(i,j)
      dudt(i,j)=dudt(i,j)+swirl*v(i,j)
    end do
  end do
end if
return
end

