! v0.2 13 October 2008
!*** numerical integration of 2-d flow, either cartesian or axisymmetric.
!*** ezan.f90 links with ezasub.f90 and blktri.f90
!*** The stress in the primitive equation part needs some work, if the
!*** viscosity varies with height.  This will be included soon .  It
!*** is currently in the improper form used in fiedler's atmosphere-
!*** ocean paper, but causes only a very small error.
!*** The streamfunction vorticity code tends to be more
!*** efficient but is untested in this version.
!*** (Apparently  it provides accuracy and stability with  
!*** lower resolution).  The primitive equation version is easier
!*** to interpret physically, because of the availability of pressure.
!*** Set useprim=1. to use the primitive equations.
!*** The input file in noslip/eza.input gives a vortex similar to 
!*** the "2 Omega= 0.7" of the 1994 A-O paper of B. Fiedler.

program ezan
implicit real*8 (a-h,o-z)
include 'eza.size'
include 'eza.par'
include 'eza.common'
!common/spare/spares(in,jn,5)
common/tends/dvordt(in,jn,3),dbdt(in,jn,3),dvdt(in,jn,3),dwdt(in,jn,3)
common/foraxi/etar(in,jn),vr(in,jn)
common/others/dt1,dt2,rmuc,rmwc,divm,turnon
dimension ap(100)
equivalence (ap(1),time)
character cn*7,filename*80
pi=4.*atan(1.)

write (*,*) " enter directory for input/output of files (e.g. 'myexp'):"
read (*,*) path 

write (*,*) "enter restart timestamp"
write (*,*) "or enter -1 to start fresh" 
read (*,*) icontin

if (icontin.lt.0) then  
  
  call change(1)
  call maketrans  
  
  open (unit=19,file=trim(path)//'/eza.hist',  status='unknown')
  ap(99)=in
  ap(100)=jn
  do n=1,100
    write (19,*) ap(n)
  enddo
  
  do  j=1,jn
    ypi=pi*(j-1.)/jm
    do  i=1,in
      xpi=pi*(i-1.)/im
      v(i,j)=0.
      u(i,j)=0.
      w(i,j)=0.
      stm(i,j)=0.
      vor(i,j)=0.
      !b(i,j)=0.01*cos(xpi)*sin(ypi)
      !if (off.gt.0.)  b(i,j)=1.2636*exp(-20.*.5*xa(i)**2)*exp(-20.*(za(j)-.5)**2)
      rbub=sqrt( (xa(i)/bdfrx)**2 + ( (za(j)-bdfz)/bdfrz)**2 )
      if  (rbub < 1.0)  then  
          bdf(i,j)=(.25/bdfrz) * (1 +cos(rbub*pi)) 
      else
          bdf(i,j)=0.
      endif
      b(i,j)=1./(1.+exp(-(za(j)-binterface)*bsmooth))
    enddo
  enddo
  
  time=0.
  
else
  write(cn,'(i6.6)') icontin 
  filename=trim(path)//'/eza'//trim(cn)//'.dat'
  open (unit=30,file=filename, status='old',form='unformatted')
  open (unit=19,file=trim(path)//'/eza.hist',  status='old',access='append') 
  read (30) ap,b,bdf,stm,vor,u,v,w,viscu,xa,za
  call change(0)
  call maketrans  
  
endif

do j=1,jn
!  zzm=za(j)-.7*hite
!  dummy=rei*(1+.5*cvisc*(1+tanh(4*zzm))) ! old way, but cvisc now undefined
  dummy=dmax1( rei,reit*2*(za(j)+.5-hite) ) ! use reit>rei to make damper 
!  print *, j,za(j),dummy
  do i=1,in
    viscu(i,j)=dummy
  enddo
enddo


write (*,*) ' time=',time,' enter tstop,tdump,tdiag'

read (*,*) tstop,tdump,tdiag

call utop(viscu,viscp)
call prepare

!for ezan bpres function:
iw1=0
do i=1,in
if (xa(i).le.0.5) iw1=i
enddo
iw2=in
do i=in,1,-1
if (xa(i).ge.1.0) iw2=i
enddo

na=1
nb=2
nc=3
it=0
dt1=1.
dt2=1.
rmuc=0.
rmwc=0.

300    continue
if (time.ge.bdfturnon) then
 turnon=1. 
else
 turnon=time/bdfturnon
endif
if (useprim.eq.0.) call stmvor(na,nb,nc,it)
if (useprim.gt.0.) call prim(na,nb,nc,it)

timel=time
time=time+dt
tprint=tprint+dt
it=it+1
dt2=dt1
dt1=dt
nd=nc
nc=nb
nb=na
na=nd

if (int(time/tdump)>int(timel/tdump).or.it.eq.1) then
  if (useprim.gt.0.) call ptou(pr,pru)
!  write(cn,'(i6.6)') int(1000.*time) ! makes time stamp too high resolution 
  write(cn,'(i6.6)') int(100.*time)*10 ! makes time stamp with 00 on end 
  filename=trim(path)//'/eza'//trim(cn)//'.dat'
  write (*,*) 'will write dump file ',trim(filename)
  open(unit=50,file=trim(filename), status='unknown',form='unformatted')
  !write (50) ap,b,stm,vor,u,v,w,viscu,xa,za
  write (50) ap,b,bdf,stm,vor,u,v,w,viscu,xa,za !ezan
  close(unit=50)
endif

if (tprint.gt.tdiag) then
  tprint=0.
  rmu=rmin(u)
  rmv=rmax(v)
  rmw=rmax(w)
  bcontrib=bpres(iw1,iw2)*froudeinv
  if (useprim.ne.0.) then
    rmpp=rmp(pr)
  else
    rmpp=0.
    divm=0.
  endif
  write (*,869) time,dt,rmu,rmv,rmw,rmpp,divm,bcontrib
  write (19,869) time,dt,rmu,rmv,rmw,rmpp,divm,bcontrib
  869    format (1x,f8.3,f8.5,4f10.4,f10.6,f10.4)
  2001    format(1x,1pd15.4,5(1x,1pe11.4))
  
endif

if (time.lt.tstop) goto 300

rewind (unit=30)

end program ezan

subroutine stmvor(na,nb,nc,it)
implicit real*8 (a-h,o-z)
include 'eza.size'
include 'eza.par'
include 'eza.common'
!common/spare/spares(in,jn,5)
common/tends/dvordt(in,jn,3),dbdt(in,jn,3), dvdt(in,jn,3),dwdt(in,jn,3)
common/foraxi/etar(in,jn),vr(in,jn)
common/others/dt1,dt2,rmuc,rmwc
call stuff(dbdt(1,1,na),in,jn,0.d0)
call stuff(dvdt(1,1,na),in,jn,0.d0)
call stuff(dvordt(1,1,na),in,jn,0.d0)

call velocity
call adtime(rmuc,rmwc)

call axi1(etar,vr)
if (order.eq.3.) then
  call advect3(dvordt(1,1,na),etar)
  if (swirl.ne.0.) call advect3(dvdt(1,1,na),vr)
  if (bactive.gt.0.) call advect3(dbdt(1,1,na),b)
  
else if (order.eq.5.) then
  call advect5(dvordt(1,1,na),etar)
  if (swirl.ne.0.) call advect5(dvdt(1,1,na),vr)
  if (bactive.gt.0.) call advect5(dbdt(1,1,na),b)
  
else if (order.eq.0.) then
  call arakawa(dvordt(1,1,na),etar)
  if (swirl.ne.0.) call arakawa(dvdt(1,1,na),vr)
  if (bactive.gt.0.) call arakawa(dbdt(1,1,na),b)
  
else if (order.lt.0.) then
  call marat(dvordt(1,1,na),etar)
  if (swirl.ne.0.) call marat(dvdt(1,1,na),vr)
  if (bactive.gt.0.) call marat(dbdt(1,1,na),b)
  
else
  stop
endif

if (axi.ne.0.) call axi2(dvordt(1,1,na),dvdt(1,1,na))
if (rei.gt.0.and.bactive.gt.0.) call diffuse(dbdt(1,1,na),b)

if (axi.eq.0.) then
  if (rei.gt.0.) call diffuse(dvdt(1,1,na),v)
  if (rei.gt.0.) call diffuse(dvordt(1,1,na),vor)
endif

if (axi.eq.1.) then
  call afdiff(dvdt(1,1,na),v)
  call afdiff(dvordt(1,1,na),vor)
endif

if (axi.ge.2.) then 
  call diffuse(dvdt(1,1,na),v)
  call diffuse(dvordt(1,1,na),vor)
  call rterm(dvordt(1,1,na),vor)
  call rterm(dvdt(1,1,na),v)
endif

call torque(dvordt(1,1,na),b)
if (swirl.ne.0.or.axi.gt.0.) call coriolis(dvordt(1,1,na),dvdt(1,1,na))
if (bactive.gt.0.and.ray.ne.0.) call strat(dbdt(1,1,na))

do i=1,in
  dbdt(i,1,na)=buoybound*dbdt(i,1,na)
  dbdt(i,jn,na)=buoybound*dbdt(i,jn,na)
  dvordt(i,1,na)=0.
  dvordt(i,jn,na)=0.
enddo

do j=1,jn
  dvordt(1,j,na)=0.
  dvordt(in,j,na)=0.
enddo

dts=cfl/dmax1(rmuc,rmwc,.000001d0)
dtdif=diff/(diffmax()+1.e-15)
dt=dmin1(dts,dtdif,dtmax)
if (dt.lt..000001) then
  write (*,*) ' dt too small, model is blowing up: now dt=',dt
  stop
endif

call abcoef(dt1/dt,(dt1+dt2)/dt)
call hop(min0(it+1,3),na,nb,nc,dt,dvordt,vor,in,jn)
if (bactive.gt.0.) call hop(min0(it+1,3),na,nb,nc,dt,dbdt,b,in,jn)
call hop(min0(it+1,3),na,nb,nc,dt,dvdt,v,in,jn)

if (pit.eq.0.) call poiblk(vor,stm,xa,rha)
if (pit.ne.0.) call iterate(vor,stm,xa,rha)

call boundary

return
end

subroutine prim(na,nb,nc,it)
implicit real*8 (a-h,o-z)
include 'eza.size'
include 'eza.par'
include 'eza.common'
common/tends/dudt(in,jn,3),dbdt(in,jn,3), dvdt(in,jn,3),dwdt(in,jn,3)
common/foraxi/etar(in,jn),vr(in,jn)
common/others/dt1,dt2,rmuc,rmwc,divm,turnon

call stuff(dbdt(1,1,na),in,jn,0.d0)
call stuff(dvdt(1,1,na),in,jn,0.d0)
call stuff(dudt(1,1,na),in,jn,0.d0)
call stuff(dwdt(1,1,na),in,jn,0.d0)
call pvelo
call adtime(rmuc,rmwc)

call axi1(etar,vr)
if (order.eq.3.) then
  call advect3(dudt(1,1,na),u)
  call advect3(dwdt(1,1,na),w)
  if (swirl.ne.0.) call advect3(dvdt(1,1,na),vr)
  if (bactive.gt.0.) call advect3(dbdt(1,1,na),b)
  
else if (order.eq.5.) then
  call advect5(dudt(1,1,na),u)
  call advect5(dwdt(1,1,na),w)
  if (swirl.ne.0.) call advect5(dvdt(1,1,na),vr)
  if (bactive.gt.0.) call advect5(dbdt(1,1,na),b)
  
else if (order.lt.0.) then
  call marat(dudt(1,1,na),u)
  call marat(dwdt(1,1,na),w)
  if (swirl.ne.0.) call marat(dvdt(1,1,na),vr)
  if (bactive.gt.0.) call marat(dbdt(1,1,na),b)
  
else
  stop
endif

if (axi.ne.0.) call axi2(etar,dvdt(1,1,na))
if (bactive.gt.0.) call diffuse(dbdt(1,1,na),b)

if (axi.eq.0.) then
  call diffuse(dvdt(1,1,na),v)
  call diffuse(dudt(1,1,na),u)
  call diffuse(dwdt(1,1,na),w)
endif

if (axi.eq.1.) then
  call afdiff(dvdt(1,1,na),v)
  call afdiff(dudt(1,1,na),u)
  call diffuse(dwdt(1,1,na),w)
endif

if (axi.ge.2.) then 
  call diffuse(dvdt(1,1,na),v)
  call diffuse(dudt(1,1,na),u)
  call diffuse(dwdt(1,1,na),w)
  call rterm(dudt(1,1,na),u)
  call rterm(dvdt(1,1,na),v)
endif

call buoyant(dwdt(1,1,na), turnon*bdf)
if (bactive>0)  call buoyant(dwdt(1,1,na),froudeinv*b) 
if (swirl.ne.0.or.axi.gt.0.) call coriprim(dudt(1,1,na),dvdt(1,1,na))
if (bactive.gt.0.and.ray.ne.0.) call strat(dbdt(1,1,na))

do i=1,in
  dbdt(i,1,na)=buoybound*dbdt(i,1,na)
  dbdt(i,jn,na)=buoybound*dbdt(i,jn,na)
  dwdt(i,1,na)=0.
  dwdt(i,jn,na)=0.
enddo

do j=1,jn
  dudt(1,j,na)=0.
  dudt(in,j,na)=0.
enddo
if (slipb.eq.0.) then
  do i=1,in
    dudt(i,1,na)=0.
    u(i,1)=0.
    dvdt(i,1,na)=0.
    v(i,1)=0.
  enddo
endif
if (slips.eq.0.) then
  do j=1,jn
    dwdt(in,j,na)=0.
    w(in,j)=0.
  enddo
endif

dts=cfl/dmax1(rmuc,rmwc,.000001d0)
dtdif=diff/(diffmax()+1.e-15)
dt=dmin1(dts,dtdif,dtmax)
if (dt.lt..00001) then
  write (*,*) ' dt too small, dt=',dt
  stop
endif

call incmp(dudt(1,1,na),dwdt(1,1,na),divm)


call abcoef(dt1/dt,(dt1+dt2)/dt)
call hop(min0(it+1,3),na,nb,nc,dt,dudt,u,in,jn)
call hop(min0(it+1,3),na,nb,nc,dt,dwdt,w,in,jn)
if (bactive.gt.0.) call hop(min0(it+1,3),na,nb,nc,dt,dbdt,b,in,jn)
if (swirl.ne.0.) call hop(min0(it+1,3),na,nb,nc,dt,dvdt,v,in,jn)

return
end


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

do j=1,jn
  do i=1,in
    wf(i,j)=w(i,j)*xb(i)*rha(j)
    uf(i,j)=u(i,j)*xb(i)*rha(j)
  enddo
enddo
return
end 


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

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

function rmp(p)
implicit real*8 (a-h,o-z)
include 'eza.size'
dimension p(im,jm)
dummy=p(im,1)
rmp=p(1,1)-dummy
do  j=1,jm
  do  i=1,im
    rmp=dmin1(rmp,p(i,j)-dummy)
  enddo
enddo
return
end

function bpres(ib1,ib2)
implicit real*8 (a-h,o-z)
include 'eza.size'
include 'eza.par'
include 'eza.common'
common/ezan/pbl(in)
do i=1,in
pbl(i)=0.
do j=jn,2,-1
pbl(i)=pbl(i)-.5*(za(j)-za(j-1))*(b(i,j)+b(i,j-1))
end do
end do
pbavr=0.
do i=ib1,ib2-1
pbavr=pbavr+0.5*(xa(i+1)-xa(i))*(pbl(i)+pbl(i+1))
end do
pbavr=pbavr/(xa(ib2)-xa(ib1))
bpres=pbl(1)-pbavr
return
end

