C  配列の大きさ: IHXは1日の観測回数(=4)、IDXは1月の日数(=31)、
C                ILXは気圧面の数(=26)、NXは地点数(=1)である。
      PARAMETER (IHX=4,IDX=31,ILX=26,NX=1)

C  配列(入力データ)の宣言:
C    ISTA (NX): 観測地点番号。館野(つくば)は47464など。
C    ICHECK (IHX,IDX,NX): データが存在する場合は1、
C                         存在しない場合は0。
C    IZ (IHX,IDX,ILX,,NX): ジオポテンシャル高度 [m]。
C    IU (IHX,IDX,ILX,,NX): 相対湿度 [%]。
C    ID (IHX,IDX,ILX,,NX): 風向 [°、北=0、東=90]。
C    P (IHX,IDX,ILX,,NX): 気圧 [hPa]。
C    T (IHX,IDX,ILX,,NX): 気温 [℃]。
C    S (IHX,IDX,ILX,,NX): 風速 [m/s]。
      INTEGER ISTA(NX),ICHECK(IHX,IDX,NX),
     +        IZ(IHX,IDX,ILX,NX),IU(IHX,IDX,ILX,NX),
     +        ID(IHX,IDX,ILX,NX)
      REAL P(IHX,IDX,ILX,NX),T(IHX,IDX,ILX,NX),
     +     S(IHX,IDX,ILX,NX)

C  配列(出力データ)の宣言:
C    TD (IHX,IDX,ILX,NX): 露点温度 [℃]。
      REAL TD(IHX,IDX,ILX,NX)

C  欠損値: データが無効のとき、IMISS、RMISSで指定された値が格納される。
C          整数型の場合はIMISS、実数型の場合はRMISSである。
      PARAMETER (IMISS=999999,RMISS=1.0E9)

C  地点番号を指定する。
C  入力を求めるメッセージを標準出力に書き出す。
      WRITE(6,*) 'Station No. ?'

C  READ文で入力値を読みこむ。
      READ (5,*) ISTA(1)

C  日時を指定する。
C  入力を求めるメッセージを標準出力に書き出す。
      WRITE(6,*) 'Year, Month, Day, Hour(3,9,15,21) ?'

C  READ文で入力値を読みこむ。
      READ (5,*) IYY,IMM,IDD,IHH
      IHH = (IHH + 3) / 6

C  サブルーチンAERDTによって、
C  指定された年、月の観測データを読むこむ。
C  ※サブルーチンAERDTは、IYY、IMMで指定された年、月について
C    データを読みこむ。
      CALL AERDT
     +  (IHX,IDX,ILX,NX,ISTA,IYY,IMM,IMISS,RMISS,
     +   ICHECK,P,IZ,T,IU,ID,S)

C  ここからDOループ11〜14が始まる。
C  各時刻の各気圧面について、露点温度を計算する。
      DO 11 N=1,NX
      DO 12 IL=1,ILX
      DO 13 IDDD=1,IDX
      DO 14 IHHH=1,IHX

C  出力データを格納する配列に、あらかじめ欠損値を入れておく。
        TD(IHHH,IDDD,IL,N) = RMISS

C  データが存在する場合だけ、
C  つまり、T(IHH,IDD,IL,N)、IU(IHH,IDD,IL,N)が欠損値でない場合だけ、
C  計算を実行する。
      IF ((T(IHHH,IDDD,IL,N).NE.RMISS)
     +    .AND.(IU(IHHH,IDDD,IL,N).NE.IMISS)) THEN

C  関数TDEWを用いて、露点温度を計算する。
C  結果をTDに代入する。
        TD(IHHH,IDDD,IL,N) = TDEW(T(IHHH,IDDD,IL,N)+273.15,
     +                            0.01*REAL(IU(IHHH,IDDD,IL,N)))
     +                       -273.15

      ENDIF

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

C  出力ファイルを開く。
C  機番は10以降の番号を指定する。
C  STATUSは書き出しの場合は'UNKNOWN'を指定する。
C  FORMはテキストファイルの場合は'FORMATTED'を指定する。
      OPEN(10,FILE='output.txt',STATUS='UNKNOWN',
     +     FORM='FORMATTED')

C  ここからDOループ21が始まる。
C  指定気圧面について、結果を出力する。
C  地上(IL=1)の結果は出力しない。
      DO 21 IL=2,ILX

C  データが存在する場合だけ、つまり、P(IHH,IDD,IL,1), T(IHH,IDD,IL,1),
C   TD(IHH,IDD,IL,1)が欠損値でない場合だけ、
C  結果を出力する。
      IF ((P(IHH,IDD,IL,1).NE.RMISS)
     +    .AND.(T(IHH,IDD,IL,1).NE.RMISS)
     +    .AND.(TD(IHH,IDD,IL,1).NE.RMISS)) THEN

C  気圧、気温、露点温度を出力する。
C  関数INTで気圧の値を整数型に変換する。INTは切り捨て、NINTは四捨五入。
        IP = INT(P(IHH,IDD,IL,1))
        WRITE(10,'(1X,I4,1X,F5.1,1X,F5.1)')
     +    IP,T(IHH,IDD,IL,1),TD(IHH,IDD,IL,1)

      ENDIF

C  ここでDOループ21が終了する。
 21   CONTINUE

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

      STOP
      END

C  主プログラムの後に、関数やサブルーチンのような副プログラムを書く。

C  露点温度を計算するための関数。
C  関数の種類は実数(REAL)、関数名はTDEWとする。
C  引数はTとHである。
      REAL FUNCTION TDEW (T,H)

C  露点温度を計算する。
C  ALOGは対数関数である。
        TDEW = 35.86
     +       + (T-35.86) / (1. - (T-35.86) / 237.3 * ALOG(H) / 17.27)

C  関数はRETURN文とEND文で終了する。
C  RETURN文が実行された時点でのTDEWの値が関数の値として返される。
      RETURN
      END