program crop
c
c --- crop depth array to eliminate small islands and unconnected water bodies
c
      include 'dimension_raw.h'
      real depth(idm,jdm),work(idm,jdm)
      integer mark(idm,jdm)
      character fmt2*9,char2*2,preambl(5)*79,util(idm*jdm+14)*2
      data fmt2/'(i4,19i4)'/
      data cutoff/15./
c
      common /conre1/ ioffp      ,spval
      common /conre4/ isizel     ,isizem     ,isizep     ,nrep       ,
     1                ncrt       ,ilab       ,nulbll     ,ioffd      ,
     2                ext        ,ioffm      ,isolid     ,nla        ,
     3                nlm        ,xlt        ,ybt        ,side
c
      ioffm=1
      nulbll=1
      isizel=0
c
      open (unit=11,
     &  file='~/your_depth_file',
     &  status='old')
      read (11,'(a79)') preambl
      write (*,'(1x,a79)') preambl
      read (11,'(3i8/(40a2))') idim,jdim,length,(util(l),l=1,length)
      call unpakk(depth,idm,idim,jdim,util,length)
      close (unit=11)
c
      do 11 j=1,jj1
      do 11 i=1,ii1
      if (depth(i,j).lt.cutoff) depth(i,j)=0.
      if (depth(i,j).gt.0.) depth(i,j)=max(depth(i,j),2.*cutoff)
 11   continue
c


c      call pakk(depth,idm,ii1,jj1,util,length)
c      open (unit=10,file='depth.1437x1678',status='unknown')
c      write (10,'(5(a79/),3i8/(40a2))') preambl,ii1,jj1,length,
c     .  (util(l),l=1,length)
c      iss=1
c      if(iss.eq.1) stop ' cut off '

      
c      call zebra(depth,idm,ii1,jj1)
      call bigrid(depth)
ccc      do 3 i=1,ii1
ccc 3    write (*,'('' i=''i5'' jfp,jlp=''5(1x,2i5)/(17x,5(1x,2i5)))') i,
ccc     . (jfp(i,l),jlp(i,l),l=1,jsp(i))
ccc      do 4 j=1,jj1
ccc 4    write (*,'('' j=''i5'' ifp,ilp=''5(1x,2i5)/(17x,5(1x,2i5)))') j,
ccc     . (ifp(j,l),ilp(j,l),l=1,isp(j))
c
      do 8 j=1,jj1
      do 8 i=1,ii1
 8    mark(i,j)=0
c
ccc      call prtmsk(ip,depth,work,idm,ii1,jj1,0.,1.,
ccc     .    'water depth - 1437-sq.pts')
c
c --- separate unwanted basins from main basin
c --- (Atlantic 0.08 degree grid)
c
c --- hudson bay
      do 21 i=69,118
 21   depth(i,417)=0.
c
c --- baltic
      do 22 j=1367,1373
 22   depth(191,j)=0.
c
c --- remove single-point islands
c
      do 14 j=1,jj1
      ja=max(  1,j-1)
      jb=min(jj1,j+1)
      do 14 i=1,ii1
      ia=max(  1,i-1)
      ib=min(ii1,i+1)
      if (depth(i,j).eq.0..and.
     .  min(depth(ia,j ),depth(ib,j ),depth(i ,ja),depth(i ,jb),
     .      depth(ia,jb),depth(ib,ja),depth(ia,ja),depth(ib,jb))
     . .gt.0.) then
      depth(i,j)=2.*cutoff
      write (*,'('' island at i,j ='',2i5,''  removed'')') i,j
      end if
 14   continue
c
c --- build up windward islands (lesser antilles)
c
      depth(874,458)=0.		!  dominica
      depth(875,458)=0.		!  dominica
      depth(884,462)=0.		!  martinique
      depth(885,462)=0.		!  martinique
      depth(897,464)=0.		!  st. lucia
      depth(898,464)=0.		!  st. lucia
      depth(907,461)=0.		!  st. vincent
      depth(908,461)=0.		!  st. vincent
      depth(913,459)=0.		!  grenadines
      depth(914,458)=0.		!  grenadines
      depth(906,481)=0.		!  barbados
      depth(907,482)=0.		!  barbados
c
      depth(825,209)=0.		!  cayman
      depth(825,210)=0.		!  cayman
      depth(850,177)=0.		!  swan
c
      depth(643,416)=0.         !  bermuda
      depth(643,417)=0.         !  bermuda
      depth(644,415)=0.         !  bermuda
      depth(644,416)=0.         !  bermuda
c      
      depth(686,215)=2.*cutoff	!  sever land bridge to fake offshore island
      depth(553,1420)=2.*cutoff	!  sever land bridge to Sicily
      depth(553,1421)=2.*cutoff	!  sever land bridge to Sicily
      depth(554,1420)=2.*cutoff	!  sever land bridge to Sicily
      depth(554,1421)=2.*cutoff	!  sever land bridge to Sicily
      depth(555,1420)=2.*cutoff	!  sever land bridge to Sicily
      depth(555,1421)=2.*cutoff	!  sever land bridge to Sicily
c
      depth(53,1311)=0.         !  isolate bay
      depth(54,1311)=0.         !  isolate bay
      mark(56,1314)=-1          !  remove bay
c
      depth(956,200)=0.         !  isolate bay
      depth(956,201)=0.         !  isolate bay
      depth(957,200)=0.         !  isolate bay
      mark(957,200)=-1          !  remove bay
c
c --- mark one point in each basin to be filled by "-1"
c --- mark one point on each island to be submerged by "+1"
c
      mark(  69, 416)=-1	!  hudson bay
      mark(  32,  83)=-1
      mark(   4, 297)=-1
      mark( 377, 338)=-1
      mark( 220, 220)=-1
      mark( 463, 241)=-1	!  great lakes
      mark(  98,1294)=-1	!  norwegian fjord
      mark( 116,1303)=-1	!  norwegian fjord
      mark( 144,1296)=-1	!  norwegian fjord
      mark(  97,1306)=-1	!  norwegian fjord
      mark( 123,1292)=-1	!  norwegian fjord
      mark(  37,1349)=-1	!  norwegian fjord
      mark(  97,1289)=-1	!  norwegian fjord
      mark( 106,1288)=-1	!  norwegian fjord
      mark( 192,1371)=-1	!  baltic
      mark( 178,1437)=-1	!  baltic
      mark( 203,1437)=-1	!  baltic
      mark( 248,1352)=-1	!  baltic
      mark( 253,1367)=-1	!  baltic
      mark(1000,   1)=-1	!  pacific
      mark( 609,1301)=-1	!  salt lake, e.algeria
      mark( 615,1310)=-1	!  salt lake, e.algeria
      mark( 808, 257)=-1	!  cuban lowlands
      mark( 946, 324)=-1	!  lake maracaibo
      mark(1104, 709)=-1	!  brazil
      mark(1105, 726)=-1	!  brazil
      mark(1058,1344)=-1	!  gabon
      mark(400,85)   =-1        !  interior lakes 
      mark(388,104)  =-1        !  interior lakes
      mark(403, 77)  =-1        !  interior lakes
      mark(390,115)  =-1        !  interior lakes
      mark(398,137)  =-1        !  interior lakes
      mark(398,152)  =-1        !  interior lakes
      mark(400,160)  =-1        !  interior lakes
      mark(406,153)  =-1        !  interior lakes
      mark(449,140)  =-1        !  interior lakes
      mark(406,149)  =-1        !  interior lakes
      mark(400,105)  =-1        !  interior lakes
      mark(395,108)  =-1        !  interior lakes
      mark(963,189)  =-1        !  interior lakes
      mark(508,311)  =-1        !  interior lakes
      mark(378,341)  =-1        !  interior lakes
      mark(950,330)  =-1        !  interior lakes
      mark(941,447)  =-1        !  interior lakes
      mark(436,413)  =-1        !  interior lakes
      mark(1091,669) =-1        !  interior lakes
      mark(760,1676) =-1        !  other lakes
      mark(700,1659) =-1        !  other lakes
      mark(694,1660) =-1        !  other lakes
      mark(573,1623) =-1        !  other lakes
      mark(480,1676) =-1        !  other lakes
      mark(424,1676) =-1        !  other lakes
      mark(512,1578) =-1        !  other lakes
      mark(552,1496) =-1        !  other lakes
      mark(552,1502) =-1        !  other lakes
      mark(545,1515) =-1        !  other lakes
      mark(537,1514) =-1        !  other lakes
      mark(561,1516) =-1        !  other lakes
      mark(551,1525) =-1        !  other lakes
      mark(556,1527) =-1        !  other lakes
      mark(189,1521) =-1        !  other lakes
      mark(179,1504) =-1        !  other lakes
      mark(471,1438) =-1        !  other lakes
      mark(442,1412) =-1        !  other lakes
      mark(263,1407) =-1        !  other lakes
      mark(261,1401) =-1        !  other lakes
      mark(191,1152) =-1        !  other lakes
      mark(617,1303) =-1        !  other lakes
      mark(678,1488) =-1        !  other lakes
      mark(675,1578) =-1        !  other lakes
      mark(674,1631) =-1        !  other lakes
      mark(652,1672) =-1        !  other lakes
      mark(179,1508) =-1        !  other lakes
      mark(532,1476) =-1        !  other lakes
      mark(436,1406) =-1        !  other lakes
      
c
      mark( 1029,605)=+1        !  off brazil
      mark( 469, 452)=+1	!  sable island neighbor
      mark( 449, 494)=+1	!  sable island neighbor
      mark( 448, 605)=+1	!  grand banks
c      mark( 521, 330)=+1	!  n.atlantic bight
c      mark( 563, 294)=+1	!  n.atlantic bight
c      mark( 569, 289)=+1	!  n.atlantic bight
c      mark( 578, 287)=+1	!  n.atlantic bight
c      mark( 587, 287)=+1	!  n.atlantic bight
c      mark( 657, 223)=+1	!  s.atlantic bight
c      mark( 667, 221)=+1	!  s.atlantic bight
c      mark( 687, 220)=+1	!  s.atlantic bight
c      mark( 618, 266)=+1	!  s.atlantic bight
c      mark( 820,  22)=+1	!  gulf of mexico
      mark( 789,  76)=+1	!  gulf of mexico
      mark( 800,  76)=+1	!  gulf of mexico
      mark( 793,  82)=+1	!  gulf of mexico
      mark( 782, 104)=+1	!  gulf of mexico
      mark( 773, 114)=+1	!  gulf of mexico
      mark( 680, 127)=+1	!  gulf of mexico
      mark( 681, 129)=+1	!  gulf of mexico
      mark( 733, 174)=+1	!  gulf of mexico
      mark( 748, 184)=+1	!  gulf of mexico
      mark( 946, 470)=+1	!  brazil
      mark( 954, 477)=+1	!  brazil
      mark( 976, 511)=+1	!  brazil
      mark( 982, 531)=+1	!  brazil
      mark( 989, 563)=+1	!  brazil
      mark( 989, 563)=+1	!  brazil
      mark(1020, 595)=+1	!  brazil
      mark(1040, 622)=+1	!  brazil
      mark(1060, 647)=+1	!  brazil
      mark(1065, 648)=+1	!  brazil
      mark(1067, 653)=+1	!  brazil
      mark(1069, 660)=+1	!  brazil
      mark(1316, 741)=+1	!  brazil
      mark(1313, 747)=+1	!  brazil
      mark(1300, 754)=+1	!  brazil
      mark(1393, 666)=+1	!  brazil
      mark(1365, 709)=+1	!  brazil
      mark( 694,1083)=+1	!  morocco
      mark( 689,1093)=+1	!  morocco
c
      iter=0
 7    num1=0
      num2=0
c
c --- reverse scanning direction after every iteration
      incr=1-2*mod(iter,2)
      if (incr.lt.0) then
        j1=jj1
        i1=ii1
      else
        j1=1
        i1=1
      end if
      do 10 j=j1,jdm-j1,incr
      ja=max(  1,j-1)
      jb=min(jj1,j+1)
      do 10 i=i1,idm-i1,incr
      ia=max(  1,i-1)
      ib=min(ii1,i+1)
c
      if (depth(i,j).gt.0..and.min(mark(i,j),
     .  mark(ia,j),mark(ib,j),mark(i,ja),mark(i,jb)).lt.0) then
ccc      write (*,'(6i5,5x,4i2)')
ccc     .   ia,i,ib,ja,j,jb,mark(ia,j),mark(ib,j),mark(i,ja),mark(i,jb)
      mark(i,j)=-1
      depth(i,j)=0.
      num1=num1+1
      end if
c
      if (depth(i,j).eq.0..and.max(mark(i,j),
     .  mark(ia,j ),mark(ib,j ),mark(i ,ja),mark(i ,jb),
     .  mark(ia,ja),mark(ib,jb),mark(ib,ja),mark(ia,jb)).gt.0) then
ccc      write (*,'(6i5,5x,4i2)')
ccc     .   ia,i,ib,ja,j,jb,mark(ia,j),mark(ib,j),mark(i,ja),mark(i,jb)
      mark(i,j)=+1
      depth(i,j)=2.*cutoff
      num2=num2+1
      end if
 10   continue
      iter=iter+1
      write (*,100) num1,num2,iter
 100  format (i7,' points removed,',i7,' points added in scan',i4)
      if (num1+num2.gt.0) go to 7
c
c      call zebra(depth,idm,ii1,jj1)
      call bigrid(depth)
ccc      do 5 i=1,ii1
ccc 5    write (*,'('' i=''i5'' jfp,jlp=''5(1x,2i5)/(17x,5(1x,2i5)))') i,
ccc     . (jfp(i,l),jlp(i,l),l=1,jsp(i))
ccc      do 6 j=1,jj1
ccc 6    write (*,'('' j=''i5'' ifp,ilp=''5(1x,2i5)/(17x,5(1x,2i5)))') j,
ccc     . (ifp(j,l),ilp(j,l),l=1,isp(j))
c
      call pakk(depth,idm,ii1,jj1,util,length)
      open (unit=10,file='your_corrected_depth_file',status='unknown')
      write (10,'(5(a79/),3i8/(40a2))') preambl,ii1,jj1,length,
     .  (util(l),l=1,length)
c
c --- highlight coastal data points to ease recognition of islands
c
      do 12 j=1,jj1
      do 12 i=1,ii1
      mark(i,j)=0
      if (depth(i,j).gt.0.) mark(i,j)=1
 12   continue
c
      do 13 j=1,jj1
      ja=max(  1,j-1)
      jb=min(jj1,j+1)
      do 13 i=1,ii1
      ia=max(  1,i-1)
      ib=min(ii1,i+1)
      work(i,j)=depth(i,j)
      if (mark(i,j).eq.0.and.mark(ia,j)+mark(ib,j)+mark(i,ja)+
     .  mark(i,jb).gt.0) work(i,j)=-9999.
 13   continue
c
c --- data are written in strips 19 points wide
      jsec=jj1/19
      write (*,'(5(a79/))') preambl
      do 9 jfrst=0,19*(jsec-1),19
      write (*,'('' columns''i5'' --''i5)') jfrst+1,jfrst+19
 9    write (*,fmt2)
     .   (i,(int(work(i,j)+.5),j=jfrst+1,jfrst+19),i=1,ii1)
      jlast=19*jsec
      if (jj1.gt.jlast) then
      write (char2,'(i2)') jj1-jlast
      fmt2(5:6)=char2
      write (*,'('' columns''i5'' --''i5)') jlast+1,jj1
      write (*,fmt2) (i,(int(work(i,j)+.5),j=jlast+1,jj1),i=1,ii1)
      end if
ccc      call prtmsk(ip,depth,work,idm,ii1,jj1,0.,1.,
ccc     .    'water depth - 1437-sq.pts')
c
c --- generate contour plot of depth array
c
      imax=(ii1+1)/2
      jmax=(jj1+1)/2
      do 4 j=1,jmax
      do 4 i=1,imax
 4    work(i,j)=depth(2*i-1,2*j-1)
      call opngks
      call set(.05,.95,.05,.95,1.,float(imax),1.,float(jmax),1)
      call conrec(work,idm,imax,jmax,25.,25.,25.,1,-1,0)
      call conrec(work,idm,imax,jmax,50.,50.,50.,1,-1,0)
      call conrec(work,idm,imax,jmax,100.,100.,100.,1,-1,0)
      call conrec(work,idm,imax,jmax,200.,200.,200.,1,-1,0)
      call conrec(work,idm,imax,jmax,400.,400.,400.,1,-1,0)
      call conrec(work,idm,imax,jmax,800.,800.,800.,1,-1,0)
      call conrec(work,idm,imax,jmax,1600.,1600.,1600.,1,-1,0)
      call conrec(work,idm,imax,jmax,3200.,3200.,3200.,1,-1,0)
      call set(.05,.95,.05,.95,1.,float(ii1),1.,float(jj1),1)
      call dashdb(1+16*(1+16*(1+16*1)))
      do 15 j=20,jdm,20
      call frstd(1.        ,float(j))
      call vectd(float(idm),float(j))
 15   continue
      do 16 i=20,idm,20
      call frstd(float(i),1.        )
      call vectd(float(i),float(idm))
 16   continue
      call dashdb(4095)
      do 17 j=100,jdm,100
      call frstd(1.        ,float(j))
      call vectd(float(idm),float(j))
 17   continue
      do 18 i=100,idm,100
      call frstd(float(i),1.        )
      call vectd(float(i),float(idm))
 18   continue
ccc      call perim(1,jdm,1,idm)
      call perim(1,1,1,1)
      call frame
      call clsgks
      stop
      end