
c
c
c =========================================================
      subroutine src1(maxmx,meqn,mbc,mx,xlower,dx,q,maux,aux,t,dt)
c =========================================================
      implicit double precision (a-h,o-z)
      dimension q(1-mbc:maxmx+mbc, meqn)
      dimension qstar(3)
      dimension f_vec(3),u_vec(3),out_vec(3)
c
c     geometric source term + Coriolis force
c
      dt2 = dt/2.d0
      f0 = 0.d0 !4.0d0

      do i = 1-mbc, mx+mbc

         xc = xlower + (i-0.5d0)*dx

c        # Stage 1
         h = q(i,1)
         u = q(i,2)/h
         v = q(i,3)/h

         rhs1 = h*u*dtan(xc)
         rhs2 = h*u*u*dtan(xc) - h*v*v*dsin(xc)*dcos(xc)
         rhs3 = 3.d0*h*u*v*dtan(xc)

         rhs2 = rhs2 - h*v*dcos(xc)*dsin(xc)*f0
         rhs3 = rhs3 + h*u*dtan(xc)*f0

         qstar(1) = q(i,1) + dt2*rhs1
         qstar(2) = q(i,2) + dt2*rhs2
         qstar(3) = q(i,3) + dt2*rhs3


c        # Stage 2
         h = qstar(1)
         u = qstar(2)/h
         v = qstar(3)/h

         rhs1 = h*u*dtan(xc)
         rhs2 = h*u*u*dtan(xc) - h*v*v*dsin(xc)*dcos(xc)
         rhs3 = 3.d0*h*u*v*dtan(xc)

         rhs2 = rhs2 - h*v*dcos(xc)*dsin(xc)*f0
         rhs3 = rhs3 + h*u*dtan(xc)*f0

         q(i,1) = q(i,1) + dt*rhs1
         q(i,2) = q(i,2) + dt*rhs2
         q(i,3) = q(i,3) + dt*rhs3
      enddo
c
c
      return
      end
