      IMPLICIT NONE
      INCLUDE 'netcdf.inc' ! netcdf ファイルを操作するのに必要
      INTEGER,PARAMETER :: NLON=144,NLAT=73,NLEV=17,NTIME=366
      INTEGER(2) :: ITEMP(NLON,NLAT,NLEV)
      INTEGER :: ISTART(4),ICOUNT(4)
      INTEGER :: IC,NCID,IDAIR,IDLON,IDLAT,IDLEV,IDTIME
      INTEGER :: ILON,ILAT,ILEV,ITIME
      INTEGER :: NY,NM,ND
      REAL :: RLON(NLON),RLAT(NLAT),RLEV(NLEV),RTIME(NTIME)
      REAL :: RTEMP(NLON,NLAT,NLEV),RTEMPM(NLAT,NLEV)      
      REAL :: RH
      REAL :: OFFSET,SCALE      
      CHARACTER :: CY*4,CM*2,CD*2

! ファイルを開く
      
      IC=NF_OPEN('air.2008.nc',0,NCID)

! 各変数の ID を取得      
      
      IC=NF_INQ_VARID(NCID,'air',IDAIR)
      IC=NF_INQ_VARID(NCID,'lon',IDLON)
      IC=NF_INQ_VARID(NCID,'lat',IDLAT)
      IC=NF_INQ_VARID(NCID,'level',IDLEV)
      IC=NF_INQ_VARID(NCID,'time',IDTIME)

! RLON, RLAT, RLEV, RTIME の値を読み込む      

      IC=NF_GET_VAR_REAL(NCID,IDLON,RLON)
      IC=NF_GET_VAR_REAL(NCID,IDLAT,RLAT)
      IC=NF_GET_VAR_REAL(NCID,IDLEV,RLEV)
      IC=NF_GET_VAR_REAL(NCID,IDTIME,RTIME)

! OFFSET と SCALE の値を読み込む

      IC=NF_GET_ATT_REAL(NCID,IDAIR,'add_offset',OFFSET)
      IC=NF_GET_ATT_REAL(NCID,IDAIR,'scale_factor',SCALE)

! dcl の描画の最初の部分はここでやっておく

      CALL SWLSET('LWAIT',.FALSE.) ! 自動的に次のページに
      CALL SWLSET('LALT',.TRUE.)   ! 裏画面で描画して変化を滑らかに
      CALL GROPN(1)

! 毎日のデータをアニメーション表示するためのループ

      DO ITIME=1,366  ! ITIME は 年頭1月1日からの1日毎のデータの番号
      
        CALL GTDATE(RTIME(ITIME),NY,NM,ND,RH)
        WRITE(CY,'(I4.4)') NY
        WRITE(CM,'(I2.2)') NM
        WRITE(CD,'(I2.2)') ND

! 指定された時刻のデータを配列 ITEMPに読み込む
! ISTART はデータの先頭位置, ICOUNT は各次元方向に読み込む数

        ISTART(1)=1
        ISTART(2)=1
        ISTART(3)=1
        ISTART(4)=ITIME      
      
        ICOUNT(1)=NLON
        ICOUNT(2)=NLAT
        ICOUNT(3)=NLEV
        ICOUNT(4)=1

        IC=NF_GET_VARA_INT2(NCID,IDAIR,ISTART,ICOUNT,ITEMP)

! 温度データをちゃんとした値に変換

        DO ILEV=1,NLEV
          DO ILAT=1,NLAT
            DO ILON=1,NLON
              RTEMP(ILON,ILAT,ILEV)= &
&               SCALE*ITEMP(ILON,ILAT,ILEV)+OFFSET
            END DO
          END DO
        END DO

! 経度平均の計算

        DO ILEV=1,NLEV
          DO ILAT=1,NLAT
            RTEMPM(ILAT,ILEV)=0
            DO ILON=1,NLON
              RTEMPM(ILAT,ILEV)=RTEMPM(ILAT,ILEV)+ &
&               RTEMP(ILON,ILAT,ILEV)
            END DO
            RTEMPM(ILAT,ILEV)=RTEMPM(ILAT,ILEV)/NLON          
          END DO
        END DO

! ここから先は dcl による描画

        CALL GRFRM
        CALL GRSWND(-90.0,90.0,1000.0,10.0)
        CALL GRSVPT(0.2,0.9,0.2,0.9)
        CALL GRSTRN(2)
        CALL GRSTRF
        CALl UWSGYA(RLEV,NLEV)
        CALl UWSGXA(RLAT,NLAT)
        CALL UEGTLA(180.0,310.0,10.0)
        CALL UETONE(RTEMPM,NLAT,NLAT,NLEV)
        CALL UXAXDV('B',10.0,30.0)
        CALL UXAXDV('T',10.0,30.0)
        CALL ULISET('IYTYPE',3)
        CALL ULYLOG('L',3,9)
        CALL ULYLOG('R',3,9)
        CALL UXSTTL('T','zonal mean temp. '//CY//'/'//CM//'/'//CD,0.0)
        CALL UXSTTL('B','Latitude(deg.)',0.0)
        CALL UYSTTL('L','Pressure(hPa)',0.0)
        CALL UDCNTR(RTEMPM,NLAT,NLAT,NLEV)

      END DO ! 毎日のループの終了

! dcl の終了

      CALL GRCLS

! ファイルを閉じる.
      
      IC=NF_CLOSE(NCID)       
      
      END
!----------------------------------------------------------------------- &
!     西暦1年1月1日0時からの時間(RTIME)から年(NY),月(NM),日(ND),時刻(RH)
!     を求めるサブルーチン      
!-----------------------------------------------------------------------
      SUBROUTINE GTDATE(RTIME,NY,NM,ND,RH)

      IMPLICIT NONE
      REAL,PARAMETER :: RT1948=17067072 ! 1948年1月1日0時の時間
      REAL :: RTIME,RH
      INTEGER :: IT,NY,NM,ND

      IT=(RTIME-RT1948)/24
      RH=(RTIME-RT1948)-IT*24
      
      CALL DATEF3(IT,1948,1,1,NY,NM,ND) ! DATEF3は DCLのサブルーチン

      END
