c
c   program to do van leer in 2d on a sphere.  Based on Dales ETA code.
c
      parameter (long$=144,latedge$=88,latcen$=latedge$+1,polebox=3) 
      real*4 const(long$,latcen$),coscen(latcen$),cosedge(latedge$),
     1       u(long$,latedge$),v(long$,latedge$),
     2       f1(long$,latcen$),f2(long$,latcen$)
     

c
c     read in stuff
c
      read (10) u
      read (10) v
      read (11) const
      read (12) coscen,cosedge
c
      pi=4.*atan(1.)
      dphi=(180-2.*polebox)/(latedge$-1.)*pi/180.
c
      dt=1.e4
      dt=1.e3
      ttot=0.
      ictout=0
      iout=50
      dt2=dt/2.
c
      write (20) ttot,dt,dphi,coscen,cosedge
      write (20) const,u,v
c
      do itime=1,1000
c
      call xadv(const,coscen,dt2,u,f1)
      call yadv(const,coscen,cosedge,dt2,v,dphi,f2,polebox)
      do j=1,latcen$
      do i=1,long$
      const(i,j)=const(i,j)+f1(i,j)+f2(i,j)
      end do
      end do
c
      call yadv(const,coscen,cosedge,dt2,v,dphi,f2,polebox)
      call xadv(const,coscen,dt2,u,f1)
      do j=1,latcen$
      do i=1,long$
      const(i,j)=const(i,j)+f1(i,j)+f2(i,j)
      end do
      end do
c
      ttot=ttot+dt
      ictout=ictout+1
c
      if (ictout .eq. iout) then
      ictout=0
      write (20) ttot
      write (20) const,f1,f2
      end if
c
c     end time loop
c
      end do
      end
c
c 
      subroutine xadv(denjk,coscen,dt,u1,f1)
c
      parameter (long$=144,latedge$=88,latcen$=latedge$+1) 
      real*4 denjk(long$,latcen$),coscen(latcen$),
     1       u1(long$,latedge$),f1(long$,latcen$)
c
      real*4 delta(long$),flux(long$),
     1       top(long$),bot(long$),
     2       dnew(long$),d(long$),dtovdxl(latcen$)
c
      a=6.37e6
      pi=4.*atan(1.)
      dist=2.*pi*a/float(long$)
c
c     cos poles = 0, pole points will be treated separately
c
      do j=2,latcen$-1
      dtovdxl(j)=dt/(dist*coscen(j))
      end do
c
c    loop over interior centers, ie not the pole
c
      do 100 j=2,latcen$-1
c
      do i=1,long$
      d(i)=denjk(i,j)
      delta(i)=0.
      end do
c 
      do i=1,long$
      im=i-1
      ip=i+1
      if (im .lt. 1) im=long$
      if (ip .gt. long$) ip=1
      top(i)=(d(i)-d(im))*(d(ip)-d(i))
      bot(i)=d(ip)-d(im)
c
c    set delta
c
      if (top(i) .gt. 0. .and. bot(i) .ne. 0.) then
      delta(i)=2.*top(i)/bot(i)
      else
      delta(i)=0.
      end if
      end do
c
c    calculate dale's fluxes,  at the edge of each box
c
c
      s=1.
      do i=1,long$
      ip=i+1
      if (ip .gt. long$) ip=1
      if (u1(i,j) .ge. 0. ) then
      flux(i) = d(i)+(.5*(1.-u1(i,j)*dtovdxl(j)))*delta(i)*s
      else
      flux(i) = d(ip)-(.5*(1.+u1(i,j)*dtovdxl(j)))*delta(ip)*s
      end if
      end do
c
c    calculate the changes
c
      do i=1,long$
      im=i-1
      ip=i+1
      if (im .lt. 1) im=long$
      if (ip .gt. long$) ip=1
      dnew(i)=-dtovdxl(j)*
     1          (u1(i,j)*flux(i)-u1(im,j)*flux(im))
      f1(i,j)=dnew(i)
      end do
  100 continue
c
c    take care of the poles
c
      do i=1,long$
      j=1
      f1(i,j)=0.
      j=latcen$
      f1(i,j)=0.
      end do
c
      return
      end
      subroutine yadv(denjk,coscen,cosedge,dt,vfld,dphi,f1,polebox)
c
      parameter (long$=144,latedge$=88,latcen$=latedge$+1) 
      real*4 denjk(long$,latcen$),coscen(latcen$),cosedge(latedge$),
     1       vfld(long$,latedge$),f1(long$,latcen$)
c
      real*4 delta(long$,latcen$),flux(long$,latedge$),
     1       top(long$,latcen$),bot(long$,latcen$),
     2       dnew(long$,latcen$),d(long$,latcen$)
c
      a=6.37e6
      dtovdxj=dt/(dphi*a)
      pi=4.*atan(1.)
      dphi2=polebox*pi/180.
      polcon=2.*dt/(a*long$*sin(dphi2))
c
      do j=1,latcen$
      do i=1,long$
      d(i,j)=denjk(i,j)
      delta(i,j)=0.
      end do
      end do
c
      do 100 i=1,long$
c
c    loop over interior centers, ie not the pole
c
      do j=2,latcen$-1
      jm=j-1
      jp=j+1
      top(i,j)=(d(i,j)-d(i,jm))*(d(i,jp)-d(i,j))
      bot(i,j)=d(i,jp)-d(i,jm)
      end do
c
c    now the poles go here.  Step 1. leave the poles out
c   
      j=1
      top(i,j)=0.
      bot(i,j)=1.
      j=latcen$
      top(i,j)=0.
      bot(i,j)=1.
c
c    end long$ loop
c
  100 continue
c
c    set delta
c
      do j=1,latcen$
      do i=1,long$
      if (top(i,j) .gt. 0. .and. bot(i,j) .ne. 0.) then
      delta(i,j)=2.*top(i,j)/bot(i,j)
      else
      delta(i,j)=0.
      end if
      end do
      end do
c
c    calculate dale's fluxes,  at the edge of each box
c
c    notation changes from dale here.
c
      s=1.
      do j=1,latedge$
      do i=1,long$
      if (vfld(i,j) .ge. 0. ) then
      flux(i,j) = d(i,j)+(.5*(1.-vfld(i,j)*dtovdxj))*delta(i,j)*s
      else
      flux(i,j) = d(i,j+1)-(.5*(1.+vfld(i,j)*dtovdxj))*delta(i,j+1)*s
      end if
      end do
      end do
c
c    calculate the changes
c
c    loop over interior centers, ie don't do the poles
c    notation changes from dale here.
c
      do 200 i=1,long$
c
      do j=2,latcen$-1
      dnew(i,j)=-(dtovdxj/coscen(j))*
     1          (vfld(i,j)*flux(i,j)*cosedge(j)-
     2           vfld(i,j-1)*flux(i,j-1)*cosedge(j-1))
      end do
c
c    do the polar boxes
c
      j=1
      dnew(i,j)=vfld(i,j)*flux(i,j)
      j=latcen$
      dnew(i,j)=vfld(i,j-1)*flux(i,j-1)
c
c    end longitude loop
c
  200 continue
c
c    now work the polar boxes
c
      sumsp=0.
      sumnp=0.
      do i=1,long$
      j=1
      sumsp=sumsp-polcon*dnew(i,j)
      j=latcen$
      sumnp=sumnp+polcon*dnew(i,j)
      end do
c
c    well mixed polar box
c
      do i=1,long$
      j=1
      dnew(i,j)=sumsp
      j=latcen$
      dnew(i,j)=sumnp
      end do
c
c    change array for main routine
c
      do j=1,latcen$
      do i=1,long$
cc      if (i .eq. 9 .or. i .eq. 19) write (6,10) dnew(i,j)
c   10 format (e15.5)
      f1(i,j)=dnew(i,j)
      end do
      end do
c
      return
      end
