EXAMPLE (GEOの編集)

計算領域	南緯30度-北緯60度，東経120度-西経80度
解像度	1 $^\circ$ $\times$ 1 $^\circ$ ，鉛直33層(下表)
海底地形	現実的な地形

level	1	2	3	4	5	6	7	8	9	10
depth(m)	0	10	20	30	50	75	100	125	150	200
level	11	12	13	14	15	16	17	18	19	20
depth(m)	250	300	400	500	600	700	800	900	1000	1100
level	21	22	23	24	25	26	27	28	29	30
depth(m)	1200	1300	1400	1500	1750	2000	2500	3000	3500	4000
level	31	32	33
depth(m)	4500	5000	5500

地形情報に現実的な海底地形を書き込むために ETOPO5 global depth dataを使用します(//lanina/E/d/datagrid/topographyに topo.grdという名前で置いてあります．コントロールファイルtopo.ctlも置いてあるので，topo.ctlを開き地形データtopo.grdの解像度やグリッド数を確認してください．またgradsで地形の様子を見るとよいでしょう)．
この地形データの解像度は $\frac{1}{12}^\circ \times \frac{1}{12}^\circ$ であり，高さ(深さ)はmで表現されています．そこでまずは解像度を1 $^\circ$ $\times$ 1 $^\circ$ に変更します．そしてGEOファイルの形式に合わせて，深さをmから層番号に変更します．
以下の作業はMATLABを使って行っています．
GEOを編集するためのMATLABプログラムは//lanina/E/d/COCO3.3/programに置いてあります．以下はプログラムMktopo1_1.mの説明です．

まずはtopo.grdを読み込みます．

fid=fopen('topo.grd','rb');

rdgrd1,count

=fread(fid,[4320 2161],'float32');

次に計算領域部分のデータを取り出します．

    i1=120;, i2=280;, j1=60;, j2=150;
    rdgrd1=rdgrd1(i1*12+1:1:i2*12,j1*12+1:1:j2*12);

領域平均を行い，解像度を1 $^\circ$ $\times$ 1 $^\circ$ に変更します．

    for j=0:(j2-j1)-1
        for i=0:(i2-i1)-1
            a=rdgrd1(12*i+1:1:12*i+12,12*j+1:1:12*j+12);
            a=reshape(a,144,1);
            newgrd(i+1,j+1)=sum(a)/144 ;
        end
    end

ここでnewgrdが南緯30度-北緯60度，東経120度-西経80度，解像度1 $^\circ$ $\times$ 1 $^\circ$ に変更された地形データとなります．
次に深さの単位をmから層番号に変更します．陸はGEOファイルで0と定義されるので，0 mより大きい値を0に置き換えます．

    newgrd(newgrd > 0)=0;

ここでnewgrdを便宜上Aと置き換えます．

    A=newgrd;

深さの単位をmから層番号に変更します(上の表に示した33層で離散化する場合)．

        i1=1; i2=160; j1=1; j2=90;
        for i=i1:i2
            for j=j1:j2
                if (A(i,j) < -5500);
                    A(i,j) = 33;
                end
                if (-5500 <= A(i,j)) & (A(i,j) < -5000);
                    A(i,j) = 32;
                end
                if (-5000 <= A(i,j)) & (A(i,j) < -4500);
                    A(i,j) = 31;
                end
                if (-4500 <= A(i,j)) & (A(i,j) < -4000);
                    A(i,j) = 30;
                end
                if (-4000 <= A(i,j)) & (A(i,j) < -3500);
                    A(i,j) = 29;
                end
                if (-3500 <= A(i,j)) & (A(i,j) < -3000);
                    A(i,j) = 28;
                end
                         $\cdots$ 省略 $\cdots$
                if (-75 <= A(i,j)) & (A(i,j) < -50);
                    A(i,j) = 5;
                end
                if (-50 <= A(i,j)) & (A(i,j) < -30);
                    A(i,j) = 4;
                end
                if (-30 <= A(i,j)) & (A(i,j) < -20);
                    A(i,j) = 3;
                end
                if (-20 <= A(i,j)) & (A(i,j) < -10);
                    A(i,j) = 2;
                end
                if (-10 <= A(i,j)) & (A(i,j) < -0);
                    A(i,j) = 1;
                end
            end
        end

GEOファイルに地形情報を書き込むときは，一番上の行が北端，下の行が南端となるように書きこむので，Aの配列を変更する必要がある．
転置して

    A=A';

行の順番をひっくり返す

    for n=1:45;
        A([91-n,n],:)=A([n,91-n],:);
    end

モデルにおいては，領域の東西端は常に周期的につながっているものとして扱われている．そこでモデルの東端を陸にする．

    A(:,160)=0;

そしてファイルGEOにテキスト形式で地形情報を書き込む．

    save GEO A -ascii

これでGEOファイルの地形情報の部分は完成である．
地形情報以外の部分はGEOファイルをemacs, viなどで開き，直接編集するとよい．
編集するときは以下の点に注意すること．

(1) KZ=1とすれば，単純なz座標モデルとなる．
(2) 東西方向の領域分割数，南北方向の領域分割数はそれぞれNX，NYの約数でなければならない．
(3) 座標系回転を行わない場合は，座標系回転の角度を0とする．

以上でGEOファイルの編集が終了した．

次にdimocn.Fを開き，NX，NY，NZ，KZの値を適切に書き込む (NX=160，NY=90，NZ=33，KZ=1).

	-- dimocn.F --
c ***	Write here the values for NX, NY, NZ, and KZ as you specified
c ***	in the geometry file
	INTEGER NX, NY, NZ
	PARAMETER(NX = 160, NY = 90, NZ = 33)	$\leftarrow$ ここ
	INTEGER KZ
	PARAMETER(KZ = 1)	$\leftarrow$ ここ
c ***	You need not change the followings
	INTEGER NXDIM, NYDIM, NZDIM
	PARAMETER(NXDIM = NX+4, NYDIM = NY+4, NZDIM = NZ+2)