C    <@(#)	 mapping.F	4.8 3/28/94>
            
      subroutine MAPPING(len,im,jm,Depth,minimum,nwater,nland,
     &                   mapi,gridi,gridj,ieast,iwest,inorth,isouth)
      implicit none
      integer len
      integer im
      integer jm
      

      real    Depth(im,jm)
      real    minimum
      integer nwater
      integer nland
      integer mapi(len)
      integer gridi(len)
      integer gridj(len)
      integer ieast(len)
      integer iwest(len)
      integer inorth(len)
      integer isouth(len)
      integer i
      integer j
      integer n_open
      integer xx(IM*JM)
      integer yy(IM*JM)
      real    test(IM*JM)

      
      
      integer j1s
      integer j1n
      
c Functions

      integer INTMIN
      integer INTMAX

c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

      print *,'# '
      print *,'#   --------------------------------------------------'
      print *,'#    MAPPING    Depth(50,50) = ', Depth(50,50)
      print *,'# '
      
      
      do i = 1 , IM * JM
       xx(i) = MOD( i-1 , IM )
       yy(i) = (i-1)/IM
      end do
  

      
      call WHENFGT(IM*JM , Depth , 1 , minimum , mapi , nwater)
      
      

      
      nland = IM
      print *,'# '
      print *,'#     nwater     = ',nwater
      print *,'#     nland      = ',nland
      
      if (nwater+nland.gt.LEN) then
        print *,'#   nwater+nland greater than LEN = ',LEN
        stop
      endif
      
      
      do i = 1,nwater
       gridi(i) = xx(mapi(i))
      enddo
      print *,'#   gridi(1)=',gridi(1),'  gridi(nwater)=',gridi(nwater)
      do i = 1,nwater
       gridj(i) = yy(mapi(i))
      enddo
      print *,'#   gridj(1)=',gridj(1),'  gridj(nwater)=',gridj(nwater)
      
      do i = nwater+1 , LEN
        mapi(i) = IM*JM
        gridi(i) = 0
        gridj(i) = 0
        ieast(i) = i
        iwest(i) = i
        inorth(i) = i
        isouth(i) = i
      end do
      
      DO i = 1 , nwater
        j = nwater + 1 + MOD( i, nland )
        ieast(i)  = j
        iwest(i)  = j
        isouth(i) = j
        inorth(i) = j
      enddo
      

      call FIND_NEIGHBORSX( nwater, 
     &             ieast,  iwest,  .TRUE. )
     
      call FIND_NEIGHBORSY( nwater, 
     &             inorth, isouth, .FALSE. )

     
    
      print *,'#ieast(1)=',ieast(1),'   ieast(nwater)=',ieast(nwater)
      print *,'#iwest(1)=',iwest(1),'   iwest(nwater)=',iwest(nwater)
      print *,'#inorth(1)=',inorth(1),' inorth(nwater)=',inorth(nwater)
      print *,'#isouth(1)=',isouth(1),' isouth(nwater)=',isouth(nwater)
  
      
      print *,'# '
      print *,'#   --------------------------------------------------'
      print *,'# '

      end
      
c#undef LAND

C--------------------------------------------------------------
      
      subroutine MAPDOWN(A)
#include "gloparam.F.h"
#include "radpth.h"
#include "mapping.h"

      real A(IM*JM)
      

      real Wk1(LEN)
      integer i

      call GATHER(nwater,Wk1,A,mapi)

      do i=1,nwater
        A(i) = Wk1(i)
      end do
      do i=nwater+1,nwater+nland
        A(i) = 0.
      end do
      
      end  
C--------------------------------------------------------------
      
      subroutine MAPCMP( nu, nc, U, C, indx)
      implicit none
      integer nu
      integer nc
      real    U(nu)     ! Uncompressed data array
      real    C(nc)     ! compressed data
      integer indx(nc)  ! mapping index
      
      call GATHER(nc, C, U, indx)

      end
           
C--------------------------------------------------------------
      subroutine MAPEXP(nc, nu,C,U,indx)
      implicit none
      integer nc
      integer nu
      real    U(nu)
      real    C(nc)
      integer indx(nc)
     
      call SCATTER(nc,U,indx,C)

      end
           
C--------------------------------------------------------------
      subroutine FIND_NEIGHBORSX( npts , 
     &             iplus, iminus , cyclic )
#include "gloparam.F.h"
#include "radpth.h"
#include "mapping.h"
      integer  iplus( npts )
      integer  iminus( npts )
      logical  cyclic
      
      integer row , i, itest, row_length
      integer indices( im ), indx, indx_test, ih, it, ip1, im1
      

 
      do row = 1,jm
      
        call WHENIEQ( npts, gridj, 1 , row-1 , indices , row_length )
        
        if ( row_length .GT. 0 ) then

           do i = 1 , row_length 
           
             indx = indices( i )
             ih   = gridi( indx )
             
             if ( cyclic ) then
                ip1 = MOD( ih + 1 , im )
                im1 = MOD( ih + im - 1 , im )
             else
                ip1 = ih + 1
                im1 = ih - 1
             endif
             
             do itest = 1 , row_length
               indx_test = indices(itest)
               it = gridi( indx_test )
               if ( it .EQ. ip1 ) iplus( indx )  = indx_test
               if ( it .EQ. im1 ) iminus( indx ) = indx_test
             enddo ! test
             
           enddo ! i
        endif ! n
        
      enddo ! row
      
      end
C--------------------------------------------------------------
      subroutine FIND_NEIGHBORSY( npts , 
     &             iplus, iminus , cyclic )
#include "gloparam.F.h"
#include "radpth.h"
#include "mapping.h"
      integer  iplus( npts )
      integer  iminus( npts )
      logical  cyclic
      
      integer row , i, itest, row_length
      integer indices( jm ), indx, indx_test, ih, it, ip1, im1
      

 
      do row = 1,im
      
        call WHENIEQ( npts, gridi, 1 , row-1 , indices , row_length )
        
        if ( row_length .GT. 0 ) then

           do i = 1 , row_length 
           
             indx = indices( i )
             ih   = gridj( indx )
             
             if ( cyclic ) then
                ip1 = MOD( ih + 1 , jm )
                im1 = MOD( ih + jm - 1 , jm )
             else
                ip1 = ih + 1
                im1 = ih - 1
             endif
             
             do itest = 1 , row_length
               indx_test = indices(itest)
               it = gridj( indx_test )
               if ( it .EQ. ip1 ) iplus( indx )  = indx_test
               if ( it .EQ. im1 ) iminus( indx ) = indx_test
             enddo ! test
             
           enddo ! i
        endif ! n
        
      enddo ! row
      
      end