program etopo5_corr c c --- Extract ocean bottom depth from the NAVY Getech 1/24-deg. data set c --- Compute grid-square averages on grid defined by user. c c --- User must define size of output array, depth(idm,jdm), and furnish c --- statement functions 'xoflat' and 'yoflon' returning (fractional) x,y c --- grid indices as functions of latitude,longitude (deg.). The inverse c --- functions 'latofx,lonofy' translating x,y into lat,lon must also be c --- supplied. Grid indices i,j and their real counterparts x,y are assumed c --- to increase southward/eastward respectively. c c --- For this routine to work satisfactorily, spatial resolution of output c --- grid must not exceed that of input grid. c c --- written by Linda Smith and Rainer Bleck, U.Miami, July 1993. c c c parameter (idm=147,jdm=324) c parameter (idm=443,jdm=971) parameter (idm=1071,jdm=1678) real depth(idm,jdm),hplot(jdm,idm),latofx,lonofy ccc . ,weight(idm,jdm) real depth2(idm,jdm),diff(idm,jdm) REAL TLABEL, HTOT, HD(3338,1562) dimension idep(20) character*2 compac(idm*jdm+14) logical skip data itest,jtest/166,1/ c data dx,dy/0.087890625,0.062500000/ c c ----------------------------user-supplied statements---------------------- data radian/57.29578/,pi/3.14159265/ ccc data grid/.9/,xpiv/96./,rlong/-97.5/ ccc data grid/.45/,xpiv/191./,rlong/-98./ c c --- output grid : rectangular, evenly spaced on Mercator projection c --- 'grid' = grid resolution in degrees longitude c --- 'xpiv' = location of equator in 'i' space c --- 'rlong' = longitude of westernmost grid point (j=1) c data grid/.08/,xpiv/1071./,rlong/-98./ c data grid/.1/,xpiv/533./,rlong/-97./ c data grid/.3/,xpiv/177./,rlong/-97./ c xoflat(alat)=xpiv-alog(tan((2.*alat/radian+pi)/4.))*radian/grid yoflon(alon)=1.+(amod(alon-rlong+cut,360.)-cut)/grid c latofx(x)=(2.*atan(exp((xpiv-x)*grid/radian))-pi/2.)*radian lonofy(y)=amod(rlong+(y-1.)*grid+720.,360.) c c --- 'cut' = diff. between rlong and longitude opposite to center of grid cut=amod(rlong-lonofy(.5*float(jdm+1))+540.,360.) c -------------------------------------------------------------------------- c dx=1./24. dy=dx do 104 j=1,jdm do 104 i=1,idm ccc weight(i,j)=0. diff(i,j)=0. 104 depth(i,j)=0. c c --- read through Navy corrected data file c READ( 20,800) HD 800 format(10f10.0) c c do 103 ilat=1,1562 alat=min(65.,max(0.,0.+float(ilat-1)*dy)) c --- get locations of N and S edge of (1/12) deg. wide strip centered c --- on present row of input points xnorth=xoflat(alat+dy/2.) xsouth=xoflat(alat-dy/2.) ccc write (*,'('' latitude''f7.2'' translates into i =''f8.3)') ccc . alat,xoflat(alat) skip=.false. c --- are we outside geographic area of interest? c if (xsouth.lt..5) skip=.true. c --- have we read all the data needed? c if (xnorth.gt.float(idm)+.5) go to 100 c c --- latitude circle cuts through output grid. start extracting data. i=xnorth+.5 if (.not.skip.and.xsouth.gt.float(i)+1.5) then write (*,'('' error -- xnorth,xsouth,i =''2f9.3,i7)') . xnorth,xsouth,i end if x=max(xnorth,min(xsouth,float(i)+.5)) c do 102 jlon=1,3338 c --- don't process unneeded records c if (skip) go to 102 c --- do the 20 values just read overlap the desired longitude interval? alone=262. alonw=262.+(jlon-1)*dx ccc if (i.eq.itest) then ccc write (*,'('' longitude''f7.2'' translates into j =''f8.3)') ccc . alone,yoflon(alone) ccc write (*,'('' longitude''f7.2'' translates into j =''f8.3)') ccc . alonw,yoflon(alonw) ccc end if c if (yoflon(alone+dx/2.).lt..5) go to 102 c if (yoflon(alonw-dx/2.).gt.float(jdm)+.5) go to 102 c c --- Depths below sea level in ETOPO5 data set are negative. value=max(-hd(jlon,ilat),0.) c print *,-hd(jlon,ilat),value if (value.eq.0.) go to 101 c alon=(262.+(jlon-1)*dx) c --- determine longitude of W and E edge of (1/12) x (1/12) deg. square c --- centered on input point ywest=yoflon(alon-dx/2.) yeast=yoflon(alon+dx/2.) c --- are we in the proper longitude interval? ccc if (i.eq.itest) write (*,'('' nrec,jrec,ywest,yeast ='' ccc . 2i5,2f9.2)') nrec,jrec,ywest,yeast if (yeast.lt..5.or.ywest.gt.float(jdm)+.5) go to 101 c c --- process this data point j=ywest+.5 if (yeast.gt.float(j)+1.5) then write (*,'('' error -- ywest,yeast,j =''2f9.3,i7)') . ywest,yeast,j end if y=max(ywest,min(yeast,float(j)+.5)) ccc if (i.eq.itest) write (*,'('' i,j,x,y,value =''2i5,3f9.2)') ccc . i,j,x,y,value c c --- Input depth value is assumed constant in (1/12) x (1/12) deg square c --- centered on given lat/lon point. This square is embedded in area c --- formed by output grid boxes (i,j),(i+1,j),(i,j+1),(i+1,j+1). c --- Divvy up depth value among those 4 boxes. c c print *,'x=',x,' y=',y, value if (i.lt.idm.and.j.lt.jdm) .depth(i+1,j+1)=depth(i+1,j+1)+value*(xsouth-x)*(yeast-y) if (i.gt.0 .and.j.lt.jdm) .depth(i ,j+1)=depth(i ,j+1)+value*(x-xnorth)*(yeast-y) if (i.lt.idm.and.j.gt.0 ) .depth(i+1,j )=depth(i+1,j )+value*(xsouth-x)*(y-ywest) if (i.gt.0 .and.j.gt.0 ) .depth(i ,j )=depth(i ,j )+value*(x-xnorth)*(y-ywest) ccc if (i.lt.idm.and.j.lt.jdm) ccc .weight(i+1,j+1)=weight(i+1,j+1)+(xsouth-x)*(yeast-y) ccc if (i.gt.0 .and.j.lt.jdm) ccc .weight(i ,j+1)=weight(i ,j+1)+(x-xnorth)*(yeast-y) ccc if (i.lt.idm.and.j.gt.0 ) ccc .weight(i+1,j )=weight(i+1,j )+(xsouth-x)*(y-ywest) ccc if (i.gt.0 .and.j.gt.0 ) ccc .weight(i ,j )=weight(i ,j )+(x-xnorth)*(y-ywest) 101 continue 102 continue 103 continue print *,'ok' c c --- User must supply statements to store output array 'depth' c 100 continue c call pakk(depth,idm,idm,jdm,compac,length) write (30,'(3i8/(40a2))') idm,jdm,length, . (compac(l),l=1,length) c call zebra(depth,idm,idm,jdm) c c c --- plot output depth array - Mercator projection call opngks call gsclip(0) call maproj('MERCATOR',0.,0.,0.) call mapstc('OUTLINE','CO') call mapstl('DOT',.FALSE.) call mapset('CORNERS',latofx(1.),lonofy(1.), . latofx(float(idm)),lonofy(float(jdm))) call mapsti('DA',34952) call mapsti('LA',0) call mapint call mapdrw call getset(xa,xb,ya,yb,xc,xd,yc,yd,linlog) call set(xa,xb,ya,yb,1.,float(jdm),1.,float(idm),1) do 4305 i=1,jdm do 4305 j=1,idm 4305 hplot(i,j)=depth(idm+1-j,i) call conrec(hplot,jdm,jdm,idm,0.,0.,500.,1,-1,0) call perim(1,jdm,1,idm) call frame call clsgks stop '(e-o-f)' end c c SUBROUTINE PAKK(ARRAY,IDIM,II,JJ,COMPAC,LENGTH) C C --- Converts the contents of -ARRAY- into an ASCII character string which C --- is stored in character*2 array -COMPAC-. COMPAC(1)...COMPAC(7) contain C --- the base value, i.e., the minimum value encountered in -ARRAY-. C --- COMPAC(8)...COMPAC(14) contain a scale factor by which the individual C --- 6-bit integers encoded as ASCII character pairs in COMPAC(8),... C --- COMPAC(LENGTH) must be multiplied before the base value is added C --- during an UNPAKKing operation. Base value and scale factor are encoded C --- in E14.7 format. C C --- The printable ASCII characters used to encode the integers include C --- the numbers 0...9, upper- and lower-case letters A...Z, a...z, plus C --- two additional characters '.' and '/' (total of 64). C C --- A packing operation fills (II*JJ+14) array elements in -COMPAC- which C --- must be dimensioned accordingly in the calling program. The total C --- number of occupied array elements is returned in -LENGTH-. In calls to C --- UNPACK, -LENGTH- is treated as input variable. C real array(idim,1) character*2 char,compac(1),comp2(14) character*14 comp14(2) equivalence (comp2,comp14) data nbits/12/ base=1.e22 do 1 i=1,ii do 1 j=1,jj 1 base=amin1(base,array(i,j)) scal=0. do 2 i=1,ii do 2 j=1,jj 2 scal=amax1(scal,array(i,j)-base) scal=scal/float(2**nbits-1) i1=0 i2=0 length=14 do 3 i=1,ii do 3 j=1,jj if (scal.eq.0.) go to 7 numb=(array(i,j)-base)/scal+.5 i1=numb/64 i2=numb-64*i1 c c --- Map 6-bit numbers onto character set consisting of numbers c --- 0...9, letters A...Z, a...z, and the two characters '.' and '/'. c --- (If mapping into the character range 32...95 -- which includes the c --- characters !"#$%&'()*+,-./:;<=>?@[\]^_ -- is deemed safe, delete c --- the next 6 lines.) c if (i1.gt.37) i1=i1+6 c if (i1.gt.11) i1=i1+7 c i1=i1+14 c if (i2.gt.37) i2=i2+6 c if (i2.gt.11) i2=i2+7 c i2=i2+14 c 7 length=length+1 compac(length)(1:1)=char(i1+32) compac(length)(2:2)=char(i2+32) 100 format (a2) 3 continue write (comp14(1),101) base write (comp14(2),101) scal 101 format (1pe14.7) do 8 l=1,14 8 compac(l)=comp2(l) c return c c ENTRY UNPAKK(ARRAY,IDIM,II,JJ,COMPAC,LENGTH) c do 9 l=1,14 9 comp2(l)=compac(l) read (comp14(1),101) base read (comp14(2),101) scal lngth=14 do 4 i=1,ii do 4 j=1,jj lngth=lngth+1 i1=ichar(compac(lngth)(1:1)) i2=ichar(compac(lngth)(2:2)) c c --- 6-bit numbers are mapped onto character set consisting of numbers c --- 0...9, letters A...Z, a...z, and the two characters '.' and '/'. c --- (If mapped into character range 32...95, delete next 6 lines) c if (i1.gt.96) i1=i1-6 c if (i1.gt.64) i1=i1-7 c i1=i1-14 c if (i2.gt.96) i2=i2-6 c if (i2.gt.64) i2=i2-7 c i2=i2-14 c 4 array(i,j)=scal*float(64*(i1-32)+(i2-32))+base if (lngth.ne.length) stop 'UNPACK' return end c c SUBROUTINE ZEBRA(ARRAY,IDIM,II,JJ) C C --- FIND NICE CONTOUR INTERVAL RESULTING IN 7 TO 10 CONTOUR LINES AND C --- DRAW CONTOURS ON LINE PRINTER THROUGH THE FOLLOWING SET OF GRID POINTS: C C ARRAY( 1, 1) . . . . . . . . . ARRAY( 1,JJ) C . . C . . (PLOT WILL APPEAR C . . ON PAPER AS SHOWN, C . . I DOWN, J ACROSS) C . . C ARRAY(II,JJ) . . . . . . . . . ARRAY(II,JJ) C C --- II MAY BE SMALLER THAN IDIM, THE FIRST (ROW) DIMENSION OF 'ARRAY' C --- IN THE CALLING PROGRAM. THUS, PLOTTING OF PARTIAL ARRAYS IS POSSIBLE. C REAL ARRAY(IDIM,1) DATA SQRT2/1.414/ C AMX=-1.E25 AMN= 1.E25 DO 1 I=1,II DO 1 J=1,JJ AMX=AMAX1(AMX,ARRAY(I,J)) 1 AMN=AMIN1(AMN,ARRAY(I,J)) C IF (AMX.GT.AMN) GO TO 2 WRITE(6,100)ARRAY(1,1) 100 FORMAT (//' FIELD TO BE CONTOURED IS CONSTANT ...'E15.5/) RETURN C 2 CONTUR=(AMX-AMN)/6. Q=10.**INT(ALOG10(CONTUR)) IF (CONTUR.LT.1.) Q=Q/10. RATIO=CONTUR/Q IF (RATIO.GT.SQRT2*5.) CONTUR=Q*10. IF (RATIO.LE.SQRT2*5.) CONTUR=Q*5. IF (RATIO.LE.SQRT2*2.) CONTUR=Q*2. IF (RATIO.LE.SQRT2) CONTUR=Q WRITE(6,101)CONTUR,AMN,AMX 101 FORMAT (' CONTOUR INTERVAL IN PLOT BELOW IS'1PE9.1, . 6X'MIN/MAX ='2E11.3/) CALL DIGPLT(ARRAY,IDIM,II,JJ,CONTUR) C RETURN END C C SUBROUTINE DIGPLT(ARRAY,IDIM,II,JJ,DEC) C C --- SIMULATE A CONTOUR LINE PLOT ON THE PRINTER C REAL ARRAY(IDIM,1) CHARACTER*1 DIGIT(130),DIG(20) DATA DIG/'0',' ','1',' ','2',' ','3',' ','4',' ', . '5',' ','6',' ','7',' ','8',' ','9',' '/ C C NCHAR = NUMBER OF CHARACTER INCREMENTS IN 'J' DIRECTION C RATIO = CHARACTER WIDTH / LINE SPACING C DATA NCHAR/74/,RATIO/.58/ XINC=FLOAT(JJ-1)/(FLOAT(NCHAR)*RATIO) YINC=FLOAT(JJ-1)/ FLOAT(NCHAR) K=FLOAT(NCHAR)*RATIO*FLOAT(II-1)/FLOAT(JJ-1)+1.00001 DO 1 I=1,K X=1.+FLOAT(I-1)*XINC IA=MIN0(II-1,INT(X)) DX=X-FLOAT(IA) DO 2 J=1,NCHAR+1 Y=1.+FLOAT(J-1)*YINC JA=MIN0(JJ-1,INT(Y)) DY=Y-FLOAT(JA) DXDY=DX*DY VALUE=ARRAY(IA,JA)*(1.-DX-DY+DXDY) . +ARRAY(IA+1,JA)*(DX-DXDY) . +ARRAY(IA,JA+1)*(DY-DXDY) . +ARRAY(IA+1,JA+1)*DXDY N=MOD(MOD(INT(2.*VALUE/DEC+SIGN(.5,VALUE)),20)+20,20)+1 2 DIGIT(J)=DIG(N) 1 WRITE(6,100)'I',' ',(DIGIT(J),J=1,NCHAR+1),' ','I' 100 FORMAT(1X,130A1) RETURN END