program hypso c Author: Dr. Peter W. Sloss c NOAA - NGDC character region*40,more*1,outfile*30 dimension clats(61),sums(100) dimension depth(61,61),bins(0:100) DOUBLE PRECISION TAREA,SUBTOT,SUBTOP,SUMS,TAMST DATA REGION /' '/ write(*,'(a)') 1 ' National Geophysical Data Center ETOPO5 Hypsometry' write(*,'(a)') 1 ' computed for 5-degree square blocks' call getet5(0,0.,0.,0.) write(*,'(a)') ' Enter OUTPUT FILE name (=>screen): ' read(*,'(a)') outfile iout = 6 if (outfile .ne. ' ') then iout = 2 open(iout,file=outfile) end if write(*,'(a,\)')' Enter a name for the specified REGION:' READ(*,'(A)') REGION write(*,'(a,\)') 1 ' How many Depth/Elevation Bins, or zero (0) for auto: ' READ *,NBIN nbins = iabs(nbin) write(*,'(a,\)') ' Enter minimum elevation (bin # 0) ' READ *,BINS(0) if (NBIN .le. 0) then write(*,'(a,\)') ' Enter bin size (m): ' read*,bsiz write(*,'(a,\)') ' Enter maximum elevation: ' read *,bmax nbins = (bmax - bins(0))/bsiz if (nbins .gt. 0) then do 35 nb = 1,nbins 35 BINS(nb) = BINS(nb-1) + bsiz print*,'The ',nbins,' bins are: ',(bins(nb),nb=1,nbins) end if else write(*,'(a)') 1 ' Enter the Bin Boundaries (max value/bin, INCREASING+):' C Initialize some constants READ *,(BINS(I),I=1,NBIN) end if DTOR = 0.0174532925 ACELL = 86.0576299 TAREA = 0.0 C Cell areas computed in square Kilometers. C Now, get down to it: 1 write(*,'(a)') ' Give LAT,LON of 5-degree block:' READ(*,*,ERR=1) SWLAT1,SWLON1 write(*,'(a,\)') 1 ' Enter # 5-deg blocks Eastward and Northward from start:' READ*,NEAST,NORTH if (NEAST.lt.1) NEAST = 1 if (NORTH.lt.1) NORTH = 1 ELAT = NORTH*5 + SWLAT1 - .05 ELON = NEAST*5 + SWLON1 - .05 DO 31 SWLON = SWLON1,ELON,5.0 DO 31 SWLAT = SWLAT1,ELAT,5.0 C Get the data CALL getet5(IFOUND,DEPTH,SWLAT,SWLON) C IF (IFOUND .NE. 1) GO TO 1 C Initialize the bins DO 10 J=1,61 AT = SWLAT + FLOAT(J-1)/12.0 10 CLATS(J) = COS(AT*DTOR) * ACELL ATS = SWLAT ATN = SWLAT + 5. ONW = SWLON ONE = SWLON + 5. C CLATS initialized to area for each latitude in sq km. I1 = 1 I2 = 60 J1 = 1 J2 = 60 IMASK = 0 if (NEAST.lt.0 .or. NORTH.lt.0) then PRINT*,'Enter 1 to mask out part of this cell, else 0:' READ (*,ERR=11) IMASK end if IF(IMASK .NE. 1) GO TO 11 4 write(*,'(a,\)') ' Latitude range within this cell (S,N)?' READ(*,*,ERR=4) ATS,ATN write(*,'(a,\)') ' Longitude range within this cell? (W,E)?' READ(*,*,ERR=4) ONW,ONE IF (ONW .LT. 0.0) ONW = ONW + 360. IF (ONE .LT. 0.0) ONE = ONE + 360. C ** Compute I,J limits for partial cells IF (AMOD(ONE,5.) .EQ. 0) ONE = ONE - 0.005 IF (AMOD(ATN,5.) .EQ. 0) ATN = ATN - 0.005 IF(ATN .LE. 0.0) ATN = ATN + 90.0 IF(ATS .LE. 0.0) ATS = ATS + 90.0 C ** Don't sweat the added latitude -- just for determining i,j limits C ** Above keeps AMOD from throwing away whole cells! I1 = AMOD(ONW,5.)*12 + 1 I2 = AMOD(ONE,5.)*12 + 1 J1 = AMOD(ATS,5.)*12 + 1 J2 = AMOD(ATN,5.)*12 + 1 11 continue if (IMASK.eq.1) then PRINT*,'Latitude (J) and Longitude (I) Limits are:' PRINT*,'J1, J2 = ',J1,J2,', and I1, I2 = ',I1,I2 end if C Now, add up the depths DO 20 J=J1,J2 DO 20 I=I1,I2 C ** Leave off top and east boundaries for no overlaps. DO 15 N=1,NBINS IF (DEPTH(I,J) .GT. BINS(N)) GO TO 15 C ** Add to total in appropriate bin SUMS(N) = SUMS(N) + CLATS(J) GO TO 16 15 CONTINUE C ** End Bin filling 16 CONTINUE TAREA = TAREA + CLATS(J) 20 CONTINUE 31 CONTINUE C ** Now, check if there are more cells to process-- write(*,'(a,\)') ' Are there more Areas? (Y/N)' READ(*,'(A)') MORE IF(MORE .EQ. 'Y' .OR. MORE .EQ. 'y') GO TO 1 C ** Do the Statistics! WRITE(iout,'(///)') write(iout,'(a,f5.0,a,f5.0,a,/,12x,f6.0,a,f6.0,a,/)') 1 ' Area covers ',SWLAT1,' to ',SWLAT1+5.*NORTH, 2 ' deg. latitude, and ',swlon1,' to ',swlon1+neast*5, 2 ' deg. longitude.' write(iout,'(a)') 1 ' The Hypsometry of your data areas is as follows:' write(iout,'(a,a)') 1 ' -----------------------------------------------------', 2 ' ---------------' write(iout,'(a,a)') 1 ' Bin # from to (m) Area (sq.km) % of Total', 2 ' Cum. %' SUBTOT = 0.0 SUBTOP = 0.0 DO 30 N=1,NBINS SUBTOT = SUBTOT + SUMS(N) SUBTOP = SUBTOT/TAREA WRITE(iout,'(I5,2F10.0,F14.3,2F15.5)') N,BINS(N-1),BINS(N), 1 SUMS(N),100.0*SUMS(N)/TAREA,100.0*SUBTOP 30 CONTINUE TAMST = TAREA - SUBTOT WRITE(iout,'(a,F10.0,10X,F14.3,F15.5)') 1 ' >',BINS(N-1),TAMST,100.0*TAMST/TAREA c PRINT*,'-----------------------------------------------------', c 1 '---------------' WRITE(iout,'('' Total area of region "'',A,/,''" = '',F14.1, 1 '' sq. km.'')') REGION, TAREA END subroutine getet5(ifound,outbuf,swlat,swlon) c ** get a 5-degree square of ETOPO5 data save init integer*2 inbuf(4320) dimension outbuf(61,61) data init /0/ if (init .eq. 0) then open(1,file='ETOPO5',status='old', 1 access='direct',form='unformatted',recl=8640,err=999) init = 1 return end if nr = 0 c ** repeated calls to this subroutine read more 5-degree blocks alat = swlat alon = swlon num1 = (90. - alat)*12 + 1 nr = nr + 1 if (num1 .le. 60) num1 = 61 num2 = num1 - 60 c ** Find first and last data points... lon1 = amod(alon+360.,360.)*12 + 1 lon2 = lon1 + 60 j = 1 do 10 num = num1,num2,-1 if (num .le. 2160) then read(1,rec=num) inbuf do 2 i = 0,60 jj = mod(i+lon1-1,4320) + 1 outbuf(i+1,j) = inbuf(jj) 2 continue else c ** south pole: elevation = 2810 m.-- insert where needed!! do 3 i = 1,61 3 outbuf(i,j) = 2810 end if nr = nr + 1 j = j + 1 10 continue print*,'Block done for lat/lon = ',alat,alon return 999 write(*,'(a)') '* * * * * * * * * * * * * * * * * * * * * *' write(*,'(a)') ' I can''t find the ETOPO5 data file! Please ' write(*,'(a)') ' type the following line:' write(*,'(a)') ' SET ETOPO5=[path to ETOPO5 data file]' write(*,'(a)') ' -- replacing [...] with the ETOPO5 file path' write(*,'(a)') ' e.g., SET ETOPO5=L:\TOPO\ETOPO5\ETOPO5.DOS' write(*,'(a)') ' and start program again' stop end