c--------( PLMAPS : Plot contour map with shading )-------- parameter (MAXSZ=3001*3001) parameter (rad=(180./3.14159)) c dimension g(MAXSZ), s(MAXSZ) character*80 wdr, flog, fnin, fnam character area*8, s4*4, yn*1, spl*3/' '/ character out*72, buf*80, sdtm*20 external cviken c call premsg('PLMAPS : Plot contour map with shading') call opnpin() ldr = lwkdir(50,wdr) + 1 flog = wdr fnin = wdr fnam = wdr call gparma('Enter LOG filename ==> ', 81-ldr, flog(ldr:80)) call gparma('Enter Input data filename ==> ', * 81-ldr, fnin(ldr:80)) open(10,file=fnin,status='old') call gparma('Enter Output PS filename ==> ', * 81-ldr, fnam(ldr:80)) if(fnam(ldr:ldr).eq.' ') fnam(1:1) = char(0) if (fnam(80:80).ne.' ') then call abendm('too long filename') endif open(9,status='scratch') c 1 read(10,'(a)',end=8) buf if(buf(1:1).eq.'#') goto 1 2 read(buf,'(a8,i4,4x,4i8)') area, ncc, iorg, korg, ispa, ispb if(ncc.ge.200) then mc = 1 nc = ncc - 200 else mc = 0 nc = ncc endif if(nc.ge.1.and.nc.le.62) then iorg = 0 ispa = 0 ispb = 0 korg = (nc-30)*360 - 180 if(nc.gt.60) korg = 0 if(nc.eq.61) iorg = +5400 if(nc.eq.62) iorg = -5400 endif read(10,*) ixs, iys, mszx, mszy, nmx, nmy, vnul xs = float(ixs)/1000. ys = float(iys)/1000. call premsg(' ') write(out,'(a,a8,a,i4,4x,4i8)') * '--- Area : ', area, ' Coord.', ncc, iorg,korg, ispa,ispb call premsg(out) write(out,'(14x,2a10,2x,2a8,2a6)') * 'xs', 'ys', 'msz-x', 'msz-y', 'mx', 'my' call premsg(out) write(out,'(16x,2f10.2,2i8,2i6)') xs, ys, mszx, mszy, nmx, nmy call premsg(out) if((nmx*nmy).gt.MAXSZ) then call abendm('too many grids') endif read(10,*) ((g(i+j),i=0,(nmx-1)*nmy,nmy),j=1,nmy) do 4 i=1,nmx*nmy if(g(i).ne.vnul) goto 5 4 continue call abendm('no valid data') 5 n = 1 gsum = g(i) gmin = g(i) gmax = g(i) do 6 j=i+1,nmx*nmy if(g(j).eq.vnul) goto 6 n = n + 1 gsum = gsum + g(j) if(g(j).lt.gmin) gmin = g(j) if(g(j).gt.gmax) gmax = g(j) 6 continue gsum = gsum / float(n) n = 0 gdev = 0. do 7 i=1,nmx*nmy if (g(i).eq.vnul) goto 7 dg = g(i) - gsum n = n + 1 gdev = gdev + dg*dg 7 continue gdev = sqrt(gdev / float(n)) write(out,'(4x,3(a,f7.1),a,f6.1)') 'Data Range : [', gmin, * ' ,', gmax, ' ] Mean :', gsum, ' St.Dev.:', gdev call premsg(out) call premsg('') c if(spl.eq.' ') then call gparmi(' Select 1(portrait) or 2(landscape) ==> ', kpl) if(kpl.eq.1) spl = 'A4P' if(kpl.eq.2) spl = 'A4L' if(spl.eq.' ') call abendm('invalid value') call psopn(fnam, spl) endif call premsg(' Specify contour interval ' * // '(minus to draw even crowded contour)') call gparmi(' ( enter 0 to skip countouring ) ==> ', istep) if(istep.eq.0) goto 3 istp = iabs(istep) lcap = 5 if(istp.eq.25.or.istp.eq.50.or.istp.eq.250) lcap = 4 if(istep.lt.0) lcap = -lcap call premsg('Specify illuminant direction (in Degrees)') call gparmf(' Azimuth angle : ', az) call gparmf(' Elevation angle : ', el) call premsg('Assign vertical magnification') call gparmf(' Value to be scaled to 1km : ', vs) fx = vs * float(mszx)/1000. fy = vs * float(mszy)/1000. cx = cos(el/rad)*cos(az/rad) cy = cos(el/rad)*sin(az/rad) cz = sin(el/rad) if(kpl.eq.1) then call gparmf(' Width of Drawing (in cm) ? ==> ', wide) fscl = float((nmy-1)*mszy) / wide high = float((nmx-1)*mszx) / fscl else call gparmf(' Height of Drawing (in cm) ? ==> ', high) fscl = float((nmx-1)*mszx) / high wide = float((nmy-1)*mszy) / fscl endif fscl = fscl * 100. ccm = fscl/100000. kscl = 0 sclr = 0. sclu = 0. klll = 0 kshl = 0 c if(fscl.lt.3333334..and.fscl.ge.20000.) then km = nint(fscl)/25000 if(km.eq.0) km = 1 if(km.gt.50) km = 50 if(km.gt.2.and.km.lt.5) km = 2 if(km.gt.5.and.km.lt.10) km = 5 if(km.gt.10.and.km.lt.20) km = 10 if(km.gt.20.and.km.lt.50) km = 20 iinc = km kinc = km + km/2 if(km.eq.50) iinc = 60 if(km.eq.1.or.km.eq.5.or.km.eq.50) kinc = iinc dscl = float(km)/ccm 11 call prompt(' Draw Scale Bar ? ') call gparma(' "y" (Yes) or "n" (No) ==> ', 1, yn) if(yn.ne.'y'.and.yn.ne.'n') goto 11 if(yn.eq.'y') then kscl = 1 write(out,'(a,f4.1,a)') * ' size :', dscl, 'cm Wide, 0.8cm High' call premsg(out) call prompt(' Lower-Left Pos. ? (right and up)') call gparmf2(' in cm ==> ', sclr, sclu) endif 12 call prompt(' Draw Lat.Lon.Lines ?') call gparma(' "y" (Yes) or "n" (No) ==> ', 1, yn) if(yn.ne.'y'.and.yn.ne.'n') goto 12 if(yn.eq.'y') klll = 1 c endif 13 call prompt(' Draw Coastlines ? ') call gparma(' "y" (Yes) or "n" (No) ==> ', 1, yn) if(yn.ne.'y'.and.yn.ne.'n') goto 13 if(yn.eq.'y') then kshl = 2 14 call prompt(' Rivers ? ') call gparma(' "y" (Yes) or "n" (No) ==> ', 1, yn) if(yn.ne.'y'.and.yn.ne.'n') goto 14 if(yn.eq.'y') kshl = 12 15 call prompt(' Pref.Boundary ?') call gparma(' "y" (Yes) or "n" (No) ==> ', 1, yn) if(yn.ne.'y'.and.yn.ne.'n') goto 15 if(yn.eq.'y') then kshl = kshl + 200 if(kshl.eq.202) kshl = 102 endif endif write(9,'(a8,i4,4x,4i8)') area, ncc, iorg, korg, ispa, ispb write(9,*) xs, ys, mszx, mszy, nmx, nmy write(9,*) istp, mval, az, el, vs, fscl, high, wide write(9,*) kscl, sclr, sclu, klll, kshl write(6,'(a)') '----------------------------------------' write(6,'(a)') 'PLMAPS : Plot contour map with shading' write(6,'(a,a)') ' Input data filename : ', fnin(1:lrtrim(fnin)) write(6,'(a,a)') ' Output PS filename : ', fnam(1:lrtrim(fnam)) write(6,'(a,a8,a,i4,4x,4i8)') * 'Areaname : ', area, ' Coord.', ncc, iorg,korg, ispa,ispb write(6,'(14x,2a10,2x,2a8,2a6)') * 'xs', 'ys', 'msz-x', 'msz-y', 'mx', 'my' write(6,'(16x,2f10.2,2i8,2i6)') xs, ys, mszx, mszy, nmx, nmy write(6,'(6x,a,i4,4x,a,i4)') 'step=', istp, 'lcapt=', lcap write(6,'(6x,a,f7.1,a,f7.1)') * 'illuminant : Azimuth', az, ', Elevation', el write(6,'(6x,a,f8.2)') ' value to be scaled to 1km :', vs write(6,'(6x,a,f8.0,4x,a3,f8.1,a1,f8.1,a1)') * 'scale 1 /', fscl, spl, high, 'H', wide, 'W' if(kscl.ne.0.or.klll.ne.0.or.kshl.ne.0) then write(6,'(a,$)') ' with ' if(kshl.ne.0) write(6,'(a,$)') ' Coast' if(mod(kshl,100)/10.ne.0) write(6,'(a,$)') '+River' if(kshl/100.ne.0) write(6,'(a,$)') '+PrefB' if(klll.ne.0) write(6,'(a,$)') ' LatLonLines' if(kscl.ne.0) then write(6,'(a,f5.1,$)') ' ScaleBar [', sclr write(6,'(a,f5.1,a,$)') ',', sclu, ' ]' endif write(6,'(a)') endif c do 50 i=1,nmx do 50 j=1,nmy n = (i-1)*nmy + j if (g(n).eq.vnul) then s(n) = 1. else if (i.eq.1.or.g(n-nmy).eq.vnul) then if (i.eq.nmx.or.g(n+nmy).eq.vnul) then dx = 0. else dx = g(n+nmy) - g(n) endif else if (i.eq.nmx.or.g(n+nmy).eq.vnul) then dx = g(n) - g(n-nmy) else dx = (g(n+nmy)-g(n-nmy)) / 2. endif endif if (j.eq.1.or.g(n-1).eq.vnul) then if (j.eq.nmy.or.g(n+1).eq.vnul) then dy = 0. else dy = g(n+1) - g(n) endif else if (j.eq.nmy.or.g(n+1).eq.vnul) then dy = g(n) - g(n-1) else dy = (g(n+1)-g(n-1)) / 2. endif endif p = dx/fx q = dy/fy r = sqrt(p*p + q*q + 1.) sn = (p*cx + q*cy - cz) / r c s(n) = (1.-sn) / 2. if (sn.lt.0.) then s(n) = 0.75 - sn/4. else s(n) = (1.-sn) * 0.75 endif endif 50 continue c call plots(1.5, 1.5) call dftone(-1) call dframe(0., 0., wide, high, nmy, nmx) call paintm(s, -1., +1., 9.9999) call wrect(0., 0., wide, high) call conts(0., 0., wide, high, nmy, nmx, 0.) call contx(g, vnul, istp, lcap, -99999, -1, 11) call contx(g, vnul, istp, lcap, 0, 99999, 10) if(kscl.ne.0) then call plot(sclr, sclu+0.1, 3) call plot(sclr+dscl, sclu+0.1, 2) call newpen(2) call plot(sclr+dscl, sclu+0.3, 3) call plot(sclr+dscl, sclu, 2) call plot(sclr, sclu, 2) call plot(sclr, sclu+0.3, 2) write(s4,'(i2,a2)') km, 'km' call lstyle('T', 0.5, 0., 0, -255) call ptext(s4, 4, sclr+dscl/2., sclu+0.5, 2) call newpen(1) endif call scisor(0., 0., wide, high) call cvinit(ncc,float(iorg),float(korg),float(ispa),float(ispb)) if(klll.ne.0.or.kshl.ne.0) then xt = xs + float((nmx-1)*mszx)/1000. yt = ys + float((nmy-1)*mszy)/1000. call cvenik(ys, xs, fi1, fk1) call cvenik(yt, xs, fi2, fk2) call cvenik(ys, xt, fi3, fk3) call cvenik(yt, xt, fi4, fk4) isf = nint(amin1(fi1,fi2) - 0.5) itf = nint(amax1(fi3,fi4) + 0.5) ksf = nint(amin1(fk1,fk3) - 0.5) ktf = nint(amax1(fk2,fk4) + 0.5) endif if(klll.ne.0) then do 200 i=isf/iinc,itf/iinc idf = i * iinc if(idf.lt.isf.or.idf.gt.itf) goto 200 call cviken(float(idf), float(ksf), y, x) call plot((y-ys)/ccm, (x-xs)/ccm, 3) do 100 kdf=ksf+1,ktf call cviken(float(idf), float(kdf), y, x) call plot((y-ys)/ccm, (x-xs)/ccm, 2) 100 continue 200 continue do 400 k=ksf/kinc,ktf/kinc kdf = k * kinc if(kdf.lt.ksf.or.kdf.gt.ktf) goto 400 call cviken(float(isf), float(kdf), y, x) call plot((y-ys)/ccm, (x-xs)/ccm, 3) do 300 idf=isf+1,itf call cviken(float(idf), float(kdf), y, x) call plot((y-ys)/ccm, (x-xs)/ccm, 2) 300 continue 400 continue endif if(kshl.ne.0) then call rshore(mc, isf/60, (itf+59)/60, ksf/60, (ktf+59)/60) call pshore(-ys/ccm, -xs/ccm, kshl, ccm, cviken) endif call scisor(0., 0., 0., 0.) call lstyle('T', 0.4, 0., 0, -255) call ptext(area, 8, -1., -1., 0) call ptext(fnin, 50, 2., -1., 0) call plote() 3 read(10,'(a)',end=8) buf if(buf(1:1).eq.'#') goto 3 20 call prompt(' Process Next set ?') call gparma(' "y" (Yes) or "n" (No) ==> ', 1, yn) if(yn.ne.'y'.and.yn.ne.'n') goto 20 if(yn.eq.'y') goto 2 c 8 close(10) endfile(9) rewind(9) call pscls() call clspin() c if(flog(ldr:ldr).ne.' ') then open(99,file=flog,access='append') write(99,'(//a)') 'PLMAPS : Plot contour map with shading' write(99,'(a,a)') * ' Input data filename : ', fnin(1:lrtrim(fnin)) write(99,'(a,a)') * ' Output PS filename : ', fnam(1:lrtrim(fnam)) 92 read(9,'(a8,i4,4x,4i8)',end=93) area, ncc, iorg,korg, ispa,ispb read(9,*) xs, ys, mszx, mszy, nmx, nmy read(9,*) istp, mval, az, el, vs, fscl, high, wide read(9,*) kscl, sclr, sclu, klll, kshl write(99,'(a,a8,a,i4,4x,4i8)') * 'Areaname : ', area, ' Coord.', ncc, iorg,korg,ispa,ispb write(99,'(14x,2a10,2x,2a8,2a6)') * 'xs', 'ys', 'msz-x', 'msz-y', 'mx', 'my' write(99,'(16x,2f10.2,2i8,2i6)') xs, ys, mszx, mszy, nmx, nmy write(99,'(6x,a,i4,4x,a,i4)') 'step=', istp, 'lcapt=', lcap write(99,'(6x,a,f7.1,a,f7.1)') * 'illuminant : Azimuth', az, ', Elevation', el write(99,'(6x,a,f8.2)') ' value to be scaled to 1km :', vs write(99,'(6x,a,f8.0,4x,a3,f8.1,a1,f8.1,a1)') * 'scale 1 /', fscl, spl, high, 'H', wide, 'W' if(kscl.ne.0.or.klll.ne.0.or.kshl.ne.0) then write(99,'(a,$)') ' with ' if(kshl.ne.0) write(99,'(a,$)') ' Coast' if(mod(kshl,100)/10.ne.0) write(99,'(a,$)') '+River' if(kshl/100.ne.0) write(99,'(a,$)') '+PrefB' if(klll.ne.0) write(99,'(a,$)') ' LatLonLines' if(kscl.ne.0) then write(99,'(a,f5.1,$)') ' ScaleBar [', sclr write(99,'(a,f5.1,a,$)') ',', sclu, ' ]' endif write(99,'(a)') endif goto 92 93 call strdtm(sdtm) write(99,'(a,a,a)') '=============== ', sdtm, '===============' close(99) endif close(9) stop end