C  '<?>'の部分を自分で埋めてプログラムを完成してください。

C  PARAMETER文で定数PIとREを宣言する。
C  定数PIは円周率を、REは地球半径を表す。
      PARAMETER (PI=3.1415927)
      PARAMETER (RE=6368.*1000.)

C  PARAMETER文で定数LMAXとMMAXを宣言する。
C  定数LMAXは経度方向の格子数を、MMAXは緯度方向の格子数を表す。
      PARAMETER (LMAX=144,MMAX=73)

C  配列U,V,XIを宣言する。
      REAL U(LMAX,MMAX),V(LMAX,MMAX),XI(LMAX,MMAX)

C  データを読みこむ。

C  東西風データのファイルを開く。
C  機番は10以降の番号を指定する。
C  STATUSは読みこみの場合は'OLD'を指定する。
C  バイナリファイルの場合は、FORMには'UNFORMATTED'を指定する。
C  さらに、ACCESSには'DIRECT'を指定し、
C  RECLにはレコード長をバイト数(整数値)で指定する。
C  各要素は4バイトなので、レコード長は4×要素数である。
      OPEN(10,FILE='U500_201104.dat',STATUS='OLD',
     +     FORM='UNFORMATTED',
     +     ACCESS='DIRECT',RECL=4*LMAX*MMAX)

C  READ文でデータを読みこむ。
C  機番10を指定する。REC=...でレコード番号を指定する。
C  データファイルには4月1日0時UTCから6時間ごとのデータが記録されていて、
C  ここでは、93番目のレコード(24日0時UTCのデータ)を読みこむ。
C  読みこんだデータを配列Uに格納する。
      READ(10,REC=93) U

C  ファイルを閉じる。
      CLOSE(10)

C  南北風データのファイルを開く。
C  機番は10以降の番号を指定する。
C  STATUSは読みこみの場合は'OLD'を指定する。
C  バイナリファイルの場合は、FORMには'UNFORMATTED'を指定する。
C  さらに、ACCESSには'DIRECT'を指定し、
C  RECLにはレコード長をバイト数(整数値)で指定する。
C  各要素は4バイトなので、レコード長は4×要素数である。
      OPEN(10,FILE='V500_201104.dat',STATUS='OLD',
     +     FORM='UNFORMATTED',
     +     ACCESS='DIRECT',RECL=4*LMAX*MMAX)

C  READ文でデータを読みこむ。
C  機番10を指定する。REC=...でレコード番号を指定する。
C  データファイルには4月1日0時UTCから6時間ごとのデータが記録されていて、
C  ここでは、93番目のレコード(24日0時UTCのデータ)を読みこむ。
C  読みこんだデータを配列Vに格納する。
      READ(10,REC=93) V

C  ファイルを閉じる。
      CLOSE(10)

C ==== ここから相対渦度の計算 ====

C  相対渦度を計算する。

C  東西、南北格子間隔の値を計算する。
C  東西間隔をDLA、南北格子間隔をDPHとする。
C  単位はラジアンである。
      DLA = 2.*PI/REAL(LMAX)
      DPH = PI/REAL(MMAX-1)

C  ここから相対渦度の計算のためのDOループが始まる。
C  北極と南極を除いて、LMAX*(MMAX-1)回だけ反復する。
      DO 11 M=2,MMAX-1

C  南北格子番号Mに対応する緯度を計算する。
C  その格子点の緯度をPH、ひとつ北の格子点の緯度をPHP、
C  ひとつ南の格子点の緯度をPHMとする。
C  M=1が北極、M=MMAXが南極である点に注意する。
        PH = 0.5*PI - PI*REAL(M-1)/REAL(MMAX-1)
        PHP = 0.5*PI - PI*REAL(M-1-1)/REAL(MMAX-1)
        PHM = 0.5*PI - PI*REAL(M+1-1)/REAL(MMAX-1)

      DO 12 L=1,LMAX

C  となりの格子における東西格子番号を計算する。
C  ひとつ東の格子における東西格子番号をLP、
C  ひとつ西の格子における東西格子番号をLMとする。
        LP = L+1
        IF (L.EQ.LMAX) LP = 1
        LM = L-1
        IF (L.EQ.1) LM = LMAX

C  ひとつ東(LP)と西(LM)の格子における南北風Vの値を求める。
C  ひとつ東の格子(LP,M)における南北風Vの値をVP、
C  ひとつ西の格子(LM,M)における南北風Vの値をVMとする。
        VP = V(LP,M)
        VM = V(LM,M)

C  ひとつ北(M-1)と南(M+1)の格子における東西風Uの値を求める。
C  ひとつ北の格子(L,M-1)における東西風Uの値をUP、
C  ひとつ南の格子(L,M+1)における東西風Uの値をUMとする。
        UP = U(L,M-1)
        UM = U(L,M+1)

C  相対渦度を計算する。
        XI(L,M) = <?>

C  ここでDOループ11と12が終了する。
 12   CONTINUE
 11   CONTINUE

C  北極における相対渦度を計算する。
        USUM = 0.
      DO 21 L=1,LMAX
        USUM = USUM + U(L,2)
 21   CONTINUE
        UAVE = USUM / REAL(LMAX)
      DO 22 L=1,LMAX
        XI(L,1) = 2. * UAVE / (RE * DPH)
 22   CONTINUE

C  南極における相対渦度を計算する。
        USUM = 0.
      DO 31 L=1,LMAX
        USUM = USUM + U(L,MMAX-1)
 31   CONTINUE
        UAVE = USUM / REAL(LMAX)
      DO 32 L=1,LMAX
        XI(L,MMAX) = - 2. * UAVE / (RE * DPH)
 32   CONTINUE

C ==== ここまで相対渦度の計算 ====

C  データを書き出す。

C  ファイルを開く。
C  機番は10以降の番号を指定する。
C  STATUSは書き出しの場合は'UNKNOWN'を指定する。
C  バイナリファイルの場合は、FORMには'UNFORMATTED'を指定する。
C  さらに、ACCESSには'DIRECT'を指定し、
C  RECLにはレコード長を整数値で指定する。
      OPEN(10,FILE='output.dat',STATUS='UNKNOWN',
     +     FORM='UNFORMATTED',
     +     ACCESS='DIRECT',RECL=4*LMAX*MMAX)

C  配列XIの各要素の値を機番10で指定されたファイルに書き出す。
C  REC=1と指定し、1番目のレコードとして出力する。
      WRITE(10,REC=1) XI

C  ファイルを閉じる。
      CLOSE(10)

      STOP
      END