#include "options.inc"




      subroutine rigid_lid
c
c=======================================================================
c     Solve for the barotropic velocities for "tau+1"
c     original author:   r.c.pacanowski       e-mail rcp@gfdl.gov
c     SPFLAME version : c.eden
c=======================================================================
c
      use spflame_module
      implicit none
      real, dimension(is_pe:ie_pe,js_pe:je_pe,-1:1,-1:1) :: coef
c     Note: coefficients could be saved and computed only once
      real, dimension(is_pe:ie_pe,js_pe:je_pe)           :: forc
      real, dimension(is_pe-1:ie_pe+1,js_pe-1:je_pe+1)   :: res
      integer   is,ie,js,je,i,j
c
c-----------------------------------------------------------------------
c     prepare for a mixing time step
c-----------------------------------------------------------------------
c
      if (eulerfore .or. eulerback1) then
        psi(:,:,2) = psi(:,:,1)
        ptd=ptd/2.
        guess=guess/2.
      endif

#ifdef detailed_timing
      call tic('tropic zu exchg')
#endif
      call border_exchg(zu(:,:,1),1,1)
      call border_exchg(zu(:,:,2),1,1)
      call set_cyclic(zu(:,:,1),1,1)
      call set_cyclic(zu(:,:,2),1,1)
#ifdef detailed_timing
      call toc('tropic zu exchg')
#endif
c
c-----------------------------------------------------------------------
c     construct the forcing for the stream function equation
c-----------------------------------------------------------------------
c
#ifdef detailed_timing
      call tic('tropic forc')
#endif
      call sfforc (forc)
c
c-----------------------------------------------------------------------
c     initialize coefficients using 9 point numerics
c-----------------------------------------------------------------------
c
      call sfc9pt (coef)
c
c-----------------------------------------------------------------------
c     linear extrapolation between the last two solution for the
c     purpose of computing an initial guess for the present solution.
c-----------------------------------------------------------------------
c
      res = ptd; ptd = (2.*ptd- guess); guess=res ; res=0.

      call border_exchg(ptd,1,1)
      call set_cyclic(ptd,1,1)
#ifdef detailed_timing
      call toc('tropic forc')
#endif
c
c-----------------------------------------------------------------------
c     solving for the "tau+1" stream function change
c-----------------------------------------------------------------------
c
#ifdef detailed_timing
      call tic('tropic congr')
#endif
      call congr (ptd, res, forc, coef, eps_solver,.true.)
#ifdef detailed_timing
      call toc('tropic congr')
#endif
c
c     Distribute new solution
c
      call border_exchg(ptd,1,1)
      call set_cyclic(ptd,1,1)
c
c     calculate first the phase velocities for psi
c
#ifdef detailed_timing
      call tic('tropic phase')
#endif
      call phase_velocities_psi
#ifdef detailed_timing
      call toc('tropic phase')
#endif
c
c-----------------------------------------------------------------------
c     update the stream function based upon the solution
c-----------------------------------------------------------------------
c
      if (eulerback2) then
        psi(:,:,1)=psi(:,:,2)+ptd
      else
        res(:,:)    = psi(:,:,2) + ptd
        psi(:,:,2)  = psi(:,:,1)
        psi(:,:,1)  = res(:,:)
      endif
c
c     apply open boundary conditions
c
#ifdef detailed_timing
      call tic('tropic obc')
#endif
      call addobcpsi
#ifdef detailed_timing
      call toc('tropic obc')
#endif
c
c     Distribute psi
c
      call border_exchg(psi(:,:,1),1,1)
      call set_cyclic(psi(:,:,1),1,1)
c
c     add external mode to the internal velocities
c
#ifdef detailed_timing
      call tic('tropic addext')
#endif
      call add_ext_mode
#ifdef detailed_timing
      call toc('tropic addext')
#endif
c
c-----------------------------------------------------------------------
c     save ptd to compute 1st guess for relaxation next timestep
c     (..note.. on 1st pass of euler backward timestep, bypass this
c            save, since it will be done on the 2nd pass)
c     (..note.. on a mixing timestep, alter ptd to be consistent with
c            normal, leap-frog stepping)
c-----------------------------------------------------------------------
c
      if( eulerfore .or. eulerback2 ) then
         ptd=2.*ptd
         guess=guess*2.
      endif

      end subroutine rigid_lid




      subroutine add_ext_mode ()
      use spflame_module
      implicit none
      integer i,j,n,k,is,ie,js,je
      real diag1,diag0,ext1,ext2

      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)

      do j=js,je
        do i=is,ie
         diag1       = psi(i+1,j+1,1) - psi(i  ,j,1)
         diag0       = psi(i  ,j+1,1) - psi(i+1,j,1)
	 ext1  = -(diag1+diag0)*dyu2r(j)*hr(i,j)
	 ext2  =  (diag1-diag0)*dxu2r(i)*hr(i,j)*csur(j)
         do k=1,km
          u(i,k,j,1,taup1) = (u(i,k,j,1,taup1)+ext1)*umask(i,k,j)
          u(i,k,j,2,taup1) = (u(i,k,j,2,taup1)+ext2)*umask(i,k,j)
         enddo
	enddo
      enddo

      do n=1,2
       call set_cyclic(u(:,:,:,n,taup1),km,2)
       if (my_blk_j == 1 .and. enable_obc_south) 
     &        u(:,:,1,n,taup1) = u(:,:,2,n,taup1)
       if (my_blk_j == n_pes_j .and. enable_obc_north) then
              u(:,:,jmt-1,n,taup1) = u(:,:,jmt-2,n,taup1)
              u(:,:,jmt  ,n,taup1) = u(:,:,jmt-2,n,taup1)
       endif
       if (my_blk_i == 1 .and. enable_obc_west) 
     &        u(1,:,:,n,taup1) = u(2,:,:,n,taup1)
       if (my_blk_i == n_pes_i .and. enable_obc_east) then
              u(imt  ,:,:,n,taup1) = u(imt-2,:,:,n,taup1)
              u(imt-1,:,:,n,taup1) = u(imt-2,:,:,n,taup1)
       endif
       call border_exchg(u(:,:,:,n,taup1),km,2)
      enddo

      end subroutine add_ext_mode


      subroutine sfforc (forc)
      use spflame_module
      implicit none
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
      real, dimension(is_pe:ie_pe,js_pe:je_pe)   :: forc
      real, dimension(is_pe-1:ie_pe,js_pe-1:je_pe) :: ustuff, vstuff
      real, dimension( 0:1, 0:1)     :: cddxu,  cddyu
      real, dimension(-1:0,-1:0)     :: cddxt,  cddyt
      real :: p5=0.5
      integer i,j,i1,j1,is,ie,js,je
c
      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)
c
c-----------------------------------------------------------------------
c     initialize the forcing
c-----------------------------------------------------------------------
      ustuff=0.; vstuff=0.; forc=0.
      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     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
      do j=js-1,je
        do i=is-1,ie
          ustuff(i,j) = zu(i,j,1)*dxu(i)*csu(j)
          vstuff(i,j) = zu(i,j,2)*dyu(j)
        end do
      end do
      do j1=-1,0
        do j=js,je
          do i=is,ie
            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
      end subroutine sfforc



      subroutine spforc (nabla_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 spflame_module
      implicit none
      real, dimension(is_pe-1:ie_pe+1,js_pe-1:je_pe+1,2) :: nabla_p
      real, dimension(is_pe:ie_pe,js_pe:je_pe)   :: forc
      real, dimension(is_pe-1:ie_pe+1,js_pe-1:je_pe+1) :: ustuff,vstuff
      real, dimension( 0:1, 0:1)     :: cddxu,  cddyu
      real, dimension(-1:0,-1:0)     :: cddxt,  cddyt
      real :: h
      real, parameter :: p5=0.5
      integer ::  i,j,i1,j1,is,ie,js,je
c
      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)

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

      forc = 0.; ustuff=0.; vstuff=0.

      do j=js,je
        do i=is,ie
#ifdef partial_cell
          h = hup(i,j)
#else
          h = zw(max(1,kmu(i,j)))*umask(i,1,j)
          
#endif
	  ustuff(i,j) = h*nabla_p(i,j,1)*dyu(j)
	  vstuff(i,j) = h*nabla_p(i,j,2)*dxu(i)*csu(j)
        end do
      end do
      call border_exchg(ustuff,1,1)
      call set_cyclic(  ustuff,1,1)
      call border_exchg(vstuff,1,1)
      call set_cyclic(  vstuff,1,1)
c
      do j1=-1,0
	do j=js,je
	  do i=is,ie
             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)
      use spflame_module
      implicit none
c
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
c=======================================================================
c
      real coef  (is_pe:ie_pe,js_pe:je_pe,-1:1,-1:1)
      real ustuff(is_pe-1:ie_pe+1,js_pe-1:je_pe+1)
      real vstuff(is_pe-1:ie_pe+1,js_pe-1:je_pe+1)
      real cddxu(0:1,0:1), cddyu(0:1,0:1)
      real cddxt(-1:0,-1:0), cddyt(-1:0,-1:0)
      real :: p5=0.5
      integer is,ie,js,je, i,j,i1,j1,i2,j2
c
      js = max(2,js_pe); je = min(jmt-1,je_pe)
      is = max(2,is_pe); ie = min(imt-1,ie_pe)
c
c-----------------------------------------------------------------------
c     initialize the work area
c-----------------------------------------------------------------------
      coef=0.
      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
      do i=max(1,is_pe-1),min(imt-1,ie_pe+1)
        do j=max(1,js_pe-1),min(jmt-1,je_pe+1)
          ustuff(i,j) = dxu(i)*csu(j)*hr(i,j) / (c2dt*dyu(j))
          vstuff(i,j) = dyu(j)*hr(i,j) / (c2dt*dxu(i)*csu(j))
        end do
      end do
c
c---------------------------------------------------------------------
c     calculate 9 point coefficients
c---------------------------------------------------------------------
c
        do j1=0,1
            do j2=-1,0
              do j=js,je
      do i1=0,1
          do i2=-1,0
                do  i=is,ie
                  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 spflame_module
      implicit none

      real coef  (is_pe:ie_pe,js_pe:je_pe,-1:1,-1:1)
      real ustuff(is_pe-1:ie_pe+1,js_pe-1:je_pe+1)
      real vstuff(is_pe-1:ie_pe+1,js_pe-1:je_pe+1)

      real :: cddxu(0:1,0:1),   cddyu(0:1,0:1)
      real :: cddxt(-1:0,-1:0), cddyt(-1:0,-1:0)
      real :: h
      real, parameter :: p5=0.5
      integer :: i,j,i1,j1,i2,j2,is,ie,js,je
c
      js = max(2,js_pe); je = min(jmt-1,je_pe)
      is = max(2,is_pe); ie = min(imt-1,ie_pe)
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 j=js,je
        do i=is,ie
#ifdef partial_cell
          h = hup(i,j)
#else
          h = zw(max(1,kmu(i,j)))*umask(i,1,j)
#endif
          ustuff(i,j) = h*dyu(j)/(dxu(i)*csu(j))
          vstuff(i,j) = h*dxu(i)*csu(j)/dyu(j)
        end do
      end do
      call border_exchg(ustuff,1,1)
      call set_cyclic(  ustuff,1,1)
      call border_exchg(vstuff,1,1)
      call set_cyclic(  vstuff,1,1)
c
c     calculate divergence = ddx (ddx (ustuff)) + ddy( ddy (vstuff))
c
      do j1=0,1
        do j2=-1,0
          do j=js,je
            do i1=0,1
              do i2=-1,0
                do  i=is,ie
                  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 



      subroutine phase_velocities_psi
      use spflame_module
      implicit none
      real var,var1
      integer i,j
c
c-----------------------------------------------------------------------
c     calculate phase velocities for psi
c     Orlanski radiation condition (passive open boundary)   
c-----------------------------------------------------------------------
c
      if (my_blk_j == 1 .and. enable_obc_south
     &     .and. .not. prescribe_psi_obc_south ) then
       var=-dyu(3)/dt
       do 4030 i=is_pe,ie_pe
         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 (my_blk_j == n_pes_j .and. enable_obc_north
     &    .and. .not. prescribe_psi_obc_north) then
       var=dyu(jmt-2)/dt
       do 4031 i=is_pe,ie_pe
         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 (my_blk_i == 1 .and. enable_obc_west
     &    .and. .not. prescribe_psi_obc_west ) then
       do 4032 j=js_pe,je_pe
         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 (my_blk_i == n_pes_i .and. enable_obc_east
     &    .and. .not. prescribe_psi_obc_east ) then
       do 4033 j=js_pe,je_pe
         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
      use spflame_module
      implicit none
      integer is,ie,js,je,i,j
      real data
c
c-----------------------------------------------------------------------
c     add open boundary values to the solution (active open boundary)
c-----------------------------------------------------------------------
c
      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)
c
c-----------------------------------------------------------------------
c     passive Orlanski radiation condition (passive open boundary)   
c-----------------------------------------------------------------------
c
      if (my_blk_j == 1 .and. enable_obc_south
     &     .and. .not. prescribe_psi_obc_south ) then
        do i=is,ie
         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 (my_blk_j == n_pes_j .and. enable_obc_north
     &    .and. .not. prescribe_psi_obc_north) then
       do i=is,ie
         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 (my_blk_i == 1 .and. enable_obc_west
     &    .and. .not. prescribe_psi_obc_west ) then
       do j=js,je
         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 (my_blk_i == n_pes_i .and. enable_obc_east
     &    .and. .not. prescribe_psi_obc_east ) then
       do j=js,je
         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-----------------------------------------------------------------------
c prescribed streamfunction at open boundary
c (which is however not exact, in a strict mathematically sense
c  since the solver assumes zero boundary conditions for the
c  time tendency of psi )
c-----------------------------------------------------------------------
c

      if (my_blk_j == 1 .and. enable_obc_south
     &     .and. prescribe_psi_obc_south ) then
        do i=is,ie
           psi(i,1:2,1) = psi_wall_south(i)
        enddo
      endif
c
      if (my_blk_j == n_pes_j .and. enable_obc_north
     &    .and. prescribe_psi_obc_north) then
        do i=is,ie
           psi(i,jmt-1:jmt,1) = psi_wall_north(i)
        enddo
      endif
c
      if (my_blk_i == 1 .and. enable_obc_west
     &    .and. prescribe_psi_obc_west ) then
        do j=js,je
          psi(1:2,j,1) = psi_wall_west(j)
        enddo
      endif

      if (my_blk_i == n_pes_i .and. enable_obc_east
     &    .and. prescribe_psi_obc_east ) then
        do j=js,je
          psi(imt-1:imt,j,1) = psi_wall_east(j)
        enddo
      endif
      call adjust_islands_obc
      end subroutine addobcpsi


      subroutine adjust_islands_obc
      use spflame_module
      use island_module
      implicit none
c
c-----------------------------------------------------------------------
c     set southern land masses and perimeter points to southern values
c     of the meridional open boundaries
c-----------------------------------------------------------------------
c
      integer i,j,isle,n,ierr,ii,jj,nn,ifound
      logical connected
      real flow
      logical :: first = .true.
      integer, allocatable :: map_big(:,:)
      integer mnisle_obc,maxipp_obc
      parameter ( mnisle_obc=500,maxipp_obc=40000)
      integer,allocatable ::  iperm_all(:),jperm_all(:)
      integer,allocatable ::  iofs_all(:),nippts_all(:)
      integer, save ::  iperm_obc(maxipp_obc),jperm_obc(maxipp_obc)
      integer, save ::  iofs_obc(mnisle_obc),nippts_obc(mnisle_obc)
      integer, save :: nisle_obc

      if (first) then
       allocate ( map_big(imt,jmt) )
       allocate(  iperm_all(maxipp_obc),jperm_all(maxipp_obc) )
       allocate(  iofs_all(mnisle_obc),nippts_all(mnisle_obc))

       call isleperim (kmt_big, map_big, 
     &    iperm_all, jperm_all, iofs_all, nippts_all, nisle_obc,
     &    imt, jmt,  mnisle_obc, maxipp_obc,1,cyclic,
     &    .false.,.false.,.false.,.false.,.false.)
c
c    Get the island perimeter points on the actual PE
c
       nn = 0
       do isle=1,nisle_obc
        iofs_obc(isle) = nn
        do n=1,nippts_all(isle)
          i = iperm_all(iofs_all(isle)+n)
          j = jperm_all(iofs_all(isle)+n)
          if(i.ge.is_pe.and.i.le.ie_pe.and.
     +       j.ge.js_pe.and.j.le.je_pe) then
            nn = nn+1
            iperm_obc(nn) = i
            jperm_obc(nn) = j
          endif
        end do
        nippts_obc(isle) = nn - iofs_obc(isle)
       end do
       deallocate( map_big,iperm_all,jperm_all, iofs_all,nippts_all)
      endif

      if (enable_obc_east) then
c      see with which islands the obc is connected to
        do isle=1,nisle_obc
         connected=.false.

         do n=1,nippts_obc(isle)
          i = iperm_obc(iofs_obc(isle)+n)
          if (i==imt) connected=.true.
         end do

         call global_lor(connected)

         if (connected) then
c         determine the boundary value for each connected island
          flow=0.
          ifound=0
          do n=1,nippts_obc(isle)
           i = iperm_obc(iofs_obc(isle)+n)
           j = jperm_obc(iofs_obc(isle)+n)
           if (i==imt) then
             ifound=1
             flow=psi(i,j,1)
           endif
          end do

          call global_sum_int(ifound)
          call global_sum(flow)

          if (ifound>0) then
            flow=flow/(1.*ifound)
          else
            print*,' cannot find a boundary value for island',isle
            call halt_stop(' in adjust_islands_obc ')
          endif
c       set all connected island points which are not the obc
c       to the boundary value
          do n=1,nippts_obc(isle)
           i = iperm_obc(iofs_obc(isle)+n)
           j = jperm_obc(iofs_obc(isle)+n)
           psi(i,j,1)=flow
          enddo
         endif  ! connected
        end do ! isle
      endif ! enable_obc_east



      if (enable_obc_west) then
c      see with which islands the obc is connected to
        do isle=1,nisle_obc
         connected=.false.

         do n=1,nippts_obc(isle)
          i = iperm_obc(iofs_obc(isle)+n)
          if (i==1) connected=.true.
         end do

         call global_lor(connected)

         if (connected) then
c         determine the boundary value for each connected island
          flow=0.
          ifound=0
          do n=1,nippts_obc(isle)
           i = iperm_obc(iofs_obc(isle)+n)
           j = jperm_obc(iofs_obc(isle)+n)
           if (i==1) then
             ifound=1
             flow=psi(i,j,1)
           endif
          end do

          call global_sum_int(ifound)
          call global_sum(flow)

          if (ifound>0) then
            flow=flow/(1.*ifound)
          else
            print*,' cannot find a boundary value for island',isle
            call halt_stop(' in adjust_islands_obc ')
          endif
c       set all connected island points which are not the obc
c       to the boundary value
          do n=1,nippts_obc(isle)
           i = iperm_obc(iofs_obc(isle)+n)
           j = jperm_obc(iofs_obc(isle)+n)
           psi(i,j,1)=flow
          enddo
         endif  ! connected
        end do ! isle
      endif ! enable_obc_west


      if ( (enable_obc_south .and. .not. cyclic).or.
     &     (enable_obc_south .and. enable_obc_adjust_islands)
     &            ) then
c      see with which islands the obc is connected to
        do isle=1,nisle_obc
         connected=.false.

         do n=1,nippts_obc(isle)
          j = jperm_obc(iofs_obc(isle)+n)
          if (j==1) connected=.true.
         end do

         call global_lor(connected)

         if (connected) then
c         determine the boundary value for each connected island
          flow=0.
          ifound=0
          do n=1,nippts_obc(isle)
           i = iperm_obc(iofs_obc(isle)+n)
           j = jperm_obc(iofs_obc(isle)+n)
           if (j==1) then
             ifound=1
             flow=psi(i,j,1)
           endif
          end do

          call global_sum_int(ifound)
          call global_sum(flow)

          if (ifound>0) then
            flow=flow/(1.*ifound)
          else
            print*,' cannot find a boundary value for island',isle
            call halt_stop(' in adjust_islands_obc ')
          endif
c       set all connected island points which are not the obc
c       to the boundary value
          do n=1,nippts_obc(isle)
           i = iperm_obc(iofs_obc(isle)+n)
           j = jperm_obc(iofs_obc(isle)+n)
           psi(i,j,1)=flow
          enddo
         endif  ! connected
        end do ! isle
      endif ! enable_obc_souht



      if ((enable_obc_north .and. .not. cyclic ) .or.
     &    (enable_obc_north .and. enable_obc_adjust_islands)
     &   ) then
c      see with which islands the obc is connected to
        do isle=1,nisle_obc
         connected=.false.

         do n=1,nippts_obc(isle)
          j = jperm_obc(iofs_obc(isle)+n)
          if (j==jmt) connected=.true.
         end do

         call global_lor(connected)

         if (connected) then
c         determine the boundary value for each connected island
          flow=0.
          ifound=0
          do n=1,nippts_obc(isle)
           i = iperm_obc(iofs_obc(isle)+n)
           j = jperm_obc(iofs_obc(isle)+n)
           if (j==jmt) then
             ifound=1
             flow=psi(i,j,1)
           endif
          end do

          call global_sum_int(ifound)
          call global_sum(flow)

          if (ifound>0) then
            flow=flow/(1.*ifound)
          else
            print*,' cannot find a boundary value for island',isle
            call halt_stop(' in adjust_islands_obc ')
          endif
c       set all connected island points which are not the obc
c       to the boundary value
          do n=1,nippts_obc(isle)
           i = iperm_obc(iofs_obc(isle)+n)
           j = jperm_obc(iofs_obc(isle)+n)
           psi(i,j,1)=flow
          enddo
         endif  ! connected
        end do ! isle
      endif ! enable_obc_north


      first=.false.
      end subroutine adjust_islands_obc





      subroutine checkerboard (solution)
c-----------------------------------------------------------------------
c     removes "checkboard" null space from an array "solution"
c     author              c.h. goldberg    e-mail => chg@gfdl.gov
c-----------------------------------------------------------------------
      use spflame_module
      implicit none
      real    :: solution(is_pe-1:ie_pe+1,js_pe-1:je_pe+1)
      real    :: sum(0:1,0:1)
      integer :: noceansum(0:1,0:1)
      real    :: correction(0:1,0:1),c,redsum,blacksum,diff,x
      integer :: nred,nblack,i1,j1,i,j
c
      sum = 0.0; noceansum = 0
c
      do i1=0,1
        do j1=0,1
          do j=js_pe+j1,je_pe,2
            do i=is_pe+i1,ie_pe,2
              sum(i1,j1) = sum(i1,j1) + solution(i,j)
	      if (map(i,j) .le. 0) then
		noceansum(i1,j1) = noceansum(i1,j1) + 1
              end if
            end do
          end do
        end do
      end do
c
      do i1=0,1
        do j1=0,1
          call global_sum(sum(i1,j1))
          x=noceansum(i1,j1)
          call global_sum(x)
          noceansum(i1,j1)=int(x)
        enddo
      enddo

      redsum   = sum(0,0) + sum(1,1)
      blacksum = sum(1,0) + sum(0,1)
      nred     = noceansum(0,0) + noceansum(1,1)
      nblack   = noceansum(1,0) + noceansum(0,1)
      diff = redsum/nred - blacksum/nblack
      c    = diff / 2.0
c
      if (my_pe==0) then
      print *, ' '
      print '(a,i6,a,i6,a,e14.7)'
     &,         '=> checkerboard: nred = ',nred, ', nblack = ',nblack
     &,         ', removing a checkerboard correction of ', c
      endif
c
      correction (0,0) = -c; correction (1,1) = -c
      correction (1,0) =  c; correction (0,1) =  c
c
      do i1=0,1
        do j1=0,1
          do j=js_pe+j1,je_pe,2
            do i=is_pe+i1,ie_pe,2
              if (map(i,j) .le. 0) then
		solution(i,j) = solution(i,j) + correction(i1,j1)
              end if
            end do
          end do
        end do
      end do
      end subroutine checkerboard



      subroutine zero_level(surfpres, variable)
      use spflame_module
      implicit none
      real :: surfpres(is_pe-1:ie_pe+1,js_pe-1:je_pe+1)
      character (len=*) :: variable
      real :: sum=0., area_ocean=0., area, surfpres0
      integer :: i,j
c
c     this does not correctly handle multiple basins
c
      do i=is_pe,ie_pe
        do j=js_pe,je_pe
	  if (map(i,j) .le. 0) then
            area = dxt(i)*cst(j)*dyt(j)
            sum = sum + surfpres(i,j)*area
            area_ocean = area_ocean + area
	  end if
        end do
      end do
      call global_sum(sum) 
      call global_sum(area_ocean)
      surfpres0 = sum / area_ocean
      do i=is_pe,ie_pe
        do j=js_pe,je_pe
	  if (map(i,j) .le. 0) then
            surfpres(i,j) = surfpres(i,j) - surfpres0
	  end if
        end do
      end do
      if (my_pe==0) then
       print '(a,e14.7,a,a/)'
     & , '=> zero_level: removing a mean of ', surfpres0, ' from '
     & , variable
      endif
      end subroutine zero_level
