*" GT2NCDF:  GTOOLファイルを netCDF ファイルに変換する。
*" [履歴]
*"  ver.0 堀之内武 1997/08/15  製作
*"        堀之内武 1997/09/05  バグフィックス
*"  ver.1 堀之内武 1997/09/12
*
*" コマンドライン引数
*"   入力ファイル名     (default: $GTTMPDIR/gtool.out)
*"   var:主変数名として使うgtoolヘッダー項目名 (default: item)
*"   out:出力ファイル名 (default: $GTTMPDIR/gtool.out)
*"   help
*"  
*" ver.0 コメント:
*"  netCDFライブラリーとインクルードファイルが必要; 通常は
*"  -L/usr/local/lib -lnetcdf -I/usr/local/include であろう。
*"  (※) 単位が udunits に合致するように変換されたかどうか確認するため
*"  に -ludunits もリンクしたいのだが(subroutine GT2NC_UNIT)、今のとこ
*"  ろうまく呼べていない(正しいはずの単位を入れてもエラーが返る)ので、
*"  使っていない。
*"  
*" ver.1 コメント(ver.0との違い. 詳しくはREADMEを見よ):
*"  変数が複数あるファイルへの対応。但し、各変数の時間進行は同時であるこ
*"  とを仮定: a(t=1),b(t=1),c(t=1),a(t=2),b(t=2),c(t=2),..という具合い.
*"  今のところ a,b,c,.. の精度は同じでなければならない
*"  軸重み変数を追加.
*"  属性の追加とデバグ.
*"  
*********************************************************************
      PROGRAM GT2NCDF
*
      include   'netcdf.inc'
*
      INTEGER    IDIM,JDIM,KDIM,IJKDIM
*
#ifdef SYS_IBMS
      INCLUDE    (GTSINC)
      INCLUDE    (GZSIZE)
#else
#include         "gtsinc.F"
#include         "gzsize.F"
#endif
*
*     [NAMELIST]
      CHARACTER  HFILE( 1 ) *(NFILN)
      DATA       HFILE / '$GTTMPDIR/gtool.out' /
      CHARACTER  VAR*(10)  !" The item in the header used as var_name
      DATA       VAR   / 'item' /
      CHARACTER  OUT    *(NFILN)
      DATA       OUT    / '$GTTMPDIR/gtool.nc' /
      LOGICAL    HELP
      DATA       HELP   / .FALSE. /
*     
      NAMELIST  /OPTION/ HFILE, VAR, OUT, HELP
     &                   
*
*     [OTHERS]
      CHARACTER  HHEAD ( NDC )*(NCC)
      REAL       GDATA ( IJKDIM )
      CHARACTER  HHEAD2 ( NDC )*(NCC)
      REAL       GDATA2 ( IJKDIM )
*
      INTEGER    IFILE, JFILE
      DATA       IFILE / 50 /
      DATA       JFILE / 60 /
*
      LOGICAL    TIMESEQ          !" time sequence or not
      LOGICAL    LSECOND
      INTEGER    NOPT, NFILE, IOS, IL, IEOD
      INTEGER    LENC,LENZ
      INTEGER    NCERR, NCID
      INTEGER    ITIME                     !" 時刻(通し)
      INTEGER    STAT,IAX
*
      INTEGER    START(4),COUNT(4)
      INTEGER    NTIME
      INTEGER    AXNUM                       !" 読み込んだデータの次元数
      INTEGER    NCVAR_TYPE       !" netCDF var type. NF_[FLOAT|DOUBLE]
      INTEGER    VARID_DATA       !" varid of the data
      INTEGER    VARID_DATA1      !" varid of the first data
      INTEGER    VARID_TIME       !" varid of time if TIMESEQ
      REAL*8     MISS_VAL
*
      LOGICAL    REINIT
*
      CALL OPTARG ( 91, 'OPTION', 'HFILE', NOPT, NFILE )
      READ (91,OPTION,IOSTAT=IOS)
      CLOSE(91)
      IF ( HELP ) THEN
         WRITE(6,*) 
         WRITE(6,*) 'Usage: gt2ncdf [in_file] [out:o_file]'//
     &              ' [var:gtool_header_item] [-help]'
         WRITE(6,*) 'Example % gt2ncdf u.gtdata out:u.nc'
         WRITE(6,*) 
     &     '        % gt2ncdf restart.gtdata out:restart.nc var:titl'
         WRITE(6,*) 
         WRITE(6,*) 'Current setting:'
         WRITE(6,OPTION)
         STOP
      ELSE IF ( IOS.NE.0 ) THEN
         WRITE(6,OPTION)
         STOP
      ENDIF
*
      CALL ENT_HEAD2VAR(VAR)
*
*     / open the input file (gtool format) /
*
      CALL GTOPEN
      CALL GTSIZE ( HHEAD, IJKDIM )
      CALL GTSIZE ( HHEAD2, IJKDIM )
*
      CALL GURNTF ( HFILE( 1 ), OUT  , '$GTTMPDIR/gtool.in' )
*
      CALL GFROPN ( IFILE , HFILE( 1 ) )
*
*     / set path for output file (netcdf format) /
*
      CALL GUNENV( OUT,'.',.FALSE. ) !" set path ($GTTMPDIR, or .)
      IL=LENC(OUT)
      WRITE (6,*) 'output='//OUT(1:IL)
*
      NCERR=NF_CREATE(OUT,NF_CLOBBER, NCID)
      IF(NCERR.NE.NF_NOERR) CALL MSGDMP('E','gt2ncd',NF_STRERROR(NCERR))
*
*     / first step : HEADER -> ncdf file initialization /
*
      CALL   GFREAD
     O        ( HHEAD , GDATA , IEOD  ,
     I          IFILE , 1               )
*
      CALL   GFREAD
     O        ( HHEAD2 , GDATA2 , IEOD  ,
     I          IFILE , 1               )
*
      IF ( IEOD .EQ.0 ) THEN
        TIMESEQ=.TRUE.  !" since the data has not ended
      ELSE
        TIMESEQ=.FALSE.
      ENDIF
*
*     / ncdf file initialization /    
*
      CALL GT2NC_INIT
     I     (NCID, HHEAD, TIMESEQ,
     O      AXNUM,NCVAR_TYPE,VARID_DATA,VARID_TIME,COUNT,MISS_VAL)
      VARID_DATA1=VARID_DATA
*
*     / output the first data /    
*
      IF(TIMESEQ) THEN
        CALL GHPGET( HHEAD, 'TIME',  ITIME  )
        STAT=NF_PUT_VAR1_REAL(NCID,VARID_TIME,1,REAL(ITIME))
      ENDIF
*
      DO IAX=1,AXNUM+1
        START(IAX)=1
        NTIME=1
      ENDDO
*
      IF (NCVAR_TYPE .EQ. NF_FLOAT) THEN
        CALL MISS_CHG_R(HHEAD,GDATA,MISS_VAL)
        STAT=NF_PUT_VARA_REAL(NCID,VARID_DATA,START,COUNT,GDATA)
      ELSE
        CALL MISS_CHG_D(HHEAD,GDATA,MISS_VAL)
        STAT=NF_PUT_VARA_DOUBLE(NCID,VARID_DATA,START,COUNT,GDATA)
      ENDIF
*
*     / succeeding data /    
*
      IF (TIMESEQ) THEN
        LSECOND=.TRUE.
 1100   CONTINUE
          IF(LSECOND) THEN  !" already read
            CALL GHCOPY(HHEAD2,HHEAD)
            CALL DTCOPY(IJKDIM,GDATA2,GDATA)
            LSECOND=.FALSE.
          ELSE
            CALL   GFREAD
     O        ( HHEAD , GDATA , IEOD  ,
     I          IFILE , 1               )
          ENDIF
          IF ( IEOD .EQ.0 ) THEN
*
            CALL GT2NC_REINIT(NCID,HHEAD,
     O           AXNUM,REINIT,VARID_DATA,COUNT)
            IF (REINIT) THEN !" initialztn for a new variable
              CALL GT2NC_INIT
     I           (NCID, HHEAD, TIMESEQ,
     O           AXNUM,NCVAR_TYPE,VARID_DATA,VARID_TIME,COUNT,MISS_VAL)
            ENDIF
*
            IF(VARID_DATA .EQ. VARID_DATA1) THEN
              NTIME=NTIME+1
*
              CALL GHPGET( HHEAD, 'TIME',  ITIME  )
              STAT=NF_PUT_VAR1_REAL(NCID,VARID_TIME,NTIME
     &                            ,REAL(ITIME))
            ENDIF
*
            DO IAX=1,AXNUM
              START(IAX)=1
            ENDDO
            START(AXNUM+1)=NTIME
*
            IF (NCVAR_TYPE .EQ. NF_FLOAT) THEN
              CALL MISS_CHG_R(HHEAD,GDATA,MISS_VAL)
              STAT=NF_PUT_VARA_REAL(NCID,VARID_DATA,START,COUNT,GDATA)
            ELSE
              CALL MISS_CHG_D(HHEAD,GDATA,MISS_VAL)
              STAT=NF_PUT_VARA_DOUBLE(NCID,VARID_DATA,START,COUNT,GDATA)
            ENDIF
            IF(STAT.NE.NF_NOERR)
     &           CALL MSGDMP('E','GT2NCD',NF_STRERROR(STAT))
*
        GOTO 1100
          ENDIF
      ENDIF
*
      STAT=NF_CLOSE(NCID)
      IF(STAT.NE.NF_NOERR) CALL MSGDMP('E','GT2NCD',NF_STRERROR(STAT))
*
**      CALL GTCLSE
      STOP
      END
*********************************************************************
      SUBROUTINE MISS_CHG_R(HHEAD,GDATA,MISS_VAL)
      CHARACTER  HHEAD  ( * )*(*)        !" ヘッダー(入力)
      REAL       GDATA  ( * )            !" データ(入出力)
      REAL*8     MISS_VAL                !" 欠損値(入力)
      INTEGER    IXSTR,IXEND,IXDIM, IYSTR,IYEND,IYDIM, IZSTR,IZEND,IZDIM 
      REAL       VMISS
      INTEGER    I
*
      CALL GUQSIZ
     I         ( HHEAD ,
     O           IXSTR , IXEND , IXDIM ,
     O           IYSTR , IYEND , IYDIM ,
     O           IZSTR , IZEND , IZDIM  )
      CALL GHPGET( HHEAD, 'MISS', VMISS )
*
      DO I=1,IXDIM*IYDIM*IZDIM
        IF( GDATA(I).EQ.VMISS ) THEN
          GDATA(I)=MISS_VAL
        ENDIF
      ENDDO
      END
*********************************************************************
      SUBROUTINE MISS_CHG_D(HHEAD,GDATA,MISS_VAL)
      CHARACTER  HHEAD  ( * )*(*)        !" ヘッダー(入力)
      REAL*8     GDATA  ( * )            !" データ(入出力)
      REAL*8     MISS_VAL                !" 欠損値(入力)
      INTEGER    IXSTR,IXEND,IXDIM, IYSTR,IYEND,IYDIM, IZSTR,IZEND,IZDIM 
      REAL       VMISS
      INTEGER    I
      LOGICAL    LREQ1
      EXTERNAL   LREQ1
*
      CALL GUQSIZ
     I         ( HHEAD ,
     O           IXSTR , IXEND , IXDIM ,
     O           IYSTR , IYEND , IYDIM ,
     O           IZSTR , IZEND , IZDIM  )
      CALL GHPGET( HHEAD, 'MISS', VMISS )
*
      DO I=1,IXDIM*IYDIM*IZDIM
        IF( LREQ1(REAL(GDATA(I)),VMISS) ) GDATA(I)=MISS_VAL
      ENDDO
      END
*********************************************************************
      SUBROUTINE DTCOPY(NUM,IN,OUT)
      INTEGER NUM,I
      REAL IN(NUM),OUT(NUM)
      DO I=1,NUM
        OUT(I)=IN(I)
      ENDDO
      END
*********************************************************************
      LOGICAL FUNCTION NH_INQ_VAREXIST(NCID,VAR)

      include   'netcdf.inc'

      CHARACTER  VAR*(*)
      CHARACTER  VN*(100)
      INTEGER    NCID,I,STAT,NVARS

      NH_INQ_VAREXIST=.FALSE.

      STAT=NF_INQ_NVARS(NCID,NVARS)
      DO I=1,NVARS
        STAT=NF_INQ_VARNAME(NCID,I,VN)
        IF(VAR.EQ.VN) THEN
          NH_INQ_VAREXIST=.TRUE.
          RETURN
        ENDIF
      ENDDO

      RETURN
      END
*********************************************************************
      LOGICAL FUNCTION NH_INQ_DIMEXIST(NCID,DIM)

      include   'netcdf.inc'

      CHARACTER  DIM*(*)
      CHARACTER  DM*(100)
      INTEGER    NCID,I,STAT,NDIMS

      NH_INQ_DIMEXIST=.FALSE.

      STAT=NF_INQ_NDIMS(NCID,NDIMS)
      DO I=1,NDIMS
        STAT=NF_INQ_DIMNAME(NCID,I,DM)
        IF(DIM.EQ.DM) THEN
          NH_INQ_DIMEXIST=.TRUE.
          RETURN
        ENDIF
      ENDDO

      RETURN
      END
*********************************************************************
      SUBROUTINE HEAD2VAR(HHEAD,VAR)
      CHARACTER  HHEAD(*)*(*) !" the first header of the input file
      CHARACTER  VAR*(*)
      CHARACTER  HIT*(10),HITE*(*)
      SAVE       HIT
      INTEGER    I,J,LN
*
      IF(HIT.NE.'TITL') THEN
        CALL GHCGET( HHEAD, HIT, VAR  )
      ELSE
        CALL GHCGTS( HHEAD, HIT, VAR  )
      ENDIF
* --- / keep only the first ward (FOR RESTART FILE) /
      LN=LEN(VAR)
      J=-999
      DO I=1,LN
        IF(VAR(I:I).EQ.' ') THEN
          J=I
          GOTO 99
        ENDIF
      ENDDO
   99 CONTINUE
      IF(J.GT.1) THEN
        VAR(J:LN)=''
      ENDIF
      RETURN
*----------------------------
      ENTRY ENT_HEAD2VAR(HITE)
      HIT=HITE
      CALL CUPPER(HIT) !" capitarize
      RETURN
      END
*********************************************************************
*     / Do we need initialization ? /
      SUBROUTINE GT2NC_REINIT
     I     (NCID, HHEAD,
     O      AXNUM,REINIT, VARID, COUNT)
*
      include   'netcdf.inc'
*
#ifdef SYS_IBMS
      INCLUDE    (GZSIZE)
#else
#include         "gzsize.F"
#endif
*
*     [INPUT]
      CHARACTER  HHEAD(NDC)*(NCC) !" the first header of the input file
      INTEGER    NCID             !" netcdf ID of the output file
*     [OUTPUT]
      INTEGER    AXNUM                       !" 読み込んだデータの次元数
      LOGICAL    REINIT
      INTEGER    VARID  !" defiend if .NOT.GT2NC_REINIT (i.e. already defined)
      INTEGER    COUNT(4)!" defiend if .NOT.GT2NC_REINIT (i.e. already defined)
*
      LOGICAL    NH_INQ_VAREXIST
      EXTERNAL   NH_INQ_VAREXIST
      CHARACTER  VARNAME*(2*NCC)             !" 識別名称
      INTEGER    STAT,DIMIDS(4),NDIMS,I
*
      CALL HEAD2VAR(HHEAD,VARNAME)
      REINIT = .NOT. NH_INQ_VAREXIST(NCID,VARNAME)
      IF(.NOT.REINIT) THEN
        STAT=NF_INQ_VARID(NCID,VARNAME,VARID)
        STAT=NF_INQ_VARNDIMS(NCID,VARID,NDIMS)
        STAT=NF_INQ_VARDIMID(NCID,VARID,DIMIDS)
        DO I=1,NDIMS-1
          STAT=NF_INQ_DIMLEN(NCID,DIMIDS(I),COUNT(I))
        ENDDO
        COUNT(NDIMS)=1  !" this must be the time dimmension
        AXNUM=NDIMS-1
      ENDIF
*
      END
*********************************************************************
      SUBROUTINE GT2NC_INIT
     I     (NCID, HHEAD, TIMESEQ,
     O      AXNUM,NCVAR_TYPE,VARID_DATA,VARID_TIME,COUNT,MISS_VAL)
*
      include   'netcdf.inc'
*
      INTEGER    NW,NAXIS !" defined in GZIWRK
      INTEGER    IDIM,JDIM,KDIM,IJKDIM
*
#ifdef SYS_IBMS
      INCLUDE    (GTSINC)
      INCLUDE    (GZSIZE)
      INCLUDE    (GZIWRK)
#else
#include         "gtsinc.F"
#include         "gzsize.F"
#include         "gziwrk.F"
#endif
*
*     [INPUT]
      CHARACTER  HHEAD(NDC)*(NCC) !" the first header of the input file
      INTEGER    NCID             !" netcdf ID of the output file
      LOGICAL    TIMESEQ          !" time sequence or not
*     [OUTPUT]
      INTEGER    AXNUM                       !" 読み込んだデータの次元数
      INTEGER    NCVAR_TYPE       !" netCDF var type. NF_[FLOAT|DOUBLE]
      INTEGER    VARID_DATA       !" varid of the data
      INTEGER    VARID_TIME       !" varid of time if TIMESEQ
      INTEGER    COUNT(4)
      REAL*8     MISS_VAL
*
*     [CONTROLING LOGICALS]
      LOGICAL    FIRST  ! Is this first call?
      DATA       FIRST/.TRUE./
      SAVE       FIRST
      LOGICAL    NEWDIM(3)
*
*     [GTOOL AXIS]
      REAL       AX(NW,3)
      REAL       AXW(NW,3)
      CHARACTER  HHEADA(NDC)*(NCC)
      CHARACTER  HDSETA(3)*(NCC)             !" データセット名
      CHARACTER  HTITLA(3)*(2*NCC)           !" タイトル
      CHARACTER  HUNITA(3)*(NCC)             !" 単位
      INTEGER    STYPA(3)                    !" スケーリングタイプ
      CHARACTER  HUNITAW(3)*(NCC)            !" 単位(重み)
*
*     [GTOOL HEADER ITEMS]
      CHARACTER  HDSET  *(NCC)             !" データセット名
      CHARACTER  HITEM  *(NCC)             !" 識別名称
      INTEGER    IFNUM                     !" ファイル番号
      INTEGER    IDNUM                     !" データ番号
      CHARACTER  HTITL  *(2*NCC)           !" タイトル
      CHARACTER  HUNIT  *(NCC)             !" 単位
      INTEGER    ITIME                     !" 時刻(通し)
      CHARACTER  HDATE  *(NCC)             !" 時刻(yyyymmddhhmmss)
      CHARACTER  HUTIM  *(NCC)             !" 時刻単位
      INTEGER    ITDUR                     !" 代表する時間
      CHARACTER  HAITM(3)*(NCC)            !" 軸1-3の格子名称
      INTEGER    IASTR(3)                  !" 軸1-3の格子始め
      INTEGER    IAEND(3)                  !" 軸1-3の格子終り
      INTEGER    IANUM(3)                  !" 軸1-3の格子数
      CHARACTER  HDFMT  *(NCC)             !" データフォーマット
      REAL       VMISS                     !" 欠損値の値
      REAL       DMIN                      !" レンジ(最小)
      REAL       DMAX                      !" レンジ(最大)
      REAL       DIVS                      !" 間隔(小)
      REAL       DIVL                      !" 間隔(大)
      INTEGER    ISTYP                     !" スケーリングタイプ
      CHARACTER  HCDATE *(NCC)             !" データ作成日付
      CHARACTER  HCSIGN *(NCC)             !" データ作成者
      CHARACTER  HMDATE *(NCC)             !" データ変更日付
      CHARACTER  HMSIGN *(NCC)             !" データ変更者
*
      CHARACTER  HETTL*(8*NCC)
*
*     [NETCDF]
      CHARACTER  VARNAME*(2*NCC)  !" 主変数名 <- HITEM
      INTEGER    DIMID(4)
      INTEGER    VARID_AX(3)
      INTEGER    VARID_AXW(3)
      CHARACTER  DIMNAME(3)*(16)   !" 次元/軸変数の名前
      CHARACTER  DIMNAMEW(3)*(19)     !" 軸重みの名前
      CHARACTER  TITLEW*(26)          !" 軸重みのタイトル
      CHARACTER  CONVENTIONS*(NCC)
      DATA       CONVENTIONS/'GTOOL.1'/ !" <-- 1-th version
      INTEGER    LEN_UDUNI
      PARAMETER  (LEN_UDUNI=80)
      CHARACTER  UDUNIT_DATA *(LEN_UDUNI)
      CHARACTER  UDUNIT_TIME *(LEN_UDUNI)
      CHARACTER  UDUNIT_AX(3) *(LEN_UDUNI)
      CHARACTER  UDUNIT_AXW(3) *(LEN_UDUNI)
      CHARACTER  FIRST_DATE
      REAL       VALID_MIN,VALID_MAX
*
*     [WORK]
      INTEGER    IENUM
      INTEGER    IAX,STAT,NVDIMS,IEOD,LEN,I
*
*     [FUNCTION]
      INTEGER    LENZ
      LOGICAL    NH_INQ_DIMEXIST
      EXTERNAL   LENZ,NH_INQ_DIMEXIST
*
*     / initialize arrays /
*
      CALL GTSIZE ( HHEADA, IJKDIM )
*
*
*     / READ GTOOL HEADER AND GET AXIS INFO /
*
      CALL      GHUPAC
     I         ( HHEAD ,
     O           HDSET , HITEM , IFNUM , IDNUM ,
     O           HTITL , HUNIT ,
     O           ITIME , HDATE , HUTIM , ITDUR ,
     O           HAITM(1), IASTR(1), IAEND(1),
     O           HAITM(2), IASTR(2), IAEND(2),
     O           HAITM(3), IASTR(3), IAEND(3),
     O           HDFMT , VMISS ,
     O           DMIN  , DMAX  , DIVS  , DIVL  , ISTYP ,
     O           HCDATE, HCSIGN, HMDATE, HMSIGN         )
*
      CALL GHQENM(HHEAD,IENUM)
      CALL GHCGTS(HHEAD,'ETTL',HETTL)
*
      AXNUM=0
      IF(LENZ(HAITM(1)).NE.0) AXNUM=1
      IF(LENZ(HAITM(2)).NE.0) AXNUM=2
      IF(LENZ(HAITM(3)).NE.0) AXNUM=3
*
      DO IAX=1,AXNUM
        IANUM(IAX)=IAEND(IAX)-IASTR(IAX)+1
*
        CALL GUQAXV
     I       ( HHEAD,  IAX, 'WGT' ,
     O         HHEADA, AXW(1,IAX) , IEOD )      
        IF(IEOD.NE.0) CALL MSGDMP('E','INIT','Failed to read:'
     &                            //HAITM(IAX)//'(weight)')
        CALL GHCGET( HHEADA, 'UNIT',  HUNITAW(IAX)  )
*
        CALL GUQAXV
     I       ( HHEAD,  IAX, 'LOC' ,
     O         HHEADA, AX(1,IAX) , IEOD )      
        IF(IEOD.NE.0) CALL MSGDMP('E','INIT','Failed to read:'
     &                            //HAITM(IAX)//'(location)')
        CALL GHCGTS( HHEADA, 'TITL',  HTITLA(IAX)  )
        CALL GHCGET( HHEADA, 'DSET',  HDSETA(IAX)  ) !" cyclic if '^C'
        CALL GHCGET( HHEADA, 'UNIT',  HUNITA(IAX)  )
        CALL GHPGET( HHEADA, 'STYP',  STYPA(IAX)   )
      ENDDO
*
*     / NAME OF DIMMENSIONS /
*
      DO IAX=1,AXNUM
*       use title
        DIMNAME(IAX)=HTITLA(IAX)
*       replace ilegal/dangerous characters
        LEN=LENZ(DIMNAME(IAX))
        CALL SED_SG_NONREG(DIMNAME(IAX)(1:LEN),' ','_')
        CALL SED_SG_NONREG(DIMNAME(IAX)(1:LEN),'.','_')
        CALL SED_SG_NONREG(DIMNAME(IAX)(1:LEN),'@','A') !" some start with @
        IF(DIMNAME(IAX)(1:9).EQ.'longitude') THEN
          DIMNAME(IAX)='lon'  !" common name
        ELSE IF(DIMNAME(IAX)(1:8).EQ.'latitude') THEN
          DIMNAME(IAX)='lat'  !" common name
        ELSE IF(DIMNAME(IAX)(1:5).EQ.'sigma') THEN
          LEN=LENZ(HAITM(IAX))
          IF(HAITM(IAX)(LEN-1:LEN).EQ.'.M') THEN
            DIMNAME(IAX)='sigm'   !" intermediate sigma level
          ELSE
            DIMNAME(IAX)='sig'
          ENDIF
        ENDIF
      ENDDO
*
* ///////// NetCDF file initialization //////////////

*     / (later than first time:) redifine mode /

      IF(.NOT.FIRST) THEN
        STAT=NF_REDEF(NCID)
        IF(STAT.NE.NF_NOERR) CALL MSGDMP('E','INIT','NF_REDEF failed')
      ENDIF

*     / Are dims set ? /

      DO IAX=1,AXNUM
        NEWDIM(IAX)=.NOT. NH_INQ_DIMEXIST(NCID,DIMNAME(IAX))
        IF(.NOT.NEWDIM(IAX)) THEN
          STAT=NF_INQ_DIMID(NCID,DIMNAME(IAX),DIMID(IAX))
        ENDIF
      ENDDO

      IF(NH_INQ_DIMEXIST(NCID,'time')) THEN
        STAT=NF_INQ_DIMID(NCID,'time',DIMID(AXNUM+1))
      ENDIF

*     / SET DIMENSIONS /

      DO IAX=1,AXNUM
        IF(NEWDIM(IAX)) THEN
          STAT=NF_DEF_DIM(NCID,DIMNAME(IAX),IANUM(IAX),DIMID(IAX))
          IF(STAT.NE.NF_NOERR) THEN
            CALL MSGDMP('E','INIT',
     &             NF_STRERROR(STAT)//'set dimension '//DIMNAME(IAX))
          ENDIF
        ENDIF
      ENDDO
*
*     time
      IF (TIMESEQ.AND.FIRST) THEN
        STAT=NF_DEF_DIM(NCID,'time',NF_UNLIMITED,DIMID(AXNUM+1))
        IF(STAT.NE.NF_NOERR) CALL MSGDMP('E','INIT'
     &       ,NF_STRERROR(STAT)//'set dimension time')
      ENDIF
*
*     / DEFINE VARIABLES /
*
*     axes
      DO IAX=1,AXNUM
        IF(NEWDIM(IAX)) THEN
          STAT=NF_DEF_VAR(NCID,DIMNAME(IAX),NF_FLOAT,1,DIMID(IAX),
     &                  VARID_AX(IAX))
          IF(STAT.NE.NF_NOERR )CALL MSGDMP('E','INIT',NF_STRERROR(STAT))
        ENDIF
      ENDDO
*
*     time
      IF (TIMESEQ.AND.FIRST) THEN
        STAT=NF_DEF_VAR(NCID,'time',NF_FLOAT,1,DIMID(AXNUM+1),
     &                  VARID_TIME)
        IF(STAT.NE.NF_NOERR) CALL MSGDMP('E','INIT',NF_STRERROR(STAT))
      ENDIF
*
*     axes weight
      DO IAX=1,AXNUM
        IF(NEWDIM(IAX)) THEN
          LEN=LENZ(DIMNAME(IAX))
          DIMNAMEW(IAX)=DIMNAME(IAX)(1:LEN)//'wgt'
          STAT=NF_DEF_VAR(NCID,DIMNAMEW(IAX),NF_FLOAT,1,DIMID(IAX),
     &                  VARID_AXW(IAX))
          IF(STAT.NE.NF_NOERR )CALL MSGDMP('E','INIT',NF_STRERROR(STAT))
        ENDIF
      ENDDO
*
*     main variable
*
      DO IAX=1,AXNUM
        COUNT(IAX)=IANUM(IAX)
      ENDDO
      IF (TIMESEQ) THEN
        NVDIMS=AXNUM+1
        COUNT(NVDIMS)=1
      ELSE
        NVDIMS=AXNUM
      ENDIF
*
      IF(HDFMT.EQ.'UR8') THEN
        NCVAR_TYPE=NF_DOUBLE
      ELSE
        NCVAR_TYPE=NF_FLOAT
      ENDIF        
*
      CALL HEAD2VAR(HHEAD,VARNAME)
      STAT=NF_DEF_VAR(NCID,VARNAME,NCVAR_TYPE,NVDIMS,DIMID,VARID_DATA)
      IF(STAT.NE.NF_NOERR) CALL MSGDMP('E','INIT',NF_STRERROR(STAT))
*
*     / UNIT CONVERSION for udunits /
*
      CALL GT2NC_UNIT(HUNIT,'dummy','dum',UDUNIT_DATA)
      DO IAX=1,AXNUM
        IF(NEWDIM(IAX)) THEN
          CALL GT2NC_UNIT
     I                 (HUNITA(IAX),HAITM(IAX),'LOC',
     O                  UDUNIT_AX(IAX))
          CALL GT2NC_UNIT
     I                 (HUNITAW(IAX),HAITM(IAX),'WGT',
     O                  UDUNIT_AXW(IAX))
        ENDIF
      ENDDO
*
      CALL DATE_TIME(HDATE,ITIME,HUTIM)  !" HDATE=HDATE-ITIME
*
      CALL GT2NC_UNIT(HUTIM,'dummy','dum',UDUNIT_TIME)
      LEN=LENZ(UDUNIT_TIME)
      UDUNIT_TIME(LEN+2:LEN+7)='since '
      CALL SED_SG_NONREG(HDATE(1:4),' ','0')
      UDUNIT_TIME(LEN+9:LEN+12)=HDATE(1:4)   !" yyyy
      UDUNIT_TIME(LEN+13:LEN+13)='-'
      UDUNIT_TIME(LEN+14:LEN+15)=HDATE(5:6)  !" mm
      UDUNIT_TIME(LEN+16:LEN+16)='-'
      UDUNIT_TIME(LEN+17:LEN+18)=HDATE(7:8)  !" dd
      UDUNIT_TIME(LEN+19:LEN+19)=' '
      UDUNIT_TIME(LEN+20:LEN+21)=HDATE(9:10)  !" hh
      UDUNIT_TIME(LEN+22:LEN+22)=':'
      UDUNIT_TIME(LEN+23:LEN+24)=HDATE(11:12) !" mm
      UDUNIT_TIME(LEN+25:LEN+25)=':'
      UDUNIT_TIME(LEN+26:LEN+27)=HDATE(13:14) !" ss
      UDUNIT_TIME(LEN+28:LEN+33)=' +0:00'     !" assume UT
*
*     / ATTRIBUTES /
*
*     global
*
      IF (FIRST) THEN
        STAT=NF_PUT_ATT_TEXT(NCID,NF_GLOBAL,'Conventions',
     &                     LENZ(CONVENTIONS),CONVENTIONS)
        IF(STAT.NE.NF_NOERR) CALL MSGDMP('E','INIT',NF_STRERROR(STAT))
*
        STAT=NF_PUT_ATT_TEXT(NCID,NF_GLOBAL,'title',
     &                     LENZ(HDSET),HDSET)
        IF(STAT.NE.NF_NOERR) CALL MSGDMP('E','INIT',NF_STRERROR(STAT))
*
        STAT=NF_PUT_ATT_TEXT(NCID,NF_GLOBAL,'history',(NCC)*4+34,
     &               'CDATE: '//HCDATE//', CSIGN: '//HCSIGN//
     &                ', MDATE: '//HMDATE//', MSIGN: '//HMSIGN)
        IF(STAT.NE.NF_NOERR) CALL MSGDMP('E','INIT',NF_STRERROR(STAT))
      ENDIF
*
*     main data
*
      STAT=NF_PUT_ATT_TEXT(NCID,VARID_DATA,'long_name',
     &                     LENZ(HTITL),HTITL)
      STAT=NF_PUT_ATT_TEXT(NCID,VARID_DATA,'units',
     &                      LENZ(UDUNIT_DATA),UDUNIT_DATA)

**      MISS_VAL=-NF_FILL_FLOAT !"\sim 1e37(ver3): normaly outside valid ranges
      MISS_VAL=-9.9D36 !"＜NF_FILL_FLOATはやたら桁数の多い数なので気持良くない
      VALID_MIN=-MIN( 0.8*ABS(NF_FILL_FLOAT), 0.9*ABS(MISS_VAL) )
      VALID_MAX= MIN( 0.8*ABS(NF_FILL_FLOAT), 0.9*ABS(MISS_VAL) )

      STAT=NF_PUT_ATT_DOUBLE(NCID,VARID_DATA,'missing_value',NF_DOUBLE,
     &                       1,MISS_VAL)
      STAT=NF_PUT_ATT_REAL(NCID,VARID_DATA,'valid_min',NF_DOUBLE,
     &                       1,VALID_MIN)
      STAT=NF_PUT_ATT_REAL(NCID,VARID_DATA,'valid_max',NF_DOUBLE,
     &                       1,VALID_MAX)

      STAT=NF_PUT_ATT_TEXT(NCID,VARID_DATA,'history',NCC*IENUM,HETTL)
*
*     axes
*
      DO IAX=1,AXNUM
        IF(NEWDIM(IAX)) THEN
          STAT=NF_PUT_ATT_TEXT(NCID,VARID_AX(IAX),'long_name',
     &                         LENZ(HTITLA(IAX)),HTITLA(IAX))
          STAT=NF_PUT_ATT_TEXT(NCID,VARID_AX(IAX),'units',
     &                         LENZ(UDUNIT_AX(IAX)),UDUNIT_AX(IAX))
**          IF(HDSETA(IAX)(1:1).EQ.'C') THEN
**            STAT=NF_PUT_ATT_TEXT(NCID,VARID_AX(IAX),'cyclic',3,'yes')
**          ELSE
**            STAT=NF_PUT_ATT_TEXT(NCID,VARID_AX(IAX),'cyclic',3,'no ')
**          ENDIF
          IF(STYPA(IAX).GE.0) THEN
            STAT=NF_PUT_ATT_TEXT(NCID,VARID_AX(IAX),'positive',
     &                           2,'up')
          ELSE
            STAT=NF_PUT_ATT_TEXT(NCID,VARID_AX(IAX),'positive',
     &                           4,'down')
          ENDIF
          IF(DIMNAME(IAX).EQ.'lon') THEN
            STAT=NF_PUT_ATT_TEXT(NCID,VARID_AX(IAX),'quantity',
     &                           9,'longitude')
            STAT=NF_PUT_ATT_REAL(NCID,VARID_AX(IAX),'modulo',
     &                           NF_DOUBLE,1,360.)
          ELSE IF(DIMNAME(IAX).EQ.'lat') THEN
            STAT=NF_PUT_ATT_TEXT(NCID,VARID_AX(IAX),'quantity',
     &                           8,'latitude')
          ELSE IF(DIMNAME(IAX)(1:3).EQ.'sig') THEN
            STAT=NF_PUT_ATT_TEXT(NCID,VARID_AX(IAX),'quantity',
     &                           6,'sigma')
          ELSE
            STAT=NF_PUT_ATT_TEXT(NCID,VARID_AX(IAX),'quantity',
     &                           LENZ(DIMNAME(IAX)),DIMNAME(IAX))
          ENDIF
          STAT=NF_PUT_ATT_TEXT(NCID,VARID_AX(IAX),'weight',
     &                           LENZ(DIMNAMEW(IAX)),DIMNAMEW(IAX))
        ENDIF
      ENDDO
*
*     axes weight
*
      DO IAX=1,AXNUM
        IF(NEWDIM(IAX)) THEN
          TITLEW='weight of '//HTITLA(IAX)
          STAT=NF_PUT_ATT_TEXT(NCID,VARID_AXW(IAX),'long_name',
     &                         LENZ(TITLEW),TITLEW)
          STAT=NF_PUT_ATT_TEXT(NCID,VARID_AXW(IAX),'units',
     &                         LENZ(UDUNIT_AXW(IAX)),UDUNIT_AXW(IAX))
        ENDIF
      ENDDO
*
*     time axis
      IF (TIMESEQ.AND.FIRST) THEN
        STAT=NF_PUT_ATT_TEXT(NCID,VARID_TIME,'long_name',
     &                      4,'time')
        STAT=NF_PUT_ATT_TEXT(NCID,VARID_TIME,'units',
     &                      LENZ(UDUNIT_TIME),UDUNIT_TIME)
        STAT=NF_PUT_ATT_TEXT(NCID,VARID_TIME,'positive',
     &                       2,'up')
        STAT=NF_PUT_ATT_TEXT(NCID,VARID_TIME,'quantity',
     &                      4,'time')
      ENDIF
*
*     / quit define mode /
*   
      STAT=NF_ENDDEF(NCID)
      IF(STAT.NE.NF_NOERR) CALL MSGDMP('E','INIT',NF_STRERROR(STAT))
*
*     / write dimension variables exept for time /
*
      DO IAX=1,AXNUM
        IF(NEWDIM(IAX)) THEN
          STAT=NF_PUT_VAR_REAL(NCID,VARID_AX(IAX),AX(1,IAX))
          IF(STAT.NE.NF_NOERR) CALL MSGDMP('E','INIT',NF_STRERROR(STAT))
          STAT=NF_PUT_VAR_REAL(NCID,VARID_AXW(IAX),AXW(1,IAX))
          IF(STAT.NE.NF_NOERR) CALL MSGDMP('E','INIT',NF_STRERROR(STAT))
        ENDIF
      ENDDO
*
      FIRST=.FALSE.
      RETURN
      END
*********************************************************************
      SUBROUTINE GT2NC_UNIT
     I     (UNIT_GT,AX_ITEM,LOC_WGT,
     O      UNIT_NC)
*
      include   'netcdf.inc'
**#include        "udunits.inc"
*
      CHARACTER UNIT_GT*(*)
      CHARACTER AX_ITEM*(*)
      CHARACTER LOC_WGT*(3) !" 'LOC' or 'WGT' or 'DUM'etc(=dummy)
      CHARACTER UNIT_NC*(*)
      LOGICAL   LCHREQ !" dcl func
      INTEGER   LENZ
      INTEGER   utInit,utScan
      INTEGER   STAT
**      PTR   UDUNI_PTR  !" PTR is defined as INTEGER in udunits.inc
*
*     longitude:
      IF( LCHREQ(AX_ITEM(2:4),'lon') .OR. LCHREQ(AX_ITEM(1:3),'lon')
     &                                                         ) THEN
        IF( LCHREQ(UNIT_GT,'deg') ) THEN
          IF(LOC_WGT.EQ.'LOC') THEN
            UNIT_NC='degree_east'
          ELSE
            UNIT_NC='angular_degree'
          ENDIF
          RETURN
        ELSE
          UNIT_NC=UNIT_GT
          CALL MSGDMP('W','UNIT',' NOT CONVERTED TO UNIDATA FORM')
          write (6,*) AX_ITEM, ' unit:', UNIT_GT
          RETURN
        ENDIF
*     latitude:
      ELSE IF(LCHREQ(AX_ITEM(2:4),'gla') .OR. 
     &    LCHREQ(AX_ITEM(2:4),'lat') .OR. LCHREQ(AX_ITEM(1:3),'lat')
     &                                                         ) THEN
        IF( LCHREQ(UNIT_GT,'deg') ) THEN
          IF(LOC_WGT.EQ.'LOC') THEN
            UNIT_NC='degree_north'
          ELSE
            UNIT_NC='angular_degree'
          ENDIF
          RETURN
        ELSE
          UNIT_NC=UNIT_GT
          CALL MSGDMP('W','UNIT',' NOT CONVERTED TO UNIDATA FORM')
          write ( 6,*) AX_ITEM, ' unit:', UNIT_GT
          RETURN
        ENDIF
*     general:
      ELSE
        IF (LENZ(UNIT_GT).EQ.0) THEN      !" nondimensional
          UNIT_NC=''
**          UNIT_NC='count'   !" should be used for integer only
          RETURN        
        ELSE
          UNIT_NC=UNIT_GT
          CALL SED_SG_NONREG(UNIT_NC,'MIN','min') !" udunits is case sensitive
          CALL SED_SG_NONREG(UNIT_NC,'SEC','sec')
          CALL SED_SG_NONREG(UNIT_NC,'HOUR','hour')
          CALL SED_SG_NONREG(UNIT_NC,'DAY','day')
          CALL SED_SG_NONREG(UNIT_NC,'mb','hPa') !"mb is not supported
          CALL SED_SG_NONREG(UNIT_NC,'g/g','') !" for specific humidity
          CALL SED_SG_NONREG(UNIT_NC,'**','^') !" though not necessary
          CALL SED_SG_NONREG(UNIT_NC,'|1"','^1') !" convert dcl style
          CALL SED_SG_NONREG(UNIT_NC,'|2"','^2')
          CALL SED_SG_NONREG(UNIT_NC,'|3"','^3')
          CALL SED_SG_NONREG(UNIT_NC,'|4"','^4')
          CALL SED_SG_NONREG(UNIT_NC,'|5"','^5')
          CALL SED_SG_NONREG(UNIT_NC,'|6"','^6')
          CALL SED_SG_NONREG(UNIT_NC,'|7"','^7')
          CALL SED_SG_NONREG(UNIT_NC,'|8"','^8')
          CALL SED_SG_NONREG(UNIT_NC,'|9"','^9')
          CALL SED_SG_NONREG(UNIT_NC,'|-1"','^-1')
          CALL SED_SG_NONREG(UNIT_NC,'|-2"','^-2')
          CALL SED_SG_NONREG(UNIT_NC,'|-3"','^-3')
          CALL SED_SG_NONREG(UNIT_NC,'|-4"','^-4')
          CALL SED_SG_NONREG(UNIT_NC,'|-5"','^-5')
          CALL SED_SG_NONREG(UNIT_NC,'|-6"','^-6')
          CALL SED_SG_NONREG(UNIT_NC,'|-7"','^-7')
          CALL SED_SG_NONREG(UNIT_NC,'|-8"','^-8')
          CALL SED_SG_NONREG(UNIT_NC,'|-9"','^-9')
        ENDIF
      ENDIF
*
***"    以下何故かうまく行かないので、udunitインクルードはとりあえずやめ.
**      STAT=UTOPEN()
**      STAT=utdec(UNIT_NC,UDUNI_PTR) !" 変換がうまく行けば認められる
**      IF(STAT.NE.0) write(6,*) '** WARNING (GT2NC_UNIT) ** ',
**     &                'DERIVED UNIT IS NOT VALID FOR UDUNITS LIBRARY'
*
      RETURN
      END
*********************************************************************
      SUBROUTINE SED_SG_NONREG(INOUT,BEFORE,AFTER)
      CHARACTER INOUT*(*),BEFORE*(*),AFTER*(*)
      CHARACTER WORK*(160)
      INTEGER   NB,NA,NIN,NLEN,I,J
      INTEGER   LENZ
*
*     CAUTION: I CAN'T HANDLE REGULAR EXPRESIONS, DO JUST NAIVE REPLACEMENT
*
      NB=LEN(BEFORE)
      NA=LEN(AFTER)
      NIN=LENZ(INOUT)
      NLEN=LEN(INOUT)
*
      IF(NIN.GT.160) CALL MSGDMP('W','SED_SG','INSUFFICIENT WORK SPACE')
      WORK=INOUT
*
      INOUT=''  !" initialization
*
      J=1
      I=1
 1100 CONTINUE
        IF(WORK(I:MIN(I+NB-1,NIN)).NE.BEFORE(1:NB)) THEN
          INOUT(J:J)=WORK(I:I)
          J=J+1
          I=I+1
        ELSE
          IF(NA.NE.0) INOUT(J:J+NA-1)=AFTER(1:NA)
          J=J+NA
          I=I+NB
        ENDIF
        IF(J.GT.NLEN+1) THEN
          CALL MSGDMP('W','SED_SG','INSUFFICIENT OUT SPACE')
          I=NIN+1 ! quit the loop
        ENDIF
      IF(I.LE.NIN) GOTO 1100
*
      END
***********************************************************************
*     DATE=DATE-TIME      
*-------------------
      SUBROUTINE DATE_TIME(DATE,TIME,TIMEUNIT)
      CHARACTER DATE*14 ! yyyymmddhhmmss
      INTEGER TIME
      CHARACTER TIMEUNIT*(*)

      INTEGER DATE1,HW,MW,SW, DAY2, DATEO
      INTEGER SDMOD1,SDMOD2,SDMODO
      LOGICAL LCHREQ
      EXTERNAL LCHREQ

      READ(DATE(1:8),'(I)') DATE1
      READ(DATE(9:10),'(I)') HW
      READ(DATE(11:12),'(I)') MW
      READ(DATE(13:14),'(I)') SW
      SDMOD1=SW*60*MW+3600*HW

      IF( LCHREQ(TIMEUNIT(1:3),'DAY') ) THEN
        DAY2=TIME
        SDMOD2=0
      ELSE IF( LCHREQ(TIMEUNIT(1:4),'HOUR') ) THEN
        DAY2=TIME/24
        SDMOD2=3600*(TIME-24*DAY2)
      ELSE IF( LCHREQ(TIMEUNIT(1:3),'MIN') ) THEN
        DAY2=TIME/1440
        SDMOD2=60*(TIME-1440*DAY2)
      ELSE IF( LCHREQ(TIMEUNIT(1:3),'SEC') ) THEN
        DAY2=TIME/86400
        SDMOD2=TIME-86400*DAY2
      ELSE
        CALL MSGDMP('E','DATE_T','INVALID UNIT')
      ENDIF

      CALL DATEF1(-DAY2,DATE1,DATEO)  !" DATEO=DATE1-DAY2
      SDMODO=SDMOD1-SDMOD2

      IF(SDMODO.LT.0) THEN
        DATE1=DATEO !" use as work area
        CALL DATEF1(-1,DATE1,DATEO)
        SDMODO=SDMODO+86400
      ENDIF

      WRITE(DATE(1:8),'(I8)') DATEO
      HW=SDMODO/3600
      MW=(SDMODO-3600*HW)/60
      SW=SDMODO-60*MW-3600*HW

      WRITE(DATE(9:10),'(I2)') HW
      WRITE(DATE(11:12),'(I2)') MW
      WRITE(DATE(13:14),'(I2)') SW

      CALL SED_SG_NONREG(DATE(1:14),' ','0')

      RETURN
      END

