c--------( GOJOIN : overlay-join multiple grid data )-------- parameter (LX=4501, LY=4501) c dimension f(LX,LY), g(LX,LY), r(LX,LY) dimension a(LX,LY), h(LX,LY), t(LX,LY) character*80 wdr, flog, fnam1, fnam2 character area*8, narea*8, area1*8, anul*10, salt*8 character buf*80, sdtm*20 c call premsg('GOJOIN : overlay-join multiple grid data') 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('') call gparma('Enter Name-label for new grid ==> ', 8, narea) call gparmi('Select Projection-Number ==> ', nccn) ncn = nccn if(nccn.ge.200) ncn = nccn - 200 ispan = 0 ispbn = 0 if(ncn.ge.1.and.ncn.le.62) then iorgn = 0 korgn = (ncn-30)*360 - 180 if(ncn.gt.60) korgn = 0 if(ncn.eq.61) iorgn = +5400 if(ncn.eq.62) iorgn = -5400 else if(ncn.eq.65) then iorgn = 0 call gparmi(' Central Longitude ? (Deg.) : ', i1) korgn = i1*60 else call gparmi2(' Origin Latitude ? (Deg. Min.) : ', i1,i2) iorgn = i1*60 + i2 call gparmi2(' Origin Longitude ? (Deg. Min.) : ', i1,i2) korgn = i1*60 + i2 if(ncn.eq.72) then call gparmi2(' St.Parallel-1 Lat. ? (Deg. Min.) : ', i1,i2) ispan = i1*60 + i2 call gparmi2(' St.Parallel-2 Lat. ? (Deg. Min.) : ', i1,i2) ispbn = i1*60 + i2 endif endif call premsg('Specify gridding parameters') call premsg(' Southwest corner Coord. in km ?') call gparmf(' Northing : ', xl) ixl = nint(xl*1000.) call gparmf(' Easting : ', yl) iyl = nint(yl*1000.) call gparmi(' Unit Mesh size in m : ', msz) if(msz.le.0) call abendm('invalid value') call gparmi(' Number of Mesh to N : ', mxn) if(mxn.gt.LX) call abendm('too large array size') call gparmi(' Number of Mesh to E : ', myn) if(myn.gt.LY) call abendm('too large array size') vnuln = 9999.9 call gparmf('Boundary zone Width ? (in km) ==> ', sb) nb = nint(sb*1000./float(msz)) mxnb = mxn + nb + nb mynb = myn + nb + nb if(mxnb.gt.LX .or. mynb.gt.LY) * call abendm('too large array size') dnorm = float(nb) nddx = nb*nb*4 call gparma('Enter new grid output filename ==> ', * 81-ldr, fnam2(ldr:80)) open(2,file=fnam2,status='new') do 100 j=1,mynb do 100 i=1,mxnb g(i,j) = vnuln h(i,j) = vnuln 100 continue write(6,'(a)') '----------------------------------------' write(6,'(a)') 'GOJOIN : overlay-join multiple grid data' write(6,'(a,a)') ' Output new file : ', fnam2 write(6,'(a,f5.1,a)') ' Boundary zone Width :', sb, 'km' write(6,1000) narea, nccn, iorgn,korgn, ispan,ispbn, * xl, yl, msz, msz, mxn, myn, ' var.' call strdtm(sdtm) open(9,status='scratch') c nsrc = 0 1 nsrc = nsrc + 1 call premsg('') call premsg('Enter Input source data filename') call gparma(' (Hit RETURN to End) ==> ', * 81-ldr, fnam1(ldr:80)) if(fnam1(ldr:ldr).eq.' ') goto 8 open(1,file=fnam1,status='old') 2 read(1,'(a)') buf if(buf(1:1).eq.'#') goto 2 read(buf,'(a8,i4,4x,4i8)') area, ncc, iorg,korg, ispa,ispb nc = ncc if(ncc.ge.200) nc = ncc - 200 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 if(ncc.ne.nccn .or. * iorg.ne.iorgn.or.korg.ne.korgn .or. * ispa.ne.ispan.or.ispb.ne.ispbn) * call abendm('Header inconsistent') read(1,'(a)') buf read(buf,*) ixs, iys, mszx, mszy, mx, my, vnul, alt write(salt,'(1x,f7.0)') alt if (alt.eq.0.) salt = ' var.' if (alt.lt.0.) salt = ' undef' if (mszx.ne.msz.or.mszy.ne.msz) call abendm('illegal mesh-size') is = (ixs-ixl)/mszx js = (iys-iyl)/mszy if (is*mszx+ixl.ne.ixs .or. * js*mszy+iyl.ne.iys) call abendm('ranges conflict') if(mx.gt.LX.or.my.gt.LY) call abendm('too large array size') xs = float(ixs)/1000. ys = float(iys)/1000. read(1,*) ((r(i,j),i=1,mx),j=1,my) if(alt.eq.0.) goto 3 if(alt.lt.0.) then call gparmf(' Altitude [in m] for this data ==> ', alt) if (alt.le.0.) call abendm('invalid value') write(salt,'(1x,f7.0)') alt endif goto 4 3 read(1,'(a)',end=4) buf if(buf(1:1).eq.'#') goto 3 read(buf,'(a8,i4,4x,4i8)') area1, ncc1, iorg1,korg1, ispa1,ispb1 nc1 = ncc1 if(ncc1.ge.200) nc1 = ncc1 - 200 if(nc1.ge.1.and.nc1.le.62) then iorg1 = 0 ispa1 = 0 ispb1 = 0 korg1 = (nc1-30)*360 - 180 if(nc1.gt.60) korg1 = 0 if(nc1.eq.61) iorg1 = +5400 if(nc1.eq.62) iorg1 = -5400 endif if(area1.ne.area.or.ncc1.ne.ncc .or. * iorg1.ne.iorg.or.korg1.ne.korg .or. * ispa1.ne.ispa.or.ispb1.ne.ispb) * call abendm('Header1 inconsistent') read(1,*) ixs1, iys1, mszx1, mszy1, mx1, my1, vnul1 if(ixs1.ne.ixs.or.iys1.ne.iys .or. * mszx1.ne.mszx.or.mszy1.ne.mszy .or. * mx1.ne.mx.or.my1.ne.my) * call abendm('Header2 inconsistent') write(anul,'(f10.3)') vnul1 if(anul(9:10).ne.'00'.or.anul(1:1).ne.' ') * call abendm('unable to handle vnul1 value') read(1,*) ((t(i,j),i=1,mx),j=1,my) 4 close(1) write(6,'(a,a,a8,a,i4,4x,4i8,a7)') '--- ', * '--- Area : ', area, ' Coord.', ncc, iorg,korg, ispa,ispb write(6,'(4x,a,2a10,2x,2a8,2a6)') * ' (mesh) ', 'xs', 'ys', 'msz-x', 'msz-y', 'mx', 'my', * 'alt' write(6,'(4x,a,2f10.2,2i8,2i6,a8)') * 'source data ', xs, ys, mszx, mszy, mx, my, salt write(9,'(a8,i4,4x,4i8)') area, ncc, iorg,korg, ispa,ispb write(9,*) xs, ys, mszx, mszy, mx, my, alt do 200 j=1,my jp = js + j + nb if (jp.le.0.or.jp.gt.mynb) goto 200 do 210 i=1,mx if (r(i,j).eq.vnul) goto 210 ip = is + i + nb if (ip.le.0.or.ip.gt.mxnb) goto 210 g(ip,jp) = r(i,j) if (alt.ne.0.) then h(ip,jp) = alt else if (t(i,j).ne.vnul1) then h(ip,jp) = t(i,j) else write(buf,'(a30,f8.1,a8,i4,a4,i4)') * ' t(i,j)=vnul where r(i,j) =', r(i,j), * ' at i=', i, ' j=', j call premsg(buf(1:58)) call abendm('Aborted') endif 210 continue 200 continue if (nsrc.eq.1) then do 300 j=1,mynb do 300 i=1,mxnb f(i,j) = g(i,j) a(i,j) = h(i,j) g(i,j) = vnuln h(i,j) = vnuln 300 continue else call dpcini('... Gojoin processing : ') do 400 j=1,mynb lf = 0 lg = 0 do 400 i=1,mxnb ndd = nddx if (lf.lt.0) then do 500 jm=-nb,+nb if ((i+nb).gt.mxnb) goto 500 if ((j+jm).le.0.or.(j+jm).gt.mynb) goto 500 if (f(i+nb,j+jm).ne.vnuln) then nddp = nb*nb + jm*jm if (nddp.lt.ndd) ndd = nddp endif 500 continue if (ndd.ne.nddx) then lf = 0 df = -sqrt(float(ndd)) / dnorm endif else if (lf.gt.0) then do 510 jm=-nb,+nb if ((i+nb).gt.mxnb) goto 510 if ((j+jm).le.0.or.(j+jm).gt.mynb) goto 510 if (f(i+nb,j+jm).eq.vnuln) then nddp = nb*nb + jm*jm if (nddp.lt.ndd) ndd = nddp endif 510 continue if (ndd.ne.nddx) then lf = 0 df = sqrt(float(ndd-1)) / dnorm endif else if (f(i,j).eq.vnuln) then do 600 jm=-nb,+nb do 600 im=-nb,+nb if ((i+im).le.0.or.(i+im).gt.mxnb) goto 600 if ((j+jm).le.0.or.(j+jm).gt.mynb) goto 600 if (f(i+im,j+jm).ne.vnuln) then nddp = im*im + jm*jm if (nddp.lt.ndd) ndd = nddp endif 600 continue if (ndd.eq.nddx) then lf = -1 df = -2. else df = -sqrt(float(ndd)) / dnorm endif else do 610 jm=-nb,+nb do 610 im=-nb,+nb if ((i+im).le.0.or.(i+im).gt.mxnb) goto 610 if ((j+jm).le.0.or.(j+jm).gt.mynb) goto 610 if (f(i+im,j+jm).eq.vnuln) then nddp = im*im + jm*jm if (nddp.lt.ndd) ndd = nddp endif 610 continue if (ndd.eq.nddx) then lf = +1 df = +2. else df = sqrt(float(ndd-1)) / dnorm endif endif endif ndd = nddx if (lg.lt.0) then do 700 jm=-nb,+nb if ((i+nb).gt.mxnb) goto 700 if ((j+jm).le.0.or.(j+jm).gt.mynb) goto 700 if (g(i+nb,j+jm).ne.vnuln) then nddp = nb*nb + jm*jm if (nddp.lt.ndd) ndd = nddp endif 700 continue if (ndd.ne.nddx) then lg = 0 dg = -sqrt(float(ndd)) / dnorm endif else if (lg.gt.0) then do 710 jm=-nb,+nb if ((i+nb).gt.mxnb) goto 710 if ((j+jm).le.0.or.(j+jm).gt.mynb) goto 710 if (g(i+nb,j+jm).eq.vnuln) then nddp = nb*nb + jm*jm if (nddp.lt.ndd) ndd = nddp endif 710 continue if (ndd.ne.nddx) then lg = 0 dg = sqrt(float(ndd-1)) / dnorm endif else if (g(i,j).eq.vnuln) then do 800 jm=-nb,+nb do 800 im=-nb,+nb if ((i+im).le.0.or.(i+im).gt.mxnb) goto 800 if ((j+jm).le.0.or.(j+jm).gt.mynb) goto 800 if (g(i+im,j+jm).ne.vnuln) then nddp = im*im + jm*jm if (nddp.lt.ndd) ndd = nddp endif 800 continue if (ndd.eq.nddx) then lg = -1 dg = -2. else dg = -sqrt(float(ndd)) / dnorm endif else do 810 jm=-nb,+nb do 810 im=-nb,+nb if ((i+im).le.0.or.(i+im).gt.mxnb) goto 810 if ((j+jm).le.0.or.(j+jm).gt.mynb) goto 810 if (g(i+im,j+jm).eq.vnuln) then nddp = im*im + jm*jm if (nddp.lt.ndd) ndd = nddp endif 810 continue if (ndd.eq.nddx) then lg = +1 dg = +2. else dg = sqrt(float(ndd-1)) / dnorm endif endif endif if (dg.lt.0.) then if (df.lt.0.) then if ((1.+df+dg).gt.0.) then avf = avrg(f,i,j,mxnb,mynb,vnuln, (1.+dg-df)*dnorm) avg = avrg(g,i,j,mxnb,mynb,vnuln, -dg*2.*dnorm) r(i,j) = (avf*(1.+df) + avg*(1.+dg)) / (2.+df+dg) ava = avrg(a,i,j,mxnb,mynb,vnuln, (1.+dg-df)*dnorm) avh = avrg(h,i,j,mxnb,mynb,vnuln, -dg*2.*dnorm) t(i,j) = (ava*(1.+df) + avh*(1.+dg)) / (2.+df+dg) else r(i,j) = vnuln t(i,j) = vnuln endif else if (df.lt.1.) then if ((1.+dg-df).gt.0.) then avf = avrg(f,i,j,mxnb,mynb,vnuln, (1.+dg-df)*dnorm) avg = avrg(g,i,j,mxnb,mynb,vnuln, -dg*2.*dnorm) r(i,j) = (avf + avg*(1.+dg-df)) / (2.+dg-df) ava = avrg(a,i,j,mxnb,mynb,vnuln, (1.+dg-df)*dnorm) avh = avrg(h,i,j,mxnb,mynb,vnuln, -dg*2.*dnorm) t(i,j) = (ava + avh*(1.+dg-df)) / (2.+dg-df) else r(i,j) = f(i,j) t(i,j) = a(i,j) endif else r(i,j) = f(i,j) t(i,j) = a(i,j) endif else if (dg.lt.1.) then if (df.lt.0.) then if ((1.+df-dg).gt.0.) then avf = avrg(f,i,j,mxnb,mynb,vnuln, (1.-df)*dnorm) r(i,j) = (avf*(1.+df-dg) + g(i,j)*(1.+dg)) / (2.+df) ava = avrg(a,i,j,mxnb,mynb,vnuln, (1.-df)*dnorm) t(i,j) = (ava*(1.+df-dg) + h(i,j)*(1.+dg)) / (2.+df) else r(i,j) = g(i,j) t(i,j) = h(i,j) endif else if (df.lt.1.) then avf = avrg(f,i,j,mxnb,mynb,vnuln, (1.-df)*dnorm) r(i,j) = (avf*(1.-dg) + g(i,j)*(1.+dg-df)) / (2.-df) ava = avrg(a,i,j,mxnb,mynb,vnuln, (1.-df)*dnorm) t(i,j) = (ava*(1.-dg) + h(i,j)*(1.+dg-df)) / (2.-df) else r(i,j) = f(i,j)*(1.-dg) + g(i,j)*dg t(i,j) = a(i,j)*(1.-dg) + h(i,j)*dg endif else r(i,j) = g(i,j) t(i,j) = h(i,j) endif call dpcent((j-1)*mxnb+i, mxnb*mynb) 400 continue do 900 j=1,mynb do 900 i=1,mxnb f(i,j) = r(i,j) a(i,j) = t(i,j) g(i,j) = vnuln h(i,j) = vnuln 900 continue endif goto 1 c 8 call clspin() endfile(9) rewind(9) write(2,'(a)') '# prog.GOJOIN output' write(2,'(a8,i4,4x,4i8)') narea, nccn, iorgn,korgn, ispan,ispbn write(2,'(2i12,4i6,1x,f7.1,1x,f7.0)') * ixl,iyl, mszx,mszy, mxn,myn, vnuln, 0. do 50 j=1,myn write(2,'((f7.1,9(1x,f7.1)))') (f(i+nb,j+nb),i=1,mxn) 50 continue write(2,'(a8,i4,4x,4i8)') narea, nccn, iorgn,korgn, ispan,ispbn write(2,'(2i12,4i6,1x,f7.1,1x,f7.0)') * ixl,iyl, mszx,mszy, mxn,myn, vnuln, -1. do 60 j=1,myn write(2,'((f7.1,9(1x,f7.1)))') (a(i+nb,j+nb),i=1,mxn) 60 continue close(2) c if(flog(ldr:ldr).ne.' ') then open(99,file=flog,access='append') write(99,'(//a)') 'GOJOIN : overlay-join multiple grid data' write(99,'(a,a)') ' Output new file : ', fnam2 write(99,'(a,f5.1,a)') ' Boundary zone Width :', sb, 'km' write(99,1000) narea, nccn, iorgn,korgn, ispan,ispbn, * xl, yl, msz, msz, mxn, myn, ' var.' write(99,'(a,a,a)') '--------------- ', sdtm, '---------------' 92 read(9,'(a8,i4,4x,4i8)',end=93) * area, ncc, iorg, korg, ispa, ispb read(9,*) xs, ys, mszx, mszy, mx, my, alt write(salt,'(1x,f7.0)') alt write(99,'(a,a,a8,a,i4,4x,4i8)') '--- ', * 'Area : ', area, ' Coord.', ncc, iorg,korg, ispa,ispb write(99,'(4x,a,2a10,2x,2a8,2a6,a7)') * ' (mesh) ', 'xs', 'ys', 'msz-x', 'msz-y', 'mx', 'my', * 'alt' write(99,'(4x,a,2f10.2,2i8,2i6,a8)') * 'source data ', xs, ys, mszx, mszy, mx, my, salt goto 92 93 call strdtm(sdtm) write(99,'(a,a,a)') '=============== ', sdtm, '===============' close(99) endif close(9) c stop 1000 format('Areaname : ', a8, 5x, 'Coord.', i4, 4x, 4i8/ * 3x, ' (mesh) ', ' xs ', ' ys ', * ' msz-x', ' msz-y', ' mx', ' my', ' alt '/ * 3x, 'gen.new grid', 2f10.2, 2i8, 2i6, a8) end c function avrg(f,i,j,mx,my,vnul, d) parameter (LX=4501, LY=4501) dimension f(LX,LY) dd = d * d n = nint(d) s = 0. k = 0 do 10 jm=-n,+n if((j+jm).le.0.or.(j+jm).gt.my) goto 10 do 20 im=-n,+n if((i+im).le.0.or.(i+im).gt.mx) goto 20 if (f(i+im,j+jm).eq.vnul) goto 20 ddp = float(im*im + jm*jm) if (ddp.gt.dd) goto 20 s = s + f(i+im,j+jm) k = k + 1 20 continue 10 continue if (k.eq.0) call abendm('Error in function subprogram (argv)') avrg = s / float(k) return end