*" 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