c
c-----------------------------------------------------------------------
c     file contains barotropic part of BARBI
c
c-----------------------------------------------------------------------
c
      subroutine tropic()
c=======================================================================
c     Integrate the barotropic velocity 
c
c                      C.Eden Feb. 2002 (mailto ceden@phys.ocean.dal.ca)
c=======================================================================
      use barbi_module
      implicit none
      real(kind=8) :: diff_fe(imt,jmt),fxa,fxb,t1,t2,diff_fn(imt,jmt)
      integer      :: i,j,n
c
c-----------------------------------------------------------------------
c     statement functions for advection and diffusion
c-----------------------------------------------------------------------
c
      real (kind=8) :: diff_ux,diff_uy,diff_metric
      DIFF_Ux(i,j) = hur(i,j)*(diff_fe(i,j)-diff_fe(i-1,j))*csudxur(i,j)
      DIFF_Uy(i,j) = hur(i,j)*(diff_fn(i,j) - diff_fn(i,j-1))*csudyur(j)
      DIFF_metric(i,j,n) = am3(j)*u(i,j,n,taum1)+am4(j,n)*dxmetr(i)
     &                    *(u(i+1,j,3-n,taum1) - u(i-1,j,3-n,taum1))
c
c-----------------------------------------------------------------------
c     calculate the gradient of G_1 to force the barotropic mode
c-----------------------------------------------------------------------
c
      grad_E=0.
      do j=1,jmt-1
        fxa  = csur(j)
        fxb  = dyu2r(j)
        do i=1,imt-1
          t1  = G(i+1,j+1,1,tau) - G(i, j,1,tau)
          t2  = G(i,j+1,1,tau) - G(i+1,j,1,tau)
          grad_E(i,j,1) = (t1-t2)*fxa*dxu2r(i)*hur(i,j)
          grad_E(i,j,2) = (t1+t2)*fxb*hur(i,j)
        enddo
      enddo
      if (cyclic) then
       grad_E(1,:,:)   = grad_E(imt-1,:,:)
       grad_E(imt,:,:) = grad_E(2,:,:)
      endif
c
c-----------------------------------------------------------------------
c     loop over the velocity components and solve for forcing
c     of the streamfunction. We are able to neglect the bottom
c     pressure here,  since it does not force the streamfunction.
c     Velocity is calculated afterwards from the stramfunction.
c     Note that velocities are vertical averages here.
c-----------------------------------------------------------------------
c

      do n=1,2

c       diffusive flux across east face of "u" cell
        do j=2,jmt-1
         do i=1,imt-1
          diff_fe(i,j)=am_csudxtr(i,j)
     &            *(u(i+1,j,n,taum1)-u(i,j,n,taum1))
     &            *(h(i+1,j)+h(i+1,j+1))/2.
         enddo
        enddo
c       diffusive flux across north face of "u" cell
        do j=1,jmt-1
         fxa = am*cst(j+1)*dytr(j+1)
         do i=2,imt-1
          diff_fn(i,j) = fxa* (u(i,j+1,n,taum1)-u(i,j,n,taum1))*
     &                          (h(i,j+1)+h(i+1,j+1))/2.
         enddo
        enddo
c
c-----------------------------------------------------------------------
c       calculate the forcing for the streamfunction
c-----------------------------------------------------------------------
c
        do j=2,jmt-1
         do i=2,imt-1
	  zu(i,j,n) = (
     &     +DIFF_Ux(i,j)+DIFF_Uy(i,j) 
     &       + DIFF_metric(i,j,n)
     &     +cori(i,j,n)*u(i,j,3-n,tau) 
     &     -grad_E(i,j,n)
     &     +wind(i,j,n)*hur(i,j))*umask(i,j)
         enddo
        enddo
        if (cyclic) then
         zu(1,:,n)   = zu(imt-1,:,n)
         zu(imt,:,n) = zu(2,:,n)
        endif
      enddo
c
c-----------------------------------------------------------------------
c     Now solve for the streamfunction and get velocites again
c-----------------------------------------------------------------------
c
      call streamfunction
      end subroutine tropic




      subroutine streamfunction
c=======================================================================
c    Construct the forcing for the streamfunction and solve for it
c
c                      C.Eden Feb. 2002 (mailto ceden@phys.ocean.dal.ca)
c=======================================================================
      use barbi_module
      implicit none
      real (kind=8) diag1,diag0, coef(imt,jmt,-1:1,-1:1)
      real (kind=8) forc(imt,jmt),res(imt,jmt)
      real (kind=8) eps,error,fxa,grad_p(imt,jmt,2)
      integer i,j,isle
      real (kind=8) c1ps(imt),c1pn(imt),c1pw(jmt),c1pe(jmt)

      coef=0.; forc=0.; res=0.
      if (eulerfore .or. eulerback1) psi(:,:,2) = psi(:,:,1)

      call sfforc(forc)
      if (cyclic) then
        forc(1,:)   = forc(imt-1,:)
        forc(imt,:) = forc(2,:)
      endif
      call sfc9pt (coef)

      fxa=1. ; if (eulerfore.or.eulerback1.or.eulerback2) fxa=.5
      res = dpsi; dpsi = fxa*(2.*dpsi- guess); guess=res   
      call congrad(coef,dpsi,forc,res,eps_solver,mscans,error,.true.)
      if (cyclic) then
        dpsi(1,:)   = dpsi(imt-1,:)
        dpsi(imt,:) = dpsi(2,:)
      endif
c
c     calculate first the phase velocities for psi
c
      call phase_velocities_psi(c1ps,c1pn,c1pw,c1pe)
c
      if (eulerback2) then
        psi(:,:,1)=psi(:,:,2)+dpsi
      else
        res(:,:)    = psi(:,:,2) + dpsi
        psi(:,:,2)  = psi(:,:,1)
        psi(:,:,1)  = res(:,:)
      endif

      call addobcpsi(c1ps,c1pn,c1pw,c1pe)

      do j=2,jmt-1
       do i=2,imt-1
        diag1       = psi(i+1,j+1,1) - psi(i  ,j,1)
        diag0       = psi(i  ,j+1,1) - psi(i+1,j,1)
	u(i,j,1,taup1)  = -(diag1+diag0)*dyu2r(j)
	u(i,j,2,taup1)  =  (diag1-diag0)*dxu2r(i)*csur(j)
       enddo
      enddo
      u(:,:,1,taup1) = u(:,:,1,taup1)*umask*hur
      u(:,:,2,taup1) = u(:,:,2,taup1)*umask*hur
      if (cyclic) then
        u(1,:,:,taup1)   = u(imt-1,:,:,taup1)
        u(imt,:,:,taup1) = u(2,:,:,taup1)
      endif
      if (enable_obc_south) u(:,1,:,taup1) = u(:,2,:,taup1)
      if (enable_obc_north) then
       u(:,jmt-1,:,taup1) = u(:,jmt-2,:,taup1)
       u(:,jmt  ,:,taup1) = u(:,jmt-2,:,taup1)
      endif
      if (enable_obc_west)  u(1,:,:,taup1) = u(2,:,:,taup1)
      if (enable_obc_east) then
        u(imt  ,:,:,taup1) = u(imt-2,:,:,taup1)
        u(imt-1,:,:,taup1) = u(imt-2,:,:,taup1)
      endif
c
c     now solve for the bottom pressure
c
      if (snapshot_time_step) then
       grad_p(:,:,1) = (U(:,:,1,taup1)-U(:,:,1,taum1))/c2dt-zu(:,:,1)
       grad_p(:,:,2) = (U(:,:,2,taup1)-U(:,:,2,taum1))/c2dt-zu(:,:,2)
       if (cyclic) then
        grad_p(1,:,:)   = grad_p(imt-1,:,:)
        grad_p(imt,:,:) = grad_p(2,:,:)
       endif

       coef=0.; forc=0.; res=0.; press=0.
       call spforc (grad_p, forc )
       if (cyclic) then
        forc(1,:)   = forc(imt-1,:)
        forc(imt,:) = forc(2,:)
       endif
       call spc9pt (coef)
       call congrad(coef,press,forc,res,eps_solver_p,mscans_p,
     &              error,.false.)
       if (cyclic) then
        press(1,:)   = press(imt-1,:)
        press(imt,:) = press(2,:)
       endif
      endif
      end subroutine streamfunction



      subroutine phase_velocities_psi(c1ps,c1pn,c1pw,c1pe)
c=======================================================================
c     calculate phase velocities for psi
c     Orlanski radiation condition (passive open boundary)   
c
c                      C.Eden Feb. 2002 (mailto ceden@phys.ocean.dal.ca)
c=======================================================================
      use barbi_module
      implicit none
      real (kind=8) var,var1,c1ps(imt),c1pn(imt),c1pw(jmt),c1pe(jmt)
      integer i,j
      if (enable_obc_south ) then
       var=-dyu(3)/dt
       do 4030 i=1,imt
         var1=psi(i,4,2)-psi(i,3,2)
         if (var1== 0.) then
           c1ps(i)=var
           goto 4030
         endif
         c1ps(i)=var*(psi(i,3,1)-psi(i,3,2))/var1
         if (c1ps(i).ge. 0.) c1ps(i)=0.
         if (c1ps(i).lt.var) c1ps(i)=var
4030   continue
      endif
c
      if (enable_obc_north) then
       var=dyu(jmt-2)/dt
       do 4031 i=1,imt
         var1=psi(i,jmt-2,2)-psi(i,jmt-3,2)
         if (var1.eq.0.) then
           c1pn(i)=var
           goto 4031
         endif
         c1pn(i)=-var*(psi(i,jmt-2,1)-psi(i,jmt-2,2))/var1
         if (c1pn(i).le. 0.) c1pn(i)=0.
         if (c1pn(i).gt.var) c1pn(i)=var
4031   continue
      endif
c
      if (enable_obc_west) then
       do 4032 j=1,jmt
         var=-dxu(3)*csu(j)/dt
         var1=psi(4,j,2)-psi(3,j,2)
         if (var1.eq.0.) then
           c1pw(j)=var
           goto 4032
         endif
         c1pw(j)=var*(psi(3,j,1)-psi(3,j,2))/var1
         if (c1pw(j).ge. 0.) c1pw(j)=0.
         if (c1pw(j).lt.var) c1pw(j)=var
4032   continue
      endif
c
      if (enable_obc_east) then
       do 4033 j=1,jmt
         var=dxu(imt-2)*csu(j)/dt
         var1=psi(imt-2,j,2)-psi(imt-3,j,2)
         if (var1.eq.0.) then
           c1pe(j)=var
           goto 4033
         endif
         c1pe(j)=-var*(psi(imt-2,j,1)-psi(imt-2,j,2))/var1
         if (c1pe(j).le. 0.) c1pe(j)=0.
         if (c1pe(j).gt.var) c1pe(j)=var
4033   continue
      endif
      end subroutine phase_velocities_psi



      subroutine addobcpsi(c1ps,c1pn,c1pw,c1pe)
c=======================================================================
c     passive Orlanski radiation condition (passive open boundary)   
c
c                      C.Eden Feb. 2002 (mailto ceden@phys.ocean.dal.ca)
c=======================================================================
      use barbi_module
      implicit none
      integer i,j
      real (kind=8) :: c1ps(imt),c1pn(imt),c1pw(jmt),c1pe(jmt)
      if (enable_obc_south) then
        do i=2,imt-1
         psi(i,2,1) =psi(i,2,2)  - c1ps(i)*dt*dyur(2)
     &                          *(psi(i,3,2)-psi(i,2,2))
         psi(i,1,1) =psi(i,2,1)
        enddo
      endif

      if (enable_obc_north) then
       do i=2,imt-1
         psi(i,jmt-1,1)=psi(i,jmt-1,2)- c1pn(i)*dt*dyur(jmt-1)
     &                          *(psi(i,jmt-1,2)- psi(i,jmt-2,2))
         psi(i,jmt,1)  =psi(i,jmt-1,1)
       enddo
      endif
c
      if (enable_obc_west ) then
       do j=2,jmt-1
         psi(2,j,1)   = psi(2,j,2)
     &                       - c1pw(j)*dt*dxur(2)*csur(j)
     &                       *(psi(3,j,2) -psi(2,j,2))
         psi(1,j,1)   = psi(2,j,1)
       enddo
      endif

      if (enable_obc_east) then
       do j=2,jmt-1
         psi(imt-1,j,1) = psi(imt-1  ,j,2)
     &                    - c1pe(j)*dt*dxur( imt-1 )*csur(j)
     &                    *(psi(imt-1,j,2)- psi(imt-2,j,2))
         psi(imt,j,1)     = psi(imt-1,j,1)
       enddo
      endif
c-----------------------------------------------------------------------
c      call adjust_islands_obc
      end subroutine addobcpsi





      subroutine sfforc (forc)
c
c=======================================================================
c           S T R E A M   F U N C T I O N   F O R C I N G
c=======================================================================
c
      use barbi_module
      implicit none
      real(kind=8), dimension( 0:1, 0:1)     :: cddxu,  cddyu
      real(kind=8), dimension(-1:0,-1:0)     :: cddxt,  cddyt
      real(kind=8), dimension(imt,jmt)       :: ustuff,vstuff,forc
      real(kind=8) p5
      integer i,j,i1,j1,i2,j2

      p5=0.5
      cddxu( 0, 0) = -p5; cddxu( 0, 1) = -p5
      cddxu( 1, 0) =  p5; cddxu( 1, 1) =  p5
c
      cddxt(-1,-1) = -p5; cddxt(-1, 0) = -p5
      cddxt( 0,-1) =  p5; cddxt( 0, 0) =  p5
c
      cddyu( 0, 0) = -p5; cddyu( 0, 1) =  p5
      cddyu( 1, 0) = -p5; cddyu( 1, 1) =  p5
c
      cddyt(-1,-1) = -p5; cddyt(-1, 0) =  p5
      cddyt( 0,-1) = -p5; cddyt( 0, 0) =  p5
c
c-----------------------------------------------------------------------
c     multiply the u eqn by dx*cos, the v eqn by dy, then subtract their 
c     partial differences to eliminate the unknown surface pressure from 
c     the resulting equation
c-----------------------------------------------------------------------
c
      ustuff=0.
      vstuff=0.
      do j=1,jmt-1
        do i=1,imt-1
          ustuff(i,j) = zu(i,j,1)*dxu(i)*csu(j)
          vstuff(i,j) = zu(i,j,2)*dyu(j)
        end do
      end do
c
      forc=0.0
      do j1=-1,0
        do j=2,jmt-1
          do i=2,imt-1
            do i1=-1,0
              forc(i,j) = forc(i,j) 
     &                     - cddyt(i1,j1)*ustuff(i+i1,j+j1)
     &                     + cddxt(i1,j1)*vstuff(i+i1,j+j1)
            end do
          end do
        end do
      end do
c
      end subroutine sfforc




      subroutine spforc (grad_p, forc)
c
c=======================================================================
c           S U R F A C E   P R E S S U R E   F O R C I N G
c     author:      c. h. goldberg      e-mail=> chg@gfdl.gov
c=======================================================================
c
      use barbi_module
      implicit none
      real(kind=8) :: grad_p(imt,jmt,2), forc(imt,jmt)
      real(kind=8), dimension(imt,jmt) :: ustuff, vstuff
      real(kind=8), dimension( 0:1, 0:1)     :: cddxu,  cddyu
      real(kind=8), dimension(-1:0,-1:0)     :: cddxt,  cddyt
      real(kind=8), parameter :: p5=0.5
      integer ::  i,j,i1,j1

      forc = 0.
c
      cddxu( 0, 0) = -p5; cddxu( 0, 1) = -p5;
      cddxu( 1, 0) =  p5; cddxu( 1, 1) =  p5
      cddxt(-1,-1) = -p5; cddxt(-1, 0) = -p5
      cddxt( 0,-1) =  p5; cddxt( 0, 0) =  p5
      cddyu( 0, 0) = -p5; cddyu( 0, 1) =  p5
      cddyu( 1, 0) = -p5; cddyu( 1, 1) =  p5
      cddyt(-1,-1) = -p5; cddyt(-1, 0) =  p5
      cddyt( 0,-1) = -p5; cddyt( 0, 0) =  p5

      do j=1,jmt-1
        do i=1,imt-1
	  ustuff(i,j) = hu(i,j)*grad_p(i,j,1)*dyu(j)
	  vstuff(i,j) = hu(i,j)*grad_p(i,j,2)*dxu(i)*csu(j)
        end do
      end do
c
      do j1=-1,0
	do j=2,jmt-1
	  do i=2,imt-1
             do i1=-1,0
	      forc(i,j) = forc(i,j)
     &                     + cddxt(i1,j1)*ustuff(i+i1,j+j1)
     &                     + cddyt(i1,j1)*vstuff(i+i1,j+j1)
	    end do
	  end do
	end do
      end do
      end subroutine spforc




      subroutine sfc9pt  (coef)
c
c=======================================================================
c     9  P T    C O E F F I C I E N T   I N I T I A L I A Z A T I O N
c     coefficient initialization for 9 point elliptic solvers
c=======================================================================
c
c-----------------------------------------------------------------------
c     generate arrays of coefficients
c     construct coefficients for partial differences. a partial
c     difference in "x" is defined as an "x" difference of a quantity 
c     which is averaged in "y". (and symmetrically for "y" differences).
c     Note that this is an x difference and NOT an x derivitive.
c     partial differences of quantities on the "t" grid are defined on
c     the "u" grid and visa versa.
c     therefore partial differences at: 
c     u/v points (i,j), involve nearby t/s points with subscripts: 
c        (i  ,j+1)    (i+1,j+1)
c        (i  ,j  )    (i+1,j  )
c     t/s points (i,j), involve nearby u/v points with subscripts:
c        (i-1,j  )    (i  ,j  )
c        (i-1,j-1)    (i  ,j-1)
c     thus if qu(i,j) is defined on u/v points, its partial
c     difference ddxqt = ddxt(qu) is defined on t/s points and has the
c     value 
c     ddxqt(i,j) = cddxt(-1,-1)*qu(i-1,j-1) + cddxt(-1,0)*qu(i-1,j+0)
c                + cddxt( 0,-1)*qu(i+0,j-1) + cddxt( 0,0)*qu(i+0,j+0)
c-----------------------------------------------------------------------
c
      use barbi_module
      implicit none
      real(kind=8) :: coef(imt,jmt,-1:1,-1:1),p5=0.5
      real(kind=8), dimension( 0:1, 0:1)     :: cddxu,  cddyu
      real(kind=8), dimension(-1:0,-1:0)     :: cddxt,  cddyt
      real(kind=8), dimension(imt,jmt)       :: ustuff,vstuff
      integer i,j,i1,j1,i2,j2
c
      cddxu( 0, 0) = -p5; cddxu( 0, 1) = -p5
      cddxu( 1, 0) =  p5; cddxu( 1, 1) =  p5
      cddxt(-1,-1) = -p5; cddxt(-1, 0) = -p5
      cddxt( 0,-1) =  p5; cddxt( 0, 0) =  p5
      cddyu( 0, 0) = -p5; cddyu( 0, 1) =  p5
      cddyu( 1, 0) = -p5; cddyu( 1, 1) =  p5
      cddyt(-1,-1) = -p5; cddyt(-1, 0) =  p5
      cddyt( 0,-1) = -p5; cddyt( 0, 0) =  p5
c
c-----------------------------------------------------------------------
c     compute coefficients for all points
c-----------------------------------------------------------------------
c
      ustuff=0.
      vstuff=0.
      do j=1,jmt-1
       do i=1,imt-1
          ustuff(i,j) = dxu(i)*csu(j) / (c2dt*dyu(j)) *hur(i,j) 
          vstuff(i,j) = dyu(j) / (c2dt*dxu(i)*csu(j)) *hur(i,j) 
        end do
      end do
c
c---------------------------------------------------------------------
c     calculate 9 point coefficients
c---------------------------------------------------------------------
c
      coef=0.0
      do j1=0,1
         do j2=-1,0
           do j=2,jmt-1
      do i1=0,1
          do i2=-1,0
                do  i=2,imt-1
                  coef(i,j,i1+i2,j1+j2) = coef(i,j,i1+i2,j1+j2) +
     &               cddyu(i1,j1)*cddyt(i2,j2)*ustuff(i+i2,j+j2)  +
     &               cddxu(i1,j1)*cddxt(i2,j2)*vstuff(i+i2,j+j2)
                end do
              end do
            end do
          end do
        end do
      end do

      end subroutine sfc9pt



      subroutine spc9pt (coef)
c
c=======================================================================
c     S U R F A C E   P R E S S U R E    C O E F F I C I E N T
c     author: c.h. goldberg         e-mail => chg@gfdl.gov
c=======================================================================
c
      use barbi_module
      implicit none

      real(kind=8) :: coef(imt,jmt,-1:1,-1:1)
      real(kind=8) :: ustuff(imt,jmt), vstuff(imt,jmt)
      real(kind=8) :: cddxu(0:1,0:1),   cddyu(0:1,0:1)
      real(kind=8) :: cddxt(-1:0,-1:0), cddyt(-1:0,-1:0)
      real(kind=8), parameter :: p5=0.5
      integer :: i,j,i1,j1,i2,j2
c
      cddxu( 0, 0) = -p5; cddxu( 0, 1) = -p5
      cddxu( 1, 0) =  p5; cddxu( 1, 1) =  p5
      cddxt(-1,-1) = -p5; cddxt(-1, 0) = -p5
      cddxt( 0,-1) =  p5; cddxt( 0, 0) =  p5
      cddyu( 0, 0) = -p5; cddyu( 0, 1) =  p5
      cddyu( 1, 0) = -p5; cddyu( 1, 1) =  p5
      cddyt(-1,-1) = -p5; cddyt(-1, 0) =  p5
      cddyt( 0,-1) = -p5; cddyt( 0, 0) =  p5
c
      coef = 0.0; ustuff = 0.0; vstuff = 0.0
      do i=1,imt-1
        do j=1,jmt-1
          ustuff(i,j) = hu(i,j)*dyu(j)/(dxu(i)*csu(j))
          vstuff(i,j) = hu(i,j)*dxu(i)*csu(j)/dyu(j)
        end do
      end do
c
c     calculate divergence = ddx (ddx (ustuff)) + ddy( ddy (vstuff))
c
      do j1=0,1
        do j2=-1,0
          do j=2,jmt-1
            do i1=0,1
              do i2=-1,0
                do  i=2,imt-1
                  coef(i,j,i1+i2,j1+j2) = coef(i,j,i1+i2,j1+j2)
     &                 + cddxu(i1,j1) * cddxt(i2,j2) * ustuff(i+i2,j+j2)
     &                 + cddyu(i1,j1) * cddyt(i2,j2) * vstuff(i+i2,j+j2)
                end do
              end do
            end do
          end do
        end do
      end do
      end subroutine spc9pt 


