
c
c-----------------------------------------------------------------------
c     Island and Island Perimeter Mapping Routines
c     now contained in a module. Original code taken from GFDL MOM2.1.
c
c     Note: this version applies changes to the kmt mask
c   
c     SPFLAME version : c.eden
c                C.Eden Feb. 2002 (mailto ceden@phys.ocean.dal.ca)
c-----------------------------------------------------------------------
c
      module island_module

      public :: isleperim,showmap

      private :: mmm, expand, qinit, qpush,
     &           qpop, qempty, jn, js, ie, iw

      logical :: cyclic,change_kmt=.false.
      private :: cyclic,change_kmt

      contains
 
      subroutine isleperim(kmt,map, iperm, jperm, iofs, 
     &                     nippts, nisle
     &,       imt, jmt, mnisle, maxipp,my_pe,cyclic_bc,
     &        obc_north,obc_south,obc_west,obc_east,change_kmt_in)

c
c=======================================================================
c         Island and Island Perimeter Mapping Routines
c         December 1993/revised February 1995
c=======================================================================
c
      common /qsize/ maxqsize
c
      integer kmt(imt,jmt), map(imt,jmt), iperm(maxipp)
      integer jperm(maxipp), nippts(mnisle), iofs(mnisle)
      integer nisle,imt,jmt,mnisle,maxipp
      logical cyclic_bc,obc_north,obc_south,obc_west,obc_east
      logical change_kmt_in
c
      parameter (maxq=10000)
c
      dimension iq(maxq)
      dimension jq(maxq)
      integer qfront, qback
      integer ocean,my_pe
c
      parameter (land=1, ocean=0)
      parameter (kmt_land=0, kmt_ocean=1)
c

      cyclic=cyclic_bc
      change_kmt = change_kmt_in

      if (my_pe==0) 
     &      print '(/,a)','Finding perimeters of all land masses'
      if (my_pe==0.and. cyclic) 
     &      print '(a,/)',' using cyclic boundary conditions'
      if (my_pe==0.and. change_kmt) 
     &      print '(a,/)',' applying changes to kmt mask    '
c
c     initialize number of changes to kmt
c
      nchanges = 0
c
    1 continue
c
c-----------------------------------------------------------------------
c     copy kmt to map changing notation
c     initially, 0 means ocean and 1 means unassigned land in map
c     as land masses are found, they are labeled 2, 3, 4, ...,
c     and their perimeter ocean cells -2, -3, -4, ...,
c     when no land points remain unassigned, land mass numbers are
c     reduced by 1 and their perimeter ocean points relabelled accordingly
c-----------------------------------------------------------------------
c
      do i=1,imt
        do j=1,jmt
          if (kmt(i,j) .gt. 0) then
            map(i,j) = ocean
          else
            map(i,j) = land
          end if
        end do
      end do

c-----------------------------------------------------------------------
c     first row is land for calculating poisson equation
c-----------------------------------------------------------------------
      if (obc_south) map(:,1)     = land
c-----------------------------------------------------------------------
c     last row is land for calculating poisson equation
c-----------------------------------------------------------------------
      if (obc_north) map(:,jmt)   = land
c-----------------------------------------------------------------------
c     first column is land for calculating poisson equation
c-----------------------------------------------------------------------
      if (obc_west)  map(1,:)     = land
c-----------------------------------------------------------------------
c     last column is land for calculating poisson equation
c-----------------------------------------------------------------------
      if (obc_east)  map(imt,:)   = land
c
c-----------------------------------------------------------------------
c     find unassigned land points and expand them to continents
c-----------------------------------------------------------------------
c
      maxqsize = 0
      call qinit (iq, jq, qfront, qback)
      label = 2
      iofs(label) = 0
      nippts(label) = 0
      nerror = 0
      jnorth = jmt
      if (cyclic) then
       iwest = 2
       ieast = imt-1
      else
       iwest = 1
       ieast = imt
      endif
      do j=jnorth,1,-1
        do i=iwest,ieast
          if (map(i,j) .eq. land) then
            call qpush (i, j, iq, jq, qfront, qback)
            call expand (map, label, iq, jq, qfront, qback, nerror
     &,                  iperm, jperm, iofs, nippts
     &,                  imt, jmt, mnisle, maxipp, kmt)
           if (my_pe==0)
     &      print '(a,i2,a,i6)',
     &        'number of island perimeter points: nippts(',label-1,')=',
     &         nippts(label)
            label = label + 1
            if (label .gt. mnisle) then
              print '(a,i3,a)','ERROR==> mnisle=',mnisle,' is too small'
              stop '==> expand'
            end if
            iofs(label) = iofs(label-1) + nippts(label-1)
            nippts(label) = 0
          end if
        end do
      end do
      nisle = label - 1
c
c-----------------------------------------------------------------------
c     relabel land masses and their ocean perimeters
c------------------------------------------------------------------------
c
      do i=iwest,ieast
        do j=1,jnorth
          if (map(i,j) .ne. 0) then
            map(i,j) = map(i,j) - sign(1, map(i,j))
          end if
        end do
      end do
      do isle=2,nisle
        iofs(isle-1) = iofs(isle)
        nippts(isle-1) = nippts(isle)
      end do
      nisle = nisle - 1

      if (cyclic) then
       do j=1,jmt
         map(1,j) = map(imt-1,j)
         map(imt,j) = map(2,j)
       end do
      endif
c
c      print '(/a)',
c     &       'The following changes fix "PERIMETER VIOLATIONS"'
c
      if (nerror .gt. 0) then
       if (change_kmt) then
        if (my_pe ==0) print *
     &, '==>Remapping land masses to see if PERIMETER VIOLATIONS remain'
        if (my_pe ==0) print *, ' '
        goto 1
       endif
       if (my_pe==0) call showmap(map,imt,jmt)
       call barrier
       call halt_stop(' in isleperim')
      end if
c
      if (my_pe==0) then
       print *,' Island perimeter statistics:'
       print *,'maximum queue size was ',maxqsize
       print *,'number of land masses is ', nisle
       print *,'number of island perimeter points is ',
     &         nippts(nisle) + iofs(nisle)
      endif
c

      return
      end subroutine isleperim


      subroutine showmap (map, imt, jmt)
      dimension map(imt,jmt)
      linewidth = 125
      print '(/,132a)',(' ',l=1,5+min(linewidth,imt)/2-13)
     &,     'Land Masses and Perimeters'
      istart = 0
      iremain = imt
      do isweep=1,imt/linewidth + 1
        iline = min(iremain, linewidth)
        iremain = iremain - iline
        if (iline .gt. 0) then
          print *, ' '
          print '(t6,32i5)', (istart+i+4,i=1,iline,5)
          do j=jmt,1,-1
            print '(i4,t6,160i1)', j,(mmm(map(istart+i,j)),i=1,iline)
          end do
          print '(t6,32i5)', (istart+i+4,i=1,iline,5)
          istart = istart + iline
        end if
      end do
      print *, ' '
      return
      end subroutine showmap



      function mmm(m)
      if (m .eq. 0) then
        mmm = 0
      else if (m .gt. 0) then
        mmm = mod(m,10)
      else
        mmm = m
      end if
      return
      end function mmm




      subroutine expand (map, label, iq, jq, qfront, qback, nerror
     &,                  iperm, jperm, iofs, nippts
     &,                  imt, jmt, mnisle, maxipp, kmt)
c
c-----------------------------------------------------------------------
c          The subroutine expand uses a "flood fill" algorithm
c          to expand one previously unmarked land
c          point to its entire connected land mass and its perimeter
c          ocean points.   Diagonally adjacent land points are
c          considered connected.  Perimeter "collisions" (i.e.,
c          ocean points that are adjacent to two unconnected
c          land masses) are detected and error messages generated.
c
c          The subroutine expand uses a queue of size maxq of
c          coordinate pairs of candidate points.  Suggested
c          size for maxq is 4*(imt+jmt).  Queue overflow stops
c          execution with a message to increase the size of maxq.
c          Similarly a map with more that maxipp island perimeter
c          points or more than mnisle land masses stops execution
c          with an appropriate error message.
c-----------------------------------------------------------------------
c
      dimension map(imt,jmt), kmt(imt,jmt)

      dimension iperm(maxipp)
      dimension jperm(maxipp)
      dimension nippts(mnisle)
      dimension iofs(mnisle)


      character * 32 problem

      parameter (maxq=10000)
      dimension iq(maxq)
      dimension jq(maxq)
      integer qfront, qback
c      logical qempty

      integer offmap, ocean
      parameter (offmap = -1)
      parameter (land = 1, ocean = 0)

      parameter (mnisle2=500)
      logical bridge_to(1:mnisle2)
c
c      print '(a,i3)', 'Exploring land mass ',label-1
c
      if (mnisle2 .lt. mnisle) then 
        print '(a,i4,a)', 
     & 'ERROR:  change parameter (mnisle2=',mnisle,') in isleperim.F'
        stop '==>isleperim'
      end if
c
      do isle=1,mnisle
        bridge_to(isle) = .false.
      end do
c
c-----------------------------------------------------------------------
c     main loop:
c        Pop a candidate point off the queue and process it.
c-----------------------------------------------------------------------
c
 1000 continue

      if (qempty (iq, jq, qfront, qback)) then
        call qinit (iq, jq, qfront, qback)
        return
      else
        call qpop (i, j, iq, jq, qfront, qback)
c
c       case: (i,j) is off the map
        if (i .eq. offmap .or. j .eq. offmap) then
          goto 1000
c
c       case: map(i,j) is already labeled for this land mass
        else if (map(i,j) .eq. label) then
          goto 1000
c
c       case: map(i,j) is an ocean perimeter point of this land mass
        else if (map(i,j) .eq. -label) then
          goto 1000
c
c       case: map(i,j) is an unassigned land point
        else if (map(i,j) .eq. land) then
          map(i,j) = label
c         print *, 'labeling ',i,j,' as ',label
          call qpush (i,         jn(j,jmt), iq, jq, qfront, qback)
          call qpush (ie(i,imt), jn(j,jmt), iq, jq, qfront, qback)
          call qpush (ie(i,imt), j,         iq, jq, qfront, qback)
          call qpush (ie(i,imt), js(j,jmt), iq, jq, qfront, qback)
          call qpush (i,         js(j,jmt), iq, jq, qfront, qback)
          call qpush (iw(i,imt), js(j,jmt), iq, jq, qfront, qback)
          call qpush (iw(i,imt), j,         iq, jq, qfront, qback)
          call qpush (iw(i,imt), jn(j,jmt), iq, jq, qfront, qback)
          goto 1000
c
c       case: map(i,j) is an ocean point adjacent to this land mass
        else if (map(i,j) .eq. ocean .or. map(i,j) .lt. 0) then

c         subcase: map(i,j) is a perimeter ocean point of another mass
          if (map(i,j) .lt. 0) then
            nerror = nerror + 1
            print '(a,a,i3,a,i3,a,a,i3,a,i3)',
     &            'PERIMETER VIOLATION==> ',
     &            'map(',i,',',j,') is in the perimeter of both ',
     &            'land masses ', -map(i,j)-1, ' and ', label-1
c           if we just quit processing this point here, problem points
c           will be flagged several times.
c           if we relabel them, then they are only flagged once, but
c           appear in both island perimeters, which causes problems in
c           island integrals.  current choice is quit processing.
c
c           only fill first common perimeter point detected.
c           after the first land bridge is built, subsequent collisions
c           are not problems.
c
            if (.not. bridge_to(-map(i,j)-1)) then
c             option 1: fill common perimeter point to make land bridge
                  kmt(i,j)= 0
                  bridge_to(-map(i,j)-1) = .true.

c              do n=1,nchanges
c                if (kmt_changes(n,1) .eq.i .and.
c     &              kmt_changes(n,2) .eq.j .and.
c     &              kmt_changes(n,4) .eq.0) then
c                  bridge_to(-map(i,j)-1) = .true.
c                end if
c              end do

            end if
c
            goto 1000
          end if

c         case: map(i,j) is a ocean point--label it for current mass
          map(i,j) = -label
          nippts(label) = nippts(label) + 1
c         print *, 'iofs(label)=',iofs(label)
c         print *, 'nippts(label)=',nippts(label)
          if (iofs(label) + nippts(label) .gt. maxipp) then
            print *, 'ERROR==>  maxipp=',maxipp,' is not large enough'
            stop '==>expand'
          end if
          iperm(iofs(label) + nippts(label)) = i
          jperm(iofs(label) + nippts(label)) = j
          goto 1000
c
c       case: map(i,j) is probably labeled for another land mass
c       ************* this case should not happen **************
        else
          nerror = nerror + 1
          print '(a,a,i3,a,i3,a,a,i3,a,i3)',
     &          'ERROR ==>  ',
     &          'map(',i,',',j,') is labeled for both ',
     &          'land masses ', map(i,j)-1,' and ',label-1
        end if
        goto 1000
c
      end if

      return
      end subroutine expand


      subroutine qinit (iq, jq, qfront, qback)
      parameter (maxq=10000)
      dimension iq(maxq)
      dimension jq(maxq)
      integer qfront, qback
c
      qfront = 1
      qback = 0
c
c     fake assignments to iq and jq to avoid "flint" warning
c
      iq(qfront) = 0
      jq(qfront) = 0
c
      return
      end subroutine qinit


      subroutine qpush (i, j, iq, jq, qfront, qback)
c
      common /qsize/ maxqsize
c
      parameter (maxq=10000)
      dimension iq(maxq)
      dimension jq(maxq)
      integer qfront, qback
c
      qback = qback + 1
      if (qback .gt. maxq) then
        if (qfront .ne. 1) then
c         shift queue left to make room
          ishift = qfront - 1
          do ip=qfront,qback-1
            iq(ip-ishift) = iq(ip)
            jq(ip-ishift) = jq(ip)
          end do
          qfront = 1
          qback = qback - ishift
        else
          print *, 'queue fault in qpush'
          stop '==>qpush'
        end if
      end if
      iq(qback) = i
      jq(qback) = j
c
      maxqsize = max(maxqsize, (qback-qfront))
c
      return
      end subroutine qpush



      subroutine qpop (i, j, iq, jq, qfront, qback)
c
      parameter (maxq=10000)
      dimension iq(maxq)
      dimension jq(maxq)
      integer qfront, qback
c
      i = iq(qfront)
      j = jq(qfront)
      qfront = qfront + 1
      return
      end subroutine qpop


      function qempty (iq, jq, qfront, qback)
c
      parameter (maxq=10000)
      dimension iq(maxq)
      dimension jq(maxq)
      integer qfront, qback
      logical qempty
c
      qempty = (qfront .gt. qback)
c      
      return
      end function qempty



      function jn(j,jmt)

c     j coordinate to the north of j

      integer offmap
      parameter (offmap = -1)

      if (j .lt. jmt) then
        jn = j + 1
      else
        jn = offmap
      end if
      return
      end function jn


      function js(j,jmt)

c     j coordinate to the south of j

      integer offmap
      parameter (offmap = -1)

      if (j .gt. 1) then
        js = j - 1
      else
        js = offmap
      end if

      return
      end function js



      function ie(i,imt)

c     i coordinate to the east of i

      integer offmap
      parameter (offmap = -1)

      if (cyclic) then
       if (i .lt. imt-1) then
         ie = i + 1
       else
         ie = (i+1) - imt + 2
       end if
      else
       if (i .lt. imt) then
        ie = i + 1
       else
        ie = offmap
       end if
      end if
      return
      end function ie



      function iw(i,imt)

c     i coordinate to the west of i

      integer offmap
      parameter (offmap = -1)

      if (cyclic) then
       if (i .gt. 2) then
        iw = i - 1
       else
        iw = (i-1) + imt - 2
       end if
      else
       if (i .gt. 1) then
        iw = i - 1
       else
        iw = offmap
       end if
      end if
      return
      end function iw


      end module island_module
