付録

 

 

 

 

 

 

 

 

 

 

 

 

1.1       標高データ作成プログラム(1600×1600

 

1.2       標高データ作成プログラム(4800×4800

 

2  WGRIBを用いた外部気象データの作成

 

 

 

 

 

1.1 標高データ作成プログラム(1600×1600

使用方法(Windows,UNIX

1.プログラム文中のファイルパスを使用に応じてかえる。

2.CD『数値地図50mメッシュ(標高)』を入れる。

3.標高地域メッシュコード一覧を参考に地域番号を入力する(4桁)。

4.確認データを作るか作らないかを選択する。

5.直接使うならファイル名を変更、4800×4800に使うならファイル名はそのままで。

***** 標高データをCReSSの使用形式に変換するプログラム(1600×1600) ******

 

      character mn1*4, mn2*2

      integer  i,j,m,n

      real a(1600,1600), e(200,200),s

      write(*,*)'地域番号(4桁入れてください。)'

      read(*,*)mn1

************************************************************************

***** 番号指定 *********************************************************

      do i=0,77

       j=i

       if(j.eq. 8.or.j.eq. 9)go to 2

       if(j.eq.18.or.j.eq.19)go to 2

       if(j.eq.28.or.j.eq.29)go to 2

       if(j.eq.38.or.j.eq.39)go to 2

       if(j.eq.48.or.j.eq.49)go to 2

       if(j.eq.58.or.j.eq.59)go to 2

       if(j.eq.68.or.j.eq.69)go to 2

               write(mn2,'(I2.2)') j

       call input(mn1,mn2,e)

       call import(a,e,i)

    2 end do

   

      do n=1,1600

       do m=1,1600

      if(a(n,m).lt.0)go to 3

      go to 4

    3  a(n,m)=-1

    4  end do

      end do   

       open(unit=7,file='c:\Fortran\'//mn1//'.bin',access='direct'

     +                  ,form='unformatted',recl=6400,err=952)

       do n=1,1600

         write(7,rec=n,err=951)(a(n,m)/10,m=1,1600)

       end do

       close(unit=7)

 

      go to 600

   

     

************************************************************************

  951 write(*,*)'出力エラー'

      go to 999

  952 write(*,*)'最後のエラー'

      go to 999

************************************************************************

  600 write(*,*)'どちらかを番号により選択してください。'

      write(*,*)'1:確認ファイルを作成してから終了'

      write(*,*)'2:すぐに終了'

      read(5,*)s

      if(s.eq.1)call kakunin(mn1)

      if(s.ne.1)go to 999

      go to 999

     

  999 end

 

***** 標高データを読み込むプログラム(input) ****************************

      subroutine input(mn1,mn2,e)

      integer   c(200),d(200)

      integer   xdeg, xmin, ydeg, ymin

      character line*1011,fname*1024

      character mn1*4 , mn2*2

      real      e(200,200)

      real      xsec, ysec

      fname = 'd:\data\'//mn1//'\'//mn1//mn2//'.mem'

      open(1, file=fname, access='direct',err=12,

     +                 form='unformatted',  recl=1011)

      read(1,rec=1) line

      read(line(766:773), '(I3,I2,F3.1)') ydeg, ymin, ysec

      read(line(774:781), '(I3,I2,F3.1)') xdeg, xmin, xsec

      do 11 m=1,200

      read(1,rec=m+1)  line(1:1009)

      read(line, '(i6,i3,200f5.0)')c(m),d(m),(e(m,n),n=1,200)

   11 continue

      go to 13

   12 do m=1,200

       do n=1,200

         e(m,n)=-1

       end do

      end do

   13 close(1)

      return

      end

************************************************************************

 

***** 標高データを配列に代入するプログラム(import) *********************

      subroutine import(a,e,n)

      real a(1600,1600), e(200,200)

      integer i,j,n,nn,x,y

       if(n.ge. 0.and.n.le. 7)nn=0

       if(n.ge.10.and.n.le.17)nn=1

       if(n.ge.20.and.n.le.27)nn=2

       if(n.ge.30.and.n.le.37)nn=3

       if(n.ge.40.and.n.le.47)nn=4

       if(n.ge.50.and.n.le.57)nn=5

       if(n.ge.60.and.n.le.67)nn=6

       if(n.ge.70.and.n.le.77)nn=7

       x=n-10*nn

       y=nn

      do i=1,200

       do j=1,200

        a(j+200*y,i+200*x)=e(201-j,i)

       end do

      end do

      return

      end

************************************************************************

 

***** 確認プログラム ***************************************************

      subroutine kakunin(mn)

      integer i, j

      real    k(1600,1600)

      character mn*4

      open(unit=77,file='c:\Fortran\'//mn//'.bin',access='direct',

     +                 form='unformatted',  recl=6400,err=996)

 

      do 53 i=1,1600

      read(77,rec=i)(k(i,j),j=1,1600)  

   53 continue

 

      close(77)

     

     

      open(unit=78,file='c:\gmt-ground\check\check.dat')

      do 51 i=1,1600

      do 52 j=1,1600

      write(78,100)j, i, k(i,j)

  100 format(' ',i4,' ',i4,' ',f10.3)

   52 continue

   51 continue

      close(78)

      go to 999

***** プログラム作成エラー *********************************************

*                                                                      *

  996 write(*,*) 'You are bad!'

      go to 999

*                                                                      *

************************************************************************

  999 write(*,*)'確認終了'

      return

      end

************************************************************************

 

 

 

 

 

 

 

 

 

 

 

1.2 標高データ作成プログラム(4800×4800

使用方法(UNIX

作成したい領域に入る1600×1600のデータを9個作る。

9個のデータを左下から順番に4桁の地域番号を入力する。

************************************************************************

***** 標高データを配列に代入するプログラム(import) *********************

c      subroutine(a,e,n)

      real    a(4800,4800),k(1600,1600)

      integer i, j, n,x,y

      character mn*4

************************************************************************

      do n=1,9

      write(*,*) '地域番号を入れてください。'

      read(*,*)mn

      open(unit=77,file=mn//'.bin',access='direct',

     +                 form='unformatted',  recl=6400)

       do i=1,1600

       read(77,rec=i)(k(i,j),j=1,1600)  

       end do

      close(77)

************************************************************************

 

 

       if(n.ge.1.and.n.le.3)y=0

       if(n.ge.4.and.n.le.6)y=1

       if(n.ge.7.and.n.le.9)y=2

       if(n.eq.1 .or. n.eq.4 .or. n.eq.7)x=0

       if(n.eq.2 .or. n.eq.5 .or. n.eq.8)x=1

       if(n.eq.3 .or. n.eq.6 .or. n.eq.9)x=2

       do i=1,1600

        do j=1,1600

        a(i+1600*y,j+1600*x)=k(i,j)

        end do

       end do

      end do

     

      open(unit=7,file='data.terrain.bin',access='direct'

     +                  ,form='unformatted',recl=19200)

       do n=1,4800

         write(7,rec=n)(a(n,m),m=1,4800)

       end do

      close(unit=7)

     

     

c      write(*,*)a(3000,4000)

************************************************************************

c      close(78)

************************************************************************

      end

 

2.WGRIBを用いた外部気象データの作成

 注):太字は任意のファイル名および値を入力する。

  1.プログラムをコンパイルする。

  コマンド・・・cc wgrib.c > wgrib

2.GRIBファイルの内容物を確認する。

コマンド・・・wgrib GRIBファイル or wgrib –s GRIBファイル

  3.地表面以外のデータを取り出す

  コマンド・・・wgrib GRIBファイル | tail –n 91 | wgrib –i GRIBファイル –o 出力ファイル名

  4.各解析値をヘッダー無しのバイナリ形式で出力する。

コマンド・・・wgrib GRIBファイル | grep “:HGT:” | wgrib –i –bin –nh GRIBファイル –o出力ファイル名

  5.足りない解析値を用意する。

 6.CReSSの各解析値を読み込む順番にcatコマンド(UNIX)でつなげる。