!--------( GGRID : gen. grid from DPAM line data or equiv.)-------- !-------- This program must be linked with C subroutine implicit none integer :: lwkdir, lrtrim real(8) :: v(30) character(80) :: wdr, flog, fnam1, fnam2 character(150) :: buf; character(20) :: sdtm; character(9) :: strim character(8) :: area; character(6) :: cu; character(3) :: cg integer :: ldr, ktyp, kg, ku, k1, k2, k3, kmx, i, lc, IER, i1, i2 integer :: nc, ncc, iorg, korg, ispa, ispb, ixl, iyl, mszx, mszy, mx, my real :: sr, xl, yl, smx, smy, fmag real(8) :: dxn, dye, vlat, vlon, tlat, tlon, talt, wlat, wlon, walt call premsg('GGRID : gen. grid from DPAM line data or equiv.') call opnpin() ldr = lwkdir(50,wdr) + 1 flog = wdr; fnam1 = wdr; fnam2 = wdr call gparma('Enter LOG filename ==> ', 81-ldr, flog(ldr:80)) call premsg('Select AM Line data type :') call premsg(' DPAM Line data (1), AMDB-GSJ Line data (2),') call premsg(' AMDB-NEDO Line data (3), StdLIN data (4),') call gparmi(' or general Random data (0) ==> ', ktyp) if ((ktyp < 0) .or. (ktyp > 4)) call abendm('invalid value') call gparma('Enter Input line/random data filename ==> ', & 81-ldr, fnam1(ldr:80)) open(1,file=fnam1,status='old') call gparmi(' Lat./Lon. ? in WGS(1) or TokyoDatum(2) ==> ', kg) if (kg == 1) then cg = 'WGS' else if (kg == 2) then cg = 'TYO' else call abendm('invalid value') endif ku = 1; cu = ' Min. ' if (ktyp == 1) then k1 = 5; k2 = 6; k3 = 9; kmx = 9; ku = 2; cu = ' Deg. ' else if (ktyp == 2) then k1 = 2; k2 = 3; k3 = 4; kmx = 4 else if (ktyp == 3) then k1 = 3; k2 = 4; k3 = 8; kmx = 8 read(1,'(a)') buf if (buf(1:8) /= '## NEDO') call abendm('seems non AMDB-NEDO data') else if (ktyp == 4) then k1 = 1; k2 = 2; k3 = 4; kmx = 4 else do read(1,'(a)') buf call premsg(buf(1:lrtrim(buf))) if ((buf(1:1) == '&') .or. (buf(1:1) == '%')) cycle if (buf(1:1) /= '#') exit enddo rewind(1) call prompt('LatLon unit ? ') call gparmi(' in Min.(1), Deg.(2) or D:Min.(3) ==> ', ku) if (ku == 1) then cu = ' Min. ' else if (ku == 2) then cu = ' Deg. ' else if (ku == 3) then cu = 'D:Min.' else call abendm('invalid value') endif call gparmi2(' data positions of Lat,Lon : ', k1,k2) if ((k1 <= 0) .or. (k2 <= 0) .or. (k2==k1)) call abendm('invalid value') call gparmi(' data position of Mag. data : ', k3) if ((k3 <= 0) .or. (k3==k1) .or. (k3==k2)) call abendm('invalid value') kmx = max(k1, k2, k3) if (ku == 3) then if ((k1+1)==k2 .or. (k1+1)==k3 .or. (k2+1)==k1 .or. (k2+1)==k3) & call abendm('invalid value') kmx = max(k1+1, k2+1, k3) endif if (kmx > 30) call abendm('too large data position') endif call premsg('Specify trim radius in km') call gparmf(' ( 0 for no-trimming ) ==> ', sr) if (sr < 0.) call abendm('invalid value') write(strim,'(f6.2,1x,a2)') sr, 'km' if (sr == 0.) strim = ' no-trim' call gparma('Enter Name-label for grid-data ==> ', 8, area) call premsg('Specify gridding parameters') call gparmi(' Coordinate Number (+800 for Old-TOKYO) : ', ncc) nc = ncc if (ncc >= 800) nc = ncc - 800 ispa = 0; ispb = 0 if ((nc >= 1) .and. (nc <= 62)) then iorg = 0 korg = (nc-30)*360 - 180 if (nc > 60) korg = 0 if (nc == 61) iorg = +5400 if (nc == 62) iorg = -5400 else if (nc == 65) then iorg = 0 call gparmi2(' Central Longitude ? (Deg. Min.) : ', i1,i2) korg = i1*60 + i2 else if (nc == 72) call abendm('nc=72 not implemented') call gparmi2(' Origin Latitude ? (Deg. Min.) : ', i1,i2) iorg = i1*60 + i2 call gparmi2(' Origin Longitude ? (Deg. Min.) : ', i1,i2) korg = i1*60 + i2 endif call premsg(' Southwest corner Coord. in km ?') call gparmf(' Northing : ', xl) ixl = nint(xl*1000.) call gparmf(' Easting : ', yl) iyl = nint(yl*1000.) call premsg(' Mesh size and numbers to N ?') call gparmi(' Unit Mesh size in m : ', mszx) if (mszx <= 0) call abendm('invalid value') smx = float(mszx) / 1000. call gparmi(' Number of Mesh to N : ', mx) if (mx <= 0) call abendm('invalid value') call premsg(' Mesh size and numbers to E ?') call gparmi(' Unit Mesh size in m : ', mszy) if (mszy <= 0) call abendm('invalid value') smy = float(mszy) / 1000. call gparmi(' Number of Mesh to E : ', my) if (my <= 0) call abendm('invalid value') call gparma('Enter grid-data output filename ==> ', 81-ldr, fnam2(ldr:80)) call clspin() call cvinit(ncc,float(iorg),float(korg),float(ispa),float(ispb)) open(2,file=fnam2,status='new') write(6,'(a)') '----------------------------------------' write(6,'(a)') 'GGRID : gen. grid from DPAM line data or equiv.' write(6,'(a,a)') ' Input line/random data file : ', fnam1(1:lrtrim(fnam1)) write(6,'(8x,i2,5a,i2,5a,i2,a)') k1, ': ', cg, '-Lat.(', cu, ') ', & k2, ': ', cg, '-Lon.(', cu, ') ', k3, ': Mag.(nT)' write(6,'(a,a)') ' trim radius :', strim write(6,'(a,a)') ' Gen. grid output filename : ', fnam2(1:lrtrim(fnam2)) write(6,'(a12,a8,6x,a6,i4,4x,4i8)') & 'Areaname : ', area, 'Coord.', ncc, iorg, korg, ispa, ispb write(6,'(4x,a)') & ' (mesh) xs ys msz-x msz-y mx my alt ' write(6,'(4x,a12,2f10.2,2i8,2i6,a8)') & 'gen. grid ', xl, yl, mszx, mszy, mx, my, ' undef' write(2,'(a)') '# output from GGRID' write(2,'(a,a)') '# input line/random data file : ', fnam1(1:lrtrim(fnam1)) write(2,'(a8,4x,i4,4i8)') area, ncc, iorg,korg, ispa,ispb write(2,'(2i12,4i6,1x,f7.1,1x,f7.0)') ixl,iyl, mszx,mszy, mx,my, 9999.9, -1. close(2) call gsurf0(xl, yl, smx, smy, mx, my) do read(1,'(a)',iostat=IER) buf if (IER < 0) exit if ((buf(1:1)=='&') .or. (buf(1:1)=='%') .or. (buf(1:1)=='#')) then call premsg(buf(1:lrtrim(buf))); cycle endif if (ktyp == 4) then i = index(buf,'N') if (i /= 0) buf(i:i) = ' ' i = index(buf,'E') if (i /= 0) buf(i:i) = ' ' i = index(buf,'m') if (i /= 0) buf(i:i) = ' ' i = index(buf,'nT') if (i /= 0) buf(i:i+1) = ' ' read(buf,*) (v(i),i=1,4) else if (ktyp == 0) then do lc = index(buf,':') if (lc == 0) exit buf(lc:lc) = ' ' enddo endif read(buf,*) (v(i),i=1,kmx) endif fmag = real(v(k3)) if (ku == 1) then vlat = v(k1); vlon = v(k2) else if (ku == 2) then vlat = v(k1) * 60.; vlon = v(k2) * 60. else vlat = v(k1)*60. + v(k1+1); vlon = v(k2)*60. + v(k2+1) endif if ((kg == 1) .and. (ncc >= 800)) then call xw84td(vlat, vlon, 0.d0, tlat, tlon, talt) call cvdiken(tlat, tlon, dye, dxn) else if ((kg == 2) .and. (ncc < 800)) then call xtw84d(vlat, vlon, 0.d0, wlat, wlon, walt) call cvdiken(wlat, wlon, dye, dxn) else call cvdiken(vlat, vlon, dye, dxn) endif call gsurf2(1, sngl(dxn), sngl(dye), fmag) enddo close(1) call gsurf2(0, 0., 0., 0.) call gsurf3(sr, 9999.9) call gsurf4(fnam2) call gsurf5() if (flog(ldr:ldr) /= ' ') then open(99,file=flog,access='append') write(99,'(//a)') 'GGRID : gen. grid from DPAM line data or equiv.' write(99,'(a,a)') & ' Input line/random data file : ', fnam1(1:lrtrim(fnam1)) write(99,'(8x,i2,5a,i2,5a,i2,a)') k1, ': ', cg, '-Lat.(', cu, ') ', & k2, ': ', cg, '-Lon.(', cu, ') ', k3, ': Mag.(nT)' write(99,'(a,a)') ' trim radius :', strim write(99,'(a,a)') & ' Surface grid output filename : ', fnam2(1:lrtrim(fnam2)) write(99,'(a12,a8,6x,a6,i4,4x,4i8)') & 'Areaname : ', area, 'Coord.', ncc, iorg, korg, ispa, ispb write(99,'(4x,a)') & ' (mesh) xs ys msz-x msz-y mx my alt ' write(99,'(4x,a12,2f10.2,2i8,2i6,a8)') & 'gen. grid ', xl, yl, mszx, mszy, mx, my, ' undef' call strdtm(sdtm) write(99,'(a,a,a)') '=============== ', sdtm, '===============' close(99) endif stop end