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

      CALL GROPN(1)
      
      DO ITIME=1,NTIME ! 描画するデータ番号の指定
      
        CALL GTDATE(RTIME(ITIME),NY,NM,ND,RH)
        WRITE(CY,'(I4.4)') NY
        WRITE(CM,'(I2.2)') NM
        WRITE(CD,'(I2.2)') ND
        WRITE(CH,'(I2.2)') INT(RH)

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

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

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

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

        DO ILAT=1,NLAT
          DO ILON=0,NLON-1
            RTEMP(ILON,ILAT)=SCALE*ITEMP(ILON,ILAT)+OFFSET
          END DO
        END DO

! 経度の範囲を [0,360] から [-180,180]に詰め替え

        DO ILAT=1,NLAT
          DO ILON=0,NLON/2
            RTEMPG(ILON,ILAT)=RTEMP(ILON,ILAT)
          END DO
          DO ILON=1,NLON/2
            RTEMPG(-ILON,ILAT)=RTEMP(NLON-ILON,ILAT)
          END DO
        END DO

        DO ILON=0,NLON/2
          RLONG(ILON)=RLON(ILON)
        END DO        
        DO ILON=1,NLON/2
          RLONG(-ILON)=RLON(NLON-ILON)-360
        END DO

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

        CALL GRFRM
        CALL SGSVPT( 0.1, 0.9, 0.1, 0.9 )
        CALL SGSSIM(0.4, 0., 0. )
        CALL SGSMPL( 0., 90., 0. )
        CALL SGSTXY(-180.0,180.,0., 90. )
        CALL SGSTRN(30)
        CALL SGSTRF
        CALL SLPWWR( 1 )
        CALL SLPVPR( 1 )
        CALl UWSGXA(RLONG,NLON+1)
        CALl UWSGYA(RLAT,NLAT)
        CALL UEGTLA(180.0,310.0,10.0)
        CALL UETONE(RTEMPG,NLON+1,NLON+1,NLAT)
        CALL UDCNTR(RTEMPG,NLON+1,NLON+1,NLAT)      
        CALL UMPMAP( 'coast_world' )
        CALL UMPGLB
        CALL SGTXR(0.5,0.95,& 
&         'temp.'//CY//'/'//CM//'/'//CD//'/'//CH//'oclock')

      END DO

! ファイルを閉じる.

      CALL GRCLS
      IC=NF_CLOSE(NCID) 

      END
!----------------------------------------------------------------------- &
!     西暦1800年1月1日0時からの時間(RTIME)から
!      年(NY),月(NM),日(ND),時刻(RH)
!     を求めるサブルーチン      
!-----------------------------------------------------------------------
      SUBROUTINE GTDATE(RTIME,NY,NM,ND,RH)

      IMPLICIT NONE
      REAL :: RTIME,RH
      INTEGER :: IT,NY,NM,ND

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

      END
