program setregionmask

!-------------------------------------------------------------------
! PURPOSE:
! o Reads POP input binaries kmt.da & grid.pop.da for a given grid,
!   and writes out a region.ieeei4 file.
! o For modern topography, generates full region set; for paleo
!   topography, generates north/south mask.
! o USER must specify boundary segment start/end points which, 
!	together with land, will completely box in a region.  A
!	boundary set for one region also serves as a boundary
!	for the region immediately adjacent.
!   AND
!	must provide a seed point (i_region,j_region) from which 
!	a given region will be filled.
!
!-------------------------------------------------------------------
!
!	Compile on dataproc with f90 setregionmask.f90
!
!-------------------------------------------------------------------

implicit none

!--- local ---

   integer :: imt, jmt, i, j, k, recl_factor,nb,nbad
   integer :: imiss=99
   logical :: paleo
   integer*4, allocatable :: kmt(:,:), mask(:,:), cnt(:,:)
   real*8, allocatable :: ulat(:,:),ulng(:,:),tlat(:,:),tlng(:,:)
   real :: atlantic(7,4),pacific(5,4),indian(1,4)
   real :: north(4),northsea(2,4),labrador(2,4)
   character(len=80) :: gridf, kmtf, maskf

!  line segments which demarc region boundaries 
! 		        /x0,  y0,  x1,  y1/  in degrees
!
   data north(1:4)      /0., 0., 360., 0. /	  !N. Hemisphere
   data atlantic(1,1:4) /300.,-32., 22., -32. /	  !S. Atlantic
   data atlantic(2,1:4) /353., 32., 353., 39. /	  !Med.
   data atlantic(3,1:4) /12., 51., 12., 60. /	  !Baltic
   data atlantic(4,1:4) /8., 60., 341., 65. /	  !Finland-Iceland
   data atlantic(5,1:4) /341., 65., 325., 72. /	  !Iceland-Greenland
   data atlantic(6,1:4) /315., 62., 315., 50. /	  !Greenland-Atlantic
   data atlantic(7,1:4) /315., 50., 300., 50. /	  !Atlantic-NovaScotia
   data pacific(1,1:4) /143.,-32., 300., -32. /	  !S. Pacific
   data pacific(2,1:4) /142.,-20., 142., -6. /	  !Torres Strait (?)
   data pacific(3,1:4) /130.,-1.5, 116., -1.5 /	  !Borneo Strait (?)
   data pacific(4,1:4) /98.5,-1.8, 112., -1.8 /	  !Java Strait (?)
   data pacific(5,1:4) /180.,66., 203., 66. /	  !Bering Strait
   data indian(1,1:4) /26.,-32., 117., -32. /	  !S. Pacific
   data northsea(1,1:4) /26.,70., 16.6, 79. / 
   data northsea(2,1:4) /16.6,79., 336., 81. / 
   data labrador(1,1:4) /286.5,61.,290.5,64./ 
   data labrador(2,1:4) /279.,73.7,275.,78./ 

   print *, 'Enter grid dimensions : imt, jmt'
   read *, imt, jmt
   print *, 'Enter name of POP grid binary:'
   read *, gridf
   print *, 'Enter name of POP kmt binary:'
   read *, kmtf
   print *, 'Enter name of POP region mask binary (output):'
   read *, maskf
   print *, 'Is this a paleo grid (T/F)?'
   read *, paleo

   allocate(kmt(imt,jmt))
   allocate(mask(imt,jmt))
   allocate(cnt(imt,jmt))
   allocate(ulat(imt,jmt))
   allocate(ulng(imt,jmt))
   allocate(tlat(imt,jmt))
   allocate(tlng(imt,jmt))

   open (1, file=kmtf, access='direct', &
     form='unformatted', recl=jmt*imt*4, status='old')
   open (2, file=gridf, access='direct', &
     form='unformatted', recl=jmt*imt*8, status='old')
   open (3, file=maskf, access='direct', &
     form='unformatted', recl=jmt*imt*4, status='unknown', &
     action='write')

! read in kmt and grid:
   read  (1, rec=1) kmt
   read  (2, rec=1) ulat 
   read  (2, rec=2) ulng
   close(1)
   close(2)

! get T grid locations
   call calc_tpoints(imt,jmt,ulng,ulat,tlng,tlat)
     
! mask 99 for ocean, 0 for land
   mask = imiss
   where(kmt == 0) mask=0
   cnt = 0

! set region boundaries (where necessary), then fill
   if (paleo) then
!	For paleo grids, presumably do something simple like
!	Northern Hemisphere & Southern Hemisphere
     nb = size(north,1)
     call set_bounds(imt,jmt,tlng,tlat,mask,nb,north,1,'north')

     call neighbor_fill(imt,jmt,mask,1,3,87,imiss)  ! fill North
     call neighbor_fill(imt,jmt,mask,2,3,40,imiss)  ! fill South
   else
     nb = size(atlantic,1)
     call set_bounds(imt,jmt,tlng,tlat,mask,nb,atlantic,6,'atlantic')
     nb = size(pacific,1)
     call set_bounds(imt,jmt,tlng,tlat,mask,nb,pacific,2,'pacific')
     nb = size(indian,1)
     call set_bounds(imt,jmt,tlng,tlat,mask,nb,indian,3,'indian')
     nb = size(northsea,1)
     call set_bounds(imt,jmt,tlng,tlat,mask,nb,northsea,9,'northsea')
     nb = size(labrador,1)
     call set_bounds(imt,jmt,tlng,tlat,mask,nb,labrador,8,'labrador')

!  Filling regions from seed points.  Note that some regions require
!	more than one fill because they are not contiguous
     call neighbor_fill(imt,jmt,mask,6,2,57,imiss) ! fill Atlantic
     call neighbor_fill(imt,jmt,mask,2,66,56,imiss) ! fill Pacific
     call neighbor_fill(imt,jmt,mask,3,29,48,imiss) ! fill Indian
     call neighbor_fill(imt,jmt,mask,1,34,18,imiss) ! fill Southern Ocean
     call neighbor_fill(imt,jmt,mask,8,92,96,imiss) ! fill Labrador
     call neighbor_fill(imt,jmt,mask,9,19,107,imiss) ! fill North Sea
     call neighbor_fill(imt,jmt,mask,10,53,96,imiss) ! fill Arctic
     call neighbor_fill(imt,jmt,mask,7,16,76,imiss) ! fill Med
     call neighbor_fill(imt,jmt,mask,7,11,76,imiss) ! fill Med
     call neighbor_fill(imt,jmt,mask,-13,20,78,imiss) ! fill Black
     call neighbor_fill(imt,jmt,mask,-14,26,79,imiss) ! fill Caspian
     call neighbor_fill(imt,jmt,mask,-5,22,70,imiss) ! fill Red
     call neighbor_fill(imt,jmt,mask,-5,21,72,imiss) ! fill Red
     call neighbor_fill(imt,jmt,mask,-5,23,67,imiss) ! fill Red
     call neighbor_fill(imt,jmt,mask,4,27,71,imiss) ! fill Persian
     call neighbor_fill(imt,jmt,mask,4,25,72,imiss) ! fill Persian
     call neighbor_fill(imt,jmt,mask,-12,19,87,imiss) ! fill Baltic
     call neighbor_fill(imt,jmt,mask,-12,20,90,imiss) ! fill Baltic
     call neighbor_fill(imt,jmt,mask,11,86,91,imiss) ! fill Hudson

   end if 

   where(mask==imiss) cnt=1
   nbad = sum(cnt)
   if (nbad /= 0) then
      write(6,*) 'Still have ',nbad,' points with unassigned value'
   end if

   write (unit=3, rec=1)mask
   write(6,*) ' '
   write(6,*) 'DONE writing ',maskf

end program setregionmask

!=================================================================

recursive SUBROUTINE neighbor_fill(imt,jmt,mask,value,i0,j0,miss)

!-----------------------------------------------------------------------
!  Recursively fills a region starting from seed point (i0,j0), and
!  filling with "value" wherever north/south/east/west shifts indicate
!  a missing value.
!-----------------------------------------------------------------------

   integer :: i0, j0, imt, jmt, miss, value
   integer :: ie,iw,jn,js
   integer, dimension(imt,jmt) :: mask

   if (mask(i0,j0) == miss) then 
      mask(i0,j0)=value
      iw = iadd(i0,-1,imt)
      call neighbor_fill(imt,jmt,mask,value,iw,j0,miss)
      ie = iadd(i0,1,imt)
      call neighbor_fill(imt,jmt,mask,value,ie,j0,miss)
      jn = iadd(j0,1,jmt)
      call neighbor_fill(imt,jmt,mask,value,i0,jn,miss)
      js = iadd(j0,-1,jmt)
      call neighbor_fill(imt,jmt,mask,value,i0,js,miss)
   end if

END SUBROUTINE neighbor_fill

!=================================================================

SUBROUTINE set_bounds(imt,jmt,tlng,tlat,mask,nb,bounds,value,rstr)

!-----------------------------------------------------------------------
!
!  Sets MASK array to VALUE along linear segments specified in the
!  BOUNDS array (there are NB segments to process for a given region).
!  The line segments are meant to follow great circles which intersect
!  the endpoints in BOUNDS, but in practice, a slope is computed in
!  grid i,j space, and this slope is traced from a mesh starting cell
!  to a mesh ending cell.  
!
!-----------------------------------------------------------------------

   integer :: i, j, k, imt, jmt, nb, value, kcnt, lcnt
   integer :: sti,stj,eni,enj,thisi,thisj
   integer, allocatable :: movei(:),movej(:)
   integer, dimension(imt,jmt) :: mask, iind, jind
   real, dimension(nb,4) :: bounds
   real :: x0,x1,y0,y1
   real*8, dimension(imt,jmt) :: tlng,tlat,dist0,dist1
   real*8  :: mindist0,mindist1
   character(len=*) :: rstr

   rstr = trim(rstr)

   do i=1,imt 
      iind(i,1:jmt)=i
   end do
   do j=1,jmt 
      jind(1:imt,j)=j
   end do

!  loop over number of boundary segments for this region
   do k=1,nb
     sti=0
     stj=0
     eni=0
     enj=0
     lcnt=0
     kcnt=0
     x0 = bounds(k,1)
     y0 = bounds(k,2)
     x1 = bounds(k,3)
     y1 = bounds(k,4)
     dist0 = sqrt((tlng-x0)**2+(tlat-y0)**2)
     dist1 = sqrt((tlng-x1)**2+(tlat-y1)**2)
     mindist0 = 1.e30
     mindist1 = 1.e30
     do j=1,jmt
     do i=1,imt
       if (dist0(i,j)<mindist0) then
         mindist0 = dist0(i,j)
       end if
       if (dist1(i,j)<mindist1) then
         mindist1 = dist1(i,j)
       end if
     end do
     end do

     do j=1,jmt
     do i=1,imt
!    write(6,*) 'i,j,dist0,dist1,mindist0,mindist1',i,j,dist0(i,j), &
! &      dist1(i,j),mindist0,mindist1
       if (dist0(i,j)==mindist0) then
         sti = i
         stj = j
       end if
       if (dist1(i,j)==mindist1) then
         eni = i
         enj = j
       end if
     end do
     end do
     if (min(sti,stj,eni,enj)==0) then
       write(6,*) 'sti,stj,eni,enj = ',sti,stj,eni,enj
       stop 'Did not locate segment endpoint on mesh'
     end if
!    write(6,*) ' '
!    write(6,*) 'processing segment ',k,' of ',rstr
!    write(6,1000) 'starting point: ',x0,y0,'  model grid: ', &
!   &    tlng(sti,stj),tlat(sti,stj)
!    write(6,1000) 'ending point: ',x1,y1,'  model grid: ', &
!   &    tlng(eni,enj),tlat(eni,enj)
1000  format(a16,f7.3,2x,f7.3,a14,f7.3,2x,f7.3)

!   get slope (in units of grid cell steps) from start to end of segment
     idx = eni-sti
     idy = enj-stj
     if (abs(idx) > imt/2) then
       if (idx > 0) idx = idx-imt
       if (idx < 0) idx = idx+imt
     end if
     nidx = abs(idx)
     nidy = abs(idy)

!    factor out common denominators in slope idy/idx
10   continue
     do i=2,min(nidx,nidy)
       if (mod(idx,i)==0 .and. mod(idy,i)==0) then
          idx = idx/i
          idy = idy/i
          goto 10
       end if
     end do

!   now that slope is reduced, determine best move sequence
!    write(6,*) 'using slope DY = ',idy,' DX = ',idx
     nmove = nidx+nidy
     allocate(movei(nmove))
     allocate(movej(nmove))
     call get_moves(nmove,idx,idy,movei,movej)
!    write(6,*) 'movei = ',(movei(i),i=1,nmove)
!    write(6,*) 'movej = ',(movej(i),i=1,nmove)

!    step from starting point on mesh (sti,stj) to adjacent cells until
!    you get to the grid cell nearest the segment end (eni,enj)
!    Set ocean point mask to "value", leave land point mask values 0
     thisi = sti
     thisj = stj
     if (mask(thisi,thisj) /= 0) mask(thisi,thisj) = value
     do
       do i=1,nmove
	thisi=iadd(thisi,movei(i),imt)
	thisj=jadd(thisj,movej(i),jmt)
        kcnt = kcnt + 1
!       write(6,*) thisi,thisj,mask(thisi,thisj)
        if (mask(thisi,thisj) /= 0) then
	   mask(thisi,thisj) = value
	   lcnt = lcnt+1
        end if
        if (thisi==eni .and. thisj==enj) exit
       end do
       if (thisi==eni .and. thisj==enj) exit
       if (kcnt > 100) stop 'not getting anywhere'
     end do
!    write(6,*) 'assigned region value to ',lcnt,' cells'

     deallocate(movei)
     deallocate(movej)
   end do	! k, loop over nb boundary segments

END SUBROUTINE set_bounds

!***********************************************************************

integer FUNCTION iadd(i,di,ni)

   implicit none

   !--- arguments ---
   integer :: i, di, ni

   if (di < 0) then
     iadd = mod(i+di-1+ni,ni)+1
   else
     iadd = mod(i,ni)+di
   endif

END FUNCTION iadd

!===============================================================================

integer FUNCTION jadd(j,dj,nj)

   implicit none

   !--- arguments ---
   integer :: j, dj, nj

   if (dj < 0) then
     jadd = max(j+dj,1)
   else
     jadd = min(j+dj,nj)
   endif

END FUNCTION jadd

!===============================================================================

SUBROUTINE get_moves(nmove,idx,idy,movei,movej)

!-----------------------------------------------------------------------
!
!     For a slope in i,j grid space given by idy & idx, this
!	routine determines a set of adjacent cell moves to
!	accomplish this slope in the most linear fashion possible.
!	The moves are stored as 1-D arrays movei, movej which
!	hold values of -1, 0, or 1.
!
!-----------------------------------------------------------------------

   integer ::  i, j , nmove, idx, idy, nidx, nidy, cnt
   integer, dimension(nmove) ::  tmp1, tmp2, movei, movej

   movei = 0
   movej = 0
   cnt = 0
   nidx = abs(idx)
   nidy = abs(idy)

!  cannot improve the linearity when either idx or idy <= 1
   if (min(nidx,nidy) <= 1) then
      do j=1,nidy
        cnt = cnt + 1
	movej(cnt)=idy/nidy
      end do
      do i=1,nidx
        cnt = cnt + 1
        movei(cnt)=idx/nidx
      end do
      return
   end if

!  
   if (nidy>nidx) then
      slp = nint(real(nidy)/real(nidx))
      nfull = min(int(nidy/slp),nidx)
      nxleft = nidx-nfull
      nyleft = nidy-(slp*nfull)
      cnt = 0
      do i=1,nfull
         cnt = cnt + 1
         movei(cnt) = idx/nidx
         do j=1,slp
           cnt = cnt + 1
           movej(cnt) = idy/nidy
         end do
      end do
      do i=1,nxleft
        cnt = cnt + 1
        movei(cnt) = idx/nidx
      end do
      do j=1,nyleft
         cnt = cnt + 1
         movej(cnt) = idy/nidy
      end do
   else			! nidx>nidy since they can't be equal
      slp = nint(real(nidx)/real(nidy))
      nfull = min(int(nidx/slp),nidy)
      nyleft = nidy-nfull
      nxleft = nidx-(slp*nfull)
      cnt = 0
      do i=1,nfull
         cnt = cnt + 1
         movej(cnt) = idy/nidy
         do j=1,slp
           cnt = cnt + 1
           movei(cnt) = idx/nidx
         end do
      end do
      do i=1,nyleft
        cnt = cnt + 1
        movej(cnt) = idy/nidy
      end do
      do j=1,nxleft
         cnt = cnt + 1
         movei(cnt) = idx/nidx
      end do
   end if

END SUBROUTINE get_moves

!***********************************************************************

SUBROUTINE calc_tpoints(imt,jmt,ux,uy,tlng,tlat)

!-----------------------------------------------------------------------
!
!     calculates lat/lon coordinates of T points from U points.
!     this computation is identical to the method used in the ccsm
!     ice model subroutine Tlatlon in module ice_gird.F as of
!     14 September 2001
!
!-----------------------------------------------------------------------

   integer ::  i, j , im1, imt, jmt
   real*8, parameter :: pi=3.14159265358979323846
   real*8 :: z1,x1,y1,z2,x2,y2,z3,x3,y3,z4,x4,y4,tx,ty,tz,da
   real*8, intent(in), dimension(imt,jmt) :: ux, uy
   real*8, intent(out), dimension(imt,jmt) :: tlat, tlng

   write(6,*) 'ux[50,50] = ',ux(50,50)*(180./pi)
   write(6,*) 'uy[50,50] = ',uy(50,50)*(180./pi)
      do j=2,jmt
        do i=1,imt

            if (i.eq.1) then
               im1=imt
            else
               im1=i-1
            endif

            z1 = cos(uy(im1,j-1))
            x1 = cos(ux(im1,j-1))*z1
            y1 = sin(ux(im1,j-1))*z1
            z1 = sin(uy(im1,j-1))

            z2 = cos(uy(i,j-1))
            x2 = cos(ux(i,j-1))*z2
            y2 = sin(ux(i,j-1))*z2
            z2 = sin(uy(i,j-1))

            z3 = cos(uy(im1,j))
            x3 = cos(ux(im1,j))*z3
            y3 = sin(ux(im1,j))*z3
            z3 = sin(uy(im1,j))

            z4 = cos(uy(i,j))
            x4 = cos(ux(i,j))*z4
            y4 = sin(ux(i,j))*z4
            z4 = sin(uy(i,j))

            tx = (x1+x2+x3+x4)/c4
            ty = (y1+y2+y3+y4)/c4
            tz = (z1+z2+z3+z4)/c4
            da = sqrt(tx**2+ty**2+tz**2)

            tz = tz/da
            ! tlng in radians
            tlng(i,j) = 0.0
            if (tx.ne.0.0.or.ty.ne.0.0) &
                 tlng(i,j) = atan2(ty,tx)
            ! tlat in radians
            tlat(i,j) = asin(tz)

        end do
      end do

      ! j=1: linear approximation
      do i=1,imt
         tlng(i,1) = tlng(i,2)
         tlat (i,1) = 2.0*tlat(i,2) - tlat(i,3)
      end do

      where (tlng > 2.*pi) tlng = tlng - 2.*pi
      where (tlng < 0.0 ) tlng = tlng + 2.*pi

! convert to degrees
    tlat = tlat*(180./pi)
    tlng = tlng*(180./pi)
   write(6,*) 'tx[50,50] = ',tlng(50,50)
   write(6,*) 'ty[50,50] = ',tlat(50,50)
!-----------------------------------------------------------------------

END SUBROUTINE calc_tpoints



