SMS模型格网转换为MIKE21的格网源代码

    program main !sms网格转换成mike21网格
    DIMENSION X(60000),Y(60000),H(60000),NDNN(110000,3),ncbd(10000)
    dimension NBS(50),NOBD(5000,50),NSED(5000,50),id0(80000),nbi(50)
    dimension x0(50000),y0(50000),h0(50000)
    character*2 a*2,b*2,gnn*3,demo*500


    open(3,file='./2013.8.30(after).grd',status='old')
    read(3,*)
    read(3,*)nel,nde
    do i=1,nde
    read(3,*)num,x(i),y(i),h(i)
    h(i)=h(i)+2.2
    enddo
    do i=1,nel
    read(3,*)num,nd,(ndnn(i,j),j=1,3)
    enddo

    do i=1,4 !
    read(3,*)
    enddo
    read(3,*)kb
    read(3,*)
    do k=1,kb
    read(3,*)nbs(k)
    do i=1,nbs(k)
    read(3,*)nobd(i,k)
    enddo
    enddo
    close(3)
c-----------------------------------------------------------
      open(4,file='滑道1.grd',status='old')
    read(4,*)
    read(4,*)nel0,nde0
    do i=1,nde0
    read(4,*)num,x0(i),y0(i),h0(i)
    enddo
    close(4)
    
    do 25 i=1,nde
    do j=1,nde0
    dis=sqrt((x(i)-x0(j))**2+(y(i)-y0(j))**2)
    if(dis.le.5)then
    h(i)=2.2-h0(j)
    goto 25
    endif
    enddo
25    continue    
     
      
      open(4,file='滑道2.grd',status='old')
    read(4,*)
    read(4,*)nel0,nde0
    do i=1,nde0
    read(4,*)num,x0(i),y0(i),h0(i)
    enddo
    close(4)
    
    do 35 i=1,nde
    do j=1,nde0
    dis=sqrt((x(i)-x0(j))**2+(y(i)-y0(j))**2)
    if(dis.le.5)then
    h(i)=2.2-h0(j)
    goto 35
    endif
    enddo
35    continue    
c----------------------------------------------------------
    id0=0


    open(18,file='closebd.txt',status='old')
    do nn=1,100
    read(18,*,end=181)nk
    write(*,*)nk
    read(18,*)(ncbd(n),n=1,nk)
    do i=1,nk
    id0(ncbd(i))=1 !+nn
    enddo
    enddo
181    continue
    close(18)
    do k=1,kb


    do i=1,nbs(k)
    id0(nobd(i,k))=k+1
    enddo
    enddo

    write(*,*)


    demo=' PROJCS["Beijing_1954_3_Degree_GK_CM_111E",GEOGCS
     *["GCS_Beijing_1954",DATUM["D_Beijing_1954",SPHEROID["Krasovsky
     *_1940",6378245.0,298.3]],PRIMEM["Greenwich",0.0],UNIT["Degree",
     *0.0174532925199433]],PROJECTION["Gauss_Kruger"],PARAMETER["False
     *_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER
     *["Central_Meridian",120.0],PARAMETER["Scale_Factor",1.0],
     *PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]'
  
    open(30,file='2013.8.30(after).mesh')
    write(30,'(I8,a)')nde,demo
    d0=0
    do i=1,nde
    write(30,'(i8,2f12.2,f10.4,I4)')i,x(i),y(i),h(i)*-1,id0(i)  !
    enddo
    write(30,*)nel,'  3  21'
    do i=1,nel
    write(30,'(5I6)')i,(ndnn(i,j),j=1,3),d0
    enddo


    close(30)


    end