function xdist,rlat,rlon,ip,rp
; geographic to map coordinate transformation for x coordinate

alpha=rp(6)*!dtor

; transformation for approximate polyconic projection

if(ip(4) eq 0) then begin
 dlat=rlat-rp(0)
 dlon=rp(1)-rlon
 xp=rp(7)*dlon+rp(8)*dlat+rp(9)*dlon*dlat+rp(10)*dlon*dlon
 yp=rp(11)*dlon+rp(12)*dlat+rp(13)*dlat*dlon+rp(14)*dlon*dlon
 xp=xp*1000.
 yp=yp*1000.
end

; transformation for lambert conformal conic projection

if(ip(4) eq 1) then begin
 r=(tan(!pi/4.-!dtor*rlat/2.))^rp(10)
 xp=rp(11)*r*sin(rp(10)*(rp(7)-!dtor*rlon))
 yp=-rp(11)*r*cos(rp(10)*(rp(7)-!dtor*rlon))
 xp=xp-rp(12)
 yp=yp-rp(13)
end

; transform to 'unprimed' system

;  first rotate

xd=xp*cos(alpha)+yp*sin(alpha)

; now translate

xd=xd+ip(2)*rp(2)
 
return,xd
end

;----------------------------------------------------------------------
function ydist,rlat,rlon,ip,rp
; geographic to map coordinate transformation for y coordinate

alpha=rp(6)*!dtor

; transformation for approximate polyconic projection

if(ip(4) eq 0) then begin
 dlat=rlat-rp(0)
 dlon=rp(1)-rlon
 xp=rp(7)*dlon+rp(8)*dlat+rp(9)*dlon*dlat+rp(10)*dlon*dlon
 yp=rp(11)*dlon+rp(12)*dlat+rp(13)*dlat*dlon+rp(14)*dlon*dlon
 xp=xp*1000.
 yp=yp*1000.
end

; transformation for lambert conformal conic projection

if(ip(4) eq 1) then begin
 r=(tan(!pi/4.-!dtor*rlat/2.))^rp(10)
 xp=rp(11)*r*sin(rp(10)*(rp(7)-!dtor*rlon))
 yp=-rp(11)*r*cos(rp(10)*(rp(7)-!dtor*rlon))
 xp=xp-rp(12)
 yp=yp-rp(13)
end

; transform to 'unprimed' system

;  first rotate

yd=yp*cos(alpha)-xp*sin(alpha)

; now translate

yd=yd+ip(3)*rp(2)

return,yd
end

;----------------------------------------------------------------------
function rlat,x,y,ip,rp  
; find lat for x,y (meters) on bathymetric grid  

alpha=rp(6)*!pi/180.  
xx=x-ip(2)  
yy=y-ip(3)  
xp=xx*cos(alpha)-yy*sin(alpha)
yp=yy*cos(alpha)+xx*sin(alpha)

; transformation for approximate polyconic projection

if(ip(4) eq 0) then begin
 xp=xp/1000.
 yp=yp/1000.
 dl=rp(19)*xp+rp(20)*yp+rp(21)*xp*yp+rp(22)*xp*xp  
 rl=rp(0)+dl  
end

; transformation for lambert conformal conic projection

if(ip(4) eq 1) then begin
 xp=xp+rp(12)
 yp=yp+rp(13)
 xx=xp/rp(11)
 yy=yp/rp(11)
 r=sqrt(xx^2+yy^2)
 rl=360.*(!pi/4.-atan(r^(1./rp(10))))/!pi
end

return,rl  
end  
  
;----------------------------------------------------------------------
function rlon,x,y,ip,rp  
; find lon for x,y (meters) on bathymetric grid  
alpha=rp(6)*!pi/180.  
xx=x-ip(2)
yy=y-ip(3)  
xp=xx*cos(alpha)-yy*sin(alpha)  
yp=yy*cos(alpha)+xx*sin(alpha)  

; transformation for approximate polyconic projection

if(ip(4) eq 0) then begin
 xp=xp/1000.
 yp=yp/1000.
 dl=rp(15)*xp+rp(16)*yp+rp(17)*xp*yp+rp(18)*xp*xp  
 rl=rp(1)-dl  
end

; transformation for lambert conformal conic projection

if(ip(4) eq 1) then begin
 xp=xp+rp(12)
 yp=yp+rp(13)
 xx=xp/rp(11)
 yy=yp/rp(11)
 rl=180.*(rp(7)-atan(xx,-yy)/rp(10))/!pi
end
return,rl  
end  
 

;----------------------------------------------------------------------
pro rgrid,fname,gname,dgrid,iparm,rparm

; Purpose:
;           To read a bathymetric grid data file
;           and return grid parameters and depths.
; Arguments:
;   On input:
;           fname - Full pathway name of bathymetric data file
;   On output:
;           gname - Lake name (from bathymetric data file)
;           dgrid - Depth grid
;           iparm - Integer parameters for grid description and
;                      coordinate conversion
;           rparm - Real parameters for grid description and
;                      coordinate conversion
;  
;       ------------------------------------------------------------------
;                      FORMAT OF BATHYMETRIC DATA FILE
;       ------------------------------------------------------------------
;                 FIELD                            FORTRAN FORMAT  COLUMNS
;       ------------------------------------------------------------------
;       RECORD 1: LAKE NAME                            40A1          1-40
;       RECORD 2: FIRST (I) DIMENSION OF DEPTH ARRAY   I5            1-5
;                 SECOND (J) DIMENSION OF DEPTH ARRAY  I5            6-10
;                 BASE LATITUDE                        F12.7        11-22
;                 BASE LONGITUDE                       F12.7        23-34
;                 GRID SIZE                            F5.0         35-39
;                 MAXIMUM DEPTH                        F5.0         40-44
;                 MINIMUM DEPTH                        F5.0         45-49
;                 BASE ROTATION (CCW FROM E-W)         F6.2         50-55
;                 I DISPLACEMENT                       I5           56-60
;                 J DISPLACEMENT                       I5           61-65
;                 ROTATION FROM BASE (CCW)             F6.2         66-71
;                 POWER OF TEN TO CONVERT DEPTHS TO 
;                                       METERS         I3           72-74
;                 MAP PROJECTION INDICATOR
;                  0=APPROXIMATE POLYCONIC
;                  1=LAMBERT CONFORMAL CONIC           I2           75-76
;       RECORDS 3-6 (FOR APPROXIMATE POLYCONIC PROJECTION):
;                 MAP PROJECTION COORDINATE CONVERSION
;                  COEFFICIENTS                       4E15.6         1-60
;       RECORD 3 (FOR LAMBERT CONFORMAL CONIC PROJECTION):
;                 CENTRAL MERIDIAN OF PROJECTION       E15.6         1-15
;                 SOUTHERNMOST STANDARD PARALLEL       E15.6        16-30
;                 NORTHERNMOST STANDARD PARALLEL       E15.6        31-45
;       FOLLOWING RECORDS:
;                 DEPTHS IN ASCENDING I, ASCENDING J
;                  SEQUENCE, 19 TO A RECORD           19F4.0         1-76
;       ------------------------------------------------------------------
;
; ON OUTPUT:
;                     D - DEPTH ARRAY. ZERO FOR LAND, AVERAGE DEPTH
;                         OF GRID BOX IN METERS FOR WATER.
;                     RPARM - ARRAY CONTAINING REAL-VALUED BATHYMETRIC
;                            GRID PARAMETERS AS FOLLOWS:
;                        1.  BASE LATITUDE
;                        2.  BASE LONGITUDE
;                        3.  GRID SIZE (M)
;                        4.  MAXIMUM DEPTH (M)
;                        5.  MINIMUM DEPTH (M)
;                        6.  BASE ROTATION (COUNTERCLOCKWISE FROM E-W)
;                        7.  ROTATION FROM BASE (COUNTERCLOCKWISE)
;
;                      FOR IPARM(45)=0 (APPROXIMATE POLYCONIC PROJECTION):
;                        8-11. GEOGRAPHIC-TO-MAP COORDINATE CONVERSION
;                             COEFFICIENTS FOR X
;                        12-15.  GEOGRAPHIC-TO-MAP COORDINATE CONVERSION
;                                COEFFICIENTS FOR Y
;                        16-19.  MAP-TO-GEOGRAPHIC COORDINATE CONVERSION
;                                COEFFICIENTS FOR LONGITUDE
;                        20-23.  MAP-TO-GEOGRAPHIC COORDINATE CONVERSION
;                                COEFFICIENTS FOR LATITUDE
;
;                      FOR IPARM(45)=1 (LAMBERT CONFORMAL CONIC PROJECTION):
;                        8.  CENTRAL MERIDIAN OF PROJECTION (RADIANS)
;                        9.  SOUTHERNMOST STANDARD PARALLEL (RADIANS)
;                        10. NORTHERNMOST STANDARD PARALLEL (RADIANS)
;                        11. LOGARITHMIC COEFFICIENT FOR TRANSFORMATIONS
;                        12. DISTANCE SCALING FACTOR FOR TRANSFORMATIONS
;                        13. X DISPLACEMENT OF BATHYMETRIC GRID ORIGIN
;                             FROM MAP PROJECTION ORIGIN
;                        14. Y DISPLACEMENT OF BATHYMETRIC GRID ORIGIN
;                             FROM MAP PROJECTION ORIGIN
;
;                    IPARM - ARRAY CONTAINING INTEGER-VALUED BATHYMETRIC
;                            GRID PARAMETERS AS FOLLOWS:
;                        1.  NUMBER OF GRID BOXES IN X DIRECTION
;                        2.  NUMBER OF GRID BOXES IN Y DIRECION
;                        3.  I DISPLACEMENT - THE NUMBER OF NEW GRID
;                            SQUARES IN THE X-DIRECTION FROM THE NEW
;                            GRID ORIGIN TO THE OLD GRID ORIGIN
;                            (USED ONLY FOR IPARM(45)=0)
;                        4.  J DISPLACEMENT - THE NUMBER OF NEW GRID
;                            SQUARES IN THE Y-DIRECTION FROM THE NEW
;                            GRID ORIGIN TO THE OLD GRID ORIGIN
;                            (USED ONLY FOR IPARM(45)=0)
;                        5. MAP PROJECTION USED FOR BATHYMETRIC GRID:
;                            0=APPROXIMATE POLYCONIC (GREAT LAKES GRIDS)
;                            1=LAMBERT CONFORMAL CONIC
gname=' '
rparm=fltarr(23)
iparm=intarr(5)
openr,1,fname
readf,1,gname
iw1=intarr(2)
iw2=intarr(2)
rw1=fltarr(6)
on_ioerror,next
readf,1,iw1,rw1,iw2,r1,idexp,iproj,  $
 format='(2I5,2F12.7,3F5.0,F6.2,2I5,F6.2,I3,I2)'
goto,next1
next:iw2=[0,0]
r1=0
idexp=0
iproj=0
next1:on_ioerror,null
iparm(0:1)=iw1
rparm(0:5)=rw1
iparm(2:3)=iw2
iparm(4)=iproj
rparm(6)=r1
if(iproj eq 0) then begin
 rw1=fltarr(16)
 readf,1,format='(4e15.6)',rw1
 rparm(7:22)=rw1
end
if(iproj eq 1) then begin
 rw1=fltarr(3)
 readf,1,format='(3e15.6)',rw1
 rparm(7:9)=rw1
end
im=iparm(0)
jm=iparm(1)
dgrid=fltarr(im,jm)
readf,1,dgrid,format='(19f4.0)'
close,1

; adjust depths

dfac=10.^idexp
rparm(4)=rparm(4)*dfac
rparm(5)=rparm(5)*dfac
dgrid=dgrid*dfac

; for lambert projection compute required constants

if(iparm(4) eq 1) then begin
 a45=atan(1.)
 rparm(7)=rparm(7)*!dtor
 rparm(8)=rparm(8)*!dtor
 rparm(9)=rparm(9)*!dtor
; alon0=rparm(7)
 a1=rparm(8)
 a2=rparm(9)
 rparm(10)=(alog(cos(a1))-alog(cos(a2)))/          $
    (alog(tan(a45-a1/2.))-alog(tan(a45-a2/2.)))

; set scale factor for n-s distance from a1 to a2

 aexp=rparm(10)
 y1=(tan(a45-a1/2.))^aexp
 y2=(tan(a45-a2/2.))^aexp
 rparm(11)=6378140.*(a2-a1)/(y1-y2)
 rparm(12)=0.
 rparm(13)=0.
 dx=xdist(rparm(0),rparm(1),iparm,rparm)
 dy=ydist(rparm(0),rparm(1),iparm,rparm)
 rparm(12)=dx
 rparm(13)=dy
end
return
end
;-------------------------------------------------------------------------
pro shore,fname,ip,rp
openr,1,fname
readf,1,nseg
for i=1,nseg do begin
 readf,1,m,n
 a=fltarr(2,n)
 readf,1,a
 x=xdist(a(1,*),a(0,*),ip,rp)
 y=ydist(a(1,*),a(0,*),ip,rp)
 plots,x/rp(2),y/rp(2)
end
close,1
return
end

;-------------------------------------------------------------------------
; main program
;-------------------------------------------------------------------------
rgrid,'mich2km.dat',gn,dgrid,ip,rp
im=ip(0)
jm=ip(1)
ptype='ps'
set_plot,ptype

;  set the plot scaling

xs=5.
device,xsize=xs,ysize=jm*xs/im,/inch,xoffset=0.5,yoffset=0.5   
!x.s=[0.,1./im]
!y.s=[0.,1./jm] 

; calculate area, volume, and mean depth

dx=rp(2)
k=where(dgrid ne 0,n)
area=dx*dx*n
volume=dx*dx*total(dgrid(k))
dmean=volume/area
print,'Area=      ',area
print,'Volume=    ',volume
print,'Mean depth=',dmean
;  
;    set array htop equal to depth array dgrid
;    set land values of scalar to special value to avoid plotting
;
spval=9999.0  
htop=dgrid
htop(where(dgrid eq 0))=spval
;  
;   contour the depth array htop 
; 
contour,htop,findgen(im)+0.5,findgen(jm)+0.5,  $   
  levels=[20.,40.,60.,80.,100.,120.,140.,160.,180.,200.,  $
	   220.,240.,260.,280.], $ 
;   levels=[50.,100.,150.,200.,250.],   $
  c_linestyle=[1,1,1,1,1,1,1,1,1,1,1,1,1,1],  $
  XTITLE=gname, $
  c_labels=[1,1,1,1,1,1,1,1,1,1,1,1,1,1],   $      
  xstyle=1,ystyle=1, $ 
  xmargin=[0,0],ymargin=[0,0] ,$  
  xrange=[0,im],yrange=[0,jm],  $   
  max_value=spval   

; plot grid outline

for i=0,im-2 do begin
 for j=0,jm-1 do begin
  if((dgrid(i,j) eq 0) and (dgrid(i+1,j) ne 0)) then begin
   plots,[i+1,i+1],[j,j+1],thick=2 
  endif
  if((dgrid(i,j) ne 0) and (dgrid(i+1,j) eq 0)) then begin
   plots,[i+1,i+1],[j,j+1],thick=2 
  endif  
 endfor
endfor

for i=0,im-1 do begin
 for j=0,jm-2 do begin  
  if((dgrid(i,j) eq 0) and (dgrid(i,j+1) ne 0)) then begin
   plots,[i,i+1],[j+1,j+1],thick=2 
  endif
  if((dgrid(i,j) ne 0) and (dgrid(i,j+1) eq 0)) then begin
   plots,[i,i+1],[j+1,j+1],thick=2 
  endif
 endfor
endfor   
  
shore,'mi_shore.dat',ip,rp
device,/close 
end